Implementation of the multiconfiguration time-dependent Hatree-Fock method for general molecules on a multi-resolution Cartesian grid

1.9 MB·

by  Ryohto Sawada

Implementation of the multiconfiguration time-dependent Hatree-Fock method for
general molecules on a multi-resolution Cartesian grid
Ryohto Sawada,1,2 Takeshi Sato,2,3 and Kenichi L. Ishikawa2,3
1Department of Applied Physics, Graduate School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
2Photon Science Center, Graduate School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
3Department of Nuclear Engineering and Management,
Graduate School of Engineering, The University of Tokyo,
6 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
1 (Dated: March 23, 2016)
0
We report a three-dimensional numerical implementation of multiconfiguration time-dependent
2
Hartree-Fock (MCTDHF) based on a multi-resolution Cartesian grid, with no need to assume any

r symmetryofmolecularstructure. Wesuccessfullycomputehigh-harmonicgeneration(HHG)ofH
a 2
and H O. The present implementation will open a way to the first-principle theoretical study of
M 2
intense-field and attosecond-pulse induced ultrafast phenomena in general molecules.

2
PACSnumbers: 33.80.Rv,31.15.A-,42.65.Ky,
2

]
h I. INTRODUCTION time. The multiconfiguration time-dependent Hartree-
p Fock (MCTDHF) [13, 14] considers all the possible elec-
- The dynamics of atoms and molecules under intense tronic configuration for a given number of spin or-
p
m (typically (cid:38) 1014W/cm2) laser pulses is of great inter- bitals. As its flexible generalizations, we have re-
est in a variety of fields such as attosecond science and cently formulated the TD complete-active-space self-
o high-field physics [1–3], with a goal to directly measure consistent field (TD-CASSCF) [8] and TD occupation-
c
. and manipulate electronic motion. Numerical simula- restricted multiple-active space (TD-ORMAS) [9] meth-
s tions of such electron dynamics are a challenging task ods. The latter is valid for general MCSCF wave func-
c
i [4]. Direct solution of the time-dependent Schr¨odinger tions with arbitrary CI spaces [4, 15] including, e.g.,
ys equation (TDSE) cannot be applied beyond He, H2, and the TD restricted-active-space self-consistent-field (TD-
h Li,duetoaprohibitivecomputationalcost. Thus,oneof RASSCF) theory developed by Miyagi and Madsen [16].
p major recent directions attracting increasing interest is Numerical implementations of MCTDHF for atoms as
[ themulticonfiguration self-consistent-field(MCSCF) ap- well as diatomic molecules have been reported for the

proach, which expresses the total wave function Ψ(t) as calculationofvalenceandcorephotoionizationcrosssec-
3
v a superposition [5–10] tions [17]. We have also implemented TD-CASSCF for
atomsbyexpandingorbitalfunctionswithsphericalhar-
9
8 |Ψ(t)(cid:105)=(cid:88)c (t)|J(cid:105), (1) monics and successfully computed high-harmonic gener-
J
3 ation and nonsequential double ionization of Be [18].
J
1 Practically all the existing implementations are in-
0
of Slater determinants |I(cid:105) built from the spin orbitals tended for atoms and diatomic molecules, exploiting the
.
1 |φ (cid:105) = |φ (cid:105)⊗|σ(cid:105), where {φ } and σ ∈ {α,β} denote underlying symmetries with either the spherical [19–23],
(i,σ) i i
0 one-electronspatialorbitalfunctionsandspineigenfunc- cylindrical [24–27], or prolate spheroidal [28–31] coordi-
6 tions, respectively. Different variants with this ansatz nates.
1
have recently been actively developed [4]. In this study, we report a three-dimensional (3D) nu-
:
v The time-dependent configuration-interaction (TDCI) merical implementation of MCTDHF based on a multi-
Xi methods take the orbital functions to be time- resolution Cartesian grid, with no need to assume any
independent and propagate only CI coefficients c (t). symmetry of molecular structure, this can in principle
r I
a Santra et al. [11] have implemented its simplest variant, be applied to any molecule. With the use of a multi-
i.e., the time-dependent configuration-interaction singles resolution finite-element representation of orbital func-
(TDCIS)methodtotreatatomichigh-fieldprocesses. In tions, we can fulfill a high degree of refinement near nu-
this method, only up to single-orbital excitation from clei and, at the same time, a simulation domain large
the Hartree-Fock (HF) ground-state is included. Bauch enough to sustain departing electrons. As demonstra-
elal. [12]haverecentlydevelopedTDgeneralized-active- tions, wesuccessfullycomputehigh-harmonicgeneration
space CI based on a general CI truncation scheme and (HHG) from H and H O. The present implementation
2 2
discusseditsnumericalimplementationforatomsanddi- will open a way to the first-principle theoretical study of
atomic molecules. intense-field and attosecond-pulse induced ultrafast phe-
In the other class of MCSCF approaches, not only nomena in general molecules.
CI coefficients but also orbital functions are varied in Thispaperisorganizedasfollows. InSec.II,webriefly
2
summarize the MCTDHF method. Section III describes with,
the multi-resolution cartesian grid. Section IV explains
the numerical procedure that we implement. In Sec. V, Pˆ =1ˆ−(cid:88)M |φ (cid:105)(cid:104)φ |, (9)
we show examples of simulation results for He, H , and j j
2
H O. Conclusions are given in Sec. VI. Atomic units are j=1
2
(cid:88)
used throughout unless otherwise stated. ρ = (cid:104)Ψ|aˆ† aˆ |Ψ(cid:105), (10)
i,j iσ jσ
σ
ρ(2) =(cid:88)(cid:104)Ψ|aˆ† aˆ† aˆ aˆ |Ψ(cid:105), (11)
jklm jσ lτ mτ kσ
II. MCTDHF στ
(cid:90) 1
g (x)= dx(cid:48)φ∗(x(cid:48)) φ (x(cid:48)), (12)
In the MCTDHF method [13, 14], the sum in Eq.(1) lm l |x−x(cid:48)| m
runs over the complete set of (cid:0)M(cid:1)(cid:0)M(cid:1) Slater determi-
Nα Nβ where 1ˆ denotes the identity operator, and aˆ† and aˆ
nants|J(cid:105)thatcanbeconstructedfromNα electronswith iσ iσ
the Fermion creation and annihilation operators, respec-
spin-projection α, N electrons with spin-projection β,
β
tively,associatedwithspatialorbitaliandspinσ. Equa-
and M spatial orbitals. Their spin-projection is conse-
tion (12) is computed by solving the Poisson equation,
quently restrited to S =(N −N )/2.
z α β
Let us consider a Hamiltonian in the length gauge, ∇2g (x)=−4πφ∗(x)φ (x). (13)
lm l m
H(t)=H (t)+H , (2) It is convenient to rewrite Eq. (8) as,
1 2
iφ˙ =Pˆ(Tφ +W (t)) (14)
i i i
H (t)=(cid:88)N (cid:32)−∇2i −(cid:88) Za +x ·E(t)(cid:33),(3) where T is kinetic energy.
1 2 |x −X | i
i a
i=1 a
III. MULTI-RESOLUTION CARTESIAN GRID
N i−1
(cid:88)(cid:88) 1
H = , (4) We discretize spatial orbital functions on a multi-
2 |x −x |
i j resolutionCartesiangrid,inspiredbytheworkofBischoff
i=1j=1
and Valeev [35] and the finite volume method [36]. Fig-
where N =N +N , X and Z are the charge and po- ure 1 (a) schematically shows how to generate it. We
α β a a
sitionofthea-thatom,respectively,andE(t)isthelaser startfromanequidistantCartesiangridcomposedofcu-
electronic field. One can derive the equations of motion bic cells. If a given cell is too large to represent orbital
for the CI coefficients c (t) and spatial orbital functions functions with sufficient accuracy, typically near the nu-
J
φ (t), resorting to the time dependent variational princi- clei, we subdivide it into eight cubic cells with half the
i
ple [32–34], side length of the original cell. We continue the subdivi-
sionuntilaccuracyrequirementsaresatisfied. Thecenter
(cid:18)(cid:90) t2 (cid:19) of each cube is taken as the grid point representing the
δ (cid:104)Ψ|H(t)−i∂ |Ψ(cid:105)dt =0 (5)
t cell.
t1 The Laplacian ∇2φ of orbital function φ(r,t) is eval-
uated at each grid point by finite difference. We first
with additional constraints for uniqueness [13],
illustrate it for a one-dimensional case for simplicity in
Fig. 1(b). One can evaluate the second derivative of a
(cid:28) ∂φ (cid:29) function f(x) at grid point xi as,
(cid:104)φ |φ (cid:105)=δ , φ | k =0. (6)
j k j,k j ∂t d2 g+−g−
f(x )≈ i i , (15)
dx2 i ∆x
i
The equations of motion are,
where ∆x is the size of cell i, and the first derivatives
i
ic˙ =(cid:88)(cid:104)J|H(t)|K(cid:105)c , (7) g± at the cell boundaries are approximated by,
J K i
K
2[f(x )−f(x )]
g+ ≈ i+1 i , (16)
and i ∆x +∆x
i i+1
  g− ≈ 2[f(xi)−f(xi−1)]. (17)
i|φ˙i(cid:105)=PˆH1(t)|φi(cid:105)+ (cid:88)(ρ−1)ijρ(j2k)lmgˆlm|φk(cid:105), i ∆xi+∆xi−1
jklm We show the extension to two dimensions in Fig. 1(c).
(8) The grid points are marked by red and blue circles. In
3
order to evaluate the second derivative with respect to at the simulation boundary ∂V is given by the multipole
the vertical direction at the center of cell i, we need the expansion
firstderivativeevaluatedatthecellboundarymarkedby
(cid:90) 1
theorangetriangle,forwhichweneed,inturn,thevalue g (x )= φ∗(x(cid:48))φ∗ (x(cid:48))dx(cid:48) (22)
of the function at the position marked by the star in cell lm bound V |xbound−x(cid:48)| l m
if(+ri1+.1)Watetahpepgrroixdimpoaitnet,thi.ies.,lathtteercebnyterthoef vthaeluceelfli+i+1 ≡1. =(cid:88)∞ (cid:90) |x|x(cid:48)|lPl(c|ols+θ1)φ∗l(x(cid:48))φ∗m(x(cid:48))dx(cid:48)
Thoughinferiorintermsofaccuracy,thisschemeismuch l=0 V bound
(23)
more advantageous in terms of computational cost over
conventional methods such as the alternating direction
for x ∈ ∂V where P (z) denotes the Legendre
implicit method [37], moving least squares [38, 39], and bound l
polynomial, and θ the angle between x(cid:48) and x . In
symmetric smoothed particle hydrodynamics [40, 41]. bound
the present study, we truncate the sum in Eq. (23) at
Then, in the 3D case, we evaluate the Laplacian
l=2 (second-order multipole expansion).
∇2φ(r) as,
Step 4: Time propagation of c and φ
(cid:88) J i
∇2φ(r )≈L φ(r )+ (cid:48)L φ(r ), (18)
a aa a ab b
b WesolvetheequationsofmotionEqs.(7)and(14)us-
ingasecond-orderexponentialintegrator[44,45]. Equa-
where a and b are cell indices, the primed sum is taken
tion (7) is integrated as,
over the cells adjacent to the a-th cell, and,
c(1)(t+∆t)=c (t)+∆t(cid:88)(cid:104)J|H|K(cid:105)c , (24)
L =−(cid:88)(cid:48)L , (19) J J K
aa ab K
b c(2)(t+∆t)=c(1)(t+∆t)
l2 2 1 J J
Lab = lab2 la+lbla (a(cid:54)=b and lb <la), (20) +∆t(cid:88)(cid:104)J(1)|H|K(1)(cid:105)c(K1)(t+∆t), (25)
2 1 K
Lab = la+lbla (a(cid:54)=b and lb ≥la), (21) c (t+∆t)= cJ(t)+c(J2)(t+∆t), (26)
J 2
with l being the side length of the a-th cell. If l < l ,
a b a
a face of a cell of side length l would contact with l2/l2 where|J(1)(cid:105)withsuperscript“(1)”denotestheSlaterde-
adjacentcellsofsidelengthl .aTheprefactorl2/l2 ofaEqb. terminantconstructedwithorbitalfunctionsφ(1) defined
b b a i
(20) takes into account the weight of each of the latter. below in Eq. (27). Equation (14) is integrated as,
φ(1) =φ (t)
IV. NUMERICAL PROCEDURE i i
1
+Pˆ [(−i∆tT)φ (t)+∆tW (t)], (27)
1+i∆tT/2 i i
We present the essential steps of MCTDHF simula-
tions using multi-resolution cartesian grid as follows:
φ(2)(t+∆t)=φ(1)
i i
Step 1: Generation of grid and Laplacian matrix
1 (cid:104) (cid:105)
+Pˆ(1) (−i∆tT)φ(1)+∆tW (t+∆t) ,
1+i∆tT/2 i i
We consider a cuboid simulation region V centered at
(28)
the origin: {x = (x,y,z) ∈ R2 | x ∈ [−x ,x ],y ∈
L L
[−y ,y ],z ∈ [−z ,z ],x > 0,y > 0,z > 0}. We
L L L L L L L
set the locations of the grid points and prepare the M
Laplacian matrix elements L using Eqs. (19)–(21). Pˆ(1) =1ˆ−(cid:88)|φ(1)(cid:105)(cid:104)φ(1)|, (29)
ab j j
These are done only once in the beginning. j=1
Step 2: Computation of ρ and ρ(2) φ(2)(t+∆t)+φ (t)
φ (t+∆t)= i i . (30)
i 2
Each time step starts with the computation of ρ and
ρ(2), using Eqs. (10) and (11), respectively. In Eqs. (27) and (28), (1+i∆tT/2)−1 is operated by the
conjugate residual method [42, 43].
Step 3: Computation of g
lm
Step 5a: Absorbing boundary (only in real time propaga-
We solve the Poisson equation (13) to obtain g , by tion)
lm
the conjugate residual method [42, 43]. The condition
4
(a)(cid:1)
divide(cid:1) divide(cid:1) Δx
i
f (cid:1)
(b)(cid:1) (c)(cid:1) l i+1
f(x)(cid:1) a cell i+1(cid:1)
i f+(cid:1)
f-(cid:1) f+(cid:1) i
i i
l
b fi(cid:1)
f(x )(cid:1) f(x )(cid:1)
i-1 i+1
Δx cell i(cid:1)
i
FIG.1. (a)Schematicofthecartesian-basedmulti-resolutiongrids. (b)Schematicofthecomputationofsecondorderdifferential
at one-dimensional irregular grids. (c) Schematic of the computation of the first differential at the surface of the grid at two-
dimensional irregular grids. Real grid points are red and blue circles.
To prevent the reflection from the grid boundaries, af- V. EXAMPLES
ter each time step, φ is multipled by a cos mask func-
i
tion M(x,y,z) that varies from 1 to 0 between the ab- A. Benchmark : HHG from helium
sorption boundary set at x = x , y = y , and z = z
0 0 0
(0 < x < x ,0 < y < y ,0 < z < z ) and the outer
0 L 0 L 0 L We simulate the HHG from a helium atom located at
boundary ∂V [46, 47]:
the origin. The side length of the cell is set to be 0.6
(cid:18)|x|−x (cid:19) (cid:18)|y|−y (cid:19) (cid:18)|z|−z (cid:19) (r > 4), 0.3 (2 < r < 4) and 0.15(r < 2) respectively,
M(x,y,z)=C 0 C 0 C 0
dependingonthedistancerofthegridpointatthecenter
x −x y −y z −z
L 0 L 0 L 0
of each cell and the origin. We also set x = 70,y =
(31) L L
z = 35 and x = 0.7x ,y = 0.7y ,z = 0.7z . The
L 0 L 0 L 0 L
where time step size ∆t is set to be 0.0025. We consider a laser
pulse linearly polarized along the x axis, whose electric
C(x)=1 (x≤0), cos(x) (x>0). (32)
field E(t) is given by,
Alternatively,onemayuse,e.g.,exteriorcomplexscaling
E(t)=E (t)sinωt, (33)
[48, 49]. env

ωt/2π (ωt<2π),

Step 5b: Rescaling of cJ and orthonormalization of φi Eenv(t)= 2−ωt/2π (2π <ωt<4π), (34)
(only in imaginary time propagation) 0 (otherwise),
We obtain the initial ground state via the imaginary withacentralwavelengthof400nmandapeakintensity
time propagation [50]. After each (imaginary) time of8×1014 W/cm2. Forsuchanultrashortpulse,thecut-
step, c is rescaled so that (cid:80) |c |2 = 1, and φ is offenergypredictedbythesemiclassicalthreestepmodel
J J J i
orthonormalized through the Gram-Schmidt algorithm. [51, 52] is I +2.07U = 49.6eV, which corresponds to
p p
the 16.0 th order where I is ionization potential and
p
Step 6: End of time step U is pondermotive energy. The harmonic spectrum is
p
obtained from the Fourier transform of the dipole accel-
We go back to Step 2 to start next time step. eration.
5
(cid:1) W/cm2, and an eight-cycle sine-squared envelope,
)
t(
ni E(t)=E sin2(ωt/16)cos(ωt). (35)
u 0
b.(
r Figure3presentstheHHGspectraforlaserpolarization
a
( parallel to the molecular axis (the x axis) [Fig. 3(a)(b)]
ty(( (cid:1) and 30 degrees from the molecular axis [Fig. 3(c)]. The
nsi u. ) cutoff energy predicted by the semiclassical three step
te a. model is 34.3 eV, which corresponds to order 22.1. One
n (
I can see that the simulation is converged with respect
to the number of orbitals [Fig. 3(a)] and grid spacing
[Fig. 3(b)]. Our multi-resolution Cartesian-grid MCT-
DHF, with no a priori assumption of symmetry, can also
FIG. 2. HHG spectrum of helium computed by multi- handle laser polarization oblique to the molecular axis
resolution MCTDHF(red-solid) and method of [18](black- [Fig. 3(c)].
dashed). The arrow represents the cutoff energy.
In Fig. 3 we can clearly see the second plateau, some-
what weaker than the first one, extending beyond the
number of largest cell side length l0 cutoff (∼ order 22.1). The second cutoff position is con-
orbital M 0.7 0.7∗ 0.6 0.55
sistentwiththevalue(53.6 eVorthe34.6-thorder)pre-
1 -1.83661 -1.83622 -1.84318 -1.84123
dictedbythethreestepmodelwiththeionizationpoten-
2 -1.85451 -1.85467 -1.86164 -1.85964
tialofH+ (34.7 eV). Hence,basedonaspeculationthat
3 -1.86218 -1.86233 -1.86925 -1.86723 2
thesecondplateauharmonicsaregeneratedfromH+pro-
6 -1.87329 -1.87342 -1.88027 -1.87756 2
duced via strong-field ionization, we have simulated the
HHGfromthismolecularionwiththesamelaserparam-
TABLE I. Ground-state energy (a.u.) of a hydrogen
molecule,obtainedbyrelaxationinimaginarytime. Theval- eters. The obtained harmonic spectrum multiplied with
uesincolumnlabeled“0.7∗”areobtainedwithgridsdisplaced theionizationprobabilityofH2 (2.4×10−4)isplottedas
parallel to the x axis by 0.025 a.u. a yellow dashed line in Fig. 3(a). The spectrum is much
weaker than the second plateau from H .
2
Presumably, the harmonic response from H+ is sub-
2
In Fig. 2 we compare the HHG spectrum calculated stantiallyenhancedbytheactionoftheoscillatingdipole
with the present implementation with that calculated formed by the recolliding first electron ejected from the
with another implementation in spherical coordinates neutral molecule and the neutral ground state. This me-
[18]. One can see that they agree with each other very chanics is similar to enhancement by an assisting har-
well. monicpulse[54–57],buttheenhancementisduetodirect
Coulomb force from the oscillating dipole, rather than
harmonics emitted from it. In the words of the semiclas-
B. HHG from a hydrogen molecule sical three-step model, the recolliding first electron vir-
tually excites H+, facilitating second ionization. Thus,
2
Next, we simulate the HHG from molecular hydrogen electron-electron interaction plays an important role in
where two hydrogen atoms are located at (±0.7,0,0), high-harmonic generation in some cases (see also [58]),
respectively. The side length of cell is set to be l whereas HHG is usually considered as a predominantly
0
(r >4),l /2(2<r <4)andl /4(r <2),respectively single-electron process.
0 0 0 0 0
where l is the side length of the largest cells (see Table
0
I for its values). We also set x = y = z = 27 and
L L L
x = 0.7x ,y = 0.7y ,z = 0.7z . The time step size C. HHG from a water molecule
0 L 0 L 0 L
∆t is set to be 0.01.
The ground-state energy, obtained through relaxation As an example of application to molecules of lower
in imaginary time, is shown in Table I where M is the symmetry, we simulate the HHG from a water molecule
numberoforbitals. Itconsistentlytendstotheliterature with its oxygen atom located at the origin and two hy-
value -1.8884 a.u. [53] with an increasing number of drogen atoms at (±1.4299,1.10718,0). The side length
orbitals. The slight dependence on M and l has only a of cell is set to be 0.6 (r > 4), 0.3 (2 < r < 4) and
0 0 0
small impact on calculated harmonic spectra, as we will 0.15(r < 2) respectively where r is the distance from
0 0
see below in Fig. 3(a,b). The values in column labeled the nearest atom. The outer boundary x ,y ,z is set
L L L
“0.7∗”areobtainedwithgridsdisplacedparalleltothex to be 60 (axis parallel to the polarization) and 30 (axis
axis by 0.025. One can see that the resulting loss of grid parpendicular to the polarization) and the absorption
symmetry with respect to the yz plane also has only a boundary x ,y ,z is set to be 0.7 times as long as the
0 0 0
small impact. outerboundary. Thetimestepsize∆tissettobe0.0025.
Let us consider a linearly polarized laser pulse with a We use the same laser pulse shape as in Sec. VA. The
centralwavelengthof800nm,apeakintensityof1×1014 cutoff energy predicted by the semiclassical three step
6
(cid:1)(cid:1) three different directions of laser polarization, demon-
b.(unit()(cid:1)unit()b.(unit() (a)(cid:1) sgcturrirdavteeMssCohbiTgtDhaiHnfleeFdxiibmwilipittlhyemMoefn=tth5aetiamonnud.lt6iO-ranelsemolcouastntioonsveeeCraltarhptaetswiatihtnhe-
arb.(ar
sity(((y(((arsity((( eoancha soitnhgelre.nTohdeeswimituhlatwtioonhewxiath-coMre=3.363toGoHkzcaX.e2o8ndparyos-
ntn cessors. In this case, the computational bottleneck was
esie
tnt thesolutionofPoisson’sequation(Step3ofSecIV).We
nen
IntI expectthatthedistributedparallelizationofthecodewill
I
(cid:1)(cid:1) substantiallyreducethecomputationaltime,andtheex-
))
t((cid:1)t( tensiontoTD-CASSCF[8]andTD-ORMAS[9]methods
b.(uniunit()b.(uni (b)(cid:1) will further extend the applicability to larger systems.
arb.(ar
sity(((y(((arsity((( VI. CONCLUSION
ntn
esie
tnt
nen We have numerically implemented the MCTDHF
ItI
In method on a multi-resolution Cartesian grid. Whereas
previous approaches have relied on the underlying sym-
metriesofthesimulatedatomsandmolecules,thepresent
(cid:1)(cid:1)
t()(cid:1)t() (c)(cid:1) implementation offers a flexible framework to describe
b.(uniunit()b.(uni smtroolencgu-lfieesl.dEaxntednsaiottnosteococnomdppurtoacteiossneasllyofmroeraelcogmenpearcatl
arb.(ar methods such as TD-CASSCF [8] and TD-ORMAS [9]
sity(((y(((arsity((( wlairlglebme oraletchuelress.traightforward and enable application to
ntn
tensite As demonstrations, we have successfully calculated
nen
IntI high-harmonic spectra from He, H2, and H2O. As
I the presence of the second plateau in Fig. 3 implies,
the present implementation will uncover yet unexplored
multi-electron, multi-channel, and multi-orbital effects,
FIG.3. Calculatedhigh-harmonicspectrafromahydrogen which only first-principles simulations can reveal.
molecule. See text for laser parameters. (a) Comparison of
the results with M = 1(blue-dotted), M = 2 (black-solid),
and M = 3(purple-solid), for laser polarization parallel to
the molecular axis. Yellow-dashed line: spectrum from H+ ACKNOWLEDGMENTS
2
multipliedwiththeionizationprobabilityofH (2.4×10−4).
2
(b) Comparison of the results with l = 0.7 (red-solid) and
0 This work was supported in part by Japan Society
0.55(black-dashed). ThecalculationwasdonewithM =1for
for the Promotion of Science (JSPS) KAKENHI Grants
polarizationparalleltothemolecularaxis. (c)Resultforlaser
No.25286064, No. 26390076, No. 26600111, and No. 26-
polarization 30 degrees from the molecular axis [polarization
direction is (cos30◦,sin30◦,0)]. M =3 was used. Red-solid: 10100. This research was also partially supported by
the Photon Frontier Network Program of the Ministry
harmonics emitted in the y direction, black-dashed: in the x
direction. Arrows in each panel indicate the cutoff positions of Education, Culture, Sports, Science and Technology
expected for H (22.1-th order) and H+ (34.6-th order). (MEXT) of Japan, the Advanced Integration Science In-
2 2
novation Education and Research Consortium Program
ofMEXT,theCenterofInnovationProgramfromJapan
model is 37.3 eV, which corresponds to the 12.0 th or- Science and Technology Agency (JST), and Core Re-
der. search for Evolutional Science and Technology, Japan
Figure 4, which presents the harmonic spectra for Science and Technology Agency (CREST, JST).
[1] F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163 Camillis, S. Anumula, F. Frassetto, L. Poletto, A. Pala-
(2009). cios, P. Decleva, J. B. Greenwood, F. Mart´ın, and
[2] R. Levis, G. Menkir, and H. Rabitz, Science 292, 709 M. Nisoli, Science 346, 336 (2014).
(2001). [4] K. Ishikawa and T. Sato, IEEE J. Sel. Topics Quantum
[3] F. Calegari, D. Ayuso, A. Trabattoni, L. Belshaw, S. D. Electron. 21, 8700916 (2015).
7
(a)(cid:1) (b)(cid:1) (c)(cid:1)
101 O(cid:1) 101 O(cid:1) 101
100 H(cid:1) H(cid:1) 100 H(cid:1) H(cid:1) 100
cut off(12.2)(cid:1)
10-1 10-1 10-1
10-2 polarization(cid:1) 10-2 polarization(cid:1) 10-2
(cid:1)
a. u. ) 1100--43 cut off(12.2)(cid:1) 1100--43 cut off(12.2)(cid:1) 1100--43 H(cid:1) O(cid:1)H(cid:1)
y ( 10-5 10-5 10-5
sit 10-6 10-6 10-6 polarization(cid:1)
n
e 10-7 10-7 10-7
nt
i 10-8 10-8 10-8
10-9 10-9 10-9
10-10 10-10 10-10
10-11 10-11 10-11
10-12 10-12 10-12
5 10 15 20 25 30 35 5 10 15 20 25 30 35 5 10 15 20 25 30 35
Harmonic order(cid:1)
FIG.4. High-harmonicspectrafromawatermolecule,calculatedwithM =5(dashed)andM =6(solid),forlaserpolarization
along (a) the x axis (b) the y axis (c) the z axis, as indicated in each panel. Laser polarization in (c) is perpendicular to the
plane of the molecule. See text for laser parameters.
[5] T.-T. Nguyen-Dang, M. Peters, S.-M. Wang, and 013407 (2005).
F. Dion, Chem. Phys. 366, 71 (2009). [22] K. L. Ishikawa and K. Ueda, Phys. Rev. Lett. 108,
[6] T.-T. Nguyen-Dang and J. Viau-Trudel, J. Chem. Phys. 033003 (2012).
139, 244102 (2013). [23] K. L. Ishikawa and K. Ueda, Appl. Sci. 3, 189 (2013).
[7] R. P. Miranda, A. J. Fisher, L. Stella, and A. P. Hors- [24] K. Harumiya, I. Kawata, H. Kono, and Y. Fujimura, J.
field, J. Chem. Phys. 134, 244101 (2011). Chem. Phys. 113, 8953 (2000).
[8] T. Sato and K. L. Ishikawa, Phys. Rev. A 88, 023402 [25] K. Harumiya, H. Kono, Y. Fujimura, I. Kawata, and
(2013). A. D. Bandrauk, Phys. Rev. A 66, 043403 (2002).
[9] T. Sato and K. L. Ishikawa, Phys. Rev. A 91, 023417 [26] S. Ohmura, T. Oyamada, T. Kato, H. Kono, and
(2015). S. Koseki, eds., Molecular Orbital Analysis of High Har-
[10] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, The monic Generation, Vol. 1 (Asia Pacific Physics Confer-
Journal of Chemical Physics 127, 154103 (2007). ence, 2014).
[11] L. Greenman, P. J. Ho, S. Pabst, E. Kamarchik, D. A. [27] S. Ohmura, H. Kono, T. Oyamada, T. Kato, K. Nakai,
Mazziotti, and R. Santra, Phys. Rev. A 82, 023406 and S. Koseki, J. Chem. Phys. 141, 114105 (2014).
(2010). [28] X.Guan,K.Bartschat, andB.I.Schneider,Phys.Rev.
[12] S.Bauch,L.K.S.rensen, andL.B.Madsen,Phys.Rev. A 82, 041404 (2010).
90, 062508 (2014). [29] L.Tao,C.W.McCurdy, andT.N.Rescigno,Phys.Rev.
[13] J. Caillat, J. Zanghellini, M. Kitzler, O. Koch, A 79, 012719 (2009).
W. Kreuzer, and A. Scrinzi, Phys. Rev. A 71, 012712 [30] L.Tao,C.W.McCurdy, andT.N.Rescigno,Phys.Rev.
(2005). A 80, 013402 (2009).
[14] T.KatoandH.Kono,Chem.Phys.Lett.392,533(2004). [31] L.Tao,C.W.McCurdy, andT.N.Rescigno,Phys.Rev.
[15] D. J. Haxton and C. W. McCurdy, Phys. Rev. A 91, A 82, 023423 (2010).
012509 (2015). [32] P.A.M.Dirac,Proc.CambridgePhil.Roy.Soc.26,376
[16] H. Miyagi and L. B. Madsen, Phys. Rev. A 87, 062511 (1930).
(2013). [33] J.Frenkel, Wave Mechanics. Advanced General Theory.,
[17] D.J.Haxton,K.V.Lawler, andC.W.McCurdy,Phys. edited by C. Press (Clarendon Press, 1934).
Rev. A 86, 013406 (2012). [34] P. Kramer and M. Saraceno, Geometry of the Time-
[18] T. Sato and K. L. Ishikawa, Unpublished. Dependent Variational Principle, edited by Springer
[19] J. Parker, K. T. Taylor, C. W. Clark, and S. Blodgett- (Springer, 1981).
Ford, Journal of Physics B: Atomic, Molecular and Op- [35] F. A. Bischoff and E. F. Valeev, J. Chem. Phys. 134,
tical Physics 29, L33 (1996). 104104 (2011).
[20] E.S.Smyth,J.S.Parker, andK.Taylor,Comput.Phys. [36] H. Versteeg and W. Malalasekera, An Introduction to
Commun. 114, 1 (1998). Computational Fluid Dynamics: The Finite Volume
[21] K. L. Ishikawa and K. Midorikawa, Phys. Rev. A 72, Method (Prentice Hall, Upper saddle river, 2007).
8
[37] D. W. Peaceman and J. H. H. Rachford, Journal of the [47] M.Beck,A.Ja¨ckle,G.Worth, andH.-D.Meyer,Physical
Society for Industrial and Applied Mathematics 3, 28 Reports 324, 1 (2000).
(1955). [48] A.ScrinziandB.Piraux,Phys.Rev.A58,1310(1998).
[38] D.Levin,MathematicsofComputation67,1517(1998). [49] A. Scrinzi, Phys. Rev. A 81, 053845 (2010).
[39] C. L. Lopreore and R. E. Wyatt, Phys. Rev. Lett. 82, [50] R. Kosloff and H. Tal-Ezer, Chem. Phys. Lett. 127, 223
5190 (1999). (1986).
[40] R.C.BatraandG.M.Zhang,ComputationalMechanics [51] P. B. Corkum, Phys. Rev. Lett. 71, 1994 (1993).
41, 527 (2007). [52] K. Kulander, K. Schafer, and J. L. Krause, in Super-
[41] C.L.Tsai,Y.L.Guan,R.C.Batra,D.C.Ohanehi,J.G. Intense Laser-Atom Physics, NATO ASI, Ser. B, Vol.
Dillard, E. Nicoli, and D. A. Dillard, Computational 316,editedbyB.Piraux,A.L’Huillier, andK.Rza¸z˙ewski
Mechanics 51, 19 (2012). (Plenum Press, New York, 1993) p. 95.
[42] E. Steiefel, Commentarii Mathematici Helvetici 29, 157 [53] A. V. Turbiner and N. L. Guevara, Collection of
(1955). CzechoslovakChemicalCommunications72,164(2007).
[43] Y. Saad, Iterative Methods for Sparse Linear Systems, [54] K. Ishikawa, Phys. Rev. Lett. 91, 043002 (2003).
edited by S. for Industrial and A. Mathematics (Society [55] K. L. Ishikawa, Phys. Rev. A 70, 013412 (2004).
for Industrial and Applied Mathematics, 2003). [56] E.J.Takahashi,T.Kanai,K.L.Ishikawa,Y.Nabekawa,
[44] S. Cox and P. Matthews, Journal of Computational andK.Midorikawa,Phys.Rev.Lett.99,053904(2007).
Physics 176, 430 (2002). [57] K. L. Ishikawa, E. J. Takahashi, and K. Midorikawa,
[45] A. Bandrauk and H. Lu, Journal of Theoretical & Com- Phys. Rev. A 80, 011807 (2009).
putational Chemistry 12, 1340001 (2013). [58] K.-J. Yuan, H. Lu, and A. D. Bandrauk, Phys. Rev. A
[46] J. L. Krause, K. J. Schafer, and K. C. Kulander, Phys. 92, 023415 (2015).
Rev. A 45, 4998 (1992).