We study the influence of nodal structures in two-dimensional quantum mechanical densities on wave packet entanglement. This is motivated by our recent study [Entropy, 25, 970 (2023)], which showed that the mutual information derived from the momentum-space probability density of a coupled two-particle system exhibits an unusual time dependence, which is not encountered if the position-space density is employed in the calculation. In studying a model density, here, we identify cases where the mutual information increases with the number of nodes in the wave function and approaches a finite value, whereas in this limit, the linear correlation vanishes. The results of the analytical model are then applied to interpret the correlation measures for coupled electron-nuclear dynamics, which are treated by numerically solving the time-dependent Schrödinger equation.

Nodes appearing in quantum mechanical wave functions are closely related to the excitation state of the system under consideration. For example, the number of nodes in the radial functions of the hydrogen atom equals the principal quantum number plus one.1 In addition, regarding normal modes of vibrations, the nodes, being directly related to the quantum numbers, determine the symmetry selection rules for infrared transitions.2 Another aspect of zeros that are present in probability densities is that they cause problems in applying classical approximations to the density, e.g., in quantum diffusion calculations.3,4

Nodal structures in wave functions influence the appearance of different observables if detected as a function of a single parameter. For example, using the “reflection principle,”5 structures appearing in absorption spectra can be traced back to the nodes of the initial vibrational wave functions,6 and the same applies to resonance Raman spectra.7 As another example, we mention the power of time-resolved photoelectron spectroscopy to detect quantum topologies in vibrational wave packets.8 In this paper, we are concerned with densities from which the correlation between two variables can be inferred. The focus is to investigate how nodes, being present in such densities, are related to the correlation. Although our numerical example treats the wave-packet entanglement between electronic and nuclear degrees of freedom, the derived results are more general, and they apply also to other types of correlation effects. This, in particular, is the case for many-body electronic systems where, due to their fermionic character, nodes arise naturally. Electron–electron correlation plays a prominent role in quantum chemistry, and it is thus important to characterize such collective effects. Electronic structure calculations are static and do not yield information on the time evolution of the correlation. Here, we treat dynamical effects that arise from the correlated motion of electrons and nuclei. Because the solution of the time-dependent Schrödinger equation for molecules is, in general, not possible without imposing restrictive approximations, a simple numerical model is employed. The results may serve as a starting point for further work on more complex systems. Within the employed model, it is possible to determine electron-nuclear wave packets in position- and momentum-space. The respective densities exhibit distinct quantum structures and, at certain times, are characterized by different node patterns. We address the following question: how do these time-dependent patterns influence particle correlations? This is of importance regarding molecules because chemical bonding and reactivity are determined by such correlations. It is then necessary to find reasonable ways to characterize the latter.

There are different measures that are used to estimate entanglement; for applications in the context of chemical reactions and electronic structure theory, see, e.g., Refs. 9–20. We consider two of such functions, namely, the mutual information (MI)21,22 and the linear correlation coefficient,23 and investigate their time dependence arising from nuclear-electron dynamics.

The mutual information has been calculated to characterize correlation effects in atoms and molecules.11,14,15,19 It is important, however, to emphasize that the MI and the correlation coefficient are no observables but measures to visualize and quantify correlation.

The paper is organized as follows: In Sec. II, we summarize the main definitions. The results are given in Sec. III, where analytical expressions for the correlation measures are derived for a parameterized model density (Sec. III A). A numerical example is presented in Sec. III B. The latter treats the adiabatic motion of an electron and a nucleus, which are coupled through screened Coulomb interactions. Using the analytical results, the behavior of the MI and the correlation coefficients derived from the coordinate-space and momentum-space densities, respectively, is interpreted. The article is closed with a summary given in Sec. IV.

Our considerations are limited to two variables x and y and a probability density ρ(x, y), which may also depend on time. Although this describes a more general situation, we will, in what follows, refer to a nucleus and an electron having nuclear (nuc) and electronic (el) degrees of freedom of x and y, respectively.

The non-linear correlation between the two particles is measured using the mutual information (MI).21,22 It is a functional that is related to the difference between the probability density ρ(x, y) and the density built from the single particles, which are defined as
ρnuc(x)=dyρ(x,y),
(1)
ρel(y)=dxρ(x,y).
(2)
Then the MI is
I=dxdyρ(x,y)lnρ(x,y)ρnuc(x)ρel(y).
(3)
Equivalently, this function can be written as
I=Snuc+SelS.
(4)
Here appear the (differential) Shannon entropies24,25 determined from the total density,
S=dxdyρ(x,y)ln[ρ(x,y)],
(5)
and the single-particle entropies,
Snuc=dxρnuc(x)ln[ρnuc(x)],Sel=dyρel(y)ln[ρel(y)].
(6)
If the Shannon entropy is regarded as a measure of information, the MI can be interpreted as the amount of information lost when comparing what is inferred from the densities ρnuc(x) and ρel(y) with the information extracted from the total density ρ(x, y).
As another measure for entanglement, we define the linear correlation coefficient as23 
corr=covσxσy,
(7)
where the nominator is the covariance function,
cov=xyxy,
(8)
and the denominator contains the variances,
σx2=x2x2,σy2=y2y2.
(9)

In the latter equations, the brackets denote the expectation values taken with respect to the density ρ(x, y). The correlation coefficient takes a value between −1 and 1, and its absolute value determines how well the density ρ(x, y) is approximated by a straight line.23 This means that, approximately, we find a linear relationship ⟨y⟩ ≈ λ0 + λ1x⟩, with λ0 and λ1 being numbers. The sign determines if this line has a negative or positive gradient. For a vanishing function, i.e., corr = 0, no linear correlation between the variables exists, but higher order correlations may be present in the MI.26 

