We present 3D Particle-in-Cell (PIC) modeling of an ultra-intense laser experiment by the Extreme Light group at the Air Force Research Laboratory using the Large Scale Plasma (LSP) PIC code. This is the first time PIC simulations have been performed in 3D for this experiment which involves an ultra-intense, short-pulse (30 fs) laser interacting with a water jet target at normal incidence. The laser-energy-to-ejected-electron-energy conversion efficiency observed in 2D(3v) simulations were comparable to the conversion efficiencies seen in the 3D simulations, but the angular distribution of ejected electrons in the 2D(3v) simulations displayed interesting differences with the 3D simulations' angular distribution; the observed differences between the 2D(3v) and 3D simulations were more noticeable for the simulations with higher intensity laser pulses. An analytic plane-wave model is discussed which provides some explanation for the angular distribution and energies of ejected electrons in the 2D(3v) simulations. We also performed a 3D simulation with circularly polarized light and found a significantly higher conversion efficiency and peak electron energy, which is promising for future experiments.

Ultra-intense laser systems offer a diverse range of potential applications and interrogate a variety of interesting physical regimes. While many research groups focus on highly relativistic laser intensities where the associated quiver velocity is large (a01), there are still interesting experiments to perform involving laser intensities near 1018Wcm2 where the electron dynamics are only moderately relativistic (a01). Importantly, at these intensities, current technology allows few-mJ laser pulses to be created at a kHz repetition rate,1,2 which is advantageous for applications and for accumulating statistics. These intensities can also exhibit high reflectivity (70%) from near-solid density targets,3 which is very unlike laser interactions at much higher intensity.4 

As discussed in Orban et al.,5 high reflectivity allows for efficient acceleration of electrons to ∼MeV energies. While the incident laser pulse is reflecting off of the target, the reflected half of the pulse interferes with the incident half of the pulse which results in the creation of a standing wave pattern. The strong electric and magnetic fields in the standing wave can launch electrons in the forward and backwards directions with a momentum p/mc similar to a0 of the incident field, where p is the momentum, m, the mass of the electron, and c the speed of light. This was first noticed by Kemp et al.6 who recognized that this can be used as a mechanism for propelling electrons forward into the compressed fuel of an inertial confinement fusion experiment. Orban et al.5 studied this phenomenon for laser intensities with a01 for the first time and focused on the backwards-directed electrons that can receive an additional energy boost when they are overtaken by the reflected laser pulse. Both the theoretical study of Orban et al.5 and the corresponding experimental investigation of Morrison et al.2 point to large numbers of MeV electrons being produced from individual few-mJ laser pulses normally incident on a water jet target. Recently, MeV energies were confirmed in a using an electron spectrometer.7 

This paper presents the first 3D Particle-in-Cell (PIC) simulations of laser-matter interactions for this experiment. In earlier theoretical studies, Orban et al.5 relied on 2D(3v) PIC simulations using the Large Scale Plasma code.8 As is well known, 2D(3v) Cartesian simulations assume symmetry along the vertical dimension, so while there can be currents in all three dimensions (hence the 3v), particles cannot move in one of these dimensions as they would in a fully 3D real world. A simple application of Gauss' law to this 2D(3v) geometry indicates that electrostatic fields decay with a different radial dependence than they would normally. So although 2D(3v) simulations should be qualitatively accurate, ultra-intense laser-matter interactions are an intrinsically 3D phenomena and must be treated as such.

While many groups have successfully performed 3D PIC simulations of low density targets (e.g., laser wakefield acceleration), it is especially computationally intensive to perform 3D PIC simulations of near-solid density plasmas as exist in the experiment considered here. There can be energy conservation problems if the Debye length, λD, of the simulated plasma is smaller than the grid resolution9 and for near solid density plasmas the Debye length is very small, of order nanometers in scale.10 As discussed in Section II, the LSP code is designed with an implicit algorithm that avoids this grid heating problem. Consequently, one can perform 3D PIC simulations without resolving the Debye length in every area of the grid and still maintain good energy conservation. We present the results of such simulations and compare to earlier 2D(3v) results throughout this paper. We also address, using analytic methods, aspects of standing wave acceleration with a considerably more sophisticated model than was used in Orban et al.5 

The initial conditions and configuration for these simulations will be discussed in Sec. II. Results will be summarized in Sec. III including a comparison to experimental measurements from Morrison et al.2 Finally, a discussion of these results will be presented in Sec. IV. The  Appendix describes an analytic model for electron acceleration in these simulations.

We performed 3D PIC simulations with the LSP code.8 The initial conditions for these simulations were three dimensional analogs to the simulations described in Orban et al.5 For this paper, we use the following Cartesian coordinate system for these simulations: the positive x-axis is the direction of the laser, the y-axis is the polarization direction, and z-axis is the axis of the water column. As shown in Fig. 1, the target is a 27.5 μm × 30 μm × 30 μm (27.5 μm along the laser axis) rectangular block of water-like plasma, consisting of free electrons, protons, and O+ ions, in proportion to make the target match water's chemical composition and to ensure charge neutrality (O+ to p+ to e ratio of 1:2:3). For 10 μm along x, the target consists of plasma at solid density with electron density of 1.0 × 1023 cm−3, and others in proportion. Beyond 10 μm in x, we included a decaying pre-plasma profile with an exponential scale length of 1.5 μm, again along x, up to the target's front edge. The initial density of the target does not vary in the transverse directions. Here, the critical electron density nc=4πmω2/e2 is 1.74 × 1021 cm−3, where m is the mass of an electron, e is elementary charge, ω=2πλ/c is angular frequency of the laser, and c is the speed of light. Thus, the solid density target is 57nc. The electron species had a starting temperature of 1 eV.

