Similarity of magnetized plasma wake channels behind relativistic laser pulses with different wavelengths

Similarity of magnetized plasma wake channels behind
relativistic laser pulses with different wavelengths

A. Bierwage T.Zh. Esirkepov J.K. Koga A.S. Pirozhkov National Institutes for Quantum and Radiological Science and Technology, Rokkasho Fusion Institute, Rokkasho, Aomori 039-3212, Japan National Institutes for Quantum and Radiological Science and Technology, Naka Fusion Institute, Naka, Ibaraki 311-0193, Japan National Institutes for Quantum and Radiological Science and Technology, Kansai Photon Science Institute, Kizugawa, Kyoto 619-0215, Japan

Using particle-in-cell simulations of relativistic laser plasma wakes in the presence of an external magnetic field, we demonstrate that there exists a parameter window where the dynamics of the magnetized wake channel are largely independent of the laser wavelength . One condition for this manifestation of “limited similarity” is that the electron density is highly subcritical, so that the plasma does not affect the laser. The freedom to choose a convenient laser wavelength can be useful in experiments and simulations. In simulations, an up-scaled wavelength (and, thus, a coarser mesh and larger time steps) reduces the computational effort, while limited similarity ensures that the overall structure and evolutionary phases of the wake channel are preserved. In our demonstrative example, we begin with a terrawattpicosecond pulse from a laser with , whose field reaches a relativistic amplitude at the center of a sub-millimeter-sized focal spot. The laser is shot into a sparse deuterium gas () in the presence of a tesla-scale magnetic field. Limited similarity is demonstrated in 2D for and is expected to extend to shorter wavelengths. Assuming that this limited similarity also holds in 3D, increasing the wavelength to enables us to simulate the after-glow dynamics of the wake channel all the way into the nanosecond regime.

Similarity scaling, Relativistic laser wakefield, Magnetized plasma, PIC simulation, Code benchmark
journal: arXiv

1 Introduction

Laser-induced wake waves Akhiezer56 () have been studied extensively during the last decades, in particular in the context of particle acceleration Tajima79 () and for the generation of compact sources of energetic photons ELIwhitebook2011 (); Corde13 (); Pirozhkov17b (). Many interesting and useful applications have emerged, and there still exists much potential for expansion BulanovSV16 ().

Beginning with the present paper, we investigate the dynamics of laser-induced wakes in a parameter regime that, to our knowledge, has not been examined before: we consider a laser pulse with relativistic intensity Mourou06 () and study the formation and long-term evolution of the wake channel that this laser produces in a magnetically confined plasma as used in current mainstream thermonuclear fusion reactor designs IterPlan2018 (); Tobita17 (). Here, we begin with relativistic laser wakes in the cold plasma limit, which is valid up to the nanosecond scale, laying the basis for future work that will include thermal motion.

Present-day laser systems that are able to accelerate electrons to relativistic speeds operate with principal wavelengths in the near infrared: Danson15 (); Polyanskiy11 (). Since the electron densities in magnetically confined fusion (MCF) plasmas produced in tokamak and stellarator devices rarely exceed Greenwald02 (), they are effectively transparent to this kind of electromagnetic radiation, so that the laser pulse passes nearly unchanged through the plasma, while plowing away all electrons in its path. The accelerated electrons will be deflected by the force, where is the electron’s charge, its velocity vector, and is the externally imposed magnetic field. The magnetic field in present-day MCF devices is typically a few tesla strong, which can be regarded as an intermediate strength in the sense that ; i.e., the Larmor frequency


is comparable to the Langmuir frequency


where is the electron’s rest mass, the Lorentz factor, and the vacuum permittivity. (We employ SI units in this work.)

The effect of an externally imposed magnetic field on laser plasmas has received increasing attention during the last decade Hosokai06 (); Hur08 (); Hosokai10 (); Bourdier12 (); Schmit12 (); BulanovSV13 (); Ghotra15 (); Wilson17 (), but these works focused primarily on magnetization effects in the vicinity of the laser pulse. Here, our motivation is also to learn how the magnetized plasma channel evolves long after the laser pulse has gone.

Even without an external magnetic field, long-lived plasma channels can be routinely seen in experiments where intense laser pulses are shot through subcritical gas. Simulations have revealed the formation of a quasi-static azimuthal magnetic field and magnetic vortices that are maintained by electron currents flowing along the channel axis BulanovSV96 (); Kuznetsov01 (); BulanovSV05 (); Nakamura08 (); Nakamura10 ().

In order to study the after-glow dynamics numerically in our parameter regime of interest, the simulations need to cover a relatively long time interval (nanoseconds), include the ion response, and cover a relatively large spatial volume (centimeters). The simulation box should encompass the entire plasma channel around the focal point, with sufficiently large buffer zones in each direction to minimize boundary effects. At the same time, a high spatio-temporal resolution (sub-micron, sub-femtosecond) is needed in order to capture the evolution of the laser’s electromagnetic field and the motion of electrons accelerated to relativistic speeds. This constitutes a multi-scale problem and poses many challenges, including unphysical boundary effects, numerical heating, and a large computational cost.

