Efficient numerical algorithm for multi-level ionization of high-atomic-number gases

Aiqi Cheng Roman Samulyak Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY 11794, United States Computational Science Initiative, Brookhaven National Laboratory, Upton, NY 11973

Abstract

An efficient numerical algorithm for the laser driven multi-level ionization of high-atomic-number gases is proposed and implemented in an electromagnetic particle-in-cell code SPACE. The algorithm is based on analytical solutions to the system of differential equations describing ionization evolution. Using analytical solutions resolves the multiscale issue of ionization due to different characteristic time scales of ionization processes and the main code time step. Algorithm efficiency is improved by using a locally-reduced system of differential equations. The effects of the orbital quantum numbers and their projections have been examined. The algorithm is applied to the study of ionization injection of electrons into laser-driven plasma wakefields.

keywords:

Laser plasma interaction , ADK model , multi-level tunneling ionization , Particle-in-cell , Ionization injection

^†^†journal: Journal of Computational Physics

1 Introduction

An accurate numerical model for ionization of gases by high-power lasers is critical for numerical simulations of laser-plasma interactions and for the laser-plasma wakefield acceleration (LPWA) of particles in particular[1, 2, 3]. Although modern lasers used for LPWA are capable of completely ionizing hydrogen gas during several periods of laser filed oscillation and within a small fraction of the laser pulse duration, the assumption of a pre-ionized plasma is not accurate for laser intensities close to the tunneling ionization threshold [4, 5]. Ionization effects such as creation of plasma through tunneling ionization, frequency up-shifting of the laser, harmonic generation, and scattering instabilities play a major role in the laser wakefield acceleration process. An accurate modeling of the laser ionization of gases significantly improves the resolution of such processes.

The large interest for multi-level ionization of high-atomic-number gases in the context of LPWA is driven by ideas of the ionization injection of electrons into plasma wakes [9, 10, 11], a technique that may significantly reduce the complexity of the LPWA method and improve the quality of resulting accelerated electron beams. While the ionization of hydrogen, a typical medium for LPWA experiments, is easy to resolve numerically, multi-level ionization of high-atomic-number gases presents a set of more complex modeling challenges. These include the multiscale nature of ionization, complex atomic states characterized by sets of orbital quantum numbers and their projections, and various software development issues such as a significant increase of the memory usage if the evolution of multiple ionization states is tracked at every computational mesh block in the simulation domain. In this paper, we propose a numerical alorithm for electromagnatic Particle-in-Cell (PIC) codes that improves the accuracy and numerical efficiency compared to previously published models [6]. PIC [7, 8] is a common method for the simulation of laser-plasma interactions and LPWA which solves the Vlason-Maxwell equations using a hybrid particle-grid method. In PIC, the electromagnetic Maxwell equations are discretized and updated on a grid while the distribution functions of charged particles is sampled by a set of Lagrangian particles, the motion of which is coupled to the mesh via currents.

The algorithm has been implemented in SPACE, a parallel, relativistic, three-dimensional particle-in-cell code developed for the simulation of electromagnetic fields, relativistic particle beams, and plasmas [13]. A distinct feature of SPACE is a software module for atomic physics transformations.

This paper is organized as follows. In Section 2, we present the main equations governing laser tunneling ionization of high-atomic-number gases. Section 3 describes details of the mathematical model and presents verification tests. We conclude this paper by Section 4 with a summary of our results and plans for the future work.

2 Laser ionization of high-Z gases

The Ammosov–Delone–Krainov (ADK) model [12, 6] for tunneling ionization in high-atomic-number gasses derives the following expression for the probability of ionization

	$\displaystyle W_{lm}=$		$\displaystyle\omega_{\alpha}\sqrt{\frac{3{n^{}}^{3}E_{L}}{\pi Q^{3}E_{a}}}\frac{Q^{2}}{2{n^{}}^{2}}\left(\frac{2e}{n^{}}\right)^{2n^{}}\exp{\left[-\frac{2E_{a}}{3E_{L}}\left(\frac{Q}{n^{*}}\right)^{3}\right]}\times$		(1)
			$\displaystyle\frac{(2l+1)(l+\|m\|)!}{2\pi n^{}2^{\|m\|}(\|m\|)!(l-\|m\|)!}\left(2\frac{E_{a}Q^{3}}{E_{L}{n^{}}^{3}}\right)^{2n^{*}-\|m\|-1}.$		(1)

