Electron Transport in Carbon Nanotubes with Adsorbed Chromium Impurities

We employ Green’s function method for describing multiband models with magnetic impurities and apply the formalism to the problem of chromium impurities adsorbed onto a carbon nanotube. Density functional theory is used to determine the bandstructure, which is then fit to a tight-binding model to allow for the subsequent Green’s function description. Electron–electron interactions, electron–phonon coupling, and disorder scattering are all taken into account (perturbatively) with a theory that involves a cluster extension of the coherent potential approximation. We show how increasing the cluster size produces more accurate results and how the final calculations converge as a function of the cluster size. We examine the spin-polarized electrical current on the nanotube generated by the magnetic impurities adsorbed onto the nanotube surface. The spin polarization increases with both increasing concentration of chromium impurities and with increasing magnetic field. Its origin arises from the strong electron correlations generated by the Cr impurities.


Introduction
The theory for disordered materials is less well developed than the theory for periodic systems. The simplest theory of a disordered system comes from the Born approximation to scattering theory for particles moving in the periodic system with isolated scatterers due to the disorder. However, this is a weak-coupling theory which works well only when the scatterer is similar to the host. This will certainly not be the case with transition metal or rare-earth impurities in conventional sp-metals. Alternatives, such as methods based on pseudopotentials [1] often fail because the nonlocal nature of the pseudopotential makes them difficult to transfer from one environment to another. This situation has changed, to some degree, with the introduction of Vanderbilt's ultrasoft pseudopotential [2,3] and the use of projector augmented waves in density functional theory, as proposed by Blohl [4,5]. This approach was further developed via the generalized gradient approximation (GGA) to density functional theory This calculation of real-time Green's functions of the disordered crystal is based on a diagrammatic technique, analogous to the technique used for a homogeneous system [32]. The set of equations of real-time Green's functions, along with expressions for the free energy and the electrical conductivity of disordered crystals, is based on the work presented in [31]. These methods are employed to obtain our final results. In the next two sections, we sketch the formalism to establish our notation and summarize the methods used.

