# Numerical modeling of laser tunneling ionization in particle-in-cell codes with a laser envelope model

Francesco Massimo, Arnaud Beck, Julien Dérouillat, Imen Zemzemi, Arnd Specka

► To cite this version:

Francesco Massimo, Arnaud Beck, Julien Dérouillat, Imen Zemzemi, Arnd Specka. Numerical modeling of laser tunneling ionization in particle-in-cell codes with a laser envelope model. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2020, 102 (3), pp.033204. 10.1103/Phys-RevE.102.033204. hal-02959553

HAL Id: hal-02959553

https://hal.science/hal-02959553

Submitted on 9 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical modeling of laser tunneling ionization in particle-in-cell codes with a laser envelope model

F. Massimo,¹,*, A. Beck,¹ J. Derouillat,² I. Zemzemi,¹ and A. Specka¹

¹Laboratoire Leprince-Ringuet-École polytechnique, CNRS-IN2P3, Palaiseau 91128, France

²Maison de la Simulation, CEA, CNRS, Université Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

<figure><img src="image_2.png" /></figure>

(Received 11 June 2020; accepted 19 August 2020; published 14 September 2020)

The resources needed for particle-in-cell simulations of laser wakefield acceleration can be greatly reduced in many cases of interest using an envelope model. However, the inclusion of tunneling ionization in this time-averaged treatment of laser-plasma acceleration is not straightforward, since the statistical features of the electron beams obtained through ionization should ideally be reproduced without resolving the high-frequency laser oscillations. In this context, an extension of an already known envelope ionization procedure is proposed, valid also for laser pulses with higher intensities, which consists in adding the initial longitudinal drift to the newly created electrons within the laser pulse ionizing the medium. The accuracy of the proposed procedure is shown with both linear and circular polarization in a simple benchmark where a nitrogen slab is ionized by a laser pulse and in a more complex benchmark of laser plasma acceleration with ionization injection in the nonlinear regime. With this addition to the envelope ionization algorithm, the main phase space properties of the bunches injected in a plasma wakefield with ionization by a laser (charge, average energy, energy spread, rms sizes, and normalized emittance) can be estimated with accuracy comparable to a nonenvelope simulation with significantly reduced resources, even in cylindrical geometry. Through this extended algorithm, preliminary studies of ionization injection in laser wakefield acceleration can be easily carried out even on a laptop.

DOI: 10.1103/PhysRevE.102.033204

# I. INTRODUCTION

In the last few decades, the limits in accelerating gradients of conventional electron accelerators based on metallic cavities prompted considerable efforts in the development of alternative electron acceleration techniques. Hitherto, the acceleration of electrons in the wake of an intense laser pulse propagating in an underdense plasma [laser wakefield acceleration (LWFA) [1–4]] has been proven particularly promising, generating smaller electron accelerators with high accelerating gradients [5–7], GeV level final energies [8,9], and femtoseconds duration accelerated beams [10]. An important role in this acceleration scheme is played by the technique used to inject relativistic electrons in the accelerating phase of the involved plasma waves [3,4]. Among the numerous demonstrated injection techniques, one of particular simplicity and often chosen is called ionization injection [11–18]. It consists in using a gas mixture to generate the required plasma. This mixture is mainly composed of a low-atomic-number Z gas, like hydrogen or helium, that is ionized very early in the laser-gas interaction. The other part of the mixture is composed by a higher Z dopant gas, like nitrogen, whose first levels of ionization are reached early in the laser-gas interaction. However, the higher ionization levels of this gas can be accessed only at higher values of the laser transverse electric field, normally near the peak of the envelope of the laser, i.e., later than the first levels of the high-Z gas and the ones of the low-Z gas. Tailoring properly the laser and

plasma parameters, these ionization levels are reached only during the short period when the laser is near its maximum focusing. It provides a reserve of electrons that can be stripped off from the high-Z ions near the peak of the laser pulse “just in time” to be trapped in the plasma wave in the wake of the laser (this case is referred to as self-truncated ionization injection [16–19]).

Particle-in-cell (PIC) [20] modeling of this interaction is an essential investigation technique to design and analyze LWFA experiments with ionization injection. The most common technique to take into account tunneling ionization in PIC simulations of these phenomena is to compute at each time step the number of electrons freed from their atoms through the Ammosov-Delone-Krainov direct current (ADK dc) ionization rate [21–23]. When a new electron is created in this way, a sufficiently accurate approximation for momentum conservation is to assign it an initial zero momentum, since the heavy ion from which it is created can be considered immobile. At later times, the electron quickly acquires a transverse momentum whose normalized magnitude is $p_{\perp}(t) = |\mathbf{A}_{\perp}(t) - \mathbf{A}_{\perp}(t_{\text{ioniz}})|$ [24–26] (see Appendix B), where $\mathbf{A}_{\perp}(t)$, $\mathbf{A}_{\perp}(t_{\text{ioniz}})$ are the instantaneous laser transverse vector potential at time t and at the ionization time $t_{\text{ioniz}}$. Given the high-frequency oscillations of the laser, this results in a quiver motion following the laser oscillations. This occurs through the instantaneous interaction force of the electron with the laser pulse, the Lorentz force, which takes into account the high-frequency laser oscillations. In the following, this kind of simulation will be referred to as standard laser simulations, which can be performed in Cartesian geometries or quasicylindrical geometry [27]. This kind of simulation,

<footnote>*Corresponding author: massimo@llr.in2p3.fr</footnote>
which needs to be performed in three dimensions in order to have physical accuracy [28], requires large amounts of resources, due to the disparity between the typical length of a LWFA accelerator, at least of the order of 1 mm, and the smallest scale to resolve, the laser pulse carrier wavelength $\lambda_0$, of the order of 1 µm. The parameter space to explore to understand and design LWFA experiments is vast and cannot be explored directly with numerous standard laser simulations in three dimensions. For these reasons, reduced models for LWFA simulations are of paramount interest for preliminary studies, because they can significantly reduce the computation time at the cost of introducing physical approximations which are reasonable in many LWFA configurations. One of the most general of these reduced models is the use of azimuthal Fourier decomposition or quasicylindrical geometry [27], which takes into account only the first azimuthal modes of the electromagnetic fields, reducing the cost of simulation with quasi-three-dimensional (3D) accuracy to the cost of approximately $N_m$ two-dimensional (2D) PIC simulations on a cylindrical grid, where $N_m$ is the number of considered azimuthal modes. For most preliminary studies of LWFA, $N_m = 2$ is sufficient, while for the reconstruction of more realistic laser pulses a higher number of modes is necessary [29].

Another reduced model of interest for LWFA is given by the envelope or ponderomotive guiding center model. In most LWFA setups, normally the ratio between the scales involved in the excitation of the required plasma waves by the ponderomotive force and the smallest scale to resolve is more than 10. Thus, an averaged formulation of the ponderomotive interaction between the laser pulse and the plasma can be used to obtain accurate results in shorter simulation times, often by one or two orders of magnitude with the same geometry, lifting the requirement of the resolution of the laser wavelength [24,30–37]. As well, preliminary studies of LWFA with this model can be performed in quasicylindrical geometry considering only one azimuthal mode representing perfect cylindrical symmetry [30, 33, 38, 39]. In that particular geometry, the savings in computational resources are even larger because a single azimuthal mode is used.

Modeling tunneling ionization in envelope simulations such as in standard laser simulations, using the ADK dc ionization rate and the zero-momentum initialization for the new electrons, does not yield accurate results, since the high-frequency laser oscillations and thus the electron motion is not well resolved. The residual momentum of the electrons stripped from the atoms or ions after the passage of the laser pulse strongly depends on the extraction field phase, which is normally poorly resolved in an envelope simulation. To circumvent this problem, in the cylindrical envelope code INF&RNO a reconstruction of the high-frequency laser oscillations near the laser pulse is performed at each time step, calculating the full ionization rate and describing the ionization-quiver dynamics of the new electrons [40]. Although very accurate, a disadvantage of this approach is that an additional grid is required to interpolate the force acting on the ions for the ionization and on the new electrons for their quiver motion. However, if the quiver motion of the electrons does not need to be reconstructed, a reduced approach would ideally reconstruct the main integrated parameters of the electron bunches obtained through ionization injection in LWFA,

i.e., charge, average energy, energy spread, rms sizes, and normalized emittance. In particular in Ref. [41] it was shown that the use of averaged ADK ionization rate at each time step of envelope simulations gives a more accurate estimate of the injected charge, and in Ref. [34] a procedure to reconstruct the residual transverse momentum spread in envelope simulations, based on the analytical results in Ref. [42], was outlined. The results in Ref. [34] have been obtained with values of the maximum normalized vector potential of the laser $a_0 = A_{max} < 1$.

For $a_0 > 1$ the same procedure computes a trapped charge in LWFA envelope PIC simulations that is lower than the one computed in standard laser PIC simulations. Therefore, in this work the reason for this disagreement is discussed, i.e., the initialization of the longitudinal momentum of the electrons created by ionization, and an extension of the tunneling ionization modeling procedure described in Ref. [34] is proposed, to obtain accurate results also for $a_0 > 1$. This extended procedure has been implemented in the open source PIC code SMILEI [43,44], used for the simulations of the paper. In an unified framework, SMILEI can run both standard laser and envelope simulations, in Cartesian geometries [36] and quasicylindrical geometry [39], the latter used for the simulations of this paper. The envelope ionization procedure proposed in this paper can be used for Cartesian geometries, but also in purely cylindrical geometry (i.e., with only the azimuthal mode $m = 0$) and at the same time take into account the initial electron momentum asymmetry intrinsic with a linear polarization for the laser, provided that the wakefields present a significant degree of cylindrical symmetry (i.e., the envelope of the laser has cylindrical symmetry). Also, using only one azimuthal mode in an envelope cylindrical simulation yields a speed-up compared to a standard laser simulation in quasicylindrical geometry, where the minimum number of azimuthal modes to use is two [39]. This is of particular interest for LWFA, where 2D Cartesian simulations fail to give accurate results [28] and at least 3D, quasicylindrical simulations [27] or envelope cylindrical simulations are necessary [39]. The proposed envelope ionization procedure cannot yield an accurate description of LWFA setups where carrier-envelope effects play an important role in ionization injection, like in LWFA with few-cycle laser pulses [45,46], where in general the envelope model cannot yield accurate results.

The paper is organized as follows: in the second section, the procedure to model ionization in PIC envelope simulations through the ADK model is presented and its characteristic elements are discussed, i.e., the use of the averaged ionization rate and the initialization of the transverse and longitudinal momentum of the newly created electrons. The authors’ original contribution is included in the initialization of the longitudinal momentum of the new electrons, a key to obtaining accurate results with $a_0 > 1$. In the third section, a basic ionization benchmark case is introduced to compare the results of a standard PIC simulation and the equivalent envelope simulation, for linear and circular polarization. In the fourth section a well-known one-dimensional (1D) model of LWFA is reviewed to hint at the importance of an accurate initialization of the longitudinal momentum of the new electrons like in the proposed ionization procedure. In the
fifth section, the comparison between a standard simulation and an envelope simulation with the proposed ionization technique for a full LWFA benchmark with ionization of N³⁺ is presented. Both benchmarks of the paper are run in regimes with $a_0 > 1$. In the fifth section it is shown how results with sufficient accuracy for preliminary studies can be obtained with the proposed ionization procedure in a very short time. In Appendix A the equations of the envelope model used in the paper simulations are briefly reviewed. In Appendix B the derivation of the initial momentum values assigned to the electrons created with the proposed ionization procedure is described.

## II. TUNNELING IONIZATION ALGORITHM WITH AN ENVELOPE MODEL

