Fully general time-dependent multiconfiguration self-consistent-field method for the electron-nuclear Dynamics

TL;DR

This paper introduces a comprehensive time-dependent multiconfiguration self-consistent-field method capable of modeling complex electron-nuclear dynamics involving multiple particle types and interactions, enabling first-principles simulations of molecular responses to intense laser fields.

Contribution

It develops a fully general, flexible framework for simulating correlated multielectron and multinucleus quantum dynamics with arbitrary configuration spaces.

Findings

Derives equations of motion for CI coefficients and spin-orbitals.

Allows modeling of complex molecular dynamics under intense laser fields.

Provides a first-principles approach for electron-nuclear interactions.

Abstract

We present the fully general time-dependent multiconfiguration self-consistent-field method to describe the dynamics of a system consisting of arbitrary 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.

Equations

$$\Psi (t)=\sum _{\mathbf{I}}{\Phi }_{\mathbf{I}}(t)C_{\mathbf{I}}(t),$$

$${\Phi }_{\mathbf{I}}(t)={\Phi }_{{\mathbf{I}}_1}^{(1)}(t)\otimes {\Phi }_{{\mathbf{I}}_2}^{(2)}(t)\otimes \cdots \otimes {\Phi }_{{\mathbf{I}}_K}^{(K)}(t),$$

$$[{\hat{c}}_{{\mu }_{\alpha }}^{(\alpha )},{\hat{c}}_{{\nu }_{\alpha }}^{(\alpha )}]=[{\hat{c}}_{{\mu }_{\alpha }}^{(\alpha )†},{\hat{c}}_{{\nu }_{\alpha }}^{(\alpha )†}]=0,[{\hat{c}}_{{\mu }_{\alpha }}^{(\alpha )},{\hat{c}}_{{\nu }_{\alpha }}^{(\alpha )†}]={\delta }_{{\nu }_{\alpha }}^{{\mu }_{\alpha }},$$

$$\{{\hat{c}}_{{\mu }_{\alpha }}^{(\alpha )},{\hat{c}}_{{\nu }_{\alpha }}^{(\alpha )}\}=\{{\hat{c}}_{{\mu }_{\alpha }}^{(\alpha )†},{\hat{c}}_{{\nu }_{\alpha }}^{(\alpha )†}\}=0,\{{\hat{c}}_{{\mu }_{\alpha }}^{(\alpha )},{\hat{c}}_{{\nu }_{\alpha }}^{(\alpha )†}\}={\delta }_{{\nu }_{\alpha }}^{{\mu }_{\alpha }},$$

$$|\Psi (t)⟩=\sum _{\mathbf{I}}{C}_{\mathbf{I}}(t)|\mathbf{I}(t)⟩,$$

$$|\mathbf{I}(t)⟩=|{\mathbf{I}}_1(t)⟩\otimes |{\mathbf{I}}_2(t)⟩\otimes \cdots \otimes |{\mathbf{I}}_K(t)⟩,$$

$$|{\mathbf{I}}_{\alpha }⟩=\frac{1}{{\prod }_{j_{\alpha }=1}^{n_{\alpha }}{I}_{\alpha ,j_{\alpha }}!}[{\hat{c}}_1^{†}{]}^{I_{\alpha ,1}}[{\hat{c}}_2^{†}{]}^{I_{\alpha ,2}}\cdots [{\hat{c}}_{\mu }^{†}{]}^{I_{\alpha ,n_{\alpha }}}|vac⟩.$$

$$|\Psi (t)⟩$$

$$C_{{\mathbf{I}}_1{\mathbf{I}}_2\cdots {\mathbf{I}}_K}(t)|{\mathbf{I}}_1(t)⟩\otimes |{\mathbf{I}}_2(t)⟩\otimes \cdots \otimes |{\mathbf{I}}_K(t)⟩.$$

$$H=H_1+H_2+\cdots +H_M,M\le N.$$

$$H_1=\sum _{\alpha =1}^K\sum _{p_{\alpha }=1}^{N_{\alpha }}{H}_{\alpha }(x_{\alpha ,p_{\alpha }},x_{\alpha ,p_{\alpha }}^{\prime }),$$

$$H_2$$

$$+\sum _{\alpha =1}^K\sum _{\beta >\alpha }^K\sum _{p_{\alpha }=1}^{N_{\alpha }}\sum _{q_{\alpha }=1}^{N_{\beta }}{H}_{\alpha \beta }(x_{\alpha ,p_{\alpha }},x_{\beta ,q_{\beta }},x_{\alpha ,p_{\alpha }}^{\prime },x_{\beta ,q_{\beta }}^{\prime }).$$