In this section, we start from a density of a particular form, and the reason for its choice will become apparent when the numerical example presented in Sec. III B is discussed. The density is
ρa,b(x,y)=2e(x2+y2)cos2a1+b2(x+by)π(1+ea2),
(10)
where a and b are real valued parameters. From the total density, the single-particle densities are determined using MATHEMATICA,27 and the results are
ρa,bnuc(x)=ea2b2+1x2ea2b2b2+1+cos2axb2+1πea2+1,
(11)
ρa,bel(y)=ea2b2b2+1y2ea2b2+1+cos2abyb2+1πea2+1.
(12)

The total density is a Gaussian multiplied by a cos2 function, and it is normalized to a value of one. The number of nodes that fall into the region where the Gaussian is located can be controlled by the variation in the parameter a, where in the limit of a → 0, a standard non-correlated Gaussian is obtained. For non-zero values of the second parameter b, the Gaussian is rotated in the (x, y) plane, and the rotation angle is given by α = tan−1[b]. The factor 1+b2, which is included in the argument of the cos2 function, ensures that b does not distort the overall shape of the density.

Examples for densities ρa,b(x, y) obtained for different choices of the parameters a and b are displayed in Fig. 1. Also included are the respective nuclear densities ρa,bnuc(x) (blue lines) and electronic densities ρa,bel(y) (red lines). The upper panels contain functions for a fixed value of a = 7 and different values of b taken from the interval b ∈ [0, 1]. As it is seen in the figure, with the variation in b, the nodal lines change their orientation with respect to the x-axis, but the frequency with which they appear is not influenced. In the case of b = 1, the nodes are parallel to the x = −y diagonal. Here, electronic and nuclear densities are nodeless and of equal shape.

FIG. 1.

Densities ρa,b(x, y), as defined in Eq. (10). The upper panels show functions for a = 7 and selected values of the parameter b. In the lower panel, b = 0.2 is used, and a is varied, as indicated. Also shown are nuclear particle densities ρa,bnuc(x) (blue lines) and electronic particle densities ρa,bel(y) (red lines).

FIG. 1.

Densities ρa,b(x, y), as defined in Eq. (10). The upper panels show functions for a = 7 and selected values of the parameter b. In the lower panel, b = 0.2 is used, and a is varied, as indicated. Also shown are nuclear particle densities ρa,bnuc(x) (blue lines) and electronic particle densities ρa,bel(y) (red lines).

Close modal

The influence of the parameter a on the density is illustrated in the lower panels of Fig. 1, where we set b = 0.2. With increasing values of a, more nodes appear in the density. It is also seen that the nuclear particle density exhibits structures that can be traced back to the nodes, whereas for the present choice of parameters, the electronic particle densities do not, and they appear to be nearly identical.

The mutual information Ia,b is calculated numerically from the densities given in Eqs. (10)(12) as a function of the two parameters a and b. Representative curves are shown in the upper panel of Fig. 2. The curves are displayed as a function of the frequency factor a and for selected values of b. Note that for b = 0, the total density factorizes and the MI vanishes so that Ia,0 = 0, which is independent of the value of a. All curves start at zero and grow monotonically. They approach an upper bound at Ia,b ≈ 0.31, and this limit is reached faster with the increasing value of b. The functions illustrate that the MI is directly related to the number of nodes in the density. Note that this number is not a well-defined concept for the model density ρa,b(x, y) because, technically, there exist an infinite number of zeros. As mentioned before, here, we mean the zeros that fall into the region where the Gaussian envelope is of non-negligible intensity. Approximately, it holds that Ia,b is proportional to this number, where the proportionality factor itself is proportional to b. The asymptotic value for a and b = 1 can be calculated analytically, see the  Appendix, and one finds that
limaIa,1=1ln[2]0.3069,
(13)
which excellently agrees with the numerical results. The latter suggests that this limit does not depend on the particular value of b.
FIG. 2.

Mutual information (upper panel) and correlation (lower panel) as a function of the parameter a for selected values of b, as indicated. The black line in the upper panel corresponds to the analytical limit of a (b = 1).

FIG. 2.

Mutual information (upper panel) and correlation (lower panel) as a function of the parameter a for selected values of b, as indicated. The black line in the upper panel corresponds to the analytical limit of a (b = 1).

Close modal
We next calculate the covariance (again using MATHEMATICA) with the result
cova,b=a2bea2+1b2+1.
(14)
The variances in the two variables are determined as
σx2=12a2ea2+1b2+1,
(15)
σy2=12a2b2ea2+1b2+1.
(16)
Both functions are bounded by a value of 0.5, and they are identical for b = 1. By inserting the obtained results into the definition of the linear correlation coefficient, one finds that
corra,b=a2bea2+1b2+1×a22ea2+1+a4b2ea2+12b2+12+1412.
(17)
In the considered parameter regime (b ϵ [0, 1]), the covariance and correlation are negative semidefinite. By considering them as a function of the parameter a, i.e., as a function of the number of nodes, they exhibit a single minimum (obtained by symbolic differentiation employing MATHEMATICA) at
amin=W1e+11.13,
(18)
where W is the Lambert W-function.28 Both the covariance and correlation are zero in the limit of vanishing parameters a and b, respectively. This is expected because in both cases, the density equals an uncorrelated Gaussian. In addition, we find that they also approach zero for a.

The functions corra,b are shown in the lower panel of Fig. 2 for different values of b. As predicted, extrema are found at amin ≈ 1.13 for all values of b. They occur when the wavelength of the oscillations is comparable to the variance of the density. This approximately is the case for a value of a = 2, as can be seen in Fig. 1. There, the Gaussian shape-function suppresses all maxima of the cos2 function except the central one. This has the consequence that the density is located in a restricted region around the line with orientation angle tan−1 (b = 0.2). Furthermore, the smaller the value of b, the smaller the value of corra,b because the orientation angle diminishes. In addition, if compared to the MI, the correlation approaches its asymptotic limit already for lower values of a.

