Leonhard Mayrhofer

Milena Grifoni

Linear and nonlinear transport across carbon nanotube quantum dots

Keywords: 73, 63, Fg Nanotubes -71, 10, Pm Fermions in reduced dimensions -73, 23, Hk Coulomb blockade; single-electron tunnelling

We present a low energy-theory for non-linear transport in finite-size interacting single-wall carbon nanotubes. It is based on a microscopic model for the interacting pz electrons and successive bosonization. We consider weak coupling to the leads and derive equations of motion for the reduced density matrix. We focus on the case of large-diameter nanotubes where exchange effects can be neglected. In this situation the energy spectrum is highly degenerate. Due to the multiple degeneracy, diagonal as well as off-diagonal (coherences) elements of the density matrix contribute to the nonlinear transport. At low bias, a four-electron periodicity with a characteristic ratio between adjacent peaks is predicted. Our results are in quantitative agreement with recent experiments.

Introduction

Single-walled carbon nanotubes (SWNTs) are graphene sheets rolled up to seamless cylinders, which possess spectacular mechanical and electrical properties [START_REF] Saito | Physical Properties of Carbon Nanotubes[END_REF][START_REF] Loiseau | Understanding Carbon Nanotubes[END_REF]. In particular, as suggested in the seminal works [3,4,5], due to the peculiar one-dimensional (1D) character of their electronic bands, metallic SWNTs are expected to exhibit Luttinger liquid behavior at low energies, reflected in power-law dependence of various quantities and spincharge separation. Later experimental observations have provided a confirmation of the theory [6,7].

In order to study the internal electronic properties, including the effects of electron -electron correlations, a quantum dot setup, where a SWNT is coupled weakly to a source and drain contact as well as capacitatively to a gate electrode (cf. Fig. 1), is a very well suited device. In fact, the current as a function of the applied bias voltage V b and gate voltage V g depends on the energy spectrum but also on the actual form of the eigenstates of the considered system, which for example determines whether some transitions are allowed or forbidden. At low enough temperatures and bias voltage, adding new particles -and therefore transport across the quantum dot -is hindered by the Coulomb repulsion [8]. Eventually the energy cost for placing an electron in a given energy level adds to the effect of Coulomb blockade. In metallic SWNTs two bands cross at the Fermi energy. Together with the spin degree of freedom this leads to the formation of electron shells, each accommodating up to four electrons.

As a result, a characteristic even-odd [9] or fourfold [10,11,12] periodicity of the Coulomb diamond size as a function of the gate voltage is found. While the Coulomb blockade can be explained merely by the ground state properties of a SWNT, the determination of the current at higher bias voltages requires the inclusion of transitions of the system to electronic excitations. For an interacting 1D system the occurrence of spin-charge separation is expected, i.e. the collective excitations can be divided into independent modes with different spectra. In the case of metallic SWNTs one finds three so called neutral modes [4], whose energies coincide with the corresponding ones of the noninteracting system. Of the three modes two can be identified as collective spin excitations and one as collective charge excitations. Additionally, there is one mode of charge excitations with energies increased compared to the neutral modes due to the repulsive electron -electron interaction. In [11,12] not only Coulomb diamonds but also low lying excitation lines have been resolved in I -V b -V g measurements. Since the associated energies are not large enough in order to identify some of the observed lines with transitions to excitations of the interactiondependent mode, the position of all the measured lines should be explainable by invoking the excitation spectrum of the system with neutral modes only.

Additionally in [12] the effect of an applied magnetic field on the transport properties of a SWNT quantum dot was examined. Clear evidence for an exchange splitting of the ground states with two electrons in the highest occupied shell was found. Exchange effects were also inferred from the excitation lines of the spectrum of samples A and B of [11], in agreement with predictions of the mean-field theory presented in [13]. On the other hand, the measured data of sample C in [11] do not show any sign of exchange effects, unless an unreasonably high exchange energy is assumed. The magnitude of the exchange effects depends on the spatial extension of the electron wave functions and scales like 1/N atoms , where N atoms is the number of carbon atoms in the particular SWNT [15]. Therefore, e.g. the ratio between the level spacing ε 0 of a noninteracting SWNT, scaling like 1/L, and the exchange related energies is roughly proportional to D, where L and D are the nanotube length and diameter respectively. On that score it is interesting to note that the experiments in [12], where the exchange splitting could be seen, was conducted with a SWNT of L = 300 nm and a comparably small diameter of 0.8 nm. Similarly for samples A and B of [11] the reported lengths and diameters were L A = 180 nm and D A = 1.1 nm and L B = 500 nm and D B = 1.3 nm. In contrast, for the SWNT of sample C in [11], showing no measurable exchange effects, L C = 750 nm and D C = 2.7 nm were determined. For comparison, in the case of armchair tubes the diameters 0.8 nm and 2.7 nm would correspond to tubes with the wrapping indices (6,6) and (20,20), respectively. In the following we concentrate on large enough nanotubes and disregard any exchange effects. A detailed discussion of the processes leading to exchange splittings and their interplay with the bosonic excitations is postponed to a forthcoming article [15]. In the absence of exchange, a large degeneracy of the energy spectrum is expected [16], which in turn can be seen in a peculiar four-electron periodicity of the stability diagrams (three equal in size Coulomb diamonds followed by a larger one) and of the Coulomb oscillation traces as discussed below.

So far excitation lines of SWNT quantum dots were only addressed in [13,14,17] in a meanfield approach. In the present work we do not only calculate the expected excitation lines but also give a quantitative calculation of the nonlinear current across a metallic SWNT quantum dot as a function of the gate and bias voltage beyond meanfield. As long as only the ground state energies are concerned, the meanfield treatment yields the right result but fails for a complete description of the excitation spectrum. Some results were already presented in abridged form in [16]. We focus on the low energy regime, which means that we consider the situation where only the gapless subbands of the SWNT are relevant. Since those bands exhibit to a good approximation a linear dispersion relation, the Tomonaga-Luttinger theory can be applied. For a typical SWNT this means a range of roughly 1 eV around the Fermi energy, which is still quite large if for example compared to the level spacing of 1.7 meV for a noninteracting SWNT of 1 µm length.

Starting from a microscopic description of noninteracting SWNTs with open boundary conditions we include the electron -electron interactions. The resulting Hamiltonian can be diagonalized using bosonization techniques [3], obtaining the spectrum and eigenstates of the isolated metallic finite size SWNT. We apply the so called constructive bosonization procedure [18] (for a detailed review we refer to [19]), coping well with the discrete energy spectrum and especially with finite electron numbers.

The electron dynamics of the quantum dot is obtained by solving the equation of motion for the reduced density matrix (RDM) of the SWNT. The metal leads are described as equilibrated Fermi gases. Due to many-fold degenerate eigenstates of the SWNT Hamiltonian we find that offdiagonal elements (coherences) of the RDM can not be generally ignored. Apart from the degeneracy of the states, this is also a consequence of the interactions as we will show. The importance of taking into account coherences for an interacting system weakly coupled to unpolarized leads was also recently discussed in [20], where transport through a metal grain is treated. The rates entering the master equation depend on the involved energies and on the matrix elements of the electron operators in the SWNT eigenstate basis. Noteworthy is a variation of the tunnelling amplitudes along the nanotube axis, also if only ground states are involved in transport, depending on the complete energy spectrum of the collective excitations.

The outline of this article is the following. In Section 2 we derive the master equation describing the dynamics of a generic quantum dot in second order of the tunnelling, gaining an expression for the current. In order to specify the explicit form of the master equation we derive the low energy Hamiltonian of metallic SWNTs in Section 3, and use the bosonization technique in order to diagonalize the SWNT Hamiltonian. The actual calculations for the current in the low and high bias regime are performed in Section 4.

Quantum dots

In this section we show how the stationary current through the system described by the Hamiltonian (1), see below, can be determined as a function of the electrochemical potentials in the leads and in the dot by calculating the dynamics of the RDM of the SWNT. Transport through the SWNT quantum dot can occur if the electrochemical potentials of source and drain are adjusted to different values by the voltages V s and V d . Since current can only flow if one transport direction is favored, the current-carrying system will have to be treated out of equilibrium beyond the linear response regime. Hence, the state of the quantum dot will in general be described by a nonequilibrium density matrix. For completeness and clarity, we show the determination of the corresponding equation of motion in some detail mainly following [START_REF] Blum | Density Matrix Theory and Applications[END_REF]. The outcomes of this and of the following section will be used in Section 4 to obtain the I-V b -V g characteristics of SWNT quantum dots.

Model Hamiltonian

In the following we examine the physics of a generic quantum dot, i.e. we consider a small metallic system weakly coupled to a source and a drain via tunnelling junctions. Moreover the electrochemical potential µ g within the dot can be controlled by a capacitatively coupled gate voltage V g . We describe the overall system by the Hamiltonian

H = H ⊙ + H s + H d + H T + H g , (1)
where

H l = σq ε q,l c † qσl c qσl , (2)
where c † qσl creates a quasi-particle with spin σ and energy ε q,l = ε q -eV s/d in lead s/d. The transfer of electrons between the leads and the central system is taken into account by the tunnelling Hamiltonian

H T = l=s,d σ d 3 r T l (r)Ψ † σ (r)Φ σl (r) + h.c. , (3)
where Ψ † σ (r) and Φ † σl (r) = q φ * q (r)c † qσl are electron creation operators in the dot and in lead l, respectively, and T l (r) describes the transparency of the tunnelling contact at lead l. Finally, H g = -µ g N c = -ecV g accounts for the gate voltage capacitively coupled to the dot, with N c counting the total electron number in the dot and c being a conversion factor, that relates the electrochemical potential to the gate voltage.

Dynamics of the reduced density matrix

Our starting point is the Liouville equation for the time evolution of the density matrix ρ(t) of the system consisting of the leads and the dot. The tunnelling Hamiltonian H T from equation ( 3) is treated as perturbation. We calculate the time dependence of ρ(t) in the interaction picture, i.e. we define

