Electronic-nuclear entanglement in Born-Oppenheimer wave functions and beyond

TL;DR

This paper investigates electronic-nuclear entanglement in molecular wave functions using various models, revealing how entanglement varies with vibrational excitation, electronic state changes, and potential energy surface features.

Contribution

It introduces a comprehensive analysis of electronic-nuclear entanglement across different molecular ansatzes, highlighting the role of vibrational states, avoided crossings, and the Born-Huang expansion.

Findings

Ground vibronic states are nearly separable and non-entangled.

Entanglement increases with vibrational excitation.

Avoided crossings significantly enhance entanglement in the adiabatic picture.

Abstract

We analyze the entanglement between electronic and nuclear motions in molecular wave functions, by using different widely used ansatzes in molecular Hamiltonian models (H$^+_2$ in 1D and the Shin-Metiu model); namely, i) Born-Oppenheimer separation in both adiabatic and ii) diabatic pictures and beyond it, iii) using a Born-Huang expansion, which involves non-adiabatic couplings in the solution. Any molecule can be considered a bipartite system in terms of electronic and nuclear halfspaces. Accordingly, the Schmidt decomposition theorem can be applied to each molecular ansatz to find a much shorter representation in terms of Schmidt basis and to estimate the entanglement content through the calculation of von Neumann entropies. Although here we justify that the ground BO vibronic state of any molecule may be regarded as an almost separable and non-entangled state, this property worsens with the vibrational excitation, which increases the entanglement monotonically. The interlace between electronic and nuclear wave functions for each vibronic state becomes even more critical for those electronic excited states whose electronic wave functions change drastically with the molecular geometry. We find that the entanglement may be quantified by the variation of the electronic wave function along the different nuclear geometries, and that the nuclear wave function indeed plays the role of a tester. In addition, the presence of avoided crossings among the potential energy curves brings about a strong enhancement of entanglement in the Born-Oppenheimer adiabatic picture, which instead is much reduced in a diabatic picture. Finally, the Born-Huang expansion uncovers some synergistic vibronic states whose entanglement content is larger than that of any Born-Oppenheimer component in the superposition.

Equations

$$\hat{H}={\hat{T}}_e+{\hat{T}}_N+\hat{U}(\mathbf{x},\mathbf{R}),$$

$${\Psi }_n^{\text{BO}}(\mathbf{x},\mathbf{R})={\varphi }_n(\mathbf{x};\mathbf{R}){\chi }_{n,m}(\mathbf{R}),$$

$$\left[{\hat{T}}_e+\hat{U}(\mathbf{x},\mathbf{R})\right]{\varphi }_n(\mathbf{x};\mathbf{R})=E_n(\mathbf{R}){\varphi }_n(\mathbf{x};\mathbf{R}),$$

$$\left[{\hat{T}}_N+E_n(\mathbf{R})\right]{\chi }_{n,m}(\mathbf{R})=W_m^n{\chi }_{n,m}(\mathbf{R}),$$

$${\Psi }^{\text{BH}}(\mathbf{x},\mathbf{R})=\sum _n{f}_n(\mathbf{R}){\varphi }_n(\mathbf{x};\mathbf{R}),$$

$$\left({\hat{T}}_N+E_n(\mathbf{R})-W_n^{\text{BH}}\right)f_n(\mathbf{R})$$

$$+\sum _{n^{\prime }}{\hat{C}}_{n,n^{\prime }}(\mathbf{R},\hat{\mathbf{P}})f_{n^{\prime }}(\mathbf{R})=0,$$

$${\hat{C}}_{n,n^{\prime }}(\mathbf{R},\hat{\mathbf{P}})=\frac{1}{2{\mu }_M}(2{\mathbf{A}}_{n,n^{\prime }}(\mathbf{R})\cdot \hat{\mathbf{P}}+B_{n,n^{\prime }}(\mathbf{R})),$$

$${\mathbf{A}}_{n,n^{\prime }}(R)=\int d\mathbf{x}{\varphi }_n(\mathbf{x};R)\frac{\partial }{\partial R}{\varphi }_{n^{\prime }}(\mathbf{x};R)$$

$$B_{n,n^{\prime }}(\mathbf{R})=\int d\mathbf{x}{\varphi }_n(\mathbf{x};R)\frac{{\partial }^2}{\partial R^2}{\varphi }_{n^{\prime }}(\mathbf{x};R).$$

$$f_n(\mathbf{R})=\sum _m{\Lambda }_{nm}{\chi }_{n,m}(\mathbf{R}),$$

$$(W_m^n-W_n^{\mathrm{BH}}){\Lambda }_{nm}+\sum _{n^{\prime }{m}^{\prime }}{C}_{n,n^{\prime },m,m^{\prime }}{\Lambda }_{n^{\prime }{m}^{\prime }}=0,$$

$${\Psi }^{\text{BH}}(\mathbf{x},\mathbf{R})=\sum _n{\varphi }_n(\mathbf{x};\mathbf{R})\sum _m{\Lambda }_{nm}{\chi }_{n,m}(\mathbf{R}).$$

$$\Psi (\mathbf{x},\mathbf{R})=\sum _{mn}{C}_{m,n}{\varphi }_n(\mathbf{x}){\chi }_m(\mathbf{R}),$$

$$|\Psi ⟩=\sum _{i=1}^{\min (M,N)}\sqrt{{\lambda }_i}|u_i⟩\otimes |v_i⟩,$$

$$S=-\text{Tr}\left[{\hat{\rho }}^x\log {\hat{\rho }}^x\right]=-\sum _i{\lambda }_i\log {\lambda }_i,$$

$$⟨{\mathbf{R}}^{\prime }\left|{\hat{\rho }}^{\mathbf{R}}\right|{\mathbf{R}}^{\prime \prime }⟩$$

$$⟨{\mathbf{x}}^{\prime }\left|{\hat{\rho }}^{\mathbf{x}}\right|{\mathbf{x}}^{\prime \prime }⟩$$

$$⟨{\mathbf{R}}^{\prime }\left|{\hat{\rho }}^{\mathbf{R}}\right|{\mathbf{R}}^{\prime \prime }⟩$$

$$⟨{\mathbf{x}}^{\prime }\left|{\hat{\rho }}^{\mathbf{x}}\right|{\mathbf{x}}^{\prime \prime }⟩$$

$$\Psi (\mathbf{x},\mathbf{R})=\sum _i^{\min (M,N)}\sqrt{{\lambda }_i}{u}_i(\mathbf{x})v_i(\mathbf{R}).$$

$$r_i=i\Delta r,$$

$$H_{ij}=\frac{2}{N}\sum _{n=1}^{(N-1)/2}\cos \left(2\pi n\frac{(i-j)}{N}\right)T_n+V(r_i){\delta }_{ij},$$

$$⟨R_i\left|{\hat{\rho }}^R\right|R_j⟩$$

$$⟨x_i\left|{\hat{\rho }}^x\right|x_j⟩$$

$$\hat{H}$$

$$-\frac{Z_{\alpha }}{\sqrt{(\hat{x}+\frac{R}{2}{)}^2+a(R)}}-\frac{Z_{\beta }}{\sqrt{(\hat{x}-\frac{R}{2}{)}^2+a(R)}},$$

$$\left[-\frac{1}{2{\mu }_N}\frac{{\partial }^2}{\partial R^2}+\frac{J(J+1)}{2{\mu }_N{R}^2}+E_n(R)-W_{J,m}^n\right]{\chi }_{n,m}^J=0,$$

$${\tilde{\chi }}_{n,m}(R_k)={\chi }_{n,m}(R_k)/\sqrt{\sum _{\ell =0}^m{\chi }_{n,m}(R_{\ell })}.$$

$${\tilde{\rho }}_{k,\ell }^{(R)}={\tilde{\chi }}_{n,m}(R_k){\tilde{\chi }}_{n,m}(R_{\ell })S_n(R_k,R_{\ell }),$$

$${\tilde{\rho }}_{k,\ell }^{(R)}={\tilde{\chi }}_{n,m}(R_k){\tilde{\chi }}_{n,m}(R_{\ell })+{ϵ}_{k,\ell }^{(n,m)},$$

$${\lambda }_{max}^{(n,m)}=1+2\sum _{k<\ell =0}^m{ϵ}_{k,\ell }^{(n,m)}{\tilde{\chi }}_{n,m}(R_k){\tilde{\chi }}_{n,m}(R_{\ell })<1.$$

$$S({\lambda }_{max}^{n,m})=-\log \left[{{\lambda }_{max}^{(n,m)}}^{{\lambda }_{max}^{(n,m)}}{\left(1-{\lambda }_{max}^{(n,m)}\right)}^{1-{\lambda }_{max}^{(n,m)}}\right].$$

Full text

††preprint: APS/123-QED

I Introduction

