Spin-dependent transport through interacting graphene armchair nanoribbons

We investigate spin effects in transport across fully interacting, finite size graphene armchair nanoribbons (ACNs) contacted to collinearly spin-polarized leads. In such systems, the presence of short ranged Coulomb interaction between bulk states and states localized at the ribbon ends leads to novel spin-dependent phenomena. Specifically, the total spin of the low energy many-body states is conserved during tunneling but that of the bulk and end states is not. As a consequence, in the single-electron regime, dominated by Coulomb blockade phenomena, we find pronounced negative differential conductance features for ACNs contacted to parallel polarized leads. These features are however absent for an anti-parallel contact configuration, which in turn leads within a certain gate and bias voltage region to a negative tunneling magneto-resistance. Moreover, we analyze the changes in the transport characteristics under the influence of an external magnetic field.

. A graphene armchair nanostripe contacted by ferromagnetic leads. At the long sides, the lattice is terminated in armchair, at the small ends in zig-zag configuration. The length of a bond between two carbon atoms is a 0 ≈ 0.14 nm. We choose the orientation of the coordinate system such that the x-axis points along the zig-zag ends, the y-axis along the long armchair edges of the stripe.

Introduction
Since the first successful separation of a one atom thick sheet of graphite by Novoselov and Geim and coworkers [1], graphene, i.e. an isolated single sheet of graphite, has been a material attracting ever raising interest. Not only a great potential for applications [2,3], but also fundamental physics issues [4] arise from the linear dispersion relation in the electronic band structure of graphene, predicted about sixty years ago [5], as confirmed by various recent experimental findings [1,6,7,8]. Increasing effort is presently put in the understanding of the electronic properties of graphene nanodevices, which can be obtained by etching or lithographic techniques and may achieve lateral dimensions of a few tenth of nm [9,10,11]. Studies on the effects of electron-electron interactions and confinement in transport across graphene nanodevices have been carried out on single-dots [10,12] and, only recently, on double-dot structures [13,14]. A desirable goal is the fabrication of clearly defined geometries, and of particular interest for applications [2] could be narrow stripes of graphene, so-called carbon nanoribbons. Conductance quantization has been observed in 30nm wide ribbons [11], while an energy gap near the charge neutrality point scaling with the inverse ribbon width was reported in [9] and theoretically [15,16] attributed to Coulomb interaction effects.
Crucial for the properties of a nanoribbon is the form of its ends. The two most regular possibilities are an armchair and a zig-zag end (see Fig. 1). At present, the shape of the ends cannot be controlled, but there is ongoing progress in developing methods to fabricate stripes with clear edges, by scanning tunneling microscope lithography [17], chemical synthesis [18] or by unrolling carbon nanotubes [19]. Moreover, there exist theoretical studies [20,21] claiming that any narrow stripe should show either the behavior of a zig-zag or of an armchair ribbon (ACN), where the names specify the form of the long side edges. A peculiar property of the zig-zag edge is the existence of localized states [22] which were indeed observed experimentally by means of scanning tunneling microscopy [23,24,25]. Due to their high degeneracy, flat-band ferromagnetism is expected from the Hubbard model, leading to a spin-polarized many particle ground state. It was suggested [26,27] to exploit this property for spintronic applications, where transport is governed by carriers in the oppositely polarized channels along the two long zig-zag edges. Recently both Hubbard and long-ranged Coulomb interaction effects have been analyzed [28] under the assumption of a filled valence and an empty conduction band (half-filling). There was a prediction of strong spin features in case of a low population of the localized mid-gap states. In contrast, in armchair ribbons the localized states are at the far apart zig-zag ends. As we showed recently [29], for narrow armchair ribbons short-ranged Coulomb interaction, i.e. local scattering events between either the two sublattices of graphene or between the extended bulk and the localized end states, gain increasing weight. They lead to an exchange coupling, inducing entanglement between the bulk and the end states. This has a decisive impact on eigenstates and transport properties.
In this work we extend the investigations of Ref. [29] to the case in which the ACN is contacted to ferromagnetic leads and/or additionally subjected to a magnetic field. Various remarkable features due to the presence of the entangled end-bulk states, e.g. a negative tunneling magneto-resistance in a fully symmetric set-up, are discussed. The paper is structured as follows: In Sec. 2 and Sec. 3, we present our low energy theory for narrow ACNs. Readers not interested in technical details can, on the basis of the results for the ACN Hamiltonian summarized in Sec. 3.3, directly start from Sec. 4, where we shortly analyze a minimal set of states relevant for the explanation of features in transport across ACNs. Specifically, we show that scattering between end and bulk electrons causes an entanglement of states with the same total spin-projection S z , but different configurations of end and bulk spins. As a consequence, neither for the spins trapped at the ribbon ends, nor for the bulk electrons, S z is a good quantum number any more, which could be a handicap for proposed quantum information applications [2]. Although the formation of states symmetric or antisymmetric under the exchange of end spins remains without a special signature in the spectrum, it in fact leaves strong fingerprints in transport, as discussed in Ref. [29] for an unpolarized set-up as well as in Sec. 5 for collinear contact magnetizations. The most prominent feature in spindependent transport for collinear lead magnetizations is a negative tunneling magnetoresistance (TMR) within a narrow region along the edges of the Coulomb diamonds for even fillings, Sec. 5.2. Finally we explain how the transport characteristic is expected to change under application of an external magnetic field, both for non-magnetic and collinearly polarized contacts, Sec. 6.