ρ I (t) = U I (t, t 0 )ρ(t 0 )U † I (t, t 0 ), (4)
with the time evolution operator U I (t, t 0 ), given by

U I (t, t 0 ) = e i (H⊙+Hs+H d )(t-t0) e -i (H⊙+Hs+H d +HT )(t-t0) , (5)
with t 0 being some reference time. Using (4) and ( 5) the equation of motion becomes

i ∂ρ I (t) ∂t = H I T (t), ρ I (t) , (6)
with

H I T (t) = e i (H⊙+Hs+H d )(t-t0) H T e -i (H⊙+Hs+H d )(t-t0
) . Equivalently we can write

ρ I (t) = ρ I (t 0 ) - i t t0 H I T (t 1 ), ρ(t 1 ) dt 1 . (7)
Reinserting ( 7) back into (6) yields

ρI (t) = - i H I T , ρ(t 0 ) + i 2 t t0 dt 1 H I T (t), H I T (t 1 ), ρ I (t 1 ) . ( 8
)
Since we are interested in the transport through the central system, it is sufficient to consider the RDM ρ I ⊙ of the dot, which can be obtained from ρ I by tracing out the lead degrees of freedom, i.e.

ρ I ⊙ = T r leads ρ I . (9)
In general the leads can be considered as large systems compared to the dot. Besides we only consider the case of weak tunnelling, such that the influence of the central system on the leads is only marginal. Thus from now on we treat the leads as reservoirs which stay in thermal equilibrium and make the following ansatz to factorise the density matrix ρ I (t) of the total system as

ρ I (t) = ρ I ⊙ (t)ρ s ρ d =: ρ I ⊙ (t)ρ leads , (10)
where ρ s and ρ d are time independent and given by the usual thermal equilibrium expression

ρ s/d = e -β(H s/d -µ s/d N s/d ) Z s/d , (11)
with β the inverse temperature. As it can be formally shown [START_REF] Blum | Density Matrix Theory and Applications[END_REF], the factorization (10) corresponds, like Fermi's Golden Rule, to a second order treatment in the perturbation H T . Furthermore, we can significantly simplify equation (8) by making the so called Markov approximation (MA). The idea is that the dependence of ρI (t) on ρ I (t ′ ) is only local in time. In more detail, ρ I (t ′ ) is replaced by ρ I (t) in (8). The MA is closely connected to the correlation time τ of the leads. In our case τ is the time after which the correlation functions

Φ † σl (r, τ )Φ σl (r ′ , 0) th := T r leads Φ † σl (r, τ )Φ σl (r ′ , 0)ρ leads
are vanishing and therefore erasing the memory of the system. If ρ I (t ′ ) is not considerably changing during τ the MA is valid. It must be noted that the MA leads to an averaging of the time evolution of ρ I (t) on timescales of the order of τ , such that details of the dynamics on short time scales are not accessible. Since we are interested in the dc current through the system, this imposes no restriction on our purpose. Finally we get, by plugging equations ( 9), ( 10) and ( 11) into (8), the following expression for the equation of motion for the RDM,

ρ I ⊙ (t) = - i T r leads H I T , ρ I ⊙ (t 0 )ρ leads + i 2 T r leads t t0 dt 1 H I T (t), H I T (t 1 ), ρ I ⊙ (t)ρ leads . (12)
The first term vanishes because Φ σl (r) th = 0. Since we are only interested in the longterm behavior of the system we send t 0 → -∞. Writing out the double commutator in (12) according to

[A, [B, C]] = ABC + CBA -ACB -BCA,
and introducing the variable

t ′ = t -t 1 one obtains ρ I ⊙ (t) = - 1 2 T r leads ∞ 0 dt ′ H I T (t)H I T (t -t ′ )ρ I ⊙ (t)ρ leads + h.c. -H I T (t)ρ I ⊙ (t)ρ leads H I T (t -t ′ ) + h.c. . (13)
Now we insert the explicit form of H T from equation (3) into (13) and perform the trace over the lead degrees of freedom. Thereby we exploit the relations

Φ σl Φ σl ′ th = Φ † σl Φ † σl ′ th = 0,
as well as

Φ σl Φ † σ ′ l ′ th = 0 for σl = σ ′ l ′ , such that ρI ⊙ (t) = - 1 2 l=s,d σ d 3 x d 3 y ∞ 0 dt ′ E σl (x, y, t ′ )Ψ † σ (x, t)Ψ σ (y, t -t ′ )ρ I ⊙ (t) + F σl (x, y, t ′ )Ψ σ (x, t)Ψ † σ (y, t -t ′ )ρ I ⊙ (t) + h.c. -F * σl (x, y, t ′ )Ψ † σ (x, t)ρ I ⊙ (t)Ψ σ (y, t -t ′ ) + E * σl (x, y, t ′ )Ψ σ (x, t)ρ I ⊙ (t)Ψ † σ (y, t -t ′ ) -h.c. . (14)
In (14) we have introduced

E σl (x, y, t ′ ) := T l (x)T * l (y) Φ σl (x)Φ † σl (y, -t ′ ) th = T l (x)T * l (y) dερ ⊖ l (ε) q|ε φ q (x)φ * q (y)e -i (ε-eV l )t ′ , ( 15
) with ρ ⊖ l (ε) = ρ l (ε)(1 -f (ε))
, where ρ l (ε) is the density of energy levels in lead l and f (ε) is the Fermi distribution. Similarly

F σl (x, y, t ′ ) := T * l (x)T l (y) Φ † σl (x)Φ σl (y, -t ′ ) th = T * l (x)T l (y) dερ ⊕ l (ε) q|ε φ * q (x)φ q (y)e i (ε-eV l )t ′ . ( 16
)
Here ρ ⊕ l (ε) = ρ l (ε)f (ε). In order to proceed, it is convenient to represent the RDM in the eigenstate basis of the dot Hamiltonian H ⊙ . Assuming that we can diagonalize the many-body Hamiltonian H ⊙ (see Section 3), this allows us to extract the t and t ′ dependence of the electron operators in (14). To proceed further, we carry out the following two major approximations: I.) We assume that matrix elements between states representing different charge states vanish (the number of electrons in the dot influences the electrostatics of the whole circuit, hence is "measured" permanently).

II.) The so called secular approximation is applied, i.e. we only retain those terms of (14), which have no oscillatory behavior in t. The latter causes that coherences between nondegenerate states are not present in the stationary solution of (14). For the dynamics this means that we can not resolve the evolution of ρ I ⊙ (t) on time scales of /(E m -E n ), where E m and E n are two distinct energy levels of H ⊙ .

In the end we can divide ρ I ⊙ (t) into block matrices ρ IEN ⊙ (t), where E N denotes an energy level of the isolated dot containing N electrons. To simplify the notation we give the resulting equations of motion in Bloch-Redfield form, ρI,EN nm (t) = -

kk ′ R EN nm kk ′ ρ I,EN kk ′ (t) + M=N ±1 E ′ kk ′ R EN E ′ M nm kk ′ ρ I,E ′ M kk ′ (t), (17)
where the indices n, m, k, k ′ refer to the eigenstates of H ⊙ .

The Redfield tensors are given by (l = s, d)

R EN nm kk ′ = l M,E ′ ,j δ mk ′ Γ (+)EN E ′ M l,njjk + δ nk Γ (-)EN E ′ M l,k ′ jjm , (18) and
R EN E ′ M nm kk ′ = l,α=± Γ (α)E ′ M EN l,k ′ mnk . ( 19
)
The quantities

Γ (α)EN E ′ M l,k ′ mnk
determine the transitions between states with N particles and energy E N to states with M particles and energy E ′ M . In detail we obtain for

transitions N → N + 1 Γ (α)EN E ′ N +1 l k ′ mnk = 1 2 σ d 3 x d 3 y× (Ψ σ (x)) EN E ′ N +1 k ′ m Ψ † σ (y) E ′ N +1 EN nk × ∞ 0 dt ′ F σl (x, y, t ′ )e -α i (E ′ N +1 -EN )t ′ . ( 20
)
For transitions N → N -1 is

Γ (α)EN EN-1 l k ′ mnk = 1 2 σ d 3 x d 3 y× Ψ † σ (x) EN E ′ N -1 k ′ m (Ψ σ (y)) E ′ N -1 EN nk × ∞ 0 dt ′ E σl (x, y, t ′ )e -α i (E ′ N -1 -EN )t ′ , ( 21
)
where we have defined the matrix elements

Ψ † σ (x) EN E ′ N +1 km := k Ψ † σ (x) m , (22)
with the states |k and |m having energy E N , E ′ N +1 and particle number N, N + 1 respectively. Equation ( 17) governs the dynamics of the SWNT electrons. In the following we deduce therefrom the current through the system.

Current

The current is essentially the net tunnelling rate in a certain direction at one of the leads. Thus the current at lead l will be of the form

I l = le N Σ N →N +1 l -Σ N →N -1 l , (23)
where we use the convention l = s/d = ±1. On the long run the currents at the two leads have to be equal, otherwise charge would accumulate on the dot which is prevented by the charging energy. The rates Σ N →N ±1 l can be obtained from the time evolution of the occupation probabilities P N = T r ρ I,N ⊙ , where ρ I,N ⊙ is the RDM for states containing N electrons. In more detail, the occupation probability of the charge state N is reduced by tunnelling events changing the number of electrons from N to N ± 1 and is increased by processes transferring the charge states N ± 1 to N, hence (17) we can easily identify the rates appearing in [START_REF] Avery | Creation and Annihilation Operators[END_REF]:

ṖN = T r ρI,N ⊙ = l -Σ N →N +1 l -Σ N →N -1 l + Σ N +1→N l + Σ N -1→N l . Using T r ρI,N ⊙ = E T r ρI,EN ⊙ together with
Σ N →N ±1 l = E,E ′ nkj Γ (+)EN E ′ N ±1 l,njjk ρ I,EN kn + ρ I,EN nk Γ (-)EN E ′ N ±1 l,kjjn . (24)
Since the relation

