Relativistic ionization dynamics for a hydrogen atom exposed to super-intense XUV laser pulses

0.22 MB·

by  Tor Kjellsson

Relativisticionizationdynamics fora hydrogen atom exposedto super-intense XUV laserpulses
Tor Kjellsson,1 Sølve Selstø,2 and Eva Lindroth1
1DepartmentofPhysics,StockholmUniversity, AlbaNovaUniversityCenter,SE-10691Stockholm,Sweden
2FacultyofTechnology,ArtandDesign,OsloandAkershusUniversityCollegeofAppliedSciences,NO-0130Oslo,Norway
Wepresent atheoretical study oftheionization dynamics of ahydrogen atomexposed toattosecond laser
pulsesintheextremeultravioletregionatveryhighintensities.Thepulsesaresuchthattheelectronisexpected
toreachrelativisticvelocities,thusnecessitatingafullyrelativistictreatment.WesolvethetimedependentDirac
equationandcompareitspredictionswiththoseofthecorrespondingnon-relativisticSchro¨dingerequation.We
7
findthatastheelectronisexpectedtoreachabout20%ofthespeedoflight,relativisticcorrectionsintroduces
1
afiniteyetsmalldecreaseintheprobabilityofionizingtheatom.
0
2

n I. INTRODUCTION full field. To handle the computational load, which is quite
a heavy even for hydrogen, highly optimized parallel applica-
J
There are a number of infrastructure projects worldwide tionshavebeendeveloped.
5
With these we investigate to what extent the ionization
2 that strive for higher laser intensities, see e.g., the review
probabilitiespredictedbytheTDDEdifferfromthoseofthe
in Ref. [1]. For lasers operating with long wavelengths,
TDSE,inaregimewhererelativisticcorrectionsaretobeex-
] the intensities have already reached the regime where mag-
h pected.
netic interactions play a crucial role and ionized electrons
p
Ofcourse,relativisticeffectsarisewhentheelectronisac-
move with relativistic velocities, see e.g., Refs. [2, 3]. With
-
m the latest generation of free electron lasers (XFEL, SACLA, celeratedtovelocitiescomparabletothespeedoflight. This
maycomeaboutintwoways;forhighlychargednucleihigh
o LCLS), an unprecedented brilliance in the extreme ultravi-
velocities may be induced by the Coulomb potential alone.
t olet (XUV) region and beyond is reached, and new tech-
a Alternativelyastrongexternalelectromagneticfieldcandrive
niques[4]tofocusthebeam,aswellaspreliminaryresults[5],
.
s promise intensities also in this wavelength region exceeding electronstowardsrelativisticspeeds. Intheformercase,rela-
ic 1020 W/cm2. The treatment of light-matter interaction in a tivisticcorrectionstotheenergystructuredoofcourseinflu-
s ence the ionizationdynamics, e.g., by modifyingthe ioniza-
relativisticframeworkisthusatimelyissue.Tothisend,some
y
tionpotential[7,9,10]. However,inthepresentworkwewill
technicalobstaclesmust be overcome– obstacles specific to
h
restrictourselvesto hydrogenandinvestigatecases in which
p the relativistic time-dependentDirac equation (TDDE). One
theexternalfieldaloneisstrongenoughtopotentiallyinduce
[ issueishowtodealwiththenegativeenergypartofthespec-
relativistic dynamics. As a “measure of relativity” we may
trumof the Dirac Hamiltonian. Specifically, the stiffnessin-
1 takethemaximumquivervelocityofaclassicalfreeelectron
ducedbythehugeenergydifferencebetweenthepositiveand
v
exposedto a homogenouselectric field of strengthE oscil-
9 negative part of the spectrum may cause severe problems in 0
4 resolving the dynamics. An even more challenging issue is, latingwithfrequencyω:
3 asitturnsout,theconsistentinclusionofhigherordermulti-
eE
7 polesofthefullelectromagneticfield. v = 0 . (1)
quiv
0 mω
InRef.[6]theTDDEforhydrogen-likesystemsexposedto
.
1 strongattosecondlaserpulseswassolvednumerically.Unfor- Asv becomescomparabletothespeedoflightc,relativis-
0 tunately,computationalconstraintsdidnotallowforcalcula- tic efqfueicvtsare expected. In orderto actually see such effects
7
tionspenetratingintothe relativistic regimeforhydrogen. It inthe ionizationprobability,the laserspecificationsmust, of
1
: wasfound,however,thatevenbelowrelativisticvelocities,the course, be such that saturation is avoided even in this limit.
v inclusionofthenegativeenergypartofthespectrumiscrucial Thus, the calculationswill involvephotonenergieswell into
Xi in orderto accountfor dynamicsbeyondthe dipoleapproxi- theXUVregion. Varioustechniqueshavebeenappliedinor-
mation. Moreover,in Ref. [7] itwas shownthatevenwithin dertosolveTDDEnumerically,e.g.,splitoperatormethods–
r
a the dipole approximation, negative energy states are crucial combinedwith Fouriertransforms[10, 11] or the method of
forionizationprocessesinvolvingmorethanonephoton. characteristics[12,13],theclosecouplingmethod[9,14]and
More recently, it was demonstrated in Ref. [8] that for a Krylov methods [15]. In the present work TDDE is solved
laser pulse with photon energy in the XUV-region and field within a spectralbasis as in Refs. [6, 7], i.e., the state is ex-
strengthsbelowtherelativisticregime,higherordermultipole pandedinasetofeigenstatesoftheHamiltonianwithoutany
effectsarewellaccountedforbyusingthesocalledenvelope externalelectromagneticfieldpresent. Thetime propagation
approximation,whichdoesnotcontainthespatialdependence isperformedusingaloworderMagnusexpansion[16]while
of the electromagnetic field in full. It was further demon- theactualmatrixexponentiationisapproximatedbyaKrylov
strated that within this approximation, the solutions of the subspaceapproachasinRef.[15].
TDDE and the non-relativistic time dependent Schro¨dinger This paper is structured as follows: The next section out-
equation(TDSE)wereinagreement. linesthetheoreticalframework,anddetailsontheimplemen-
Inthisworkwegofurtherandsolvetheequationsofmotion tation are provided in Sec. III. Our results and findings are
intherelativisticregimeincludingmultipoleeffectsfromthe presentedanddiscussedinSec.IV,whileourconclusionsare
2
drawn in Sec. V. Atomic units are used throughout the text
unlessexplicitlystatedotherwise.
H (t)=cα A(x,t) . (9)
I z
Asthistermdependsonbothtimeandspace,tandx,adi-
II. THEORY
rectcalculationofthex-dependentcouplingsinducedbythis
interactionwouldhavetobeperformedateachandeverytime
Our starting point is the time-dependent Dirac equation
stepinordertorepresentH(t)numerically.Thiscumbersome
(TDDE):
feature may be removed by writing the vector potential as a
d sumoftermswithapurelytime-dependentandaspatiallyde-
i~ Ψ=H(t)Ψ, (2) pendentpart,
dt
withtheHamiltonian ntrunc
A(η) c T (t)X (x) . (10)
n n n
H(t)=cα [p+eA(η)]+V(r) 4+mc2β ≈ nX=0
· 1
=H0+ecα A . (3) Such separationsmay be achievedby, e.g., a Fourier expan-
·
sion in η or a Taylor expansion around η = t. In Ref. [6]
Fortherepresentationofα,thePaulimatricesareused,
both these approaches were followed. In the Fourier imple-
0 σ mentation,thenumberoftermswasminimizedintwoways:
α=(cid:18)σ 0 (cid:19) , (4) First, by taking A to have the pulse length T as period and,
second,byneglectingthespatialdependenceoftheenvelope
and f(η),c.f.,Eq.(8).Bothoftheseapproacheshavesevereshort-
comings.Theformerobviouslyintroducesanerroneousperi-
β = 12 0 . (5) odicity,while for the latter ithasbeen shownthatin general
(cid:18) 0 2 (cid:19) it is the spatial dependence of the envelope, not the carrier,
−1
that provides the dominant correction to the dipole approxi-
Thefour-componentwavefunctioncanbewrittenas
mation[8]. InviewofthisweresorttoaTaylorexpansionin
thepresentwork:
Ψ (r,t)
Ψ(r,t)= F , (6)
(cid:18)ΨG(r,t) (cid:19) ntrunc 1 x n ntrunc
A(η) A(n)(t) = a (t)xn (11)
n
whereΨF andΨG aretwo-componentspinors. Thepotential ≈ nX=0 n! (cid:16)−c(cid:17) nX=0
V(r)issimplytheCoulombpotentialofapointnucleus,i.e.,
weneglectretardationeffectsintheelectron-nucleusinterac- withoutneglectingspatialdependenceinneithertheenvelope
tion and take the nuclear mass to be infinite, thus allowing northecarrier.
forseparationbetweentheelectronicandthenucleardegrees The electric dipole approximation, which is generally not
of freedom. The mass energy term, i.e., mc2β, introduces applicableforthecasesofinteresthere,see, e.g.,thediscus-
a 2mc2 gap in the spectrum, dividing it into the aforemen- sion by Reiss [17], consists in substituting η with the time
tionednegativeandpositiveparts. Sincechangesinthepop- t, i.e., neglecting the spatial dependence of the laser pulse
ulationofnegativeenergystates isinterpretedastheappear- completely. In going beyond the dipole approximation, i.e.,
anceofpositrons(throughpair-creation),onemightarguethat assigning a value larger than zero to n in Eq. (11), it is
trunc
the negative spectrum should be excluded in simulating the however important that the spatial dependence is introduced
dynamics induced between an electron and an external field consistently in the equations of motion. This is not guaran-
withstrengthfarbelowthelimitofpairproduction. Thiswas teedbyjustchoosingn =1,whichmaybeexplainedby
trunc
disproved,however,inRefs.[6,7]andwewillbrieflyreturn firstlookingintothenon-relativisticinteraction.
tothisissuealsointhiswork.
WechoosetoworkinCoulombgauge, A=0,withthe
externalvectorpotentialAlinearlypolariz∇ed· alongthez-axis A. Non-relativisticinteraction
andpropagatingalongthex-axis;
Thetime-dependentSchro¨dingerequationisgivenby:
E
A(η)= 0f(η)sin(ωη+ϕ)zˆ , (7)
ω d
i~ Ψ =H (t)Ψ (12)
NR NR NR
dt
where η = t x/c. The envelope function is chosen to be
−
sine-squared; withtheHamiltonian
sin2 πη , 0<η <T p2 e e2A2
f(η)= T . (8) HNR(t)= +V(r)+ p A+ .
(cid:26)0, (cid:0) (cid:1) otherwise (cid:20)2m m · 2m (cid:21)
With the pulse being linearly polarized in the z-direction, InsolvingtheTDSE,Eq.(12),ithasbeenfoundthattheA2-
thetime-dependentpartoftheHamiltonianbecomes term,thesocalleddiamagneticterm, providespracticallyall
3
correctionstothedipoleapproximation[18,19].Moreover,in Specifically,theA2-termseeninthenon-relativisticHamil-
Ref.[19]itwasfoundsufficienttoincludeonlythefirstorder tonian of Eq. (12) reappears. Thus, we see that the A2-
correctionin the diamagnetic term. Here it is worth empha- termisimplicitlypresentintheDiracequationandentersthe
sizingthatinthecaseoftheSchro¨dingerequationitiscrucial equation for the large componentof the wave function, Ψ ,
F
thattheHamiltonian,thusnotthevectorpotentialitself,isex- throughthe small component, Ψ . Whenever we try to for-
G
panded consistently in x. Our resulting non-relativistic first mulatetheinteractionwith theelectromagneticfield through
orderHamiltonianisthen; operatorsthatareanti-diagonalwithrespecttothesmalland
largecomponent,itreappears. Thisimplicitoccurencecom-
H p2 +V(r)+ e p A(t) e2xA(t)A(1)(t) , (13) plicatesaconsistentinclusionofspatialeffectsintherelativis-
NR z
≈ 2m m −m c ticinteraction,aswillbedemonstratedinSec.IV.
where the purely time-dependent A(t)2-term has been re-
movedbythetrivialgaugetransformation:
C. Propagation
Ψ˜(r,t)=e~i R0t 2em2A(ωt′)2dt′Ψ(r,t). (14)
Bothintherelativisticandthenon-relativisticcase,thestate
In the non-relativistic regime the TDSE and the TDDE vectorΨ(t)ispropagatedbymeansofasecondorderMagnus
should,ofcourse,agree.Thisconditionwillnowproveuseful propagator,
in understanding the consistent incorporation of effects be-
yond the dipole approximation in the Dirac Hamiltonian by Ψ(t+τ)=exp[ iτH(t+τ/2)]Ψ(t)+ (τ3) . (19)
− O
studyingitsnon-relativisticlimit.
OneofthemajoradvantagesofaMagnus-typepropagatorfor
the Schro¨dinger equation is stated clearly in Ref. [21]: “In
B. Thenon-relativiticlimitofthelight-matterinteraction contrast to standard integrators, the error does not depend
onhighertimederivativesofthesolution,whichisingeneral
highlyoscillatory”. Dueto thestiffnessinherentin themass
While the time dependentSchro¨dingerequation, Eq.(12),
energyterm,i.e.,themc2β-termoftheHamiltonianinEq.(3),
hasbotha linearand a quadratictermin A, the time depen-
thisbecomesevenmoreadvantageousfortheDiracequation
dent Dirac equation, Eqs. (2,3), is only linear in the vector
than for the Schro¨dinger equation. In fact, the accuracy of
potential.Wecanstudythenon-relativisticlimitoftheTDDE
manytime-propagationschemes,suchasCrank-Nicolsonand
by using the form of the wave functiongiven in Eq. (6) and
Runge-Kutta, suffer greatly from the 2mc2 energy splitting.
rewriteEq.(2):
InRef.[22],e.g.,itisstatedthat“Themajordrawbackofthe
dΨ Diractreatmentisthetemporalstepsize∆t.~/Erequired,
VΨ +cσ (eA+p) Ψ =i~ F (15)
F · G dt which has to be significantly smaller than for Schro¨dinger
dΨ treatments, because of the large rest mass energy mc2 that
cσ (eA+p)ΨF + V 2mc2 ΨG =i~ G , (16) iscontainedintheparticle’stotalenergyE.”Indeed,several
· − dt
(cid:0) (cid:1) works dealing numerically with the TDDE apply time-steps
with ΨF and ΨG being the upper and lower component, re- of the order 10−5 a.u. or smaller, see, e.g., [9, 10, 13, 23].
spectively,c.f. Eq.(6). Forpositiveenergystates, Obviously,sucharestrictionrendersthedescriptionofalaser
pulseofarealisticdurationratherinfeasible.Thisproblemis,
ΨF c ΨG , however, circumventedby the application of Magnuspropa-
| |∼ | |
gators[21]. Itshouldbenotedthatextremelysmalltimesteps
soacomparisonwithΨ shouldbedictatedbyΨ .
NR F are,ofcourse,requiredwhenresolvingphenomenawhichre-
By assuming that the Coulomb potential is negligible in
allydotakeplaceatsuchshorttimescales, suchasZitterbe-
comparison with the mass energy term, V 2mc2, and
wegung[15].
≪
thatthetimevariationofthesmallcomponentΨ ismodest,
G Weemphasizethatinatime-dependentcontextitiscrucial
Eq.(16)yields
tokeepinmindthatatimedependentHamiltonianleadstoa
time-dependentdistinctionbetweenpositiveandnegativeen-
1
ΨG σ (eA+p)ΨF , (17) ergy states; the Dirac sea is not calm, so to speak. As was
≈ 2mc ·
discussedinRef.[6],negativeenergysolutions,asdefinedby
whichinsertedintoEq.(15)provides thetime-independentHamiltonian,H ,areessentialtocalcu-
0
late effects beyondthe dipole approximations, while the dy-
p2 +V + e p A+ e2A2 + e~ σ B Ψ =i~dΨF. namicnegativeenergystates,asdefinedbyH(t),shouldonly
(cid:20)2m m · 2m 2m · (cid:21) F dt comeintoplaywhenpair-productionstartstoplayarole. For
(18) fieldsfarawayfromthatlimitthepropagatorofEq.(19)may
Thatis,theTDSE,Eq.(12),isreproduced–withanadditional inprinciplebemodifiedto
term corresponding to the interaction between the spin and
the magnetic field. The same can be achieved via a Foldy- Ψ(t+τ)= (t+τ/2)exp[ iτH(t+τ/2)]Ψ(t)+ (t3) ,
P − O
Wouthuysentypeoftransformationonthewavefunction[20]. (20)
4
where (t) projects the state onto the time-dependent sub- energy index n can attain per spin-angular symmetry. The
P
space spanned by the positive spectrum of the Hamiltonian boundaryconditionsappliedonthecomponentsare:
H(t), i.e., the negative energy states are blocked. For the
fields used here the population of the negative energy states P (0) = P (R )=0 (26)
n,κ n,κ max
oftheHamiltonianH(t)isso tinythattheprojectedandthe Q (0) = Q (R )=0 . (27)
n,κ n,κ max
unprojectedpropagatorgivesthesameresults,asexplainedin
Sec.IV.
Weincludeallspin-orbitalswithangularmomentumupto
When,ultimately,thefieldsareincreasedevenfurtherand
l = 30(asdefinedforthelargecomponent)andkeepall
max
pair-productiondo startto playa role, this has, of course, to
the associated magnetic quantumnumbersm . To speed up
j
behandledthroughfieldtheory,wherethedistinctionbetween
the propagation without compromising the results, high en-
positive and negative energy states is inherent. Then transi-
ergycomponentshavebeenfiltered outleaving, forthistyp-
tionsintonegativeenergystates(annihilation)arecoupledto
ical choice of parameters a final number of 1902594 states
excitations out of them (pair-production)and the number of
in ourbasis. Similarly, the non-relativisticspectralbasishas
particlesisnolongerconserved.
eigenfunctionsoftheform
P (r)
III. IMPLEMENTATION Φ(r)= nl Yl(θ,φ) , (28)
r m
We expandour wave function in eigenstates of the unper- where, as in the relativistic case, the radial component is
turbedHamiltonianH , expanded in B-splines while the angular part is analytically
0
known. Weusethesamelinearknotsequenceasintherela-
Ψ(t)= cn,j,m,κ(t)ψn,j,m,κ(r) , (21) tivisticcasewith500B-splines,k = 7andRmax = 150a.u.
n,Xj,m,κ andtheboundaryconditions,
with
P (0)=P (R )=0 . (29)
n,l n,l max
F (r)
ψn,j,m,κ(r)=(cid:18)Gnn,,jj,,mm,,κκ(r)(cid:19) , (22) For convergence we needed to keep all orbital angular mo-
mentauptol =40,withallassociatedm -values. Justas
max l
where intherelativisticcasewefilteredouthighenergycomponents
that do not affect the dynamics, leaving a total of 827351
F (r) 1 P (r)X (Ω)
n,j,m,κ = n,κ κ,j,m . (23) statesinthenon-relativisticbasis.
(cid:18)Gn,j,m,κ(r)(cid:19) r (cid:18)iQn,κ(r)X-κ,j,m(Ω)(cid:19)
Here κ = l for j = l 1/2 and κ = (l + 1) for j =
l + 1/2, and X rep−resents the spin-a−ngular part which A. Computingthematrixelements
κ,j,m
hastheanalyticalform
WiththevectorfieldgivenbyEq.(11),thelight-matterin-
X = l ,m ;s,m j,m Ylκ(θ,φ)χ , (24) teraction,Eq.(9),givesrisetomatrixelementsbetweeneigen-
κ,j,m hκ l s| i ml ms
mXs,ml states to the time-independent Hamiltonian. Labelling two
such states k = nκjm and k˜ = n˜κ˜˜jm˜ , the matrix
where Ymlκl(θ,φ) is a spherical harmonic and χms is an elementconn|eicting|themisigiven|byiasum| ofterimsincluding
eigenspinor. TheradialcomponentsPn,κ(r)andQn,κ(r)are higherandhigherpowersofthespatialcoordinate:
expandedinB-splines[24];
ntrunc
Pn,κ(r)= aiBik1(r) , Qn,κ(r)= bjBjk2(r) . Hkk˜(t)= ai(t)hnκjm|αzxγ|n˜κ˜˜jm˜i. (30)
Xi Xj γX=0
(25)
In Ref. [25] it is demonstrated that specific choices of B- SinceweareworkingintheeigenbasisofH0, wemaycom-
splineordersk andk controltheoccurrenceofthesocalled putethesecouplingsusingangularmomentumtheorybyfirst
1 2
spuriousstates, which are knownto appearin the numerical expressing the operators in terms of spherical tensor opera-
spectrumafterdiscretizationoftheDiracHamiltonian.While tors[26]. Forthepowersofx, thesearesphericalharmonics
the choicek = k contaminatesthe spectrumwith suchin- Yλ;
1 2 µ
correctstates,thechoicesk =k 1arereportedtobestable
1 2
±
combinationsthatdonotproducethem. γ λ
In this work we use k = 7,k = 8. Converged re- xγ = cλµrγYλ , (31)
1 2 γ µ
sultswereobtainedusingalinearknotsequencewith500B- λX=0µX=−λ
splinesforthelargecomponentand501forthesmalloneup
to R = 150 a.u.. This gives a total of 1001 bound and while σ = σ1 is a componentof a rank-1 spherical tensor
max z 0
pseudocontinuum(bothpositiveandnegative)statesthatthe operator.ThecouplingsinEq.(30)arenowgivenbysumming
5
uptermsfactoredinaradialandspin-angularpart: where,justasinEq.(35),theWigner-Eckarttheoremhasbeen
applied. Note that this enables an efficient memory storage
nκjmα rγYλ n˜κ˜˜jm˜ =
h | z µ| i since all dependence on the projection numbers may be ac-
countedforbymultiplicationofsinglescalars. Thisfacthas
i P∗ (r)rγQ (r) κjmσ1Yλ -κ˜˜jm˜
Zrh nκ n˜κ˜ h | 0 µ| i− alsobeenusedforcache-efficientmemoryusageintheprop-
agationmethodused.
Q∗ (r)rγP (r) -κjmσ1Yλ κ˜˜jm˜ dr (32)
nκ n˜κ˜ h | 0 µ| ii
where -κ implies the spin-angular part of the small compo-
B. Propagation
nent. With the radial components expressed in B-splines,
theintegralsarecomputedtomachineaccuracyusingGauss-
Asthe useof a Magnuspropagator,Eq.(19), involvesex-
Legendrequadrature. Toobtainthespin-angularpart,theop-
ponentiatingatime-dependentHamiltonianmatrix,fulldiag-
eratorproductcan be expressedin a coupledtensor operator
onalization at each time step is called for. This would also
basis:
be necessary in order to make the distinction between (dy-
λ+1 namical) positive and negative energy states needed for the
σ1Yλ = 10;λµKQ σσσ1Yλ K , Q=µ .
0 µ h | i Q projection in Eq. (20). Incidentally, in the work of Ref. [6],
K=X|λ−1| (cid:8) (cid:9) imposingthisprojectionwasactuallynecessaryduetotheuse
ofthecomplexscalingmethod,whichwasthenakeytomini-
Withthischoicethespin-angularpartmaybecomputedas:
mizethebasissize. Sincenegativeenergystatesattainaposi-
λ+1 tiveimaginaryenergycomponentundercomplexscaling,they
hκjm|σ01Yµλ|κ˜˜jm˜i= (−1)λ−1−Q√2K+1 cannotbepropagatedforwardintimeandhastobeseparated
K=X|λ−1| fromthepositiveenergystates.
×(cid:18)01 µλ Kµ(cid:19)hκjm| σσσ1Yλ KQ |κ˜˜jm˜i . (33) cisNioenedmleastrsixtoosfasyi,zerep(e1a0t6e,d1d0i6a)goinsaelixztaretimoneloyfeaxpdeonusbivleeparned-
− (cid:8) (cid:9) completelyoutofthequestion.Anyefficientwayofdoingthe
TheWigner-Eckarttheoremcanbeappliedtothematrixele- exponentiationapproximativelyisthuswelcome.Tothisend,
mentofthecombinedoperator σσσ1Yλ K: Krylovsubspacemethodsturnouttobequiteuseful[15,27].
Q
(cid:8) (cid:9) Such methodsprovideaccurateapproximationsto the action
κjm σσσ1Yλ K κ˜˜jm˜ = oftheexponentialofanoperatoronaspecificvector. More-
h | Q | i over, their numericalimplementation can be made very effi-
(-1)j−m(cid:8) j (cid:9)K ˜j j σσσ1Yλ K ˜j (34) cient.AteachtimesteptheKrylovsubspaceofdimensionm,
(cid:18) m Q m˜(cid:19)×h || || i Km(t+τ/2),whichisspannedbythesetofstatesobtained
− (cid:8) (cid:9) byiterativelymultiplyingthestatewiththeHamiltonian,
withthereducedmatrixelementgivenby:
[H(t+τ/2)]kΨ(t), k =0,...,(m 1) , (38)
j σσσ1Yλ K ˜j = (2j+1)(2K+1)(2˜j+1) is constructed using the Arnoldialgorithm.−The exponential
h || || i q
(cid:8) (cid:9) l ˜l 1 isnowprojectedontoKm(t+τ/2)andexponentiatedwithin
×hs||σσσ1||s˜ihl||Yλ||˜lisj Js˜˜ Kλ . (35) othfisτs≈ubs1p0a−c3e.a.uT.ypiincaolulyr cmalc≈ulat5i0onws.heOnnucesincognavetrigmeencseteips
achieved,a backtransformationto the originalHilbertspace
With this scheme the couplings induced by α xγ, for γ = is performed to give Ψ(t + τ). The computationally heavy
z
0,1,...,n in Eq. (30) are readily computed to repre- partinthisapproachistherepeatedmatrix-vectorproductsin
trunc
sentthelight-matterinteractionterm,cf.Eq.(9),intheDirac Eq.(38). DespiteH(t+τ/2)beingquitesparse,thenumber
Hamiltonian, Eq. (3). For the Schro¨dinger Hamiltonian in ofnon-zeroelementsintherelativisticcaseisapproximately
Eq.(13)theneededcouplingsarep -couplings; 3.9 1011, whichisquiteachallengetohandle. Thisisdone
z ·
bystoringtheprojectionfactorsseparatefromtherestofthe
couplings, as discussed in Sec. IIIA, such that the memory
n l m p n l m =
a a a z b b b requirementfor our parameters is reduced by roughly a fac-
h | | i
( 1)la−ma la 1 lb n l p1 n l , (36) tor ∼ 100 while high computational throughput is achieved
− (cid:18) ma 0 mb(cid:19)×h a a|| || b bi byperformingallmultiplicationscorrespondingtotransitions
−
betweenstatessharingallotherquantumnumberssimultane-
andx-couplings;
ously.
In order to prevent reflection at the computational box
n l m xn l m =
a a a b b b
h | | i boundary, a complex absorbing potential (CAP) is added to
2π ( 1)la−ma la 1 lb la 1 lb H(t)duringthetime-propagation–forboththeSchro¨dinger
r 3 · − (cid:18)(cid:18) ma 1 mb(cid:19)−(cid:18) ma 1 mb(cid:19)(cid:19) andtheDiracequation,
− − −
×Z Pn∗ala(r)rPnblb(r)drhnala||Y1||nblbi , (37) VCAP =(cid:26)−η(r0−,r0)2, rr >rr0 . (39)
0
≤
6
TypicalCAP values are η = 0.05iand r = 110.0a.u., the andtheTDDEdoindeedagreebelowandaroundthisregime.
0
latterbeinglargeenoughtoensurethatanyfluxreachingthis However,theconvergencepatternswithrespecttothespatial
distance will really represent the ionization current. In that dependence in Eq. (11) of the calculations are actually very
case the absorptionof the flux beyondr doesnotaffectthe different. We will now study this convergence behavior in
0
ionization dynamics – provided that reflection is negligible. somedetail.
For the TDDE the complex absorbing potential is found to
havethesameeffectonallthecomponentsofthewavefunc-
tion. A. Therepresentationofthevectorfieldbeyondthedipole
approximation
IV. RESULTSANDDISCUSSION Aspreviouslymentioned,Ref.[19]providesstrongsupport
fortheclaimthatbelowtherelativisticregion,thefirstorder
The probatility of ionizing a hydrogen atom from the term in Eq. (11) alone providespractically all correctionsto
ground state, P , by exposing it to laser pulses of various the dipole prediction. At first glance there does not seem to
ion
intensitieshasbeenstudied.Thepulseischaracterizedbythe beanyapriorireasonwhythisconclusionshouldnotapplyto
followingparameters,cf. Eqs.(7,8): therelativistictreatmentaswell. Yet,asshowninFig.2,we
typicallyhavetoresorttoathirdorderexpansiontoreproduce
ω =3.5a.u., φ=0, T =2πNωc a.u. and Nc =15. (40) thenon-relativisticresultswhensolvingtheDiracequation.
P hasbeencalculatedwithpeakelectricfieldstrengths,E ,
ion 0
sohighthattheelectronquivervelocity,vquiv,correspondsto ω=3.5 a.u., NCycles=15
almost 20 % of c, cf. Eq. (1). But before we discuss such
0.7 TDSE dip
extremeconditionswewillconsidertheionizationprobability TDSE BYD
inthenon-relativisticregime. 0.6 BYD1 v≈0.1c
Figure 1 shows the convergedP calculated both within BYD2
thedipoleapproximationandbeyonidon(BYD)forE0 70a.u. eld0.5 BYD3
t(icvoerlryeslopwonfideinldgsttoreIngt∼hst1h.e7i·on1i0z2a0tioWn/pcrmo2b)a.bilFitoyrisc≤oinmcpreaaras-- tion yi0.4
a
ing monotonouslywith increasing field strength, but around z0.3
ni
E 10a.u. theso-calledstabilizationsetsin[28]. Foreven o
hi0gh≈er field strengths the ionization starts to increase again. I0.2
At E 30 a.u. the dipole approximation breaks down, a
0 ≈ 0.1
behaviourthathasbeendiscussed,e.g.,inRef.[29].
0.0
0 10 20 30 40 50 60
Peak electric field stength (a.u.)
ω=3.5 a.u., NCycles=15
0.8
TDSE dip
FIG.2: (coloronline). ComparisonbetweenTDDEresultswithin-
0.7 TDDE dip
creasingorderofxincludedintheexpansionofthevectorpotential.
TDSE BYD
0.6 Thenumberstrailingtheacronym“BYD”isthenumericalvalueof
TDDE BYD
n in Eq. (11). BYD1 starts deviating from TDSE already at
d trunc
el0.5 E ≈ 12 a.u. and from there on consistently overestimates P .
n yi BY0D2, on the other hand, agrees with TDSE up to E0 = 40 aio.nu.
tio0.4 where it starts to give a lower value for Pion. BYD3 then pushes
niza0.3 PionupandagreeswiththeTDSEresultsuptoE0≈60a.u..
o
I
0.2 TheunderlyingproblemistheimplicitinclusionoftheA2-
termintheDirac-equation,discussedinSec.IIB.Whenonly
0.1
v≈0.1c firstordercorrectionsto Aareincluded,someofthe second
0.0 ordercorrectionsintheimplicitA2-termwillstillbeincluded,
0 10 20 30 40 50 60 70
Peak electric field stength (a.u.) namely the (A(1))2x2 contribution, while the other con-
∼
tribution, namely the A(2)A(0)x2-term, require the inclu-
∼
FIG.1:(coloronline).ComparisonofTDSEandTDDEcalculations sionalsoofsecondordercorrectionstoA. FørreandSimon-
withinthedipoleapproximationandbeyond(BYD).Theverticalline sen [19] foundthat in the non-relativistic limit there are im-
indicates a maximum electric field strength E0 corresponding to a portantcancellationsbetweenthesex2-terms;itsneteffectis
maximumquivervelocityvquiv =0.1c. vanishinglysmallwhileseparatecontributionswouldshiftre-
sultsdramatically–andartificially.
The verticalline in Fig. 1 marksthe electric field strength It is thus natural to assume that these cancellations are
E that corresponds to a maximum electron quiver velocity a key-issue and that higher-order correction terms to A are
0
v = 0.1c. As expected, the predictions by the TDSE neededtoachievethemwhenimplementingasolutionofthe
quiv
v≈0.1c
7
Dirac equation [30]. Of course, inclusion of second order
correction terms in A, in turn, introduceseffective third and ω=3.5 a.u., NCycles=15
fourth order terms in an equally inconsistent manner and so 0.8 TDSE BYD
forth. However,sincethemagnitudeofthesecorrectionsde- BYD3
BYD4
crease with n , cf. Eq. (11), the problem should for a
trunc 0.7 BYD5
givenfieldstrengthdiminishwithincreasingorders.Thiscon- d
vergence behavior is observable in Fig. (2); TDDE BYD1, yiel
i.e., n = 1, severely overestimates the ionization yield n 0.6
trunc o
evenbeforethebreakdownofthedipoleapproximation,while ti
a
z
TDDEBYD2, i.e., ntrunc = 2, givesinitialagreementupto ni0.5 v≈0.1c v≈0.2c
o
E0 = 40a.u. butisthenseentounderestimatePion (aspre- I
dicted by the TDSE). Inclusion of an additional third order
0.4
terminEq.(11)increasesP ascomparedtoTDDEBYD2,
ion
givingcloseagreementwiththeTDSEuptoE =60a.u..
0
0.3
Beforewegoon,acommentisnecessary. Theconclusion 50 60 70 80 90 100
concerningtheinconsistentrepresentationoftheimplicitA2- Peak electric field stength (a.u.)
termmayseemtocontradicttheresultspresentedinRef.[8],
where the agreement was found between TDSE and TDDE FIG.3:(coloronline).ContinuationofFig.2withntrunc=3,4and
usingafirstorderexpansioninxofAutilizingthesocalled 5 for the TDDE. At the highest fieldstrengths the relativisticPion
starts to show a decreasing value compared to the non-relativistic
envelopeapproximation. Withthisapproximationthespatial
prediction.
dependenceinthecarrier,i.e.,sin(ωη+ϕ)inEq.(7),isdis-
garded,whichfortheTDDEinfactremovestheinconsistency
problemtoa largeextent: Theremainingtermin thederiva-
field,issmallerbyafactorv /c,i.e.,stillfarfromrelativis-
tive of the vector potential, A(1)(t), now only consists of a quiv
tic.
small contribution from the envelope, cf. Eq. (8), and thus
the problematic (A(1))2x2 contribution from the implicit
∼
A2 termin Eq.(18) is moreor lessinsignificant–atleast in
ω=3.5 a.u., NCycles=15
the non-relativistic regime. Within the envelope approxima-
0.02
Dip
tionitisthussufficienttouseonlythefirstorderterminthe
BYD5
Taylor-expansionofA, butitisnotwhenalsoitsfullspatial
dependenceisconsidered. 0.01
Inordertofindrelativisticeffectswe nowconsiderhigher
valuesforE0. ual
d 0.00
si
e
R
B. Relativisticeffects
−0.01
Figure3showsP forelectricfieldstrengthsuptoE =
ion 0
90a.u. (correspondingtoI 2.8 1020 W/cm2). Forhigher
∼ · −0.02
field strengths than those seen in Fig. 1, higher values of
50 60 70 80 90 100
ntrunc arenecessary,asshowninFig.3. Forthehighestfield Peak electric field stength (a.u.)
strengthsconsideredhere, relativistic effectsstart to surface.
Inordertohighlightthesewepresentthedifferencebetween FIG.4: (coloronline). ResidualplotforthefinalconvergedTDDE
therelativisticandthenon-relativisticpredictionsinFig.(4). calculations with respect to TDSE. Both within and beyond the
Both within and beyond the dipole approximation the rela- dipoleapproximationtherelativisticcorrectionsshowadecreaseof
tivisticP seemstodecreasesteadilycomparedtothecorre- Pionasthequivervelocityvquivapproaches0.2c.
ion
spondingnon-relativisticvalue. Albeitsmall,thediscrepancy
iswellwithintheaccuracyofourcalculations.
ItisclearfromFig.4thatrelativisticeffectsreducetheion-
It is interesting to note from Fig. 4 that the differencebe- izationprobability. Itwouldseem reasonableto assume that
tweentherelativisticandthenon-relativistictreatmentismore this reduction is related to the increased inertia of the elec-
orlessthesamewithinthedipoleapproximationasbeyondit. tron. Inordertotestthisassumptionwehaveperformednon-
Thisindicatesthat,althoughthevelocityoftheelectroninthe relativisticcalculationsinthedipoleapproximationinwhich
polarizationdirection,vquiv,cf.Eq.(1),inducedbytheelec- the electron mass has been substituted with the relativistic
tricfield,issubjecttorelativisticcorrections,themagneticin- massofaclassicalfreeelectron,
teractionisessentiallyunaffected.Thisisnotsurprising,how-
ever, consideringthe magnitude of the relativistic correction
hereandthefactthatthevelocityin thedirectionperpendic- m e
m where v(t)= A(t) . (41)
ulartothedirectionofpolarization,inducedbythemagnetic → 1 (v(t)/c)2 m
−
p
v≈0.2c
v≈0.1c
8
The difference between the corresponding ionization proba-
bility and the one obtained without mass shift is shown in 0.5 ω=3.5 a.u., NCycles=15
Fig. (5) – along with the TDDE results. Indeedwe find that TDSE BYD
the “relativistic substitution”, Eq. (41), shifts the ionization 0.4 BYD1
BYD2
probabilitydownwards. Moreover,althoughourmodelover-
BYD3
estimatesthediscrepancysomewhat,theionizationprobabil-
0.3 BYD4
ityobtainedbythemodelbehavesinamannerverysimilarto al BYD5
thetrulyrelativisticcalculations. u
d 0.2
si
e
R
0.1
ω=3.5 a.u., NCycles=15
0.02
TDDE Dip
0.0
TDSE Dip mass correction
0.01 −0.1
0 20 40 60 80 100
Peak electric field stength (a.u.)
al
u
sid 0.00 FIG. 6: (color online). Difference between Pion computed within
e
R and beyond the dipole approximation for TDSE and TDDE in re-
spectivecase.
−0.01
−0.02
50 60 70 80 90 100
ω=3.5 a.u., NCycles=15
Peak electric field stength (a.u.) 0.7
TDSE dip
FIG.5:(coloronline).Residualplotsfortherelativisticdipolecalcu- TDSE BYD
lationandthenon-relativisticdipolecalculationwiththemasssub- 0.6 Dip NoNES
BYD1 NoNES
stitutedasinEq.(41).
eld BYD2 NoNES
This is a strong indication that the dominating relativistic n yi0.5 BYD3 NoNES v≈0.1c
o
effectsforthefieldsconsideredhereoriginatefromdipolein- ti
a
z
teraction. In supportof this, Fig. 6 shows the differencebe- ni
o0.4
tween calculations performed within and beyond the dipole I
approximation–bothfortheTDSEandtheTDDE.Itisseen
that the corrections to the dipole approximation in the two
0.3
frameworksdoinfactcoincide. Althoughitis hardto judge
fromFig.3aloneifntrunc = 5issufficientforE0 = 90a.u., 5 10 15 20 25 30 35 40 45 50
theagreementinFig.6providesstrongsupportthereof. Peak electric field stength (a.u.)
FIG.7: (coloronline). AcomparisonofP fortheTDSEandthe
ion
TDDE calculated with and without the time-independent negative-
C. Dealingwiththenegativeenergystates
energy statesin therelativisticbasis. Withthislimitation, none of
thesimulationsareabletoprovideanycorrectiontothedipoleap-
In Ref. [6] it was shown that exclusion of the time- proximation.
independentnegative-energystates from the propagationba-
sisremovedalleffectsbeyondthedipoleapproximation.Here
we showthatthisconclusionseemstoholdregardlessofthe
valueofn inEq.(11). InFig.7wepresentresultsfrom
trunc
solvingtheTDDEwithoutthetimeindependentnegativeen- evaluated within a Krylov subspace, cf. Eq. (38), since the
ergy states, i.e., with those eigenstates of H corresponding obtainedspectrumisnotthesameasthetruetimedependent
0
to negative eigenenergies excluded from the basis set, for spectrumofH(t). Totestthis,twosimulationsfortheTDDE
n 3. It is seen that all these different Hamiltonians BYD5, withE = 70a.u. and80a.u.,wereperformedboth
trunc ≤ 0
nowpredictthesameP . with and without the projection imposed within the Krylov
ion
In Sec. IIC it was argued that the projection onto time- subspace, cf. Eqs. (20, 38). As it turns out, the only differ-
dependent (dynamic) positive energy states, i.e., eigenstates ence found between these calculations was slightly different
of the dynamical Hamiltonian H(t), was adequate as long convergencepropertiesin the time step, which indicatesthat
as the fields are not within the pair productionregime. This even in the approximate propagation method the arguments
should howevernotbe taken for grantedif the propagatoris madeinSec.IICshouldhold.
v≈0.2c
9
V. CONCLUSION mand of increasingly higher order corrections. A better ap-
proach is to look for a transformationof the Hamiltonian to
Wehavesolvedthetime-dependentDiracequationforahy- a form where the consistent inclusion of higher order multi-
drogenatomexposedtoextremelaserpulses. Uponcompari- polesoftheelectromagneticfieldiseasiertoachieve. Sucha
sonwiththenon-relativisticcounterpart,i.e.,theSchro¨dinger transformationwillbepresentedinaforthcomingpaper.
equation,itwasfoundthateffectsbeyondthedipoleapprox-
imation are morecomplicatedto incorporatecorrectlyin the
relativistic framework. Whereas first order space-dependent
Acknowledgments
corrections to the vector field are sufficient for the TDSE,
the TDDE demands an expansion to at least third order to
evenreproducethenon-relativisticresultsbelowtherelativis- The authors would like to thank prof. M. Førre for fruit-
tic regime. With increasing field strengths the demand for fuldiscussionsandforprovidinguswithdataforcomparison.
higherordercorrectionsincreasesevenfurtherwith theneed SignificantdevelopmentoftheTDSEandTDDEcodeswere
of a fifth-order space-dependent correction when the quiver done at resources provided by the Swedish National Infras-
velocityisv 0.19c. tructure for Computing (SNIC) HPC2N and NSC. We espe-
quiv
Emergingrelat≈ivisticcorrectionsarefoundintheionization cially would like to thank A˚ke Sandgren at HPC2N for as-
yield P for the test cases starting at v 0.17c. Both sistanceinthisdevelopment. TheNorwegianMetacenterfor
ion quiv
≈
withinandbeyondthedipoleapproximationtherelativisticef- ComputationalScience(UNINETTSigma2)isacknowledged
fectsgivealowerP suggestingincreasedstabilizingeffect forprovidingcomputationalresourcesforcarryingoutthecal-
ion
against ionization. It was demonstrated that this shift could culations(AccountnumberNN9417K).Financialsupportby
beexplainedbytheelectron’sincreasedrelativisticinertiade- theSwedishResearchCouncil(VR),GrantNo. 2012-3668is
scribedalreadywithinthedipoleapproximation. gratefully acknowledged. We also acknowledge support for
Toincreasethefieldstrengthfurtherandrevealstrongerrel- ourcollaborationthroughtheNordicInstituteforTheoretical
ativistic effects, possibly also relativistic corrections beyond Physics(Nordita)andfromtheresearchgroupMathematical
thedipoleapproximation,thepresentformofthelight-matter Modelling at Oslo and Akershus University College of Ap-
interactionintheDiracequationisnotsuitableduetothede- pliedSciences.
[1] A.DiPiazza,C.Mu¨ller,K.Z.Hatsagortsyan,andC.H.Keitel, gan,Phys.Rev.A85,032701(2012).
Rev.Mod.Phys.84,1177(2012). [15] R. Beerwerthand H.Bauke, Comput. Phys. Comm. 188, 189
[2] C.I.Moore,A.Ting,S.J.McNaught,J.Qiu,H.R.Burris,and (2015).
P.Sprangle,Phys.Rev.Lett.82,1688(1999). [16] S.Blanes, F.Casas, J.Oteo,andJ.Ros,PhysicsReports470,
[3] A. D. DiChiara, I. Ghebregziabher, R. Sauer, J. Waesche, 151 (2009).
S.Palaniyappan,B.L.Wen,andB.C.Walker,Phys.Rev.Lett. [17] H.R.Reiss,Phys.Rev.A63,013409(2000).
101,173002(2008). [18] K.J.Meharg,J.S.Parker,andK.T.Taylor,J.Phys.B38,237
[4] H.Yumoto, H.Mimura, T.Koyama, S.Matsuyama, K.Tono, (2005).
T.Togashi, Y. Inubushi, T.Sato, T.Tanaka, T.Kimura, et al., [19] M.FørreandA.S.Simonsen,Phys.Rev.A90,053411(2014).
Nat.Photon7,43(2013). [20] L.L.FoldyandS.A.Wouthuysen,Phys.Rev.78,29(1950).
[5] H.Yoneda,Y.Inubushi,M.Yabashi,T.Katayama,T.Ishikawa, [21] M.HochbruckandC.Lubich, SIAMJ.Numer.Anal.41, 945
H.Ohashi,H.Yumoto,K.Yamauchi,H.Mimura,andH.Kita- (2003).
mura,Nat.Commun5,5080(2014). [22] Y. I. Salamin, S. Hu, K. Z. Hatsagortsyan, and C. H. Keitel,
[6] S. Selstø, E. Lindroth, and J. Bengtsson, Phys. Rev. A 79, Phys.Rep.427,41 (2006).
043418(2009). [23] G. R. Mocken and C. H. Keitel, Phys. Rev. Lett. 91, 173202
[7] Y.V.VanneandA.Saenz,Phys.Rev.A85,033411(2012). (2003).
[8] A. S. Simonsen, T. Kjellsson, M. Førre, E. Lindroth, and [24] C. de Boor, A Practical Guide to Splines (Springer-Verlag,
S.Selstø,Phys.Rev.A93,9(2016). 1978),revisededition2001ed.
[9] M.S.Pindzola,J.A.Ludlow,andJ.Colgan,Phys.Rev.A81, [25] C. F. Fischer and O. Zatsarinny, Comput. Phys. Comm. 180,
063431(2010). 879(2009).
[10] H.Bauke, H.G.Hetzheim, G.R.Mocken, M.Ruf,andC.H. [26] I. Lindgren and J. Morrison, Atomic Many-Body Theory, Se-
Keitel,Phys.Rev.A83,063414(2011). riesonAtomsandPlasmas(Springer-Verlag,NewYorkBerlin
[11] G. R. Mocken and C. H. Keitel, Comput. Phys. Comm. 178, Heidelberg,1986),2nded.
868 (2008). [27] M.HochbruckandC.Lubich,SIAMJ.Numer.Anal.34,1911
[12] F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk, Comput. (1997).
Phys.Comm.183,1403 (2012). [28] M.Gavrila,J.Phys.B35,R147(2002).
[13] F.Fillion-Gourdeau, E.Lorin, and A. D. Bandrauk, J. Comp. [29] A.S.SimonsenandM.Førre,Phys.Rev.A92,013405(2015).
Phys.272,559 (2014). [30] M.Førre,privatecommunication(2015).
[14] M.S.Pindzola,S.A.Abdel-Naby,F.Robicheaux,andJ.Col-