Physical entanglement between localized orbitals

The notion of locality is the starting point for the theory of entanglement. We consider two distant parties, Alice and Bob, each possessing a local quantum system with a finite-dimensional Hilbert space denoted by $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively. The two systems are local in the sense that an action of one party on their subsystem has no effect on the other. Mathematically speaking, this separation is described by the tensor product structure of both the total Hilbert space $\mathcal{H} = \mathcal{H}_{A} \otimes \mathcal{H}_{B}$ and the total algebra of observables $\mathcal{A} = \mathcal{A}_{A} \otimes \mathcal{A}_{B}$.

With the stage set, a quantum state ρ_AB shared by Alice and Bob is said to be separable, if and only if it can be prepared by local operations and classical communication [41]. Local operations allow Alice and Bob to prepare any product state $\rho_{A}\otimes \rho_{B}$, and classical communication can be used to generate convex mixtures of those states. This then gives rise to the convex set of bipartite separable states $D_\mathrm{sep}$: (bipartite separability [41]).

Definition 2.1 A bipartite state ρ_AB on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ is separable if and only if it can be expressed as

$$\begin{equation} \rho_{AB} = \sum_i \,p_i\, \rho_{A}^{\left(i\right)} \otimes \rho_{B}^{\left(i\right)} \end{equation} \tag{ 1 }$$

with $p_i \unicode{x2A7E} 0$ and $\sum_i p_i = 1$. The set of separable states is denoted as $D_\mathrm{sep}$. Any state that is not in $D_\mathrm{sep}$ is necessarily entangled.

A non-separable, or entangled state contains valuable entanglement resource for QIP tasks [23–28]. This resourcefulness of entanglement is distinct from the mutual information, which contains both quantum and classical correlation [42]. Therefore the detection and quantification of entanglement in a bipartite system is of crucial interest. For a pure state $\rho_{AB} = |\Psi_{AB}\rangle\langle\Psi_{AB}|$, it is easy to verify whether or not it belongs to $D_\mathrm{sep}$. It suffices to check if it is a pure product state. Or equivalently one can check the purity of the reduced state $\rho_{A/B} = \mathrm{Tr}_{B/A}[|\Psi_{AB}\rangle\langle\Psi_{AB}|]$. The degree of mixedness of $\rho_{A/B}$ thus accounts for the entanglement in ρ_AB . This leads to a well-defined entanglement measure for pure states

$$\begin{equation} E\left(|\Psi_{AB}\rangle\langle\Psi_{AB}|\right) = S\left(\rho_{A/B}\right), \end{equation} \tag{ 2 }$$

where $S(\rho) = -\mathrm{Tr}[\rho\log(\rho)]$ is the von Neumann entropy. For mixed states, however, there is no such straightforward analytic check for separability, except for qubit–qubit and qubit–qutrit systems [43, 44]. In fact, even the problem of determining the separability of a general mixed quantum state is non-deterministic polynomial-tim (NP)-hard [45, 46], and in practice can only be approximately solved with significant computational efforts using symmetric extension and semidefinite programming [47, 48]. On the formal level, the degree of entanglement of a mixed state can be defined geometrically as the distance to the closest separable state (here and in the following we occasionally skip the index AB of the joint state)

$$\begin{equation} E\left(\rho\right) = \min_{\sigma \in D_\mathrm{sep}} S\left(\rho\|\sigma\right), \end{equation} \tag{ 3 }$$

where the generalized distance function is chosen to be the relative entropy

$$\begin{equation} S\left(\rho\|\sigma\right) = \mathrm{Tr}\left[\rho\left(\log\left(\rho\right)-\log\left(\sigma\right)\right)\right] \end{equation} \tag{ 4 }$$

for its rich quantum informational meaning. To be more precise, the quantum relative entropy $S(\rho\|\sigma)$ exponentially suppresses the 'difficulty' of telling ρ and σ apart [49], in the spirit of the so-called quantum Sanov's theorem [49]. The measure (3) is called the relative entropy of entanglement [50], and will be used as the standard entanglement measure throughout this paper. When applied to pure states, (3) reduces to the pure state measure via the von Neumann entropy (2) [49].

In the previous section, we explained that the notion of entanglement necessitates two ingredients of locality: (1) a tensor product structure of the Hilbert space, and (2) a tensor product structure of the algebra of observables. In the context of fermionic systems, this raises the following questions. Does there also exist a tensor product structure in a fully-antisymmetrized N-fermion Hilbert space, e.g. for a hydrogen molecule with two electrons? And how can one identify two commuting algebras of observables when fermionic creation/annihilation operators anticommute?

The fermionic analog of the bipartition into distant laboratories corresponds to a choice of splitting the one particle Hilbert space $\mathcal{H}_1$ (or a subspace of it) into two orthogonal subspaces. For instance, $\mathcal{H}_1$ can be realized by an array of Wannier orbitals, or an arrangement of molecular spin orbitals. From here on we shall refer to the elements of an orthonormal basis of $\mathcal{H}_1$ as orbitals. A partition of an orbital basis $\{|\varphi_i\rangle\}_{i = 1}^{d}$ into two complementing subsets A and B dissolves the one particle Hilbert space into a direct sum $\mathcal{H}_1 = \mathcal{H}_1^{(A)} \oplus \mathcal{H}_1^{(B)}$. This for example can correspond to separating a lattice into the left and right parts at some cut. Each part generates a local Fock space and together they form a tensor product of the total Fock space $\mathcal{F}(\mathcal{H}_1) \cong \mathcal{F}(\mathcal{H}_1^{(A)}) \otimes \mathcal{F}(\mathcal{H}_1^{(B)})$ through the following basis identification

$$\begin{equation} |n_1,n_2,\ldots,n_d\rangle = |n_1, \ldots, n_{d_{A}} \rangle \otimes |n_{d_{A}+1},\ldots ,n_d\rangle, \end{equation} \tag{ 5 }$$

where $n_1,\ldots,n_d \in \{0,1\}$,

$$\begin{equation} \begin{split} |n_1,\ldots,n_d\rangle \equiv \left(f^{\,\dagger}_{1}\right)^{n_1}\cdots\left(f^{\,\dagger}_{d}\right)^{n_d}|0\rangle, \end{split} \end{equation} \tag{ 6 }$$

and $f^{\,(\dagger)}_{i}$ denotes the operator that annihilates (creates) a fermion in state $|\varphi_i\rangle$.

