First-principles calculations of steady-state voltage-controlled magnetism: Application to x-ray absorption spectroscopy experiments

Recent x-ray absorption experiments have demonstrated the possibility to accurately monitor the magnetism of metallic heterostructures controlled via a time-independent perturbation caused, for example, by a static electric ﬁeld. Using a ﬁrst-principles, nonequilibrium Green’s function scheme, we show how the measured dichroic signal for the corresponding steady-state situation can be related to the underlying electronic structure and its response to the external stimulus. The suggested approach works from the inﬁnitesimal limit of linear response to the regime of strong electric ﬁeld effects, which is realized in present experimental high-sensitivity investigations.

The creation, control and detection of spin polarized currents lie at the heart of spintronics.Accordingly, a number of magneto-optical experiments have been proposed over the years to explore related magnetic phenomena and corresponding materials.Obviously, Xray-based measurement techniques seem to be especially suited for this purpose because they provide a high signalto-noise ratio, precise targeting of individual chemical species and tunable sensitivity to either surface or bulklike properties of a sample, depending on photon energy [1,2].
In particular the time-dependence of magnetic properties has been studied over a broad range of time scales by means of magneto-optical techniques, addressing for example ultra-fast demagnetization [3][4][5][6][7][8][9] or optical manipulation of the magnetic order [10][11][12][13].In general, the most comprehensive information can be acquired by corresponding experiments using pump-probe techniques.Early experiments in the field, such as time-dependent MOKE (magneto-optical Kerr effect) or XMCD (X-ray magnetic circular dichroism) experiments [3,4] have been performed by controlling the investigated system via an external magnetic field and using different delay times for the read-out via x-ray pulses.
Another type of interesting experiments exploits the XMCD to study the impact of an external static electric field on the magnetization.For the resulting outof-equilibrium situation one may have pure charge rearrangements with a continuous flow of electric charge prevented by applying the electric field across an interface to vacuum or an insulating layer, i.e. having a capacitor like experimental setup [14,15].For the combination of conducting subsystems, on the other hand, a steadystate out-of-equilibrium situation will be created with a constant electric current flowing [16,17].In this case one may focus on the electric field induced change of the magnetization longitudinal [16] or transverse [17] with respect to the electric field.The electric field-induced electric current will in general be accompanied by a spin current that might be used for example for switching the magnetization via the spin transfer torque (STT), spinorbit torque (SOT) or the spin Hall effect (SHE) [18].Accordingly, the question concerning the connection between the observed XMCD and the induced spin current arises in a natural way.
In the following, a theoretical description of the electric field induced XMCD for the case of a conducting system is given, with a focus on the longitudinal setup investigated by Kukreja et al. [16].These authors investigated the bilayer system Co/Cu by XMCD measurements at the Cu L 3 -edge.It was found that a voltage applied across the layer system changes the XMCD spectra primarily in the vicinity of the Fermi level.By switching the sign of the voltage it was possible to separate in a reliable way the electric field induced contribution from that due to the so-called proximity effect (see below).This allowed in particular to demonstrate the linear dependence of the induced changes for the XMCD spectra on the applied voltage.Here, we present an ab-initio description of the observed phenomenon, that accounts in particular for the out-of-equilibrium situation when a finite voltage is applied to a conducting system.To deal with the XMCD in this case, an appropriate expression for the X-ray absorption coefficient together with a corresponding extension of the XMCD sum rules [19][20][21][22] are suggested.This allows to relate the XMCD spectra of an atom to its spin and orbital moments under the new out-of-equilibrium conditions.
Considering the electronic structure of a system out-ofequilibrium, the occupied and unoccupied states can be represented in an appropriate way in terms of the lesser and greater Green functions, respectively [23].Here we deal with the steady state situation encountered for a layered system exposed to a constant electric field across the layers, or equivalently a layered system connected to left and right leads with a corresponding voltage drop Φ in between.In this case, it is most convenient to consider the lesser and greater Green functions, G < Φ (E) and G > Φ (E), respectively, as functions of the energy E and dependent on the applied voltage Φ.For the setup of a layered system as sketched schematically in Fig. 1 it was shown by several authors [24][25][26] in terms of the retarded and advanced Green functions, G r Φ (E) and G a Φ (E), respectively, that in turn can be calculated using standard techniques [30].In this context, the finite voltage Φ determines the finite and fixed difference in the chemical potentials of the left and right leads, µ L = E F + Φ/2 and µ R = E F − Φ/2, respectively, where E F is the Fermi level of the unperturbed reference system with Φ = 0 V.To deal with the gradual voltage drop across the region between the leads, denoted interaction zone in Fig. 1, a corresponding model suggested in the literature [31,32] was adopted.(See Appendix A for more details).
The scheme to calculate the lesser and greater Green functions, briefly sketched above, was implemented within the framework of local spin density approximation making use of the Korringa Kohn Rostoker (KKR) Green function technique as suggested by Ogura and Akai [33].As the study aimed among others at the calculation of spin-orbit induced XMCD spectra, a fully relativistic formulation of the scheme was used on the basis of the four-component Dirac equation, i.e. all Green functions are 4 × 4-matrix functions [34].
With the lesser and greater Green functions G ≶ Φ (E) available, many electronic properties of interest can be calculated straightforwardly.For example the density of states (DOS) related to the occupied (<) and unoccupied (>) states is given by n where the angular-momentum representation used for the Green function [34] allows in a simple way to project out the contribution due to d-electrons as indicated by the optional index d.In a similar manner one has for the spin magnetic moment m where β is one of the standard Dirac matrices, σ z a 4 × 4-Pauli matrix [35] and the integration regime has been restricted according to the maximum of µ L and µ R (See Appendix B for more details.) For the steady-state situation considered here, an appropriate expression for the X-ray absorption coefficient µ q,λ Φ (ω) can straightforwardly be derived by replacing in the standard expression for an unperturbed ground state [36] the product of the retarded Green function and the Fermi distribution at given temperature, i.e.G r (E) [1 − f T (E)], by the greater Green function G > Φ (E) representing the unoccupied final states.This leads to: where Ψ i represents an occupied initial core state i while X q,λ stands for the interaction of the electrons with the photons with wave-vector q, energy ω and polarization λ [36].Dealing with X q,λ we make use of the dipole approximation.Accordingly, the index q can be omitted in the following.Finally, for the case Φ = 0 V the rather general expression for µ λ Φ (ω) in Eq. ( 3) is of course fully equivalent to the standard one in terms of G r (E) [36].This is demonstrated by a numerical example in Appendix C.
An extremely attractive tool to interpret experimental as well as theoretical XMCD data, i.e. the difference in absorption of left or right circularly polarized x-rays, is provided by the XMCD sum rules [19][20][21][22].Focusing on L 2,3 -spectra, the spin-related sum rule is given by where m d is the spin magnetic moment associated with d-states, T d z denotes the magnetic dipole or asphericity  4) was initially derived for the situation that no external voltage is applied to the system.Nevertheless, it can be straightforwardly adapted to the situation of a finite applied voltage Φ by calculating the absorption coefficients µ λ Φ (ω) by means of Eq. ( 3) and using the greater Green function G > Φ (E).The quantities m d , T d z and N d h , on the other hand, have to be calculated accordingly using the lesser Green function G < Φ (E).Finding the number of d holes N d h requires in particular an integration of the DOS n <d Φ (E) related to G <d Φ (E) up to the energy where it goes to zero.This upper threshold energy lies typically above the original Fermi energy E F for the case Φ = 0 V. On the other hand, the lower threshold energy from which on unoccupied states may contribute to the X-ray absorption via To avoid numerical problems and to ensure pronounced field-induced changes of the magnetic properties the bilayer system Co/Pd has been considered in this work, as Pd has a much higher spin polarizability than Cu.Otherwise, the set-up sketched in Fig. 1 followed essentially the work of Kukreja et al. [16] on Co/Cu; i.e. we investigated a fcc (001) textured Co/Pd bilayer system consisting of 8 layers of Co with their magnetization oriented perpendicular to the layers (out-of plane) and 8 layers of Pd.On both sides 4 buffer layers of Cu have been added to allow for a smooth connection to the fixed Cu leads with their respective chemical potential shifted by the voltage drop Φ.
The top panel of Fig. 2 shows the resulting profile of the spin magnetic moment m n Φ for the Co/Pd bilayer for Φ = 0 and ±0.34 V, where n is the layer index.For Φ = 0 V one finds an enhancement for the Co moments at the Co/Pd interface [37], while for the Pd layers there is an appreciable induced moment that decays rapidly with the distance from the interface (proximity effect).These well known features of the magnetization profile at the Co/Pd interface are obviously essentially unchanged when a finite voltage Φ is applied, that modifies the magnetic moments depending on the sign of Φ.The corresponding electric field induced contribution to the magnetization profile ∆m n Φ = m n Φ − m n 0 is shown for Φ = +0.34V in the bottom panel of Fig. 2. As one can see, the changes of the Co and Pd moments are of the same order of magnitude but interestingly of different sign.In addition one notes that ∆m n Φ has nearly exclusively to be ascribed to the part of the d-electrons (∆m d n Φ ).The layer resolved induced magnetic moments, and with these also the averaged induced moments of Co and Pd, scale fairly well linearly with the applied voltage Φ (see below).
These observed properties of the induced magnetic moments suggest an alternative description of the phenomenon of an electric field induced magnetization by means of Kubo's linear response theory.The corresponding Kubo-Bastin like expression for the non-local response coefficient p nn ′ zz that gives for layer n the induced magnetization along the z-axis due to an electric field in layer n ′ along the same direction can be found in Appendix D. Corresponding numerical results for the local layer dependent polarization coefficient p n zz defined by the sum n ′ p nn ′ zz are given in Fig. 3 for the investigated Co/Pd bilayer.Obviously, the profile of p n zz is in full accordance with the profile of the induced magnetic moment ∆m n Φ shown in Fig. 2. In particular, the different sign for the induced moment for Co and Pd is fully confirmed by the linear response calculations.In addition, we find also that p n zz is by far dominated by its part stemming from the d-electrons (p d n zz ).Here it should be stressed that a perfect one-to-one correspondence of ∆m in Fig. 2 and p (d) n zz in Fig. 3, respectively, cannot be expected as the evaluation of the induced magnetic moment requires the non-local response coefficient p nn ′ zz together with the local electric field E n ′ z for the various layers.However, the later one will depend in a non-trivial way on the layer index n ′ as it is suggested by the local electrical conductivity σ n zz that varies in an appreciable way from layer to layer for the Pd subsystem as can be seen in Fig. 3 (orange circles).Nevertheless, the com-  and 3, respectively, clearly shows that the electric field induced magnetization is first of all a manifestation of the Edelstein effect [38,39].Concerning this, it should be stressed that the lack of inversion symmetry as the central precondition for the occurrence of the Edelstein effect [38].This is of course given for the Co/Pd interface region.
The proximity effect mentioned above, i.e. the occurrence of an induced spin magnetic magnetization in an non-magnetic metal at the interface to a ferromagnetic metal, could be demonstrated in the past for several systems by means of XMCD measurements on the nonmagnetic metal.A prominent examples for this is Pt in Co/Pt [40] but also Cu in Co/Cu [41] with the high and small, respectively, spin susceptibility of the nonmagnetic metal reflected by the correspondingly high and small induced spin magnetic moments that differ by one to two orders of magnitude (see also Fig. 2).
In an earlier study we predicted that the magnetization induced by an external magnetic field should give rise to a corresponding XMCD signal that scales in a one-to-one manner with the spin and orbital susceptibilities [42].This could indeed be demonstrated experimentally for the non-magnetic metals Pt and Pd [43] but also for Au [44].In a completely analogous way, one can expect that the magnetization induced by an external electric field via the Edelstein effect leads also to the occurrence of an XMCD in case of an otherwise non-magnetic system.For an element in a system with a spontaneous or induced magnetization, on the other hand, this mechanism will alter the already present XMCD spectrum accordingly.This was indeed found when calculating the Pd L 2,3 -spectra for the investigated Co/Pd bilayer system when a finite voltage drop Φ is applied (see Appendix E).To deduce the relationship of the additional XMCD signal and spin magnetic moments induced by the electric field, the L 2,3 XMCD spectra for the individual Pd layers have been calculated in a first step for the voltages ±Φ.Taking for each layer the difference leads accordingly to the field induced contribution of the XMCD spectra, that have been analyzed by making use of the modified version of the spin sum rule given in Eq. ( 4).This leads finally to the field induced spin magnetic moments ∆m XMCD d n Φ of the d-electrons as deduced from the XMCD spectra.As a corresponding experiment does not allow to distinguish the XMCD spectra of the individual Pd layers n, the average over all these have been taken.This average value ∆m XMCD d Φ is given in Fig. 4 4)) (full circles).
approximations used to derive the XMCD sum rules [36] the agreement between the directly and XMCD-derived moments is rather satisfying.This finding is very similar to the results of studies on systems with the magnetization occurring spontaneously or induced by the proximity effect [45].
The close connection between the spin magnetic moment ∆m d Φ calculated directly for the Pd layers and ∆m XMCD d Φ deduced from the spectroscopic data clearly shows that the electric field induced XMCD indeed reflects in a one-to-one manner the induced spin magnetic moments for a steady-state situation.Accordingly, we conclude that the experimental findings of Kukreja et al. [16] for Co/Cu have to be interpreted as well this way; i.e. that the observed additional XMCD signal due to the applied voltage has to be seen as a rather direct measure for the electric field induced spin magnetic moment and as a manifestation of the Edelstein effect.The fact that the calculated as well as the observed field induced changes of the XMCD spectra occur primarily in the vicinity of the Fermi level supports this conclusion and interpretation of the experimental findings.Nevertheless, it should be mentioned that the observed XMCD spectra may have additional contributions due to spin accumulation caused by the specific features of the experimental set-up not accounted for within the present theoretical study.In summary, a numerical study on the electric field induced changes in the electronic, magnetic and spectroscopic properties of the ferromagnet/non-magnet bilayer system Co/Pd has been presented.The field induced magnetic moments were found to scale essentially linearly with the applied voltage.This finding as well as additional linear response calculations allow to ascribe the induced magnetic moments to the Edelstein effect.To make contact with recent XMCD investigations on Co/Cu an appropriate expression for the calculation of the X-ray absorption coefficient for a steady-state outof-equilibrium situation was suggested.Together with a corresponding extension of the XMCD sum rules the one-to-one corresponding relation between the electric field induced contributions to the spin magnetic moment and XMCD spectra could be demonstrated.From this it could be concluded that the results of a recent experimental XMCD study on the bilayer system Co/Cu first of all reflect the field induced magnetic moment in the vicinity of the Co/Cu interface.
We gratefully acknowledge insightful discussions with Prof. Hisazumi Akai, Institute for Solid State Physics (ISSP), Tokyo.This work was supported by the Deutsche Forschungsgemeinschaft grant: DFG EB 154/35.