The Born-Oppenheimer (BO) approximation in terms of the small (large) ratio between electron and nuclear masses (velocities) has been the preamble Born1927 for the widely used adiabatic separation between electron and nuclear motion in the quantum molecular problem Born1954; Longuet-Higgins1961; Bransden2003. The corresponding BO wave function (WF) is an approximation to the exact molecular WF, namely $\Phi_{\mathrm{Mol}}({\bf x},{\bf R})\sim\phi^{\mathrm{BO}}({\bf x};{\bf R})\chi^{\mathrm{BO}}({\bf R})$, where the electronic WF $\phi^{BO}({\bf x};{\bf R})$ is obtained as an eigenstate of the electronic Hamiltonian for a chosen clamped nuclear geometry ${\bf R}$. This separation has been championed by quantum chemists for decades because it allows the concept of a nuclear motion (with a WF $\chi^{BO}({\bf R})$) driven by unambiguous potential energy surfaces provided by the electronic solution made a priori.

Other separations are possible, even an exact factorization in a similar form; $\Phi_{\mathrm{Mol}}({\bf x},{\bf R})=\phi({\bf x}|{\bf R})\chi({\bf R})$, but as a conditional probability amplitude for the electronic part, as stated by the Bayes theorem Hunter1975. This proposal has been further investigated by other authors Bishop1975; Czub1978; Gidopoulos2014; Cederbaum2013 to find that the reduction to conditional probability amplitude requires knowledge of the total WF, that the equations governing the new two factorized WFs are much more complicated to solve, and that the non-adiabatic potentials depend upon each rovibrational state. To our knowledge, the time-independent exact factorization has been applied to a limited number of soluble model systems Cederbaum2013, and the prospect of its usability in other complex molecular systems remains doubtful. Whatever the case, despite an extant exact factorization with a difficult manipulation, the molecular state $\phi({\bf x}|{\bf R})\chi({\bf R})$ still remains entangled. However, this kind of factorized solution, unfamiliar to most quantum chemists, is not within the scope of the present work.

The BO separation also provides a particular protocol for approaching the full solution of the total Hamiltonian, by proposing the BO electronic states as a complete basis (see Sutcliffe2010; Sutcliffe2012 for a historical account and a critical viewpoint on solutions based on the BO adiabatic approximation) in the so-called Born-Huang (BH) expansion Born1954. This BH expansion is not an strictly exact solution and it introduces non-adiabatic couplings between the BO adiabatic electronic states, which partially dilutes the concept of nuclei moving on a single BO adiabatic potential energy surface. The accurate computation of non-adiabatic couplings near avoided crossings or conical intersections (at geometries where the BO fails) is usually a remarkably complicated task that is avoided by the proposal of unitary transformations to a diabatic picture Smith1969. The diabatization procedure eliminates non-adiabatic couplings in favor of softer electrostatic couplings, since the diabatic WFs are no longer eigenstates of the electronic Hamiltonian.

In this work, we deal with approximate BO WFs in the adiabatic picture and their transformation into diabatic WFs, as well as beyond the BO approximation, with a BH expansion in terms of adiabatic BO WFs and finally by using a direct variational ansatz to solve the total molecular Hamiltonian. The entanglement between electron and nuclear motions has already been analyzed by several authors Bouvrie2014; McKemmish2015; Izmaylov2017; Sanz2017 and also from a time-dependent perspective Vatasescu2013, including the use of the time-dependent exact factorization ansatz Agostini2019. In our previous work Sanz2017, we computed the electro-nuclear entanglement of the simplest molecule H${}^{+}_{2}$ from first principles, without considering BO or Born-Huang ansatzes. Without bias, we diagonalize the total molecular Hamiltonian represented with a basis of tensor products of crude electron and nuclear basis functions $\Phi({\bf x},{\bf R})=\sum C_{nm}\phi_{n}({\bf x})\chi_{m}({\bf R})$. With this method of solution beyond BO, we obtained a variational solution for all vibronic states up to a given energy, devoid of any reference to BO potential energy curves. This kind of solution requires large expansions and is computationally very demanding for large molecules. However, this approach can be fairly implemented for 1D molecular models. The total variational WF is devoid of any semiclassical BO clamped nuclei approximation and much closer to the exact total WF, thus serving us as a reference.

The advantage of using molecular models of reduced dimensionality (for H${}^{+}_{2}$ and the Shin-Metiu molecular model) is that both electronic and nuclear WFs are accessible and manageable in computational terms, and can be easily visualized. Also, non-adiabatic couplings at avoided crossings can be readily computed with almost exact variational electronic WFs. Therefore, although we analyze electro-nuclear entanglement in these simple models, they include all the ingredients already present in the molecular dynamics involving excited states of complex molecules.

The correlation, coherence, and entanglement between particles within a whole molecular system are ongoing continuous objects of study from different perspectives. Correlation among particles are fully imprinted in the Hamiltonian, and its proper description depends upon the chosen ansatz for the WFs, the latter carrying the characters of coherence and entanglement. For example, the particular electron-nuclei entanglement in the photoionization of H${}_{2}^{+}$ subject to two delayed attosecond pulses has recently been theoretically and experimentally investigated Vrakking2021; vrakking2022. A desirable future goal is to understand how to take control on the content of entanglement between particles or motions within complex systems, and to take advantage of it in building efficient entanglement-based quantum-key distribution systems.

The paper is organized as follows. We introduce the theory in Section II: BO states, BH expansion, non-adiabatic couplings, the application of the Schmidt theorem, and the calculation of von Neumann entropies. Both electronic and nuclear problems are solved numerically in a 1D grid using a Fourier grid Hamiltonian method. In Section III we first analyze the molecular case of 1D H${}^{+}_{2}$ (without avoided crossings between the lowest excited states), and introduce a simplified procedure to assess the entanglement using a much reduced representation of the nuclear density matrix. Secondly, the Shin-Metiu model serves as an excellent platform to analyze entanglement with BO adiabatic vs. diabatic pictures, and with BH expansion, all compared to the entanglement content in the variational WF of the total molecular Hamiltonian. We end up with conclusions and perspectives. Atomic units are used unless otherwise stated.

II Theory

II.1 Born-Oppenheimer approximation in molecules

In molecular systems, composed of both electrons and nuclei, a natural bipartite arises between electronic and nuclear degrees of freedom. The total molecular Hamiltonian can be written as

[TABLE]

where $\hat{T}_{e}$ and $\hat{T}_{N}$ are the kinetic energy operators for electrons and nuclei, respectively, and $\hat{U}$ encompasses all Coulomb interactions: electron-electron, nuclei-nuclei, and electron-nuclei. The collective electronic coordinates are denoted by $\mathbf{x}$ and the nuclear coordinates by $\mathbf{R}$. The BO approximation is based on the observation that nuclei are much heavier and thus move more slowly than electrons. As a result, one can treat the nuclei with a fixed geometry while solving the electronic problem. This allows the total molecular WF to be approximated as a product of electronic and nuclear components, that reads

[TABLE]

where $\phi_{n}(\mathbf{x};\mathbf{R})$ is the electronic WF parameterized by the nuclear coordinates and $\chi_{n,m}(\mathbf{R})$ the nuclear WF for the $n$-th electronic state and $m$-th nuclear state. To determine the electronic WF and energy, one solves the electronic Schrödinger equation in each fixed nuclear geometry $\mathbf{R}$, i.e.

[TABLE]

where $E_{n}(\mathbf{R})$ defines the potential energy surface (PES) associated with the $n$-th electronic state. These surfaces serve as effective potentials for the nuclear motion, governed by the respective nuclear Schrödinger equation

[TABLE]

with $W^{n}_{m}$ being the total vibronic energy of the $m$-th nuclear state bound by the $n$-th PES.

The BO approximation is widely used in molecular physics and quantum chemistry because of its computational simplicity and its ability to yield accurate results when non-adiabatic effects can be safely neglected. Moreover, this approximation provides a valuable insight and serves as a foundation for more general solutions, as described in the next section.

II.2 Beyond BO: the Born-Huang ansatz and the variational ansatz.

To account for non-adiabatic effects and going beyond the standard BO approximation, one can use a more general ansatz for the total electro-nuclear WF, given by

[TABLE]

named Born-Huang (BH) expansion, where $\phi_{n}(\mathbf{x};\mathbf{R})$ is the complete set of adiabatic BO electronic eigenfunctions in Eq. (3), and $f_{n}(\mathbf{R})$ are the nuclear coefficients that carry the $\mathbf{R}$-dependence in the total state Born1954; Sutcliffe2010; Barford2013. Inserting this BH ansatz into the total molecular Hamiltonian [Eq. (1)] and projecting with $\int d\mathbf{x}\,\phi_{n^{\prime}}^{*}$, it yields the following set of coupled equations for the nuclear coefficients

[TABLE]

with $W^{\mathrm{BH}}_{n}$ being the total BH energy, and $\hat{C}_{n,n^{\prime}}(\mathbf{R},\hat{\mathbf{P}})$ corresponds to the non-adiabatic coupling operator between electronic states $n$ and $n^{\prime}$, This coupling operator is Hermitian and takes the well-known form

[TABLE]

where $\mathbf{A}_{n,n^{\prime}}$ and $B_{n,n^{\prime}}$ include the first- and second-order coupling terms, and $\mu_{M}$ is the reduced mass of the nuclei. These electronic non-adiabatic couplings are represented, for example, for diatomic molecules as