The need for large computational resources is connected with the necessity to perform these simulations in three dimensions (3D), so as to capture the intrinsically 3D motion of charged particles in an external magnetic field . Although 2D simulations can be used to identify interesting parameter regimes and trends, the reduced degrees of freedom and enforced symmetries imply that great care is needed when trying to extract physical insight from 2D results. The ability to perform 3D simulations is important for identifying verifiable aspects of the 2D results and obtaining accurate qualitative and quantitative information that can be validated against experiments.

If the dynamical structures of interest are larger than the laser wavelength , the computational cost (memory, time) is determined by the value of , which determines the spatial resolution, the time step and (via its focal length) the required box size. It can be shown (A) that even a single long-time 3D simulation (not to speak of parameter scans and convergence tests) would be prohibitively expensive if it was to be performed with realistic laser wavelengths in the range . Under some conditions, however, this problem can be circumvented by using an artificially increased wavelength.

In this paper, we identify a parameter regime, where such an artificial scaling of the laser wavelength can be justified on the basis of so-called limited similarity or limited scaling Block67 (); Esirkepov12 (). In Section 2 we describe the simulation setup, where a relativistically strong laser pulse is shot into a sparse deuterium gas with highly subcritical electron density and a moderate external magnetic field. In Section 3, the properties of the magnetized plasma wake channel produced by a relativistic laser pulse with wavelength are illustrated in 2D. In Section 4, 2D simulations are used to demonstrate the limited similarity of the magnetized wake channel for laser wavelengths in the range . Using a scaled wavelength , we are able to perform long-time 3D simulations, and first results are reported in Section 5. A summary and conclusions are given in Section 6.

2 Simulation scenario, model and codes

2.1 Reference case

As a starting point and reference case, we consider a short picosecond-scale pulse from a terrawatt-scale gas laser with wavelength . The laser pulse travels through a quasi-neutral deuterium plasma with number density in the presence of a magnetic field, as is typical for present-day tokamaks.

We follow the dynamics of the laser-induced wake channel for up to 1 nanosecond. This means that the plasma can be safely assumed to be ideal (collisionless). Thermal motion will be ignored for simplicity, although it should be noted that, under typical tokamak conditions with thermal energies around , electrons travel about and deuterons about per nanosecond. The effect of thermal motion is briefly discussed in Section 6 and will be examined in a separate paper, which will be dedicated to the physics of the laser-induced wake channels. The present paper lays the conceptual and computational foundations.

A right-handed Cartesian coordinate system in the laboratory frame of reference is used and the laser travels in the positive direction. The simulation domain is


where is the trailing edge of the simulation window along the laser path and the time is denoted by . All simulations begin with . In simulations with a moving window, advances with the speed of light after a suitable delay.

The initial electron density is uniform, except for a narrow boundary layer, where it is ramped up or down. In mathematical form, it may be written as


with and


The thickness of the boundary layer where is ramped up or down linearly is . The thickness of the surrounding vacuum region is also , so that (same for ) in the transverse direction. Along the laser path we use and in the case of a moving window, or for a stationary window.

Figure 1: Spatial dependence of (a) the amplitude radius and (b) the normalized amplitude of a Gaussian pulse propagating through vacuum in the paraxial wave approximation. The red curves are for a laser wavelength of and blue for . In panel (b), solid lines are used for the 3D case and dashed lines for 2D.

2.2 Laser model

The Gaussian laser pulse is linearly polarized, with its electric and magnetic field vectors and oscillating along the and the axis, respectively. The laser wavelength, wavenumber and angular frequency are related as


We are interested in laser-plasma interactions near the threshold of the relativistic regime Mourou06 (), so we focus the laser such that its normalized amplitude


reaches at the center of the focal spot, which is located at .

Since the electron density in a tokamak is highly subcritical (transparent) for near-infrared light, the laser pulse propagates as if in vacuum, where its electromagnetic field can be well approximated by the solution of the paraxial wave equation Boyd (). The complex-valued amplitude of the electric field can then be written as


where () for 2D and for 3D. The auxiliary coordinate


has its origin at the focal point, where the laser field takes the form of a plane wave. The first line of Eq. (10) describes a harmonic wave. The second line is the Gaussian temporal envelope of the laser pulse with amplitude width parameter . The functions


describe the shape of the pulse envelope with curvature radius . The corresponding deformation of the optical wave front with respect to the - plane is captured by the complex exponents on the last two lines of Eq. (10), where


measures the phase shift of the wave front relative to the plane wave at . Also on the last two lines, there are factors describing the transverse profile of the laser spot and the variation of the field amplitude with the spot size, where


is the amplitude radius of the Gaussian spot in the - plane. Note that the last line of Eq. (10) is present only in the 3D case (), where . In the 2D case (, ) we have and the phase shift is .

Laser wavelength
Norm. amplitude
Spot size () ()
Pulse length () ()
Table 1: Laser parameters. The wavelength is scanned with respect to reference case (10 microns) as with scaling factors . The pulse has a Gaussian profile in time and space, and it is linearly polarized, consisting of and field components.

The Rayleigh length appearing in Eqs. (13)–(15) can be written in terms of the amplitude waist radius as


At a distance from the focal point , the mean laser intensity in 3D is 2 times lower than at the focus ( times lower in 2D).