As for the second necessary ingredient, i.e. commuting local algebras of observables, it seems at first sight to be in direct contradiction with the nature of fermionic operators. Suppose Alice and Bob each possess a fermionic mode at the opposite ends of the Universe. Due to the anticommutation relation between the creation operators, Alice's action could in principle still influence Bob's mode, and accordingly realize superluminal signaling [51, 52]. This manifest violation of special relativity is in reality not possible, due to nature's restriction on the local operations by Alice and Bob, which can be formulated as so-called SSRs [29, 30]. The relevant SSRs depend on the conserved observables in the theory or even practical limitations. Nonetheless, the parity SSR [29] (P-SSR) is ubiquitous: one cannot superpose quantum states with different local fermionic parity. Equivalently, physical local observables must always commute with the local parity operator [32, 33, 53, 54] $\hat{\Pi} \equiv \exp(i\pi\sum_l \hat{n}_l) = \hat{P}_+ - \hat{P}_-$ where the sum is over all local orbital degrees of freedom, $\hat{P}_\pm$ is the orthogonal projection onto the even(+)/odd(−) subspace, and $\hat{n}_l \equiv f_l^\dagger f_l$ denotes the particle number operator for mode $|\varphi_l\rangle$. The local observables $\hat{O}$ fulfilling this requirement are polynomials of $f^{\,\dagger}_i f^{\,\dagger}_j$, $f^{\phantom{\dagger}}_i f^{\phantom{\dagger}}_j$, and $f^{{\dagger}}_if^{\phantom{\dagger}}_j$, and satisfy

$$\begin{equation} \hat{O} = \mathcal{G}_{\Pi}\left[\hat{O}\right]\equiv \hat{P}_+ \hat{O} \hat{P}_+ + \hat{P}_- \hat{O} \hat{P}_-, \end{equation} \tag{ 7 }$$

where $\mathcal{G}_\Pi$ denotes the projective map onto the P-SSR compatible sub-algebra. Then, it follows directly that any two observables $\hat{A}$ and $\hat{B}$ referring to two distinct sets of fermionic modes and satisfying locally (7) necessarily commute with each other. In other words, algebras of parity preserving observables on two distinct fermionic subsystems do exhibit a tensor product structure. More generally, any physical quantum channel $\mathcal{E}$ (completely positive and trace-preserving map [55]) on the reduced state must satisfy

$$\begin{equation} \mathcal{E} = \mathcal{G}_\Pi \circ \mathcal{E} \circ \mathcal{G}_\Pi. \end{equation} \tag{ 8 }$$

Nature's restriction on the implementable local observables has a profound implication on the entanglement in fermionic states. From the perspective of state tomography, coherent terms between different local parity states can never be observed under P-SSR, for the observable needed for such measurement is simply not available by the laws of physics. Therefore, the entanglement arising from SSR-violating coherence, also known as the 'fluffy bunny entanglement' [32], is not accessible. This reduction of entanglement due to SSR is proven rigorously in [33] on an abstract mathematical level. Here we demonstrate it in concrete terms with a swapping protocol.

In figure 1 we depict schematically an entanglement swapping protocol between two fermionic orbitals $(\varphi_{A},\varphi_{B})$ in the state ρ_AB , and two qubits $(q_{A},q_{B})$ in the state $|00\rangle$, both shared by Alice and Bob. We shall denote the orbital Fock spaces as $\mathcal{F}_{\varphi_{A/B}}$ and the qubit Hilbert spaces as $\mathcal{H}_{q_{A/B}}$. For both the Fock spaces and the qubit Hilbert spaces, we denote the local basis as $\mathcal{B} = \{|0\rangle,|1\rangle\}$. One should keep in mind that for the qubit space $\mathcal{B}$ is the computational basis, while for the Fock space, $\mathcal{B}$ refers to the occupational basis. The protocol in figure 1 extracts the entanglement between the fermionic modes onto the qubit registers, on which this entanglement could be used to perform more flexible quantum information protocols.

Zoom In Zoom Out Reset image size

The protocol consists of Alice swapping the two states on ϕ_A and q_A , and Bob doing the same on ϕ_B and q_B . Specifically, the SWAP channel $\mathcal{S}$ on a bipartite system is defined element-wise as

$$\begin{equation} \mathcal{S}\!: \,\,|i\rangle \langle \, j | \otimes |k\rangle\langle l | \mapsto |k\rangle \langle l | \otimes |i\rangle\langle \, j |, \end{equation} \tag{ 9 }$$

for all $i,j,k,l\in\{0,1\}$. In a world without SSR, after the SWAP operation, the final state of the two qubits would be exactly ρ_AB and it would accordingly contain exactly the same amount of entanglement as the original two-orbital state. However, the SWAP operation clearly does not preserve the local parities, and hence cannot be realized with physical fermionic operators. Instead, Alice/Bob can only perform the superselected channel

$$\begin{equation} \tilde{\mathcal{S}}_{A/B} = \mathcal{G}_{\Pi_{A/B}} \circ \mathcal{S}_{A/B} \circ \mathcal{G}_{\Pi_{A/B}}. \end{equation} \tag{ 10 }$$

To see the difference between $\mathcal{S}_{A/B}$ and $\tilde{\mathcal{S}}_{A/B}$, let us first observe the action of the SWAP channel on a pure state $|\psi_{A}\rangle = |+\rangle\otimes|+\rangle \in \mathcal{F}_{\varphi_{A}} \otimes \mathcal{H}_{q_{A}}$ where $|+\rangle = \frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$, which clearly violates the P-SSR. In a world free of SSR, $|\psi_{A}\rangle$ would be invariant under the SWAP channel. Under the P-SSR, however, the action of the physically implementable $\tilde{\mathcal{S}}_{A}$ consists of three steps: (1) $\mathcal{G}_{\Pi_{A}}$ turns $|\psi_{A}\rangle$ into a mixed state $\mathcal{G}_{\Pi_{A}} [\rho_{A}] = \frac{\unicode{x1D7D9}}{2} \otimes |+\rangle\langle+|$. Namely, the coherent terms between different fermionic parities are erased. (2) $\mathcal{S}_{A}$ swaps the two tensor states to $\mathcal{S}_{A} \circ \mathcal{G}_{\Pi_{A}} [\rho_{A}] = |+\rangle\langle+| \otimes \frac{\unicode{x1D7D9}}{2}$, which again contains P-SSR violating terms. (3) $\mathcal{G}_{\Pi_{A}}$ then eliminates said terms and the final state is $\mathcal{G}_{\Pi_{A}} \circ \mathcal{S}_{A} \circ \mathcal{G}_{\Pi_{A}} [\rho_A] = \frac{\unicode{x1D7D9}}{2}\otimes\frac{\unicode{x1D7D9}}{2}$. In this example we see evidently that SWAP does not leave $|\psi_{A}\rangle$ invariant, despite its apparent symmetry. Instead, P-SSR filters out any coherence between different parity sectors, and only the superselected state $\tilde{\rho}^\mathrm{P}_{A} \equiv \mathcal{G}_{\Pi_{A}}[|\psi_{A}\rangle\langle\psi_{A}|]$ can be transferred onto the qubit register. Moreover, the final state on the fermionic mode ϕ_A is also the superselected variant of the original qubit state.