[TABLE]

[TABLE]

The vector coupling matrix elements ${\bf A}_{n,n^{\prime}}$ are usually computed pseudo-analytically or numerically within some molecular computational toolkits, but the second scalar coupling $B_{n,n^{\prime}}$ is usually ignored. This BH approach provides a convenient representation for evaluating observables, population transfer, and coherence effects in a fully coupled electro–nuclear picture. As such, it forms the foundation for many modern techniques in the treatment of photo-excitation, proton-coupled electron transfer, and other non-adiabatic collisional phenomena Nelson2020; Song2021; Bransden1992.

II.2.1 Matrix representation in a nuclear basis

To solve the coupled Eqs. (II.2) numerically, one can expand the nuclear coefficients in terms of the BO nuclear eigenfunctions $\{\chi_{n,m}(\mathbf{R})\}$,

[TABLE]

where $\Lambda_{nm}$ are the expansion coefficients. This leads to the following matrix eigenvalue problem,

[TABLE]

being $C_{n,n^{\prime},m,m^{\prime}}$ the vibronic non-adiabatic coupling matrix elements. Finally, the total electro-nuclear WF is then expanded as

[TABLE]

In case all non-adiabatic couplings are neglected, the expansion reduces to a single BO term.

II.2.2 Variational ansatz for the total wave function.

A direct variational solution of the total molecular Hamiltonian in Eq.  (1) can also be implemented by using a configuration interaction tensor product, with separate orthogonal basis functions for electrons and nuclei to expand each molecular state, i.e.,

[TABLE]

where the expansion coefficients $C_{m,n}$ are obtained computationally after diagonalization. Of course, this approach can be implemented in practice for Hamiltonians with a reduced number of degrees of freedom.

II.3 Schmidt theorem for bipartite systems

The structure of molecular WFs under the BO, the BH or the full variational ansatz naturally leads to a bipartite quantum system for the electronic and nuclear halfspaces. For this case, the singular value decomposition (SVD) provides a complete representation in terms of two subsystems, $U$ and $V$. This SVD-based representation is mathematically equivalent to the Schmidt theorem Gerry2004; Sanz2017; Blavier2022. Therefore, employing the Schmidt decomposition with all singular values provides an exact description of the quantum state, although in practice only a limited number of separable terms are often required to achieve an accurate approximation.

II.3.1 Schmidt theorem for a pure quantum state

Given a pure state $\left|\Psi\right\rangle$ of a bipartite system ($U$,$V$) with Hilbert space $\mathcal{H}=\mathcal{H}_{U}\oplus\mathcal{H}_{V}$, the Schmidt theorem confirms the existence of new orthonormal sets $\left\{u_{1},\dots,u_{M}\right\}\in\mathcal{H}_{U}$ and $\left\{v_{1},\dots,v_{N}\right\}\in\mathcal{H}_{V}$ such that

[TABLE]

with $\left|u_{i}\right>$ and $\left|v_{i}\right>$ being the reduced Schmidt bases, obtained as eigenstates of the reduced density matrices in the half-spaces, namely $\hat{\rho}^{U}=\text{Tr}_{V}\left[\hat{\rho}\right]$ and $\hat{\rho}^{V}=\text{Tr}_{U}\left[\hat{\rho}\right]$. It should be noted that both reduced operators share the same eigenvalue spectrum $\{\lambda_{i}\}$. These eigenvalues satisfy $0\leq\lambda_{i}\leq 1$ and $\sum_{i=1}^{\min(M,N)}\lambda_{i}=1$, as expected for a normalized pure state. For the indices $\min(M,N)<i\leq\max(M,N)$, the Schmidt eigenvalues $\lambda_{i}$ and the expansion coefficients $\sqrt{\lambda_{i}}$ are identically zero. In this work, we are using a coordinate grid representation for both the electronic and nuclear half-spaces. Then $\{x_{1},...,x_{N}\}\in\mathcal{H}_{el}$ and $\{R_{1},...,R_{M}\}\in\mathcal{H}_{nuc}$ and the Schmidt basis of the eigenstates are thus represented as $u_{i}(x_{j})=\langle x_{j}|u_{i}\rangle$ and $v_{i}(R_{j})=\langle R_{j}|v_{i}\rangle$.

II.3.2 Entanglement measure

The degree of entanglement can be quantified by the von Neumann entropy, which is defined as

[TABLE]

where $x$ stands for $U$ or $V$. When the expansion coefficient in Eq. (14) satisfies $\sqrt{\lambda_{i}}=\delta_{ij}$ the expansion reduces to a single product state term with a single eigenvalue $\lambda=1$ (that is, without entanglement) and the entropy $S$ is clearly zero. As the entanglement increases, so does the entropy, with larger values indicating substantial quantum correlations between subsystems in this particular half-space separation.

II.3.3 Wave functions in the Schmidt form

The BO, BH and the variational total wave functions may eventually be separable in electronic and nuclear coordinates, although they cannot be recast into a product state in the form $\Psi(\mathbf{x},\mathbf{R})=\phi(\mathbf{x})\chi(\mathbf{R})$. However, some degree of separability can be examined via the Schmidt decomposition. For example, in the case of a BO WF in Eq. (2), the position matrix representation of the reduced density matrices reads

[TABLE]

Instead, for a BH wave function, Eq. (12), the reduced density matrix elements now are

[TABLE]

and an analogous form for the variational total wave function. The Schmidt WF, obtained after diagonalizing the reduced nuclear and electronic density operators, is now represented as a sum of products of Schmidt eigenfunctions

[TABLE]

This expansion (which works for the BO wave function, the BH expansion and the variational ansatz) always contains some degree of entanglement between electronic and nuclear halfspaces and reaches a maximum entanglement entropy for $\sqrt{\lambda_{i}}=1/\sqrt{\mathrm{min}(M,N)}$, according to Eq. (15).

II.4 Discrete representation of the molecular Hamiltonian

To implement numerically the method of solution for both the electronic and the nuclear motions, several methods based on finite grids are available. We choose the Fourier grid Hamiltonian (FGH) method Marston1989, which belongs to the family of discrete variable representation (DVR) approaches Light1985. In the FGH method, a continuous spatial variable $r$ is discretized as

[TABLE]

with spacing $\Delta r$ for a total number of points $N$. The kinetic energy is represented in momentum space with a spacing given by $\Delta k=2\pi/(N\Delta r)$. The latter spacing in $k$-space defines the longest wave length (smallest frequency) that can be resolved. The Hamiltonian matrix for a single degree of freedom on this grid is then constructed as

[TABLE]

where $T_{n}=\frac{1}{2m}(n\Delta k)^{2}$ is the discretized kinetic energy term, and $V(r_{i})$ is the potential energy at each grid point. For a full derivation, see Ref. Marston1989. The eigenvectors resulting from this FGH are composed directly of the WF amplitudes at each grid point. Then $\psi=\sum_{i=1}^{N}\psi(r_{i})$ and the norm $\int dr\left|\psi(r)\right|^{2}=1$ in the FGH method is given by $\sum_{i=1}^{N}\left|\psi(r_{i})\right|^{2}=1$.

II.4.1 FGH for a molecular problem in one dimension

In this work (we believe without much loss of generality) we solve the electro-nuclear molecular problem for 1D systems in which both electrons and nuclei move in the same axis. This means that we neglect the rotation for the nuclear problem, and we propose simplified BO WFs in the form $\phi_{n}(x;R)\chi_{m,n}(R)$, where $R$ is just the vibrational coordinate.

Electronic Hamiltonian on a Grid.

The electronic Hamiltonian [Eq. (3)] (now in 1D) is solved for each fixed internuclear distance $R$, thus diagonalizing the FGH matrix, Eq. (22), in the electronic grid $\{x_{j}\}$. Note that the potential energy operator $\hat{V}$ depends parametrically on $R$. We are then provided with a set of electronic eigenfunctions $\phi_{n}(x_{j};R)$, along with the corresponding eigenvalues $E_{n}(R)$, which become the potential energy curves (PEC) for the 1D vibrational problem.

Nuclear Hamiltonian on a Grid.

The nuclear motion is solved with the nuclear Hamiltonian in Eq. (4) for 1D. We provide a grid also for the internuclear distances $\{R_{j}\}$ and the nuclear FGH Hamiltonian is diagonalized for the PEC of each electronic state $n$. Then we have at our disposal a set of vibrational eigenfunctions $\chi_{n,m}(R_{j})$ represented on the nuclear grid, along with their corresponding total vibronic energies $W_{m}^{n}$.

II.4.2 Reduced density matrices in the FGH discrete basis

In the FGH representation, the matrix elements for the electronic and nuclear reduced density matrices are now discretized within the 1D grids $\{x_{k}\}$ and $\{R_{i}\}$. For instance, for the case of BO WFs they are

[TABLE]

and similarly for the Born-Huang expansion and the total variational ansatz. Here we choose for the ongoing analysis the nuclear reduced density matrix in Eq. (23), which allows for a clearer interpretation.

III Results