FIG. 1.

Initial conditions of the 3D simulations in the left pane and of the 2D(3v) simulations in the right pane, specifically the electron density. The target consists of solid density, colored dark blue, that is 10 μm deep in the laser axis direction (x) and is 30 μm along the transverse dimensions. The critical density, which is 1.74 × 1021 cm−3, is highlighted in gray on both plots. An exponential density of electrons is placed on the surface facing the laser with a scale length of 1.5 μm. The laser comes to focus at (x,y,z)=(0,0,0), which is in the solid density of the target for all simulations.

FIG. 1.

Initial conditions of the 3D simulations in the left pane and of the 2D(3v) simulations in the right pane, specifically the electron density. The target consists of solid density, colored dark blue, that is 10 μm deep in the laser axis direction (x) and is 30 μm along the transverse dimensions. The critical density, which is 1.74 × 1021 cm−3, is highlighted in gray on both plots. An exponential density of electrons is placed on the surface facing the laser with a scale length of 1.5 μm. The laser comes to focus at (x,y,z)=(0,0,0), which is in the solid density of the target for all simulations.

Close modal

In these simulations, an 800 nm wavelength laser pulse is normally incident onto the target, well off-focus with the peak focus being at x=y=z=0 in Fig. 1 which is as in Orban et al.5 The pulse has a sine-squared envelope with a period of 60 fs, or a FWHM of 30 fs. Among the simulations executed, two laser intensities were used, 5.4×1017Wcm2 and 3×1018Wcm2. Also, these studies used a linearly polarized laser for both intensities, with an additional 3D simulation with a 5.4×1017Wcm2, circularly polarized pulse. We note that the 3×1018Wcm2 and 5.4×1017Wcm2 beams had corresponding a0 values at peak focus of 0.5 and 1.2, respectively, where a0=eE0/ωmc, where E0 is the peak electric field of the laser. The incident laser profiles had 2.26 μm (2.6 μm FWHM) and 2.174 μm (2.5 μm FWHM) spot sizes for the 5.4×1017Wcm2 and 3×1018Wcm2 simulations, respectively. Time steps of 0.1 fs, or about 1/30th of the laser oscillation period were used for all simulations. Along each dimension, the spatial resolution was 8 cells per wavelength (λ/8×λ/8×λ/8) for the 3D simulations and 32 cells per wavelength (λ/32×λ/32) for the correspondent 2D(3v) simulations. While these simulations do not resolve the Debye length in every cell (since there are cells in the simulation with near-solid densities and nanometer-scale Debye lengths), the phenomena of interest are electron acceleration in the underdense pre-plasma extending from the target where the Debye length is much larger and more easily resolved by the simulation. The implicit algorithm in LSP is designed to avoid grid-heating issues so that the near-solid density regions do not ruin the overall energy conservation of the simulation. All simulations had 27 macro-particles per cell of each species (free electrons, protons, and O+ ions). As mentioned in Orban et al.,5 O+ ions in this simulation are capable of being further ionized by an ionization model in LSP which follows the Ammosov-Delone-Krainov rate.11 Moreover, electrons scatter by a Monte-Carlo algorithm as in Kemp, Pfund, and Meyer-Ter-Vehn12 with a scattering rate determined by a Spitzer model.13 

The 3D simulations were run for 185 fs or longer; specifically, the 5.4×1017Wcm2 and 3×1018Wcm2 simulations were run for 185 fs and 209 fs, respectively, and the corresponding 2D(3v) simulations were run for 250 fs each. By the end of both 3D simulations, the cumulative number of ejected, high energy (> 120 keV) electrons had reached a plateau, suggesting that at that point in the evolution of that simulation, most accelerated electrons had already escaped the target, and so either 3D simulation would not have needed to have been run any further. An additional 3D simulation was performed with circular polarization but with other parameters being the same as the 5.4×1017Wcm2 simulation with linear polarization.

See Fig. 2 for an analysis of the angle and energy distributions of ejected electrons from the linear polarization simulations run. The bins in these plots have a radial (energy) spacing of ΔE=50 keV for the 5.4×1017Wcm2 simulations and ΔE=100 keV for the 3×1018Wcm2 simulations. For all plots, angular bins of Δϕ=2° were used. The ejection angle for each macroparticle was determined from the momentum of the macroparticle as it left the simulation volume. To compare results between the various simulations which warranted different energy spacings, the charge histogram weights are divided by the spacing per cell (reflected in the unit labels above the color bar). Given that the macro-particles in the 2D(3v) simulations are in actuality line charges, the ejected macro-particles in 2D(3v) simulations cannot be compared directly to ejected macroparticles from 3D simulations. Thus, to obtain units of charge instead of linear charge density for the 2D(3v) plots in the first column of Fig. 2, every ejected line charge was multiplied by the spotsize of the respective simulation, 2.26 μm and 2.174 μm for the 3×1018Wcm2 (top) and 5.4×1017Wcm2 (bottom), respectively.