Hamiltonian of an Electron-Phonon System for a Disordered Crystal
The many-body Hamiltonian of a disordered crystal consists of a single-particle Hamiltonian of the electrons in the external potential of the (disordered) ionic cores, the potential energy of the electron-electron interaction, the quadratic Hamiltonian for the phonons, the contribution from the electron-phonon interaction, and the anharmonic phonon potential terms.
We represent the Hamiltonian in the basis of free neutral atoms. In the Wannier representation, the system Hamiltonian is [31] H = H 0 + H int (1) where the zeroth-order Hamiltonian consists of the single-particle Hamiltonian of the electrons in the field of the ionic cores niγ,n i γ a + niγ a n i γ (3) and the harmonic phonon Hamiltonian for the motion of the ion cores niα,n i α u niα u n i α .
Here, the ion cores are located on a periodic lattice (i.e., the unperturbed system is periodically ordered and has no disorder). The symbol n denotes the unit cell, i denotes the ith basis vector in the nth unit cell, and γ denotes all of the other quantum numbers for the orbital, including spin. Disorder will enter for the species of ion at a particular lattice site, which need not be periodic via a perturbed Hamiltonian term (see below). The symbol h (0) denotes the "hopping integral" that connects the respective orbitals. For the phonon Hamiltonian, n and i are the same as before, namely the unit cell and basis site within the unit cell, while α is a spatial direction (x, y, or z). P is the ionic momentum, M is the ionic mass, u is the deviation of the ion from the equilibrium position of the lattice site, and Φ (0) is the corresponding spring-constant matrix.
The interaction Hamiltonian in Equation (1) is the perturbation of the system due to all of the effects we will be including. It is composed of six pieces: δΦ is the modification of the ion-core-ion-core Coulomb interaction due to the disordered ions added to the system; it is the difference between the original ion-ion repulsion Hamiltonian and the new one. The single-particle electronic Hamiltonian is modified by the change in the ion core and the extra term: H ei = ∑ niγ n i γ w niγ,n i γ a + niγ a n i γ , (6) which is the difference between the new hopping Hamiltonian and the original periodic one. The electron-phonon interaction is given by H eph = ∑ niγ n i γ v niγ,n i γ a + niγ a n i γ .
This is described in more detail below. The screened Coulomb interaction between electrons is given by the different multiband interaction terms, including density-density interactions and exchange interactions (2)niγ,n i γ n i γ ,n i γ a + niγ a + n i γ a n i γ a n i γ .
The modification of the interaction of the phonons with the impurity ion cores is given by niα,n i α P niα P n i α + 1 2 ∑ niα n i α ∆Φ niα,n i α u niα u n i α (9) where ∆M −1 niα,n i α = 1 M ni − 1 M i δ nn δ ii δ αα (10) ∆Φ niα,n i α = Φ niα,n i α − Φ (0) niα,n i α (11) and M ni and M i are the masses of the atoms at site (ni) for disordered and ordered alloys, respectively. We also include the cubic anharmonic potential terms for the phonons (under the assumption that they remain small and can be treated perturbatively) via niα,n i α ,n i α u niα u n i α u n i α .
The strategy for determining the matrix elements of the Hamiltonian is to employ the conventional density functional theory within the generalized gradient approximation to solve the Kohn-Sham equations in the presence of a Hartree plus exchange-correlation potential. We then employ those wavefunctions as the basis for expanding the full Hamiltonian, where we include explicitly the Coulomb repulsion of the electrons, which avoids the double-counting problem usually associated with trying to add correlations on top of density functional theory (which includes some correlation effects at the mean-field level).
More concretely, the strategy is as follows: We first determine the electronic wavefunctions from the standard density functional theory approach. The operators a + niγ , a niγ create and destroy electrons in the state described by Vane's function φ niγ (ξ) = ξ|niγ , where ξ = (r, σ) are the spatial and z-component of spin coordinates of the wavefunction [33]. These wavefunctions (of an electron in the field of a free neutral atom species λ located at site (ni) are obtained from the Kohn-Sham equation in density functional theory [10]: where γ is a superindex which incorporates the quantum numbers for the principle energy eigenvalueε, the standard angular momentum quantum numbers l and m, and the z-component of spin σ. To reduce the length of Equation (13), the relative coordinate (r − r ni ) is denoted by the simplification (r). In Equation (13), V λ ext (r) is the Hartree potential energy for an electron in the atomic core Coulomb field of type λ at basis site i, which is given by the simple expression In Equation (14), the electron density is summed over both spin components: n λi (r) = n λiσ (r) + n λi−σ (r) .
The electron density for each component of spin is given by where Z λ niγ is the occupation number electrons in the state denoted by γ (of type λ at basis site i and with spin component σ).
The exchange-correlation potential is more complicated. In the MGGA [8,10], it is represented by where is the GGA contribution, and µ XC,σ (r) =  (1) are denoted by ϕ niγ (r). They are found by solving the Kohn-Sham equations in Equation (13) and are expressed in a factorized form for the radial and angular components via where R λ iεlσ (|r − r ni |) is the radial component and we employ the spherical functions Y lm (θ, ϕ) for the angular component.
The potential-energy operator of an electron in the field of the different ionic cores is given by where r is the electron position vector, r ni = r n + ρ i is the position vector for the ion core at site (ni) in equilibrium, and u s ni ionic core's static displacement from equilibrium due to phonons, and u ni is the dynamic ion core displacement due to phonons. The total potentials of the ionic core v ni (r − r ni ) are found from Equation (20) and require a summation over all electronic states.
The matrix element of the electron-ion interaction Hamiltonian in Equation (6) is given by where w λn i niγ,n i γ = v λn i niγ,n i γ + ∆v λn i niγ,n i γ − v λ i n i niγ,n i γ (23) with λ i the type of ion at (n i ). Here, c λ ni are random numbers, taking the values 1 or 0, depending on whether the atom of type λ is at site (ni) or not. The symbol ν is a matrix element of the potential of the ionic core v ni (r − r ni ). The symbol ∆v will be defined next.
The expression for the electron-phonon interaction in Equation (7) is found through derivatives of the potential energy of the electrons in the ion core field due to a displacement of the atom by the vector u ni . In Equation (7), the value of v niγ,n i γ is given by v niγ,n i γ = ∑ n i α v n i α niγ,n i γ u n i α (24) where v with v λn i α niγ,n i γ the matrix elements of the following operator where e n i = r − r n i |r − r n i | .
The term ∆v λn i niγ,n i γ in Equation (21) describes electron scattering on the static displacement of the atoms and is defined by the equation The matrix of the force constants arising from the direct Coulomb interaction of the ionic cores has the form where Z ni is the valence of the ion cores at (ni). This matrix satisfies the following constraint: The force constants with the (0) superscript are defined in the same fashion, but correspond to the force constants of the initial periodic system with no disorder.
The matrix elements v (2)niγ,n i γ n i γ ,n i γ in Equation (8) are calculated by integrating over the corresponding angular variables. Integrals of the product of three spherical functions (a so-called Gaunt integral) are found by using Clebsch-Gordan coefficients [34,35]. This yields v (2)εlm,ε l m ε l m ,ε l m = e 2 ∑ |l − l | l 1 l + l |l − l | l 1 l + l l + l + l 1 = 2k, k ∈ N l + l + l 1 = 2k 1 , k 1 ∈ N 1 2l 1 + 1 × (2l 1 + 1)(2l + 1)(2l 1 + 1)(2l + 1) (2l + 1)(2l + 1)  where l and m are the standard angular momentum quantum numbers, c(l l l; m , m ) are the standard Clebsch-Gordan coefficients [34], R˜ε l (r) is the radial part of wave function, andε is the principle quantum number. There is a further simplification that we invoke when we treat the system with Gaussian orbitals. Thus, the matrix elements in the real wave function basis v (2)niγ,n i γ n i γ ,n i γ for different sites (ni) are approximately represented in a form similar to that in Equation (19). When the radial part is a Gaussian orbital, as is done in the method of molecular orbitals via linear combinations of the atomic orbitals 12 , multicenter integrals v (2)niγ,n i γ n i γ ,n i γ have the form of single-center integrals, because the product of two Gaussian orbitals that are localized at different centers can be reduced to the product of orbitals that are localized about a common center.