Γ (+)EN E ′ N ±1 l,njjk = Γ (-)EN E ′ N ±1 l,kjjn * holds,
we arrive, by inserting ( 24) into [START_REF] Avery | Creation and Annihilation Operators[END_REF], at the following expression,

I l = l2eRe N,E,E ′ Γ (+)EN E ′ N +1 l,njjk -Γ (+)EN E ′ N -1 l,njjk ρ I,EN kn . ( 25
)
For the actual calculations we are going to replace ρ I,EN in ( 25) by the stationary solution ρ I,EN st of ( 14), because we are only interested in the longterm behavior of the system.

Until now our treatment has been quite general, having made no assumptions about the nature of the dot so far. But the actual transport properties depend on the microscopic structure of the system. As seen from equations ( 20) and ( 21), we shall have to determine the spectrum of H ⊙ and the matrix elements (Ψ σ (x))

EN E ′ N ±1 kk ′
in order to solve master equation (17). This is the subject of the next section. Moreover the leads and the geometry of the tunnelling contacts will influence the system via the quantities F σl and E σl as discussed in Section 4.1 and Appendix A.

Low energy description of metallic finite size SWNTs

In the following we derive the Hamiltonian H ⊙ describing the interacting electrons in a metallic finite size SWNT at low energies. Typically (depending on the diameter of the nanotube) "low energies" mean a range of 1 eV around the Fermi energy, such that the dispersion relation for the noninteracting electrons is linear [3] around the Fermi points. Since we are interested in the transport across a finite size system, open boundary conditions (OBCs) for the wave functions must be imposed and we can deduce the wave functions and the energy spectrum of the noninteracting system. Switching on interactions we still can diagonalize the Hamiltonian by using the so called constructive bosonization procedure [18,19], which we also employ in order to rewrite the electron operators in a suitable form for the determination of the matrix elements (Ψ σ (x))

EN E ′ N +1 kk ′ .

Noninteracting SWNTs

The bandstructure of SWNTs is conveniently derived from the one of the unbound 2p z electrons in graphene sheets [START_REF] Saito | Physical Properties of Carbon Nanotubes[END_REF]. Since each unit cell of the graphene lattice contains two atoms p = 1, 2, the band structure consists of a valence and a conduction band which touch at the corner points of the first Brillouin zone. Two of those Fermi points, F = ±K 0 , are independent, i.e. don't differ by a reciprocal lattice vector. The energy dispersion around the Fermi points is very well linearly approximated with a Fermi velocity v F = 8.1 • 10 5 m/s. Since SWNTs are essentially graphene sheets rolled up to a cylinder of a certain diameter, we obtain the SWNT band structure by imposing periodic boundary conditions around the circumference L ⊥ of the tube, leading to the quantization of the allowed transverse wave vectors

OBC 0 L R L R L ~R K -K 0
k ⊥ = 2π L ⊥ m, m = 0, ±1, ±2 . . . , (26)
and thus to the formation of subbands labelled by m. For metallic SWNTs the Fermi points F will satisfy condition [START_REF] Kim | [END_REF], hence the corresponding valence and conduction subbands will have no gap. In the following the focus is on armchair SWNTs. Then only the gapless sub-bands with linear energy dispersion nearby the Fermi points are relevant [3,4]. At each Fermi point there are two different branches r = R/L associated to right and left moving electrons. The corresponding Bloch waves are of the form

ϕ R/L,F,κ (r) = e iκx ϕ R/L,F (r), (27)
where κ measures the distance from the Fermi points ±K 0 (Fig. 2

left).

In more detail, the Bloch waves ϕ R/L,F (r) at the Fermi points are given by a superposition of wave functions living either on sublattice p = 1 or p = 2,

ϕ R/L,F (r) = 1 √ N L R,p e iF •R f prF χ(r -R -τ p ), ( 28
)
where N L is the number of lattice sites which are identified by the lattice vector R. The position of the atoms p = 1, 2 in the unit cell is given by τ p and the functions χ are the p z orbitals. The value of the coefficients f prF depends on the considered SWNT type. For simplicity we concentrate on armchair SWNTs for which we have

f 1rF = 1/ √ 2, f 2rF = -1/ √ 2sgn(rF )
, where we use the convention that R/L = ±1. In this article, we are interested in finite size effects and therefore we impose OBCs instead of periodic boundary conditions on the single electron wave functions. Generalizing [START_REF] Fabrizio | [END_REF] to the case of SWNTs we introduce standing waves which fulfil the OBCs (Fig. 2 right):

ϕ OBC R/ L,κ (r) = 1 √ 2 ϕ R/L,K0,κ (r) -ϕ L/R,-K0,-κ (r) , (29)
with quantization condition

κ = π(m κ + ∆)/L, m κ = 0, ±1, ±2, . . . ,
where L is the SWNT length. The offset parameter ∆ occurs if there is no integer n with K 0 = (πn/L)ê (where ê is the unit vector along the tube axis) and is responsible for a possible energy mismatch between the states of the R and L branches, defined by relation (29). From the linear energy dispersion relation around the Fermi points, the energies of the standing waves follow as

ε rκ = sgn(r) v F κ. ( 30
)
Including the spin degree of freedom, the electron operator reads

Ψ (r) = r= R, L κ,σ ϕ OBC rκ (r)c rσκ =: σ Ψ σ (r), (31)
where the operator c rσκ annihilates the state ϕ OBC rκ |σ . Since the Bloch waves can be divided into a slowly and a fast oscillating part (cf. ( 27)) we can as well split off a slowly varying part from the operators Ψ σ (r), i.e. we introduce the 1D operators

ψ rσF (x) = 1 √ 2L κ e isgn(F )κx c rσκ , (32)
in terms of which the 3D electron operators read

Ψ σ (r) = √ L F sgn(F )× ϕ sgn(F )R,F (r)ψ RσF (x) + ϕ sgn(F )L,F (r)ψ LσF (x) . ( 33
)
Later on, the bosonizability of ψ rσF (x) will be of quite some significance for the calculation of the transport properties. From (30) the appropriate Hamiltonian for the noninteracting system is easily derived. It reads

H ⊙,0 = v F rσ sgn(r) κ κc † rσκ c rσκ . (34)

Interacting SWNTs

For the inclusion of the electron -electron interactions we have to add the following term to the SWNT Hamiltonian,

V ⊙ = 1 2 σσ ′ d 3 r d 3 r ′ Ψ † σ (r)Ψ † σ ′ (r ′ )V (r -r ′ )Ψ σ ′ (r ′ )Ψ σ (r), (35)
where V (rr ′ ) is the possibly screened Coulomb potential. Upon inserting (33) into (35), integration over the coordinates perpendicular to the tube axis yields the interacting Hamiltonian expressed in terms of 1D operators and an effective 1D interaction V eff (x, x ′ ), see (37) below.

As in [3] we expand the Bloch waves into their sublattice contributions (28) and ignore the difference between the intra and inter lattice correlations, which are only of relevance at the length scale of the next neighbor spacing of the carbon atoms. The final form of V ⊙ we obtain by ignoring all fast oscillating terms in the interaction, or in other words by keeping only the so called forward scattering processes associated to the effective potential V eff as depicted in Fig. 3. This leads to an expression of V ⊙ in terms of the 1D densities ρ rσF (x) = ψ † rσF (x)ψ rσF (x),

V ⊙ = 1 2 rr ′ F F ′ σσ ′ L 0 dxdx ′ ρ rσF (x)V eff (x, x ′ )ρ r′ σ ′ F ′ (x ′ ). (36)
The approximations that have been made for getting from (35) to (36) mainly mean that we neglect, as already mentioned in the introduction, any kind of exchange effects and thus are only valid for large enough SWNTs. Using (28) and (31) the effective potential V eff is obtained from an integral over the coordinates perpendicular to the tube axis of the 3D Coulomb interaction weighted by the p z orbitals, i.e.

V eff (x, x ′ ) = L 2 N 2 L × R,R ′ d 2 r ⊥ d 2 r ′ ⊥ |χ(r -R)| 2 V (r-r ′ ) |χ(r ′ -R ′ )| 2 .
(37)

Bosonization

Here we show how the introduction of bosonic excitations enables us to diagonalize the SWNT Hamiltonian H ⊙ = H ⊙,0 + V ⊙ by recasting it into a sum of a fermionic and a bosonic part. Besides we give the identity that expresses the electron operators in terms of the boson operators. The connection between fermionic and bosonic operators can be obtained from the Fourier coefficients of the electron density operators. To be more specific, we Fourier expand the electron density operator,

ρ rσF (x) = 1 2L q e isgn(F )qx ρ rσq , (38)
with q = π L n q , n q ∈ Z and define the operators b σsgn(r)q = 1 √ n q ρ rσsgn(r)q , q > 0.

(39)

The operators b σq fulfil the canonical bosonic commutation relations b σq , b † σ ′ q ′ = δ σσ ′ δ qq ′ as it can be shown [19] by using the explicit expression b σsgn(r)q = 1 √ n q κ c † rσκ c rσκ+sgn(r)q , q > 0.

(40)

From the previous equation we can see that the b, b † operators annihilate or create collective particle hole excitations within the branch (rσ).

Diagonalization of the SWNT Hamiltonian

As for example shown in [19] for a generic 1D system, the noninteracting part of the Hamiltonian is already diagonal in the bosonic operators. For the SWNTs in particular, we obtain from (34) and ( 40),

H ⊙,0 = ε 0 σ   q =0 |n q | b † σq b σq + r 1 2 N 2 rσ + ∆sgn(r)N rσ   , (41)
where

ε 0 = v F π L
is the level spacing of the noninteracting nanotube. The first term in (41) describes collective particle hole excitations, whereas the second term represents the energy cost of the shell filling due to Pauli's principle. Specifically, N rσ = κ c † rσκ c rσκ counts the number of electrons N rσ in the (rσ)-branch. Without loss of generality we have assumed that we only have states with a positive number of particles in each branch.