FIG. 2.

Ejected electron distributions are presented for both 2D(3v) (left column) and 3D (middle and right column) simulations for 3×1018Wcm2 (upper panels) and 5.4×1017Wcm2 (lower panels) intensities with linear polarization. For each plot, distance from the origin indicates the electron kinetic energy, angle from the origin indicates the angle of ejection from the target, and color represents the charge per angular-energy bin in pC MeV−1 rad−1 units. The 3D simulation results show ejected electrons projected both onto the polarization plane (center column) and the perpendicular (B-k̂) plane (right column). A gray line on each plot shows ±25° of the laser back-reflection axis. Ejected electrons that leave within the conical section defined by these angles are counted in the conversion efficiency measurement (see Table I). Electrons with kinetic energies below 120 keV are not counted in this measurement (as discussed in the text) which is why these lines stop at 120 keV in all plots. To directly compare 2D(3v) and 3D simulations, the charges in the 2D(3v) simulation (left column) were multiplied by the spot-size of the laser.

FIG. 2.

Ejected electron distributions are presented for both 2D(3v) (left column) and 3D (middle and right column) simulations for 3×1018Wcm2 (upper panels) and 5.4×1017Wcm2 (lower panels) intensities with linear polarization. For each plot, distance from the origin indicates the electron kinetic energy, angle from the origin indicates the angle of ejection from the target, and color represents the charge per angular-energy bin in pC MeV−1 rad−1 units. The 3D simulation results show ejected electrons projected both onto the polarization plane (center column) and the perpendicular (B-k̂) plane (right column). A gray line on each plot shows ±25° of the laser back-reflection axis. Ejected electrons that leave within the conical section defined by these angles are counted in the conversion efficiency measurement (see Table I). Electrons with kinetic energies below 120 keV are not counted in this measurement (as discussed in the text) which is why these lines stop at 120 keV in all plots. To directly compare 2D(3v) and 3D simulations, the charges in the 2D(3v) simulation (left column) were multiplied by the spot-size of the laser.

Close modal

Finally, a 3D simulation with the same laser parameters as the 5.4×1017Wcm2 simulations was performed with circular (instead of linear) polarization. The ejected electron results for this simulation are shown in Fig. 3 and compared to the results from the simulation with linear polarization.

FIG. 3.

Ejected electron distributions for linearly polarized and circularly polarized lasers with intensity 5.4×1017Wcm2. Both plots are results of 3D simulations and show projections onto the xy plane, which in the linearly polarized case is the polarization plane. The plot layout and color mapping are the same as for Fig. 2.

FIG. 3.

Ejected electron distributions for linearly polarized and circularly polarized lasers with intensity 5.4×1017Wcm2. Both plots are results of 3D simulations and show projections onto the xy plane, which in the linearly polarized case is the polarization plane. The plot layout and color mapping are the same as for Fig. 2.

Close modal

As discussed in Ref. 2, the final focusing mirror in the experiment, the so-called Off-Axis Parabola (OAP), was used as a Faraday cup to measure the number of electrons ejected from the target within the solid angle subtended by the mirror. In one experimental configuration, a 100 μm thick transparent glass slide was placed between the target and the OAP, allowing the laser pulse to pass through but blocking many of the lower energy electrons from reaching the OAP. As discussed in Ref. 2, the energy threshold for electrons to pass through the slide is estimated to be about 120 keV. Using this value as a minimum energy for the electrons that are still detected at the OAP gives a robust lower limit on the conversion efficiency from laser energy to electrons with kinetic energies above 120 keV. This lower bound conversion efficiency from experiment can be compared to measurements from simulations. In our simulations, the conversion efficiency is inferred in an analogous way to the experiment, counting only the electrons with kinetic energies above 120 keV that are ejected within the angle subtended by the 50.5° opening angle (0.60 sr solid angle) of the OAP. Experimental and simulated conversion efficiencies are summarized in Table I. The experimental results indicate a larger total amount of ejected electrons and a higher conversion efficiency than that observed in these simulations. We discuss possible reasons for this in Sec. IV. For the 3×1018Wcm2 case, the discrepancy is less severe and the conversion efficiency from simulations is more similar to the lower bound on the conversion efficiency inferred from the experiment.

TABLE I.

Summary of the simulation results presented here and experimental results from Morrison et al.2 Conversion efficiencies can be calculated for 2D(3v) PIC simulations but because of the geometry of these simulations an exact value for the charge ejected cannot be determined as discussed in the text. “Charge Accelerated Backwards” is defined as summed charge of electrons ejected from the target with energy >120 keV and within an angle of 25° from the incoming laser axis. “Laser Conversion Efficiency” is defined as the total kinetic energy of these electrons divided by the energy of the incoming laser.

