Alternative gauge for the description of the light-matter interaction in a relativistic framework

TL;DR

This paper introduces a generalized velocity gauge for relativistic light-matter interaction that improves numerical convergence and multipole expansion accuracy, demonstrated through simulations of hydrogen under intense ultrashort laser pulses.

Contribution

A new velocity gauge formulation for relativistic laser-matter interaction that enhances numerical and multipole expansion convergence compared to traditional minimal coupling.

Findings

Superior convergence in numerical solutions of the Dirac equation.

Improved accuracy in multipole expansion of laser fields.

Effective modeling of hydrogen dynamics under intense ultrashort pulses.

Abstract

We present a generalized velocity gauge form of the relativistic laser-matter interaction. In comparison with the (equivalent) regular minimal coupling description, this new form of the light-matter interaction results in superior convergence properties for the numerical solution of the time-dependent Dirac equation. This applies both to the numerical treatment and, more importantly, to the multipole expansion of the laser field. The advantages of the alternative gauge is demonstrated in hydrogen by studies of the dynamics following the impact of superintense laser pulses of extreme ultraviolet wavelengths and sub-femtosecond duration.

Equations

$$i\hslash \frac{\mathrm{d}}{\mathrm{d}t}{\Psi }_{\mathrm{NR}}=H_{\mathrm{NR}}(t){\Psi }_{\mathrm{NR}},$$

$$H_{\mathrm{NR}}(t)=\left[\frac{{\mathbf{p}}^2}{2m}-e\phi +\frac{e}{m}\mathbf{p}\cdot \mathbf{A}+\frac{e^2{A}^2}{2m}\right].$$

$$\mathbf{A}(\eta )=A(\eta )\hat{\mathbf{z}}=\frac{E_0}{\omega }f(\eta )\sin (\eta +\varphi )\hat{\mathbf{z}},$$

$$f(\eta)=\left\{\begin{array}[]{lc}\sin^{2}\left(\frac{\pi\eta}{\omega T}\right),&0<\eta<\omega T\\ 0,&\text{otherwise}\end{array}\right..$$

$$\displaystyle A(\eta)\approx\sum_{n=0}^{n_{\mathrm{trunc}}}\frac{1}{n!}\frac{d^{n}A(\eta)}{d\eta^{n}}\bigg{\rvert}_{\eta=\omega t}\left(-\frac{\omega x}{c}\right)^{n}.$$

$$H(t)=c\mathbf{\alpha }\cdot \left[\mathbf{p}+e\mathbf{A}(\eta )\right]-e\phi (r)\text{mathbbm}{1}_4+m{c}^2\beta ,$$

$$\boldsymbol{\alpha}=\left(\begin{array}[]{cc}0&\boldsymbol{\sigma}\\ \boldsymbol{\sigma}&0\end{array}\right).$$

$$\beta=\left(\begin{array}[]{cc}\mathbbm{1}_{2}&0\\ 0&-\mathbbm{1}_{2}\end{array}\right).$$

$$i\hslash \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\Psi }=H(t)\tilde{\Psi },$$

$$\tilde{\Psi}({\bf r},t)=\left(\begin{array}[]{c}\tilde{\Psi}_{F}({\bf r},t)\\ \tilde{\Psi}_{G}({\bf r},t)\end{array}\right),$$

$$\tilde{\Psi }(\mathbf{r},t)=\Psi (\mathbf{r},t)e^{-im{c}^{2}t}.$$

$$i\hslash \frac{\mathrm{d}}{\mathrm{d}t}\Psi (\mathbf{r},t)=\left(H(t)-m{c}^2\right)\Psi (\mathbf{r},t).$$

$$-e\phi {\Psi }_F+c\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right){\Psi }_G=i\hslash \frac{\mathrm{d}{\Psi }_F}{\mathrm{d}t}$$

$$c\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right){\Psi }_F+\left(-e\phi -2m{c}^2\right){\Psi }_G=i\hslash \frac{\mathrm{d}{\Psi }_G}{\mathrm{d}t}.$$

$${\Psi }_G\approx \frac{1}{2mc}\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right){\Psi }_F.$$

$$-e\phi {\Psi }_F+\frac{1}{2m}{\left(\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right)\right)}^2{\Psi }_F=i\hslash \frac{\mathrm{d}{\Psi }_F}{\mathrm{d}t},$$

