A fully general time-dependent multiconfiguration self-consistent-field method for the electron–nuclear dynamics

Abstract

We present a fully general time-dependent multiconfiguration self-consistent-field method to describe the dynamics of a system consisting of arbitrarily different kinds and numbers of interacting fermions and bosons. The total wave function is expressed as a superposition of different configurations constructed from time-dependent spin-orbitals prepared for each particle kind. We derive equations of motion followed by configuration–interaction (CI) coefficients and spin-orbitals for general, not restricted to full-CI, configuration spaces. The present method provides a flexible framework for the first-principles theoretical study of, e.g., correlated multielectron and multinucleus quantum dynamics in general molecules induced by intense laser fields and attosecond light pulses.

---

1 Introduction

We are now witnessing rapid progress in ultrashort intense light sources in different spectral ranges such as terahertz radiation, optical-parametric-chirped-pulse-amplification mid-infrared lasers, high-harmonic extreme-ultraviolet (XUV) pulses, and XUV/X-ray free-electron lasers. These technological advances have triggered various research activities, including attosecond science,^1–3 with a goal to directly measure and, ultimately, control electron and nuclear motion in atoms and molecules.

Ab initio simulations of the electronic and nuclear dynamics in atoms and molecules remain a challenge. A multiconfiguration time-dependent Hartree–Fock (MCTDHF) method^4,5 has been developed for the investigation of multielectron dynamics in strong and/or ultrashort laser fields.⁶ In this approach, the time-dependent total electronic wave function Ψ(t) is expressed as a superposition of different Slater determinants Φ_I(t),

where C_I(t) is the configuration–interaction (CI) coefficients. Both {C_I(t)} and the spin-orbitals constituting {Φ_I(t)} are allowed to vary in time. In the community of high-field phenomena and attosecond physics, the term MCTDHF is conventionally used for the full-CI case, in which the sum in eqn (1) runs over all the possible ways to distribute the electrons among a given number of spin-orbitals. On the other hand, also under active development are variants without the restriction to the full-CI expansion, generically referred to as the time-dependent multiconfiguration self-consistent-field (TD-MCSCF) methods hereafter. Representative examples include the time-dependent complete-active-space self-consistent-field,^7,8 the time-dependent restricted-active-space self-consistent-field,⁹ and the time-dependent occupation-restricted multiple active-space (TD-ORMAS)¹⁰ methods. These allow a compact and computationally less demanding description of the multielectron dynamics, without sacrificing the accuracy. In particular, the TD-ORMAS method can treat arbitrary CI expansions of the form eqn (1) in principle.

Among successful approaches for nuclear dynamics is the multiconfiguration time-dependent Hartree (MCTDH) method.¹¹ Developed for systems consisting of distinguishable particles, this method expresses the time-dependent total nuclear wave function as a superposition similar to eqn (1) but that of Hartree products. The other way around, the MCTDHF method can be viewed as an extension of the MCTDH method to fermions. By hybridizing the MCTDHF method for electrons and the MCTDH method for nuclei, one can construct a multiconfiguration electron–nuclear dynamics (MCEND) method^12,13 to describe the non-Born–Oppenheimer coupled dynamics. Nuclei forming molecules are, however, indistinguishable particles, either fermions or bosons. Alon et al. have explored methods for systems consisting of identical particles,^14–17e.g., MCTDH methods for mixtures consisting of two¹⁴ and three¹⁷ kinds of identical particles, and inclusion of particle conversion.¹⁶ Kato and Yamanouchi have extended the MCTDHF theory to molecules composed of electrons, (fermionic) protons, and two heavy (either classical or quantum) nuclei.¹⁸

In this paper, further stepping forward in this direction, we present a fully general TD-MCSCF method for a system comprising arbitrarily different kinds and numbers of interacting fermions and bosons (without particle conversion¹⁶). Treating all the constituent particles on an equal footing, we expand the total wave function in terms of configurations of the whole system [see eqn (5) below]. Thus, based on the time-dependent variational principle, we derive the equations of motion (EOM) of CI coefficients and spin-orbitals for general configuration spaces, not restricted to full-CI. As a simple example, we apply the TD-MCSCF method to high-harmonic generation (HHG) from a one-dimensional (1D) model hydrogen molecular ion H₂⁺ induced by an intense near-infrared (NIR) laser pulse, for which numerically exact solution is available.