Simplified models can provide valuable insight into the behavior of more complex molecular systems. A common approach involves reducing the dimensionality of the real system by focusing on the most relevant or active reaction coordinates. In this work, we study two illustrative low-dimensional models; the first is the hydrogen molecular ion in 1D and the second is the Shin-Metiu Hamiltonian model Shin1996. Both systems must be solved following the guidelines in Section II, and support solutions in BO form $\Psi^{\text{BO}}(x,R)$ and in BH form $\Psi^{\text{BH}}(x,R)$ in Eqs. (2) and (12), respectively. Electronic states in H${}^{+}_{2}$ hardly show avoided crossings within the manifold of the lowest excited states of the same symmetry (Wigner-von Neumann non-crossing rule). However, it is worth noting that a series of extremely broad avoided crossings in H${}^{+}_{2}$ among particular orbitals passed unnoticed for many years rost1989, but this will be of no concern here, since we limit ourselves to the two lowest electronic states for simplicity.

III.1 One-dimensional H${}_{2}^{+}$

Our first system considers the simplest molecule H${}_{2}^{+}$ in reduced dimensionality. Although the 3D case admits an exact solution in confocal elliptic coordinates Power1973, the exact analytical solution for 1D has not yet been worked out to our knowledge. However, this simplified model has been widely used for decades, for example, in the study of nonlinear multiphoton processes in molecules Kulander1996. Also, the plots and analysis of the electronic and nuclear WFs in 1D are much more simplified than in higher dimensions.

In this reduced model, only two coordinates, the internuclear separation $R$ and the electronic coordinate $x$, measured with respect to the center of mass of the nuclei, are taken into account along a single axis [see Fig. 1(a)].

In the associated Hamiltonian,

[TABLE]

$\mu$ is the reduced electron mass and $M_{p}$ is the proton mass. A softening parameter $a(R)$ is used for the Coulomb singularity. It is chosen to produce the exact three-dimensional $1s\sigma_{g}$ BO PEC Yue2014; Madsen2012. In this bipartite system, with three particles (nuclei with charges $Z_{\alpha}$ and $Z_{\beta}$ and one electron), the vibration and the electronic motion span the two Hilbert spaces $\mathcal{H}_{R}$ and $\mathcal{H}_{x}$.

The solution of the Schrödinger equation associated with the Hamiltonian in Eq. (25), follows the method described in Section II.4.1. We obtain electronic WFs in different nuclear configurations $R_{i}$ for the bonding orbitals $\phi_{1\sigma_{g}}(x;R)$ and $\phi_{2\sigma_{g}}(x;R)$, along with their associated vibrational states $\chi_{1\sigma_{g}/2\sigma_{g}v}(R)$, as illustrated in Fig. 2. The PEC for the electronic state $1\sigma_{g}$ has 18 vibrational bound states up to the dissociation limit at $E=-0.5$ a.u. while the PEC for the state $2\sigma_{g}$ supports 28 vibrational bound states up to $E=-0.23$ a.u. (note that $a(R)$ is only adjusted to produce the exact dissociation limit of the ground state).

For each electro-nuclear BO state, we have computed the entanglement content through the von Neumann entropy, either using Eq. (23) or Eq. (24) for the discrete FGH representation of the reduced density matrices and Eq. (15), with the obtained set of eigenvalues $\{\lambda_{i}\}$. The high amount of oscillations in the vibrational grid makes that in practice the number of points in the FGH representation of the nuclear grid be larger than the one in the electronic grid ($M>N$). In Fig. 4(a), we include the von Neumann entropy associated with the vibronic bound states of the electronic states $1\sigma_{g}$ and $2\sigma_{g}$. Both cases exhibit a monotonic increase in the amount of entanglement $S$ with the vibrational excitation $v$ and satisfy $S_{1\sigma_{g},v}<S_{2\sigma_{g},v}$. The latter effect is in principle related to the spatial extension of the vibrational WF. In general, the more spreading of the WF in the configuration space, the more delocalization and then its related entropy increases. In fact, above the dissociation threshold the vibronic WFs pertain to the nuclear continuum (states with infinite extension) and the entanglement must diverge for the exact continuum states. The entanglement for the electronic state $1\sigma_{g}$ has a convex curvature, whereas that for $2\sigma_{g}$ has a concave curvature. A more detailed explanation is given below using a simplified approach for the reduced densities in Eq. (24).

Contribution of molecular rotation to entanglement.

At this point, we may wonder how molecular rotation may affect the electro-vibrational entanglement. For simplicity, we may assume a 3D nuclear WF in the form $\chi^{J,M}_{n,m}({\bf R})=\chi^{J}_{n,m}(R)\mathcal{Y}^{M}_{J}(\theta,\phi)$, where $\mathcal{Y}^{M}_{J}(\theta,\phi)$ is a spherical harmonic function. For a rotating diatomic molecule, the reduced nuclear Schrödinger equation for the radial part $\chi^{J}_{n,m}(R)$ reads

[TABLE]

where $\mu_{N}$ denotes the nuclear reduced mass. The entanglement for the resulting rovibrational states is shown in Fig. 4(b), for rotational quantum numbers from $J=0$ to $J=15$. The blue dots represent all bound states of the $1\sigma_{g}$ PEC in the 1D system. The horizontal gray line connects states with the same vibration number $m$ but different rotational number $J$. Although there are shifts in energy induced by the additional centrifugal potential, the von Neumann entropy remains nearly unaffected for the lowest vibrational levels. In contrast, for high vibrational excitation, the entanglement contributed by rotation within the same vibrational level seems to also follow a convex pattern with increasing rotational excitation. In conclusion, our 1D model for H${}^{+}_{2}$ already captures the leading vibrational part of the entanglement.

III.1.1 Model with a simplified density matrix

We begin by simplifying the vibrational WF $\chi_{n,m}(R)$, defined on a discrete vibrational grid with points $\{R_{i}\}_{i=1}^{M}\in{\mathbf{F}}$, as obtained via the FGH method. All vibrational WFs exhibit oscillatory behavior, and we propose to reduce the representation of the FGH with a subset ${\mathbf{S}}\subset{\mathbf{F}}$ of grid points that correspond only to the location of maximum amplitudes in the WF $\{R_{k}\}_{k=0}^{m}\in{\mathbf{S}}$. The new amplitudes $\tilde{\chi}_{n,m}(R_{k})$ at each point $R_{k}$ are conveniently scaled to produce a new normalized state. The number of extrema (maxima and minima) in a vibrational wave function $\chi_{n,m}(R)$ with $m$ nodes is determined by the vibrational quantum number $m$, therefore with $m+1$ grid points in the new reduced set. Consequently, the simplified (normalized) vibrational WF $\tilde{\chi}_{n,m}(R)$ at the grid point $R_{k}$ reads

[TABLE]

Within the FGH grid, this construction ensures that the trace of the density matrix remains equal to unity, due to the renormalization of the quantum states. Fig. 3 shows the lowest three vibrational WFs ($\chi_{1\sigma_{g},m}$ and $\chi_{2\sigma_{g},m}$ for $m=0,1,2$) for the 1D H${}^{+}_{2}$ molecular system, along with their counterparts in our simplified model ($\tilde{\chi}_{1\sigma_{g},m}$ and $\tilde{\chi}_{2\sigma_{g},m}$ for $m=0,1,2$). Using Eq. (27), we obtain a much simplified expression for the density matrix, given by

[TABLE]

where $S_{n}(R_{k},R_{\ell})$ denotes the overlap integral between electronic WFs at nuclear distances $R_{k}$ and $R_{\ell}$, i.e. $S_{n}(R_{k},R_{\ell})=\int{dx\phi_{n}\left(x;R_{k}\right)\phi_{n}^{*}\left(x;R_{\ell}\right)}$. In this discrete representation, the dimension of the simplified density matrix corresponds to the number of extrema in the vibrational state. For example, the vibrational state $\tilde{\chi}_{n,0}$ leads to a density matrix of dimension unity, with zero entanglement by construction, while the state $\tilde{\chi}_{n,17}$ results in a dimension $18\times 18$ for the matrix.

The overlap integrals, $0\leq S_{n}(R_{k},R_{\ell})\leq 1$, thus introduce weights in the elements of the reduced density matrix. In the limit of perfect overlap, when $S_{n}(R_{k},R_{\ell})=1$, for all $R_{k}$ and $R_{\ell}$, it implies that the electronic WF is the same for all internuclear distances and the BO state is then fully separable and not entangled between electronic and nuclear motions. This is consistent with the fact that the reduced density matrix $\tilde{\rho}_{k,\ell}^{(R)}=\tilde{\chi}_{n,m}(R_{k})\tilde{\chi}_{n,m}(R_{\ell})$ and its eigenvalues are determined solely by the vibrational functions. This density matrix is a rank-one matrix, which means that the matrix is generated by the outer product of a vector with itself, ${\mathbf{M}}={\mathbf{x}}{\mathbf{x}}^{T}$, and has at most one non-zero eigenvalue $\lambda={\mathbf{x}}^{T}{\mathbf{x}}$ equal to the trace. Since vibrational states are normalized within the FGH grid, then $\text{Tr}\left[\tilde{\rho}^{(R)}\right]=1$. Consequently, we also have $\text{Tr}\left[\left(\tilde{\rho}^{(R)}\right)^{2}\right]=1$, the von Neumann entropy vanishes, and the fully separable BO state is not entangled. However, any molecule departs from this ideal situation.

