Introduction. The term localization carries many different meanings across physics—from Anderson localization in disordered systems to Wannier localization in periodic solids and spatial confinement in quantum dots. Here we focus on a distinct notion: the propensity of electrons to stay together or apart, as encoded in two‑particle correlations within regions associated with bonds, atomic shells, and lone pairs.

Real‑space indicators of electron localization can be viewed as modern counterparts to Lewis structures [1] to understand the organization of electrons in atoms, molecules, and materials. They are also pervasive in both density‑functional and wave‑function electronic‑structure methods. Among them, the electron localization function (ELF) is arguably the most widely used. Initially introduced by Becke and Edgecombe through the curvature of the like‑spin conditional pair density at the Hartree–Fock level [2], and later reinterpreted as a Pauli (Kohn–Sham) kinetic‑energy excess relative to its bosonic von Weizsäcker limit [3], [4], the ELF has proved invaluable for visualizing bonds and shells.

Yet, despite its success [5]–[8], its construction includes empirical ingredients, and its semi‑local form [9] can limit its ability to capture intrinsically nonlocal electronic features [10]. This empirical content is often accepted on the grounds that ELF analyses are largely qualitative. However, ELF is also widely used to carefully analyze and compare bonds across different systems and environments. Elements of the ELF also enter modern density‑functional approximations [11]–[13], further motivating the need for a more transparent theoretical foundation to support quantitative objectives. Interpretative challenges for the ELF have been reported as well [14]. Significant efforts have been made over the years to place the ELF on safer theoretical grounds [15]–[17], and Kohout’s ELI family of descriptors [18] avoids ad-hoc scaling and normalization by adopting a space-partitioning approach.

Quantum‑information approaches have meanwhile shown that electronic entanglement is a powerful descriptor of bonding, reactivity, and correlation in molecules [19]–[27]. Yet most of these measures are orbital‑resolved and can vary under orbital rotations [28], [29], limiting their suitability as universal real-space indicators.

Along parallel developments, spin entanglement has been investigated in the non‑interacting electron gas [30]–[33], suggesting that a real‑space, basis‑independent localization indicator may be derived from spin entanglement [32]. This viewpoint offered an alternative rationalization of the internal workings of the ELF. However, the same analysis did not overcome the empirical elements in the ELF construction.

In the following, we address the longstanding challenge of formulating electron localization through a well-defined spin-entanglement indicator. We achieve this by employing a measure of entanglement derived from the underlying two‑spin correlations. For each pair of spatial points, we construct a properly normalized mixed two‑spin state and quantify its entanglement via the concurrence function [34]. As we show below, this procedure provides a physically grounded measure of localization in correlated many‑electron states and can be applied straightforwardly, yielding consistent and reliable results.