Appendix A: Accounting for finite voltage effects
In order to include effects from a non-zero applied voltage, the initially equal chemical potential of the left and right leads of the considered trilayer system is rigidly shifted by ± Φ 2 with respect to the Fermi level E F coming from selfconsistent calculations performed without any external perturbation, i.e., in terms of the retarded Green function (GF) alone.This translates into a constant shift for the potentials which describe atoms of the left or right semi-infinite leads.
Different levels of self-consistency can be deployed to study the effects of such boundary conditions on the quantities of interest within the interaction zone, i.e., the middle part of the considered trilayer system.In this work we rely on a "ohmic conductor" assumption, that makes use of Ohm's law V = R × I to simply extend the same manipulation of the potentials also to the interaction zone.The layer-dependent voltage offset is applied gradually in a piecewise constant way, so that the whole interval − Φ 2 , + Φ 2 is covered accross the interaction zone.The difference between the offsets for adjacent layers is, in principle, proportional to the local resistivity, which undergoes only relatively minor variations across the interaction zone.

Appendix B: Practical calculation of observables
When evaluating the induced spin moment via Eq.( 2) in the main text, the lesser GF goes by definition to zero for unoccupied states.The choice of the upper integration limit is hence not crucial, provided that it samples a large enough energy range up to which the applied voltage can shift spectral weight, with respect to the original Fermi level for Φ = 0 V.This interval has a width roughly proportional to the magnitude of the external perturbation: − Φ 2 , + Φ 2 , with a cutoff that becomes smoother with larger electronic temperature (if applied) in the Fermi-Dirac distribution that enforces the different definition of lesser and greater GFs in the steadystate non-equilibrium Green function (NEGF) scheme of Refs.[26,33].
Obtaining the observables within the NEGF formalism requires that the energy integrations [as in Eq. ( 3) of the main paper] are performed on the real axis.One cannot make use of the complex contour integration as it is common when dealing with G r or G a .This makes both the Brillouin zone and energy integrations numerically quite demanding.One has to address numerical problems arising from the need to evaluate the relatively small differences between the majority and minority integrated density of states (IDOS).
This difficulty can be addressed by exploiting the fact that changes in the electronic structure due to the NEGF self-energy are typically significant only close to the Fermi level.Following proposals in the literature [26], we split the energy integration for the spin projected IDOS over two sub-intervals, one of them covering most of the energy range up to a small value below E F and the other one covering the rest.Within the first interval (lying lower in energy) the lesser GF can be safely replaced by the retarded GF, computed for the same local potentials, as discussed in Appendix A above.This allows to resort to well-known Gaussian-Legendre quadrature in the upper complex semiplane and invest more computational resources for the adjacent energy window at higher energy, where we sample the lesser GF along a contour parallel and close to the real energy axis.Although not adopted for the present investigations, analytical continuation schemes for each individual term of the underlying KKR-NEGF expression also appear viable, after adaptation from the retarded GF case [46].The expression for the absorption coefficient µ λ Φ (ω) given in terms of the greater GF in Eq. (3) in the main text is applicable also for the equilibrium case (Φ = 0).This offers the opportunity to check the numerical results based on Eq. ( 3) with results based on the standard expression using the retarded GF [36].Fig. 5 shows a corresponding comparison for the L 2,3 -edge of pure Cu, that confirms the correctness of our generalized expression for calculating the absorption coefficient.

