2.1 Simulation set-up

In our 2-D investigation of the autoresonant PBWA, the two co-propagating laser beams are chosen to have parallel linear polarization and identical intensities. They have the shape of a sixth-order super-Gaussian temporal profile  $\propto \exp [-(t/T_{\rm pulse})^6]$, and $\propto \exp [-(y/w_0)^6]$ in the transverse ($y$) direction, where the width of the laser beam $w_0$ is set to be $2.4{\rm \pi} k_{\rm pe}^{-1}$. Here, $k_{\rm pe}=\omega _{\rm pe}/c$ is the plasma wavenumber, $\omega _{\rm pe}=\sqrt {n_{e} e^2/\epsilon _0 m_e}$ is the electron plasma frequency, $n_{e}$ is the background electron density, $-e$ and $m_e$ are the electron charge and mass, and $c$ is the speed of light in vacuum.

The laser and plasma parameters are chosen based on our previous 1-D study to reach the wave-breaking limit. Namely, the amplitudes of the two lasers (indicated by subscripts $1$ and $2$) are $a_1=a_2= 0.2$ in terms of the normalized vector potential $a = e A/m_ec$. The ratio of homogeneous plasma density $n_{e}$ and the critical density $n_{cr}$ is $0.0004$, $n_{cr}=\omega _{1}^2m_e\epsilon _0/e^2$ is the critical density corresponding to a laser frequency $\omega _1$, with $\omega _1$ the central frequency of the beam one (1). Before entering the homogeneous plasma, a linear density ramp with the gradient length $L_{\rm ramp}k_{\rm pe}=40{\rm \pi}$ is applied. A chirp rate $\alpha = -0.0014$ is applied to the first laser beam, giving the frequency difference of the two laser beams $\Delta \omega =\omega _{\rm pe}[1+\alpha (t-t_0)\omega _{\rm pe}]$. At time $t_0=22.5{\rm \pi} /\omega _{\rm pe}$, the frequency difference between the two lasers matches the nominal plasma frequency, $\Delta \omega =\omega _1-\omega _2=\omega _{\rm pe}$. Here, the corresponding laser frequencies are $\omega _1/\omega _{\rm pe}=50$ and $\omega _2/\omega _{\rm pe}=49$. The laser duration $T_{\rm pulse}\omega _{\rm pe}=64{\rm \pi}$ is chosen to promote the autoresonant growth of the plasma wave until the wave-breaking limit. Further details on the choice of these parameters are given in Appendix A. The ions are set to be immobile, as they are not essential for the dynamics on the time scale of the laser pulse passing by (Mora et al. Reference Mora, Pesme, Héron, Laval and Silvestre1988). The longitudinal and transverse spatial resolutions are ${{\rm d} \kern0.7pt x}=0.008k_{\rm pe}^{-1}$ and ${{\rm d} y}=2\,{\rm d}\kern0.7pt x$, respectively. Table 1 shows the experimentally relevant parameters for a standard chirped pulse amplified $800\,{\rm nm}$ Ti:sapphire laser.

Table 1. Simulation parameters.

2.2 Autoresonant plasma wave excitation in two dimensions

First, we consider the growth phase of the autoresonant beat-wave excitation, comparing 1-D and 2-D results. Figure 1 shows the longitudinal electrostatic field $E_L$ in the 2-D simulation at time $t \omega _{\rm pe}\approx 250$. The electric field is normalized to the Rosenbluth–Liu limit $E_{\rm RL}=(16a_1a_2/3)^{1/3}E_0$, where $E_0=m_ec\omega _{pe}/e$ is the cold, non-relativistic wave-breaking field. The longitudinal coordinate is the co-moving position $\xi =\omega _{\rm pe}(x/c-t)$. Figures 1(a) and 1(b) correspond to simulations without chirp and with a chirp rate $\alpha =-0.0014$, respectively, that was deemed optimal for field enhancement with these parameters. The field enhancement with respect to the zero chirp case is evident in figure 1(a,b). The maximum value of the electric field in the chirp case is close to the cold, non-relativistic wave-breaking threshold $E_0$ that takes the value $E_0/E_{\rm RL}\approx 1.67$.

Figure 1. Longitudinal electric field of the plasma wave normalized to the RL field, $E_L / E_{\rm RL}$, in the $2 \mathrm {D}$ simulations, evaluated at $t \omega _{\mathrm {pe}} \approx 250$. The duration of the laser pulses is $T_{\text {pulse}} \omega _{\text {pe}} \approx 64 {\rm \pi}$. (a) No chirp. (b) With chirp rate $\alpha =-0.0014$. (c) Electric field at the symmetry axis ($y=0$). The blue (red) line corresponds to simulations without (with) chirp in two dimensions. Corresponding 1-D simulation data are shown by the light-blue (grey) lines without (with) chirp. The black dashed line indicates the wave-breaking field.