$$\left(\frac{{\mathbf{p}}^2}{2m}-e\phi +\frac{e}{m}\mathbf{p}\cdot \mathbf{A}+\frac{e^2{A}^2}{2m}+\frac{e\hslash }{2m}\mathbf{\sigma }\cdot \mathbf{B}\right){\Psi }_F=i\hslash \frac{\mathrm{d}{\Psi }_F}{\mathrm{d}t},$$

$$\mathbf{A}\to \mathbf{A}+\nabla \zeta$$

$$\phi \to \phi -\frac{\partial \zeta }{\partial t},$$

$$H=c\mathbf{\alpha }\cdot \left[\mathbf{p}+e\mathbf{A}(\eta )+e\nabla \zeta \right]+\left[e\frac{\partial \zeta }{\partial t}-e\phi (r)\right]\text{mathbbm}{1}_4+m{c}^2\beta .$$

$$p_k(\eta )=\frac{mc}{2}{\left(\frac{eA(\eta )}{mc}\right)}^2,$$

$$\zeta \left(\eta \right)$$

$$e\nabla \zeta =-\frac{m{c}^2}{\omega }\hat{\mathbf{x}}\frac{\partial \eta }{\partial x}\frac{d}{d\eta }{\int }_{-\infty }^{\eta }d{\eta }^{\prime }\frac{1}{2}{\left(\frac{eA({\eta }^{\prime })}{mc}\right)}^2\aleph \left({\eta }^{\prime }\right)=$$

$$+\hat{\mathbf{x}}k\frac{m{c}^2}{\omega }\frac{1}{2}{\left(\frac{eA(\eta )}{mc}\right)}^2\aleph \left(\eta \right)=\hat{\mathbf{x}}mc\frac{1}{2}{\left(\frac{eA(\eta )}{mc}\right)}^2\aleph \left(\eta \right),$$

$$e\frac{\partial \zeta }{\partial t}=-\frac{m{c}^2}{2}{\left(\frac{eA(\eta )}{mc}\right)}^2\aleph \left(\eta \right),$$

$$H_{\mathrm{PG}}=c\mathbf{\alpha }\cdot \left[\mathbf{p}+e\mathbf{A}(\eta )\right]-e\phi (r)\text{mathbbm}{1}_4+m{c}^2\beta$$

$$+\frac{e^2{A}^2(\eta )}{2m}\aleph \left(\eta \right)\left({\alpha }_x-\text{mathbbm}{1}_4\right),$$

$$-\left(\phi +\frac{e^2{A}^2}{2m}\aleph \right){\Psi }_F+$$

$$\left(c\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right)+{\sigma }_x\frac{e^2{A}^2}{2m}\aleph \right){\Psi }_G=i\hslash \frac{\mathrm{d}{\Psi }_F}{\mathrm{d}t}$$

$$\left(c\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right)+{\sigma }_x\frac{e^2{A}^2}{2m}\aleph \right){\Psi }_F$$

$$-\left(e\phi +2m{c}^2+\frac{e^2{A}^2}{2m}\aleph \right){\Psi }_G=i\hslash \frac{\mathrm{d}{\Psi }_G}{\mathrm{d}t}.$$

$${\Psi }_G\approx \frac{1}{2mc}\left(\mathbf{\sigma }\cdot \left(\mathbf{p}+e\mathbf{A}\right)+{\sigma }_x\frac{e^2{A}^2}{2mc}\aleph \right){\Psi }_F.$$

$$H_{\mathrm{PG}}^{\mathrm{NR}}=\frac{{\mathbf{p}}^2}{2m}+\frac{e}{m}\mathbf{p}\cdot \mathbf{A}+\frac{e^2{A}^2}{2m}-e\varphi +\frac{e\hslash }{2m}\mathbf{\sigma }\cdot \mathbf{B}$$

$$+\frac{1}{2mc}\left\{\frac{e^2{A}^2}{2m}\aleph ,p_x\right\}-\frac{e^2{A}^2}{2m}\left(\aleph -{\aleph }^2\frac{e^2{A}^2}{4m^2{c}^2}\right),$$

Full text

I Introduction