1018Wcm25.4·1017Wcm23·1018Wcm2
IntensityExperiment (linear Pol.)2D linear Pol.3D linear Pol.3D circular Pol.2D linear Pol.3D linear Pol.
Run time … 209 fs 185 fs 
Charge accelerated backwards (>120 keV) ∼300 pC … 20.6 pC 39.2 pC … 88.5 pC 
Conversion efficiency (>120 keV) >1.5% 0.71% 0.56% 1.82% 1.42% 1.01% 
1018Wcm25.4·1017Wcm23·1018Wcm2
IntensityExperiment (linear Pol.)2D linear Pol.3D linear Pol.3D circular Pol.2D linear Pol.3D linear Pol.
Run time … 209 fs 185 fs 
Charge accelerated backwards (>120 keV) ∼300 pC … 20.6 pC 39.2 pC … 88.5 pC 
Conversion efficiency (>120 keV) >1.5% 0.71% 0.56% 1.82% 1.42% 1.01% 

All simulations presented here demonstrate significant numbers of super-ponderomotive electrons ejected backwards in the direction of the reflection of the normally incident pulse, where the ponderomotive potential for the linearly polarized simulations is 32 keV for 5.4×1017Wcm2 and 180 keV for 3×1018Wcm2. This extends and confirms the results of earlier 2D(3v) simulations presented in Ref. 5, and it corroborates recent experimental measurements with an electron spectrometer indicating that significant numbers of MeV-scale electrons are ejected from the target.7 

Table I summarizes the results for the conversion efficiencies. It is interesting to compare the 2D(3v) and 3D simulation results. For the linear polarized simulations, the 2D(3v) and 3D simulations have conversion efficiencies that agree within 30% for both simulated laser intensities. We had expected the 3D simulations to exhibit a significantly lower conversion efficiency from the simple consideration that 3D simulations include a third spatial dimension and electrons are no longer confined to the plane of the laser polarization. Thus, we are currently without a satisfactory explanation for this result. We summarize the qualitative differences between 2D(3v) and 3D simulations in Sec. IV C.

Although the conversion efficiencies are comparable, the angle-energy-spectra show interesting discrepancies between the 2D(3v) simulations and the analogous 3D simulations. For the 5.4×1017Wcm2 simulations, note that the lower row, first and second columns of Fig. 2, shows an increase in the angular spread, manifesting as a large amount of charge beyond the angular shadow of the OAP, while for the 3D simulations, most of the spectrum is focused on ejection angles directly backwards, almost precisely within the opening angle of the OAP. Meanwhile, for 3×1018Wcm2 results, the 2D(3v) simulations show more charges backwards ejected at an angle from 180° above 2 MeV, and almost no charge ejected directly backwards, whereas backwards is the significant ejection angle for charges ejected from the 3D simulations.

Fig. 3 presents 3D simulation results for the circularly polarized and linearly polarized simulations which were performed with the 5.4×1017Wcm2 intensity and the same spot size. There is a stark difference between the circularly polarized and linearly polarized results. For the circular polarized simulation, ejected electrons range from low energies until 3 MeV, while for the linearly polarized simulation, ejected electron energies only reach ∼1 MeV (Fig. 3). Interestingly, the electrons from the circular polarized simulation are preferentially ejected directly backwards from the target, rather than exhibiting preferred angles away from the back direction. This is likely due to the nature of circularly polarized laser pulses, which exhibit no preferred direction for the laser electric fields. By contrast, linearly polarized laser pulses do have a specified polarization angle and, as discussed in Section IV B and in the  Appendix, linearly polarized plane waves will preferentially direct electrons towards specific angles, depending on the electron energy. This does not occur for circularly polarized laser pulses and we are still working to incorporate circularly polarized light into the framework described in Section IV B.

One possible explanation for this is that the reflectivity of circularly polarized light is greater than that of linearly polarized light (both at normal incidence) which could result in back-directed electrons experiencing higher peak electric fields in the circularly polarized case. However, we measured the reflectivity in the 5.4×1017Wcm2 simulations and found that the reflectivity was similar in the circularly polarized and linearly polarized cases. The conversion efficiency and peak energy is very important for applications and we plan to continue to pursue a theoretical explanation for the enhanced electron acceleration in future work.

To help understand the dynamics of electrons in simulation with linearly polarized laser pulses and to interpret earlier results, we developed an analytic model for single particle motion that can provide an estimate for both the energy of the ejected electrons and the direction of ejection. This model significantly improves upon5 which only provided an order-of-magnitude estimate for the ejected electron energy.

The model is as follows: First, consider if an electron is ejected from the target in the region near the critical density by the standing-wave mechanism,5,6 possibly with an extra boost from the electrostatic fields due to the evacuation of charges due to the ionization of the target by the incident pulse. For electrons ejected into the reflected pulse, we can consider, for simplicity, if the electrons then propagate according to the classical motion in a electromagnetic pulsed plane wave in vacuum, ignoring other effects such as quasi-static electric and magnetic fields, and the focusing (or defocusing) of the reflected laser pulse. The motion of a charged particle in a plane-wave in vacuum is a well-known result from Landau and Lifshitz,14 in which electrons move with the quiver velocity in the polarization direction in response to the laser electric fields and in the forward direction in response to the laser magnetic fields and the Lorentz force. The well-known, parametric Landau and Lifshitz result is often employed in the non-relativistic case when a0<1. In the relativistic case in which the particle has an arbitrary initial velocity in any direction, and with an arbitrary time envelope, we derive the momenta and energies of charged particles in the  Appendix cast in a form that is parametric in the instantaneous phase of the laser as observed by the charged particle, η=ωtk·x. We then choose a temporal pulse shape of a sine squared envelope, which matches the temporal envelope of our incident beam in our simulations. An important parameter in the model is the ratio of the pulse period (FWFM) to the laser period, which we denote as α=60fs/2.67fs22.5.