Notice that, as clearly visible in figure 1(a), for the case $\alpha =0$, the plasma wave amplitude is larger off-axis than closer to the symmetry axis, in the region $\xi \approx 100$–$150$. This is caused by the finite size of the laser field in the transverse direction and the related intensity variation. As a consequence, the growth and saturation time of the plasma wave excitation varies transversally. Specifically, the saturation time in the standard PBWA scheme is approximately $t_{\rm sa}\omega _{\rm pe}=3.68(\omega _{\rm pe}\tau _{r})^{4/7}(a_1a_2)^{-2/7}$ (Forslund et al. Reference Forslund, Kindel, Mori, Joshi and Dawson1985), where $\tau _{r}$ is the rise-time of the laser beam's temporal profile, introduced in § 2.1, leading to slower saturation of the plasma wave off-axis. Eventually, the RL limit is also reached at the off-axis regions, as visible in the plot.

Figure 1(c) shows the electric field variation along the symmetry axis, i.e. $E_L(\xi, y=0)$. The blue (red) line corresponds to the case without (with) chirp in two dimensions. For reference, we also present the 1-D simulation results, shown by the light-blue (grey) line without (with) chirp. The no-chirp results in one dimension and two dimensions overlap, and show that the RL value is indeed the limit of the achieved electric field; the maximum of the blue curves is at $E_L/E_{RL}=1$. When the optimal chirp value of $\alpha =-0.0014$ is employed, the autoresonant growth of the plasma wave is achieved and $E_L/E_{RL}=1$ is exceeded, but the overlap of the 2-D and 1-D results (red and grey, respectively) extends only for $\xi \gtrsim 125$. Below this value, the 2-D results start to deviate from the 1-D ones. The reasons for this are discussed in detail in § 3.

2.3 Energetic electron dynamics

After concluding that autoresonant plasma wave excitation is possible in two dimensions, we will now turn to the electron acceleration aspect of the process. Figure 2(a,b) show the electron phase space in terms of the co-moving coordinate $\xi$ and the kinetic energy reached by the electrons $U_{{\rm kin}}$ (the data are integrated over the transverse coordinate). Self-injected electrons can be effectively accelerated by the electric field structure of the plasma wave and, as expected, electrons can reach higher energies – in this case, $200\,{\rm eV}$ at $t\omega _{\rm pe}\approx 750$ – when a finite chirp of $\alpha =-0.0014$ is employed, allowing for autoresonant field growth; compare figure 2(b) with figure 2(a), where $\alpha =0$. Notice that in the chirped case, electrons are mainly accelerated in the region of $100< \xi < 150$, while in the region of $50 < \xi < 100$, we observe energetic particles but at lower energies. This was also observed in the analogous 1-D kinetic simulation (Luo et al. Reference Luo, Riconda, Pusztai, Grassi, Wurtele and Fülöp2024), and can be understood in the following way. Because of autoresonance, the plasma wake field grows with decreasing $\xi$. In the fluid model for the parameters chosen here (see Appendix A), it can become larger than $E_0$, but in kinetic simulation, once it reaches $E_0$, it starts to trap very efficiently particles that in turn stop the autoresonance. As a result, the field is highest in the range of $100< \xi < 150$ with value slightly below $E_0$. This is clearly visible in figure 1(c) for both the 2-D and 1-D chirped simulations (grey and red curves), the decrease of the field for $\xi < 100$ being more pronounced in the 2-D case.

Figure 2. Kinetic energy spectrum of electrons as a function of the co-moving coordinate $\xi$, integrated over $y$, at $t \omega _{\mathrm {pe}} \approx 750$ are plotted (a) without chirp and (b) with chirp rate $\alpha =-0.0014$. (c) Time traces of the highest electron kinetic energy for a variety of cases. The red ($\alpha =-0.0014$) line and blue ($\alpha =0$) line are in two dimensions, the grey ($\alpha =-0.0014$) line and light blue ($\alpha =0$) line are in one dimension. The black dashed and dotted lines indicate the kinetic energies of an electron that was accelerated by a constant electric fields at the wave-breaking ($E_0$) and the RL limit ($E_{\rm RL}$). (d) Temporal evolution of the maximum electric field $E_{L, \max }$ is shown in magenta (right axis). If an electron was accelerated with this time-varying field, it would reach energies as shown by the green curve (left axis). For reference, the dashed and the red lines of panel (c) are repeated.

We now focus on the region $\xi > 100$. To assess the efficiency of the acceleration process, we track the kinetic energy of electrons reaching the highest energies (note that the acceleration starts at slightly different times in the various simulations). In figure 2(c), we compare the temporal evolution to estimate the acceleration efficiency corresponding to the acceleration of relativistic electrons by a constant electric field at the RL limit, $E_{\rm RL}$ (dotted line), and at the wave breaking value, $E_0$ (dashed line). As expected, we find that in the no-chirp cases, the kinetic energy mostly increases with a steepness dictated by the RL field limit. The 1-D simulation (light blue) stably accelerates at this rate up to the end of the considered time range, while we observe a reducing slope in two dimensions (dark blue). In the simulations with chirp, the electron energy initially increases closely following the wave-breaking field estimate. At later times, the slope changes slightly in one dimension and more significantly in the 2-D simulation. The 1-D and 2-D autoresonant results (grey and red line, respectively) quite closely follow each other up to $t\omega _{\rm pe}\approx 450$, indicating a similar performance in the early stage of the acceleration.