This paper is organized as follows. Section 2 introduces our TD-MCSCF ansatz for many-particle systems composed of different kinds of fermions and bosons, and also defines the target Hamiltonian considered in this work. In Section 3, we derive the general equations of motion, based on the time-dependent variational principle. Explicit working equations for a molecule interacting with an external laser field are shown in Section 4. We present numerical examples in Section 5. Concluding remarks are given in Section 6.

2 Definition of the problem

2.1 TD-MCSCF ansatz

We consider a quantum mechanical many-body system with K kinds of fermions or bosons. The subsystem of kind α consists of N_α identical particles. Thus, there are particles as a whole. For notational brevity, we call such a system an N-particle system, where the array of integers N = (N₁N₂⋯N_K) carries the information of both particle kinds and the number of particles in each kind.

Let us define, for each kind of particles, the complete orthonormal set of spin-orbitals , which spans the one-particle Hilbert space Ω_α, is time-dependent in general. Then the N-particle Hilbert space is spanned by

where is a determinant (or permanent) of α-kind fermions (or bosons), consisting of N_α spin-orbitals chosen from . We call Φ_I(t) the I’s configuration, where I = I₁I₂⋯I_K is considered, at the moment, to collectively label the chosen spin-orbitals. The objective of this paper is to formulate the TD-MCSCF theory of the N-particle system within the ansatz of total wavefunction analogous to that for the electronic system, eqn (1), but using the configurations of eqn (2).

For rigorous and compact presentation of theory, we resort to the second quantization formulation by introducing creation and annihilation operators associated with . These operators obey the (anti-)commutation relations of bosons (fermions),

for bosons, where [â, b̂] = âb̂ − b̂â, andfor fermions, where {â, b̂} = âb̂ + b̂â.

Within the TD-MCSCF ansatz, the complete set of spin-orbitals is split into n_α (≥N_α for fermions) occupied spin-orbitals and remaining virtual spin-orbitals . We call the subspace of Ω_α spanned by occupied spin-orbitals the occupied spin-orbital space Ω^occ_α, and that spanned by virtual spin-orbitals the virtual spin-orbital space Ω^vir_α, where Ω_α = Ω^occ_α ⊕ Ω^vir_α. The total state Ψ(t) is expressed as a superposition of configurations Φ_I(t) of eqn (2), but constructed from occupied spin-orbitals only. Thus we write

where C_I(t) is the CI coefficient, and |I(t)〉 is the occupation number representation of the configuration Φ_I,

|I(t)〉 = |I₁(t)〉 ⊗ |I₂(t)〉 ⊗⋯⊗ |I_K(t)〉, (6)

Now I_α = I_α,1I_α,2⋯I_α,nα is (rigorously) reinterpreted as an integer array, satisfying . Note that I_α,iα ∈ {0,1} for fermions. Here and in what follows, we use indices i_α, j_α, k_α,… for occupied (Ω^occ_α), a_α, b_α, c_α,… for virtual (Ω^vir_α), and μ_α, ν_α, κ_α, τ_α,… for general (Ω_α) spin-orbitals of kind α. The indices p_α, q_α will be used for numbering the coordinates.

It should be noted that we do not restrict the expansion eqn (5) to the full-CI one. It should also be noticed that occupied configurations are specified in terms of the whole system rather than in terms of each particle kind as,¹⁷

with summation over I₁, I₂,…,I_K taken independently. Eqn (8) implies that given the set of configurations {I_α} for each kind, the set of configurations for the whole system {I} must include all the direct products of {{I_α};α = 1, 2,…,K}. Our approach based on eqn (5), including eqn (8) as a special case, allows a highly flexible choice of CI space, e.g., including up to double excitations,¹⁰ regardless of particle kind, from a reference configuration of the whole system, thereby enabling a proper account of correlation between different kinds of particles while reducing the computational cost.

2.2 Target Hamiltonian

