Extracting the Single-Particle Gap in Carbon Nanotubes with Lattice Quantum Monte Carlo

We show how lattice Quantum Monte Carlo simulations can be used to calculate 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 and compare our results to empirical data of the Anti-Ferromagnetic Mott Insulating gap in large diameter tubes.


Introduction
Lattice Quantum Monte Carlo methods form an ideal approach for studying non-perturbative phenomena of strongly correlated electrons in low-dimensional crystal systems. Here the ions of the lattice define the physical, spatial lattice, with the lattice spacing a given by the distance between ions. As the lattice spacing is typically in the ångström range (a = 1.42 Å for graphene), the scale is set and no spatial continuum limit is required. From naïve dimensional arguments alone, the low-dimensionality of these systems enhances non-perturbative effects.
In these proceedings we report on progress made in investigating the effects of strong electron correlations in a (quasi) 1-D lattice: carbon nanotubes. In particular, we employ standard lattice QCD methods to present calculations of the Anti-Ferromagnetic Mott Insulating (AFMI) gap for tubes of various diameters. As these tubes have effectively reduced dimensionality (they can be viewed as rolled-up sheets of graphene), the effects of strong correlations are pronounced and both theoretical findings [1] and empirical data confirm this finding [2][3][4]. We compare our results to empirical data.

Nanotube system
The honeycomb (hexagonal) lattice, i.e. graphene, serves as the fundamental structure from which one can construct nanotubes of different types and shapes. This lattice can be viewed as consisting of two underlying triangular lattices, as shown in the left panel of Figure 1. Depending on the geometry of the tube, various properties intrinsic to graphene are inherited by the tube.

Geometry
Single-wall nanotubes can be viewed as rolled-up sheets of a 2-D honeycomb lattice of carbon ions. Each tube is characterized by two integers (n, m) ("chirality"), which define the vector C h = n a 1 + m a 2 that wraps around the circumference of the tube, and a perpendicular vector T whose length gives the minimal tube unit along its axis [5]. The left panel of Figure 1 shows these vectors for the 2-D lattice and the right panel the corresponding tube for (3,3) chirality. For all calculations presented here, we work with n=m ("armchair") geometries and values n ≥ 10 to minimize systematic errors due to the curvature of the tube.

Hamiltonian
Because of the localized nature of the p z -orbital electron at each ion, the dynamics of electrons on the lattice is well approximated by nearest neighbor hopping (i.e. tight-binding), and the effects of electron self-interactions are captured via an onsite interaction, providing the basis for the Hubbard model, where is the charge operator at ion site i, κ is the hopping parameter, U/κ is the Hubbard ratio, a † x (a x ) is a spin-↑ electron creation (annihilation) operator at location x, b † x (b x ) is a spin-↓ hole creation (annihilation) operator at location x and x, y indicates sums over all nearest neighbors. That electrons repel one another provides a basis for the assumption U ≥ 0. This Hamiltonian represents the electrically neutral, half-filling system, where the average number of electrons per ion site is one. Figure 2. The non-interacting dispersion for a (10, 10) armchair tube. The black lines correspond to a tube of infinite length, while red dots are for a tube of finite length L = 15| T |. Positive energies correspond to electrons, while negative energies correspond to holes. Note the two Dirac points where the dispersion for electrons and holes coincide.

Non-interacting dispersion
In the non-interacting case of U = 0, the dispersion relation for a single particle can be determined analytically [5,6]. Because of the finite length around the circumference of the tube, the dispersion exhibits discrete bands. If the tube is assumed to be infinitely long, then the bands themselves are continuous lines. In practical calculations, the tube is of finite length with periodic boundary conditions applied at its ends (see right panel of Figure 1). The dispersion relation in this case will consist of discrete points along these bands. Figure 2 shows a sample dispersion for the (10, 10) tube. All armchair configurations (n, n) exhibit two isolated points within the first Brillioun zone, called Dirac points, where the dispersion for filled (holes) and unfilled (electrons) bands intersect at E = 0 1 , indicating metallic behavior. Around these Dirac points the dispersion is linear and therefore relativistic, though the characteristic speed is not that of light. In the presence of electron interactions (i.e. U 0) a gapped anti-ferromagnetic Mott insulating (AFMI) state may be formed, separating the bands and removing the Dirac points. The Dirac points occur at non-trivial momenta | T · k| = ± 2π 3 , 2π √ 3 and appropriate momentum projection is required to access information there.