Two reasons can be responsible for the difference between the autoresonant simulations and the estimate indicated by the dashed line. First, the maximum electric field might not be sustained at the constant level of $E_0$ and, second, the particles are not necessarily positioned to be optimally accelerated by the maximum electric field. We find that the first effect is the one mainly at play. Indeed, the kinetic energy of the most energetic particle from the simulations $U_{{\rm kin}, {\rm sim}}$ as function of time should follow $U_{{\rm kin}, {\rm max}}$ defined by

(2.1)\begin{equation} U_{{\rm kin}, {\rm max}}(t)=\int_{t_0}^t eE_{L,{\rm max}}(t')c\,{\rm d}t', \end{equation}

where $E_{L, {\rm max}}$ is the maximum instantaneous electric field. In figure 2(d), the time evolution of $E_{L, {\rm max}}$ is shown (magenta, right axis) along with the corresponding $U_{{\rm kin}, {\rm max}}$ (green) and $U_{{\rm kin}, {\rm sim}}$ (red). While $E_{L, {\rm max}}$ remains in the vicinity of $E_0$, it is somewhat lower most of the time. Thus, the expected particle energy (green curve) departs downward from the theoretical upper bound (dashed). The actual electron energy (red curve), closely follows $U_{{\rm kin}, {\rm max}}$ until $t\omega _{\rm pe}\approx 400$, and then starts to deviate. This marks the beginning of a 2-D effect of the system that is strong enough to degrade the efficiency of acceleration. The responsible physical processes are to be detailed in § 3. We note that while the simulations presented here do not include a plasma–vacuum interface, additional simulations confirmed that, while the lower energy part of the electron distribution is pulled back by space charge forces at such a boundary, it has a negligible effect on the most energetic electron population.

Next, we consider the transverse distribution of high-energy particles. Figure 3 shows the transverse distribution of the most energetic electrons as a function of time. At each time, electrons in the uppermost $10\,\%$ of the currently spanned energy range are identified and their transverse distribution is colour plotted. The transverse integral of this colour plot is shown as a red curve (right axis). More specifically, defining the ‘high-energy’ electron density by $n_{\rm HE}(\xi,y,t)=\int _{0.9 U_{{\rm kin},{\rm max}}(t)}^{U_{{\rm kin},{\rm max}}(t)}f_e(U_{{\rm kin}},\xi,y,t)\,{\rm d}U_{{\rm kin}}$, with the electron distribution function $f_e$, the colour plot shows $N_{\rm HE,y}(y,t)=\int \,{\rm d}\xi \, n_{\rm HE}(\xi,y,t)/n_e$, with the nominal electron density $n_e$, while the red line is $N_{\rm HE}(t)=\int \,{{\rm d} y} \,N_{\rm HE,y}(y,t)$; note that here, the $y$ coordinate is normalized to $k_{\rm pe}^{-1}$ (similarly to $\xi$). When interpreting this image, we must keep in mind that the considered energy range corresponds to continuously increasing energies, since $U_{{\rm kin},{\rm max}}=U_{{\rm kin},{\rm max}}(t)$. The highest number of energetic electrons is reached at $t\omega _{\rm pe}\simeq 400$, then it rapidly drops and stays around $20\,\%$. After this time, the electron energy keeps increasing but for a smaller population than at earlier times. Up to $t\omega _{\rm pe}=350$, the fast electrons are mostly localized at the centre, while we observe a rapid widening of the transverse distribution after that. Even later, at approximately $t\omega _{\rm pe}=550$, the most energetic electrons are found off-axis. Note that the highest energy electrons have been continuously accelerated, while those that fail to be accelerated to $90\,\%$ of the highest energy are not present in this visualization.

Figure 3. Colour plot: the high-energy electron distribution ($N_{{\rm HE},y}(y,t)$, integrating in $\xi$ over the box length) is plotted. Red curve (right axis) shows $N_{{\rm HE}}(t)$ (integrating over the $y$-axis). High-energy electrons are defined, for the plot, as those in the uppermost $10\,\%$ of the instantaneous energy range are considered. The two vertical dashed lines mark the times considered in figure 4, when the energetic electrons are localized on-axis and off-axis, respectively.

We now consider two snapshots, at the times indicated by the white dashed vertical lines in figure 3, and we show in figure 4 the logarithm of the relative number of energetic electrons as function of $\xi$ and $y$. At $t\omega _{\rm pe}=350$, shown in figure 4(a), the maximum kinetic energy is $U_{{\rm kin},{\rm max}}=86\,{\rm MeV}$ and, as expected, the energetic electron bunches are close to the symmetry axis. At the later time, $t\omega _{\rm pe}=550$, shown in figure 4(b), the distribution peaks off-axis and the instantaneous maximum energy is $U_{{\rm kin},{\rm max}}=160\,{\rm MeV}$. The bunches are localized at the minimum of the electric field where the acceleration is the strongest, indicated by the red curves in the same figures. At the end of the simulation, the cumulative charge can be estimated to be $15\,{\rm pC}$, assuming that the electron bunch is axisymmetric within a region of length $15\,\mathrm {\mu }{\rm m}$ and a transverse extent $24\,\mathrm {\mu }{\rm m}$, given a laser wavelength of $800 \,{\rm nm}$.

Figure 4. Spatial distribution of $n_{{\rm HE}}(\xi,y)$ at (a) $t\omega _{\rm pe}\approx 350$ and (b) $t\omega _{\rm pe}\approx 550$. These times are marked with dashed lines in figure 3. The longitudinal electric field at the axis, ${E_L(\xi, y=0)/E_0}$, is indicated by the red curves (right axes).

We may quantify the angular spread of the electron beam by $\bar {\theta } =[\int \,{\rm d}^3\boldsymbol {p}(p_y/p_x)^2 f_e(\boldsymbol {p})$ / $\int \,{\rm d}^3\boldsymbol {p}f_e(\boldsymbol {p})]^{1/2}$, with the integrals covering only the uppermost $10\,\%$ of the electron energies. We find that, in spite of the rich transverse dynamics, $\bar {\theta }$ only reaches values of ${\sim }0.02$ throughout the simulation. In the following section, we will delve into the reasons for the transverse dynamics, which shall also shed light on why the acceleration efficiency is reduced when moving from one to two dimensions.

3.1 The weakly nonlinear stage

In the early phase of the evolution of the plasma wave – apart from a downward shift of its wavenumber as the plasma wave amplitude increases – the wavenumber spectrum is affected by forward Raman scattering at a slight angle with respect to the laser propagation. This leads to the appearance of a double signal in the electromagnetic spectrum, mainly around $k_L/k_{\rm pe} = 50$ (Forslund et al. Reference Forslund, Kindel, Mori, Joshi and Dawson1985), and a broadening of both the electromagnetic and the electron plasma wave spectrum in the transverse direction. Moreover, Stokes and anti-Stokes scattering results in harmonics of the electromagnetic wave. In Appendix B.1, we plot the Fourier spectrum of the waves, and all these components are clearly visible. At this stage, the only new effect compared with the 1-D case, discussed by Luo et al. (Reference Luo, Riconda, Pusztai, Grassi, Wurtele and Fülöp2024), is the transverse broadening and the double signal around $k_s/k_{\rm pe} = 50$, while the Stokes and anti-Stokes components and the longitudinal widening were already present in the 1-D simulations. Overall, the autoresonant excitation and the particle acceleration appear robust with respect to the presence of the transverse Raman scattering.

Another 2-D effect that can be identified is due to the transverse variation of the laser intensity and the corresponding nonlinear correction to the wavelength. This leads to the formation of slightly curved plasma wave fronts, while the plasma wave structure remains coherent, as was also observed in some stimulated Raman scattering simulations (Rose & Yin Reference Rose and Yin2008; Yin et al. Reference Yin, Albright, Bowers, Daughton and Rose2008; Masson-Laborde et al. Reference Masson-Laborde, Rozmus, Peng, Pesme, Hüller, Casanova, Bychenkov, Chapman and Loiseau2010). More details on this wavefront bowing are given in Appendix B.2. Similarly to the mentioned parametric wave coupling processes, the wavefront bowing effect has no appreciable negative effect on the efficiency of the autoresonant PBWA.

3.2 Transverse magnetic field growth due to a Weibel-like instability

We now focus our attention on processes that become important around the time the acceleration efficiency starts to degrade compared with the 1-D simulations. The first of these effects is the growth of a Weibel-like instability (Fried Reference Fried1959; Weibel Reference Weibel1959), which provides a seed for subsequent plasma wave filamentation (Masson-Laborde et al. Reference Masson-Laborde, Rozmus, Peng, Pesme, Hüller, Casanova, Bychenkov, Chapman and Loiseau2010).

In figure 5, the longitudinally averaged transverse magnetic field $\bar {B}_z$ is shown as a function of the transverse coordinate $y$, normalized as $\bar {H}_z=e\bar {B}_z/m_e\omega _1$. The times shown illustrate the development of transverse magnetic field filaments (time increases from panel a to panel g). The average in the longitudinal direction is performed over the $\xi$-range indicated by the dashed rectangle in figure 6(a). In this region, the plasma wave has a significant amplitude, but wave breaking has not yet been reached, and only a few particles are trapped. Indeed, the $\xi - p_x$ phase space (not shown here) has a very regular structure with longitudinal velocity up to $\simeq c$. Stronger particle acceleration is instead observed in the range $100< \xi < 150$, leading to a more disordered magnetic field structure. In the early stages shown in figure 5(a,b), the magnetic field starts growing at the edges of the laser beams ($y\sim \pm 2.4{\rm \pi}$) because of the large-scale current generated by the density and velocity perturbation in the plasma wave (with a homogeneous component plus one varying in the longitudinal direction with wavenumber $\sim 2 k_{\rm pe}$), as discussed by Gorbunov, Mora & Antonsen (Reference Gorbunov, Mora and Antonsen1996, Reference Gorbunov, Mora and Antonsen1997) and Sheng, Meyer-ter Vehn & Pukhov (Reference Sheng, Meyer-ter Vehn and Pukhov1998). The plasma current comes from the collective electron motion induced by the longitudinal electric field of the plasma wave. Note that in the central region of the laser beams, no magnetic field is generated, as the displacement current compensates for the plasma current. Some filaments start to appear in the central region and become more pronounced by the time shown in figure 5(c). As expected for a Weibel-like instability, the filaments grow with a typical wavenumber of the order of $k_{\rm pe}$. Both wavelength and field amplitude increase with time (figure 5d–g) up to saturation, happening at $\omega _{\rm pe}t \simeq 430$.

Figure 5. Longitudinally averaged transverse profiles of the $z$-component of the magnetic field, taken at different times, normalized as $\bar {H}_z=e\bar {B}_z/m_e\omega _1$. The averaging is performed over the same $\xi$-range as covered by the dashed rectangle in figure 6(a).

Figure 6. (a) Spatial variation of the magnetic field $H_z=eB_z/(m_e\omega _1)$ at time $t \omega _{\rm pe}\approx 350$. The filamentation of the magnetic field is clearly visible inside the magenta dashed box, where the plasma wave is driven to significant amplitudes while wave breaking is not reached. (b) Zoom-in of the magnetic field, taken from the region highlighted in panel (a). (c) Temperature anisotropy $A$ in the same region as panel (b). Note that for the visibility, the range of $A$ is capped at $1$, while $A$ can reach values of $\approx 3$ but only in very restricted regions.

To understand the origin of these magnetic field filaments, figure 6(a) shows the magnetic field $H_z=eB_z/(m_e\omega _1)$ as a function of $\xi$ and $y$ at time $t\omega _{\rm pe}\approx 350$ (corresponding to the longitudinally averaged result of figure 5c). The transverse limits of the dashed magenta rectangle are located at $y=\pm {\rm \pi}$ ($y$ is normalized by the plasma wavenumber $k_{\rm pe}$) and are chosen such that the rectangle covers the flat-top region of the transverse profile of the laser beams. Inside this marked area, magnetic field filaments are clearly visible, as highlighted by figure 6(b) which shows a zoom-in on the same data.

To identify the drive of the instability and estimate its growth rate, we define the longitudinal and transverse thermal speeds $v_{tx}$ and $v_{ty}$, respectively, as follows: $v_{tj}^2(\xi,y)=n^{-1}\int (v_{j}-V_j)^2 f(\xi,y,\boldsymbol {v})\,{\rm d}^3v$, where $n=\int f(\xi,y,\boldsymbol {v})\,{\rm d}^3v$, $V_j=n^{-1}\int v_j f(\xi,y,\boldsymbol {v})\,{\rm d}^3v$ and $j=\{x,y\}$. These definitions allow us to introduce an effective temperature anisotropy parameter $A(\xi,y)=(v_{tx}^2/v_{ty}^2)-1$, adapting previous definitions (Ruyer et al. Reference Ruyer, Gremillet, Debayle and Bonnaud2015; Silva et al. Reference Silva, Schoeffler, Vieira, Hoshino, Fonseca and Silva2020; Silva, Afeyan & Silva Reference Silva, Afeyan and Silva2021) to our conditions. Its values are plotted in figure 6(c) at $t\omega _{\rm pe}\approx 350$. Since most of the acceleration happens longitudinally, the longitudinal spread of the distribution function and hence the thermal speed is higher than the transverse one, yielding values $A>0$. Significant values of $A$ up to $\approx 3$ are reached in patches, as shown in figure 6(c). The anisotropy in the distribution makes the plasma unstable to a Weibel-like instability, leading to the growth of the magnetic filaments.Footnote ¹

The time evolution of the maximum value of $|H_z|$ and the magnetic energy content inside the marked rectangular region of figure 6(a), $\mathcal {E}_{B_z}=\int \int {|eB_z/(m_e\omega _1)|^2}\,{\rm d}\xi \,{{\rm d} y}$, are shown in figure 7, with solid and dash-dotted lines, respectively. We choose to focus only on the magnetic field inside the marked region, where the Weibel-like instability is the dominant effect that generates the magnetic filaments. We hence neglect the field growing at the intensity gradients of the laser beams and the more complex structures in the rear region ($100<\xi < 140$), where efficient acceleration occurs near the axis, modifying the plasma wave structures and the phase space in a non-trivial way. This allows us to identify the typical growth of the field generated by the Weibel-like instability, as discussed below.

Figure 7. Maximum of the $z$-component of the magnetic field (solid black curve, $10 |H_z|_{\rm max}$) and energy contained in the $z$-component of the magnetic field (black dash-dotted, $\mathcal {E}_{B_z}$); both within the region indicated by the dashed magenta rectangle in figure 6(a). The solid magenta line indicates an exponential growth with the estimated growth rate of the Weibel-like instability. The red dashed line is the quasi-longitudinal electric field energy $\mathcal {E}_{\rm in}$ and the blue dotted line is the electric field energy scattered outside the quasi-longitudinal wavenumber region $\mathcal {E}_{\rm out}$.

The general trends of $\mathcal {E}_{B_z}$ and $|H_z|_{\rm max}$ are similar, and reflect the observations from figure 5. A first early rapid increase and saturation is observed at $t \omega _{\rm pe}\approx 190$. This initial phase of magnetic field generation is caused by the development of the self-generated quasi-static magnetic field as stated by Gorbunov et al. (Reference Gorbunov, Mora and Antonsen1996, Reference Gorbunov, Mora and Antonsen1997) and Sheng et al. (Reference Sheng, Meyer-ter Vehn and Pukhov1998). This is followed by another significant growth until $t\omega _{\rm pe}\approx 400$, after which the magnetic field saturates around a higher level, with some oscillations. This second growth phase corresponds to the Weibel-like instability. The maximum growth rate of this Weibel-like instability, $\gamma _{\rm wb}$, in terms of the anisotropy parameter above defined, is approximately given by (Okada, Yabe & Niu Reference Okada, Yabe and Niu1978; Satou & Okada Reference Satou and Okada1997; Sugie, Ogawa & Okada Reference Sugie, Ogawa and Okada2006; Okada & Ogawa Reference Okada and Ogawa2007; Zhou et al. Reference Zhou, Zheng, Liu and Cao2022)

(3.1)\begin{equation} \gamma_{\rm wb}/\omega_{\rm pe}=\sqrt{\frac{8}{27{\rm \pi}}}\frac{v_{ty}A^{3/2}}{A+1}. \end{equation}

Considering values extracted from the simulation for the transverse thermal speed $v_{ty}\approx 0.1c$ and temperature anisotropy $A \approx 2.5$ – where the latter is close to the upper end of the observed range of $A$ values – yields a growth rate $\gamma _{\rm wb}/\omega _{\rm pe}\approx 0.035$. Such an exponential growth, $\propto \exp (\gamma _{\rm wb}t)$, is indicated by the magenta line in figure 7, showing a satisfactory agreement with the growth of $|H_z|_{\rm max}$ (solid black curve). We note that a similar comparison between the magnetic field growth and the predicted temperature anisotropy-driven Weibel growth rate at a lower laser intensity, $a=0.12$, showed similarly good agreement.

In the weakly nonlinear stage, discussed in detail in Appendix B.1, we already observed that forward/side Raman scattering and wave bending will tend to deform the plasma wave wavefront and introduce a transverse component to its wavenumber. The Weibel-like instability also affects the transverse wavenumber of the plasma wave, $k_{Ly}$. To quantify the effect of the Weibel-like instability on the transverse wavenumber of the plasma wave $k_{Ly}$, we can take advantage of the fact that it will induce shorter transverse lengths (larger $k_{Ly}$) than that already identified in the weakly nonlinear stage. The longitudinal wavenumber $k_{Lx}$ instead is given by the matching condition, roughly $k_{Lx}\sim k_{\rm pe}$, so we may consider a $k_{Lx}$ range that provides a comfortable margin around this,Footnote ² such as $k_{Lx}\in [0.5,1.5]k_{\rm pe}$.

The characteristic transverse broadening of the wavenumber spectrum due to the Raman process can be estimated as $\delta k_{Ly}=\sqrt {k_{\rm pe}^2-k_{\rm np}^2}$, where $k_{\rm np}$ is the nonlinear plasma wavenumber introduced in Appendix B.1 and discussed in figure 12 of Luo et al. (Reference Luo, Riconda, Pusztai, Grassi, Wurtele and Fülöp2024): $k_{\rm np}\approx k_{\rm pe} -3 k_{\rm pe} (E_L/E_0)^2/16$. If $|k_{Ly}|<\delta k_{Ly}$, we expect Raman scattering and wavenumber nonlinearity to dominate, with a corresponding ‘quasi-longitudinal’ electric energy content $\mathcal {E}_{\rm in}$. If $|k_{Ly}|>\delta k_{Ly}$, we expect additional processes, such as the Weibel-like instability to dominate, with electric energy content $\mathcal {E}_{\rm out}$. More specifically,

(3.2)\begin{equation} \left. \begin{aligned} \mathcal{E}_{\rm in} & ={\int_{0.5k_{\rm pe}}^{1.5k_{\rm pe}}\,{\rm d}k_{Lx}\int_{|k_{Ly}|<\delta k_{Ly}}}{{\rm d}k_{Ly} |E_L(k_{Lx},k_{Ly})|^2}, \\ \mathcal{E}_{\rm out} & ={\int_{0.5k_{\rm pe}}^{1.5k_{\rm pe}}\,{\rm d}k_{Lx}\int_{|k_{Ly}|>\delta k_{Ly}}}{{\rm d}k_{Ly} |E_L(k_{Lx},k_{Ly})|^2}. \end{aligned} \right\} \end{equation}