Formulation.  Given an \(N\)‑electron state \(\Psi\), the corresponding one‑body reduced density matrix (1RDM) is \[\begin{align} \gamma_1({\mathbf{x}}_1;{\mathbf{x}}_1') = &N \!\! \int \! d{\mathbf{x}}_{2}\dots \! \int \! d{\mathbf{x}}_{N} \nonumber \\ &\Psi({\mathbf{x}}_1,{\mathbf{x}}_2 \dots {\mathbf{x}}_N) \Psi^*({\mathbf{x}}_1',{\mathbf{x}}_2 \dots {\mathbf{x}}_N)\, \label{eqn:1RDM} \end{align}\tag{1}\] with natural orbitals \(\psi_i\) and occupations \(0\le n_i\le 1\) defined by \[\begin{align} \int d{\mathbf{x}}_{1}'\,\gamma_1({\mathbf{x}}_1;{\mathbf{x}}_1')\,\psi_i({\mathbf{x}}_1') = n_i\,\psi_i({\mathbf{x}}_1). \label{nt} \end{align}\tag{2}\] From \(\gamma_1\), we may form the anti-symmetrized two-body correlation function \[\begin{align} \tilde{\gamma}_{2}({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathbf{x}}_1',{\mathbf{x}}_2') = \left[\gamma_{1} \wedge \gamma_{1}\right]({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathbf{x}}_1',{\mathbf{x}}_2') \label{eqn:dm2dm1} \end{align}\tag{3}\] where “\(\wedge\)” stands for the wedge product ¹.

The correlation function \(\tilde{\gamma}_{2}\) is significant for three main reasons: (i) it naturally generalizes the conventional Hartree and Fock terms to the correlated \(\gamma_1\) [35]; (ii) as the seminal works by Bader [36], [37], Becke [2], Savin [3], and collaborators have demonstrated, \(\tilde{\gamma}_{2}\) encodes key elements of the physics that underlies electron localization; and (iii) the auxiliary object \(\tilde{\gamma}_{2}\) can encode singlet–triplet mixing that is not visible when the concurrence is extracted directly from the pure-state 2RDM \(\gamma_2\) ².

Here we exploit the singlet–triplet mixing in \(\tilde{\gamma}_2\) to introduce a measure of electron localization based on spin entanglement. The underlying physical picture is straightforward: when two electrons share a sufficiently small region of space, the antisymmetry of the total state forces them into a singlet—and thus maximally entangled—spin configuration. As the two electrons move apart, the weights of the triplet components tend to increase, and the two-spin state becomes less entangled. Thus, \(\tilde{\gamma}_2\) provides a connection between electron localization and spin entanglement that can be used to analyze both.

We stress that \(\tilde{\gamma}_{2}\) is not used here as an approximation to \(\gamma_2\), but as an auxiliary quantity: \(\gamma_2\) and \(\tilde{\gamma}_{2}\) encode different information, with \(\gamma_2\) determining the two-particle probability of the state under consideration and \(\tilde{\gamma}_{2}\) carrying the local singlet–triplet structure relevant for the localization analysis developed here. We also note that the triplet (antisymmetric) component of the 2RDM has been used previously in localization-related analyses [38]. The focus is thus on the state described by \(\tilde{\gamma}_{2}\), while the full two-body cumulant [39]—i.e., \(\left(\gamma_2 - \tilde{\gamma}_{2}\right)\)—is neither required nor employed.

The properties of \(\tilde{\gamma}_2\) deserve further comments. When \(\gamma_1\) is idempotent, as in Restricted Hartree-Fock (RHF)[40], \(\tilde{\gamma}_2\) coincides with \(\gamma_2\), while when \(\gamma_1\) carries fractional occupation numbers—as for a multi-configurational state—\(\tilde{\gamma}_2\) need not satisfy pure-state \(N\)-representability conditions. In order to stress this point and avoid any confusion, in this paper \(\gamma_2\) denotes the pure-state 2RDM defined above.

To estimate the entanglement of the spin pairs, we consider the spatially-dependent \(4\times4\) spin matrix \(\tilde{\boldsymbol{\gamma}}_2\), whose “hyper-diagonal" \((\boldsymbol{r}_1=\boldsymbol{r}_1',\boldsymbol{r}_2=\boldsymbol{r}_2')\) provides a bona‑fide two‑qubit state: \[\begin{align} \tilde{\boldsymbol{\Lambda}}(\boldsymbol{r}_1,\boldsymbol{r}_2) =\frac{\tilde{\boldsymbol{\gamma}}_2(\boldsymbol{r}_1,\boldsymbol{r}_2)}{\mathrm{Tr}\,[\tilde{\boldsymbol{\gamma}}_2(\boldsymbol{r}_1,\boldsymbol{r}_2)]}, \label{eqn:D} \end{align}\tag{4}\] where \(\mathrm{Tr}\) denotes a trace in the two‑spin space. This step involves a non‑empirical trace normalization and yields a proper density operator (unit trace, positive semidefinite). In systems with spin‑rotational symmetry, \(\tilde{\boldsymbol{\Lambda}}\) assumes the Werner form [41] \[\tilde{\boldsymbol{\Lambda}}(\boldsymbol{r}_1,\boldsymbol{r}_2) = p(\boldsymbol{r}_1,\boldsymbol{r}_2)\,|\Psi^{-}\rangle\langle\Psi^{-}| +\bigl[1-p(\boldsymbol{r}_1,\boldsymbol{r}_2)\bigr]\;\mathcal{I}/4. \label{werner}\tag{5}\] Here, \(|\Psi^{-}\rangle\) is the singlet state, \(\mathcal{I}\) is the two‑spin identity, \(p(\boldsymbol{r}_1,\boldsymbol{r}_2) = |\gamma_1(\boldsymbol{r}_1,\boldsymbol{r}_2)|^2 / \left[2 n (\boldsymbol{r}_1)n(\boldsymbol{r}_2)- |\gamma_1(\boldsymbol{r}_1,\boldsymbol{r}_2)|^2 \right]\) is the excess singlet probability determined by the spin coherence in \(\tilde{\boldsymbol{\gamma}}_2\), and \(n(\boldsymbol{r})=\gamma_1(\boldsymbol{r},\boldsymbol{r})\) is the electron density.

We quantify the spin entanglement encoded in \(\tilde{\boldsymbol{\Lambda}}\) through the concurrence, \[\begin{align} \widetilde{{\mathcal{C}}}(\boldsymbol{r}_1,\boldsymbol{r}_2) =\max\!\left\{0,\;2\tilde{\lambda}_{\max}(\boldsymbol{r}_1,\boldsymbol{r}_2)-\sum_{i=1}^4\tilde{\lambda}_i(\boldsymbol{r}_1,\boldsymbol{r}_2)\right\}, \label{eqn:C1} \end{align}\tag{6}\] where \(\{\tilde{\lambda}_i\}\) are the eigenvalues of the Wootters matrix \(R=\sqrt{\sqrt{\tilde{\boldsymbol{\Lambda}}}\,{\tilde{\boldsymbol{\Lambda}}_{\rm SF}}\,\sqrt{\tilde{\boldsymbol{\Lambda}}}}\), with spin-flipped (SF) density matrix \({\tilde{\boldsymbol{\Lambda}}_{\rm SF}}=(\sigma_y\!\otimes\!\sigma_y)\tilde{\boldsymbol{\Lambda}}^*(\sigma_y\!\otimes\!\sigma_y)\) [34]. For the states considered here \({\tilde{\boldsymbol{\Lambda}}_{\rm SF}}=\tilde{\boldsymbol{\Lambda}}\), so that \(R=\tilde{\boldsymbol{\Lambda}}\) and \(\{\tilde{\lambda}_i\}\) reduce to the eigenvalues of \(\tilde{\boldsymbol{\Lambda}}\) itself. For a Werner state this gives \[\widetilde{{\mathcal{C}}}(\boldsymbol{r}_1,\boldsymbol{r}_2)=\max\!\left\{\frac{3p(\boldsymbol{r}_1,\boldsymbol{r}_2)-1}{2},\,0\right\}, \label{conc}\tag{7}\] so that \(\widetilde{{\mathcal{C}}}=0\) for \(p\le 1/3\) and \(\widetilde{{\mathcal{C}}}=1\) when \(\tilde{\boldsymbol{\Lambda}}\) coincides with a pure singlet state (\(p=1\)). The level at which \(\gamma_1\) is determined (e.g., RHF, or beyond) propagates through Eqs. 3 –7 and thus influences the degree of entanglement given by \(\tilde{\mathcal{C}}\). Again, the notation underlines that the concurrence computed for estimating the localization uses \(\tilde{\gamma}_{2}\) [see Eq. 3 ] rather than the pure-state representable \(\gamma_2\). Because \(\tilde{\gamma}_{2}\) is generally not pure-state \(N\)-representable, \(\tilde{\boldsymbol{\Lambda}}\) naturally incorporates singlet–triplet mixing that emerges in strongly correlated regimes.

This construction has three immediate consequences: (i) it associates maximum localization with singlet configurations rather than with fermions endowed with a bosonic kinetic-energy form [3], [4], [42] — an interpretation which, if taken too far, would contradict the fermionic origin of chemical bonding, underscoring that full nonlocality is conceptually essential; (ii) consistent with this, it predicts \(\widetilde{{\mathcal{C}}}= 0\) whenever the exchange term is small compared to the Hartree term in Eq. 3 [see Eq. 7 ]; and (iii) \(\widetilde{{\mathcal{C}}}\) is rigorously bounded in \([0,1]\) with a direct quantum-information meaning, requiring neither heuristic normalization nor ad hoc non-linear remapping (see below).

Increasing internuclear separation in molecules enhances static correlation and leads to a gradual disentanglement of the spin pairs (see below).

Next, it is instructive to see how the conventional ELF is related to \(\tilde{\mathcal{C}}\). To recover the ELF from our formulation, four reductions must be imposed: (i) One first replaces \(\tilde{\boldsymbol{\gamma}}_2(\boldsymbol{r},\boldsymbol{r}+{\boldsymbol{u}})\) with its angular average over \({\boldsymbol{u}}\), thereby discarding directional information; (ii) one then expands the concurrence for small inter-electronic separations \(u\), keeping the leading term, \(\widetilde{{\mathcal{C}}}(\boldsymbol{r},u) \approx \max\!\left\{0,\, 1 - \left[u/l_E(\boldsymbol{r})\right]^2 + \ldots\right\},\) where (assuming real-valued natural orbitals) \(l_E(\boldsymbol{r}) \propto \left\{ \tau(\boldsymbol{r})/n(\boldsymbol{r}) -|\nabla n(\boldsymbol{r})|^2/4n^2(\boldsymbol{r}) \right\}^{-1/2}\) is a characteristic entanglement length, expressed in terms of the particle density \(n(\boldsymbol{r}) = \sum_i n_i |\psi_i(\boldsymbol{r})|^2\) and semi-local quantities: the kinetic-energy density \(\tau(\boldsymbol{r}) = \nabla \cdot \nabla' \gamma_1(\boldsymbol{r},\boldsymbol{r}')|_{\boldsymbol{r}'=\boldsymbol{r}} = \sum_i n_i |\nabla \psi_i(\boldsymbol{r})|^2\) and the density gradient \(\nabla n(\boldsymbol{r})\); (iii) one may normalize \(l_E(\boldsymbol{r})\) to its homogeneous‑electron‑gas value \(l_E^{\rm unif} (\boldsymbol{r})\propto n(\boldsymbol{r})^{-1/3}\); and (iv) finally, one may apply a nonlinear logistic‑like mapping that is not dictated by any known physical principle, \({\rm ELF}(\boldsymbol{r}) = 1 / \,\{1 + [\,l_E^{\rm unif}(\boldsymbol{r})/l_E(\boldsymbol{r})\,]^2\}\,\), which reproduces the Becke–Edgecombe ELF when \(\tilde{\gamma}_{2}\) is evaluated at the RHF level. By contrast, the concurrence itself retains full two-point information and requires no heuristic normalization or empirical logistic‑like mapping.

The above analysis provides a unified perspective to view both the traditional RHF‑based ELF and its correlated extensions [18], [42], [43] as different instances of the same reduction. Furthermore, it highlights the empirical elements in the previous ELF constructions.

We emphasize that the reduction from \(\widetilde{{\mathcal{C}}}(\boldsymbol{r}_1,\boldsymbol{r}_2)\) to the \({\rm ELF}(\boldsymbol{r})\) yields \(l_E(\boldsymbol{r})\), a natural three-dimensional indicator that assigns a quantum-information meaning to a key ingredient of the ELF—the Pauli kinetic energy excess (here divided by the particle density). This observation was already made at the HF/Kohn–Sham level in Ref. [32]; the present work extends it to correlated 1RDMs and identifies the full six-dimensional \(\widetilde{{\mathcal{C}}}(\boldsymbol{r}_1,\boldsymbol{r}_2)\) as the underlying nonlocal parent quantity encoding electron localization from the outset. Conversely, this perspective also suggests that \(\widetilde{{\mathcal{C}}}(\boldsymbol{r}_1,\boldsymbol{r}_2)\) can serve as a starting point for constructing other, quantum-information based, lower-dimensional projections.

Applications.  We examine the dissociation of closed‑shell diatomic molecules into open‑shell atomic fragments, including both covalent and ionic bonds as well as a representative singlet excitation. We perform RHF and complete active space self-consistent field (CASSCF) [44] calculations ³. In the following, CASSCF(\(n\),\(m\)) indicates a CASSCF calculation with \(n\) active electrons distributed within \(m\) active orbitals.

All figures show a projection \(\widetilde{{\mathcal{C}}}(z_1,z_2)\) onto the internuclear \(z\)‑axis, with the two atoms appearing in opposite corners of each panel, except in Fig. 4 where a wider window is shown. Each value \(\widetilde{{\mathcal{C}}}(z_1,z_2)\) gives the entanglement of the two‑spin state extracted from \(\tilde{\gamma}_{2}\) [see Eq. 4 ]. Regions with \(\widetilde{{\mathcal{C}}}>0\) contain entangled pairs, with larger values indicating stronger singlet character. Along the diagonal (\(z_1=z_2\)), antisymmetry enforces \(\widetilde{{\mathcal{C}}}=1\), while along the anti‑diagonal (\(z_1=-z_2\)) the electrons occupy opposite sides of the bond midpoint.

The ground state of H\(_2\) is a pure singlet with only two electrons. Hence, the concurrence computed directly from the corresponding pure‑state 2RDM equals one everywhere, independently of whether the underlying wavefunction is RHF, CASSCF(2,2), or exact. Thus, a visualization based solely on \(\gamma_2\) cannot reveal whether an approximate wavefunction correctly describes dissociation.

Figure 1 shows that this limitation is resolved by resorting to the auxiliary \(\tilde{\gamma}_{2}\). At the RHF level, \(\widetilde{{\mathcal{C}}}\) remains equal to one, consistent with the fact that RHF does not describe bond breaking. At the CASSCF level, however, triplet components mix in as the bond is stretched, and \(\widetilde{{\mathcal{C}}}\) decreases accordingly. Near equilibrium (left), the entanglement is only mildly reduced upon moving the electrons toward opposite atoms; beyond the Coulson–Fischer point [46] (right), the collapse of \(\widetilde{{\mathcal{C}}}\) clearly renders the idea of a stretched bond and the emergence of two separate fragments.

For multi‑electron systems, F\(_2\) in Fig. 2 and N\(_2\) in Fig. 3, the concurrence maps encode both bonding and shells. Bonding information resides in the regions where \(z_1\) and \(z_2\) lie on opposite sides of the diagonal: a central bright island connecting the two atoms indicates inter‑atomic singlet entanglement, while its collapse into a thin bridge (a size to be intended relative to the internuclear distance) signals dissociation. Thus, the effect of static correlations as compared to a mean-field solution is made apparent. We notice that the narrow dip at the equilibrium distance, previously reported in ELF studies of F\(_2\) [47], is also visible in the fully non‑local concurrence. Bright rectangular blocks close to each nucleus—where both electrons reside on the same atom—reflect the atomic shell structure: within a given shell, \(\widetilde{{\mathcal{C}}}\) remains high, whereas dark channels mark the shell boundaries (see Figs. 2 and 3).

Energetically, RHF is an acceptable approximation near equilibrium but fails in stretched geometries, where static correlation becomes essential. Consistently, at the RHF level the concurrence spuriously increases away from the bond center in stretched regimes, artificially suggesting a reinforcement of the bond (upper right panel of Figs. 2 and 3). At the CASSCF level, by contrast, the collapse of \(\widetilde{{\mathcal{C}}}\) between the bright atomic blocks (bottom right panel) is the real‑space signature of bond breaking and the emergence of two independent fragments.

Figure 4 provides a zoomed‑out visualization for N\(_2\) in a slightly stretched geometry, making lone pairs visible through \(\widetilde{{\mathcal{C}}}\). It is also clear that, at this distance, correlation effects are minor, and yet not entirely negligible along the stretched bond.

We next examine LiF, which passing through an avoided crossing undergoes an ionic–to–neutral charge transfer upon dissociation—a regime beyond the capabilities of RHF. Equal‑weight SA‑CASSCF(2,2) calculations for the \(S_0\) and \(S_1\) singlets correctly resolve the avoided crossing near \(R\approx 4.3\) Å. Figure 5 shows \(\widetilde{{\mathcal{C}}}\) at internuclear distances about \(0.5\) Å before and \(0.5\) Å after this crossing. Before the crossing (top row), \(S_0\) is ionic, with \(\widetilde{{\mathcal{C}}}\) displaying two disconnected blocks centered on the atoms, while \(S_1\) is covalent, featuring a bright diagonal island (\(z_1\!\approx\! z_2\)) between the atomic regions. After the crossing (bottom row), \(S_0\) becomes covalent and its \(\widetilde{{\mathcal{C}}}\) develops a clear bridging island, while \(S_1\) becomes ionic and its \(\widetilde{{\mathcal{C}}}\) shows separated blocks that merge into a broader region. The key diagnostic feature is thus the presence (covalent) or absence (ionic) of a distinct diagonal bridging island; after the crossing, the covalent island shrinks and merges into the wider ionic pattern.

Conclusions.  We have introduced a quantum-information-based measure of electron localization derived from the spin entanglement encoded in a normalized two‑spin state built from the \(\gamma_1\)-derived Fock–Dirac two-body object. By extracting a spatially resolved two‑spin state and evaluating its concurrence, the method captures both local and non‑local aspects of electronic structure, dispensing with the short‑range assumptions and empirical elements of traditional localization functions. This establishes a rigorous quantum‑information foundation for next‑generation localization indicators relevant to both density‑functional and wave‑function methodologies.

Notably, the calculation involves the one‑body reduced density matrix while still reflecting non‑local correlations. This feature enables a computationally efficient measure of quantum correlations that can be evaluated with existing software and readily integrated into data-driven workflows [48], [49] for enhanced electronic‑structure modeling.

Acknowledgments.  The authors acknowledge financial support from: the Ministero dell’Università e della Ricerca (MUR) under the Project PRIN 2022 number 2022W9W423 (SP, FT, CA) and the PNRR Project PE0000023-NQSTI (FT); an Australian Research Council (ARC) Discovery Project DP200100033 (TG) and an ARC Future Fellowship FT210100663 (TG).

References↩︎