Green's Functions of Electrons and Phonons
We employ a Green's function-based formalism to perform the calculations. Ultimately, we need the real-time retarded G AB r (t, t ) and advanced G AB a (t, t ) Green's functions, which are each defined as follows [35]: Here, the operators are expressed in the Heisenberg representatioñ whereh is Planck's constant, H = H − µ e N e , µ e is the chemical potential of the electronic subsystem, and N e is the electron number operator given by In addition, the commutator or anticommutator is defined via where the commutator is used for Bose operators (-), and the anticommutator is used for Fermi operators (+). The symbol θ(t) is Heaviside's unit step function. The angle brackets ... denote the thermal averaging with respect to the density matrix ρ where Ω is the thermodynamic potential of the system given by exp(Ω/Θ) = Tr exp(−H/Θ) and Θ = k b T, with k b Boltzmann's constant and T the temperature. Note that, even though the real-time Green's functions appear to depend on two different times, because of the time-translational invariance for equilibrium systems, they actually depend only on the time difference t − t . Our procedure for calculating the real-time Green's functions follows the standard one-we first determine the thermal Green's functions (defined below) and then analytically continue using them as real-time functions using the conventional spectral relations.
In addition, the time-ordering operator satisfies where the plus sign is used for Bose operators and the minus sign for Fermi operators. We next go to the interaction representation by introducing the operator with Differentiating the expression for σ(τ) in Equation (40) with respect to τ and then integrating from 0 with the boundary conditionσ(0) = 1, we obtain where H int (τ) = e H 0 τ H int e −H 0 τ . Employing this result yields with A(τ) in the Heisenberg representation with respect to the noninteracting Hamiltonian. Substituting these results into the definition of the thermal Green's function creates the alternate interaction-representation form for the Green's function, given by where all time dependence is with respect to the noninteracting Hamiltonian and the trace over all states is with respect to the noninteracting states This last result forms the starting point for the perturbative expansion employed here. The diagrammatic method is generated by expanding σ(τ) in a power series in terms of H int (τ) and then using Wick's theorem to evaluate the resulting operator averages (since the noninteracting Hamiltonian is quadratic [31]). This technique then generalizes the approach used for the homogeneous system [31]. The denominator in Equation (43) cancels all disconnected diagrams in the expansion, as usual. Therefore, the thermal Green's function are expanded in terms of connected diagrams. Using the standard relations between the spectral representations of thermal and real-time Green's functions [31], we obtain the following Dyson equation for the electronic Green's function in the real frequency domain (hereinafter the dependence on r is suppressed) [31]: , and G Pu (ε) are the real-frequency representations of the single-particle Green's function of the electrons, the coordinate-coordinate, momentum-momentum, coordinate-momentum, and momentum-coordinate Green's functions of the phonons, respectively; Σ eph (ε), Σ phe (ε), Σ ee (ε), and Σ phph (ε) are the corresponding self-energies (mass operators) for the electron-phonon, phonon-electron, electron-electron, and phonon-phonon interactions. The real-time and real-frequency Green's functions are related by standard Fourier transform relations given by and The thermal Green's functions are periodic (bosons) or antiperiodic (fermions) on the interval −1/Θ τ < 1/Θ, and hence have a Fourier series representation in terms of their Matsubara frequencies, as follows: and where the Matsubara frequencies satisfy ω n = 2nπ Θ for Bose particles (2n + 1)π Θ for Fermi particles n = 0, ±1, ±2, ... .
Note that the thermal Green's functions are equal to the retarded Green's functions evaluated at the Matsubara frequencies for positive Matusbara frequencies and are equal to the advanced Green's functions evaluated at the Matsubara frequencies for negative Matsubara frequencies.
The electronic Green's functions are infinite matrices with indices given by the lattice site n, the basis site i, and the other quantum numbers γ. Similarly, the phonon Green's functions are also infinite matrices with the same lattice and basis site dependence plus a dependence on the spatial coordinate direction α. This produces some simple equations for the noninteracting Green's functions, namely [31] and Here, the double lines denote a matrix. When the perturbations are small, given by then the solution of the system of equations in Equation (45) becomes where These solutions only include the first-order corrections due to the small terms in Equation (56). Ref. [31] summarizes the explicit formulas for the corresponding self-energies, which we do not reproduce here.
The electronic self-energy due to the electron-phonon interaction Σ eph (τ, τ ) is described by the diagram in Figure 1. The solid lines correspond to electronic propagators G aa + niγ,n i γ (τ, τ ) and the dashed lines correspond to phonon propagators G uu niα,n i α (τ, τ ). The vertex correction Γ n 2 i 2 α 2 niγ,n 1 i 1 γ 1 (τ 2 , τ, τ 1 ) is given by the diagrams in Figure 2. Note that the unshaded triangle corresponds to the equation In Figures 1 and 2, the internal summations forñγ imply both a summation over niγ and an integration over the internal time τ. Each diagram has an overall sign determined by (−1) n+F , where N is the order of the diagram (number of vertices Γ 0 ), and F is the number of electronic Green's function lines G aa + . Explicitly, the electron-phonon self-energy becomes where repeated indices are summed over.
The self-energy of the phonon due to the phonon-electron interaction is shown in Figure 3. Evaluating the diagram using the above rules yields f (ε) is the so-called Fermi-Dirac distribution function. The electron-electron scattering contribution to the electronic self-energy Σ ee (τ, τ ) is shown in Figure 4, while the vertex part Γ n 2 i 2 γ 2 ,n 1 i 1 γ 1 niγ,n i γ (τ 2 , τ 1 τ, τ ) is given in Figure 5.
Note that the unshaded square in Figure 5 corresponds to the equation Using this vertex function, then yields the contribution to the electron self-energy from the electron-electron interaction: A similar result for the contribution to the phonon self-energy Σ phph (ε) from phonon-phonon coupling is given in [31]. In deriving the expressions in Equations (60), (61), and (64), we employed the standard resummation techniques for any function ϕ(z) that is analytic in the region covered by the contour C, which encloses all of the Matsubara frequencies. Namely, we have for the Bosonic case, and for the Fermionic case, withf We comment that, for the many-body Green's functions described here, it is customary to have the chemical potential situated at zero frequency, as we do here.
The Fermi level ε F ≡ µ e of the system is determined by the equation where Z is the average number of electrons per atom, and g e (ε) is the many-body electronic density of states, which satisfies Here, ... c denotes configurational averaging over the disorder, N is the number of primitive lattice cells, and ν is the number of atoms per primitive cell. We drop the letter c on the configurational averaging for simplicity. In Equation (71), Z is the average number of electrons per atom.
It should be noted that the first term in the electron self-energy due to electron-electron interactions, Σ ee niγ,n i γ (ε), which is caused by corrections beyond Hartree-Fock, describes the effects of electron correlations. As opposed to the procedures used in [13][14][15][16][17][18][19], the long-range Coulomb interaction of electrons located at different lattices sites of the crystal is described by taking into account an arbitrary number of energy bands.
The expression for the Green's function in Equation (54) differs from the corresponding expressions for the Green's function of a single-particle Hamiltonian of a disordered system only from the different self-energy contributions. Hence, we solve for the Green's function using the well-known methods of disordered systems theory [28].