Now that we understand the action of (10), we can proceed to write down the action of the protocol $\tilde{S} \equiv \tilde{\mathcal{S}}_{A} \otimes \tilde{\mathcal{S}}_{B}$ on the composite Hilbert space $\mathcal{F}_{\varphi_{A}} \otimes \mathcal{H}_{q_{A}} \otimes \mathcal{F}_{\varphi_{B}} \otimes \mathcal{H}_{q_{B}}$: for any states ρ_AB on $\mathcal{F}_{\varphi_{A}} \otimes \mathcal{F}_{\varphi_{B}}$ and σ_AB on $\mathcal{H}_{q_{A}} \otimes \mathcal{H}_{q_{B}}$, we have

$$\begin{equation} \begin{split} \tilde{\mathcal{S}}\left[\rho_{AB} \otimes \sigma_{AB}\right] = \tilde{\sigma}^{\mathrm{P}}_{AB} \otimes \tilde{\rho}^{\mathrm{P}}_{AB} \end{split}, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} \begin{split} \tilde{\rho}^\mathrm{P}_{AB} &\equiv \mathcal{G}_{\Pi_{A}} \otimes \mathcal{G}_{\Pi_{B}}\left[\rho_{AB}\right] \\ & = \sum_{s,s^{\prime} = +,-} \hat{P}^{\left(A\right)}_s \otimes \hat{P}^{\left(B\right)}_{s^{\prime}} \rho_{AB}\hat{P}^{\left(A\right)}_s \otimes \hat{P}^{\left(B\right)}_{s^{\prime}} \end{split} \end{equation} \tag{ 12 }$$

and the same for $\tilde{\sigma}^\mathrm{P}_{AB}$. The final state $\tilde{\rho}^{\mathrm{P}}_{AB}$ on the qubit registers contains the same information as the original two-orbital state ρ_AB , except it is rid of the coherence between the different local parity sectors, which does not contribute to the physical entanglement. The decohering map $\mathcal{G}_{\Pi_{A}} \otimes \mathcal{G}_{\Pi_{B}}$ thus acts as an entanglement filter as depicted in figure 1, and effectively reduces the available entanglement in the fermionic state. Only the entanglement within each local parity sector can be extracted for operational purposes.

In the previous section, we concluded that the physical entanglement in a bipartite fermionic state ρ_AB in the presence of P-SSR, is precisely the entanglement in its superselected variant $\tilde{\rho}^{\mathrm{P}}_{AB} = \mathcal{G}_{\Pi_{A}} \otimes \mathcal{G}_{\Pi_{B}}[\rho_{AB}]$ quantified in the usual manner without P-SSR. In this section, we will briefly summarize the analytic formula for physical entanglement between two spatial orbitals in electronic systems, derived in [40].

Consider two orbitals A and B with internal spin-$\frac{1}{2}$ degrees of freedom. This system is described by a state ρ_AB acting on the Fock space $\mathcal{F}_A \otimes \mathcal{F}_B$, spanned by $\{|\alpha\rangle_A \otimes |\beta\rangle_B\}$, where $\alpha,\beta\in\{0,\uparrow,\downarrow,\uparrow\downarrow\}$. The mixedness of ρ_AB accounts for the physically relevant scenario where A and B are in contact with other orbitals. A key step in quantifying the entanglement $E(\tilde{\rho}^\mathrm{P}_{AB})$ is understanding its restriction to several smaller subspaces of the total Fock space, after imposing several common symmetries exhibited by realistic condensed matter and quantum chemical systems.

The relevant systems in our case, the tight-binding model and the hydrogen ring, enjoy many symmetries, including (1) particle number, (2) magnetization and (3) reflection symmetry between orbitals A and B (e.g. two orthogonalized 1s orbitals in a hydrogen ring). Under these symmetries, there are only four entangled pure states compatible with local P-SSR:

$$\begin{equation} \begin{split} |\Phi_{\pm}\rangle & = \frac{|0\rangle_A \otimes |\!\uparrow\downarrow\rangle_B \pm |\!\uparrow\downarrow\rangle_A \otimes |0\rangle_B}{\sqrt{2}}, \\ |\Psi_{\pm}\rangle & = \frac{|\!\uparrow\rangle_A \otimes |\!\downarrow\rangle_B \pm |\!\downarrow\rangle_A \otimes |\!\uparrow\rangle_B}{\sqrt{2}}. \end{split} \end{equation} \tag{ 13 }$$

Naturally, if $\tilde{\rho}^{\mathrm{P}}_{AB}$ enjoys the symmetries (1), (2) and (3), then $|\Phi_\pm\rangle$ and $|\Psi_\pm\rangle$ are eigenstates of $\tilde{\rho}^{\mathrm{P}}_{AB}$, with eigenvalues $p_{\pm}$ and $q_{\pm}$, respectively. Even so, the exact P-SSR orbital–orbital entanglement, though entirely analytic, has a heavily involved form, which can be found in the appendix of [40]. Instead, we highlight another related entanglement formula from [40] in the case where the particle number SSR (N-SSR) applies. In this setting, no superposition between different local particle number states is possible and therefore even $|\Phi_\pm\rangle$ are forbidden. The N-SSR compatible two-orbital state is obtained with the same reasoning via the superselection map

$$\begin{equation} \begin{split} \tilde{\rho}^\mathrm{N}_{AB} &\equiv \mathcal{G}_{N_A} \otimes \mathcal{G}_{N_B}\left[\rho_{AB}\right] \\ & = \sum_{n_A,n_B = 0}^{2} \hat{P}_{n_A,n_B} \rho_{AB} \hat{P}_{n_A,n_B}, \end{split} \end{equation} \tag{ 14 }$$

where $\hat{P}_{n_A,n_B} = \hat{P}_{n_A} \otimes \hat{P}_{n_B}$ and $\hat{P}_{n_{A/B}}$ is the projection onto the sector of local particle number $n_{A/B}$ on subsystem $A/B$. The N-SSR entanglement of $\tilde{\rho}^\mathrm{N}_{AB}$ can then be derived, after considerable efforts, to be

