Formation and evolution of post-solitons following a high intensity laser-plasma interaction with a low-density foam target

The formation and evolution of post-solitons has been discussed for quite some time both analytically and through the use of particle-in-cell (PIC) codes. It is however only recently that they have been directly observed in laser-plasma experiments. Relativistic electromagnetic (EM) solitons are localised structures that can occur in collisionless plasmas. They consist of a low-frequency EM wave trapped in a low electron number-density cavity surrounded by a shell with a higher electron number-density. Here we describe the results of an experiment in which a 100 TW Ti:sapphire laser (30 fs, 800 nm) irradiates a 0.03 gcm−3 TMPTA foam target with a focused intensity Il=9.5×1017 Wcm−2 . A third harmonic ( λprobe≃266 nm) probe is employed to diagnose plasma motion for 25 ps after the main pulse interaction via Doppler-Spectroscopy. Both radiation-hydrodynamics and 2D PIC simulations are performed to aid in the interpretation of the experimental results. We show that the rapid motion of the probe critical-surface observed in the experiment might be a signature of post-soliton wall motion.


Introduction
The hydrodynamic evolution of over-dense foams has been the subject of considerable study in the longer pulse regime (pulse duration >1 ns) and with intensities between 10 13-15 W cm −2 (see [1][2][3][4][5][6]), the production of shock-like perturbations by ultra-short ultra-intense laser pulses is less well understood. Such phenomena could have impacts on the possible use of foams in fast particle and x-ray production [7,8].
In this letter we present a pump(ω)-probe(3ω) Doppler spectroscopy experiment on low-density foam targets using the UPHILL 100TW Ti-Sapphire laser system at the Tata Institute for Fundamental Research (TIFR) in Mumbai. The density of the foam target, at 30 × 10 −3 g cm −3 , is larger than the pump critical density of 5.5 × 10 −3 g cm −3 , whilst still smaller than the probe critical density of 50 × 10 −3 g cm −3 , so that n crit < n foam < n 3ω .
Relativistic EM solitons are localised structures that can occur in collisionless plasmas. They consist of a lowfrequency EM wave trapped in a low electron density cavity surrounded by a shell with a higher electron density. The relativistic electron motion is strongly coupled with the trapped electromagnetic mode, so that both particles and fields contribute equally to the total energy of the soliton [19].
Post-solitons form as a result of the formation of relativistic EM solitons in a plasma. Relativistic EM solitons form as a result of dispersion effects [9,15,16] due to the finite inertia of the electron and the relativistic increase in the effective mass of the electron. PIC simulations [18,28] show EM solitons are generated in the wake [29] of an EM pulse and propagate at a velocity significantly lower than the speed of light, and with a frequency significantly lower than the local plasma frequency.
A laser pulse travelling through a region of low density plasma is susceptible to a range of instabilities which either generate reflections, or alter the frequency of the pulse itself. If these phenomena generate EM waves with frequencies corresponding to the local plasma frequency then they can become trapped. If these trapped EM waves are of a large enough amplitude the ponderomotive force pushes electrons out of regions of high field, lowering the local electron density at the centre of the pulse and raising the density around the edges of the pulse. Ion motion can be neglected during the formation of a soliton as the oscillation period of the field inside the soliton is of the order [18] of 2πω −1 pe .
The net effect of the density distortion is that the dielectric function of the plasma in the central region is significantly increased, whilst it is decreased around the edges. The decrease in the dielectric function is not only a result of the increased electron density, but also due to the increase in the electrons' Lorentz factor γ, in a similar manner to the relativistic self-focusing effect [27,30].
The ions are pushed out of this region at a much lower speed by the charge separation effect of the electrons that have already exited the region. The charge separation can be significant enough to result in ion acceleration due to a coulomb explosion [20]. On the much longer timescale of ion dynamics, the soliton ceases to be a soliton, but a low frequency EM wave packet can still remain confined in the slowly expanding plasma cavity. Over time this low frequency EM structure can leak low frequency radiation [18,20,28] (though higher than the oscillation frequency of the EM field in the soliton) and so gradually lose energy.
The experiments described and modelled here provide some tentative further evidence of the formation and existence of such solitons as a result of high-intensity laser-plasma interaction, but also provide what may be the first direct experimental measurements of post-soliton wall ion velocities via pump-probe Doppler spectroscopy.

