Pairwise entanglement and the Mott transition for correlated electrons in nanochains

Pairwise entanglement, calculated separately for charge and spin degrees of freedom, is proposed as a ground-state signature of the Mott transition in correlated nanoscopic systems. Utilizing the exact diagonalization—ab initio, for chains containing N ≤ 16 hydrogenic-like atoms (at the half filling), we find that the vanishing of the nearest-neighbor charge concurrence indicates the crossover from a partly-localized quantum liquid to the Mott insulator. Spin concurrence remains nonzero at the insulating phase, showing that the decopling of spin and charge degrees of freedom may manifest itself by wavefunctions entangled in spin, but separable in charge coordinates. At the quarter filling, the analysis for N ≤ 20 shows that spin concurrence vanishes immediately when the charge-energy gap obtained from the scaling with 1 / N → 0 vanishes, constituting a finite-system version of the Mott transition. Analytic derivations of the formulas expressing either charge or spin concurrence in terms of ground-state correlation functions are also provided.


Introduction
In 90 years after the seminal work by Schrödinger [1] analytical properties of quantum-mechanical wavefunctions describing hydrogenic-like atoms (or ions) continue to surprise. For instance, recent derivation of the Wallis formula for π by Friedmann and Hagen [2] refers to the variational principle and the correspondence principle, but also employs acounterintuitive fact that some variational approaches may truncate particular excited states of atoms with the same or better accuracy than the ground state (GS). This feature has direct analogs in many-electron hydrogenic systems and motivated the proposal of the exact diagonalization-ab initio method (EDABI) [3], which were recently used to discuss the influence of electron correlations on the metallization of solid hydrogen in arigorous manner [4,5]. Also, as amethod putting equal footing on single-and multiparticle aspects of the correlated quantum states, EDABI seems to be apromising candidate for the theoretical tool capable of giving abetter insight into the superconductivity mechanism in sulfur hydrides [6,7] or into the magnetic moment formation in functionalized graphenes [8].
When considering ageneric second-quantized Hamiltonian with both spin and charge degrees of freedom, afew numerical techniques can be regarded as exact ones; i.e., giving the desired correlation functions within the accuracy limited (in principle) by the machine precision only. These includes: exact diagonalization (ED) for relatively small systems [9], density matrix renormalization group (DMRG) for low-dimensional systems [10], and quantum Monte Carlo (QMC) for non-frustrated systems [11]. Aseparate class is outlined by variational approaches designed to treat one-dimensional (1D) atomic chains in the insulating phase [12]. Even though each of these methods provides us with detailed information about the closed-system GS, it is usually achallenging task to determine whether GS is insulating, metallic, or of amore complex nature [13] 1 . This is primarily because standard GS signatures of the metal-insulator transition, such as avanishing charge gap [14], are absent at any finite system size (quantified by the number of lattice sites N), but may appear only after the so-called finite-size scaling (with N 1 0  ), aprocedure introducing systematic errors difficult to estimate in some cases. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
For this reason, specially-designed GS correlation functions, not only signaling the Mott transition in the thermodynamic limit, but also showing fast convergence (with N 1 0  ) for finite systems close to the metalinsulator boundary, may constitute avaluable complement of the existing numerical techniques.
Sacramento et al [44] used entanglement measures to complement the long-lasting discussion of decoupling of charge and spin degrees of freedom in 1D Hubbard model [45] and its relation to the charge-spin separation phenomenon in Luttinger liquids [46,47]. In particular, it is pointed out in [44] that charge-charge fluctuations, when quantified by properly defined correlation functions, show substantially different asymptotic behavior than spin-spin fluctuations in both the metallic and the insulating phases.
Here we focus on linear chains containing up to N=20 equally-spaced hydrogenic-like atoms, each one containing asingle valence orbital (see figure 1). Although such chains are not directly observable due to the well-known Peierls instability, they have been intensively studied to benchmark different ab initio approaches to strongly-correlated systems [48]. In this paper, two distinct physical regimes are considered: at the half filling the system shows a crossover behavior, with the increasing interatomic distance R, from apartly localized quantum liquid to the Mott insulating phase [3,49]. At the quarter filling, the Mott transition is reconstructed with . (The above-mentioned findings are consistent with more recent results of QMC simulations for Hubbard chains with long-range interactions [50].) It is worth to mention that the two regimes are also significantly different from amore fundamental point of view: in the strong-coupling (i.e., large R) limit the charge fluctuations are suppressed and the half-filled chain can be effectively described as 1D Heisenberg antifferomagnet, for which long-range order is absent due to the Mermin-Wagner theorem [51]. For the quarter-filling similar reasoning cannot be applied in the presence of long-range Coulomb interactions, and the charge-density wave phase is predicted to appear [3,[52][53][54]. In other words, we consider the two complementary model cases allowing one to test entanglement-based phase-transition indicators. We further discuss the concurrence [22], defined separately for charge and spin degrees of freedom, and employ it to recognize the decoupling of charge and spin degrees of freedom accompanying afinite-system version of the Mott transition.
The paper is organized as follows. In section 2 we present the system Hamiltonian and summarize the EDABI method. In section 3 the findings of [3], concerning the finite-size scaling of charge-energy gap, are revisited after taking larger systems under consideration. In section 4 we determine the local and pairwise entanglement, the latter for charge and spin degrees of freedom, and discuss their evolution with the interatomic distance and the system size. The conclusions are given in section 5. A system studied numerically (schematic). Hydrogenic-like atoms, containing 1s-type orbitals with their radii 1 a -, are arranged linearly with the interatomic distance R. The Hamiltonian parameters, including the atomic energy a  , the intraatomic Coulomb integral U, the hopping integrals t ij , and interatomic Coulomb integrals K ij , are determined for each value of α and R (see section 2 for the details).

The exact diagonalization-ab initio method (EDABI)
The EDABI method, together with its application to the system of figure 1, have been presented in several works [3][4][5]49]. Here we briefly recall the main findings, to which we refer in the remaining parts of the paper.