In the next subsections, the elements of the envelope ionization procedure implemented in SMILEI are outlined, for both laser and envelope simulations. In the following equations of the paper, unless specified, normalized units will be used. Charges are normalized to the unit charge e, velocities to the speed of light c, masses to the electron mass $m_e$, lengths to the inverse of the laser carrier wave number $\lambda_0/2\pi$, and momenta by $m_e c$.

### A. Tunneling ionization rate with an envelope

Following the notation in Ref. [47], the ADK dc tunneling ionization rate of an atom or ion under the effect of a constant electric field of magnitude $|E|$ is, in atomic units (4.134×10¹⁶ Hz) [21–23,47]:

$$W_{\text{ADK, dc}} = A_{n^*, l^*} B_{l, |m|} I_p \left[ \frac{2(2I_p)^{3/2}}{|E|} \right]^{2n^* - |m| - 1} \times \exp \left[ -\frac{2}{3} \frac{(2I_p)^{3/2}}{|E|} \right], \quad (1)$$

where the coefficients $A_{n^*, l^*}$, $B_{l, |m|}$ are given by

$$A_{n^*, l^*} = \frac{2^{2n^*}}{n^* \Gamma(n^* + l^* + 1) \Gamma(n^* - l^*)}, \\ \times B_{l, |m|} = \frac{(2l + 1)(l + |m|)!}{2^{|m|} |m|! (l - |m|)!} \quad (2)$$

In the previous equations, $I_p$ is the ionization potential for the Z + 1 level of ionization normalized to 27.2116 eV, $|E|$ is normalized to 0.514224 TV/m, $\Gamma(x)$ is the gamma function, $n^* = Z/\sqrt{2I_p}$, $l^* = n^* - 1$, and $l$ and $m$ are the angular momentum and its projection on the laser polarization vector, respectively. In Ref. [22] it is shown that the ionization rate for $m = 0$ is dominant, thus only this one is taken into account in SMILEI.

In a standard laser simulation with single-level ionization, the probability of ionization at each time step for each atom or ion is computed as $1 - \exp(-W_{\text{ADK, dc}} \Delta t)$, where $\Delta t$ is the integration time step of the simulation. The extension to multiple-level ionization is described in Ref. [47]. The use of the ADK dc rate in a standard PIC is physically justified from the fact that normally the field does not vary significantly within an interval $\Delta t$, chosen to resolve the laser oscillation frequency. As discussed in Ref. [41], this approximation is no

longer valid in an envelope simulation, where the integration time step can potentially contain multiple laser oscillations. Thus, as recommended in the same reference, the averaged version of the ionization rate in alternate current $W_{\text{ADK, ac}}$ must be used with an envelope model for accurate results. For circular polarization, since the laser electric field magnitude does not change within a laser oscillation, $W_{\text{ADK, ac}} = W_{\text{ADK, dc}}$. For linear polarization, the averaging of $W_{\text{ADK, dc}}$ yields [22,41]

$$W_{\text{ADK, ac}} = \left[ \frac{3}{\pi} |E_{\text{envelope}}| (2I_p)^{-3/2} \right]^{1/2} W_{\text{ADK, dc}}(|E_{\text{envelope}}|). \quad (3)$$

For the computation of the ac ionization rate, the magnitude of the laser envelope electric field $|E_{\text{envelope}}|$, including the longitudinal and transverse field components (see Appendix A for their computation), must be used instead of the instantaneous field $|E|$ used in Eq. (1). Note that the ponderomotive force does not change the ionization rate [48].

The use of $W_{\text{ADK, ac}}$ in envelope simulations ensures more correct computations of the total amount of electrons created by tunneling ionization [22,34]. However, as discussed in the next subsections, to correctly model LWFA, a correct initialization of the transverse and longitudinal momenta of these electrons also is needed with high laser intensities.

### B. Transverse momentum initialization

In standard laser simulations, to ensure the conservation of momentum, the momentum assigned to the new electrons created by ionization is normally the same of the atom or ion from which they originated. In most LWFA setups, this momentum is initially zero. The new electrons almost instantly acquire a quickly oscillating component in the transverse momentum in their interaction with the laser (see Appendix B), which in general can be nonzero and depends on the position or phase of the electron. Thus, these electrons initially have a certain spread in the transverse momentum, quantified in Ref. [42]. Neglecting the high-frequency oscillations, their averaged dynamics is then determined by the ponderomotive force of the laser and in the case also by the averaged plasma wakefield in LWFA.

The aim of the proposed ionization procedure for an envelope simulation is not necessarily to obtain injected electron bunches identical in the phase space to those in an equivalent standard laser simulation, but to obtain residual statistical properties that are at least similar. Here the term “residual” denotes the values of the physical quantities computed through the averaging over the laser oscillations, but taking into account their initial momentum spread. To satisfy this requirement, in SMILEI the same procedure for transverse momentum initialization of the new electrons in envelope simulations as in the PIC code ALADyn [35,49] and the hybrid PIC-fluid quasistatic code QFluid [34] was implemented, for both linear and circular polarizations. This procedure is briefly reviewed in the following and discussed in detail in Appendix B.

To preliminary define the notation used in this work (the same notation from Ref. [24], used also in the Appendixes), the transverse vector potential of a laser that
can be described by a laser envelope is denoted with $\hat{\mathbf{A}}_{\perp}(\mathbf{x}, t) = \text{Re}[\tilde{\mathbf{A}}_{\perp}(\mathbf{x}, t)e^{i(x-t)}]$, where $\tilde{\mathbf{A}}_{\perp}$ is the complex envelope of the laser transverse vector potential. This envelope, whose evolution is described by the envelope equation (see Appendix A), takes the form $\tilde{\mathbf{A}}_{\perp} = \mathbf{e}_y \tilde{\mathbf{A}}$ for a laser linearly polarized in the y direction and $\tilde{\mathbf{A}}_{\perp} = \frac{\tilde{\mathbf{A}}}{\sqrt{2}}(\mathbf{e}_y \pm i\mathbf{e}_z)$ for a circularly polarized laser. In SMILEI, at a given time step the envelope equation is solved for $\tilde{\mathbf{A}}$ for linear polarization and for $\tilde{\mathbf{A}}/\sqrt{2}$ for circular polarization (see Appendix A).

Electrons stripped from an atom or ion by a laser with linear polarization have an initial Gaussian distribution in the transverse momentum $p_{\text{pol}}$ in the polarization direction, with rms width [42]

$$\sigma_{p_{\text{pol}}} = |\tilde{\mathbf{A}}|(\frac{3}{2}|\tilde{\mathbf{E}}_{\text{envelope}}|)^{1/2}(2I_p)^{-3/4}, \quad (4)$$

where $|\tilde{\mathbf{A}}|$ is the magnitude of the complex envelope of the laser vector potential in the polarization direction. Thus, for linear polarization simulations, the initial transverse momentum in the polarization direction of the new electron $p_{\text{pol},0}$ is assigned as a pseudorandom number drawn from a Gaussian distribution of rms width $\sigma_{p_{\text{pol}}}$, computed from Eq. (4) using the $|\tilde{\mathbf{A}}|$ interpolated at the new electron position (i.e., the same position of the atom/ion from which it originated). Since the laser is propagating in the x direction, the initial momentum in the direction perpendicular to x and to the polarization direction is assigned as zero. Thus, in linear polarization the magnitude of the initial transverse momentum of the new electrons is $|\mathbf{p}_{\perp,0}| = p_{\text{pol},0}$. We remark that, using the notation in Ref. [42], Eq. (4) is valid when

$$\Delta = \left(\frac{3}{2}|\tilde{\mathbf{E}}_{\text{envelope}}|\right)^{1/2}(2I_p)^{-3/4} \ll 1. \quad (5)$$

For example, the normalized $I_p$ corresponding to the ionization of the sixth (seventh) level of nitrogen (used in the benchmarks of this work) are equal to 552 eV/13.6 eV = 40.6 [667 eV/13.6 eV = 49], thus $\Delta \ll 1$ as long as $|\tilde{\mathbf{E}}_{\text{envelope}}| \ll (2/3)(2I_p)^{3/2} = 28$ [34] TV/m, corresponding to values of the normalized maximum vector potential $a_0 \ll 7$ [8.5]. These conditions are satisfied in the benchmarks of this work.

For circular polarization, since the magnitude of the laser transverse vector potential does not change within a time step, the transverse momentum of the new electron is assigned as $|\mathbf{p}_{\perp,0}| = |\hat{\mathbf{p}}_{\perp}| = |\hat{\mathbf{A}}_{\perp}| = |\tilde{\mathbf{A}}|/\sqrt{2}$, with an azimuthal angle randomly extracted from 0 to $2\pi$. Again $|\tilde{\mathbf{A}}|$ is interpolated at the new electron position at the iteration when ionization occurs.

### C. Longitudinal momentum initialization

In standard laser simulations, also the longitudinal momentum assigned to the new electrons created by ionization is typically the same as the atom or ion which they originated from, in general zero for LWFA setups. As in the transverse case, the new electrons almost instantly acquire a quickly oscillating component in the longitudinal momentum in their interaction with the laser (see Appendix B), which depends on the vector potential to which the electron is subject. Considering that the transverse momentum also depends on the vector potential, the two momentum components are not independent. Electrons subject to a higher $|\hat{\mathbf{A}}_{\perp}|$ should have both

higher transverse and longitudinal momentum at the same time. Also, since the residual transverse momentum spread is not zero, the initial longitudinal momentum spread will be in general nonzero. To the authors' knowledge, in the literature there is no analytical result that estimates this spread. Additionally, the calculation of the average residual $p_x$ implies the averaging over the square of a locally sinusoidal function (see Appendix B), hence it can be inferred that in an envelope simulation a residual positive average longitudinal momentum (a drift in the positive x direction) should be present.

For values of $a_0 < 1$, the use of the ADK ac ionization rate and the initialization of the transverse momentum of the new electrons described in the previous section are sufficient to have a statistically accurate description of the electrons obtained from ionization in an envelope ionization [34]. The ionization process $Ar^{8+} \rightarrow Ar^{9+}$ requires $a_0 = 0.4$ for $\lambda_0 = 0.4$ µm, for example, and can be correctly modeled with this approach [34]. However, to ionize some dopants commonly used in LWFA with ionization injection, higher values of $a_0$ are necessary. In the next sections it is shown that in these cases an initialization of the longitudinal momentum of the new electrons also is necessary. For example, the commonly chosen process $N^{5+} \rightarrow N^{6+}$ requires $a_0 > 1.8$ for $\lambda_0 = 0.8$ µm. This process is used in order to benchmark the momentum initialization that is proposed in this work, as shown in the next sections.

The initial $p_{x,0}$ of the new electrons can be chosen to be initialized as

$$p_{x,0} = \begin{cases} \frac{|\tilde{\mathbf{A}}|^2}{4} + \frac{p_{\text{pol},0}^2}{2} & \text{for linear polarization} \\ \frac{|\tilde{\mathbf{A}}|^2}{2} & \text{for circular polarization.} \end{cases} \quad (6)$$

This choice, motivated in Appendix B, yields statistical characteristics in the longitudinal momentum that are similar to those of the electrons created through tunneling ionization within a time step $\Delta t$, as is shown in the following sections. A worthwhile remark is that the maximum values assigned to $p_{x,0}$ scale as $a_0^2/2$ or smaller, thus they start to become relevant only for values of $a_0$ approaching 1. In the following sections it is shown also that for $a_0 > 1$ this choice has effects also on the transverse momentum evolution, due to the contribution of $p_x$ to the relativistic inertia of the electrons, quantified by their Lorentz factor. In LWFA, this choice of $p_{x,0}$ also allows the new electrons to have an averaged momentum in the x direction high enough to be trapped and efficiently accelerated by the plasma wave behind the laser (see Secs. IV and V).