Optimally, the energy is mostly contained in the quasi-longitudinal field, while non-ideal 2-D effects may scatter energy outside this region, and the ratio $\mathcal {E}_\textrm {out}/\mathcal {E}_\textrm {in}$ thus quantifies the importance of these 2-D effects that reduce the transverse coherence of the plasma wave. In figure 7, the electric energies $\mathcal {E}_{\textrm {in}}$ and $\mathcal {E}_{\textrm {out}}$ are also shown by red dashed and blue dotted lines, respectively. Most of the energy is contained in the quasi-longitudinal field that saturates around $t \omega _\textrm {pe}\approx 190$. The energy scattered outside this wavenumber region is approximately following the same trends as the $\mathcal {E}_{B_z}$ and $|H_z|_\textrm {max}$ curves, with a first saturation around $t \omega _\textrm {pe}\approx 190$. A more pronounced growth phase occurs during $t\omega _\textrm {pe}\approx 300$–$400$, when the magnetic fields grow due to the Weibel-like instability, followed by a second saturation. It is noteworthy that around the time of the onset of the Weibel-like instability, $\mathcal {E}_{\textrm {in}}$ also starts to decrease, indicating a reduced coherence of the plasma wave, even if the ratio $\mathcal {E}_{\textrm {in}}/\mathcal {E}_{\textrm {out}}$ is always larger than one.