With high laser intensities, available already now or in the near future Moore et al. (1999); DiChiara et al. (2008); Yumoto et al. (2013); Yoneda et al. (2014), and the interesting possibilities then opening, discussed for example in Ref. Di Piazza et al. (2012), the description of the light-matter interaction in a relativistic framework is of growing importance. The ionization dynamics initiated with few-cycle laser pulses calls further for a time-dependent treatment. Several attempts have consequently been made to solve the time-dependent Dirac equation (TDDE), see, e.g., Refs. Selstø et al. (2009); Pindzola et al. (2010); Bauke et al. (2011); Vanne and Saenz (2012), but it has been proven a hard task to explore the truly relativistic region while simultaneously accounting for the spatial dependence of the electromagnetic field and the full dimensionality of the problem. Recently, however, a numerical study was made where high orders of multipole interaction terms were successfully accounted for Kjellsson et al. (2017). Field intensities up to the strength where electrons are expected to reach quiver velocities, $v_{\mathrm{quiv}}\approx eE_{0}/m\omega$, of around 20% of the speed of light were treated, and emerging relativistic effects could be detected. Still, the study also underlined some severe problems appearing when one is tackling the TDDE, concerning in particular the inclusion of magnetic effects.

When electrons are driven to high velocities by laser fields the magnetic part of the electromagnetic field inevitably becomes increasingly important. A qualitatively new effect emerging is then the force imposed on the particle in the propagation direction of the light. Simulations in the low- or medium-intensity regimes are usually made within the dipole approximation, where the spatial dependence of the vector potential of the pulse is neglected completely. Since this approximation implies a neglect of all magnetic effects it is rather pointless in the high intensity regime Reiss (2000). To understand the importance of different types of effects beyond the dipole approximation it is illustrative to look at the studies within this regime that have been done with the non-relativistic time-dependent Schrödinger equation (TDSE). In that case the spatial dependence of the vector potential may conveniently be treated through a Taylor expansion Førre and Simonsen (2014); Simonsen and Førre (2015) and the lowest order contribution has been shown to dominate the dynamics beyond the dipole completely – at least up to intensities that drive the electron to velocities just above ten percent of the speed of light Førre and Simonsen (2014). Surprisingly enough, as shown in Ref. Kjellsson et al. (2017), when the same approach is used with the TDDE, the lowest order spatial contribution from the Taylor expansion gives results that deviates significantly from the non-relativistic results already at modest intensities, far below the relativistic regime. This can be corrected by adding the next term in the expansion, but when the intensity is increased further the situation is repeated and one is forced to include also the following term and so on. This behavior can be analyzed and understood in the non-relativistic limit, as shown in Ref. Kjellsson et al. (2017) and also discussed in Sec. II.1 below. The problem stems from contributions that are known to cancel (approximately), but which enter in different formal orders with respect to the Taylor expansion when it is applied to the Dirac equation. Wherever the expansion is truncated, there will be unbalanced contributions which at some intensity will play a significant role. This imbalance is inherent to the four-component Dirac equation, and if the TDDE is to be solved for strong relativistic pulses the Taylor expansion approach in the regular minimal coupling Hamiltonian quickly leads to an intractable problem.

Recently a generalized velocity gauge form of the non-relativistic light-matter interaction was presented Førre and Simonsen (2016); Simonsen and Førre (2016). Within this gauge the dipole contribution is given as in velocity gauge, while the so called diamagnetic term disappears and instead new terms appear. Of these, the leading order ones depend explicitly on the momentum in the direction of the laser propagation and the new gauge was consequently coined the propagation gauge. It was further shown that it is possible to use a series of gauge transformations to successively remove all field-dependent terms that do not depend explicitly on the momentum. The lowest order interaction within this gauge was further tested Førre and Simonsen (2016) and compared to simulations performed with the traditional minimal coupling Hamiltonian. Impressive numerical advantages were then demonstrated. One manifestation of this was the evolution of the momentum expectation value along the direction of propagation of the light. In the minimal coupling description it showed a strong oscillating behavior, but in the new gauge it was replaced by a smooth curve that could be sampled with much larger time steps. Moreover, for a wave function expanded in spherical harmonics, the ionized wave packet could be described with considerably fewer angular momenta. In the following we will show that a corresponding gauge choice for the TDDE is even more advantageous. It requires just a single gauge transformation, it takes a simpler form, and it is a promising candidate for studying strong relativistic multipole interactions.

The paper is structured as follows: In the next section we outline the theoretical framework. Brief details on the implementation are provided in Sec. III, while the results are presented and discussed in Sec. IV. Finally, we present our conclusions in Sec. V. Atomic units are used throughout the text unless explicitly stated otherwise.