### III. BASIC CASE STUDY: IONIZATION OF NITROGEN SLAB

LWFA with ionization injection with high $a_0$ is rich in nonlinear phenomena, which take place at the same time with the further ionization of the high-Z gas: relativistic plasma wave excitation, relativistic self-focusing [50], and trapping and acceleration of electrons. To show more in detail the comparison in the momentum space obtained between the two kinds of simulations, standard laser PIC simulation and envelope simulation, a benchmark is presented in this section, where the interaction of the new electrons with the laser is limited, no plasma wave is present to accelerate electrons,
<figure><img src="image_3.png" /><figcaption>FIG. 1. Initial state of the N$^{5+}$ ionization benchmark in the envelope simulations along the propagation axis x. The envelope module |$\tilde{A}$| and the charge density of the N$^{5+}$ slab are normalized to 1.</figcaption></figure>

and the laser pulse evolution is negligible compared to that
in vacuum.

This simple benchmark consists in initializing an immobile
cylindrical slab of N⁵⁺ ions and their neutralizing electrons
(those obtained ionizing the first five atomic levels) in vacuum
and let the ions be further ionized by the passage of an intense
Gaussian laser pulse with a₀ = 2 in the linear polarization
simulation and with a₀ = 2√2 in the circular polarization
simulation (see Fig. 1). The nitrogen slab is a cylinder of
radius R = 100 λ₀/2π and length L = 120 λ₀/2π, whose axis
coincides with the laser propagation axis (the x axis). The
longitudinal density profile of the nitrogen slab is a plateau
with density n₀ = 0.0000056 n_c, where n_c is the critical den-
sity for the chosen laser wavelength λ₀, with an upramp and a
downramp of zero length.

The Gaussian laser pulse, propagating along the positive x direction with carrier wavelength $\lambda_0 = 0.8$ µm (as for a Ti:Sa laser system), has a waist size $w_0 = 100 \lambda_0/2\pi$ and full width half maximum duration in field $L_{FWHM} = 92.5 \lambda_0/2\pi c$. The laser pulse initial position at the beginning of the simulation is at a distance $875 \lambda_0/2\pi$ from the nitrogen slab along the x axis, and its focal plane is placed at the beginning of the nitrogen slab. In the linear polarization simulations, the polarization direction of the laser is the y direction.

The laser pulse passes through the N⁵⁺ slab triggering
injection, and the new electrons are left to move subject to
the initial ponderomotive force of the laser and to the electric
field that is created from the charge imbalance progressively
created in the slab after the ionization and the movement of
the electrons.

After a time $T = 1277.2 \lambda_0/2\pi c$, when the laser is far from the N⁵⁺ slab, the momentum distributions of the newly created electrons in the envelope and standard laser simulations are compared. It is important to remark that this comparison can show an agreement only after the interaction with the laser is finished. Indeed, the transverse momentum distribution in the standard laser simulation when the laser is still in the nitrogen slab would show two peaks due to the quivering

motion, which of course would be completely absent in the
envelope simulation. In other words, the momenta in the enve-
lope simulations represent the slowly varying part **p** of the real
momenta **p** = **p**̃ + **p**̇ (see Appendix A). After the interaction
with the laser, the quickly oscillating part of the momenta **p**̇
is negligible, and the momenta of the envelope simulations **p**̃
can be compared to the total momenta of the standard laser
simulation **p**. For this reason, the bar over the momenta in the
envelope simulations results is dropped in the following for
the sake of brevity.

Both simulation types, standard laser and envelope, are
run in quasicylindrical geometry [27], with transverse reso-
lution Δr = 2 λ₀/2π. The standard laser simulations model
the laser-plasma interaction using two azimuthal modes
(modes m = 0, 1), while the envelope simulations does it
with only one azimuthal mode (m = 0, representing perfect
cylindrical symmetry). Although in the simulation with lin-
ear polarization the cylindrical symmetry is not respected
in the momentum initialization of the envelope simulation,
the resulting asymmetries do not seem to lead to signifi-
cant statistical differences in the momentum distribution of
the electrons, since the magnitude of the currents and field
with m = 1 is negligible compared to those with m = 0.
The standard laser simulations have a longitudinal resolution
Δxlaser = 0.125 λ₀/2π and integration time step Δt lase r =
0.124 λ₀/2πc, while in the envelope simulations these are
Δx envelope = 0.8 λ₀/2π and Δt envelope = 0.79825Δx envelope /c.
With this choice of integration time step, the envelope solver
is still stable [36] and Δt envelope = (515/10), thus the results
of the two kinds of simulation can be easily compared at the
same time, provided that the ratio between the iterations is the
same for the standard laser and envelope simulations.

The laser envelope is initialized already in the simulation
domain, while the laser enters from the left border of the
simulation domain in the standard laser simulation, through
Silver-Müller boundary conditions [51,52]. Despite this dif-
ference, it was ensured that the laser focal plane position
and the initial distance between the laser temporal center and
the N⁵⁺ slab were the same for envelope and standard laser
simulations.

The total number of macroparticles per cell is the same
in all simulations, Nₓ × Nᵣ × Nθ = 32, placed regularly along
the three directions. In the laser simulation, the distribution
of the macroparticles is [Nₓ, Nᵣ, Nθ] = [1, 4, 8], where Nₓ
and Nᵣ are the number of particles in the x and r directions
and Nθ is the number of particles evenly distributed in the
azimuthal angle interval between 0 and 2π. For the envelope
simulations, the macroparticle distribution is [Nₓ, Nᵣ, Nθ] =
[8, 4, 1]. Since only one azimuthal mode is used in the enve-
lope simulations, i.e., the physical phenomena are assumed
to be cylindrically symmetric, in the envelope simulations
Nθ = 1 can be used.

The N⁵⁺ ions are immobile, i.e., although they project
charge on the grid and the electromagnetic field and the laser
envelope are interpolated at their positions for the ionization
operations, their positions and momenta are not changed. This
is a reasonable approximation in most LWFA regimes, since
the timescales of ion motion are much longer than those of
electron motion, and avoids adding a significant computation
time for operations which do not impact the final results.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Comparison of (a) the momentum distributions and (b) the position distributions of the electrons created by ionization in the N⁵⁺ benchmark at time T = 1277.2 λ₀/2πc, with linear polarization along the y axis.</figcaption></figure>

To show the effects of the $p_x$ initialization for the new electrons in the envelope simulations presented in Sec. II C, in the top panels of Figs. 2 and 3 the comparison of the electron momentum distributions at time T (i.e., 2000 iterations for the envelope simulations, 10300 iterations for the standard laser simulations) is reported. The figures also show the results obtained without longitudinal momentum initialization, with linear polarization along the y axis and circular polarization. Without this initialization, the longitudinal momentum distribution results are completely different in the two simulations, and the transverse momenta distributions display differences as well, although the shape of these distributions is the same as with longitudinal momentum distribution. The electrons in the envelope simulation without $p_x$ initialization seem to be much slower than in the standard laser simulations, pushed on the left by the field created by the charge separation that progressively forms and the laser ponderomotive force. For LWFA, this difference is particularly critical because if the electrons obtained from the dopant do not possess an averaged longitudinal momentum high enough to be trapped in the plasma wave in the wake of the laser, the injected charge will be significantly lower in the envelope simulation (see Secs. IV and V). Thus, as will be shown in the next sections, agreement can be expected in the injected charge with the proposed $p_x$ initialization procedure in a LWFA simulation with $a_0 > 1$. Although the transverse momenta initialization is the same for the two envelope simulations (the one described in Sec. II B), the resulting distribution of the transverse momenta at time T is quantitatively different, because with the $p_x$ initialization a

portion of macroparticles starts to become nearly relativistic ($1 + |\mathbf{p}|^2 = 1 + p_x^2 + |\mathbf{p}_\perp|^2 \gtrsim 2$), changing their inertia towards the ponderomotive force and the force created by the charge separation.

Even in this simple benchmark the difference in momenta at the moment of creation of the new electrons in the envelope simulations results in completely different position distributions at later times. In the bottom panel of Figs. 2 and 3 the comparison of the electron position distribution at time T is reported. The positions are referred to the center of the N⁵⁺ slab. As was evident from the top panels, in the envelope simulation without $p_x$ initialization the drift in the x direction is lower, thus the macroparticles on average are moving towards the negative x direction, while in the standard laser simulation and in the envelope simulation with $p_x$ initialization they are moving in the opposite direction. The transverse position distributions are different as well.

#### IV. ELECTRON TRAPPING FROM LWFA WITH IONIZATION INJECTION

Before showing the results of an envelope LWFA simulation, a brief review of a well-known 1D model of the dynamics of an electron in a plasma wave [3,14,53,54] can hint at the importance of an accurate $p_x$ initialization. In this model a plasma of constant density $n_0$ is modeled as a cold relativistic fluid with immobile ions. The driver of the plasma oscillations is a linearly polarized laser pulse with carrier wavelength $\lambda_0$, in the quasistatic approximation [30]
<figure><img src="image_5.png" /><figcaption>FIG. 3. Comparison of (a) the momentum distributions and (b) the position distributions of the electrons created by ionization in the N$^{5+}$ benchmark at time $T = 1277.2 \lambda_0/2\pi c$, with circular polarization.</figcaption></figure>

described by a transverse vector potential $A(\xi = x - v_g t) = a_0 \exp(-\xi^2/\sigma) \cos(\xi)$, where $\sigma = 0.5(L_{FWHM})^2/\ln 2$ ($L_{FWHM}$ is the laser pulse FWHM duration in field). The plasma wave phase velocity is assumed to coincide with the driver laser pulse group velocity in the plasma $v_g$. The laser pulse excites a wakefield with electrostatic potential $\Psi$ and longitudinal electric field $E_x = -\partial_\xi \Psi$ described by the differential equation [55–57]

$$ \partial_{\xi}^{2} \Psi = \left\{ \beta_{p} \left[ 1 - \frac{(1+A^{2})}{\gamma_{p}^{2}(1+\Psi)^{2}} \right]^{-1/2} - 1 \right\}, \quad (7) $$

where $\gamma_p = (1 - \beta_p^2)^{-1/2} = (n_0)^{-1/2}$ is the Lorentz factor associated to the plasma wave phase or group normalized velocity $\beta_p$. Once Eq. (7) is solved assuming $\Psi|_{\xi \to +\infty} = \partial_\xi \Psi|_{\xi \to +\infty} = 0$, through a second-type generating function a conserved Hamiltonian function $H$ can be defined to describe the dynamics of a test electron in this system [53,54]:

$$ \begin{aligned} H &= \sqrt{1 + p_x^2 + p_{\perp}^2} - \Psi - \beta_p p_x \\ &= \sqrt{1 + p_x^2 + |A - A(\xi_0)|^2} - \Psi - \beta_p p_x. \end{aligned} \quad (8) $$

The last identity was obtained using conservation of transverse canonical momentum $p_\perp - A = p_\perp(\xi_0) - A(\xi_0)$ (since the Hamiltonian is not explicitly dependent on the transverse coordinate $x_\perp$; see also Appendix B) and $p_\perp(\xi_0) = 0$. For an electron created at rest through ionization, $\xi_0$ is the point where it is released from its atom or ion and $A(\xi_0) = A_{\text{ioniz}} \to p_\perp = A - A_{\text{ioniz}}$. Free background electrons at rest far from

