Unidirectional magnetoresistance and spin-orbit torque in NiMnSb

Spin-dependent transport phenomena due to relativistic spin-orbit coupling and broken space-inversion symmetry are often difficult to interpret microscopically, in particular when occurring at surfaces or interfaces. Here we present a theoretical and experimental study of spin-orbit torque and unidirectional magnetoresistance in a model room-temperature ferromagnet NiMnSb with inversion asymmetry in the bulk of this half-heusler crystal. Besides the angular dependence on magnetization, the competition of Rashba and Dresselhaus-like spin-orbit couplings results in the dependence of these effects on the crystal direction of the applied electric field. The phenomenology that we observe highlights potential inapplicability of commonly considered approaches for interpreting experiments. We point out that, in general, there is no direct link between the current-induced non-equilibrium spin polarization inferred from the measured spin-orbit torque and the unidirectional magnetiresistance. We also emphasize that the unidirectional magnetoresistance has not only longitudinal but also transverse components in the electric field -- current indices which complicates its separation from the thermoelectric contributions to the detected signals in common experimental techniques. We use the theoretical results to analyze our measurements of the on-resonance and off-resonance mixing signals in microbar devices fabricated from an epitaxial NiMnSb film along different crystal directions. Based on the analysis we extract an experimental estimate of the unidirectional magnetoresistance in NiMnSb.


I. INTRODUCTION
Anisotropic magnetoresistance (AMR) is an example of a relativistic transport phenomenon in ferromagnets with a long history. McGuire and Potter 1 provided a phenomenological model that fully explained the angular dependence of AMR. The model was based on a general argument that the magnetisation dependent conductivity tensor σ ij (m), like any other observable, 2 reflects the symmetry of the ferromagnetic crystal. This means that if a symmetry operation belongs to the crystallographic point group, the conductivity tensor is left invariant under the same symmetry operation. As an example, when an electrical current flows in an isotropic (polycrystalline) ferromagnet, the magnetisation direction defines the only axis of rotational symmetry, which results in a cos(2θ) angular dependence of AMR, with θ being the angle between the magnetisation and the current.
While AMR can be observed in any magnetic system, ferromagnetic conductors with an inversion asymmetric crystal structure have been more recently identified as a fruitful platform for discovering and utilizing a range of new relativistic spintronics phenomena. [3][4][5][6][7][8] Apart from the exchange field B ex , which splits the electronic structure in majority and minority spin bands, the inversion symmetry breaking in the lattice together with spinorbit coupling introduces an additional splitting ∆H SO = B SO (k) · s. Here B SO (k) is an effective magnetic field that depends on the crystal momentum k, while s represents the spinpolarisation. The key property of B SO (k) is that it is odd in k and thus it results in opposite spin splitting for opposite k. In inversion asymmetric strained zinc-blende semiconductors like GaAs, B SO (k) is a combination of Rashba and Dresselhaus symmetry terms, B R SO (k) ∝ (k y , −k x ) and B D SO (k) ∝ (k x , −k y ). A direct manifestation of this splitting is found in the spin-orbit torque (SOT) 9 in, e.g., (Ga,Mn)As that emerge when an electrical current is applied to this inversion-asymmetric ferromagnetic semiconductor. [3][4][5][6] In metallic systems such as half-heusler NiMnSb studied here, B SO (k) is not described by the simple linear-ink Rashba and Dresselhaus form. SOT, nevertheless, contains Rashba and Dreselhaus-like terms analogous to those observed in the zinc-blende semiconductors as those reflect the crystal symmetry which is the same in zinc-blende and half-heusler crystals. 8 SOT in bulk non-centrosymmetric crystals is associated with a current-induced nonequilibrium spin polarisation δs SO , an effect often referred to as the Edelstein effect or the inverse spin-galvanic effect (iSGE). [9][10][11][12][13][14] When δs SO is perpendicular to the magnetization, it will exert a torque on it via the exchange coupling. The component of δs SO that is parallel to B ex does not lead to a magnetization torque and is, therefore, transparent to any experimental method that relies on driving the magnetisation out of equilibrium by SOT. 5, 9,15 Besides SOT, B SO (k) together with B ex can also lead to magneto-transport terms that are second order in the applied electric field and, unlike the first-order AMR, are odd under the magnetization. The unidirectional magnetoresistance (UMR) is an example of such a secondorder magneto-transport effect that was previously reported in non-centrosymmetric ferromagnet/paramagnet bilayers or bulk ferromagnets. 16,17 The origin of UMR was considered to be closely connected to the phenomenology of the giant magnetoresistance (GMR). [17][18][19][20] While in the GMR multilayer, the fixed reference ferromagnetic layer acts as an external source of spin δs, in UMR this is replaced with δs SO generated internally by the spin-orbit coupling. As for GMR, spin-dependent scattering within the ferromagnet results in this scenario in a different resistance depending on the orientation of δs SO relative to the magnetization of the probed ferromagnet. Moreover, the accumulated spin can introduce a proportional change in the exchange splitting of the bands, further affecting the conductivity by influencing spin transmission of minority and majority spins. Another UMR mechanism considers that the magnons' population of the ferromagnet is increased or decreased depending on the orientation of δs SO relative to the magnetisation. This leads to a change in the average magnetisation, which also results in a net change of the magnetoresistance.
Although it is possible to experimentally distinguish these different contributions for their particular dependence on current density and magnetic field, 20,21 they all share a common origin in the component of δs SO collinear with the ferromagnet's magnetisation.
In Sec. II we report our symmetry analysis and ab initio calculations of current induced spin polarization and UMR in NiMnSb. We use this model system to highlight potential misconceptions when inferring these quantities from experiment. In Sec. III we then discuss our measurements in NiMnSb microbars patterned along different crystal directions, and in Sec. IV we summarize our main results.