In this article, we consider the Hamiltonian of an N-particle system composed of up to M-body terms,

H = H₁ + H₂ +⋯+ H_M, M ≤ N. (9)

The Hamiltonian is explicitly time-dependent in general, but the time argument t is dropped in this section for simplicity. Here, the m-body Hamiltonian is assumed to be given explicitly in terms of the coordinates (and momenta, see below) in a general sense characterizing m particles (or degrees of freedom), and symmetric under exchange of coordinates among particles of the same kind. One-particle Hamiltonian, e.g., is written aswhere the non-local form allows us to describe the momentum dependence of the Hamiltonian, and two-body interaction is generally given by¹⁶The reasons why we here consider the (non-local) higher-than-two body terms, which will not actually be used in Section 4, are (1) that such a form is used in the multiconfiguration Hartree (MCH) method for distinguishable particles, and (2) their possible appearance upon coordinate transformations, or in the effort of removing translational and rotational degrees of freedom.^19,20

The Hamiltonian is equivalently expressed in the second quantization formalism as

where the net m-body Hamiltonian is further classified into those contributions Ĥ_m, hereafter called the m-body Hamiltonian, involving m_α particles of the kind α, (0 ≤ m_α ≤ N_α, ),where μ = (μ₁μ₂⋯μ_K), and μ_α = (μ_α,1μ_α,2⋯μ_α,mα) indexes the set of spin-orbitals to represent m_α particles in the Hamiltonian. Ê^μ_ν is the m-particle replacement operator , withand (H_m)^m_ν is given bywhere x = (x₁x₂,…,x_K), x_α = (x_α,1x_α,2⋯x_α,mα) is the set of m_α coordinates of particle α, andFor the later discussion, we define the m-body reduced density matrix (RDM) as

(ρ_m)^μ_ν = 〈Ψ|Ê^ν_μ|Ψ〉. (17)

One- and two-particle RDMs are also denoted aswith β ≠ α.

3 Equations of motion

In this section, we derive the EOMs for the CI coefficients and spin-orbitals by imposing the time-dependent variational principle^21–23 on our TD-MCSCF ansatz. We require the action integralto be stationary, δS = 0, with respect to the variation of the total wavefunction δΨ within our TD-MCSCF ansatz eqn (5), subject to the boundary conditions δΨ(t₀) = δΨ(t₁) = 0. To this end, let us introduce anti-Hermitian matrices Δ_α and X_α as,(Recall that indices μ_α, ν_α refer to both occupied and virtual spin-orbitals.) We also define,with which orthonormality-conserving spin-orbital variations and time derivatives can be written asThen, the variation and time derivative of total state are compactly given by,^7,10,24and their Hermitian conjugates areIt follows from eqn (19) that,Substituting eqn (23) and (24) into this equation, after some algebraic manipulation,^7,24 we obtain,where denotes the projector onto the CI space, i.e., the subspace of N-electron Hilbert space spanned by the configurations included in eqn (5). The action functional S should be made stationary with respect to all independent variations; {δC_I, δC_I*} for CI coefficients and for spin-orbitals.

First, the EOM for CI coefficients are obtained from δS/δC_I* = 0,

Requiring δS/δC_I = 0 derives the complex conjugate of eqn (27). Next from , one obtainswhere [capital Pi, Greek, macron] = 1 − Π. Eqn (28) is to be solved for , thus determines the time derivative of spin-orbitals. We now take a closer look at eqn (28) for the following two distinct cases:

Case 1: (μ_α, ν_α) = (i_α, j_α)

In this case we focus on the components of the spin-orbital variations within the subspace spanned by the occupied spin-orbitals. Since and in general, one needs to directly work with eqn (28) within the occupied spin-orbital spaceIn the full-CI case, where , , eqn (29) reduces to an identity 0 = 0. Therefore, the corresponding may be arbitrary anti-Hermitian matrix elements, of which the simplest choice is .

Case 2: (μ_α, ν_α) = (i_α, a_α)