Noninteracting electrons in a finite size graphene armchair nanoribbon
Being crucial for the subsequent analysis, we review here the end and bulk states of narrow, noninteracting ACNs. The carbon atoms in the graphene lattice are arranged in a hexagonal honeycomb lattice. There are two atoms per unit cell such that we can define two different sublattices p = ± as shown in Fig. 2. Hybridization of the 2s-orbital with the 2p x -and 2p y -orbitals leads to strong σ-bonds in the lattice plane. The electrons in the remaining 2p z orbitals form π-bands which determine the electronic properties at low energies. Characteristic for the structure of the π-band are the valence and conduction bands that touch at the corner points of the first Brillouin zone, also called Dirac points. Since there is one p z -electron per carbon atom, in isolated graphene the valence band is completely filled whereas the conduction band is empty. In the vicinity of the Dirac points the band structure exhibits a linear dispersion relation, Fig. 3a, resembling, up to a reduced propagation velocity of v = 8.1 · 10 5 m/s, the one of massless relativistic particles. From now on we focus on the region of linear dispersion. A description in terms of a Dirac equation for the p z -electrons capturing the essential features of the π-band close to the Dirac points can be obtained by a nearest neighbor tight binding calculation [30]. The six corner points of the first Brillouin zone can be decomposed into two subsets of equivalent Dirac points. As particular representatives we choose where a 0 ≈ 0.14 nm is the nearest neighbor distance. Restricting the discussion to the vicinity of K ± , the π-electrons are described by Bloch waves where r, R ∈ R 2 , N L is the number of sites of the considered lattice, α = ± denotes the conduction and valence band, respectively, and χ R p ( r) is the p z orbital on sublattice p at lattice site R. Finally, κ is the wave vector relative to the Dirac point K F . Defining the spinors it is found that they fulfill the Dirac equation where reflects the linear dispersion relation and σ x , σ y are the Pauli matrices. From (4) it follows that for κ y = ±iκ x it holds the relation A specific solution of Eq. (6) we will use in the remaining of the article is given by

Boundary effects in narrow ribbons
So far no boundary effects have been included. However, we wish to discuss the electronic properties of finite size carbon nanoribbons, which implies a wave function which vanishes all along the ribbon edges (open boundary conditions). The geometry we want to study is depicted in Fig. 1. In particular we assume the long edges of the ribbon, along the y direction, in an armchair configuration, while the narrow terminations of the ribbon, along the x axis, have zig-zag character. We are interested in quasi onedimensional ribbons and thus we restrict our discussion to stripes with a large aspect ratio L y ≫ L x where L x and L y are the extensions in x and y direction, respectively. From Fig. 2 we can easily see that the two zig-zag ends each consist of atoms either sitting on sublattice p = + or p = −. We will use the convention that the "left" end at y = 0 is formed by atoms living on sublattice p = − whereas on the other end (at y = L y ) we have atoms from sublattice p = + only. To fulfill the appropriate boundary condition for the zig-zag ends, we have to construct new wave functionsφ F α ( r, κ) from linear combinations of ϕ F α ( r, κ), in such a way that they vanish on the "missing" atoms at the ends, namely on sublattice p = + on the left end and p = − on the other end. A lengthy but straightforward calculation leads tõ with a normalization constant ‡ C zz (F κ x , κ y ) ∈ C and the quantization condition [31] Additionally, the wave function we are looking for must also vanish at the armchair edges. In contrast to the zig-zag end, the terminating atoms where the wave function ‡ The normalization constant is non-trivial as the functions given in Eq. (2) are non-orthogonal. To verify that there is the dependence on F κ x and κ y , but not on α, Eq. (2) must be consulted in combination with Eqs. (1) and (6).
is required to go to zero are from both sublattices. In order to build up suited linear combinations,we have to mix states of Eq. (7a) which belong to different Dirac points. Then, the resulting wave function ϕ α ( r, κ) vanishes on the lattice sites with R x = 0 and with C ac (κ x , κ y ) ∈ C another normalization constant, and the quantization [31] , the two conditions in Eq. (8b) are equivalent.

