Two-dimensional particle-in-cell simulation of the expansion of a plasma into a rarefied medium

The expansion of a dense plasma through a more rarefied ionized medium has been studied by means of two-dimensional particle-in-cell simulations. The initial conditions involve a density jump by a factor of 100, located in the middle of an otherwise equally dense electron–proton plasma with uniform proton and electron temperatures of 10 eV and 1 keV, respectively. Simulations show the creation of a purely electrostatic collisionless shock together with an ion-acoustic soliton tied to its downstream region. The shock front is seen to evolve in filamentary structures consistently with the onset of the ion–ion instability. Meanwhile, an un-magnetized drift instability is triggered in the core part of the dense plasma. Such results explain recent experimental laser–plasma experiments, carried out in similar conditions, and are of intrinsic relevance to non-relativistic shock scenarios in the solar and astrophysical systems.

4 a shock. The rapid formation time and propagation speed of such a shock in the simulation matched that of its experimental counterpart and the electrostatic field amplitudes were also in good agreement. The rarefaction wave was destabilized in the simulation by an electrostatic drift instability [26,27], which triggered the formation of plasma structures that resembled the ion-acoustic waves observed in the experiment.
This good agreement between the experimental and simulation results motivated us to further study the thermal expansion of a plasma with the initial conditions used in [22], but with the insertion of a transverse dimension. This is possible by the rapid formation of the shock observed in [4,22]. Our main results are as follows. We demonstrate the creation of an ion-acoustic soliton, trailing the main collisionless shock, which traps a population of protons counter-streaming from the dilute plasma. Such a structure presents an electric field distribution that closely resembles the one experimentally detected in [4]. Such a soliton was not triggered in the 1D case, possibly due to a poorer plasma representation. We also show that a protonbeam-driven instability is triggered in a hot thermal electron population [28], resulting in an oblique structuring of the plasma protons on a Debye length scale just ahead of the shock. This instability is fed by electric fields rather than by magnetic fields and is thus compatible with the experimental observations in [4].
The absence of magnetic fields, which appears to be in contradiction to the results reported in [18], can be explained by considering that the electrons of the ambient plasma, which were not considered in [18], are accelerated by the rarefaction wave and form a super-thermal population behind the expanding ion front. Super-thermal electron populations can stabilize the electron distribution against the Weibel instability [29] and may explain why we could not observe it here. Moreover, simulations show a very modest cooling of the electrons along the plasma propagation direction, another factor against the onset of the Weibel instability.
It is worth mentioning here that no significant magnetic fields were detected in the 1D case [22] either. However, a 1D simulation allows only for a magnetic field growth perpendicular to the simulation axis. The Weibel instability requires a wavevector component orthogonal to the expansion direction and could therefore have been artificially suppressed in a 1D geometry. A 2D setup, instead, allows for the growth of in-plane magnetic fields and more complex current systems that can also lead to more realistic out-of-plane magnetic fields. The lack of magnetic fields in the present simulation results is therefore a more reliable indication of the suppression of the Weibel instability in this plasma scenario.
The structure of this paper can be summarized as follows: the PIC method and the initial conditions will be discussed in section 2 and the simulation results will be presented in section 3. Section 4 will instead be entirely devoted to the filamentary structures at the shock front and a discussion of the main reported results will be given in section 5. Finally, a conclusive summary will be presented in section 6.