As an example, Fig. 1 shows the -dependence of the pulse radius and the normalized amplitude for two wavelengths, and . The Rayleigh length is indicated by dotted lines at .

In the simulation, the laser is injected by applying an oscillating electromagnetic source on the left-hand boundary () of the simulation domain. The electric component of the source is chosen to be the sine (imaginary) part of in Eq. (10) evaluated at ,


The amplitude factor , the phase and the phase front delay are given by


where in 3D () and in 2D ().

Box size [mm] Resolution
() () ()
Table 2: Numerical parameters for 2D EPOCH simulations with moving window and weakly reflecting “open” or “simple outflow” boundaries. The length of the window is fixed, while its height is varied with the wavelength scaling factor , where corresponds to the reference case with . is the number of simulation particles (electrons) per cell. Ions are immobile. The values in parentheses on the bottom row were used in convergence tests.
Case Dims. Laser wavelength Focusing time Box size Resolution
(i) 2D ()
(ii) 2D ()
(iii) 2D () ()
(iv) 3D ()
Table 3: Numerical parameters for the 2D and 3D simulations with stationary window and “open” (EPOCH) or absorbing (REMP) boundaries. 2D simulations are performed for three different wavelengths, with scaling factors . The transverse box size is fixed while varying its length as . The rationale behind the chosen box size is discussed in the first part of A. The results shown in this paper were obtained using 1 PPC (simulation particles per cell) per species (electrons and deuterons). The case in 2D was also simulated with 10 PPC per species (in parentheses for case (iii)), which gave very similar result.

2.3 Parameters

The laser parameters for the 10 micron reference case are given in Table 1, where the focal spot size and pulse length are specified in terms of the respective full-width-half-maximum (FWHM) diameters of the intensity,


whose relation to and follows from


The values in Table 1 are chosen such that, even for the largest wavelength () in our parameter scans, it is ensured that and that the pulse still contains wave cycles. The constant time delay in Eq. (20) is used to ensure that the jump of the laser amplitude at the head of the laser pulse is sufficiently small, such that the associated numerical artifacts are negligible. Here, we choose for , where the increase of for larger is needed to compensate for the larger curvature of the wave front (second term in Eq. (20)).

Note that the FWHM pulse length of our Gaussian laser pulse corresponds to the amplitude length . This is shorter than the theoretically predicted optimum for the generation of a large-amplitude wakefield with a flat-top pulse: for BulanovSV16 (), which is for our parameters. Our choice of a shorter pulse length has a practical motivation as it reduces the energy required to reach the relativistic regime. A lower pulse energy, in turn, allows to reduce the cost of the laser or to increase its repetition rate.

We will focus primarily on the case where an external magnetic field with strength is applied in the -direction (); i.e., along the laser axis. Only a few results for the unmagnetized case and with a perpendicular field () are presented for comparison (Section 3.1, Fig. 2).

Short-time simulations with a moving window, which cover less than are performed with immobile ions. The ion response becomes important after a few , and we take it into account in simulations with a stationary window, where we follow the after-glow dynamics of the plasma wake channel for up to . In those cases, the ions are chosen to have the mass and charge of deuterium, so that initially.

The numerical parameters are summarized in Tables 2 and 3 for simulations with moving and stationary simulation windows, respectively. In the 2D simulations with a moving window, we fix the focal distance and the window length , while varying its height in proportion to the wavelength-dependent initial spot radius in Eq. (15) via the wavelength scaling factor . In the 2D and 3D simulations with stationary windows, we fix and, hence, use the same transverse box size in all cases, while varying the box length as in proportion to the wavelength-dependent focal distance, , which decreases with increasing wavelength .

For ease of comparison, snapshot times will be given in terms of


i.e., relative to the instant where the laser reaches its (theoretical) focal point . In the simulations with stationary windows, varies with and the values can be found in Table 3.

All results reported in this paper were obtained using only 1 particle per cell (PPC). Convergence tests in 2D indicate that this is sufficient for the scenario and phenomena discussed in this paper. The physical reason is that the plasma is collisionless and has a highly subcritical density, so individual particles do not interact with each other and the laser is not altered by the plasma’s presence. Moreover, we are presently concerned only with the overall large-scale structure of the plasma wake and its associated electromagnetic field.111Singular structures and associated phenomena are not well-resolved, as will be briefly discussed in Section 4.2 using Fig. 10.

We have performed convergence tests with respect to the spatial resolution and the number of simulation particles in 2D. The range of values scanned is shown in parentheses in Table 2 (bottom row) and Table 3 ( column). The results of these convergence tests are not presented here, since they are similar and readily reproducible using the EPOCH code Arber15 (). The spatial resolution for the 3D simulation (case (iv) in Table 3) was shown to be sufficient in 2D, so the 3D results are considered to be reliable for the low-density plasma at hand.

2.4 Codes

For the purpose of verification, we have performed our simulations with two relativistic particle-in-cell (PIC) codes, EPOCH Arber15 () and REMP Esirkepov01 (). Details about the equations solved and numerical methods used can be found in the respective references for each code.