Foam target pump-probe experiment
The experiment was performed at the Tata Institute of Fundamental Research with a chirped-pulse-amplified 100 TW Ti:sapphire laser system, which can provide 30 fs laser pulses at a central wavelength of 800 nm. The p-polarised pump pulse has a 500 ps contrast of ∼10 7 , and is focused to a spot size of 14 µm diameter using an f/3 off-axis parabolic mirror at an angle of incidence of 40 • . The resulting laser intensity is The probe pulse is created by splitting off a small amount of the pump beam and up-shifting it to the third harmonic (λ probe ≃ 266 nm). This probe pulse is then fed through a motorised delay-stage which can delay the pulse with micrometer resolution. The probe is directed onto the target with a spot size of ∼ 50 µm, at approximately 4 • from the surface normal. The probe is centered at, and completely overlaps the pump interaction region.
The foams are formed by polymerisation in-situ [7] from a polymer TriMethylol Propane TriAcrylate (TMPTA, C 18 H 20 O 6 ). The foam has a density of 0.03 g cm −3 , which, if fully ionised (Z = 3.85), gives an electron density of 9.64 × 10 21 cm −3 (equivalent to 5.5n crit ). The fully ionised electron density is smaller than the probe-critical density of n 3ω = 1.6 × 10 22 cm −3 , which would allow the probe pulse to pass through the foam with little or no reflection if the foam were An example of several reflected probe spectra including fits used to determine the peak wavelength. The top left plot shows a reference spectrum of the probe reflecting off the frame, with a Gaussian fit (red) and vertical dashed lines indicating the peak wavelength and fitting errors. The bottom left shows the spectrum of a 10 ps delayed probe, with a gaussian fit and peak wavelength with errors. The right plot shows an example of a 4 ps delayed reflected probe spectrum with multiple peaks, the green indicates the summed peak fits, with the individual peak fits shown as dashed lines coloured red, magenta, and blue. of a uniform density. The foam has a microscopic structure (characterised using the process described in [31]) consisting of pores ∼900 nm in size with a distance between pore centres of approximately 1400 nm.
If some of the foam were to remain intact through the prepulse then it is possible that the probe could interact with a plasma up to four times higher than the macroscopic density (i.e. higher than the pre-interaction average density). However, this scenario is unlikely if there is a significant compressive shock formed by the pre-pulse.
Sample results from the experiment are shown in figure 1. The probe spectrum (seen in the top left plot) has a very narrow spectral width ∆λ = 0.9 ± 0.2 nm. A typical reflected spectrum of a pump-probe shot can be seen in the bottom left of figure 1. This shows a significantly wider and red-shifted peak, with an additional wider base on the shorter wavelength side of the peak. On several occasions, most notably during the experiments with a 4 ps and 10 ps delay setting, more than one peak is observed: with a wider base peak; a sharp narrow peak and a highly red-shifted smaller peak. An example of results with multiple peaked reflections is shown in the right plot of figure 1.
To obtain a value for the spectral location of each of the spectral peaks from the probe's reflected-spectra, first a noise baseline with an amplitude of A b is fitted to the noise floor. Once the noise floor is found the data is fitted to a function of the form: where λ n is the spectral location of the nth peak, A n is the amplitude of the nth peak, and σ n is the spectral width of the nth peak. An example of a fit with multiple peaks can be seen in figure 1. At least four reflected spectra are gathered for each delay. For each foam-reflected spectrum gathered, a reflection of the The results are subdivided into three sets, the blue set indicate the peaks with the least deviation from the probe wavelength, green have a large red-shift, and the magenta set have a very large blue shift. The calculated probe-surface velocity is denoted on the right hand axis of both plots. Here the positive velocity direction was chosen so that it represents material flowing towards the target surface and away from the spectrometer, while material ablating away from the target surface would have a negative velocity.
probe onto the frame housing the foams is also gathered as a reference spectrum.
In the spectra with multiple peaks the peaks can be broadly separated into three groups. The first group with a very large positive Doppler shift (12 to 18 nm, i.e. a large red-shift) corresponding to the smaller highly shifted peaks seen in figure 1. A second group with a moderately large negative Doppler shift (−0.76 to −2.5 nm, i.e blue-shifted), corresponding to the smaller broad response seen in figure 1 to the left of the sharper peak. The final group which appears to follow a similar trend to earlier experiments with smaller Doppler shifts [32], corresponding to the sharp peaks seen in figure 1. To calculate the red/blue-shift of a probe-foam reflection, the spectral location of the probe-frame reflection is used as the initial wavelength, so that ∆λ = λ foam − λ frame . All the ∆λ values for each delay setting (separated into groups as appropriate) are used to calculate an average red/blue shift for each group at each delay setting. These average values are shown in figure 2.
There is a significant amount of variation in the Doppler shift of the probe reflection from shot-to-shot for the same delay setting. This shot-to-shot variation is significantly larger than both the spectral width of the peaks from the probe-reflected spectra and the error from fitting, and so the standard error from the shot-to-shot variation is used to calculate error bars for the resulting time evolution of the reflected probe Doppler-shift. Whilst the error bars shown in figure 2 are very large, there is still a visible trend: an early time outward motion of approximately 40 × 10 6 cm s −1 ; followed by an inwards motion between 2 and 6 ps (though this could vary significantly given the large variation in the data at 4 ps); and then finally an outwards motion is observed gradually falling back down towards zero velocity between 10 and 25 ps.
An additional delay setting is also used in an attempt to characterise the pre-plasma conditions, the delay stage is set to −20 ps, i.e. 20 ps before the impact of the main pulse. The measured Doppler-shift of the negatively delayed probe reflection is measured at (0.39 ± 0.54) nm with a corresponding velocity of (22 ± 30) × 10 6 cm s −1 .
The results in blue in figure 2 are possibly Doppler-shifted as a result of a shock-like phenomena travelling through an over-dense plasma similar to that seen in solid density targets [33,34]. However, the larger wavelength shifted results are unlikely to be the result of Doppler-shifting from reflections off moving over-critical plasma, simply due to the much larger Doppler-shifting in both cases. Whilst there is likely to be material travelling at high speeds away from the target in the blow-off region, this material is unlikely to be dense enough to affect the probe-pulse. However, depending on the scale-length and collisionality of the pre-plasma, there maybe other phenomena that can account for the observed shift in wavelength.