The particle-in-cell (PIC) code and the initial conditions
A PIC simulation code solves the Maxwell's equations on a grid together with the relativistic Lorentz force equation for an ensemble of computational particles. Each computational particle (CP) can have a charge and a mass that are different from those of a plasma particle, but the ratios are the same. The ensemble properties of CPs are a good approximation to those of a collisionless plasma, provided that an adequate number of CPs are used and that the resolution of the simulation grid is high enough. The currents and the charges of individual CPs are 5 interpolated to the simulation grid, and a global current J(x) and charge ρ(x) are obtained from a summation over all CPs. These global distributions are used to update the electromagnetic fields E(x) and B(x) with the help of Maxwell's equations. The updated electromagnetic fields are then interpolated to the position x j of each CP with index j and its momentum p j is updated with the help of the relativistic Lorentz force equation. Its position is advanced with the product of the simulation time step t and the CP's velocity v j . These steps are repeated n times, until the interval N t of interest is covered. Maxwell's equations and the Lorentz force equation for the jth CP are The numerical scheme our code uses is described in detail in [30]. For a more general description of the PIC method, see [31]. We take initial conditions for the plasma which are identical to those in [22], but shorten the simulation box along the plasma expansion direction and introduce instead a second spatial dimension. The box length along the plasma expansion direction x is L x and the box length in the y-direction is L y . We subdivide the box along x into two halves. The interval −L x /2 < x < 0 contains the dense plasma cloud and the interval 0 < x < L x /2 the tenuous one. The electron and proton number densities of the tenuous cloud are each n 0 = 10 15 cm −3 . The electron and proton number densities of the dense cloud are both 100n 0 . The electron (proton) temperature of both clouds is 1 keV (10 eV).
The simulation domain covers the spatial interval L x × L y = 540λ D × 60λ D , where λ D = v e /ω pt is the Debye length of the dilute plasma. An electron thermal speed of v e ≈ 1.3 × 10 7 m s −1 and the plasma frequency of the diluted plasma ω pt = (n 0 e 2 /m e 0 ) 1/2 = 1.8 × 10 12 s −1 are realistic for the ambient plasma in laser-plasma experiments. The dimension of the interval L x /2 × L y , which is initially occupied by the dilute plasma, is then ≈ 2 mm × 0.44 mm. It is comparable in size to that displayed in figure 1 of [4]. Our simulation time N t ω pt = 775 (N = 2.9 × 10 5 ) corresponds to approximately half a nanosecond, which is comparable to the time at which the snapshot in figure 1 of [22] was taken. Our simulation should, in principle, capture the source process which is resulting in the structures visible in this picture.
The full simulation domain is resolved by a grid with 4000 × 440 cells and all boundary conditions are periodic. The periodic boundary conditions imply that we have two rarefaction waves: one at x ≈ 0 and a second practically identical one at the boundary between x ≈ L x /2 and x ≈ −L x /2. We will focus our attention only on the central one. We use 360 CPs per simulation cell to represent the electrons of the dense cloud, and we use the same number for the protons. The electrons and protons of the tenuous cloud are represented each by 160 particles per cell. The CPs of the dense cloud thus have a higher numerical weight in that their charge and mass are higher than those of the CPs of the tenuous cloud by a factor of ≈ 45.
It is worth noting here that, in the dense plasma region, the electron Debye length is resolved by 0.75 grid cell. Generally speaking, such an under-resolved electron population may lead to numerical instabilities via the interaction of alias (i.e. unphysical) phase velocities with the particle distribution [32]. An artificial heating of the plasma might occur [33], whose most relevant signatures are an increment of the total internal energy and a thermal broadening of the particles phase space. However, neither of these effects is seen in our simulation. By the end of the simulation, the internal energy has increased by only a factor of 1.0005 (see figure 1), which can be considered as entirely negligible in simulating the plasma dynamics. Moreover, the electron phase space does not show any sensible thermal broadening (see movie 1, available from stacks.iop.org/NJP/13/073023/mmedia). This can be understood by considering that the alias phase velocity can interact with either the flat part (implying that the instability is confined to long wavelengths and is very weak) or the steep part of the particle distribution function (resulting in strong artificial instabilities). In the case of linear weighting, such as the one used in our simulation, significant self-heating should occur only at λ D / x 1/π ≈ 0.3 [32] ( x being the grid spacing). For our simulation parameters, λ D / x ≈ 0.75, justifying why no significant self-heating is observed. It is therefore reasonable to assume that the plasma dynamics observed in the simulation is exclusively driven by realistic, physical processes.