In this case we deal with the components of the spin-orbital variations outside the occupied spin-orbital space. Since and , eqn (28) becomes,However, the matrix element in the left-hand side of the above equation survives only when β = α and κ_α = a_α ∈ Ω^vir_α, namely . Thus eqn (30) is simplified toThe m-body Hamiltonian contribution to the RHS of eqn (31) is evaluated as follows;In the second line of the above equation, we note that the matrix element survives when one and only one of the m_α creation operators in refers to a_α ∈ Ω^vir_α, and all the others to the occupied spin-orbitals. All such cases [μ_α,p = a_α, μ_α,q≠p ∈ Ω^occ_α;1 ≤ p ≤ m_α] give the same contribution since the phase (∓)^p−1 [+ (−) sign for bosons (fermions), arising in (anti-)commuting the creation operators] is canceled by shifting the corresponding annihilation operator ν_α,p, and the Hamiltonian is symmetric for interchange of particles of the same kind. The third line is thus obtained after renaming summation variables, where k^[α] = (k₁,…,k_α,…,k_K) is the array of m − 1 indices with k_α = (k_α,2,…,k_α,mα) and k_β = (k_β,1,…,k_β,mβ) for β ≠ α, with l^[α] defined similarly. The fourth line introduces the short-hand notation for the array of m indices, μ_αk^[α] = (k₁,…,μ_αk_α,…,k_K) with μ_αk_α = (μ_α,k_α,2,…,k_α,mα) (μ_αl^[α] is defined similarly), and the fifth line uses the definition of the m-body RDM, eqn (17).

Now the RHS of eqn (31) is given by the sum over m,

where is the effective one-particle operator,and is given in the coordinate representation aswhere y^[α] = (y₁,…,y_α,…,y_K) is the set of m − 1 coordinates with y_α = (y_α,2,…,y_α,mα) and y_β = (y_β,1,…,y_β,mβ) for β ≠ α, and x_αy^[α] = (y₁,…,x_αy_α,…,y_K) is the array of m coordinates. z^[α] and x_α′z^[α] are defined similarly.

Finally, gathering the occupied and virtual components of the time derivative completes the derivation of the EOM for spin-orbitals

where is the spin-orbital projection operator onto the occupied spin-orbital space, with which the virtual space is referenced as a whole, thus avoiding the explicit use of virtual spin-orbitals. in the first term is to be obtained by solving eqn (29), and, as discussed above, can be set zero in the full-CI case. Eqn (27) for CI coefficients and eqn (36) for spin-orbitals form fully general TD-MCSCF equations of motion, not restricted to full CI, for a system composed of any arbitrary kinds and numbers of fermions and bosons.

4 Molecules interacting with an external laser field

In this section we present the working equations for a molecule subject to an external laser field. Let the molecule consist of electrons and K_n different kinds of nuclei treated quantum mechanically (the kind does not necessarily corresponds to the nuclear species, see discussion below), and N_cl nuclei treated as a classical point charge. For clarity and notational simplicity, we assign the electrons to the first kind of particle (α = 1), and kinds α = 2, 3,…,K represent quantum nuclei with K = 1 + K_n. The numbers of identical particles are, as before, denoted by {N_α}. Then the number of electrons is N₁, the number of quantum nuclei is , and the total number of atoms is N_atom = N_n + N_cl. We use atomic units in this section.

The spin-independent molecular Hamiltonian in the coordinate representation is given by

where U_αβ(r) = Z_αZ_β/r is the Coulomb interaction with Z_α being the electric charge, andis the one-particle Hamiltonian composed of the kinetic energy [the first term with m_α being the mass (not to be confused with the number of particles)], Coulomb interaction with classical nuclei with the charges {Z_A} located at {R_A} (the second term), and the time-dependent laser-particle interaction V^ext_α, given, e.g., within the dipole approximation either in the length gauge (LG) or in the velocity gauge (VG), by

V^ext_α,LG(r,r′,t) = −δ(r − r′)Z_αE(t)r, (39)

where E(t) is the laser electric field, and is the vector potential.

The general formulation of Section 3 is readily applicable to the molecular Hamiltonian of eqn (37). The CI EOM reads

wherewith x_α = (r_α,σ_α) being the composite spatial- and spin-coordinates, and the EOM for spin-orbitals is given bywhere the one-body contribution to the second term of eqn (36) is extracted to lead to by noting , and