Radiation-hydrodynamic modelling of pre-pulse absorption
The HYADES radiation hydrodynamics simulation code is used [35]. A foam with a macroscopic electron density n m (considering the foam to be fully ionised) consists of pores which are either empty or have a density n pore ≪ n m and strands with a density n strand > n m . The size of the pores and the fractal-like structure of the foam dictates the density of the strands. A more fractal-like structure will lead to narrower, denser strands. The very early time pre-pulse foam interaction will form a highly non-equilibrium plasma consisting of ablating strands and pores filled with hot low-density plasma [4]. The time taken for this plasma to equilibriate will depend on the transfer of energy from the pre-pulse through the forming plasma and how quickly a shock-like feature forms. If the plasma is under-dense for the majority of this interaction then the foam could heat volumetrically to form a large nonequilibrium plasma [2,36].
If the target ionises quickly to an over-critical density then the heating process is likely to be dominated by the formation of a shock-like feature forming around the critical density region, similar to the behaviour seen in homogeneous targets.
At higher intensities (10 13-14 W cm −2 ) the absorption of laser light into foam targets is likely to be significantly altered, with much higher absorption than solid targets [3]. At later times, in cases with high intensity, greater penetration of the laser would lead to the formation of a cavity inside the foam [4]. It is not clear if these effects are likely to take place at the lower intensity seen in the pre-pulse (∼10 11 W cm −2 ) and at timescales shorter than 1 ns. With a pore-like structure on the front surface (scaled at ∼900 nm) that has a similar scale to the pump wavelength (800 nm) it is also possible that surface effects could play a part in laser absorption [37,38].
While there are theoretical models that describe the hydrodynamic motion of foams [2] under intense laser illumination, these are not generally implementable in standard fluid codes due to the fact they require a two-scale laser absorption model. Two such codes that have had some success at this include the 2D ALE hydro-code PALE [5] and the 1D hydro-code MULTI [6].
Since the target is ∼50% hydrogen atoms by number and the target has a macroscopic density ∼5n crit if fully ionised, it is possible that the target ionises quickly enough to a density > 2n crit that it might be possible for the foam-like structure to behave in a similar way to a homogenous low density target.
An average-atom LTE ionisation model is used along with the polystyrene equation of state from record 7592 of the SES-AME library. The multi-group radiation package is used with 100 photon groups from 10 −3 -10 −1 eV up to 10 keV, with the correct atomic ratios for the TMPTA foam. The laser solver uses parameters of λ = 800 nm, θ = 40 • , and plane polarisation, along with a pre-pulse profile consistent with experimental setup. The flux-limited thermal conductivity model is used with a flux-limiter value of 0.3. The extra transport allowed by a higher flux limiter may offer a better description of the foam target dynamics given the possibility of greater penetration by the pre-pulse laser [2].
In order to attempt modelling the foam structure a structured mesh using alternating densities between 6 × 10 18 cm −1 for 1.2 µm and then 4 × 10 18 cm −3 for 0.8 µm was tested against a uniform grid with a homogeneous density of 5.4 × 10 18 cm −3 . However both simulations returned almost identical simulation outputs after the very short period where the profile with alternating densities was being ionised.
To further investigate the interaction of the high-intensity pulse with the front surface pre-plasma, two separate density profiles are obtained from the HYADES output. The first density profile is obtained directly from the HYADES simulation initialised using the modulated profile at a time t = 490 ps. The density profile is taken from the HYADES results 10 ps prior to the onset of the main pulse so that only the effect of the prepulse is taken into account. The second density profile uses a steepened version of this profile. This approach is taken as a shorter length-scale pre-plasma may be a better approximation of the density profile in the experiment. In reality the preplasma may have a more spherical expansion profile which is better approximated by a shorter 1D profile. This profile is calculated by making a fit to the existing density profile and reducing the characteristic length-scale L. The best fit was found using the function: where a is a minimum density, b is the maximum density, d is the location of the maximum density, and L is the characteristic length-scale. This function was chosen as the much longer decay of the density from the front surface was not well represented when using a Gaussian or super-Gaussian function.
The best fit to the HYADES profile gives a scale length of 10.0 ± 0.2 µm, the reduced profile is set to 7 µm.