Here $\omega_{a}=\alpha^{3}c/r_{e}=4.13\times 10^{16}s^{-1}$ is the atomic unit frequency, $E_{a}=510\ GV/m$ is the atomic unit of Electric field, and $E_{L}\ (GV/m)$ represents the strength of the local laser electric field. $l$ and $m$ are the electron’s orbital quantum number and projection of the orbital quantum number respectively, $n^{*}=Q\sqrt{U_{H}/U_{ion}}$ is the effective principle quantum number, and $U_{H}$ and $U_{ion}$ are the ionization potential of hydrogen and ion of the material being ionized, respectively. $Q$ is the charge number after ionization.

Denoting the number density of atoms in the neutral ground state as $n_{0}$ and the number density of ions ionized to the state $j+$ as $n^{j+}$ , the system of equation governing the ionization dynamics is as follows

$\displaystyle\frac{dn_{0}}{dt}$	$\displaystyle=-W_{0}*n_{0}(t),$
$\displaystyle\frac{dn^{+}}{dt}$	$\displaystyle=W_{0}n_{0}(t)-W_{1}n^{+}(t),$
$\displaystyle\frac{dn^{(Z-1)+}}{dt}$	$\displaystyle=W_{Z-2}n^{(Z-2)+}(t)-W_{Z-1}n^{(Z-1)+}(t),$	(2)
$\displaystyle\frac{dn^{Z+}}{dt}$	$\displaystyle=W_{Z-1}*n^{(Z-1)+},$

where $Z$ denotes the gas atomic number. We would like to note that because overall time scales of LPWA simulations are short compared to recombination processes in gases, the recombination terms were not included in the governing system of equations.

This system of equations must be solved independently at each numerical grid cell of the computational domain which presents a set of challenging numerical issues. First, the system includes multiple time scales associated with short characteristic times for ionization dynamics of each individual level and the overall characteristic time of the electromagnetic problem being studied. Second, recording the evolution dynamics of all individual ionization levels in each computational mesh cell requires a significant increase of the total memory. In the next Section, we propose an efficient numerical algorithm for multi-level ionization of high-atomic-number gases that resolves these and other issues, as well as a set of verification tests for the algorithm.

3 Numerical model and its verification

We start with writing equations (2) in the following matrix form:

\begin{bmatrix}\dot{n_{0}}\\ \dot{n}^{+}\\ \dot{n}^{2+}\\ ...\\ \dot{n}^{(Z-1)+}\\ \dot{n}^{Z+}\end{bmatrix}=\begin{bmatrix}-W_{0}&0&0&...&0&0\\ W_{0}&-W_{1}&0&...&0&0\\ 0&W_{1}&-W_{2}&...&0&0\\ \vdots&...&...&\ddots&\vdots&\vdots\\ \vdots&...&...&W_{Z-2}&-W_{Z-1}&0\\ 0&...&...&0&W_{Z-1}&0\end{bmatrix}\begin{bmatrix}n_{0}\\ n^{+}\\ n^{2+}\\ ...\\ n^{(Z-1)+}\\ n^{Z+}\end{bmatrix}

(3)

Since the ionization probabilities contained in the matrix depend on the local electric field values, our model is based on the natural assumption that ionization probabilities do not change within each time step of the overall electromagnetic simulation which updates states of the electromagnetic field. Assuming that the ionization probabilities are updated at the beginning of each time step of the overall simulation denoted at $t$ , we compute the analytical solution to (3) at the end of the time step $t+dt$ in terms of the eigenvalues $\lambda_{i}=-W_{i}(t)$ and eigenvectors $\vec{V}_{j}$ of the corresponding matrix. The general solution is:

n^{i+}(t+dt)=\sum\limits_{j=0}^{Z}C_{j}\vec{V}_{j}(i)e^{-W_{j}(t)dt}

(4)

The initial conditions at each time step $n^{j+}(t)$ , $j=0,\ldots,Z$ are used to compute coefficients $C_{j}$ .

\begin{bmatrix}\vec{V}_{0}&|&\vec{V}_{1}&...&|&\vec{V}_{Z}\end{bmatrix}\begin{bmatrix}C_{0}\\ C_{1}\\ \vdots\\ C_{Z}\end{bmatrix}=\begin{bmatrix}n_{0}(t)\\ n^{+}(t)\\ \vdots\\ n^{Z+}(t)\end{bmatrix}

(5)

Analytic solutions are computed independently in every computational cell and the population of ionization levels $n^{i+}$ is updated. Using analytical solution efficiently resolves the multiscale nature of ionization evolution and it eliminates the need for fractional time stepping for a numerical solver for the ODE system (2), as implemented in some electromagnetic codes. Explicit expressions for the analytical solutions subject to another simplifying approximations that reduces the total storage will be presented later in this Section.

The second challenge is associated with strong dependence of ionization probabilities on the orbital quantum numbers and their projections for electrons being ionized. A schematic of the electron structure of krypton, the medium of interest for the two-color injection in laser plasma wakefield acceleration, is shown in Figure 1.

Refer to caption — Figure 1: Schematic the electron structure of krypton

In the present implementation, we assume that electrons are ionized starting from out-most orbitals and their ionization probabilities are computed according to the quantum numbers and their projections as indicated in Figure 1. However, we also investigated other scenarios in which after removing the first two electrons from the $4p$ orbital, the unoccupied p-orbitals coupled to produce five-fold degenerate singlet- $d$ , nondegenerate singlet- $s$ , and triplet- $s$ states, and similar transformations in the lower orbitals. Since we observed negligibly small changes in the overall ionization dynamics, we compute probabilities based on quantum numbers presented in Figure 1.

The final challenge is associated with a significant increase of memory allocation for computing ionization processes. To store all ionization levels in krypton, 36 floating point numbers must be added to each computational cell. Of course, the total number of ionization levels could be reduced if ionization beyond some $N+$ , $N<Z$ , is unlikely at given laser parameters. This, however, would require a preliminary study for each new simulation setting, and fine-tuning of the ionization model for every gaseous material involved in the laser-matter interaction.

We performed detailed studies of the ionization dynamics of various ionization levels for different laser settings and concluded that only 2-3 intermediate levels are populated at any moment of time, with the population of lower and higher levels effectively approaching zero. Figure 2(a) shows the evolution of ionization levels in krypton interacting with a $CO_{2}$ laser pulse. A 4J laser pulse with 2ps duration was used and the laser wavelength was 9.2 $\mu m$ . The data was collected from a point on the laser propagation axis. The highly oscillatory gray line shows, in normalized units, the laser amplitude evolution at the given point. We observe that the density of neutrals reduces to zero within 0.35 ps after the laser pulse arrival. Due to a high ionization probability, this happens at a low laser field amplitude. The density of $Kr+$ rapidly grows and then reduces to zero at about 0.5 ps. Higher ionization levels undergo a similar dynamics. As the laser amplitude reduces from its maximum after 2.5 ps, the population of the ionization level $Kr^{10+}$ slowly reduces and the levels $Kr^{11+}$ and $Kr^{12+}$ slowly grow but their population remains at the order of magnitude of 1%. Similar dynamics occurs in other locations of the computational domain, off the laser axis, except that the maximum ionization level is lower due to lower laser field intensities. We observe that at any time, the effectively non-zero population is present for only 2-3 ionization levels.