3.3 Strong nonlinearity and plasma wave filamentation

We have established the influence of Weibel-like instability on the transverse dynamics of the plasma waves: it gives rise to transverse magnetic perturbations and a scattering of electric field energy beyond that caused by Raman scattering, potentially leading to a reduction in the coherence of plasma waves. However, the direct impact of Weibel-like instability on the laser evolution in the plasma remains negligible as the instability quickly saturates. In this section, we explore how the later-stage evolution of the laser propagation in the strong-nonlinear phase, where density modulations and filamentation arise.

Figure 8 shows the wavenumber spectra of the electric field of the plasma wave $E_L$ (left column) and that of the electromagnetic wave $E_s$ (right column) at different time instances, with the corresponding wavenumbers denoted by $k_L$ and $k_s$, respectively. Here, we focus on times $t \omega _\textrm {pe} \ge 350$. At early times ($t \omega _\textrm {pe}=350$), the main mode spectrum of the plasma wave assumes an arc shape between the circle given by mode matching (dashed line) and that with the nonlinear shift accounted for (solid line). By the time $t \omega _\textrm {pe} = 630$, the arc-shaped spectrum is not clearly seen any longer and a stronger transverse broadening is observed. In addition, in the vicinity of $k_{s,x}\approx 50$ of the electromagnetic wave, the slightly downshifted and upshifted components start to evolve differently, the downshifted one spreading more in the transverse direction. The Stokes/anti-Stokes scattering is barely affected.