EPOCH particle-in-cell modelling
The 1D version of the EPOCH PIC code is used [39] to model the main pulse interaction. Several simplifications are assumed for all of the PIC code simulations performed in studying the main-pulse foam-target dynamics. The first simplification is that ionisation is assumed not to have a significant effect, because the foam pre-plasma is seen in HYADES simulations, to be almost completely ionised up to the densest part of the pre-plasma shock front. The second simplification is that of using an average ion (with an atomic mass of 7.2m p and charge of 3.85e) rather than simulating multiple species of ions. This was necessary as it significantly reduced the number of particles required to achieve a convergent result. A cell size of 5 nm is used with 4000 particles per cell evenly divided between ions and electrons. (The particles referred to here are in fact macro-particles which represent 'clouds' of real particles. Each macro-particle being statistically weighted according to the density profile in the initial setup). Collisions are calculated using a Coulomb logarithm which is a function of the local density and temperature. The laser is set to an intensity of 10 18 W cm −2 and has a gaussian temporal profile with a FWHM of 30 fs. Results from simulations using both the HYADES density profile and the shortened density profile can be seen in figure 3. Both of these plots are obtained at ∼200 fs after the reflection of the main pulse.
There are two groups of velocity perturbations visible in the longer scale-length pre-plasma (see figure 3(b)): a large perturbation around the n e ≃ n crit /4 and a much smaller perturbation around n e ≃ n crit . Around the same locations as both sets of velocity perturbations, are temperature peaks. Leading up to the larger of the two temperature peaks is a region of heating, which appears to be a part of the phenomenon leading up to the larger of the two velocity peaks, and will be explored in the next section. The larger of the two sets of perturbations is visible in both the longer scale-length pre-plasma calculation (seen in figures 3(a)-(c)) and the shorter scale-length pre-plasma (seen in figures 3(d)-(f)). The smaller perturbation, located at a denser region in the longer scale-length pre-plasma, does not appear to exist in the shorter pre-plasma results. Interestingly this smaller perturbation does not appear when these simulations are run without collisions.
A hydrodynamic calculation can now be initialised by combining outputs from the EPOCH simulation above, with T i and u f outputs from the original HYADES calculation 10 ps prior to the main pulse. To test whether the smaller velocity perturbation (visible in the inset in figure 3(b)) is related to the phenomena observed in the experimental results, only the region between −10 µm and 40 µm is simulated initially. While there is a shock formed from this perturbation, it does not travel far enough into the denser region of the pre-plasma to be seen by the n 3ω probe. Also the shock front travels at a significantly smaller velocity (∼4 × 10 6 cm s −1 ) than any of the experimentally measured phenomena.
The outward expansion of the pre-plasma at densities greater than n crit , which initially starts at 4 × 10 6 cm s −1 travelling into the target, and then gradually turns into an outward There are significant perturbations leading up to a region near n crit /4, at which a large temperature spike is observable with concurring disturbances in the density and velocity plots. The inset plot highlights a comparatively small velocity perturbation of 20 × 10 6 cm s −1 close to the critical density which is not observed in the shorter pre-plasma simulation. These profiles plotted here are used as the initial conditions for the hydrodynamic simulations of the shock-like perturbation after the main pulse interaction. motion of 2.5 × 10 6 cm s −1 at later times. The motion of the dense region of the target also does not appear to relate to phenomena observed in the experiment, and so further investigation into the lower density behaviour is required.