The 2D results reported in Sections 3 and 4 were obtained with EPOCH using “open” or “simple outflow” boundaries, which gave similar results: the radially expanding (defocusing) laser pulse was partially reflected, but with a negligible effect as its amplitude was attenuated by a factor when it passed again through .222For reflected pulses, see Fig. 12 in Section 4.4. Preliminary studies were performed using EPOCH versions 4.9.2 and 4.11.0, which gave similar results to version 4.14.4 that was used here. EPOCH is freely available and all simulations should be easily reproducible with the input.deck files provided in this paper’s meta data. Most of the 2D simulations, especially with laser wavelengths , can be completed within a few hours using a few 100 CPU cores.

The 3D results presented in Section 5 were obtained with REMP using open boundaries supplemented with an absorbing layer of thickness in all directions, which damps all electromagnetic fluctuations with a super-Gaussian factor after the laser has entered the simulation domain.

REMP and EPOCH were benchmarked against each other and gave similar results in long-time 2D simulations (up to ) and in short-time 3D simulations (). At later times (after a few ), we encountered discrepancies between the 3D results of REMP and our adaptation of EPOCH. The reason for these discrepancies has not been found yet, so we chose to report only the 3D results of our in-house code REMP in this paper. Note that the standard EPOCH code loads the quasi-particles in a uniform but random fashion. For the present study, we have modified EPOCH such that the quasi-particles are loaded “uniformly” (i.e., at equidistant positions in regions where the density is constant), as we usually do with REMP. This makes some discretization effects clearly visible and was helpful when making code-to-code comparisons. EPOCH was run with the default 1st-order b-spline (triangular) shape function,333In the present scenario, 2D EPOCH simulations using a 3rd-order b-spline shape function took 30% longer and gave results that were very similar to those obtained with the default 1st-order b-spline. whereas REMP uses a quadratic form.

Figure 2: Snapshots of the electron density in the wake of a laser pulse with and in 2D EPOCH simulations with moving window and immobile ions (parameters in Tables 1 and 2 with scaling factor ). Results are shown for cases (a) without external magnetic field, (b) with an axial magnetic field , and (c) with a perpendicular out-of-plane magnetic field . The laser has entered from the left and propagates to the right. The location and width ( radius of the amplitude) of the laser pulse is indicated by a dotted yellow ellipse. The black dotted curves indicate the local spot diameter of the laser (cf. Fig. 1(a)), and the vertical dashed lines indicate the location of the focal point (). Note that these snapshots show only a portion of the simulation domain and the contour plots were averaged over grid points, so that the density spikes are smoothed out on that scale.

3 Evolution of a magnetized laser plasma wake channel at low density in 2D

3.1 Effect of magnetization in the vicinity of the laser

We begin with a short introduction of the laser-driven electron dynamics in the parameter regime of interest. For illustration, we use the 2D simulation results in Figs. 2 and 3, which show snapshots of the perturbed electron density behind a 2 picosecond laser pulse with a wavelength of 10 microns that propagates through a sparse gas with low electron density . The laser parameters are given in Table 1 (with scaling factor ) and the numerical parameters are listed in Table 2. The laser was followed for about with a moving window covering () along the direction. Parallelized over 1,280 CPU cores, each simulation took about 13 hours on the supercomputer JFRS-1 jfrs1 ().

Figure 2 shows the structure of the laser wake in three cases: (a) without external magnetic field, , (b) with an axial field directed along the laser path, and (c) with a perpendicular out-of-plane field . The snapshots were taken at a time where the laser had traveled about past its focal point.

In Fig. 2(a) one can see the structure of the unmagnetized laser wake in the relativistic regime, which consists of an electron-free positively charged cavity and a detached bow wave Esirkepov08 (). Immediately after the laser has passed, the weakly perturbed electrons near the cavity boundary perform approximately radial plasma oscillations with angular frequency . In agreement with theoretical predictions for Akhiezer56 (); Tajima79 (); BulanovSV16 (), this produces a wake with characteristic wavelength


behind the laser pulse moving away with velocity . However, for the present combination of parameters, only the first half cycle is clearly visible in the form of a cavity because transverse wave breaking BulanovSV97 () dominates and overshadows the subsequent wakes, which have a D-shaped form but are not clearly visible in Fig. 2(a).

Panels (b) and (c) of Fig. 2 show that adding an external magnetic field with a strength of (and otherwise identical parameters) gives rise to a force that is sufficiently strong to confine the bow waves near the wake channel and cause repeated transverse wave breaking, so that the laser pulse leaves behind a trail of multiple wakes and cavities. In the case with an axial magnetic field, Fig. 2(b) shows that the oscillations possess nearly perfect up-down symmetry. The oscillation frequency is consistent with that predicted theoretically for an extraordinary electromagnetic wave near the upper hybrid (UH) resonance BulanovSV13 (),


which produces a magnetized wake with characteristic wavelength


behind the laser pulse moving away with velocity . When the magnetic field is perpendicular to the laser axis, the two modes of oscillation — and — exist independently, which produces asymmetric flow patterns, as can be seen in Fig. 2(c) for an out-of-plane magnetic field that causes counter-clockwise electron gyration. Since the snapshots in Figs. 2(b) and 2(c) were taken during the relativistic phase, the oscillation period is somewhat obscured by the multi-flow structures associated with the repeated transverse breaking of wake and bow waves. The wave cycle can be seen more clearly in weakly perturbed (nonrelativistic) regions away from the focal point, as will be briefly discussed in Section 3.2 below for the case with axial magnetic field (for instance, see Fig. 3(a)).