The eigenstates of metallic ACNs
In total we obtain with Eqs. (7a) and (8a) the following expression for the eigenstates of noninteracting electrons in finite size ACNs, where the generally complex number C(κ x , κ y ) = C zz (κ x , κ y )C ac (κ x , κ y ) guarantees that ϕ α ( r, (κ x , κ y )) is normalized to 1.
We want to investigate now the solutions of the Dirac equation fulfilling our quantization conditions Eqs. (7b) and (8b). It has been shown [22,20] that the presence of zig-zag ends leads to the formation of localized states characterized by a purely imaginary κ y giving rise to an enhanced density of states around the Dirac energy. Though those states are localized at the zig-zag ends with an exponential decay in y direction, they are of decisive relevance for the transport properties of ACNs.
At first, however, let us focus on the bulk states, where both κ x and κ y are real numbers. Since L y ≫ L x , the quantization condition Eq. (8b) leads to the formation of sub-bands assigned to different κ x . For a metallic ACN, the sub-band must cross the Dirac points F = K ± , which demands κ x = 0, or equivalently n x = Lx π K 0 , n x ∈ N, which follows from Eq. (8b). The width L x depends on the number M of hexagons in a row parallel tox like L x = √ 3a 0 (M + 1/2) (compare Fig. 1). Hence n x = 1 3 (4M + 2) (cf. Eq. (1)). Obviously, this means that the geometrical condition for having metallic ACNs with gap-less sub-bands reads M mod 3 = 1. Focussing on such ribbons, κ ∝k y as κ x = 0, which corresponds to a cut through the Dirac cone as shown in Fig. 3a. The corresponding states are characterized by the band index α and κ y , see Eq. (8a). With κ x = 0, Eq. (7b) yields as allowed values of κ y : For solutions with imaginary κ y the cone opens along ik y with a slope determined by v. As for the eigenstates it must hold κ y = ±iκ x , the corresponding dispersion relation is obtained from the intersection of the cone with the two planes k x = ±ik y . There is only tangency along two straight lines within the k x -ik y -plane, resulting in a dispersion which is identically zero.
Since by definition, Eq. (9), ϕ α ( r, (0, κ y )) = −ϕ α ( r, (0, −κ y )), we can further restrict for each α our analysis to either κ y > 0 or κ y < 0. Thus it is an allowed choice to just consider states with sgn(κ y ) = sgn(α), which we define as ϕ b κy ( r) := ϕ α=sgn(κy) ( r, (0, κ y )) . Doing so, we select the positive slope of the two branches of the dispersion relation Fig 3a. Bearing in mind the form of the graphene Bloch waves, Eq. (2), we can of course express the states ϕ b κy in terms of the sublattice wave functions ϕ F p := ( where up to a complex prefactor the coefficients f F pr are given by Note that the index r here denotes right (r = +) and left (r = −) moving waves [compare also to Fig. 3a].
Now we turn to the end states, emerging for purely imaginary κ y [30], which are allowed by both the Dirac equation Eq. (4) and the quantization condition Eq. (7b). In more detail, there exist two imaginary solutions for each κ x > 1/L y , which holds in ACNs for all κ x = nπ/L x , n ∈ N. Besides, the relation L x ≪ L y causes that to a very good approximation κ y = ±iκ x satisfies Eq. (4) and Eq. (7b). The corresponding dispersion relation is given in Fig. 3b. Notice that Eq. (6) is not applicable to κ y = ±iκ x (as explicitly exempted before), and instead the spinors fulfilling the Dirac equation Eq. (4) are given by η F αp (κ x , ±iκ x ) = δ p,±F . Using this in Eq. (2) and following the steps leading to Eq. (9), we obtain instead after straightforward insertions The corresponding ACN eigenstates can be chosen such that they live on one sublattice p = ± only: whereC(κ x ) is a normalization constant. It is evident that the decay length of ϕ e pκx ( r) from one of the ends to the interior of a specific ACN is L x /(n x π), which is always much shorter than the length of the ribbon. That is, one finds localized end states.
From the dispersion relation Eq. (5) it is easy to see that the energy of the end states is zero. Consequently, they will be unpopulated below half filling, but as soon as the Dirac point is reached, one electron will get trapped at each end (in an interacting system, Coulomb repulsion will hinder a second electron to enter). So we can conclude that at energies around the Dirac energy, not only the extended states with κ x = 0, but also the localized end states can be of importance.

Electron and Hamilton operator of the metallic ACN
All in all, the appropriate operator describing an electron with spin σ at position r reads in the low energy regimê whereĉ σk ,d σpk are the annihilation operators for electrons of momentum k and spin σ in the bulk or end states, respectively. The one-dimensional (1D) character of ACNs at low energies becomes evident by defining the slowly varying electron operatorŝ ψ rσ (y) := κy=(Z+0.5)π/Ly ϕ b κy ( r)ĉ σκy = 1 2L y κy=(Z+0.5)π/Ly e irκyĉ σκy (13) such that we obtain with ϕ F p ( r) := ϕ F p ( r, κ = (0, 0)).
From the dispersion relations Fig. 3, it is easy to give the Hamilton operator of the noninteracting metallic ACN, with the Fermi velocity v = 8.1 × 10 5 m/s corresponding to the absolute value of the slopes of the linear branches in Fig. 3a. There is no contribution of the end states as those have zero energy, see Fig. 3b. With the allowed values for κ y , Eq. (10), there results obviously a level spacing

The interaction Hamiltonian
This section is dedicated to the determination of the eigenstates of an interacting ACN. The cornerstones of the derivation and an analysis of the resulting spectrum were already presented in our recent short manuscript [29]. We provide here merely a compact outline of the technical steps and refer the interested readers for details to the appendices. A discussion of the main features needed to understand the transport results is provided in Sec. 4.
In the following we concentrate both on interaction effects regarding the extended bulk states in ACNs as well as on correlations between end and bulk states. Ignoring exchange effects, the many-body end states can be spin-or edge-degenerate. Above half filling, within a reasonable energy range, both end states can be assumed to be populated with one single electron only. This is because their charging energy exceeds, due to the strong localization in position space, the charging energy of the extended states by far: A simple estimation modelling the Coulomb repulsion between two electrons localized at the same ribbon end via the Ohno potential given below, Eq. (19), yields energies of the order of 0.1 eV for ribbon width around 10 nm. In contrast, typical charging energies for the bulk states of such ACNs range around 1 − 10 meV. The only relevant scattering processes between bulk and end states mediated by the Coulomb interaction are thus of the formV All other processes should be strongly suppressed for energy reasons. For the bulk-bulk interaction, our scattering potential is described by the expression In both Eqs. (17) and (18), the function U( r − r ′ ) models the screened 3D Coulomb interaction potential. For our calculations we use the Ohno potential [32], with U 0 = 15 eV [33] and ǫ ≃ 1.4 − 2.4 [34] the dielectric constant of graphene.
We proceed now with a discussion of the two different expressions. We start with the bulk-bulk processes, where the analysis follows largely the lines of an earlier work on interaction effects in metallic single wall carbon nanotubes (SWCNTs) [35].

Bulk-bulk interaction
With the help of the reformulation of the 3D electron operator in terms of the 1D operatorsψ rσ (y), Eq. (14), we obtain after integrating over the coordinates perpendicular to the ribbon axis an effectively 1D expression for the interaction, where {[I]} denotes the sum over all possible quadruples [I 1 , I 2 , I 3 , I 4 ], in the former case for the band index I = r. The spin-independent 1D interaction potential U [r] (y, y ′ ) reads with ⊥ indicating that the integration has to be performed over the coordinates perpendicular to y, y ′ (i.e. x, x ′ , z, z ′ ). Exploiting the explicit form, Eq. (11), of the coefficients f F pr we can easily identify those scattering processes which are indeed mediated by U [r] (y, y ′ ). In detail, where we have defined the potentials In total this means that we can rewrite Eq. (20) as § In the case of umklapp-and back-scattering with respect to r, the potential U [r] (y, y ′ ) is proportional to the difference of the inter-and intra-lattice interaction potential. Since the latter potentials differ only on the length scale of the lattice spacing a 0 ≈ 0.14 nm, this means that in the case of S r = u, b the effective 1D potential U [r] (y, y ′ ) can be considered as point-like. Introducing the coupling constant we can set in good approximation U [r] (y, y ′ ) = u δ(y − y ′ ) and write the short-ranged interaction processes aŝ Since u is derived from a short-ranged interaction it scales inversely with the size of the underlying ribbon. We find typical values [36] of uL x /ε 0 = 0.1 nm for a level spacing ε 0 , Eq. (16). The process V b−b uf + vanishes identically, because it involves the operator productψ −rσ (y)ψ −rσ (y) = 0. Only the forward scattering processes Since easily diagonalizable by bosonization, it is convenient to identify the densitydensity processes among the relevant bulk-bulk interaction processes, such with the density-density and non-densitydensity parts given bŷ

End-bulk interaction
In ACNs, we additionally have to consider scattering between the electrons living in the bulk of the ribbon and the electrons trapped in the end state existing at both zigzag terminations of the stripe. Below half-filling, the end states are unpopulated and thus all terms discussed in the following would be zero a priori. The range we want to concentrate on is the low energy regime above half-filling, where exactly one electron will permanently occupy each end state, so that we have in total two end electrons interacting with our bulk electrons.
We can certainly impose that the wave functions of the localized p z -orbitals have nonvanishing overlap for electrons on the same sublattice only. Moreover we demand that both end electron operators have to belong to the same end, and thus to the same sublattice, in order to give a nonzero contribution . If we insert then into Eq. (17) the decomposition Eq. (14) for the bulk electron operator and ψ ẽ pσ ( r) = κx ϕ ẽ pκx ( r)d σpκx for the end electron operator, and set in the coefficients from Eq. (11), we obtain: (24) Thereby, symmetry arguments have lead to the demand r = r ′ for the part containing the interaction potential related to densities, which, again for symmetry reasons, fulfils further We define a density-density part of the end-bulk interaction correspondingly aŝ where we exploited that each end state is populated with exactly one electron. The second potential, , can be considered point-like due to the localization of the end states at yp =− = 0 or yp =+ = L y , and hence simplifies to , with a short-range coupling constant which is independent ofp for symmetry reasons, For an ACN of width L x ranging from 5 to 25 nm, one finds [36] t κx nρ ≈ t nρ =: t, with tL x /ε 0 ≈ 0.55 nm, practically independent of κ x . This means that the short-ranged endbulk scattering is comparable in strength to the exchange interactions induced by the Operators acting on different ends would change the population of each end state, which is not allowed.
bipartite sublattice structure, and consequently we have to account for a non-negligible contributionV

Diagonalization of the ACN Hamiltonian
We can diagonalize the HamiltonianĤ ρρ by bosonization and subsequently express the total ACN Hamiltonian, in the eigenbasis ofĤ 0 +V ρρ . A numerical diagonalization of the so constructed total Hamiltonian, however, yields reliable results only away from half-filling: as the eigenbasis needs to be truncated for the calculation, it is crucial thatV nρρ is just a perturbation toĤ 0 +V ρρ in the sense that it only mixes states close in energy to each other. In the direct vicinity of the Dirac points, the processV b−b uf − breaks this demand, while it vanishes away from half-filling, as it will become clear in the course of this section.
Diagonalization of the density-density part Diagonalization ofĤ 0 +V ρρ can be achieved by bosonization. The end result of the procedure on which more details are given in Appendix A isĤ The first three purely fermionic contributions in Eq. (28) account for charging and shell filling. The last term counts the bosonic excitations of the system, created/annihilated by the operatorsâ † jq /â jq . The two channels j = c, s are associated to charge (c) and spin (s) excitations. The excitation energies are The energies of the charge channel are dominated by the long-ranged interactions via the coefficients Ly 0 dy Ly 0 dy ′ U intra (y, y ′ ) + U inter (y, y ′ ) cos(qy) cos(qy ′ ). (30) Finally, we can give the eigenbasis of where | N, σ e , 0 has no bosonic excitation. The fermionic configuration N = (N ↑ , N ↓ ) defines the number of electrons in each spin band. The occupation of the end states is determined by σ e = (σ e + , σ e − ), where '−' relates to y − = 0 and '+' to y + = L y as before. Below half filling the end states are empty, such that there is only one possible configuration: σ e + = 0 = σ e − . Above half filling, exactly one electron occupies each end state and thus σ e + , σ e − ∈ {↑, ↓}. Finally, m = ( m s , m c ), with m jq = ( m j ) q containing the information how many bosonic excitations are present in level q for mode j = c, s.
Non-density-density interaction In the following we use the states from Eq. (31) as basis to examine the effect ofV nρρ . For this purpose we evaluate the matrix elements of the potentialsV b−b nρρ andV e−b nρρ , using the bosonization identity for the 1D electron operators.
Generally,V b−b nρρ does not conserve the quantity m, while it must neither mix states with different electron configurations N , nor with different end spin configurations σ e : the Coulomb interaction between bulk electrons cannot change the quantity S z = 1 2 (N ↑ − N ↓ ), and it cannot touch the end states. Further, we already know nρρ are effectively local interactions, see Eqs. (23a) and (23b), such that the matrix elements of the non-diagonal bulk-bulk interaction scale with the exchange-coupling parameter u, The derivation of this expression is given in Appendix B, as well as the definitions of F (λ, m, m ′ ) andλ jq [r][σ] (y), Eq. (B.14) and Eq. (B.12), respectively. In fact it turns out that in comparison to the end-bulk non-diagonal interaction, the bulk-bulk nondiagonal interaction has only minor impact on spectrum and transport properties of narrow ACNs.
For the non-diagonal end-bulk interaction, we find, as explicitly evaluated in Appendix C, The action of this scattering process is to flip the localized spin at the p end, σ ′ p → σ p ! = −σ ′ p , while at the same time a bulk spin must be inverted to preserve the spin-projection The localized spin at the other end, i.e., thep ≡ −p end, must be conserved: σp = σ ′p . As a result, the degeneracy between the lowest lying states ofĤ 0 +V ρρ is removed; moreover, the spin-charge separation gets smeared [29].

Minimal model for the lowest lying states
The low energy spectrum resulting from the numerical diagonalization ofĤ ⊙ , Eq. (27), is discussed in Ref. [29]. Here to explain how the lowest-lying states in the truncated eigenbasis Eq. (31) transform under the influence ofV nρρ it is sufficient to restrict to a minimal set of states: For the even fillings, N c = 2n, n ∈ N, we have to take into account twelve states, allowing up to one fermionic excitation. The reason is that without an unpaired bulk spin no mixing can take place ¶. For the odd fillings, N c = 2n + 1, it is enough to include the eightfold degenerate ground state ofĤ 0 +V ρρ . To preserve lucidity, bosonically excited states are left out from our analysis, as they do not change qualitatively the mechanisms behind the observed effects. The N and m c conserving bulk-bulk scatteringV b−b nρρ can also be disregarded: besides for slight shifts in energy, forming linear combinations of states differing just in m s but identical in N, m c has not much impact. The restriction to the minimal model applies only for the subsequent analytics. Concerning the numerical calculations shown in Sec. 5, an energy cutoff of 1.9ǫ 0 above the ground states was used, including every energetically allowed bosonic or fermionic excitation.

Even electron fillings
For even electron fillings, N c = 2n, we want to restrict to the following subset of the states described by Eq. Notice that the second and the third case describes fermionically excited states. Building now all possible combinations for our minimal set of states for even fillings, we get the following set of possibilities: ¶ For transport, this is fatal, as introducing the end spins degree of freedom a priori yields four identical, completely decoupled channels. This makes the kinetic equations ill-defined. Only the endbulk interaction induced mixing guarantees a unique solution for the transport problem. There are four degenerate ground states |a , |b , |d σσ with energy E (0) 2n , while the remaining singly fermionic excited states |f σσ , |g σσ , |σ, σσ, σ have an energy E (f ) 2n . The end-bulk interaction can only mix states with same spin-projection S z . To the highest values, S z = ±2 , it belongs only one state and thus no mixing occurs. In contrast, the three states |d σσ , |f σσ , |g σσ with S z = sgn(σ) get coupled to each other. Also the four states |a , |b , |c ± with S z = 0 can in principle transform into one another under the influence of end-bulk scattering. With help of Eq. (33) we can set up the corresponding blocks of the full HamiltonianĤ ⊙ in the truncated eigenbasis (assuming n even + ): Diagonalization leads to the eigenstates and eigenenergies listed in Tab. 1. We employ there the abbreviation with α, α ′ ∈ {±1}. Obviously, ξ ++ (γ) > ξ +− (γ) , ξ −+ (γ) > ξ −− (γ), and as in our context γ ≃ t, and hence γ ≪ E 2n holds, we can rely on the relations   The resulting energy landscape is sketched on the left side of Fig. 4, where we used differently colored and shaped symbols to indicate the composition of states. In our simple model, we find then from Tab. 1 that the interaction has hardly lifted the degeneracies between the various states. It can be verified with Eq. (34a) that there is a slight splitting in the ground states, such that their energy grows from |g 1 to |g 3 . Also the excited states, of which only the two labelled ones turn out to be important for later explanations, are listed increasing in energy. Crucial is the mixing of states with different bulk and end configurations. We can single out linear combinations which are either symmetric or antisymmetric under exchange of the end spins. For the features we will observe in transport, the decisive entanglement is the one between the four S z = 0 states |a , |b , |c + , |c − , leading to two states containing the antisymmetric combination |a − |b : A ground state |g 1 with small contribution of the fermionically excited states |c + + |c − , and an excited state |e 1 where those dominate [as found with Eq. (34b)].
In the case of S z = ±3 /2 there is only one state each. For the three states with S z = ± /2, from Eqs. (28) and (33) the following mixing matrix is found (still n is assumed even): Table 2. Lowest lying eigenstates of an ACN filled with an odd number of electrons (σ ∈ {↑, ↓}). Due to spin-degeneracy, the total number of possible states is eight. The two states with |S z | = 3 /2 are marked by grey diamonds. Red and olive boxes stand for the symmetric components |a σ + |b σ and |c σ , respectively. Blue boxes label the antisymmetric combinations |a σ −|b σ . Notice that all states behave purely symmetric or antisymmetric with respect to end spin exchange.
The matrix is easily diagonalizable and yields eigenstates according to Tab. 2 at three distinct eigenenergies (compare also to Fig. 4, right). Notice that for the odd filling all emerging states are purely symmetric or antisymmetric under exchange of the two end spins. Thereby, the symmetric states |t 1 σ and |t 3 σ essentially have the same tunneling properties, because they only differ by the sign in front of |c σ . It is of crucial importance that the state |t 2 σ is formed by the antisymmetric combination |a σ − |b σ . Comparing the definition of the 2n states |a , |b and |a σ , |b σ , we see that from the 2n ground state |g 3 , a tunneling event can exclusively lead to one of the 2n + 1 states |t 1 σ or |t 3 σ . Via their |c σ components, these connect to |c + + |c − as well as to |d σσ , and thus to all the other labelled 2n states from Tab. 1, but the link to |g 1 , |e 2 σσ is weak due to Eq. (34b). This will be the key ingredient for the explanation of the stability diagrams in Sec. 5.

Spin-dependent transport across quantum-dot ACNs
Looking solely to the spectrum, there is no demand for distinguishing between symmetry or antisymmetry of a certain state under the exchange of end spins. In transport, however, this property turns out to lead to tremendous effects. The case of an unpolarized set-up was discussed in Ref. [29], where it was found that end-bulk entanglement leads to pronounced negative differential conductance (NDC) lines occurring for a completely symmetric setup. In this work the focus is on spin-dependent transport for collinear lead magnetizations. Strikingly, all the NDC features observed for the unpolarized set-up vanish for anti-parallel contact magnetization while they persist for the parallel case. As a consequence we predict negative tunneling magneto-resistance (TMR) within a narrow region along the edges of the Coulomb diamonds for even fillings.
Unless specified differently, we have employed the following parameters for all viewed plots: Energy cutoff E max 1.9ε 0 , Thermal energy k B T 0.01 meV, Charging energy * W 0 2.31 meV, Ribbon length L y 572 nm, Level spacing ǫ 0 2.93 meV, Ribbon width L x 7.8 nm, Bulk-bulk exchange u 0.036 meV, End-bulk exchange t 0.21 meV. Polarization strength P 0.8 The values of charging energy, bulk-bulk and end-bulk exchange-coupling were numerically verified for ribbon widths ranging from 5 − 20 nm.
Throughout this work we assume that the coupling between the electronic reservoirs, i.e. the contacts, and the ACN is weak. Under such condition, the total system is described bŷ with the ACN-HamiltonianĤ ⊙ given in Eq. (27). Further, the contacts are described byĤ leads = lq σ (ǫ q −µ l )ĉ † lσqĉ lσq , withĉ lσq annihilating an electron in lead l of kinetic energy ǫ q . The chemical potential µ l differs for the left and right contact by eV bias , with V bias the applied bias voltage. Next, H T = lσ d 3 r T l ( r)ψ † σ ( r)φ lσ ( r) + h.c. describes tunneling between ACN and contacts, where T l ( r) is the in general position dependent tunneling coupling andψ σ ( r) the ACN bulk electron operator as given in Eq. (14). The lead electron operator isφ lσ l ( r) = q ϕ lq ( r)ĉ lσ l q with ϕ lq ( r) denoting the wave function of the contacts. Finally, the potential term describes the influence of a capacitively applied gate voltage (0 ≤ α ≤ 1). Due to the condition that the coupling between ACN and the contacts is weak, we can calculate the stationary current by solving a master equation for the reduced density matrix to second order in the tunneling coupling. As this is a standard procedure, we refer to previous works [37, 38,39] for details about the method.
Besides the energy spectrum, the other system specific input required for transport calculations are the tunneling matrix elements of the ACN bulk electron operator, with the function F (λ, m, m ′ ) and the parameter λ jq rσ (y) as given in Appendix B in Eq. (B.14) respectively Eq. (B.11). We omit here the calculation of this identity, because the fermionic contribution follows straightforwardly from Eqs. (B.1)-(B.3), while the bosonic contribution emerges in the same manner as for the more complicated matrix elements involving more than one electron operator evaluated in Appendix B. Moreover, the detailed derivation of the corresponding expression for SWCNTs can be found in Ref. [40].

Collinearly spin-polarized transport without magnetic field
In Fig. 5, left, (a) and (b) we show the stability diagrams obtained for an ACN coupled to leads polarized in parallel and in anti-parallel, respectively. Due to electron-hole symmetry around the 2n filling, the transport characteristics are mirror-symmetric with respect to the central diamond. As anticipated at the beginning of this section, end states leave various signatures in the parallel case, Fig. 5, left, (a) which are absent in the anti-parallel configuration, Fig. 5, left, (b). In the following we give an explanation for all effects indicated in Fig. 5, left, (a) based on the minimal set of states discussed in Sec. 4.1, Tabs. 1, 2. Let us shortly recall their properties and alongside explain in which points we have to expect discrepancy with respect to the full set of real eigenstates which was employed for all calculations: Firstly, the eight states from Tab. 2, |t 1 σ , |t 2 σ , |t 3 σ and |σ, σ, σ , will in the following frequently be called the lowest lying 2n + 1 states. They occur at only three distinct energies ± √ 2t, 0.0. Inclusion of excitations within an energy cutoff of 1.9 ε 0 slightly lifts the degeneracy of |t 2 σ and |σ, σ, σ , and introduces eight high lying excited states which are almost degenerate. Secondly, for what concerns the even fillings, we refer to |g 1 , |g 2 σσ , |g 3 as 2n ground states. The fact that they are almost, but not perfectly degenerate, and that their energy grows from |g 1 to |g 3 , is not changed upon the inclusion of further excitations and plays some role in the following. Moreover, mixing between the states from Tab. 1 and bosonically excited 2n states takes place in general, but actually preserves the types of linear combinations occurring in Tab. 1, which is the relevant point for our explanations. In summary, the main effect of inclusion of excitations within an energy cutoff of 1.9 ε 0 is the lifting of the degeneracies among the excited 2n states. In fact, the lowest lying excited state will be of the same nature as |e 1 , and the separation to the state corresponding to |e 2 σσ exceeds √ 2t, thus it is well resolvable. With this additional information to Tabs. 1, 2 we can now start to explain the features marked by the different arrows in Fig. 5 (left). The dashed red arrow points towards a triple of three parallel lines which are split by √ 2t. Those mark transitions from the 2n ground states to the 2n + 1 lowest lying states. Hereby, the antisymmetric state |t 2 σ , associated to the second line of the triple, is special, because it is the only one strongly connected to the 2n state |g 1 . The first line of the triple is the |g 3 →|t 1 σ ground state transition line.
The blue dotted arrow marks the lines around the tip of the Coulomb diamond. Those appear in four clearly distinct positions, separated by about √ 2t . The lower lying triple of lines arises from transitions of the lowest lying 2n + 1 states to the 2n + 2 ground states. By coincidence of parameters, the highest line, which is split, follows also in a distance of about √ 2t and marks transitions from 2n ground states to the aforementioned higher lying excited 2n + 1 states arising upon inclusion of bosonic excitations. The solid green arrow highlights the negative differential conductance features (I), (II) and (III). The former two originate from trapping in the state |g 1 , while the latter occurs due to depletion of the transport channel |t 2 σ ↔|g 1 . The mechanisms work as follows: The NDCs (I) and (II) mark the opening of the 2n + 1 → 2n back-transition channels |t 1 σ →|e 2 σσ and |t 1 σ →|e 1 , respectively. The situation is sketched in Fig. 6. Once they get populated, from both of these excited 2n states the system can decay into any of the lowest lying 2n + 1 states, and in particular there is a chance to populate the antisymmetric state |t 2 σ . This state is strongly connected to the 2n ground state |g 1 , which contains a large contribution of the antisymmetric combination |a − |b . But in the region where the NDCs occurs, the forward channel |g 1 →|t 2 σ is not yet within the bias window such that |g 1 serves as a trapping state. Fig. 5a (right) confirms this explanation: the population of the state |g 1 is strongly enhanced in the concerned region where the back-transitions |t 1 σ →|e 2 σσ and |t 1 σ →|e 1 can take place, while the forward transition |g 1 →|t 2 σ is still forbidden. NDC (III) belongs to the back-transition |t 2 σ →|e 1 , which is a weak channel because |t 2 σ is a purely antisymmetric state, while the antisymmetric contribution in |e 1 is rather small. From time to time, nevertheless the transition will take place, and once it happens the system is unlikely to fall back to |t 2 σ , but will rather change to a symmetric 2n + 1 state. Thus the state |t 2 σ is depleted, and with it the transport channel |t 2 σ ↔|g 1 , which leads to NDC. The statement can also be verified from the plot of the occupation probability for |g 1 , Fig. 5 (right): a pronounced dark region of decreased population follows upon the NDC transition.
Major changes in the stability diagram of the ACN are observed for anti-parallel contact configuration, Fig. 5b (left): Compared to Fig. 5a (left), various transitions lines are suppressed and the NDC features have vanished. The reason is that an anti-parallel contact configuration as drawn in Fig. 1 favors (g 2 ) ↓↓ as 2n ground state, because intunneling of ↓ -electrons and subsequent out-tunneling of ↑ -electrons is preferred. As a consequence, all transport channels related to |g 1 and |g 3 are of minor relevance, which weakens various transition lines and in particular destroys the NDCs effects: in the anti-parallel configuration, the occupation of the trapping state |g 1 is significantly lowered, as seen in Fig. 5b (right). In detail, starting from the 2n ground state |g 3 = 1 √ 2 (|a + |b ), in-tunneling of a majority (↓ -) electron from the source takes the system to (t 1 ) ↓ = 1 √ 2 (|a ↓ + |b ↓ ) + |c ↓ . Via the |c ↓ (= | ↓, ↑, ↓ ) -component of this state, it is possible to tunnel out with a majority (↑ -) electron of the drain, yielding a transition to (g 2 ) ↓↓ . Similarly, also starting from |g 1 in-tunneling of a ↓ -and subsequent out-tunneling of an ↑ -electron changes the 2n ground state to (g 2 ) ↓↓ . Depending on the bias voltage, transport is either carried by ↑ -electron via the ground state channel (g 2 ) ↓↓ ↔ (t 1 ) ↓ , or by ↓ -electrons via (g 2 ) ↓↓ ↔ | ↓, ↓, ↓ , where | ↓, ↓, ↓ forms a blocking state unless a back-transition to the 2n excited state |↓, ↓↓, ↓ is energetically allowed.

Tunneling magneto-resistance
In Fig. 7 we have plotted the tunneling magneto-resistance (TMR), which is a measure for the ratio of the current in the parallel configuration, I P A to the current in the anti-parallel configuration, I P A . Along the edge of the 2n Coulomb diamond, the TMR acquires a negative value, i.e. I AP exceeds I P A . This is unusual: for lowest order calculations without Zeeman splitting between the spin species typically strictly positive TMR is observed [41,38,39]. For the ACN, however, the effect originates from a reduced feeding of the |g 1 trapping state. This statement can be confirmed by comparing its occupation probability for the parallel, Fig. 5a (right), and anti-parallel polarized case, Fig. 5b (right), in the concerned region. Namely, |g 1 lies slightly lower in energy than |g 3 -for the values we chose, the energy difference amounts, according to Tab. 1 and Eq. (34a), to about 6 k B T . Hence it can serve as a perfect trapping state within a narrow region along the edge of the 2n Coulomb diamond: here, the bias is high enough to allow the ground state transition |g 3 →|t 1 σ , but not |g 1 →|t 1 σ . Though the latter channel is weak in any case, it nevertheless provides a nonzero escape rate from |g 1 . That is why in Fig. 5a (right), in the region where the NDC mechanism Fig. 6 (II) can populate |g 1 , the occupation probability approaches 1 only straight along the edge of the Coulomb diamond, and a value of 0.6 − 0.8 further apart from it. In contrast, we observe no comparable increase of the |g 1 population in Fig. 5b (right), because for the anti-parallel configuration, as explained above, the transition channels involved in the NDC mechanisms are strongly disfavored compared to (g 2 ) ↓↓ ↔ (t 1 ) ↓ . For this reason, no trapping occurs and I AP can exceed I P A , leading to the negative TMR.

Magnetic field sweep
Finally we study, both for non-magnetic and collinearly polarized contacts (see Fig.  1), the transport behavior under the influence of an external magnetic field. In its presence, formerly degenerate states with different spin projections S z components become Zeeman split. This means, at a fixed gate voltage one half of the transitions occur at a higher, one half at a lower bias compared to the situation without magnetic field. In detail, simple thoughts can confirm that the forward transitions involving ↑electrons, i.e. 2n For three distinct values of the Zeeman splitting, Fig. 8 shows stability diagram zooming on the region between the 2n and 2n + 1 Coulomb blockade diamonds. Those plots complement Fig. 9a, where we show the differential conductance versus bias and Zeeman splitting, at a fixed gate voltage eαV gate ≈ 2.0 meV (marked in Fig. 8 with the dashed line). In turn, the three values of the Zeeman splitting considered in Fig. 8 are marked in Fig. 9a by dashed white lines. At first we focus on Fig. 9a. In Fig. 9b the case of polarized leads is considered. A small Zeeman splitting will actually select (g 2 ) ↑↑ as 2n ground state. The 2n → 2n+ 1 ground-state-to-ground-state transition is then (g 2 ) ↑↑ →(t 1 ) ↑ , as indicated left of the figure. It is the first line of a triple marking transitions to the lowest lying 2n + 1 states. Upon introducing a Zeeman energy, the spin-degeneracies of those are lifted, but only the two excited lines split in "V"-like manner, while the ground state-to-ground state transition (g 2 ) ↑↑ →(t 1 ) ↑ has only one rightwards slanted branch (i.e. raises in energy with increasing field). The reason is that (g 2 ) ↑↑ is connected to the |c ↑ (= | ↑, ↓, ↑ )component of the energetically favored state (t 1 ) ↑ = 1 √ 2 (|a ↑ +|b ↑ )+|c ↑ by in-tunneling of ↓ -electrons. There is no possibility for a transition with ↑ -electrons, hence a left branch does not exist. At a Zeeman splitting of √ 2t/2, the process (g 2 ) ↑↑ ↑ → | ↑, ↑, ↑ becomes the ground state-to-ground state transition. The crossover is marked with (P) in Fig. 9a. Due to the in-tunneling of ↑ -electrons, this resonance is continuously lowered in bias upon increasing the Zeeman energy further. At a Zeeman splitting of about 0.4 meV, we are exactly at resonance. As seen in the middle plot of Fig. 8, a line triple has clearly separated from the ground state transition line. Upon comparison with Fig. 9a it is immediately understood that it belongs to ↓electron transitions to the lowest lying 2n + 1 states. Concerning the corresponding ↑ -electron transitions, something interesting happens: In the point (P'), the left branch of the second "V"-shaped pattern, which belongs to the 2n + 1 state (t 3 ) ↑ , ends. The reason is that (t 3 ) ↑ consists of the same components as (t 1 ) ↑ . By in-tunneling of ↑electrons, it can thus not be connected to (g 2 ) ↑↑ , but rather to |g 3 [= 1 √ 2 (|a + |b )], see Fig. 9, sketch (P'). This 2n state is, compared to (g 2 ) ↑↑ , lifted by the Zeeman energy and can only be populated by back-transitions from (t 1 ) ↑ . The state (t 1 ) ↑ , however, is not available below the transition (g 2 ) ↑↑ → (t 1 ) ↑ [ Fig. 9, sketch (P'), dashed arrow]. This explains why the point (P') is positioned at the crossing with the resonance line marking this transition. Going on to a value of 0.8 meV of the Zeeman splitting, where the rightmost plot in Fig.  8 is taken, we reside at low bias voltages within the 2n + 1 Coulomb blockade diamond; the ground state transition is now the out-tunneling process | ↑, ↑, ↑ ↑ → (g 2 ) ↑↑ , thus raising in bias for an increasing magnetic field. The behavior reverts again in the point (P"), where the Zeeman splitting has lowered the excited 2n state |↑, ↑↑, ↑ enough to change the ground state transition to |↑, ↑, ↑ →|↑, ↑↑, ↑ , which involves out-tunneling of ↓ -electrons.
Finally, the Fig. 9b shows the data obtained if the calculation yielding Fig. 9a is performed for ferromagnetic leads, polarized in parallel to the applied field. The only thing changing is the intensity of the lines. In particular, several of them are transformed into negative differential conductance lines. Such an effect is expected for any type of single electron transistor with parallel polarized contacts in magnetic field: upon opening a channel to a state from which the system can only escape via a weak (in our case ↓ -mediated) transition, NDC occurs as such slow processes hinder the current flow. The only exception are the ground state-to-ground state transitions, i.e. the edges of the Coulomb diamonds, where current starts to flow: to those, of course always positive differential conductance (PDC) lines belong. An obvious example is the | ↑, ↑, ↑ −↓ → | ↑, ↑↑, ↑ transition, which turns from NDC to PDC beyond the point (P").

Conclusion
We have studied the transport characteristics of fully interacting graphene armchair nanoribbons (ACNs) attached to ferromagnetic contacts. Short-ranged Coulomb interactions play an essential role in such systems, leading to an entanglement between bulk states and the ones localized at the zig-zag ends of the ribbons, thereby lifting degeneracies between various states. Importantly, the entanglement breaks the otherwise strict conservation of the bulk spin-S z component, which leaves strong fingerprints in transport. The stability diagrams predicted for ACNs possess a two-electron periodicity and show already for a completely symmetric, unpolarized setup in zero magnetic field unique features like a characteristic transition line triple and pronounced negative differential conductance [29]. These effects, originating from the interaction-induced lifting and formation of states symmetric or anti-symmetric under exchange of the ribbon ends, are preserved for a parallel contact polarization. For an anti-parallel contact polarization, absence of various transition lines is observed due to spin-blockade effects and also the NDC features have vanished. The reason is that transition channels feeding the trapping state are disfavored, which leads even to a negative tunneling magneto-resistance. We have further investigated the transport behavior in magnetic field, for unpolarized as well as for in parallel polarized contacts. A change of the odd filling ground state with S z = /2 to one with S z = 3 /2 is observed at a Zeeman splitting of √ 2t/2, such that the value of the end-bulk exchange coupling t can directly be read off. Upon imposing a parallel contact magnetization, at several transition lines the differential conductance changes from the positive to the negative regime, because all ↓ -mediated transport channels become weak.
All in all, we found that short-ranged Coulomb interactions yield a strong influence of localized end states on the properties of ACNs. In particular, exchange makes the isolated bulk and end spin-S z components a bad quantum number: only the sum of both is a conserved quantity. Due to this fact, ACNs might not be as ideal candidates for certain spintronic devices as previously regarded. On the other hand, the entanglement is a rich source of ACN specific features in transport. Recent achievements in fabrication of carbon nanostripes with defined geometries [18,19] raise the hope of an experimental confirmation of our predictions within the near future. In fact, any term of the form Eq. (A.5) takes can be absorbed in the quadratic part of the Hamiltonian without any relevant impact on the spectrum, and we remain withĤ 0 +V b−b ρρ , which can be diagonalized in a standard way [42] by a Bogoliubov transformation [43]: One introduces new bosonic operatorsâ jq andâ † jq which relate to the old bosonic operators,b σq viâ The transformation coefficients B jq and D jq can be expressed in terms of X jq and A jq , which were introduced in Sec. 3.3: With our values, Eq. (29), for X jq and A jq we find approximately and for the transformation coefficients to the spin mode The transformation coefficients for the charge modes depend, as for SWCNTs [35,40], on the ratio g q = ε 0q /ε cq : Exploiting these relations yields then the diagonalized Hamiltonian Eq. (28).
The operatorη σ is the so called Klein factor, which annihilates an electron in the σ-branch and thereby takes care of the right sign as required from the fermionic anticommutation relations; in detail, K rσ (y) yields a phase factor depending on the number of electrons of spin σ, Finally, we have the boson fields iφ rσ (y), The fermionic part is then given by and the bosonic part reads where Q N [r]σ (y) = exp i π Ly N σ sgn(r 4 − r 1 ) − N −σ sgn(r 3 − r 2 ) + sgn(r 4 +r 3 −r 2 −r 1 ) 2 y . Hence, for S r = u, [r] u = [r, r, −r, −r], we obtain which is oscillating fast with N c = N ↑ + N ↓ and thus completely suppresses the S r = u contribution away from half-filling. The only remaining term in V b−b nρρ is consequently S r S σ = bf − , for which we get with [r] b = [r, −r, r, −r] We can now restrict our further analysis to the bosonic part M where we have introduced The function F (λ, m sq , m ′ sq ) = m s e −λ * â † sq e λâsq m ′ s is given by [40,35]  For q m sq − m ′ sq ≤ 1, the evaluation of this expression is problematic as the divergence arising from 1/[4 sin 2 (πy/L y )] remains uncompensated. Hence, the evaluation of the corresponding matrix elements needs special care. The origin of this divergence lies in the fact that, if no bosonic excitations are present, the N conserving processes depend on the total number of electrons in the single branches [compare to the fermionic contributions toĤ 0 +V ρρ in Eq. (28)]. Since the bosonization approach requires the assumption of an infinitely deep Fermi sea [42], this leads, without the correct regularization, necessarily to divergences. These findings are in complete analogy to the theory for SWCNTs [35]. In the following we exemplify the proper calculation where O(sin 2 ) contains only terms q (sin(qy)) tq with q t q ≥ 2 and which 'cure' the 1/ sin 2 (πy/L y ) divergence appearing in M