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

Current extraction

To extract the in-plane sheet current IC flowing inside the sample from the measured terahertz signal STHz, we first measure our set-up response function HSE 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 Eref, HSE 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 HSE (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 equation16.

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

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


Here, −e is the electron charge and the sample impedance Z(ω) is given by Z0/[1 + nsub + Z0(ω)] with the free-space impedance Z0, the substrate refractive index nsub ≈ 2 (ref. 56) and the metal-stack thickness d = dFM + dPM. 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 IC 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, SiO2 (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 dW < 10 nm, with an α-phase content that grows with dW and dominates for dW > 10 nm (ref. 9).

Estimate of electronic temperatures

We calculate the electronic temperature increase ΔTe0 upon pump-pulse absorption by

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


Here, T0 = 300 K is the ambient temperature, Fl is the absorbed fluence in the respective layer l = FM or PM (Supplementary Table 1), d is the layer thickness and γTe is the specific electronic heat capacity at electronic temperature Te 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}}}}$$


with the total absorbed fluence Ftot 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 tmax of the transient charge current IC(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 IC(tmax) (Fig. 4c) and the relative area of IC(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 IC(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 ΔLk0 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 ΔLk(z, t) = ΔLk0 δ(z − vkzt), where vkz is the z component of the wave-packet group velocity. Note that we restrict ourselves to k values with non-negative vkz 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).$$


Assuming that ΔLk0 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 vkz ≥ 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)$$


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})$$


is the L weight of the z-axis group velocity vz. In equation (7), we phenomenologically account for the relaxation of the ballistic current with time constant τ by introducing the factor et/τ. 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).$$


To determine a plausible shape of w(v), we note that ΔLk0 is non-zero within several 0.1 eV around the Fermi energy EF owing to the width of the photoexcited and rapidly relaxing electron distribution33. Assuming a spherical Fermi surface and isotropic ΔLk0, we have ΔLk0 δ(Ek − EF) and vkz = vFcosθ, where Ek is the band structure, vF 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),$$


where Θ is the Heaviside step function. In other words, all velocities \({v}_{z}={v}_{{\rm{F}}}{{\cos }}\theta\) from 0 to vF 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)$$


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, vL is their group velocity and τ is their velocity relaxation time. To determine the current density, we apply Fick’s law23, jL = −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).$$


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 code2, which implements the full-potential linearly augmented plane-wave (FLAPW) method59. The exchange–correlation effect is included in the scheme of the generalized gradient approximation by using the Perdew–Burke–Ernzerhof functional60. The lattice parameter of the cubic unit cell is set to 5.96a0, where a0 is the Bohr radius. For the muffin-tin potential, we set RMT = 2.5a0 for the radius and lmax = 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 code61. 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} }},$$


where vα and Lβ are the α component of the velocity and the β component of the orbital-angular-momentum operators, respectively. The |nk〉 is the eigenstate of the Hamiltonian with energy Enk 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}}}}},$$


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}}},$$


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 kB 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 })$$


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 }$$


where V is the volume of the system and Γ = 25 meV is a phenomenological broadening parameter. The LCC 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 },$$


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 LCC 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.