Effects of simulation dimensionality on laser-driven electron acceleration and photon emission in hollow micro-channel targets

The interaction of high-power ultra-intense lasers and structured (both nanostructured Rajeev et al. [2004]; Nishikawa et al. [2004]; Andreev et al. [2011]; Margarone et al. [2012]; Blanco et al. [2017]; Moreau et al. [2019]; Vallières et al. [2021] and microstructured Naumova et al. [2004]; Klimo et al. [2011]; Xiao et al. [2016]; Gozhev et al. [2020]; Feng et al. [2018]; Jiang et al. [2016]; Ji et al. [2016]; Stark et al. [2016]; Yu et al. [2018]; Snyder et al. [2019]; Bailly-Grandvaux et al. [2020]) targets has been a topic of great interest for its capability of enhancing the laser energy conversion efficiency Rajeev et al. [2004], high-order harmonics generation Lavocat-Dubuis and Matte [2010]; Cantono et al. [2018], charged particles (relativistic electrons and ions Andreev et al. [2011]; Margarone et al. [2012]; Bailly-Grandvaux et al. [2020]) acceleration and the production of X-ray Rajeev et al. [2004]; Yi et al. [2016]; Lècz and Andreev [2017]; Rolles [2018] to $\gamma$-ray Stark et al. [2016]; Yu et al. [2018]; Wang et al. [2018] radiation. The produced charged particle and photon beams have a wide range of applications from medical ion therapy Malka et al. [2004]; Robson et al. [2007], nuclear physics Bychenkov et al. [1999]; Bertozzi et al. [2008] to photon-photon pair production Yu et al. [2019]; Wang and Arefiev [2020]; Wang et al. [2020]. The microstructured targets with characteristic size of surface modulation comparable to laser wavelength are able to extensively absorb laser energy through various processes, including surface plasmon resonance excitation Rajeev et al. [2003], multipass stochastic heating Breizman et al. [2005] in dense clusters, prolonged acceleration distance in hollow channels Xiao et al. [2016] and microwires Gozhev et al. [2020] and relativistic transparency in prefilled channels Stark et al. [2016]. In this paper, we examine the regime involving hollow micro-channels (see Fig. 1).

When ultra-intense lasers irradiate hollow micro-channels, strong laser fields directly act on electrons, dragging them into the channel and forming periodic electron bunches which then surf along with the laser pulse, gaining energy from the laser. Additionally, the presence of a channel guides the propagation of the electromagnetic fields, confines the electron motion and as a result, leads to a well collimated photon emission. This setup can serve as a promising electron source to further stimulate ion acceleration Bailly-Grandvaux et al. [2020] as long as the ion expansion does not significantly impact electron acceleration Wang et al. [2019]. However, to successfully apply such an electron source in experiments, careful numerical investigations are needed in order to determine where the electron energy peaks, i.e. the location to cut off the channel and collect an electron source with an optimal spectrum.

To carry out such numerical studies, one can choose between 2D3V and 3D3V Particle-In-Cell (PIC) simulations. Both 2D Klimo et al. [2011]; Xiao et al. [2016]; Naumova et al. [2004]; Zou et al. [2017]; Gong et al. [2019] and 3D Ji et al. [2016]; Jiang et al. [2016, 2014]; Snyder et al. [2019]; Serebryakov et al. [2019] numerical simulations have been widely used to characterize laser interactions with structured targets. The appeal of 2D simulations is that they require significantly less computational resources than 3D simulations, so one is able to perform extensive parameter scans using 2D simulations. However, the field topology differs between 2D and 3D setups and it is not immediately clear how the differences impact the particle dynamics. It is then important to evaluate the dimensionality effects on a case-by-case basis. A few publications Stark et al. [2017]; Liu et al. [2013]; Xiao et al. [2018]; Blanco et al. [2017] have discussed the influences of simulation dimensionality on ion acceleration with various target geometries. But, to our knowledge, nobody has examined and explained the physics of dimensional effects on laser-irradiated hollow micro-channel targets.