Instability analysis
As the laser pulse travels through the under-dense plasma it suffers some significant degradation, possibly as the result of at least one plasma instability. To ascertain the nature of the instability, spectral analysis is required. Whilst both the longer and shorter scale-length simulations show significant disturbances, the longer scale-length pre-plasma is the simplest to analyse, as the instability occurs over a longer distance. The occurrence of many types of plasma instability are affected by collisions and so the analysis is performed with and without collisions switched on.
Using standard Fourier analysis techniques, pulse chirping, where the frequency or wavelength becomes a function of x or t, cannot be clearly resolved. By performing a Wigner-Ville transform on the laser pulse, rather than a Fourier transform, a determination of the wave-number k as a function of x is possible. The amplitude of a particular point, a(x, k), in the Wigner-Ville transform can be understood as a quasiprobability that a wave, with a wave-number k, is travelling through the point at position x. The 'quasi' nature of the quasiprobability comes from the fact that the amplitude a is also proportional to the amplitude of the original wave at this point. The sacrifice made in this type of analysis is that the output can suffer from some interference patterns. Figure 4 shows Wigner-Ville transforms of two different simulations with two different laser pulses, the first with t pulse = 30 fs, the second with a shorter. t pulse = 8 fs. The first simulation ( figure 4, left) shows a pulse as it passes through plasma where n e ∼ n crit /4. The Wigner-Ville transform of the laser pulse, 220 fs later in time, show significant disturbance of the pulse from instabilities in the plasma. The wave matching criteria for forward Raman scattering requires that the wave-number k s = k l ± k pe , where the subscript s refers to a scattered wave, l refers to the incoming laser pulse, and pe refers to a co-propagating Langmuir wave. The criteria for Raman Backwards scattering requires a wave travelling in the opposing direction and has wave matching criteria of k s = k pe − k l . Distinguishing between these two phenomena is not always possible.
In order to determine whether forward or backward Raman scattering is dominating the interaction the forward and backwards travelling EM waves can be analysed separately by normalising to the fields E p = m e ω p c/e and B p = m e ω p /e and computing the forward propagating EM wave profile a f = (E y + B z ), and a b = (E y − B z ). The results from this can be seen in figure 4.
The laser pulse, when travelling through a vacuum, has a wavelength of 800 nm and so has a wave-number of k = 7.85 × 10 6 m −1 . So with the forward scattering wavematching criteria there should be scattered light wave-numbers of |k f+ | ∼ 12 × 10 6 m −1 and |k f− | ∼ 4 × 10 6 m −1 , which are drawn on figure 4. The criteria for backwards scattering give a Despite the larger growth rate of backwards Raman scattering, the peak from the backwards travelling EM wave signal shown in figure 4 around k ≃ 4 × 10 6 m −1 is smaller than the forward scattering peak. This discrepancy is likely due to the fact that the growth rate for backwards Raman scattering does not take into account very short pulse lengths, while the forward Raman scattered wave travels with the laser pulse and so continues to grow.
Whilst it is possible that Raman scattering is the cause of the significant red-shifting of part of the main pulse, it may not be the sole cause. To eliminate the SRS instability a shorter pulse, shortened from 30 fs to 8 fs, can be used to see if the pulse still endures a significant amount of distortion despite insufficient time for SRS instabilities to grow. Figure 4 shows the Wigner-Ville transform of the shortened 8 fs pulse at the same point in time (220 fs) and along the same density profile as the previous pulse. The pulse undergoes a significant chirp with a large amount of the pulse red-shifted to k l = 5 × 10 6 m −1 (or λ ∼ 1300 nm), the shift in wavelength of the affected pulse equates to a critical density n shift ∼ n crit /4. However, in this simulation the downshifted light simply reflects at the point where n e = n crit /4 rather than becoming trapped in the plasma. The anomalous red-shift of the laser pulse in this case is strongly indicative of scattering due to photon acceleration [40] (or deceleration in this case) which has been previously observed in laser wake-field settings [41].