Simulation results
In this section we will discuss the simulated plasma dynamics following the evolution in time of the electron and proton phase spaces and the related electric field distributions. It has to be noted that the gross behavior of the system is seen to be predominantly confined to the x-direction (i.e. parallel to the density gradient induced as the initial condition); therefore, unless specifically stated, all the physical quantities shown in the following plots will be the result of an average along the y-direction. This averaging operation offers the advantage of notably increasing the phase-space resolution of the simulation results. Compared to the 1D case [22], the average over the y-direction induces a statistical improvement of a factor 1000. Movie 1 (available from stacks.iop.org/NJP/13/073023/mmedia) shows the temporal evolution of the electron (a) and proton (b) phase spaces, averaged over the y-coordinate, along the x-coordinate for a temporal window of 0 < t < 900. In the following, we select several snapshots of this movie and discuss them in detail. Unless specifically stated, throughout the rest of the paper the spatial and temporal coordinates will be normalized in terms of the Debye length (≈ 7 µm) and the inverse of the plasma frequency (≈ 0.5 ps) of the rarefied plasma. The densities will be normalized in terms of the initial electron density of the rarefied plasma as well. On the other hand, in order to make our results more easily comparable with experiments, the electric and magnetic fields will always be presented in physical units. However, due to the normalization system used, the normalized electric field can be easily obtained by dividing the electric field in physical units by (m e λ D ω 2 pt )/e ≈ 10 8 . As discussed in the introduction, the initial conditions of our simulation involve a sharp density jump located at x = 0. The electrons of the higher-density plasma will immediately start to diffuse, guided by their own pressure, into the more tenuous plasma, roughly preserving, at least in proximity to the boundary, a Maxwellian distribution (figure 2(a)). The more massive protons will not be able to keep up with this motion and a net charge imbalance will be set close to the plasmas interface (located approximately at 0 x 2 (figure 2(c)). This imbalance sets a net positive electrostatic field (figure 2(d)) which is almost uniform in the range 0 < x < 1 and presents a peak at x ≈ 2, in correspondence to the front of the expanding plasma ( figure 2(b)). This field accelerates the protons of both plasmas inducing a bending in the proton phase space, the first signature of the onset of a shock, and a downstream rarefaction wave. This initial evolution of the plasma closely resembles the one obtained in the 1D case (see [22]) and this can be explained by taking into account the planar geometry of the plasmas boundary. As we proceed in time, the electrons of the dense plasma keep on streaming towards the dilute plasma, and the onset of a return current of electrons, counter-streaming from the dilute to the dense plasma, starts to become visible (the narrow beam at the bottom of the phase space in figure 3(a)). At this time (t 2 = 150, corresponding in physical units to ≈ 88 ps) the rarefaction wave has almost completely settled, ranging from x ≈ −4 to x ≈ 15, disrupted in the middle by a bending in the phase space, a signature of the onset of a shock (figure 3(b)). Such a structure is associated with a positive hump in the electric field distribution (x ≈ 6 in figure 3(d)), the result of a charge imbalance (see figure 3(c)). A secondary peak in the electric field distribution is present at x ≈ 14 and it is associated with a charge imbalance at the tip of the forward accelerated proton beam which, at this time, as already reached a peak velocity of the order of approximately 60 times the initial proton thermal speed (in agreement with the 1D results reported in [22]). The downstream part of the rarefaction wave is instead associated with an almost uniform electric field distribution (−4 x 4 in figure 3(d)). A tiny downstream region separates this part of the rarefaction wave from the shock (4 x 6 in figure 3(b)). In this region, two different proton populations can be distinguished in the proton phase space. The first has a higher density and a higher velocity and is associated with the proper downstream part of the rarefaction wave, created by the protons of the dense plasma. The second is slower and less dense and is associated with the protons of the dilute plasma. These two proton populations therefore locally produce a density plateau forcing the electrostatic field to locally tend to zero (see x ≈ 5 in figure 3(d)). This tiny structure was also visible in the 1D simulations [22]. It is interesting to note that many more protons are present in this low density region compared to the 1D simulation results (see, e.g., figure 4(b) in [22]). Such a difference can be explained just in terms of dimensionality effects: a purely collisionless and un-magnetized plasma, such as the one modeled here, does not allow crossing of the phase-space trajectories. In the 1D case, this implies that trajectories are always separated in the v x space. If an additional dimension is inserted, the particles acquire a further degree of freedom. This implies that, in a reduced phase space resulting from integration along y, the particles are allowed to mix. Again, the integration of the additional dimension y allows a much larger dynamical range for the densities along x, due to the better statistical representation. This is true thanks to the weak dependence of the shock upon the transverse dimension (quasi-planarity). At later times (t 3 = 300 or, in physical units, t 3 ≈ 176 ps) the shock has moved up to x ≈ 12 (figure 4(b)), meaning an average speed of 0.04 times the electron thermal velocity. This velocity lies in the ion-acoustic range. Generally speaking, the ratio between the ion-acoustic speed (c s ) and the electron thermal velocity (v e ) in a Maxwellian plasma is c s /v e = (γ Z m e /m i ) 1/2 . In our simulation m e /m i = 1/1836, Z = 1 and the adiabatic index γ = (N + 2)/N = 2. This is concluded on the basis that each particle has only two translational degrees of freedom (N ). The ion-acoustic speed is then 0.03 times the electron thermal speed. This structure is associated with a narrow, monopolar, electrostatic peak with an amplitude of 6 × 10 7 V m −1 and a full-width at half-maximum of the order of the Debye length of the dilute plasma. The spatial overlapping of the slow and fast protons in the downstream region is still visible, located at 8 x 12 ( figure 4(b)). In this region, in fact, the fast protons form, consistent with the rarefaction wave onset, a narrow beam with an average speed of the order of 20, whereas the colder protons have an average speed centered at approximately 10. It is this latter proton population that induces a kink in the electric field distribution in the wake of the shock (figure 4(d)), as discussed above. This structure is a distinct result of the presence of a background medium as it is clear by separating the contributions of the protons of the two plasmas in the proton density distribution of figure 4(c) (see figure 5). In this graph the overall proton density distribution is plotted (black line) together with the dense and dilute plasma ones. The protons of the dense expanding plasma (blue line) are seen to decrease their density, in an exponential manner (dashed line), from the core to x ≈ 10. This implies a smooth decrease of the density without the appearance of any shock-like structure. This is the expected behavior of the proton density distribution in the case of a freely expanding plasma into a vacuum (see, e.g., for instance [12]). However, when the density of the denser plasma becomes comparable to the rarefied one (i.e. in the region around x ≈ 12), a sharp density decrease is seen, which is associated with the shock. On the other hand, the proton density of the rarefied plasma is seen to be roughly constant for x > 15 and piles up only in correspondence to the shock. This latter proton density is then negligible before the shock, thus not perturbing the rarefaction wave of the dense plasma. The presence of a rarefied, ambient plasma therefore significantly modifies the idealized behavior of a freely expanding plasma in mainly two ways. The first apparent one is that the proton density distribution of the expanding plasma is forced to converge down to the density of the rarefied plasma. The second, less intuitive one is that a population of protons, streaming from the rarefied plasma, peaks up in correspondence to the downstream region of the shock (peak of the red curve located approximately at x = 10 in figure 5). It is this proton population that can also be easily spotted in the proton phase-space distribution (figure 4(b)), which contributes to the shock creation and, as will be seen in the following, is responsible for creating an ion-acoustic soliton in the downstream region of the shock. Meanwhile, nothing relevant appears to happen in the electron phase space except for a density jump in correspondence to the shock and a broadening of the return current of cold electrons. This is an expected consequence of the much higher mobility of the electrons; except for the first preliminary stage of the simulation (compare figure 2), in which the electrons trigger the dense plasma expansion, the electrons are too fast to suitably respond to the ion motion and can thus be treated simply as a neutralizing negative fluid.
Later in time (t 4 = 600 or, in physical units, 350 ps), the shock has moved up to x ≈ 26 (figure 6(b)), thus preserving an average propagation velocity in the ion-acoustic range. The associated electrostatic peak has decreased its amplitude down to ≈ 4 × 10 7 V m −1 and broadened its width in order to keep approximately constant the associated potential (figure 6(c)). The low-energy proton population at the downstream part of the shock has further developed into a soliton-like structure with an associated bipolar electric field (x ≈ 20 in figure 6(c)). This structure, which is clearly related to the mixing in space but not in velocity of the two plasmas, could not be detected in the 1D case where only a few slow protons were present at the downstream part of the shock [22]. It is interesting to note that the y-component of the electric field distribution (figure 6(d)) does not show any significant peak in correspondence to both the shock (x ≈ 26) and the soliton (x ≈ 20), suggesting the 1D structure of their electric field distributions. A broad peak in the electric field is in fact present only in the 'foreshock' region (i.e. that spatial interval, just in front of the shock, in which the ambient plasma is perturbed by the presence of the shock, approximately at x = 28); such an electric field is related, as will be extensively discussed in the following section, to the onset of an instability, which is transverse to the shock front, driven by the two counter-propagating proton beams. The absence of such E y components, in correspondence to the shock and the soliton, is also a good indication that the plasma flow direction (i.e. the direction along which the strongest currents are expected to be) is mostly parallel to E x and therefore those structures are purely electrostatic.
Moreover, in the very high density part of the rarefaction wave, an electrostatic ionic instability is set (see the region −10 x −7 in figures 6(b) and (c)). Such instability induces a strong modulation in the proton phase space associated with a bipolar electrostatic field with a maximum amplitude of the order of 4 × 10 7 V m −1 . The long wavelength at which this instability is seen to occur is a good indication of it being related to the so-called drift instability [26,27]. The relatively small drift velocity (v drift 0.9v te ) permits us in fact to rule out the onset of the Buneman instability [34]. The drift instability occurs whenever a gradient in a particular plasma parameter (say temperature, density or magnetic fields) is set on a scale that is comparable to the related plasma scale length [26]. Nevertheless, the region at which this instability is seen to be triggered shows neither significant temperature gradients nor magnetic fields (see figure 14 and the following discussion) and we can thus assume that only density gradients are responsible for this. A small fluctuation in proton density, due to an oscillation in proton drift speed, generates in fact a modulation of the electric field, which further amplifies the density gradient setting an instability loop. Generally speaking, such fluctuations can occur throughout the entire extension of the rarefaction wave but the instability is seen to occur just at the foot of it. This can be explained by considering that the growth rates of instabilities in un-magnetized plasmas scale with the plasma frequency of the relevant species or, equivalently, to the square root of the species density. As this instability is purely ionic, the maximum growth rate is expected to be seen where the proton density is maximum, i.e. at the foot of the rarefaction wave (see, e.g., figure 6(a)). Therefore, even if the instability can, in general, be triggered throughout the entire rarefaction wave, we see it first at the foot of the wave, because it is here that the growth rate is maximized.
Detailed theoretical works on this instability usually deal with the presence of a magnetic field (see, e.g., [27]), but, to the knowledge of the present authors, a clear treatment in the case of only electrostatic phenomena has not been reported in the literature to date and will be the subject of a further work. It is worth mentioning here that, at the foot of the rarefaction wave, the grid size is smaller than the electron Debye length (see section 2); this instability could thus be, in principle, just a consequence of the onset of numerical instability brought about by the artificial grid heating that this implies. Even if the simulation results provide several indications that significant self-heating should not occur (see section 2), the grids could, in principle, modify the fluctuation spectrum through aliasing. Analytical expressions of fluctuation spectra in plasma gradients are extremely difficult to derive, even in the absence of aliasing [35]. They are thus beyond the scope of the present paper and will be the subject of a future work.
Finally, the plasma state at t 5 = 900 is displayed in figure 7. The main shock structure does not show any significant difference in its velocity, while the amplitude of its electric field has decreased and the width has broadened in order to keep the associated potential roughly constant. The y-component of the electric field still presents a pronounced peak just in the foreshock region. Again, this is related to an instability set by the two counter-propagating proton beams, which will be discussed in detail in the next section. On the other hand, the population of slow protons trailing the shock has finally developed a solitary structure in the phase space. In its comoving reference frame, this density hole presents a width in space on the Debye length scale and a width in velocity of the order of twice the ion-acoustic speed of the dilute plasma. This latter structure, once formed, is seen to propagate, almost unperturbed, closely tied to the shock. A zoom-in of such a phase-space region is shown in figure 8(a). The related electric field presents a bipolar distribution with a maximum amplitude of ≈ 10 7 V m −1 with which is associated a monopolar positive potential profile ( figure 8(b)) with a maximum amplitude of approximately 70 V. All these characteristics corroborate the thesis that this structure is an ion-acoustic soliton. It is worth noting that this structure might be the one detected experimentally by Romagnani and collaborators by the proton imaging technique [4]. These authors report, in fact, the detection of a laminar structure, almost spherically symmetric, ahead of a warm expanding plasma, as resulting from laser solid ablation, through a colder and more rarefied one. The electric field associated with this front experimentally shows a bipolar distribution with a maximum amplitude of 3 × 10 7 V m −1 . Also in the experimental case, such a structure is seen to have a width comparable to the Debye length of the rarefied plasma and is seen to propagate with a constant velocity of the order of the ion-acoustic speed of the tenuous plasma. It might be possible that the shock, which the simulations show to be ahead of the soliton, has also been detected in this experiment. Simulations show that the shock should be associated with a monopolar electric field peak with maximum amplitude of the order of 2-3 × 10 7 V m −1 . Such an electric field should induce a deflection of the probing proton beam that is consistent with the one detected ahead of the soliton (see region IV in figure 1(a) of [4]). Finally, also the chaotic deflection pattern detected in correspondence to the solid foil (see region I in figure 1(a) of [4]) might be explained by the simulations. The drift instability mentioned above triggers in fact a strongly modulated electric field distribution, in correspondence to the core of the dense plasma (see −20 x 10 in figure 7(a)), which reaches a maximum amplitude of 2 × 10 7 V m −1 , i.e. strong enough to induce the proton deflections detected in the experiment.

Filamentation onset at the shock front
An interesting feature directly arising from the analysis of the simulation results is the onset of filamentary structures in the proton distribution at the foreshock region (x = 27 in figure 9). Such structures are seen to have a clear periodical pattern with a spatial wavelength of the order of twice the Debye length of the dilute plasma and to be almost perpendicular to the shock front (average deviation of θ = 10 • ). They are associated with an electrostatic field, polarized along the y-direction, with a maximum amplitude of about 5 × 10 7 V m −1 ( figure 10(b)), suggesting the strong directionality, longitudinal to the shock propagation axis, of these structures. The electrostatic field along the x-direction (figure 10(a)) does not show in fact a similar pattern but is simply constituted by a localized positive peak at x ≈ 28, associated with the shock. In order to obtain a better signal-to-noise ratio, these two figures have been obtained by performing a low-pass filtering of the electric field distribution. First both components of the electric fields have been Fourier transformed only along the x-direction, and subsequently, all the spectral components with a wavenumber modulus larger than a reference value |k (x,m) | have been forcefully set equal to zero. The chosen value of |k (x,m) | is 30 × 2π/30, where 2π/30 represents the minimum wavenumber resolvable in this Fourier transform (since the transform has been performed in a spatial window of x = 30; see figure 10). This implies a minimum modulation wavelength of λ (x,m) ≈ 1. An inverse Fourier transform was then performed to obtain the lowpass real-valued electric field distributions; the noise in figure 10 has indeed a typical spatial wavelength of the order of λ (x,m) , whereas the electric field associated with this instability has a spatial period of the order of 2λ (x,m) . Moreover, the electric field modulations appear to be much longer along x (which is the filtered direction) than the noise fluctuations. No filtering was in fact performed along the y-direction. These two arguments represent a good indication of the physical nature of the striations observed. Indeed, similar electric field structures were also seen in the unfiltered ion distribution.
An interesting feature of these striations is that slow and fast protons (i.e. protons that are slower or faster than the reference velocity represented by the black line in figure 9(a)) appear to be clearly separated spatially (see figure 11) and to propagate in time maintaining, independently, a roughly constant distribution (for the temporal evolution see movie 2, available from stacks.iop.org/NJP/13/073023/mmedia).
Such a filamentation of the proton beam that, for obvious dimensionality arguments, was not visible in the 1D simulations [22], might have been induced by the onset of ion-ion instability [36,37].
Such an instability can be seen as the ionic counterpart of the well-known electron twostream instability as first discussed by Pierce [38] and can occur when a proton beam is streaming through an un-magnetized, collisionless plasma (refer to figure 14 and the discussion to see that no coherent magnetic fields are generated during the entire simulation time). The counter-streaming proton beam is provided by the existence of the shock-reflected ions. The electric field distribution generated is expected to be polarized in the y-direction (see figure 10), whereas the density striations are oriented along the plasma expansion axis (see figures 9 and 11). The typical phase velocities of the waves generated are much smaller than those of the electrons, which can therefore be seen simply as a neutralizing negative fluid (see figure 6(a), where no relevant modulations of the electron population can be seen in correspondence to the filaments). The dispersion relation of such instability was first derived in [36] and can be expressed as [39] where k i is the instability wavevector, k D is the Debye wavevector of the ambient protons, V i = 2c s is the drift velocity of the shock-reflected ions and v thi is the thermal speed of the  figure 9) and (b) the density of slow protons (i.e. slower than 20 times the proton thermal speed; see the black line in figure 9).
ion beam (which the simulation shows to be approximately equal to the thermal speed of the ambient protons). Here Z represents the plasma dispersion function [40], defined as Inserting equation (7) into equation (5) reveals the maximum linear growth rate of this instability to be = (ω) ≈ 5.6 × 10 −3 ω pt , which implies a characteristic time, for the instability to grow, of t i ≈ 178. This time scale is comparable to the time at which the instability is first seen in the simulation (t 2 = 150). Moreover, the instability wavevector is found to scale as k i ≈ 60ω pt /c, which gives a spatial period of λ i ≈ 2.4λ De , in good agreement with the simulation results. The absence of a magnetic field (see figure 14 and the discussion) is already a good indication that such filamentary structures are not the consequence of the Weibel instability which is instead triggered if the plasma expands into a vacuum [18]. The Weibel instability is also invoked to explain part of the numerical results reported in [39]: a comparison of these results with our numerical results is given in section 5. The Weibel instability requires in fact the electrons of the dense plasma to be reflected from the shock front, which moves at the ionacoustic speed. The electrons should then significantly lose kinetic energy along the ions' front. A thermally anisotropic electron velocity distribution should then develop which is responsible for the Weibel instability [15][16][17] and thus for a magnetized plasma front [18]. This picture, however, is significantly changed when a rarefied medium is introduced. The electrons in the dense plasma are in fact seen to not develop significant deviations from the initial Maxwellian distribution, as is clear from figure 12. This figure in fact shows that only symmetric deviations of the order of 5% will be present by the end of the simulation. This is because the heated, super-thermal electron population present in the ambient medium will act as a suppressor for consistent deviations in the Maxwellian distribution. Such deviations are indeed insufficient to develop detectable magnetic fields, as figure 14 confirms.
It is interesting to note that this instability might explain the filamented front of the shock waves reported in [4]. In that paper, modulations of the probing proton beam are detected in the foreshock region (see region III in figure 1 of [4]) generated during the expansion of a dense plasma, resulting from laser-solid ablation, through a rarefied ionized medium. The striations detected are purely electrostatic, with a typical width of 2-3 Debye lengths of the rarefied plasma and appear to be almost perpendicular to the shock fronts, all elements in agreement with our simulation results.