In this paper, we show that the chosen dimensionality has a considerable effect on electron acceleration and the associated photon emission. First, we provide theoretical analysis to demonstrate that the dephasing rate between the accelerated electron and laser wave-fonts strongly depends on simulation geometry. Later we show numerical evidence to demonstrate that the dephasing rate differs with simulation dimensionality and such a difference is the key reason for the distinguished observation in terms of electron and photon beam generation. Through collectively evaluating generated particles and detailed particle tracking, we show that the occurrence of high phase velocity in the 3D setup terminates electron acceleration process early in space and time, and leads to a reduction of photon emission.

Previous studiesXiao et al. [2016]; Wang et al. [2019]; Gong et al. [2019] have shown that, in a hollow channel without significant ion expansion, the dominant contributor to electron acceleration is the electric field in the direction of the laser propagation, i.e. the longitudinal field $E_{x}$. The corresponding structure of propagating electromagnetic fields can be characterized as TM (transverse magnetic) modes. The wave equation of TM modes can be written as

$$\bigg{(}\dfrac{\partial^{2}}{\partial y^{2}}+\dfrac{\partial^{2}}{\partial z^{2}}\bigg{)}E_{x}+\bigg{(}\dfrac{\omega^{2}}{c^{2}}-k^{2}\bigg{)}E_{x}=0,$$

(1)

where $\omega$ is the wave frequency, $k$ is the wavenumber and $c$ is the speed of light. For a two-dimensional waveguide (shown in Fig. 1a), $E_{x}$ is a function of $y$ and $\partial E_{x}/\partial z=0$. Equation 1 then becomes ${\partial^{2}}E_{x}/{\partial y^{2}}+(\omega^{2}/c^{2}-k^{2})E_{x}=0$. We choose $E_{x}=E_{0}sin(\pi y/R)$ as the TM-mode solution that matches the field structure in the incoming pulse (i.e. $E_{x}=0$ on axis) and the boundary conditions which requires $E_{x}(y=\pm R)=0$. Here $R$ represents the radius of the plasma channel. The dispersion equation in the 2D waveguide is then given by

$$\dfrac{\omega^{2}}{c^{2}}=k^{2}+\dfrac{\pi^{2}}{R^{2}}.$$

(2)

For the 3D cylindrical waveguide (shown in Fig. 1b), it is convenient to rewrite Eq. 1 in cylindrical form as

$$\dfrac{1}{r}\dfrac{\partial}{\partial r}\bigg{(}r\dfrac{\partial E_{x}}{\partial r}\bigg{)}+\dfrac{1}{r^{2}}\dfrac{\partial^{2}E_{x}}{\partial\psi^{2}}+\bigg{(}\dfrac{\omega^{2}}{c^{2}}-k^{2}\bigg{)}E_{x}=0,$$

(3)

where $r$ is the axial distance and $\psi$ is the azimuth. Assuming a solution in the form of $E_{x}=f\sin(\psi)$, the resulting equation for $f$ then reads

$$s^{2}\dfrac{\partial^{2}f}{\partial s^{2}}+s\dfrac{\partial f}{\partial s}+(s^{2}-1)=f.$$

(4)

Here $s\equiv\beta r$ and $\beta^{2}=\omega^{2}/c^{2}-k^{2}$. The solution for $f$ is given by $f=E_{\parallel}J_{1}(s)=E_{\parallel}J_{1}(\beta r)$ where $J_{1}(x)$ is a Bessel function of the first kind and $E_{\parallel}$ is the amplitude of the longitudinal field. At the boundary of the cylindrical waveguide, $f(r=R)=0$ which yields $\beta R\approx 3.8$. The dispersion relation of the 3D cylindrical waveguide is then written as

$$\dfrac{\omega^{2}}{c^{2}}=k^{2}+\dfrac{14.7}{R^{2}}.$$

(5)

It is convenient to derive a general expression for $v_{ph}$ of a propagating wave with the dispersion relations given in Eqs. 2 and 5.

$$u=\dfrac{v_{ph}}{c}=\dfrac{\omega}{kc}=\sqrt{1+\dfrac{\alpha^{2}}{R^{2}k^{2}}}\approx\sqrt{1+\dfrac{\alpha^{2}}{R^{2}k_{0}^{2}}},$$

(6)

where we set $k\approx k_{0}=\omega/c$, $\alpha^{2}=\pi^{2}\approx 9.9$ for the 2D channel and $\alpha^{2}=14.7$ for the 3D cylindrical channel. Subtracting the $v_{ph}$ from $c$, in the limit of $\alpha^{2}\ll R^{2}k_{0}^{2}$ we have

