Quantum-ionic features in the absorption spectra of homonuclear diatomic molecules

We show that additional features can emerge in the linear absorption spectra of homonuclear diatomic molecules when the ions are described quantum mechanically. In particular, the widths and energies of the peaks in the optical spectra change with the initial configuration, mass, and charge of the molecule. We introduce a model that can describe these features and we provide a quantitative analysis of the resulting peak energy shifts and width broadenings as a function of the mass.


I. INTRODUCTION
Molecular spectroscopy deals with the response of a molecule interacting with an external electromagnetic field. The development of attosecond sources [1-3] allows one to probe in real time coupled electron-ion dynamics after photoionization processes. Two types of processes are seen in these experiments. The motion of the ions is associated with chemical transformations such as dissociation [4] in the femtosecond domain. The motion of the electrons is associated with electronic rearrangement processes such as charge redistribution [5,6], localization [7,8] as well as ionization processes such as tunneling [9] in the attosecond domain.
Modeling coupled electronic-ionic dynamics in photoionization processes is a formidable challenge for most systems. For this reason, previous studies have been limited to one (H + 2 ) and two (H 2 ) electron benchmark systems [8,[10][11][12][13]. A full coupled electronic-ionic 3D treatment has only been achieved for the one electron system H + 2 , where the ionic motion is confined to the direction of the laser's polarization [10,14]. For a full quantum mechanical treatment of two electron two ion systems (H 2 ), it is necessary to confine both the electronic and ionic motion to the laser's polarization direction. This is a reasonable semiclassical approximation, as the electronic and ionic motion should be predominantly along this direction [10]. Therefore, for most molecules, any quantum-ionic features are typically neglected by instead using classical approximations, e.g., the Born Oppenheimer approximation (BOA) and Ehrenfest dynamics (ED). These approaches rely on a weak coupling between the electronic and ionic wave functions. However, the validity of such approximations breaks down for light atoms, when hybridization between the electronic and ionic wave functions must be included. A quantum versus classical treatment of the ions has been previously used to investigate the localization [8], nonsequential double ionization [10] and harmonic generation [11] of H 2 , as well as the dissociation [12] and proton kinetic energies [13] for H + 2 . * alisonc1986@gmail.com † duncan.mowbray@gmail.com ‡ david.cardamone@gmail.com The aim of this paper is a comparison between a quantum mechanical (QMI) and classical (BOA/ED) treatment of the ionic motion to describe coupled electronic and ionic processes [15]. In particular, we consider three and four body systems of electrons and ions for which a fully quantum mechanical treatment of the coupled electron-ion system is feasible. This comparison with respect to the QMI solution is performed both for the static spectra and for the time dependent linear response spectra. In fact, we find significant differences between the QMI and BOA/ED spectra. These features can be quantitatively analyzed using a simple two-level twoparameter model based on the BOA electronic energy levels and the electron-ion mass ratio. The results of our work will help us to determine the domain of applicability of the simplified BOA and ED approaches to interpret coupled electron-ion experiments for more complicated systems.
The paper is organized as follows: in Sec. II, we introduce the theoretical methods and models employed to simulate the coupled electronic and ionic processes; in Sec. III, we explain the methodology to obtain both the ground state and time dependent linear response spectra, as well as the computational details of our calculations; in Sec. IV, we show our results for both, the H + 2 and H 2 molecules, which we then analyze according to the model we provide; and finally, in Sec. V we summarize the main conclusions and relevant results of our work. Atomic units a.u. ( = m e = e = a 0 = 1) are used throughout, unless stated otherwise.

A. Quantum electron-ion approach
A many-body system composed of N ions and n electrons, where both the electrons and ions are treated quantum mechanically (QMI), is described by the total electron-ion timedependent Hamiltonian H(t) =T I +T e +V II +V Ie +V ee +V ext (t), whereT I andT e are the ionic and electronic kinetic energy operators, respectively, andV II ,V Ie ,V ee , andV ext (t) are the ion-ion, ion-electron, electron-electron, and external potential energy operators, respectively. The kinetic energy operators take the formT where M α is the mass of ion α, and The interaction between the ions is given bŷ where Q α , Q β , R α and R β are the corresponding charges and positions of ion α and β. Similarly, the electron-electron repulsion isV where r i and r j are the positions of electrons i and j, while the interaction between electrons and ions iŝ Finally,V ext (t) describes the interaction of the system of electrons and ions with an external electromagnetic timedependent field, defined explicitly in Sec. III B. The time-dependent Schrödinger equation takes the form where ψ is the time-dependent electron-ion wavefunction. This depends on the positions R α and r i and on the spin coordinates S α and s i of ion α and electron i, respectively. For time-independent problems (V ext (t) = 0), the general solution of the time-dependent Schrödinger equation can be written as where ε k and ψ k are the k th eigenvalue and eigenstate of the electron-ion stationary Schrödinger equation withĤ =T I +T e +V II +V Ie +V ee .
We will next focus on the time-independent solution until introducing an external field in Sec. III B.
Solving the QMI problem is very demanding computationally for many-body systems. In fact, it quickly becomes unfeasible for systems with more than three independent variables. For this reason, we restrict consideration herein to one or two-electron diatomic molecules whose motion is confined to one direction (see Sec. II D). By applying an appropriate coordinate transformation, such systems may be modeled with only two or three independent variables (see Appendix A). In Secs. II B and II C, we introduce two of the most widely used approximations to simplify the general many-body electronion problem.

B. Born-Oppenheimer approximation
Within the Born-Oppenheimer approximation (BOA) [16], the total electronic-ionic wavefunction ψ is assumed to be separable into an ionic χ and electronic ϕ part. As the electrons move much faster than the ions, we assume that the kinetic energy of the ions does not cause the excitation of the electrons to another electronic state, i.e., an adiabatic approximation. Such an approximation is valid so long as the ratio of vibrational to electronic energies, E vib to E elec , which goes as the root of the electron-ion mass ratio, i.e., E vib /E elec ≈ √ m e /M, is small [16](see Appendix B for details). Since for a proton M p ≈ 1836m e and E vib /E elec ∼ 0.02, the BOA is expected to work quite well for our molecules. We thus may neglect T I = 0 from Eq. (10), although the electrons still feel the static field of the ions (V eI ,V II 0).
The separable BOA solution ψ of the electron-ion stationary Schrödinger equation (8) is given by [12] ψ = χ(R 1 S 1 , R 2 S 2 , ...R N S N )ϕ (R 1 ,R 2 ,...R N ) (r 1 s 1 , r 2 s 2 ...r n s n ), (11) where χ depends on the ionic coordinates only and ϕ depends on both the electronic coordinates and on the ionic coordinates which, however, only enter into the electronic wavefunctions as parameters. As shown in Ref. [12], this may be done for the Hamiltonian of Eq. 1 without loss of generality.
If we insert Eq. (11) directly into Eq. (9), we obtain the general coupled electron-ion BOA problem However, one normally separates Eq. (12) into an electronic problem only in ϕ and an ionic problem only in χ. To do so, one first solves the electronic-only BOA frozen ion Schrödinger equation, where the ionic coordinates R α only enter as fixed parameters in ϕ: whereĤ e =T e +V Ie +V ee +V II .
In this way, one may find the so-called i th potential energy surfaces E i (R 1 , R 2 , ...R N ) (PES). These are representations of the electronic energy landscape as a function of the ionic coordinates.
In the next step, the ionic BOA Schrödinger equation is solved by adding the previously neglected kinetic energy of the ions to the potential energy surfaces obtained from the frozen ion Schrödinger equation and the ionic excitations j depend on the electronic excitations, i. Comparing Eq. (12) with Eq. (16) we realize that the second and third terms of Eq. (12) are neglected in the BOA. This is because we assume that the kinetic energy of the ions is not affecting the electronic part ϕ, i.e., ∇ α ϕ ≈ 0.

C. Ehrenfest dynamics
Within the Ehrenfest dynamics (ED) scheme [17], we solve the coupled evolution of the electrons and ions. The electrons evolve quantum mechanically, whereas the ions evolve classically on a mean time-dependent PES ϕ i (t) weighted by the different BOA PES ϕ i in Eq. (13) The ions are evolved according to Newton's equation of motion which satisfies the following potential energy derivative condition where Eq. (17) and the Hellmann-Feynman theorem have been employed. The Ehrenfest electron-ion scheme consists of the time propagation of the coupled Eqs. (13) and (19).

D. Model systems: Initial configurations and Hamiltonians
We model the positively charged one electron H + 2 and neutral two electron H 2 homonuclear diatomic molecules, assuming their motion is confined to one direction. Such a model should provide a reasonable description of a molecule excited by a laser field, where the electronic and ionic motion are confined to the polarization axis of the laser field [12]. In this case the QMI problem described in Sec. II A, where both electrons and ions are treated quantum mechanically, can be solved exactly. Furthermore, by working in center of mass coordinates, the computational effort required to solve Eq. (7) is significantly reduced.
However, the singularity in the bare Coulomb interaction of Eqs. (4), (5), and (6) in 1D makes the direct numerical solution of the Schrödinger equation (7) unfeasible. Instead, one employs the so-called "soft Coulomb interaction" [12,18]. For two particles i and j with charges Q i and Q j , the soft Coulomb interaction V int has the general form where s is the separation between the two charges and ∆ is the soft Coulomb parameter [18]. Typically, ∆ = a 0 , although other values can also be used [10]. In essence, the soft Coulomb interaction amounts to a displacement of the trajectories of the two particles in an orthogonal direction. So for a hydrogen atom, a soft Coulomb interaction of is equivalent to having a bare Coulomb interaction with the electron and proton trajectories required to be parallel, with a minimum separation of a 0 . This is a quite reasonable assumption, as the most probable electron-proton separation in a hydrogen atom is the Bohr radius a 0 . Soft Coulomb parameters correspond to the separation between the 1D trajectories that each electron and ion will move along in 3D with a bare Coulomb interaction. We may directly map the 1D soft Coulomb problem to a bare Coulomb problem in 3D where the electrons and ions are separated by the soft Coulomb parameter distances shown in Fig. 1, while their motion is confined in one direction. In this way, one clearly sees that constraint rotations of the molecules are possible in 3D, while their motion is still confined to 1D. As discussed in Sec. III B, the molecules will be perturbed by a kick confined in one direction.
In Fig. 1 we show the various configurations we have employed to model an H + 2 or H 2 molecule whose electronic and ionic motion is confined to one direction. These configurations are specified by the soft Coulomb parameters between the ions ∆ II , the electrons ∆ ee and the ions and electrons ∆ Ie . One such configuration has been used previously [10] to study the dynamics of a one-dimensional H 2 model molecule in strong laser fields by means of QMI.
respectively. Here, M is the ion mass; V 1 , V 2 , X 1 and X 2 are the ionic velocities and positions for both molecules along the direction of motion along the direction of motion. The first three and four terms of Eqs. (22) and (23) are the kinetic energies of the electrons and ions in the molecules, as explained above. The remaining terms correspond to the attractive and repulsive electrostatic potential energy terms between such electrons and ions.
The spatial configuration of positively charged or neutral homogeneous diatomic molecules in 1D does not change if the particle positions are translated uniformly. This reduces our three-and four-body coordinate problems into two-and three-body ones, respectively.
We rewrite the classical energies in Eqs. (22) and (23) in terms of the center-of-mass transformation [20] (see Appendix A) to obtain the following two-body (X,ξ) and threebody (X,x,ξ) Hamiltonianŝ and for positively charged and neutral homogeneous diatomic molecules, respectively, after removing the center of mass term. Here X and x are the ionic and electronic separations and ξ is the separation between ionic and electronic centers of mass along the direction in which their motion is confined.
Although the electrons are treated quantum mechanically along their direction of motion, the confinement of their motion and position along one direction is inherently classical. For this reason, our treatment herein is essentially semiclassical: quantum mechanical along the direction of motion, and classical perpendicular to the direction of motion. This has important repercussions for the H 2 configurations shown in Fig. 1(g) and (h). For these cases, the Hamiltonian of Eq. (25) is no longer symmetric under ion or electron exchange, since ∆ Ie is either 1 3 a 0 or 2 3 a 0 . This reflects the limitations of such a semiclassical treatment. So although the Hamiltonian of Eq. (25) is still symmetric under ion or electron exchange for the H 2 configurations shown in Fig. 1(e), (f), and (i), the confinement of the electron's position perpendicular to its motion may still have an important impact for these configurations.
Our aim here is to assess the accuracy of the approximations introduced in Secs. II B and II C, in which the ions are treated classically. To accomplish this, we vary the ionic mass M in our homogeneous diatomic molecules for the many-body problem, while fixing the the ionic charge Q = e. We only consider Q = e because the repulsion between the ions of more massive homonuclear diatomic molecules with a single electron would be so large that the molecules would be unstable [19]. Furthermore, this allows us to directly compare absorption spectra between these model systems for a fixed interaction potential.

E. Symmetries of the many-body wavefunction
Since for our positively charged homogeneous diatomic molecules there are one electron and two ions, the antisymmetry of the many-body wavefunction must be enforced for the ions only as for the triplet and for the singlet. For our neutral homogeneous diatomic molecules there are two electrons and two ions. Therefore, the antisymmetry of the many-body wavefunction must be enforced both for the ions and the electrons as for the ionic triplet, (29) for the ionic singlet, for the electronic triplet, and for the electronic singlet, respectively. Therefore, due to the exchange symmetry of the many-body wavefunction, in order to have a total antisymmetric manybody wavefunction, the spatial part of the ionic and electronic wavefunction must be odd for the triplet and even for the singlet under the exchange of two identical particles. Consequently, we will only be concerned with the spatial part of the wavefunction, with the spin part already being separated off due to the exchange symmetry of the many-body wavefunction.

A. Ground state
The QMI eigenvalues are obtained by inserting Eqs. (24) and (25) into Eq. (9) for the H + 2 and H 2 molecules, respectively.
To obtain the PES within the BOA and ED we would insert Eqs. (24) and (25), neglecting the first term, into Eq. (13) for the H + 2 and H 2 molecules, respectively. For the BOA and ED ground state electron-ion level, we do not compute Eq. (15). Instead, we fit the ground state PES around its minimum energy at the inter-ionic distance X eq using a harmonic approximation E gs (X eq ) + 1 2 k 1 (X − X eq ) 2 , where k 1 = ω 2 I µ p is the harmonic constant, ω I is the harmonic oscillator vibrational frequency and µ p is the ionic reduced mass defined in Eq. (A5).
From ω I , we obtain the ground state electron-ion eigenvalue of a harmonic oscillator ε BOA/ED gs = E gs (X eq ) + 1 2 ω I in the BOA and ED PES picture.
The inversion symmetry with respect to the inter-ionic X coordinate of the potential in Eqs. (24) and (25), leads to a doubly-degenerate solution ε k for each state ψ k in Eq. (9), for sufficiently bound global ground state potentials. The inversion symmetry with respect to the inter-electronic x coordinate of the potential in Eqs. (24) and (25) is not related to the statistics of the ions, but to the symmetry of the electronic molecular orbital.

B. Time dependent linear response spectra
To obtain the linear response photoabsorption spectra we apply an initial impulsive perturbation, or "kick" [21] to the ground state wavefunctions ψ gs of our H + 2 and H 2 molecules, respectively, for the BOA and QMI approaches. K is a measure of the strength of the kick. We employ a converged kick strength of K = 0.001, for which the linear response spectra does not change if it is decreased further. Using the center of mass coordinates defined in Appendix A, the terms in Eq. (32) become where X CM 2 is the global center of mass coordinate, and ξ is the separation between the ionic and electronic centers of mass along the direction in which their motion is confined. The perturbative kick K will only induce polarization on the coordinates ξ defined for H + 2 and H 2 in Eqs. (A3) and (A8) for the BOA and QMI methods.
In linear response, we expand Eq. (33) in terms of K, neglecting higher order terms For the ED approach one should follow the same procedure starting from Eq. (32), but substituting the electronic coordinates x for − 2M+2 2M+1 ξ ≈ −ξ when M ≫ 1 for H + 2 and x 2 + x 1 for −2ξ for H 2 . During the time propagation the ions are not kicked, but evolve as parameters according to Eq. (19). In this case, the electron is kicked relative to the center of mass of the ions for the H + 2 molecule, and the two electrons are kicked relative to their distance to the ions for the H 2 molecule. However, the linear response absorption spectra does not depend on uniform translations of the ions and electrons.
The enforced time-reversal symmetry evolution operator [22] we apply to propagate our equations after this external perturbation has been applied is given by where the Hamiltonian H(t + ∆t) is calculated from and the kicked initial state we propagate is where ψ gs is the ground state eigenstate of the time independent Hamiltonian H 0 of Eqs. (24) and (25) for H + 2 and H 2 , respectively.
The time dependent Hamiltonian H(t) is then obtained by time propagation at each time step self consistently according to Eq. (36), starting from the kicked initial state given in Eq. (37). The expectation value of the dipole moment d at time t is If we assume that the Hamiltonian does not evolve in time and we insert Eq. (34) into Eq. (37), using the completeness relation k |ψ k ψ k | = 1 and Eq. (8) we get for the H + 2 molecule and for the H 2 molecule. From Eqs. (39) and (40), we see that only the ψ gs to odd k dipole moment matrix elements are non-zero by symmetry, i.e parity, sinceξ is an odd operator. Eq. (38) can be written as for the H + 2 molecules using Eq. (39) and for the H 2 molecules using Eq. (40), where d(t) depends linearly on K and ω k ≡ ε k − ε gs . However, in the ED approach the ionic coordinates are updated at each time step. This makesξ and the Hamiltonian time dependent. For this reason, the dipole moment from the ED approach does not necessarily have the form of Eqs. (41) and (42). As we will show in Sec. IV A, the time-dependent effects of zero-point motion within ED, which are incorporated into the coordinate ξ(t) = x − X CM 1 , have an important impact on the spectra. The optical photoabsorption cross section spectra σ abs is obtained by performing a discrete Fourier transform of d(t) [23]. More precisely, is a Gaussian damping applied to improve the resolution of the photoabsorption peaks, ω is the frequency of the oscillations of d(t), α is the fine structure constant, T = 1000 is the total propagation time, and ∆t = 0.01 is the time step.

C. Computational details
All numerical calculations have been performed using the real space electronic structure code Octopus [24]. We discretize the configuration space of the H + 2 and H 2 molecules, using a finite set of values (i.e. a so-called grid) for the coordinates X, x and ξ in the box intervals X ∈ [−L X , L X ], x ∈ [−L x , L x ] and ξ ∈ [−L ξ , L ξ ]. These are discretized as using N X , N x and N ξ equally spaced points, respectively. The spacing between two adjacent points in the X, x and ξ direc- Convergence is achieved when a decrease in ∆X, ∆x, ∆ξ and an increase in L X , L x , L ξ does not change the electron-ion static and time propagation linear response spectra.
Within the BOA and ED, the X coordinate does not need to be discretized quantum mechanically. It is either fixed as a parameter in BOA or it changes according to the dynamic equations in ED. As a consequence, the two and three variable bare Coulomb QMI problems confined to 1D trajectories for the H + 2 and H 2 molecules respectively, become one and two variable BOA and ED problems. These are easier to compute numerically, thus providing a more attractive alternative.

A. H + 2 and H 2 results
In Fig. 2, we show how the H + 2 and H 2 BOA ground state PES change as a function of the ionic separation X 2 + ∆ 2 II for each configuration shown in Fig. 1. The PES fitted minimum energies at X eq , E 0 (X eq ) and positions X 2 eq + ∆ 2 II are shown in Table I for the H + 2 and H 2 molecules with the configurations shown in Fig. 1.
The ground state PESs (dotted black lines in Fig. 2 and taken from Ref. 25) have been obtained by solving the stationary Schrödinger equation in 3D using basis sets within the BOA. Here, the electronic and ionic positions were allowed to vary in all spatial directions. The overall shape of these 3D PES is reproduced qualitatively by configurations (b) and (c) for H + 2 and (g) for H 2 from Fig. 1 Fig. 1.
The Hamiltonian for configuration (g) for H 2 in Fig. 1 (b) is not invariant under electron exchange. Yet, we still consider  (X eq ) and positions X 2 eq + ∆ 2 II obtained from a harmonic fit around X eq for the configurations shown in Fig. 1.  Fig. 1(e,f,i)), yield PES that differ qualitatively from the 3D PES, as shown in Fig. 2 (b). The ions sometimes undergo a strong inter-ionic repulsion for larger and small X, depending on the initial configuration. For the strongly repulsive configurations for small X, the ions are farther apart because the repulsion between the ions is stronger than the attraction between the ions and electrons. For the strongly attractive configurations for larger X, the ions are closer together because the repulsion between the ions is weaker than the attraction between the ions and electrons. For the latter configurations, more energy is required to dissociate the molecule.
For H + 2 , when ∆ Ie = a 0 , the potential becomes less repulsive for small X as ∆ II increases. However, when ∆ Ie = 0.5a 0 , the potential becomes strongly attractive for larger X.
For H 2 , the potential becomes strongly repulsive for small X for the linear configuration (g). When ∆ II = a 0 and ∆ Ie a 0 , configuration (d) in Fig. 1, the ground state PES is unbound. As the electrons are necessarily very close to each other (∆ ee = 1 3 a 0 ) when the molecule is bound, their repulsion forces the dissociation of the H 2 molecule into two isolated stabler H atoms. For this reason we will disregard this configuration from hereon.
In Fig. 3 we compare the H + 2 and H 2 optical spectra we obtain by classically fixing and letting the ions evolve according to ED in time from X eq . Essentially, including the classical movement of the ions hardly changes the spectra. However, new peaks appear before the first electronic excitation for both the H + 2 and H 2 molecule, at 1 and 12 eV, respectively. For H + 2 , the new peak corresponds to the frequency of the ionic zeropoint motion around X eq , which vanishes for large masses (M p × 10 4 ) because heavy ions hardly move around X eq . On the other hand, H 2 's higher energy peak does not vanish for large M. Due to its width, the ED and fixed ion spectra do not overlap completely. We explain the origin of this peak in Sec. IV C.
The inset of Fig. 3(a) illustrates the time-dependent effects of zero-point motion. In the BOA the ionic center of mass X CM 1 is fixed, so the electron can only oscillate about it. ED (and QMI), however, allow ionic motion, so long as the global center of mass X CM 2 is conserved. The increasing difference ∆d between the ED and BOA dipole moments thus demon-strates the ions move, e.g., ∆d(24fs) ∼ 0.7mD.
In Fig. 4, we show how a quantum mechanical treatment of the ions (QMI) affects the optical absorption spectra for H + 2 and H 2 molecules in the configurations of ∆ II , ∆ ee and ∆ Ie shown in Fig. 1. We see that new features emerge in the spectra when the ions are treated quantum mechanically instead of classically. The peaks are broadened, become asymmetric, and their amplitudes and energies change as a function of the initial configuration and charge of the molecule. In particular, comparing the BOA and QMI spectra shown in Fig. 4, we find that each peak splits into a lower and higher energy contribution. Depending on the energy shift and amplitude of each contribution, these can appear as separate peaks or shoulders in the spectra. The shoulders are giving rise to an asymmetry that can be seen for almost every peak. These quantum features are not as strong for the neutral H 2 homonuclear diatomic molecule, regardless of the initial configuration. With a classical description of the ions, we do not obtain these quantum mechanical features in the optical spectra.
Generally, we find treating the ions quantum mechanically substantially affects both the peak positions and widths in the absorption spectra for most of the configurations considered. For inter-ionic potentials which are less repulsive ( Fig. 4(b) and (e)), the line shape of the QMI peaks is narrowed, and approaches the fixed-ion at X eq limit. For potentials which are attractive for larger X (Fig. 4(d) and (f)), all the QMI peaks are blue shifted with respect to the fixed-ion at X eq spectra. In this case, the peak excitation energies are larger because more energy is required to excite these transitions.