Localized Magnetic Moments
Since we will be working with magnetic moments for the remainder of the paper, we now slightly modify our notation so that the symbol γ now refers to all other quantum numbers except for spin, and we will introduce the spin quantum number σ explicitly in all of the following equations. We will be employing the approximate results for the electron-electron self-energy in a self-consistent fashion to allow correlation effects to modify the electronic bands. This requires a heterogeneous distribution of the electron density. We assume that this electron density distribution corresponds to the minimum of the free energy. The electron-electron self-energy in Equation (64) (66), where the total electronic density of states g e (ε) is replaced by the partial density of states g λm λi niγσ (ε) for energy band γ and spin projection σ to allow for magnetic solutions. The occupation numbers Z λm λi niγσ and the partial density of states g λm λi niγσ (ε) then satisfy Note that the disorder averaging is done under the assumption that an atom of type λ is located at the site (ni), and its projection of the localized magnetic moment onto the z-axis is equal to m λi . The probability of this configuration is P In this fashion, we allow for localized magnetic moments which are inhomogeneously distributed through the crystal lattice and correspond to static magnetization fluctuations.
The total charge and magnetization for each orbital on a site are given by and by respectively. We need to sum over all other quantum numbers to get the totals: Next, we calculate the phonon Green's function in the coordinate basis G uu (ε) by solving Equation (57). Here, we employ the following procedure to take into account the heterogeneity: First, we work with the homogeneous system, which is pure and has no disorder. Then we introduce the disorder and compute its effects via a cluster expansion related to cluster expansions from the theory of alloys. Therefore, the zeroth approximation for the Green's function is the Green's function of the pure system given by G aa + (ε), which we call the effective medium Green's function. Since the system is pure for the effective medium, we compute the effective medium Green's function via a Fourier transformatioñ where N is the number of primitive unit cells, and σ e (k, ε) is the self-energy of the effective medium, also called the coherent potential. The coherent potentials are determined via the coherent potential approximation, which is described in detail below, with the coherent potential given in Equation (90). We do a similar procedure for the effective medium phonon Green's function, which satisfies where we have and the coherent potential satisfies Equation (91). Note that the wavevector k varies within the first Brillouin zone. Furthermore,Σ eph (k, ε) is the Fourier transformation of the matrix Σ eph niγ,n i γ (ε) given in Equation (60) for which the terms v n 1 i 1 α 1 niγ, n 3 i 3 γ 3 are replaced by the values for a pure crystal and the corresponding Green's functions are those of the effective medium. The other self-energies given byΣ ee (k, ε),Σ phe (k, ε), andΣ phph (k, ε) are defined similarly.
In Equation (82), Φ (0) (k) is the Fourier transform of the matrix Φ (0) niα,n i α , which describes the atomic nucleus repulsion. The self-energyΣ phe (k, ε) describes the attractive interaction between the atomic nuclei and the electrons.
The Green's functions in Equation (57)

