## 3 PIC Simulation Results
The velocity is directly proportional to the laser intensity in the HB model, and thus it increases on the rising edge of the pulse up to a maximum value (\(t\sim\frac{40}{\sqrt{2}}\) fs) and decreases thereafter. In the simulation, whereas the increasing velocity is driven by the same increase in intensity, the decreasing velocity after the laser peak occurs becausethe thermal pressure starts to compete again with the radiation pressure when the intensity drops. For this relatively thick target, the compressed electron layer does not reach the target rear side before the laser radiation pressure decays and thus the dynamics are fully accounted for by HB. At the other extreme, for the case \(l=10\) nm, the compressed plasma electron layer quickly reaches the target rear. In this case, as shown in Figure 4b, the PCDS velocity results from the simulation are explained by the LS model, until RSIT occurs (at \(\sim\)56 fs), at which point the laser pressure becomes less efficient and the electron acceleration due to radiation pressure ceases. Consequently, relativistic transparency limits the maximum velocity of the PCDS. We note that the agreement is with the LS model when the Gaussian temporal profile of the laser pulse is accounted for (as opposed to a constant intensity profile). The two target thicknesses in between correspond to hybrid scenarios in which the dynamics are explained by the HB model up to a point in time at which the compressed electron layer reaches the target rear and the LS model thereafter. The corresponding time is \(\sim\)20 fs and \(\sim\)38 fs for \(l=40\) nm and \(l=100\) nm, respectively. All of these simulation results correspond to the velocity of the PCDS along the laser propagation axis, where the radiation pressure is highest. The off-axis velocities are plotted in Figure 3b,d,f,h, which shows a clear dependency on the target thickness. The thickest target (\(l=200\) nm) shows that the highest velocities are extended over a broad distribution around the interaction region, whereas in \(l=10\) nm target, they are distributed within a narrow region around the laser axis, and moreover,  
Figure 3: (**a**) Combined 2D plot of electron density and laser intensity, from a 2D particle-in-cell (PIC) simulation with \(l=10\) nm, at the time step corresponding to the moment that the target becomes relativistically transparent; (**b**) 2D PIC simulation of the plasma critical density surface (PCDS) velocity as a function of time and position along the Z-axis, again for \(l=10\) nm. The red shaded region corresponds to the parameter range for which the target has become relativistically underdense and therefore the PCDS no longer exists. The green dashed curve is the laser pulse intensity profile as a function of time, for reference. Corresponding plots for: (**c,d**) \(l=40\) nm; (**e,f**) \(l=100\) nm; and (**g,h**) \(l=200\) nm. Insets in (**b**) and (**h**) show the relativistic critical density velocity as a function of time and position as calculated using the LS and HB models, respectively (for a Gaussian laser temporal profile). on the falling edge of the laser pulse, just prior to the onset of RSIT. The grid-shading corresponds to the parameter space for which RSIT has occurred. The inset in Figure 3b,h shows the velocity of the PCDS, in time and space, along the Z direction, predicted by the LS and HB models, respectively. These model calculations are in very good agreement with the simulation results over the full laser focal spot. Transitional characteristics between these two extremes are observed for \(l\) = 40 nm and \(l\) = 100 nm. Figure 4: **(a)** 2D PIC simulation results showing the evolution of the electron density and normalized perpendicular component of the laser electric field squared, \(F_{y}^{2}\), along the Z = 0 axis as a function of position in X and time, for \(l\) = 10 nm.