B. Mass dependency
To provide a quantitative analysis of the differences between a classical (BOA/ED) or quantum (QMI) treatment of the ions, we will compare the total ground state energies and the peak positions and widths in the absorption spectra as we vary the ionic mass over seven orders of magnitude.
The accuracy of the static BOA and ED calculations can be understood from a perturbation theory argument in terms of the small parameter κ = (m e /M) 1 ⁄4 , defined as the ratio between the ionic and electronic displacement [29], where X = X eq + κζ. We illustrate this in detail in Appendix B.
To test the accuracy of the BOA and ED approximations, we first compare the BOA and ED ground state electron-ion eigenvalues to those obtained from QMI. We expect that the BOA and ED should be accurate around the minimum of the ground state PES, as the exact eigenvalues of the electron-ion problem can be interpreted in terms of the ionic vibrational levels for the electronic ground state PES. As discussed in Sec. III A, the ionic contribution comes from the ground state level of a quantum harmonic oscillator where the mass is included via ω I .
The ground state electron-ion eigenvalue for H + 2 and H 2 molecules whose motion is confined in one direction is given by  where ε (0) and ε (2) correspond to the electronic and ionic motion eigenvalues, and ε (1) and ε (3) are equal to zero by symmetry (see Appendix B). The first order correction to the ground state energy for the full electron-ion problem using the BOA/ED is the term of fourth order in κ.
To check the dependence of the static ground state eigenvalue accuracy of the BOA and ED approaches, on the electron-ion mass ratio, we use the following power law relation Note that most of the molecules used in this analysis are fictitious because we do not change the charge of the ions as explained in Sec. II D, except for H, D, and T, as these have a positive electric charge of e.
Fitting the ground state error from Eq. (47) to the data in Tables II and III, we obtain a power law of b ≈ 0.92(2) and b ≈ 1.05(2) for the BOA/ED approaches, as shown in Fig. 5. This means the BOA and ED energy expression gives the correct total ground state energy of the full electron-ion problem up  Figs. 1(b) and (c), and the H 2 configuration shown in Fig. 1(g).
In Fig. 6(b) we see that in the large mass limit (M ≈ 10 4 × M p ), the QMI spectra exhibits even to odd transitions which are allowed by the symmetry of the electronic wavefunctions ϕ i (ξ), shown as insets. For every allowed transition in Figs. 6 and 7, we have a red shifted and a blue shifted contribution.
The position of the first, second and fourth peaks in Fig. 6(a) and (b) (ω 1 , ω 2 , and ω 4 ) are red-shifted and the third and fifth peaks (ω 3 and ω 5 ) are blue-shifted with respect to the fixed-ion at X eq spectra. As the mass increases, all peaks tend towards the fixed-ion at X eq limit. In Fig. 7 all the peaks are red-shifted, although the second peak is also a classical peak as shown in Fig. 4(h), which disappears for smaller masses.
In Fig. 8, we show the symmetry of the occupied electronic wavefunction ϕ 0 (ξ, x) and the first ten unoccupied electronic wavefunctions ϕ i (ξ, x) for H 2 . Only transitions to unoccupied electronic wavefunctions that are even functions of x and odd functions of ξ should contribute to the absorption spectra by symmetry, i.e. ϕ i (ξ, x)|ξ|ϕ 0 (ξ, x) > 0. However, Fig. 4(h) shows there is an absorption peak in the BOA spectra for the ϕ 0 → ϕ 3 transition, despite ϕ X eq 3 (ξ, x) being an odd function of x, as shown in Fig. 8. This is because the Hamiltonian for the configuration ∆ II = 1 3 a 0 , ∆ ee = a 0 , and ∆ Ie = 1 3 ; 2 3 a 0 is not invariant under electron exchange.
From the ED spectra in Fig. 3(b), we also have an additional peak at a lower energy. When the ions are fixed, this peak is less intense than when they are allowed to evolve. In this case, the peaks' energy is given by the first excited transition (ϕ 0 → ϕ 1 ) as seen from the energy of the vertical dotted frozen ion lines in Fig. 4(f). However, the unoccupied wavefunction ϕ 1 (ξ, x) is even with respect to ξ and odd with respect to x, as shown in Fig. 8. This suggests such a transition should initially be parity forbidden.
To calculate the dipole moment for H 2 and different M, we only kick our molecules along ξ, as shown in Sec. III B.   When calculating the spectra in the BOA, the ionic coordinate is frozen at X eq and cannot evolve in time. The electronic coordinate x forms part of the integral, but can evolve in time.
When we apply a kick along ξ, the distribution of the charge in the molecule will change with time. The electrons and ions will feel the charge distribution of the other particles. Thus, the electronic coordinate x can evolve in time, although this effect is not taken into account when the dipole moment is calculated. Essentially, the (ξ, x) basis is rotated by the kick to a (ξ ′ , x ′ ) basis. As ξ is time dependent within ED, the ϕ 1 (ξ ′ , x ′ ) rotated basis has a mixture of even and odd components in both x and ξ, removing the parity constraint on the ϕ 0 → ϕ 1 transition.
As shown in Fig. 3(b), this extra parity forbidden peak is not as intense as the other peaks which are allowed by symmetry. This peak is rather weak because the electrons and ions are close to each other, but not on the same plane, for the configuration shown in Fig. 1(i). Additionally, as the ϕ 0 → ϕ 1 transition is initially forbidden by symmetry for both x and ξ, the rotated contribution with even and odd symmetry in x and ξ is small.
Overall, heavier ions have narrower peaks as we approach the classical limit. Figure 9 presents this effect in the time domain. Energy transfer between the excited electrons and the ionic system is already clearly seen after a few fs. Even in the large-mass limit (M p × 10 4 ) energy transfer is clearly evident. Oscillatory behavior, including beat frequencies, is still present 24 fs after the initial kick.
When the ions evolve quantum mechanically, the electrons can transfer part of their dipole moment to the ions. The amplitude of the dipole moment thus decreases at different rates depending on the ionic mass. This process will take longer as the mass of the ions increases and it becomes more difficult to displace the ions. For very large ion masses, the interaction with the electronic motion becomes nearly elastic. This allows the electrons to oscillate back and forth without the influence of any external ionic displacements. Since the widths of the absorption peaks are proportional to the energy transfer from the electrons to the ions, we expect the widths to scale as the electron-ion mass ratio to the one fourth, as discussed in Section II B and Appendix B.