Figure 8. (a,b) Wavenumber spectra of the longitudinal electric field of the (a i,b i) plasma wave $E_L$ and the (a ii,b ii) electromagnetic wave $E_s$ at different times. The green curves correspond to the solutions of (a i,b i) (B1) and (a ii,b ii) (B2). In panel (a i,b i), the solid (dashed) curves include (do not include) corrections for the nonlinear wavenumber shift. In panel (a ii,b ii), different green curves represent different scattering order. (c) Temporal evolution of the purely longitudinal ($k_{s,y}=0$) wavenumber component of the main modes of the two laser beams. Vertical red lines indicate the maximum bandwidth of the first laser beam $\Delta \omega /\omega _\textrm {pe}\approx 0.28$.

In figure 8(c), the temporal evolution of the wavenumber spectrum $k_{s,x}$ of the two dominant laser beams is shown for the purely longitudinal mode ($k_{s,y}=0$). The two vertical red lines indicate the maximum bandwidth of the first, chirped, laser beam – located at $k_{s,x}/k_\textrm {pe}=50$ – that can be estimated as $\Delta \omega /\omega _\textrm {pe}\approx |\alpha |T_\textrm {pulse}=0.28$. Without the interaction of the lasers and the autoresonant creation of the plasma wave, the wavenumber spectrum of this chirped pulse would be nearly constant within this region. Instead, one observes that the signal in the vicinity of $k_{s,x}\approx 50$ becomes very weak, leaving two bands on the sides, partly related to laser depletion due to the plasma wave excitation, also present in one dimension and observed in our previous simulations, and partly related to scattering from $k_{s,y}=0$ (plotted) to $k_{s,y}\ne 0$, as evident in figure 8(b ii).