satisfy a Dyson equation that can be expressed in terms of a T-matrix via
where the T-matrix T is represented by a series, in which each term describes the scattering of clusters with different numbers of nodes expressed schematically by Here, we have where t n 1 i 1 is the on-site scattering operator given by The self-energy employed in Equation (86), Σ n 1 i 1 e (ε), satisfies for the electrons. For the phonons, we have The coherent potential approximation requires that which yields a system of coupled equations [31] for the electronic and phononic coherent potentials. Now we are ready to tackle the effects of disorder. Starting from Equation (83), we can obtain a cluster decomposition for the Green's function of electrons and phonons. The electronic and phononic density of states, the free energy, and the electrical conductivity can all be expanded in an infinite series. Each term in the series describes the scattering on clusters with different numbers of atoms. It turns out that the strength of the scattering within a cluster decreases with an increasing number of atoms in the cluster and can be represented by the following small parameter: where r is the total number of energy bands included in the calculation. We have shown previously [36,37] that this parameter remains small when many parameters of the system are changed, except possibly for narrow energy intervals near the band edges. Next, we employ Equations (72), (74), and (83) to perform an average over the distribution of different types of atoms and different projections of the localized magnetic moments on the sites of the crystal lattice. During the averaging process, we neglect the contribution of electron scattering in clusters consisting of three or more atoms, since they are guaranteed to be small due to the smallness of the parameter in Equation (92). After performing the averaging, the electronic density of states becomes Similarly, averaging of the phonon Green's function G uu (ε) yields the phononic density of states: , whereG =G uu (ε). In these formulas, the single-center scattering operator t λn 1 i 1 is given by Equation (80). According to Equations (6), (60), (64), and (87), the self-energy denoted by Σ λn 1 i 1 e (ε) describes the electron scattering: The value ofZ n 1 i 1 γ 1 σ 1 ,n 2 i 2 γ 2 σ 2 in Equation (95) In Equation (93), is the conditional probability of finding an atom of type λ at site (lj) with magnetic moment m λ j , provided that the sites in the unit cell at the origin (0i) have an atom of type λ with a magnetic moment m λi . Here, t λ m λi ni is the value of the matrix element of a single-center operator for scattering in the case where an atom of type λ is located at site (ni) and has a magnetic moment m λi .
When the system is disordered, we need to consider a random arrangement of the disordered atomic sites. Hence, in Equation (94), the probability of an atom of type λ to be at site (0i) is given by where c λ ni is a discrete binary random number taking the values of 1 or 0, depending on whether an atom of type λ is at site (ni) or not, respectively. The joint probabilities in Equations (86) and (94) are defined by the following: The probabilities are determined by the interatomic pair correlations ε B B lj 0i , ε [23,27] where δ is the Kronecker delta function. Note that the interatomic pair correlations also satisfy The notations P m λi 0i , P m λ j /m λi lj 0i denote the probabilities of the static fluctuations of the magnetization.
As an example, when we have a binary alloy, consisting of two sublattices, and two types of atoms A and B, we obtain for the first sublattice and for the second sublattice, with Here, we have ν = ν 1 + ν 2 is the total number of (given by more atoms of one type on Sublattice 1 and vice versa on Sublattice 2).
We assume that the projections of the localized magnetic moment onto the z axis are given by two values m λi = µ + λi , µ − λi . The probability P m λi 0i is connected with the long-range magnetic parameter η m via the expressions for Sublattice 1 and P Here, are equal to the relative number of lattice sites with localized magnetic moment projections µ + λi and µ − λi , respectively. The value x µ + λ = x µ − λ = 0.5 when the external magnetic field vanishes H = 0, corresponding to a paramagnetic state.
So far, we have described how one performs one iteration of the self-consistent calculation. Once one iteration has been completed, the values of the occupation numbers of the electron states in Equation (16) are determined by Equation (73): The iterations continue until the densities of states have converged. Once convergence has been reached, we can then employ the Green's functions to calculate observables, which is described in the next two sections.