Plugging the Fourier expansion of the density operators (38) into (36) and taking into account the definition of the b-operators yields the bosonized form of V ⊙ ,

V ⊙ = 1 2 E c N 2 c + 1 2 q>0 n q W qq × σσ ′ rr ′ b σsgn(r)q + b † σsgn(r)q b σ ′ sgn(r ′ )q + b † σ ′ sgn(r ′ )q . ( 42
)
Here E c = W 00 is the SWNT charging energy responsible for the Coulomb blockade and the effective interaction potential is absorbed into

W qq = 1 L 2 L 0 dx L 0 dx ′ V eff (x, x ′ ) cos(qx) cos(qx ′ ). ( 43
)
Since the bosonic operators appear only quadratically in (41) and (42), the Bogoliubov transformation [START_REF] Avery | Creation and Annihilation Operators[END_REF] can be applied to diagonalize the bosonic part of the total SWNT Hamiltonian. Here we only give the result, namely

H ⊙ = 1 2 E c N 2 c + ε 0 rσ N 2 rσ 2 + ∆sgn(r)N rσ + q>0 j=c,s δ=± ε jδq a † jδq a jδq . (44)
The first line of (44) describes the energy cost to add new particles to the system. The excitations are created by the bosonic operators a † jδq . Four channels are associated to total (jδ = c+, s+) and relative (jδ = c-, s-) (with respect to the occupation of the R and L branch) charge and spin excitations. To be more precise, the new operators a jδq are related to the old operators b σsgn(r)q via the Bogoliubov transformation:

b σsgn(r)q = jδ Λ jδ rσ S jδq a jδq + C jδq a † jδq , (45)
where Λ jδ σr is given by

Λ jδ rσ = 1 2    1 1 1 1 1 1 -1 -1 1 -1 1 -1 1 -1 -1 1    , jδ = c+, c-, s+, s rσ = R ↑, R ↓, L ↑, L ↓ .
(46) For the transformation indices S jδq and C jδq we get in the case of the three modes jδ = c-, s+, s-, S jδq = 1 and C jδq = 0.

(47)

Only for jδ = c+ there is an interaction dependence,

S c+q = 1 2 ε 0q ε cq + ε cq ε 0q , C c+q = 1 2 ε 0q ε cq - ε cq ε 0q .
(48) Generalized spin-charge separation occurs, since for the three interaction independent channels jδ = c-, s+, sthe energy dispersion is the same as for the noninteracting system,

ε jδq = v F q = v F π L n q =: ε 0q , n q = 1, 2, . . . ,
but the energies of the c+ channel are enhanced by the repulsive interaction,

ε c+ q = ε 0q (1 + 8W qq /ε 0 ) 1/2 . ( 49
)
The ratio between ε c+ q and ε 0 q for a (20,20)-armchair SWNT is shown in Fig. (4). For the actual calculation of ε c+ q we have substituted |χ(r -R)| 2 by δ(r -R) in (37). The extension of the p z orbitals has been modeled by introducing the average distance a ⊥ ≈ 0.15 nm of the electrons from their nuclei [3], i.e. we used for the 3D interaction potential the following expression The ratio εc+ q /ε0 q as a function of ε0 q for a (20,20) armchair SWNT of 980 nm length. Here we show the decay of εc+ q /ε0 q for an unscreened and a screened (screening length 25 nm) Coulomb interaction. In both cases, a dielectric constant ε of 1.4, see Ref. [3], is assumed.

V (R -R ′ ) = e 2 4πε 0 ε e -|R-R ′ |/Lscreen (R -R ′ ) 2 + a 2 ⊥ . 0
Here L screen is the screening length of the potential and ε 0 ε is the dielectric constant. Due to the finite range of the interaction potential we find a decay of ε c+ q /ε 0 q with increasing q, which should be taken into account if higher excitations in the c+ channel are involved in transport. In the Luttinger liquid theory ε c+ q /ε 0 q is usually assumed to be equal to the constant 1/g, where

g = (1 + 8W qq /ε 0 ) -1/2 , q = π L . (50)
The eigenstates of the SWNT Hamiltonian in (44) are

|N , m := q>0,jδ (m jδq !) -1/2 a † jδq m jδq |N , 0 , (51)
where |N , 0 has no bosonic excitations and the vector N defines the number of electrons in each of the four branches (rσ).

Bosonization of the electron operators

The determination of the transport properties through the SWNT quantum dot involves the calculation of the matrix elements of the electron operators between the eigenstates (51). Therefor we give the so called bosonization identity for the operators ψ rσF (x),

ψ rσF (x) = 1 1 -e -α π L η rσ K rσF (x)e iφ † rσF (x)+iφrσF (x) .
(52) The way (52) is derived can be found e.g. in [19]. Here α is a convergence factor needed to ensure that only excitations with physically sensible energies are considered. It will be set to zero at the end of the calculation. The operator η rσ is the so called Klein factor. Its main effect if acting on a state |N , m is to annihilate a particle in the (rσ)-branch, more precisely

η rσ |N , m = (-1) rσ-1 j=1 Nj |N -êrσ , m ,
where we use the convention j = R ↑, R ↓, L ↑, L ↓ = 1, 2, 3, 4. The fact that the Klein factor annihilates fermions after all, expresses itself in the factor (-1) rσ-1 j=1 Nj . K rσF (x) yields a phase depending on the number of electrons in the (rσ) branch, in our case,

K rσF (x) = 1 √ 2L e i π L sgn(F )(sgn(r)Nrσ +∆)x .
At last we have the fields φ rσF (x) given in terms of the bosonic operators b σq ,

iφ rσF (x) = q>0 e -αq/2
√ n q e isgn(rF )qx b σsgn(r)q .

The matrix elements of the electron operators

In the following we are going to determine the matrix elements (Ψ σ (x))

EN E ′ N ±1 kk ′
, appearing in the transition rates ( 20) and ( 21). As we know from equation (33) the 3D electron operators Ψ σ (x) can be written in terms of the 1D operators ψ rσF (x) yielding