The Hamiltonian for alinear chain
The analysis starts from the second-quantized Hamiltonian, which can be written in the form where 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 t ij ºif i and j are nearest neighbors, otherwise t 0 ij » ), U is the intrasite Coulomb repulsion, K ij is the intersite Coulomb repulsion, and V e R ) expresses the repulsion of infinite-mass ions. The single-particle and interaction parameters of the Hamiltonian R ,  â ( )(1), also marked in figure 1, can be defined as follows 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 where the Slater 1s function r r R exp , with α being the inverse orbital size, here taken as variational parameter, and R i being the position of ith ion. The coefficients ij b in equation (4) can be uniquely defined by imposing that

The range of Coulomb interactions in the Hamiltonian
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 condensedmatter 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] 2 , 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 20  atoms.

Single-particle basis optimization
Next, each Slater function is approximated as follows where p is the number of Gaussian functions truncating j y , and B , q q G { }, q = 1,K, p, are adjustable parameters chosen to minimize the atomic energy for 1 a = and agiven value of p. Here we set p=3, for which the deviation from the exact energy for s 1 function is lower that 1% 3 . Subsequently, the parameters a  , t ij , U, and K ij are calculated from equations (2) and (3) as functions of α and R. The Hamiltonian R ,  â ( )(1) is diagonalized numerically in the Fock space, using the Lanczos algorithm, and the orbital size is optimized to find min a a = corresponding to the minimal GS energy E G (N) for each R. Defining the efective atomic energy as one finds that min a and other parameters converges rapidly with N. This observation allows us to speed up computations by using the values obtained for smaller systems (N 10  ) to perform extrapolation with

Finite-size scaling for the charge-energy gap
Astandard numerical approach [14,50], allowing one to determine whether acorrelated system described by the Hamiltonian such as given by equation (1) is metallic or insulating in the GS, involves calculating the chargeenergy gap according to where N el is the number of electrons. Next, the limit of . In particular, the largest considered subspace dimension corresponds to N=20, N 1 11 el + = , and is equal to 600 935 040. Moreover, we impose the periodic (t t i j N ij , table 2).  (6); the remaining symbols are t t ij º -    , the finite-size scaling results can be rationalized with the empirical function with the best-fitted parameters (numbers in parentheses are standard deviations for the last digit) For asake of conciseness, our discussion of finite-size scaling estimates of the Mott transition is limited to the charge-energy gap. Amore detailed analysis, presented in [3,49,55] and involving calculations of the electron momentum distribution, the Drude weight, as well as the so-called modern theory of polarization [63], justify the transition appearance in the N N 2 el = case, which coincides with the results for related parametrized model studies [50,52,53,57]. In the N N el = case, the situation is of aslightly more complex nature: apart from anonzero gap for any R, following from the finite-size scaling, several GS and dynamical characteristics (in particular-the Drude weight, see [3]) exhibit, for any finite N, acrossover behavior between apartly localized quantum liquid, appearing for small R, and afully-reconstructed Mott insulator, typically appearing for R a 4 0  . These findings coincides with more recent variational Monte-Carlo studies of hydrogenic chains (see [12]), and can be attributed to the increasing (with growing R) role of electron correlations in either the exact or variational GS wavefunction.
It is worth to mention here that some surprisingly efficient variational methods, such as presented in [12], become significantly less efficient when treating open-shell configurations in attempt to calculate E C D from equation (7), or generally when utilized away from the half filling (N N el ¹ ). Also for this reason, GS correlation functions signaling afinite-system version of the Mott transition for N small enough to be treated with some more flexible numerical techniques, are desired.