Figure 3: Evolution of the electron density in the wake of a laser pulse in the same 2D case as in Fig. 2(b) above, with and an axial magnetic field, . Each snapshot is arranged as in Fig. 2, with the addition of a red dashed rectangle indicating the estimated length of the region, where the laser is sufficiently intense to carve out an electron-free cavity (cf. Eq. (31)). In the present case, we have , which is only partially visible in each panel. Snapshot (e) is the same as in Fig. 2(b).

One physical consequence of the smaller radial excursion of the magnetized bow wave is that the induced electric field is weaker, while the induced current density and associated magnetic perturbations are stronger than in the unmagnetized case. From the viewpoint of numerical modeling, the magnetic confinement of the transversely accelerated electrons near the wake channel is advantageous, because it can prevent the bow wave from reaching the unphysical boundaries of the simulation box in the transverse direction, as happened in the unmagnetized case in Fig. 2(a). The minimal transverse size of the simulation box is then determined only by the initial spot size of the unfocused pulse at () (cf. Fig. 1(a)).

Figure 4: Early phase () of the evolution of the magnetized plasma wake in three cases: (i) , (ii) , and (iii) (i.e., ). For each case, three snapshots (a)–(c) of the electron density are shown in close succession, apart, while the laser pulse (yellow ellipse) propagates away from its focal point (). The black dotted curves indicate the local spot radius of the laser pulse, and the red dashed rectangle indicates the cavitation length . For snapshot (a) of each case, panel (d) shows the transverse profile of the electron density measured at . (2D EPOCH simulations with stationary window and mobile ions, using the parameters in Tables 1 and 3.)

As we have already noted in the introduction, 2D simulations cannot fully capture the true dynamics of laser plasma wakes in the presence of an external magnetic field. At first glance, the choice of an external field along the direction as in Fig. 2(c) may seem to be more easily justified for 2D simulations than the use of an in-plane field, such as in Fig. 2(b). However, with an out-of-plane field, the particle motion is unmagnetized in the third dimension, so that, in reality, electrons would be pulled into the wake channel along the -axis and perform plasma oscillations when the laser pulse has a finite extent in that direction. Since this effect is absent in 2D simulations with an out-of-plane magnetic field, such a setup is not useful for our purpose of investigating the long-time evolution of the plasma wake channel, where parallel streaming of charged particles along the magnetic field is expected to be a key factor that determines the channel’s evolution and life time.

Figure 5: Long-time evolution () of the electron density (left) and the deuteron density (right) in the three cases (i)–(iii) shown in Fig. 4 above.

Thus, in the following, we will consider only cases with an axial magnetic field as in Fig. 2(b) with , which will serve us as a reference case. The structures seen in Fig. 2(b) may be viewed as a distorted planar projection of the 3D dynamics, where the charged particles gyrate around the -axis in a helical manner. The 2D code evolves the particle momenta in all three directions, but the motion along is not executed, and we have , where is the fluctuating part of the magnetic field. In any case, regardless of where is pointing, the results of 2D simulations must always be interpreted with great care.

3.2 Length of the electron-free cavity

For the same case as in Fig. 2(b), the series of snapshots in Fig. 3 shows how the structure of the magnetized wake behind the laser pulse changes as the laser focuses and defocuses. In snapshots (a,b,g,h,i), where the laser intensity is low, one can clearly see the theoretically predicted wavelength of the wake given by Eq. (27). In the other snapshots (c–f), where the laser intensity is relativistic, the wave period is obscured by multi-flow structures as we already saw in Fig. 2.

Figure 6: Long-time evolution of the transverse profiles of the electron density (left) and the deuteron density (right) at the focal point () in each of the three cases (i)–(iii) shown in Fig. 5. For noise reduction, the profiles were spatially averaged over a distance of .

Figure 3 shows that the laser is able to carve out an electron-free cavity in a region extending approximately on both sides of the focal point . According to the predicted evolution of the 10 micron laser pulse in Fig. 1 above, this corresponds to the region where in 2D () or in 3D (). Based on this observation, the length of the domain where a cavity forms may be estimated by substituting the condition


into Eq. (10), which yields




for . We will refer to this quantity as the “cavitation length”. The exponent in our definition of in Eq. (28) implies that the cavitation length is the same in 2D and 3D when . This empirical choice seems to be valid at least for the scenario studied here, as will be confirmed later by the 3D simulation reported in Section 5. The empirical value on the right-hand side of Eq. (28) is somewhat arbitrary and may be replaced by another value of similar magnitude. In our reference case with and , where , the formula in Eq. (30) yields the cavitation length


which is indicated in Fig. 3 by a red dashed rectangle with boundaries at , which is only partially visible due to its large size.

Note that the length of the cavity can be shorter than the cavitation length , as in the unmagnetized case shown in Fig. 2(a). Furthermore, note that the length of the plasma channel evolves in time and its boundaries will deviate from our estimate as particles enter the cavity by streaming along the magnetic field.

