Chirality-Induced Spin Selectivity (CISS) Effect: Magnetocurrent–Voltage Characteristics with Coulomb Interactions I

One of the manifestations of chirality-induced spin selectivity (CISS) is the appearance of a magnetocurrent. Magnetocurrent is the observation that the charge currents at finite bias in a two terminal device for opposite magnetizations of one of the leads differ. Magnetocurrents can only occur in the presence of interactions of the electrons either with vibrational modes or among themselves through the Coulomb interaction. In experiments on chiral molecules assembled in monolayers, the magnetocurrent seems to be dominantly cubic (odd) in bias voltage while theory finds a dominantly even bias voltage dependence. Thus far, theoretical work has predicted a magnetocurrent which is even bias. Here we analyze the bias voltage dependence of the magnetocurrent numerically and analytically involving the spin–orbit and Coulomb interactions (through the Hartree–Fock and Hubbard One approximations). For both approximations it is found that for strong Coulomb interactions the magnetocurrent is dominantly odd in bias voltage, confirming the symmetry observed in experiment.


INTRODUCTION
Chirality-induced spin selectivity (CISS) is a term that classifies a collection of experimental observations on chiral molecules. These observations were made in photoemission, 1−4 Halltype, 5,6 and transport experiments 7−15 (for an extensive overview, see ref 16). Photoemission experiments show that a layer of chiral molecules has a different transmission probability for spin up and down electrons, i.e., the transmission probability is spin dependent. Hall-type experiments show that selfassembled monolayers of chiral molecules magnetize when placed on a substrate. This magnetization changes with the chirality of the molecules, 5 and it decreases over time. 6 In two terminal transport experiments, CISS manifests itself as the appearance of magnetocurrent (MC). MC is the observation that the currents for non-zero bias differ for opposite magnetizations of the lead. Theory initially has mostly focused on the spin dependence of the transmission. The spin−orbit coupling of the molecule's constituents in combination with the chirality of the molecule induce an asymmetry in the transmission probability for spin up and down electrons of the order 10 −5 % when no decoherence is considered. 17 It has been shown that the chirality of the molecule in combination with the spin−orbit coupling of the substrate can induce an asymmetry in the transmission probability for spin up and down electrons of a few percent, 18,19 consistent with theoretical work on 20,21 and findings in photoemission experiments. A spin dependent transmission does however not imply a MC. 19,22,23 MC simply is not possible in a fully coherent, noninteracting particle picture according to Buẗtiker's reciprocity theorem for two terminal systems 24 so that modeling beyond this simplified picture is necessary. 22 Thus, in addition to the chirality of the molecule and spin− orbit coupling, interactions need to be present 25 to obtain a magnetocurrent; in other words, interactions are a necessary ingredient for translating a spin dependent transmission into a non-zero MC. Several authors have made attempts to reveal MC in a theoretical description in chiral structures by including electron−phonon interactions, 26 electron−electron interactions, 27 or a generic decoherence probe. 19,23 Some authors 28,29 have proposed alternative explanations based on chiralityinduced interface magnetization. In describing two terminal transport measurements of the CISS effect there are three points that need to be addressed by theory. The first point is the Onsager−Casimir reciprocity, which states that the linear conductance terms at equilibrium are equal for opposite magnetizations: G 1 (m) = G 1 (−m). 22,30−32 Deviations from Onsager−Casimir reciprocity are not expected and theories should therefore reproduce this or provide a strong reason for their violation. 33 The second point is the odd/even behavior of the magnetocurrent ΔI in bias voltage. From Onsager−Casimir reciprocity, it follows that the magnetocurrent ΔI is nonlinear in bias voltage. 19, find that the MC is dominantly odd in bias voltage indicating that a cubic dependence dominates (ΔI ∝ V 3 ).
The theoretical work of Yang, van der Wal, and van Wees 23 modeled interactions with the vibrational modes via an extra node which is placed between the molecule and one of the leads, forcing the electrons to move through this node and thereby fully decohere. They found that the magnetocurrent is dominantly even in bias voltage (ΔI ∝ V 2 ). An analysis from our group that modeled interactions with the vibrational modes via the Buẗtiker voltage probe method and realistic parameters for the electronic structure came to the same conclusion. 19 Theoretical work of refs 25 and 34 for mesoscopic metallic samples (quantum Hall bar, chaotic cavity) finds that ΔI is dominantly even in bias voltage. In ref 35, ΔI is odd, but it seems to violate Onsager−Casimir reciprocity since ΔI is linear in bias voltage. The discrepancy between the odd ΔI−V characteristics of the experiment and the even ones of the theory remains a puzzling problem and it is this discrepancy which is the main topic of this work. The third point entails the size of the effect. In experiments, it can reach values of 1−80%, while from theoretical model calculations for realistic parameters, 19,36 values of less than 1% have been reported. Fransson 27 found that for Coulomb interactions in the Hubbard One approximation the size of the CISS effect can reach values on the order 10%. In this study, we will consider Coulomb interactions using the Hartree−Fock approximation (HFA) and the Hubbard One approximation (HIA) for more or less realistic parameters, focusing on the bias voltage dependence of the magnetocurrent. The article is structured as follows: in section 2, we describe the scattering region; in section 3, we present our numerical results; in section 4, we give an explanation of our results; and we present our main conclusions in section 5.