This observation allowed us to replace the full system of equations for all ionization levels by a reduced-order system containing only 3 or 4 ionization levels which is specified independently for each computational cell. The active ionization levels are shifted to the higher ones independently in each computational cell when the population of the lowest level drops below a prescribed threshold. Figure 2 presents a comparison of the evolution of ionization levels along the laser axis computed by analytic solutions to the full system of ODE’s (3). The reduced 4-level system of equations describes this process very accurately. We do not plot the corresponding solution obtained with the 4-level system since there is no visible difference with the plot 2(a) obtained with the full system of ODE’s. Instead, we plot the evolution of solution errors for the reduced 4-level system in Figure 2(b) and observe that errors are of the order $10^{-4}$ for the intermediate states and for the final dominant ionization state $Kr^{10+}$ . Errors of the partially ionized states $Kr^{11+}$ and $Kr^{12+}$ are on the order of $10^{-3}$ , but these states contribute very small quantity of electrons to the final superposition state. Since errors in the final ionization states partially compensate each other, the relative error of the electron density, depicted in Figure 2(c), demonstrates high overall accuracy of the method.

The analytical solution for the reduced 4-level system is:

$\displaystyle n_{p}(t+dt)$	$\displaystyle=n_{p}(t)e^{-W_{p}dt}$	(6)
$\displaystyle n_{p+1}(t+dt)$	$\displaystyle=n_{p}(t)\frac{W_{p}}{W_{p}-W_{p+1}}(e^{-W_{p+1}dt}-e^{-W_{p}dt})+n_{p+1}(t)*e^{W_{p+1}dt}$	(7)
$\displaystyle n_{p+2}(t+dt)$	$\displaystyle=n_{p}(t)*\left(\frac{W_{p}W_{p+1}e^{-W_{p+2}dt}}{(W_{p}-W_{p+2})(W_{p+1}-W_{p+2})}-\frac{W_{p}W_{p+1}e^{-W_{p+1}dt}}{(W_{p}-W_{p+1})(W_{p+1}-W_{p+2})}\right.$
	$\displaystyle\left.+\frac{W_{p}W_{p+1}e^{-W_{p}dt}}{(W_{p}-W_{p+1})(W_{p}-W_{p+2})}\right)$
	$\displaystyle+n_{p+1}(t)*\left(\frac{(W_{p}W_{p+1}-W_{p+1}W_{p+2})e^{-W_{p+2}dt}}{(W_{p}-W_{p+2})(W_{p+1}-W_{p+2})}-\frac{(W_{p}W_{p+1}-W_{p+1}^{2})e^{-W_{p+1}dt}}{(W_{p}-W_{p+1})(W_{p+1}-W_{p+2})}\right)$
	$\displaystyle+n_{p+2}(t)*\left(\frac{(W_{p+2}^{2}+W_{p}W_{p+1}-W_{p}W_{p+2}-W_{p+1}W_{p+2})e^{-W_{p+2}dt}}{(W_{p}-W_{p+2})(W_{p+1}-W_{p+2})}\right)$	(8)
$\displaystyle n_{p+3}$	$\displaystyle=n_{p}(t)*\left(1+\frac{W_{p}W_{p+1}e^{-W_{p+2}dt}}{(W_{p}-W_{p+2})(W_{p+1}-W_{p+2})}+\frac{W_{p}W_{p+2}e^{-W_{p+1}dt}}{(W_{p}-W_{p+1})(W_{p+1}-W_{p+2})}\right.$
	$\displaystyle\left.-\frac{W_{p+1}W_{p+2}e^{-W_{p}dt}}{(W_{p}-W_{p+1})(W_{p}-W_{p+2})}\right)$
	$\displaystyle+n_{p+1}(t)*\left(1+\frac{(W_{p}W_{p+1}-W_{p+1}W_{p+2})e^{-W_{p+2}dt}}{(W_{p}-W_{p+2})(W_{p+1}-W_{p+2})}+\frac{(W_{p}W_{p+2}-W_{p+1}W_{p+2})e^{-W_{p+1}dt}}{(W_{p}-W_{p+1})(W_{p+1}-W_{p+2})}\right)$
	$\displaystyle+n_{p+2}(t)*\left(1+\frac{(W_{p+2}^{2}+W_{p}W_{p+1}-W_{p}W_{p+2}-W_{p+1}W_{p+2})e^{-W_{p+2}dt}}{(W_{p}-W_{p+2})(W_{p+1}-W_{p+2})}\right)+n_{p+3}(t),$	(9)