II Theory

In a non-relativistic framework, the evolution of a wave packet representing a particle of mass $m$ and charge $-e$ in the scalar field $\varphi$ and vector potential ${\bf A}$ is governed by the TDSE,

[TABLE]

with the Hamiltonian

[TABLE]

Here the potential ${\bf A}$ has been taken to fulfill the Coulomb gauge condition, $\nabla\cdot{\bf A}=0$. Letting the electromagnetic pulse be defined in terms of the vector potential ${\bf A}$, and assuming the field to be linearly polarized along the $z$ axis and propagating along the $x$ axis, the pulse may be written

[TABLE]

where $\eta=\omega t-kx$ and $k=\omega/c$. The envelope function is chosen to be sine squared:

[TABLE]

In the dipole approximation, when the spatial dependence of the vector potential is neglected, the $A^{2}$ term can be removed by a gauge transformation and consequently does not affect the dynamics. On the other hand, this diamagnetic term is known to give the leading contribution beyond the dipole approximation in the high-intensity limit. This has, e.g., been shown in Ref. Førre and Simonsen (2014), where the spatial dependence of the vector potential was examined with the help of a Taylor series expansion:

[TABLE]

In Ref. Førre and Simonsen (2014) $n_{\mathrm{trunc}}\leq 2$ was considered. A point worth noticing here is that the diamagnetic term is second order in $\mathbf{A}$, and thus an expansion of $A$ to a particular order in $(\omega x/c)^{n}$ does not imply an expansion of the Hamiltonian to the same order. The contributions to the Hamiltonian with, for example, $n=2$ come from the square of the $n=1$ term in Eq. (5) and from the cross term between the $n=0$ and $n=2$ terms. Furthermore, there should be considerable cancellations between these terms; their sum oscillates with twice the frequency of the light while each of them have a constant sign. Large cancellations were indeed found in Ref. Førre and Simonsen (2014), and it was concluded that it is decisive to include all terms that contribute to the Hamiltonian to a given order. When the corresponding time-dependent Dirac equation is solved using the minimal coupling Hamiltonian Kjellsson et al. (2017), the diamagnetic term, which is then only implicitly included, causes severe convergence problems in terms of the multipoles of the external field. This is connected to an effective blocking of the aforementioned cancellations as will be clear in the following.

Turning now to the TDDE, the first step will be to consider the minimal coupling Dirac Hamiltonian;

[TABLE]

and

[TABLE]

As usual, $\boldsymbol{\sigma}$ is given by the Pauli matrices, and

[TABLE]

We set out to solve the TDDE:

[TABLE]

where the four-component wave function $\tilde{\Psi}$ can be written as

[TABLE]

with $\tilde{\Psi}_{F}$ and $\tilde{\Psi}_{G}$ being two-component spinors, often called the large and small component, respectively. The potential $\varphi(r)$ is for the present purposes simply the Coulomb potential from a point nucleus, i.e., we neglect retardation effects in the electron-nucleus interaction and take the nuclear mass to be infinite.

II.1 The non-relativistic limit of the relativistic minimal coupling Hamiltonian

In order to understand the origin of the problems encountered with the TDDE expressed in terms of the orignal minimal coupling Hamiltonian, it is important to study its non-relativistic limit. Since we are aiming for a solution to the TDDE which describes a positive energy state, we may write

[TABLE]

Eq. (9) can then be rewritten as:

[TABLE]

Using the form of the wave function given in Eq. (10) we can write Eq. (12) as two coupled differential equations:

[TABLE]

If only the dominating terms on the second line is retained (i.e. assuming that the mass-energy term is large both compared to the potential energy, $2mc^{2}\gg e\varphi$, and to the time variation of the small component) it is possible to write the small component as:

[TABLE]

When inserting this into the first line of Eq. (13) we get

[TABLE]

and with some operator algebra, detailed in Ref. Kjellsson et al. (2017), this expression can be rewritten as

[TABLE]