3.3 Effect of the laser wavelength

It is clear that the cavitation length varies with the laser wavelength via the Rayleigh length , which appears as a factor in Eq. (30). In the following, we compare the magnetized wake dynamics induced by relativistic laser pulses with different wavelengths. All other parameters are the same as in Table 1.

For each of the three cases (i)–(iii) in Table 3, Fig. 4 shows three snapshots (a)–(c) of the electron density during the first after the laser pulse has passed its focal point. The laser’s amplitude radius is indicated by a yellow dotted ellipse. While the contour plots provide an image of the overall structure of the density field, more quantitative information about the density perturbations can be gleaned from panel (d) in Fig. 4, which shows the transverse profile of the electron density at for snapshot (a) taken at in each case.

The results in Fig. 4 show that laser pulses with different wavelength, here (), produce similar wake channels. One can also see that the starting point of the strongly perturbed region is consistent with the cavitation length (red dashed rectangles) estimated by Eq. (30).

About - after the formation of the electron-free cavity, the ion response becomes visible in the contour plots of shown in Fig. 5(right) and in the corresponding transverse profiles in Fig. 6(right). Owing to their large inertia, the ion motion becomes the dominant factor determining the further evolution of the plasma channel thereafter. During the interval shown in Figs. 5 and 6, the ion density also develops a pronounced cavity in the original cavitation region of the electron density, (red dashed rectangle). At this stage, the electrons in Fig. 5(left) seem to mostly follow the relatively slow ion motion.

The plasma channel grows in length along the axis, far beyond the initial cavity boundaries, first towards the left () and with a small delay also towards the right (). The diameter of the elongated channel remains similar to the diameter of the primary cavity near the focal point, which is a clear indication of the fact that the axial elongation is due to the particles flowing along the magnetic field.444If the axial elongation of the channel was a direct but retarded effect of the laser pulse, we would expect it to become wider with increasing distance from the focal point, but this is not the case here.

4 Limited similarity of magnetized wakes in 2D

The results in Figs. 46 indicate that the overall structure and temporal evolution of the magnetized plasma wake channel in the parameter regime considered here is largely independent of the laser wavelength . We interpret this as a manifestation of so-called limited similarity or limited scaling (see Ref. Esirkepov12 () for a review of these concepts). In this section, we examine in more detail to what extent the 2D dynamics are similar.

Figure 7: Evolution of the laser pulse in 2D with wavelengths scaled by factors , , , and . (a): Peak amplitude of the laser’s magnetic field multiplied by , which is proportional to the normalized amplitude in Eq. (9). (b): spot radius. (c): FWHM pulse length. All curves are plotted as functions of the position of the Gaussian pulse center. The symbols are the simulation results (2D EPOCH with moving window) and the solid lines are the analytical solutions (10) and (15) of the paraxial wave equation, with according to Eq. (9).

4.1 Preliminary considerations

Two systems are guaranteed to behave identically if their dynamics are governed by the same equations and if all dimensionless parameters have the same values. In contrast, limited scaling requires only that the same phenomena dominate in the simulation “and in nature, i.e. dimensionless quantities in nature which are small compared to unity should be small in the model, but not necessarily by the same order of magnitude” Block67 (). In our case, the small key parameter is the ratio of the initial electron density to the critical density,




is the density for which the laser frequency equals the electron plasma frequency given by Eq. (2).

Figure 8: Contour plots of (I) the electron density field (left) and (II) the out-of-plane magnetic field fluctuations (right) in 2D simulations with axial magnetic field and different wavelengths (a)–(e). The laser pulse has entered from the left-hand side and all snapshots are taken at the time where the pulse center is located at (yellow ellipses), a few millimeters beyond its focal point (). The black dotted curves indicate the local spot diameter of the laser (cf. Fig. 7(b)). In (d) and (e), the boundaries of the cavitation domain, (cf. Eq. (31)), lie within the plotted range and are indicated by a dashed red rectangle. In (e), that rectangle was shifted by to match the actual location of the focal point inferred from Fig. 7(b). The magnetic fluctuation amplitude was clipped at , so the magnetic field in the vicinity of the laser appears as a saturated dark area in the contour plots on the right-hand side. The vertical stripes in the contours seem to be harmless post-processing artifacts of the EPOCH code and coincide with the MPI domain boundaries along . (2D EPOCH simulations with moving window and immobile ions using the parameters in Tables 1 and 2.)

As explained in Sections 1 and 2 above, we intend to apply our simulations to magnetically confined fusion plasmas, in particular tokamaks. Such plasmas tend to have a relatively flat density profile and are constrained by an empirical limit known as the “Greenwald density limit” Greenwald02 (), which implies that the electron density near the plasma surface cannot be much higher than . The reference density of in our simulations was chosen to lie somewhat below this limit. From Eq. (33) one can readily see that this is at least 5 orders of magnitude smaller than the critical density for typical high power lasers, whose wavelengths lie between (for solid-state Ti:sapphire lasers) and (for gas lasers). This means that the condition in Eq. (32) will be satisfied even if we increase the laser wavelength in the simulation far above the original reference value by a large scaling factor as