the laser pulse, at $\xi_0 \to +\infty$, have an initial condition $A(\xi_0) = 0 \to p_\perp = A$. The curves of constant $H$ in the phase space $\xi - p_\perp$ describe the evolution of the test electron momentum and give insightful information on its trapping state. Given a value $H_0$ of the Hamiltonian, the evolution of the test electron momentum for each value of $\xi$ can be found inverting Eq. (8). For a background electron, i.e., $A(\xi_0 \to +\infty) = 0 \to p_\perp = A$:

$$ p_x = \beta_p \gamma_p^2 (H_0 + \Psi) \pm \gamma_p \sqrt{\gamma_p^2 (H_0 + \Psi)^2 - (1 + A^2)^2}. \quad (9) $$

Electrons with $p_x|_{\xi \to +\infty} = 0$ are associated to fluid orbits with $H_0 = H_{\text{fluid}} = 1$, describing an electron indefinitely drifting in the negative $\xi$ direction and oscillating with the plasma wave (blue line in Fig. 4). In this sense, the fluid orbit is an untrapped orbit. A trajectory which remains in the plasma wave bucket behind the laser pulse is referred to as a trapped orbit, which characterizes electrons that can be accelerated by the plasma wave. In this system a separatrix curve (red line in Fig. 4) with $H_{\text{sep}} = \sqrt{1 + \min(A^2)/\gamma_p - \min(\Psi)}$ separates the untrapped ($H_0 > H_{\text{sep}}$) and trapped orbits ($H_0 < H_{\text{sep}}$).

This simple model can give some insight on the ionization injection process, as discussed thoroughly in Refs. [14,54]. If a test electron is considered as initially at rest and stripped through ionization from a high-Z gas atom or ion within the laser pulse at position $\xi_0$, in particular near a peak of the electric field [where $A(\xi_0) \approx 0 \to p_\perp = A$ since $p_\perp - A = \text{const}$], the same Hamiltonian of Eq. (8) can be used to
<figure><img src="image_6.png" /><figcaption>FIG. 4. Adapted from Figs. 2 and 5 of Ref. [54] with the almost the same laser and plasma parameters: $a_0 = 2$, $n_0 = 0.44\% n_c$, $L_{FWHM} = 44 \lambda_0 / 2\pi$. (a) Phase space trajectories of test electrons, including a fluid orbit, the separatrix, and the trajectory of an electron stripped from its atom or ion through ionization; (b) laser vector potential A, electrostatic potential $\Psi$, and longitudinal electric field $E_x$.</figcaption></figure>

describe its motion, with value $H_e = 1 - \Psi(\xi_0)$. The ionization injection scheme of LWFA relies on creating through ionization enough electrons in trapped orbits ($H_e < H_{sep}$) to accelerate them. Figure 4 shows the trajectory (cyan line) of such a test electron created through ionization ($\xi_0 \approx 0$) in a trapped orbit of a nonlinear wakefield ($a_0 = 2$). Although its initial longitudinal momentum is zero, as discussed in the previous section it acquires a momentum under the effect of the laser, and it is trapped in the plasma wave then to be accelerated. In a standard laser simulation of this phenomenon it is accurate to set the $p_{x,0}$ of these electrons as zero, since most of them are created near the peaks of the electric field, where the vector potential is near zero [54]. However, although the electron $p_x$ starts oscillating, the average value of these oscillations is nonzero. Thus, in an envelope simulation where the real laser field varies significantly in an integration time step it is necessary to initialize a nonzero averaged $p_{x,0}$ to be physically accurate. An averaged $p_{x,0} = 0$ would make the electron move outside the separatrix curve, preventing its trapping.

Additionally, as was shown in the previous section, in a multidimensional system the correct initialization of $p_{x,0}$ with $a_0 > 1$ ensures an accurate evaluation of the electron averaged Lorentz factor after its interaction with the laser, which determines its inertia towards the laser ponderomotive force, the

forces present in the plasma, and the correct description of the transverse coordinates' evolution.

V. BENCHMARK CASE STUDY: LASER WAKEFIELD ACCELERATION WITH IONIZATION INJECTION

In this section, a full LWFA simulation with ionization injection is presented as benchmark. A Gaussian laser pulse of wavelength $\lambda_0 = 0.8$ µm, propagating in the positive x direction and linearly polarized in the y direction with $a_0 = 2.5$, is focused at the beginning of a target made by an already ionized mixture of 99% He and 1% N$^{5+}$, exciting a wakefield and further ionizing the nitrogen. Some electrons extracted from the last two ionization levels of nitrogen are then trapped and accelerated in the laser wakefield. The results of two envelope simulations are shown, simulations whose only difference is that one of them uses the $p_x$ initialization described in Sec. II C. They are benchmarked against a standard laser simulation. The ADK ac ionization rate and the initialization of the transverse momentum of the electrons created by ionization in the envelope simulations are those reported in Secs. II A and II B. As in the simulations of Sec. III, the ions are immobile in all simulations.

The laser pulse has a waist $w_0 = 18.7$ µm and a full width half maximum duration in intensity $L_{FWHM} = 33$ fs. The N$^{5+}$
<table><caption>TABLE I. Numerical parameters for the LWFA benchmark simulations.</caption><thead><tr><th></th><th>Standard laser simulation</th><th>Envelope simulations</th></tr></thead><tbody><tr><td>Δx [λ<sub>0</sub>/2π]</td><td>0.125</td><td>1</td></tr><tr><td>cΔt/Δx</td><td>0.993</td><td>0.82749</td></tr><tr><td>N<sub>modes</sub></td><td>2</td><td>1</td></tr><tr><td>N<sub>cells,x</sub></td><td>3328</td><td>416</td></tr><tr><td>N<sub>macroparticles per cell</sub> [N<sub>x</sub>, N<sub>r</sub>, N<sub>θ</sub>]</td><td>1,4,8</td><td>8,4,1</td></tr><tr><td></td><td colspan="2">All simulations</td></tr><tr><td>L<sub>x</sub> = N<sub>cells,x</sub>Δx [λ<sub>0</sub>/2π]</td><td colspan="2">416</td></tr><tr><td>Δr [λ<sub>0</sub>/2π]</td><td colspan="2">2</td></tr><tr><td>N<sub>cells,r</sub> (half plane)</td><td colspan="2">192</td></tr><tr><td>L<sub>r</sub> = N<sub>cells,r</sub>Δr [λ<sub>0</sub>/2π]</td><td colspan="2">384</td></tr></tbody></table>

target longitudinal profile density is $n_0 = 3.4 \times 10^{18} \text{ cm}^{-3}$, with a linear upramp of length $100 \lambda_0/2\pi$. The target has a radius $R = 360 \lambda_0/2\pi$ and constant radial density profile. The laser focal plane is placed at the beginning of the target.

For the reader's convenience, the numerical parameters of the laser and envelope simulations are reported in Table I. As in the previous section benchmark, all the simulations have the same transverse resolution. The moving window physical size is the same, but $\Delta x_{envelope} = 8 \Delta x_{laser}$ and $\Delta t_{envelope} = 0.82749 \Delta x_{envelope}/c$, respectively. The time-step choice ensures that $\Delta t_{envelope}/\Delta t_{laser} \approx 20/3$, allowing one to easily compare the results at approximately the same time, provided that the number of iterations of the two kind of simulations has the same ratio. The same macroparticle distribution used for the benchmark in Sec. III was chosen. Since $\Delta x_{envelope} = 8 \Delta x_{laser}$, the standard laser and envelope simulations have the same spatial sampling $N_x \Delta x \times N_r \Delta r$ in the x-r plane.

In Fig. 5 the evolution in time of the total charge of the electrons created by ionization and with $p_x > 50 m_e c$ is compared. A very good agreement between the standard laser and envelope simulation with $p_x$ initialization is found. In the envelope simulation without $p_x$ initialization, the charge is significantly lower (23 pC instead of $\approx 225$ pC at 800 $\mu$m, for example). As expected, in this simulation the electrons created through ionization do not have an average longitudinal momentum high enough to be trapped in the wakefield behind the laser. In other words, referring, for example, to Fig. 4 and the treatment in Refs. [14,54], their initial phase and average momentum in the wakefield lie within the region of untrapped orbits in the phase space, i.e., outside the separatrix curve. Thus, in an envelope simulation, to correctly model these conditions, the averaged initial momentum of the new-born electrons must be properly initialized, otherwise too many electrons would start from a point in the phase space outside the separatrix curve and would not be trapped. The evolution of the average energy of the same electrons is reported in the bottom panel of Fig. 4. A good agreement is found between the standard laser and envelope simulation with $p_x$ initialization. The average accelerating gradient and thus the final energy in the envelope simulation without the $p_x$ initialization are higher. With lower charge trapped in the wakefield bucket behind the laser, the beam-loading effect [58] on the longitudinal electric field in this simulation is significantly lower and a higher $E_x$ field accelerates the macroparticles. To better show this

phenomenon, in Fig. 6 the longitudinal electric field and the electron charge density on the propagation axis are reported for the three simulations after 800 $\mu$m of propagation. The beam loading of the longitudinal electric field and the density perturbation in the bunch zone of the standard laser simulation are well reproduced by the envelope simulation with $p_x$ initialization, ensuring a more accurate estimate of the final energy ($\approx 90$ MeV), while without the $p_x$ initialization its value is overestimated ($\approx 120$ MeV).

In the standard laser and envelope simulation with $p_x$ initialization sudden changes of the charge and energy evolution curve slopes at $\approx 550 \mu$m and $\approx 600 \mu$m of propagation occur (see bottom panels of Fig. 5). These changes are caused by a complex interplay of self-focusing, injection of a second bunch immediately after the first one (see right panel of Fig. 6), and consequent change in beam loading. The study of the charge and energy evolution of this particular setup is

<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Comparison of the evolution of the laser peak transverse electric field in the LWFA benchmark simulations. For the other two panels, all the electrons created by ionization, present in the moving window and with longitudinal momentum $p_x > 50 m_e c$, are considered. The electrons injected in the plasma wave bucket behind the laser are a subset of these electrons. (b) Comparison of the electron total charge evolution; (c) comparison of the average electron energy.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison of (a) the longitudinal electric field and (b) the electron normalized electron density on the propagation axis after 800 µm, computed with the standard laser and envelope simulations in the LWFA benchmark. The maximum value shown for the charge density is 0.01, to highlight the zone near the injected electron beam.</figcaption></figure>

beyond the scope of this work. The important aspect for its scope is that, with the $p_x$ initialization, the envelope simulation can reproduce this behavior present in the standard laser simulation.

As can be seen in Fig. 5, the peak laser electric field value is increasing due to relativistic-self focusing, arriving even at values of $5 m_e\omega_0c/e$. This shows that the accuracy of the ionization procedure with $p_x$ initialization seems robust even at these values of the ionizing field.

In Fig. 7 a snapshot of the electron density after 800 µm of propagation is shown, comparing the results of the standard laser and envelope simulation with $p_x$ initialization. Although more accurate comparisons of the electron density should be done in one dimension (as in Fig. 6, right panel) without the saturation of a color map in two dimensions, the shape of the injected bunch appears very similar, apart from a fish-bone-like shape in the standard laser simulation. This phenomenon

could be caused by the nonlinear mapping in the phase space of the electron discrete injection by laser ionization described in Ref. [59]. Further investigation, however, is necessary to verify this hypothesis. If the mechanism in Ref. [59] is indeed the cause of this bunched structure, therefore it cannot be correctly described by an envelope simulation with our model, since the injection does not take place in a discrete way near the peaks of the laser field. However, the bunch parameters do not seem too different between the two simulations (see Table II), hence this phenomenon does not appear to have significant consequences at this distance.

