Ground-State Properties of the Hydrogen Chain: Dimerization, Insulator-to-Metal Transition, and Magnetic Phases

(Simons Collaboration on the Many-Electron Problem) IBM Almaden Research Center, San Jose, California 95120, USA Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA SISSA-International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy The Center for Advanced Quantum Studies and Department of Physics, Beijing Normal University, Beijing 100875, China Department of Physics and Astronomy, University of California, Irvine, California 92697-4575, USA Institute for Theoretical Physics, University of Amsterdam, Science Park 904 Postbus 94485, 1090 GL Amsterdam, The Netherlands Department of Chemistry, Wesleyan University, Middletown, Connecticut 06459, USA Center for Computational Quantum Physics, Flatiron Institute, New York, New York 10010, USA Department of Physics, Columbia University, New York, New York 10027, USA Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany Istituto per i Processi Chimico Fisici del CNR (IPCF-CNR), Via G. Moruzzi 1, 56124 Pisa, Italy Democritos Simulation Center CNR-IOM Istituto Officina dei Materiali, Via Bonomea 265, 34136 Trieste, Italy Department of Physics, College of William and Mary, Williamsburg, Virginia 23187-8795, USA


I. INTRODUCTION
The study of interacting many-electron systems is profoundly important to many areas of science, including condensed-matter physics and quantum chemistry. Materials properties result from a delicate interplay of competing factors, including atomic geometry and structure, quantum mechanical delocalization of electrons from atoms, the entanglement implied by quantum statistics, and electronelectron interaction. Capturing such effects accurately is essential for understanding properties, for predictive computations, and for the realization of materials genome or materials-by-design initiatives, but it requires a sufficiently accurate solution of the interacting electron problem.
Over the years, many approaches to the problem of solving the many-electron Schrödinger equation have been formulated. Density functional theory (DFT) obtains the ground-state energy of the many-body system in terms of the solution of an auxiliary one-electron problem plus a self-consistency condition. DFT-based methods have had an enormous impact on materials science and condensedmatter physics, but the approaches become less reliable, or even break down, in the presence of strong electronic correlation effects, including magnetic, structural, conductive, and superconductive phase transitions [1][2][3][4][5][6][7]. Alternative approaches in terms of model systems, such as the single-orbital Hubbard model, have contributed very significantly to our understanding of electronic correlation physics but involve dramatic simplifications, including the neglect of many of the terms in the physical Coulomb interaction and truncation to a small fixed orbital basis, whose consequences are not fully understood. The need to establish systematic approaches that are chemically realistic and fundamentally many-body is thus intensely investigated.
A linear chain of hydrogen atoms (N protons equispaced along a line, with N electrons) [8][9][10][11][12] embodies many central themes of modern condensed-matter physics while retaining a certain degree of computational tractability. It features a periodic atomic potential, incorporates realistic Coulomb interactions, requires the construction of accurate and compact basis sets, and yet maintains a connection with the fundamental Hubbard model, which has been a hallmark of the theory of interacting fermions.
Here, we deploy a combination of the most advanced and accurate many-body methods to investigate the groundstate phase diagram of the H chain. Our methods have been benchmarked in a previous collaborative effort to obtain the ground-state equation of state [12]. In this work, we have introduced methodological advances to compute properties beyond the total energy. There have also been several previous studies of this system [8][9][10][11]. However, they were restricted either to small basis sets or finite system sizes, which prevented a realistic description of the H chain, or by their accuracy and capabilities, which prevented reliable resolution of the delicate scales and discovery of all the phases. The synergistic application of complementary methods distinguishes the present work and allows us to make robust predictions with a multimessenger approach in this challenging problem.