Fig. 4 shows that the von Neumann entropy computed using the simplified density matrix closely matches the results obtained from the full density matrix; specifically, it also satisfies the asymptotic behavior at the dissociation thresholds. The subset with pairs of maxima and minima $\{R_{k},R_{\ell}\}\in{\mathbf{S}}$ in the vibrational WFs suffices to monitor the variation of the electronic WFs in those different nuclear geometries. For example, the lowest three vibrational WFs $\chi_{1\sigma_{g},n}(R)$ are mainly localized in the interval $1.7<R<2.7$, for which the electronic WFs barely change, as shown in Fig. 3(a)-(b) (the electronic overlaps at the grid points $\{R_{k},R_{\ell}\}$ remain close to unity).

For higher vibrational excitations, the WFs increase their spatial extension, and their extrema (or their nodes) are more separated from each other. Electronic overlaps noticeably depart from unity for many distant pairs $\{R_{k},R_{\ell}\}$ and these overlaps reduce the magnitude of the off diagonal density matrix elements. This effect is even more pronounced in the electronic excited state $2\sigma_{g}$ for which the distributions of pairs $\{R_{k},R_{\ell}\}$ are even more separated than in $1s\sigma_{g}$ and the overlap factor assesses very different electronic densities at two distant points.

In this way, nuclear WFs act as testers of the similarity of the electronic WFs at different geometries. The most relevant matrix elements contributing to the eigenvalue spectrum of the reduced density matrix are those that pertain to bands close to the diagonal; for a fixed row $k$, close to the diagonal, the electronic wave function is monitored with $S_{n}(R_{k},R_{\ell})$ at adjacent nuclear points for $R_{\ell}$, $\{R_{k-1},R_{k},R_{k+1}\}$. The worse the overlap in this tridiagonal matrix, the larger the entanglement since the eigenspectrum contains several sizeable eigenvalues $0<\lambda_{i}<1$.

Beyond the previous analysis, for large overlap values, $S_{n}(R_{k},R_{\ell})\sim 1$, we can introduce a perturbation to the reduced density rank-one matrix defined solely by the vibrational WFs. Then, the simplified density matrix from Eq. (28) can be recast as

[TABLE]

where $\epsilon_{l,\ell}^{(n,m)}$ is a perturbation defined as $\epsilon_{k,\ell}^{(n,m)}$ = $-\tilde{\chi}_{n,m}(R_{k})\tilde{\chi}_{n,m}(R_{\ell})\left[1-S_{n}(R_{k},R_{l})\right]$. In the case of perfect overlap, the factor $1-S_{n}(R_{k},R_{\ell})$ vanishes and $\tilde{\rho}^{(R)}$ reduces to the rank-one matrix. For non-overlapping electronic WFs, the largest eigenvalue of $\tilde{\rho}^{(R)}$ takes the form

[TABLE]

In such a case, for moderately entangled states, it is often sufficient to consider only the two largest eigenvalues of the reduced density matrix to provide a good estimation of the entanglement content. By just considering the two largest eigenvalues, the von Neumann entropy takes a simple form in comparison with the one listed in Eq. (15), given by

[TABLE]

This approximate von Neuman entanglement is included in Fig. 4(a) with black dotted lines. It works fine for entropies below 0.5 (for all bound vibrational states in $1s\sigma_{g}$ and up to the 8th vibrational state in curve $2\sigma_{g}$). From this perturbative approach, we conclude that even slight changes in the electronic overlap may induce (up/down) concavity in the entropy curves.

Compact vibrational states bounded by narrow PEC are expected to display similar electronic overlaps in the region between the two turning points. In this case, we find that the largest eigenvalue in the $2s\sigma_{g}$ state departs from 1 faster but still linearly with the vibrational excitation energy $E$, $\lambda_{max}^{(2\sigma_{g},E)}\sim 1-aE$ with a constant $a$. On the other hand, for compact vibrational states bounded by the narrow $1s\sigma_{g}$ PEC, we find a slower and nonlinear departure from unity with the vibrational excitation energy, which we can model with a nonlinear term $\lambda_{max}^{(1\sigma_{g},E)}\sim 1-aE^{2}$. These simple trends for $\lambda_{max}^{(n,m)}$, introduced in Eq. (31), explain the different concavities of the von Neuman entropies with increasing energy, as shown in Fig. 4.

Fig. 5(a) shows two BO wave functions for two vibrational states (the lowest $m=0$ and the last bound $m=17$) in the ground state $1\sigma_{g}$, $\phi_{1\sigma_{g}}(x,R)\chi_{1\sigma_{g},0}(R)$ and $\phi_{1\sigma_{g}}(x,R)\chi_{1\sigma_{g},17}(R)$, while Fig. 5(b) shows two BO wave functions also for the lowest and the upper vibrational states, $\phi_{2\sigma_{g}}(x,R)\chi_{2\sigma_{g},0}(R)$ and $\phi_{2\sigma_{g}}(x,R)\chi_{2\sigma_{g},27}(R)$ (see also Fig. 2). Here we use the whole set of grid points $\{R_{i}\}_{i=1}^{M}\in{\bf F}$. However, the full BO wave function can also be exactly reproduced with the complete set of Schmidt basis in Eq. (14) (eigenstates of the reduced density matrices). For the vibrational state $v=0$ in state $1\sigma_{g}$, the BO wave function is nearly separable since the first Schmidt coefficient $\lambda_{1}=0.998$ dominates, and the product of the lowest Schmidt electronic and nuclear eigenfunctions, $u_{1}(x)$ and $v_{1}(R)$, are sufficient to accurately describe the quantum state. A similar situation occurs for the lowest vibrational state in state $2\sigma_{g}$ although $\lambda_{1}$ now departs more from unity. All this reflects an almost complete separability for the ground state, thus a very low electron-nuclear entanglement and an almost vanishing von Neumann entropy. It should be noted that separability is exact and the entanglement entropy is zero for the same BO state represented in the simplified grid $\{R_{k}\}_{k=0}^{m}\in{\bf S}$ described above, because the ground state is represented by only a single grid point. In practice, the BO wave function reproduced with the Schmidt basis (in the electronic and nuclear halfspaces) using only the two highest eigenvalues is visually indistinguishable from those included in Fig. 5.

III.2 The Shin-Metiu model

In the previous section, we analyzed how changes in the configuration space affect the degree of entanglement in a molecule for which BO separability is appropriate. We now turn to examine how non-adiabatic effects, particularly those arising from avoided crossings between PECs, influence the entanglement content. The Shin-Metiu model has already been used in different scenarios, and it allows us to study the structure and dynamics through different kinds of PECs and avoided crossings, let them be sharp or broad.

The Shin-Metiu molecular model consists of three ions and a single electron in 1D [see Fig. 1(b)]. The ions labeled $Z_{\alpha}$ and $Z_{\beta}$ are fixed at a distance of $L=10\ \text{\AA }\approx 18.9$ a.u. The third ion, $Z_{\gamma}$, with nuclear mass $M$ is allowed to move along the $x$-axis between the two fixed ions, while the electron with mass $m_{e}$ moves freely along the same axis. The degrees of freedom are the position of the electron $x$, and the position $R$ of the ion with charge $Z_{\gamma}$ Seokmin1995; Shin1996, with respect to the origin in $x=0$. The Shin-Metiu Hamiltonian reads

[TABLE]

This model includes two Coulomb repulsion terms between the nuclei and three attractive electron-nucleus interactions. The latter are represented using screened Coulomb potentials in the form $V_{Ne}=Z_{i}\text{erf}\left(\hat{x}/R^{i}_{c}\right)/r_{i}$, where $r_{i}$ is the distance between the electron and the nucleus $i$, and $R_{c}^{i}$ is a cut-off parameter.

This model supports several stable nuclear configurations, corresponding to minima in the three adiabatic PECs. Adjusting charges and cut-off parameters $R_{c}^{i}$ allows us to adjust the electron’s ability to localize the ion with charge $Z_{i}$; different parameters can produce dramatic changes in adiabatic PECs and, in particular, in the strength and width of the avoided crossings (see Seokmin1995). Throughout this work, we use a linear chain of protons $Z_{\alpha}=Z_{\beta}=Z_{\gamma}=1$ and we select the cut-off parameters to be $R_{c}^{\alpha}=3.00$ a.u. $R_{c}^{\beta}=2.20$ a.u. and $R_{c}^{\gamma}=4.00$ a.u. This choice produces the adiabatic PECs $V^{A}_{1}(R)$, $V^{A}_{2}(R)$ and $V^{A}_{3}(R)$, shown in Fig. 6(a).

III.2.1 Adiabatic and Diabatic states