where $\mathbf{B}=\nabla\times\mathbf{A}$. Apart from the spin-dependent term, the operators on the left-hand side are the same as those that appeared in Eq. (2). In particular, we note the diamagnetic contribution, which apparently is implicitly included in the TDDE through the coupling between the small and large component of the wave function. Hence the advice from Ref. Førre and Simonsen (2014) regarding the consistent inclusion of the $x^{n}$ terms from the Taylor expansion of $\mathbf{A}$ is not easy to follow for the TDDE. Since the Dirac equation is linear in the vector potential, truncation after a particular $n$ in Eq. (5) will result in an implicitly included diamagnetic term that contains the $x^{2n}$ contributions from the square of the $x^{n}$ term, but not the cross terms between higher- and lower-order terms that are also $x^{2n}$ contributions. As a consequence, the solution of the TDDE with only the lowest order spatial correction ($n=1$) to the vector potential generally gives meaningless results, as demonstrated in Ref. Kjellsson et al. (2017). Furthermore, the convergence of the dynamics with respect to $n_{\mathrm{trunc}}$ in Eq. (5) was shown to be very slow once the laser pulse parameters started to enter the relativistic regime. In the case where the electron was accelerated to a quiver velocity of $v_{\mathrm{quiv}}\sim 0.2c$, a fifth order expansion was necessary for converged results. It is natural to assume that an even higher order expansion would be necessary further into the relativistic regime, and with each additional term $x^{n}$ in Eq. (5), the computational demand quickly turns this into an intractable problem. It is thus highly relevant to instead seek an alternative route less prone to grow so complex when the dynamics become increasingly relativistic.

II.2 The relativistic propagation gauge

Since it is the actual way the diamagnetic contribution resurfaces in the TDDE that causes the convergence problems, it might be possible to find an alternative form where it is easier to balance the terms included in Eq. (5). We are for instance free to make a gauge transformation to change the scalar field and vector potential as;

[TABLE]

which will yield a transformed Hamiltonian:

[TABLE]

In Ref. Førre and Simonsen (2016) it was shown that by choosing a gauge that followed the classical electron momentum in the direction of the light propagation, $p_{k}$, the diamagnetic term in the Schrödinger equation could be removed and replaced by operators that showed superior convergence properties. In the relativistic case, as shown in Refs. Sarachik and Schappert (1970); Salamin and Faisal (1996), a free classical particle that is initially at rest will acquire the momentum

[TABLE]

when exposed to the electromagnetic field given by $\bf A(\eta)$. It is natural to assume that a suitable gauge can be found if $\zeta$ is defined using Eq. (19), but we start by defining it with an additional operator, $\aleph(\eta)$, that remains to be determined:

[TABLE]

The introduction of $\aleph$ is related to the distinction between the relativistic and the non-relativistic version of the gauge transformation leading to the propagation gauge formulation Førre and Simonsen (2016); Simonsen and Førre (2016). We will return to its specific form in the following.

With the vector potential polarized along the $z$ axis and the field propagating along the $x$ axis we obtain

[TABLE]

where $k=\omega/c$ has been used in the last step.

Equation (II.2) is a vector operator in the propagation direction of the field. Further, with

[TABLE]

we may now write down the propagation gauge Dirac Hamiltonian, $H_{\rm PG}$:

[TABLE]

where the first line is just the minimal coupling Dirac Hamiltonian from Eq. (6). The second line, on the other hand, displays one operator proportional to $\alpha_{x}$, the relativistic velocity operator in the direction of the propagation of the light, and one counter term. As we will see, this counter term cancels the implicit diamagnetic term contributed by the first line when the equation is examined in the non-relativistic limit.

II.3 The non-relativistic limit of the relativistic propagation gauge

Starting again from Eq. (12) but now adding the new terms from the second line in Eq. (II.2), we will instead of Eq. (13) find

[TABLE]

Following the derivation preceding Eq. (14), and assuming in addition that $2mc^{2}$ dominates also over $(e^{2}A^{2}/2m)\aleph$, we obtain a new approximate relation between the large and small component:

[TABLE]

Inserting this expression for $\Psi_{G}$ into Eq. (24) we find the propagation gauge Hamiltonian in the non-relativistic limit (cf. the expression for the minimal coupling Hamiltonian on the left-hand side of Eq. (16));

[TABLE]

where $\left\{a,b\right\}$ denotes an anticommutator. In addition to the original terms in Eq. (16), two new terms have appeared on the last line of Eq. (II.3). It is evident that if we put $\aleph=1$ the diamagnetic term is cancelled. However, another possibility is to require

[TABLE]

and thus get rid also of the term proportional to $A^{4}$. If Eq. (28) is regarded as the defining equation for $\aleph$, we can readily write down its expression as