In Fig. 8 the energy spectrum of the electrons created by ionization with $p_x > 50 m_e c$ after 800 µm of propagation is reported for the standard laser and envelope simulations. A very good agreement with the laser simulation is found between the laser and envelope simulation with $p_x$ initialization. Without the $p_x$ initialization, as previously discussed, the

<figure><img src="image_9.png" /><figcaption>FIG. 7. Comparison of the electron normalized density on the x-y plane after 800 µm of propagation computed with (a) the envelope simulation and (b) the standard laser simulation in the LWFA benchmark.</figcaption></figure>

<table><caption>TABLE II. Electron beam parameters at 800 µm of propagation: charge $Q$, rms sizes $2\sigma_i$ ($i = x, y, z$), normalized emittances $\epsilon_{n,i}$ ($i = y, z$), mean energy $E$, rms energy spread $\sigma_E/E$. First column: beam parameters at the beginning of the simulation. Second and third columns: beam parameters after 800 µm of propagation in the standard laser simulation and envelope simulation. The electron beam is defined as all the electrons within $2.5\Delta E_{rms}$ of the spectrum peak (see Fig. 8), where $\Delta E_{rms}$ is the rms width of the electrons within a FWHM width in the spectrum around the energy peak.</caption><thead><tr><th></th><th>Standard laser</th><th>Envelope with p<sub>x</sub> initialization</th></tr></thead><tbody><tr><td>Q [pC]</td><td>175</td><td>182</td></tr><tr><td>2σ<sub>x</sub> [µm]</td><td>3.4</td><td>3.5</td></tr><tr><td>2σ<sub>y</sub> [µm]</td><td>2.3</td><td>2.3</td></tr><tr><td>2σ<sub>z</sub> [µm]</td><td>1.1</td><td>1.1</td></tr><tr><td>ε<sub>n,y</sub> [mm-mrad]</td><td>3.9</td><td>4.0</td></tr><tr><td>ε<sub>n,z</sub> [mm-mrad]</td><td>1.2</td><td>1.2</td></tr><tr><td>E<sub>avg</sub> [MeV]</td><td>90.2</td><td>89.7</td></tr><tr><td>σ<sub>E</sub>/E [rms, %]</td><td>11.91</td><td>11.93</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Comparison of the energy spectrum of all the electrons created by ionization in the moving window after 800 µm of propagation computed with the standard laser and envelope simulations in the LWFA benchmark.</figcaption></figure>

trapped charge is underestimated and the average energy is overestimated.

Table II summarizes the bunch parameters at the same distance. The bunch is defined as all the electrons within $2.5\Delta E_{rms}$ of the spectrum peak, where $\Delta E_{rms}$ is the rms width of the electrons within a FWHM width in the spectrum around the energy peak. Note that the bunch electrons constitute only a subset of the electrons with $p_x > 50 \text{ mec}$ considered in Fig. 5.

It is important to highlight that the envelope simulation with $p_x$ initialization needed a significantly smaller amount of resources to run (102 cpu-h) compared to the standard laser simulation (9.3 kcpu-h), a factor 91 of difference. This speed-up comes from considering only one azimuthal mode and using larger $\Delta x$ and $\Delta t$ in the envelope simulation.

Thus, with the envelope ionization technique proposed in this work, preliminary envelope LWFA simulations with ionization injection, even with relativistic values for $a_0$, can become affordable with a small cluster. For even quicker simulations, a hybrid fluid-PIC approach [33, 35, 38, 60] can be envisaged in weakly nonlinear regimes, modeling the immobile ions with macroparticles, the background electrons as a relativistic cold fluid, and the electrons created with ionization with macroparticles. In the next section it is shown how the envelope simulation can yield results that are quantitatively accurate enough for preliminary studies with even less resources.

## VI. EFFECTS OF REDUCING THE NUMBER OF PARTICLES

In the previous section it was shown that using the same spatial sampling by the macroparticles in the x-r plane and the same number of macroparticles per cell a LWFA envelope simulation with $p_x$ initialization can yield results that are very similar to those of a standard laser simulation. This agreement

was obtained at a propagation distance of 800 µm, which is sufficient for a preliminary parametric study, e.g., to design the laser and plasma parameters of an experiment. Since the parameter space to explore is vast, having an estimate of the injected charge and the bunch energy with quick preliminary simulations can greatly speed up the design process, and this is exactly the purpose of reduced models. More cumbersome simulations with nonreduced models can then investigate further a region of parameters of interest found with a coarse study made with reduced models. In this section it is shown that even degrading the accuracy of the envelope simulation by reducing the number of macroparticles the accuracy of the results remains acceptable, especially considering the reduction of the needed computing resources. The possibility to use only one azimuthal mode and the absence of high-frequency oscillations significantly relaxes the sampling requirements in a cylindrical envelope simulation.

The case study for the “degraded” envelope simulations is the same LWFA simulation of the previous section, with the same physical and numerical parameters, except for the distribution of macroparticles per cell. The results of two envelope simulations are reported, with the regular macroparticle distributions [$N_x$, $N_r$, $N_\theta$] = [4, 2, 1], [1, 1, 1] respectively. Therefore, these “degraded” envelope simulations have, respectively, four times and 32 times less macroparticles per cell than the envelope simulation of the previous section. In the following figures and table the results are compared to the standard laser simulation of the previous section, which are taken as a reference.

Figure 9 compares the evolution in time of the peak transverse electric field and of the total charge in the moving

<figure><img src="image_11.png" /><figcaption>FIG. 9. (a) Comparison of the evolution of the laser peak transverse electric field in the LWFA benchmark simulations. For the other two panels, all the electrons created by ionization, present in the moving window and with longitudinal momentum $p_x > 50 \text{ mec}$, are considered. The electrons injected in the plasma wave bucket behind the laser are a subset of these electrons. (b) Comparison of the electron total charge evolution; (c) comparison of the average electron energy. The results of the standard laser simulation and of the envelope simulations with $[N_x, N_r, N_\theta] = [4, 2, 1]$, [1, 1, 1] are reported.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Comparison of (a) the longitudinal electric field and (b) the electron normalized density on the propagation axis after 800 µm, computed with the standard laser and envelope simulations in the LWFA benchmark. The maximum value shown for the charge density is 0.01, to highlight the zone near the injected electron beam.</figcaption></figure>

window computed from the electrons created by ionization with $p_x > 50 m_ec$ in the standard laser simulation and envelope simulations. The evolution of the average energy of the same electrons is reported as well. Reducing the number of particles per cell does not seem to significantly influence the integrated charge and average energy of the considered electron population.

Figure 10 reports a comparison of the longitudinal electric field and of the electron charge density on axis computed with the standard laser and envelope simulations. As expected, reducing the number of particles increases the noise in these grid quantities, but the zone near the injected electron beam, where physical phenomena of interest for LWFA occur, displays a high degree of agreement with the standard laser simulation.

Figure 11 compares the spectrum of the electrons obtained with the standard laser simulation and with the envelope simulations. Again, reducing the number of particles increases the level of noise, but the main features of the spectrum are well reproduced even with only one particle per cell. To delve into the details of the accelerated bunch, Table III compares the beam parameters of the envelope simulations

and of the standard laser simulation. Even in the most noisy envelope simulation the beam parameters have a high degree of agreement with the standard laser simulation. The computing time needed for this envelope simulation ($[N_x, N_r, N_\theta] = [1, 1, 1]$) was 35 minutes, without using MPI or OpenMP. Therefore, with the ionization procedure presented in this work applied to the cylindrical geometry, envelope simulations with results reasonably accurate for preliminary studies can be carried out even from a laptop.

VII. CONCLUSIONS

In the context of an existing ionization algorithm for LWFA simulations with an envelope model, an extension of this algorithm was presented, showing a good agreement with standard laser LWFA simulations also for $a_0 > 1$. This feature proves useful for simulations involving high-Z dopant gases like nitrogen where the last ionization levels are accessed with lasers driving highly nonlinear wakefields. The introduced feature of the proposed algorithm is the initialization of the longitudinal momentum $p_x$ of the new electrons created by ionization,

<table><caption>TABLE III. Comparison of standard laser simulation and envelope simulations with $p_x$ initialization and with decreasing number of macroparticles per cell at 800 µm of propagation. The electron beam parameters: charge $Q$, rms sizes $2\sigma_i$ ($i = x, y, z$), normalized emittances $\epsilon_{n,i}$ ($i = y, z$), mean energy $E$, and rms energy spread $\sigma_E/E$. First column: beam parameters at the beginning of the simulation. Second and third columns: beam parameters after 800 µm of propagation in the standard laser simulation and envelope simulation. The electron beam is defined as all the electrons within $2.5\Delta E_{rms}$ of the spectrum peak (see Fig. 8), where $\Delta E_{rms}$ is the rms width of the electrons within a FWHM width in the spectrum around the energy peak.</caption><thead><tr><th></th><th>Standard laser</th><th>Envelope [$N_x, N_r, N_\theta$] = [4, 2, 1]</th><th>Envelope [$N_x, N_r, N_\theta$] = [1, 1, 1]</th></tr></thead><tbody><tr><td>$Q$ [pC]</td><td>175</td><td>182</td><td>179</td></tr><tr><td>2$\sigma_x$ [µm]</td><td>3.4</td><td>3.5</td><td>3.5</td></tr><tr><td>2$\sigma_y$ [µm]</td><td>2.3</td><td>2.3</td><td>2.4</td></tr><tr><td>2$\sigma_z$ [µm]</td><td>1.1</td><td>1.1</td><td>1.2</td></tr><tr><td>$\epsilon_{n,y}$ [mm-mrad]</td><td>3.9</td><td>4.0</td><td>4.0</td></tr><tr><td>$\epsilon_{n,z}$ [mm-mrad]</td><td>1.2</td><td>1.1</td><td>1.2</td></tr><tr><td>$E_{avg}$ [MeV]</td><td>90.2</td><td>89.6</td><td>89.6</td></tr><tr><td>$\sigma_E/E$ [rms, %]</td><td>11.91</td><td>11.95</td><td>11.52</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Comparison of the energy spectrum of all the electrons created by ionization in the moving window after 800 µm of propagation computed with the standard laser and envelope simulations with $[N_x, N_r, N_\theta] = [4, 2, 1], [1, 1, 1]$ in the LWFA benchmark.</figcaption></figure>

reproducing the initial electron drift in the $x$ direction, which becomes significant with high-intensity lasers. This extended ionization procedure does not need an additional finer grid to reproduce the quiver motion of electrons and can be used in Cartesian and cylindrical geometries as well. Two benchmarks have been presented: the ionization of a slab of $N^{5+}$ and a full LWFA simulation with ionization injection with plasma containing $N^{5+}$ as dopant. In the first benchmark it was shown how this $p_x$ initialization allows us to accurately reproduce the momenta and evolved positions of the electrons created by ionization, with both linear and circular polarization. In the nonlinear LWFA benchmark the $p_x$ initialization has been proven to be essential to accurately compute the total charge trapped in the wakefield bucket behind the laser driver and the other statistical parameters of the trapped electron bunch like average energy, emittance, energy spread, energy spectrum, and rms sizes. It was shown that an estimate of these parameters sufficiently accurate for preliminary studies can be obtained even reducing the number of macroparticles per cell to one, a case in which the LWFA benchmark can run on a laptop in less than 1 h. Considering the degree of agreement obtained for these parameters after 800 µm of propagation with significantly smaller amounts of resources, the proposed extended envelope ionization algorithm can pave the way for quick preliminary studies of LWFA with ionization injection. Future studies should check the effects of the envelope ionization procedure on the accelerated bunches at longer distances of propagation, where long-term 3D effects due to the asymmetry of the linear polarization also can become significant.