In the following subsections 4.24.4, we compare in detail the magnetized plasma dynamics induced by relativistic laser pulses with wavelength and scaling factors


All other parameters are the same as in Table 1. For these values, the chosen pulse length is sufficient to accommodate a reasonable number of wave cycles, which is for .

Figure 9: Contour plots of the and components of the induced current density (top) and the electric field strength (bottom) in the cases with laser wavelengths from Fig. 8. (Results of 2D EPOCH simulations with moving window and immobile ions.)

Before examining the plasma dynamics, it is useful to take a glance at the evolution of the laser pulse in each case, which is shown in Fig. 7 as a function of the position of the Gaussian pulse center. Figure 7(a) shows the peak amplitude of the laser’s magnetic field scaled by , which is proportional to the normalized amplitude in Eq. (9). Figure 7(b) shows the spot radius (defined in terms of the amplitude), which drops to at the focal point in all cases. The simulation results (symbols) are close to the respective analytical solutions (solid lines) of the paraxial wave equation for a Gaussian laser pulse propagating through vacuum; namely, Eqs. (10) and (15). There is, however, a small but noticeable difference between the theoretically estimated focal point (defined to be at ) and the simulated one: in the simulation, the laser reaches its focus slightly before , and that difference seems to increase with increasing wavelength .

Figure 7(c) serves as an additional test for numerical accuracy as it shows the constancy of the pulse duration throughout the simulation. Note that the pulse length plotted in Fig. 7(c) was simply measured from the profile of the oscillatory field magnitude, , so it tends to be underestimated and fluctuates by some value of magnitude .

4.2 Demonstration of limited similarity in 2D

We consider the time slice where the laser has traveled about past the focal point located at . Figure 8 shows color-shaded contour plots of the spatial structure of (I) the electron density in the left column and (II) the fluctuating component of the out-of-plane magnetic field in the right column (here since is directed along ). One can see that the structure of the electron density fluctuations and of the induced magnetic perturbations is similar in the cases with (). The same can be said about the structure of the induced current density and electric field , whose and components are shown in Fig. 9.

The case in row (d) of Figs. 8 and 9 is marginal because the cavitation length — which is according to Eq. (31) and indicated by a dashed red rectangle in Figs. 8(d) and 9(d) — can still accommodate nearly two cycles of the magnetized wake, whose characteristic length is approximately (Eq. (27)). In contrast, in the case in Figs. 8(e) and 9(e), the cavitation length has shrunk to roughly the same size as the wake length or shorter, , so that similarity is lost for the present choice of parameters.555As is clear from Eqs. (26) and (27), the wake length may be reduced by increasing the strength of the external field or increasing the density .

Note that the data shown in rows (a)–(c) of Figs. 8 and 9 have been filtered and down-sampled to the same resolution as in the case in row (d). Originally, the electron density spikes and cusps in the case are much narrower and about 5 times higher () than in the case ().

Figure 10 shows a zoom-up of the density spike located near the laser pulse, where the bow wave detaches from the wake. Recently, this structure has attracted attention since it may serve as a practically useful source of coherent energetic photons via a mechanism dubbed Burst Intensification by Singularity Emitting Radiation (BISER) Pirozhkov17b (); Pirozhkov12 (); Pirozhkov14 (). This structure is a true (albeit integrable) singularity, which means that its width in the simulation varies with the spatial resolution used. This is demonstrated in Fig. 10, where panels (a) and (c) show the singularities for the and cases with the default resolution of ( points per wavelength), and panel (b) shows the case simulated with 4 times higher resolution (); i.e., with the same grid spacing as used in the case. One can see that the singularity in (b) is sharper than that in (c) and it is located a little closer to the laser, whose amplitude radius is indicated by the yellow dotted curve in the bottom right portion of each panel.

Note that there are differences between the equally resolved singularities in panels (a) and (b). For instance, one can observe a modulation of both the wake and the bow wave on the scale of the laser wavelength, and the amplitude of this modulation is larger in Fig. 10(b) than in (a). This observation implies that singularity-related phenomena such as harmonic generation by the BISER mechanism Pirozhkov17b (); Pirozhkov12 (); Pirozhkov14 () depend at least quantitatively on the laser wavelength . Meanwhile, the required resolution and, hence, the computational effort required to study these phenomena is not determined by the laser wavelength but by the singular structures, which should be resolved in as much detail as possible.

Next, let us inspect the electron energy distribution , where is the kinetic energy of an electron with rest mass and Lorentz factor . Figure 11(a) shows the spatially averaged energy distribution for the same snapshot time as in Fig. 8. One can see that the form and amplitude of is similar for wavelengths and energies in the range .