Finally, eqn (29) is formulated as the linear system of equations,

where , , withIn order for eqn (47) to be solvable (with non-singular coefficient matrix A), one needs a systematic method of constructing non-full-CI space analogous to the TD-ORMAS method¹⁰ for electrons. We shall discuss this issue in the future publication.

Equations of motions (41) and (44), with the matrix eqn (47), define the general TD-MCSCF method, not restricted to full CI, for molecules interacting with an external field. Our formulation is very flexible; it includes as special cases both the electron dynamics at the classical-nuclei approximation (N_n = 0, N_cl = N_atom) and the full quantum molecular dynamics (N_n = N_atom, N_cl = 0). Furthermore, it allows various approaches to the same physical problem; e.g., the same nuclear species in the molecule can be treated either as identical particles or distinguishable ones to investigate the physical outcomes of the particle statistics during the course of laser-molecule interactions.

5 Numerical examples

In order to demonstrate the utility of the TD-MCSCF method for an exactly solvable system, in this section, we consider a one-dimensional model hydrogen molecular ion with the singlet protonic spin configuration (1D para-H₂⁺), driven by an intense near-infrared laser pulse. The Hamiltonian within the dipole approximation in the velocity gauge is given in terms of electron coordinate x and proton coordinates X₁, X₂ bywhere M denotes the proton mass, and , and are the electron–proton and proton–proton soft Coulomb potentials, respectively. Following ref. 25, we set c_ne = 1 and c_nn = 0.03, and transform the coordinates as z = x − (X₁ + X₂)/2, X = (X₁ + X₂)/2, and R = X₁ − X₂. Neglecting the terms involving X, ∂/∂X and transforming into the length gauge lead to the reduced Hamiltonian,²⁵with μ_e = 2M/(2M + 1) and q_eff = −(2M + 2)/(2M + 1). This Hamiltonian was used, e.g., for investigating charge-resonance enhanced ionization,²⁵ but application of the multiconfiguration method to this system in the context of highly nonlinear laser-driven dynamics such as high-harmonic generation has not previously been reported to the best of our knowledge.

Let us compare the following treatments: (1) numerically exact solution of the time-dependent Schrödinger equation (TDSE) for the spatial part Ψ_exact(z,R,t) of the 1D × 1D wavefunction by discretizing electron coordinate z and nuclear distance R with constant grid spacing δz = 0.4 and δR = 0.2, respectively, within the simulation volume |z| < 640 and |R| < 80, (2) the electron only dynamics by solving for electronic wavefunction Ψ_fixed(z,t) given the first line of eqn (51) with a fixed (near to equilibrium) internuclear distance R₀ = 2.6 [fixed (classical point) nuclei approximation], and (3) TD-MCSCF method with the total wavefunction given by

using n occupied spatial orbitals both for the electron and the R degree-of-freedom of the protons. Note that n = 1 corresponds, in this example, to the time-dependent Hartree approximation, Ψ_TDH = Ψ₁ = χ^e₁(z,t)χⁿ₁(R,t). (4) Based on this TDH method we also consider the frozen (quantum) nuclei approximation, where the protonic orbital is fixed to the initial form,

Ψ_frozen(z,R,t) = χ^e₁(z,t)χⁿ₁(R,0). (53)

The same spatial discretization as that in the first treatment is applied to the latter three. For each method described above we first obtain the ground state by imaginary time propagation of the respective EOM, and then investigate the molecular dynamics induced by the NIR laser pulse with a peak intensity of 10¹⁴ W cm⁻² and a central wavelength of 800 nm, with a sin² envelop of the foot-to-foot pulse duration of 12 optical cycles. The TDH ground state orbitals are used as initial orbitals in eqn (53).