It is well known that there is no net energy gain for an electron in a plane wave in vacuum.15 Moreover, it remains true that an electron overtaken by a time-pulsed plane wave in the absence of focusing will receive no net energy gain. However, if we assume that the injection of the electron into the reflected pulse is due to standing wave fields, a net energy gain will occur because the electron will be launched into the mid-way into the reflected pulse and not the beginning or ahead of the pulse.

The analytic solution from the  Appendix yields expressions for the kinetic energy (γ1)mc2 where γ is given by (A6) and the ejected angle ϕ in the polarization plane using py/px=tanϕ, where the momentum from (A5) is utilized. From these two expressions, we can create parametric plots of energy vs. ejection angle that can be compared with earlier plots from the LSP simulations. Fig. 4 shows this comparison for the 2D(3v) simulations with 5.4×1017Wcm2 (top) and 3×1018Wcm2 (bottom). In the model predictions in Fig. 4, we choose the incident momentum into the reflected pulse as up to a0 in the traverse directions and 2a0 in the backwards (longitudinal) direction, as discussed in the  Appendix. Although it is a parametrized model, there appears to be reasonable choices that produce preferred angles away from the backwards direction that are similar to the preferred angles seen in the simulation results. The model also provides a more-precise explanation for how electrons are accelerated to super-ponderomotive energies than was presented in Ref. 5.

FIG. 4.

Analytic, pulsed plane wave solutions from the  Appendix overlaid with the 2D(3v) simulation angle-energy spectra, the first column of Fig. 2. The curved lines are plots of kinetic energy vs. ejection angle (here, in the polarization plane) with each color representing a particular injection momentum. Here, the injection momentum is scaled to a0. Specifically, each line corresponds to an injection momentum p/mc=γ0β0=a0ρ on injection into the reflected pulse. Each line is labelled by its ρ=(ρx,ρy) value (as in the legend inset). The lines are created by varying the injection phase η0 from πα to 2 πα, corresponding to phase of the laser at the point and time of injection. For more specific details, refer to the  Appendix.

FIG. 4.

Analytic, pulsed plane wave solutions from the  Appendix overlaid with the 2D(3v) simulation angle-energy spectra, the first column of Fig. 2. The curved lines are plots of kinetic energy vs. ejection angle (here, in the polarization plane) with each color representing a particular injection momentum. Here, the injection momentum is scaled to a0. Specifically, each line corresponds to an injection momentum p/mc=γ0β0=a0ρ on injection into the reflected pulse. Each line is labelled by its ρ=(ρx,ρy) value (as in the legend inset). The lines are created by varying the injection phase η0 from πα to 2 πα, corresponding to phase of the laser at the point and time of injection. For more specific details, refer to the  Appendix.

Close modal

Although the plane-wave model gives a possible explanation for some of the features seen in the angular-energy spectra, it does not explain the differences between the 2D(3v) and 3D simulations. Note that the form of Eqs. (A5), (A6) implies that for electrons in a plane wave, the maximum kinetic energy occurs when the transverse momentum is also maximized. Thus, one expects very few energetic electrons to be moving directly backwards (βy=0) according to this model. So while this model may provide some insight into the 2D(3v) simulations, its utility may be limited in interpreting the 3D simulation results, which do not appear to exhibit preferred angles away from the backwards direction.

As mentioned earlier, the differences between 2D(3v) and 3D simulations were more pronounced for the 3×1018Wcm2 case. We do not yet have a satisfying explanation for this, but it is worth outlining the known differences between 2D(3v) and 3D PIC simulations since there is likely multiple causes for the observed difference. Qualitative differences between 2D3(v) and 3D simulations include:

  1. Electro-static fields decay as 1/r in 2D(3v) simulations instead of 1/r2 as would be the case in a truly 3D world.

  2. Particles are confined to the plane of the laser polarization in 2D(3v) instead of being free to move in all three spatial directions.

  3. Laser light comes to focus differently in 2D(3v) than would occur in 3D.

Regarding the first point, Orban et al.5 concluded that charge-separation effects caused by the overlap of the forward-going and reflected laser pulse will produce quasi-static electric fields that can enhance electron acceleration. These quasi-static electric fields develop in a process known as “ponderomotive steepening,”16 which is so named because these fields will ultimately cause a steepening of the ion and electron density profiles. But before this occurs, quasi-static electric fields can provide an extra boost to some electrons that ultimately reach higher energies than they would without this effect.5 Since the quasi-static electric fields decay with a different radial dependence in 3D and 2D(3v) simulations, this will lead to a different acceleration. With this in mind, in future work we plan to measure the precise energy gain from quasi-static electric fields in the simulations from analyzing particle trajectories.