MODEL DESCRIPTION
The Hamiltonian of a molecular transport junction is given by where H os is the on-site Hamiltonian, H T is the hopping Hamiltonian, H SOC is the hopping Hamiltonian due to spin− orbit coupling, H U describes the Coulomb interactions, H lead−molecule describes the coupling of the molecule to the leads, and H leads describes the Hamiltonian of the leads. The onsite Hamiltonian is given by H os = ∑ k ϵ k n̂k, the on-site energy will be set to zero (ϵ k = 0) throughout this paper. The hopping Hamiltonian is given by where t is the hopping parameter and h.c. denotes the Hermitian conjugate. In this work, the sites are arranged in a helix with radius a and pitch c. N is the number of sites within one winding and M is the number of windings such that MN is the total number of sites in the molecule. For the hopping Hamiltonian due to spin−orbit coupling, we use the model of Fransson 27 which couples next- where λ is the spin−orbit coupling parameter, the components of σ⃗ are the Pauli matrices, v⃗ k = d ⃗ k+1 × d ⃗ k+2 and d ⃗ k+n = (r⃗ k − r⃗ k+n )/|r⃗ k − r⃗ k+n |, with r⃗ k the coordinates of site k on a helix: We take M = 1, N = 8 and a = 1, c = 1. Note that due to the spin-dependent hopping term the lattice is non-bipartite. H U contains the Coulomb interactions, we take those to be on-site: with U the Coulomb interaction strength. We model the leads using the wide-band limit, meaning that the self-energies are purely imaginary and independent of energy. The diagonal matrix elements of lead α that are coupled to the molecule are given by γ(1 + p z α σ z ) and are zero otherwise. Here γ is the coupling strength and p z α ∈ [−1, 1] is the magnetic polarization of lead α. For the right lead, we take p z R = 0, and we couple the left lead to the two leftmost sites and the right lead to the two right most sites. We aim at using realistic parameters corresponding to a molecule consisting of carbon atoms. We take the hopping parameter t = 2.4 eV. 37 Due to the image-charge effect, 38 the effective on-site Coulomb interaction of carbon, which for an isolated molecule is U C = 10.06 eV, 39 will be lowered to an extent which sensitively depends on the molecule−lead separation. To investigate the effect of U on the bias dependence of the magnetocurrent, we vary U to a maximum value of 4.8 eV. The spin−orbit coupling parameter of helicene is λ = 6 meV 17 (therefore, λ/t ≈ 10 −3 ). To also investigate its effect on the bias dependence of the magnetocurrent, we will vary λ between 10 −3 t and 10 −1 t. Furthermore, we take T = 300 K, the coupling strength the lead is taken as γ = 0.5 eV 40 and p z L = 0.5. The Coulomb interactions cause a shift in the on-site energies of the Hamiltonian of U/2 causing the molecular spectrum to be symmetric around U/2 for bipartite lattices. At zero bias voltage the chemical potentials of the left and right lead are given by the Fermi energy E F . In that case, for E U F 2 = , the molecule is charge neutral and E F lies precisely between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) energy. However, in molecular junctions, the molecule is rarely charge neutral due to charge transfer, which corresponds to E U F 2 (lying closer to either the HOMO or LUMO energy). Therefore, we also vary the Fermi energy between the energy of the HOMO and LUMO level. In appendix A, we state the retarded and advanced two-point Green's functions in the Hartree−Fock and the Hubbard One approximation derived with the equation of motion technique (analogous to chapter 12 of ref 41). These Green's functions are expressed in terms of the average electron densities for site k with spin s, ⟨n ks ⟩, which we determine self-consistently from the Green's functions; see eq (eq 11). Every iteration m has an input and an output electron density, and as convergence criterion for the mth iteration, we use: |⟨n ks in,m ⟩ − ⟨n ks out,m ⟩| < 10 −5 . The Hamiltonian without interactions (U = 0) is defined as H 0 = H os + H T + H SOC and is constructed with the Kwant code 42 and the Qsymm code. 43 We have implemented a nonequilibrium transport code which can be found in https://github.com/ khhuisman/CISS_CoulombInteraction. In this code, we determine the electron density as follows. Suppose we want to calculate the electron density for the decreasing or increasing bias voltages {0, V 1 , V 2 , ...}, (|V i+1 | > |V i |). First of all, we start our self-consistent calculation at zero bias voltage where we expect that every site is approximately half filled; therefore we take this as an initial guess ( )