$${\hat{H}}_m=\sum _{m_1,\cdots ,m_K}{\hat{H}}_{m_1,\cdots ,m_K}=\sum _{\mathbf{m}}{\hat{H}}_{\mathbf{m}},$$

$${\hat{H}}_{\mathbf{m}}$$

$$=\sum _{\mathbf{\mu }\mathbf{\nu }}(H_{\mathbf{m}}{)}_{\mathbf{\nu }}^{\mathbf{\mu }}{\hat{E}}_{\mathbf{\nu }}^{\mathbf{\mu }},$$

$$({\hat{E}}_{\alpha }{)}_{{\mathbf{\nu }}_{\alpha }}^{{\mathbf{\mu }}_{\alpha }}={\hat{c}}_{{\mu }_{\alpha ,1}}^{(\alpha )†}{\hat{c}}_{{\mu }_{\alpha ,2}}^{(\alpha )†}\cdots {\hat{c}}_{{\mu }_{\alpha ,m_{\alpha }}}^{(\alpha )†}{\hat{c}}_{{\nu }_{\alpha ,m_{\alpha }}}^{(\alpha )}\cdots {\hat{c}}_{{\nu }_{\alpha ,2}}^{(\alpha )}{\hat{c}}_{{\nu }_{\alpha ,1}}^{(\alpha )},$$

$$(H_{\mathbf{m}}{)}_{\mathbf{\nu }}^{\mathbf{\mu }}=\frac{1}{{\prod }_{\alpha =1}^K{m}_{\alpha }!}\sum _{\mathbf{\mu }\mathbf{\nu }}\int \mathrm{d}\mathbf{x}\mathrm{d}{\mathbf{x}}^{\prime }{\phi }_{\mathbf{\mu }}^{\ast }(\mathbf{x})H_{\mathbf{m}}(\mathbf{x},{\mathbf{x}}^{\prime }){\phi }_{\mathbf{\nu }}({\mathbf{x}}^{\prime }),$$

$${\phi }_{\mathbf{\mu }}(\mathbf{x})$$

$$=\prod _{\alpha =1}^K{\chi }_{{\mu }_{\alpha ,1}}^{(\alpha )}(x_{\alpha ,1}){\chi }_{{\mu }_{\alpha ,2}}^{(\alpha )}(x_{\alpha ,2})\cdots {\chi }_{{\mu }_{\alpha ,m_{\alpha }}}^{(\alpha )}(x_{\alpha ,m_{\alpha }}).$$

$$({\rho }_{\mathbf{m}}{)}_{\mathbf{\nu }}^{\mathbf{\mu }}=⟨\Psi |{\hat{E}}_{\mathbf{\mu }}^{\mathbf{\nu }}|\Psi ⟩.$$

$$({\rho }_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }}=⟨\Psi |({\hat{E}}_{\alpha }{)}_{{\mu }_{\alpha }}^{{\nu }_{\alpha }}|\Psi ⟩=({\rho }_{0_1\cdots 1_{\alpha }\cdots 0_K}{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }},$$

$$({\rho }_{\alpha \alpha }{)}_{{\nu }_{\alpha }{\lambda }_{\alpha }}^{{\mu }_{\alpha }{\gamma }_{\alpha }}=⟨\Psi |({\hat{E}}_{\alpha }{)}_{{\mu }_{\alpha }{\gamma }_{\alpha }}^{{\nu }_{\alpha }{\lambda }_{\alpha }}|\Psi ⟩=({\rho }_{0_1\cdots 2_{\alpha }\cdots 0_K}{)}_{{\nu }_{\alpha }{\lambda }_{\alpha }}^{{\mu }_{\alpha }{\gamma }_{\alpha }},$$

$$({\rho }_{\alpha \beta }{)}_{{\nu }_{\alpha }{\lambda }_{\beta }}^{{\mu }_{\alpha }{\gamma }_{\beta }}=⟨\Psi |({\hat{E}}_{\alpha }{)}_{{\mu }_{\alpha }}^{{\nu }_{\alpha }}({\hat{E}}_{\alpha }{)}_{{\gamma }_{\beta }}^{{\lambda }_{\beta }}|\Psi ⟩=({\rho }_{0_1\cdots 1_{\alpha }\cdot \cdot 1_{\beta }\cdots 0_K}{)}_{{\nu }_{\alpha }{\lambda }_{\beta }}^{{\mu }_{\alpha }{\gamma }_{\beta }},$$

$$S={\int }_{t_0}^{t_1}\mathrm{d}t⟨\Psi |(\hat{H}-\mathrm{i}{\partial }_t)|\Psi ⟩,$$