$$\delta u=u-1\approx\dfrac{1}{2}\dfrac{\alpha^{2}}{R^{2}k_{0}^{2}}.$$

(7)

$\delta u$ is a dimensionless parameter used to quantify the degree of supuerluminosity Robinson et al. [2015]; it can be understood as a dephasing rate, since it illustrates how quickly the local laser wave-front outpaces the electron in question. The ratio of the dephasing rate for the case in 2D to that in 3D is

$$\dfrac{\delta u_{3D}}{\delta u_{2D}}\approx 1.5.$$

(8)

In the following sections we are going to demonstrate via numerical simulations that the lower dephasing rate in 2D leads to an overestimate of electron acceleration and $\gamma$-ray emission.

We model the electron acceleration through fully relativistic PIC simulations using the EPOCH codeArber et al. [2015] which is available in both 2D and 3D. The target geometries implemented in the simulations are depicted in Fig. 1. In 2D, the target is a straight empty channel enclosed by two uniform plasma slabs. In 3D, the channel is enclosed by a cylindrical wall of plasma, with a channel diameter such that a slice along the target axis would be identical to the 2D setup. The 2D and 3D simulations share the same plasma composition, plasma density, laser intensity, temporal profile, and laser spot size. The laser intensity is set as $1.37\times 10^{22}\ \rm{W/cm}^{2}$, corresponding to $a_{0}=100$. Here $a_{0}\equiv|e|E_{0}/(m_{e}c\omega)$ is the normalized laser amplitude, where $E_{0}$ is the peak amplitude of the electric field in the incoming laser pulse, and $m_{e}$ and $e$ are the electron mass and charge. The laser pulse is always focused at the channel entrance. We choose gold as the original target material with a density of 1.5 g/cm³. According to the field ionization model, the considered laser pulse is capable of ionizing gold atoms to the level of $Z=69$. In the simulations, the target is pre-ionized accordingly to a plasma composed of Au⁺⁶⁹ (density 4 $n_{cr}$) and $e^{-}$ (density 276 $n_{cr}$) where $n_{cr}\equiv m_{e}\omega/(4\pi e^{2})$ is the critical density. A detailed comparison of parameters used in the 2D and 3D simulations is listed Table 1. Note that although the target lengths are set differently, this is done such that in different dimensionalities the target is long enough for the electrons to reach their first energy peak.

Figure 2(a) and (b) illustrate the time history of the electron energy spectrum observed in both 2D and 3D simulations. There exists more than one energy peak in both Fig. 2(a) and Fig. 2(b) and our focus is on the first energy peak which happens before any deceleration takes effect. It is clear that the maximum electron energy gain achieved in 2D significantly exceeds the gain observed in the 3D case. In fact, the first energy peak in 2D occurs at $t=500$ fs with $\varepsilon_{e}=1920$ MeV while in 3D the peak occurs at $t=255$ fs with $\varepsilon_{e}=1240$ MeV. Figure 2(c) gives a direct comparison of the peak spectra observed in the 2D and 3D simulations. In order to facilitate a comparison between the two simulations, we added a gray star to Fig. 2(a) that represents the energy and time of the first peak from Fig. 2(b). Evidently, the acceleration in 2D lasts longer, which results in a more energetic electron spectrum. Note that to make a quantitative comparison with the 3D simulations, we assume a uniform third dimension with a length of $2~{}{\upmu\mathrm{m}}$ for the 2D simulation.

To further understand the electron acceleration process, we tracked energetic electrons in both simulations. Figure 2(d) and Figure 2(e) illustrate the trajectories of representative electrons selected from both the 2D and 3D simulations. Plotted in the space of the transverse location ($y$ or $r$) and time, the trajectories are in agreement with the time history of electron spectrum given in Fig. 2(a-b); the electron energies peak at the corresponding moments. Regardless of the dimensionality, the electrons surf along the channel wall while getting accelerated by longitudinal electric fields until reaching their first energy peak. However, the horizontal surfing of electrons in the 3D simulation terminates earlier due to its higher dephasing rate, which will be elaborated in the next section. It is worth noting that the second energy peak observed in Fig. 2(a-b) is correlated with the electron motion of crossing the central axis, indicating an involvement of the transverse electric field in electron acceleration.

