# Time-domain observation of ballistic orbital-angular-momentum currents with giant relaxation length in tungsten

### Current extraction

To extract the in-plane sheet current *I*_{C} flowing inside the sample from the measured terahertz signal *S*_{THz}, we first measure our set-up response function *H*_{SE} by having a reference electro-optic emitter (50 μm GaP on a 500 μm glass substrate) at the same position as the sample, which yields a reference terahertz signal \(S_{\mathrm{THz}}^{\mathrm{ref}}\). By calculating the emitted terahertz electric field from that reference emitter *E*^{ref}, *H*_{SE} is determined by solving the convolution \({{{S}}}_{{\rm{THz}}}^{{\rm{ref}}}\left({{t}}\right)=\left({{{H}}}_{{{SE}}} * {{{E}}}^{{\rm{ref}}}\right)\left({{t}}\right)\) (equation (2)) for *H*_{SE} (ref. ^{43}). Further measured inputs for this calculation are the excitation spot size with a full-width at half-maximum of 22 μm, the excitation pulse energy of 1.9 nJ and a Fourier-transform-limited pump pulse with a spectrum centred at 800 nm and a full-width at half-maximum of 110 nm. We perform the deconvolution directly in the time domain by recasting it as a matrix equation^{16}.

Next, the electric field *E* directly behind the sample is obtained from the recorded terahertz signal *S*_{THz} with the help of the derived function *H*_{SE} by solving the analogous equation *S*_{THz}(*t*) = (*H*_{SE} * *E*)(*t*) for *E*. Finally, the sheet charge current (Supplementary Table 1) as shown in Fig. 3 is derived from a generalized Ohm’s law^{28}, which in the frequency domain at frequency *ω*/2π reads

$$E\left(\omega \right)={eZ}\left(\omega \right){I}_{{{C}}}\left(\omega \right).$$

(3)

Here, −*e* is the electron charge and the sample impedance *Z*(*ω*) is given by *Z*_{0}/[1 + *n*_{sub} + *Z*_{0}*dσ*(*ω*)] with the free-space impedance *Z*_{0}, the substrate refractive index *n*_{sub} ≈ 2 (ref. ^{56}) and the metal-stack thickness *d* = *d*_{FM} + *d*_{PM}. The measured mean sample conductivity *σ* (Supplementary Table 1) is approximately frequency-independent due to the large Drude scattering rate (Supplementary Fig. 2). To enable a comparison of terahertz currents from different samples, we normalize *I*_{C} by the absorbed fluence in the FM. The data shown in Fig. 3 were obtained in a dry-air atmosphere.

### Sample preparation

The FM|PM samples (FM = Ni and Py, PM = Pt, Ti, Cu and W) are fabricated on glass substrates (thickness 500 μm) or thermally oxidized Si substrates (625 μm) by radio-frequency magnetron sputtering under an Ar atmosphere of 6N purity. The sample structure and thickness are described in Supplementary Table 1. For the sputtering, the base pressure in the chamber is lower than 5 × 10^{−7} Pa. To avoid oxidation, SiO_{2} (thickness 4 nm) is sputtered on the surface of the films. All sputtering processes are performed at room temperature. The W films are predominantly in the β-phase for *d*_{W} < 10 nm, with an α-phase content that grows with *d*_{W} and dominates for *d*_{W} > 10 nm (ref. ^{9}).

### Estimate of electronic temperatures

We calculate the electronic temperature increase Δ*T*_{e0} upon pump-pulse absorption by

$$\Delta {T}_{{\rm{e}}0}=\sqrt{{T}_{0}^{2}+\frac{2{F}_{l}}{\gamma d}}-{T}_{0}.$$

(4)