What we have shown so far is that, for a model Gaussian density ρa,b(x, y), which is modulated with a periodic function, the MI and correlation coefficient behave rather differently as a function of the number of occurring nodes. The correlation coefficient shows a single extremum and quickly reaches a value of zero. This feature does not depend on the orientation of the nodal lines in the (x, y)-plane. On the other hand, the MI is a monotonous function of the number of nodes, and it converges, for all values of b ≠ 0, to a finite limit of ln[e/2]. This limit is reached for much larger values of a than the asymptotic value of zero reached by corra,b.

We now treat coupled electron-nuclear wave packet dynamics as a non-trivial example, where the analytical results of Sec. III A can be used for interpretation. The model consists of an electron (coordinate r) and a proton (coordinate R) moving in one dimension.29,30 They interact with each other and two fixed nuclei at positions of R1,2 = ±5 Å through screened Coulomb interactions. The interaction potential is (in atomic units)
V(r,R)=1|R1R|+1|R2R|erf[|R1r|/Rf]|R1r|erf[|Rr|/Rc]|Rr|erf[|R2r|/Rf]|R2r|+Δ,
(19)
where erf denotes the error function. For the screening parameters, we use Rf = 1.5 Å and Rc = 1.0 Å, and Δ is the energy shift adjusted such that the potential minimum occurs at zero energy in the region of our spatial grid.
The time-dependent Schrödinger equation for the coupled motion reads
itΨ(r,R,t)=ĤΨ(r,R,t),
(20)
with the Hamiltonian
Ĥ=p̂R22M+p̂r22me+V(r,R).
(21)
The momentum operators for the electron and the nucleus are denoted as p̂r and p̂R, respectively, the mass of the proton is M, and me is the electron mass.
We solve Eq. (20) numerically with the split-operator method.31 The grid ranges from −12 to +12 Å and from −6 to +6 Å for the coordinates r and R, respectively. The number of grid points is taken as 512 for both coordinates, and a time step of Δt = 0.0024 fs is employed. The initial condition is chosen as
Ψ(r,R,0)=N0eβ02(RR0)2φ0(r;R),
(22)
with the normalization factor N0, and the Gaussian is parameterized with R0 = −3.5 Å and β0 = 7.14 Å−2. The function φ0(rR) is the electronic ground state wave function, which parametrically depends on the nuclear coordinate, and it is the solution of the electronic Schrödinger equation
pr̂22me+V(r,R)φ0(r,R)=V0(R)φ0(r;R),
(23)
where V0(R) is the ground state potential curve. Equation (23) is solved numerically using imaginary time propagation.32 

In Sec. II, the mutual information and the correlation are defined with respect to the variables y and x. We first identify these variables with the coordinates r and R, respectively, and the density ρ(r, R, t) = |Ψ(r, R, t)|2 in the coordinate space is obtained from the time-dependent wave function. Taking the Fourier transform of the latter, one arrives at the momentum-space wave function Ψ(k, K, t), where y and x now correspond to the electronic (k) and nuclear (K) momenta, respectively, and the density is ρ(k, K, t) = |Ψ(k, K, t)|2. From the two densities, the MI in coordinate space (Iξ(t)) and in momentum space (Iπ(t)) is readily calculated. Likewise, the correlations corrξ(t) and corrπ(t) are determined.

Instead of using either position-space or momentum-space densities, one can think of using the Wigner function to obtain the MI from a phase-space distribution. This is not straightforward because it may assume negative values33 so that the logarithm cannot be calculated. There have been attempts to circumvent the associated difficulties. We will not follow this route and refer the reader to the literature.34,35

For the given parameterization of the interaction potential and the chosen initial wave function, the electronic ground state potential V0(R) is separated by a large energy gap from the excited state potentials.36 Accordingly, the Born–Oppenheimer (BO) approximation applies, and the wave function can be written as a product of φ0(rR) and the nuclear wave function ψ(R, t) depending on time. In a previous work on the present model,37,38 we found that information theoretical properties derived from coordinate dependent densities can be predicted analytically if the numerically propagated wave function is approximated as
Ψ(r,R,t)=Nteβt2(RRt)2eγ2(rR)2,
(24)
which is of the BO form. Here appear the normalization factor Nt, the Gaussian width parameter βt, and the center of the Gaussian Rt, where the subscripts indicate that these parameters depend on time. However, significant deviations were found in the momentum space MI38 so that its time dependence cannot be explained within the ansatz of Eq. (24). One reason is that the BO Gaussian function does not incorporate nodal structures of the nuclear-electron wave packet. In what follows, we will take up this point and address their importance.
The time-dependent mutual information obtained for times up to 300 fs is displayed in the upper two panels of Fig. 3. It is seen that the spatial MI exhibits an overall increase, superposed by oscillations that are related to the vibrational wave packet dynamics. Adopting the BO ansatz of Eq. (24), it can be shown37 that the coordinate space MI takes the form
Iξappr(t)=12ln1+2γσR2(t),
(25)
where the time-dependent variance σR2(t) enters the nuclear coordinate. This curve is also included in Fig. 3. The values for σR2(t) are extracted from the numerical propagation, and an average value of γ = 0.733 Å−2 is used, which is determined from the electronic wave function φ0(rR). It is seen that Iξappr(t) reproduces the numerically obtained MI accurately, emphasizing the quality of the Gaussian approximation.
FIG. 3.

Mutual information (upper two panels) and correlation coefficients (lower two panels) calculated from the coordinate-space (subscript ξ) and momentum-space densities (subscript π). The numerically derived curves (num) are compared to the ones obtained from the simplified ansatz for the BO wave function [Eqs. (25), (26), (28) and (29)]. In addition, the MI determined from Eq. (27) is shown in the two upper panels (Gauss). The red crosses correspond to functions derived within the harmonic ansatz (HO) for the nuclear density, as described in the text.