(Ψ σ (x)) EN E ′ N ±1 kk ′ = √ L F sgn(F )× ϕ sgn(F )R,F (r) (ψ RσF (x)) EN E ′ N ±1 kk ′ + ϕ sgn(F )L,F (r) (ψ LσF (x)) EN E ′ N ±1 kk ′ . (53
N , m| ψ rσF (x) |N ′ , m ′ = δ N +êrσ ,N ′ (-1) rσ-1 j=1 Nj × K rF (N ′ )rσ (x) e -1 2 q>0 e -αq jδ |λ jδq rσF (x)| 2 1 -e -α π L =:A(x) × q>0 jδ F (λ jδq rσF (x), m jδq , m ′ jδq ). ( 54
)
The parameters λ jδq rσF (x) (cf. Appendix C) for the three neutral modes jδ = c-, s+, s-are given by

λ jδq rσF (x) = e isgn(F r)qx √ n q Λ jδ rσ , (55)
and for the c+ mode we have The function F (λ, m, m ′ ), describing the dependence of the matrix elements on the bosonic excitations, can be expressed in terms of the Laguerre polynomials L m n as we show explicitly in Appendix B, see equations ( 93) and (94) there. It is interesting to note that the interaction leads to the formation of a non-oscillatory x dependence of the matrix elements k |ψ rσF (x)| k ′ . For matrix elements between states with no bosonic c+ excitations this effect is described solely by the function A(x). From expressions (55) and (56) for the parameters λ jδq rσF we find

λ c+q rσF (x) = 1 2 √ n q ε cq ε 0q cos(qx) + i ε 0q ε cq sgn(F r) sin(qx) . ( 56
)
A(x) = e -q>0 e -αq 8nq
εcq ε 0q cos 2 (qx)+ ε 0q εcq sin 2 (qx)-1 . ( 57
)
Since we are considering a finite ranged interaction, the ratio ε cq /ε 0q goes to 1 for large q and so the sum in (57) converges even for α = 0 (using a finite α in order to disregard unphysical states outside of the low energy regime has only little effect on A(x)). In Fig. 5 we show A(x) for a (20,20) armchair SWNT with screened and unscreened Coulomb potential, using the values of ε cq /ε 0q as shown in Fig. (4). It is evident that A(x) (and hence the tunnelling amplitude at low energies) is smallest at the tube ends and increases towards the middle of the nanotube. For an attractive interaction the opposite behavior would be found and for the noninteracting system A(x) ≡ 1. Hence A(x) is closely related to the dispersion relation for the excitation energies ε cq . However, in this article the contact geometry is fixed and so the actual form of A(x) doesn't influence the transport properties qualitatively, as we will show below when discussing the role of the tunnelling contacts.

Current through a SWNT quantum dot

Knowing the spectrum of H ⊙ and the expressions for the electron operator matrix elements we can return to Section 2.2 and perform the actual calculation of the current through a SWNT quantum dot.

Influence of the tunnelling contacts

Due to the integration d 3 x d 3 y, the rates from ( 20) and ( 21) depend on the geometry of the tunnelling contact. At low enough energies (i.e. as long as no degenerate states with very different bosonic excitations or fermionic configurations are considered), the matrix elements of the 1D operators are slowly varying compared to the extension of the tunnelling contacts described by the transparencies T l (r). Then the influence of the tunnelling contacts at the tube ends can be taken into account by introducing the parameters (cf. Appendix A)

Φ lr r′ (ε) = d 3 r d 3 r ′ T * l (r)T l (r ′ ) q|ε φ * lq (r)φ lq (r ′ )A(x)A(x ′ ) × F F ′ sgn(F F ′ )ϕ sgn(F )r,F (r)ϕ * sgn(F ′ )r ′ ,F ′ (r ′ )ν lF F ′ (∆), (58)
where the phase factor ν lF F ′ (∆) = e iδ l,d πsgn(F -F ′ )∆ accounts for the band mismatch ∆. With the new parameters Φ lr r′ (ε), equation ( 20) can be rewritten as

Γ (α)EN E ′ N +1 l k ′ mnk = 1 2 rr ′ σ dερ ⊕ l (ε)Φ lr r′ (ε) (ψ rσl ) EN E ′ N +1 k ′ m × ψ † r′ σl E ′ N +1 EN nk ∞ 0 dt ′ e α i (ε-eVl-(E ′ N +1 -EN ))t ′ (59)
with ρ ⊕ l (ε) = ρ l (ε)f (ε) as defined below equation ( 16). Additionally we have introduced the notations ψ rσl := ψ rσK0 (x l )/A(x l ) and x l = 0, L for l = s, d. For [START_REF] Blum | Density Matrix Theory and Applications[END_REF] the parameterization can also be performed yielding,

Γ (α)EN EN-1 l k ′ mnk = 1 2 rr ′ σ dερ ⊖ l (ε)Φ * lr r′ (ε) ψ † rσl EN E ′ N -1 k ′ m × (ψ r′ σl ) E ′ N -1 EN nk ∞ 0 dt ′ e α i (-ε+eVl-(E ′ N -1 -EN ))t ′ , ( 60
)
where

ρ ⊖ l (ε) = ρ l (ε)(1 -f (ε))
. The integrals over ε and t ′ in (59) and (60) can be carried out by using

dεg(ε) ∞ 0 dt ′ e α i (ε-E)t ′ = π g(E) + αi P g(ε) ε -E dε,
(61) where P denotes the Cauchy principal value. The first term on the right hand side of (61) corresponds to processes which conserve the energy. Additionally, we have the principal value terms which can be attributed to so called virtual transitions since they cancel in the expression for the current but nevertheless can affect the transport properties indirectly via the time evolution of the RDM.

From now on we assume furthermore that the leads are properly described by a 3D electron gas (think of gold leads for example). In Appendix A we then find for a realistic range of the lead electron wavenumbers q,

Φ lr r′ (ε) = δ rr ′ Φ l (ε), (62)
with

Φ l ∼ R,p |T l (R + τ p )| 2 A 2 (R x ), (63)
simplifying equations ( 59) and (60) further. We should note that the case

Φ lr r′ (ε) = δ rr ′ Φ l (ε) (64)
corresponds to leads which are polarized with respect to the band degree of freedom r, in complete analogy to spin polarized leads. Hence it will be interesting to discuss in the future for which type of contacts, (64) eventually could be fulfilled.

Excitation lines

If certain transitions are possible depends among other things on the available energy and hence on the applied gate and bias voltages. From ( 59) and (60) together with the evaluation of the appearing ∞ 0 dt ′ . . . integral according to (61), we know the resonance condition for tunnelling in/out of lead l. In detail, at temperature T = 0, transitions from a state with N electrons and eigenenergy E N to a state with M = N ± 1 electrons and eigenenergy E M are possible under the condition

eV l ≤ E N -E N +1 + µ g , for M = N + 1, eV l ≥ E N -1 -E N + µ g , for M = N -1. ( 65
)
We expect that for low temperatures the current can only change considerably at the position of the corresponding lines in the eV -µ g plane. But we should mention that in principle the so called virtual transitions (cf. (61)) are sensitive to the values of the bias and gate voltage outside of the excitation lines, if coherences of the RDM influence transport. However, in our case of unpolarized leads we do not find a significant change of the current between two excitation lines. Furthermore not all of the resonance conditions lead to a considerable change of the current at the corresponding lines in the eV b -µ g plane. Relevant transitions are those between states with considerable overlap matrix elements N , m| ψ rσF (x) |N ′ , m ′ and for which the occupation probability of the initial state is large enough. In order to add another electron to a N particle ground state the extra energy δµ

N = E 0 N +1 -E 0 N -E 0 N -E 0 N -1
has to be payed. Here E 0 M denotes the ground state energies. From the SWNT Hamiltonian H ⊙ , equation (44), we deduce

δµ 4m+1 = δµ 4m+3 = E c , δµ 4m+2 = E c + 2|∆|ε 0 , δµ 4m = E c + (1 -2|∆|)ε 0 . (66)
The energy δµ N is a direct measure of the height and the width of the Coulomb diamonds in the eV b -µ g plane. Thus a repeated pattern of one large Coulomb diamond followed by three smaller ones is expected for ∆ = 0, 1/2. Otherwise the pattern will consist of a large diamond followed by a small, a medium and again a small one. In the experiments of [11] sample C showed the first pattern repeating very regularly, whereas samples A and B revealed the second pattern.

For sample C of [11] not only the Coulomb diamonds but also a bunch of excitation lines could be resolved. In [11] the positions of the Coulomb diamonds and of the lowest lying excitation lines were determined. The mean field theory of [13] was used for comparison. Apart from the height of the large Coulomb diamonds, all the lines from the experimental I -V characteristics of sample C could be reproduced by an appropriate choice of five mean field parameters E c , ε 0 , ∆, J and dU . The first three parameters also appear in our theory. The mean field parameters J and dU are the exchange energies with respect to the spin and band degree of freedom. For sample C the choice of [11] was E c = 6.6 meV, ε 0 = 8.7 meV, 2∆ε 0 = J = 2.9 meV and dU = 0 meV. With our theory exactly the same Coulomb diamonds and excitation lines are recovered by choosing E c = 9.5 meV, ε 0 = 2.9 meV and ∆ = 0. We think that in this case our choice of parameters is much more realistic than the one made in [11], with an unreasonably high J of 2.9 meV (compare this to J = 0.5 meV in [12], where a considerably smaller SWNT was examined). We therefore conclude that our treatment of the interaction, where only forward scattering events are considered and exchange contributions are not present, is valid here. In Section 4.4, Figure 9 shows the numerical result for the corresponding current.

Low Bias Regime

In the following we consider the low bias regime, i.e. the bias voltages and the temperature are low enough that only ground states with N and N + 1 particles and energies E 0 N , E 0 N +1 can have a considerable occupation probability. In this case the importance of taking into account the off-diagonal elements of the RDM in ( 17) depends crucially on the parameters Φ lrr ′ (ε) from (58). From the expressions (59) and (60) for the rates it is evident that for our assumption of unpolarized leads, i.e. under condition (62), the time evolution of the RDM elements between states with the same band filling vector N is decoupled from elements between states with different N . Since the current only depends on the time derivative of the diagonal elements of the RDM, the elements mixing states with different N will have no influence on the current and therefore can be ignored. Because all considered ground states do have a different N , in the low bias regime we only have to take into account diagonal matrix elements in the master equation. We will refer to this kind of master equations as "commonly used master equations" (CMEs).

In the low bias regime the occupation probabilities ρ

I,E 0 M N N (1x) (4x) (6x) (4x) D=0:
D=1/8: of the groundstates containing M particles will all be the same in the stationary solution. Therefore, we introduce the probability for finding the system in charge state M by

(1x) (2x) (1x) (2x) L R ~L R ~L R ~L R Ñ=4m N=4m+1 N=4m+2 N=4m+3
P M = d M ρ I,E 0 M N N (t)
, where d M is the degeneracy of the corresponding ground states. Using (17), the general expression for the master equation, we find the following CME for P M , ṖM

d M = -R E 0 M N N N N P M d M + N ′ R E 0 M E 0 M ′ N N N ′ N ′ P M ′ d M ′ . ( 67
)
With equations ( 18) and ( 19) for the Redfield tensors we obtain

ṖM = - l Σ M→M ′ l -Σ M ′ →M l = -2 l N ′ Γ E 0 M E 0 M ′ l N N ′ P M - d M d M ′ Γ E 0 M ′ E 0 M l N ′ N P M ′ , ( 68
)
where we have defined

N ′ Γ E 0 M E 0 M ′ l,N N ′ := N ′ Γ (+)E 0 M E 0 M ′ l N N ′ N ′ N . ( 69
)
Using the relation

d M d M ′ N ′ Γ E 0 M ′ E 0 M l N ′ N = N ′ Γ E 0 M ′ E 0 M l N N ′ (70)
the stationary solution of (68) is easily obtained,

P M = N ′ Γ E 0 M ′ E 0 M l N N ′ Γ tot , (71)
with

Γ tot = l N ′ Γ E 0 N E 0 N +1 l N N ′ + Γ E 0 N +1 E 0 N l N N ′ . ( 72
)
The current is obtained from [START_REF] Avery | Creation and Annihilation Operators[END_REF]. Evaluated at the source for example it yields,

I N,N +1 = 2e N ′ Γ E 0 N E 0 N +1 s N N ′ Γ E 0 N +1 E 0 N d N N ′ -Γ E 0 N +1 E 0 N s N N ′ Γ E 0 N E 0 N +1 d N N ′ Γ tot . ( 73
)
The rates Γ

E 0 M E 0 M ′ l,N N ′
are obtained by using equations ( 59) and (60),

Γ E 0 N E 0 N +1 l N N ′ = π Φ l ρ l f (ε l ) rσ (ψ rσl ) E 0 N E 0 N +1 N N ′ ψ † rσl E 0 N +1 E 0 N N ′ N ,
and

Γ E 0 N +1 E 0 N l N N ′ = π Φ l ρ l (1 -f (ε l )) rσ ψ † rσl E 0 N +1 E 0 N N N ′ (ψ rσl ) E 0 N E 0 N +1 N ′ N .
Here we have assumed that Φ l and ρ l are constant in the relevant energy range. Furthermore we have defined ε l = eV l -∆E and ∆E = E 0 N -E 0 N +1 . From (54) we obtain (remember that ψ rσl := ψ rσK0 (x l )/A(x l )),

rσ N ′ ψ † rσl E 0 M E 0 M ′ N N ′ (ψ rσl ) E 0 M E 0 M ′ N ′ N = 1 2L C M,M ′ ,
where C M,M ′ is the number of ground states with M ′ particles that differ from a given band filling vector N for one of the ground state with M particles only by a unit vector (see Fig. 6). Therefore we get 2

N ′ Γ E 0 N E 0 N +1 l N N ′ = γ l f (ε l )C N,N +1 (74)
as well as 2

N ′ Γ E 0 N +1 E 0 N l N N ′ = γ l (1 -f (ε l )) C N +1,N , (75)
with γ l = π L Φ l ρ l . Inserting the rates (74) and (75) into expression (73) for the current results in

I N,N +1 = e C N,N +1 C N +1,N γ s γ d [f (ε s ) -f (ε d )] l γ l [f (ε l )C N,N +1 + (1 -f (ε l ))C N +1,N ]
.

(76)

Linear conductance

In the regime |eV | ≪ kT ≪ ε 0 we can further simplify (76) by linearising I N N +1 in the bias voltage. We choose -eV s = eV d =: e V b 2 and obtain

I N,N +1 = e 2 βC N,N +1 C N +1,N γ s γ d l γ l (f 0 C N,N +1 + (1 -f 0 )C N +1,N ) e -∆E (e -∆E + 1) 2 V b ,
with f 0 = f (-∆E). Unlike one could expect, the maxima of the conductance G N,N +1 = I N,N +1 /V b are not at ∆E = 0, but at

∆E max = 1 2β ln C N +1,N C N,N +1 ,
which is only zero for C N +1,N = C N,N +1 . The height of the conductance peaks is

G max N,N +1 = γ s γ d C N,N +1 C N +1,N (γ s + γ d ) C N,N +1 + C N +1,N + 2 C N,N +1 C N +1,N e 2 β. (77)
We still have to determine the values of C N,N +1 and C N +1,N , which depend on the mismatch between the R and the L band. If both bands are aligned (∆ = 0, 1/2) one finds from Fig. 6, C N,N +1 = 4, 3, 2, 1 and C N +1,N = 1, 2, 3, 4 for N = 4m, 4m + 1, 4m + 2, 4m + 3. Then the conductance G N,N +1 shows fourfold electron periodicity with two equally high central peaks for N = 4m + 1, 4m + 2 and two smaller ones for N = 4m, 4m + 3 (cf. Fig. 7 a)).