The second point mentioned—namely, the confining of particles to the plane of the polarization—affects the realism of the simulation in a number of ways. For example, the real laser pulse will produce gradients in the electron density both in and out of the plane of the polarization of the laser. There is also the possibility that a charge near the laser axis will experience larger electric and magnetic fields than it would otherwise because it is confined to the plane of polarization and cannot move out of the laser spot in the direction perpendicular to this plane. This latter consideration is another reason for surprise that the conversion efficiencies of the 3D simulations did not turn out to be significantly lower than the 2D(3v) simulations.

For an illustration of the third point, see Fig. 5. This figure shows the spatial dependence of the laser intensity along the laser axis for the laser in vacuum as a function of the distance from peak focus (x = 0). This dependence is different in 2D(3v) simulations than it is in 3D because the 2D(3v) geometry can be thought of as a waveguide. This causes the light to come to focus less rapidly than a 3D laser pulse with the same spot size and intensity at peak focus. As a result, the laser intensity near the critical density is substantially larger in the 2D(3v) simulations than in the 3D simulations. This could artificially enhance the conversion efficiency of 2D(3v) simulations relative to 3D, yet, as mentioned earlier, the conversion efficiencies of 2D(3v) and 3D simulations are overall quite similar. Moreover, as we discuss in the  Appendix, the spatially escaping electrons could be considered as electrons escaping a foreshortened time pulse with a smaller value of α, as long as they remain in the pulse for more two or three laser periods (α>3). For a tighter focus, more electrons will be free to escape the pulse, after which they will no longer propagate according to the pulsed wave model and other effects can become important. Regardless, this remains a puzzling result that we will return to in future work.

FIG. 5.

Identically defined Gaussian lasers evolve to the same peak intensity in different ways between 2D(3v) and 3D Cartesian. In 2D(3v) PIC simulations, a Gaussian laser propagates as a 2D Gaussian “wedge” pulse while it propagates as a true Gaussian laser in 3D PIC simulations. This leads to intensity differences before the focus between 2D(3v) and 3D, even though both pulses reach the same peak intensity. The laser axis intensity of a Gaussian in 2D PIC (blue) and 3D PIC (red) in vacuum is plotted as a function propagation distance along k (x). The vertical dashed line marks the location of the critical density surface in our PIC simulations.

FIG. 5.

Identically defined Gaussian lasers evolve to the same peak intensity in different ways between 2D(3v) and 3D Cartesian. In 2D(3v) PIC simulations, a Gaussian laser propagates as a 2D Gaussian “wedge” pulse while it propagates as a true Gaussian laser in 3D PIC simulations. This leads to intensity differences before the focus between 2D(3v) and 3D, even though both pulses reach the same peak intensity. The laser axis intensity of a Gaussian in 2D PIC (blue) and 3D PIC (red) in vacuum is plotted as a function propagation distance along k (x). The vertical dashed line marks the location of the critical density surface in our PIC simulations.

Close modal

The results of 2D(3v) and 3D PIC simulations using the LSP code were compared in the context of electron acceleration using a high intensity, linearly polarized laser pulse interacting with a water jet target. Orban et al.5 presented the first 2D(3v) PIC simulations of this experiment (Morrison et al.2 and Feister et al.7) and commented on the mechanisms of electron acceleration. In this paper, we present 3D PIC simulations of this experiment for the first time and compare with 2D(3v) PIC simulations. These simulations were performed at intensities 5.4×1017Wcm2 and 3×1018Wcm2 which bracket the experimental intensity (1018 W cm−2).

Comparing 2D(3v) and 3D results, there are some differences in the angular distribution of ejected electrons, but both 2D(3v) and 3D simulations confirm significant back-directed electrons with super-ponderomotive energies, in agreement with the experiment. More quantitatively, while the laser-to-ejected-electron conversion efficiencies from 2D(3v) and 3D simulations with the same intensity and spot size were similar, the precise values were somewhat lower than the estimated lower bound conversion efficiency from experiment of 1.5% for electrons with kinetic energies greater than 120 keV.2 

In an effort to understand the angular and energy distribution of ejected electrons, we developed a parameterized analytic model that considers the dynamics of electrons that are accelerated initially by the standing wave fields and quasi-static electric fields present in the laser-interaction region. Later, these electrons experience the reflected laser pulse which we model as a plane wave. We find that this model can describe a number of features observed in the angular-energy spectra of 2D(3v) simulations.

We also performed a 3D PIC simulation using circular polarized laser light with the same parameters as the 5.4×1017Wcm2 simulation. Remarkably, electron energies were observed up to ∼3 MeV and the conversion efficiency increased to 1.8%. This result indicates that experiments with circularly polarized light should prove to be more effective than the experiments that have been conducted with linearly polarized light. The nature of this enhanced electron acceleration will be considered in future work.

This research was sponsored by the Quantum and Non-Equilibrium Processes Division of the Air Force Office of Scientific Research, under the management of Dr. Enrique Parra, Program Manager. The authors acknowledge significant support from the Department of Defense High Performance Computing Modernization Program (DOD HPCMP) Internship Program. Supercomputer time was used on the DOD HPC Spirit and Garnet supercomputers. Resources were also used at the Ohio Supercomputer Center.