Discussion
A general summary of the plasma state evolution, and all the features discussed in the text, can be extracted by looking at the electric field distributions in both the xand y-directions throughout the entire simulation time ( figure 13). It is important to note that in the y-direction, the average of the energy density must have been performed. Directly averaging the electric field itself would have given trivially zero, owing to the oscillatory behavior of this component of the electric field along this direction. In the x-polarized electric field distribution ( figure 13(a)), a strong signal, located around (x, t) = (0, 0), is a signature of the initial charge imbalance that triggers the forward ion acceleration and shock creation (see the description of the plasma state at t 1 = 25). The shock itself is then represented by an intense line of electric field that moves, towards positive x, with an approximately constant velocity of the order of the ion-acoustic speed of the rarefied plasma and an amplitude in the range of a few tens of megavolts. As discussed previously, such an electric field progressively decreases its amplitude and slightly broadens its width. However, it is seen to neither decrease its propagation velocity nor vanish, a good indication of the stability of such a structure, at least for time scales of relevance to laser-plasma experiments [4]. On the other hand, the linear color scale of this graph does not allow us to resolve the presence of an ion-acoustic soliton tied to the downstream part of the shock. A further strong electric field region is visible, starting approximately at (x, t) = (−15, 400) and moving towards negative x. This is the field associated with the onset of the drift instability. Such a field is seen to increase in time until a saturation value is reached at approximately t = 650 and then progressively decreases and broadens into a much weaker signal. This weakening of the electric field is due to a turbulent stage of such an instability in which the electric field is transferred to the y-direction (compare the same temporal window in figure 13(b)). The plateau of constant electric field present between the drift instability and the shock is associated with the part of the rarefaction wave that is not affected by the presence of an ambient medium. It is interesting to note the clear dip in such an electric field distribution which is located just at the left of the shock and starting approximately at t = 300. This is related to the trapping of the protons counter-streaming from the dilute plasma which develops in the ion-acoustic soliton discussed above.
Secondary electric field structures are present ahead of the shock, moving approximately three times faster than the ion-acoustic speed of the rarefied plasma or, equivalently, with a speed equal to 60 times the proton thermal speed. Such structures are associated with the tip of the rarefaction wave (see, for instance, the peak at x ≈ 14 in figure 3(d)). This electric field has not been detected in experiments carried out in similar plasma conditions (such as [4]), possibly suggesting that its amplitude lies below the threshold of the proton imaging technique (which in the case of [4] was of the order of a few times 10 6 V m −1 ). This result thus seems to be in disagreement with experimental findings. However, it has to be noted that the rarefaction wave in the simulation is launched from a sharp density jump, whereas in the experiment it will be the consequence of a much smaller and extended density gradient. It is thus reasonable to conclude that, in this latter case, the tip of the rarefaction wave will be associated with a much smaller electric field, possibly lying below the detection threshold of the proton probing technique. Further numerical work is currently under way in order to clarify this disagreement. It is interesting to compare figures 13(a) and (b), i.e. the xand y-components of the electric field distribution along the x-direction. In figure 13(b), the shock position (which at t = 900 is highlighted by the black line in figure 13(a)) is evidenced by the overdrawn black line. It is thus clear that the filamented structure of the shock front originates only at t = 600 and it is associated with a strong electric field, always located at the foreshock region, which is only polarized along the y-direction (no similar electric field distribution can in fact be seen in the x-component of the electric field). The onset of the drift instability, starting at t ≈ 400 and located at y < 0, is also visible.
A peculiar result of these simulations is the evidence that significant magnetic fields never grow in the plasma for the entire simulation time. Figure 14 shows the magnetic field distributions, at the different simulation times discussed in the text, averaged over the y-direction. As can be clearly seen, the amplitude of the magnetic field remains within the noise induced by the PIC representation of the plasma. In particular, the main consequence of such a result is that shocks and ion-acoustic solitons generated during the expansion of a warm and dense plasma into a colder and more rarefied one are purely electrostatic modes. Also the drift instability appears to be driven exclusively by density gradients.
The absence of a magnetized shock front appears to be in contradiction to recent theoretical results dealing with plasma expansion into a vacuum [18]. This work shows that, during the expansion, the electrons will be mostly cooled down along the expansion axis, due to the transfer of energy from the electrons to the expanding protons. This sets an anisotropy in the electron velocity distribution that induces magnetic fields of the order of the MegaGauss that feed the Weibel instability at the plasma front. It is shown that the cold population of electrons is eventually able to inhibit the instability only in a preliminary stage but does not play any significant role in the second expansion stage [18]. Magnetic fields are therefore always created regardless of the presence of a cold electron population in the expanding plasma. These results significantly differ from what was deduced by our simulations, where a rarefied plasma is introduced as an ambient medium. In this case a superthermal population is created, which will drastically reduce the growth rate of the Weibel instability, possibly hampering its onset.
Also the unstable shock front can be explained by taking into account only electrostatic phenomena. In particular, it appears to be consistent with the onset of the ion-ion instability as first numerically observed by Karimabadi and collaborators [41]. However, the topological distribution of the filaments obtained in our simulations closely resembles that obtained by Tako and Takabe [39] while simulating the creation of a collisionless shock in a neutral plasma with an ion-to-electron temperature ratio of about 1/9. That paper indicates that collisionless, ion-acoustic shock fronts are unstable to the ion-ion instability and, at a later time, to the Weibel instability as well. These authors claim that such an unstable front causes the vanishing of the shock at late times. However, the disruption of the shock reported in [39] might be caused by their different simulation setup. The ions (and electrons) are reflected at one boundary of the simulation box in [39] by assuming boundary conditions which correspond to a reflecting wall. Two ion beams are present close to the wall. The incoming ion beam moves with the speed v b to the wall and the reflected one moves with the speed −v b away from the wall. As long as no shock has formed, we find a spatial interval in which two counter-streaming ion beams are thus present. The ions of the reflected ion beam are increasingly modulated by their interaction through instabilities with the incoming electrons and ions. This modulation is strongest close to the front of the reflected ions as PIC simulations of colliding plasma slabs in a parallel magnetic field geometry have demonstrated [42]. These modulations eventually evolve into a shock. The shock forms far from the reflecting wall if a high ion-to-electron mass ratio is used, which is in agreement with the simulation results in [39]. Once the shock has formed, the incoming ions are thermalized in correspondence to the shock and they can no longer reach the reflecting wall in the form of a beam. This, in turn, implies that the reflected ion beam is cut off. However, the shock has formed out of the collision between the incoming and the reflected ion beam and it must collapse once it is no longer fed by the reflected ion beam. The shock in figure 3 of [39] forms at a time of 1000 inverse plasma frequencies and its collapse becomes apparent before a time equal to 2000 inverse plasma frequencies. The shock forms at a distance of 12 electron skin depths away from the reflecting wall and the incoming plasma has the flow speed of 0.1c.
Ions flowing at 0.1c cross 24 electron skin depths in 240 inverse plasma frequencies. The shock may not collapse immediately. It is thus possible that it is the cut-off of the incoming ions by the shock formation that lets the shock in [42] collapse, and not the ion-ion instability.