[TABLE]

It is clear from Eq. (29) that its range of validity is restricted to the region where

[TABLE]

which is consistent with the approximation made to obtain Eq. (26). In this case we may also expand Eq. (29) and find

[TABLE]

which in fact is the series that was found in Refs. Førre and Simonsen (2016); Simonsen and Førre (2016), i.e.,

[TABLE]

With $\aleph$ defined this way we may write Eq. (II.3) as

[TABLE]

which, apart from the spin-dependent term, is identical to the propagation gauge Hamiltonian obtained directly from the TDSE in Refs. Førre and Simonsen (2016); Simonsen and Førre (2016).

II.4 The long-wavelength approximation

While a vector potential without spatial dependence does not introduce any magnetic interaction in the ordinary minimal coupling Dirac Hamiltonian, Eq. (6), a purely time-dependent ${\bf A}$ does provide an additional dynamical term in Eq. (II.2); the term proportional to $\alpha_{x}$. Again this is in agreement with the findings in Refs. Førre and Simonsen (2016); Simonsen and Førre (2016); in the propagation gauge the radiation pressure is accounted for through a velocity gauge-like operator acting along the propagation direction of the laser in spite of a spatially independent $\mathbf{A}$. The effective Dirac Hamiltonian in this long-wavelength approximation (LWA) is given by

[TABLE]

where the terms that lack spatial dependence altogether have been removed since they do not affect the dynamics. In Sec. IV we will show that the Hamiltonian Eq. (II.4) can account fully for the dominating effects beyond the dipole approximation for a wide range of electromagnetic pulses. It gives in fact excellent agreement with the much more demanding fifth order expansion of the Hamiltonian Eq. (6), as applied in Ref. Kjellsson et al. (2017).

In principle we are free to choose either $\aleph=1$, to follow the relativistic momentum in the direction of the propagation of the laser light, or as given in Eq. (29) to allow for a more straight forward comparison with the non-relativistic treatment. For high enough fields there will of course be differences for any non-exact implementation, as will be demonstrated in Sec. IV. Lastly, although the properties of the LWA-Hamiltonian Eq. (II.4) are very promising, we want to emphasize that a practical implementation of Eq. (II.2) is by no means restricted to the LWA-approximation. It is indeed possible to go further and introduce spatial dependence in $\mathbf{A}$, which should become important for large enough field strengths $E_{0}$ and/or in the limit of very high laser frequencies. For illustrative purposes we present a first order beyond the LWA Hamiltonian in the next section and later demonstrate that it gives negligible contributions for the laser pulses considered in this article, which accelerate the electron to quiver velocities $v_{\rm quiv}$ up to about $0.2c$.

II.5 Beyond the long-wavelength approximation

We introduce a spatial dependence in $\mathbf{A}$ by using $n_{\rm{trunc}}=1$ in Eq. (5). Then, using $\aleph=1$, the first order beyond the long-wavelength approximation (BYLWA1) Hamiltonian can be written as

[TABLE]

It may seem odd that the two terms proportional to $A^{2}$ in Eq. (II.2) have been expanded differently. However, according to the discussion in Sec. II.1, this is indeed the proper way of expanding the field-dependent terms as this minimizes the problem with inconsistent terms appearing in the corresponding non-relativistic Hamiltonian. Note also that terms lacking spatial dependence altogether have been removed from Eq. (II.5).

In the continuation, we will demonstrate that the Hamiltonian $H_{\rm LWA}$ in Eq. (II.4) provides practically all dynamics for fields penetrating into the relativistic region, simply by comparing its results with the corresponding results obtained with $H_{\rm BYLWA1}$ as defined above. For clarity we emphasize that $\aleph=1$ has been used in both Hamiltonians for a just comparison. Before presenting our results we briefly describe our numerical implementation.

III Implementation

We expand the wave function in eigenstates of the time-independent Hamiltonian, i.e., Eq. (6) without $\mathbf{A}$, giving

[TABLE]

with

[TABLE]

where

[TABLE]

Here $\kappa=l$ for $j=l-1/2$ and $\kappa=-(l+1)$ for $j=l+1/2$. $X_{\kappa,j,m}$ represents the spin-angular part which has the analytical form

[TABLE]