FIG. 3.

Mutual information (upper two panels) and correlation coefficients (lower two panels) calculated from the coordinate-space (subscript ξ) and momentum-space densities (subscript π). The numerically derived curves (num) are compared to the ones obtained from the simplified ansatz for the BO wave function [Eqs. (25), (26), (28) and (29)]. In addition, the MI determined from Eq. (27) is shown in the two upper panels (Gauss). The red crosses correspond to functions derived within the harmonic ansatz (HO) for the nuclear density, as described in the text.

Close modal
The function corrξ(t) is shown in Fig. 3, second panel from below. It increases as a function of time and converges to a value of 0.9 after about 100 fs. The curve is compared to the analytical function derived from the BO ansatz, which reads37 
corrξappr(t)=σR2(t)σR2(t)+12γ.
(26)
From the figure, it emerges that this formula presents an excellent approximation to the numerical result.
There is a relationship between the linear correlation coefficient and the mutual information for two-dimensional Gaussians,24,39 which is
IGauss=12ln1corr2.
(27)
The curve for IξGauss(t), determined from the numerically derived function corrξ(t) and using Eq. (27), is also included in Fig. 3. Inserting Eq. (26) in Eq. (27) yields the result of Eq. (25). However, the numerical curves show minor deviations, which are due to the Gaussian approximation. These findings mean that the coordinate-space MI is mainly determined by linear correlations between the two moving particles. In particular, because the wave function in Eq. (24) does not include nodal patterns but is just the product of two correlated Gaussians, such nodes do not influence the spatial MI and the spatial correlation coefficient.
The numerically calculated momentum space MI is displayed in Fig. 3 (second panel from the top). In addition, we show the function that is derived from the analytical ansatz, which reads38 
Iπappr(t)=12ln1γ2σK2(t),
(28)
with the nuclear variance in momentum space being σK2(t). The latter is determined numerically from the propagated wave function Ψ(k, K, t), and γ is chosen as specified above. The curves agree only for short times (this is not seen in the figure because of their low magnitude), and afterward, no agreement is found. The MI seems to saturate for later times at values slightly above a value of Iπ ≈ 0.30, where the appearing fluctuations are not due to numerical errors. We note that this value is close to what is found for the MI derived from the model density (Sec. III A) in the limit of a large number of nodes.
The momentum-space correlation is shown in the lower panel of Fig. 3. It exhibits some oscillations related to the nuclear dynamics and, for longer times, remains at a value of zero. The approximate expression for the correlation coefficient, which is38,
corrπappr(t)=γ2σK2(t),
(29)
is able to reproduce the main features of the numerical results. Thus, here, the assumption of an unstructured Gaussian wave packet is appropriate. We emphasize that the values of the various functions obtained for longer times do not necessarily remain constant over time. In particular, if the level structure of the quantum system allows for a wave-packet revival,40–44 the initial values will be recovered.

From the discussion of the results presented in Fig. 3, it emerges that the time-evolution of the position-space MI and correlation as well as the momentum-space correlation is not sensitive to an eventually existing node pattern in the probability densities. Only the momentum-space MI cannot be explained within an ansatz of nodeless Gaussians. In analyzing the density ρ(k, K, t), one finds that the deviations between the numerically and approximately calculated MI are accompanied by the appearance of nodal patterns. In particular, the maxima seen in Iπ(t) occur when clear nodal lines are detectable. To illustrate these patterns, we show in Fig. 4 the coordinate- and momentum-space densities for two selected times. It is seen that in position space, the nodes are oriented vertically (visualized by a green line) and the density is oriented along the diagonal (r = R), as is indicated by the red line. Due to wave-packet dispersion, the density becomes more de-localized over time, and more structures are seen. In momentum space, the nodes appear parallel to the line, defined by k = −K, and the density is located along the horizontal axis with k = 0. As is the case in position space, a richer node pattern is encountered at a later time, and the density becomes more extended along the horizontal axis. This explains why the correlation coefficient is small (see Fig. 3).

FIG. 4.

Coordinate-space (upper panels) and momentum-space densities (lower panels) of the electron-nuclear wave packets shown for times of t = 41 fs and t = 201 fs. As a guide to the eye, the red and green lines indicate the overall orientation of the densities and the nodal lines, respectively.

FIG. 4.

Coordinate-space (upper panels) and momentum-space densities (lower panels) of the electron-nuclear wave packets shown for times of t = 41 fs and t = 201 fs. As a guide to the eye, the red and green lines indicate the overall orientation of the densities and the nodal lines, respectively.

Close modal

The densities shown in Fig. 4 exhibit characteristics that we encountered in those of the model system studied in Sec. III A. Regarding the density ρa,b(x, y)) with the parameter b = 0, the nodal lines are oriented vertically (see Fig. 1), as is also seen in the spatial densities of Fig. 4, and for b ≠ 0, they are rotated in the plane, as is found for the momentum densities. In the model, the b = 0 case leads to an uncorrelated wave function so that both, the MI and the correlation coefficient, are zero. In our numerical example, the orientation of the position space density is along the diagonal and not along one of the axis, which can be traced back to the particular form of the BO wave function. This leads to a non-zero MI and correlation, with the latter approaching the value of one, which are not (or only weakly) influenced by the vertical nodes.