II. METHODS
We consider a system of N protons at fixed, equally spaced positions along a line, with N electrons:Ĥ where ðr 1 ; …; r N Þ are the electron coordinates, and R a ¼ aRe z is the coordinate of the ath proton. This Hamiltonian has been studied in the finite basis or finite systems (or both) and has generated increasing interest [8,10,11,[13][14][15][16]. We use atomic units throughout, in which energies are measured in Hartrees (me 4 =ℏ 2 ) and lengths in units of the Bohr radius a B ¼ ℏ 2 =ðme 2 Þ. In the thermodynamic limit (TDL) of infinite system size at zero temperature, which is our primary focus, our system is characterized by only one parameter, R.
We access the ground-state properties using multiple first-principles many-body methods, including variational Monte Carlo (VMC), standard and sliced-basis density-matrix renormalization group (DMRG, sb-DMRG) [17][18][19][20], auxiliary-field (AFQMC) [21], and latticeregularized diffusion Monte Carlo (DMC) methods [22]. For reference, independent-electron calculations are also performed, including restricted and unrestricted Hartree-Fock (RHF and UHF, respectively) [23] and DFT [24]. While it is not practically possible with any single method to simultaneously converge all aspects of the electronicstructure calculations (such as the basis set limit, manyelectron correlation, and the thermodynamic limit) across different regimes of the ground-state phase diagram, we draw our conclusions based on the convergence of multiple approaches to a consistent physical picture.
In this work, we introduce several technical advances that enabled the calculations on the hydrogen chain presented here. Within AFQMC, we implement techniques for the study of metal-insulator transitions based on the modern theory of polarization and the calculation of structure factors in chemical simulations. A procedure is developed for generating basis sets for confined systems (as in the planes perpendicular to the chain) from downfolded virtual Kohn-Sham orbitals [25]. Within DMC, we implement a scheme for computing unbiased estimates of correlation functions and develop basis sets for the short-R regime, as well as a more efficient parametrization of the trial wave-function ansatz. Within DMRG, we introduce new basis-set techniques building on the idea of Gausslets [26], and we implement the computation of entanglement entropy and dimerization in electronic density. To approach the continuum limit for small values of R, we devise and extensively compare basis sets including plane waves, continuum coordinate space, and problem-specific Gaussians with diffuse and polarized functions, as well as the new basis sets mentioned above. To approach the TDL in our numerical computations, we study increasingly long chains with open boundary conditions (OBC) or supercells with periodic Born-von-Karman boundary conditions (PBC). Further descriptions of our methods and discussions of technical advances and issues can be found in the Supplemental Material (SM) [27].

III. RESULTS
In this section, we present our results as a function of R. We start by describing the large-R regime, where we find that the chain features quasi-long-range antiferromagnetic correlations, followed by dimerization in several properties, namely, electronic density, kinetic energy, and entanglement entropy as R is reduced. We then proceed to describe the short-R regime, distinguished by a metal-insulator transition occurring at a critical bond length R MIT , so that the chain is insulating for R > R MIT and metallic for R < R MIT . We then discuss the physical properties of the metallic phase and show that the MIT arises from a selfdoping mechanism.
A. Large-R regime

Antiferromagnetic correlations
We begin our study of the H-chain phase diagram at large proton-proton separation. Here, to a first approximation, the system is a collection of isolated H atoms, each with a single electron in the atomic 1s orbital. This system is very similar to the half-filled Hubbard model in the large coupling (U=t) limit. However, the weakly bound nature of H − and its very diffuse orbitals can create excitons with strong binding in the H chain.
At large R, the correlations in the H chain can be characterized in terms of a spin-1 2 Heisenberg chain. The Heisenberg chain is a critical system, with power-law decay of AFM spin-spin correlations, hŜ 0 ·Ŝ i i. In response to a local perturbation (such as a local magnetic field), one observes local ordering (such as local Néel order), which decays as a power law in the distance from the perturbation. To probe AFM correlations in the H chain, in Fig. 1 we present the quantity C i ¼ hn 0↑ni↓ i − hn 0↑ ihn i↓ i computed with DMRG in the minimal (STO-6G) basis, wheren iσ denotes the number of electrons occupying (orthogonalized) atomic orbital i with spin polarization σ along z. Here, C i oscillates with wavelength λ ¼ 2R, corresponding to the Néel vector of two-sublattice antiferromagnetism; the wave vector may also be thought of as twice the Fermi wave vector q ¼ ð2π=λÞ ¼ 2k 0 F of a paramagnetic 1D ideal Fermi gas of density ρ ¼ ð1=RÞ.
In Fig. 1, we further observe that the oscillations in C i decay with a power-law envelope, indicates weakening of AFM as the bond length is reduced. The power-law decay of the AFM correlation is consistent with quasi-long-range order in 1D. The fitted exponent of η ≃ 1.11ð1Þ, likely affected by finite-size effects, is slightly higher than the prediction from conformal field theory, which gives hS i ·S 0 i ∝ ð−1Þ i ffiffiffiffiffiffi lni p =i η , with η ¼ 1 for systems within the same universality class as the 1D Heisenberg chain [28,29].

Dimerization
We observe dimerization as R is reduced in the large-R regime. With PBC, dimerization can be probed by densitydensity correlations. With OBC, or with a local perturbation, density dimerization can be measured by the electronic density, for example, integrated along transverse slabs (i.e., over x and y for a finite δz), nðzÞ. The upper portion of Fig. 2(a) shows nðzÞ versus z, for a segment at the center of the chain under OBC, computed with AFQMC. The density has maxima at the proton sites and minima halfway between. We define the dimerization measure for an N-atom chain, Δ N , by the difference between the two adjacent local minima of nðzÞ in the center of the chain as illustrated in the lower portion of Fig. 2(a). When dimerization is present, we find that its amplitude increases as R is decreased.
In Fig. 2(b), we investigate the dimerization in the H chain at R ¼ 2.0a B , using AFQMC, sb-DMRG, DMC, and VMC. We find that, in the middle of the chain, dimerization decays with chain length as Δ N ∝ N −d , with exponents d ¼ 0.57ð3Þ, 0.563(4), 0.58 (11), and 0.90(26) from AFQMC, sb-DMRG, DMC, and VMC, respectively. The subtle power-law physics of 1D correlated systems is not easy to capture accurately with numerical methods. It is encouraging that a quantitative agreement is seen between our very different many-body methods. HF is qualitatively incorrect here, giving either no order or true long-range order; DFT results are similar.
It is illuminating to relate the dimerization order observed in the H chain to a kind of dimerization order that appears in the Heisenberg model, as seen in nearestneighbor spin correlations, hðŜ 0 ·Ŝ 1 ÞðŜ i ·Ŝ iþ1 Þi under PBC. In the case of OBC, the open ends act as a local perturbation, which couples to the dimerization, making odd-numbered bonds stronger than even ones, with an amplitude that decays as a power law away from the ends. Criticality can thus be observed by studying the strength of the dimerization in the center of a chain as a function of N. In the Heisenberg chain, the dimerization profile as a function of site index i is described by Δ N ðiÞ ∝ ½N sinðπi=NÞ −d , where d ¼ 1=2. This case is consistent with our observations for the H chain. In the half-filled Hubbard model, we observe dimerization in the correlation function of the nearest-neighbor hopping, i.e., the kinetic energy operator,ĥ i ¼ P σ ðâ † i;σâiþ1;σ þ H:c:Þ. As shown in the SM [27], the correlation function Dimerization manifests itself not only in the electron density, but also in other observables such as the kinetic energy and entanglement entropy. Results for a 10-atom H chain are shown in Fig. 2(c) from sb-DMRG, where the kinetic energy is measured as hopping (integrated in x and y) between slices, and the entanglement entropy is measured across planes that bisect a bond. At R ¼ 2.0 a B , where dimerization is consistently observed in larger systems from the electronic density, S vN ðzÞ shows five entangled dimers (bottom of lower panel). The entanglement exhibits a striking asymmetry on the two sides of each proton along the chain, with a strength discrepancy much larger than in the electron density.
B. Short-R regime

Insulator-to-metal transition
Correlated electron materials often exhibit metalinsulator transitions as parameters such as temperature, pressure, or crystal structure are varied [3,6]. From the perspective of the 1D one-band Hubbard model, which, as we have seen, captures the universal aspects of the physics at large R-and more generally, from the perspective of one-band models with second-order Umklapp processesno MIT should occur here [30]. With multiple methods and multiple probes, we show below conclusive evidence that a MIT occurs in the H chain, and we provide a characterization of the physical origin and properties of the transition.
The concept of macroscopic localization [11,31,32] provides a direct wave-function-based characterization of system properties. For periodic systems, one defines the complex polarization Z N ¼ hΨ N je ið2π=LÞ P iẑ i jΨ N i, where jΨ N i is the ground state of the N electrons in a supercell of size L ¼ NR along the chain direction. The electron localization length Λ ¼ ð ffiffiffiffi D p =2πρÞ is related to the complex polarization by D ¼ − lim N→∞ N log jZ N j 2 . In localized systems, Z ≡ lim N→∞ jZ N j ¼ 1, and Λ is finite; in metallic systems, Z ¼ 0, and Λ → ∞.
In Fig. 3(a), we establish the MIT by computing Z N with AFQMC, DMC, and VMC. (Additional support for a metal insulator transition is provided by energy gaps as included in the SM [27].) For small R, all methods give Z N equal to or statistically compatible with 0 across a wide range of system sizes N, indicating a metallic many-body ground state. For large R, all methods yield a nonzero Z N . The many-body methods point to a transition point located approximately at R MIT ∼ 1.70ð5Þ a B . While there is some uncertainty in the critical value, because of computational limitations, it is remarkable that two methods working with completely different basis sets and projection algorithms yield results in quantitative agreement.
In particular, the values of jZ N j in the insulating phase, which fall further from unity at smaller R, are sensitive to finite-size effects. At large N, we expect jZ N j ¼ 1 −gðξ=LÞ in the insulating phase, whereg is a scaling function and ξ the MIT correlation length. The smooth decrease of jZ N j as R MIT is approached suggests a secondorder transition (related to gap closure). We further discuss finite-size corrections and their implication on ξ in the SM [27].
Before exploring the origin of the MIT, we briefly discuss magnetic correlations in the metallic phase. The evolution of the magnetic moments within DFT is shown in Fig. 3(b), and the spin density at R ¼ 0.9 a B is plotted in Fig. 3(c) as a point of reference. To account for spin correlations, the DFT solution breaks translational symmetry to create antiferromagnetic domains of varying periods, often associated with very diffuse orbitals as seen in Fig. 3(c). Of course, DFT observations are independent-electron in nature (and somewhat functional dependent, see SM [27]). In the many-body solution, in particular, translational symmetry is restored, and two-body correlation functions are needed to probe magnetic correlations.

Origin of the MIT and properties of the metallic phase
It is theoretically established that the ground state of a one-band model with commensurate filling is insulating. Our calculations reveal that the MIT arises from a selfdoping mechanism in which the one-band picture breaks down. We illustrate the basic idea in Fig. 4(a), using a bandtheory-based cartoon of the electronic structure. The isolated H atom has multiple states, including the occupied 1s and excited states 2s, 2p, etc. At large R, the bandwidths are small compared to the energy gaps, and a one-band approximation is reasonable. As R is decreased, the bands broaden, band overlap occurs, and metallic behavior results. Quantitative calculations involve a correlation problem that requires treating interactions in a multiband situation, effects that are entirely absent in the standard Hubbard model.
We compute the spin structure factor, defined as SðqÞ ¼ ð1=NÞhΨjρ † qρq jΨi, whereρ q ¼ P iŜz;i e iq·r i is the Fourier transform of the spin density at q ¼ ð0; 0; qÞ, witĥ S z;i andr i ¼ ðx i ;ŷ i ;ẑ i Þ denoting the spin-z and position operators of electron i, respectively. The result is shown in  Fig. 4(b). At R ¼ 2.5, only one peak is seen at q ¼ 2k 0 F , signaling power-law AFM order as discussed earlier. To probe the nature of the metallic phase, we focus on a representative case of R ¼ 0.9 a B , away from the vicinity of the MIT transition. Note that SðqÞ is shown from two different QMC calculations, in supercells of N ¼ 48 atoms averaging over 11 twist angles, in addition to DMRG calculations in a supercell of N ¼ 24 atoms. In a metallic system, we expect peaks at q ¼ 2k F , where k F is one of the Fermi wave vectors of the system. Two cusps are seen in each result at locations q 1 and q 2 , in precise agreement among the different calculations. We interpret the larger wave vector as arising from the 2k F process in the 1s-dominated lower band. The position q 2 ¼ 2ð1 − xÞk 0 F then gives the doping x of this band. In a simple two-band picture, the lower wave vector is given by q 1 ¼ 2xk 0 F =g, where g gives the degeneracy of the "upper band," which is occupied. At R ¼ 0.9 a B , the locations of q 1 and q 2 satisfy gq 1 þ q 2 ¼ π=R in all the results, with g ¼ 2. This case is consistent with a doubly degenerate upper band (e.g., 2p x;y ), although multiple bands strongly contribute to the DMRG SðqÞ, as we discuss below.
Within DFT, one may characterize the metallic state by the nature of the occupied Kohn-Sham bands. Here, we use the Perdew-Burke-Ernzerhof (PBE) exchangecorrelation functional in DFT calculations. As R decreases, DFT-PBE solutions indicate a band of mainly atomic 1s origin with a relatively large (but not integer) occupancy and a Fermi level near k 0 F , and bands of mostly 2p x , 2p y character (i.e., π bands) becoming occupied, with Fermi levels near 0.
We probe the many-body metallic state by computing the one-body density matrix, which we then diagonalize. Each eigenvalue gives the occupancy of the corresponding eigenstate ("natural orbital," NO). The occupancies computed by two methods are shown in Fig. 5(a). For the N ¼ 24 supercell (red symbols), a sharp transition is seen in the orbital occupancy, reminiscent of a Fermi surface. We label the characters of each NO, which, similar to Kohn-Sham orbitals, can be analyzed via the single-particle basis or projection onto atomic orbitals. In the N ¼ 24 system, the last two eigenvalues before the sharp drop (11th and 12th points from left, shown as filled symbols) are twofold degenerate, consistent with the notion that the metallic phase corresponds to shifting electrons from 1s-like to 2p x;y -like orbitals. However, the many-body results display variations in the occupancy and the nature of the orbitals, as illustrated in Figs. 5(b) and 5(c). We find that the energy differences associated with different occupancies, such as 2s instead of or together with 2p, are very small, on the scale of sub-mHa/atom, and the precise energy ordering is affected by boundary conditions, supercell sizes, and basis sets, as well as the choice of the correlation algorithm. This low-energy scale is associated with the very diffuse nature of "upper band" NOs.
The occupation and nature of these NOs shed light onto magnetic correlations. Electrons associated with the "lower band" are bound to the chain axis, with a range similar to that of the atom. Note that SðqÞ does not show a peak at q 2 but only a kink, indicating short-range correlations from these electrons. The electrons in the "upper bands" are diffuse, with NOs of transverse size of about 5 a B [see Figs. 5(b) and 5(c)]. They display a tendency for magnetic ordering driven by exchange, as indicated by the kink at q 1 . Their average interparticle separation along the chain, ξ ∼ π=ð2q 1 Þ, is large (and becomes larger as R MIT is approached), with only a small number of electrons in a large supercell. These features make it challenging to determine the precise behavior of the correlation over a sufficiently long range, which is further magnified by the sensitivity of the nature of the upper bands discussed earlier. This picture suggests a fascinating model of a quasi-onedimensional electron wire for probing both magnetic and charge correlations under the effect of long-range Coulomb interactions [33].

IV. CONCLUSIONS AND OUTLOOK
While much has been understood about model systems, most notably the Hubbard model, understanding correlated electron physics and reliably performing predictive computations in real materials remain central challenges in condensed-matter physics and materials science. Significant progress has been made from several fronts, for example, with the combination of DFT and the GW and Bethe-Salpeter equation formalisms [34][35][36], approaches based on dynamical mean field theory [5,[37][38][39], as well as several of the methods included here.
In this paper, we have focused on pushing state-of-the-art many-body ground-state methods towards a comprehensive theory of correlated materials. Individually, advances are made within each method to address the technical challenges for systematically accurate calculations of correlated systems in the continuum and thermodynamic limit. Collectively, with a combination of multiple methods, we are able to achieve a high level of robustness in our characterization and predictions of the ground-state properties.
The hydrogen chain provides a good compromise between physical and chemical realism and tractability. We studied the balance and competition between bandstructure, electron-interaction, and multiorbital physics as a function of a single parameter, the interproton distance R. At large R, the system exhibits a Mott insulating phase with quasi-long-range AFM order, as well as dimerization effects visible in the density, kinetic energy, and entanglement spectrum. Upon decreasing R, it transitions via a selfdoping mechanism into a metallic phase, which exhibits magnetic correlations. At short bond length, the magnetic correlation is modulated by an incommensurate wave vector. Near the MIT, the occupancy of the "higher band" is very low, and the size of the magnetic domain grows, potentially allowing for the possibility to develop global polarization. We find that the higher natural orbitals (occupied in the metallic phase and important as virtual excitations in the large-R insulating phase) have an extremely diffuse and occupation-dependent structure. Thus, despite its apparent simplicity, the H chain embodies a surprisingly rich set of the major themes in contemporary quantum materials physics, including metal-insulator transitions, magnetism, and intrasite and intersite electroncorrelation physics.
It is interesting to place our results in historical context. The study of the MIT was pioneered by Nevill Mott. In 1949, Mott noted [1] that many materials were insulating when band theory considerations indicated that they should be metallic. He argued that one should view the physics in terms of localized atomiclike orbitals rather than plane waves, and he observed that, in this picture, conductivity required moving an electron from one site to another, creating a particle-hole pair. He argued that this process could be blocked if the charging energy on a given lattice site was large enough to prevent an additional electron from moving onto it. This physics is famously instantiated in the extensively studied Hubbard model [40], and our results in the large-R regime are consistent with it. In subsequent work [2,3], Mott noted that, in the insulating phase, the unscreened Coulomb interaction would bind the particle and hole to a neutral object, which would not conduct; however, in a metal, the self-consistent screening of the Coulomb interaction would lead to unbound pairs, so the MIT occurs via a first-order process. Our results demonstrate a different route to the MIT: The H chain becomes metallic via a self-doping process, which is (within our resolution) second order. The carriers in the second band are both very delocalized and at very different k points from the carriers in the dominantly occupied band, so excitonic effects are either absent or unobservably weak. Our work underscores the importance of the interplay of correlation physics with the real materials effects of band crossing and electronic structure.
Experimentally, the physics we have described could potentially be realized in ultracold atoms in optical lattices, arrays of Rydberg atoms [41], or perhaps with carbon nanotubes, in which carbon chains have been isolated [42]. Experimental realization would provide an in silico insulator-to-metal transition, which would be a significant contribution to the intensely investigated topic of metallicity in hydrogen [43][44][45][46]. We also hope that our study will stimulate experimental research on the characterization of hydrogen at extreme density [46][47][48][49][50], in connection with the theoretical description of astrophysical bodies [51] and the discovery and characterization of exotic magnetic phases in low-dimensional condensed-matter systems [46,52].
More generally, our study establishes the H chain as an important model system and demonstrates the power of the synergistic application of independent numerical methodologies [12,53,54]. This mode of attack has the potential to address other challenging problems in condensed-matter physics, including the systematic investigation of real materials with immediate technological applications by first-principles many-body methods.