ACKNOWLEDGMENTS

F.M. was supported by P2IO LabEx (ANR-10-LABX-0038) in the framework “Investissements d’Avenir” (ANR-11-IDEX-0003-01) managed by the Agence Nationale de la

Recherche (ANR, France) as well as by the Laboratoire de physique des deux infinis Irène Joliot-Curie (IJCLab IN3P3-PALLAS project). The authors are grateful to Kevin Cassou and the IJCLab team for fruitful discussions. Computer time has been granted by GENCI (project 2019-A0050510062) and by the École polytechnique through the LLR-LSI project. The authors are grateful to the TGCC and CINES engineers for their support. The authors thank the engineers of the LLR HPC clusters for resources and help. The authors are grateful to D. Terzani and Gilles Maynard for fruitful discussions.

APPENDIX A: EQUATIONS OF THE ENVELOPE MODEL IN SMILEI

Under the hypothesis of a laser pulse propagating in the positive $x$ direction, with complex envelope of the transverse component of the vector potential $\tilde{\mathbf{A}}_\perp(\mathbf{x}, t)$ slowly varying along the $x$ and transverse directions (e.g., a Gaussian laser pulse with large waist and long duration) compared to the laser wavelength $\lambda_0$, a perturbative treatment can be formulated. This derivation leads to averaged equations which describe some phenomena in the laser-plasma interaction in terms of $\tilde{\mathbf{A}}_\perp$. In many cases of interest for LWFA, this reduced formulation of laser-plasma interaction can accurately describe the relevant physical phenomena that are involved, e.g., self-focusing of the laser pulse, self-injection of electrons, and the radiation pressure or ponderomotive force of the laser acting on the particles. The details of this perturbative treatment leading to the following equations can be found in many references, e.g., Refs. [24,30,35,61]. In this theoretical framework, upon which the envelope simulation model of this work is based, the physical quantities as macroparticles momenta $\mathbf{p}$ can be written as a quickly oscillating part $\hat{\mathbf{p}}$ (denoted with a hat) added to a slowly varying part $\tilde{\mathbf{p}}$ (denoted with a bar). The slowly varying positions $\tilde{\mathbf{x}}$, momenta $\tilde{\mathbf{p}}$, and the so-called ponderomotive Lorentz factor $\tilde{\gamma}$ of the electrons in an envelope simulation follow the ponderomotive equations of motion [24,35]:

$$ \frac{d\tilde{\mathbf{x}}}{dt} = \frac{\tilde{\mathbf{p}}}{\tilde{\gamma}}, \quad \frac{d\tilde{\mathbf{p}}}{dt} = -\left(\tilde{\mathbf{E}} + \frac{\tilde{\mathbf{p}}}{\tilde{\gamma}} \times \tilde{\mathbf{B}}\right) - \frac{1}{2\tilde{\gamma}} \nabla\Phi, \\ \tilde{\gamma} = \sqrt{1 + |\tilde{\mathbf{p}}|^2 + \Phi}, \qquad (A1) $$

where $\Phi = |\tilde{\mathbf{A}}_\perp|^2/2$ is the ponderomotive potential.

From d'Alembert's inhomogeneous wave equation for $\tilde{\mathbf{A}}_\perp = \text{Re}[\tilde{\mathbf{A}}_\perp e^{i(x-t)}]$, the envelope evolution equation in non-comoving coordinates can be derived [35]:

$$ [\nabla^2 + 2i(\partial_x + \partial_t) - \partial_t^2]\tilde{\mathbf{A}}_\perp = \chi \tilde{\mathbf{A}}_\perp, \\ \chi = \sum_{p=1}^{N_{\text{macroparticles}}} \frac{S(\tilde{\mathbf{x}} - \tilde{\mathbf{x}}_p)}{\tilde{\gamma}_p}, \qquad (A2) $$

where $\tilde{\mathbf{x}}_p$ is the particle $p$ slowly varying part of the position, $S(\tilde{\mathbf{x}})$ the shape function of the particle, and $\tilde{\gamma}_p$ its ponderomotive Lorentz factor. In Refs. [24,30–32,37] a reduced version of Eq. (A2) is used, but in SMILEI the full form of Eq. (A2) is solved as in Ref. [35].

It is worth noting that, once the susceptibility $\chi$ is known at a given time step, the envelope equation is linear with respect
to $\tilde{\mathbf{A}}_{\perp}$. Thus at each time step, instead of a vector equation, only a scalar equation can be solved explicitly for a laser of given polarization, determined by the parameter $\varepsilon$:

$$ \tilde{\mathbf{A}}_{\perp} = \mathbf{e}_y \varepsilon \tilde{\mathbf{A}} + i \mathbf{e}_z (1 - \varepsilon^2)^{1/2} \tilde{\mathbf{A}}, \quad (A3) $$

where $\varepsilon = 0, 1$ for linear polarization or $\varepsilon = \pm 1/\sqrt{2}$ for circular polarization. The envelope equation can be solved for a nonzero component of the vector potential at each time step, provided that the ponderomotive potential is defined as $\Phi = |\tilde{\mathbf{A}}_{\perp}|^2/2 = |\tilde{\mathbf{A}}|^2/2$ for all polarizations. Thus, for example, the envelope equation could be solved for $\tilde{\mathbf{A}}/\sqrt{2}$ in circular polarization and for $\tilde{\mathbf{A}}$ in linear polarization (and storing these values would be practical for an easy comparison with a standard laser simulation), but the ponderomotive potential would still be $\Phi = |\tilde{\mathbf{A}}_{\perp}|^2/2 = |\tilde{\mathbf{A}}|^2/2$ for both polarizations. This method is used in SMILEI to facilitate the comparison with the laser electric field computed with standard laser simulations. The numerical schemes and the steps used in SMILEI to solve Eqs. (A1) and (A2) are detailed in Refs. [35,36,39].

To easily compare the electric field components between envelope simulations and standard laser simulations and in the specific context of this work to correctly compute the total electric field of the laser envelope that might trigger injection (see Sec. II A), it is necessary to express the envelopes of the electric field components in terms of the envelope of the transverse vector potential $\tilde{\mathbf{A}}_{\perp}$. Without loss of generality, here they are derived in case of linear polarization along the y axis. Since the real field is related to the envelope field through the identity $\hat{\mathbf{A}}_y = \text{Re}[\tilde{\mathbf{A}}_y e^{i(x-t)}]$, the y component of the electric field can be written, similarly:

$$ \hat{\mathbf{E}}_y = -\partial_t \hat{\mathbf{A}}_y = \text{Re} \left[ -(\partial_t - i) \tilde{\mathbf{A}}_y e^{i(x-t)} \right] = \text{Re} \left[ \tilde{\mathbf{E}}_y e^{i(x-t)} \right]. \quad (A4) $$

Thus, the envelope of the electric field along y can be defined as

$$ \tilde{\mathbf{E}}_y = -(\partial_t - i)\tilde{\mathbf{A}}_y. \quad (A5) $$

For the envelope of the electric field longitudinal component $\tilde{\mathbf{E}}_x$, it is useful to change variables, using $(\xi = x - t, \tau = t) \rightarrow (\partial x = \partial_{\xi}, \partial_t = \partial_{\tau} - \partial_{\xi})$. In the context of the envelope perturbative treatment presented in Refs. [24,30,35,61], $\partial_{\xi} \gg \partial_{\tau}$, thus, using also the Coulomb gauge $\partial_x \hat{\mathbf{A}}_x + \partial_y \hat{\mathbf{A}}_y = \partial_{\xi} \hat{\mathbf{A}}_x + \partial_y \hat{\mathbf{A}}_y = 0$, the longitudinal component of the electric field can be found:

$$ \begin{aligned} \hat{\mathbf{E}}_x = \text{Re}[\tilde{\mathbf{E}}_x e^{i(x-t)}] &= -\partial_t \hat{\mathbf{A}}_x = \partial_t \hat{\mathbf{A}}_x + \partial_{\xi} \hat{\mathbf{A}}_x \approx \partial_{\xi} \hat{\mathbf{A}}_x \\ &= -\partial_y \hat{\mathbf{A}}_y = \text{Re} [(-\partial_y \tilde{\mathbf{A}}_y) e^{i(x-t)}]. \end{aligned} \quad (A6) $$

Thus, the envelope of the longitudinal electric field can be defined as

$$ \tilde{\mathbf{E}}_x = -\partial_y \tilde{\mathbf{A}}_y. \quad (A7) $$

Once the envelopes of the electric field components are computed, the total envelope field defined as $|\tilde{\mathbf{E}}_{\text{envelope}}| = \sqrt{|\tilde{\mathbf{E}}_x|^2 + |\tilde{\mathbf{E}}_y|^2}$ can be used in the calculation of the ADK ac ionization rate in Eq. (3).

## APPENDIX B: INITIAL MOMENTUM OF THE ELECTRONS FROM ENVELOPE IONIZATION

The derivation of the initial momentum for an electron obtained from envelope ionization is discussed. The macroparticles representing such electrons are initialized in an envelope PIC simulation with this momentum, which then evolves under the effect of the ponderomotive force and of the averaged electromagnetic fields. The following derivation assumes that the mentioned electrons are extracted when the force due to the laser field, locally a plane wave, constitutes the dominating force acting on them. Indeed, near the peak of the laser pulse the amplitude of the electromagnetic field due to the wake charge separation is much lower than the amplitude of the transverse electromagnetic field of the laser.

Given a test electron in a plane wave, the relations between the electron momentum and the laser vector potential can be analytically found and are briefly reviewed in this Appendix, following the derivations in Refs. [25,26]. For coherence the hat and bar notation of the previous Appendix for the quickly and slowly oscillating quantities are maintained. The ansatz of the following derivation is a plane wave propagating along the positive x direction, i.e., described by a vector potential with zero component on the x direction and form $\mathbf{A}_{\perp}(x-t)$. From this laser vector potential, the laser transverse electromagnetic fields can be expressed as $\mathbf{E}_{\perp} = -\partial_t \mathbf{A}_{\perp}$, $\mathbf{B}_{\perp} = \mathbf{e}_x \times \partial_x \mathbf{A}_{\perp}$. The momentum evolution equation for an electron in this plane wave can be rewritten as

$$ \frac{d\mathbf{p}_{\perp}}{dt} = -\mathbf{E}_{\perp} - (\mathbf{v} \times \mathbf{B}_{\perp})_{\perp} = (\partial_t + v_x \partial_x)\mathbf{A}_{\perp} = \frac{d\mathbf{A}_{\perp}}{dt}, \quad (B1) $$

which implies $\mathbf{p}_{\perp} - \mathbf{A}_{\perp} = \text{const}$. In LWFA the atom or ion from which the electrons are extracted through tunneling ionization can be considered at rest, hence for conservation of momentum the extracted electrons can be considered at rest ($p_{x,t_{\text{ioniz}}} = |\mathbf{p}_{\perp},t_{\text{ioniz}}| = 0$). Therefore, from Eq. (B1) it can be inferred that

$$ \begin{aligned} \mathbf{p}_{\perp} - \mathbf{A}_{\perp} &= \mathbf{p}_{\perp,t_{\text{ioniz}}} - \mathbf{A}_{\perp,t_{\text{ioniz}}} \\ &= -\mathbf{A}_{\perp,t_{\text{ioniz}}} \rightarrow \mathbf{p}_{\perp} = \mathbf{A}_{\perp} - \mathbf{A}_{\perp,t_{\text{ioniz}}}, \end{aligned} \quad (B2) $$