Soliton-like disturbance formation and evolution
The large disturbances visible in figure 3 at around the n e ∼ n e /4 regions, in both simulations are re-plotted at later times in figure 5, along with the local electric and magnetic fields. The phenomena observable in the longer pre-plasma simulation shows a deep depression in both the electron and ion density surrounded by denser plasma. Within the depression there is a large-amplitude oscillating electromagnetic field with the same polarisation as the laser pulse. The phenomena in the shorter pre-plasma is similar, though the magnetic field does not appear to correlate, this could be related to the expansion of post-soliton walls. This particular type of feature appears to be very similar to electromagnetic post-solitons observed in similar conditions [9,18,19,25,27]. The low density material in the central depression is surrounded by a higher density region of ejected electrons (presumably ponderomotively pushed out of the higher field region) and ions. The surrounding material is travelling at extremely high cell-averaged ion velocities; with the wall on the low density side having a velocity of 400 × 10 6 cm s −1 out away from the target, and the high density wall having a velocity of 200 × 10 6 cm s −1 towards the target bulk.
Calculations of the longer term evolution of these structures has been attempted using a Lagrangian hydrodynamics code, in this case simply using the PIC code ion density n i , ion temperature T i , and a cell averaged ion velocity u i . The density at the shock fronts, n front , varies significantly, starting at n front ⩽ n crit and then rising as high as n front ∼ 4n crit . The velocity of the walls is compared to experimental measurements in figure 6.
Whether a probe, with a critical density of n 3ω ≃ 9n crit , will interact with these structures depends on the evolution of the density in the shock fronts formed at the walls of the postsolitons. The shock front density depends largely on the surrounding material, as there are several smaller features that produce smaller shocks which can then merge into the postsoliton wall-shocks. The large compression of these wallshocks is due to the initial conditions set in the fluid code from EPOCH.
The initial formation of the walls, due to ponderomotive expulsion of electrons from a region of high electric field, is likely to be close to collisionless, due to the extremely high electron temperature in this region (∼10 5 eV).
The hydrodynamic calculation above assumes that the motion of the wall edges is a collisional hydrodynamic process rather than a collisionless one. This may not be a realistic assumption if the walls are driven by the electric field inside the soliton for a long (more than 0.5 ps) time after they are formed. In the simulations above the electric field continues oscillating inside the post-soliton cavity until the end of the PIC simulation, and so it is unclear whether the walls are driven apart by the ponderomotive force for a significant length of time or not.