In this setup, the non-adiabatic couplings between the three lowest electronic states are non-negligible. As shown in Fig. 6(a), the system exhibits two different avoided crossings: a sharp one between the curves $V_{1}^{A}$ and $V_{2}^{A}$ located at $R^{(1,2)}=-3.0$ a.u., and a broader one between the curves $V_{2}^{A}$ and $V_{3}^{A}$ at $R^{(2,3)}=1.2$ a.u. At avoided crossings it is known that the BO approximation fails and the adiabatic states become coupled through non-adiabatic couplings. To force a BO approximation implies neglecting the non-adiabatic couplings. Strong non-adiabatic couplings cause character exchange between the states involved. For example, the electronic states $\phi^{A}_{1}$ and $\phi^{A}_{2}$ fully exchange their character after the sharp avoided crossing, and it can be observed in Fig. 7(b) for the narrow avoided crossing located at $R_{c}=-3.0$ a.u. in Fig. 6(a). However, the broad avoided crossing at a distance $R_{c}=1.25$ a.u. (see Fig. 6(a)) between the curves $V_{2}^{A}$ and $V_{3}^{A}$ produces little character exchange after the avoided crossing, as illustrated in Fig. 7(a). The presence of avoided crossings, non-adiabatic couplings and character exchanges indicates that the WF becomes a superposition of BO vibronic eigenstates, which ends up with a change in the electro-vibrational entanglement. For a briefing on the calculation of non-adiabatic couplings, see Appendix B.

A widely used alternative is the diabatic picture, brought about by a unitary transformation by which the non-adiabatic couplings vanish with a diabatic basis of states. The latter diabatic bases are not eigenstates of the electronic Hamiltonian, so that they are still coupled by electrostatic couplings, smaller and softer than non-adiabatic couplings. There is no single prescription to furnish a diabatic representation. Our approach to obtain diabatic states is described in more detail in the Appendix  A. The resulting diabatic electronic curves are shown in Fig. 6(b) (the avoided crossings become real crossings).

The electronic PECs and their ensuing vibrational states for the diabatic picture version of the Shin-Metiu system are illustrated in Fig. 6(b). After diabatization, the potential curves cross to each other and the vibrational eigenstates must be obtained from the nuclear equation that includes the new diabatic PECs. Then, above the energy of any energy crossing, the adiabatic and diabatic vibrational states differ. Consequently, in the diabatic picture the character exchange is now reduced to a minimum [see Fig. 8(b)] but electronic diabatic states may differ from adiabatic ones [compare Fig. 8(a) with Fig. 7(a)].

III.2.2 Electro-nuclear entanglement in the Shin-Metiu model

Adiabatic and diabatic pictures

The von Neumann entropy for the lowest three potential energy curves of the Shin-Metiu model, in both adiabatic and diabatic pictures, is shown in Fig. 9. These entanglement entropies are computed using the eigenvalues of both the full reduced density matrix and the reduced density obtained with the simplified model introduced in the previous section for H${}^{+}_{2}$. The dimension of the simplified reduced density matrix is $m+1$, with $m$ being the label of the vibrational state or its number of nodes.

From Fig. 9, we realize that the entanglement of the BO vibronic states in the adiabatic picture for PECs $V^{A}_{1}$ and $V^{A}_{3}$ follows the same pattern already observed in H${}^{+}_{2}$ for the $1\sigma_{g}$ and $2\sigma_{g}$ states, respectively. The states in the PEC $V^{A}_{1}$ are not coupled below the lowest sharp avoided crossing at $E=-0.30$ a.u. and those in the PEC $V^{A}_{3}$ are barely coupled due to a wider avoided crossing with a weak coupling. The electronic WF of the electronic ground state $\phi_{1}(x;R)$ varies little from $R=-2.5$ to $R=6.0$ a.u., so that the two-point electronic overlap $S_{1}(R_{k},R_{\ell})$ remains close to unity for any pair. However, a sudden increase in the entanglement for $V^{A}_{1}$ occurs in the energy region around the lowest avoided crossing at distance $R_{c}^{(1,2)}=-3.0$ a.u. and energy $E\sim-0.3$ a.u. In H${}^{+}_{2}$ an asymptotic increase at $E=-0.5$ a.u. occurs at the dissociation threshold and the ensuing delocalization of the nuclear WF.

However, in the Shin-Metiu model, the entanglement has a smooth behavior except for the regions of avoided crossings, where it shows the form of a step function. Note that in the avoided crossing zone there is a character exchange in the electronic WF; also the configuration space for the vibration changes. Both effects contribute to increasing the entanglement content by up to $0.6$ for BO vibronic states $\phi_{1}(x;R)\chi_{1,m}(R)$ in the PEC $V^{A}_{1}$ at the lowest avoided crossing, and from $0.6$ to $1.1$ for states $\phi_{2}(x;R)\chi_{2,m}(R)$ in the PEC $V^{A}_{2}$. Above the energy $E=-0.30$ a.u. the entanglement for $V^{A}_{1}$ increases monotonically but very slowly. Taking into account the character exchange in the electronic WF, the overlap $S_{1}(R_{k},R_{\ell})$ now departs from unity for pairs of nuclear distances $\{R_{k},R_{\ell}\}$ that meet $R_{k}<R_{c}^{(1,2)}<R_{\ell}$ or $R_{\ell}<R_{c}^{(1,2)}<R_{k}$. Apart from this exchange character, the electronic WFs change little when both distances $\{R_{k},R_{\ell}\}$ lie above $R_{c}^{(1,2)}$ or below it. Then the entanglement above the avoided crossing remains similar in magnitude, as shown in Fig. 9(a). The BO states $\phi_{2}(x;R)\chi_{2,m}(R)$ in the PEC $V^{A}_{2}$ lying between the two avoided crossings [$E=-0.30$ and $E=-0.25$ a.u. approximately in Fig. 6(a)] increase rapidly their entanglement for the same reason. However, although the electronic WF $\phi_{2}(x;R)$ has a character exchange in $R_{c}^{(1,2)}$ the configurational space accessible to vibration that evaluates different overlaps $S_{2}(R_{k},R_{\ell})$ is much reduced.

This case study indicates that the compactness in space (localization) of wave functions is not a sufficient condition to reduce the entanglement. Strong electron-nuclear correlations at a given geometry may produce a dramatic change in entanglement. In fact, even though BO PECs $V^{A}_{1}$ and $V^{A}_{3}$ are similar in shape, their lowest tens of vibronic states strongly differ in their entanglement, the former with a small electron-nuclear correlation that facilitates separability and the latter with an exchange of character across an avoided crossing.

Our simplified model for the reduced density matrix reproduces very well the entanglement behavior for the BO vibronic states $\phi_{2}(x;R)\chi_{2,m}(R)$ and $\phi_{3}(x;R)\chi_{3,m}(R)$, but for $\phi_{1}(x;R)\chi_{1,m}(R)$ the simplification fails above the lowest avoided crossing at $E=-0.30$ a.u. We attribute this failure to insufficient sampling of nuclear points $\{R_{k},R_{\ell}\}\in{\bf S}$, when the crossings are close to the turning points in the PECs. In this region nuclear WFs are characterized by the last wide oscillation with the highest amplitude. A single grid point located at the maximum cannot capture the abrupt changes in the electronic WF caused by a nearby avoided crossing. However, the simplified model nicely follows the trend of the entanglement behavior calculated with the full grid.

In the diabatic picture [see Fig. 9(b)], the magnitude of the von Neumman entanglement for the vibronic states in the PEC $V_{1}^{D}$ and $V_{2}^{D}$ is much reduced and with two small steps starting at the potential crossings. The diabatization procedure generates a smooth continuation of the electronic WFs $\phi_{1}(x,R)$ and $\phi_{2}(x,R)$ across the sharp avoided crossing $R_{c}^{(1,2)}$ and therefore the electronic overlaps $S_{n}(R_{k},R_{\ell})$ become closer to unity. In addition, WFs $\phi_{2}(x,R)$ and $\phi_{3}(x,R)$ are even smoother across the avoided crossing $R_{c}^{(2,3)}$. The entanglement of vibronic states in the diabatic PEC $V_{3}^{D}$ is not much affected and, indeed, only slightly reduced compared to the adiabatic picture. The only explanation comes from the fact that the diabatic electronic WF $\phi^{D}_{3}$ varies with nuclear geometry more than $\phi^{D}_{1}$ and $\phi^{D}_{2}$. This can be visualized in Fig. 8(a)-(b); while the diabatic electronic WFs $\phi^{D}_{1}$ and $\phi^{D}_{2}$ remain localized and similar at two separate points $R_{c}^{1,2}=-3.0$ and $R_{c}^{2,3}=1.2$, the WF $\phi^{D}_{3}$ changes noticeably. A similar pattern was found in H${}^{+}_{2}$ for state $2\sigma_{g}$.

Born-Huang expansion

The Born–Huang expansion is believed to offer a systematic way to improve the total molecular state beyond the adiabatic BO approximation by proposing a total wave function as a superposition of electronic BO states in Eq. (12). Non-adiabatic couplings between the different electronic states in the expansion are now explicitly included. The diagonalization of Eq. (II.2), in the FGH grid provides the total vibronic energies and states. Having access to Born-Huang solutions in a realistic molecular system may become prohibitive. This Shin-Metiu model allows for the accurate computation of all ingredients (electronic and vibrational wave functions and the whole set of non-adiabatic couplings).