YSIS OF EXPERIMENTS
A. Current-induced spin polarization and spin-orbit torque The tetragonal distortion of the non-centrosymmetric cubic unit cell of NiMnSb, induced by the lattice mismatch with the substrate, lowers its symmetry to a −42m symmetry point group and results in a Dresselhaus-like k-linear term of the spin-orbit coupling. Experiments in NiMnSb epilayers show an additional Rashba-like k-linear term of the spin-orbit coupling which we model by introducing a shear strain. 8 This lowers the symmetry further to a point group mm2. When an electrical current is passed in the plane perpendicular to the growth direction, carriers acquire a non-equilibrium spin polarization δs SO , which in general can be decomposed into a component that is parallel to the in-plane magnetisation of the NiMnSb film, δs SO , and a component that is perpendicular to it, δs ⊥ SO . This includes both intrinsic and extrinsic contributions.
ization as a function of the directions of magnetization and applied electric field (see Appendix C for the description of the numerical method). We consider here only the in-plane components of current-induced spin-polarization since this is the component relevant for the UMR and the field-like torque. For these calculations we used the relaxation time τ chosen so that the theoretical conductivity matches the experimental value of 3.3 × 10 4 Scm −1 . The tetragonal and shear strains are chosen to make the Rashba and Dresselhaus contributions of comparable strength, 8 resulting in a non-trivial dependence of δs SO on the crystal direction of the applied electric field. In Fig. 1 we plot the perpendicular (δs ⊥ SO ) and parallel (δs SO ) components of the current-induced non-equilibrium spin polarization. The dependence of δs ⊥ SO ∼ cos(θ − θ SO ) on the magnetization angle θ corresponds to what we would expect to find for a field-like SOT generated by a magnetisation-independent effective spin-orbit field, Remarkably, when we include δs SO , total δs SO becomes magnetization-dependent (see Figs. 1). We emphasize that this dependence of δs SO (h SO ) on the magnetization angle would be invisible when measuring SOT, which only depends on δs ⊥ SO . SOT can be reliably obtained from on-resonance mixing-signal measurements since the detected signal can be directly attributed to the precessing magnetization driven by the SOT. Our calculations in Fig. 1 demonstrate, however, that a simple harmonic dependence of δs ⊥ SO on the magnetization angle, when inferred from the SOT measurement, does not imply that the total δs SO is constant and that its parallel component δs SO is a 90 • -phase shifted replica of the perpendicular component. Therefore, if considering δs SO as the driving mechanism behind UMR, SOT measurements cannot be, in general, used to quantify experimentally δs SO in a given structure.