Appendix D: Edelstein effect calculations
In the linear response regime, the static spin polarization induced by an external electric field applied to a noncentrosymmetric solid is described by the so-called Edelstein or inverse spingalvanic effect [38,47].Here spinorbit coupling (SOC) and absence of inversion symmetry lead to an induced spin polarization of the conduction electrons.
Using Kubo's linear response formalism the Edelstein response tensor p that gives the relation between spin polarization s and electric field E, s = pE, can be given by an expression analogous to the Kubo-Bastin formula for the conductivity tensor [48], where f (E) is the Fermi-Dirac distribution, evaluated here for the same chemical potential µ = E F at T = 0 K throughout the trilayer system, G r (G a ) is the retarded (advanced) single-particle Green function at energy E (arguments have been dropped for the sake of readability), Pµ is the spin-polarization operator and ĵν is the current density operator (see below).The presence of f (E) implies that in the temperature limit T → 0 K the first term in Eq. (D1) above has to be evaluated only for the Fermi energy E F (Fermi surface term), while the second one requires an integration over the occupied part of the valence band (Fermi sea term), which has been omitted based on the experience with the anomalous and spin Hall coefficient [49].The product of two Green functions in Eq. (D1) implies that the response in an atomic cell n is due to the perturbation in all cells n ′ .This leads in a natural way to the site resolved response function p nn ′ µν mentioned in the main text.The current (density) operator ĵν = −|e|cα ν in Eq. (D1) represents the perturbation due to the electric field component E ν .Adopting a fully relativistic formulation to account coherently for the impact of SOC, ĵν is expressed by the corresponding velocity operator vν = cα ν , where c is the speed of light and α ν is one of the standard 4 × 4 Dirac matrices [50].The spinpolarization operator Pµ on the other hand represents the spin-magnetic moment along the µ axis induced by the electric field, and is therefore most conveniently expressed in a relativistic four-component Dirac notation by: with σ µ being one of the standard 2 × 2 Pauli matrices [50].Figure 3 in the main text shows results obtained by adding an imaginary part of 10 −5 Ry to the Fermi energy and assuming T = 0 K.
Appendix E: Evaluating the sum rule use of the modified sum rules given in Eq. ( 4) in the main text, similar considerations concerning the energy path and integration as outlined above also apply to the greater GF, which is used to compute the XAS/XMCD spectra in the non-equilibrium situation.Far above the Fermi level, results remain unmodified with respect to a retarded GF calculation, obtained under direct inclusion of the layer-to-layer voltage drop within each atomic potential.This allows to invest more computational effort for the low-lying contributions from the greater GF.
This work deals with layered systems via the tightbinding KKR (TB-KKR) formalism, meaning that for spectroscopy application, further care is required in choosing a suitably large (TB-KKR) screening parameter [51].Namely, the greater GF needed to calculate the spectra has to be evaluated for higher energies than the lesser GF needed to calculate the ground state observables [compare Eqs. ( 2) and (3) in the main text].The use of the TB-KKR scheme and considering energies well above the Fermi energy also implies a larger angular momentum cutoff, as already known from studies based on the retarded GF alone.Throughout this work, XAS/XMCD calculations were performed by using a TB-KKR screening parameter of 12 Ry, and spherical harmonics cut-off at ℓ max = 4.
For the practical use of the sum rule expression, one finally has to determine the energy cutoff parameter E cut up to which 10 d electrons are entirely accounted for.Similarly as in Ref. [52,53], we find this energy by requiring that the density of the d states accounts exactly for 10 electrons when integrated up to E cut .For transition metal systems as considered in Ref. [16], we are typically in a situation where this IDOS value from the sum of lesser and greater GF lies fairly above the E F − Φ 2 , E F + Φ 2 energy window, in which the NEGF self-energy is operative.Noting the complementarity in the definition of lesser and greater Green functions [26,33], we can exploit the DOS relationship: n < (E) + n > (E) = n r (E), where all quantities are computed under the same assumed voltage drop.It is then possible to estimate the upper limit for absorption spectra integrations E cut in terms of the IDOS from the right hand side term alone.
FIG.1.Sketch of the layered system considered: a Co/Pd layer system is connected to the left and to the right to Cu leads.The Co and Pd layers, together with few neighboring Cu layers, form the central so-called interaction zone.
FIG. 2. Top: layer (n) resolved profile of the spin magnetic moment m n Φ (in µB) for the investigated Co/Pd bilayer system for three applied voltages: Φ = −0.34(red bars), 0 (black bars) and +0.34 V (blue bars).Bottom: electric field induced contribution to the profile ∆m n Φ = m n Φ − m n 0 (full bars) for the voltage drop Φ = +0.34V with ∆m d n Φ its part due to the d-electrons (shaded bars).

