Time-dependent quantum chemistry of laser driven many-electron molecules.

Time-dependent quantum chemistry of laser driven many-electron molecules Thanh-Tung Nguyen-Dang, Étienne Couture-Bienvenue, Jérémy Viau-Trudel, and Amaury Sainjon Citation: The Journal of Chemical Physics 141, 244116 (2014); doi: 10.1063/1.4904102 View online: http://dx.doi.org/10.1063/1.4904102 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/24?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Explicitly time-dependent coupled cluster singles doubles calculations of laser-driven many-electron dynamics J. Chem. Phys. 134, 054113 (2011); 10.1063/1.3530807 Time-dependent multiconfiguration theory for electronic dynamics of molecules in intense laser fields: A description in terms of numerical orbital functions J. Chem. Phys. 128, 184102 (2008); 10.1063/1.2912066 Photoionization-induced dynamics of ammonia: Ab initio potential energy surfaces and time-dependent wave packet calculations for the ammonia cation J. Chem. Phys. 124, 214306 (2006); 10.1063/1.2202316 Time-dependent configuration-interaction calculations of laser-pulse-driven many-electron dynamics: Controlled dipole switching in lithium cyanide J. Chem. Phys. 123, 074105 (2005); 10.1063/1.1999636 Electron propagator theory calculations of molecular photoionization cross sections: The first-row hydrides J. Chem. Phys. 121, 4143 (2004); 10.1063/1.1773135

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.237.165.40 On: Sun, 30 Aug 2015 06:26:43

THE JOURNAL OF CHEMICAL PHYSICS 141, 244116 (2014)

Time-dependent quantum chemistry of laser driven many-electron molecules Thanh-Tung Nguyen-Dang, Étienne Couture-Bienvenue, Jérémy Viau-Trudel, and Amaury Sainjon Département de Chimie, Université Laval, Québec, Québec G1V 0A6, Canada

(Received 25 September 2014; accepted 1 December 2014; published online 31 December 2014) A Time-Dependent Configuration Interaction approach using multiple Feshbach partitionings, corresponding to multiple ionization stages of a laser-driven molecule, has recently been proposed [T.-T. Nguyen-Dang and J. Viau-Trudel, J. Chem. Phys. 139, 244102 (2013)]. To complete this development toward a fully ab-initio method for the calculation of time-dependent electronic wavefunctions of an N-electron molecule, we describe how tools of multiconfiguration quantum chemistry such as the management of the configuration expansion space using Graphical Unitary Group Approach concepts can be profitably adapted to the new context, that of time-resolved electronic dynamics, as opposed to stationary electronic structure. The method is applied to calculate the detailed, sub-cycle electronic dynamics of BeH2, treated in a 3–21G bound-orbital basis augmented by a set of orthogonalized plane-waves representing continuum-type orbitals, including its ionization under an intense λ = 800 nm or λ = 80 nm continuous-wave laser field. The dynamics is strongly non-linear at the field-intensity considered (I ≃ 1015 W/cm2), featuring important ionization of an inner-shell electron and strong post-ionization bound-electron dynamics. C 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4904102]

I. INTRODUCTION

Solving the time-dependent Schrödinger equation (TDSE) for a many-electron atom or molecule driven by an external time-dependent force field remains an important challenge in computational physics and chemistry. Two major difficulties are encountered as one attempts to describe accurately this laser-driven electron dynamics by numerical simulations. The dimensionality of the problem is one of these difficulties. Grid-based methods, for example, are limited to systems of low or reduced dimensionality,1,2 as the increased number of electrons, mutually in interaction, described in full threedimensional space, implies memory and CPU requirements that rapidly become prohibitive. Large-amplitude electron motion in strong-field ionization is another one: In grid-based methods, it implies large grids, while in approaches using basis functions, such as B-Splines,3,4 Sturmians,5 or those using L 2 Gaussian primitives, and adapted from stationary-state abinitio quantum chemistry,6–10 this large-amplitude electronic excursion calls for extensive use of polarization functions and basis functions not centered at the actual nuclei of the molecule, and designed to cover the expected large spatial extension of the electronic wavepackets. In a preceding paper,11 partitioning ideas have been proposed to address this difficulty within a quantum chemical Time-Dependent Configuration Interaction (TDCI) type approach, toward the development of a fully ab-initio method for the calculations of N-electron wavepackets, describing, in a time-resolved manner, the electronic dynamics of a molecule in an intense laser-field, including multiple ionizations. The ideas begin with realizing that, in both grid-based and 0021-9606/2014/141(24)/244116/17/$30.00

functions-based approaches, the computation costs become prohibitive only if one maintains the same level of treatment of electron correlation throughout, i.e., for ionized as well as for bound electrons. On the other hand, a qualitative understanding of attosecond electron dynamics has been made possible only thanks to a number of approximations. Among these, the strong-field approximation (SFA) considers that any ionized electron would feel the external field only, and not the Coulomb forces of the ionic core, while for bound electrons, the latter are the dominant forces.12,13 The idea of a separate treatment of bound and free electrons, suggested by this strong-field approximation, has been exploited in grid-based schemes using R-matrix theory,14–16 and a partitioning of real, electron-position (coordinate) space into internal and external regions, corresponding to bound and ionized states of an electron. In basis-function-based descriptions, such as the one used here, it is a Feshbach partitioning of the many-electron Hilbert space (state space) and Hamiltonian that is required to express the same idea.11,17–21 A unitary solution to the partitioned TDSE is obtained in Ref. 11 first by regrouping the various components of the Feshbach-partitioned Hamiltonian into (block-)diagonal and non-diagonal parts, a procedure called Adams repartitioning.22 The total time-evolution operator is then factorized into a product of corresponding propagators, the one associated with the non-diagonal bound-free interaction Hamiltonian in this repartitioning being represented, in each elementary (short) timeinterval, by a Cayley-Crank-Nicholson type propagator.23,24 Extension to many successive ionization stages, as generically appears in any many-electron problem, makes use of successive two-subspace partitionings of the many-electron state

141, 244116-1

© 2014 AIP Publishing LLC


244116-2

Nguyen-Dang et al.

space. The partitioning scheme defines orthogonal subspaces of the N-electron Hilbert space corresponding to different ionization stages. Typically, the dynamics within the bound neutral molecule and cation can be treated in a fully correlated (for example, TDCI) level, while correlation effects involving the freed electrons can be treated separately, at a mean-field level for instance. These partitioning ideas can be expressed in any Nelectron basis. Our preference went to the basis of the socalled configuration-state functions (CSF) as these are readily constructed and thus known, as soon as the molecular orbital basis is given, whereas large-scale, high-level multiconfiguration electronic-structure calculations, such as found in Refs. 25 and 26 would be needed in a preparatory step if eigenstates of the field-free Hamiltonian are used instead. The multiple Feshbach partitioning operators, expressed in the N-electron CSF basis, are recalled in Sec. II, using a second-quantized occupation-number representation of the CSF and restricting to the first two ionization steps, which are the most important and pertinent ones in current research in attosecond physics. In the previous paper,11 the methodology was illustrated on a model diatomic molecule with just two active electrons distributed among two bound molecular orbitals, and a basis of orthogonalized plane-waves. The expressions of various matrices appearing in the implementation of the propagation scheme are relatively simple in this two-electron, two-boundmolecular orbital (MO) model. How are these matrices expressed and systematically calculated in a general N-electron, Norb-orbital system (N, Norb > 2) is the type of question we will address in Sec. II, using a four-electron, five-(active)orbital singlet model system for explicit illustrations. It is in this more general case that many tools of ab-initio quantum chemistry take on their full significance and will prove to be most useful for time-dependent electron dynamics. Among these, is the so-called direct-row table (DRT),27,28 a concept derived from the Graphical Unitary Group Approach (GUGA),29,30 which was used in the construction, identification, and handling of the very large number of CSF encountered in large CI (configuration interaction) and MCSCF (multi-configuration self-consistent-field) calculations. These GUGA concepts are reviewed in Appendix A. Normally, if the N-electron system is treated in a full TDCI approach including all (singlet) CSF that one can construct out of an extended orbital basis containing both bound and a large number of continuum-type orbitals, then the number of CSF can become so large, and the DRT so long that the problem would become intractable. Would this dimensionality problem be overcome or avoided with the Feshbach partitioning scheme alluded above? By considering matrix elements describing dipole transitions between the bound-state subspace and subspaces associated with one- and two-electron ionizations, we found that a reduced DRT is sufficient, in which at most two continuum-type orbitals are required to be added to the finite set of bound orbitals, to represent symbolically all the oneand two-electron ionization continua. This kind of reduction, illustrated in Sec. II, and demonstrated in more details in Appendix B, has been noted and used in an earlier work dealing with photoionization treated in a time-independent scattering theoretical context.31

J. Chem. Phys. 141, 244116 (2014)

We illustrate, in Sec. II, how general propagation equations for the multicomponent, many-electron wavepacket (given in full in Appendix B, for single and double ionizations), together with GUGA concepts, are explicitly implemented in numerical calculations of the (single) ionization of the six-electron BeH2 molecule in an intense laser field. The calculations used a rather modest (3–21G) bound orbital basis, augmented by a basis of plane-waves defined in a two-dimensional reciprocal space. With the molecular core orbital considered frozen (we have thus 4 active electrons instead of 6), and an active space of five bound molecular orbitals, the model comprises 50 bound CSF for the neutral parent molecule and 40 for the singly ionized cation. The results of calculations on this model system, in particular channel-resolved ionization data, defined in Sec. III, are given and discussed in Sec. IV, for a continuous-wave (CW), λ = 800 nm and λ = 80 nm, field polarized either parallel or perpendicular to the internuclear axis of the linear molecule. Very detailed dynamical effects observed there are given a thorough interpretation using selection rules and strong-field non-linear mixing mechanisms. Some non-linear processes are seen to favor states that are normally inaccessible by a one-photon dipole transition, over doorway states that are directly dipole-coupled to the initial state. This predominance even leads to a dramatic quenching of these doorway states in the high-frequency (λ = 80 nm) case. Another non-linear effect found in this study concerns the ionization of the molecule under the λ = 800 nm field. Removal of an innershell electron is favored over that of the valence one and the ionization to the corresponding ionic channel indicates the involvement of many competing pathways. The simulations on this simple four-active-electron molecule show important excitations of bound electrons accompanying the ionization processes, a dimension obviously absent in simulations using the so-called single-active electron (SAE) approximation,32 and which is currently of much interest.33–35

II. MANY-ELECTRON LASER-DRIVEN DYNAMICS AND AB-INITIO METHODS A. Second-quantized Hamiltonian