B. Unidirectional magnetoresistance
When writing the total current density j up to the second order in the applied electric filed E as, UMR has been associated with the longitudinal component of the ξ ijk transport coefficient. 16,17 Formally, ξ ijk is described by the second order quantum mechanical Kubo formula. However, finding accurate solutions of the formula is a major challenge, in particular in the presence of electron scattering. Here we analyze ξ ijk using a semiclassical Boltzmann approximation, and v n (k) = 1 ∂ n ∂k , f 0 is the Fermi-Dirac equilibrium distribution function, and n (k) is the band energy (see Appendix A for the derivation of this formula).
Since the group velocity v(k) is odd under space inversion, the second-order term will vanish in crystals that have inversion symmetry. Moreover, it will also vanish in nonmagnetic crystals since v(k) is odd under time-reversal. Furthermore, similarly to the anomalous Hall effect in coplanar magnetic systems, the second-order term will vanish in the absence of spin-orbit coupling, as the system will then be invariant under combined spin rotation and time-reversal symmetry. 22 ξ ijk will thus be present in magnetic crystals with broken inversion symmetry as a consequence of the spin-orbit coupling. These are the same symmetry requirements as for the existence of SOT.
In Fig. 2  When ferromagnetic resonance is excited, rectification between the microwave current and the oscillating AMR results in a resonance in V dc . 5 The resonance is clearly visible in the 2D plots in Fig. 3b. It shows V dc as a function of the external magnetic field and its direction θ with respect to the current direction, for bars parallel to the [110] and [100] axes. From the analysis of its line-shape (see Ref. 6 for details) we are able to quantify h ⊥ SO . In agreement with previous studies on different systems, 5,6,15 we experimentally identify the θ-dependence of h ⊥ SO for the different crystal directions in which the microbars are patterned, as shown in Fig. 4 In Fig. 3b we notice that the resonance is sitting on a sinusoidally varying background, which we refer to as V BG . This is also shown in We note that these amplitude ratios of V BG , as well as the crystal direction dependent phase shifts, are reproducible in different physical sets of samples patterned along the four crystal directions and do not depend on applied power, as shown in Appendix D.
To interpret the measured V BG we now consider UMR and the thermoelectric contribution, namely the anomalous Nernst effect (ANE). When exciting the system by the applied ac current, an out-of-plane temperature gradient due to Joule heating can result in an electrical signal detected in the sample plane due to ANE. Since the heat deposited by the current scales with the square of the current density, ANE is a second-order effect in the electric field, just like the second order term of the conductivity. Based on a recent experimental measurement 24 of ANE in NiMnSb and our numerical simulation of the heat gradient, we estimate that the contribution to V BG due to ANE is ∼ 0.01 − 0.1 µV per current density of 10 10 Am −2 (see Appendix E). While the magnitude is similar to that of the measured signals in Fig.5, we do not expect a strong dependence of the ANE contribution to V BG on the crystal direction of the applied current. This is because ANE requires only the timereversal symmetry breaking by the magnetization while broken spatial symmetries of the crystal only lead to additional, higher order corrections.
UMR, on the other hand, is generated by the inversion symmetry breaking which is of the combined Rashba-Dresselhaus-like form in our NiMnSb samples. As discussed in the theory section, this leads to a strong crystal direction dependence of the UMR. Since theory suggests that the amplitude ratios for the different crystal directions of SOT and UMR are similar, we can use this as a constraint when fitting the measured V BG data. The results of the fitting shown in Fig. 5 were obtained by assuming fixed amplitude ratios of UMR, corresponding to the measured amplitude ratios of the SOT fields h ⊥ SO , plus an ANE contribution with an amplitude which is independent of the crystal direction. The extracted experimental ANE component is of the same order of magnitude as the above estimate. The fitted UMR contribution is an order of magnitude larger than our theoretical value which we attribute to the crude semiclassical Boltzmann approximation used in the UMR calculations. In the future, more elaborate Kubo formula calculations seem necessary to capture UMR in NiMnSb on the quantitative level. They will also allow for verifying the correspondence between the SOT and UMR amplitude ratios, and by this for more firmly establishing the fitting method we used to separate the UMR and thermoelectric contributions in the measured data.