RESULTS
We now turn to a two terminal system with Coulomb interactions. The transmission then depends on the bias voltage V through the Coulomb potential, and in the HFA and HIA, this is expressed in terms of the electron densities, that depend on the bias voltage applied to the molecule. This means that the transmission becomes voltage dependent: T LR (m) → T LR (m, V) = T LR (m, ⟨n 1↑ (m, V)⟩, ⟨n 1↓ (m, V)⟩, ..., ⟨n ks (m, V)⟩). The current into the left lead is then given by where f α = f(E, μ α , β) is the Fermi−Dirac distribution of the lead α with chemical potential μ α at Using eq 3, we can write the magnetocurrent as In Figure 1, ΔI(m, V) is plotted as a function of bias voltage in the HFA and HIA respectively. In Figure 1a, the magnetocurrent is plotted for E U F 2 in the HFA, and we see that for small U/t (upper panels) the magnetocurrent is dominantly even in bias voltage and for large U/t (lower panels) the magnetocurrent is dominantly odd in bias voltage for E U F 2 . In Figure 1b, the

DISCUSSION
In this section, we provide theoretical arguments to qualitatively explain our numerical results. We can expand the magnetocurrent (eq 4) in bias voltage as ΔI(m, Onsager−Casimir reciprocity, 30 where ∂ ks = ∂ ⟨nd ks (m,V=0)⟩ . Equation 7 can be understood as follows: the Green's function is time-reversal symmetric except for the electron densities. As the zeroth order term Δn ks (m, V = 0) = 0, ΔG 2 (m) can only scale with the first order derivative of Δn ks (m, V) at V = 0. Intuitively this makes sense: to what extent TRS is broken out of equilibrium (V ≠ 0) scales with the occupation difference: Δn ks (m, V) = ⟨n ks (m, V)⟩ − ⟨n ks ̅ (−m, V)⟩. Deviations from TRS thus manifest themselves through Δn ks (m, V). If this deviates from zero, ΔG 2 (m) will too. If TRS is present out of equilibrium for every voltage (⟨n ks (m, V)⟩ = ⟨n ks ̅ (−m, V)⟩) then ΔG 2 (m) ∝ ∂ V Δn ks (m, V = 0) = 0 and there is no magnetocurrent as expected due to Buẗtiker's reciprocity theorem for two terminal systems. 24 Furthermore, if U = 0, then ΔI(m, V) = 0, since every derivative with respect to bias V yields an electron density multiplied by U. This shows the importance of going beyond the noninteracting particle picture. From eq 5 (SI section 1.1.3), we have for ΔG 3 (m): . This is consistent with Figure 1a since ΔI(m, V) changes from even to odd with increasing U, and we see that for large U, the cubic term in the magnetocurrent tends to dominate. In the HIA we see similar behavior only for much smaller U than considered in Figure 1b and V = 0, the system is not exactly half-filled which is more the rule than the exception in molecular junctions due to charge transfer to the molecule. Due to the large Coulomb interactions this will result in a magnetocurrent which is odd in bias voltage.

CONCLUSION
In this work, we studied the voltage dependence of the magnetocurrent for a system with Coulomb interactions (in the HFA and HIA). The system we studied has next-nearest neighbor, spin-dependent hopping that causes the lattice to be non-bipartite. Our numerical results show that the magnetocurrent is odd in bias voltage in both the HFA and HIA for strong Coulomb interactions (U > t) in agreement with experiments. 8−15 Furthermore, we verified that the Onsager− Casimir reciprocity is satisfied, as expected. For a large spin− orbit coupling parameter (λ/t = 0.1), we found that the size of the effect is on the order of 0.1% which is of the same order as our previous work on Buẗtiker voltage probes. 19 How a bipartite lattice with spin−orbit coupling affects the voltage dependence of the magnetocurrent will be considered in a separate work.

■ APPENDIX A: ELECTRON GREEN'S FUNCTION
A derivation of the Green's functions is done in Chapter 12 of Haug and Jauho. 41 In this section, we simply give the Green's functions for a system with spin−orbit coupling.  (10) where Σ is the retarded self-energy which in the wide band limit for a magnetized left lead is given by and n = ∑ ks ⟨n ks ̅ ⟩nk s is a diagonal matrix with the electron densities on the diagonal. For these Green's functions, the density of states are symmetric around the Fermi energy E U F 2 = if the lattice is bipartite. The electron density for site k with spin s is given by In eq 9, the term ∑ ks ⟨n ks ̅ ⟩nk s is the result of a Wick contraction of the term n k↑ n k↓ in eq 2. This contraction in principle also allows for the term 45,46 −⟨c k↑ † c k↓ ⟩c k↓ † c k↑ + h.c., which we will call the noncollinear Hubbard model. In that case, we obtain (SI, section 2) the following noncollinear (NC) Hartree−Fock Green's function: where we defined ρ = ∑ ks ⟨c ks ̅ † c ks ⟩c ks † c ks ̅ . The expectation value ⟨c ks ̅ † c ks ⟩ vanishes in the absence of spin−orbit coupling, and it is expected that it does not change the result of our calculations. Indeed, when we include this term in our calculations the bias dependence of the magnetocurrent does not change with respect to ⟨c ks ̅ † c ks ⟩ = 0 for all of the considered values of the spin−orbit coupling, Coulomb interaction strength, and Fermi energy. Only in the HFA for 2 U t = , 0.1 t = , and E U F 2 the size of the effect increases significantly to a few percent for low bias voltage. However, this increase in the size of the effect happens for a particular choice Fermi energy, and it is not clear to us why that happens; however, we suspect this is due to a numeric instability.