ˆ In a fixed, finite MO basis, {φ s }, the Hamiltonian, H(t), of an N-electron molecular system driven by an external laser field can be expressed in terms of second-quantized operators associated with these orbitals in the form29,36 ˆ H(t) = Hˆ 1(t) + Hˆ 2,  Hˆ 1(t) = hr s (t) Eˆr s , r

(1a) (1b)

s

1  gr su v eˆr su v . Hˆ 2(t) = 2 r s u v

(1c)

In the one-electron part Hˆ 1, ˆ r ,t)|φ s (⃗r )⟩ hr s (t) = ⟨φr (⃗r )| h(⃗

(2a)

are one-electron integrals, where


244116-3

Nguyen-Dang et al.

ˆ r i ,t) = h(⃗

pî2  −Zα ⃗ i) + + Vînt (⃗r i ,t), (⃗ pi = −i ∇ 2 ⃗α | ri − R α |⃗

J. Chem. Phys. 141, 244116 (2014)

(2b)

contains the field-electron, spin-conserving interaction poten⃗ ⃗ = E(t) ⃗ tial Vînt (⃗r ,t) = d⃗ · F(t), with d⃗ = ⃗r , F(t) (the electric ⃗ ⃗ ⃗ (the vector field) in the length gauge, and d = p⃗, F(t) = A(t) potential) in the velocity gauge.37,41 The operator Eˆr s = aˆr†α aˆ sα + aˆr† β aˆ s β

(3)

is called one-electron generator. It describes the excitation of an electron (in any spin state) from a first orbital φ s to a second orbital φr , when r , s, while Eˆr r is the rth orbital occupation number operator. In the two-electron part Hˆ 2 of the Hamiltonian, 1 |φ s (⃗r 1)φ v (⃗r 2)⟩ (4) r 12 are two-electron integrals, and the two-electron generators are gr su v (t) = ⟨φr (⃗r 1)φu (⃗r 2)|

eˆr su v = Eˆr s Eû v − Eˆr v δ su .

(5)

† The fermion creation and annihilation operators aˆrν , aˆrν associated with the spin orbital φrν , ν = α, β (corresponding to m s = ±1/2) satisfy well-known anticommutation relations, which imply the following commutation properties of the Eˆr s :

[ Eˆr s , Eû v ] = Eˆr v δu s − Eû s δr v .

(6)

FIG. 1. Shavitt graph (see Appendix A for definition) representing the DRT of the Be H2 model used in the calculations [panel (a)]. The numbers on the leftmost vertical column refer to active orbitals, with 0 denoting the vacuum state. The graph represents all the singlet CSF that can be obtained by distributing four electrons in n ac t = 5 active bound orbitals, (arcs in black thick lines), and allowing a single occupation of a symbolic continuum orbital, (MO 6). Walks leading to this orbital and terminating with arcs traced in red thick lines represent singlet CSF of the Be H2++e system, (|J +, κ⟩ in the text). Walks following arcs traced in thin lines represent the same situation, but with the ionized electron ending in another continuum-type orbital. They are never used in the calculations of matrix elements of one- and two-electron operators. The active bound MO are 1 ↔ 2σ g , 2 ↔ 1σ u , 3 ↔ 1π u(x), 4 ↔ 1π u(y), and 5 ↔ 3σ g . Their relative disposition on an orbital energy scale together with their shape and nodal structure is shown in panel (b). Note that the 1σ u is frozen, and does not appear in the DRT, or Shavitt graph of panel (a).

B. Feshbach partitionings

The orbital basis {φ s } generally comprises bound MOs, νb of which are considered active, i.e., their populations may change during the laser-driven dynamics, and are designated by ϕi , i = 1,2,...,νb , and νc orbitals representing a discretized one-electron continuum, χκ , κ = 1,2,...,νc , preorthogonalized with respect to the ϕi ’s, so that νb νc {φr } = {ϕi }i=1 ∪ { χκ }κ=1 ,

corresponding to a partition of the one-electron Hilbert space ν b H (1) by orthogonal projection operators qˆ = i=1 |ϕi ⟩⟨ϕi |,  | χκ ⟩⟨ χκ | = 1 − q. ˆ pˆ = κ

For the specific case of the BeH2 molecule considered in Secs. III and IV, ϕi are field-free MOs calculated (using Columbus42), at the HF-SCF level using the 3–21G basis, in its singlet ground-state and at its ground-state equilibrium linear geometry. Figure 1(b) shows the six first MOs of the molecule calculated at this level. The orbitals are designated by the usual notations for symmetric linear molecule (exact point group D∞h , although the calculations were made in C1). The core orbital 1σg ≃ 1s Be is considered frozen, and active orbitals are the HOMO (highest-occupied molecular orbital, 1σu ), HOMO-1 (2σg ), the LUMO (lowest unoccupied MO, 1πu = 2p⊥Be , those are in fact doubly degenerate, one is along the x axis, another is in the y direction), LUMO+1, (3σg ), giving 5 active MOs and 4 active electrons. The time-dependent state of the N-electron molecule can be expressed in terms of a basis of CSF. These are Nelectron spin-adapted antisymmetrized states corresponding to a definite occupation scheme of the MO, with a given spin

multiplicity defined by a given value of the total spin quantum number S. A CSF of a given multiplicity S, denoted |ΦSI ⟩ = |n1,n2,n3,...,nν b ,···nν b +k ···nν b +ν c ⟩      nB

(7)

nC

= |n B;nC ⟩, where nr is the occupation number of orbital φr which is a linear combination of Slater determinants |φ1α φ1 β φ2α φ2 β ···| = |n1α ,n1β ,n2α ,n2β ,...⟩ † n 1α † n 1β † n 2α † n 2β = (aˆ1α ) (aˆ1β ) (aˆ2α ) (aˆ2β ) ···|0⟩,

(|0⟩ designates the vacuum state in occupation number representation), with nr = nr α + nr β , and appropriate spincoupling coefficients to produce an eigenvector of the total spin operator Sˆ 2 with eigenvalue S(S + 1). The partition of the one-electron Hilbert space by qˆ and pˆ implies a hierarchy of partitions of the N-electron state space, as shown in Ref. 11. Restricting to two ionization steps, and using the present notations, we thus define   |n B;0C ⟩⟨n B;0C | = |I⟩⟨I| (8) Qˆ 0 = nB

I

as the projector onto N-electron bound-states (i.e., zeroionized-electron states, as indicated by nC = 0C in the CSF used in Eq. (8)). Its complement, Pˆ0 = 1 − Qˆ 0, is then the projector onto states with at least one ionized electron. In


244116-4

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

turn, it is decomposable into Pˆ0 = Qˆ 1 + Pˆ1 = Qˆ 1 + Qˆ 2 + Pˆ2 = ···,   Qˆ 1 = |n B;1κ ⟩⟨n B;1κ | = Qˆ 2 =

nB

κ

I+

κ

    nB

=

(9a)

|n B;1κ ,1κ ′⟩⟨n B;1κ ,1κ ′|

κ,κ ′

  I ++

|I +,κ⟩⟨I +,κ|,

|I ++,κ,κ ′⟩⟨I ++,κ,κ ′|.

(9b)

κ,κ ′

As defined, Qˆ n , (n = 0,1,2), denotes the projector onto states with exactly n ionized electrons, and Pˆn projects onto states with at least (n + 1) ionized electrons. The corresponding Nelectron subspace will be referred to by the same name as the associated projector, without the caret. The dimension of a Q i subspace will henceforth be called NQ i . An alternative system of notations is evoked in Eqs. (8)-(9b) and uses an index I l+ to identify bound CFSs of the neutral molecule (l = 0, not explicitized) and singly and doubly charged cations (l = 1,2). This index refers to a lexical ordering of the CSF in their identification with distinct paths (walks) in a Shavitt graph, the graphic representation of the so-called DRT of the GUGA,27–30 which has largely been used in multiconfiguration ab-initio theories of quantum chemistry.36 A review of the concepts of DRT and Shavitt graph, and of the way CSF are defined and used in this GUGA representation, is given in Appendix A.

2. Coupling between subspaces

The matrices HQ 0Q 1(= H†Q 1Q 0), L±0 , represent, in the CSF ± basis, the operators H¯ Q 0Q 1(= H¯ Q 1Q 0) and Lˆ Q of Eqs. (B4) 0 and (B3), respectively. As shown in Appendix B, if the couplings between subspaces are restricted to the radiative ones, which are contained in the one-electron part Hˆ 1 of the Hamiltonian, then the calculation of the rectangular matrix H¯ Q 0Q 1 requires the knowledge of the one-electron integrals t ˆ ′)dt ′, and that of the hrκ = ⟨ϕr |h(t n )| χκ ⟩, h(t n ) = t nn+1 h(t matrix element of a one-electron generator Eˆr κ for a single value of κ, κ = νb + 1, say 

 HQ 0 Q 1

I, {J +,κ}

=

νb 

( ) hrκ ⟨I| Eˆr,κ | J +,κ⟩ . κ=ν b +1 r =1

(12)

Thus, the mapping, induced by Hˆ 1, of aQ1 component of the   N-electron wavefunction (Qˆ 1|Ψ⟩ = J + κ γJ +,κ | J +,κ⟩), onto a vector of the Q0 subspace, does not use matrix elements involving all the CSF | J +,κ⟩, with κ covering its full range of values. It can be carried out in two separate steps (i) a sum over the cation CSF index J + only, for each orbital pair (r,κ), yielding ) (  , (13a) γJ +,κ ⟨I| Eˆr,κ | J +,κ⟩ m I (r,κ) = κ=ν b +1 J+ followed by (ii) a convolution of the result with the continuum → bound orbital transformation matrix element hr κ (variable κ in this one-electron integral)     hrκ .m I (r,κ). (13b) HQ 0 Q 1 γ = I

r

κ

C. N-electron wavepacket propagation 1. Propagation equations

How the TDSE is to be solved under this multiple partitioning scheme has been shown in Ref. 11. The working equations for the calculation, over a short time-slice, t ∈ [t n ,t n+1], of the time-dependent N-electron wavefunction |Ψ(t)⟩, including two ionization steps are recalled in full in Appendix B. We give here the shorter version for the propagation in time of the CI wavefunction (singlet state S = 0) of the BeH2 including only the first ionization of the molecule    |Ψ,t⟩ = cI (t)|I⟩ + γJ +,κ | J +,κ⟩, (10) I

J+ κ

for which, the general algorithm of Eqs. (B1a)-(B1c), rewritten in terms of the vectors c(t), γ(t) of CI coefficients cI and γJ +,κ reduces to  c I (t) = UQ 0Q 0 L−0 [L+0 ]−1c I (t n−1)  −i[L+0 ]−1HQ 0Q 1γ(t n−1) , (11a)   (1 γ(t) = UQ 1Q 1 γ(t n−1) − HQ 1Q 0 [L+0 ]−1HQ 0Q 1γ(t n−1) 2 ) +i[L+0 ]−1c I (t n−1) , (11b) for t ∈ [t n−1,t n ].

3. Reduced DRT and Shavitt graph in GUGA

The fact that the Q0 ↔ Q1 transitions require, for each bound orbital index r, only one (one-electron) generator Er κ associated with a unique, arbitrary continuum orbital χκ means that, in constructing the DRT or the Shavitt graph for the system, the ionization continuum can be represented by but a single, arbitrary (symbolic) orbital, which was designated above by (κ = νb + 1). This kind of reduction has been noted and used in an earlier work dealing with photoionization treated in a time-independent scattering theoretical context.31 Figure 1, panel (a), illustrates this reduction for the above model of BeH2, with four active electrons and νb = 5 active bound orbitals, showing also how the Shavitt graph would look like if νc > 1 continuum MOs are included. Within the lexical ordering convention defined in Ref. 30 the CSF (the walks in the Shavitt graph of Figure 1, panel (a)), are numbered by an index m, such that m ≤ NQ 0 corresponds to the CSF of the bound neutral molecule, and the bound CSF of the cation are numbered from NQ 0 + 1 to NQ 0 + Nb+, Nb+ being the number of bound CSF of the cation formed from the nb bound orbitals, with spin-coupling to the ionized electron so as to give a total spin S, (NQ 1 = Nb+ × νc ). The indices I, J + used above in the alternative notation for the CSF are thus defined by I = m ≤ NQ 0, J + = m − NQ 0, m ≥ NQ 0 + 1. Thus, in the specific case of the above model for BeH2,


244116-5

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

the bound CSF of the neutral molecule are numbered from m = 1 to 50 (thus I = 1 to NQ 0 = 50), and those of the BeH2+ cation from 51 to 90 (i.e., J + = 1–40, and Nb+ = 40). The electronic configurations of those singlet CSF that are most populated during the action of the field are listed in detail in Table I.

4. Dynamics within each subspace: UQi Qi matrices

Using the matrices HQ 0Q 1(= H†Q 1Q 0), L±0 , calculated at the beginning of the [t n−1,t n ] time interval, and the values of cI and γJ,κ known at t n , the quantities in the curly brackets in Eqs. (11a) and (11b) are first computed, to be followed by the application of UQ 0Q 0 = e−iH UQ 1 Q 1 = e

M (t−t ) n

,

−iH M +(t−t n )

(14a) ⊗ u f (t, t n ),

(14b)

H M , H M+ designating the Hamiltonian matrices of the neutral molecule (M) and of the cation (M +) expressed in the basis of the bound CSF of these species. The second line of the above is the implementation of the factorization in Eq. (B11) of Appendix B, the first factor causing a transformation of the cation bound CSF’s indices J +, while the freed electron propagator u f (t, t n ) transforms the continuum orbitals (indices κ) among themselves. In the calculations reported below, it is expressed in the well-known SFA, which considers the ionized electron to be driven by the action of the external strong field only, neglecting the Coulomb forces. In the plane wave representation, u f (t, t n ) is thus represented by the Volkov propagator [u f (t, t n )]k⃗′, k⃗ = e−iΦ(t, t n ) exp[−i⃗α (t, t n ).⃗k ′]   (k ′)2 ⃗ t n )]), (15a) × exp −i δt δ(⃗k ′ − [⃗k + A(t, 2  t  t ′⃗ ′ ⃗ ⃗ ′,t n ), A(t, t n ) = dt E(t ), ⃗α (t, t n ) = dt ′ A(t tn

Φ(t, t n ) =

1 2



tn t

dt ′ A2(t ′, t n ).

III. COMPUTATIONAL DETAILS A. Input from quantum chemical electronic structure calculations

Field-free MOs, and one- and two-electron integrals hr s , gr su v , (and among them the dipole matrix elements, d r s , (d = x, y,z)), for BeH2 are computed in prior electronic-structure calculations with Columbus42 and its Argos integral package, at the HF-SCF level using the 3–21G basis. These integrals are then transformed into corresponding integrals in the MO basis using the LCAO coefficients of the MOs returned in the HF-SCF run. All these data are read in as input items of the general set of codes, (named MEDYS for Many-Electron DYnamicS), we are developing. Another important input item is the definition of the active space, i.e., the identification of active orbitals with their initial population, from which a DRT routine constructs the bound CSF of the neutral molecule and its cations. We have described this in Sec. II, with the Shavitt graph representing the DRT of the model used for single-ionization of BeH2 shown in Figure 1(a). The DRT routine also calculates matrix elements of the one- and two-electron generators Eˆr s , eˆr st u needed for the matrix representation of various components, (HQ i,Q j ), of the partitioned electronic Hamiltonian. B. Continuum orbital basis and bound-free transition dipole integrals

The present calculations assume that the one-electron ionization discretized continuum is represented in a plane ⃗ wave basis {⟨⃗r |⃗k⟩ ∝ ei k .⃗r } with wavevector ⃗k discretized on a three-dimensional grid. The continuum orbitals χκ in fact are plane waves orthogonalized to all bound (active) orbitals, and the index κ is in correspondence with the wavevector ⃗k, i.e., κ = κ(⃗k), and | χκ(k) ˆ ⃗k⟩. ⃗ ⟩ ↔ (1 − q)|

(16)

Using well-know analytical expressions for the Fourier transform of Cartesian Gaussian primitive functions43

(15b)

⃗

⃗

2

c) G(nlRm (⃗r ;ζ) = Nnl m x n y l z m e−ζ(⃗r − R c ) ,

tn

TABLE I. Sample CSF, and corresponding configurations, of the Be H2 molecule and singly charged cation as defined by the DRT represented by the Shavitt graph of Figure 1(a). The CSF numbering derives from their lexical ordering in the DRT: |I = m⟩, m ≤ 50 are bound CSF of the parent molecule, |J +⟩, J + = m − 50, 50 < m ≤ 90 are bound CSF of the cation. These actually come with a continuum orbital χ κ singly occupied, as indicated on the first row of the table. The CSF that are most populated during the action of the field are preceded by a (∗) sign. Q space CSF CSF

Configuration

P space CSF CSF

Configuration

Q space CSF CSF

Configuration

∗|1⟩

1σ g2 2σ g2 1σ u2

∗|1+⟩

1σ g2 2σ g2 1σ u1 [χ κ1 ]

|7⟩

1 1σ g2 2σ g2 1σ u1 1π u(y)

∗|2+⟩

1 1σ g2 2σ g2 1σ u1 1π u(x)

|2⟩

1σ g2 2σ g1 1σ u2

∗|8⟩

|3+⟩

1 1σ g2 2σ g1 1σ u2 1π u(x)

1 1σ g2 2σ g2 1σ u0 1π u(y)

∗|15⟩

2 1σ g2 2σ g2 1σ u0 1π u(x)

1 1σ g2 2σ g1 1σ u1 1π u(y)

|17⟩

2 1σ g2 2σ g0 1σ u2 1π u(x)

1 1σ g2 2σ g2 1σ u0 1π u(x)

∗|21⟩ ∗|22⟩

|3⟩

1 1σ g2 2σ g1 1σ u2 1π u(y)

∗|4⟩

1σ g2 2σ g2 1σ u0 1π 2y

|4+⟩

∗|5⟩

1σ g2 2σ g1 1σ u1 1π 2y

|9+⟩

∗|6⟩

2 1σ g2 2σ g0 1σ u2 1π u(y)

|10+⟩

1 1σ g2 2σ g1 1σ u1 1π u(x)

P space CSF CSF |11+⟩ |16+⟩ ∗|23+⟩

Configuration 1 1σ g2 2σ g0 1σ u2 1π u(y) [χ κ1 ] 1 1 1σ g2 2σ g1 1σ u0 1π u(y) 1π u(x)

1σ g2 2σ g1 1σ u2 3σ g1

∗|31+⟩

1σ g2 2σ g1 1σ u1 3σ g1

1σ g2 2σ g2 1σ u1 3σ g1

∗|21+⟩

1σ g2 2σ g2 1σ u0 3σ g1

1σ g2 2σ g1 1σ u2 3σ g1

∗|22+⟩

1σ g2 2σ g1 1σ u1 3σ g1


244116-6

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

as they are defined and used in most quantum chemistry packages, knowing the contraction coefficients, exponents of these primitives in the constitution of the atomic orbitals (basis functions), we calculate dipole matrix elements ⃗ c) ⟨⃗k|d α |G(nlRm ⟩,

α = x, y,z.

⃗ ⟨ χκ(k) ˆ α |ϕr ⟩ ⃗ |d α |ϕr ⟩ = ⟨ k|(1 − q)d νb 

⟨ϕ s |d α |ϕr ⟩.

(18)

while the channel-resolved ionization probabilities (i.e., the populations of the cation’s states), are defined by   ΓJ + (t) = |γJ +,κ (t)|2 (19) κ

Then, with the LCAO coefficient matrix, we calculate from this the corresponding dipole transition matrix element between a bound MO ϕr and a plane wave, and that between ϕr and a continuum orbital χκ

= ⟨⃗k|d α |ϕr ⟩ −

PI (t) = |cI (t)|2,

(17)

s=1

Together with the matrix elements of the one-electron generators Eˆr s , these are used to construct the matrix elements of Hˆ Q 0Q 1 (in matrix notations, HQ 0Q 1), according to Eq. (B5) (and of Hˆ Q 1 P1 according to Eq. (B6), if the second ionization step is included). From these, we calculate L±0 ,(L±1 ), the matrix ± ± ). , (and Lˆ Q representation of Lˆ Q 1 0 C. Choice of gauge and field

The expression of the Volkov propagator given above (Eq. (15a)) implies the use of the length gauge. We indeed consider a linearly polarized CW field, represented by an electric field vector of the form ⃗ = ϵ⃗ E0 sin(ωt + δ), E(t) ϵ⃗ is an arbitrary time-independent unit vector which defines the polarization of the field. Results will be shown for a polarization parallel to the internuclear axis of the BeH2 molecule (z polarization) and one perpendicular to this axis (called x polarization). E0 is the peak amplitude of the field, and is related to field intensity (in W/cm2) by I(W/cm2) = 3.50934 × 1016 E02(a.u.), E0 is fixed to the value of 0.16 a.u., corresponding to 8.9 × 1014 W/cm2. Results will be shown for calculations using two values of the field frequency ω: ω = 0.055 a.u. corresponds to a λ = 800 nm (Ti:Sa) laser field, and ω = 0.55 a.u., the tenth harmonic of the (Ti:Sa) field. The absolute phase δ is set to zero in the present calculations. D. Observables

Two classes of observables were considered: (i) Populations. The time-dependent populations of various channels, corresponding to the various CSF of the neutral molecule, are defined simply by

=

2 d 3⃗k|γJ +,κ(k) ⃗ (t)| .

(ii) Channel resolved photoelectron spectra. We will be showing the photoelectron momentum distribution, in the (k x ,k y ) plane, for different channels, corresponding to different CSF (labeled J +) of the cation. These are defined simply by 2 f J +(⃗k) = |γJ +,κ(k) ⃗ (t)| .

(20)

E. Grid parameters

In all calculations, the time step defining the time intervals [t n−1,t n ] above is systematically taken to be δt = TL /1000, where TL = 2π/ω is the period of the field.44 This choice was made after a thorough exploration of the convergence of the calculations with respect to δt. We found that with δt = TL /1000, the norm is conserved up to numerical errors of order 10−14, and the population profiles PI (t), ΓI (t) are invariant with further refinement of the temporal grid. The initial time is taken to be t = 0, when the field is zero, and the time slices into which the propagation time is divided are defined by t n = nδt. The numerical calculations use plane waves with wavevectors ⃗k restricted to the (k x ,k z ) plane. (Recall that the molecule’s internuclear axis is the z axis.) Systematically, the momentum grid in the direction parallel to the field (be it x ∥ = −2π or z) uses δk ∥ = π/300(a.u.), and extends from k min ∥ to k ma x = +2π, comprising 1200 gridpoints. The grid in the transverse direction uses δk ⊥ = π/50(a.u.), and extends from ⊥ ⊥ = −π to k ma k min x = +π and comprises 100 gridpoints. F. Initial state and channel-ionization potentials

All the calculations start with the neutral molecule in its ground state. This is obtained by diagonalizing the matrix of the field-free Hamiltonian in the basis of the 50 bound CFS of BeH2. This state, of energy E0Q = −19.0373 a.u., is dominated by the first CSF of the Q-subspace, corresponding to the configuration |1⟩ in Table I (1σg2 2σg2 1σu2 ). Small contributions of other CSF |m⟩ of the same symmetry, (Σg or Ag in D2h ), are found for m = 4, 6, 15, 17. Similarly, diagonalising the field-free Hamiltonian matrix restricted to the 40 remaining CSF defined by the reduced DRT of Figure 1 [those are identified with the bound states of the (N −1)-electron cation, with the N t h electron in a zero-energy orbital and non-interacting with the rest], we also obtain the ground-state of the cation and its energy, E0P = −18.5937 a.u., from which we estimate the ionisation potential to reach this channel IP (channel 1) = 0.443 a.u. (vs. 0.461 a.u. by the Koopman theorem). The first excited state of the cation


244116-7

Nguyen-Dang et al.

is at E1P = −18.5618 a.u., corresponding to IP (channel 2) = 0.475 a.u. (vs. 0.496 a.u. by the Koopman theorem). IV. RESULTS AND DISCUSSIONS

Since selection rules for bound-bound transitions and for channel-resolved ionizations vary radically with the field polarization, we will present the results for each polarization in a separate part, considering together the λ = 800 nm field and its tenth harmonic (λ = 80 nm) case. A. Parallel (z) polarization 1. Dynamics in Q-subspace

Figure 2 shows the variations, during the first cycle of the field, of the populations of the major CSF of the Q subspace. Results for λ = 800 nm are shown in the left column and those for λ = 80 nm in the right column. This organization will be kept throughout the rest of the paper, i.e., in all the figures to follow. In each column of Figure 2, panels (a) and (b) show the populations of the CSF which contribute to the initial state,

J. Chem. Phys. 141, 244116 (2014)

the main one being the CSF |1⟩ d 1σg2 2σg2 1σu2 , the time evolution of which is shown in panel (a), while panel (b) shows that of states |4⟩ − |6⟩. Two states, |21⟩ and |22⟩, are found to play an important role in the dynamics with a field polarized along the z direction. The curves showing the time-evolution of their populations are shown together in panel (c) of Figure 2. a. Case of the λ = 800 nm field. In this case, the left side of Figure 2, panel (a) depicts a fast depopulation of the CSF |1⟩ during the optical cycle, with a small, but noticeable, revival at mid-cycle. As for the other CSF |m⟩ that contribute to the initial state, m = 4, 6, 15, 17, panel (b) shows that starting with an initial value that is two orders of magnitudes smaller than that of channel |6⟩ d 1σg2 2σg0 1σu2 1π 2y (blue solid line), the population of channel |4⟩ d 1σg2 2σg2 1σu0 1π 2y (red dashed line) acquires, at about mid-cycle, an important value seemingly through an (oscillatory) population exchange with channel |6⟩. We note, however, that the latter is directly fieldcoupled to channel |5⟩ d 1σg2 2σg1 1σu1 1π 2y , rather, not to |4⟩. Indeed, for z polarization, the selection rules only allow the following radiative transitions among the bound MOs of the

FIG. 2. Variations, during the first cycle of the z-polarized field, of (i) the bound CSF |1⟩ in panel (a), (ii) bound states |4⟩ (red dashed line), |5⟩ (green dotted line), |6⟩ (blue solid line) in panel (b), (iii) bound states |21⟩ (red dashed line) and |22⟩ (blue solid line) in panel (c). The results for λ = 800 nm are shown in the left column, those for λ = 80 nm in the right column.


244116-8

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

neutral molecule or its cation:

two-photon one. Yet, it is far more important than the onephoton resonant one. This inversion between states |21⟩ and |22⟩ is reminiscent of the highly non-perturbative, strongfield effect37 that gives rise to Above-Threshold Ionization and/or Dissociation (ATI and/or ATD) spectra.38,39 As found in these processes, non-linear interactions with the field might favor the two-photon pathway (here leading to state |22⟩) over the (resonant) one-photon one (leading to state |21⟩) at some field-intensity. The effect would tend to revert (to a more normal distribution) as the intensity is decreased. To check this, we have reconsidered these populations as a function of the field-intensity. It is found that, by lowering E0 down to 0.05 a.u. (corresponding to dividing the field-intensity by 9), the (cycle-averaged) branching ratio P22/P21 drops from 104 to 103 (with P22 of order 10−4), in this λ = 80 nm case, as expected. At E0 = 0.16 a.u., the two-photon (downward) transition |6⟩ → |4⟩ is favored over the one-photon |6⟩ → |5⟩ one. These two transitions are much weaker, and appears to be of equal strength at E0 = 0.05 a.u. In the λ = 800 nm case, at least 11 photons are needed for excitations to both states |21⟩, |22⟩ and/or for ionization, and the ATI/ATD-like non-linear effects (favoring the excitation of state |22⟩ over that of state |21⟩), involve now higher-order multiphoton processes. The situation should then approach a quasi-static dynamical regime in which ionization is a tunnelling process, and non-linear interactions of bound states (and orbitals), with the field would be better described by continuous (adiabatic, field-following) distortions of these states,40 distortions that are all the more important, the higher-energy the state is. Actually, what matters is the way these field-following states are used in the time-dependent dynamics, i.e., transition from one CSF to another depends on how these CSF are involved in these distorted states. Also, just as tunnel ionization occurs mainly in the vicinity of the maximum instantaneous field-intensity, it is the projections of the initial state on the distorted states defined at these moments which govern the bound-state dynamics. The fine oscillations found in Figure 2, panels (b) and (c), in the λ = 800 nm reflect this highly multiphoton dynamics approaching a quasi-static regime. More precisely, they reflect the difference in the dynamical phases of the distorted states correlating, in the zero-field limit, to states |4⟩ and |6⟩ for the case of panel (b), and to states |21⟩ and |22⟩ for the case of panel (c). No such oscillation is found in the λ = 80 nm case. These oscillations are thus mainly field-driven as also seen by the fact that they set in only after the first maximum in the field amplitude. Note that (time-parametrized) distortions of the electron distribution in these field-following states will necessarily be accompanied by fluctuations at the level of electron interaction forces, modulating the correlation energy of these states, which in turn enters the definition of their dynamical phases.

1σu ↔ 3σg ,

1σu ↔ 2σg .

(21)

Thus, the population acquired by channel |4⟩ can only be explained by an indirect transfer from the channel |6⟩ (initially the dominant one among the 3 CSF in consideration), through the allowed transitions     |6⟩ d 1σg2 2σg0 1σu2 1π 2y ↔ |5⟩ d 1σg2 2σg1 1σu1 1π 2y   ↔ |4⟩ d 1σg2 2σg2 1σu0 1π 2y . The intermediate state |5⟩ does acquire a measurable population, as shown in Figure 2, panel (b) (curve in green dotted line). The same type of analysis holds for the observed transitions (not shown here) between |15⟩, |16⟩, and |17⟩ which correspond exactly to the same type of configurations than the three just discussed, but for the replacement of the π y (1b2u in D2h ) by π x (1b3u in D2h ). Given the above selection rules for radiative dipole transitions between the bound orbitals, it would not be surprising that the field induces a strong transition from the CSF |1⟩, (1σg2 2σg2 1σu2 ), which constitutes mostly (99.9%) the initial state, to the CSF |21⟩(d 1σg2 2σg2 1σu1 3σg1 ). This transition is clearly seen in panel (c) of Figure 2 (curve in red dashed line). But what is surprising is that the population of the bound CSF |22⟩(d 1σg2 2σg1 1σu2 3σg1 ) of BeH2, which is not directly coupled to the ground CSF by a dipole allowed transition, acquires very important values, and even at an earlier time than that of the state |21⟩, as shown in the same panel (c) of Figure 2, (curve in blue solid line). This is the manifestation of the same strongly non-linear dynamics within the system of three dipole-coupled states     |1⟩ d 1σg2 2σg2 1σu2 ↔ |21⟩ d 1σg2 2σg2 1σu1 3σg1   ↔ |22⟩ d 1σg2 2σg1 1σu2 3σg1 , (22) as in the case of the population exchange between the |4⟩ − |6⟩ states just discussed, but this time with the first state, |1⟩, initially much more populated. b. Case of the λ = 80 nm field. The last effect is much more dramatic in the case of the λ = 80 nm field, represented by the right column of Figure 2. Indeed, in panel (c), the population of state |21⟩ hardly exceeds 10−6 , while that of state |22⟩ reaches values four orders of magnitude higher. In panel (b), the same type of dynamics, but involving the states |4⟩ − |6⟩ seems to occur, with the amplitude of the intermediate state |5⟩ being much weaker as compared to the 800 nm case (shown on the left side). We now attempt to understand these observations, focusing first on the behaviour of the relative populations of states |21⟩, |22⟩. The main difference between the λ = 80 nm case and that of the λ = 800 nm field is that the transition |1⟩ → |21⟩ is resonant (transition energy ∆ϵ ≃ 0.554 a.u. ≃ ω) at λ = 80 nm. In fact, given the relative proximity of the HOMO and HOMO-1, even the transition |1⟩ → |22⟩ is almost resonant [∆ϵ(|1⟩ → |22⟩) ≃ 0.57 a.u.], but this is symmetry forbidden as a one-photon transition and is strongly non-resonant as a

2. Channel-resolved ionization dynamics

Figure 3 shows the time-evolution of the populations, ΓJ +, Eq. (19), of the major ionic channels, J + = 1,2 in panel (b), J + = 21, 22, 23, and 31 in panel (c). These are also the probabilities of ionization into the J +-th ionic channels. Panel


244116-9

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

FIG. 3. Variations, during the first cycle of the z-polarized field, of (i) the total ionization probability [panel (a)], (ii) the ionization probabilities into channels |1+⟩ (red dashed line) and |2+⟩ (blue solid line) [panel (b)], (iii) the ionization probabilities into channels |21+⟩ (red solid line), |22+⟩ (blue solid line), |23+⟩ (green dashed line), and |31+⟩ (black dotted line) [panel (c)]. The results for λ = 800 nm are shown in the left column, those for λ = 80 nm in the right column.

(a) shows the total ionization probability as a function of time. a. For λ = 800 nm field. The dominant ionization channel is the J + = 2 channel (ionic CFS |2+⟩). In fact, the first two channels |1+⟩(d 1σg2 2σg2 1σu1 ), |2+⟩(d 1σg2 2σg1 1σu2 ), appear to be populated equally at the start of the field, up to almost t = TL /4 (quarter of a period), where the Γ2(t) curve shoots up, with a clear discontinuity at that time. This seems to coincide with the sharp drop in the population of the bound CSF |22⟩ d 1σg2 2σg1 1σu2 3σg1 , seen in Figure 2, left side of panel (c). This discontinuity thus indicates the onset of a second ionization pathway leading to the J + = 2 channel, specifically the ionization of the excited 3σg electron of |22⟩ bound state BeH2(1σg2 2σg1 1σu2 3σg1 ) −→ BeH2+(1σg2 2σg1 1σu2 ) + e(⃗k) (23) that adds to the direct removal of one of the two electrons occupying the 2σg orbital in the ground state configuration, CSF |1⟩. Similarly, the ionization channel J + = 1 (ionic |1+⟩ CFS) can also be reached by the ionization (removal) of the single 3σg electron of the bound CSF |21⟩ d 1σg2 2σg2 1σu1 3σg1 . However, given that the signal of this parent molecule bound state is weaker, and comes later than that of

|22⟩, its addition to the signal of direct ionization from state |1⟩ onto this ionic channel produces a rise in its population that comes later (at a t ≃ 0.4TL ⟩TL /4) and is smaller than the rise of the population of the |2+⟩ channel. Ionization of the parent molecule pre-excited to the CFSs |21⟩, |22⟩ can also reach the ionic channels J + = 21,22,23, and 31 by the removal of either a 3σg electron or a 2σg electron. Removing the single 2σg electron from the bound state |22⟩ yields the ionic state |23+⟩(d 1σg2 2σg0 1σu1 3σg1 ). Its population, Γ23(t), is shown in green dashed line in panel (c) of Figure 3. Similarly, removing the 1σu electron from the bound state |21⟩ can yield the corresponding ionic state |21+⟩(d 1σg2 2σg2 1σu0 3σg1 ). The time-dependent population of this state, Γ21(t), shown in red solid line, seems to follow that of state |23+⟩. Disregarding the fine oscillations, the two populations are comparable in amplitude. Now, given that the population of the bound |22⟩ CSF is, in Figure 2 [panel (c)], one order of magnitude larger than that of the bound state |21⟩, one expects the signal of the resulting ionic channel |23+⟩ to be stronger than that of state |21+⟩. To explain why this expectation is not met, we note that with a furtherdipole allowed transition 1σu ↔ 2σg , we can obtain state |22+⟩(d 1σg2 2σg1 1σu1 3σg1 ) [Γ22(t) is shown in blue solid line in Figure 2, panel (c)] from |23+⟩. Then, from |22+⟩ one can again get |21+⟩ by the same dipole allowed transition


244116-10

Nguyen-Dang et al.

between the bound MOs. This is to say that the three ionic states |21+⟩, |22+⟩, and |23+⟩ are coherently linked together and the time-profiles of their populations do reflect this. We saw already that a large part of the population of the bound state |22⟩ is used to produce the ion in channel J + = 2 (second pathway discussed above). Now that population in turn gives, by an allowed (radiative) transition 1σu → 3σg , an excitation of the ionic state |31+⟩(d 1σg2 2σg1 1σu1 3σg1 ). Note that this shares the same formal configuration with state |22+⟩, but differs from it by the partial spin S (N −1) of the cation [(N − 1)-electron system], (and in the way this is coupled to the spin of the freed electron). The stepwise variation of the population Γ31(t), of this state in time, shown in black dotted line in panel (c) of Figure 3, appears indeed after the sharp increase in Γ2 at t ≃ TL /4, attributed previously to the branching of the second pathway to the ionization channel J + = 2, Eq. (23). It tracks relatively well the variations of Γ2 past that point; in particular, its increase in the second half of the optical cycle corresponds to the drop of Γ2 at these times. b. For λ = 80 nm field. In this case, the thresholds for ionization into channels J + = 1 and J + = 2 lie below the photon energy ω = 0.55 a.u. [IP (channel 1) = 0.443 a.u.; IP (channel 2) = 0.475 a.u.], and the direct formation of states |1+⟩ and |2+⟩ from the bound CSF |1⟩, |2⟩, by removal of an

J. Chem. Phys. 141, 244116 (2014)

electron in the HOMO or HOMO-1, is favored. The population of the bound state |22⟩, though still much larger than that of state |21⟩, does not exceed 6%, while in the λ = 800 nm, it reaches up to 20%. Transition from |22⟩ to |2+⟩, i.e., the second pathway to produce the cation in the J + = 2 channel is thus less probable, and the time profile of Γ2(t) shown in Figure 3, panel (b), right side, does not exhibit a discontinuity as in the case of the λ = 800 nm field, and state |1+⟩ defines the dominant ionization channel under the 80 nm field. Again, the populations of three ionic states |21+⟩, |22+⟩, and |23+⟩ are coherently linked together, as seen on the right side of panel (c) of Figure 3. The population of state |31+⟩ is no longer the largest among this group of ionic state. Its relation to state |2+⟩ is still clear from the drop in Γ2(t) accompanying the increase in Γ31(t) in the last quarter of the field period (t ≥ 3TL /4). At those times, the population of the bound state |22⟩ in this λ = 80 nm case remains appreciable, as seen in Figure 2, panel (c) for this case, and the |31+⟩ can also be obtained by removal (ionization) of one of the two 1σu electrons of state |22⟩. 3. Photoelectron spectra

Figure 4 shows, in the (k x ,k z ) plane, the photoelectron spectra associated with the two main ionization channels J + = 1, panel (a), and J + = 2, panel (b), and taken at t = TL ,

FIG. 4. Photoelectron spectra, in the (k x, k z ) plane, associated with the two main ionization channels J + = 1, |1+⟩ [panel (a)] and J + = 2, |2+⟩ [panel (b)], and taken at t = T L , the end of the first period of the λ = 800 nm (left) or λ = 80 nm (right) z-polarized field.


244116-11

Nguyen-Dang et al.

the end of the first period of the λ = 800 nm (left) or λ = 80 nm (right) field. The spectra for the shorter wavelength reflect faithfully the distributions of the transition dipole moment |⟨⃗k|z|ϕi ⟩|, where |ϕi ⟩ is the HOMO 1σu , in the case of channel J + = 1, panel (a), and the HOMO-1, 2σg , in the case of channel J + = 2, panel (b). This is not surprising as, in the previous work,11 we have seen how a few-cycle pulse of wavelengths in the XUV spectral region, corresponding to the n t h harmonics of the Ti:Sa (λ = 800 nm) field is perceived by the molecule practically as an impulsive perturbation, and give photoelectron spectra characteristic of such a perturbation. From simple symmetry considerations, it can be easily shown that ⟨⃗k|z|ϕi ⟩ is symmetric with respect to inversion of both k x and k z for the HOMO, while it is symmetric along k x but antisymmetric along k z for the HOMO-1. A contour plot of |⟨⃗k|z|1σu ⟩|2 in the (k x ,k z ) plane (not shown here) actually gives a distribution with an elliptically shaped main structure peaked at the origin, delimited by a nodal surface outside of which are found, on the two sides of this central structure, two minor structures (lobes) symmetrically placed with respect to the k z = 0 axis. The photoelectron spectrum of panel (a) for the λ = 80 nm case does have this form, the

J. Chem. Phys. 141, 244116 (2014)

central structure dominating the rest which are not seen at the scale of the plot. The one for channel J + = 2, panel (b), reflects the nodal structure of ⟨⃗k|z|2σg ⟩ (with the nodal plane at k z = 0). It is interesting to note that for λ = 80 nm, the spectrum for channel J + = 31, not shown here, appears to be the sum of the two types of spectrum as shown in panels (a) and (b), i.e., it has the central structure characteristic of ionization out of the 1σu orbital, and the two side lobes characteristic of an ionization out of the 2σg orbital. This is consistent with the discussion given above of the origin of the population of the |31+⟩ ionic state in this case. At the longer wavelength, many ionization events are possible within the field cycle, leading to photoelectron signals born at different times, and carrying nontrivial relative phase relationship between them, producing the interference patterns shown on the left side of each panel of Figure 4. There is also a symmetry breaking along the k z direction, which is normal, as the molecule now sees the oscillating field reducing its symmetry from D∞h to C∞v , and the ionization process depends on the absolute phase [δ is fixed to zero here]. This becomes the CEP (carrier-envelope phase) in the case of a pulsed field, a parameter now acknowledged

FIG. 5. Variations, during the first cycle of the x-polarized field, of the populations of (i) the bound CSF |1⟩ [panel (a)], (ii) bound states |8⟩ (red dotted line) and |22⟩ (blue solid line) [panel (b)], and (iii) ionic states |1+⟩ (red dotted line) and |2+⟩ (blue solid line) [panel (c)]. The results for λ = 800 nm are shown in the left column, those for λ = 80 nm in the right column.


244116-12

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

FIG. 6. Same as for Figure 4, but for x-polarized fields.

as playing an important role in attosecond processes. The photoelectron spectrum for ionization into the J + = 1 channel under the λ = 800 nm field is the one on the right of it in panel (a), but broken into many islands by interference effects, and is polarized toward the negative k z direction. The spectrum for the J + = 2 channel, panel (b), is also strongly polarized in this direction. However, it is more difficult to correlate this spectrum with its counterpart in the short wavelength case shown on its right, partly because signals contributing to this can come from removal of either a 2σg electron or a 3σg one, and these contributions will in turn be subjected to the interference and symmetry breaking effects mentioned above. B. Perpendicular (x ) polarization 1. Bound and ionic channel populations

To highlight the main differences between the case of a perpendicular (here x-) polarization and that of the parallel polarization discussed in detail above, we gathered in Figure 5 the variations, during the first optical cycle of an x-polarized field, of the populations of the bound configuration state |1⟩ in panel (a), those of states |8⟩ (red dotted line) and |22⟩ (blue solid line) in panel (b). Panel (c) shows the time-dependent probability of ionization into ionic channels |1+⟩ (red dotted line) and |2+⟩ (blue solid line). The total ionization probability is essentially the sum of

these two curves, although ionic states |21+⟩, |22+⟩, |23+⟩, and |31+⟩ do contribute and exhibit the same coupled dynamical behaviours as found in the case of the z− polarized field, and shown in panel (c) of Figure 3. The fundamental differences between this case of a perpendicular polarization and the case of the z− polarization arises from different selection rules. A field polarized in x can only couple σ to π orbitals, with the further selection g ↔ u rule. Since the only π orbitals included in the active space, i.e., the model system, are the (non-bonding) Beryllium 2px(y) ≡ 1πu orbitals, we can have only transitions that use the following interactions: 2σg ↔ 1πu ,

3σg ↔ πu .

(24)

Starting from the bound configuration state |1⟩ (the major component of the initial state), we would thus have tran1 sitions to state |8⟩ d 1σg2 2σg1 1σu2 1πu(x) (by a 2σg → 1πu(x) transfer) and from this to state |22⟩(d 1σg2 2σg1 1σu2 3σg1 ) (by a 1πu(x) → 3σg transfer), and both channels |8⟩ and |22⟩ acquire appreciable populations during the optical cycle, in the case of the λ = 800 nm field. Like the situation for state |21⟩ in the parallel polarization case [cf. Figure 2, panel (c)], the signal of state |8⟩ is strongly quenched in the λ = 80 nm. Again, we have here a two-photon process favored over the one-photon one by non-linear interactions with the intense field. An interesting observation can be made regarding the (field-driven) oscillations found, at λ = 800 nm, in P8(t) and


244116-13

Nguyen-Dang et al.

P22(t), and those seen in the ionization probability Γ2(t) for channel |2+⟩, panel (c) of Figure 5. The obvious relation between the oscillations found in these three different observables can be understood by realizing that there are three pathways toward ionization into channel |2+⟩, namely, (i) + |1⟩ |2+⟩, [2σg ↔ χκ(k) ⃗ ], (ii) |8⟩ |2 ⟩, [1πu ↔ χ κ(k) ⃗ ], and (iii) |22⟩ |2+⟩, [3σg ↔ χκ(k) ]. We emphasize that these are ⃗ † reversible, i.e., the action of Hˆ PQ(= Hˆ QP ) can reverse that of Hˆ QP , returning populations from the ionized channel to the bound one(s), and that the ionization into channel J + = 2 can be direct (from bound state |1⟩) or indirect (from bound state |8⟩ or |22⟩). 2. Photoelectron spectra

Figure 6 shows the photoelectron spectra associated with the two ionization channels J + = 1, panel (a), and J + = 2, panel (b). As in the case of the parallel polarization, the spectra for the λ = 80 nm case reflect faithfully the distributions of the transition dipole moment ⟨⃗k|x|ϕi ⟩ : This is antisymmetric with respect to inversion of both k x and k z for ϕi = 3σu , and antisymmetric with respect to the inversion of k x , symmetric with respect to that k z , for ϕi = 2σg . These symmetry characters give the nodal structures seen in the photoelectron spectra for channel |1+⟩ and |2+⟩ shown on the right column of the figure. In the left column, the spectra found for the λ = 800 nm exhibit, along the z-direction, i.e., in the direction normal to the polarization of the field, the nodal structure of the transition dipole moments in that direction (one nodal plane, z = 0, in the case ϕi = 3σu , the initial orbital for the ionization into the channel J + = 1, no nodal plane in the case of ionization into channel J + = 2), while in the opposite direction, they exhibit interference patterns characteristic of a non-impulsive ionization dynamics, and a symmetry breaking that washes out the nodal structure of the ⟨⃗k|x|ϕi ⟩ transition dipole moments. With a perpendicular field polarization, the symmetry is reduced from D∞h to C2v and the commutativity of the conserved symmetry operations, (those or C2v ), with the timeevolution operator implies the conservation of the characters of the initial orbital with respect to these symmetry operations. For example, the 3σu orbital is of symmetry b1 in C2v , and so is the asymptotic wavefunction describing the photoelectron originating form this MO, giving the nodal structure observed in panel (a). This symmetry conservation principle is at the root of the reading and inversion of the Laser-Induced Electron Diffraction (LIED) data of Ref. 45 and is reflected in the nodal structure of the photoelectron spectra shown on the left side of Figure 6. V. CONCLUSION AND PERSPECTIVES

The results of Sec. IV illustrate the quality of the very detailed dynamical informations that we have been able to obtain on the response of this simple, four active electrons model for the BeH2 molecule in an intense laser field. These results highlight many non-linear effects of the field: The predominance of higher-order processes is one of these, as reflected by the predominance of the population of a normally

J. Chem. Phys. 141, 244116 (2014)

inaccessible bound state (|22⟩) over the one of a doorway, in principle directly accessible state (|21⟩ in the parallel field case, |8⟩ in the perpendicular field case). The contribution of the higher-order pathway, shown in of Eq. (23), toward ionization into channel |2+⟩ in the case λ = 800 nm, is another non-linear dynamical effect of the intense laser field. The simulations show important excitations of bound electrons accompanying the ionization process, excitation that cannot be described by simulations using the so-called SAE approximation,32 and which are currently of much interest.33–35 We have shown results pertaining to the electron dynamics within a single cycle of the CW field. Calculations over a ten-cycle pulse of the λ = 800 nm, made in an extension of the present work to explore gauge invariance conditions, show that much of the subcycle details discussed in Sec. IV remain at longer times. For example, the non-monotonous ionization profile shown in panel (a) of Figure 3, is found to be repeated in all its details, with amplitude scaling that follows the modulation of the electric field oscillations by the pulse envelope, in each successive cycle of the pulse. The bases used in the present model for the expansion of bound molecular orbitals (the rather small 3–21G basis set) and for that of the N-, (N − 1)- electron bound states (as defined by the restricted active space) are obviously not complete. One thus expects that quantitative details of the dynamics will change as these bases are enlarged. However, the above main conclusions, concerning the general traits of the many-electron dynamics, are expected to remain valid. A series of test calculations including the next LUMO (the 2σu LUMO + 2 orbital) in the active space, generating a much larger CSF expansion space (105 CSF in Q0 and 70 in Q1), give an ionization probability profile with the same characteristic shape as in Figure 3, in particular, with the same abrupt rise at t ≃ TL /4, for a parallel λ = 800 nm field. The same strong quenching of state |21⟩ to the benefice of state |22⟩, as shown in Figure 2(c), is found in the case of a parallel field at λ = 80 nm. To these are added new excitation pathways, corresponding to the promotion of an electron, either from the HOMO or the HOMO-1 to the new active orbital, and exhibiting the same non-linear field-induced interactions as observed previously between states |21⟩ (or |8⟩) and |22⟩, for a parallel (or perpendicular) field, respectively. The set of codes that we are developing for many electron laser driven dynamics (MEDYS) is being improved and extended in many respects to become ready for distribution, hopefully in a near future. We are working, in particular, on the inclusion of bound-free electron repulsion, using a meanfield approach, starting from the expression of this interaction such as that given in Eq. (B7b) or a similar one but involving CSF of different subspaces. This is important for the study of electron correlation effects in the non-sequential double ionization (NSDI) process,33 for example. Double ionization is not considered in the present calculations on BeH2, but is an available option, implementing the more general propagation Eqs. (B1a)-(B2b) of Appendix B. However, the restriction of the couplings between subspaces (defined by the Feshbach partitionings of Sec. II) to radiative interactions, as evoked in Appendix B, implies that only sequential double ionization can presently be calculated. Inclusion of electron repulsion in these


244116-14

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

couplings, at least partially, is needed to describe NSDI, and this is yet to be done. The use of a basis of orthogonalized plane waves to represent continuum orbitals amounts to using a grid-based, rather costly and restrictive approach for the description of the freed electrons’ dynamics. The present calculations used a two-dimensional momentum grid for this representation of the ionized electron’s motion, and yet they have reached a limit in memory usage when calculations using an active space larger than that obtained by including the LUMO+2 orbital are attempted. The two-dimensional photoelectron momentum distributions of the type shown in Figures 4 and 6 are the observables to be analyzed in LIED-based dynamical (geometrical and structural) imaging of a molecule.45 The precision of the imaging results depends on the quality of this spectrum, i.e., on that of the momentum grid, or of the continuum orbital basis employed. Other representations of these continuum orbitals are thus being explored. Two possibilities are presently considered: The Gaussian representation of one-electron scattering wavefunctions,46 the B-spline representation of atomic and diatomic continuum orbitals.4,47 The question of the choice of gauge is also an important one. Presently, the codes implement integral reading and/or construction for calculations in the length or velocity gauge. Provisions for calculation in the Kramers-Henneberger representation,48 (also referred to as the acceleration gauge) are to be worked out, as this would require computations of the nuclear attraction potential displaced by a time-dependent quantity, i.e., a different integral reading procedure. An ongoing work explores the quality of the calculations in velocity and length gauge, with respect to gauge invariance and Ehrenfest theorem obedience, both important in proper calculations of high-order harmonic spectra.41,49 The last criterion is of particular importance as calculations of the dipole acceleration, (from which the high-order harmonic spectrum derives) by direct timedifferentiation of the time-dependent induced dipole moment, or of that of the average electron momentum, would be easier (more accessible) within our approach than by calculating the time-dependent average of the one-electron force operator. Finally, implementing a time-dependent polarization is a straightforward extension. This will permit the electron dynamics under a circularly or elliptically polarized field to be studied with the level of details illustrated here.

1. Paldus tableau

ACKNOWLEDGMENTS

3. Shavitt graph

Financial support of this research by the Natural Sciences and Engineering Research Council of Canada (NSERC, Grant No. 42395) and by the Fond de Recherche Nature et Technologie du Québec (FRQNT, Grant No. 167238) are gratefully acknowledged by T.T.N.D.

A convenient graphical representation of the DRT is provided by the Shavitt graph. In this graph, the values of the pair (a,b), as appearing in the DRT, are set on the horizontal scale, while on the vertical one are placed the orbital indices, with zero corresponding to the vacuum state. A vertex at the intersection of a horizontal line drawn at the level of the ith orbital (ith level in the DRT), and a vertical line drawn at a given pair (a,b), corresponds to a row of the DRT. The lowest vertex is called the tail, the highest one the head of the graph. The line (arc) joining two vertices at two adjacent levels indicates, by the slant it makes, how electrons are added in going from the lower level to the higher one: A vertical arc

APPENDIX A: CSF AND DRT

The concept of DRT derives from that of Paldus tableau, which is a convenient way to represent an N-electron CSF constructed out of, say n orbitals, with spin coupling to give a resultant N-electron spin quantum number S.

To each orbital, numbered i in some pre-determined order, the Paldus tableau gives a set of three non-negative integers ai , bi , ci , such that ai + bi + ci = i, 2ai + bi = Ni ,

bi = 2Si ,

(A1)

Ni being the number of electrons used in filling orbitals up to the ith one, Si is the cumulative spin quantum number up to this point, i.e., it is the spin of the said Ni -electron. The Paldus tableau is then of the form i(MO number)

ai

bi

ci

n n−1 .. . 1

an an−1 .. . a1 0

bn bn−1 .. . b1 0

cn cn−1 .. . c1 0

(A2)

and is to be read upward. The first (bottom) row, with a,b,c all equal to zero, represents the vacuum state. The next indicates how the first orbital is filled: only one of the three numbers can be nonzero (and equal to one). If a1 = 1, we have a doubly occupied orbital φ1, (so that S1 = 0), while if b1 = 1 , the orbital will be occupied only by a single electron (S1 = 1/2), and c1 = 0 indicates an empty φ1. Taking into account of the filling of the previous orbitals, conveyed by the informations from rows 1 ≤ j < i, the ai ,bi ,ci values of the ith row uniquely specify how orbital φi is to be filled. 2. DRT

For a given orbital set and total electron number N and spin S, one can regroup all the Paldus tableaux into a table called a distinct row table, DRT, from which all possible Nelectron CSF with total spin S can be read. For each orbital, we now have a group of distinct rows, and a CSF is specified by a Paldus tableau that uses one specific row of each group. The rows of the ith group (level) are related to those of the (i − 1)th group by four possible values of a chaining index, corresponding to addition (to go from the lower to the higher level) of two paired electrons, of one electron with spin-raising or spin-lowering coupling, or of zero electron.


244116-15

Nguyen-Dang et al.

corresponds to adding zero electron, a slant closest to the vertical corresponds to adding one electron with spin-raising (Si = Si−1 + 1/2) coupling, a next-to-closest to the vertical slant corresponds to adding an electron with spin-lowering (Si = Si−1 − 1/2) coupling. Finally, the strongest slant from the vertical corresponds to addition of two (paired) electrons (∆S = 0). 4. CSF and CSF lexical ordering

A CSF (or Paldus tableau) is represented in the Shavitt graph by a walk from the tail to the head of the graph, following a set of arcs joining vertices at subsequent levels. The slant of these arcs dictate the type of spin coupling coefficients to be used to construct the CSF from Slater determinants. The lexical numbering of the CSF is defined by the ordering convention that tableau m ′ (CSF |m ′⟩) precedes tableau m (CSF |m⟩), i.e., m ′ < m if at the highest level i where the two walks differ, the vertex in the walk representing |m ′⟩ is on the left of the one in the walk representing |m⟩. 5. Matrix elements of one- and two-electron generators

The value of the matrix element of a one-electron generators Eˆr s , (or of a two-electron one eˆr st u ), between two CSF |m⟩, |m ′⟩, depends on the composition of the CSF in terms of Slater determinants, and is entirely determined by the shape of the loop made up by the walks representing |m⟩ and |m ′⟩. Simple rules and relations concerning the evaluation of these matrix elements can be found tabulated in detail in the literature, for example, in Ref. 30.

APPENDIX B: GENERAL WAVEPACKET PROPAGATION EQUATIONS

Restricting to two ionization steps, the following equations govern the propagation of the |ΨQ 0⟩, |ΨQ 1⟩, |ΨP1⟩ = |ΨQ 2⟩ components of the time-dependent N-electron wavefunction |Ψ(t)⟩ over a short time-slice, t ∈ [t n ,t n+1], starting from their values supposed known at the initial time t n . First    −1 − ˆ+ |ΨQ 0(t n )⟩ |ΨQ 0(t)⟩ = Uˆ Q 0Q 0(t, t n ) Lˆ Q L Q0 0  −1   + ¯ Q P |ΨP (t n )⟩ , − i Lˆ Q H (B1a) 0 0 0 0 then, with |ΨP0(t)⟩ = |ΨQ 1(t)⟩ + |ΨP1(t)⟩, one has    −1 − ˆ+  Q (t n )⟩ |ΨQ 1(t)⟩ = Uˆ Q 1Q 1(t, t n ) LQ L |Φ Q1 1 1   −1  +  P (t n )⟩ , − i Lˆ Q H¯ Q 1 P1|Φ (B1b) 1 1 |ΨP1(t)⟩ (   −1 ) 1 +  P (t n )⟩ ¯ Q P |Φ = Uˆ P1 P1(t, t n ) 1 P1 − H¯ P1Q 1 Lˆ Q H 1 1 1 1 2   −1  +  Q (t n )⟩ , −i H¯ P1Q 1 Lˆ Q |Φ (B1c) 1 1 where   −1 +  Q (t n )⟩ = |ΨQ (t n )⟩ − i H¯ Q Q Lˆ Q |Φ |ΨQ 0(t n )⟩ 1 1 1 0 0

J. Chem. Phys. 141, 244116 (2014)

  −1  1 + − H¯ Q 1Q 0 Lˆ Q H¯ Q 0Q 1|ΨQ 1(t n )⟩ 0 2  + H¯ Q 0 P1|ΨP1(t n )⟩ , (B2a)   −1 +  P (t n )⟩ = |ΨP (t n )⟩ −i H¯ P Q Lˆ Q |Φ |ΨQ 0(t n )⟩ 1 1 1 0 0    −1 1 + H¯ Q 0Q 1|ΨQ 1(t n )⟩ − H¯ P1Q 0 Lˆ Q 0 2  + H¯ Q 0 P1|ΨP1(t n )⟩ . (B2b) In the above   1 ± Lˆ Q = 1ˆ Q i ± H¯ Q i Q i+1 H¯ Q i+1Q i , i 4

(B3)

1ˆ Q i is the identity operator restricted to the Q i subspace (note that here Qˆ 2 = Pˆ1), and H¯ Q i Q j , i, j = 0,1,2, is defined by  t ¯ HQ i Q j = dt ′Hˆ Q i Q j (t). (B4) tn

Two types of terms appear in these equations: 1. Radiative coupling between subspaces

ˆ If one restricts these to the radiative interaction Vint of h(t), ˆ Eq. (2b), then HQ 0 P1 = 0, and using Eqs. (1c) and (8)-(9b)    H Q 0Q 1 = |nb ;0C ⟩⟨nb ;0C | κ

n b n′ b

×

 r

hr s Eˆr s |nb′ ;1κ ⟩⟨nb′ ;1κ |

s

    * = |nb ;0C ⟩ hr κ ⟨nb′ ;1κ |+ n b n′ , κ r b ′ × ⟨n ;0 | Eˆr,ν +1|n ;1ν +1⟩. b

C

b

b

b

(B5) The following facts were used in going from the first to the second line: (i) Eˆr s excites an electron from the (occupied) orbital φ s , (here, φ s must be identified with χκ , otherwise the result would be zero), to orbital φr , (ii) ⟨nb ;0C | Eˆrκ |nb′ ;1κ ⟩ is independent of the value of κ. Similarly, one has, for the Q1 P1 component of Hˆ 1  H Q 1 P1 = J + K ++ r

    * | J +,κ⟩hr κ ′⟨K ++,κ,κ ′|+ , κ κ′ ′ + ++ ⟨J ,κ| Eˆr,κ ′|K ,κ,κ ⟩,

(B6)

where κ, κ ′ are two distinct, fixed, arbitrary continuum orbital indices, (for example, κ = νb + 1, κ ′ = νb + 2), representing the two symbolic orbitals that are needed to accommodate the ionization of up to two electrons and that represent two ionization continua in the reduced DRT or Shavitt graph. A consequence of Eq. (B6) is that the square matrix repre± senting the operator Hˆ Q 1 P1 Hˆ P1Q 1, and hence that of Lˆ Q (see 1 + Eq. (B3)), in the basis of the CSF | J ,κ⟩ is diagonal with respect to the continuum orbital index κ. Thanks to this, the inversion


244116-16

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

dynamics within the M + + e and M ++ + 2e, respectively, a separation between the correlated motion of the N − 1 (N − 2) electrons which remain bound in the cation and that of the liberated electron(s) can be made, if the electron correlation effects involving the ionized electrons is neglected. This rests on the fact that among the one-electron generators, those with a single continuum index, Eˆr κ , r ≤ νb , κ ≥ νb + 1, do not operate within the subspace Q i , i = 1 or 2, and among the two-electron generators eˆr st u , those with an odd number of continuum orbital indices do not operate within the subspace Q i , i = 1 or 2. One can thus write

± of Lˆ Q needed in calculating double ionization amplitudes 1 seemingly a huge matrix [of dimension (NQ 1 × NQ 1)], reduces to that of the much smaller (Nb+ × Nb+) matrix.

2. Dynamics within each subspace: Representation of Uˆ Qi ,Qi , i = 0,1,2

Among the intra-block propagators Uˆ Q i,Q i appearing in Eqs. (B1a)-(B1c), Uˆ Q 0,Q 0 can be straightforwardly implemented as a (NQ 0 × NQ 0) bound-state TDCI approach. As to Uˆ Q 1,Q 1 and Uˆ P1, P1 = Uˆ Q 2,Q 2, which denote the post-ionization

Qˆ i Hˆ 1Qˆ i =



 ( ( )  ) hr s Qˆ i Eˆr sQˆ i + hκκ ′ Qˆ i Eˆκκ ′Qˆ i ,

r, s ≤ν b

κ ′,κ



 





(Qˆ i Hˆ 1Qˆ i )f

(Qˆ i Hˆ 1Qˆ i )b

Qˆ i Hˆ 2Qˆ i =



(B7a)

( ) gr st u Qˆ i eˆr st u Qˆ i

 

(B7b)

r s ≤ν b t,u ≤ν b







(Qˆ i Hˆ 2Qˆ i )b

+

 

( ) (gr sκκ ′ + gκκ ′r s + gr κ sκ ′ + gκr κ ′s + gr κκ ′s + gκr sκ ′) Qˆ i eˆr sκκ ′Qˆ i .

r s ≤ν b κκ ′







(Qˆ i Hˆ 2Qˆ i )f

6P.

Adding Eq. (B7b) to Eq. (B7a) yields ˆ Qˆ i = (Qˆ i Hˆ t ot (t)Qˆ i ) + (Qˆ i Hˆ 1Qˆ i ) Qˆ i H(t) + (Qˆ i Hˆ 2Qˆ i )f b

f

(B8)

with (Qˆ i Hˆ t ot (t)Qˆ i )b = (Qˆ i Hˆ 1(t)Qˆ i )b + (Qˆ i Hˆ 2Qˆ i )b.

(B9)

The first term of Eq. (B8) describes the fully correlated dynamics of the (N − l), (l = 1,2) bound electrons of the cation. It commutes with the second term of Eq. (B8), which denotes the motion of the ionized electron(s) in the combined effects of the Coulomb attraction of the nuclei (bare or screened) and the external time-dependent field,   (Qˆ i Hˆ t ot (t)Qˆ i )b,(Qˆ i Hˆ 1Qˆ i )f = 0. (B10) Thus, if the last term in Eq. (B8) is neglected, then Uˆ Q i,Q i (t, t n ) can be factorized into the product of a propagator for the internal dynamics of the bound electrons of the cation with one for the motion of the ionized electrons Uˆ Q i,Q i (t, t n ) ≃ Uˆ Qb i Ht o t Q i (t, t n ) ⊗ Uˆ Qf i H1Q i (t, t n ). (B11) 1A. D. Bandrauk and H. Lu, J. Phys. B: At., Mol. Opt. Phys. 38, 2529 (2005). 2C.

Lefebvre, H. Z. Lu, S. Chelkowski, and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys. 89, 023403 (2014). 3Y. V. Vanne and A. Saenz, J. Mod. Opt. 55, 2665 (2008). 4S. Barmaki, K. Guessaf, and S. Laulan, Can. J. Phys. 89, 703 (2011). 5J. M. N. Djiokap, S. X. Hu, W.-C. Jiang, L.-Y. Peng, and A. F. Starace, New J. Phys. 14, 095010 (2012).

Krause, T. Klamroth, and P. Saalfrank, J. Chem. Phys. 123, 074105 (2005). 7N. Rohringer, A. Gordon, and R. Santra, Phys. Rev. A: At., Mol., Opt. Phys. 74, 043420 (2006). 8H. B. Schlegel, S. M. Smith, and X. Li, J. Chem. Phys. 126, 244110 (2007). 9T.-T. Nguyen-Dang, M. Peters, S.-M. Wang, E. Sinelnikov, and F. Dion, J. Chem. Phys. 127, 174107 (2007). 10T. Kato and H. Kono, Chem. Phys. Lett. 392, 533 (2004). 11T.-T. Nguyen-Dang and J. Viau-Trudel, J. Chem. Phys. 139, 244102 (2013). 12M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, Phys. Rev. A: At., Mol., Opt. Phys. 49, 2117 (1994). 13P. B. Corkum, Phys. Rev. Lett. 71, 1994 (1993). 14L. R. Moore, M. A. Lysaght, L. A. A. Nikolopoulos, J. S. Parker, H. W. van der Hart, and K. T. Taylor, J. Mod. Opt. 58, 1132 (2011). 15L. Torlina, M. Ivanov, Z. B. Walters, and O. Smirnova, Phys. Rev. A: At., Mol., Opt. Phys. 86, 043409 (2012). 16O. Zatsarinny and K. Bartschat, J. Phys. B: At., Mol. Opt. Phys. 46, 112001 (2013). 17C. M. Granados-Castro and J. L. Sanz-Vicario, J. Phys. B: At., Mol. Opt. Phys. 46, 055601 (2013). 18R. E. F. Silva, P. Rivière, and F. Martín, Phys. Rev. A: At., Mol., Opt. Phys. 85, 063414 (2012). 19M. Spanner and S. Patchkovskii, Phys. Rev. A: At., Mol., Opt. Phys. 80, 063411 (2009). 20B. Mignolet, R. D. Levine, and F. Remacle, Phys. Rev. A: At., Mol., Opt. Phys. 86, 053429 (2012). 21T.-T. Nguyen-Dang, M. Peters, S.-M. Wang, and F. Dion, Chem. Phys. 366, 71 (2009). 22W. H. Adams, J. Chem. Phys. 45, 3422 (1966). 23J. Crank and P. Nicolson, Math. Proc. Cambridge Philos. Soc. 43, 50 (1947). 24H. Kono, A. Kita, Y. Ohtsuki, and Y. Fujimura, J. Comput. Phys. 130, 148 (1997). 25B. Mignolet, A. Gijsbertsen, M. J. J. Vrakking, R. D. Levine, and F. Remacle, Phys. Chem. Chem. Phys. 13, 8331 (2011). 26B. Mignolet, R. D. Levine, and F. Remacle, J. Phys. B: At., Mol. Opt. Phys. 47, 124011 (2014). 27I. Shavitt, Int. J. Quantum Chem. 12, 131 (1977).


244116-17

Nguyen-Dang et al.

J. Chem. Phys. 141, 244116 (2014)

28I.

41A. D. Bandrauk, F. Fillion-Gourdeau, and E. Lorin, J. Phys. B: At., Mol. Opt.

29J.

Phys. 46, 153001 (2013). Lischka, R. Shepard, I. Shavitt, R. M. Pitzer, M. Dallos, T. Müller, P. G. Szalay, F. B. Brown, R. Ahlrichs, H. J. Böhm, A. Chang, D. C. Comeau, R. Gdanitz, H. Dachsel, C. Ehrhardt, M. Ernzerhof, P. Höchtl, S. Irle, G. Kedziora, T. Kovar, V. Parasuk, M. J. M. Pepper, P. Scharf, H. Schiffer, M. Schindler, M. Schüler, M. Seth, E. A. Stahlberg, J.-G. Zhao, S. Yabushita, Z. Zhang, M. Barbatti, S. Matsika, M. Schuurmann, D. R. Yarkony, S. R. Brozell, E. V. Beck, and J.-P. Blaudeau, COLUMBUS, An ab initio electronic structure program, release 5.9.1, 2006. 43R. Colle, A. Fortunelli, and S. Simonucci, Il Nuovo Cimento D 9, 969 (1987). 44For λ = 800 nm, T /1000 is about a hundredth of the time (period) associated L with a typical electronic transition, for example the transition from the ground CSF |1⟩ to the excited CSF |21⟩ (∆ϵ = 0.554 a.u. ≃ 10ω L ). 45M. Peters, T.-T. Nguyen-Dang, C. Cornaggia, S. Saugout, E. Charron, A. Keller, and O. Atabek, Phys. Rev. A: At., Mol., Opt. Phys. 83, 051403 (2011). 46A. Faure, J. D. Gorfinkiel, L. A. Morgan, and J. Tennyson, Comput. Phys. Commun. 144, 224 (2002). 47C. Marante, L. Argenti, and F. Martín, Phys. Rev. A 90, 012506 (2014). 48W. C. Henneberger, Phys. Rev. Lett. 21, 838 (1968). 49C. Granados and L. Plaja, Phys. Rev. A: At., Mol., Opt. Phys. 85, 053403 (2012).

Shavitt, Int. J. Quantum Chem. 14, 5 (1978). Paldus, in The Unitary Group, Lecture Notes in Chemistry 22, edited by J. Hinze (Springer-Verlag, Heidelberg, 1981), p. 1. 30I. Shavitt, in The Unitary Group, Lecture Notes in Chemistry 22, edited by J. Hinze (Springer-Verlag, Heidelberg, 1981), p. 51. 31R. E. Stratmann and R. R. Lucchese, J. Chem. Phys. 102, 8493 (1995). 32M. Awasthi, Y. V. Vanne, A. Saenz, A. Castro, and P. Decleva, Phys. Rev. A: At., Mol., Opt. Phys. 77, 063403 (2008). 33E. P. Månsson, D. Guénot, C. L. Arnold, D. Kroon, S. Kasper, J. M. Dahlström, E. Lindroth, A. S. Kheifets, A. L’Huillier, S. L. Sorensen, and M. Gisselbrecht, Nat. Phys. 10, 207 (2014). 34S. Sukiasyan, S. Patchkovskii, O. Smirnova, T. Brabec, and M. Y. Ivanov, Phys. Rev. A: At., Mol., Opt. Phys. 82, 043414 (2010). 35E. Pisanty and M. Ivanov, Phys. Rev. A 89, 043416 (2014). 36R. Sheppard, in Ab-initio Methods in Quantum Chemistry, Adv. Chem. Phys. LXIX, edited by K. P. Lawley (Wiley, New-York, 1987), p. 63. 37A. D. Bandrauk, E. E. Aubanel, and J.-M. Gauthier, in Molecules in Laser Fields, edited by A. D. Bandrauk (M. Dekker Publisher, New-York, 1994). 38A. Giusti-Suzor, X. He, O. Atabek, and F. H. Mies, Phys. Rev. Lett. 64, 515 (1990). 39A. Zavriyev, P. H. Bucksbaum, H. G. Muller, and D. W. Schumacher, Phys. Rev. A: At., Mol., Opt. Phys. 42, 5500 (1990). 40I. Kawata, H. Kono, Y. Fujimura, and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys. 62, 031401 (2000).

42H.


Quantum chemistry structures and properties of 134 kilo molecules.

Adiabatic quantum simulation of quantum chemistry.

Covalent Chemistry beyond Molecules.

Quantum Chemistry Meets Spectroscopy for Astrochemistry: Increasing Complexity toward Prebiotic Molecules.

Quantum chemistry in drug research.

Exploring laser-driven quantum phenomena from a time-frequency analysis perspective: a comprehensive study.

Manipulating molecules with quantum light.

Magnetically driven quantum heat engine.

Plasmon-driven photoregeneration of cofactor molecules.

Planck's quantum-driven integer quantum Hall effect in chaos.

Quantum-CORMs: quantum dot sensitized CO releasing molecules.

Thermalization of field driven quantum systems.

Quantum chemistry in dataflow: Density-Fitting MP2.

Early Experiences with Computational Quantum Chemistry.

Recent advances in quantum dot surface chemistry.

Environment-governed dynamics in driven quantum systems.

Laser optical separation of chiral molecules.

Organic chemistry. Nanomole-scale high-throughput chemistry for the synthesis of complex molecules.

Quantum chemistry investigation of fluorinated polymer systems of industrial interest.

Interval prediction of molecular properties in parametrized quantum chemistry.

Laser-driven 3-kA x-ray source.

Polyatomic molecules under intense femtosecond laser irradiation.

High efficiency quantum cascade laser frequency comb.

Reversible light-driven polymerization of polyoxometalate tethered with coumarin molecules.