The relative height between central and outer peaks is G max 4m+1,4m+2 /G max 4m,4m+1 = 27/(10 + 4 √ 6) ≈ 1.36. Note that this ratio is independent of a possible asymmetry γ s = γ d in the lead contacts. In addition, the conductance is symmetric under an exchange of the sign of the bias voltage.

In the case of an energy mismatch between the R and the L band exceeding well the thermal energy, the degeneracy of the ground states is either 1 or 2 and we get C N,N +1 = 2, 1, 2, 1 and C N,N +1 = 1, 2, 1, 2 for N = 4m, 4m + 1, 4m + 2, 4m + 3. Therefore, according to (77), all the conductance peaks have the same height as depicted in Fig. 8 a) .

Low bias regime well outside of the Coulomb diamonds

Now we examine the regime where still only ground states are occupied but where we are well outside the region of Coulomb blockade, i.e. ε 0 ≫ |eV l ± ∆E| ≫ kT. Then we have |f (ε s )f (ε d )| = 1. If e.g. eV s -∆E < 0 and eV d -∆E > 0, such that electrons tunnel in from the source and tunnel out at the drain, we have f (ε s ) = 1 and f (ε d ) = 0 such that (76) becomes For mismatched bands and γ s = γ d we find like for the conductance peaks that all current maxima are of the same size. In the case of asymmetric tunnelling contacts, γ s = γ d , a pattern of alternating small and large peaks is found.

I N,N +1 = e C N,N +1 C N +1,N γ s γ d γ s C N,N +1 + γ d C N +1,N . (78)
If we invert the sign of the bias voltage, the current is obtained by flipping its direction and exchanging γ s with γ d in (78). Then, if γ s = γ d , the current does not only change its sign but also changes its magnitude because of C N,N +1 = C N +1,N .

High bias regime

In the bias regime eV > ε 0 not only the ground states will contribute to transport but also states with bosonic excitations and band filling configurations N different from the ground state configurations. We refer to the latter type of excitations as fermionic excitations. Since the number of relevant states increases rapidly with increasing bias voltage, an analytical treatment is not possible any more and we have to resort to numerical methods in order to calculate the stationary solution of the master equation (17) and the respective current. From (65) we know that at low temperatures the current only changes considerably near the excitation lines given therein. Therefore we can reduce drastically the number of (eV, µ g ) points for which we actually perform the numerical calculations, saving computing time. In Figs. 9 a) and 10, the current as a function of the applied bias voltage and the electrochemical potential in the dot is depicted. The chosen parameters for E c , ε 0 and ∆ are the ones we have obtained for fitting the data of sample C and sample A of [11], respectively. Hence Figs. 9 and 10 show the current for SWNTs with the R and L band being aligned (∆ = 0) and mismatched (∆ ≈ 0.17). In both cases a symmetric coupling to the leads, i.e. γ s = γ d , is assumed. In the transport calculations the lowest lying c+ excitations are taken into account using g = 0.21 for the Luttinger parameter (cf. equation (50) for the definition of g).

In addition we have also determined the current using the CME, hence ignoring any coherences in the RDM. For the current corresponding to Fig. 9, the quantitative difference between the calculations with and without coherences is considerable in the region of intermediate bias voltage as we show in Fig. 9 b). On the other hand the deviation of the CME result from the calculation including coherences is by far less pronounced for the parameter choice of Fig. 10. The crucial point here is that the charging energy E c is smaller than ε 0 , the level spacing of the neutral system. Then the subsequent considerations are not strictly valid. Now we explain why coherences can't be generally ignored if considering interacting electrons in a SWNT. Fig. 10. Current as a function of the bias voltage and the electrochemical potential in the SWNT for symmetric contacts. The parameters here are chosen to fit the positions of the Coulomb diamonds of sample A in [11], i.e. Ec = 4.7 meV, 2∆ε0 = 2.8 meV, ε0 = 8.2 meV.

Why and when are coherences needed?

As in the low bias regime we assume to have unpolarized leads, i.e. we use condition (62), such that we can ignore coherences between states with different fermionic configurations N . Unlike in the low bias regime, we are not only left with diagonal elements of the RDM but there still might be coherences between degenerate states which have the same N but different bosonic excitations m. For the importance of these kind of coherences it is illuminating to discuss our system without electron-electron interactions, i.e. for the moment let us assume that an eigenbasis of H ⊙ is given by the Slater determinants of the single electron states ϕ OBC rκ |σ . Furthermore, we concentrate without loss of generality on the case ∆ = 0, and we assume that the charging energy E c exceeds ε 0 . Each of the Slater determinants can be denoted by the occupation n of the single electron states. In the case condition (62) holds, it is again easy to show that coherences vanish in the stationary solution of the master equation (17), if the RDM is expressed in the |n basis. But of course we still could use the states |N , m from (51) as eigenbasis, now with four neutral modes c+, c-, s+, s -. In the |N , m basis it is crucial to include the off diagonal elements in order to get the right stationary solution as we show in the following example.

We adjust the voltages such that only transitions are possible from the ground state with 4m particles and energy E 0 4m to ground and first excited states with 4m + 1 particles and energies E 0 4m+1 and E 1 4m+1 = E 0 4m+1 +ε 0 , respectively. From the 4m + 1 ground states only transitions to the 4m ground state shall be allowed as depicted in Fig. 11 a). Note that this situation is only stable, because we have made the choice E c > ε 0 . The master equation expressed in the |n basis reveals that only four of the 16 states with the lowest particle hole excitation are indeed occupied in the stationary limit (cf. side b) of Fig. 11, because not all of the corresponding energetically allowed transitions from the N = 4m ground state can be mediated by one-electron tunnelling processes. Whereas in the |N , m basis, all 16 states with the energetically lowest bosonic excitations are equally populated. Since the degenerate states of the two bases are connected by a unitary transformation, the same must be true for the corresponding matrix representations of ρ Switching on the electron -electron interactions, the |n states are no longer an eigenbasis of the SWNT Hamiltonian and hence we must work in the |N , m basis. But then, as we know from the discussion above, the coherences are expected to be of importance.

Negative differential conductance

Spin charge separation and therefore non-Fermi liquid behavior could also manifest itself in the occurrence of negative differential conductance (NDC) at certain excitation lines involving transitions to states with fermionic excitations, as was predicted for a spinful Luttinger liquid quantum dot [START_REF] Cavaliere | [END_REF] with asymmetric contacts. We also find this effect for the nonequilibrium treatment of the SWNT quantum dot. Since in our case only the energy spectrum of the c+ mode depends on the interaction and the other three modes have the same energies as the neutral system, rather large asymmetries are needed in order to observe NDC. In Fig. 12 we show the current across the first excitation line for transitions from N = 4m + 1 to N = 4m in the ∆ = 0 case. The corresponding trace in the µ b -V b plane is indicated in the inset of Fig. 12 a). Here the origin of the NDC is that some states with fermionic excitations have lower transition rates than nonexcited states, since due to the increased energy of the c+ modes less channels are available for transport. In Fig. 12 b) we show some of the excited states with N = 4m electrons which are responsible for the NDC, because their transition amplitudes to states with N = 4m + 1 electrons are reduced compared with the one of the 4m ground state. Apart from the asymmetries all other parameters are chosen as for Fig. 9. Only for asymmetries a = γ d /γ s larger than around 45, clear NDC features are seen.

Conclusions