2D EPOCH soliton formation
A 2D density profile can be generated by taking a 1D profile in the y-direction and copying it along the x-axis to form a shelf-like, or cliff-like, edge. While this flat-faced type of density gradient is not likely to exactly match the experimental conditions, it can be used to test whether the laser pulse will create similar structures to those seen in the 1D case, if it enters the density gradient at an angle. Two of these profiles are generated, the first profile uses the reduced scale length of 7 µm, the second using a further reduced scale length of 5 µm.
Two EPOCH simulations are run using the Stampede 2 cluster [42]. The first simulation consists of a 1200 × 2400 grid of cells for a box with dimensions 60 × 120 µm 2 (giving a square cell with side 50 nm) and is run with an electron species and an average ion species (A = 7.2m p and Z = 3.85), 100 particles per species per cell are used and the default collision model in EPOCH is used. The laser is set with a λ = 800 nm and entering the simulation box at an angle of 40 • . The second simulation uses a slightly smaller box size of 784 × 1568 maintaining the previous cell dimensions. The aforementioned density profile is used and an initial ion and electron temperature is set to zero. Initially runs with temperature profiles generated in a similar manner to the density profile are attempted, however the dynamics of the laserplasma interaction do not change significantly and the lower signal-to-noise ratio achieved with a cold plasma are preferred. While the cell dimensions are small enough to resolve the laser wavelength adequately they will not resolve the Debye length for the higher density regions. To account for the under resolution at higher densities current smoothing and a third order spline function is used for the particle shape. No numerical heating is observable in the simulations.  figure 7 shows the ion density in the 2D EPOCH calculation approximately 200 fs after the reflection of the main pulse. The pulse reflects off a layer of plasma at y = 20 µm, which is the n crit layer of this density gradient. In front of the critical density layer, at around y = 15 µm are three empty areas surrounded by denser plasma which are very similar to the post-solitons observed in the 1D calculation. However, the electric field quickly leaks out of the depressions and no trapped electromagnetic radiation is observable in the soliton depressions ∼50 fs after they have formed. The ion temperature is displayed in the bottom most plot in figure 7, where it can be seen that the highest temperatures occur around the empty regions seen in the density plot at x = 15 µm, y = 22.5 µm.
To check whether the soliton expansion shown here would give a similar velocity to those observed experimentally, a 1D Lagrangian hydrodynamics simulation is set up using a slice of data from the 2D EPOCH calculation. For the first simulation a slice is taken from the region denoted in a green box in figure 7 and an average across the y-direction is taken. The original input into the 2D EPOCH simulation did not include the prepulse shock front, and so there is no n 3ω density surface for a probe pulse to reflect from at t = 0. The results from this calculation can be seen in figure 8. At early times, 0-2 ps, there is no n 3ω critical surface and so the routine used to calculate the Doppler shift returns zero velocity. Between 2 and 10 ps the wall-shock formed at X = 16 µm (see figure 8) gains significant density and so is seen by the n 3ω probe. The increase in density appears to be due to two of the shock fronts briefly merging. After this time the shock front formed by the soliton wall slows down and loses density to the point where n front < n 3ω . Later in time, between 17ps and 22ps, the slower shock merges with another soliton wall (X = 19 µm figure 8), and is again dense enough to be seen by the n 3ω probe. Also displayed in figure 8 are the outermost and innermost wall-shock velocities. While the innermost wallshock (the wall-shock closest to, and travelling towards, the denser region of the plasma) only briefly gains enough density to be reflected by the 266 nm probe the velocity is consistent at earlier times with several of the experimentally measured points. The outermost wall-shock (the wall-shock furthest from, and travelling away from, the denser region of the plasma) shows similar velocities to the experimental results at later times but doe not merge with any other shocks and so does not gain enough density to be reflected by the probe pulse.