$$⟨{\chi }_{{\mu }_{\alpha }}^{(\alpha )}|\delta {\chi }_{{\nu }_{\alpha }}^{(\alpha )}⟩=({\Delta }_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }},⟨{\chi }_{{\mu }_{\alpha }}^{(\alpha )}|{\dot{\chi }}_{{\nu }_{\alpha }}^{(\alpha )}⟩=(X_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }}.$$

$$\hat{\Delta }=\sum _{\alpha }\sum _{{\mu }_{\alpha }{\nu }_{\alpha }}({\Delta }_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }}({\hat{E}}_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }},\hat{X}=\sum _{\alpha }\sum _{{\mu }_{\alpha }{\nu }_{\alpha }}(X_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }}({\hat{E}}_{\alpha }{)}_{{\nu }_{\alpha }}^{{\mu }_{\alpha }},$$

$$|\delta {\chi }_{{\mu }_{\alpha }}^{(\alpha )}⟩=\hat{\Delta }|{\chi }_{{\mu }_{\alpha }}^{(\alpha )}⟩,|{\dot{\chi }}_{{\mu }_{\alpha }}^{(\alpha )}⟩=\hat{X}|{\chi }_{{\mu }_{\alpha }}^{(\alpha )}⟩.$$

$$|\delta \Psi ⟩=\sum _{\mathbf{I}}\delta C_{\mathbf{I}}|\mathbf{I}⟩+\hat{\Delta }|\Psi ⟩,|\dot{\Psi }⟩=\sum _{\mathbf{I}}{\dot{C}}_{\mathbf{I}}|\mathbf{I}⟩+\hat{X}|\Psi ⟩,$$

$$⟨\delta \Psi |=\sum _{\mathbf{I}}\delta C_{\mathbf{I}}^{\ast }⟨\mathbf{I}|-⟨\Psi |\hat{\Delta },⟨\dot{\Psi }|=\sum _{\mathbf{I}}{\dot{C}}_{\mathbf{I}}^{\ast }⟨\mathbf{I}|-⟨\Psi |\hat{X}.$$

$$\delta S$$

$$={\int }_{t_0}^{t_1}\mathrm{d}t⟨\delta \Psi |\left[H|\Psi ⟩-\mathrm{i}|\dot{\Psi }⟩\right]+\left[⟨\Psi |H+\mathrm{i}⟨\dot{\Psi }\right]|\delta \Psi ⟩.$$

$$\delta S$$

$$-\sum _{\mathbf{I}}\left\{⟨\Psi |H-\mathrm{i}\hat{X}|\mathbf{I}⟩+\mathrm{i}{\dot{C}}_{\mathbf{I}}^{\ast }\right\}\delta C_{\mathbf{I}}$$

Full text

††preprint: APS/123-QED

I 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 Agostini and DiMauro (2004); Krausz and Ivanov (2009); Gallmann et al. (2013), 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. The multiconfiguration time-dependent Hartree-Fock (MCTDHF) method Kato and Kono (2004); Caillat et al. (2005) has been developed for the investigation of multielectron dynamics in strong and/or ultrashort laser fields Ishikawa and Sato (2015). In this approach, the time-dependent total electronic wave function $\Psi(t)$ is expressed as a superposition of different Slater determinants $\Phi_{\boldsymbol{I}}(t)$,

[TABLE]