With the onset of autoresonance, considerable pump depletion, on account of the plasma wave excitation, is expected to take place. In terms of frequency, the width of the depleted gap can be estimated by $\Delta \omega _\textrm {eff}/\omega _\textrm {pe}\approx 3[1-(E_\textrm {RL}/E_0)^2]/16\approx 0.1$. This is approximately a third of the overall bandwidth of the chirped beam, consistent with the observation. This also allows us to estimate the length over which autoresonance occurs, $L_\textrm {AR}\approx \Delta \omega _\textrm {eff}/\omega _\textrm {pe}/|\alpha |=75$. This distance $L_\textrm {AR}$ is indeed representative of the spatial extent over which effective autoresonance takes place; see the $\xi \approx 125$–$200$ range of figure 1(c).

Figure 9(a–c) captures snapshots of the plasma wave and the electromagnetic field, respectively, at time $t \omega _\textrm {pe}\approx 550$, the time moment between those in figures 8(a) and 8(b). In figure 9(a), the electric field carried by the plasma wave $E_L$ over the range of $\xi \in [100,250]$ is shown, where the transverse modulation is observed, particularly occurring over the range of $\xi \in [150,200]$, as exhibited in figure 9(b). The transverse regions where the beating patterns of the laser are weak, such as the gap at $|y|\approx 4$ in figure 9(c),Footnote ³ are correlated with a reduced plasma wave amplitude; see the reduced amplitude transverse bands in figure 9(b).