IV. SUMMARY
Based on our study of ferromagnetic NiMnSb with non-centrosymmetric bulk crystal structure we make the following observations regarding the explored non-equilibrium spin-   directions and the ANE magnitude.

Appendix A: Second-order Boltzmann formula derivation
Here we derive the second-order Boltzmann formula (2). The general form of the Boltzmann formula for a distribution function g(t, r, k) under the assumptions of a stationary and spatially homogeneous g and no magnetic field 25 can be expressed as: where e is the (positive) elementary charge, E is the electric field and dg dt col is the change of the distribution function due to scattering.
We will assume the constant relaxation-time approximation: where τ is the relaxation time and g 0 is the equilibrium distribution function, which for electrons is the Fermi-Dirac distribution function: Within the constant relaxation-time approximation the Boltzmann formula has the following form: To find a solution for g up to second order in E we expand g in powers of E: and insert it into the Boltzmanm formula (A4). Keeping only the terms up to E 2 we find: Since this equation must hold for all E, the coefficients for the E and E 2 terms on both sides of the equation must be the same. Therefore Taking into account the relation we have Electrical current is then given by: We note that this integral is done over the first Brillouin zone or any other unit cell in the reciprocal space. The first order contribution to the current is given by and the second order contribution: Note that in a multi-band system, these expression give a contribution from each individual band and the total current will be a sum over all bands.
The expression for the second order current can be further simplified. We first define a second order conductivity tensor ξ ijl : Alternatively, using Eqs. (A7) and (A8), ξ can be written in the form Using integration by parts here Γ signifies integral over the unit cell boundary and ν is the outward unit normal vector to the boundary. The first term in this relation in fact vanishes. To see that this is the case it is useful to use for the integration a unit cell of the reciprocal space spanned by the reciprocal lattice vectors b 1 , b 2 , b 3 (see Fig. (6)), instead of the first Brillouin zone. The first term in (A22) is then given by sum over 6 surfaces. Since g 0 is a periodic function of k, also ∂g 0 ∂ must be periodic. This means that at the opposite boundaries of the reciprocal unit cell ∂g 0 ∂ is the same. However, since the outward unit normal vector ν is opposite for the opposite boundaries the whole term vanishes.
Combining Eqs. (A21) and (A22) we have To simplify this further we need to show that ξ ijl is symmetric under interchanging any two indices. This is clearly satisfied for ξ a ijl , however, it is less obvious for ξ b ijl . From Eq. (A17), we see that ξ ijl = ξ ilj must hold and thus also ξ b ijl = ξ 2,b ilj should be satisfied. This can be explicitly verified from 1 ∂v l ∂k j = ∂ 2 ∂k l ∂k j , which is symmetric as long as is sufficiently smooth function. To show that the tensors are also symmetric under interchanging the other indices, we consider: Then it must also hold that the ξ ijl and ξ b ijl tensors are symmetric under interchanging indices i and l: We can thus rewrite Eq. (A23) as Therefore, ξ b ijl = −ξ b ijl /2 and we finally find: the Boltzmann formula that we use in our calculations. It is not clear whether this holds for the ξ ijk in general, thus the symmetry analysis here should be taken to refer specifically to the Boltzmann contribution. The method for the symmetry analysis is analogous to the one used in Ref. 26. We have implemented the second order symmetry analysis in the open source code Symmetr. 27 All the results shown here are given in a cartesian coordinate system described in Fig. 6(b). Since in the experiment the magnetization always lies in the [001] plane, we consider here only magnetization in this plane. In Table I we give the general shape of the ξ ijk for general direction of magnetization within this plane as well as for the [110] and [1-10] directions where the symmetry is higher. To describe the dependence of the second-order currents on magnetization, it is useful to expand the ξ ijk tensor in powers of the magnetization. We consider only the lowest order Here,Ŝ is the spin operator, u n (k) are the Bloch functions of a band n, k is the Bloch wave vector, ε n (k) is the band energy, E F is the Fermi energy,v is the velocity operator and Γ is a quasiparticle broadening parameter that describes the disorder strength. This The DFT calculations utilized 12x12x12 k-points and the GGA-PBE potential. The current-induced spin-polarization and the second-order conductivity response calculations use a 400x400x400 k-mesh, which we have confirmed to be sufficient for good convergence.
To estimate the value of Γ we calculate the first order conductivity using a conductivity formula analogous to Eq. (C1) and choose the Γ so that the conductivity matches the experimental conductivity. This corresponds to Γ ≈ 0.05 eV or alternatively to τ ≈ 6.6 fs.
In previous calculations of SOT in NiMnSb, a Γ = 0.036 eV was used, 8 because the samples used in those experiments had somewhat larger conductivity. We note that our calculations of the current-induced spin-polarization were done for the Γ = 0.036 eV value and have been afterwords rescaled to Γ = 0.05 eV, assuming 1/Γ scaling. This is quite accurate since for theses value of Γ the formula is very close to the Boltzmann formula and thus has 1/Γ scaling.