C. Spectral lineshape
To quantify the width and energy of the peaks in the spectra, we have employed both Gaussian and Lorentzian functions. Here I i is the intensity, ω i the position, σ i the standard deviation, and Γ i the full width at half maximum of the first three peaks of the QMI spectra. From Fig. 10, in which we show the QMI spectra for H + 2 , we clearly see that the tails of the peaks of the QMI spectra are Gaussian. Moreover, the three peaks can only be fitted simultaneously with Gaussian functions, as the Lorentzian fit to the first peak decays so slowly that the second and third peaks are completely obscured. Furthermore, the ionic wave packet on the ground state PES is a solution of a harmonic eigenvalue problem and thus should have a Gaussian line shape. This means the spectral line shape arises from the shape of the PES, rather than the coupling between ionic vibrations of the molecule.
Note that the width of the fixed-ion at X eq spectra in Fig. 3 is due to the artificial damping introduced in the spectra. The electronic transitions should be delta-like functions, but are convoluted with a Gaussian function to plot the spectra (see Eqs. (43) and (44)). However, the widths in the QMI spectra are physical, and the Gaussian line shape is due to the electron-electron coupling via the ionic displacements. To ensure that the optical spectra we obtain is only affected by the external perturbation K, we have used a symmeterized initial wavefunction Thus, we always excite from a ground state which is symmetric in the ionic coordinate. In Fig. 10, we show that symmetrizing the wavefunction does not change the calculated optical absorption spectrum. This means that we already obtain a nearly symmetric ground state starting configuration from the stationary Schrödinger equation (9). However, the data shown in Fig. 6(b) has been symmetrized for every M.