Figure 9. At $t\omega _\textrm {pe}\approx 550$, (a) the electric field of the plasma wave $E_L$, over the range $\xi \in [100,250]$, (b) $E_L$ in the range $\xi \in [150,200]$ and (c) the electric field of the electromagnetic wave $\tilde {a}$, normalized by the initial total laser amplitude $a_1+a_2$ in the range $\xi \in [150,200]$.

At this stage of the evolution, we conjecture that the appearance of transverse structure of the plasma wave and of the laser pulses will self-sustain and enhance the transverse modulation. To justify this, we performed an almost identical simulation, adding a co-propagating probe laser with a linear polarization orthogonal to that of the dominant laser pulses, having only $5\,\%$ of their amplitudes; we denote the normalized electric field of this probe beam by $\tilde {a}_{p}$ and its initial value by $a_{p}$ with $a_{p}=5\,\%a_{1,2}$.

The density perturbation of the plasma wave ($\delta n_e+n_e)/n_e$ is shown in figure 10(a), along with the electric field of the probe laser in figure 10(b), at the same time instance as in figure 9. To emphasize the impact of the density perturbation on the laser field, in figure 10(b), we only plot the fields of the probe laser whenever $|\tilde {a}_{p}|/a_{p}\geq 1$ (these appear as blue and red spikes). The density perturbation modifies the index of refraction $\eta _r$, thereby affecting the transverse propagation of the electromagnetic wave. We may approximate the index of refraction by

(3.3)\begin{equation} \eta_r \approx 1-\frac{\omega_{\rm pe}^2}{2\omega_0^2}\frac{\delta n_e+n_e}{n_e\gamma}, \end{equation}

with the Lorentz factor $\gamma$ being comparable in magnitude to $\gamma _{\perp }=(1+a^2)^{1/2}$. To show the spatial variation of the refractive index, we draw the contours $\eta _r=0.9996$ (white) and $n_r=0.9998$ (black), calculated with the density variation shown in figure 10(a). As a result of the spatial variation of $\eta _r$, the laser beams undergo fragmentation, with the various segments exhibiting different characteristics. Specifically, the intensity of the beamlets follow closely the density contours. Both the laser beams and the plasma have developed transverse modulation – a feature not present at earlier stages – marking the transition into the strongly nonlinear stage of the process. The evolution of the plasma wave in this stage have now a non-negligible impact on particle acceleration. The resulting reduced acceleration efficiency is apparent in figure 2(c), and a non-trivial transverse re-arrangement of the acceleration process is reflected in figures 3 and 4(b).

Figure 10. A low power ($a_{p}=5\,\%a_{1,2}$) probe laser with orthogonal polarization to the lasers that create the plasma wave co-propagates together with the two dominant laser beams. At the time $t\omega _\textrm {pe}=550$ restricted to one plasma oscillation in $\xi$, (a) the plasma wave density profile and, in panel (b), blue (red) corresponds to large negative (positive) values of the probe laser field, and for the purpose of visibility, values of the probe field less in magnitude than the initial peak probe field are set to zero. The black and white contours correspond to refractive index values of $\eta _r=0.9998$ and $\eta _r=0.9996$. See text for details.