To find an appropriate approximation to what is found numerically and explicitly address the role of the nodes, we proceed as follows: The ground state electronic eigenfunction, which enters the Born–Oppenheimer ansatz [Eq. (24)], is, to a good approximation a Gaussian for all values of R,37 and the dynamics are described by temporal changes in the nuclear wave packet ψ(R, t). Let us assume that the motion of this wave packet is harmonic where the harmonic eigenfunctions are denoted as χm(R). Expansion of the nuclear density then yields
ρnuc(R,t)=mcmeiEmtχm(R)2,
(30)
with coefficients cm and harmonic eigenenergies Em. In order to get a hand on the number of nodes in the nuclear density, we assume that at a time ts, the density is dominated by a single term in Eq. (30) (with quantum number m = n) so that it can be written as
ρnuc(R,ts)|χn(R)|2.
(31)
The normalized density then is
ρnHO(r,R,ts)=2nβtsγπn!eβtsR2γ(Rr)2Hn(βtsR)2,
(32)
with the Hermite polynomials Hn. Examples for these spatial densities are shown in the upper panels of Fig. 5 for different values of the quantum number n. Here, the values of β = γ = 1 Å−2 are employed. It is seen that the displayed functions have the same properties as the numerically determined ones shown in Fig. 4. That is to say, they are oriented along the diagonal and show vertical nodes. From the findings presented in Sec. III A, it is expected that both Iξ(t) and corrξ(t) are not influenced by such oriented nodes. In addition, the densities become more extended along the diagonal, which implies a linear relationship between r and R, and thus, the MI and the correlation here contain the same information (see below).
FIG. 5.

Representative Born–Oppenheimer densities involving the harmonic approximation for the nuclear motion. The upper and lower panels depict the coordinate-space and momentum-space densities as defined in Eqs. (32) and (33), respectively. Shown are functions for different quantum numbers n. The nuclear (blue lines) and electronic (red lines) particle densities are also included.

FIG. 5.

Representative Born–Oppenheimer densities involving the harmonic approximation for the nuclear motion. The upper and lower panels depict the coordinate-space and momentum-space densities as defined in Eqs. (32) and (33), respectively. Shown are functions for different quantum numbers n. The nuclear (blue lines) and electronic (red lines) particle densities are also included.

Close modal
To arrive at the momentum space densities, the Fourier transformation of the coordinate space wave function is calculated, which then leads to the density as follows:
ρHO(k,K,ts)=2nπn!βtsγe1γk212βts(K+k)2×Hn1βts(K+k)2.
(33)

The relation between the electronic and nuclear degrees of freedom is different from what is encountered in coordinate space. In particular, the Hermite polynomials, which are responsible for the nodes, are dependent on the sum of the electron and nuclear momenta.

Momentum space densities are shown in the bottom panels of Fig. 5. They very much reproduce the characteristics of the numerical densities depicted in Fig. 4. The orientation of the overall density and of the nodal lines are the same. With the increasing value of the quantum number n, the density spreads along the k = 0 axis, whereas the width in the direction of the electron momentum is constant. This has a consequence that the correlation becomes smaller with increasing values of n. On the other hand, because the nodal lines follow the direction of the line k = −K, they are not visible in the single particle densities so that this information encoded in the total density is lost and the MI (which measures the loss of information) is non-zero. This is in line with the discussion about the properties of the model density in Sec. III A.

Although the harmonic-like densities reproduce features of the numerically obtained electron-nuclear densities, it is unclear if they are useful in describing the correlation measures we are interested in. To answer this question, we choose a set of times t = ts and estimate the number of nodes n appearing in the spatial density ρ(r, R, ts). This is possible only approximately because the densities do not always show clearly defined nodes (in particular, in momentum-space). Rather, more or less pronounced minima are seen. However, the MI and correlations are not too sensitive with respect to the exact number of nodes. The variances of the position- and momentum-space densities, as defined in Eqs. (32) and (33), can be calculated analytically, and they are
σR2(n)=2n+12βs,
(34)
σK2(n)=2n+12βts+γ2.
(35)

Inserting the numerically derived variances into the latter equations, i.e., σR2(n)=σR2(ts) and σK2(n)=σK(ts), leads to two different values of the parameter β (namely, βR and βK), which are then used in the equations for the approximate densities. As before, a value of γ = 0.733 Å−2 is chosen. In this way, we obtain approximate densities in position- and in momentum-space, which have the widths and the number of nodes similar to the numerical ones. The thus constructed functions are used to calculate the MI and correlation coefficients, and the results are included in Fig. 3 as red crosses. In what follows, they are referred to as IξHO,IπHO,corrξHO and corrπHO, respectively. It is seen that in position-space, the approximate MI (IξHO) and correlation (corrξHO) agree well with the numerical curves. The same applies to the momentum-space correlation (corrπHO). In all three cases, the nodes of the wave functions do not play a role because the Gaussian ansatz, which neglects the node structure, yields about the same results. The only exception is the function Iπ(t). There, the assumption of a Gaussian form of the wave function fails completely. Including the nodes as described above yields a curve that deviates at shorter times but gives a reasonable prediction of the overall time dependence. The deviations occur at times when the wave packet is fairly localized, so the representation of its R-dependence by a single harmonic wave function is not very accurate. The curve increases and, for longer times, converges to a constant. This trend is in accordance with the results derived for the model density shown in Sec. III A. Thus, we have found that the momentum-space MI reflects non-linear electron-nuclear correlations, which are strongly influenced by the nodal structure of the underlying momentum density.