Here, *T*_{0} = 300 K is the ambient temperature, *F*_{l} is the absorbed fluence in the respective layer *l* = FM or PM (Supplementary Table 1), *d* is the layer thickness and *γT*_{e} is the specific electronic heat capacity at electronic temperature *T*_{e} with *γ* = 300 J m^{−3} K^{−2} for W, 320 J m^{−3} K^{−2} for Ni, 330 J m^{−3} K^{−2} for Ti and 90 J m^{−3} K^{−2} for Pt (ref. ^{57}).

To obtain the absorbed fluences in each layer, we note that the pump electric field is almost constant throughout the sample (Supplementary Fig. 9). Therefore, the local pump absorption scales solely with the imaginary part Im*ε* of the dielectric function *ε* at a wavelength of 800 nm, which equals 22.07 for Ni, 9.31 for Pt, 19.41 for Ti and 19.71 for W (ref. ^{58}). Consequently, the absorbed fluence is determined by

$${F}_{l}={F}_{{\rm{tot}}}\frac{{d}_{l}{\rm{Im}}{\varepsilon }_{l}}{{d}_{{\rm{FM}}}{\rm{Im}}{\varepsilon }_{{\rm{FM}}}+{d}_{{\rm{PM}}}{\rm{Im}}{\varepsilon }_{{\rm{PM}}}}$$

(5)

with the total absorbed fluence *F*_{tot} that is obtained from the absorbed pump power (Supplementary Table 1) and the beam size on the sample (as described previously).

### Experimental error estimation

The error bars for the maximum position *t*_{max} of the transient charge current *I*_{C}(*t*) (Fig. 4b) are estimated as ±20% of the read-off delay value (circular markers in Fig. 4a), but no less than 5 fs. The uncertainty in the relative amplitude of *I*_{C}(*t*_{max}) (Fig. 4c) and the relative area of *I*_{C}(*t*) (Fig. 4e) is, respectively, estimated as ±20% and ±10%, both reflecting the typical signal-to-noise ratio of the extracted current traces (Fig. 4a). The error bars for the full-width of *I*_{C}(*t*) at half-maximum (Fig. 4d) are obtained from the uncertainty of the delay (Fig. 4b) with subsequent multiplication by \(\sqrt{2}\), which accounts for the error propagation of a difference of two quantities.

### Model of *L* transport

To model the ballistic current in the PM, we assume that a δ(*t*)-like transient *L* accumulation in the FM generates an electronic wave packet, which has orbital angular momentum Δ*L*_{k0} along the direction of the FM magnetization **M** and mean wave vector **k** in the PM right behind the FM/PM interface, that is, at *z* = 0^{+} (Fig. 5a).

In the case of purely ballistic transport, this wave packet propagates into the PM bulk according to Δ*L*_{k}(*z*, *t*) = Δ*L*_{k0} δ(*z* − *v*_{kz}*t*), where *v*_{kz} is the *z* component of the wave-packet group velocity. Note that we restrict ourselves to **k** values with non-negative *v*_{kz} values. The total pump-induced *L* current density flowing into the depth of the PM is for *z* > 0 given by the sum

$${r}_{z}\left(t\right)=\mathop{\sum }\limits_{{\bf{k}},{v}_{{\bf{k}}z}\ge 0}\Delta {L}_{{\bf{k}}0}{v}_{{\bf{k}}z}\updelta\left(z-{v}_{{\bf{k}}z}t\right).$$

(6)

Assuming that Δ*L*_{k0} arises from states not too far from the Fermi energy, the summation of equation (6) is approximately proportional to an integration over the Fermi surface parts with *v*_{kz} ≥ 0. One obtains

$${r}_{z}(t)={\mathrm{{e}}}^{-t/\tau }{\int }_{0}^{\infty }{\rm{d}}{v}_{z}w({v}_{z}){v}_{z}\updelta (z-{v}_{z}t)$$

(7)

where *z* > 0, and

$$w({v}_{z})=\sum _{{\bf{k}},{v}_{{\bf{k}}z}\ge 0}\Delta {L}_{{\bf{k}}0}\updelta ({v}_{{\bf{k}}z}-{v}_{z})$$