Monte Carlo formulation
The treatment of Eq. (1) using the path integral formalism with inverse temperature β has been presented elsewhere [7][8][9][10]. Here we give a cursory description, concentrating on novel aspects relevant to our work.

Setting up the Path Integral
We consider the expression e −βH split into N t "time slices" according to We insert a complete set of fermionic states ψ (for electrons) and η (for holes) between each exponential to obtain the partition function with the condition that ψ N t = −ψ 0 and η N t = −η 0 (i.e. anti-periodic temporal boundary conditions).
Within each matrix element we then introduce an auxiliary field via the Hubbard-Stratonovich transformation, giving where we have introduced the dimensionless variables κ ≡ δκ,Ũ ≡ δU,φ ≡ δφ .
Once standard rules for evaluating matrix elements of normal-ordered operators with fermionic fields are implemented, the integration over fermionic degrees of freedom can be formally done to give where Dφ is a shorthand notation for N t −1 x,t=0 dφ x,t . The functional fermion matrix M is where δ x,y is unity if x and y are nearest-neighbor sites, and zero otherwise. Equation (7) shows a discretization where the forward differencing in time was done at all lattice sites. Backward differencing would be equally valid. In our calculations, we employ both forward and backward differencing on the different sublattices [10]. The probability measure P φ is positive definite, therefore standard Monte Carlo techniques can be employed to generate an ensemble Φ of the auxiliary fieldφ with probability distribution P φ . For this work we utilize pseudo-fermion fields to exponentiate the determinant and employ the hybrid Monte Carlo algorithm [11] to generate such distributions [9,10]. We generate ensembles for different temporal discretizations and tube characteristics.

Momentum wall sources
To extract the quasi-particle dispersion, we analyze the long-time behavior of the single-particle correlator, where N cfg is the number of configurations in Φ. The correlator above contains information of the spectrum at all allowed momentum points. To extract the spectrum at a particular momentum k ≡ (k x , k y ), we project onto specific momenta at both the source and sink. We employ momentum wall sources (for each time slice τ) to perform the momentum projection at the source by solving where N i is the total number of lattice sites and x i refer to sites on one of the sublattices only. It is easy to see that the solution to M is Finally, we Fourier transform the remaining coordinates to momentum project at the sink We note that there are subtle aspects involved in performing the momentum projection because of the two underlying sublattices. These aspects are discussed in detail in Ref. [10]. It remains, however, that for τ − τ ≡ t 0 (but t β), the correlator M −1 k, t;φ ∝ e −E k t , where E k is the desired quasi-particle energy at momentum k. To obtain the value of the gap, we project the correlator at one of the Dirac momenta | T · k D | = ± 2π 3 , 2π √ 3 , and analyze the large time behavior to extract E k D . The gap is given by ∆ = 2E k D .
Because of the low dimensionality of our problems, we calculate wall sources for all possible time slices at the source. This provides us with the means to equivalently calculate the momentum projection of all-to-all correlators.

Results
In Table 1 we enumerate the types of systems and their parameters that we have investigated. In principle, our calculations with different L = N u | T | (tube length) allow us to perform an infinite length extrapolation, while our multiple N t runs allow us to perform a continuum limit extrapolation in time. For weak U/κ 1, we found that our results are nearly constant in L, suggesting that N u is large enough such that our results are already near their infinite volume limits. For larger couplings we expect the dependence on length to be more pronounced, and we are currently investigating this behavior.   Regardless of coupling, all our results showed strong dependence on N t , and thus a temporal continuum extrapolation was essential. Figure 3 shows these findings. Figure 4 shows explicit examples of the dependence of the gap ∆ on the timestep N t , and our extrapolated continuum results using the function ∆(δ) = ∆ 0 + α δ , where δ = β/N t and ∆ 0 is continuum result. Our correlation functions are accurate to O(δ 2 ) [15], and so the calculated effective masses should scale as O(δ), which motivates the use of Eq. (12). We combine our extrapolated points in the right panel of Figure 4 and show the dependence of ∆/κ as a function of Hubbard ratio U/κ for different tube geometries. Also shown in the figure for comparison is the analytical result for the 1-D Hubbard model obtained with the Bethe ansatz [13,14]. The 1-D result is the limit in which the tubes have vanishing radii.  [13,14]. Also indicated in the right plot is the location of the critical coupling (black star) for the 2-D hexagonal Hubbard model [12]. The red points of the two plots coincide with each other.