Let us now discuss the limit of longer times in the regarded interval. The numerically calculated spatial- and momentum-densities for a time of tn = 281 fs are displayed in the left hand panels of Fig. 6. Counting the number of nodes in the coordinate space density yields an approximate value of n = 65. The respective densities obtained within the harmonic approximation [i.e., Eqs. (32) and (33)] are also shown in the figure (right hand panels). They are obtained for the parameters βR = 13.63 Å−2, βK = 3.01 au−2, and γ = 0.733 Å−2. As can be seen in the figure, they reproduce the main characteristics of the numerical densities. This opens up the possibility to directly study the influence of the number of nodes on the MI and the correlation coefficients at the considered fixed time. These functions are calculated by keeping the chosen parameters fixed (i.e., βR, βK, and γ) and varying the quantum number n of the harmonic functions. In the upper two panels of Fig. 7, the MI in position- and momentum-space is shown. Both curves grow monotonically in the displayed range of quantum numbers, while the gradient decreases slowly. They start at the same value because for zero nodes, the correlation of the model is linear. The momentum-space MI increases with a large gradient at lower quantum numbers and then levels slowly. This is exactly what is found in the b = 1 case, which corresponds to the orientation encountered here along the line k = −K of the model density ρa,b(x, y) shown in Sec. III A. However, it does not converge in the presented regime. Taking the estimated value of n = 65, marked in Fig. 7 as a vertical line, we find the following numbers, Iξ(281) = 1.07, IξHO=0.95, Iπ(281) = 0.31, and IπHO = 0.28, so that good agreement is found. Figure 7 also contains the MI as obtained from Eq. (27) and using the correlations depicted in the lower panels of the figure. In position space, this is a good approximation, but it completely fails in momentum space. There, the curve approaches zero because the linear correlation vanishes. This again emphasizes that in calculating the function Iπ, the nodes of the quantum wave functions are of great importance.

FIG. 6.

Comparison of the propagated coordinate-space and momentum-space densities (left hand panels) and their harmonic approximations (right hand panels) at a time of t = ts = 281 fs.

FIG. 6.

Comparison of the propagated coordinate-space and momentum-space densities (left hand panels) and their harmonic approximations (right hand panels) at a time of t = ts = 281 fs.

Close modal
FIG. 7.

Mutual information (upper two panels) and linear correlation coefficients in position- and momentum space (lower two panels). The curves (HOn, interpolated by straight lines) are calculated using the densities ρξHO and ρπHO constructed for a time of ts = 281 fs, as described in the text. They are shown for fixed parameters as a function of the quantum number n. For comparison, the MI determined from Eq. (27) using the correlations depicted in the two lower panels is also included (red curves). The horizontal line indicates the quantum number n = 65, which is determined from the numerical derived density at t = 281 fs.

FIG. 7.

Mutual information (upper two panels) and linear correlation coefficients in position- and momentum space (lower two panels). The curves (HOn, interpolated by straight lines) are calculated using the densities ρξHO and ρπHO constructed for a time of ts = 281 fs, as described in the text. They are shown for fixed parameters as a function of the quantum number n. For comparison, the MI determined from Eq. (27) using the correlations depicted in the two lower panels is also included (red curves). The horizontal line indicates the quantum number n = 65, which is determined from the numerical derived density at t = 281 fs.

Close modal

Note that the correlations corrξ(t) and corrπ(t) do not exhibit extrema, as is found for the model density ρa,b(x, y) (see Fig. 2). This is because the variance of the densities is larger than the wavelength of the oscillations for all quantum numbers. Thus, here, we are in a regime of values above the critical number of amin [Eq. (13)].

The focus of this work is on the analysis of how nodal structures in quantum wave functions influence particle correlations. Therefore, we consider two measures of such an entanglement. One is the mutual information, and the other is the linear correlation coefficient. These functions are derived from probability densities for systems with two degrees of freedom.

We first analyze model densities with a Gaussian shape function, where nodes are introduced by a cosine function. Here, two parameters (a and b) enter. The first one (a) determines the frequency with which the nodes appear, and the other (b) defines the orientation of the nodal lines in the plane. It is shown that the MI behaves as a monotonously increasing function converging to the same limit for all values of b. The latter determines the rate with which the asymptotic value is reached. On the other hand, the linear correlation coefficient shows an extremum for all values of b and approaches zero for larger numbers of the frequency. This means that, opposite to the behavior of the MI, the latter is only sensitive to the number of nodes within a region where the wavelength of the oscillations in the density is larger or comparable to the variance of the density.

Having the analytical results at hand, we apply them to the motion in a system consisting of an electron and a nucleus, which show the typical dynamics encountered in the Born–Oppenheimer regime. In integrating the time-dependent Schrödinger equation for the coupled motion, the MI and the correlation coefficients are determined as a function of time. The numerical results are compared to those derived using a Gaussian ansatz for the BO wave function, which does not account for the node patterns in the probability density. It is found that the correlation measures derived from the coordinate-space density show excellent agreement with the numerically exact results. This hints at the fact that, here, the nodal structure of the density does not play an important role. As is inferred from the analytical considerations, the reason is that nodal lines present in the spatial densities are oriented perpendicular to the electronic coordinate axis. This is related to the BO wave function where the dynamics and thus the node structure are contained in the nuclear wave function whereas the electronic wave function remains of approximately constant shape.

Another picture evolves from the momentum-space dynamics. There, the quantum density exhibits nodal lines, which are oriented and non-parallel to the nuclear momentum axis. Here, the particular form of the BO wave function puts these lines at an angle of −π/4 to the axis. This has the consequence that the MI strongly depends on the number of nodes. Because the linear correlation coefficient does not show this dependence, we identify non-linear contributions to the particle correlation which, however, are only seen if the momentum-space density is used in the calculation. To explicitly address the influence of the number of nodes, we calculate the correlation measures, assuming harmonic wave functions for the nuclear degree of freedom at different times. The then obtained densities in position- and momentum-space are a good approximation to the numerically derived ones at selected times. By increasing the number of zeros in the nuclear harmonic density, it is documented that, for the regarded BO dynamics where the electron adiabatically adapts to the nuclear geometry, it is only the momentum-space mutual information that is sensitive to structures in the quantum mechanical wave function. We thus conclude that these nodal structures increase the wave packet entanglement in a non-linear fashion, which can be visualized by calculating the mutual information between the particles.