Fig. 1 shows the time evolution of the dipole acceleration a(t) = d²〈Ψ(t)|z|Ψ(t)〉/dt² computed using the Ehrenfest theorem. As shown in Fig. 1(a), both fixed-nuclei (light blue curve) and frozen-nuclei (red) approximations strongly underestimate the nonlinear response seen in the TDSE result (black), which suggests the importance of the dynamical quantum treatment of nuclei. The TD-MCSCF method with n = 1, or TDH [Fig. 1(b)] provides improved description of the nonlinear response, but still fails to quantitatively reproduce the TDSE result. The TD-MCSCF method with n = 8 [Fig. 1(c)], on the other hand, produces the dipole acceleration that agrees with the TDSE result nearly perfectly, on the scale of the figure. The HHG spectra calculated as the modulus squared of the Fourier transform of the dipole acceleration are shown in Fig. 2. The fixed and frozen nuclei approximations underestimate the intensity of the first few harmonics, while overestimating the higher plateau (the cutoff position is estimated to be the 39th harmonic from the Lewenstein model²⁶) [Fig. 2(a)], both suggesting an important contribution from (more polarizable and loosely bound) the larger |R| region in the TDSE result. The TD-MCSCF spectrum with increasing n [shown for n = 1 (Fig. 2(b)) and n = 8 (Fig. 2(c))] shows increasingly better agreement with the TDSE one, and the TD-MCSCF method with n = 8 nicely reproduces the overall spectrum obtained using the TDSE simulations, especially for low-order harmonics, albeit with a still remaining slight overestimation of the plateau intensity which implies a strong electron–proton correlation.

Fig. 1 Time evolution of the dipole acceleration of one-dimensional H₂⁺ exposed to a laser pulse with a wavelength of 800 nm and an intensity of 1.0 × 10¹⁴ W cm⁻². Comparison of the results of fixed nuclei and frozen nuclei approximations (a), TD-MCSCF method with n = 1 (b), and TD-MCSCF method with n = 8 (c) with numerically exact TDSE results.

Fig. 2 The HHG spectra of one-dimensional H₂⁺ for the same laser pulse and simulation methods as those in Fig. 1. Comparison of the results of fixed nuclei and frozen nuclei approximations (a), TD-MCSCF method with n = 1 (b), and TD-MCSCF method with n = 8 (c) with numerically exact TDSE results. Also shown in (b) is the result of the TD-MCSCF method directly applied for the Hamiltonian of eqn (50) with one electronic orbital and two protonic orbitals.

Finally, also plotted in Fig. 2(b) with a red line is the result of TD-MCSCF method directly applied to the Hamiltonian of the original coordinate, eqn (50), with one electronic orbital and two protonic orbitals to expand the three-particle wavefunction Ψ(x,X₁,X₂), using the method described in ref. 19 to eliminate the translational degree of freedom. We find a remarkable agreement between the results from the original Hamiltonian eqn (50) and the reduced Hamiltonian eqn (51). The former approach can be applied to more complex systems in general.

6 Summary

We have developed a fully general ab initio TD-MCSCF approach to describe the dynamics of a many-body system that is a mixture of any arbitrary kinds and numbers of fermions and bosons subject to an external field. In this approach, the total wave function is expanded in terms of configurations constructed from time-dependent single-particle spin-orbitals. The expansion is not limited to the full-CI one, and the configurations used in the expansion can be specified in terms of the whole mixture. The equations of motion for the CI coefficients and spin-orbitals have been derived, based on the time-dependent variational principle. Furthermore, we have presented the working equations applicable to the investigation of the ultrafast dynamics in a molecule irradiated by intense laser fields and/or ultrashort XUV pulses.

The present framework is highly flexible. For example, we can treat identical nuclei in spatially separated subdomains of a molecule as different particle kinds. We can also treat heavy nuclei and incident projectiles as classical particles¹⁸ instead of quantum ones. The latter may be handled as an external field as well. This flexibility will be useful for unraveling the physical mechanisms underlying the phenomena under investigation.²⁷ As a future prospect, an introduction of a CI space that adaptively changes following the dynamics as well as an inclusion of particle conversion¹⁶ would enable even more efficient and flexible simulations.