$$\begin{equation} \begin{split} E\left(\tilde{\rho}^\mathrm{N}_{AB}\right) = \begin{cases} r\log\left(\frac{2r}{r+t}\right) + t\log\left(\frac{2t}{r+t}\right), \quad &r < t, \\ 0, \quad &r \unicode{x2A7E} t, \end{cases} \end{split} \end{equation} \tag{ 15 }$$

where

$$\begin{equation} \begin{split} &t \equiv \max\left\{q_\pm\right\}, \quad r \equiv \mathrm{Tr}\left[\hat{P}_{1,1} \rho_{AB}\right] - t. \end{split} \end{equation} \tag{ 16 }$$

Note that r < t, or explicitly

$$\begin{equation} \mathrm{Tr}\left[\hat{P}_{1,1} \rho_{AB}\right] < 2\max\left\{q_\pm\right\}, \end{equation} \tag{ 17 }$$

is the exact entanglement criterion. Below we will see that the P- and N-SSR orbital–orbital long-distance entanglement on a localized orbital chain are virtually indistinguishable. The simple form of the latter thus allows us to derive analytically asymptotic properties of the physical entanglement between orbitals for both P- and N-SSR.

We remark that (15) is only applicable for states with the desired symmetry. Applying our formula to symmetry-broken states could lead to erroneous results. For example, according to (15), the entanglement of $\frac{1}{\sqrt{2}}(|\!\!\uparrow\uparrow\rangle + |\!\!\downarrow\downarrow\rangle)$ would be zero, while this state is clearly entangled. However, one could apply a local spin flip (which preserves the entanglement), and arrive at the symmetric $|\Psi_+\rangle$ state, for which the entanglement formula (15) is applicable.

In the reminder of our paper we will apply the entanglement measure (15) and the corresponding one for the weaker P-SSR to two model systems in order to investigate their orbital–orbital entanglement as a function of, e.g. the electron density and the distance between both localized orbitals. This will provide new insights which may guide the utilization of the abundant entanglement resource for QIP tasks performed in realistic materials by tuning the particle density and interaction strength.

A general particle number conserving non-interacting fermionic Hamiltonian can be written as a sum of one-body terms

$$\begin{equation} \hat{H}_\text{free} = \sum_{ij} h_{ij} f^{\,\dagger}_i f^{\phantom{\dagger}}_j, \quad h_{ij} = h_{ji}^\ast. \end{equation} \tag{ 18 }$$

The spectrum of H_free is analytically solved by finding a unitary matrix U such that $\mathbf{d} = \mathbf{U} \mathbf{h} \mathbf{U}^\dagger$ is diagonal. We then rewrite $\hat{H}_\text{free}$ as

$$\begin{equation} \hat{H}_\text{free} = \sum_i d_{ii} c^{\dagger}_i c^{\phantom{\dagger}}_i \equiv \sum_i d_{ii} \hat{n}_i, \quad c^\dagger_i = \sum_j U_{ij} f^{\,\dagger}_j. \end{equation} \tag{ 19 }$$

The N-particle ground states of $\hat{H}_\text{free}$ is simply a Slater determinant state of the form

$$\begin{equation} |\Psi_N\rangle = c^\dagger_1c^\dagger_2\cdots c^\dagger_N|0\rangle \end{equation} \tag{ 20 }$$

with ground state energy $\sum_{i = 1}^N d_{ii}$, assuming the d_ii 's are ordered increasingly.

To study the entanglement within the ground state $|\Psi_N\rangle$, a natural object to consider is the reduced density matrix ρ_A of the ground state on a sub-lattice A. Its entropy quantifies the entanglement between the subsystem A and its complement, and serves as a widely used analytic tool [36–38, 61, 62]. Note that the lattice basis where such partition is performed is not necessarily the one in which h is diagonalized. For the ground states of Hamiltonians of the form (18), according to Wick's theorem the expectation value of any observable is a function of the two-point correlator $\langle \Psi_N| f^{\,\dagger}_i f^{\phantom{\dagger}}_j|\Psi_N \rangle \equiv \gamma_{ji}$. The latter defines the one-particle reduced density matrix (1RDM) γ of the ground state. As ρ_A is uniquely defined via the expectations of all physical observables on the subsystem A [63], one can then in principle express ρ_A using γ. While the explicit expression can be cumbersome, spectra of ρ_A and γ can be easily related via the remarkable shortcut put forward by Peschel [39]. The shortcut makes use of two insights: (1) the reduced state ρ_A is a Gaussian state $\rho_A = K \exp(\sum_{ij}C_{ij}f^{\,\dagger}_i f^{\phantom{\dagger}}_j)$, where K is the normalization factor. (2) The defining matrix $\mathbf{C} = (C_{ij})$ and the 1RDM restricted to the sub-algebra on subsystem A, γ_A , can be simultaneously diagonalized. Combining these two ingredients, Peschel discovered that the two matrices are related via $\mathbf{C} = \ln[(\unicode{x1D7D9}-\gamma_A)/\gamma_A]$.

This shortcut greatly simplifies the complexity of computing the spectrum of ρ_A : one simply needs to find the restricted 1RDM and its eigenvalues. As the von Neumann entropy is only a function of the spectrum, $S(\rho_A)$ can also be readily calculated, without explicit construction of C or ρ_A . In general, however, not all correlation quantities can be determined with the spectrum of the reduced states alone. In particular, suppose A contains two orthonormal orbitals. The entanglement between the two orbitals is not invariant under rotations between them and cannot be retrieved using the same method for $S(\rho_A)$. In the next section, we will break the current paradigm by calculating explicitly the orbital–orbital entanglement in the tight-binding model, using the analytic formula (15).

We consider the following tight-binding Hamiltonian for electrons on a periodic chain with L lattice sites

$$\begin{equation} \hat{H}_\mathrm{tb} = -\frac{1}{2}\sum_{\sigma = \uparrow, \downarrow}\sum_{l = 1}^{L} f^{\,\dagger}_{l\sigma} f^{\phantom{\dagger}}_{l+1,\sigma} + h.c. \end{equation} \tag{ 21 }$$

with periodic boundary condition imposed. Here, $f^{(\dagger)}_{l\sigma}$ annihilates(creates) an electron with spin σ at the Wannier orbital/lattice site labeled by l. Each Wannier orbital can host up to two electrons with different spin. The discrete translation invariance of the system allows us to diagonalize the Hamiltonian via a Fourier transform

$$\begin{equation} c^\dagger_{k\sigma} = \frac{1}{\sqrt{L}} \sum_{l = 1}^L \mathrm{e}^{-\frac{2\pi i}{L}kl} f^{\,\dagger}_{l\sigma}. \end{equation} \tag{ 22 }$$