Free Energy
We first focus on the Gibbs free energy (also called the thermodynamic potential) of the system which satisfies [32]: The Hamiltonian H is defined in Equation (1). To perform the trace, we need to sum over all of the band states, but we also need to take into account the disorder averaging. The latter is commonly handled via a configurational average [31]. As for the correlations, we employ Equation (40) for the interaction picture, which allows us to re-express the Gibbs free energy as where Ω (0) e and Ω (0) ph are the thermodynamic potentials for the electrons and the phonons in the field of the ionic cores, respectively. As before, the equilibrium ion core positions are chosen to be the same as those of the crystal lattice for the pure ordered crystal, even when we introduce disorder and change the type of some of the atoms. The symbol Ω is the component of the thermodynamic potential that is caused by the mutual scattering of electrons and phonons; it is defined by with σ given in Equation (40) for the interaction picture.
In addition, S c = − ln P c is the configurational entropy, where P c denotes the distribution function for atoms with a specific z-component of the magnetic moment on a given lattice site. The angular brackets ... denote the configurational averaging over different disorder configurations for a given density of disorder.
Next, we use the "integration over the coupling constant" method to simplify the results further. By replacing the interacting Hamiltonian H int (defined in Equation (5)) by H int (λ) = λH int , differentiating the expression for the piece of the thermodynamic potential Ω (λ) in Equation (111) with respect to λ, and then integrating them (with the boundary conditions Ω (0) = 0, Ω (1) = Ω ), we obtain the following after a long derivation: This expression can be immediately evaluated, because we know all the Green's functions and self-energies.
The contribution to the thermodynamic potential from the electrons (in the field of the ionic cores) is also simple to find. It is given by Similarly, the contribution to the thermodynamic potential from the phonons (in the field of the ionic cores) is given by The densities of states g ph (ε) in Equations (113) and (114) are given by the results in Equations (93) and (94); note that the Green's functions G aa + (ε) and G uu (ε) are replaced by Green's functions G aa + 0 (ε) and G uu 0 (ε) for the zeroth-order approximation.
Finally, the configurational entropy can be represented as [31] S c = − ∑ λ , m λi ,ni P λm λi ni ln P λm λi ni Ultimately, we are interested in determining the Helmholz free energy, F, as a function of the volume V, the temperature T, the number of electrons N e , and the parameters of the interatomic correlations (ε n 1 i 1 ,n 2 i 2 , η). The Helmholz free energy can be found directly from the thermodynamic potential. Namely, it satisfies F = Ω + µ e N e . The free energy per atom can be approximated by [31] where Ω e , Ω ph are given by Equations (113) and (114), but with g ph (ε) replaced by g e (ε) and g ph (ε) (see Equations (93) and (94)), in the situation where the electron scattering is weak.