(8)

is the *L* weight of the *z*-axis group velocity *v*_{z}. In equation (7), we phenomenologically account for the relaxation of the ballistic current with time constant *τ* by introducing the factor e^{−t/τ}. Performing the integration of equation (7) yields

$${r}_{z}\left(t\right)=\frac{{{\rm{e}}}^{-t/\tau }}{t}\frac{z}{t}w\left(\frac{z}{t}\right).$$

(9)

To determine a plausible shape of *w*(*v*), we note that Δ*L*_{k0} is non-zero within several 0.1 eV around the Fermi energy *E*_{F} owing to the width of the photoexcited and rapidly relaxing electron distribution^{33}. Assuming a spherical Fermi surface and isotropic Δ*L*_{k0}, we have Δ*L*_{k0} ∝ δ(*E*_{k} − *E*_{F}) and *v*_{kz} = *v*_{F}cos*θ*, where *E*_{k} is the band structure, *v*_{F} is the Fermi velocity and *θ* is the angle between **k** and the *z* axis. Therefore, after turning equation (8) into an integral, the integrand \(\Delta {L}_{{\bf{k}}0}\updelta\left({v}_{{\rm{F}}}{{\cos }}\theta -{v}_{z}\right)\) \({{\rm{d}}^{3}{\bf{k}}}\) becomes proportional to \(\updelta\left({v}_{{\rm{F}}}{{\cos }}\theta -{v}_{z}\right)\) \({\rm{d}}{{\cos }}\theta\) in spherical coordinates, leading to

$$w\left({v}_{z}\right)\propto {\int }_{0}^{1}{\rm{d}}{{\cos }}\theta\ \updelta\left({v}_{{\rm{F}}}{{\cos }}\theta -{v}_{z}\right)\propto \varTheta \left({v}_{{\rm{F}}}-{v}_{z}\right),$$

(10)

where Θ is the Heaviside step function. In other words, all velocities \({v}_{z}={v}_{{\rm{F}}}{{\cos }}\theta\) from 0 to *v*_{F} have equal weight.

In the case of purely diffusive transport, we use the *L* diffusion equation for *μ*_{L} (ref. ^{23}). With a localized accumulation *μ*_{L}(*z*, *t*) ∝ δ(*z*) at time *t* ≈ 0, the accumulation disperses according to the well-known solution

$${\mu }_{L}\left(z,t\right)\propto \frac{1}{\sqrt{Dt}}{{\exp }}\left(-\frac{{z}^{2}}{4Dt}\right)$$

(11)

for *z* > 0 and *t* > 0. Here, *D* is the diffusion coefficient that equals \({v}_{L}^{2}\tau /3\) in the case of a spherical **k**-space surface carrying the *L* wave packets, *v*_{L} is their group velocity and *τ* is their velocity relaxation time. To determine the current density, we apply Fick’s law^{23}, *j*_{L} = −*D*∂*μ*_{L}/∂*z* to equation (11) and obtain

$${r}_{z}\left(t\right)\propto \varTheta \left(t\right)\frac{z}{t}\frac{1}{\sqrt{Dt}}{{\exp }}\left(-\frac{{z}^{2}}{4Dt}\right).$$

(12)

### Ab initio estimate of the *L* velocity

The orbital velocity is estimated for the bulk W in a body-centred cubic structure. The ab initio self-consistent calculation of the electronic states is performed within density functional theory by using the FLEUR code^{2}, which implements the full-potential linearly augmented plane-wave (FLAPW) method^{59}. The exchange–correlation effect is included in the scheme of the generalized gradient approximation by using the Perdew–Burke–Ernzerhof functional^{60}. The lattice parameter of the cubic unit cell is set to 5.96*a*_{0}, where *a*_{0} is the Bohr radius. For the muffin-tin potential, we set *R*_{MT} = 2.5*a*_{0} for the radius and *l*_{max} = 12 for the maximum of the harmonic expansion. Further, we set the plane-wave cut-offs for the interstitial region to \(4.0{a}_{0}^{-1}\), \(10.1{a}_{0}^{-1}\) and \(12.2{a}_{0}^{-1}\) for the basis set, the exchange–correlation functional and the charge density, respectively. For the **k**-points, a 16 × 16 × 16 Monkhorst–Pack mesh is defined.