In this article we have analyzed the linear and nonlinear current as a function of the gate and bias voltage across metallic SWNT quantum dots. The properties of the SWNT itself were derived from a microscopic model including electron -electron interactions. Exchange and related effects, which become relevant for small diameter SWNTs, have not been taken into account. The energy spectrum of a metallic SWNT, which can be obtained with the help of bosonization, turned out to be highly degenerate as a consequence of both fermionic and bosonic excitations. In the linear bias regime, the degeneracy of the groundstates leads to a characteristic pattern of conductance peaks depending on whether the two branches of the dispersion relation are aligned or not. Leaving the linear regime, asymmetry effects become relevant. Thus measurements of the current at low bias voltages, in the linear and nonlinear regime, in principle allow the separate determination of the source and drain tunnelling resistances as a function of the gate voltage. At higher bias voltages also excited states become relevant. The correct calculation of the nonequilibrium dynamics of the system then requires the inclusion of coherences in the reduced density matrix between degenerate states with bosonic excitations. At intermediate bias voltages there is a considerable deviation between the transport calculations with and without coherences. We emphasize that for a noninteracting system with unpolarized leads, coherences do Only for asymmetries larger than around 45 negative differential conductance occurs. All parameters here are chosen as for Fig. 9, except for kBT a value of 0.026 meV was used. b) The transitions responsible for NDC. i) In region A of the inset from a) only transitions between ground states are possible at low enough temperature. The transition rate from the N = 4m ground state to the N = 4m + 1 ground states is given by 4γs (see equation ( 74)). ii) In region B additionally excited states become occupied. For some of the states with N = 4m electrons and fermionic excitations the transition rates to states with N = 4m + 1 electrons and neutral bosonic excitations is decreased compared to the ground state rate as a consequence of the larger energies of c+ excitations, which are not available for transport yet. Notice that we only show the most important types of transitions from N = 4m to N = 4m + 1 that take place in region B. Since we are considering a nonequilibrium situation other types of transitions are possible in principle.

not have to be considered if expressing the reduced density matrix in terms of Slater determinants, formed by the one electron wave functions of the noninteracting system. Another consequence of the electron correlations is the formation of a non-oscillatory spatial dependence of the tunnelling amplitudes along the nanotube axis. For transitions between states with energetically low excitations we find a strong suppression of the tunnelling amplitudes near the SWNT ends. Furthermore we have addressed the influence of the tunnelling contacts on the transport. We have shown that extended contacts described as 3D Fermi gas do not lead to a polarization of the contacts with respect to the two branches of the dispersion relation. We think that a further investigation of this point for other types of contacts is worthwhile.

*** Useful discussions with S. Sapmaz and support by the DFG under the program GRK 638 are acknowledged.

A Description of the tunnelling contacts

In this appendix we discuss the dependence of the rates Γ

(α)EM E M ′ l,k ′ mnk
on the properties of the tunnelling contacts. In Section 4 we have introduced the parameters Φ lrr ′ (ε) with equation ( 58). Here we show the calculation of these parameters assuming a 3D electron gas in the leads. Furthermore, we require that the tunnelling region extends over several sites of the SWNT lattice. Starting from (20) together with ( 16) and (33) we can characterize the part of Γ

(α)EN E ′ N +1 l,k ′ mnk
that depends on the tunnelling contact l by the following expression,

T l,nkk ′ m (ε) := F F ′ sgn(F F ′ ) d 3 r d 3 r ′ T * l (r)T l (r ′ )× q|ε φ * lq (r)φ lq (r ′ ) rr ′ σ φsgn(F )r F (r) φ * sgn(F ′ )r ′ F ′ (r ′ )× (ψ rF σ (x)) nk ψ † r′ F ′ σ (x ′ ) k ′ m . (79)
The wave functions in the leads are denoted φ lq and the sum q|ε extends over all q values that correspond to the energy ε. At low enough energies, oscillations of the product (ψ rF σ (x)) nk ψ † r′ F ′ σ (y)

k ′ m can be ignored along the length of the tunnelling interfaces. Using (54) we thus can rewrite (79) as

T l,nkk ′ m (ε) = L F F ′ sgn(F F ′ )ν lF F ′ (∆)× d 3 r d 3 r ′ T * l (r)T l (r ′ )A(x)A(x ′ ) q|ε φ * lq (r)φ lq (r ′ )× rr ′ σ ϕ sgn(F )r F (r)ϕ * sgn(F ′ )r ′ F ′ (r ′ ) (ψ rσl ) nk ψ † r′ σl k ′ m . (80)
Here the factor ν lF F ′ (∆) = e iδ l,d πsgn(F -F ′ )∆ takes into account a possible phase shift due to the band mismatch ∆ and we have used the notation ψ rσl = ψ rK0σ (x l )/A(x l ), where x l = 0, L for l = s, d. Now the parameters Φ lr r′ (ε) from (58) are recovered by the relation

T l,nkk ′ m (ε) = rr ′ Φ lr r′ (ε) σ (ψ rσl ) nk ψ † r′ σl k ′ m .
In the next step we exploit that the Bloch waves ϕ sgn(F )r F from equation (28) are only nonvanishing around the positions of the carbon atoms in the SWNT lattice. On the length scale of the extension of the p z orbitals all other quantities in Φ lr r′ are slowly varying. Hence we can rewrite the integrals over x and y as a sum over the positions of the carbon atoms

Φ lrr ′ (ε) = C L N L F F ′ sgn(F F ′ )ν lF F ′ (∆)× p,p ′ f psgn(F )r F f * p ′ sgn(F ′ )r ′ F ′ R,R ′ T * l (x R,p )T l (x R ′ ,p ′ )× A(R x )A(R ′ x ) q|ε φ * lq (x R,p )φ lq (x R ′ ,p ′ )e iF R e -iF R ′ , (81)
where we have defined x R,p := R + τ p . The constant C results from the integration over the p z orbitals. For the 3D electron gas the wave functions φ lq (x) are simply given by plane waves,

φ lq (x) = 1 √ V l e iqx .
Therefore we can easily perform the sum over the wave numbers q associated with the energy ε,

q|ε φ * lq (x)φ lq (y) ≈ 2π 0 dϕ 1 -1 d cos θe iq|ε|x-y| cos θ = 4π sin(q| ε |x -y|) q| ε |x -y| .
The previous expression is peaked around x = y. For q| ε > 1/ |x -y| the correlations between x and y drop off fast. It means that if q| ε is larger than 1/a 0 , where a 0 is the nearest neighbor distance on the SWNT lattice, correlations between different carbon atom sites are suppressed such that we arrive at the following approximation:

q|ε φ * lq (x R,p )φ lq (x R ′ ,p ′ ) ≈ 4πδ R,R ′ δ p,p ′ , which yields Φ lrr ′ (ε) = 4πC L N L × F F ′ sgn(F F ′ )ν lF F ′ (∆) p f psgn(F )r F f * psgn(F ′ )r ′ F ′ × R |T l (x R,p )| 2 e i(F -F ′ )R A 2 (R x ). ( 82
)
Since we assume an extended tunnelling region, the fast oscillating terms with F = -F ′ are supposed to cancel. Furthermore we assume that both sublattices are equally well coupled to the contacts, such that the sum over R in (82) should approximately give the same result for p = 1 and p = 2. Hence we can separate the sum over p from the rest. Since the Bloch waves φsgn(F )r F and φsgn(F )r ′ F are orthogonal to each other for r = r ′ , the relation p f psgn(F )r F f * psgn(F )r ′ F = δ rr ′ must hold, as it can be deduced from the explicit expression for the Bloch wave, equation (28), and we arrive at

Φ lr r′ (ε) = δ rr ′ 4πC L N L R |T l (x R,p )| 2 A 2 (R x ). ( 83
)
We should mention that the treatment of the contact geometry here relies on several assumptions that might not be fulfilled under all circumstances and so it should be seen as a first estimate. For example we assume that the transmission functions T l do not depend on the q vector of the incoming wave. But such a dependence would increase the correlation between different atom sites. Furthermore q| ε will not be much larger than 1/a 0 for all kind of contacts (for gold the Fermi wave number is about 1.5/a 0 ), nor can the lead electrons always be described by a 3D electron gas.

Here we want to emphasize, the tunnelling contacts itself play an important role for transport through a SWNT quantum dot. Especially it will be interesting for which kind of contacts the condition Φ lrr ′ (ε) = δ rr ′ Φ l (ε), of having leads polarized with respect to the band degree of freedom r, might be fulfilled.

B The matrix elements of the electron operators

In this appendix we calculate the expressions for the matrix elements N , m| ψ rσF (x) |N ′ , m ′ using the bosonization identity (52) of the electron operators ψ rσF (x), namely e -αq/2 √ n q e isgn(F r)qx C jδq a † jδq + S jδq a jδq -e -isgn(F r)qx S jδq a † jδq + C jδq a jδq (96)

By comparing (96) to (87) the values of the λs can be read off:

λ jδq rσF (x) = 1 √ n q Λ jδ
rσ -e -isgn(F r)qx C jδq + e isgn(F r)qx S jδq .

With the coefficients C jδq , S jδq and Λ jδ σr from (46), ( 47) and (48) we get for jδ = c-, s+, s-, λ jδq rσF (x) = e isgn(F r)qx √ n q Λ jδ rσ (97)

and for the c+ mode we have

λ c+q rσ (x) = 1 2 √ n q
ε cq ε 0q cos(qx) + i ε 0q ε cq sgn(F r) sin(qx) . (98)

Fig. 1 .

 1 Fig. 1. SWNT quantum dot setup. A metallic SWNT is coupled via tunnelling junctions to a source and a drain contact. The electrochemical potential in the dot is adjusted by a gate voltage.

Fig. 2 .

 2 Fig. 2. On the right, energy spectrum of a SWNT with open boundary conditions (OBCs) described in terms of left ( L) and right ( R) branches. It is constructed from suitable combinations of travelling waves (cf. equation (29)) whose spectrum is shown on the left side.

Fig. 3 .

 3 Fig. 3. Scheme of forward scattering terms. All those processes conserve the number of electrons in each branch.

  Fig.4. The ratio εc+ q /ε0 q as a function of ε0 q for a(20,20) armchair SWNT of 980 nm length. Here we show the decay of εc+ q /ε0 q for an unscreened and a screened (screening length 25 nm) Coulomb interaction. In both cases, a dielectric constant ε of 1.4, see Ref.[3], is assumed.

)

  Now the bosonization procedure used for the diagonalization of H ⊙ and the redefinition of ψ rσF (x) in terms of boson operators pays off. Denoting the SWNT eigenstates |k and |k ′ as |N , m and |N ′ , m ′ according to equation (51), we get, as shown in Appendix B,

