Carbon nanotubes as excitonic insulators

Fifty years ago Walter Kohn speculated that a zero-gap semiconductor might be unstable against the spontaneous generation of excitons–electron–hole pairs bound together by Coulomb attraction. The reconstructed ground state would then open a gap breaking the symmetry of the underlying lattice, a genuine consequence of electronic correlations. Here we show that this excitonic insulator is realized in zero-gap carbon nanotubes by performing first-principles calculations through many-body perturbation theory as well as quantum Monte Carlo. The excitonic order modulates the charge between the two carbon sublattices opening an experimentally observable gap, which scales as the inverse of the tube radius and weakly depends on the axial magnetic field. Our findings call into question the Luttinger liquid paradigm for nanotubes and provide tests to experimentally discriminate between excitonic and Mott insulators.

detailed in Supplementary Note 7. The envelope F τ η is a pseudospinor with respect to valley and sublattice indices, F ≡ (F KA , F KB , F K A , F K B ) T . In the valley-sublattice product space, F obeys the Dirac equation of graphene: γ σ x ⊗ 1 τkx + σ y ⊗ τ zky F (r) = ε F (r). ( Here σ x and σ y are 2 × 2 Pauli matrices acting on the sublattice pseudospin, τ z and the 2 × 2 identity matrix 1 τ act on the valley pseudospin,k x = −i∂/∂x is is the wave vector operator along the circumference direction x andk y = −i∂/∂y acts on the tube axis coordinate y, γ is graphene's band parameter, and ε is the single-particle energy. Furthermore, F obeys the boundary condition along the tube circumference: where L is the chiral vector in the circumference direction of the tube and |L| = L = 2πR is the circumference. A magnetic field may or may not be applied along the tube axis, with ϕ = φ/φ 0 being the ratio of the magnetic flux φ piercing the tube cross section to the magnetic flux quantum φ 0 = ch/e. Supplementary Eq. (2) depends on the reference frame. Note that in our effective-mass treatment the x and y directions are parallel to the circumference and axis of the tube, respectively, as shown in Supplementary Fig. 1a, whereas in the main text as well as in the first-principles treatment the z axis is parallel to the tube.
The energy bands are specified by the valley index τ , the valence index α = c, v denoting either the conduction (α = c) or the valence band (α = v), and the wave vector k in the axis direction. The wave functions in K and K valleys are respectively F ≡ (F K αk (r), 0) T and (0, F K αk (r)) T , with F τ αk (r) ≡ (F τ A αk , F τ B αk ) T being a plane-wave pseudospinor in the sublattice space, where A is the tube length and the wave function ξ τ αk (x) for the motion along the circumference direction is The constant pseudospinor F τ αk is a unit vector with a k-dependent phase between the two sublattice components, where and s α = ±1 for conduction and valence bands, respectively. In Supplementary Eqs. (5) and (7) the transverse wave vector k ⊥ is proportional to the magnetic flux ϕ, In each valley, the energy is where the origin of the k axis is located at the Dirac point K (K ). K point, is obtained by specular reflection. The DFT / GW location of the Dirac point is K = 0.289 (2π/a), whereas the effective-mass estimate is K = 1/3 (2π/a) (the discrepancy between DFT and tight-binding predictions is well documented in the literature 4,5 ). As seen in Supplementary   Fig. 2a, the GW bands are approximately linear in an energy range of at least ± 0.4 eV around the Dirac point, which validates the effective-mass model at low energy.
Note that, in the absence of the magnetic field, electron states have a well defined chirality 6-8 , which is one of the two projections, C, of the sublattice pseudospin onto the momentum direction, expressed as the eigenvalues C = ±1 of the operator σ y ⊗ τ z . The chirality index is highlighted by red (C = +1) and black (C = −1) colour in Supplementary Fig. 2a.

Supplementary Note 2
Electron-electron interaction: Effective-mass vs first-principles description Within the effective-mass framework, the Coulomb interaction v between two electrons on the carbon nanotube cylindrical surface located at r ≡ (x, y) and r ≡ (x , y ), respectively, is 2 where κ r is a static dielectric constant that takes into account polarization effects due to the electrons not included in the effective-mass description plus the contribution of the dielectric back-ground. The interaction matrix element between single-particle states is 9,10 V (τ,α,k+q),(τ ,β ,k );(τ ,α ,k +q)(τ,β,k) where the one-dimensional effective interaction resolved in momentum space, is modulated by a form factor given by overlap terms between sublattice pseudospinors [I 0 (z) and K 0 (z) are the modified Bessel functions of the first and second kind, respectively 11 ]. The effect of screening due to the polarization of those electrons that are treated within the effective-mass approximation is considered by replacing v(q) with in the matrix element (11), where ε(q) is the static dielectric function (to discriminate between screened and unscreened matrix elements we use respectively 'w' and 'v' letters throughout the  Fig. 3 Effect of chirality on Coulomb interaction matrix elements. a Energy bands and chiralities of electron states in armchair carbon nanotubes in the absence of the magnetic field.
Solid and dashed lines highlight chirality C = ±1, respectively. b Allowed scattering processes induced by long-range Coulomb interaction. The indices τ = K, K and α = c, v label valleys and bands, respectively. The chirality is conserved at each vertex of diagrams.
Effect of chiral symmetry. The chiralities of electron states, which is illustrated in Supplementary Fig. 3a (solid and dashed lines label C = +1 and C = −1, respectively), signficantly affects Coulomb interaction matrix elements. This occurs through the form factors of the type F † ·F appearing in Supplementary Eq. (11), which are overlap terms between sublattice pseudospinors.
As apparent from their analytical structure, the chiral symmetry of the states is conserved at each vertex of Coulomb scattering diagrams (see Supplementary Fig. 3b), hence initial and final states scattered within the same band must have the same momentum direction. This significantly affects the Bethe-Salpeter equation for excitons, as we show below. We are especially interested in the dominant long-range Coulomb matrix element 12 that binds electrons and holes: This term scatters electron-hole pairs from the initial pair state (c, k)(v, k) to the final state (c, k + q)(v, k + q) within the same valley τ . Throughout this Supplementary Information we use the tilde symbol for quantities whose dimension is an energy multiplied by a length, like V =Ṽ /A.
In the first instance we neglect screening, since for low momentum transfer, q → 0, polarization is suppressed hence ε(q) → 1. In this limit Coulomb interaction diverges logarithmically, but this is harmless to the Bethe-Salpeter equation, since v(q) occurs only in the kernel of the scattering term, hence it is integrated over q for macroscopic lengths A, which removes the divergence. Note that, throughout this Supplementary  Hence, the regularized matrix element, integrated over the mesh, is and N = 900. The two matrix elements agree almost quantitatively, as they both exhibit: (i) zero or very small values in the second and fourth quadrants, i.e., k > K and k < K or k < K and k > K (ii) a logarithmic spike on the domain diagonal, i.e., k → k. This behavior has a simple interpretation in terms of exciton scattering, as an electron-hole pair with zero center-of-mass momentum, (c, k)(v, k), has a well-defined chirality with respect to the noninteracting ground state, i.e., ∆C = +2 = 1 − (−1) for k > K (∆C = −2 for k < K). The chirality of the e-h pair is conserved during Coulomb scattering, i.e., as the pair changes its relative momentum from Effect of electronic polarization. In order to appreciate the minor differences between V (k, k ) and W DFT (k, k ) it is convenient to compare the cuts of the maps of Supplementary Fig. 4 along a line k = k 0 , as shown in Supplementary Fig. 5 for k 0 = 0.289(2π)/a (panel a) and 0.28(2π)/a (panel b), respectively. For small momentum transfer, q = k − k 0 ≈ 0, V (k, k 0 ) (squares) exhibits a sharper spike than W DFT (k, k 0 ) (filled circles). This is an effect of the regularization of the singularitity occurring in the DFT approach, as in the first-principles calculation the tube is actually three-dimensional. As |q| increases, V is systematically blushifted with respect to W DFT since it does not take into account the effect of the RPA polarization, Π(q), which acquires a finite value. Within the effective-mass approximation, Π(q) enters the dressed matrix element W through the dielectric function 13 , Here we use the simple ansatz as this choice makes the dressed Coulomb interaction scale like the three-dimensional bare Coulomb potential for large q (i.e., at short distances), W ∼ 1/q 2 . In Supplementary Fig. 5a, b the dressed matrix element W [empty circles, A ansatz = 50/(πγ), γ/a = 1.783 eV] quantitatively agrees with its ab initio counterpart, W DFT (filled circles), in the whole range of k in which electrons are massless (cf. Supplementary Fig. 2). Note that for k > K = 0.289(2π)/a the effective-mass potentials are exactly zero whereas W DFT shows some numerical noise.
Effect of the magnetic field. The magnetic field along the tube axis adds an Aharonov-Bohm phase to the transverse momentum, k ⊥ . This breaks the chiral symmetry C of single-particle

Supplementary Note 3
Effective mass: Bethe-Salpeter equation In this Note we detail the calculation of low-lying excitons of armchair carbon nanotubes, |u , within the effective mass theory. The analysis of the first-principles exciton wave function for the (3,3) tube shows that the lowest conduction and highest valence bands contribute more than 99.98% to the spectral weight of excitons. Therefore, according to conventional taxonomy, these excitons are of the M 00 type. Within the effective-mass approximation, |u is written as where |0 is the noninteracting ground state with all valence states filled and conduction states empty, and the operatorĉ τ + k,σ (v τ + k,σ ) creates an electron in the conduction (valence) band labeled by wave vector k, spin σ, valley τ . The exciton |u is a coherent superposition of electron-hole pairs having zero center-of-mass momentum and amplitude ψ τ (k). The latter may be regarded as the exciton wave function in k space. The 2 × 2 spin matrix χ σσ is the identity for singlet excitons, χ = 1 s , whereas for triplet excitons χ = σ s ·n, where n is the arbitrary direction of the spin polarization (|n| = 1) and σ s is a vector made of the three Pauli matrices. Throughout this work we ignore the small Zeeman term coupling the magnetic field with electron spin, hence triplet excitons exhibit three-fold degeneracy. Here we use the same notation, |u , for both singlet and triplet excitons, as its meaning is clear from the context.

The Bethe-Salpeter equation for the triplet exciton is
The diagonal term E eh (τ, k) is the energy cost to create a free electron-hole pair (τ, c, k)(τ, v, k), including the sum of self-energy corrections to electron and hole energies, Σ τ (k), which may be evaluated e.g. within the GW approximation. This self-energy, which describes the dressing of electrons by means of the interaction with the other electrons present in the tube, is responsible for the small asymmetry of the Dirac cone close to K, as shown by the GW dispersion of Supple- Since this asymmetry appears already at the DFT level of theory and is similar to the one predicted for the Dirac cones of graphene 14 , it necessarily originates from mean-field electron-electron interaction and it does not depend on R. We take into account the effect of Σ τ (k) onto E eh (τ, k) by explicitly considering different velocities (slopes of the linear dispersions) for respectively left-and right-moving fermions, according to: We infer the actual values of γ and slope mismatch parameter α sl from the linear fit to the firstprinciples GW dispersion (in Supplementary Fig. 2b the solid lines are the fits and the dots the GW data), which provides γ = 5.449 eV·Å and α sl = 0.05929.
where Ω 0 = ( √ 3/2)a 2 is the area of graphene unit cell and w 2 > 0 is the characteristic energy associated with short-range Coulomb interaction. We reasonably reproduce first-principles results taking w 2 = 2.6 eV-this would be a plane located at 0.24 meV in Supplementary Fig. 8. This estimate is not far from Ando's prediction w 2 = 4 eV. Note that the previous theory proposed by one of us 3 relies on the scenario W KK > W , which is ruled out by the present study.
The Bethe-Salpeter equation for the singlet exciton is obtained from Supplementary Eq. (22) by simply adding to the kernel the bare exchange term where w 1 > 0 is a characteristic exchange energy 3,9 . From first-principles results we estimate w 1 = 4.33 eV, whose magnitude is again comparable to that predicted by Ando 9

. Supplementary
Eq. (22), with or without the exchange term, is solved numerically by means of standard linear algebra routines.
Minimal Bethe-Salpeter equation. The minimal Bethe-Salpeter equation illustrated in the main text includes only one valley (with α sl = 0) and long-range Coulomb interaction. Within the effective-mass approximation, the Dirac cone indefinitely extends in momentum space, hence one has to introduce a cutoff onto allowed momenta, |k| ≤ k c . Supplementary Fig. 9a shows the convergence of the lowest-exciton energy, ε u , as a function of k c . Reassuringly, ε u smoothly converges well within the range in which GW bands are linear. This is especially true for the screened interaction W (black circles), whereas the convergence is slower for the unscreened interaction V (red circles), as it is obvious since W (q) dies faster with increasing q. This behavior implies that the energy scale associated with the exciton is intrinsic to the tube and unrelated to the cutoff, as we further discuss below.
In the reported calculations we took k c = 0.05(2π)/a as a good compromise between accuracy and computational burden (we expect that the maximum absolute error on ε u is less than 0.1 meV). This corresponds to an energy cutoff of 1.4 eV for e-h pair excitations. Whereas for these calculations, as well as for the data of Supplementary Fig. 9a, the mesh ∆k in momentum space is Supplementary Fig. 9b shows the convergence of ε u as a function of the mesh, ∆k. Interestingly, ε u smoothly decreases with ∆k only for a very fine mesh, whereas for larger values of ∆k the energy exhibits a non-monotonic behaviour. This is a consequence of the logarithmic spike of the Colulomb potential at vanishing momentum, which requires a very fine mesh to be dealt with accurately.
that the wave vector k appearing as an argument of the interaction is always multiplied by R. This is important for the exciton scaling behaviour.
Neglecting the small corrections to the exciton binding energy due to intervalley scattering (w 2 = 0) and cone asymmetry (α sl = 0), the dimensionless Bethe-Salpeter equation for armchair tubes in the absence of a magnetic flux becomes where α graph = e 2 /(κ r γ) is graphene fine-structure constant, the scaled exciton wave function must satisfy the scale invariant normalization requirement, τ dκ |ξ τ (κ)| 2 = 1, and the dielectric function entering Ω takes the dimensionless form The only scaling length leaving Supplementary Eqs. (27) and (28) with E 0 being calculated once for all for the (3,3) tube radius, R = 2Å. The same conclusion holds for finite cone asymmetry α sl and dimensionless magnetic flux ϕ. Note that, for a fixed value of ϕ, the possible values of the magnetic field B scale like 1/R 2 .
The above demonstration relies on the assumption that the parameters κ r and A ansatz , which control the screening behavior of the carbon nanotube, do not depend significantly on R. On the other hand, one might expect to recover, for large R, the screening properties of graphene. This in turn would imply that κ r would tend to smaller values and hence ε u would decay slower than 1/R.
The first-principles investigation of this issue is left to future work.
Supplementary Note 4 Self-consistent mean-field theory of the excitonic insulator The ground-state wave function of the excitonic insulator, |Ψ EI , exhibits a BCS-like form, where η is the arbitrary phase of the condensate, the e-h pairsĉ τ + The symbolW τ τ (k, k ) in Supplementary Eq.

(32) is a shorthand for both intra and intervalley
Coulomb interaction matrix elements. For the spin-triplet EI [χ σσ = (σ s ·n) σσ ], which is the absolute ground state, one has, for τ = τ , the long-range intravalley term,W τ τ (k, k ) =W τ (k, k ), and, for τ = τ , the short-range intervalley term,W τ τ (k, k ). For the spin singlet (χ σσ = δ σσ ), the unscreened direct term must be subtracted from the dressed interaction,W τ τ (k, k ) →W τ τ (k, k ) − with the pseudo exciton wave function defined as This shows that, at the onset of the EI phase, when ∆(τ k) is infinitesimal-at the critical magnetic field-the exciton wave function for ε u = 0 is the same as ϕ apart from a constant, ϕ(τ k) ∼ where T is the temperature and k B is Boltzmann constant.
The quasiparticles of the EI are the free electrons and holes. For example, in the simplest case of the spin-singlet EI (χ σσ = δ σσ ), the electron quasiparticle wave function Ψ τ k↑ EI differs from the ground state |Ψ EI as the conduction electron state labeled by (τ, k, ↑) is occupied with probability one as well as the corresponding valence state: Here the symbol means that the dummy indices τ k take all values but τ k. The quasiparticle energy dispersion is with the reference chemical potential being zero, as for the noninteracting undoped ground state.
E(τ k) is increased quadratically by the amount |∆(τ k)| with respect to the noninteracting energy, ε(τ k) = E eh (τ, k)/2. This extra energy cost is a collective effect reminescent of the exciton binding energy, since now the exciton condensate must be ionized to unbind one e-h pair and hence have a free electron and hole.

Supplementary Note 5 Inversion symmetry breaking in the excitonic insulator phase
Carbon nanotubes inherit from graphene fundamental symmetries such as time reversal and spatial inversion. Time reversalT swaps K and K valleys whereas the inversionÎ is a π rotation around an axis perpendicular to the tube surface and located in the origin of one of the frames shown in Supplementary Fig. 1. This swaps the valleys as well as the A and B sublattices. Whereas the noninteracting ground state |0 is invariant under both inversion and time reversal,T |0 = |0 and I |0 = |0 , the EI ground state breaks the inversion symmetry 16 . Here we consider a spin-singlet exciton condensate (χ σσ = δ σσ ) withT |Ψ EI = |Ψ EI , hence the excitonic order parameter is real, η = 0, π (otherwise the EI ground state would exhibit transverse orbital currents).
To see that the inversion symmetry of the EI ground state is broken we use the following transformations (whose details are given in Supplementary Note 7): where the shorthand −τ labels the valley other than τ . The transformed ground state iŝ where we have used the fact that u τ k = u * τ k = u −τ −k and v τ k = v * τ k = v −τ −k , as a consequence of time reversal symmetry. The original and transformed ground states are orthogonal in the thermo-dynamic limit, since u 2 − v 2 < 1. On the contrary, 0|Î |0 = 1. Therefore, the symmetry of the EI ground state is lower than that of the noninteracting ground state so the EI phase has broken inversion symmetry, i.e., charge is displaced from A to B sublattice or vice versa.

Supplementary Note 6 Charge displacement between A and B sublattices
In this section we compute the charge displacement between A and B carbon sublattices in the EI ground state. To this aim we must average over the ground state the space-resolved charge density where the sum runs over all electrons in the Dirac valleys. The explicit form of the charge density, in second quantization, iŝ We recall that the states of conduction (α = c) and valence (α = v) bands appearing in Supplementary Eq. (42), ϕ ατ k (r), are products of the envelope functions F times the Bloch states ψ τ at Brillouin zone corners τ = K, K , where ψ τ A (r) [ψ τ B (r)] is the component on the A (B) sublattice. Neglecting products of functions localized on different sublattices, like ψ * τ A ψ τ B , as well as products of operators non diagonal in τ and k indices, which are immaterial when averaging over the ground state, one obtains: The first and second line on the right hand side of Supplementary Eq. (44) are respectively the intra and interband contribution to the charge density. Only the intraband contribution survives when averagingˆ over |0 , providing the noninteracting system with the uniform background charge density 0 (r), with k 1 = A/a. Since |ψ KA (r)| = |ψ K A (r)| = |ψ A (r)|, and similarly for B, this expression may be further simplified as It is clear from this equation that 0 is obtained by localizing the two π-band electrons uniformly on each sublattice site. When averagingˆ over |Ψ EI , the charge density (r) exhibits an additional interband contribution, which is proportional to τ k u τ k v τ k and hence related to the EI order parameter. This term, whose origin is similar to that of the transition density shown in Fig. 3d of main text, as it takes into account the polarization charge fluctuation between |0 and a state with one or more e-h pairs excited, is driven from the long-range excitonic correlations. Importantly, the charge displacement is uniform among all sites of a given sublattice and changes sign with sublattice, the sign depending on the phase of the exciton condensate, η. The charge displacement per electron, ∆e/e, on-say- which is the same as Eq. (3) of main text. In order to evaluate numerically ∆e/e, for the sake of simplicity we neglect the exchange terms splitting the triplet and singlet order parameters (i.e., we assume w 1 = 0). The quantum Monte Carlo order parameter AB defined in the main text is, in absolute value, twice |∆e/e| as there are two relevant electrons per site.

Supplementary Note 7
Reference frame and symmetry operations The reference frame of the armchair carbon nanotube shown in Supplementary Fig. 1a is obtained by rigidly translating the frame used by Ando in a series of papers 2,9,13 , recalled in Supplementary In the frame of Supplementary Fig. 1a used throughout this Supplementary Information, the vectors locating the sites of A and B sublattices are respectively and where (n a , n b ) is a couple of integers and R A 0 (R B 0 ) is the basis vector pointing to the origin of the A (B) sublattice. Besides, one has a ≡ a( where N is the number of sublattice sites, φ π (r) is the 2p z carbon orbital perpendicular to the graphene plane, normalized as in Secchi & Rontani 10 , and ω = exp (i2π/3).
The relative phase between the two sublattice components of Bloch states within each valley, shown in Supplementary Eq. (51), is determined by symmetry considerations 18 . Specifically, the sublattice pseudospinor transforms as a valley-specific irreducible representation of the symmetry point group of the triangle, C 3v : In addition, the relative phase between Bloch states of different valleys is fixed by exploiting the additional C 2 symmetry. The latter consists of a rotation of a π angle around the axis perpendicular to the graphene plane and intercepting the frame origin. This rotation, which in the xy space is equivalent to the inversionÎ, swaps K and K valleys as well as A and B sublattices. With the choice of phases explicited in Supplementary Eq. (51) the inversion operatorÎ takes the form whereR is the inversion operator in the xy space. In contrast, the time-reversal operatorT swaps valleys but not sublattices,T whereK is the complex-conjugation operator. The orthogonal time-reversal of Supplementary Eq. (54) should not be confused with the symplectic transformation 19 , which does not exchange valleys.
The magnetic field along the tube axis breaks bothÎ andT symmetries. However, the reflection symmetry y → −y along the tube axis still swaps the valleys (but not sublattices), as it may be easily seen from a judicious choice of K and K Dirac points. This protects the degeneracy of states belonging to different valleys in the presence of a magnetic field.

Supplementary Discussion
Effects of Dirac cone asymmetry and magnetic field on the exciton wave function The origin of the asymmetry of the exciton wave function in k space, illustrated by Fig. 3b of main text, may be understood within the effective mass model applied to a single Dirac valley-say K.
In the presence of a vanishing gap, electrons (and excitons) acquire a chiral quantum number, C, which was defined above. With reference to the noninteracting ground state, |0 , the e-h pairŝ have chiral quantum number ∆C = +2 for positive k and ∆C = −2 for negative k.
Since long-range Coulomb interaction conserves chirality, we expect the wave function of a chiral exciton to live only on one semi-axis in k space, either ψ K (k) = 0 for k < 0 and ∆C = +2, or In the symmetric case ( Supplementary Fig. 10a) ψ K (k) is even in k since nothing prevents the numerical diagonalization routine from mixing the two degenerate amplitude distributions with ∆C = ±2. Hovever, as the Dirac cone symmetry under axis inversion, k → −k, is lifted by energetically favoring e-h pairs with ∆C = −2 ( Supplementary Fig. 10b), the wave function weight collapses on the negative side of the axis. Therefore, the asymmetry of the exciton wave function is explained by the combined effects of chiral symmetry and cone distortion. The EI mean-field wave function as specialization of the QMC variational wave function The QMC variational wave function, |Ψ QMC , is the zero-gap state, |0 , multiplied by the Jastrow factor, J = J 1 J 2 , which accounts for one-and two-body correlations encoding the variational degrees of freedom. In this section we focus on a relevant specialization of the pair Jastrow factor, showing that a proper choice of the two-body term u(r, r ) allows to recover the mean-field EI wave function to first order in u, i.e., |Ψ QMC takes the form Note that the first-order restriction is consistent with the range of validity of EI mean-field theory 22 .
Throughout this section we take J 1 = 1 and suppress spin and valley indices, as they may be included straightforwardly in the derivation, as well as we assume positive order parameter for the sake of clarity (η = 0).
To first order in the two-body factor u, the QMC wave function is where N e is the number of electrons. The Slater determinant Φ 0 in real space is obtained by where |vac is the vacuum with no electrons present. The Fermi field annihilation operatorψ is spanned by the basis of conduction and valence band operators, withψ where the explicit effective-mass form of Bloch states ϕ ck and ϕ vk was given in Supplementary Note 1.
Similarly, we work out the form of Ψ EI in real space, Ψ EI (r 1 , r 2 , . . . , r Ne ) = vac|ψ(r Ne )ψ(r Ne−1 ) . . .ψ(r 1 ) where the valence band states k 1 , k 2 , . . ., k Ne , are filled up to the Dirac point in |0 and we defined where B = k u k is a constant. After expanding the field operatorsψ in the second row onto the basis ofv andĉ [cf. (60) and (59)], we observe that the only non-vanishing contributions consist in products of N e − 1 operatorsv times a single operatorĉ k . Sinceĉ k occurs N e times in theψ(r i )'s, with i = 1, . . . , N e , we may write To make progress, we consider the generic operator identitŷ ψ(r)ψ † (r ) +ψ † (r )ψ(r) = δ(r − r ).
Since electrons are mainly localized at honeycomb lattice sites R and there is-on the averageone electron per site (N e = 2N ), this identity may approximately be expressed aŝ which provides a useful representation of the identity operatorÎ for any position of the ith electron: occurring in the second row of (68) depends on r − r only, which allows to decouple the sums over r i and r i − r j , respectively. Then Supplementary Eq. (68) may be rearranged as where, among all addenda of the mixed sum over momenta k and index i, the only non-vanishing contributions are those permutating the annihilation operators applied to |0 that belong to the set The final result is Ψ EI (r 1 , r 2 , . . . , r Ne ) = B 1 +

=1
Φ exc (r ) Φ 0 (r 1 , r 2 , . . . , r Ne ) + (contact term), with the exciton wave function Φ exc being defined as where the sum over k is limited to those levels that are filled in |0 and r is the electron-hole of electrons are in contact, QMC and mean-field EI wave functions coincide apart from a normalization factor, provided that u(r, r ) = 2Φ exc (r − r )/N . When two electrons touch, say r i = r j , a discrepancy arises, which is expected since the QMC wave function enforces the cusp condition whereas the mean-field ansatz does not.

Detection of Peierls charge density wave through the order parameter Transl
The QMC analysis of main text introduces the order parameter Transl as a measure of the charge displacement between adjacent unitary cells along the tube axis. If the ground state is a charge density wave (CDW) with period 2a (the characteristic wave vector is q = π/a), then the quantum average of Transl extrapolated to the thermodynamic limit is finite. In this section we discuss whether the order parameter Transl may also detect a Peierls CDW with nesting vector q = 2k F , the Fermi wave vector being located at Dirac point K.
A first issue is the commensurability of the QMC supercell with respect to the period of Peierls CDW. According to DFT calculation k F = 0.289(2π)/a, hence q = 2k F = 0.422(2π)/a (folded back to first Brillouin zone) and the period is 2.37 a (Supplementary Fig. 12b). This implies that the size of the commensurate supercell exceeds our computational capability. On the other hand, the size of a smaller supercell may approximately match a multiple of the Peierls CDW  Fig. 12 Model charge density wave ground state. a-c Model charge density along the axis, n q (z) − n 0 , vs axial coordinate, z. The density is reference from its average value, n 0 , and the blue (red) colour stands for positive (negative) charge deviation. The wave vector q identifies the period of the charge density wave as (2π)/ |q|: undistorted structure, q = (2π)/a (panel a); charge density waveà la Peierls, q = 2k F , with k F ≡ K (panel b); dimerized charge density wave, q = π/a (panel c). d Square order parameter, model period. This is the case e.g. of a supercell made of seven units, whose length compares with three times the period, 7.11 a.
The key issue is the finite-size scaling of Transl averaged over the Peierls CDW ground state.
To gain a better understanding, we introduce a simple model for a generic CDW. The charge density profile, n q (z), is a sinusoidal modulation of wave vector q along the axis z, n q (z) = n mod sin qz + n 0 , where n mod is the modulation amplitude, n 0 is the homogeneous background, and we ignore the relaxation of the ground state occurring in a finite-size supercell. The order parameter model Transl that fits to the model (73) is where N cell is a number of unitary cells such that N cell a is approximately commensurate with the CDW period, 2π/q. If N cell is even, then, except for a prefactor, model Transl is equivalent to Transl as defined in the main text.
The extrapolated value of model Transl in the thermodynamic limit, N cell → ∞, is trivial in two cases. For the undistorted structure, Transl = 0 as the integral of n q − n 0 over the unitary cell vanishes. This is illustrated in Supplementary Fig. 12a, where the blue (red) colour stands for positive (negative) charge deviation, n q (z) − n 0 . Second, for the dimerized CDW of period 2a, which is discussed in the main text, any cell with N cell even is commensurate and hence model Transl = 2an mod /π, the integral of n q − n 0 over the unitary cell exhibiting alternate sign between adjacent cells ( Supplementary Fig. 12c). In the following, in order to compare with the VQMC extrapolated order parameter AB discussed in the main text, we take n mod a/π = AB /2 = 0.00825.
We now focus on the Peierls case of nesting vector q = 2k F (Supplementary Fig. 12b)