Electrical Conductivity
In this section, we discuss how to calculate the electrical conductivity. We assume the system will not be driven too far from equilibrium, so we employ the linear response formalism of Kubo for the electrical conductivity tensor [38,39], which is given by In this equation, J α is the current operator along the α spatial direction. The real part of the conductivity, called the optical conductivity, can then be represented in terms of the imaginary part of the retarded response function, or equivalently as in terms of the retarded and advanced response functions. The current operator is just the number operator for the electrons, multiplied by their velocity and the electric charge, and then summed over all states. It is compactly represented via where Ψ + (ξ, t) and Ψ (ξ, t) are the field operators for the creation and annihilation of electrons, respectively, ν α is the operator of the α component of the band velocity, and e is the electron charge. The integration over ξ sums over all states.
To get the retarded response function on the real frequency axis, we must analytically continue the thermal response functions. The thermal current-current response function is defined to be where V 1 is the volume of the primitive unit cell and the two-particle thermal Green's function is given by the following time-ordered expectation value: G n 1 τ , n 2 τ, n 3 τ , n 4 τ = T τ a n 1 τ a n 2 (τ) a + n 3 τ a + n 4 (τ) (121) The two-particle Green's function consists of a bare direct, a bare exchange, and a vertex-corrected piece, which are illustrated schematically in Figure 6. The numbers at the vertices are a shortcut for all of the relevant quantum numbers and imaginary time, e.g., 1 corresponds to (n 1 i 1 γ 1 τ 1 ).
n ≡ niγσ. The electron velocity satisfies the conventional definition The interaction piece is evaluated approximately via ∆G I I αβ (ε 1 ; ε 2 ) ≈ ∆G I I αβ (ε 1 ; ε 2 ), where ∆G I I αβ (ε 1 ; ε 2 ) is given in Equation (124) but replacing G aa + (ε) withG aa + (ε). The above derivation yields the dc conductivity in the presence of a static electric field. In this work, we are also interested in the effects of a weak external magnetic field. Within a nonrelativistic approximation, an external magnetic field is introduced into the kinetic energy of an electron via the matrix element h (0) n 1 i 1 γ 1 ,n 2 i 2 γ 2 (see Equation (3)). This is done using the minimal substitutionp 2 /2m → p − e c A 2 /2m, with A the vector potential of the electromagnetic field and c the speed of light. In addition, we need to include a term in the Hamiltonian that corresponds to the interaction energy of an intrinsic magnetic moment with the external magnetic field: where µ B is the Bohr magneton, H is the external magnetic field, and σ is the projection of the magnetic moment (spin) onto the direction of the magnetic field, and m is the angular momentum quantum number.
Since we are only interested in a weak external magnetic field, we only include the spin projection term [in Equation (126)] in the Hamiltonian when a magnetic field is present. By properly treating the local magnetic moment (spin) dependence, we can then compute spin-dependent transport. Note that, since these calculations are approximate, the accuracy is determined both by the strength of the vertex corrections that are ignored and by the smallness of the parameter that governs the cluster expansion.