FIG. 3 .
FIG. 3. Layer n resolved Edelstein coefficient p n zz (full bars) together with its part p d,n zz due to the d-electrons (shaded bars).In addition the layer resolved charge conductivity σ n zz

FIG. 4 .
FIG. 4. Electric field induced spin magnetic moment of the d states of Pd, as a function of the applied voltage Φ, averaged over the Pd layers: ∆m d Φ calculated directly from the lesser Green function G < Φ (E) via Eq.(2) (open diamonds) and ∆m XMCD d Φ deduced from the XMCD spectra at the L2,3-edges via the modified sum rules (see Eq. (4)) (full circles).
Appendix C: Absorption spectra via retarded GF and via greater GF

FIG. 5 .
FIG.5.Calculated XAS spectra for the L2,3-edge of pure Cu based on the standard expression for µ λ (ω) shown using the retarded GF (full black curve) and the generalized expression in Eq. (3) in the main text shown using the greater GF (dashed red curve).
F and μL 2(3) (ω) stands for the polarization averaged absorption coefficient for the L 2 (L 3 ) edge.All energy integrals in Eq. (4) run up to the cutoff energy E cutoff for which the integrated d-like DOS covers 10 dstates; i.e. ´Ecutoff dE is the number of d-like holes derived from the d-like DOS n d (E) above the Fermi level E −∞ n d (E) dE = 10.Equation (