where ${Y}^{l_{\kappa}}_{m_{l}}(\theta,\phi)$ is a spherical harmonic and $\chi_{m_{s}}$ is an eigenspinor. The radial components $P_{n,\kappa}(r)$ and $Q_{n,\kappa}(r)$ are expanded in B-spline functions de Boor (1978);

[TABLE]

Just as in Ref. Kjellsson et al. (2017) we use B-spline functions of orders $k_{1}=7$ and $k_{2}=8$ for the small and large components, respectively. As has been shown by Froese Fischer and Zatsarinny Froese Fischer and Zatsarinny (2009), the use of different $k$ for the two components effectively removes the so called spurious states which are known to appear when the Dirac equation is solved within a finite basis set. We use a linear knot sequence with 500 B-spline functions for the large component and 501 for the small component up to $R_{\mathrm{max}}=150$ a.u.. To avoid unphysical reflections at the box boundary we have used a complex absorbing potential starting from $r=110$ a.u.. We include all spin-orbitals with angular momenta up to a certain $l_{\mathrm{max}}$ (as defined for the large component) and keep all the associated magnetic quantum numbers $m$, cf. Eqs. (38, 39). To speed up the propagation without compromising the results, high energy components have been filtered out from the basis.

In Sec. IV we present converged data for the energy distribution, the expectation value of the momentum operator along the pulse propagation direction and, finally, the total ionization yield from the hydrogen ground state exposed to a $15$ cycle, $95$ eV ($\omega=3.5$ a.u.) laser field of intensity $7\times 10^{19}$ W/cm2 ($E_{0}=45$ eV). For the two former quantities, converged data were obtained with $l_{\rm max}=30$ for the propagation gauge LWA, cf. Eq (II.4), corresponding to 1,902,594 states and about $2.32\times 10^{10}$ non-zero matrix elements. In order to arrive at the same result with the minimal coupling Hamiltonian, Eq. (6), it was necessary to include $l_{\rm max}=50$ for the case with $n_{\mathrm{trunc}}=5$ in Eq. (5). From now on we will refer to this level of approximation as fifth order beyond dipole (BYD5). The BYD5 simulation required 5,125,954 states and about $1.64\times 10^{12}$ non-zero matrix elements, i.e., roughly 70 times more than our converged propagation gauge simulations.

For the ionization yield, which was systematically investigated for both lower and higher values of $E_{0}$, convergence was always achieved with $l_{\rm max}=40$ for both $H_{\rm LWA}$, Eq. (II.4), and $H_{\rm BYLWA1}$, Eq. (II.5). For further details on the implementation, such as how interaction matrix elements are computed and which numerical schemes that are applied, readers are referred to Ref. Kjellsson et al. (2017).

IV Results

We will first present results for the following scenario: A hydrogen atom with the electron initially prepared in the ground state is exposed to a laser pulse, defined in Eqs. (3 - 4), with the parameters

[TABLE]

The pulse parameters are such that the electron’s quiver velocity is expected to reach about $v_{\mathrm{quiv}}\sim 0.1c$ and have been chosen primarily to demonstrate the convergence property of the relativistic propagation gauge – not to reveal relativistic effects per se. To show the convergence properties, the lowest order interaction in the propagation gauge, LWA, cf. Eq. (II.4), has been compared to the minimal coupling Hamiltonians ranging from BYD1 to BYD5, that is with $n_{\mathrm{trunc}}=1-5$ in Eq. (5). Figure 1 shows a comparison of the energy distribution of the ionized electron after interaction with the pulse. A somewhat typical convergence pattern for the minimal coupling simulations can be seen, where each successive interaction type pushes the distribution to either side of the fully converged result.

Figure 2 also shows the energy distribution but now only for the minimal coupling BYD5 result and the propagation gauge LWA result. The energy grid has been extended to include the three first ionization peaks and a logarithmic scale is used to better resolve the data. The coincidence between the LWA result in the propagation gauge, which only involves purely time-dependent fields, and the result using a fifth order Taylor expansion of the vector field within the minimal coupling formulation Eq. (6), is evident.

In the non-relativistic version of the propagation gauge Førre and Simonsen (2016) an important demonstration of its computational advantages was the smooth evolution of the expectation value of the momentum in the propagation direction of the laser field. The same behaviour is found also for the corresponding relativistic results in Fig. 3, where violent oscillations seen using the minimal coupling Hamiltonian are transformed into a smooth development in the propagation gauge. Interestingly, to the naked eye, $\langle p_{x}\rangle$ seems to be converged already at BYD3 while the probability distribution clearly requires at least BYD5, as seen in Fig. 1.