Through evaluating the electron energy spectrum and tracking individual electron motion, we have shown that the dimensionality of simulations significantly impacts electron acceleration. The 2D simulations tend to extend the electron acceleration process, leading to an overestimate of the maximum electron energy when compared to more realistic 3D simulations.

An electron travelling in the laser fields gains energy only while staying in the accelerating phase of the electric field. The electron can gain energy from both longitudinal ($E_{\parallel}=E_{x}$) and transverse ($E_{\perp}=E_{y},\ E_{z}$) electric fields. The total work done by the electric fields on a given electron can be expressed as

$$W_{tot}=W_{\parallel}+W_{\perp}=-|e|\int_{-\infty}^{t}(E_{\parallel}v_{\parallel}+E_{\perp}v_{\perp})\,dt{{}^{\prime}}.$$

(9)

Figure 3 shows the contributions of $W_{\parallel}$ and $W_{tot}$ to electron relativistic energy at the first peak by binning all electrons according to their energies. In the 2D simulation 95% of the energy of energetic electrons ($\varepsilon_{e}>500$ MeV) comes from the work done by the longitudinal electric fields and in the 3D case the quantity is 86%. It is then reasonable to approximate

$$W_{tot}\approx W_{\parallel}=-|e|\int_{-\infty}^{t}E_{\parallel}v_{\parallel}\,dt{{}^{\prime}}.$$

(10)

This simply allows us to narrow down the investigation of electron acceleration to a single component of the electric fields.

At ultra-high laser intensities, the dephasing rate $\delta u$ is approximately $(v_{ph}-c)/c$, since for a relativistic electron co-propagating with laser, $c-v_{x}\ll v_{ph}-c$. Figure 4(a-b) describes the temporal profiles of the transverse electric fields $E_{y}$ recorded in a moving window. Note that $E_{y}$ and $E_{x}$ share the same wave mode and phase velocity. By tracking a fixed field segment, we find that $v_{ph}\approx 1.0031\ c$ in the 2D simulation, and $v_{ph}\approx 1.0063\ c$ in the 3D simulation. The dephasing rate in 2D is nearly two times lower than that in 3D, explaining the distinction between the electron energy gain seen in the two simulations. The duration for electrons staying in an accelerating phase of electric fields can be estimated by

$$\delta t\approx 0.5\lambda_{L}/(\delta u\cdot c),$$

(11)

where $0.5\lambda_{L}$ is the width of the accelerating phase. We find that $\delta t\approx$ 540 fs and 260 fs in 2D and 3D, matching well with the times found from the simulations, demonstrating the accuracy of the approximation used to calculate the phase velocity. The previously derived Eq. 7 however gives the analytical values of the phase velocity as $v_{ph}=1.0078\ c$ in 2D, and $v_{ph}\approx 1.0116\ c$ in 3D. Despite the difference in numerical value, the same trend of an increase in phase velocity is consistent in both the analytical treatment and the numerical simulations. There are a number of factors that can lead to this discrepancy between the analytical calculations and the numerical simulations, for example: the analytical calculation assumes a perfectly conducting boundary and so requires the transverse fields to be zero at the channel edgesGong et al. [2019]; sinusoidal waves are not a perfect match (though close) for the fields observed in the hollow channels. Our additional simulations with different target size $R$ show that as phase velocity varies with $R$, the ratio of $\delta u$ in 2D to $\delta u$ in 3D is preserved.

To further check how the phase velocity influences electron acceleration, we track representative electrons with respect to $E_{x}$ fields sampled by them, as illustrated in Fig. 4(c-d). It is clear that the electrons remain accelerated when they stay in the favourable phase of $E_{x}$ fields (the negative fields colored in blue). After exiting the accelerating phase, the electron energy declines. As can be seen from the time scale in Fig. 4(c-d) electrons in the 2D simulation are subject to far longer periods of acceleration than those in 3D, due to the lower phase velocity in 2D. The evolution of energy and longitudinal work of the chosen electrons is shown in Fig. 4(e). At early moments (up to $\sim$260 fs), the curves of $W_{\parallel_{\ {\rm{2D}}}}$ and $W_{\parallel_{\ {\rm{3D}}}}$ almost overlap with each other implying that the amplitudes of accelerating fields in 2D and 3D are similar. It is then clear that the major cause for the smaller electron energy observed in the 3D simulation is the early termination of the acceleration process due to the higher dephasing rate.