Although the expressions for the full density matrix in this basis [Eqs. (18) and (19)] are more elaborate than those derived for pure BO states [Eqs. (23) and (24)], its application is straightforward. Within the BH method the vibronic eigenstates cannot be attributed to a single potential curve (beyond BO there are no unambiguous PECs). Despite that, Fig. 9(c) illustrates that the entanglement for most BH vibronic states aligns closely along one of the BO entropy values for $V^{A}_{1}$, $V^{A}_{2}$ and $V^{A}_{3}$. This happens when the mixing component of one BO wave function $\phi^{A}_{n}$ dominates. Nevertheless, there is always a BO superposition, so that the BH entanglement is always slightly above any BO counterpart. Incidentally, above $E=-0.30$ eV a couple of two nearly degenerated BO vibronic states from $V^{A}_{1}$ and $V^{A}_{2}$ with energies $E_{1,m_{1}}>E_{2,m_{2}}$, interact with a (approximately) constant non-adiabatic coupling to produce a splitting in the energies and a new coupled BH vibronic eigenstate where the leading component is the $V^{A}_{1}$ component, then inherited by the entanglement content. The opposite also occurs for $E_{1,m_{1}}<E_{2,m_{2}}$ with a dominant contribution from $V^{A}_{2}$.

However, Fig. 9(c) also displays many other BH vibronic states whose entanglement deviates from the BO results and they correspond to states with significant mixing between multiple BO states, leading to an entropy value larger than that of the individual contributing BO states. These BH states typically appear in nearly degenerate doublets or triplets, a feature that reflects their origin in strongly coupled sets of BO vibronic states through the non-adiabatic couplings. As an illustration, we highlight the case of BH vibronic states 95 and 96 (in the energy region of the second avoided crossing), which exhibit the highest mixture between two distinct BO vibronic states associated with different PECs $V^{A}_{1}$ and $V^{A}_{2}$. In Fig. 10(a)-(b), we plot the probability densities for the two BH vibronic states separated by a small energy splitting and with the superposition $0.4880\phi_{1}\chi_{1,72}$ + $0.5118\phi_{2}\chi_{2,23}$ for the lower BH vibronic state and $0.5117\phi_{1}\chi_{1,72}$ + $0.4881\phi_{2}\chi_{2,23}$, for the upper BH vibronic state. Both states have almost the maximum mixing 0.5/0.5 and their entropies are almost identical. These BH states are built from two energy-degenerated BO vibronic states, whose densities are plotted in Fig. 10(c)-(d), clearly inherited by the BH states. Above $E=-0.25$ a.u. there is a plethora of BH vibronic states with an entanglement that clearly departs from the BO discipline, and they correspond to other complex superpositions of three BO states, this time $\phi_{2}$ and $\phi_{3}$ having the leading role.

Entanglement in the total molecular eigenfunction.

The total energies and WFs were also obtained by direct exact diagonalization of the Shin-Metiu Hamiltonian in Eq. (32). We use the same FGH grid basis mentioned above for the electronic and nuclear coordinates, in tensor product form. The purpose of this calculation is to provide a reference to be compared to those previous solutions obtained within the adiabatic (both the BO and BH ansatzes) and diabatic frameworks.

With the set of screening parameters used in this work, the Shin–Metiu model allows for the analysis of two distinct non-adiabatic regimes within the BO picture: a sharp avoided crossing (strong coupling) and a broad avoided crossing (weak coupling). Fig. 9(d) presents the von Neumann entropy of the variational vibronic states (with colored dots). The colored solid lines correspond to the entanglement computed for the adiabatic and diabatic states discussed previously [note that we have also included the entanglement corresponding to the adiabatic 4th state with PEC $V_{4}^{A}(R)$]. To clarify the observed trend, we performed the diabatization procedure only for the lowest sharp avoided crossing $V_{1}-V_{2}$ (entanglement of vibronic states in $V^{D}_{1}$ with red line and of vibronic states in PEC $V^{D}_{2}(r)$ with cyan line), but not for the broad avoided crossing $V_{2}-V_{3}$ (vibronic states in $V^{D}_{2}$ with blue line and vibronic states in PEC $V^{D}_{3}(R)$ with green line).

To facilitate the analysis we plot the total nuclear propability densities for each total vibronic state $n$. They are obtained from the Schmidt decomposition (20) of the total variational wave function (13)

[TABLE]

Fig. 11 displays the nuclear probability densities of all non-BO variational vibronic states with energies below $E=-0.15$ a.u. (higher states are omitted for clarity). We realize that below the lowest avoided crossing at $E\sim-0.30$ a.u. the nuclear density of vibronic states follow the envelope of the adiabatic PEC $V^{A}_{1}(R)$, as expected, so do their related entropies of entanglement. Some near degenerate vibronic states around this crossing exhibit strong mixing. Above this sharp avoided crossing, the nuclear densities mainly follow the diabatic picture with vibronic densities alternating between PECs $V^{D}_{1}(R)$ (red curve) and $V^{D}_{2}(R)$ (cyan curve). In contrast, the upper avoided crossing centered around $E~\sim-0.25$ a.u.largely preserves the adiabatic following of the nuclear densities within each respective PEC $V^{A}_{2}(R)$ and $V^{A}_{3}(R)$, although some degree of mixing always remains. The densities in solid gray lines in Fig. 11 correspond to vibronic states that follow the adiabatic curve $V_{4}^{A}(R)$, which barely exhibit non-adiabatic mixing, so does the entanglement in Fig. 9.

In Fig. 9(d), the color of each point corresponds to the nuclear density displayed in Fig. 11. Although the total states are not purely adiabatic or diabatic, their classical turning points follow either adiabatic or diabatic PECs. The nuclear densities associated with red points mainly follow the diabatic curve $V_{1}^{D}(R)$, the cyan ones follow $V_{2}^{D}(R)$, and the green and gray points correspond to adiabatic curves $V_{3}^{A}(R)$ and $V_{4}^{A}(R)$, respectively. Appendix C illustrates the comparison between two quasi-degenerated total variational eigenstates (labelled 82 and 83 in Fig. 11) with the corresponding (closer in energy) adiabatic and diabatic states. The better comparison of the eigenstates with the diabatic states endorses our entanglement findings and supports the chemical intuition for a diabatization procedure when necessary.

Entanglement seems to be a very subtle property of multipartite WFs and it provides indeed a witness to their accuracy. Our total variational eigenfunction (and its entanglement content) represents the molecular Shin-Metiu model at best, and provides a guidance on the quality of other approximated WFs. In terms of entanglement, the adiabatic BO approximation works fine for states far from avoided crossings (low mixing). In this way entanglement content for vibronic states within the ground state below the crossing, and for all vibronic states within the 4th state is well represented by BO WFs. Instead, our entanglement results dictate that the vibronic states within the 1st and 2nd electronic states must be treated diabatically above the lowest avoided crossing, while in the region of the second avoided crossing and above, the vibronic states must be better treated adiabatically. It is worth noting that the BH expansion based on the BO states does not provide the correct diabatic trend when needed, thus indicating a slow convergence for the proper representation of the total WF.

IV Conclusions

We have investigated the content of electro-nuclear entanglement in molecular wave functions. The quantum correlations between particles in a complex system like a molecule are already encoded within the Hamiltonian. A good representation of the correlations depends on the choice of a variational ansatz to solve the problem. Each variational molecular WF is then subjected to an alternative analysis in terms of quantum entanglement, which is not equivalent to correlation.

Without loss of generality, we have analyzed the quantum entanglement between the electronic and nuclear motions present in variational WFs of molecular models in 1D. Although the BO approximation is sometimes invoked to claim that the electronic motion is separated from the nuclear one, the BO WF is indeed a kind of conditional product state $\Psi^{BO}({\bf r},{\bf R})$ = $\phi^{BO}({\bf r};{\bf R})\chi^{BO}({\bf R})$ but not a fully separable product state $\phi({\bf r})\chi({\bf R})$, so it is rigorously a quantum entangled state. In this work, we have evaluated the electro-nuclear entanglement of BO WFs in H${}^{+}_{2}$ for two lowest electronic states of symmetry $\sigma_{g}$ and their associated vibrational states. From this example, we learn that, contrary to the intuitive expectation that a larger configurational space (delocalization of the WF) leads to a larger entropy of entanglement, the variation of the electronic (nuclear) WF along the nuclear (electronic) distances plays a central role in the discrimination of molecular entanglement. Interestingly, our present results for the entanglement in 1D H${}_{2}^{+}$ compare very well with those previously obtained with a non-BO ansatz to solve the 3D H${}_{2}^{+}$ (states $1s\sigma_{g}$ and $3d\sigma_{g}$) Sanz2017. Our results also show that the inclusion of molecular rotation does not change the entanglement content dramatically, indicating that the leading contribution comes in fact from the vibrational degrees of freedom.

