Pairwise entanglement and the Mott transition for correlated electrons in nanochains

The analysis starts from the second-quantized Hamiltonian, which can be written in the form

$$\begin{eqnarray}&&\hat{{ \mathcal H }}(\alpha ,R)={\epsilon }_{a}\displaystyle \sum _{j}{\hat{n}}_{j}+\displaystyle \sum _{i\lt j,\sigma }{t}_{{ij}}({\hat{c}}_{i\sigma }^{\dagger }{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{\dagger }{\hat{c}}_{i\sigma })+U\displaystyle \sum _{i}{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }+\displaystyle \sum _{i\lt j}{K}_{{ij}}{\hat{n}}_{i}{\hat{n}}_{j}+{V}_{\mathrm{ion}-\mathrm{ion}},\end{eqnarray} \tag{ 1 }$$

where ${\epsilon }_{a}$ is the atomic energy (same for all N sites upon applying periodic boundary conditions), t_ij is the hopping integral between ith and jth site (we further set ${t}_{{ij}}\equiv -t$ if i and j are nearest neighbors, otherwise ${t}_{{ij}}\approx 0$), U is the intrasite Coulomb repulsion, K_ij is the intersite Coulomb repulsion, and ${V}_{\mathrm{ion}-\mathrm{ion}}={\sum }_{i\lt j}{e}^{2}/{R}_{{ij}}$ (with the distance ${R}_{{ij}}\equiv R\min (| i-j| ,N-| i-j| )$) expresses the repulsion of infinite-mass ions.

The single-particle and interaction parameters of the Hamiltonian $\hat{{ \mathcal H }}(\alpha ,R)$ (1), also marked in figure 1, can be defined as follows

$$\begin{eqnarray}&&\langle {w}_{i}| T| {w}_{i}\rangle =\displaystyle \int {{\rm{d}}}^{3}r\,{w}_{i}^{\star }({\bf{r}})T({\bf{r}}){w}_{j}({\bf{r}})\equiv {\epsilon }_{a}{\delta }_{{ij}}+{t}_{{ij}}(1-{\delta }_{{ij}}),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\langle {w}_{i}{w}_{j}| V| {w}_{i}{w}_{j}\rangle =\displaystyle \int {{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}r^{\prime} \,| w({\bf{r}}){| }^{2}V({\bf{r}}-{\bf{r}}^{\prime} )| w({\bf{r}}^{\prime} ){| }^{2}\equiv U{\delta }_{{ij}}+{K}_{{ij}}(1-{\delta }_{{ij}}),\end{eqnarray} \tag{ 3 }$$

where T is the single particle Hamiltonian describing an electron in the medium of periodically arranged ions, and V is the Coulomb repulsive interaction of two electrons. The Wannier functions are defined via atomic (Slater-type) functions, namely

$$\begin{eqnarray}&&{w}_{i}({\bf{r}})=\displaystyle \sum _{j}{\beta }_{{ij}}{\psi }_{j}({\bf{r}}),\end{eqnarray} \tag{ 4 }$$

where the Slater 1s function ${\psi }_{i}({\bf{r}})={({\alpha }^{3}/\pi )}^{1/2}\exp (-\alpha | {\bf{r}}-{{\bf{R}}}_{i}| )$, with α being the inverse orbital size, here taken as variational parameter, and ${{\bf{R}}}_{i}$ being the position of ith ion. The coefficients ${\beta }_{{ij}}$ in equation (4) can be uniquely defined by imposing that (i) $\langle {w}_{i}| {w}_{k}\rangle ={\delta }_{{ik}}$ and that (ii) $\langle {w}_{i}| {\psi }_{i}\rangle$ is maximal [3, 55].

The range of Coulomb interactions in the Hamiltonian $\hat{{ \mathcal H }}(\alpha ,R)$ (1), quantified by parameters ${K}_{{ij}}\approx {e}^{2}/{R}_{{ij}}$ for ${R}_{{ij}}\gg R$, is a priori limited only by the chain length $L\equiv {NR}$. Previous studies on linear chains, both ab initio ones [12, 48] and these starting from second-quantized model Hamiltonians [50, 52], usually have imposed some form of charge screening reducing the range of such interactions (leaving the on-site Coulomb repulsion only in the extreme case [53]). This can be partly justified by possible influence of the enviroment in condensed-matter realizations such as quantum wires [54] or self-organized chains [56]. A slightly different situation occurrs for cold atom systems, where long-range dipole–dipole interactions may be relevant [57]. Apart from these experiment-related premises, long-range interactions usually lead to a noticeable slowdown of the convergence for numerical methods such as DMRG or QMC [58, 59]. This is not the case for the EDABI method, as the determination of parameters K_ij corresponds to a relatively small portion of the overall computation time [60]² , and the convergence of subsequent diagonalization in the Fock space (see the next subsection) is unaffected by the fact whether or not long-range interactions are included. Also, as charge screening becomes systematically less effective when the system dimensionality is reduced [61], one can expect it to be insignificant in a hypotetical realization of chains containing $N\leqslant 20$ atoms.

Next, each Slater function is approximated as follows

$$\begin{eqnarray}&&{\psi }_{i}({\bf{r}})\approx {\alpha }^{3/2}\displaystyle \sum _{q=1}^{p}{B}_{q}{\left(\displaystyle \frac{2{{\rm{\Gamma }}}_{q}^{2}}{\pi }\right)}^{3/4}{{\rm{e}}}^{-{\alpha }^{2}{{\rm{\Gamma }}}_{q}^{2}| {\bf{r}}-{{\bf{R}}}_{i}{| }^{2}},\end{eqnarray} \tag{ 5 }$$

where p is the number of Gaussian functions truncating ${\psi }_{j}$, and $\{{B}_{q},{{\rm{\Gamma }}}_{q}\}$, q = 1,..., p, are adjustable parameters chosen to minimize the atomic energy for $\alpha =1$ and a given value of p. Here we set p = 3, for which the deviation from the exact energy for $1s$ function is lower that 1%³ . Subsequently, the parameters ${\epsilon }_{a}$, t_ij, U, and K_ij are calculated from equations (2) and (3) as functions of α and R. The Hamiltonian $\hat{{ \mathcal H }}(\alpha ,R)$ (1) is diagonalized numerically in the Fock space, using the Lanczos algorithm, and the orbital size is optimized to find $\alpha ={\alpha }_{\min }$ corresponding to the minimal GS energy E_G(N) for each R. Defining the efective atomic energy as

$$\begin{eqnarray}&&{\epsilon }_{a}^{\mathrm{eff}}={\epsilon }_{a}+\displaystyle \frac{1}{N}\left(\displaystyle \sum _{i\lt j}{K}_{{ij}}+{V}_{\mathrm{ion}-\mathrm{ion}}\right),\end{eqnarray} \tag{ 6 }$$

one finds that ${\alpha }_{\min }$ and other parameters converges rapidly with N. This observation allows us to speed up computations by using the values obtained for smaller systems ($N\geqslant 10$) to perform extrapolation with $1/N\to 0$. (The results are listed in table 1.)

Table 1.  Microscopic parameters (specified in eV—unless stated otherwise) of the Hamiltonian $\hat{{ \mathcal H }}(\alpha ,R)$ (1) with $\alpha ={\alpha }_{\min }$ minimizing the ground-state (GS) energy E_G. The Bohr radius is ${a}_{0}=0.529\,$Å. The effective atomic level ${\epsilon }_{a}^{\mathrm{eff}}$ is defined by equation (6); the remaining symbols are $t\equiv -{t}_{{ij}}$ for $j=i\pm 1$ (mod N), ${K}_{m}\equiv {K}_{| i-j| }$ (with $m=1,2,3$), and the correlated hopping integral $V\equiv \langle {w}_{i}{w}_{i}| V| {w}_{i}{w}_{i\pm 1}\rangle$, quantifying the largest neglected term in $\hat{{ \mathcal H }}(\alpha ,R)$. The numerical extrapolation with $1/N\to 0$ is performed for all parameters.

$\ R/{a}_{0}\$	$\ {\alpha }_{\min }{a}_{0}\$	$\ {\epsilon }_{a}^{\mathrm{eff}}\$	$\ t\$	$\ U\$	$\ {K}_{1}\$	$\ {K}_{2}\$	$\ {K}_{3}\$	$\ V\$	$\ {E}_{G}/N\$
1.5	1.363	1.36	11.31	27.95	15.85	9.08	6.08	−0.597	−10.19
2.0	1.220	−7.48	6.02	23.58	12.40	6.82	4.54	−0.324	−12.65
2.5	1.122	10.85	3.60	20.83	10.20	5.46	3.63	−0.203	−13.32
3.0	1.062	12.27	2.32	19.14	8.69	4.54	3.02	−0.150	−13.48
3.5	1.031	12.90	1.57	18.16	7.58	3.89	2.60	−0.128	−13.50
4.0	1.013	13.20	1.08	17.57	6.71	3.40	2.27	−0.119	−13.50
5.0	1.004	13.43	0.51	17.12	5.43	2.72	1.81	−0.096	−13.50
6.0	1.001	13.48	0.23	16.99	4.53	2.27	1.51	−0.058	−13.50
7.0	1.000	13.49	0.10	16.97	3.89	1.99	1.29	−0.027	−13.50

Let us consider a general example at first: the quantum system, which can be divided into two distinct subsystems A and B. A pure state $| {\rm{\Psi }}\rangle$ can be represented as follows

$$\begin{eqnarray}&&| {\rm{\Psi }}\rangle =\displaystyle \sum _{\alpha \beta }{{\rm{\Psi }}}_{\alpha \beta }| \alpha \rangle \otimes | \beta \rangle ,\end{eqnarray} \tag{ 11 }$$

where ${{\rm{\Psi }}}_{\alpha \beta }$ denotes a complex probability amplitude corresponding to a basis state of the full system $| \alpha \rangle \otimes | \beta \rangle$, while $\{| \alpha \rangle \}$ and $\{| \beta \rangle \}$ are complete basis sets for the subsystems A and B (respectively). The reduced density operator ${\hat{\rho }}_{A}$ is defined as:

$$\begin{eqnarray}{\hat{\rho }}_{A} & = & {{\rm{Tr}}}_{A}| {\rm{\Psi }}\rangle \langle {\rm{\Psi }}| \equiv \displaystyle \sum _{\beta }\langle \beta | {\rm{\Psi }}\rangle \langle {\rm{\Psi }}| \beta \rangle \\ & = & \displaystyle \sum _{\beta ,\alpha ,\alpha ^{\prime} }{{\rm{\Psi }}}_{\alpha \beta }{{\rm{\Psi }}}_{\alpha ^{\prime} \beta }^{\star }| \alpha \rangle \langle \alpha ^{\prime} | ,\end{eqnarray} \tag{ 12 }$$

where the last equality follows from equation (11) and ${{\rm{\Psi }}}_{\alpha \beta }^{\star }$ denotes the complex conjugate of ${{\rm{\Psi }}}_{\alpha \beta }$. Subsequently, the matrix elements of ${\hat{\rho }}_{A}$ are given by

$$\begin{eqnarray}&&{\rho }_{\alpha ,\alpha ^{\prime} }\equiv \langle \alpha | {\hat{\rho }}_{A}| \alpha ^{\prime} \rangle =\displaystyle \sum _{\beta }{{\rm{\Psi }}}_{\alpha \beta }{{\rm{\Psi }}}_{\alpha ^{\prime} ,\beta }^{\star }.\end{eqnarray} \tag{ 13 }$$

We define now projection operators P_α and P_β, associated with the basis states $| \alpha \rangle$ and $| \beta \rangle$, via

$$\begin{eqnarray}&&{P}_{\alpha }| \alpha ^{\prime} \rangle ={\delta }_{\alpha \alpha ^{\prime} }| \alpha \rangle ,\ \ \ {P}_{\beta }| \beta ^{\prime} \rangle ={\delta }_{\beta \beta ^{\prime} }| \beta \rangle .\end{eqnarray} \tag{ 14 }$$

Useful properties, directly following from equation (14) are

$$\begin{eqnarray}&&{P}_{\alpha }^{2}={P}_{\alpha },\ \ \ {P}_{\beta }^{2}={P}_{\beta },\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{\alpha }{P}_{\alpha }=\displaystyle \sum _{\beta }{P}_{\beta }=1,\end{eqnarray} \tag{ 16 }$$

where equation (16) also employs the completeness of the basis sets $\{| \alpha \rangle \}$ and $\{| \beta \rangle \}$. Next, we define the transfer operator

$$\begin{eqnarray}&&{T}_{\alpha \alpha ^{\prime} }| \alpha ^{\prime} \rangle =| \alpha \rangle ,\ \ \ \langle \alpha ^{\prime} | {T}_{\alpha \alpha ^{\prime} }^{\dagger }=\langle \alpha | ,\end{eqnarray} \tag{ 17 }$$

which is unitary, i.e.

$$\begin{eqnarray}&&{T}_{\alpha \alpha ^{\prime} }^{\dagger }={({T}_{\alpha \alpha ^{\prime} })}^{-1}\equiv {T}_{\alpha ^{\prime} \alpha }.\end{eqnarray} \tag{ 18 }$$

Explicit forms of P_α, P_β, and ${T}_{\alpha \alpha ^{\prime} }$, corresponding to the particular splittings of the system into A and B subsystems, are to be specified later in terms of the operators $\{{\hat{c}}_{i\sigma },{\hat{c}}_{i\sigma }^{\dagger }\}$ for lattice spin-1/2 fermions. Here we notice that as the indices $(\alpha ,\alpha ^{\prime} )$ refer the subsystem A, whereas the index β refers to the subsystem B, we have

$$\begin{eqnarray}&&[{P}_{\alpha },{P}_{\beta }]=0,\ \ \ [{P}_{\beta },{T}_{\alpha \alpha ^{\prime} }]=0,\end{eqnarray} \tag{ 19 }$$

for all $(\alpha ,\alpha ^{\prime} )$ and β.

With the help of operators ${P}_{\alpha }$ and ${P}_{\beta }$ (14) one can write down, for $| {\rm{\Psi }}\rangle$ represented according to equation (11),

$$\begin{eqnarray}&&{P}_{\alpha }{P}_{\beta }| {\rm{\Psi }}\rangle ={{\rm{\Psi }}}_{\alpha \beta }\,| \alpha \rangle \otimes | \beta \rangle ,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\langle {\rm{\Psi }}| {P}_{\alpha ^{\prime} }{P}_{\beta }=\langle \alpha ^{\prime} | \otimes \langle \beta | \,{{\rm{\Psi }}}_{\alpha ^{\prime} \beta }^{\star }.\end{eqnarray} \tag{ 21 }$$

In turn, the reduced density matrix, as given by the rightmost equality in equation (13), can be rewritten as

$$\begin{eqnarray}{\rho }_{\alpha \alpha ^{\prime} } & = & \displaystyle \sum _{\beta }\langle {\rm{\Psi }}| {P}_{\alpha ^{\prime} }{P}_{\beta }{T}_{\alpha ^{\prime} \alpha }{P}_{\alpha }{P}_{\beta }| {\rm{\Psi }}\rangle \\ & = & \langle {\rm{\Psi }}| {P}_{\alpha ^{\prime} }{T}_{\alpha ^{\prime} \alpha }{P}_{\alpha }| {\rm{\Psi }}\rangle ,\end{eqnarray} \tag{ 22 }$$

with the last equality following from equations (15), (16), and (19). Finally, using equation (18), we find

$$\begin{eqnarray}&&{({\rho }_{\alpha \alpha ^{\prime} })}^{\star }=\langle {\rm{\Psi }}| {P}_{\alpha }{T}_{\alpha \alpha ^{\prime} }{P}_{\alpha ^{\prime} }| {\rm{\Psi }}\rangle .\end{eqnarray} \tag{ 23 }$$

Remarkably, the reduced density matrix elements in equation (23) are expressed by pure-state expectation values (correlation functions) of the operators acting only on the subsystem A. This feature is further explored in the remaining parts of this section.

The local entanglement [27] exhibits quantum correlations between the local state of a selected jth site (subsystem A) and the rest of the system (B). As the basis set $\{| \alpha \rangle \}$ is simply given by equation (10), one can write down the corresponding projection operators

$$\begin{eqnarray}&&({P}_{\alpha })={[(1-{n}_{\uparrow })(1-{n}_{\downarrow }),{n}_{\uparrow }(1-{n}_{\downarrow }),(1-{n}_{\uparrow }){n}_{\downarrow },{n}_{\uparrow }{n}_{\downarrow }]}^{T},\end{eqnarray} \tag{ 24 }$$

and the transfer operators

$$\begin{eqnarray}({T}_{\alpha \alpha ^{\prime} })=\left(\begin{array}{cccc}1 & {c}_{\uparrow } & {c}_{\downarrow } & {c}_{\downarrow }{c}_{\uparrow }\\ {c}_{\uparrow }^{\dagger } & 1 & {c}_{\uparrow }^{\dagger }{c}_{\downarrow } & -{c}_{\downarrow }\\ {c}_{\downarrow }^{\dagger } & {c}_{\downarrow }^{\dagger }{c}_{\uparrow } & 1 & {c}_{\uparrow }\\ {c}_{\uparrow }^{\dagger }{c}_{\downarrow }^{\dagger } & -{c}_{\downarrow }^{\dagger } & {c}_{\uparrow }^{\dagger } & 1\end{array}\right),\end{eqnarray} \tag{ 25 }$$

where we have omitted the site index j for brevity. For the system with the total spin and charge conservation, such as our linear chain (see equation (1)), we have

$$\begin{eqnarray}&&\langle {c}_{\sigma }\rangle =\langle {c}_{\sigma }^{\dagger }\rangle =0,\ \ \ \langle {c}_{\sigma }{c}_{\bar{\sigma }}\rangle =\langle {c}_{\sigma }^{\dagger }{c}_{\bar{\sigma }}^{\dagger }\rangle =0,\\ &&\quad \mathrm{and}\ \ \ \ \langle {c}_{\sigma }^{\dagger }{c}_{\bar{\sigma }}\rangle =\langle {c}_{\bar{\sigma }}^{\dagger }{c}_{\sigma }\rangle =0,\end{eqnarray} \tag{ 26 }$$

(with $\overline{\sigma }$ being the spin index opposite to σ) providing that the averaging takes place over an eigenstate $| {\rm{\Psi }}\rangle$ of the system Hamiltonian (not necessarily the GS). Substituting the expressions for P_α (24) and ${T}_{\alpha \alpha ^{\prime} }$ (25) into equation (23) and taking equation (26) into account we get, after a straightforward algebra,

$$\begin{eqnarray}&&({\rho }_{\alpha \alpha ^{\prime} })={\rm{diag}}({u}_{+},{w}_{1},{w}_{2},{u}_{-})\end{eqnarray} \tag{ 27 }$$

with

$$\begin{eqnarray}&&{u}_{+}=\langle (1-{n}_{j\uparrow })(1-{n}_{j\downarrow })\rangle ,\ \ \ {w}_{1}=\langle {n}_{j\uparrow }(1-{n}_{j\downarrow })\rangle ,\\ &&\quad {w}_{2}=\langle (1-{n}_{j\uparrow }){n}_{j\downarrow }\rangle ,\ \ \ {u}_{-}=\langle {n}_{j\uparrow }{n}_{j\downarrow }\rangle .\end{eqnarray} \tag{ 28 }$$

Equation (27) corresponds to the reduced density operator of the well-known form

$$\begin{eqnarray}&&{\hat{\rho }}_{A}=\displaystyle \sum _{\alpha \alpha ^{\prime} }{\rho }_{\alpha \alpha ^{\prime} }| \alpha \rangle \langle \alpha ^{\prime} | ={u}_{+}| 0\rangle \langle 0| +{w}_{1}| \uparrow \rangle \langle \uparrow | +{w}_{2}| \downarrow \rangle \langle \downarrow | +{u}_{-}| \uparrow \downarrow \rangle \langle \uparrow \downarrow | ,\end{eqnarray} \tag{ 29 }$$

where we have used the basis given by equation (10).

A quantitative measure of the entanglement between the state of jth site and that of the remaining $N-1$ sites is given by the von Neumann entropy

$$\begin{eqnarray}&&{E}_{v}=-{u}_{+}{\mathrm{log}}_{2}{u}_{+}-{w}_{1}{\mathrm{log}}_{2}{w}_{1}-{w}_{2}{\mathrm{log}}_{2}{w}_{2}-{u}_{-}{\mathrm{log}}_{2}{u}_{-}.\end{eqnarray} \tag{ 30 }$$

If the translational invariance of the system is imposed, one can define the particle density n and the average number of double occupancies d, following

$$\begin{eqnarray}&&\langle {n}_{j\uparrow }\rangle =\langle {n}_{j\uparrow }\rangle =\displaystyle \frac{{N}_{\mathrm{el}}}{2N}\equiv \displaystyle \frac{n}{2},\ \ \ \ \ \langle {n}_{j\uparrow }{n}_{j\uparrow }\rangle \equiv d,\end{eqnarray} \tag{ 31 }$$

where we have further assumed that the averaging takes place over an eigenstate characterized by the total zth component of spin S_z = 0. Subsequently, equation (28) can be rewritten as

$$\begin{eqnarray}&&{u}_{+}=1-n+d,\ \ \ {w}_{1}={w}_{2}=\displaystyle \frac{n}{2}-d,\ \ \ {u}_{-}=d.\end{eqnarray} \tag{ 32 }$$

In particular, for the GS of the Hamiltonian $\hat{{ \mathcal H }}(\alpha ,R)$ (1), for which the effective interaction between electrons can roughly be estimated by $\sim U-{K}_{1}$ and is always repulsive, one can show that

$$\begin{eqnarray}&&\displaystyle \frac{| n-1| +n-1}{2}\leqslant d\leqslant \displaystyle \frac{{n}^{2}}{4}\ \ \ \ \mathrm{for}\ \ \ \ 0\leqslant n\leqslant 2,\end{eqnarray} \tag{ 33 }$$

where the lower bound for d corresponds to the so-called strong-correlations limit approached for $R\gg {a}_{0}$ (as the bandwidth $W\equiv 4t\ll U-{K}_{1}$, see table 1), whereas the upper bound for d corresponds to the free-electrons limit approached for $R\lesssim 2{a}_{0}$ ($W\gg U-{K}_{1}$). As the charge fluctuations are given by ${\rm{Var}}\{{n}_{j}\}=\langle {n}_{j}^{2}\rangle -\langle {n}_{j}{\rangle }^{2}\equiv n-{n}^{2}+2d$, the above-mentioned limits coincide (respectively) with the minimal and the maximal fluctuations for a given n. Moreover, the bounds for d in equation (33) can be mapped, via equation (32), onto

$$\begin{eqnarray}&&{E}_{v}^{\mathrm{corr}}(n)\leqslant {E}_{v}\leqslant {E}_{v}^{\mathrm{free}}(n),\end{eqnarray} \tag{ 34 }$$

with

$$\begin{eqnarray}&&{E}_{v}^{\mathrm{corr}}(n)=-| 1-n| {\mathrm{log}}_{2}| 1-n| -(1-| 1-n| ){\mathrm{log}}_{2}\left(\displaystyle \frac{1-| 1-n| }{2}\right)\end{eqnarray} \tag{ 35 }$$

for the strong-correlations limit, and

$$\begin{eqnarray}&&{E}_{v}^{\mathrm{free}}(n)=-n{\mathrm{log}}_{2}\displaystyle \frac{n}{2}-(2-n){\mathrm{log}}_{2}\left(1-\displaystyle \frac{n}{2}\right)\end{eqnarray} \tag{ 36 }$$

for the free-electrons limit. For instance, equations (34)–(36) lead to

$$\begin{eqnarray}&&1\leqslant {E}_{v}\leqslant 2\ \ \ \mathrm{for}\ \ \ {N}_{\mathrm{el}}=N,\end{eqnarray} \tag{ 37 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{3}{2}\leqslant {E}_{v}\leqslant 4-\displaystyle \frac{3}{2}{\mathrm{log}}_{2}3\approx 1.623\ \ \ \mathrm{for}\ \ \ {N}_{\mathrm{el}}=\displaystyle \frac{N}{2}.\end{eqnarray} \tag{ 38 }$$

Remarkably, the strong-correlations limit corresponds to the minimal E_v, whereas the free-electrons limit corresponds to the maximal E_v (for any $0\leqslant n\leqslant 2$). This is because when the free-electrons limit is approached, ground-state wavefunction can be truncated by a single Slater determinant composed of fully delocalized single-particle functions (Bloch functions), resulting in the maximal entanglement between jth site and the remaining $N-1$ sites [12, 36]. Repulsive interactions between electrons lead to the correlation-induced suppression of ${\rm{Var}}\{{n}_{j}\}$, and to the decreasing E_v with growing the interaction-to-bandwidth ratio.

General findings, presented briefly in the above, are now illustrated with the numerical examples for ${N}_{\mathrm{el}}=N$ and ${N}_{\mathrm{el}}=N/2$ (see figures 4 and 5). In both cases, the values of E_v are close to the upper bounds given by equations (37) and (38) for the smallest considered value of $R/{a}_{0}=1.5$, and systematically decrease with growing R, gradually approaching the lower bounds in equations (37) and (38). Also, a remarkably fast convergence of E_v with growing N is observed for any R, making it necessary to apply vertical shifts to the datasets in figures 4 and 5. We further notice that ${E}_{v}\approx 3/2$ starting from $R\gtrsim 5\,{a}_{0}$ for ${N}_{\mathrm{el}}=N/2$, whereas for ${N}_{\mathrm{el}}=N$ we still have ${E}_{v}\gt 1$ in this range. However, the smooth evolution of E_v with R is observed for both ${N}_{\mathrm{el}}=N$ and ${N}_{\mathrm{el}}=N/2$, making it difficult to consider the local entanglement as an estimate of the Mott transition. Pairwise entanglement is discussed next as the other candidate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We consider now the subsystem A consists of two spatially-separate quantum bits (qubits), one of which is associated with ith lattice site and the other with jth site. Each individual qubit can be realized employing charge or spin degrees of freedom.

For instance, if charge qubits are under consideration, one can choose the basis set for A in a fixed-spin sector $\sigma =\uparrow$, namely

$$\begin{eqnarray}&&{\{| \alpha \rangle \}}^{c}=\{\,| 0{\rangle }_{i}\otimes | 0{\rangle }_{j},| 0{\rangle }_{i}\otimes {| \uparrow \rangle }_{j},{| \uparrow \rangle }_{i}\otimes | 0{\rangle }_{j},{| \uparrow \rangle }_{i}\otimes {| \uparrow \rangle }_{j}\,\},\end{eqnarray} \tag{ 39 }$$

whereas for spin qubits we have

$$\begin{eqnarray}&&{\{| \alpha \rangle \}}^{s}=\{\,{| \uparrow \rangle }_{i}\otimes {| \uparrow \rangle }_{j},{| \uparrow \rangle }_{i}\otimes {| \downarrow \rangle }_{j},{| \downarrow \rangle }_{i}\otimes {| \uparrow \rangle }_{j},{| \downarrow \rangle }_{i}\otimes {| \downarrow \rangle }_{j}\,\}.\end{eqnarray} \tag{ 40 }$$

The corresponding projection operators are

$$\begin{eqnarray}&&({P}_{\alpha }^{\,c\,})={[(1-{n}_{i\uparrow })(1-{n}_{j\uparrow }),(1-{n}_{i\uparrow }){n}_{j\uparrow },{n}_{i\uparrow }(1-{n}_{j\uparrow }),{n}_{i\uparrow }{n}_{j\uparrow }]}^{T},\end{eqnarray} \tag{ 41 }$$

and

$$\begin{eqnarray}({P}_{\alpha }^{\,s\,}) & = & [{n}_{i\uparrow }(1-{n}_{i\downarrow }){n}_{j\uparrow }(1-{n}_{j\downarrow }),\quad {n}_{i\uparrow }(1-{n}_{i\downarrow })(1-{n}_{j\uparrow }){n}_{j\downarrow },\\ \quad & & {(1-{n}_{i\uparrow }){n}_{i\downarrow }{n}_{j\uparrow }(1-{n}_{j\downarrow }),(1-{n}_{i\uparrow }){n}_{i\downarrow }(1-{n}_{j\uparrow }){n}_{j\downarrow }]}^{T},\end{eqnarray} \tag{ 42 }$$

with the upper indices ($c,s$) referring to charge and spin qubits (respectively). Subsequently, the transfer operators are given by

$$\begin{eqnarray}({T}_{\alpha \alpha ^{\prime} }^{\,c})=\left(\begin{array}{cccc}1 & {c}_{j\uparrow } & {c}_{i\uparrow } & {c}_{j\uparrow }{c}_{i\uparrow }\\ {c}_{j\uparrow }^{\dagger } & 1 & {c}_{j\uparrow }^{\dagger }{c}_{i\uparrow } & {c}_{i\uparrow }\\ {c}_{i\uparrow }^{\dagger } & {c}_{i\uparrow }^{\dagger }{c}_{j\uparrow } & 1 & -{c}_{j\uparrow }\\ {c}_{i\uparrow }^{\dagger }{c}_{j\uparrow }^{\dagger } & {c}_{i\uparrow }^{\dagger } & -{c}_{j\uparrow }^{\dagger } & 1\end{array}\right),\end{eqnarray} \tag{ 43 }$$

and

$$\begin{eqnarray}({T}_{\alpha \alpha ^{\prime} }^{\,s})=\left(\begin{array}{cccc}1 & {S}_{j}^{+} & {S}_{i}^{+} & {S}_{j}^{+}{S}_{i}^{+}\\ {S}_{j}^{-} & 1 & {S}_{j}^{-}{S}_{i}^{+} & {S}_{i}^{+}\\ {S}_{i}^{-} & {S}_{i}^{-}{S}_{j}^{+} & 1 & {S}_{j}^{+}\\ {S}_{i}^{-}{S}_{j}^{-} & {S}_{i}^{-} & {S}_{j}^{-} & 1\end{array}\right),\end{eqnarray} \tag{ 44 }$$

with the spin operators ${S}_{i}^{+}={c}_{i\uparrow }^{\dagger }{c}_{i\downarrow }$, ${S}_{i}^{-}={c}_{i\downarrow }^{\dagger }{c}_{i\uparrow }$. Substituting the above expressions into equation (23) we get

$$\begin{eqnarray}({\rho }_{\alpha \alpha ^{\prime} }^{X})=\left(\begin{array}{cccc}{u}_{+}^{X} & 0 & 0 & 0\\ 0 & {w}_{1}^{X} & {z}^{X} & 0\\ 0 & {({z}^{X})}^{\star } & {w}_{2}^{X} & 0\\ 0 & 0 & 0 & {u}_{-}^{X}\end{array}\right)\ \ \ \ \mathrm{for}\ X=c,s,\end{eqnarray} \tag{ 45 }$$

where the total spin and charge conservation is imposed. The nonzero matrix elements in equation (45) are given by

$$\begin{eqnarray}&&{u}_{+}^{c}=\langle (1-{n}_{i\uparrow })(1-{n}_{j\uparrow })\rangle ,\ \ \ {w}_{1}^{c}=\langle {n}_{i\uparrow }(1-{n}_{j\uparrow })\rangle ,\\ &&\quad {z}^{c}=\langle {c}_{j\uparrow }^{\dagger }{c}_{i\uparrow }\rangle ,\ \ \ {w}_{2}^{c}=\langle (1-{n}_{i\uparrow }){n}_{j\uparrow }\rangle ,{u}_{-}^{c}=\langle {n}_{i\uparrow }{n}_{j\uparrow }\rangle ,\end{eqnarray} \tag{ 46 }$$

for charge qubits, or by

$$\begin{eqnarray}{u}_{+}^{s} & = & \langle {n}_{i\uparrow }(1-{n}_{i\downarrow }){n}_{j\uparrow }(1-{n}_{j\downarrow })\rangle ,\\ {w}_{1}^{s} & = & \langle {n}_{i\uparrow }(1-{n}_{i\downarrow })(1-{n}_{j\uparrow }){n}_{j\downarrow }\rangle ,\\ {z}^{s} & = & \langle {S}_{j}^{+}{S}_{i}^{-}\rangle =\langle {c}_{j\uparrow }^{\dagger }{c}_{j\downarrow }{c}_{j\downarrow }^{\dagger }{c}_{j\uparrow }\rangle ,\\ {w}_{2}^{s} & = & \langle (1-{n}_{i\uparrow }){n}_{i\downarrow }{n}_{j\uparrow }(1-{n}_{j\downarrow })\rangle ,\\ {u}_{-}^{s} & = & \langle (1-{n}_{i\uparrow }){n}_{i\downarrow }(1-{n}_{j\uparrow }){n}_{j\downarrow }\rangle ,\end{eqnarray} \tag{ 47 }$$

for spin qubits.

We use now the concurrence ${ \mathcal C }$, as a quantitative measure of quantum entanglement in the two-qubit subsystem A. The closed-form expression was derived by Wootters [22] and reads

$$\begin{eqnarray}&&{ \mathcal C }=\max \{0,\sqrt{{\lambda }_{1}}-\sqrt{{\lambda }_{2}}-\sqrt{{\lambda }_{3}}-\sqrt{{\lambda }_{4}}\},\end{eqnarray} \tag{ 48 }$$

where ${\lambda }_{1}\geqslant {\lambda }_{2}\geqslant {\lambda }_{3}\geqslant {\lambda }_{4}$ are eigenvalues of the matrix product

$$\begin{eqnarray}&&\rho \cdot ({\sigma }_{i}^{y}\otimes {\sigma }_{j}^{y}){\rho }^{* }({\sigma }_{i}^{y}\otimes {\sigma }_{j}^{y})\end{eqnarray} \tag{ 49 }$$

with $\rho =(\,{\rho }_{\alpha \alpha ^{\prime} }^{X})$ given by equations (45)–(47), and ${\sigma }_{i}^{y}$ (${\sigma }_{j}^{y}$) being the second Pauli matrix acting on the qubit associated with ith (jth) lattice site. It was also shown in [22] that the entanglement of formation for a pair of qubits was uniquely determined by ${ \mathcal C }$, namely

$$\begin{eqnarray}&&{E}_{f}=-\xi {\mathrm{log}}_{2}\xi -(1-\xi ){\mathrm{log}}_{2}(1-\xi ),\end{eqnarray} \tag{ 50 }$$

where $\xi =(1+\sqrt{1-{{ \mathcal C }}^{2}})/2$. In fact, ${ \mathcal C }$ can be interpreted as a distance between a given quantum state and the nearest separable state [64]. For these reasons, ${ \mathcal C }$ can be used to quantify the entanglement of two qubits instead of E_f, and is hereinafter called pairwise entanglement.

The eigenvalues of the matrix defined by equation (49) can be written as

$$\begin{eqnarray}{\tilde{\lambda }}_{1} & = & {\tilde{\lambda }}_{4}={u}_{+}^{X}{u}_{-}^{X},\\ {\tilde{\lambda }}_{\mathrm{2,3}} & = & {(\sqrt{{w}_{1}^{X}{w}_{2}^{X}}\pm | {z}^{X}| )}^{2},\end{eqnarray} \tag{ 51 }$$

where ${\tilde{\lambda }}_{1},\ldots ,{\tilde{\lambda }}_{4}$ are yet unsorted. After some straightforward steps, equation (48) leads to

$$\begin{eqnarray}&&{ \mathcal C }=2\max \{0,| {z}^{X}| -\sqrt{{u}_{+}^{X}{u}_{-}^{X}}\},\end{eqnarray} \tag{ 52 }$$

where the relevant correlation functions are given by equation (46) for X = c, or by equation (47) for X = s. The correlation functions w₁^X and w₂^X are absent in equation (52) as they contribute to $\sqrt{{\lambda }_{1}}-\sqrt{{\lambda }_{2}}-\sqrt{{\lambda }_{3}}-\sqrt{{\lambda }_{4}}$ if and only if $\sqrt{{\lambda }_{1}}-\sqrt{{\lambda }_{2}}-\sqrt{{\lambda }_{3}}-\sqrt{{\lambda }_{4}}\lt 0$.

Generalizations of ${ \mathcal C }$ for multifermionic states, appearing when larger subsystem A is under consideration, are also possible [65]. Such generalizations are, however, beyond the scope of this paper. We focus here on the pairwise entanglement, mainly because the form of equation (52) makes it possible to be determined with no significant computational costs once standard (spin or charge) correlation functions are calculated.

Our numerical results for the pairwise entanglement are presented in figures 6 and 7, constituting the central findings of this paper. The presentation is limited to the cases when sites i and j are the nearest neighbors in a linear chain (i.e., $j=i\pm 1\,$ mod $\,N$), as we have found that ${ \mathcal C }=0$ for more distant neighbors⁵ . Also, correlation functions in equation (52) are averaged over the reference site number i, to reduce the finite-precision effects.

Zoom In Zoom Out Reset image size

For ${N}_{\mathrm{el}}=N$, we have ${ \mathcal C }=0$ for the charge degrees of freedom (so-called charge concurrence) provided that $R/{a}_{0}\gtrsim 4$ (see top panel in figure 6). In such a range, the separability of a quantum state in the position representation, appearing if a bipartite subsystem (A) is selected in a fixed-spin sector, justifies the above-used notion of a fully-reconstructed Mott insulator for a finite N (see section 3 and [3, 12]). For smaller R, stronger charge fluctuations manifest itself via nonzero charge concurrence, providing an insight into the nature of a partly-localized quantum liquid. Remarkably, spin concurrence (see bottom panel in figure 6) is positive for any considered N and $R/{a}_{0}$, reaching the maximum in the crossover range of $4\lesssim R/{a}_{0}\lesssim 5$. For this reason, spin concurrence still can be considered as a prospective signature of the Mott transition.

For ${N}_{\mathrm{el}}=N/2$ the evolution of the pairwise entanglement with $R/{a}_{0}$ (see figure 7) is significantly different than for the ${N}_{\mathrm{el}}=N$ case. For each N, we have a nonzero charge concurrence for any considered $R/{a}_{0}$ (top panel in figure 7), whereas the spin concurrence (bottom panel) indicates a finite-system version of the Mott transition⁶ . In brief, for spin degrees of freedom ${ \mathcal C }=0$ for $R\gt {R}_{\star }(N)$, where the nodal value ${R}_{\star }(N)$ corresponds to $| {z}^{s}| -\sqrt{{u}_{+}^{s}{u}_{-}^{s}}=0$, and is determined numerically via least-squares fitting of a line to the actual datapoints for a given N (see the last column in table 2). For larger N, the values of ${R}_{\star }(N)$ are close to ${R}_{c}\approx 4.2\,{a}_{0}$, determined in section 3 via the finite-size scaling for ${\rm{\Delta }}{E}_{C}$.

Zoom In Zoom Out Reset image size

Additionally, we observe that the nearest-neighbor pairwise entanglement may indicate the decoupling of charge and spin degrees of freedom in correlated systems. This happens for $R/{a}_{0}\gtrsim 4$ in the ${N}_{\mathrm{el}}=N$ case (as well as for $R\gt {R}_{\star }(N)$ in the ${N}_{\mathrm{el}}=N/2$ case), where nonzero spin (charge) concurrence is accompanied by vanishing charge (spin) concurrence. In the above-mentioned ranges, quantum states of a bipartite subsystem (A) are entangled when spin, but separable when charge degrees of freedom are under consideration (${N}_{\mathrm{el}}=N$), or vice versa (${N}_{\mathrm{el}}=N/2$). Such a decoupling should be distinguished from charge–spin separation predicted to appear in 1D systems showing the Luttinger liquid phase [44–47]. Our numerical examples also show there is no one-to-one relation between the decoupling and the Mott transition, as the format can be identified either for the system showing the crossover behavior only or undergoing the Mott transition (in the large-N limit). Finally, the pairwise-entanglement analysis presented above allows one to identify noticeably different behavior of charge and spin correlation functions without referring to their asymptotic behavior, as previously proposed in the literature (see the first paper in [12, 44, 66]).