where $\mathbf{A}_{\perp,t_{\text{ioniz}}}$ is the transverse vector potential acting on the electron initially, at the ionization time $t_{\text{ioniz}}$. Similarly, the equation for the variation of the electron longitudinal momentum and energy can be rewritten as

$$ \frac{dp_x}{dt} = -(\mathbf{v} \times \mathbf{B}_{\perp})_x = -(v_y \partial_x A_z + v_z \partial_x A_y), \quad (B3) $$

$$ \frac{dy}{dt} = -(\mathbf{v} \cdot \mathbf{E}_{\perp}) = (v_y \partial_t A_y + v_z \partial_t A_z). \quad (B4) $$

Subtraction of the previous equations and the use of $(\partial_t - \partial_x)\mathbf{A}_{\perp} = 0$ yields $\frac{d(p_x - \gamma)}{dt} = 0$, implying $p_x - \gamma = \text{const}$. Coherently with the hypothesis of electron initially at rest ($p_{x,t_{\text{ioniz}}} = 0, \gamma_{t_{\text{ioniz}}} = 1$), and using the identity $\gamma^2 = 1 + p_x^2 + |\mathbf{p}_{\perp}|^2$, a quadratic relation between the longitudinal momentum and the vector potential amplitude can be found:

$$ p_x = \frac{|\mathbf{p}_{\perp}|^2}{2} = \frac{1}{2} |\mathbf{A}_{\perp} - \mathbf{A}_{\perp,t_{\text{ioniz}}}|^2. \quad (B5) $$

An envelope model describes the dynamics of the macroparticles using their averaged positions and momenta,
thus an average over the optical cycles is then computed on the momentum components represented by Eqs. (B2) and (B5):

$$ \bar{\mathbf{p}}_{\perp} = (\overline{\mathbf{A}_{\perp} - \mathbf{A}_{\perp, t_{\text{ioniz}}}}) = -\mathbf{A}_{\perp, t_{\text{ioniz}}}, \quad (B6) $$

$$ \begin{aligned} \bar{p}_x &= \frac{1}{2} (|\mathbf{A}_{\perp}|^2 + |\mathbf{A}_{\perp, t_{\text{ioniz}}}|^2 - 2 \mathbf{A}_{\perp} \cdot \mathbf{A}_{\perp, t_{\text{ioniz}}}) \\ &= \frac{1}{2} (|\mathbf{A}_{\perp}|^2 + |\mathbf{A}_{\perp, t_{\text{ioniz}}}|^2) = \frac{1}{2} |\bar{\mathbf{A}}_{\perp}|^2 + \frac{1}{2} |\mathbf{A}_{\perp, t_{\text{ioniz}}}|^2. \end{aligned} \quad (B7) $$

To derive the last two identities the constancy of $\mathbf{A}_{\perp, t_{\text{ioniz}}}$ for a given electron over the optical cycles was used, and a plane wave of the form $\mathbf{A}_{\perp} = \bar{\mathbf{A}}_{\perp} = \text{Re}[\bar{\mathbf{A}}_{\perp} e^{i(x-t)}]$ with $\bar{\mathbf{A}}_{\perp}$ described by Eq. (A3) was assumed.

Using the same plane wave definition, from Eqs. (B6) and (B7) the initial momentum conditions for an electron created through envelope ionization described in Secs. II B and II C can be motivated.

In linear polarization, the rms thermal spread $\sigma_{p_{\perp}}$ of the transverse momentum is due to the extraction of electrons at different phases of the laser, which yield different $\mathbf{A}_{\perp, t_{\text{ioniz}}}$. This thermal spread is quantified by Eq. (4), derived in Ref. [42]. Indeed, Eq. (4) is derived averaging the square of $\bar{\mathbf{p}}_{\perp}$ (which is obtained through an average over the optical cycles denoted with the bar) over the possible ionization phases, an operation which will be denoted as $\langle |\bar{\mathbf{p}}_{\perp}|^2 \rangle = \langle |\bar{\mathbf{p}}_{\perp}|^2 \rangle_{\text{ioniz phases}}$. With this notation, $\sigma_{p_{\perp}} = ((|\bar{\mathbf{p}}_{\perp}|^2))^{1/2} = ((|\mathbf{A}_{\perp, t_{\text{ioniz}}}|^2))^{1/2} = \Delta|\bar{\mathbf{A}}|$ [Eq. (9) of Ref. [42]], where $\Delta$ is defined in Eq. (5) and Ref. [42]. In an envelope simulation, the typical resolution does not allow us to resolve the fast phase oscillations of $\mathbf{A}_{\perp, t_{\text{ioniz}}}$ in linear polarization. Thus, to recreate this rms thermal spread, the initial transverse momentum of electrons in the polarization direction is initialized with a number $p_{\text{pol},0}$ drawn from a centered Gaussian distribution with the mentioned rms spread $\sigma_{p_{\perp}}$, i.e., $N(0, \sigma_{p_{\perp}})$. For circular polarization, the amplitude of the wave does not change with the phase oscillations, thus the value $|\mathbf{A}_{\perp, t_{\text{ioniz}}}| = |\bar{\mathbf{A}}|/\sqrt{2}$ can be used for the initial electron transverse momentum, randomly assigning the transverse direction.

The quadratic dependence of the longitudinal momentum with respect to the transverse momentum, which is an oscillating function, suggests that the initial longitudinal momentum assigned to the electrons created by envelope ionization $\bar{p}_{x,0}$ should be a non-negative value. For circular polarization, the

plane wave definition that was chosen gives, for Eq. (B7), $\bar{p}_{x,0} = |\bar{\mathbf{A}}|^2/2$. From Eq. (B7) it can be inferred that for linear polarization, the average $\langle \bar{p}_x \rangle$ over the possible ionization phases of electrons subject to a vector potential with envelope $\bar{\mathbf{A}}$ is dominated by $|\bar{\mathbf{A}}|^2/4$, i.e., $(1/2)|\bar{\mathbf{A}}_{\perp}|^2$ with $|\bar{\mathbf{A}}_{\perp}| = |\bar{\mathbf{A}}| \sin(x-t)$. This represents the first term of the last identity of Eq. (B7). Indeed, the average over the ionization phases of the second term gives $(1/2)\sigma_{p_{\perp}}^2$, which scales as $\frac{1}{2}[(1.5)^{1/2}|\bar{\mathbf{A}}|^{3/2}(2I_p)^{-3/4}]^2$. For $|\bar{\mathbf{A}}| = 2.5$ and the ionization of the sixth level of nitrogen ($I_p = 40.6$), its ratio over the first term is 0.01. To assign an initial $p_{x,0}$ to an electron created by envelope ionization in linear polarization, in the procedure presented in this work this second small term is neglected. Instead, coherently with the identity $\bar{p}_x = (1/2)|\bar{\mathbf{p}}_{\perp}|^2$, it is reasonable to assign a greater longitudinal momentum to electrons with a greater transverse momentum, thus a term $(1/2)p_{\text{pol},0}^2$ was chosen to be added to $|\bar{\mathbf{A}}|^2/4$ in the initial longitudinal momentum. This choice, although heuristic, also allows to give a stochastic spread in the initial longitudinal momentum of the new electrons, for which no analytical results about its distribution are available at the present moment to the authors' knowledge. Further developments of this procedure include analytical calculations that will allow to initialize the longitudinal momenta with a more formally derived stochastic component.

To resume, dropping the bar notation, in the envelope ionization procedure of this work the magnitude of the initial transverse momentum of the new electrons $|\mathbf{p}_{\perp,0}|$, subject to an envelope $\bar{\mathbf{A}}$, is assigned as in Ref. [34]:

$$ |\mathbf{p}_{\perp,0}| = \begin{cases} p_{\text{pol},0} \leftarrow N(0, \sigma_{p_{\perp}}) & \text{for linear polarization} \\ |\bar{\mathbf{A}}|/\sqrt{2} & \text{for circular polarization.} \end{cases} \quad (B8) $$

The direction of the initial momentum is the polarization direction for linear polarization and randomly assigned between 0 and $2\pi$ for circular polarization.

The initial longitudinal momentum of the new electron is assigned with the following choice:

$$ p_{x,0} = \begin{cases} |\bar{\mathbf{A}}|^2/4 + p_{\text{pol},0}/2 & \text{for linear polarization} \\ |\bar{\mathbf{A}}|^2/2 & \text{for circular polarization.} \end{cases} \quad (B9) $$

[1] T. Tajima and J. M. Dawson, Laser Electron Accelerator, Phys. Rev. Lett. 43, 267 (1979).

[2] V. Malka, S. Fritzler, E. Lefebvre, M.-M. Aleonard, F. Burgy, J.-P. Chambaret, J.-F. Chemin, K. Krushelnick, G. Malka, S. P. D. Mangles et al., Electron acceleration by a wake field forced by an intense ultrashort laser pulse, Science 298, 1596 (2002).

[3] E. Esarey, C. B. Schroeder, and W. P. Leemans, Physics of laser-driven plasma-based electron accelerators, Rev. Mod. Phys. 81, 1229 (2009).

[4] V. Malka, Laser plasma accelerators, Phys. Plasmas 19, 055501 (2012).

[5] S. P. D. Mangles, C. D. Murphy, Z. Najmudin, A. G. R. Thomas, J. L. Collier, A. E. Dangor, E. J. Divall, P. S. Foster,

J. G. Gallacher, C. J. Hooker et al., Monoenergetic beams of relativistic electrons from intense laser-plasma interactions, Nature (London) 431, 535 (2004).

[6] J. Faure, Y. Glinec, A. Pukhov, S. Kiselev, S. Gordienko, E. Lefebvre, J. P. Rousseau, F. Burgy, and V. Malka, A laser-plasma accelerator producing monoenergetic electron beams, Nature (London) 431, 541 (2004).

[7] C. G. R. Geddes, C. Toth, J. van Tilborg, E. Esarey, C. B. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, and W. P. Leemans, High-quality electron beams from a laser wakefield accelerator using plasma-channel guiding, Nature (London) 431, 538 (2004).

[8] W. P. Leemans, A. J. Gonsalves, H.-S. Mao, K. Nakamura, C. Benedetti, C. B. Schroeder, C. Tóth, J. Daniels, D. E.
Mittelberger, S. S. Bulanov et al., Multi-GeV Electron Beams from Capillary-Discharge-Guided Subpetawatt Laser Pulses in the Self-Trapping Regime, Phys. Rev. Lett. **113**, 245002 (2014).

[9] A. J. Gonsalves, K. Nakamura, J. Daniels, C. Benedetti, C. Pieronek, T. C. H. de Raadt, S. Steinke, J. H. Bin, S. S. Bulanov, J. van Tilborg et al., Petawatt Laser Guiding and Electron Beam Acceleration to 8 GeV in a Laser-Heated Capillary Discharge Waveguide, Phys. Rev. Lett. **122**, 084801 (2019).

[10] O. Lundh, J. Lim, C. Rechatin, L. Ammoura, A. Ben-Ismail, X. Davoine, G. Gallot, J-P. Goddet, E. Lefebvre, V. Malka, and J. Faure, Few femtosecond, few kiloampere electron bunch produced by a laser-plasma accelerator, Nat. Phys. **7**, 219 (2011).

[11] C. McGuffey, A. G. R. Thomas, W. Schumaker, T. Matsuoka, V. Chvykov, F. J. Dollar, G. Kalintchenko, V. Yanovsky, A. Maksimchuk, K. Krushelnick et al., Ionization Induced Trapping in a Laser Wakefield Accelerator, Phys. Rev. Lett. **104**, 025004 (2010).

[12] A. Pak, K. A. Marsh, S. F. Martins, W. Lu, W. B. Mori, and C. Joshi, Injection and Trapping of Tunnel-Ionized Electrons into Laser-Produced Wakes, Phys. Rev. Lett. **104**, 025003 (2010).