where $C_{\boldsymbol{I}}(t)$ is the configuration-interaction (CI) coefficients. Both $\{C_{\boldsymbol{I}}(t)\}$ and the spin-orbitals constituting $\{\Phi_{\boldsymbol{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 Eq. (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. The representative examples include the time-dependent complete-active-space self-consistent-field Sato and Ishikawa (2013); Sato et al. (2016), the time-dependent restricted-active-space self-consistent-field Miyagi and Madsen (2014), and the time-dependent occupation-restricted multiple active-space (TD-ORMAS) Sato and Ishikawa (2015) methods. These allow a compact and computationally less demanding description of the multielectron dynamics, without sacrificing accuracy. In particular, the TD-ORMAS can treat arbitrary CI expansions of the form Eq. (1) in principle.

Among successful approaches for nuclear dynamics is the multiconfiguration time-dependent Hartree (MCTDH) method Meyer et al. (1990). Developed for systems consisting of distinguishable particles, this method expresses the time-dependent total nuclear wave function as a superposition similar to Eq. (1) but that of Hartree products. The other way around, the MCTDHF can be viewed as an extension of the MCTDH to fermions. By hybridizing the MCTDHF for electrons and the MCTDH for nuclei, one can construct a multiconfiguration electron-nuclear dynamics (MCEND) method Ulusoy and Nest (2012) to describe the non-Born-Oppenheimer coupled dynamics. Nuclei forming molecules are, however, indistinguishable particles, either fermions or bosons.

In this Paper, further stepping forward in this direction, we present a fully general TD-MCSCF method for a system comprising of arbitrary different kinds and numbers of interacting fermions and bosons. Treating all the constituent particles on an equal footing, we expand the total wave function in terms of configurations of the whole system [see Eq. (5) below], rather than considering configurations of each particle kind separately as in Ref. Alon et al. (2012). 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.

This paper is organized as follows. Section II 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 Sec. III, 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 Sec. IV. Concluding remarks are given in Sec. V.

II Definition of the problem

II.1 TD-MCSCF ansatz

We consider a quantum mechanical many-body system with $K$ kinds of fermions or bosons. The subsystem of kind $\alpha$ consists of $N_{\alpha}$ identical particles. Thus, there are $N=\sum_{\alpha=1}^{K}N_{\alpha}$ particles in whole. For notational brevity, we call such a system an $\boldsymbol{N}$-particle system, where the array of integers $\boldsymbol{N}=(N_{1}N_{2}\cdots N_{K})$ carries information of both particle kinds and number of particles in each kind.

Let us define, for each kind of particles, the complete orthonormal set of spin-orbitals $\{\chi^{(\alpha)}_{\mu_{\alpha}}(t):\mu_{\alpha}\in\Omega_{\alpha}\}$, which spans the one-particle Hilbert space $\Omega_{\alpha}$, and are time-dependent in general. Then the $\boldsymbol{N}$-particle Hilbert space is spanned by

[TABLE]

where $\Phi^{(\alpha)}_{\boldsymbol{I}_{\alpha}}(t)$ is a determinant (or parmanent) of $\alpha$-kind fermions (or bosons), consisting of $N_{\alpha}$ spin-orbitals chosen from $\{\chi^{(\alpha)}_{\mu_{\alpha}}\}$. We call $\Phi_{\boldsymbol{I}}(t)$ the $\boldsymbol{I}$’s configuration, where $\boldsymbol{I}=\boldsymbol{I}_{1}\boldsymbol{I}_{2}\cdots\boldsymbol{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 $\boldsymbol{N}$-particle system within the ansatz of total wavefunction analogous to that for electronic system, Eq. (1), but using the configurations of Eq. (2).

For rigorous and compact presentation of theory, we resort to the second quantization formulation by introducing creation and annihilation operators $\{\hat{c}^{(\alpha)\dagger}_{\mu_{\alpha}},\hat{c}^{(\alpha)}_{\mu_{\alpha}}\}$ associated to $\{\chi^{(\alpha)}_{\mu_{\alpha}}\}$. These operators obey the (anti-)commutation relations of bosons (fermions),

[TABLE]

for bosons, where $[\hat{a},\hat{b}]=\hat{a}\hat{b}-\hat{b}\hat{a}$, and

[TABLE]

for fermions, where $\{\hat{a},\hat{b}\}=\hat{a}\hat{b}+\hat{b}\hat{a}$.

Within the TD-MCSCF ansatz, the complete set of spin-orbitals $\{\chi_{\mu_{\alpha}}^{(\alpha)}(t)\}$ is split into $n_{\alpha}$ ($\geq N_{\alpha}$) occupied spin-orbitals $\{\chi_{i_{\alpha}}^{(\alpha)}(t):i_{\alpha}=1,2,\cdots,n_{\alpha}\}$ and remaining virtual spin-orbitals $\{\chi_{a_{\alpha}}^{(\alpha)}(t):a_{\alpha}=n+1,n+2\cdots\}$. We call the subspace of $\Omega_{\alpha}$ spanned by occupied spin-orbitals the occupied spin-orbital space $\Omega^{occ}_{\alpha}$, and that spanned by virtual spin-orbitals the virtual spin-orbital space $\Omega^{vir}_{\alpha}$, where $\Omega_{\alpha}=\Omega^{occ}_{\alpha}\oplus\Omega^{vir}_{\alpha}$. The total state $\Psi(t)$ is expressed as a superposition of configurations $\Phi_{\boldsymbol{I}}(t)$ of Eq. (2), but constructed from occupied spin-orbitals only. Thus we write

[TABLE]

where $C_{\boldsymbol{I}}(t)$ is the CI coefficient, and $|\boldsymbol{I}(t)\rangle$ is the occupation number representation of the configuration $\Phi_{\boldsymbol{I}}$,

[TABLE]

[TABLE]

Now $\boldsymbol{I}_{\alpha}=\boldsymbol{I}_{\alpha,1}\boldsymbol{I}_{\alpha,2}\cdots\boldsymbol{I}_{\alpha,n_{\alpha}}$ is (rigorously) reinterpretted as an integer array, satisfying $\sum_{i_{\alpha}=1}^{n_{\alpha}}I_{\alpha,i_{\alpha}}=N_{\alpha}$. Note that $I_{\alpha,i_{\alpha}}\in\{0,1\}$ for fermions. Here and in what follows, we use indices $i_{\alpha},j_{\alpha},k_{\alpha},...$ for occupied ($\Omega^{occ}_{\alpha}$), $a_{\alpha},b_{\alpha},c_{\alpha},...$ for virtual ($\Omega^{vir}_{\alpha}$), and $\mu_{\alpha},\nu_{\alpha},\kappa_{\alpha},\tau_{\alpha},...$ for general ($\Omega_{\alpha}$) spin-orbitals of kind $\alpha$. The indices $p_{\alpha},q_{\alpha}$ will be used for numbering the coordinates.

It should be noted that we do not restrict the expansion Eq. (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 separately as Alon et al. (2012),

[TABLE]

with $\boldsymbol{I}_{\alpha},\boldsymbol{I}_{\beta},\boldsymbol{I}_{\gamma},\cdots$ being the configuration of particle kind $\alpha,\beta,\gamma,\cdots$, respectively, and $C_{\boldsymbol{I}_{\alpha}\boldsymbol{I}_{\beta}\boldsymbol{I}_{\gamma}\cdots}(t)$ the CI coefficient. Our approach allows a highly flexible choice of CI space, e.g., including up to double excitation Sato and Ishikawa (2015) regardless of particle kind, thereby enabling proper account of correlation between different kinds of particles while suppressing computational cost.

II.2 Target Hamiltonian

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

[TABLE]

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 as

[TABLE]

where the non-local form allows to describe the momentum dependence of the Hamiltonian, and two-body interaction is generally given by

[TABLE]

The reasons why we here consider the (non-local) higher-than-two body terms, which will not actually be used in Sec. IV, are (1) that such form is used in 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 freedomNakai et al. (2005a, b).

The Hamiltonian is equivalently expressed in the second quantization formalism as

[TABLE]

where the net $m$-body Hamiltonian is further classified into those contributions $\hat{H}_{\boldsymbol{m}}$, hereafter called $\boldsymbol{m}$-body Hamiltonian, involving $m_{\alpha}$ particles of the kind $\alpha$, ($0\leq m_{\alpha}\leq N_{\alpha}$, $\sum_{\alpha=1}^{K}m_{\alpha}=m$),

[TABLE]

where $\boldsymbol{\mu}=(\boldsymbol{\mu}_{1}\boldsymbol{\mu}_{2}\cdots\boldsymbol{\mu}_{K})$, and $\boldsymbol{\mu}_{\alpha}=(\mu_{\alpha,1}\mu_{\alpha,2}\cdots\mu_{\alpha,m_{\alpha}})$ indexes the set of spin-orbitals to represent $m_{\alpha}$ particles in the Hamiltonian. $\hat{E}^{\boldsymbol{\mu}}_{\boldsymbol{\nu}}$ is the $\boldsymbol{m}$-particle replacement operator $\hat{E}^{\boldsymbol{\mu}}_{\boldsymbol{\nu}}=(\hat{E}_{1})^{\boldsymbol{\mu}_{1}}_{\boldsymbol{\nu}_{1}}\cdots(\hat{E}_{K})^{\boldsymbol{\mu}_{K}}_{\boldsymbol{\nu}_{K}}$, with

[TABLE]

and $(H_{\boldsymbol{m}})^{\boldsymbol{\mu}}_{\boldsymbol{\nu}}$ is given by

[TABLE]

where $\boldsymbol{x}=(\boldsymbol{x}_{1}\boldsymbol{x}_{2}\cdots\boldsymbol{x}_{K})$, $\boldsymbol{x}_{\alpha}=(x_{\alpha,1}x_{\alpha,2}\cdots x_{\alpha,m_{\alpha}})$ is the set of $m_{\alpha}$ coordinates of particle $\alpha$, and

[TABLE]

For the later discussion, we define the $\boldsymbol{m}$-body reduced density matrix (RDM) as

[TABLE]

One- and two-particle RDMs are also denoted as

[TABLE]

with $\beta\neq\alpha$.

III Equations of Motion

In this section, we derive the EOMs for the CI coefficients and spin-orbitals by imposing the time-dependent variational principleDalgarno and Victor (1966); Löwdin and Mukherjee (1972); Moccia (1973) on our TD-MCSCF ansatz. We require the action integral

[TABLE]

to be stationary, $\delta S=0$, with respect to the variation of the total wavefunction $\delta\Psi$ within our TD-MCSCF ansatz Eq. (5), subject to the boundary conditions $\delta\Psi(t_{0})=\delta\Psi(t_{1})=0$. To this end, let us introduce anti-Hermitian matrices $\Delta_{\alpha}$ and $X_{\alpha}$ as,

[TABLE]

(Recall that indices $\mu_{\alpha},\nu_{\alpha}$ refer to both occupied and virtual spin-orbitals.) We also define,

[TABLE]

with which orthonormality-conserving spin-orbital variations and time derivatives can be written as

[TABLE]

Then, the variation and time derivative of total state are compactly given byMiranda et al. (2011); Sato and Ishikawa (2013, 2015),

[TABLE]

and their Hermitian conjugate are

[TABLE]

It follows from Eq.(19) that,

[TABLE]

Substituting Eqs.(23) and (24) into this equation, after some algebraic manipulation Miranda et al. (2011); Sato and Ishikawa (2013), we obtain,

[TABLE]

where $\Pi=\sum_{\boldsymbol{I}}|\boldsymbol{I}\rangle\langle\boldsymbol{I}|$ denotes the projector onto the CI space, i.e., the subspace of $\boldsymbol{N}$-electron Hilbert space spanned by the configurations included in Eq. (5). The action functional $S$ should be made stationary with respect to all independent variations; $\{\delta{C}_{\boldsymbol{I}},\delta{C}^{*}_{\boldsymbol{I}}\}$ for CI coefficiens and $\{(\Delta_{\alpha})^{\mu_{\alpha}}_{\nu_{\alpha}}\}$ for spin-orbitals.

First, the EOM for CI coefficients are obtained from $\delta S/\delta C^{*}_{\boldsymbol{I}}=0$,

[TABLE]

Requiring $\delta S/\delta C_{\boldsymbol{I}}=0$ derives the complex conjugate of Eq. (27). Next from $\delta S/\delta(\Delta_{\alpha})^{\mu_{\alpha}}_{\nu_{\alpha}}=0$, one obtains

[TABLE]

where $\bar{\Pi}=1-\Pi$. Equation (III) is to be solved for $(X_{\alpha})^{\mu_{\alpha}}_{\nu_{\alpha}}=\langle\chi^{(\alpha)}_{\mu_{\alpha}}|\dot{\chi}^{(\alpha)}_{\nu_{\alpha}}\rangle$, thus determines the time derivative of spin-orbitals. We now take a closer look at Eq. (III) for the following two distinct cases:

Case 1: $(\mu_{\alpha},\nu_{\alpha})=(i_{\alpha},j_{\alpha})$. In this case we focus on the components of the spin-orbital variations within the subspace spanned by the occupied spin-orbitals. Since $\bar{\Pi}(\hat{E}_{\alpha})^{i_{\alpha}}_{j_{\alpha}}|\boldsymbol{I}\rangle\neq 0$ and $\langle\boldsymbol{I}|(\hat{E}_{\alpha})^{i_{\alpha}}_{j_{\alpha}}\bar{\Pi}\neq 0$ in general, one needs to directly work with Eq. (III) within the occupied spin-orbital space

[TABLE]

In the full-CI case, where $\bar{\Pi}(\hat{E}_{\alpha})^{i_{\alpha}}_{j_{\alpha}}|\Psi\rangle=0$, $\langle\Psi|(\hat{E}_{\alpha})^{i_{\alpha}}_{j_{\alpha}}\bar{\Pi}=0$, Eq. (29) reduces to an identity $0=0$. Therefore, the corresponding $(X_{\alpha})^{i_{\alpha}}_{j_{\alpha}}$ may be arbitrary anti-Hermitian matrix elements, of which the simplest choice is $(X_{\alpha})^{i_{\alpha}}_{j_{\alpha}}=0$.

Case 2: $(\mu_{\alpha},\nu_{\alpha})=(i_{\alpha},a_{\alpha})$. In this case we deal with the components of the spin-orbital variations outside the occupied spin-orbital space. Since $\langle\Psi|(\hat{E}_{\alpha})^{i_{\alpha}}_{a_{\alpha}}\bar{\Pi}=\langle\Psi|(\hat{E}_{\alpha})^{i_{\alpha}}_{a_{\alpha}}$ and $(\hat{E}_{\alpha})^{i_{\alpha}}_{a_{\alpha}}|\Psi\rangle=0$, Eq. (III) becomes,

[TABLE]

However, the matrix element in the left-hand side of the above equation survives only when $\beta=\alpha$ and $\kappa_{\alpha}=a_{\alpha}\in\Omega^{vir}_{\alpha}$, namely $\langle\Psi|(\hat{E}_{\alpha})^{i_{\alpha}}_{a_{\alpha}}(\hat{E}_{\beta})^{\kappa_{\beta}}_{j_{\beta}}|\Psi\rangle=\delta^{\alpha}_{\beta}\delta^{a_{\alpha}}_{\kappa_{\alpha}}\langle\Psi|(\hat{E}_{\alpha})^{i_{\alpha}}_{j_{\alpha}}|\Psi\rangle=\delta^{\alpha}_{\beta}\delta^{a_{\alpha}}_{\kappa_{\alpha}}(\rho_{\alpha})_{i_{\alpha}}^{j_{\alpha}}$. Thus Eq. (30) is simplified to

[TABLE]

The $\boldsymbol{m}$-body Hamiltonian contribution to the RHS of Eq. (31) is evaluated as follows;

[TABLE]

In the second line of the above equation, we note that the matrix element survives when one and only one of the $m_{\alpha}$ creation operators in $(\hat{E}_{\alpha})^{\boldsymbol{\mu}_{\alpha}}_{\boldsymbol{\nu}_{\alpha}}$ refers to $a_{\alpha}\in\Omega^{vir}_{\alpha}$, and all the others to the occupied spin-orbitals. All such cases [$\mu_{\alpha,p}=a_{\alpha},\mu_{\alpha,q\neq p}\in\Omega^{occ}_{\alpha};1\leq p\leq m_{\alpha}$] give the same contribution since the phase $(\mp)^{p-1}$ [$+$ ($-$) sign for bosons (fermions), arising in (anti-)commuting the creation operators] is canceled by shifting the corresponding annihilation operator $\nu_{\alpha,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 $\boldsymbol{k}^{[\alpha]}=(\boldsymbol{k}_{1},\cdots\boldsymbol{k}_{\alpha},\cdots\boldsymbol{k}_{K})$ is the array of $m-1$ indices with $\boldsymbol{k}_{\alpha}=(k_{\alpha,2}\cdots,k_{\alpha,m_{\alpha}})$ and $\boldsymbol{k}_{\beta}=(k_{\beta,1}\cdots,k_{\beta,m_{\beta}})$ for $\beta\neq\alpha$, with $\boldsymbol{l}^{[\alpha]}$ defined similarly. The fourth line introduces the short-hand notation for the array of $m$ indices, $\mu_{\alpha}\boldsymbol{k}^{[\alpha]}=(\boldsymbol{k}_{1},\cdots\mu_{\alpha}\boldsymbol{k}_{\alpha},\cdots\boldsymbol{k}_{K})$ ($\mu_{\alpha}\boldsymbol{l}^{[\alpha]}$ is defined similarly), and the fifth line uses the definition of the $\boldsymbol{m}$-body RDM, Eq. (17).

Now the RHS of Eq. (31) is given by the sum over $\boldsymbol{m}$,

[TABLE]

where $(\mathbf{H}_{\alpha})_{i_{\alpha}}^{j_{\alpha}}$ is the effective one-particle operator,

[TABLE]

and $(W_{\boldsymbol{m}})^{\boldsymbol{k}^{[\alpha]}}_{\boldsymbol{l}^{[\alpha]}}$ is given in the coordinate representation as

[TABLE]

where $\boldsymbol{y}^{[\alpha]}=(\boldsymbol{y}_{1},\cdots\boldsymbol{y}_{\alpha},\cdots\boldsymbol{y}_{K})$ is the set of $m-1$ coordinates with $\boldsymbol{y}_{\alpha}=(y_{\alpha,2}\cdots,y_{\alpha,m_{\alpha}})$ and $\boldsymbol{y}_{\beta}=(y_{\beta,1}\cdots,y_{\beta,m_{\beta}})$ for $\beta\neq\alpha$, and $x_{\alpha}\boldsymbol{y}^{[\alpha]}=(\boldsymbol{y}_{1},\cdots x_{\alpha}\boldsymbol{y}_{\alpha},\cdots\boldsymbol{y}_{K})$ is the array of $m$ coordinates. $\boldsymbol{z}^{[\alpha]}$ and $x^{\prime}_{\alpha}\boldsymbol{z}^{[\alpha]}$ are defined similarly.

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

[TABLE]

where $\hat{P}_{\alpha}=\sum_{i_{\alpha}}|\chi_{i_{\alpha}}\rangle\langle\chi_{i_{\alpha}}|$ is the spin-orbital projection operator onto the occupied spin-orbital space, with which the virtual space $\sum_{a_{\alpha}}|\chi^{(\alpha)}_{a_{\alpha}}\rangle\langle\chi^{(\alpha)}_{a_{\alpha}}|=1-\hat{P}_{\alpha}$ is referenced as a whole, thus avoiding explicit use of virtual spin-orbitals. $(X_{\alpha})^{j_{\alpha}}_{i_{\alpha}}$ in the first term is to be obtained by solving Eq. (29), and, as discussed above, can be set zero in the full-CI case. Equation (27) for CI coefficients and Eq. (III) 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.

IV 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_{\rm n}$ different kinds of nuclei treated quantum mechanically (the kind does not necessarily corresponds to the nuclear species, see discussion below), and $N_{\rm cl}$ nuclei treated as a classical point charge. For clarity and notational simplicity, we assign the electrons to the first kind of particle ($\alpha=1$), and kinds $\alpha=2,3,\cdots,K$ represent quantum nuclei with $K=1+K_{\rm n}$. The numbers of identical particles are, as before, denoted by $\{N_{\alpha}\}$. Then the number of electrons is $N_{1}$, the number of quantum nuclei is $N_{\rm n}=\sum_{\alpha=2}^{K}$, and the total number of atoms is $N_{\rm atom}=N_{\rm n}+N_{\rm cl}$. We use atomic units in this section.

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

[TABLE]

where $U_{\alpha\beta}(r)=Z_{\alpha}Z_{\beta}/r$ is the Coulomb interaction with $Z_{\alpha}$ being the electric charge, and

[TABLE]

is the one-particle Hamiltonian composed of the kinetic energy [the first term with $m_{\alpha}$ being the mass (not to be confused with the number of particles)], Coulomb interaction with classical nuclei with the charges $\{Z_{A}\}$ located at $\{\boldsymbol{R}_{A}\}$ (the second term), and the time-dependent laser-particle interaction $V^{\rm ext}_{\alpha}$, given, e.g., within the dipole approximation either in the length gauge (LG) or in the velocity gauge (VG), by

[TABLE]

where $\boldsymbol{E}(t)$ is the laser electric field, and $\boldsymbol{A}(t)=-\int\boldsymbol{E}(t)\mathrm{d}t$ is the vector potential.

The general formulation of Sec. III is readily applicable to the molecular Hamiltonian of Eq. (IV). The CI EOM reads

[TABLE]

where

[TABLE]

[TABLE]

with $x_{\alpha}=(\boldsymbol{r}_{\alpha},\sigma_{\alpha})$ being the composite spatial- and spin-coordinates, and the EOM for spin-orbitals is given by

[TABLE]

where the one-body contribution to the second term of Eq. (III) is extracted to lead to $h_{\alpha}|\chi^{(\alpha)}_{i_{\alpha}}\rangle$ by noting $[(\rho_{\alpha})(\rho_{\alpha})^{-1}]^{i_{\alpha}}_{j_{\alpha}}=\delta^{i_{\alpha}}_{j_{\alpha}}$, and

[TABLE]

[TABLE]

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

[TABLE]

where $(A_{\alpha\beta})^{i_{\alpha}k_{\beta}}_{j_{\alpha}l_{\beta}}=(\bar{A}_{\alpha\beta})^{i_{\alpha}k_{\beta}}_{j_{\alpha}l_{\beta}}-(\bar{A}_{\beta\alpha})^{k_{\beta}i_{\alpha}}_{l_{\beta}j_{\alpha}}$, $(B_{\alpha})^{i_{\alpha}}_{j_{\alpha}}=(F_{\alpha})^{i_{\alpha}}_{j_{\alpha}}-(F_{\alpha})^{j_{\alpha}*}_{i_{\alpha}}$, with

[TABLE]

[TABLE]

In order for Eq. (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 Sato and Ishikawa (2015) for electrons. We shall discuss this issue in the future publication.

Equations of motions (IV) and (45), with the matrix equation (47) defines 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_{\mathrm{n}}=0,N_{\rm cl}=N_{\rm atom}$) and the full quantum molecular dynamics ($N_{\mathrm{n}}=N_{\rm atom},N_{\rm 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 interaction.

V 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 rather than each particle kind separately. 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 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.

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.

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, No. 26390076, No. 26600111, and No. 16H03881) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan and also by the Photon Frontier Network Program of MEXT. 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).