where $p$ denoted the lowest ionization level of the system. Initially, $p=0$ uniformly in the computational domain, and the value of $p$ increases non-uniformly after the cut-off criteria is met in a specific computational cell.

The corresponding reduced 3-level system of ODE’s has also been tested. Since numerical simulations with the 3-level system showed increased errors in certain parts of the computational domain, the 4-level system was selected as the most optimal solution that optimizes the accuracy and numerical cost.

4 Implementation in SPACE and numerical examples

The numerical algorithm for multi-level ionization of high-atomic-number gases was implemented as a module in SPACE [13], a parallel, relativistic, three-dimensional particle-in-cell code developed for the simulation of electromagnetic fields, relativistic particle beams, and plasmas. The algorithm proposed in this paper extends a set of atomic physics algorithms previously implemented in SPACE that resolve single-electron ionization, recombination, and electron attachment to dopants in dense neutral gases. These algorithms were used for longer-time-scale processes compared to LPWA and can be illustrated by SPACE simulations, performed in support of the experimental program on the high-pressure hydrogen gas-filled radio frequency cavity in the Mucool Test Area at Fermilab [14, 15]. The previous SPACE implementation also employed the ADK ionization model for hydrogen use for the simulation of LPWA driven by a CO2 laser [17, 18].

The new multi-level ionization model enabled simulation studies of ionization injection of electrons into wakes created by a $CO_{2}$ laser pulse in high-atomic-number materials such as krypton. After the plasma wakes are excited by the $CO_{2}$ laser, a separate high-intensity near-infrared laser with a low normalized vector potential will inject electrons through ionization inside the wakes. Detailed simulation study of the interaction of both lasers with krypton and ionization injection of electrons will be published in a forthcoming paper. Here we illustrate the use of the multi-level ionization algorithm by the simulation of a single 4J, 2 ps $CO_{2}$ laser pulse creating plasma wakes in krypton gas with neutral density of $2\times 10^{17}$ 1/cc. Figure 3 shows a distribution of electrons on a 2D plane along the laser propagation and polarization directions.

5 Conclusions and future work

In this paper, we developed an efficient multi-level laser tunneling ionization algorithm in high-atomic-number gases. The algorithm features several improvements compared to the previously published methods based on numerical solutions to the system of ODE’s describing ionization dynamics via fractional time stepping. By using analytical solutions to the system of differential equations describing ionization evolution, our method effectively resolves the multiscale nature of ionization processes which occur at time scales significantly different compared to the time step of the main code. The method also eliminates the need to store all ionization states in every cell of the numerical mesh by using a locally-reduced system of equations, thus significantly reducing the memory requirement. It also accounts for electron quantum numbers and their projections while computing ionization probabilities.

The multi-level ionization algorithm was implemented and verified in the code SPACE. The algorithm developed in this paper extends a set of atomic physics algorithms previously implemented in SPACE that resolved single-level ionization, electron-ion recombination, and electron attachment to dopants in dense neutral gases. SPACE has been extensively used in a variety of applications in the area of plasma physics and accelerator design, in particular for simulations of coherent electron cooling of relativistic ion beams [16] and the laser-plasma wakefield acceleration [17, 18]. The multi-level ionization module is currently being applied to numerical studies of the ionization injection and acceleration of electrons into plasma wakes within the Brookhaven National Laboratory program on $CO_{2}$ laser driven plasma wakefield accelerator.

Acknowledgments

This work was supported by Grant No. DE-SC0014043 funded by the U.S. Department of Energy, Office of Science, High Energy Physics.

References