The focus of the discussion in the previous sections was the impact of simulation dimensionality on the laser-driven electron acceleration and generation of energetic electrons with energies of hundreds of MeV to a few GeV. These high-energy electrons are subject to emitting energetic photons (in x-ray and even up to the $\gamma$-ray range) while accelerating in laser and channel fields. In this section, we investigate how the photon emission changes with the dimensionality of simulations. The emitted power $P_{\gamma}$ of synchrotron emission is determined by the electron acceleration in an instantaneous rest frame. This acceleration is proportional to a dimensionless parameterLandau and Lifshitz [1987] $\eta$,

$$\eta\equiv\frac{\gamma_{e}}{E_{S}}\sqrt{\bigg{(}\bf{E}+{\rm{\frac{1}{c}}}[v\times B]\bigg{)}{\rm{{}^{2}}}-{\rm{\frac{1}{c^{2}}}}(E\cdot v){\rm{{}^{2}}}},$$

(12)

where $\bf{E}$ and $\bf{B}$ are the electric and magnetic fields acting on the electron, $\gamma_{e}$ and $\bf{v}$ are the relativistic factor and the velocity of the electron, and $E_{S}\approx 1.3\times 10^{18}\ \rm{V/m}$ is the Schwinger field. The emitted power from an electron scales as $P_{\gamma}\propto\eta^{2}$. We are interested in photons with energy above 100 keV, a threshold shown to be critical for photon-photon pair production Wang et al. [2020].

Collecting all the forward-emitted photons over the duration of the simulations, Fig. 5 compares the spatial distributions of photons produced in the 2D and 3D simulations. Noting that in the simulations performed here, once a photon is generated, the emission location is marked as the photon location for the duration of the simulation. Corresponding to the electron acceleration, the first peak of emission in 3D takes place at $x\sim 60\ \mu$m, well ahead of that in 2D which is located at $x\sim 140\ \mu$m [see Fig. 5(a)]. Up to the first peaks, 38% of laser energy in 2D is transferred to particles (photons, electrons and ions) compared to 41% in 3D. Of the transferred energies 0.47% is converted into photons ($\varepsilon_{\gamma}>100\ \rm{keV}$) in the 2D simulation in contrast to 0.11% in the 3D simulation. As shown in Eq. 12, $\eta$ is directly proportional to electron’s relativistic factor, i.e. $\eta\propto\gamma_{e}$. In the previous section it is noted that the field amplitudes acting on the electrons are similar in both 2D and 3D. Thus the primary cause of the lower conversion rate in 3D is due to the lack of electrons at very high energies. Figure 5(b-c) illustrates the distribution of photons in ($x$,$y$) space. The common feature of both distributions is that most photons are generated close to the channel boundary, corresponding to the surfing motion of energetic electrons depicted in Fig. 2(d-e). The second peak of emission in both simulations becomes weaker in amplitude and accumulates less photon yield following the first peak. The distribution of the emission in the 3D simulation on the transverse plane ($y$,$z$) is given in Fig. 5(d). It is clear that the emission is concentrated around $z=0\ \mu\rm{m}$ plane with two populated lobes formed near the channel boundaries, meaning that the emission pattern on the transverse plane is correlated with the direction of laser polarization.

Figure 6 compares the angular distribution and energy spectrum of 2D and 3D simulations. The photon emission in 3D is projected onto a sphere, as illustrated Fig. 6(a). Though the target is cylindrical, the emission pattern does not preserve the same symmetry. From Fig. 6(b), the energy distribution along the azimuthal angle ${\rm{d}}E_{\gamma}/{\rm{d}}\theta_{\gamma}$ demonstrates a divergence of $10^{\circ}$ (i.e. the FWHM of the energy distribution curve), which is more than two times wider than that in the polar angle direction. The photon beam in the 2D case is found to be more collimated with a divergence of $6^{\circ}$ in ${\rm{d}}E_{\gamma}/{\rm{d}}\theta_{\gamma}$. Distributing emission over $\theta_{\gamma}$ and photon energy, Fig. 6(c-d) manifest the photon beam in 2D is better collimated throughout the whole energy spectrum. To further evaluate the quality of the two photon beams, we compare the beam brilliance. The source size of the 3D target is easily decided by its radius while in 2D the source size is set as 4$\times$2 $\mu$m². It is found that the brilliance of $\gamma$-ray beams is 2.9$\times$$10^{21}$ and 5.8$\times$$10^{20}$ photons/s mm² mrad² 0.1%BW (at 1 MeV) for 2D and 3D respectively. In Fig. 6(e) starting from 1.5 MeV, the number of $\gamma$-ray photons in 2D surpasses that in 3D, projecting a larger brilliance in 2D for energy level above 1 MeV.