From the converged electronic structure, we obtain the maximally localized Wannier functions by using the WANNIER90 code^{61}. We use 18 Wannier states with \(s,{p}_{x},{p}_{y},{p}_{z},{d}_{{z}^{2}},{d}_{{x}^{2}-{y}^{2}},{d}_{xy},{d}_{yz},{d}_{zx}\) symmetries for spin up and down as the initial guess. The maximum of the inner (frozen) energy window is set 5 eV above the Fermi energy for the disentanglement, and the outer energy window is defined by the minimum and maximum energies of the 36 valence states obtained from the FLAPW calculation. The Hamiltonian, position, orbital-angular-momentum operators and spin-angular-momentum operators are transformed from the FLAPW basis into the maximally localized Wannier function basis. The resulting electronic band structure and texture of the orbital-angular-momentum operator **L** are displayed in Supplementary Fig. 9.

From this realistic tight-binding model, the orbital-momentum-weighted velocity averaged over the Fermi surface (FS) is calculated by

$${\langle {v}_{\alpha }{L}_{\beta }\rangle }_{{\rm{FS}}}=\frac{{\sum }_{n{\bf{k}}}{\,f}_{n{\bf{k}}}^{\,{\prime} }\langle n{\bf{k}}|({v}_{\alpha }{L}_{\beta }+{L}_{\beta }{v}_{\alpha })/2|n{\bf{k}}\rangle }{{\sum }_{n{\bf{k}}}{\,f}_{n{\bf{k}}}^{\,{\prime} }},$$

(13)

where *v*_{α} and *L*_{β} are the *α* component of the velocity and the *β* component of the orbital-angular-momentum operators, respectively. The |*n***k〉** is the eigenstate of the Hamiltonian with energy *E*_{n}_{k} and band index *n*, and \({f}_{n{\bf{k}}}^{\,{\prime} }\) is the energy derivative of the Fermi–Dirac distribution function. To polarize *L*_{β}, we add a small orbital Zeeman coupling along the *β* direction to the bare Hamiltonian. We confirm that the result of equation (13) changes by less than 1% when the orbital Zeeman splitting is increased from 10 meV to 30 meV. The **k**-space integrals in equation (13) are performed on a 256 × 256 × 256 mesh.

The orbital velocity is estimated by

$${\left\langle {v}_{\alpha }^{{L}_{\beta }}\right\rangle }_{{\rm{FS}}}=\frac{{\left\langle {v}_{\alpha }{L}_{\beta }\right\rangle }_{{\rm{FS}}}}{\sqrt{{\left\langle {L}_{\beta }^{2}\right\rangle }_{{\rm{FS}}}}},$$

(14)

where \({\langle {L}_{\beta }^{2}\rangle }_{{\rm{FS}}}\) is obtained by equation (13), but with *v*_{α} replaced by *L*_{β}. The result is shown in Supplementary Fig. 10.

### Details of ab initio LCC calculations

A thin W stack of 19 body-centred cubic (110) atomic layers is calculated by the self-consistent ab initio method using the same FLAPW parameters as for the bulk calculation given previously, except for the mesh of **k**-points for which we use a 24 × 24 Monkhorst–Pack mesh. For Wannierization, we obtain 342 maximally localized Wannier functions, starting from Wannier states with \(s,{p}_{x},{p}_{y},{p}_{z},{d}_{{z}^{2}},{d}_{{x}^{2}-{y}^{2}},{d}_{xy},{d}_{yz},{d}_{zx}\) symmetries for spin up and down as the initial guess. We define the maximum of the inner window at 2 eV above the Fermi energy.