Consider an electron in a linearly polarized plane wave where k=kx̂ and ŷ is the polarization direction, and in which the electric field has the form E(x,t)=E0f(η)ŷ where η=ω(tx/c), ω is angular frequency of the plane wave in the lab frame, and c is the speed of light. The Lorentz force can be written

1ωddt(γβ)=a0f(η)(ŷ(1βx)+βyx̂),
(A1)

where we have expanded the β×ẑ term and identified a0=eE0/mωc, where e is the charge of the electron and m is the electron's mass. a0 is often referred to as the quiver velocity or the normalized vector potential. From the x̂ component of (A1) and the power equation, one can obtain

ddt(γβx)=a0ωf(η)βy=ddt(γβ)·β=ddtγ,
(A2)

which implies that the Doppler shift factor γγβx=γdη/dt is a conserved quantity.

Given that the force depends on the Lorentz invariant, but dynamic phase η, this motivates the change of variable from the frame time to the phase (tη). Noting that dη/dt=ω(1βx) and using chain rule yield the equation of motion

ddη(γβ)=a0f(η)(ŷ+x̂γβyγγβx).
(A3)

We can directly integrate the ŷ component of (A3) and then substitute the result of that integration into the x̂ component, while noting that the denominator is constant due to (A2). To facilitate this derivation, we define the first structure function

f1(η0,η)=η0ηdηf(η),
(A4)

and the second structure function

f2(η0,η)=η0ηdηf(η)f1(η)=df1(η)f1(η)=f1(η0,η0)2/2,

where we simplified the second structure function and defined it in terms of the f1. We thus obtain the following solution for γβ as a function the phase η

γβ(η)=γ0β0ŷa0f1(η0,η)+x̂γγβx(a022f1(η0,η)2γ0βy0a0f1(η0,η)),
(A5)

where γ0 and β0 are the γ-factor and speed at an arbitrary starting phase η0. We will use η0 as the phase at injection into the reflected pulse. For the case of a monochromatic sinusoidal plane wave, we have that f(η)=cosη,η0=0, which implies that f1(η0=0,η)=sinη. Substituting this result into (A5) yields the solution given in Landau and Lifshitz.14 

Given that γ0γ0βx0=γγβ, we obtain the energy for the electron

γ(η)=γ0+1γγβxa022sin2ηγ0βy0γγβxa0sinη.
(A6)

Finally, note that

β=1cdxdt=dkxdη(1βx)
kx=0ηdηγβ(η)γγβx,

which, from (A5), yields the (normalized) position as a function of η

kx(η)=kx0+γ0β0γγβxηŷa0f1i(η0,η)γγβx+x̂(γγβx)2(a02f2i(η0,η)γ0βy0a0f1i(η0,η)),
(A7)

where f1i=dηf1 and f2i=dηf2=12dηf12. Finally, using η=ωtkx, we obtain that ωt(η0,η)=η+kx(η0,η), where we use (A7). This gives all quantities as a function of frame time ωt(η) parametric in η.

Using (A5), we define a pulsed plane wave that has a sine squared shape corresponding to the time envelope of our simulations. To do this, we set

f(η)=sinηsinη2αH(η)H(2παη),
(A8)

where H(x) is the Heaviside-step function where H(x<0)=0 and H(x>0)=1. Here, α is the ratio of the pulse time to the laser period. For our simulations, the pulse is 60 fs long (30 fs full-width at half-maximum), so α=22.48. We omitted a pulse phase shift for simplicity. Using this choice for f(η) yields the following structure functions and integrals

f1(η0,η)=sin[η2α12α]sin[η02α12α](2α1)/αsin[η2α+12α]sin[η02α+12α](2α+1)/α,
(A9)
f2(η0,η)=f1(η0,η)2/2,
(A10)
f1i(η0,η)=cos[η2α12α]cos[η02α12α](2α1)2/2α2+cos[η2α+12α]cos[η02α+12α](2α+1)2/2α2f1(0,η0)(ηη0),
(A11)

and

f2i(η0,η)=α22(2α1)2[ηη02α2sin[η2α1α]sin[η02α1α]2α1]+α22(2α+1)2[ηη02α2sin[η2α+1α]sin[η02α+1α]2α+1]+α24(4α21)[sin2ηsin2η02α(sinηαsinη0α)]+f1(0,η0)[f1i(η0,η)12f1(0,η0)(ηη0)].
(A12)

For Eqs. (A9)–(A12), η0,η is restricted between 0 and 2πα due to the step functions in (A8) that make the pulse finite. Given the mechanics of a standing wave, we expect that an electron injected from the standing wave into the reflected pulse will be injected for a given phase η0 in the middle of the pulse and remain in the pulse until the pulse overtakes the electron. The end of the pulse as observed by the electron corresponds to the phase η=2πα. This leads to a net energy gain for electrons in the reflected pulse. We note that for the case in which a pulse simply passes over an electron in free space, corresponding to η0=0,η=2πα, there is no net energy gain from (A5) except for moments when the electron is in the laser field.15 Instead, superponderomotive energies are achieved because the standing wave injects electrons mid-way into the reflected pulse and with a significant momentum.