A comparison with non-relativistic simulations is also in order. As mentioned, the chosen pulse parameters result in an expected quiver velocity of $v_{\mathrm{quiv}}\sim 0.1c$ and only small relativistic corrections, if any, are expected. Figure 4 shows a comparison of the relativistic and non-relativistic probability distributions, as obtained by solving both the TDSE in the propagation gauge LWA, cf. Ref. Førre and Simonsen (2016), and the TDDE with the propagation gauge LWA Hamiltonain, Eq. (II.4), using an equivalent basis set in both cases. Indeed, there are no relativistic effects displayed in Fig. 4, and one only expects these to appear at even higher intensities.

In order to search for possible relativistic effects we now systematically increase the field strength up to about $E_{0}=100$ a.u., corresponding to $I\sim 3.5\times 10^{20}\mbox{ W/cm}^{2}$ and $v_{\mathrm{quiv}}\sim 0.2c$. Figure 5 shows the resulting ionization yield as a function of the electric field strength. The minimal coupling results from Ref. Kjellsson et al. (2017) (TDDE BYD5), the relativistic propagation gauge results obtained both within the LWA (Eq. (II.4)) and beyond (Eq. (II.5)), as well as the corresponding TDSE result, are shown for comparison. Again, the relativistic corrections seem to be very small. Nevertheless, a tiny relativistic shift manifested as a decrease in the ionization yield, is displayed as the quiver velocity approaches $v_{\rm quiv}\sim 0.2c$. Furthermore, Fig. 5 shows that the favorable behavior of the LWA propagation gauge Hamiltonian persists over a wide range of intensities, up to the onset of the relativistic regime.

Finally, going even further into the relativistic regime the important question of how to properly incorporate the full spatial dependence of the field in the propagation gauge Dirac Hamiltonian needs to be addressed. In an exact calculation, the choice of $\aleph$, be it simply $\aleph=1$ or as in Eq. (29) or, equivalently, Eq. (32), does of course not matter. However, with a truncated representation of the field, cf. Eq. (5), this choice may be of crucial importance – and increasingly so for increasing field strengths. This dependence has been investigated and the first results are shown in Fig. 6. As in Fig. 5, ionization probabilities are shown for both the TDSE and the TDDE, but here for field strengths $E_{0}$ in the interval $80-105$ a.u. and only within the LWA. In going from $\aleph$ as defined in Eq. (29) to simply choosing $\aleph=1$, i.e., truncation at the first term in Eq. (32), a small shift downwards is introduced – consistently both in the relativistic and non-relativistic treatments, respectively. Based on the present results, it is still unclear which of the choices for $\aleph$ represent the best approximation, nor where the LWA approximation breaks down. However, the correspondence between the gauge transformation and the relativistic momentum in the propagation direction for a corresponding free, classical electron moving in the field, cf. Eq. (19), suggests that simply $\aleph=1$ should be the best choice. Although we here leave these open questions for future research, simply due to the computational complexity of the problem, it should be noted that they could all be studied within the current computational framework.

V Conclusion

We have presented a generalized velocity gauge form of the relativistic light-matter interaction and demonstrated its superior convergence properties compared to the regular minimum coupling Hamiltonian. As in the non-relativistic case, the alternative relativistic gauge relaxes the requirement on the maximum angular momentum needed during the time propagation. However, the major advantage goes even beyond that. While the usual minimal coupling formulation is numerically tough for high-intensity fields treated in a non-relativistic framework, it constitutes an intractable problem in the relativistic case due to inherent imbalance in the Dirac equation. The propagation gauge to a large extent removes this imbalance and opens up for calculations on atoms subjected to electromagnetic pulses in the truly relativistic regime.

Acknowledgments

Simulations have been performed partly on resources provided by the Swedish National Infrastructure for Computing (SNIC) at TRIOLITH and KEBNEKAISE and partly on resources provided by the Norwegian Metacenter for Computational Science (UNINETT Sigma2) at HEXAGON (Account number NN9417K). Financial support by the Swedish Research Council (VR), Grant No. 2016-03789, is gratefully acknowledged. We also acknowledge support for our collaboration through the Nordic Institute for Theoretical Physics (Nordita) and from the research group Mathematical Modelling at Oslo and Akershus University College of Applied Sciences.