Fig. 5 .

 5 Fig. 5. The non-oscillatory position dependence of the matrix elements k |ψrσF (x)| k ′ is described by A(x), if no c+ excitations are involved. For a repelling interaction, A(x) is strongly suppressed at the SWNT ends. Here we used the values of εcq/ε0q of a (20,20) armchair SWNT (L = 980 nm) as shown in Fig. (4).

Fig. 6 .

 6 Fig. 6. Ground states of the SWNT for different numbers of electrons. For each charge state we show one representative ground state, on top for aligned bands and on the bottom for mismatched bands. In brackets the corresponding degeneracy is given. Full black and orange arrows represent occupied states. The bold (orange) arrows indicate electrons which can contribute to transitions N → N -1 in the low bias regime. Incoming electrons for the transitions N → N + 1 can be accommodated by the states represented by the dashed arrows. So the number of bold and dashed arrows is equal to CN,N-1 and CN,N+1 respectively.

Fig. 7 .

 7 Fig. 7. Gate traces for ∆ = 0, i.e. for aligned R and L bands, at low bias voltage. Parameters are Ec = 9.5 meV, kBT = 0.10 meV, ε0 = 2.9 meV and the chosen asymmetry is γs = 5γ d = 4.9 • 10 10 s -1 . a) Conductance in the linear regime eV b ≪ kBT ≪ ε0. Despite asymmetric contacts, we find the repeating pattern of two small outer peaks and two large central peaks. b) Current in the regime kBT ≪ |eV l | ≪ ε0. Asymmetry effects appear. The fourfold periodicity is retained. The upper and lower pattern correspond to opposite values of the bias voltage.

Fig. 8 .

 8 Fig. 8. Gate traces for ∆ = 0.2, i.e. for mismatched R and L bands, at low bias voltage. Other parameters are Ec = 9.5 meV, kBT = 0.10 meV, ε0 = 2.9 meV and the chosen asymmetry is γs = 5γ d = 4.9 • 10 10 s -1 . a) Conductance in the linear regime eV b ≪ kBT ≪ ε0. The conductance peaks are all of the same height. b) Current in the regime kBT ≪ |eV l | ≪ ε0.The contact asymmetry leads to alternating large and small current maxima.

Fig. 9 .

 9 Fig. 9. a) Current in a bias voltage -electrochemical potential plane for the symmetric contacts case. b) Difference plot of the current with and without coherences. Here kBT = 0.01 meV and the interaction parameter is g ≈ 0.21. Other parameters are as for Fig. 7.

Fig. 11 .

 11 Fig.11. a) Scheme of relevant energy levels, if transitions from the ground state with N = 4m particles to the first excited state with N = 4m + 1 particles are energetically allowed, while no transitions from the N = 4m + 1 ground states to excited states with N = 4m electrons are possible. b) Possible transitions between the N = 4m and N = 4m + 1 electron states of the noninteracting system, which are energetically allowed in situation a). The degeneracy of the eigenstates is given in brackets.

I,E 1 4m+1⊙.

 1 From the representation of the RDM in the |n basis, we know that the rank of ρ I,E 1 4m+1 ⊙ must be equal to 4. Because an unitary transformation does not change the rank of a matrix, the stationary solution in the |N , m basis can maximally have 4 linearly independent columns. Since all diagonal elements are nonvanishing this is only possible if there are also nonvanishing coherences.

Fig. 12 .

 12 Fig.12. a) Differential conductance as a function of the bias voltage for different asymmetries a = γ d /γs of the coupling to the leads. The trace in the µ b -V b plane across the boundary between regions A and B is indicated in the inset. Only for asymmetries larger than around 45 negative differential conductance occurs. All parameters here are chosen as for Fig.9, except for kBT a value of 0.026 meV was used. b) The transitions responsible for NDC. i) In region A of the inset from a) only transitions between ground states are possible at low enough temperature. The transition rate from the N = 4m ground state to the N = 4m + 1 ground states is given by 4γs (see equation (74)). ii) In region B additionally excited states become occupied. For some of the states with N = 4m electrons and fermionic excitations the transition rates to states with N = 4m + 1 electrons and neutral bosonic excitations is decreased compared to the ground state rate as a consequence of the larger energies of c+ excitations, which are not available for transport yet. Notice that we only show the most important types of transitions from N = 4m to N = 4m + 1 that take place in region B. Since we are considering a nonequilibrium situation other types of transitions are possible in principle.

ψ 1 )K 1 j=1m 1 j=1

 111 rσF (x) = η rσ K rF Nrσ 1e -α π L e iφ † rσF (x)+iφrσF (x) .(84)As a reminder, η rσ is the Klein factor reducing the electron number by one, the operator K rF Nrσ is given byK rF Nrσ (x) := 1 √ 2L e i π L sgnF (sgnr•Nrσ +δ)x(85)and yields a phase factor depending on the filling of the band (rσ). The bosonic fields φ rF σ (x) are given in terms of the bosonic annihilation operators b σq ,iφ rσF (x) = q>0 e isgn(F r)qx-αq/2 √ n q b σsgn(r)q . (86)The action of ψ rF σ (x) on the states |N , m is conveniently determined by rewriting the operators b σsgn(r)q in terms of the operators a jδq and a † jδq which diagonalize the SWNT Hamiltonian H ⊙ . As we know from the Bogoliubov transformation, the b operators depend linearly on the a operators and hence we can write in generaliφ † rσF (x) + iφ rσF (x) = = q>0 e -αq/2 jδ λ jδq rσF (x)a jδq + λ ⊕jδq rσF (x)a † jδq , (87) where λ ⊕jδq rσF (x) =λ jδq rσF (x) * , since iφ † rσF (x)+iφ rσF (x) is anti-hermitian. The actual values of the λs are calculated below in this appendix. Plugging the bosonization identity (84) together with (87) into N , m| ψ rσF(x) |N ′ , m ′ , yields N , m| ψ rσF (x) |N ′ , m ′ = δ N +erσ,N ′ (-rF (N ′ )rσ (x)× m| e q>0 e -αq/2 jδ (λ jδq rσF (x)a jδq +λ ⊕jδq rσF (x)a † jδq ) |m ′ .Remember that the factor (-1) rσ-Nj stems from the Klein factor. The term m| e q>0 e -αq/2 jδ (λ jδq rσF (x)a jδq +λ ⊕jδq rσF (x)a † jδq ) |m ′ (88) does not depend on the fermionic configuration and therefore we have dropped the N index in |N , m . Now we exploit that in our finite size SWNTs there will always be only a finite number of bosonic excitations. Hence, by using the Baker-Hausdorff formula, e A+B = e A e B e -[A,B]/2 , we commute all annihilation operators a jσq to the right and all creation operators a † jδq to the left in (88): m| e q>0 e -αq/2 jδ (λ jδq rσF (x)a jδq +λ ⊕jδq rσF (x)a † jδq ) |m ′ = jδq | e λ ⊕jδq rσF (x)a † jδq e λ jδq rσF (x)a jδq m ′ jδq , (89) with |m jδκ = (m jδκ !) -1/2 a † jδκ m jδκ |0 . Expanding the exponentials e λ jδq rσF (x)a jδq in (89), all terms which are of higher order than m ′ jδq will vanish. Analogously, all terms in the expansion of e λ ⊕jδq rσF (x)a † jδq of higher order than m jδq can be ignored. Hence we get m| e λ ⊕ a † e λa |m ′ =: F (λ, m, m ′ ) = favor of readability all indices have been dropped. In (90) only the terms with mi = m ′j survive. For m ′ > m we conveniently express j in terms of i and get:F (λ, m, m ′ ) = λ m ′ -m m! m ′ ! m i=0 (λ ⊕ λ) i i!(i + m ′m)! m ′ ! (mi)!(91) and for m ′ < m we write i in terms of j and getF (λ, m, m ′ ) = λ ⊕ m-m ′ (-λ * ) m-m ′ m ′ ! m! m ′ j=0 (λ ⊕ λ) j j!(j + mm ′ )! m! (m ′j)! .(92)Defining m max/min := max / min(m, m ′ ), we can summarize (91) and (92) asF (λ, m, m ′ ) = = Θ(m ′m)λ m ′ -m + Θ(mm ′ ) (-λ * ) i + m maxm min )! m max ! (m mini)! ,(93)or in terms of the Laguerre polynomials[START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF],F (λ, m, m ′ ) = m min ! m max ! L mmax-mmin mmin (|λ| 2 )× Θ(m ′m)λ m ′ -m + Θ(mm ′ ) (-λ * ) m-m ′ .(94)Similar expressions for F (λ, m, m ′ ) were found by Kim et al.[START_REF] Kim | [END_REF] when investigating charge plasmons in a spinless Luttinger liquid quantum dot. In the end we obtain for the matrix elements (ψ rσF (x)) nm the following expression:N , m| ψ rσF (x) |N ′ , m ′ = δ N +erσN ′ (-1) rσ-Nj K rF (N ′ )rσ (x)A(x)× q>0 jδ F (λ jδq rσF (x), m jδq , m ′ jδq )We still have to determine the values of the parameters λ jδq rσF (x) for our case. Using equation (86) we can express the sum iφ † rσF (x) + iφ rσF (x) in terms of the bosonic operators b σq , iφ † rσF (x) + iφ rσF (x) = q>0 e -αq/2 √ n q × -e -isgn(F r)qx b † σsgn(r)q + e isgn(F r)qx b σsgn(r)q , (95) which in turn are related to the operators a jδq and a † jδq via equation (45). Hence iφ † rσF (x) + iφ rσF (x) = jδ Λ jδ σr q>0