Numerical modeling of laser tunneling ionization in ... - Semantic Scholar
Journal of Computational Physics 236 (2013) 220–228
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Numerical modeling of laser tunneling ionization in explicit particle-in-cell codes M. Chen a,⇑, E. Cormier-Michel b, C.G.R. Geddes a, D.L. Bruhwiler b, L.L. Yu a, E. Esarey a, C.B. Schroeder a, W.P. Leemans a a b
Lawrence Berkeley National Laboratory, 1 Cyclotron Rd., Berkeley, CA 94720, USA Tech-X Corporation, 5621 Arapahoe Ave., Suite A, Boulder, CO 80303, USA
a r t i c l e
i n f o
Article history: Received 28 March 2012 Received in revised form 16 September 2012 Accepted 26 November 2012 Available online 7 December 2012 Keywords: Particle in cell Tunneling ionization Laser wakefield acceleration Laser plasma Ionization injection
a b s t r a c t Methods for the calculation of laser tunneling ionization in explicit particle-in-cell codes used for modeling laser–plasma interactions are compared and validated against theoretical predictions. Improved accuracy is obtained by using the direct current form for the ionization rate. Multi level ionization in a single time step and energy conservation have been considered during the ionization process. The effects of grid resolution and number of macro-particles per cell are examined. Implementation of the ionization algorithm in two different particle-in-cell codes is compared for the case of ionization-based electron injection in a laser–plasma accelerator. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction In numerical studies of ultra-intense laser–plasma interactions, the plasma production via the laser ionization is often omitted due to the fact that the ionization process occurs before the main body of the pulse arrives. For example, a laser– plasma accelerator [1], driven by a short-pulse (tens of fs) laser, typically uses low-Z gas (e.g., H or He) targets that are easily ionized by the ps laser pre-pulse present in the laser temporal profile. For relativistically intense (I J 1018 W/cm3) lasers interacting with low-Z gases, ionization effects on plasma wave generation are negligible [2], however pre-pulse effects may still strongly influence the main pulse propagation [3,4]. For laser intensities close to the tunneling ionization threshold, or when a high-Z gas target is used, ionization processes can significantly influence the laser–plasma interaction. Ionization effects in laser–plasma interactions, including laser frequency blue shifting, harmonic generation, ionizationinduced diffraction, and instabilities, have been studied theoretically and experimentally for laser pulses interacting with neutral targets [5–16]. For high-field laser–gas interactions the ionization process is dominated by tunneling ionization. In laser–solid interactions, due to the overdense target, the laser pulse cannot penetrate into the target, and ionization process inside the target is mainly due to impact ionization [17]. Although the impact ionization process is important to correctly describe the laser–solid interaction process, in this work we will focus on the modeling of laser tunneling ionization. Accurate modeling of the ionization process is critical to several recently proposed methods of radiation and particle beam generation. For example, ionization is the key process in ionization-current-based THz radiation generation in laser–gas interactions [18–20] and ionization-induced trapping in laser–plasma accelerators [21–27]. Laser-ionization-based trapping ⇑ Corresponding author. Tel.: +1 510 486 6632; fax: +1 510 486 7981. E-mail address: [email protected] (M. Chen). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.11.029
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
221
in plasma accelerators relies on the use of a high-Z gas target to generate electrons at the appropriate phase of the plasma wakefield. For this application, an accurate description of the ionization process is critical for predicting the laser–plasma accelerated beam properties. Particle-in-cell (PIC) [28,29] is a common numerical method used to study short-pulse laser–plasma interactions. The PIC numerical method models the Vlasov–Maxwell equations using a particle-grid method, such that the continuous phasespace fluid is sampled by a finite number of (Lagrangian) macro-particles and the Maxwell equations are discretized on a grid. Laser ionization has been implemented and studied in PIC codes, and several types of ionization formulae have been used in the literature (see, for example, Refs. [30,17]). In this paper, we focus on the numerical implementation in a PIC code of optical tunneling ionization based on the Ammosov–Delone–Krainov (ADK) model [31]. We compare the alternating current (AC) and direct current (DC) models and show that the DC model must be used to give accurate results in explicit PIC simulations. Multi level ionization in a single simulation step and energy conservation during ionization are considered. The effects of numerical parameters on the modeling of the laser ionization process are presented. We consider the effects of the resolution of the grid and the average number of numerical macro-particles per cell representing the phase-space distribution. A comparison of ionization injection studies using two PIC frameworks (VLPL [32] and VORPAL [33]) is also made.
2. Tunneling ionization formulae Tunneling ionization from the ground state of a Hydrogen atom in an electrostatic (DC) field was first solved schematically by Landau [34] for the dependence of the ionization rate on the applied external field. When the intensity of the external electric field satisfies Edc Ea (where eEa =me c2 ¼ a4 =re or Ea ’ 5:1 GV=cm), the ionization rate for a Hydrogen-like system is given by the Landau–Lifshitz DC tunneling ionization formula [34]:
W LL ¼ 4xa
" 5=2 3=2 # Ea U ion 2 Ea U ion ; exp 3 Edc U H Edc U H
ð1Þ
where xa ¼ a3 c=re ¼ 4:13 1016 s1 is the atomic unit frequency, U ion is the ionization potential, and U H ¼ 13:6 eV is the ionization potential of Hydrogen at the fundamental state. Here me is the electron mass, c is the speed of light, a ¼ 1=137 is the fine structure constant, and r e ¼ e2 =me c2 is the classical electron radius. For complex atoms in a static field (Edc ) the DC tunneling ionization formula is [31]
W DClm ¼ x
2 a C n l
ð2l þ 1Þðl þ jmjÞ!
U ion 2U H 2jmj ðjmjÞ!ðl jmjÞ!
" " 3=2 #2n jmj1 3=2 # Ea U ion 2 Ea U ion ; 2 exp 3 Edc U H Edc U H
ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where l; m are the electron’s orbital quantum number and its projection, respectively, n ¼ Z U H =U ion is the effective prin cipal quantum number, l ¼ n0 1 is the effective value of the orbital number (n0 is the effective principal quantum number of the ground state.), and Z is the ion charge number after ionization. The coefficients C n l can be found calculated using a semi-classical approximation (see Eq. (20) of Ref. [31]):
C 2n l ¼
n l þ1=2 22n 1 4e2 n l : ’ n Cðn þ l þ 1ÞCðn l Þ 2pn n2 l2 n þ l
ð3Þ
In the limit l n , Eq. (3) reduces to
C 2n l ¼
2n 1 2e : 2pn n
ð4Þ
The ionization rate of a hydrogen atom in an alternating (AC) electric field, e.g., a laser field [EL / cosðxtÞ, where x is the laser frequency], can be treated, in the appropriate limit, as the DC ionization rate averaged over a single period of the electric field oscillation. This model was proposed by Keldysh [35] and then later expanded to hydrogen-like atoms by Ammosov et al. [31] (referred to as the ADK model). The DC ionization rate is related to the AC ionization rate by
"
W AC
3=2 #1=2 3 E UH ¼ W DC : p Ea U ion
ð5Þ
Combining Eqs. (2), (4), and (5), yields the one cycle averaged ionization rate for a complex atom in an AC field. For a linearly polarized laser pulse EL , the ionization rate is [31]
W AClm ¼ xa
" !2n jmj1 3 1=2 2 2n 3 # 3n EL Z 2e 2 Ea Z ð2l þ 1Þðl þ jmjÞ! Ea Z 3 2 exp : 3 EL n EL n3 2n2 n pZ 3 Ea 2pn 2jmj ðjmjÞ!ðl jmjÞ!
ð6Þ
222
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
For the fundamental state l ¼ m ¼ 0, we obtain the ionization rate
pffiffiffi e 3=2 Z 2 Ea Z 3 W AC00 ¼ xa 3 4e 9=2 p EL n4 n
!2n 3=2
" exp
3 # 2 Ea Z : 3 EL n
ð7Þ
In the original ADK paper a typographic error gave an ionization rate (see Eq. (13a) of Ref. [31]) of p=e times the correct rate. This and other corrections, including errors in Eq. (20) of Ref. [31], have been pointed out by several authors [36,37]. However, the original formula with the typographic error has been used in several studies (e.g., Refs. [38,17]), and is pointed out here to prevent future confusion. In numerical modeling of laser–plasma interactions via PIC it is essential to use the correct ionization model. If the laser pulse is modeled as an envelope and the time step is much larger than the laser period, an ionization rate based on a cycleaveraged AC model should be used. However, in an explicit PIC simulation, the sub-wavelength scale is well-resolved and the simulation time step is much smaller than the laser period (dt T 0 , where dt is the simulation time step and T 0 ¼ c=k0 is the laser period), and within each simulation time step, the laser field can be viewed as a DC field. In this case the DC ionization rate model should be employed. To illustrate the differences in the ionization models, Fig. 1 shows numerical calculations comparing the DC model Eq. (2) and the AC model Eq. (6) in a calculation that resolves the laser period. Fig. 1(a) shows the ionization rate from N5þ to N6þ versus normalized electric field intensity a ¼ eEL =me xc. As shown in Fig. 1(a), the AC ionization model (as used in Ref. [17]) shows the same trend as the DC model; however, it always provides a lower ionization rate, as anticipated from Eq. (5). The ionization probability may also be considered within a finite time. Fig. 1(b) shows the ionization probability (from N5+ to N6+) after a single laser cycle (a ¼j a0 sinð2ptÞ j, with 0 6 t 6 1) versus laser intensity: The red curve shows the result by using the theoretical calculation of 1 exp½Wðt p ÞT 0 , where tp represents the time at the electric field peak. The black curve (overlapped by the red curve) shows the result from the DC integration model, and the blue curve shows the result from the AC integration model. As shown in Fig. 1(b), use of the AC model in this case has underestimated the ionization probability. Fig. 1(c) shows the difference of the two models for a temporally Gaussian pulse, with full-width-half-maximum (FWHM) length LFWHM ¼ 14:89cT 0 . Again, the AC model gives a lower ionization probability than the DC model. From this, we expect that, while results of studies using explicit codes with the AC ionization model (or using the uncorrected ADK coefficient) may give qualitatively reasonable results, the use of the corrected DC model in explicit simulations is essential to obtain accurate quantitative results.
(a)
(b)
(c)
Fig. 1. (a) Ionization rate from N5+ to N6+ versus electric field intensity a ¼ eEL =me cx. The laser frequency is x ¼ 2pc=k0 with k0 ¼ 0:8 lm. The solid black solid curve is based on the DC model and the dashed blue curve is based on the AC averaged model. (b) Ionization probability after one laser cycle calculated using different models: solid black curve integrates over W DClm , the dashed red curve uses the one-cycle-averaged formula 1 exp½Wðtp ÞT 0 , the blue dashed curve integrates over W AClm . (c) Ionization probability after a Gaussian pulse: solid black curve is integration over W DClm integration, blue dashed curve integrates over W AClm . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
223
3. Tunneling ionization model implementation in PIC codes In an explicit PIC simulation, the ionization process is generally implemented with the following prescription. At each time step, for each un-ionized particle, the local electric field and particle ionization parameters (such as ionization potential and charge state) are used to calculate the ionization rate WðtÞ. The particle ionization probability is then calculated by PðtÞ ¼ 1 exp½WðtÞdt, which is approximated as PðtÞ ’ WðtÞdt when PðtÞ 1. The code generates a random number p with a uniform distribution between 0 and 1; once p < PðtÞ, the particle will be set to be ionized, otherwise the particle will not be ionized and will be evaluated during the next time step. In the PIC code the particles are macro-particles representing the phase space fluid. The ionization probability cannot be accumulated to the next time step if the particle is not ionized; otherwise an over-counted ionization probability results. Although for most of our simulation cases, the time step of PIC simulations can satisfy PðtÞ 1, however, if the ionization potential is too low or a relatively high laser intensity is used, the above approximation can not be satisfied and high level ionization in a single PIC simulation time step may happen. A subroutine will be used to solve this issue. Let us consider an ion particle whose initial state is Zi and its highest charge is z. The evolution of the state population (N j ) is as following: i
dN =dt ¼ W i Ni dN
iþ1
ð8Þ
=dt ¼ W i Ni W iþ1 Niþ1
ð9Þ ð10Þ
... dN
z1
=dt ¼ W z2 Nz2 W z1 Nz1
ð11Þ
z
dN =dt ¼ W z1 Nz1 i
ð12Þ j
and the initial state is: N ðt ¼ 0Þ ¼ 1 and N ðt ¼ 0Þ ¼ 0 when j > i. Obviously during the whole ionization process (t 2 ½0; dt), P j the total number of all kinds of ions is constant: j¼z j¼i N ¼ 1. So only z i equations need to be solved. Depends on the specific dependance of W j on t different approximations can be used to simplify these equations. Generally a subroutine can be installed in the PIC code to solve this equation array by using the typical fourth order Runge–Kutta method. Once one knows the state population (N j ) at the end of time step t ¼ dt, a random number p with a uniform distribution between 0 and 1 will P P j k be generated to decide the final ionization state. When j¼k1 N j < p < j¼k j¼i j¼i N , the particle will be set to be ionized to Z and k i electrons will be generated. Here we neglect the different response of the different valence electrons in such a single time step, which is reasonable when the time step dt of the PIC simulation is small enough. Besides the high level ionization process, ionization induced energy loss (Eloss ion ) is another issue to be considered. Although in some simulation cases such as laser short–gas interaction, the energy loss due to ionization can be neglected compared with the high laser energy. However, in some special cases such as long distance laser propagation in neutral gas or laser high density high Z neutral matter interaction, energy loss due to ionization should be considered. Currently, to retain the energy conservation, an approach using artificial ionization current (~ jion ) is suggested. Such method has been discussed by Kemp et al. [17] before. The artificial ionization current is directed along the local electric field ~ E, and its magnitude is set ~2 to satisfy the energy conservation law: ~ jion ¼ ~ ERk Dnek ðEloss ionk þ eki Þ=dtjEj . Here k represents the different ionization state, Dnek th represents the number density of the k ionized electrons and eki represents their initial energies just after ionization which is usually set to be zero. This current is then added to the physical current in each cell, after charge and current are computed from particle positions and velocities and before the integration of Maxwell’s equations. Such process should be done for every time step and each cell once the ionization process happens. The above analysis also indicates that, to accurately model the ionization process, random number generation, selection of the macro-particle number representing the ionizing atoms, and the choice of time step size are all critical numerical parameters. In the following we examine the effect of these numerical parameters on the ionization modeling using the VLPL PIC code [32]. In the VLPL code, the ionization rate formula of the DC ionization model Eq. (2) is implemented. A library source random number generator is used to produce random numbers uniformly distributed in the region ½0; 1. To study the effects of time step and macro-particles per cell on the ionization probability calculation, we consider the ionization of a Nitrogen gas by a linearly polarized laser pulse. To exclude other physical effects, the number density of the initial neutral Nitrogen atoms is 105 nc , where nc ¼ pme c2 =e2 k20 is the critical density (nc ’ 1:7 1021 cm3 for a laser of wavelength of k0 ¼ 0:8 lm) and the simulation is in a one-dimensional (1D) geometry. The normalized laser intensity is a0 ¼ eEL =me xc ¼ 1:8 and the pulse length is L0 ¼ 10cT 0 , with the temporal shape of electric field defined as / expðt2 =L20 Þ. The neutral Nitrogen atoms are initially uniformly distributed within the region of 5k0 —35k0 . No other particle elements are used. For such a laser pulse the highest ionization rate for Nitrogen N5+ is 4:68 1014 s1 at the peak of the laser pulse, as shown in Fig. 1(a). Fig. 2(a) and (b) shows the distribution of ionization probability of N5+ (fraction of ions in charge states N6+ or N7+) along the laser propagation direction. The laser propagates from the left to the right and the grey curve shows the laser field (with peak at x ¼ 15k0 ). The red dashed curve shows the theoretical ionization calculated using Eq. (2) and assuming N5+ is uniformly pre-ionized before the laser arrives. The assumption is reasonable since Nitrogen is easily ionized to the fifth valence
224
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
δ
(a)
δ δ
(b)
Fig. 2. Ionization probability distribution of N5+ along the laser propagation direction x=k0 . (a) The grey curve is the normalized laser electric field with peak field at x ¼ 15k0 . The red dashed curve shows the theoretical result, the green dotted curve is calculate from the PIC code with dt ¼ 0:1T 0 and 40 macroparticles per cell. (b) The black solid curve uses dt ¼ 0:01T 0 and 4 macro-particles per cell, the light blue solid curve used dt ¼ 0:005T 0 and 1 macro-particle per cell. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 3. Example of the typical laser–plasma accelerator wake structure and ionization injection scheme: laser electric field (grey curve) along its propagation direction (x), the electric field of the wake (red curve), the momentum of the background plasma electron fluid (dashed blue curve), the trapping momentum threshold of the electrons in the wakefield (dotted green curve), and the momentum gap between the trapping threshold and the background electron momentum (black curve). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
state at the front of the laser pulse. The green dotted curve in Fig. 2 shows the sum density of the N6+ and N7+ normalized by the initial Nitrogen density N 0 , which represents the total ionization probability of N5+, calculated from the PIC code with a time resolution of dt ¼ 0:1T 0 ’ 0:267 fs and 40 macro-particles per cell. The ionization probability has been underestimated in the PIC code by 21.4%. At the peak of the laser pulse, Wdt ’ 0:125 and the approximation of 1 expðWdtÞ Wdt gives a relative error of +6.38%, i.e., the error in the ionization calculation does not come from the small exponent approximation, rather it is largely due to the sampling of the electric fields.
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
225
Since the ionization mainly happens at the peak of each cycle, if the resolution is too course, the sampling of fields (i.e., the resolution) may not be sufficiently accurate. To study the effect of resolution, we considered increased temporal resolution with the total macro-particle number within a fixed length constant. Increased temporal resolution is equivalent to increased longitudinal spatial resolution, since in the 1D PIC code dx=k0 ¼ dt=T 0 . The black solid curve in Fig. 2(b) shows the ionization probability calculated from the PIC code with dt ¼ 0:01T 0 and 4 macro-particles per cell. Better agreement with the theoretical curve (dashed red curve) is achieved with the increased resolution. The increase in noise is due to the fewer macro-particles per cell (i.e., larger weight given to each macro-particle). For example, if the resolution is increased to dt ¼ 0:005T 0 with 1 macro-particle per cell, increased simulation noise is observed (light blue solid curve in Fig. 2(b)). The ionization probability calculation is sensitive to the resolution, and accurate modeling of the ionization process the time step of each should be sufficiently small. In general, the resolution criterion is more strict than that due to the ionization rate’s sensitivity to the electric field strength WðtÞdt 1. The number of macro-particle per cell affects the simulation noise level. The larger the number per cell, the smaller weight of each macro-particle, and the simulation noise will be less. When modeling ionization with PIC, evaluating the convergence of the resolution is essential for an accurate description of ionization process.
(a)
(b)
(c)
Fig. 4. (a) Evolution of trapped (with momentum larger than 20 MeV/c) electron number versus laser propagation distance: VLPL using DC model (solid black squares), VORPAL using DC model (solid red circles), VLPL using AC model (open black squares), and VLPL using DC model and 4 macro-particles per cell (blue solid triangles). (b) Evolution of beam mean energy (open points) and rms energy spread (solid points) versus laser propagation distance using VLPL (black square) and VORPAL (red circle). (c) Beam transverse momentum spread evolution versus laser propagation distance using VLPL (black square) and VORPAL (red circle). In these simulations, 1 macro-particle per cell for the pre-ionized electrons and 6 macro-particles per cell for Nitrogen atoms are used. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
226
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
4. Comparison of ionization injection studies The implementation of the DC ionization model in PIC codes described above has been carried out in two PIC codes VLPL [32] and VORPAL [33]. In this section we compare results obtained in both PIC codes using the same DC ionization model. The case studied is for ionization injection into a laser–plasma accelerator, and the final simulated beam characteristics are compared. High-Z gas ionization injection has recently been employed experimentally as an effective electron injection method in the laser–plasma accelerator [22–26]. The detailed physics of ionization injection has been describe in Ref. [27]. Fig. 3 shows a typical wake structure excited by an intense laser and the injection criterion based on 1D nonlinear cold plasma wake theory at a fixed time. A laser pulse (grey curve) is used to generate a wake in a mixed plasma of pre-ionized electrons and neutral Nitrogen atoms. The red solid curve shows the longitudinal electric field of the wake when only the initial pre-ionized electron plasma is included. Since the concentration of Nitrogen is small (typically 20 MeV/c) as calculated in the two PIC codes using the DC ionization model. Using the AC model gives a lower injection number, which is approximately half that of the DC model result for these parameters. The calculated injected beam mean energy and root-mean-square (r.m.s.) energy spread are also similar for both codes when the DC ionization model is used (shown in Fig. 4(b)). For ionization-induced injection, the beam transverse momentum spread is mainly due to the residual transverse momentum as a result of the ionization process. In the 1D case, canonical momentum conservation ensures ~ p? ¼ ~ aðt i Þ, where t i is the ionization time. Therefore, the electrons have a transverse momentum equal to the laser vector potential at ionization. After the injection, the transverse momentum evolution depends on the transverse wakefield structure. Fig. 4(c) shows the FWHM transverse momentum distribution as calculated by the two codes, with differences less than 10%. This level of agreement is reasonable, considering the numerical noise generated by the ionization (i.e., the low macro-particles per cell). Increasing the number of macro-particles per cell for the background plasma electrons (from 1 to 4 macro-particles per cell), improves the agreement, indicating the source of errors is numerical noise in the PIC code. It also shows the electrons’ transverse momentum spread not only depends on the correct simulation of ionization process, it depends on the correct wake calculation too. Typical beam momentum distributions are shown in Fig. 5 (where px ; py are the longitudinal and transverse beam momenta, respectively). Differences in the beam transverse momenta distribution are the result of differences in the transverse wakefield structure calculated by the two codes. The above comparison provides a benchmark of the ionization algorithms implemented in VLPL and VORPAL for parameters relevant to studies of ionization injection in laser–plasma accelerators. 5. Summary and conclusions Implementation of laser tunneling ionization rate formulae in explicit PIC codes has been examined. The originally published ionization rate formulae were reviewed, and corrections noted. It was shown, by comparison to theoretical predictions, that the DC form of the ionization rate provides improved accuracy. A subroutine is used to deal with the high level ionization in a single time step and an ionization current is used to consider energy loss due to ionization. The effect of PIC numerical parameters, such as grid resolution and number of macro-particles per cell, on the calculated ionization probability was studied. It was found that the ionization probability is sensitive to the grid resolution, and increased number of macro-particles per cell can reduce the simulation noise. Modeling the ionization injection in a laser–plasma accelerator was used to benchmark the ionization algorithm implemented in two PIC codes. Acknowledgments The authors thank Carlo Benedetti and Stepan Bulanov for useful discussions. This work was supported by the Director, Office of Science, Office of High Energy Physics, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and by a Scientific Discovery through Advanced Computing Project ‘‘Advanced Computing for 21st Century Accelerator Science and Technology’’. Computational resources of the National Energy Research Scientific Computing Center were used to perform the simulations. References [1] E. Esarey, C.B. Schroeder, W.P. Leemans, Physics of laser-driven plasma-based electron accelerators, Rev. Mod. Phys. 81 (2009) 1229–1285. [2] C.G.R. Geddes, C. Tóth, J. van Tilborg, E. Esarey, C.B. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, W.P. Leemans, Production of high-quality electron bunches by dephasing and beam loading in channeled and unchanneled laser plasma accelerators, Phys. Plasmas 12 (2005) 056709. [3] T. Hosokai, K. Kinoshita, A. Zhidkov, K. Nakamura, T. Watanabe, T. Ueda, H. Kotaki, M. Kando, K. Nakajima, M. Uesaka, Effect of a laser prepulse on a narrow-cone ejection of mev electrons from a gas jet irradiated by an ultrashort laser pulse, Phys. Rev. E 67 (3) (2003) 036407. [4] T. Hosokai, K. Kinoshita, T. Ohkubo, A. Maekawa, M. Uesaka, A. Zhidkov, A. Yamazaki, H. Kotaki, M. Kando, K. Nakajima, S.V. Bulanov, P. Tomassini, A. Giulietti, D. Giulietti, Observation of strong correlation between quasimonoenergetic electron beam generation by laser wakefield and laser guiding inside a preplasma cavity, Phys. Rev. E 73 (3) (2006) 036407. [5] E. Esarey, A. Ting, P. Sprangle, Frequency shifts induced in laser pulses by plasma waves, Phys. Rev. A 42 (1990) 3526. [6] E. Esarey, G. Joyce, P. Sprangle, Frequency upshifting of laser pulses by copropagating ionization fronts, Phys. Rev. A 44 (1991) 3908. [7] W.P. Leemans, C.E. Clayton, W.B. Mori, K.A. Marsh, P.K. Kaw, A. Dyson, C. Joshi, J.M. Wallace, Experiments and simulations of tunnel-ionized plasmas, Phys. Rev. A 46 (1992) 1091–1105. [8] W.P. Leemans, C.E. Clayton, W.B. Mori, K.A. Marsh, A. Dyson, C. Joshi, Plasma physics aspects of tunnel-ionized gases, Phys. Rev. Lett. 68 (1992) 321– 324. [9] R.L. Savage, C. Joshi, W.B. Mori, Frequency upconversion of electromagnetic radiation upon transmission into an ionization front, Phys. Rev. Lett. 68 (1992) 946–949, http://dx.doi.org/10.1103/PhysRevLett.68.946. URL: . [10] B. Cros, X. Xu, T. Tsukada, N. Yugami, Y. Nishida, G. Matthieussent, Ionization-induced frequency up-shift of a high-power microwave interacting with a plasma, Phys. Plasmas 4 (8) (1997) 2837–2844, http://dx.doi.org/10.1063/1.872416. URL: .
228
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
[11] P. Sprangle, E. Esarey, J. Krall, Self-guiding and stability of intense optical beams in gases undergoing ionization, Phys. Rev. E 54 (4, Part b) (1996) 4211– 4232, http://dx.doi.org/10.1103/PhysRevE.54.4211. [12] E. Esarey, P. Sprangle, J. Krall, A. Ting, Self-focusing and guiding of short laser pulses in ionizing gases and plasmas, IEEE J. Quantum Electron. 33 (11) (1997) 1879–1914, http://dx.doi.org/10.1109/3.641305. [13] P. Sprangle, E. Esarey, B. Hafizi, Propagation and stability of intense laser pulses in partially stripped plasmas, Phys. Rev. E 56 (5, Part b) (1997) 5894– 5907, http://dx.doi.org/10.1103/PhysRevE.56.5894. [14] P. Sprangle, E. Esarey, B. Hafizi, Intense laser pulse propagation and stability in partially stripped plasmas, Phys. Rev. Lett. 79 (6) (1997) 1046–1049, http://dx.doi.org/10.1103/PhysRevLett.79.1046. [15] T. Antonsen, Z. Bian, Ionization induced scattering of short intense laser pulses, Phys. Rev. Lett. 82 (18) (1999) 3617–3620, http://dx.doi.org/10.1103/ PhysRevLett.82.3617. [16] Z. Bian, T. Antonsen, Ionization instabilities of an electromagnetic wave propagating in a tenuous gas, Phys. Plasmas 8 (7) (2001) 3183–3194, http:// dx.doi.org/10.1063/1.1377612. [17] A.J. Kemp, R.E.W. Pfund, J.M. ter Vehn, Modeling ultrafast laser-driven ionization dynamics with monte carlo collisional particle-in-cell simulations, Phys. Plasmas 11 (12) (2004) 5648–5657. [18] M. Kre, T. Löffler, M.D. Thomson, R. Dörner, H. Gimpel, K. Zrost, T. Ergler, R. Moshammer, U. Morgner, J. Ullrich, H.G. Roskos, Determination of the carrier-envelope phase of few-cycle laser pulses with terahertz-emission spectroscopy, Nat. Phys. 2 (2006) 327–331. [19] K.-Y. Kim, J.H. Glownia, A.J. Taylor, G. Rodriguez, Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields, Opt. Express 15 (8) (2007) 4577–4584. [20] M. Chen, A. Pukhov, X.-Y. Peng, O. Willi, Theoretical analysis and simulations of strong terahertz radiation from the interaction of ultrashort laser pulses with gases, Phys. Rev. E 78 (2008) 046406. [21] M. Chen, Z.-M. Sheng, Y.-Y. Ma, J. Zhang, Electron injection and trapping in a laser wakefield by field ionization to high-charge states of gases, J. Appl. Phys. 99 (5) (2006) 056109. [22] A. Pak, K.A. Marsh, S.F. Martins, W. Lu, W.B. Mori, C. Joshi, Injection and trapping of tunnel-ionized electrons into laser-produced wakes, Phys. Rev. Lett. 104 (2) (2010) 025003. [23] C. McGuffey, A.G.R. Thomas, W. Schumaker, T. Matsuoka, V. Chvykov, F.J. Dollar, G. Kalintchenko, V. Yanovsky, A. Maksimchuk, K. Krushelnick, V.Y. Bychenkov, I.V. Glazyrin, A.V. Karpeev, Ionization induced trapping in a laser wakefield accelerator, Phys. Rev. Lett. 104 (2) (2010) 025004. [24] C.E. Clayton, J.E. Ralph, F. Albert, R.A. Fonseca, S.H. Glenzer, C. Joshi, W. Lu, K.A. Marsh, S.F. Martins, W.B. Mori, A. Pak, F.S. Tsung, B.B. Pollock, J.S. Ross, L.O. Silva, D.H. Froula, Self-guided laser wakefield acceleration beyond 1 GeV using ionization-induced injection, Phys. Rev. Lett. 105 (2010) 105003. [25] B.B. Pollock, C.E. Clayton, J.E. Ralph, F. Albert, A. Davidson, L. Divol, C. Filip, S.H. Glenzer, K. Herpoldt, W. Lu, K.A. Marsh, J. Meinecke, W.B. Mori, A. Pak, T.C. Rensink, J.S. Ross, J. Shaw, G.R. Tynan, C. Joshi, D.H. Froula, Demonstration of a narrow energy spread, 0:5GeV electron beam from a two-stage laser wakefield accelerator, Phys. Rev. Lett. 107 (2011) 045001, http://dx.doi.org/10.1103/PhysRevLett.107.045001. URL: . [26] J.S. Liu, C.Q. Xia, W.T. Wang, H.Y. Lu, C. Wang, A.H. Deng, W.T. Li, H. Zhang, X.Y. Liang, Y.X. Leng, X.M. Lu, C. Wang, J.Z. Wang, K. Nakajima, R.X. Li, Z.Z. Xu, All-optical cascaded laser wakefield accelerator using ionization-induced injection, Phys. Rev. Lett. 107 (2011) 035001, http://dx.doi.org/10.1103/ PhysRevLett.107.035001. URL: . [27] M. Chen, E. Esarey, C.B. Schroeder, C.G.R. Geddes, W.P. Leemans, Theory of ionization-induced trapping in laser–plasma accelerators, Phys. Plasmas 19 (3) (2012) 033101. [28] C.K. Birdsall, A.B. Langdon, V. Vehedi, J.P. Verboncoeur, Plasma Physics via Computer Simulations, Adam Hilger, Bristol, 1991. [29] R. Hockney, J. Eastwood, Computer Simulation using Particles, Taylor & Francis Group, New York, 1988. [30] D.L. Bruhwiler, D.A. Dimitrov, J.R. Cary, E. Esarey, W. Leemans, R.E. Giacone, Particle-in-cell simulations of tunneling ionization effects in plasma-based accelerators, Phys. Plasmas 10 (5) (2003) 2022–2030. [31] M.V. Ammosov, N.B. Delone, V.P. Krainov, Sov. Phys. JETP 64 (1986) 1191. [32] A. Pukhov, Three-dimensional electromagnetic relativistic particle-in-cell code VLPL (virtual laser plasma lab), J. Plasma Phys. 61 (1999) 425–433. [33] C. Nieter, J. Cary, VORPAL: a versatile plasma simulation code, J. Comput. Phys. 196 (2004) 448. [34] L. Landau, E. Lifshitz, Quantum Mechanics, Butterworth-Heinemann, New York, 1965. [35] L.V. Keldysh, J. Exp. Theor. Phys. (USSR) 47 (1964) 1945. [36] X. Lu et al, Acta Opt. Sin. 8 (1999) 1014. [37] F.A. Ilkov, J.E. Decker, S.L. Chin, Ionization of atoms in the tunnelling regime with experimental evidence using hg atoms, J. Phys. B At. Mol. Opt. Phys. 25 (19) (1992) 4005–4020. [38] B.M. Penetrante, J.N. Bardsley, Residual energy in plasmas produced by intense subpicosecond lasers, Phys. Rev. A 43 (1991) 3100–3113. [39] C.B. Schroeder, E. Esarey, B.A. Shadwick, W.P. Leemans, Trapping, dark current, and wavebreaking in nonlinear plasma waves, Phys. Plasmas 13 (3) (2006) 033103.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Numerical modeling of laser tunneling ionization in explicit particle-in-cell codes M. Chen a,⇑, E. Cormier-Michel b, C.G.R. Geddes a, D.L. Bruhwiler b, L.L. Yu a, E. Esarey a, C.B. Schroeder a, W.P. Leemans a a b
Lawrence Berkeley National Laboratory, 1 Cyclotron Rd., Berkeley, CA 94720, USA Tech-X Corporation, 5621 Arapahoe Ave., Suite A, Boulder, CO 80303, USA
a r t i c l e
i n f o
Article history: Received 28 March 2012 Received in revised form 16 September 2012 Accepted 26 November 2012 Available online 7 December 2012 Keywords: Particle in cell Tunneling ionization Laser wakefield acceleration Laser plasma Ionization injection
a b s t r a c t Methods for the calculation of laser tunneling ionization in explicit particle-in-cell codes used for modeling laser–plasma interactions are compared and validated against theoretical predictions. Improved accuracy is obtained by using the direct current form for the ionization rate. Multi level ionization in a single time step and energy conservation have been considered during the ionization process. The effects of grid resolution and number of macro-particles per cell are examined. Implementation of the ionization algorithm in two different particle-in-cell codes is compared for the case of ionization-based electron injection in a laser–plasma accelerator. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction In numerical studies of ultra-intense laser–plasma interactions, the plasma production via the laser ionization is often omitted due to the fact that the ionization process occurs before the main body of the pulse arrives. For example, a laser– plasma accelerator [1], driven by a short-pulse (tens of fs) laser, typically uses low-Z gas (e.g., H or He) targets that are easily ionized by the ps laser pre-pulse present in the laser temporal profile. For relativistically intense (I J 1018 W/cm3) lasers interacting with low-Z gases, ionization effects on plasma wave generation are negligible [2], however pre-pulse effects may still strongly influence the main pulse propagation [3,4]. For laser intensities close to the tunneling ionization threshold, or when a high-Z gas target is used, ionization processes can significantly influence the laser–plasma interaction. Ionization effects in laser–plasma interactions, including laser frequency blue shifting, harmonic generation, ionizationinduced diffraction, and instabilities, have been studied theoretically and experimentally for laser pulses interacting with neutral targets [5–16]. For high-field laser–gas interactions the ionization process is dominated by tunneling ionization. In laser–solid interactions, due to the overdense target, the laser pulse cannot penetrate into the target, and ionization process inside the target is mainly due to impact ionization [17]. Although the impact ionization process is important to correctly describe the laser–solid interaction process, in this work we will focus on the modeling of laser tunneling ionization. Accurate modeling of the ionization process is critical to several recently proposed methods of radiation and particle beam generation. For example, ionization is the key process in ionization-current-based THz radiation generation in laser–gas interactions [18–20] and ionization-induced trapping in laser–plasma accelerators [21–27]. Laser-ionization-based trapping ⇑ Corresponding author. Tel.: +1 510 486 6632; fax: +1 510 486 7981. E-mail address: [email protected] (M. Chen). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.11.029
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
221
in plasma accelerators relies on the use of a high-Z gas target to generate electrons at the appropriate phase of the plasma wakefield. For this application, an accurate description of the ionization process is critical for predicting the laser–plasma accelerated beam properties. Particle-in-cell (PIC) [28,29] is a common numerical method used to study short-pulse laser–plasma interactions. The PIC numerical method models the Vlasov–Maxwell equations using a particle-grid method, such that the continuous phasespace fluid is sampled by a finite number of (Lagrangian) macro-particles and the Maxwell equations are discretized on a grid. Laser ionization has been implemented and studied in PIC codes, and several types of ionization formulae have been used in the literature (see, for example, Refs. [30,17]). In this paper, we focus on the numerical implementation in a PIC code of optical tunneling ionization based on the Ammosov–Delone–Krainov (ADK) model [31]. We compare the alternating current (AC) and direct current (DC) models and show that the DC model must be used to give accurate results in explicit PIC simulations. Multi level ionization in a single simulation step and energy conservation during ionization are considered. The effects of numerical parameters on the modeling of the laser ionization process are presented. We consider the effects of the resolution of the grid and the average number of numerical macro-particles per cell representing the phase-space distribution. A comparison of ionization injection studies using two PIC frameworks (VLPL [32] and VORPAL [33]) is also made.
2. Tunneling ionization formulae Tunneling ionization from the ground state of a Hydrogen atom in an electrostatic (DC) field was first solved schematically by Landau [34] for the dependence of the ionization rate on the applied external field. When the intensity of the external electric field satisfies Edc Ea (where eEa =me c2 ¼ a4 =re or Ea ’ 5:1 GV=cm), the ionization rate for a Hydrogen-like system is given by the Landau–Lifshitz DC tunneling ionization formula [34]:
W LL ¼ 4xa
" 5=2 3=2 # Ea U ion 2 Ea U ion ; exp 3 Edc U H Edc U H
ð1Þ
where xa ¼ a3 c=re ¼ 4:13 1016 s1 is the atomic unit frequency, U ion is the ionization potential, and U H ¼ 13:6 eV is the ionization potential of Hydrogen at the fundamental state. Here me is the electron mass, c is the speed of light, a ¼ 1=137 is the fine structure constant, and r e ¼ e2 =me c2 is the classical electron radius. For complex atoms in a static field (Edc ) the DC tunneling ionization formula is [31]
W DClm ¼ x
2 a C n l
ð2l þ 1Þðl þ jmjÞ!
U ion 2U H 2jmj ðjmjÞ!ðl jmjÞ!
" " 3=2 #2n jmj1 3=2 # Ea U ion 2 Ea U ion ; 2 exp 3 Edc U H Edc U H
ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where l; m are the electron’s orbital quantum number and its projection, respectively, n ¼ Z U H =U ion is the effective prin cipal quantum number, l ¼ n0 1 is the effective value of the orbital number (n0 is the effective principal quantum number of the ground state.), and Z is the ion charge number after ionization. The coefficients C n l can be found calculated using a semi-classical approximation (see Eq. (20) of Ref. [31]):
C 2n l ¼
n l þ1=2 22n 1 4e2 n l : ’ n Cðn þ l þ 1ÞCðn l Þ 2pn n2 l2 n þ l
ð3Þ
In the limit l n , Eq. (3) reduces to
C 2n l ¼
2n 1 2e : 2pn n
ð4Þ
The ionization rate of a hydrogen atom in an alternating (AC) electric field, e.g., a laser field [EL / cosðxtÞ, where x is the laser frequency], can be treated, in the appropriate limit, as the DC ionization rate averaged over a single period of the electric field oscillation. This model was proposed by Keldysh [35] and then later expanded to hydrogen-like atoms by Ammosov et al. [31] (referred to as the ADK model). The DC ionization rate is related to the AC ionization rate by
"
W AC
3=2 #1=2 3 E UH ¼ W DC : p Ea U ion
ð5Þ
Combining Eqs. (2), (4), and (5), yields the one cycle averaged ionization rate for a complex atom in an AC field. For a linearly polarized laser pulse EL , the ionization rate is [31]
W AClm ¼ xa
" !2n jmj1 3 1=2 2 2n 3 # 3n EL Z 2e 2 Ea Z ð2l þ 1Þðl þ jmjÞ! Ea Z 3 2 exp : 3 EL n EL n3 2n2 n pZ 3 Ea 2pn 2jmj ðjmjÞ!ðl jmjÞ!
ð6Þ
222
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
For the fundamental state l ¼ m ¼ 0, we obtain the ionization rate
pffiffiffi e 3=2 Z 2 Ea Z 3 W AC00 ¼ xa 3 4e 9=2 p EL n4 n
!2n 3=2
" exp
3 # 2 Ea Z : 3 EL n
ð7Þ
In the original ADK paper a typographic error gave an ionization rate (see Eq. (13a) of Ref. [31]) of p=e times the correct rate. This and other corrections, including errors in Eq. (20) of Ref. [31], have been pointed out by several authors [36,37]. However, the original formula with the typographic error has been used in several studies (e.g., Refs. [38,17]), and is pointed out here to prevent future confusion. In numerical modeling of laser–plasma interactions via PIC it is essential to use the correct ionization model. If the laser pulse is modeled as an envelope and the time step is much larger than the laser period, an ionization rate based on a cycleaveraged AC model should be used. However, in an explicit PIC simulation, the sub-wavelength scale is well-resolved and the simulation time step is much smaller than the laser period (dt T 0 , where dt is the simulation time step and T 0 ¼ c=k0 is the laser period), and within each simulation time step, the laser field can be viewed as a DC field. In this case the DC ionization rate model should be employed. To illustrate the differences in the ionization models, Fig. 1 shows numerical calculations comparing the DC model Eq. (2) and the AC model Eq. (6) in a calculation that resolves the laser period. Fig. 1(a) shows the ionization rate from N5þ to N6þ versus normalized electric field intensity a ¼ eEL =me xc. As shown in Fig. 1(a), the AC ionization model (as used in Ref. [17]) shows the same trend as the DC model; however, it always provides a lower ionization rate, as anticipated from Eq. (5). The ionization probability may also be considered within a finite time. Fig. 1(b) shows the ionization probability (from N5+ to N6+) after a single laser cycle (a ¼j a0 sinð2ptÞ j, with 0 6 t 6 1) versus laser intensity: The red curve shows the result by using the theoretical calculation of 1 exp½Wðt p ÞT 0 , where tp represents the time at the electric field peak. The black curve (overlapped by the red curve) shows the result from the DC integration model, and the blue curve shows the result from the AC integration model. As shown in Fig. 1(b), use of the AC model in this case has underestimated the ionization probability. Fig. 1(c) shows the difference of the two models for a temporally Gaussian pulse, with full-width-half-maximum (FWHM) length LFWHM ¼ 14:89cT 0 . Again, the AC model gives a lower ionization probability than the DC model. From this, we expect that, while results of studies using explicit codes with the AC ionization model (or using the uncorrected ADK coefficient) may give qualitatively reasonable results, the use of the corrected DC model in explicit simulations is essential to obtain accurate quantitative results.
(a)
(b)
(c)
Fig. 1. (a) Ionization rate from N5+ to N6+ versus electric field intensity a ¼ eEL =me cx. The laser frequency is x ¼ 2pc=k0 with k0 ¼ 0:8 lm. The solid black solid curve is based on the DC model and the dashed blue curve is based on the AC averaged model. (b) Ionization probability after one laser cycle calculated using different models: solid black curve integrates over W DClm , the dashed red curve uses the one-cycle-averaged formula 1 exp½Wðtp ÞT 0 , the blue dashed curve integrates over W AClm . (c) Ionization probability after a Gaussian pulse: solid black curve is integration over W DClm integration, blue dashed curve integrates over W AClm . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
223
3. Tunneling ionization model implementation in PIC codes In an explicit PIC simulation, the ionization process is generally implemented with the following prescription. At each time step, for each un-ionized particle, the local electric field and particle ionization parameters (such as ionization potential and charge state) are used to calculate the ionization rate WðtÞ. The particle ionization probability is then calculated by PðtÞ ¼ 1 exp½WðtÞdt, which is approximated as PðtÞ ’ WðtÞdt when PðtÞ 1. The code generates a random number p with a uniform distribution between 0 and 1; once p < PðtÞ, the particle will be set to be ionized, otherwise the particle will not be ionized and will be evaluated during the next time step. In the PIC code the particles are macro-particles representing the phase space fluid. The ionization probability cannot be accumulated to the next time step if the particle is not ionized; otherwise an over-counted ionization probability results. Although for most of our simulation cases, the time step of PIC simulations can satisfy PðtÞ 1, however, if the ionization potential is too low or a relatively high laser intensity is used, the above approximation can not be satisfied and high level ionization in a single PIC simulation time step may happen. A subroutine will be used to solve this issue. Let us consider an ion particle whose initial state is Zi and its highest charge is z. The evolution of the state population (N j ) is as following: i
dN =dt ¼ W i Ni dN
iþ1
ð8Þ
=dt ¼ W i Ni W iþ1 Niþ1
ð9Þ ð10Þ
... dN
z1
=dt ¼ W z2 Nz2 W z1 Nz1
ð11Þ
z
dN =dt ¼ W z1 Nz1 i
ð12Þ j
and the initial state is: N ðt ¼ 0Þ ¼ 1 and N ðt ¼ 0Þ ¼ 0 when j > i. Obviously during the whole ionization process (t 2 ½0; dt), P j the total number of all kinds of ions is constant: j¼z j¼i N ¼ 1. So only z i equations need to be solved. Depends on the specific dependance of W j on t different approximations can be used to simplify these equations. Generally a subroutine can be installed in the PIC code to solve this equation array by using the typical fourth order Runge–Kutta method. Once one knows the state population (N j ) at the end of time step t ¼ dt, a random number p with a uniform distribution between 0 and 1 will P P j k be generated to decide the final ionization state. When j¼k1 N j < p < j¼k j¼i j¼i N , the particle will be set to be ionized to Z and k i electrons will be generated. Here we neglect the different response of the different valence electrons in such a single time step, which is reasonable when the time step dt of the PIC simulation is small enough. Besides the high level ionization process, ionization induced energy loss (Eloss ion ) is another issue to be considered. Although in some simulation cases such as laser short–gas interaction, the energy loss due to ionization can be neglected compared with the high laser energy. However, in some special cases such as long distance laser propagation in neutral gas or laser high density high Z neutral matter interaction, energy loss due to ionization should be considered. Currently, to retain the energy conservation, an approach using artificial ionization current (~ jion ) is suggested. Such method has been discussed by Kemp et al. [17] before. The artificial ionization current is directed along the local electric field ~ E, and its magnitude is set ~2 to satisfy the energy conservation law: ~ jion ¼ ~ ERk Dnek ðEloss ionk þ eki Þ=dtjEj . Here k represents the different ionization state, Dnek th represents the number density of the k ionized electrons and eki represents their initial energies just after ionization which is usually set to be zero. This current is then added to the physical current in each cell, after charge and current are computed from particle positions and velocities and before the integration of Maxwell’s equations. Such process should be done for every time step and each cell once the ionization process happens. The above analysis also indicates that, to accurately model the ionization process, random number generation, selection of the macro-particle number representing the ionizing atoms, and the choice of time step size are all critical numerical parameters. In the following we examine the effect of these numerical parameters on the ionization modeling using the VLPL PIC code [32]. In the VLPL code, the ionization rate formula of the DC ionization model Eq. (2) is implemented. A library source random number generator is used to produce random numbers uniformly distributed in the region ½0; 1. To study the effects of time step and macro-particles per cell on the ionization probability calculation, we consider the ionization of a Nitrogen gas by a linearly polarized laser pulse. To exclude other physical effects, the number density of the initial neutral Nitrogen atoms is 105 nc , where nc ¼ pme c2 =e2 k20 is the critical density (nc ’ 1:7 1021 cm3 for a laser of wavelength of k0 ¼ 0:8 lm) and the simulation is in a one-dimensional (1D) geometry. The normalized laser intensity is a0 ¼ eEL =me xc ¼ 1:8 and the pulse length is L0 ¼ 10cT 0 , with the temporal shape of electric field defined as / expðt2 =L20 Þ. The neutral Nitrogen atoms are initially uniformly distributed within the region of 5k0 —35k0 . No other particle elements are used. For such a laser pulse the highest ionization rate for Nitrogen N5+ is 4:68 1014 s1 at the peak of the laser pulse, as shown in Fig. 1(a). Fig. 2(a) and (b) shows the distribution of ionization probability of N5+ (fraction of ions in charge states N6+ or N7+) along the laser propagation direction. The laser propagates from the left to the right and the grey curve shows the laser field (with peak at x ¼ 15k0 ). The red dashed curve shows the theoretical ionization calculated using Eq. (2) and assuming N5+ is uniformly pre-ionized before the laser arrives. The assumption is reasonable since Nitrogen is easily ionized to the fifth valence
224
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
δ
(a)
δ δ
(b)
Fig. 2. Ionization probability distribution of N5+ along the laser propagation direction x=k0 . (a) The grey curve is the normalized laser electric field with peak field at x ¼ 15k0 . The red dashed curve shows the theoretical result, the green dotted curve is calculate from the PIC code with dt ¼ 0:1T 0 and 40 macroparticles per cell. (b) The black solid curve uses dt ¼ 0:01T 0 and 4 macro-particles per cell, the light blue solid curve used dt ¼ 0:005T 0 and 1 macro-particle per cell. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 3. Example of the typical laser–plasma accelerator wake structure and ionization injection scheme: laser electric field (grey curve) along its propagation direction (x), the electric field of the wake (red curve), the momentum of the background plasma electron fluid (dashed blue curve), the trapping momentum threshold of the electrons in the wakefield (dotted green curve), and the momentum gap between the trapping threshold and the background electron momentum (black curve). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
state at the front of the laser pulse. The green dotted curve in Fig. 2 shows the sum density of the N6+ and N7+ normalized by the initial Nitrogen density N 0 , which represents the total ionization probability of N5+, calculated from the PIC code with a time resolution of dt ¼ 0:1T 0 ’ 0:267 fs and 40 macro-particles per cell. The ionization probability has been underestimated in the PIC code by 21.4%. At the peak of the laser pulse, Wdt ’ 0:125 and the approximation of 1 expðWdtÞ Wdt gives a relative error of +6.38%, i.e., the error in the ionization calculation does not come from the small exponent approximation, rather it is largely due to the sampling of the electric fields.
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
225
Since the ionization mainly happens at the peak of each cycle, if the resolution is too course, the sampling of fields (i.e., the resolution) may not be sufficiently accurate. To study the effect of resolution, we considered increased temporal resolution with the total macro-particle number within a fixed length constant. Increased temporal resolution is equivalent to increased longitudinal spatial resolution, since in the 1D PIC code dx=k0 ¼ dt=T 0 . The black solid curve in Fig. 2(b) shows the ionization probability calculated from the PIC code with dt ¼ 0:01T 0 and 4 macro-particles per cell. Better agreement with the theoretical curve (dashed red curve) is achieved with the increased resolution. The increase in noise is due to the fewer macro-particles per cell (i.e., larger weight given to each macro-particle). For example, if the resolution is increased to dt ¼ 0:005T 0 with 1 macro-particle per cell, increased simulation noise is observed (light blue solid curve in Fig. 2(b)). The ionization probability calculation is sensitive to the resolution, and accurate modeling of the ionization process the time step of each should be sufficiently small. In general, the resolution criterion is more strict than that due to the ionization rate’s sensitivity to the electric field strength WðtÞdt 1. The number of macro-particle per cell affects the simulation noise level. The larger the number per cell, the smaller weight of each macro-particle, and the simulation noise will be less. When modeling ionization with PIC, evaluating the convergence of the resolution is essential for an accurate description of ionization process.
(a)
(b)
(c)
Fig. 4. (a) Evolution of trapped (with momentum larger than 20 MeV/c) electron number versus laser propagation distance: VLPL using DC model (solid black squares), VORPAL using DC model (solid red circles), VLPL using AC model (open black squares), and VLPL using DC model and 4 macro-particles per cell (blue solid triangles). (b) Evolution of beam mean energy (open points) and rms energy spread (solid points) versus laser propagation distance using VLPL (black square) and VORPAL (red circle). (c) Beam transverse momentum spread evolution versus laser propagation distance using VLPL (black square) and VORPAL (red circle). In these simulations, 1 macro-particle per cell for the pre-ionized electrons and 6 macro-particles per cell for Nitrogen atoms are used. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
226
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
4. Comparison of ionization injection studies The implementation of the DC ionization model in PIC codes described above has been carried out in two PIC codes VLPL [32] and VORPAL [33]. In this section we compare results obtained in both PIC codes using the same DC ionization model. The case studied is for ionization injection into a laser–plasma accelerator, and the final simulated beam characteristics are compared. High-Z gas ionization injection has recently been employed experimentally as an effective electron injection method in the laser–plasma accelerator [22–26]. The detailed physics of ionization injection has been describe in Ref. [27]. Fig. 3 shows a typical wake structure excited by an intense laser and the injection criterion based on 1D nonlinear cold plasma wake theory at a fixed time. A laser pulse (grey curve) is used to generate a wake in a mixed plasma of pre-ionized electrons and neutral Nitrogen atoms. The red solid curve shows the longitudinal electric field of the wake when only the initial pre-ionized electron plasma is included. Since the concentration of Nitrogen is small (typically 20 MeV/c) as calculated in the two PIC codes using the DC ionization model. Using the AC model gives a lower injection number, which is approximately half that of the DC model result for these parameters. The calculated injected beam mean energy and root-mean-square (r.m.s.) energy spread are also similar for both codes when the DC ionization model is used (shown in Fig. 4(b)). For ionization-induced injection, the beam transverse momentum spread is mainly due to the residual transverse momentum as a result of the ionization process. In the 1D case, canonical momentum conservation ensures ~ p? ¼ ~ aðt i Þ, where t i is the ionization time. Therefore, the electrons have a transverse momentum equal to the laser vector potential at ionization. After the injection, the transverse momentum evolution depends on the transverse wakefield structure. Fig. 4(c) shows the FWHM transverse momentum distribution as calculated by the two codes, with differences less than 10%. This level of agreement is reasonable, considering the numerical noise generated by the ionization (i.e., the low macro-particles per cell). Increasing the number of macro-particles per cell for the background plasma electrons (from 1 to 4 macro-particles per cell), improves the agreement, indicating the source of errors is numerical noise in the PIC code. It also shows the electrons’ transverse momentum spread not only depends on the correct simulation of ionization process, it depends on the correct wake calculation too. Typical beam momentum distributions are shown in Fig. 5 (where px ; py are the longitudinal and transverse beam momenta, respectively). Differences in the beam transverse momenta distribution are the result of differences in the transverse wakefield structure calculated by the two codes. The above comparison provides a benchmark of the ionization algorithms implemented in VLPL and VORPAL for parameters relevant to studies of ionization injection in laser–plasma accelerators. 5. Summary and conclusions Implementation of laser tunneling ionization rate formulae in explicit PIC codes has been examined. The originally published ionization rate formulae were reviewed, and corrections noted. It was shown, by comparison to theoretical predictions, that the DC form of the ionization rate provides improved accuracy. A subroutine is used to deal with the high level ionization in a single time step and an ionization current is used to consider energy loss due to ionization. The effect of PIC numerical parameters, such as grid resolution and number of macro-particles per cell, on the calculated ionization probability was studied. It was found that the ionization probability is sensitive to the grid resolution, and increased number of macro-particles per cell can reduce the simulation noise. Modeling the ionization injection in a laser–plasma accelerator was used to benchmark the ionization algorithm implemented in two PIC codes. Acknowledgments The authors thank Carlo Benedetti and Stepan Bulanov for useful discussions. This work was supported by the Director, Office of Science, Office of High Energy Physics, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and by a Scientific Discovery through Advanced Computing Project ‘‘Advanced Computing for 21st Century Accelerator Science and Technology’’. Computational resources of the National Energy Research Scientific Computing Center were used to perform the simulations. References [1] E. Esarey, C.B. Schroeder, W.P. Leemans, Physics of laser-driven plasma-based electron accelerators, Rev. Mod. Phys. 81 (2009) 1229–1285. [2] C.G.R. Geddes, C. Tóth, J. van Tilborg, E. Esarey, C.B. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, W.P. Leemans, Production of high-quality electron bunches by dephasing and beam loading in channeled and unchanneled laser plasma accelerators, Phys. Plasmas 12 (2005) 056709. [3] T. Hosokai, K. Kinoshita, A. Zhidkov, K. Nakamura, T. Watanabe, T. Ueda, H. Kotaki, M. Kando, K. Nakajima, M. Uesaka, Effect of a laser prepulse on a narrow-cone ejection of mev electrons from a gas jet irradiated by an ultrashort laser pulse, Phys. Rev. E 67 (3) (2003) 036407. [4] T. Hosokai, K. Kinoshita, T. Ohkubo, A. Maekawa, M. Uesaka, A. Zhidkov, A. Yamazaki, H. Kotaki, M. Kando, K. Nakajima, S.V. Bulanov, P. Tomassini, A. Giulietti, D. Giulietti, Observation of strong correlation between quasimonoenergetic electron beam generation by laser wakefield and laser guiding inside a preplasma cavity, Phys. Rev. E 73 (3) (2006) 036407. [5] E. Esarey, A. Ting, P. Sprangle, Frequency shifts induced in laser pulses by plasma waves, Phys. Rev. A 42 (1990) 3526. [6] E. Esarey, G. Joyce, P. Sprangle, Frequency upshifting of laser pulses by copropagating ionization fronts, Phys. Rev. A 44 (1991) 3908. [7] W.P. Leemans, C.E. Clayton, W.B. Mori, K.A. Marsh, P.K. Kaw, A. Dyson, C. Joshi, J.M. Wallace, Experiments and simulations of tunnel-ionized plasmas, Phys. Rev. A 46 (1992) 1091–1105. [8] W.P. Leemans, C.E. Clayton, W.B. Mori, K.A. Marsh, A. Dyson, C. Joshi, Plasma physics aspects of tunnel-ionized gases, Phys. Rev. Lett. 68 (1992) 321– 324. [9] R.L. Savage, C. Joshi, W.B. Mori, Frequency upconversion of electromagnetic radiation upon transmission into an ionization front, Phys. Rev. Lett. 68 (1992) 946–949, http://dx.doi.org/10.1103/PhysRevLett.68.946. URL: . [10] B. Cros, X. Xu, T. Tsukada, N. Yugami, Y. Nishida, G. Matthieussent, Ionization-induced frequency up-shift of a high-power microwave interacting with a plasma, Phys. Plasmas 4 (8) (1997) 2837–2844, http://dx.doi.org/10.1063/1.872416. URL: .
228
M. Chen et al. / Journal of Computational Physics 236 (2013) 220–228
[11] P. Sprangle, E. Esarey, J. Krall, Self-guiding and stability of intense optical beams in gases undergoing ionization, Phys. Rev. E 54 (4, Part b) (1996) 4211– 4232, http://dx.doi.org/10.1103/PhysRevE.54.4211. [12] E. Esarey, P. Sprangle, J. Krall, A. Ting, Self-focusing and guiding of short laser pulses in ionizing gases and plasmas, IEEE J. Quantum Electron. 33 (11) (1997) 1879–1914, http://dx.doi.org/10.1109/3.641305. [13] P. Sprangle, E. Esarey, B. Hafizi, Propagation and stability of intense laser pulses in partially stripped plasmas, Phys. Rev. E 56 (5, Part b) (1997) 5894– 5907, http://dx.doi.org/10.1103/PhysRevE.56.5894. [14] P. Sprangle, E. Esarey, B. Hafizi, Intense laser pulse propagation and stability in partially stripped plasmas, Phys. Rev. Lett. 79 (6) (1997) 1046–1049, http://dx.doi.org/10.1103/PhysRevLett.79.1046. [15] T. Antonsen, Z. Bian, Ionization induced scattering of short intense laser pulses, Phys. Rev. Lett. 82 (18) (1999) 3617–3620, http://dx.doi.org/10.1103/ PhysRevLett.82.3617. [16] Z. Bian, T. Antonsen, Ionization instabilities of an electromagnetic wave propagating in a tenuous gas, Phys. Plasmas 8 (7) (2001) 3183–3194, http:// dx.doi.org/10.1063/1.1377612. [17] A.J. Kemp, R.E.W. Pfund, J.M. ter Vehn, Modeling ultrafast laser-driven ionization dynamics with monte carlo collisional particle-in-cell simulations, Phys. Plasmas 11 (12) (2004) 5648–5657. [18] M. Kre, T. Löffler, M.D. Thomson, R. Dörner, H. Gimpel, K. Zrost, T. Ergler, R. Moshammer, U. Morgner, J. Ullrich, H.G. Roskos, Determination of the carrier-envelope phase of few-cycle laser pulses with terahertz-emission spectroscopy, Nat. Phys. 2 (2006) 327–331. [19] K.-Y. Kim, J.H. Glownia, A.J. Taylor, G. Rodriguez, Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields, Opt. Express 15 (8) (2007) 4577–4584. [20] M. Chen, A. Pukhov, X.-Y. Peng, O. Willi, Theoretical analysis and simulations of strong terahertz radiation from the interaction of ultrashort laser pulses with gases, Phys. Rev. E 78 (2008) 046406. [21] M. Chen, Z.-M. Sheng, Y.-Y. Ma, J. Zhang, Electron injection and trapping in a laser wakefield by field ionization to high-charge states of gases, J. Appl. Phys. 99 (5) (2006) 056109. [22] A. Pak, K.A. Marsh, S.F. Martins, W. Lu, W.B. Mori, C. Joshi, Injection and trapping of tunnel-ionized electrons into laser-produced wakes, Phys. Rev. Lett. 104 (2) (2010) 025003. [23] C. McGuffey, A.G.R. Thomas, W. Schumaker, T. Matsuoka, V. Chvykov, F.J. Dollar, G. Kalintchenko, V. Yanovsky, A. Maksimchuk, K. Krushelnick, V.Y. Bychenkov, I.V. Glazyrin, A.V. Karpeev, Ionization induced trapping in a laser wakefield accelerator, Phys. Rev. Lett. 104 (2) (2010) 025004. [24] C.E. Clayton, J.E. Ralph, F. Albert, R.A. Fonseca, S.H. Glenzer, C. Joshi, W. Lu, K.A. Marsh, S.F. Martins, W.B. Mori, A. Pak, F.S. Tsung, B.B. Pollock, J.S. Ross, L.O. Silva, D.H. Froula, Self-guided laser wakefield acceleration beyond 1 GeV using ionization-induced injection, Phys. Rev. Lett. 105 (2010) 105003. [25] B.B. Pollock, C.E. Clayton, J.E. Ralph, F. Albert, A. Davidson, L. Divol, C. Filip, S.H. Glenzer, K. Herpoldt, W. Lu, K.A. Marsh, J. Meinecke, W.B. Mori, A. Pak, T.C. Rensink, J.S. Ross, J. Shaw, G.R. Tynan, C. Joshi, D.H. Froula, Demonstration of a narrow energy spread, 0:5GeV electron beam from a two-stage laser wakefield accelerator, Phys. Rev. Lett. 107 (2011) 045001, http://dx.doi.org/10.1103/PhysRevLett.107.045001. URL: . [26] J.S. Liu, C.Q. Xia, W.T. Wang, H.Y. Lu, C. Wang, A.H. Deng, W.T. Li, H. Zhang, X.Y. Liang, Y.X. Leng, X.M. Lu, C. Wang, J.Z. Wang, K. Nakajima, R.X. Li, Z.Z. Xu, All-optical cascaded laser wakefield accelerator using ionization-induced injection, Phys. Rev. Lett. 107 (2011) 035001, http://dx.doi.org/10.1103/ PhysRevLett.107.035001. URL: . [27] M. Chen, E. Esarey, C.B. Schroeder, C.G.R. Geddes, W.P. Leemans, Theory of ionization-induced trapping in laser–plasma accelerators, Phys. Plasmas 19 (3) (2012) 033101. [28] C.K. Birdsall, A.B. Langdon, V. Vehedi, J.P. Verboncoeur, Plasma Physics via Computer Simulations, Adam Hilger, Bristol, 1991. [29] R. Hockney, J. Eastwood, Computer Simulation using Particles, Taylor & Francis Group, New York, 1988. [30] D.L. Bruhwiler, D.A. Dimitrov, J.R. Cary, E. Esarey, W. Leemans, R.E. Giacone, Particle-in-cell simulations of tunneling ionization effects in plasma-based accelerators, Phys. Plasmas 10 (5) (2003) 2022–2030. [31] M.V. Ammosov, N.B. Delone, V.P. Krainov, Sov. Phys. JETP 64 (1986) 1191. [32] A. Pukhov, Three-dimensional electromagnetic relativistic particle-in-cell code VLPL (virtual laser plasma lab), J. Plasma Phys. 61 (1999) 425–433. [33] C. Nieter, J. Cary, VORPAL: a versatile plasma simulation code, J. Comput. Phys. 196 (2004) 448. [34] L. Landau, E. Lifshitz, Quantum Mechanics, Butterworth-Heinemann, New York, 1965. [35] L.V. Keldysh, J. Exp. Theor. Phys. (USSR) 47 (1964) 1945. [36] X. Lu et al, Acta Opt. Sin. 8 (1999) 1014. [37] F.A. Ilkov, J.E. Decker, S.L. Chin, Ionization of atoms in the tunnelling regime with experimental evidence using hg atoms, J. Phys. B At. Mol. Opt. Phys. 25 (19) (1992) 4005–4020. [38] B.M. Penetrante, J.N. Bardsley, Residual energy in plasmas produced by intense subpicosecond lasers, Phys. Rev. A 43 (1991) 3100–3113. [39] C.B. Schroeder, E. Esarey, B.A. Shadwick, W.P. Leemans, Trapping, dark current, and wavebreaking in nonlinear plasma waves, Phys. Plasmas 13 (3) (2006) 033103.