Along this direction, we demonstrate that it is fair to consider the BO WF for the vibronic ground state to be almost separable and unentangled. However, the electro-nuclear entanglement increases with vibrational excitation. An extended excited-vibrational WF checks out more nuclear configurations in which differences in the electronic WF may appear. Since this proves to be a quantitative origin of entanglement, we have proposed a simplified model for the nuclear reduced density matrix in which the information of the vibrational WF is encoded by retaining only the points with maximum amplitude, and we show that these nuclear positions are fairly enough to assess the variation of electronic WF. This simplification provides a computationally inexpensive method for estimating the entanglement as well as a simple perturbative approach for understanding the curvature behavior of entanglement at low vibrational excitations. We also apply this framework to study abrupt changes in entanglement near avoided crossings in the PECs. The Shin-Metiu molecular model in 1D contains many of the expected ingredients in the photodynamics of molecules involving excited states: non-adiabatic transitions, adiabatic to diabatic pictures, etc. In the BO adiabatic picture, sharp avoided crossings produce rapid changes in entanglement as a result of abrupt modifications in the electronic WF (character exchange). In contrast, broad avoided crossings lead to smoother increases in entropy, as the electronic WF evolves smoothly across nuclear configurations.

In the diabatic picture (where non-adiabatic couplings vanish) the change of entanglement across the real crossings is very much reduced. This result can be puzzling, but the entanglement of a bipartite system is indeed not unique; it depends on how the halfspaces are initially separated in the total WF. For example, in our previous work Sanz2017 we explored the entanglement of two coupled oscillators, a system that admits an exact solution in terms of normal modes. In the latter case, the total WF is a non-entangled direct product. In other local coordinates (which are related to the normal modes by a simple rotation), the WF becomes entangled. Similarly, the adiabatic and diabatic pictures are related by a basis transformation. Thus, the diabatic picture provides an ansatz for which electrons and nuclei display a much reduced entanglement, and we find that it eventually produces the correct representation for sharp avoided crossing (when compared to the total eigenfunction).

If we move beyond the BO approximation, a proposal is to expand the total WF in terms of the BO adiabatic electronic states, and to include explicitly the non-adiabatic couplings (BH expansion). In general, the entropy of entanglement of a BH vibronic state closely follows that of the dominant contributing BO electronic state in the expansion. However, in certain cases, significant mixing between BO vibronic states (localized on different PECs and quasi-degenerated in energy) leads to a noticeable increase in the entanglement of the BH vibronic state, beyond that of any BO contribution. In perspective, in the search for a maximum entanglement in vibronic states beyond the BO approximation, the entanglement content for each vibronic state might be analyzed individually.

For a deep understanding of the results obtained within the BO approximation, we computed the entanglement with a variational eigenfunction of the total Hamiltonian. As mentioned above, from the entanglement results we conclude that the BO approximation and the BH expansion in terms of BO WFs do not adequately describe the correct WF in the presence of sharp avoided crossings, in favor of a diabatic representation. In contrast, broad avoided crossings with feeble couplings are reasonably described within the BO approximation and BH expansions. Our study here is performed on systems of reduced dimensionality, but we believe that our methods and conclusions remain valid for realistic molecules of higher dimensionality.

Appendix A Adiabatic and diabatic states

There is no unique way to describe the diabatic representation of a coupled system. One approach involves the direct computation of the non-adiabatic coupling matrix elements (e.g., via the Hellmann–Feynman theorem), while another option is to perform successive unitary rotations in the vicinity of each avoided crossing Smith1969; Alden1982; Lecomte1979; Triana2022. These two methods yield similar results and are theoretically equivalent Shu2022. In the present work, we adopt the latter approach for convenience.

In the case of three adiabatic PECs, the corresponding matrix obtained by solving the Schrödinger equation [Eq. (1)]. Within the BO approximation, it is block diagonal and it reads

[TABLE]

Since this matrix is already diagonal, there exists a unitary matrix $\mathbf{U}_{T}$ such that $\mathbf{V}^{D}=\mathbf{U}_{T}^{\dagger}\mathbf{V}^{A}\mathbf{U}_{T}$, where $\mathbf{V}^{D}$ contains the diabatic PECs. In the system studied here (with two avoided crossings located at $R_{c}^{(1,2)}$ and $R_{c}^{(2,3)}$), the diabatization procedure consists of a unitary transformation, in which a total rotation matrix is built from two consecutive rotations: $\mathbf{U}_{T}=\mathbf{U}_{1}\mathbf{U}_{2}$. The first rotation, characterized by an angle $\theta(\mathbf{R})$, couples $V_{1}^{A}$ and $V_{2}^{A}$ near $R_{c}^{(1,2)}$, while the second rotation, characterized by an angle $\varphi(\mathbf{R})$, couples $V_{2}^{A}$ and $V_{3}^{A}$ near $R_{c}^{(2,3)}$. Independent rotations are justified because of the spatial separation of the avoided crossings in our Shin-Metiu model.

The diagonal elements of the diabatic matrix are then given by

[TABLE]

The off-diagonal elements (electrostatic couplings, since the diabatic states are not eigenstates of the electronic Hamiltonian) are

[TABLE]

The angles $\theta(R)$ and $\varphi(R)$ can be modeled using error functions of the form

[TABLE]

where $K$, $K_{0}$, and $\Gamma$ are parameters optimized for each angle to minimize the deviation from the conditions defined by Eq. (35). The rotation angles must also satisfy the conditions below

[TABLE]

Appendix B Non-adiabatic couplings

B.1 First order non-adiabatic couplings

A direct method for computing the non-adiabatic couplings $\mathbf{A}_{n,n^{\prime}}$, as defined in Eq. (7), relies on the off-diagonal Hellmann-Feynman theorem Balasubramanian1990, that gives

[TABLE]

for $n\neq n^{\prime}$, and we take $\mathbf{A}_{n,n}(R)=0$. However, this method can become numerically unstable near degeneracies between two PECs. To overcome such issues, several alternative approaches have been proposed Mandal2018.

In our approach, we avoid direct computation of couplings by instead deriving an analytical expression based on the error function introduced above. This function governs the rotation angles $\theta(R)$ and $\varphi(R)$ between adiabatic and diabatic states as described in Appendix A and connects with the non-adiabatic couplings as

[TABLE]

Since the error function has an analytical derivative, the corresponding expressions are given by

[TABLE]

Fig. 12 shows the agreement between the solution via Hellmann-Feynman and the couplings derived through angle rotations near the avoided crossings. Minor discrepancies are observed in $A_{1,2}(R)$ near the second avoided crossing at $R_{c}^{(2,3)}$, as seen in Fig. 12(b). However, it becomes irrelevant since the vibrational wave functions $\chi_{2,m}(R)$ vanish for $R>0$ (see Fig. 6). Thus, this model can be effectively used for systems with well-separated avoided crossings, where the adiabatic to diabatic transformation can be represented by successive, independent rotations.

B.2 Second order non-adiabatic couplings

The second-order couplings $B_{n,n^{\prime}}(R)$, involving second derivatives with respect to $R$, are typically more difficult to compute and often neglected. From the structure of Eq. (7), the second-order term can be calculated using a completeness relation for the electronic states, with

[TABLE]

Since we treat each avoided crossing independently, it suffices to construct two $2\times 2$ matrices, to deal with the non-adiabatic coupling in the vicinity of $R_{c}^{(1,2)}$ and $R_{c}^{(2,3)}$, respectively.

For the sharp avoided crossing at $R_{c}^{(1,2)}$, the matrices take the form

[TABLE]

and

[TABLE]

An equivalent representation holds for the second avoided crossing at $R_{c}^{(2,3)}$ using instead $\varphi(R)$.

Appendix C Total eigenfunctions vs. adiabatic and diabatic WFs in the Shin-Metiu model

Here we include plots for the two quasi-degenerated variational states of the Shin-Metiu Hamiltonian (states 82 and 83 with energies $E=-0.2558$ and -0.2556 a.u.), with the corresponding (closest in energy) adiabatic states $\phi^{A}_{1}(x;R)\chi^{A}_{1,66}(R)$ ($E=-0.2561$ a.u.) and $\phi^{A}_{2}(x;R)\chi^{A}_{2,16}(R)$ ($E=-0.2553$ a.u.), respectively, and the closer diabatic states $\phi^{D}_{1}(x;R)\chi^{D}_{1,61}(R)$ ($E=-0.2558$ a.u.) and $\phi^{D}_{2}(x;R)\chi^{D}_{2,21}(R)$ ($E=-0.2556$ a.u.), respectively. In this energy region above the sharp avoided crossing, the WFs in the diabatic picture are those that better reproduce the total variational eigenfuntions. This WF similarity is then translated into their respective electron-nuclei entanglement.

IV DATA AVAILABITY

All the numerical results shown in this work can be obtained straightforwardly from the formulas and parameters described in the text. However, the data are available from the authors upon reasonable request.

Acknowledgements.

J.F.P.M and J.L.S-V acknowledge financial support from Vicerrectoría de Investigación at Universidad de Antioquia, CODI, Programática Project 2022-5357 and Estrategia de Sostenibilidad.