If the electron is injected at a phase before the middle of the pulse, i.e., η0<πα, then the symmetry of the pulse in time across the pulse's peak time will average out momentum gains across half the pulse. This warrants parametrizing η0 from πα to 2πα. To create the theoretical model predictions shown in Fig. 4, we take the ratio of γβy/γβx=tanϕ from (A5) with this range of η0 and use (A9) for our model for f1 to obtain the ejection angle of an accelerated electron in the polarization.

The model just described still requires an initial momentum, γ0β0, which we assume is on the order of the normalized electric field, a0. It is well known that a0 is the impulse done during one cycle by the transverse electric field of a simple plane wave. More importantly, a0 is an estimate for the longitudinal momentum from a standing wave as found by Ref. 6. According to Ref. 6, electrons that are ejected longitudinally from a plane wave can have a momentum significantly in excess of a0, sometimes with as much as 2a0. For this reason, Fig. 4 shows this model's prediction for initial momenta scaled to a0. Specifically, we consider initial momenta of the form p/mc=a0ρ, where we vary ρx from 0 to 2 and ρy between −0.5 and 0.5.

We note that a weakness of the model is that for the intense case of 3×1018Wcm2 with large transverse injection momentum (ρy0.5), the model predicts transverse motion that would take the electron outside of the width of the laser pulse. This is not the case for the 5.4×1017Wcm2 intensity. We show predictions that include ρy=0.5 in Fig. 4 for both intensities for completeness even though the ρy=±0.5 cases with 3×1018Wcm2 are suspect. Conceivably, the energies and angles for ρy=±0.5 and 3×1018Wcm2 may still be accurate for the simple reason that our results do not strongly depend on the precise value of α. Electrons that travel outside of the laser width before the end of the pulse are similar to electrons that experience a pulsed plane wave with a smaller value of α. Thus, the model predictions for ρy=±0.5 and 3×1018Wcm2 may still be qualitatively accurate.

1.
A. G.
Mordovanakis
,
J.
Easter
,
N.
Naumova
,
K.
Popov
,
P.-E.
Masson-Laborde
,
B.
Hou
,
I.
Sokolov
,
G.
Mourou
,
I. V.
Glazyrin
,
W.
Rozmus
,
V.
Bychenkov
,
J.
Nees
, and
K.
Krushelnick
,
Phys. Rev. Lett.
103
,
235001
(
2009
).
2.
J. T.
Morrison
,
E. A.
Chowdhury
,
K. D.
Frische
,
S.
Feister
,
V. M.
Ovchinnikov
,
J. A.
Nees
,
C.
Orban
,
R. R.
Freeman
, and
W. M.
Roquemore
,
Phys. Plasmas
22
,
043101
(
2015
).
3.
D.
Panasenko
,
A. J.
Shu
,
A.
Gonsalves
,
K.
Nakamura
,
N. H.
Matlis
,
C.
Toth
, and
W. P.
Leemans
,
J. Appl. Phys.
108
,
044913
(
2010
).
4.
M. C.
Levy
,
S. C.
Wilks
,
M.
Tabak
, and
M. G.
Baring
,
Phys. Plasmas
20
,
103101
(
2013
).
5.
C.
Orban
,
J. T.
Morison
,
E. D.
Chowdhury
,
J. A.
Nees
,
K.
Frische
, and
W. M.
Roquemore
,
Phys. Plasmas
22
,
023110
(
2015
).
6.
A. J.
Kemp
,
Y.
Sentoku
, and
M.
Tabak
,
Phys. Rev. E
79
,
066406
(
2009
).
7.
S.
Feister
,
D. R.
Austin
,
J. T.
Morrison
,
K. D.
Frische
,
C.
Orban
,
G.
Ngirmang
,
E. A.
Chowdhury
,
R. R.
Freeman
, and
W. M.
Roquemore
, “
Super-ponderomotive electron spectra from efficient, high-intensity, kHz laser-water interactions
,” (unpublished).
8.
D. R.
Welch
,
D. V.
Rose
,
R. E.
Clark
,
T. C.
Genoni
, and
T. P.
Hughes
,
Comput. Phys. Commun.
164
,
183
(
2004
).
9.
C.
Birdsall
and
A.
Langdon
,
Plasma Physics Via Computer Simulation
, Series in Plasma Physics (
Taylor & Francis
,
2004
).
10.
J. D.
Huba
, Naval Research Laboratory (
2013
), available at http://www.nrl.navy.mil/ppd/sites/www.nrl.navy.mil.ppd/files/pdfs/NRL_FORMULARY_13.pdf.
11.
V. M.
Ammosov
,
N. B.
Delone
, and
V. P.
Krainov
,
Sov. Phys. JETP
64
,
1191
(
1986
).
12.
A. J.
Kemp
,
R. E. W.
Pfund
, and
J.
Meyer-Ter-Vehn
,
Phys. Plasmas
11
,
5648
(
2004
).
13.
L.
Spitzer
,
Am. J. Phys.
31
,
890
(
1963
).
14.
L. D.
Landau
and
E.
Lifshitz
,
The Theory of Classical Fields
(
Elsevier
,
Oxford
,
1975
).
15.
E.
Esarey
,
P.
Sprangle
, and
J.
Krall
,
Phys. Rev. E
52
,
5443
(
1995
).
16.
K.
Estabrook
and
W. L.
Kruer
,
Phys. Fluids
26
,
1888
(
1983
).