2D EPOCH calculation with 5 µm scale-length preplasma.
The second 2D EPOCH simulation, run with a shorter pre-plasma scale length of 5 µm, shows significantly altered structures at a slightly denser region. Instead of three smaller soliton-like features, a single similar feature is observable, see figure 9 at x = 0 µm, y = −25 µm. Similarly to the previous, longer-scale length, case the ion temperature is high in a very localised region where the soliton-like depression occurs. Profiles for a 1D hydrodynamic calculation are obtained from the 2D simulation in the same way as for the previous 2D EPOCH calculation, and the results are plotted in figure 10 . The first observation to note is that there is only one large depression visible in this simulation, while there are three in the longer scale-length case. The hydrodynamic simulation for this case did not show any density surface dense enough to reflect the probe pulse. The velocity of both the innermost and outermost wall shocks did however match some of the experimentally measured velocities.

1D modelling of main pulse interaction
A large amount of energy appears to go into laser-plasma instabilities observable in the PIC code results. Both of the observed instabilities, the possible stimulated Raman scattering instability, and the photon acceleration or selffocusing red-shifting of the laser pulse, generate post-solitonlike phenomena. These post-solitons generate large shocks at their edges in the low-density regions of the pre-plasma. These wall-shocks progress at velocities comparable to those  figure 9. The blue dashed line with error-bars shows experimental data. The red dashed line is the velocity of the shock closest to the dense region of the target, and the green broken line shows the velocity of the shock furthest away from the dense region of the target. measured in the experiment (see figure 6). However, hydrodynamic calculations of the progress of these shocks show that they do not progress deep enough into the target to reach the densities observable by the probe. An alternative method of modelling the pre-pulse interaction may give a different shaped pre-plasma which may enable these shocks to progress further into the target and thus still offer a possible explanation of the cause of the Doppler shifts seen in the experiment.

2D modelling of main pulse interaction
The laser-plasma interaction in the 2D case generates more post-soliton like structures around the n crit /4 region. The hydrodynamic motion of the wall-shocks surrounding the post-solitons are similar, in both timing and velocity, to the inward velocities measured in the experiment. However they only reach appropriate densities that will reflect a 3ω probe pulse under specific circumstances. The circumstance require that there is more than one post-soliton like depression so that the wall-shocks from one post-soliton can interact with another, and so boost the density to a high enough level to reflect the probe pulse.
No evidence of motion at sufficient velocities to explain the experimental results were observable in the very dense region (≫n crit ) of the target. This leaves only the dynamics at lower densities able to provide an explanation, even though it is a tentative one.

Conclusions
The extra reflected peaks recorded in the experimental results, (see figure 2) do not seem to match any phenomena occurring in either the HYADES or the EPOCH calculations and so no explanation can be made about the possible origins of these points. However, a possible explanation of the central larger Doppler shifted peaks can be made using the modelling shown here.
Post-soliton wall-shock motion might explain the inward motion into the target measured in the experiment (the positive velocities observed in figure 2). Electron heating of the prepulse formed shock, which explains the initial outward motion in the high contrast silicate target experiment, does not appear to explain the outward motion in this case. Looking at figure 6 it is possible that the initial outward expansion could come from an outward travelling wall-shock from a post-soliton-like phenomena, if the soliton can be formed in a region with a high enough initial density. However at later times, these outwards travelling wall-shocks would travel into less dense regions and so would be unlikely to account for later experimental results. These measurements not only support the earlier proton radiographic [24][25][26] and optical measurements [27] in providing evidence of soliton formation in intense laser-plasma interactions, they also provide what may be the first direct measurement of ion-velocities in such post-soliton wall plasmas.

Data availability statement
The data that supports the findings of this study is available from the corresponding author upon reasonable request.