Spin-Dependent Transport of Carbon Nanotubes with Chromium Atoms
In this section, we present the application of the above formalism to the problem of chromium impurities doped onto a carbon nanotube. In these final calculations, we neglect the vertex corrections in Equations (60), (61), and (64), and we neglect the static displacements of the atoms. We employ a 2s, 2p-states wave-function basis for the neutral carbon atoms and a 3d, 4s-states wave-function basis for neutral chromium atoms. The initial ion core valence of C and Cr atoms Z λ i are 4 and 6, respectively. In addition, to reduce the total computational time, we only performed one iteration of the self-consistent iterative procedure. In Equation (16), we used Z λ iδσ = 1 for the occupied electronic states. The off-diagonal matrix elements of the Hamiltonian in Equation (1) were calculated by including the first three coordination spheres. The Green's function calculation in Equations (79) and (81) employed 10 3 points in the Brillouin zone. All calculations were performed at T = 300 K.
The calculations start with a (3,0) chirality carbon nanotube that is doped with Cr atoms. The geometry for the crystal structure is then optimized by minimizing the free energy F, defined in Equation (116). Note that the carbon nanotube (doped with Cr) has a one-dimensional crystal structure. The primitive cell contains 18 nonequivalent atomic positions. Carbon atoms are located at 12 positions on the surface of the inner cylinder. The distance between the carbon atoms is 0.142 nm. Cr atoms are located at the 6 positions on the outer surface of the cylinder opposite the center of a hexagon, the vertices of which are carbon atoms. The distance between carbon atoms and neighboring Cr atoms is 0.22 nm. a cross-sectional view of the crystal structure of the (3,0) chirality carbon nanotube with adsorbed Cr impurities is plotted in Figure 7.
It turns out that the free energy is minimized by a random arrangement of Cr atoms on the surface of the nanotube. Figure 8 plots the dependence of the free energy on the pair correlations of Cr impurities with ε BB = ε BB lj0i in Equation (100) and for the first coordination sphere. The symbol B denotes an atom of Cr. The dependence of the free energy on the pair correlations is shown in the region of the free energy minimum.  As shown in Figure 8, the free energy F has its minimum at ε BB = 0. This implies that the Cr atoms are randomly located on the surface of the carbon nanotube (for this density of impurities). The relative positions of the carbon atoms and the chromium impurities are similar to those found elsewhere for transition-metal dopants on carbon nanotubes using ultrasoft pseudopotentials [40,41]. The value of the localized magneticmoment of a Cr impurity and the induced localized magnetic moment of a C atom in the direction of the magnetic field increases, as expected, with the size of the field. For carbon nanotubes that have five Cr atoms per primitive unit cell, the projection of the Cr magnetic moment varies in the range m Cr = (1.02; 2.24)µ B , while the induced C magnetic moment lies in the range m C = (0.0036; 0.02)µ B as the external magnetic field increases from zero to H = 200 A/m. In this calculation, the magnetic field is oriented along the axis of the carbon nanotube. We also find that the pair correlations for the orientation of the localized magnetic moments in the first coordination sphere satisfies ε m = 0.235 when the external field vanishes. This pair correlation nearly vanishes for the second and third coordination spheres. A positive value of ε m for the first coordination sphere indicates that the induced magnetic moment on a C atom is oriented in the same direction as the magnetic moment of the nearest Cr atom, as one might have expected. Figure 9 plots the partial g eσ (ε) = 1 v ∑ i,γ,λ P λ 0i g λ 0iγσ (ε) and total g e (ε) = ∑ σ g eσ (ε) densities of states for the electrons on a carbon nanotube with an adsorption of Cr (and vanishing external magnetic field). In this case, we have a paramagnetic phase, so g 1/2 (ε) = g −1/2 (ε). The vertical line shows the Fermi level ε F . Figure 10 plots the partial g eσ (ε) and total g e (ε) densities of states for a carbon nanotube with five atoms of Cr per primitive unit cell and in an external magnetic field of strength H = 100 A/m and oriented along the tube axis. The plot is zoomed into the region near the Fermi level.
As shown in Figure 10, the partial density of states g eσ (ε) for spin σ = 1/2 is shifted relative to those for spin σ = −1/2. These results are also qualitatively consistent with results obtained by a different method in [40]. However, the quantitative results differ because this work examines a different chirality for the nanotube than [40] does ((3,0) versus (9,0)).
In Figure 11, the dependence of the spin polarized electrical conductivity ∆σ/σ = (σ 1/2 − σ −1/2 )/σ of a (3,0) chirality carbon nanotube with five atoms of Cr per primitive unit cell versus the magnitude of the external magnetic field is plotted for T = 300 K.

Conclusions
In this work, previous methods [6][7][8][9][10][13][14][15][16][17][18][19] that describe pure ordered crystals and molecules have been generalized to include disorder effects. The method employed involves a diagrammatic expansion for the electron correlations (under the assumption they are small) along with a cluster-based method to treat the disorder effects (truncated to a small cluster). This method employs Green's functions but is rooted in density functional theory.
The theory is applied to a particular case of Cr dopants added to a carbon nanotube. In particular, a (3,0) chiral nanotube has five Cr atoms per primitive unit cell added to the system. We find that the resulting spin-dependent electron transport derives from strong electron correlations caused by the presence of chromium atoms. The magnitude of the spin polarized current stems primarily from the difference of the partial densities of states (see Figure 10) with opposite spin projection at the Fermi level. However, it is also affected by the difference between the relaxation times that arise from the different occupation numbers of the single-electron states Z λ niγσ of C and Cr (see Equation (95)). The spin polarization of the electric current increases with the concentration of Cr atoms and with the magnitude of the external magnetic field. The results presented here generically agree with those calculated with ultra-soft pseudopotentials [41,42]. The main difference is that, in the present work, we do not find that a gap opens up when the spin polarization becomes large enough.
Author Contributions: S.R. did the calculations.; I.V. constructed and optimized the computational models; Y.N. project administration; S.K. responsible of methodology and conceptualization; S.B. supervised the work and did main calculation. All authors contributed to manuscript preparation.
Funding: S.B. acknowledges support from the NATO Partnership and Cooperative Security Committee, in the framework of the Science for Peace and Security (SPS) Programme, project SPS G5351-"Nanocomposite Based Photonic Crystal Sensors of Biological and Chemical Agents." S.K. acknowledges the support from the National Academy of Sciences of Ukraine (project No. 0116U002067).