D. Model
To explain why a quantum treatment of the ions has such a strong effect on the absorption spectra for H + 2 , we propose a simple two level model. Using this model, we will show how the observed QMI spectral peaks and widths can be extracted from the electronic BOA eigenenergies ε i at equilibrium X eq of the ground state through the electron-ion mass ratio m e /M.
When an external kick is applied, a charge separation is induced in the molecule which will oscillate back and forth with time. As discussed in Sec. III B, the applied kick is simply a transformation of the ground state wavefunctions, through the application of a phase factor e −iKξ , to eigenstates of the system with momentum K along the direction of motion. This leads to a transition dipole moment between an initial and a final electronic state.
The electronic wavefunctions with even indices ϕ 2i are even functions of ξ, while the electronic wavefunctions with odd indices ϕ 2i+1 are odd functions of ξ. This parity of the electronic wavefunctions means that the transition dipole moment is zero for transitions from the ground state to even unoccupied states, i.e., ϕ 2i+2 |ξ|ϕ 0 = 0. Essentially, optical transitions ϕ 0 → ϕ 2i+2 are forbidden so long as ϕ 2i+2 is an even function of ξ.
This is the case when the ions are treated classically. In fact, Figs. 4(a-d) clearly show that the peaks in the absorption spectra obtained from a classical BOA treatment are always aligned with the energies of odd-parity unoccupied electronic levels ε 2i+1 , i.e., the Franck-Condon transitions ϕ 0 → ϕ 2i+1 .
When the ions are treated quantum mechanically, every allowed transition is split into red and blue shifted contributions, with the shifts increasing as the mass decreases. The level splitting we observe in Fig. 6 is reminiscent of level hybridization.
This motivates us to employ a simple two-level model [27,28] to describe the energies and widths of the QMI peaks.
To do so, for each odd-parity unoccupied electronic level ϕ 2i+1 at ε 2i+1 , we artificially introduce a level at ε 2i+1 to which it couples.
As mentioned in Section II B and Appendix B, the ratio between the vibrational and electronic energies, E vib /E elec , scales as the square of the ratio between the ionic and electronic displacement (δ/a 0 ) 2 [29]. This means the ionic displacement scales as the electron-ion mass ratio to the one fourth δ ≈ (m e /M) 1 ⁄4 . Since the probability of coupling is directly related to the quantum ionic displacement, we expect the coupling between the energy levels to scale as δ ≈ (m e /M) 1 ⁄4 . Further, the width of the peaks in the absorption spectra should also be related to the ionic displacement. We thus assume the coupling between the energy levels is proportional to the ionic displacement δ ∼ M − 1 /4 , i.e, αM − 1 /4 where α is the constant of proportionality.
The resulting two-level Hamiltonian solution for the ε 2i+1 allowed peak's energy of This two-level model yields two peaks that are lower and higher in energy, through level repulsion. As the mass decreases, ionic displacements become larger. This leads to a greater coupling. The coupling between the energy levels will be larger when the mass decreases and the energy separation between the two coupling energy levels will also increase.
In Figs. 11 and 12 we use the two-level model to fit the calculated QMI peaks for various ionic masses M. Specifically, we employ Eqs. (51) to fit ω 1 , ω 2 , and ω 3 for H + 2 and H 2 . In general, the calculated peak positions are within 0.1 eV of the two-level model fit, which is also the expected accuracy of such calculations. In each case, we find the coupling between the transitions has a constant of proportionality of α ≈ 11 eV. Furthermore, the artificial level ε 3 ≈ 23.5 eV for both configurations of H + 2 . For H + 2 , we find the first peak depends on M − 1 ⁄2 , since α 2 /M 1 ⁄2 ≪ ε 1 − ε 1 ≈ 9.7 eV. In other words,  where ε 1 ≈ 20 eV.
For the second and third peaks, so long as ε 3 − ε 3 ≪ 2α/M 1 /4 , we may further approximate the peaks by As we see in Fig. 11, this is indeed the case for the second and third peaks in the absorption spectra, ω 2 and ω 3 , of H + 2 in the configuration ∆ Ie = a 0 and ∆ II = 1 2 a 0 , as ε 3 − ε 3 ≈ 0.4 eV. Essentially, all the peak positions in the QMI spectra are fit using only two parameters, the artificial level's energy ε 2i+1 , and the coupling constant α. However, we are not always able to decouple the two peaks' red and blue shifted contributions 2σ i ) to the first and third peaks of the absorption spectra for a positively charged homonuclear diatomic molecule with ionic mass M of µ, H, D, T, Li, Na, or K versus the fourth root of the electron-ion mass ratio (m e /M) 1 ⁄4 in the configuration ∆ Ie = a 0 and ∆ II = 1 2 a 0 ( ) or ∆ II = a 0 ( ). Black lines are linear fits to each peak for both configurations. Gray regions denote a ±0.1 eV estimated accuracy. because of their overlap due to their finite width. This is particularly true for H 2 . As a result, we have fewer data points for the H 2 peaks, as shown in Fig. 12, reducing the reliability of the fit to Eq. 51.
In Fig. 13, we also show that the width of the first and third peaks for the configurations shown in Fig. 1 (b) and (c) scale as the electron-ion mass ratio to the one fourth, i.e. FWHM ≈ (m e /M) 1 ⁄4 , as expected from our model. The FWHM for the first peak has a larger constant of proportionality than the third peak, but the widths for both configurations may be fit simultaneously. Altogether, this demonstrates the predictive power of the simple two-level model for describing the QMI spectra as a function of the ionic mass.

V. CONCLUSIONS
We have shown that additional features may appear in the linear response spectra of charged H + 2 and neutral H 2 homonuclear diatomic molecules when the ionic motion is described quantum mechanically. Such features are strongly dependent on the molecules' configuration, i.e., the shape of the PES. The widely used classical ionic motion BOA and ED approaches fail to describe such features. We also demonstrate that these features may be understood using a predictive twolevel model. These results demonstrate how for light atoms, the quantum nature of the ions may play an important role when describing absorption processes. respect to the eigenvalue ε (0) at X eq is zero. More explicitly, From Eq. (B8) we obtain the second-order correction to the energy: where the first order correction to the wavefunction is obtained from Eqs. (B7) and (B10): Here ε (0) n and |ψ (0) n are the n th electronic eigenvalue and eigenstate of the HamiltonianĤ (0) .
We now choose χ(ζ) (see Eq. (B10)) to be the lowest eigenfunction of the harmonic oscillator problem. We can then express ε (2) in the form where is the harmonic oscillator constant. Second order corrections to the energy thus correspond to the ionic vibrations.