Appendix D: Power dependence and reproducibility
In Fig. 7 we show that, as expected, the background voltage V BG scales linearly with power and that the angular dependence of V BG is independent of the power. In Fig. 8 we show that the results are reproducible between different bars patterned on the same chip. is shown in Figure 9 and consists of the NiMnSb wire (thickness = 37 nm, width = 1 µm) deposited on a (In,Ga)As mesa (thickness = 200 nm, width = 1 µm). There is also an InP substrate (thickness = 1 µm, width= 2 µm), a thin MgO capping layer (thickness = 5 nm, width = 1 µm), and a He atmosphere included in the simulation. The reference temperature, T ref , is set to 300 K (before the Joule heating takes place).
The current density applied to NiMnSb is 10 10 or 10 11 Am −2 . At the boundary of the simulated area the temperature is fixed to T ref or a thermal insulation is assumed (except the bottom boundary again fixed to T ref ). These two types of boundary condition correspond to very efficient cooling (transfer of heat to the surroundings) or to extremely poor cooling, respectively. The real system would fall between these two limiting cases. Figures 10 and 11 show the simulated steady-state temperature distribution for the case with boundaries fixed to T ref and with insulating boundaries, respectively.
The former case (Fig. 10) is simulated with applied current density of 10 10 and 10 11 Am −2 but we show the temperature profile only for 10 10 Am −2 . The 10 times larger current density results visually in the same temperature profile but the maximum temperature increase, T-T ref , is 100 times larger.
The latter case (Fig. 11) is simulated only with applied current density of 10 10 Am −2 and all the heat is dissipated only via the bottom boundary -through the substrate. The outof-plane temperature gradient is evaluated along a vertical cut-line (along the z-coordinate) running through the middle of the mesa. Figure 12 shows the temperature increase, T-T ref , along the cut-line within the NiMnSb layer (1200 nm to 1237 nm) for the three cases simulated.
As expected, the case with insulating boundaries shows a larger increase in temperature but the dependence on the z-coordinate is the same. The case with 10 times larger applied current shows a 100 times larger increase of temperature which is due to the Joule heating scaling with the square of the applied current density. Finally, Figure 13 shows the out-of-plane temperature gradient generating ANE in the NiMnSb wire. It is evaluated simply as a numerical derivative of the temperature given in Figure 12 for the three cases. Note that the gradient decreases linearly in the NiMnSb layer towards the top surface. The efficiency of the cooling (insulating boundaries or fixed temperature) does not affect the gradient so the approximation of the cooling mechanism assumed in our model should not significantly compromise the validity of our numerical results. The main result of the FEM simulation relevant to the experimental device is that the average out-of-plane temperature gradient is of the order of 10 4 Km −1 for current density of 10 10 Am −2 . Our FEM results depend on the material parameters used. We have measured the electrical conductivity of the individual layers at room temperature. The room temperature thermal conductivity and heat capacity parameters of our films are estimated based on literature as listed in Table II. In case of NiMnSb they are estimated based on related materials. 35,36 The thermal conductivity is underestimated, considering the metallic character of the film, so we report an upper estimate of the thermal gradient. The results are obtained for a steady state so there is no significant dependence on the heat capacity or mass density. (We estimate the mass density based on the molar masses and unit cell volume.)