The differences at energies below may be explained as follows. The different values of the first sample () in Fig. 11(a) can be attributed to the fact that, for a fixed focal spot size , lasers with longer wavelengths perturb a larger volume of the gas in regions away from the focal point, where their spot size is larger (cf. Fig. 7(b)). The perturbation in those regions is, of course, very weak, which is why the resulting differences in are seen only at the lowest energies. To some extent, the differences at low energies are also due to the fact that the height of the simulation box and, thus, the number of weakly perturbed low-energy electrons is varied in proportion to the laser wavelength as while keeping . (cf. Table 2). Here, this has a noticeable effect in the case with for energies . As shown in Fig. 11(a), the differences are reduced when the box height is increased from to .666As explained in A, a larger transverse box size may be needed to reduce the accumulation of artifacts in the electromagnetic fields; in particular, for shorter wavelengths , because the laser converges only slowly due to its longer focal length, so that its outer portions interact with the artificial boundaries for a longer time. While this is not critical for the simulations with moving window discussed in the present section, more care is required when choosing the transverse size of the simulation box in long-time simulations with stationary window.

Figure 10: Zoomed-up view of the density singularity near the laser pulse. Panels (a) and (c) show data from the same simulations as Figs. 8(b) and 8(d) for the and cases, respectively, albeit at a slightly earlier time where the laser is closer to its focal point (). Panel (b) shows the results of the case simulated with 4 times higher resolution in both and ; i.e., the same resolution as was used by default in the case. The white dotted curves indicate the local spot radius of the laser (cf. Fig. 7(b)) and the partially visible dotted orange ellipse is the laser’s amplitude radius. (Results of 2D EPOCH simulations with moving window and immobile ions.)
Figure 11: Electron energy distribution in arbitrary units (a.u.). Panel (a) shows the energy distributions for the cases with at the same snapshot time as in Fig. 8. For (blue circles) results are shown for two different heights of the simulation box: (dashed) and (solid). The inset at the bottom of panel (a) shows the structure of the (noisy) high-energy tails, . Panels (b) and (c) show zoom-ups of the energy distribution in the case simulated with the default grid points per wavelength (solid) and with 4 times higher resolution, (dashed) as in Fig. 10. Notable resolution-dependent differences in were found only near and . (Results of 2D EPOCH simulations with moving window and immobile ions.)

Since we know from Fig. 10 that there are underresolved singular structures in the present simulations, we have also checked whether the spatial resolution affects the energy distribution in Fig. 11 in the case. Noteworthy differences were found only near and , and zoom-ups around these points are shown in panels (b) and (c) of Fig. 11: with 4 times higher resolution (dashed line), the “knee” near in (b) is slightly shifted towards higher energies, and in (c) one can see that the high-energy “tail” around is somewhat longer. Since these resolution-dependent differences are small, we have not tried to identify their reasons. They may be related to the structural differences at small scales, such as those seen in Fig. 10(b) and (c).

Finally, we discuss the high-energy tail. At (), the energy distribution in Fig. 11(a) drops abruptly by about one order of magnitude. This cut-off lies well below the theoretical maximum for the acceleration of a single electron in the 1D plane wave limit, which is


for acceleration due to the laser field (e.g., see p. 327 in Ref. Mourou06 () and p. 231 in Ref. BulanovSV_Book01 ()). Nevertheless, the fact that the present cut-off in Fig. 11 does not depend on the wavelength of the laser suggests that this is, in fact, the effective maximal electron energy that can be achieved with the present paraxially focused Gaussian pulse with FWHM spot diameter and pulse length . Indeed, it is reasonable to assume that some electrons can be accelerated to full speed in all our cases, because even with the Rayleigh length is an order of magnitude larger than the theoretically predicted distance over which full acceleration takes place, which is given by (again in the 1D limit)777Equations (37) and (38) determine the distance that an electron travels along the laser path until it reaches its maximal perpendicular and parallel momentum, and , respectively, assuming that this acceleration takes place while the electron experiences half a period of the laser field.


and becomes for .

In other words, our in Fig. 11 is independent of because electrons can be accelerated to full speed well within the time interval (or distance), where the laser pulse is maximally focused. We consider this to be one of the necessary conditions for limited similarity, which are discussed in Section 4.3 below.

A clear dependence of on can be seen in Fig. 11(a) only beyond the cut-off, where a noisy high-energy tail extends to higher energies, reaching nearly in the case. The inset panel in Fig. 11(a) shows the structure of these noisy tails for three cases with . We have not investigated this further, because the number of electrons in the noisy tail is very small, so we do not expect them to have any significant influence on the overall dynamics of the plasma wake channel.

4.3 Conditions for insensitivity w.r.t. laser wavelength

The results in Figs. 8, 9 and 11 discussed in the previous section show that the overall structure and dynamics of the magnetized plasma wakes in the vicinity of the laser’s focal point are fairly independent of the laser wavelength for . It is also reasonable to assume that this limited similarity holds for shorter wavelengths below . Of course, the parameter regime where such similarity can be found has some constraints, besides the key condition in Eq. (32) that was already pointed out in Section 4.1.

One constraint that limits the range of values for the wavelength scaling factor in Eq. (34) arises from the fact that, according to Eq. (30), the cavitation length scales with the wavelength as , so that the initial length of the plasma channel becomes shorter with increasing . However, the channel should not be shorter than the typical size of the structures that form in the magnetized plasma wake of the laser (cf. Eq. (27)). In our case, this condition was violated for .

In order to ensure similarity of the bow waves as well as the electron energy distributions , we require that some electrons must be accelerated to the maximally possible kinetic energy (Eq. (36) in 1D) while the pulse is fully focused. This mea