In this basis, the Hamiltonian (21) is diagonal, and the one-electron spectrum is given by the dispersion relation

$$\begin{equation} E_{k} = -\cos\left(\frac{2\pi k}{L}\right), \quad k = 0, \pm1, \ldots,\pm \Bigl\lceil\frac{L-1}{2}\Bigr\rceil, \end{equation} \tag{ 23 }$$

where in case of L even only one of the signs $\pm \bigl\lceil\frac{L-1}{2}\bigr\rceil$ should be chosen. As the system is non-interacting, the N-electron ground state is the configuration where the N lowest energy levels are filled. For simplicity, we shall assume from now on that the number of electrons is either N = 0, or $N = 4 k_\mathrm{max} + 2$ where $k_\mathrm{max} \unicode{x2A7E} 0$ denotes the highest occupied momentum. In the former case, the ground state is simply the vacuum. In the latter, the N-particle ground state $|\Psi_N\rangle$ is uniquely characterized by the Fermi level

$$\begin{equation} \langle \Psi_N | c^\dagger_{k\sigma} c^{\phantom{\dagger}}_{k^{^{\prime}}\sigma^{^{\prime}}} |\Psi_N\rangle = \delta_{kk^{^{\prime}}}\delta_{\sigma\sigma^{^{\prime}}}\Theta\left(k_\mathrm{max}-|k|\right), \end{equation} \tag{ 24 }$$

where Θ is the Heaviside step function. With (22) and (24) at hand we are ready to compute the matrix elements of the 1RDM in the spatial basis as

$$\begin{equation} \begin{split} \langle \Psi_N|f^{\,\dagger}_{l\sigma} f^{\phantom{\dagger}}_{l^{\prime} \sigma^{\prime}} |\Psi_N\rangle = \frac{1}{L}\delta_{\sigma, \sigma^{\prime}} \sum_{k = -k_\mathrm{max}}^{k_\mathrm{max}} \mathrm{e}^{\frac{2\pi i\left(l-l^{\prime}\right)k}{L} }. \end{split} \end{equation} \tag{ 25 }$$

When $l = l^{\prime}$, all summands in (25) become 1, and thus we obtain

$$\begin{equation} \begin{split} \langle \Psi_N|f^{\,\dagger}_{l\sigma} f^{\phantom{\dagger}}_{l \sigma^{\prime}} |\Psi_N\rangle & = \delta_{\sigma\sigma^{\prime}}\frac{2k_\mathrm{max}+1}{L} \\ & = \delta_{\sigma\sigma^{\prime}}\frac{N}{2L} \equiv \delta_{\sigma\sigma^{\prime}}\eta, \end{split} \end{equation} \tag{ 26 }$$

where $\eta = N/(2L)$ is the particle filling fraction. This agrees with the physical intuition that the electrons move freely along the chain and appear at every site with uniform probability. When $l\neq l^{\prime}$, by introducing $\omega \equiv \frac{2\pi(l-l^{\prime})}{L}$, we arrive at

$$\begin{equation} \begin{split} \langle \Psi_N|f^{\,\dagger}_{l\sigma} f^{\phantom{\dagger}}_{l^{\prime} \sigma^{\prime}} |\Psi_N\rangle & = \frac{1}{L} \delta_{\sigma\sigma^{\prime}}\frac{\mathrm{e}^{\mathrm{i}\omega k_\mathrm{max}}-\mathrm{e}^{-\mathrm{i}\omega k_\mathrm{max}}}{1-\mathrm{e}^{\mathrm{i}w}} \\ & = \frac{1}{L}\delta_{\sigma\sigma^{\prime}} \frac{\sin\left[\omega \left(k_\mathrm{max}+\frac{1}{2}\right)\right]}{\sin\left(\frac{\omega}{2}\right)} \\ & = \frac{1}{L}\delta_{\sigma\sigma^{\prime}} \frac{\sin\left(\frac{\omega N}{4} \right)}{\sin\left(\frac{\omega}{2}\right)}. \end{split} \end{equation} \tag{ 27 }$$

In the thermodynamic limit where $N,L\!\rightarrow\!\infty$ while $\eta\! = \!N/(2L)$ is kept fixed, the off-diagonal elements of the 1RDM become

$$\begin{equation} \begin{split} \lim_{L\rightarrow\infty}\langle \Psi_N|f^{\,\dagger}_{l\sigma} f^{\phantom{\dagger}}_{l^{\prime} \sigma} |\Psi_N\rangle = \frac{\sin\left( \pi d \eta\right)}{\pi d} \equiv W\left(d,\eta\right), \end{split} \end{equation} \tag{ 28 }$$

where $d\! = \!|l-l^{\prime}|$ is the orbital separation distance.

Next, we use Wick's theorem to calculate the quantities needed for the entanglement between the orbitals l and lʹ from the elements of the 1RDM derived above. The relevant operators for calculating r and t in the analytic formula (15) are the orthogonal projection $\hat P_{ - - }^{l{l^\prime }}{\mkern 1mu}$ onto the local odd parity subspace, and the orthogonal projections $\hat{P}^{\,ll^{\prime}}_{\Psi_{\pm}}$ onto $|\Psi_\pm\rangle$, all defined on the orbitals l and lʹ