In this section, we compared the photon yield with regard to the simulation dimensionality. Since photon emission is a direct result of electron acceleration, the emission demonstrates correlated features as observed for electrons in Section III. The first peak of photon emission in 3D arrives early in space due to the high phase velocity. The lack of high-energy electrons makes the conversion of laser energy to $\gamma$-ray photons less efficient in 3D. Besides the impact on photon number yield, the dimensionality also affects the collimation of the photon beam. With a large divergence angle and a small number of high-energy photons, the 3D $\gamma$-ray beam is less bright than the 2D beam.

We have demonstrated the effects of simulation dimensionality on electron acceleration and $\gamma$-ray production. There are significant distinctions between the results obtained in 2D and 3D setups from both analytical consideration and numerical calculation. Though in numerical values $v_{ph}$ is close to the speed of light, the dephasing rate ($v_{ph}/c-1$) which determines the energy gain varies considerably with simulation geometry. The higher dephasing rate observable in the 3D setup terminates electron acceleration process early in space and time and leads to a reduction of photon emission when compared to 2D.

A 2D channel target presents a planar geometry while a 3D target has a cylindrical symmetry. Analytically we have shown that the phase velocity of laser fields propagating inside the targets closely depends on the target geometry; the phase velocity in 2D is smaller than that in 3D. To be more specific, the dephasing rate in a 2D setup is derived to be 1.5 times slower when controlling for the channel radius and laser wavelength. Through numerical simulations, it is found that the absolute majority of work done on energetic electrons comes from the longitudinal electric field, enabling the investigation of laser-driven particle acceleration based only on one single component of the electric fields. By tracking a fixed segment of laser fields the phase velocity in 2D is again shown to be smaller than in 3D, matching the correct trend shown in the analytical derivation. It is clear that in 2D simulations electrons surf for a longer period in an accelerating phase compared to 3D simulations. As a result, electrons in 2D have an elongated acceleration distance and present a more energetic spectrum, leading to an overestimate of maximum electron energy and the number of high-energy electrons when compared to more realistic 3D simulations.

Similar to electron acceleration, photon emission is strongly impacted by simulation dimensionality. Due to a lack of high-energy electrons, in 3D the photon ($\varepsilon_{\gamma}>100$ keV) conversion rate is more than 4 times smaller and the emission falls behind on the generation of energetic photons ($\varepsilon_{\gamma}>1.5$ MeV). In 3D, the emission is accumulated close to the $z=0$ plane with a narrow divergence along the polar angle. However the emission along the azimuthal angle is far more diverged in 3D compared to that in a 2D simulation. As a result, the photon beam produced in a 3D setup is found to be less bright. It is also worth noting that the subsequent peaks of photon emission drop sharply in amplitude and duration.

Though 2D numerical simulations are widely applied in the research of laser and micro-channel interactions, one should not ignore the overestimate and inaccuracy of results caused by the low phase velocity in 2D simulations. In particular, when carrying out numerical simulations to predict and optimize the output for laser micro-channel experiments, it matters to accurately know the exact location of peaked electron spectrum and therefore 3D simulations are indispensable. We conclude that the 2D simulations are capable of qualitatively reproducing the features of 3D simulations, but for quantitative evaluations and reliable predictions, 3D modelling is strongly recommended.

T.W. and A.A. were supported by AFOSR (Grant No. FA9550-17-1-0382). D.B. was supported by NSF (Grant No. PHY 1903098). K.C. was supported by NSF (Grant No. 1821944). Simulations were performed with EPOCH (developed under UK EPSRC Grants No. EPG0549401, No. EPG0551651, and No. EPG0568031) using HPC resources provided by TACC at the University of Texas. This work used XSEDE, supported by NSF grant number ACI-1548562.