{'fraction_non_alphanumeric': 0.04039629134427811, 'fraction_numerical': 0.011375841132761985, 'mean_word_length': 4.5281715739731005, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Monte Carlo Calculations for Carbon Nanotubes

We show how lattice Quantum Monte Carlo can be applied to the electronic properties of carbon nanotubes in the presence of strong electron-electron correlations. We employ the path-integral formalism and use methods developed within the lattice QCD community for our numerical work. Our lattice Hamiltonian is closely related to the hexagonal Hubbard model augmented by a long-range electron-electron interaction. We apply our method to the single-quasiparticle spectrum of the (3,3) armchair nanotube configuration, and consider the effects of strong electron-electron correlations. Our approach is equally applicable to other nanotubes, as well as to other carbon nanostructures. We benchmark our Monte Carlo calculations against the two- and four-site Hubbard models, where a direct numerical solution is feasible.

Carbon nanotubes have proven to be a prime testing ground of our knowledge of quantum many-body physics [1][2][3][4]. Viewed as "rolled-up" sheets of its "parent material" graphene [5,6], their electronic properties are closely related to those of graphene [7,8], and depend on how the graphene sheet has been compactified. The allowed momentum modes in a carbon nanotube, for example, are quantized within the two-dimensional Brillouin zone of the graphene sheet (with appropriate use of zone folding). In the absence of electronelectron interactions, graphene exhibits a linear dispersion in the vicinity of the "Dirac points" which are characterized by a Fermi velocity of v F ≃ c/300, where c is the speed of light in vacuum [9,10]. Depending on its geometry, a nanotube can also inherit these Dirac points within its dispersion. The remarkable electronic properties of nanotubes, coupled with the their excellent mechanical and thermal properties, has spurred interest in using them as a replacement for silicon in future electronic applications.
The low dimensionality of graphene (2D), and particularly nanotubes (quasi-1D), provides a good environment for investigating strong-interaction phenomena. For example, the enhanced electron correlation and interaction effects in 1D systems has motivated the Luttinger liquid description of the electronic ground state of nanotubes, where the low-energy excitations consist of bosonic waves of charge and spin [11,12]. In contrast, the properties of 3D metals can often be well described in terms of a Fermi liquid of weakly interacting quasiparticles similar to non-interacting electrons. The possibility of an interaction-induced Mott gap at the Dirac points [13][14][15], particularly in the case of nanotubes, opens the possibility of using these systems as field-effect transistors. Many other phenomena due to strong electron-electron correlations in graphene and nanotubes have been predicted [16][17][18][19][20][21].
Because of electron screening due to underlying substrates and/or surrounding gates, the empirical observation of interaction-driven phenomena in these systems has been surprisingly difficult and for the vast field of applications inspired by these systems, the non-interacting, or tight-binding, picture has proven sufficient. However, experiments with "cleaner" environments (e.g. "suspended" graphene) provide a growing body of empirical evidence for strong electron-electron correlations [22][23][24][25][26][27][28][29][30][31] including, to our knowledge, the only tentative evidence for an interaction-induced gap in the absence of an external magnetic field [32].
In [33], for example, gaps were observed and measured by means of transport spectroscopy in "ultra-clean" samples of nanotubes. Such gaps could not be attributed to curvature effects, and therefore the ground states of nominally metallic carbon nanotubes were identified to be Mott insulators with induced gaps of 10 − 100 meV, with the largest diameter tubes exhibiting the smallest energy gaps. Just as interesting, bound "trions" were observed in doped nanotubes in [34]. In all these cases, the non-perturbative effects of electron-electron correlations cannot be ignored and, at the very least, must be placed on equal footing with other electronic couplings [35].
Monte Carlo methods are well suited for strongly interacting quantum mechanical manybody problems, as exemplified by lattice Quantum Chromodynamics (LQCD). The great advantage offered by the Monte Carlo treatment of the path-integral formalism is that quantum mechanical and thermal fluctuations are fully accounted for, without the need for uncontrolled or ad hoc approximations. For a given Lagrangian or Hamiltonian theory, the Monte Carlo results are regarded as fully ab initio. The systematical errors in any such calculation are due to discretization (non-zero spatial or temporal lattice spacing) or finite volume effects (when studying an infinite system). These errors can be systematically reduced by use of multiple lattice spacings and volumes, and by means of "improved" lattice operators.
Monte Carlo methods have been applied to graphene, using either a "quasi-relativistic" low-energy theory of Dirac fermions valid near the Dirac K points for monolayer [36][37][38][39][40][41] and bilayer [42,43] systems, or applied directly to "tight-binding" models formulated on the physical, underlying honeycomb lattice of graphene, supplemented by a long-range Coulomb interaction which may or may not be screened at short distances [44][45][46]. The former approach is attractive in the sense of being independent of the details of the tight-binding approximation, while the latter appears more amenable to connect with the framework of applied graphene research, and is furthermore closely related to the hexagonal Hubbard model, of which many lattice Monte Carlo studies exist [47][48][49][50]. Notably, in the absence of short-range screening of the Coulomb interaction, both methods predict the opening of a mass gap around a graphene fine-structure constant of α g ≃ 300α ≃ 1, which may be attainable in suspended graphene, unaffected by a supporting dielectric substrate. The AC and DC conductivities of graphene [51,52], its dispersion relation [53], and the effects of an external magnetic field and strain [54][55][56] have also been studied using lattice Monte Carlo methods.
As opposed to graphene, where a gap opens if the coupling α g is above some critical value, for 1D nanotubes it is expected that a gap is induced for any positive value of the coupling (at half-filling). Therefore a non-perturbative Monte Carlo method for nanotubes is quite appropriate. This motivates our introduction of the Monte Carlo method for carbon nanotubes, where we consider (in this paper) the spectrum of a single quasiparticle. While our method is completely applicable to any nanotube configuration (and in principle to other carbon nanostructures as well), we benchmark it for the " (3,3) armchair" tube, which does not exhibit an energy gap in the non-interacting limit. We model our electron-electron in- Previous studies of carbon nanotubes using Density Functional Theory (DFT) [57][58][59] have shown that curvature can significantly distort the band structure of small-radius carbon nanotubes, including changing the electronic properties from semiconducting to insulating and vice versa (for a recent review, see [60]). Such effects become significant for tube radii < 10 Angstroms, which includes the (3,3) nanotube we consider here. However, for armchair nanotube configurations, their symmetry protects them from developing a bandgap due to curvature effects. Moreover, our ultimate objective is to describe large-radius nanotubes, where a Mott insulating state has been experimentally identified [33]. Nevertheless, we discuss how curvature effects can be included into the tight-binding Hamiltonian.
The rest of our paper is structured as follows: In Section II, we summarize the mathematical description of a nanotube, with emphasis on the (3,3) armchair configuration.
The path-integral formalism for nanotubes is given in Section III, along with the lattice formulation that we use for our Monte Carlo calculations. In Section IV we discuss the noninteracting (tight-binding) solution in the context of our path-integral formalism, and its various approximations in discretized form. Section V provides details of our implementa-tion of the long-ranged Coulomb interactions, and the consequences of using this interaction within small dimensions. It also describes the momentum projection method that we use to extract the dispersion energies. As this paper serves as an initial description of the Monte Carlo method applied to nanotubes, we invest significant time in its description in Sections III, IV, and V. The reader interested instead in the results could skip to Section VI, where we present our results for the dispersion relation of the (3,3) armchair nanotube. Here we also discuss our analysis techniques and demonstrate in detail our continuum-limit and infinite-volume extrapolations for the Dirac point energy. We conclude with a recapitulation of our methods and results, and comment on possible future applications. We provide benchmark results of our code in Appendix B.

II. NANOTUBE GEOMETRY
We shall first review the construction of nanotubes from a planar hexagonal lattice, with emphasis on the "(3,3) nanotube" which we shall later use in our lattice Monte Carlo calculations. The geometry of the (3,3) nanotube can be obtained by first considering a planar graphene (honeycomb) lattice, shown in the left panel of Fig. 1. Each point on the graphene lattice can be obtained by integer combinations of the unit vectors where a = 1.42Å is the physical lattice spacing (lattice constant) of graphene. We find for the reciprocal lattice vectors. The hexagonal lattice can also be described in terms of two triangular lattices (labeled A and B), separated by the vector a ≡ ( a 1 + a 2 )/3 as shown in Fig. 1. Such a description of the graphene lattice will be useful for our path integral formulation in Section III.
A general nanotube of "chirality" (n, m) is given in terms of the "chiral vector" C h , where n, m are integers with 0 ≤ |m| ≤ n. The "translation vector" T perpendicular to the chiral vector C h is defined as with where d R ≡ gcd(2m + n, 2n + m) (greatest common divisor). These vectors are shown in Fig. 1 for the case of (n, m) = (3,3).
In order to construct a (3, 3) nanotube, we cut from the graphene lattice the rectangle formed by the chiral and translation vectors. Next, we roll the rectangle along the C h  vector, in order to form a nanotube. Thus, we identify C h as the vector that points along the circumferential direction of the tube, while the vector T points along the longitudinal direction of the tube. This construction represents one "unit cell" of a nanotube of length | T |. The number of hexagons N within this nanotube unit is and for the (3, 3) tube, this gives N U = 6 and | T | = √ 3a.
The length of the tube can be increased by adding additional unit cells to its ends. We denote by N L the number of unit cells along the longitudinal direction, giving an overall tube length of N L | T | and a total number of hexagons N L N U . In our lattice Monte Carlo studies of the (3, 3) nanotube, we use N L = 3, 6, and 9. In the right panel of Fig. 1, we show a (3,3) tube with N L = 3 unit cells. In Table I we summarize the other properties of the nanotubes under consideration.

III. PATH INTEGRAL FORMALISM
We note that detailed treatments of the path integral formalism for a graphene monolayer in the tight-binding description have already been given in Refs. [46,61]. Hence, our main objectives are to give a cursory overview intended to introduce notation, and to highlight the differences encountered in the application to carbon nanotubes. The Hamiltonian H of the carbon nanotube system consists of the tight binding Hamiltonian H tb that describes the interaction of the electrons with the carbon ions, and of the interaction Hamiltonian H I , responsible for electron-electron correlations. We write this in the form where x and y denote sites on the honeycomb lattice, κ ≃ 2.7 eV is the nearest-neighbor hopping amplitude for electrons in graphene, and V x,y is the electron-electron potential matrix (see Section V A). Further, x, y denotes summation over nearest neighbors, and s assumes the values (↑, "spin up") or (↓, "spin down"). Also, q i ≡ a † i,↑ a i,↑ + a † i,↓ a i,↓ − 1 is the charge operator at position i, shifted by (−1) to ensure overall neutrality ("halffilling"). In contrast to Ref. [46], we do not introduce a "staggered mass" term to our Hamiltonian (see Eqn. 10 of Ref. [46]). As mentioned in the introduction, we do not include the effects of curvature in our Hamiltonian, which induces geometrical tilting of π orbitals and hybridization of σ bonds [57]. We note that the effects of curvature can be incorporated into our calculations by using hopping parameters that are dependent on the direction of the three nearest neighbor bonds relative to the tube and azimuthal directions, i.e., κ i for i = 1, 2, 3, as described in [62].
In order to recast the Hamiltonian in a form more amenable to Quantum Monte Carlo calculations, we define the "hole" operators for spin ↓ electrons, and similarly for spin ↑ electrons. In terms of these new operators, Eqn. (10) becomes with the charge operator Finally, we flip the sign of the operators b and b † on one of the sublattices. This impacts only the nearest-neighbor hopping term of the Hamiltonian and leaves the dynamics of the system invariant, as the anticommutation relations of the hole operators remain unchanged. However, this last step is essential in ensuring a positive definite probability measure for our Monte Carlo calculations, as we shall discuss below. The Hamiltonian now becomes where the superfluous spin indices have been dropped.
The basis of our Monte Carlo calculations is Eqn. (13), and we are interested in calculating expectation values of operators O (or time-ordered products of operators), where Z ≡ Tr e −βH is the partition function. The Grassmann-valued fields ψ and η represent electrons and holes, respectively. Their products and sums (denoted by α) are over all fermionic degrees of freedom. Here, β is an inverse temperature and is identified with the temporal extent of our system.
If we now divide e −βH into N t "time slices" according to where δ ≡ β/N t , we may insert a complete set of fermionic coherent states, between each of the factors on the RHS of Eqn. (15). One then arrives at the following expression for the partition function, which depends on the Grassmann fields only. In order to account for the minus sign in the Grassmann fields generated by the trace in Eqn. (14), we identify ψ Nt = −ψ 0 and η Nt = −η 0 , which corresponds to anti-periodic boundary conditions in the temporal dimension.
We now introduce an "auxiliary field" φ by means of a Hubbard-Stratonovich (HS) transformation in the matrix element on the RHS of Eqn. (16), where we have introduced the dimensionless variables κ ≡ δκ,Ṽ ≡ δV,φ ≡ δφ, and we note that Eqn. (17) is valid up to an irrelevant overall constant and rescaling. We note that the stability of this transformation in a Monte Carlo calculation relies on V −1 x,y being positive definite.
We now apply the identity [63] ψ| exp where A x,y is a matrix of c-numbers, to the interaction term. We then obtain [45,52] where we have introduced a "time index" t for the auxiliary field φ x,t . If we insert this expression into Eqn. (16), we find x,y,t [Ṽ ] −1 x,yφx,tφy,t exp κ x,y ,t where Dφ is a shorthand notation for Nt−1 x,t=0 dφ x,t (and similarly for the other fields). The motivation for the HS transformation is now clear: Only quadratic powers of the fermion fields appear in the argument of the exponent (without the HS transformation, quartic powers would also appear). We are now in a position to perform the Gaussian-type integrals over the fermion fields. Up to irrelevant overall factors, the partition function becomes x,y,t=0 where the fermion matrix M is a functional ofφ, where δ x,y equals unity if x and y are nearest-neighbor sites, and zero otherwise. This is referred to as the "compact formulation" of the path integral for the interacting, hexagonal tight-binding system.
The feasibility of a Monte Carlo evaluation of the path integral relies on the generation of configurations ofφ that follow the probability distribution which is positive definite as long as We use global Hybrid Monte Carlo (HMC) lattice updates in order to generate the necessary ensembles of configurations, which we denote by {φ}. For a thorough discussion of the HMC algorithm and related issues, see for example Refs. [46,64]. Given an ensemble {φ}, the Monte Carlo estimate of the expectation value of any operator O is given by whereφ i ∈ {φ} and N cf is the number of configurations within the ensemble. Each such estimate carries with it an associated uncertainty which (in principle) can be arbitrarily reduced with increased statistics (i.e. by taking N cf → ∞). In this first study, we are interested in computing the single quasi-particle spectrum, which can be accessed by taking and by analyzing the temporal behavior of the resulting correlator.
We finally note that the fermion fields can be recast in terms of two-component fields, with one component for the underlying A sublattice and the other one for the B sublattice.
For instance, the electron fields can be written as where x in this case represents the location of a given hexagonal unit cell. In this manner, the ion on site A associated with this particular hexagonal unit cell is located at position x, while the ion on site B is located at x + a. An analogous definition can be made for the auxiliary HS field, and the matrix M(x, t ′ ; y, t) acting on the two-component fermion field is now given by where the coordinates x and y now represent locations of hexagonal unit cells (and not of the ions themselves), so that the definition of δ x,y must be slightly modified to account for all pairs of unit cell locations x and y that share nearest neighbor ions. While we stress that the matrix notation for M(x, y; t) in Eqn. (29) is equivalent to Eqn. (22), the underlying A/B sublattice structure has now been made explicit. We find this representation convenient in analyzing the non-interacting limit of our theory, as discussed Section IV, and also in our zero-mode analysis in Section V C and Appendix A.

IV. NON-INTERACTING SYSTEM
Before we present results of calculations that include electron-electron correlations, it is highly instructive to recall the non-interacting (tight-binding) theory and to compare with the results of our path-integral calculations in this regime. Not only does this exercise allow us to emphasize some salient features of our formalism, but it also allows us to find an accurate way of representing temporal finite differences on the lattice in a way which avoids the infamous "doubling problem", where spurious high-momentum modes contribute in the continuum limit.
A. Zero-temperature continuum limit The non-interacting case is obtained by settingφ = 0 in our expressions for the path integral. In the δ → 0 (continuous time) limit, it is straightforward to show that Eqn. (29) becomes when expressed in terms of dimensionful quantities. We shall first consider the zerotemperature limit, followed by the case of finite temperature. We now move to Fourier space in the β → ∞ (zero temperature) limit, expressing Eqn. (30) as where and following Ref. [65]. In Eqn. (31) we have introduced the momentum variables k || and k ⊥,i which satisfy where T (Eqn. (6)) is parallel to the tube axis. Since we assume that the tube is infinitely long, k || is continuous within an interval of length 2π/| T | 1 . However, the momentum k ⊥ is discrete due to the finite circumference of the tube. These discrete momenta k ⊥,j are given by [65] where the t i are translation vector components and b i the reciprocal lattice vectors, as discussed in Section II. Also, N U is given by Eqn. (9) and j ∈ [0, To determine the zero-temperature dispersion relation for a single quasiparticle in the non-interacting limit, it suffices to study the pole structure of M −1 . This is equivalent to finding simultaneous values of ω and k that satisfy the quantization condition which admits the solution the discrete dispersion to access the Dirac point [65]. In other words, the discrete momentum . This condition is the reason why we focus on tubes with N L = 3, 6, and 9 unit cells.

C. Finite temperature
In addition to calculations with a finite tube length, the path integral formalism requires the introduction of a finite temporal extent β (as discussed in Section III), which in turn can be viewed as an inverse temperature. This implies that the frequency integral should be replaced by the summation where are the Matsubara frequencies. We note that the expression for the correlator can be evaluated analytically using straightforward (though tedious) algebra. In the range 0 < τ < β, we find where and ω( k i ) is given by the positive solution in Eqn. (36). The form of Eqn. (41) is due to the underlying A/B sublattice structure 2 , and admits two linearly independent correlator solutions (see for instance Ref. [66]), which for t ≪ β behave as which shows that the "leading" exponential behavior of these correlators provides access to the (non-interacting) spectrum of the theory. As we show in Section VI, we use this aspect of the correlators when we compute the spectrum in the presence of electron-electron correlations.

D. Discretization of time
We now consider the case where the temporal dimension is also discretized. Given a temporal extent β divided into N t time steps of equal width δ = β/N t , the allowed Matsubara frequencies ω n = 2π T (n + 1/2) are those that fall within the first Brillouin zone [−π/δ, π/δ), which corresponds to −N t /2 ≤ n < N t /2 3 . The time derivative in Eqn. (30) should now be approximated using these discrete steps. As we show below, analytic expressions are still obtainable for the non-interacting case. In what follows, we make use of the representation, where t i = iδ and t j = jδ are lattice time sites (i and j are integers with 0 ≤ i, j < N t ).

Forward difference
We first consider the case of forward discretization to approximate the derivative where f (t) is an arbitrary function on the lattice. Under this differencing scheme, the matrix M in Eqn. (30) in the momentum-frequency domain becomes for which the quantization condition gives the solution for small δ. Hence, we expect our energies computed in this discretized scheme to be shifted by O(δ) from the result in the (temporal) continuum limit. We note that for a backward time difference an analogous derivation gives similar results, provided that the replacement iω n → −iω n is made in Eqs. (51) and (53).

Mixed difference
Given our results for the forward and backward differences, a natural choice would be to consider a symmetric differencing scheme to approximate the time derivative according to although it is well know that this admits spurious high-energy solutions that have no analog in the continuum limit (see for instance the discussion on the "doubling problem" in Ref. [64]) 4 . Instead, we employ a "mixed" differencing scheme where we use a forward difference on A sites and a backward difference on B sites. We are free to do this, since the mixed scheme has the correct continuum limit. This idea was first pointed out in Ref. [61], and we shall use it here for our non-interacting system. With mixed differencing, our fermion matrix becomes and the quantization condition gives which is "O(δ) improved", in comparison with Eqn. (53).
We can visualize this improvement by direct inspection of the correlators. In Fig. 3, we compare the exact analytic correlator at the Γ point to its discretized form, using the forward, backward, and mixed differencing schemes, noting that the time dependence of M is given by Eqn. (39). As can be seen from the left panel of Fig. 3 and we note that our use of this expression is motivated by the improved performance of the mixed-differencing scheme in the non-interacting case.
Since the conclusions of this section were obtained for the non-interacting system, it is not this initial study, we assume O(δ 2 ) scaling to perform our continuum limit extrapolations.
Future calculations with additional values of δ should be able to clarify this scaling with certainty.

V. INTERACTING SYSTEM
Having considered the non-interacting system in some detail, we now turn to the case with electron-electron interactions. In Monte Carlo calculations of graphene, the electrons and holes propagate on the plane defined by the hexagonal graphene sheet, and thus the interaction between the particles is constructed to reflect this geometry. Furthermore, the spatial extent of the system in graphene calculations is typically much larger. In the case of a nanotube, interactions between particles can occur when they are, for example, on opposite sides of the tube wall. Thus, the interaction is not confined to a plane, and the construction of the potential matrix V xy depends on the chirality (n, m) and length of the tube. We now turn to the construction of the potential.

A. Screened Coulomb potential
Our screened Coulomb interaction uses the results of RPA calculations performed by Wehling et al. [44] for the onsite interaction U 00 , nearest neighbor U 01 , next-to-nearest neighbor U 02 , and next-to-next-to-nearest neightbor U 03 interaction. This interaction takes into account the short-distance screening due to the σ-band electrons (which are not dynamic in our calculations i.e. do not hop). We couple this interaction with a potential parameterized as in [46] that ensures the potential approaches the bare Coulomb potential at asymptotic distances. Translational invariance of the potential is maintained by employing a procedure similar to one described in [46]: For any two points x and y on the nanotube, we determine within the tube the shortest distance r between these two points with the ends of the tube identified by periodic boundary conditions. We then assign V (r) as the potential matrix element V xy between these two points. Note that because of periodic boundary conditions at the ends of the tube, there will be cases when r 2 < ( x − y) 2 since the largest value or r || , the component of r parallel to the tube direction, is a| T |N L /2. Due to the finite length of our tube calculations, the infrared divergence of the Coulomb potential is avoided. In Fig. 4 we show the matrix elements of the screened Coulomb potential used in our calculation of the (3, 3) tube with nine unit cells (points). For comparison, we also show the bare Coulomb potential evaluated at the same distances (triangles).

B. Momentum projection
Unlike the non-interacting case, where the quasiparticle spectrum can be directly determined by analyzing the quantization conditions given by the determinant in Eqn. (35), the spectrum of the interacting system must be determined by analyzing the temporal behavior of the appropriate correlator. To access the spectrum at a particular momentum, we must first project our correlator to the corresponding momentum. Such a procedure is routinely performed in lattice QCD calculations. However, we discuss the formalism as it is applied to our system, in order to point out the specific differences to other lattice methods.
where N is the number of unit cells (and not the number of ions). Given a function f ( x i ) of the unit cell coordinates, the above relations can be used to define its Fourier and inverse Fourier transforms In addition to the unit cell locations x i , each unit cell also includes a basis vector a One can make an analogous definition for the hole operator b † ± ( x). In momentum space, the electron correlators are given by and by inserting Eqn. (63) into Eqn. (64), we find where with similar expressions for the other components of the correlator. Note the similarity of this expression to Eqn. (43).
Finally, we emphasize a significant difference with respect to lattice QCD calculations.
In a discretized cubic box of length L with N 3 lattice points x i ≡ a(n x , n y , n z ), where a is the lattice spacing and n i ∈ [0, N) integer, the conjugate momenta are k i = 2π a (l x , l y , l z ) with l i ∈ [− N 2 , N 2 ). The triplet of numbers (l x , l y , l z ) are independent of each other, and thus the momenta in different spatial dimensions can be treated independently. This is in stark contrast to the nanotube case, where for a general tube chirality, the conjugate momenta in the different tube and azimuthal directions cannot be treated independently.

C. Zero-mode analysis
Even though the infrared behavior of the Coulomb interaction is regulated by the finite length of the nanotube, the long-distance nature of the interaction coupled with the small physical dimensions of our calculations provide a setting in which the zero momentum modes of the auxiliary field Φ can introduce non-perturbative contributions. The effects of such "zero-modes" have been investigated in the context of lattice QCD calculations with longrange electromagnetic interactions [67]).
In Table II we give the values of (v AA 0 +v AB 0 )/(2κN) for the different systems we consider in this paper. Though these values were obtained assuming a low-temperature quenched approximation, they nonetheless provide a scale of the expected size of the Gaussian term in our correlators. As we show in the next section, we indeed observe linear behavior in our calculated effective masses which we attribute to the zero-modes of our theory. However, the slopes of the linear terms do not agree with those shown in Table II, and in principle depend on the momentum state of the electron that we are considering. This is to be expected as our numerical simulations are fully dynamical (i.e. they include the determinant terms in Eqn. (21)).

VI. RESULTS
For our Monte Carlo calculations of the (3, 3) nanotube, we consider three different tube lengths, N L = 3, 6, and 9 units, which (in principle) allows us to perform an "infinite volume" (infinite tube length) extrapolation. For each tube length, we generated configurations with N t = 64, 80, and 96, which allowed us to perform a (temporal) continuum limit extrapolation. corresponds to an electron temperature of 0.25 eV. We use the PARDISO package [68][69][70] to perform inversions of sparse matrices within our HMC algorithm.
For the purpose of presentation only, we make use of "effective mass plots", defined by where ∆ is an integer parameter used for statistical analysis. Such a plot provides visual information on the argument of the exponential of the correlator G − (τ ). For example, in  [71]. We also bin our data in order to reduce systematic errors due to autocorrelations. For the results presented below, we bin our data every 100 HMC trajectories.
We have also benchmarked our code to cases where analytic solutions are known, or where solutions can be obtained via direct numerical diagonalization, specifically the twoand four-site Hubbard models. We discuss these benchmark calculations in Appendix B.

A. The Dirac K point
We now describe in detail our analysis for the Dirac K point. In the left panel of Fig. 6 of N t . In the right panel of Fig. 6, we show the corresponding effective mass using ∆ = 2.
Also shown is the expected linear behavior of the effective mass for these calculations in the quenched approximation (lines), as discussed in Section V C. Our dynamical calculations exhibit the same qualitative features as the quenched approximation. In particular, there is a clear linear behavior for the effective mass, particularly for smaller N t . However, the slopes are not as steep as in the quenched approximation. Indeed, for N t = 96, a linear contribution is hardly discernible in the dynamical case for the shown range of time steps, but is nevertheless statistically significant.
As we do not know the analytic form of the slope in the dynamical case, we perform simultaneous fits of both the leading exponential term (which provides the energy) and Gaussian term (which is responsible for the slope) directly to our correlators within a specific time window to extract our spectrum. We stress that we do not perform fits to the effective mass points (Eqn. In Tabs. III-V, we give the extracted energy at the Dirac point for each of our Monte Carlo calculations, as obtained from data with momentum | T |(|k || |, |k ⊥,i |) = 2π 3 , 2π √ 3 . For each value of N L , we perform a continuum (δ → 0) extrapolation using the functional form E(δ) = E 0 + C t δ 2 . The extrapolation is determined by multiple fits of the data, where each data point is sampled according to a normal distribution given by the combined statistical and systematic errors reported in Tabs. III-V. This distribution of fits is then used to estimate the uncertainty of the extrapolation. The results in the continuum limit, and their associated uncertainties, are given in the last columns of Tabs. III-V. In Fig. 8, we show this extrapolation for the Dirac point. As can be seen from Fig. 8, the data points at each value of δ are mutually consistent within uncertainties. This prevents us from determining with certainty that the discrete-time scaling is quadratic in δ. Increased statistics, in addition to calculations at smaller values of δ, would be needed.
Using the continuum-limit results at each N L , we finally perform an infinite-volume, N L → ∞ (infinitely long tube) extrapolation. We find that our data extrapolates well with a simple linear dependence on 1/N L , and therefore we use the following functional form to perform our extrapolation: E(N L ) = E ∞ + C L /N L . Quoted uncertainties of our infinite-volume extrapolation are determined in a similar fashion as our continuum-limit extrapolations. The extrapolation is shown in the bottom right panel of Fig. 8. We find that the energy at the Dirac point is E K /κ = .551 (46). We note that the true volume dependence of our calculations may be something other than linear (see [72], for example, for a discussion of finite-volume scaling within low-dimension systems). However, our three points, and their associated uncertainties, are not sufficient to discern anything that deviates from linear dependence. Future studies, with larger values or N L , should provide answers to this question.

B. Spectrum of the (3,3) carbon nanotube
With our analysis formalism described in the preceding section, we now present the results of the remaining spectrum points for the (3,3) tube. In Fig. 9 we show the effective masses (with Gaussian term subtracted) for all accessible momenta for each of our calculations.
Note that the N L = 3 case has less effective mass points than N L = 6, which in turn has less effective mass points than N L = 9. This is due to the fact that the number of accessible momenta increases as N L increases. In generating these figures, only the G − ( k, τ ) electron correlators were analyzed since statistics for the G + ( k, τ ) correlators were too poor for analysis. Also, correlators with degenerate energies were combined to increase statistics.
A close-up of the effective mass plot for the Γ point, as well as the fitted energy, is shown in the right panel of Fig. 7 for the N L = 9, N t = 96 case.
In tabs. III-V we list all the extracted energies for each N L , N t , and momentum point.
The continuum-limit extrapolation is given in the last column of these tables. Figure 10 shows another example of this extrapolation but this time at the Γ point. As opposed to the Dirac point, the data points at different δ are statistically distinct, but still not sufficient to discern linear or quadratic in δ scaling. To be consistent with the Dirac point analysis, we assume a quadratic scaling for the Γ point as well as all other points on the dispersion.
Again, future studies with smaller values of δ should tell whether such an assumption is valid.
In Fig. 11 we show all continuum-limit results for each N L system, compared to the non-interacting dispersion (dashed lines). As we only analyze the G − (τ ) correlators, and combine degenerate correlators when possible to increase statistics, we only show the upperright portion of the dispersions in this figure (compare with Fig. 2).
To perform an infinite-volume extrapolation, we must use momenta which are common in all N L = 3, 6, and 9 cases, which as can be seen from Fig. 11 occurs for points that have | T ||k || | = 0 and 2π/3. We tabulate these points in Table VI. In the bottom right panel of Fig. 11 we plot these points along with the non-interacting dispersion. The effects of strongly-correlated electrons is clearly seen in the calculated spectrum of this system, and in general lifts all points above their non-interacting values.  Extracted energies for N L = 3 at different N t points. All energies are given in units of the hopping parameter κ. The first uncertainty is statistical, the second systematic. The last column gives the continuum-limit extrapolation assuming a δ 2 scaling. Statistical and systematic errors were combined in quadrature to perform the extrapolation. (15) (7) Table VI.   (28)  found by Ref. [44]. Apart from the requirement that the potential matrix V xy be positive definite, we stress that our formalism is not dependent on any particular parametrization of the electron-electron interaction.
As opposed to previous lattice Monte Carlo simulations, we extracted single quasi-particle energies by direct analysis of the momentum-projected one-body correlators, a method commonly used by LQCD calculations. This allowed us to not only extract the spectrum at the Dirac point, but at all allowed momentum modes. As the nanotube systems studied were relatively small, we were able to perform calculations at multiple time steps N t and multiple tube lengths N L . The former allowed us to perform a continuum limit extrapolation, and the latter allowed us to consider nanotubes of infinite length. In all cases, we found that the non-interacting spectrum is strongly modified by electron-electron correlation effects, in general raising the energies at all points in the Brillouin zone above their non-interacting To account for this effect in our analysis of the large-time behavior of our correlators, we performed simultaneous fits of both Gaussian and (leading) exponential terms.
A possible application of our method would be to consider the energy gap at the Dirac point as a function of nanotube radius, for which considerable experimental data is available.
Such calculations would allow for a direct test of different models for the electron-electron interaction, which in turn could provide additional input to the problem of a possible Mott insulating state in suspended graphene. Also, while we have so far only considered the singlequasiparticle dispersion relation, our formalism can easily be extended to the spectrum of multi-particle states. For example, the interacting J = 0 electron-hole system can be repre- , where the operators a † ± and b † ± are defined in Eqn. (63). The spectrum of such a system could be ascertained by analyzing the temporal behavior of the two-particle correlator O(x ′ , y ′ )O † (x, y) .
In light of the results found in [34], such studies would be very interesting.
Our relatively large uncertainties can be traced back to the non-trivial contribution of zero-modes to our correlators and the fact that our system is physically very small. Both contributions can be alleviated by increasing the volume (length) of our system. As in LQCD, we anticipate a suppression of uncertainties that scale as 1/V 3/2 , where V is the volume of the system [73]. This would correspond to uncertainties that scale as 1/N 1/2 L for our system. Such suppression is already evident when comparing the uncertainties between our N L =3, 6, and 9 calculations (Tables III-V). However, the dimensions of such calculations would scale linearly in N L . We expect reduced uncertainties with larger diameter tube calculations as well. The dimension of a (14,14) tube calculation with N L = 9, for example, is ≃ 20 times larger than the (3,3) system with N L = 9 and same number of time steps.
Such a calculation would require resources beyond what we have committed to this paper, i.e. desktop workstation (we are currently modifying our codes to run on larger computer clusters), and would be ideally suited for GPUs [46,61,74].
In conclusion, we emphasize that our Monte Carlo method is completely general and can be applied to other carbon nanostructures, such as graphene single-and multi-layers, multiwall nanotubes and carbon nano-ribbons. The most significant restriction of our method is the requirement of a positive definite probability measure, the availability of which has to be assessed on a case-by-case basis. In addition to periodic boundary conditions, our method also allows for arbitrary boundary conditions (such as open or twisted boundary conditions). For instance, the latter choice could prove useful in studies of carbon nanotubes in an external magnetic field. where and where m(iω n , k α , Φ x 0 (t 0 )) is the following 2×2 matrix 6 The matrix m(iω n , k α , Φ x 0 (t 0 ) contains no poles. The frequencies ω ± are We first concentrate on the frequency sum in Eqn. (71), which using Eqns. (75)-(77) can be written as Assuming that 0 ≤ τ ≤ β, one can use the Matsubara weighting function h(ω) = (1 + exp(−βω)) −1 and standard finite-temperature integration techniques [75,76] to show that the sum in Eqn. (78) is equal to We concentrate on the small time dependence, τ ≪ β, of our expression and perform a small temperature, large β expansion of the Matsubara regulator, To leading order in this expansion, we have where the function F can be determined by comparing the second and third lines of the equation above. The argument structure of F is written in such a manner as to stress the fact that its depencence on the auxiliary fields is through the difference φ A x 0 (t 0 ) − φ B x 0 (t 0 ), which can be easily verified by analyzing Eqns. (75)-(77).
We now expand our auxiliary fields in momentum-frequency space, where in the second line we explicitly separate the zero-mode contribution. Note that the frequency sum is over bosonic frequencies, ω n = 2πn/β, and N is the number of hexagons in our calculation. Further, it is convenient to define the fieldsφ ± 0 through the following linear combinations of the fields in momentum-frequency space, e −iωnτ 0 e −i kα· x 0 (φ A kα,ωn ±φ B kα,ωn ) 2 (83) where ∆ ± [Φ kα,ωn ] contains sums over terms that have non-zero momentum or frequency modes.
The action in Eqn. (73) can be cast in momentum-frequency space, wherev kα is the discrete fourier transform of the screened Coulomb potential and we have again separated out the zero-mode contribution and defined the remainder as S[Φ k,ω ]. In terms ofφ A,B 0 we have that and a is the basis unit vector. Finally we factor out the zero mode measures in the integration measure, The Jacobian from the change of variables in the last expression is unity. Combining We can now perform the integral overφ + 0 explicitly. Up to an irrelvant multiplicative factor, the result is which shows the Gaussian dependence in τ . The remaining functional integrals overφ − 0 and DΦ (kα,ω) =(0,0) produce exponential depencence in the spectrum E i ( k α ) of the system.
Thus for low temperatures, small time regime, and quenched approximation, our correlator behaves as where i, j denotes nearest neighbor summation, c † i,σ (c i,σ ) is the creation (annihilation) operator for an electron of spin σ at site i, n i,σ = c † i,σ c i,σ is the number operator for spin σ at site i, and U 00 is the onsite repulsive interaction parameter. We note that the relevant dimensionless parameter in this model is simply the ratio λ = U 00 /κ. The eigenvalues of the system can be obtained by direct diagonalization, and the singleelectron correlation function can be obtained using the expression where the sum is over eigenstates |i of the system, E i is the eigenvalue for state |i , and For a given β, U, and κ, we perform our lattice calculations and compare our calculated correlators with those derived analytically. In Fig. 12 we compare our lattice results with analytic results for the case when U/κ = 2. Details of lattice calculation are given in the caption.
The relevant correlators to extract energies are given by Momentum projection on G ± is given by where the sum is over the unit cells (not ion sites). In figs. 15 and 16 we show calculations compared to exact results (determined via diagonalization) for the different momentum projections, using U 00 = 9.3 eV and β = 6.4 eV −1 . Calculations were done with N t =128 and 256, and shows definitive convergence. Also shown are the non-interacting (NI) solutions.
In addition to onsite interactions, we have also benchmarked our codes to 2-and 4-site systems with onsite, nearest-neighbor, and next-to-nearest neighbor interactions. Though we do not show results of these studies here, we find our code gives equally good agreement with analytic and direct diagonalization methods. We note that for systems that have only onsite U 00 and nearest neighbor U 01 interactions, our Monte Carlo code fails due to instability of the Hubbard transformation.