$$\begin{equation} \begin{array}{*{20}{l}} {\hat P_{ - - }^{l{l^\prime }} = \sum\limits_{i,j = l,{l^\prime }} {\sum\limits_{\sigma {\sigma ^\prime }} {{{\hat n}_{i\sigma }}} } {{\hat n}_{j{\sigma ^\prime }}}\left( {1 - {{\hat n}_{i\bar \sigma }}} \right)\left( {1 - {{\hat n}_{j\bar \sigma \prime }}} \right)}\\ {{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \hat P_{{\Psi _ \pm }}^{l{l^\prime }} = |\Psi _ \pm ^{l{l^\prime }}\rangle \langle \Psi _ \pm ^{l{l^\prime }}|}\\ {{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} = {{\hat n}_{l \uparrow }}{{\hat n}_{{l^\prime } \downarrow }}\left( {1 - {{\hat n}_{l \downarrow }}} \right)\left( {1 - {{\hat n}_{{l^\prime } \uparrow }}} \right) \pm f_{l \uparrow }^\dagger f_{{l^\prime } \downarrow }^\dagger f_{{l^\prime } \uparrow }^{}f_{l \downarrow }^{}.} \end{array} \end{equation} \tag{ 29 }$$

Here, $\hat{n}_{l\sigma}\! = \! f^{\,\dagger}_{l\sigma} f^{\phantom{\dagger}}_{l\sigma}$ denotes the spin-σ particle number operator on orbital l. After some straightforward but lengthy calculation, we obtain the following expressions for the parameters r and t, required in equation (15),

$$\begin{equation} \begin{array}{l} t = \max \left\{ {\langle {\Psi _N}|\hat P_{{\Psi _ \pm }}^{{\kern 1pt} l{l^\prime }}|{\Psi _N}\rangle } \right\} = A + B,\\ r = \langle {\Psi _N}|\hat P_{ - - }^{l{l^\prime }}|{\Psi _N}\rangle - t = 3A - 3B, \end{array} \end{equation} \tag{ 30 }$$

where

$$\begin{equation} A \equiv \left[\eta^2 - \eta - W\left(d,\eta\right)^2 \right]^2, \quad B \equiv W\left(d,\eta\right)^2. \end{equation} \tag{ 31 }$$

With this, the N-SSR entanglement can be neatly expressed as

$$\begin{equation} \begin{split} E\left(\tilde{\rho}^\mathrm{N}_{ll^{\prime}}\right) & = \left(A+B\right)\log\left(\frac{A+B}{2A-B}\right) \\ & \quad+ \left(3A-3B\right)\log\left(\frac{3A-3B}{2A-B}\right) \end{split} \end{equation} \tag{ 32 }$$

provided the entanglement criteria (17) (here it translates to $A \lt 2B$) is met, and otherwise it is zero. We present in figure 2 the orbital–orbital entanglement as a function of the filling fraction η, at various orbital separations d, together with the entanglement when only P-SSR is assumed, also obtained analytically (formula can be found in [40]).

Zoom In Zoom Out Reset image size

In this section we discuss the orbital–orbital entanglement in the tight-binding model based on (32). Figure 2 reveals the rich entanglement structure of the Slater determinant ground state. At different particle filling fraction η, the orbital–orbital entanglement behaves in qualitatively distinct ways, which we will now account for in detail.

First, we notice that the graph of orbital–orbital entanglement is symmetric around the half-filling point $\eta\! = \!1/2$. This is due to the unbroken particle–hole symmetry of the Hamiltonian (21) in the ground states. Moreover, using (32), it is easy to prove that the entanglement as a function of η is symmetric about $\eta\! = \!1/2$. Because of this, it suffices to consider only $0 \! \unicode{x2A7D} \! \eta \! \unicode{x2A7D} \! 1/2$, and treat it as a particle or hole filling fraction. In order to illustrate more details of this η-regime, we use a log–log scale in the right panel of figure 2.

Second, the analytic properties of the orbital–orbital entanglement highly depend on the combination of orbital separation d and the filling fraction η. The common feature here is the vanishing of orbital–orbital entanglement for all d at both, $\eta\! = \!0$ and $\eta\! = \!1$. In these two cases the ground state is invariant under any orbital rotation and contains zero correlation in any orbital basis. For two neighboring orbitals ($d\! = \!1$), their entanglement is maximized exactly when the chain is half filled $(\eta\! = \!1/2)$, for both P- and N-SSR entanglement. Around half-filling, the P-SSR entanglement is significantly higher than the N-SSR one. This is not surprising, as N-SSR imposes a stronger restriction on the accessible entanglement. Yet, at small particle or hole filling fractions, the two become essentially the same. The reason for this is that in such regimes the weight in the $P_{1,1}$ sector dominates over the weight in the $P_{2,0} + P_{0,2}$ sector, as the particles/holes tend to avoid sitting on the same orbital due to the repulsive interaction. This weight domination effect is stronger when the orbitals are further apart. For $d\!\unicode{x2A7E}\!2$, the P- and N-SSR entanglement is practically indistinguishable. Thus we can restrict the investigation of analytic properties to the latter by referring to (32) for all practical purposes here, and also in the following section. In this case, two maxima occur which appear symmetrically on the η-axis about $\eta\! = \!1/2$ (see the left panel of figure 2) due to the particle–hole symmetry of (32). At lower η, more and farther separated orbital pairs become entangled, though this entanglement grows vanishingly small with the filling fraction.

Third, we find that entanglement across extremely long distances can be engineered by tuning the filling fraction to be infinitesimally small. In the right panel of figure 2, we see more clearly on the log–log scale that for each separation d, there exists a critical value $\eta(d)$, below which the two orbitals are always entangled. At the same time, for any filling fraction η, there is a minimal disentangling distance $d_\mathrm{\min}(\eta)$, beyond which all orbital pairs are in a separable state. This phenomenon of completely vanishing entanglement is called the sudden death of entanglement [64]. Geometrically, it means nothing else than that the two-orbital reduced quantum state evolves along a trajectory which at some point enters the convex set of separable states. We shall explore this phenomenon in more details in the next section.

In this section we investigate the asymptotic behavior of the orbital–orbital entanglement in various limiting regimes. Previously we observed that interesting phenomena such as long-distance entanglement and sudden death of entanglement occur when the particle filling η or hole filling $1\!-\!\eta$ is small. Moreover, the P- and N-SSR entanglement in this region are practically identical. Therefore we can focus on the latter in our analysis.

For a fixed orbital separation d, when $\eta \ll d^{-1}$, we get from (28) $W(d,\eta) = \eta + \mathcal{O}(\eta^3)$. The N-SSR entanglement between two orbitals separated by d becomes

$$\begin{equation} \begin{split} E\left(\tilde{\rho}^\mathrm{N}_{ij}\right) = 2\log\left(2\right) \eta^2 + \mathcal{O}\left(\eta^5\right), \quad \eta \ll d^{-1}, \end{split} \end{equation} \tag{ 33 }$$

i.e. the N-SSR orbital–orbital entanglement is quadratic in η for small η. In log–log scale this relation becomes linear

$$\begin{equation} \log\left(E\left(\tilde{\rho}^\mathrm{N}_{ij}\right)\right) = 2\log\left(\eta\right) + \mathcal{O}\left(1\right), \quad \eta \ll d^{-1}. \end{equation} \tag{ 34 }$$

For the sake of completeness, we also verify that the P-SSR entanglement indeed reduces to the N-SSR one, i.e.

$$\begin{equation} \log\left(E\left(\tilde{\rho}^\mathrm{P}_{ij}\right)\right) = 2 \log\left(\eta\right) + \mathcal{O}\left(1\right), \quad \eta \ll d^{-1}. \end{equation} \tag{ 35 }$$

Equations (34) and (35) explain the asymptotic linear behavior in figure 2 (right panel) of both the P-SSR (red) and N-SSR (blue) entanglement on the log–log scale, with slope 2 valid for any separation d. This relation is more robust for small separation d, as for large separation the orbital–orbital entanglement quickly deviates from the linear asymptote and plunges into sudden death as η increases. On the other hand, for small η and $d \lt d_\mathrm{min}(\eta)$, the orbital–orbital entanglement is only a function of η, but not of d. All orbital pairs below the critical separation are equally entangled with one another. Due to the particle–hole symmetry, the orbital–orbital entanglement at small hole fraction $1\!-\!\eta \ll d^{-1}$ is also quadratic in $1\!-\!\eta$.

Above the critical separation $d_\mathrm{min}(\eta)$, all orbital pairs are disentangled. Based on figure 3, $d_\mathrm{min}(\eta)$ seems to diverge as η approaches 0. To understand the asymptotic behavior of $d_\mathrm{min}(\eta)$ at small η, we analyze the entanglement criteria (17). Recall from the previous section that the two orbitals are in a separable state if A and B in (31) satisfy the relation $A \unicode{x2A7E} 2B$, which in explicit terms reads

$$\begin{equation} \begin{split} \eta^2-\eta &\unicode{x2A7D} W\left(d,\eta\right)^2 - \sqrt{2} |W\left(d,\eta\right)|. \end{split} \end{equation} \tag{ 36 }$$

Although (36) cannot be solved analytically, in the limit of large d using the expansion $W(d,\eta) \sim (\pi d)^{-1}$ we can derive an analytic estimation for the minimal disentangling separation $d_\mathrm{min}(\eta)$

$$\begin{equation} d_{\min}\left(\eta\right) = \frac{\sqrt{2}}{\pi} \frac{1}{\eta\left(1\!-\!\eta\right)} + \mathcal{O}\left(\eta^0\right). \end{equation} \tag{ 37 }$$

In log–log scale, this relation becomes linear

$$\begin{equation} \begin{split} \log\left(d_\mathrm{min}\left(\eta\right)\right) = \begin{cases} -\frac{\sqrt{2}}{\pi}\log\left(\eta\right) + \mathcal{O}\left(\eta\right),&\quad\:\:\eta \ll 1, \\ -\frac{\sqrt{2}}{\pi}\log\left(1\!-\!\eta\right) + \mathcal{O}\left(1\!-\!\eta\right),&\,1\!-\!\eta \ll 1. \end{cases} \end{split} \end{equation} \tag{ 38 }$$

In figure 3, we plot the exact minimal disentangling distance obtained by numerically locating the sudden death of entanglement using (36), together with the analytic estimation (37). We can see the agreement is excellent, even up to the order 10⁻¹ in η.

Zoom In Zoom Out Reset image size

To summarize section 3, building on [40], we derived an analytic expression for the entanglement between any two spatial orbitals in the tight-binding model (21), as a function of both inter-orbital separation d and particle filling fraction η. In particular, equation (15) revealed the existence of long-distance entanglement when the filling fraction is low or close to 1. Moreover, our asymptotic analysis for small η revealed that this entanglement between orbitals is almost independent of the orbital separation d, and quadratic in the particle filling η for small η (and quadratic in hole filling $1-\eta$ for small $1-\eta$). For a fixed filling η, the entangling range is still finite, as we also observed the sudden death of entanglement when the two orbitals are separated beyond a critical distance. Utilizing our analytic expression, the leading order of the minimal disentangling distance $d_\mathrm{min}$ is extracted to be diverging as η⁻¹ when η is small, and $(1\!-\!\eta)^{-1}$ when η is close to 1.

We consider a finite periodic hydrogen ring, defined by the uniform nearest internuclear distance R, with the minimal STO-3G basis set. This basis contains one 1s orbital at each atomic center and it is the minimal setting for describing the dissociation of hydrogen molecules.

The electronic Hamiltonian contains three terms: the nuclear potential, the kinetic energy, and the Coulomb interaction. The corresponding Hamiltonian contains one- and two-electron terms $\hat{h}$ and $\hat{V}$, respectively, and takes the form

$$\begin{equation} \begin{split} &\hat{H} = \hat{h} + \hat{V} \\ &\:\:\: = \sum_{ij}\sum_\sigma {h}_{ij} f^{\,\dagger}_{i\sigma} f^{\phantom{\dagger}}_{j\sigma} + \sum_{ijkl} \sum_{\sigma\sigma^{\prime}} V_{ijkl} f^{\,\dagger}_{i\sigma} f^{\,\dagger}_{j\sigma^{\prime}} f^{\phantom{\dagger}}_{l\sigma^{\prime}}f^{\phantom{\dagger}}_{k\sigma}, \end{split} \end{equation} \tag{ 39 }$$

where

$$\begin{equation} \begin{split} h_{ij} & = \int \mathrm{d}x^3 \phi_{i}^\ast\left(x\right) \left( -\frac{\hbar}{2m_e} \nabla^2 + \sum_{m = 0}^{M-1} \frac{Z_m}{|x-\mathbf{X}_m|} \right) \phi_j\left(x\right), \\ V_{ijkl} & = \frac{1}{2}\int \mathrm{d}x^3 \mathrm{d}x^{^{\prime} 3} \phi_{i}^\ast\left(x\right) \phi_{j}^\ast\left(x^{\prime}\right) \frac{1}{|x-x^{\prime}|} \phi_{l}\left(x^{\prime}\right) \phi_k\left(x\right). \end{split} \end{equation} \tag{ 40 }$$

Here, φ_i 's are orthonormal molecular orbitals which form a complete basis of the one-particle Hilbert space spanned by the non-orthogonal atomic orbitals. M is the total number of nuclear centers with charges Z_m . In the case of $\mathrm{H}_{16}$, $Z_m\! = \!1$ for all $m\! = \!1,2,\ldots,16$ and X_m are defined as

$$\begin{equation} \begin{split} \mathbf{X}_{m} = \frac{R}{2\sin\left(\frac{\theta}{2}\right)} \left(\cos\left(m\theta\right),\sin\left(m \theta\right),0\right), \quad \theta = \frac{2\pi}{16}, \end{split} \end{equation} \tag{ 41 }$$

so that the 16 hydrogen atoms are arranged periodically on a ring, with the straight-line distance between two nearest neighbor atoms equal to R. From here on, we shall use atomic units (denoted by a.u.) for both distance and energy, which correspond to Bohr and Hartree, respectively.

The matrix elements of the one- and two-electron Hamiltonians are explicitly orbital basis dependent. For the purpose of QIP, we choose a set of localized orthonormal orbitals, obtained by symmetrically orthogonalizing [76] the local 1s atomic orbitals. In figure 4, we plot these localized orbitals for different nearest neighbor distances R. Locally they resemble the atomic orbitals, but they necessarily have finite contributions from atomic orbitals on other nuclear centers to ensure orthogonality. At $R\! = \!1$, each localized orbital at the corresponding atomic center still contains significant weight of the atomic orbitals on the nearest neighbors, while for $R \!\unicode{x2A7E}\! 3$, the localized orbitals are virtually single-centered.

Zoom In Zoom Out Reset image size

Do these localized orbitals exhibit a local interaction structure similar to the tight-binding chain? To answer this question, we plot the hopping terms h_ij for $i\!\neq\! j$ in figure 5. Recall that the hydrogen atoms are arranged on a ring, and therefore orbital 1 and 16 are in fact nearest neighbors, despite appearing on the opposite ends of the axes in figure 5. At small separation, e.g. $R\! = \!1$ (a.u.), hopping between distant localized orbitals is possible, but the hopping strength decays exponentially as internuclear distance increases. As we stretch the hydrogen ring, the hopping strength decays faster with the orbital separation. Finally at $R\! = \!5$ (a.u.), the hopping is strictly between nearest neighbors, just as in the tight-binding model (21).

Zoom In Zoom Out Reset image size

The calculation of the ground state of the hydrogen ring is done in three steps: (1) using the atomic wavefunction overlap matrix provided by Molpro [77–79], we symmetrically orthogonalize [76] the atomic orbitals into the localized orbitals, for nearest neighbor separation $R \! = \! 1,2,3,5$ (a.u.). (2) With a preceding Hartree–Fock calculation, we then transform the Hartree–Fock canonical orbitals into the localized orbitals as a post Hartree–Fock step, and compute the matrix elements of the one- and two-electron Hamiltonians followed by outputting them in FCIDUMP [80] file format, all within the framework of Molpro. (3) With the FCIDUMP file, quantum chemistry density matrix renormalization group (DMRG) calculations are performed to solve for the ground states, using the SyTen [81, 82] package, originally created by Claudius Hubig. The calculations are carried out for various electron numbers, initialized with random configuration states. The DMRG convergence is particularly challenging, because of both the strong intrinsic correlation in the system and the additional entanglement due to the localized orbital basis representation. For the data we use, the maximal bond dimension is set at $m\! = \!1000$. The final energies are chemically accurate, verified by the energy improvement below the order 10⁻³ (a.u.) in the DMRG calculations performed with bond dimension 2000.

In this section we present our main findings on the entanglement structure of the hydrogen ring $\mathrm{H}_{16}$. We tabulate in figure 6 the P-SSR entanglement between two localized orbitals for various nearest neighbor separations R and different numbers of electrons N in the ground state. For readability we exploit the discrete rotational symmetry of the system and only plot the entanglement associated to the localized orbital labeled 1 (represented as filled red circle). Entanglement involving other localized orbitals can be understood by relabeling them. Exact values of the P- and N-SSR entanglement can be found in table 1 in appendix. With these data, we observe three consequential features of the orbital–orbital entanglement between the localized orbitals in $\mathrm{H}_{16}$.

Zoom In Zoom Out Reset image size

First, we observe that at half-filling ($N\! = \!16$), the orbital–orbital entanglement is in complete agreement with the tight-binding model (21), namely only neighboring localized orbitals are entangled. In fact, this agreement extends beyond half-filling, as there is a region of filling fraction around $1/2$ where all orbital pairs separated by more than one lattice constant are disentangled, exactly as in figure 2. This is particularly surprising for systems with small nearest neighbor distance R, since there the hopping strength beyond nearest neighbors is rather considerable.

Second, at low or high filling, the hydrogen ring also displays long-distance entanglement similar to the tight-binding model. This phenomenon here differs from the latter, however, by preferring to entangle localized orbitals that are farther apart. Recall that in the tight-binding model, the orbital–orbital entanglement is uniform when long-distance entanglement is realized. This departure from the non-interacting scenario is best exemplified by the $N\! = \!2$ column in figure 6. In the ground state, due to the strong Coulomb repulsion, the two electrons would like to situate themselves as far apart as possible, forming singlet pairs at distant localized orbitals. Singlet states of such type, usually observed in a dissociated molecule, contain significant amount of entanglement. Singlet states realized on orbitals located closer to each other are of course possible, but they contribute to the total ground state with smaller weights. These two observations put together qualitatively explain the system's preference for long-distance entanglement at low particle filling. Furthermore, we notice that despite the lack of particle–hole symmetry in the system, the orbital entanglement in the $N\! = \!30$ column in figure 6 displays a similar behavior as in $N\! = \!2$. This can be qualitatively understood with the observation that at very low (or very high) filling, the particle–particle (or hole–hole) interaction becomes negligible compared to the one-particle terms. In the extreme case where the system hosts only one particle or one hole, the corresponding Hamiltonians are related to each other by a minus sign. Consequently, the two ground states are identical up to a global phase, and therefore share the same orbital entanglement pattern. For orbitals separated by large distances, this near particle–hole symmetry is even more apparent.

Third, the long-distance entanglement in the hydrogen ring is far more potent than that in the tight-binding model. In equation (33), we see that the long-distance entanglement in the tight-binding chain is quadratic in η for small η, which is obviously not the case here. This is best demonstrated by the entanglement between orbitals 1 and 9, which form the farthest separated pair, at $N\! = \!2$. The entanglement between them is of the order 10⁻², which is the same order as the nearest neighbor entanglement. In the tight-binding model, the available entanglement at this filling fraction would yield at the order 10⁻³. This yet again validates that interacting electronic systems contain significant entanglement, and the entanglement is shown to be naturally distributed across distant regions, making it ideal for realizing quantum information protocols.

To summarize the present section, we took the interaction between electrons into account by studying the hydrogen ring, and analyzed the orbital–orbital entanglement therein. We discovered that the entanglement between the localized orbitals on the hydrogen ring can to a large degree be rationalized using the much simpler tight-binding model: close to half-filling only nearest neighbor orbitals can be entangled, but away from half-filling, long-distance entanglement appears. Yet, in contrast to the vanishingly weak long-distance entanglement in the tight-binding model, the entanglement between two localized orbitals sitting on far-separated nuclei is much more potent, at an order comparable to the nearest neighbor entanglement. This surprising feature, conjectured to be a result of Coulomb interaction, suggests that chain-like molecules can indeed be considered as reservoirs of entanglement.