From the Hamiltonian of the W thin film, the **k**-space orbital-angular-momentum texture at the Fermi surface is obtained by

$${\langle {{\bf{L}}}_{{\rm{top}}}\rangle }_{{\rm{FS}}}({\bf{k}})=-4{k}_{{\rm{B}}}T\sum _{n}{f}_{n{\bf{k}}}^{\,{\prime} }{\langle {{\bf{L}}}_{{\rm{top}}}\rangle }_{n{\bf{k}}},$$

(15)

where \({\langle {{\bf{L}}}_{{\rm{top}}}\rangle }_{n{\bf{k}}}=\langle n{\bf{k}}|{{\bf{L}}}_{{\rm{top}}}|n{\bf{k}}\rangle\) is the expectation value of the orbital angular momentum for the two atoms on the top surface, *T* = 300 K is temperature and *k*_{B} is the Boltzmann constant.

To calculate the charge current due to LCC, we consider the orbital-dependent chemical potential

$$\frac{{\varepsilon }_{\beta \gamma }}{2}({r}_{\beta }{L}_{\gamma }+{L}_{\gamma }{r}_{\beta })$$

(16)

as a perturbation. Here, *L*_{γ} is the *γ* component of the orbital-angular-momentum operator, and *r*_{β} is the *β* component of the position operator, which is well-defined along the *z* axis (Fig. 1), and \({\varepsilon }_{\beta \gamma }\) can be interpreted as an orbital-dependent electric field. The charge-current density along the *α* direction is given by the Kubo formula

$$\langle {\,j}_{\alpha }\rangle =-\frac{e}{V}\sum _{{\bf{k}}nn^{\prime} }({\,f}_{n{\boldsymbol{k}}}-{f}_{n^{\prime} {\bf{k}}}){\rm{Re}}\frac{\langle n{\bf{k}}|{v}_{\alpha }|n^{\prime} {\bf{k}}\rangle \langle n^{\prime} {\bf{k}}|V|n^{\prime} {\bf{k}}\rangle }{{E}_{n{\bf{k}}}-{E}_{n^{\prime} {\bf{k}}}+{\rm{i}}\varGamma }$$

(17)

where *V* is the volume of the system and *Γ* = 25 meV is a phenomenological broadening parameter. The *LC*C response is characterized by the tensor

$${\sigma }_{{{{LC}}}{\rm{C}},\alpha \beta }^{{L}_{\gamma }}=-\frac{e}{2V}\mathop{\sum }\limits_{{\bf{k}}n{n}^{{\prime} }}\left({\,f}_{n{\bf{k}}}-{f}_{{n}^{{\prime} }{\bf{k}}}\right){\rm{Re}}\frac{\left\langle n{\bf{k}},|,{v}_{\alpha },|,{n}^{{\prime} }{\bf{k}}\right\rangle \left\langle {n}^{{\prime} }{\bf{k}},|,\left({r}_{\beta }{L}_{\gamma }+{L}_{\gamma }{r}_{\beta }\right),|,{n}^{{\prime} }{\bf{k}}\right\rangle }{{E}_{n{\bf{k}}}-{E}_{{n}^{{\prime} }{\bf{k}}}+{\rm{i}}\varGamma },$$

(18)

which relates \(\langle {j}_{a}\rangle\) and *ε*_{βγ} by \(\langle {j}_{\alpha }\rangle ={\sigma }_{{LC\mathrm{C}},\alpha \beta }^{{L}_{\gamma }} \ \epsilon_{\beta \gamma}\). The **k**-space integral is performed on a 400 × 400 mesh. The *z*-resolved *LC*C response is shown in Supplementary Fig. 11.

### Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.