Reduced density matrix and quantum entanglement
In this section we derive explicit representations of the reduced density operator Tr where Tr A stands for tracing over all degrees of freedom except from these characterizing a selected subsystem A) relevant when discussing the local entanglement, the pairwise entanglement for charge degrees of freedom, and the pairwise entanglement for spin degrees of freedom. The derivations remain valid for ageneric spin-1/2 fermionic system at any pure state Yñ | , which can be represented assuming four basis states for each site j, namely: n n 0 , , ,

General considerations
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 Yñ | can be represented as follows where Y ab denotes acomplex probability amplitude corresponding to abasis state of the full system a b ñ Ä ñ | | , while añ {| }and bñ {| }are complete basis sets for the subsystems A and B (respectively). The reduced density operator A r is defined as: where the last equality follows from equation (11) and  Y ab denotes the complex conjugate of Y ab . Subsequently, the matrix elements of A r are given by We define now projection operators P α and P β , associated with the basis states añ | and bñ | , via Useful properties, directly following from equation (14) are where equation (16) In turn, the reduced density matrix, as given by the rightmost equality in equation (13), can be rewritten as

Entanglement entropy
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 añ {| }is simply given by equation (10), one can write down the corresponding projection operators P n n n n n n nn

and the transfer operators
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) Equation (27) corresponds to the reduced density operator of the well-known form 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 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 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 Remarkably, the strong-correlations limit corresponds to the minimal E v , whereas the free-electrons limit corresponds to the maximal E v (for any n 0 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 singleparticle 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 n Var 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 N el = and N N 2 el = (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 1.5 0 = , 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 , whereas for N N el = we still have E 1 v > in this range. However, the smooth evolution of E v with R is observed for both N N el = and N N 2 el = , making it difficult to consider the local entanglement as an estimate of the Mott transition. Pairwise entanglement is discussed next as the other candidate.

The fermionic concurrence
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.
and P n n n n n n n n n n n n n n n n 1 1 , 1 1 , with the upper indices (c s , ) referring to charge and spin qubits (respectively). Subsequently, the transfer operators are given by for spin qubits.
We use now the concurrence , as aquantitative measure of quantum entanglement in the two-qubit subsystem A. The closed-form expression was derived by Wootters [22] and reads For N N el = , we have 0  = for the charge degrees of freedom (so-called charge concurrence) provided that R a 4 0  (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 R a 4 5 0   . For this reason, spin concurrence still can be considered as aprospective signature of the Mott transition. For N N 2 el = the evolution of the pairwise entanglement with R a 0 (see figure 7) is significantly different than for the N N el = 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 6 . In brief, for spin degrees of freedom 0  = for R R N 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 4 0  in the N N el = case (as well as for R R N  > ( )in the N N 2 el = 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 N el = ), or vice versa (N N 2 el = ). Such adecoupling should be distinguished from charge-spin separation predicted to appear in 1D systems showing the Luttinger liquid phase [44][45][46][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]).

Conclusions
In this study we have demonstrated that pairwise entanglement, quantified by the fermionic concurrence determined separately for charge and spin degrees of freedom, can serve as aconvenient indicator for afinitesystem version of the Mott transition. In particular, standard finite-size scaling estimates of the Mott transition for linear chains of hydrogenic-like atoms are revisited utilizing the EDABI at the half and the quarter electronic filling. In the latter case, we find that not merely the charge gap indicates the transition for N 1 0  (with N being the number of atoms) at the interatomic distance R a 4.2 c 0 » (where a 0 denotes the Bohr radius), but also the spin concurrence vanishes for afinite N at R N  ( )taking the values relatively close to R c . Charge concurrence remains nonzero in both the metallic and the insulating phases, signaling adecoupling between charge and spin degrees of freedom.
At the half filling the Mott transition is not observed. Instead, the crossover from apartly localized quantum liquid to afully-reconstructed Mott insulator occurs. The charge concurrence vanishes at the crossover point, where the spin concurrence shows abroad peak; the latter remains nonzero in the entire parameter range.
It is worth to stress here that calculations of the fermionic concurrence, employed in this paper, generate essentially no extra computational costs, as the concurrence is determined solely via pairwise ground-state correlation functions for charge or spin degrees of freedom, without referring to the dynamical properties such as the optical or dc conductivity. The analytic relations between entanglement and the correlation functions (see equation (52)) are also derived in this paper.
Astriking feature is that the analysis holds true regardless how the exact (or approximated) GS is obtained, making it possible to employ the entanglement-based signatures of the Mott transition when discussing ageneric correlated-electron system. What is more, equation (52) also applies when correlation functions on its right-hand side are determined using more powerful numerical techniques, such as DMRG or QMC. For these reasons, we believe the approach we propose may shed new light on various open problems in the field, such as the long-standing metallization of solid hydrogen [4], or the recently-raised phase diagram of the repulsive Hubbard model on ahoneycomb lattice ( [13], see footnote 1).