Finally, we point out that even though the nodal structures in our example influence only the momentum-space MI, one can find examples where the nodes strongly influence both the momentum- and the position-space MI. However, the results presented here illustrate the general effect, which, to the best of our knowledge, has not been discussed in detail before.

The authors have no conflicts to disclose.

Peter Schürger: Conceptualization (equal); Formal analysis (lead); Software (lead); Visualization (lead); Writing – original draft (equal). Volker Engel: Conceptualization (equal); Formal analysis (supporting); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

In what follows, we provide a proof of Eq. (13), which gives the value of the MI in the limit of large values of the parameter a. This is equivalent to a large number of nodes in the model density appearing in Eq. (10). Here, we treat the special case where the parameter b assumes the value of b = 1. The density then reads
ρa(x,y)=2ex2y2cos2(ax+ay)πe2a2+1,
(A1)
where we have, for simplicity of notation, the replacement aa2.
In the case of b = 1, because the nodes are aligned along the line x = −y, the electronic and nuclear densities have the same functional dependence on their respective coordinates. Regarding (exemplarily) the nuclear density and using symbolic integration, its analytic expression is determined as
ρanuc(x)=ea2x2ea2+cos(2ax)πe2a2+1.
(A2)
We next apply a coordinate transformation as (x, y) → (xy, y) so that the density takes the form
ρa(x,y)2e(xy)2y2cos2(ax)πe2a2+1.
(A3)
This simplifies the integration over the variable y, which is performed with MATHEMATICA and yields
dxf(x,a)=dxdyρa(x,y)lnρa(x,y)=dxe2a2x22cos2(ax)2πe2a2+1×4a22lnπe2a2+1lnsec4(ax)x21+ln(4).
(A4)
In order to obtain the MI, the remaining integrals over the variable x are to be solved, i.e.,
Ia=dx2ρanuc(x)lnρanuc(x)f(x,a).
(A5)
The integrand appearing in the latter equation can be decomposed into eight different terms Tn(x, a) as
2ρanuc(x)lnρanuc(x)f(x,a)=n=18Tn(x,a),
(A6)
with
T1(x,a)=22πa2e2a2x22cos2(ax)e2a2+1,
(A7)
T2(x,a)=e2a2x22cos2(ax)2πe2a2+1,
(A8)
T3(x,a)=x2e2a2x22cos2(ax)2πe2a2+1,
(A9)
T4(x,a)=ln(4)e2a2x22cos2(ax)2πe2a2+1,
(A10)
T5(x,a)=2πe2a2x22cos2(ax)e2a2+1lnπe2a2+1,
(A11)
T6(x,a)=2e2a2x2πe2a2+1lnea2x2ea2+cos(2ax)πe2a2+1,
(A12)
T7(x,a)=2ea2x2cos(2ax)πe2a2+1lnea2x2ea2+cos(2ax)πe2a2+1,
(A13)
T8(x,a)=e2a2x22cos2(ax)lnsec4(ax)2πe2a2+1.
(A14)
Direct symbolic integration can be performed for the terms T1(x, a), T2(x, a), T3(x, a), T4(x, a), and T5(x, a), and one finds
dxT1(x,a)=2a2,
(A15)
dxT2(x,a)=12,
(A16)
dxT3(x,a)=a2a2tanha212,
(A17)
dxT4(x,a)=ln(2),
(A18)
dxT5(x,a)=lnπe2a2+1.
(A19)
The terms T6(x, a) and T7(x, a) are obtained by first using
limadxTn(x,a)=dxlimaTn(x,a).
(A20)
We therefore calculate these limits, which are
limaT6(x,a)=ex22x2+ln(π)π,
(A21)
limaT7(x,a)=0.
(A22)
Thus, T7(x, a) vanishes and does not contribute to the MI in the considered limit. The term of Eq. (A21) can be integrated, which yields
dxlimaT6(x,a)=1+lnπ.
(A23)
For the last term T8(x, a), MATHEMATICA does not provide a solution, so we need to analyze it in more detail. It can be written as
T8(x,a)=e2a2e2a2+1ex222πcos2(ax)lncos4(ax).
(A24)
First we note that
limae2a2e2a2+1=1.
(A25)
The remaining part being present in T8(x, y) consists of a product of a Gaussian and an oscillating function having a period of P = π/a. Thus, for large values of a, the oscillations become arbitrarily fast if compared to the scale on which the Gaussian changes, and, effectively, the Gaussian only sees an average value, which is determined as
1P0P=π/adxcos2(ax)lncos4(ax)=ln(4)1.
(A26)
Inserting this average value and the limit Eq. (A25) into the integral containing T8(x, a) yields
dxlimaT8(x,a)=dxex222π[1ln(4)]=1ln(4).
(A27)
We now sum over all contributions to obtain
n=18dxTn(x,a)=2a212+a2a2tanha212+ln(2)lnπe2a2+1+1+ln(π)+1ln(4),
(A28)
which simplifies to
n=18dxTn(x,a)=3a2ln2e2a2+2a2tanha2+1.
(A29)
This equation indeed converges in the limit of a, and the asymptotic value of the MI is calculated as
limaIa=liman=18dx Tn(x,a)=1ln(2)=lne20.306852.
(A30)
1.
E.
Merzbacher
,
Quantum Mechanics
(
Wiley
,
New York
,
1998
).
2.
S.
Califano
,
Vibrational States
(
Wiley
,
New York
,
1976
).
3.
P. J.
Reynolds
,
D. M.
Ceperley
,
B. J.
Alder
, and
W. A.
Lester
,
J. Chem. Phys.
77
,
5593
(
1982
).
4.
J. B.
Anderson
, in
Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy
, edited by
S. R.
Langhoff
(
Springer
,
Netherlands, Dordrecht
,
1995
), pp.
1
45
.
5.
J.
Tellinghuisen
, in
Photodissociation and Photoionization
, 3rd ed., edited by
K. P.
Lawley
(
Wiley
,
New York
,
1985
).
6.
R.
Schinke
,
Photodissociation Dynamics
(
Cambridge University Press
,
Cambridge
,
1993
).
7.
M.
Eichelsdörfer
and
V.
Engel
,
Chem. Phys. Lett.
263
,
640
(
1996
).
8.
H.
Katsuki
,
H.
Chiba
,
B.
Girard
,
C.
Meier
, and
K.
Ohmori
,
Science
311
,
1589
(
2006
).
9.
W.
Kutzelnigg
,
G.
Del Re
, and
G.
Berthier
,
Phys. Rev.
172
,
49
(
1968
).
10.
R. F. W.
Bader
and
M. E.
Stephens
,
J. Am. Chem. Soc.
97
,
7391
(
1975
).
11.
R.
Grobe
,
K.
Rzazewski
, and
J. H.
Eberly
,
J. Phys. B: At., Mol. Opt. Phys.
27
,
L503
(
1994
).
12.
P.
Ziesche
,
Int. J. Quantum Chem.
56
,
363
(
1995
).
13.
R. P.
Sagar
and
N. L.
Guevara
,
J. Chem. Phys.
123
,
044108
(
2005
).
14.
R. P.
Sagar
and
N. L.
Guevara
,
J. Chem. Phys.
124
,
134101
(
2006
).
15.
R. F.
Nalewajski
,
Found. Chem.
16
,
27
(
2014
).
16.
E. V.
Ludeña
,
F. J.
Torres
,
M.
Becerra
,
L.
Rincón
, and
S.
Liu
,
J. Phys. Chem. A
124
,
386
(
2020
).
17.
S.
López-Rosa
,
R. O.
Esquivel
,
J. C.
Angulo
,
J.
Antolín
,
J. S.
Dehesa
, and
N.
Flores-Gallegos
,
J. Chem. Theory Comput.
6
,
145
(
2010
).
18.
S.
López-Rosa
,
A. L.
Martín
,
J.
Antolín
, and
J. C.
Angulo
,
Int. J. Quantum Chem.
119
,
e25861
(
2019
).
19.
J. C.
Angulo
and
S.
López-Rosa
,
Entropy
24
,
233
(
2022
).
20.
R. O.
Esquivel
,
M.
Molina-Espíritu
, and
S.
López-Rosa
,
J. Phys. Chem. A
127
,
6159
(
2023
).
21.
C. E.
Shannon
,
Bell Syst. Tech. J.
27
,
379
(
1948
).
22.
S. J. C.
Salazar
,
H. G.
Laguna
, and
R. P.
Sagar
,
Eur. Phys. J. Plus
137
,
19
(
2022
).
23.
E.
Mandel
and
E.
Wolf
,
Optical Coherence and Quantum Optics
(
Cambridge University Press
,
Cambridge
,
1995
).
24.
T. M.
Cover
and
J. A.
Thomas
,
Elements of information theory
,
Wiley Series in Telecommunications and Signal Processing
, 2nd ed. (
Wiley-Interscience
,
2006
).
25.
A.
Hertz
and
N. J.
Cerf
,
J. Phys. A: Math. Theor.
52
,
173001
(
2019
).
26.
W.
Li
,
J. Stat. Phys.
60
,
823
(
1990
).
27.
W. R.
Inc
,
Mathematica, Version 12.3
,
Champaign, IL
, p.
2021
.
28.
R. M.
Corless
,
G. H.
Gonnet
,
D. E. G.
Hare
,
D. J.
Jeffrey
, and
D. E.
Knuth
,
Adv. Comput. Math.
5
,
329
(
1996
).
29.
S.
Shin
and
H.
Metiu
,
J. Phys. Chem.
100
,
7867
(
1996
).
30.
S.
Shin
and
H.
Metiu
,
J. Chem. Phys.
102
,
9285
(
1995
).
31.
M. D.
Feit
,
J. A.
Fleck
, and
A.
Steiger
,
J. Comput. Phys.
47
,
412
(
1982
).
32.
R.
Kosloff
and
H.
Tal-Ezer
,
Chem. Phys. Lett.
127
,
223
(
1986
).
33.
W. P.
Schleich
,
Quantum Optics in Phase Space
(
Wiley VCH
,
Berlin
,
2001
).
34.
H. G.
Laguna
and
R. P.
Sagar
,
Entropy
15
,
1516
(
2013
).
35.
J. F.
Santos
,
C. H.
Vieira
, and
P. R.
Dieguez
,
Physica A
579
,
125937
(
2021
).
36.
P.
Schürger
and
V.
Engel
,
J. Phys. Chem. Lett.
14
,
334
(
2023
).
37.
P.
Schürger
and
V.
Engel
,
Phys. Chem. Chem. Phys.
25
,
28373
(
2023
).
38.
P.
Schürger
and
V.
Engel
,
Entropy
25
,
970
(
2023
).
39.
D.-J.
Jwo
,
T.-S.
Cho
, and
A.
Biswal
,
Entropy
25
,
1177
(
2023
).
40.
J. H.
Eberly
,
N. B.
Narozhny
, and
J. J.
Sanchez-Mondragon
,
Phys. Rev. Lett.
44
,
1323
(
1980
).
41.
I. S.
Averbukh
and
N. F.
Perelman
,
Phys. Lett. A
139
,
449
(
1989
).
42.
R. W.
Robinett
,
Phys. Rep.
392
,
1
(
2004
).
43.
O.
Knospe
and
R.
Schmidt
,
Phys. Rev. A
54
,
1154
(
1996
).
44.
C.
Leichtle
,
I. S.
Averbukh
, and
W. P.
Schleich
,
Phys. Rev. A
54
,
5299
(
1996
).