Whereas our original motivation lies in ab initio simulations of the electron–nuclear dynamics in molecules driven by a laser pulse, our method will be applicable to a wide variety of problems far beyond. Especially, the Hamiltonian can contain non-local terms and involve many-body (more than two-body) interactions. Thus, it may also find applications in cold-atom/cold-molecule physics and nuclear physics.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We thank Joachim Burgdörfer for helpful discussions. This research was supported in part by a Grant-in-Aid for Scientific Research (Grants No. 25286064, 26390076, 26600111, 16H03881, and 17K05070) from Japan Society for the Promotion of Science and also by the Photon Frontier Network Program of the Ministry of Education, Culture, Sports, Science and Technology of Japan. This research was also partially supported by the Center of Innovation Program from the Japan Science and Technology Agency, JST, and by CREST (Grant No. JPMJCR15N1), JST. R. A. gratefully acknowledges support from the Graduate School of Engineering, The University of Tokyo, Doctoral Student Special Incentives Program (SEUT Fellowship).

References

1. P. Agostini and L. F. DiMauro, Rep. Prog. Phys., 2004, 67, 1563.
2. F. Krausz and M. Ivanov, Rev. Mod. Phys., 2009, 81, 163–234.
3. L. Gallmann, C. Cirelli and U. Keller, Annu. Rev. Phys. Chem., 2013, 63, 447.
4. T. Kato and H. Kono, Chem. Phys. Lett., 2004, 392, 533–540.
5. J. Caillat, J. Zanghellini, M. Kitzler, O. Koch, W. Kreuzer and A. Scrinzi, Phys. Rev. A: At., Mol., Opt. Phys., 2005, 71, 012712.
6. K. L. Ishikawa and T. Sato, IEEE J. Sel. Top. Quantum Electron., 2015, 21, 1–16.
7. T. Sato and K. L. Ishikawa, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 023402.
8. T. Sato, K. L. Ishikawa, I. Březinová, F. Lackner, S. Nagele and J. Burgdörfer, Phys. Rev. A, 2016, 94, 023405.
9. H. Miyagi and L. B. Madsen, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 063416.
10. T. Sato and K. L. Ishikawa, Phys. Rev. A: At., Mol., Opt. Phys., 2015, 91, 023417.
11. H.-D. Meyer, U. Manthe and L. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
12. I. S. Ulusoy and M. Nest, J. Chem. Phys., 2012, 136, 054112.
13. M. Nest, Chem. Phys. Lett., 2009, 472, 171–174.
14. O. E. Alon, A. I. Streltsov and L. S. Cederbaum, Phys. Rev. A, 2007, 76, 062501.
15. O. E. Alon, A. I. Streltsov and L. S. Cederbaum, J. Chem. Phys., 2007, 127, 154103.
16. O. E. Alon, A. I. Streltsov and L. S. Cederbaum, Phys. Rev. A: At., Mol., Opt. Phys., 2009, 79, 022503.
17. O. E. Alon, A. I. Streltsov, K. Sakmann, A. U. Lode, J. Grond and L. S. Cederbaum, Chem. Phys., 2012, 401, 2–14.
18. T. Kato and K. Yamanouchi, J. Chem. Phys., 2009, 131, 164118.
19. H. Nakai, M. Hoshino, K. Miyamoto and S. Hyodo, J. Chem. Phys., 2005, 122, 164101.
20. H. Nakai, M. Hoshino, K. Miyamoto and S. aki Hyodo, J. Chem. Phys., 2005, 123, 237102.
21. A. Dalgarno and G. Victor, Proc. R. Soc. London, Ser. A, 1966, 291–295.
22. P.-O. Löwdin and P. Mukherjee, Chem. Phys. Lett., 1972, 14, 1–7.
23. R. Moccia, Int. J. Quantum Chem., 1973, 7, 779–783.
24. R. P. Miranda, A. J. Fisher, L. Stella and A. P. Horsfield, J. Chem. Phys., 2011, 134, 244101.
25. E. Khosravi, A. Abedi and N. T. Maitra, Phys. Rev. Lett., 2015, 115, 263002.
26. M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L'Huillier and P. B. Corkum, Phys. Rev. A: At., Mol., Opt. Phys., 1994, 49, 2117.
27. I. Tikhomirov, T. Sato and K. L. Ishikawa, Phys. Rev. Lett., 2017, 118, 203202.

---