Conclusions
We have carried out a 2D PIC simulation aimed to study the propagation of a dense plasma into a more rarefied medium. This study represents the extension of a preliminary 1D study with similar initial conditions. The plasma expansion is first governed by the electrons, which, thanks to the favorable charge-to-mass ratio, start to diffuse into the rarefied plasma, leaving the protons behind. A charge imbalance will thus be set that is responsible for the acceleration of the protons of the tenuous plasma and the onset of a rarefaction wave. Such a wave is then disrupted by the onset of a collisionless shock structure which moves at a constant velocity of the order of the ion-acoustic speed of the tenuous plasma. The shock front and the foreshock region are seen to be subject to electrostatic ion-ion instability, which generates longitudinal filamentary structures in the upstream part of the shock. In the downstream part of such a shock, an ionacoustic soliton is formed, generated by the trapped protons of the tenuous plasma. The purely unmagnetized nature of the plasma, which is evident from the simulation results, indicates that no Weibel instability is triggered and that these filaments might be caused by the onset of an electrostatic ion-ion instability. The absence of any coherent magnetic field structures in the plasma is clear evidence of the purely electrostatic nature of shocks and solitons created in such conditions. Meanwhile, in the core part of the dense, expanding plasma, the onset of an electrostatic drift instability is observed. Such a plasma dynamics significantly differs from analytical and numerical treatments of plasma expansion into a vacuum. In addition to the relevance of this work for the understanding of non-relativistic shocks in solar and astrophysical plasma scenarios, these results allow a clear interpretation of experimental data obtained in similar conditions.