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.

## I. INTRODUCTION

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.

## II. THEORY

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.

^{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

^{24,25}determined from the total density,

*ρ*

^{nuc}(

*x*) and

*ρ*

^{el}(

*y*) with the information extracted from the total density

*ρ*(

*x*,

*y*).

^{23}

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} + *λ*_{1}⟨*x*⟩, 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}

## III. RESULTS

### A. Analytical considerations

*a*and

*b*are real valued parameters. From the total density, the single-particle densities are determined using

*MATHEMATICA*,

^{27}and the results are

The total density is a Gaussian multiplied by a cos^{2} 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 cos^{2} 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 $\rho a,bnuc(x)$ (blue lines) and electronic densities $\rho 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.

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.

*I*

_{a,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

*I*

_{a,0}= 0, which is independent of the value of

*a*. All curves start at zero and grow monotonically. They approach an upper bound at

*I*

_{a,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

*I*

_{a,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

*b*.

*MATHEMATICA*) with the result

*b*= 1. By inserting the obtained results into the definition of the linear correlation coefficient, one finds that

*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

*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 *corr*_{a,b} are shown in the lower panel of Fig. 2 for different values of *b*. As predicted, extrema are found at *a*_{min} ≈ 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 cos^{2} 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 *corr*_{a,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 *corr*_{a,b}.

### B. A numerical example: Electron-nuclear dynamics in the Born–Oppenheimer case

*r*) and a proton (coordinate

*R*) moving in one dimension.

^{29,30}They interact with each other and two fixed nuclei at positions of

*R*

_{1,2}= ±5 Å through screened Coulomb interactions. The interaction potential is (in atomic units)

*erf*denotes the error function. For the screening parameters, we use

*R*

_{f}= 1.5 Å and

*R*

_{c}= 1.0 Å, and Δ is the energy shift adjusted such that the potential minimum occurs at zero energy in the region of our spatial grid.

*M*, and

*m*

_{e}is the electron mass.

^{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

*N*

_{0}, and the Gaussian is parameterized with

*R*

_{0}= −3.5 Å and

*β*

_{0}= 7.14 Å

^{−2}. The function

*φ*

_{0}(

*r*;

*R*) 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

*V*

_{0}(

*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 values^{33} 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}

*V*

_{0}(

*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}(

*r*;

*R*) 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

*N*

_{t}, the Gaussian width parameter

*β*

_{t}, and the center of the Gaussian

*R*

_{t}, where the subscripts indicate that these parameters depend on time. However, significant deviations were found in the momentum space MI

^{38}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.

^{37}that the coordinate space MI takes the form

*γ*= 0.733 Å

^{−2}is used, which is determined from the electronic wave function

*φ*

_{0}(

*r*;

*R*). It is seen that $I\xi appr(t)$ reproduces the numerically obtained MI accurately, emphasizing the quality of the Gaussian approximation.

*corr*

_{ξ}(

*t*) is shown in Fig. 3, second panel from below. It increases as a function of time and converges to a value of $\u22480.9$ after about 100 fs. The curve is compared to the analytical function derived from the BO ansatz, which reads

^{37}

^{24,39}which is

*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.

^{38}

*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.

^{38}

^{,}

^{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).

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.

*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

*c*

_{m}and harmonic eigenenergies

*E*

_{m}. In order to get a hand on the number of nodes in the nuclear density, we assume that at a time

*t*

_{s}, the density is dominated by a single term in Eq. (30) (with quantum number

*m*=

*n*) so that it can be written as

*H*

_{n}. 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).

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.

*t*=

*t*

_{s}and estimate the number of nodes

*n*appearing in the spatial density

*ρ*(

*r*,

*R*,

*t*

_{s}). 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

Inserting the numerically derived variances into the latter equations, i.e., $\sigma R2(n)=\sigma R2(ts)$ and $\sigma K2(n)=\sigma 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\xi HO,I\pi HO,corr\xi HO$ and $corr\pi HO$, respectively. It is seen that in position-space, the approximate MI $(I\xi HO)$ and correlation $(corr\xi HO)$ agree well with the numerical curves. The same applies to the momentum-space correlation $(corr\pi 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 *t*_{n} = 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\xi HO=0.95$, *I*_{π}(281) = 0.31, and $I\pi 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.

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 *a*_{min} [Eq. (13)].

## IV. SUMMARY

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: LIMIT OF A LARGE NUMBER OF NODES

*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

*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

*x*,

*y*) → (

*x*−

*y*,

*y*) so that the density takes the form

*y*, which is performed with MATHEMATICA and yields

*x*are to be solved, i.e.,

*T*

_{n}(

*x*,

*a*) as

*T*

_{1}(

*x*,

*a*),

*T*

_{2}(

*x*,

*a*),

*T*

_{3}(

*x*,

*a*),

*T*

_{4}(

*x*,

*a*), and

*T*

_{5}(

*x*,

*a*), and one finds

*T*

_{6}(

*x*,

*a*) and

*T*

_{7}(

*x*,

*a*) are obtained by first using

*T*

_{7}(

*x*,

*a*) vanishes and does not contribute to the MI in the considered limit. The term of Eq. (A21) can be integrated, which yields

*T*

_{8}(

*x*,

*a*), MATHEMATICA does not provide a solution, so we need to analyze it in more detail. It can be written as

*T*

_{8}(

*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

*T*

_{8}(

*x*,

*a*) yields

*a*→

*∞*, and the asymptotic value of the MI is calculated as

## REFERENCES

*Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy*

*Photodissociation and Photoionization*

*Photodissociation Dynamics*

*Optical Coherence and Quantum Optics*

*Elements of information theory*

*Wiley Series in Telecommunications and Signal Processing*

*Quantum Optics in Phase Space*