[1] Tajima, T. and Dawson, J. M., Laser electron accelerator, Phys. Rev. Lett. 43, 267 (1979).
[2] E. Esarey, C.B. Schroeder, W.P. Leemans, Physics of laser-driven plasma-based electron accelerators. Rev. Mod. Phys. 81, 1229 (2009).
[3] A.J. Gonsalves et al., Petawatt laser guiding and electron beam acceleration to 8 GeV in a laser-heated capillary discharge waveguide, Phys. Rev. Lett. 122 (2019) 084801.
[4] D. L. Bruhwiler, D. A. Dimitrov, J. R. Gary, E. Esarey, W. Leemans, and R. E. Giacone, Particle-in-cell simulations of tunneling ionization effects in plasmabased accelerators, Phys. Plasmas 10 (2003) 2022–2030.
[5] Kumar, P., Yu, K., Zgadzaj, R., Amorim, L., Downer, M. C., Welch, J., Litvinenko, V. N., Vafaei-Najafabadi, N. and Samulyak R. Simulation study of co2 laser-plasma interactions and self-modulated wakefield acceleration. Phys. Plasmas 26, 083106 (2019).
[6] M. Chen, E. Cormier-Michel, C.G.R. Geddes, D.L. Bruhwiler, L.L. Yu, E. Esarey, C.B. Schroeder, W.P. Leemans, Numerical modeling of laser tunneling ionization in explicit particle-in-cell codes, Journal of Computational Physics, 236 (2013), 220–228.
[7] C.K. Birdsall, A.B. Langdon, V. Vehedi, J.P. Verboncoeur, Plasma Physics via Computer Simulations, Adam Hilger, Bristol, 1991.
[8] R. Hockney, J. Eastwood, Computer Simulation using Particles, Taylor & Francis Group, New York, 1988.
[9] L.-L. Yu, E. Esarey, C.B. Schroeder, J.-L. Vay, C. Benedetti, C.G.R. Geddes, M. Chen, W.P. Leemans, Two-color laser-ionization injection, Phys. Rev. Lett., 112 (2014) 125001.
[10] C.B. Schroeder, J.-L. Vay, E. Esarey, S.S. Bulanov, C. Benedetti, L.-L. Yu, M. Chen, C.G¿ Geddes, W. Leemans, Thermal emittance from ionization-induced trapping in plasma accelerators, PRAB, 17 (2014) 101301.
[11] 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, Journal of Applied Physics, 99 (2006) 056109.
[12] M. V. Ammosov, P. A. Golovinsky, I. Yu. Kiyan, V. P. Krainov, Tunneling ionization of atoms and atomic ions in an intense laser field with a nonhomogeneous space–time distribution Sov. Phys. JETP 64 (1986) 1191
[13] K. Yu, P. Kumar, S. Yuang, A. Cheng, R. Samuyak, SPACE: 3D Parallel Solvers for Vlasov-Maxwell and Vlasov-Poisson Equations for Relativistic Plasmas with Atomic Transformations, submitted to Comp. Phys. Comm., 2022, arXiv:2111.04139.
[14] K. Yu, R. Samulyak, K. Yonehara, B. Freemire, Simulation of beam-induced plasma in gas-filled rf cavities, Physical Review Accelerators and Beams 20 (3) (2017) 032002.
[15] K. Yu, R. Samulyak, K. Yonehara, B. Freemire, Simulation of plasma loading of high-pressure RF cavities, JINST 13 (2018) P01008.
[16] G. W. J. Ma, R. Samulyak, V. Litvinenko, X. Wang, K. Yu, Simulation studies of modulator for coherent electron cooling, Physical Review Accelerators and Beams 21 (2018) 111001.
[17] P. Kumar, K. Yu, R. Zgadzaj, L. Amorim, M. C. Downer, J. Welch, V. N. Litvinenko, N. Vafaei-Najafabadi, R. Samulyak, Simulation study of co2 laser-plasma interactions and self-modulated wakefield acceleration, Physics of Plasmas 26 (8) (2019) 083106.
[18] P. Kumar, K. Yu, R. Zgadzaj, M. C. Downer, I. Petrushina, R. Samulyak, V. N. Litvinenko, N. Vafaei-Najafabadi, Evolution of the self-injection process in long wavelength infrared laser driven lwfa, Physics of Plasmas 28 (1) (2021) 013102.