[13] B. B. Pollock, C. E. Clayton, J. E. Ralph, F. Albert, A. Davidson, L. Divol, C. Filip, S. H. Glenzer, K. Herpoldt, W. Lu et al., Demonstration of a narrow Energy Spread, 0.5 GeV Electron Beam from a Two-Stage Laser Wakefield Accelerator, Phys. Rev. Lett. **107**, 045001 (2011).

[14] M. Chen, E. Esarey, C. B. Schroeder, C. G. R. Geddes, and W. P. Leemans, Theory of ionization-induced trapping in laser-plasma accelerators, Phys. Plasmas **19**, 033101 (2012).

[15] G. Golovin, S. Chen, N. Powers, C. Liu, S. Banerjee, J. Zhang, M. Zeng, Z. Sheng, and D. Umstadter, Tunable monoenergetic electron beams from independently controllable laser-wakefield acceleration and injection, Phys. Rev. Accel. Beams **18**, 011301 (2015).

[16] M. Mirzaie, S. Li, M. Zeng, N. A. M. Hafz, M. Chen, G. Y. Li, Q. J. Zhu, H. Liao, T. Sokollik, F. Liu et al., Demonstration of self-truncated ionization injection for GeV electron beams, Sci. Rep. **5**, 14659 (2015).

[17] J. P. Couperus, R. Pausch, A. Köhler, O. Zarini, J. M. Krämer, M. Garten, A. Huebl, R. Gebhardt, U. Helbig, S. Bock et al., Demonstration of a beam loaded nanocoulomb-class laser wakefield accelerator, Nat. Commun. **8**, 487 (2017).

[18] A. Irman, J. P. Couperus, A. Debus, A. Köhler, J. M. Krämer, R. Pausch, O. Zarini, and U. Schramm, Improved performance of laser wakefield acceleration by tailored self-truncated ionization injection, Plasma Phys. Controlled Fusion **60**, 044015 (2018).

[19] M. Zeng, M. Chen, Z.-M. Sheng, W. B. Mori, and J. Zhang, Self-truncated ionization injection and consequent monoenergetic electron bunches in laser wakefield acceleration, Phys. Plasmas **21**, 030701 (2014).

[20] C. K. Birdsall and A. B. Langdon, *Plasma Physics via Computer Simulation* (Taylor and Francis, Philadelphia, 2004).

[21] A. M. Perelomov, V. S. Popov, and M. V. Terent'ev, Ionization of atoms in an alternating electric field, Sov. J. Exp. Theor. Phys. **23**, 924 (1966).

[22] M. V. Ammosov, N. B. Delone, and V. P. Krainov, Tunnel ionization of complex atoms and atomic ions in electromagnetic field, in *High Intensity Laser Processes*, Vol. 0664, edited

by John A. Alcock (SPIE, Bellingham, Washington, 1986), pp. 138–141.

[23] G. L. Yudin and M. Yu. Ivanov, Nonadiabatic tunnel ionization: Looking inside a laser cycle, Phys. Rev. A **64**, 013409 (2001).

[24] B. M. Cowan, D. L. Bruhwiler, E. Cormier-Michel, E. Esarey, C. G. R. Geddes, P. Messmer, and K. M. Paul, Characteristics of an envelope model for laser-plasma accelerator simulation, J. Comput. Phys. **230**, 61 (2011).

[25] P. Gibbon, *Short Pulse Laser Interactions with Matter* (World Scientific, Singapore, 2005).

[26] A. Macchi, *A Superintense Laser-Plasma Interaction Theory Primer* (Springer, Dordrecht, Netherlands, 2013).

[27] A. Lifschitz, X. Davoine, E. Lefebvre, J. Faure, C. Rechatin, and V. Malka, Particle-in-Cell modelling of laser-plasma interaction using Fourier decomposition, J. Comput. Phys. **228**, 1803 (2008).

[28] X. Davoine, E. Lefebvre, J. Faure, C. Rechatin, A. Lifschitz, and V. Malka, Simulation of quasimonoenergetic electron beams produced by colliding pulse wakefield acceleration, Phys. Plasmas **15**, 113102 (2008).

[29] I. Zemzemi, F. Massimo, and A. Beck, Azimuthal decomposition study of a realistic laser profile for efficient modeling of Laser Wakefield Acceleration, in *Proceedings of the 4th European Advanced Accelerator Concepts Workshop* (Elba Island, 2019).

[30] P. Mora and T. M. Antonsen, Jr., Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas, Phys. Plasmas **4**, 217 (1997).

[31] D. F. Gordon, W. B. Mori, and T. M. Antonsen, A ponderomotive guiding center particle-in-cell code for efficient modeling of laser-plasma interactions, IEEE Trans. Plasma Sci. **28**, 1135 (2000).

[32] C. Huang, V. K. Decyk, C. Ren, M. Zhou, W. Lu, W. B. Mori, J. H. Cooley, T. M. Antonsen, and T. Katsouleas, QUICKPIC: A highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas, J. Comput. Phys. **217**, 658 (2006).

[33] C. Benedetti, C. B. Schroeder, E. Esarey, C. G. R. Geddes, and W. P. Leemans, Efficient modeling of laser-plasma accelerators with INF&RNO, AIP Conf. Proc. **1299**, 250 (2010).

[34] P. Tomassini, S. D. Nicola, L. Labate, P. Londrillo, R. Fedele, D. Terzani, and L. A. Gizzi, The resonant multi-pulse ionization injection, Phys. Plasmas **24**, 103120 (2017).

[35] D. Terzani and P. Londrillo, A fast and accurate numerical implementation of the envelope model for laser-plasma dynamics, Comput. Phys. Commun. **242**, 49 (2019).

[36] F. Massimo, A. Beck, J. Derouillat, M. Grech, M. Lobet, F. Pérez, I. Zemzemi, and A. Specka, Efficient start-to-end 3D envelope modeling for two-stage laser wakefield acceleration experiments, Plasma Phys. Controlled Fusion **61**, 124001 (2019).

[37] T. Silva, A. Helm, J. Vieira, R. Fonseca, and L. O. Silva, On the use of the envelope model for down-ramp injection in laser-plasma accelerators, Plasma Phys. Controlled Fusion **62**, 024001 (2019).

[38] P. Tomassini and A. R. Rossi, Matching strategies for a plasma booster, Plasma Phys. Controlled Fusion **58**, 034001 (2015).

[39] F. Massimo, I. Zemzemi, A. Beck, J. Dérouillat, and A. Specka, Efficient cylindrical envelope modeling for laser wakefield acceleration, in *Proceedings of the 4th European Advanced Accelerator Concepts Workshop* (Elba Island, 2019).
[40] C. Benedetti et al., Presentation for the Advanced Accelerator Concepts Workshop 2016, https://indico.syntek.org/event/4/contributions/23/attachments/59/76/Benedetti_AAC2016.pdf (slide 12).

[41] M. Chen, E. Cormier-Michel, C. G. R. Geddes, D. L. Bruhwiler, L. L. Yu, E. Esarey, C. B. Schroeder, and W. P. Leemans, Numerical modeling of laser tunneling ionization in explicit particle-in-cell codes, J. Comput. Phys. 236, 220 (2013).

[42] C. B. Schroeder, J.-L. Vay, E. Esarey, S. S. Bulanov, C. Benedetti, L.-L. Yu, M. Chen, C. G. R. Geddes, and W. P. Leemans, Thermal emittance from ionization-induced trapping in plasma accelerators, Phys. Rev. Accel. Beams 17, 101301 (2014).

[43] J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai et al., Smilei: A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation, Comput. Phys. Commun. 222, 351 (2018).

[44] A. Beck, J. Derouillat, M. Lobet, A. Farjallah, F. Massimo, I. Zemzemi, F. Perez, T. Vinci, and M. Grech, Adaptive SIMD optimizations in particle-in-cell codes with fine-grain particle sorting, Comput. Phys. Commun. 244, 246 (2019).

[45] A. F. Lifschitz and V. Malka, Optical phase effects in electron wakefield acceleration using few-cycle laser pulses, New J. Phys. 14, 053045 (2012).

[46] J. Faure, D. Gustas, D. Guénot, A. Vernier, F. Böhle, M. Ouillé, S. Haessler, R. Lopez-Martens, and A. Lifschitz, A review of recent progress on laser-plasma acceleration at kHz repetition rate, Plasma Phys. Controlled Fusion 61, 014012 (2018).

[47] R. Nuter, L. Gremillet, E. Lefebvre, A. Lévy, T. Ceccotti, and P. Martin, Field ionization model implemented in particle in cell code and applied to laser-accelerated carbon ions, Phys. Plasmas 18, 033107 (2011).

[48] N. B. Delone and V. P. Krainov, Tunneling and barrier-suppression ionization of atoms and ions in a laser radiation field, Phys. Usp. 41, 469 (1998).

[49] C. Benedetti, A. Sgattoni, G. Turchetti, and P. Londrillo, Aladyn: A high-accuracy PIC code for the Maxwell-Vlasov equations, IEEE Trans. Plasma Sci. 36, 1790 (2008).

[50] G. Z. Sun, E. Ott, Y. C. Lee, and P. Guzdar, Self focusing of short intense pulses in plasmas, Phys. Fluids 30, 526 (1987).

[51] S. He and V. H. Weston, Wave-splitting and absorbing boundary condition for Maxwell's equations on a curved surface, Math. Comput. Sim. 50, 435 (1999).

[52] H. Barucq and B. Hanouzet, Asymptotic behavior of solutions to Maxwell's system in bounded domains with absorbing silver-Mueller's condition on the exterior boundary, Asympt. Anal. 15, 25 (1997).

[53] E. Esarey and M. Pilloff, Trapping and acceleration in nonlinear plasma waves, Phys. Plasmas 2, 1432 (1995).

[54] J. Faure, Plasma injection schemes for laser-plasma accelerators, CERN Yellow Rep. 1, 143 (2016).

[55] V. I. Berezhiani and I. G. Murusidze, Interaction of highly relativistic short laser pulses with plasmas and nonlinear wake-field generation, Phys. Scr. 45, 87 (1992).

[56] E. Esarey, P. Sprangle, J. Krall, A. Ting, and G. Joyce, Optically guided laser wake-field acceleration, Phys. Fluids B: Plasma Phys. 5, 2690 (1993).

[57] D. Teychenné, G. Bonnaud, and J.-L. Bobin, Wave-breaking limit to the wake-field effect in an underdense plasma, Phys. Rev. E 48, R3248 (1993).

[58] C. Rechatin, X. Davoine, A. Lifschitz, A. Ben Ismail, J. Lim, E. Lefebvre, J. Faure, and V. Malka, Observation of Beam Loading in a Laser-Plasma Accelerator, Phys. Rev. Lett. 103, 194804 (2009).

[59] X. L. Xu, C.-H. Pai, C. J. Zhang, F. Li, Y. Wan, Y. P. Wu, J. F. Hua, W. Lu, W. An, P. Yu et al., Nanoscale Electron Bunching in Laser-Triggered Ionization Injection in Plasma Accelerators, Phys. Rev. Lett. 117, 034801 (2016).

[60] F. Massimo, S. Atzeni, and A. Marocchino, Comparisons of time explicit hybrid kinetic-fluid code architect for plasma wakefield acceleration with a full PIC code, J. Comput. Phys. 327, 841 (2016).

[61] B. Quesnel and P. Mora, Theory and simulation of the interaction of ultraintense laser pulses with electrons in vacuum, Phys. Rev. E 58, 3719 (1998).