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. ⁴³). 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¹⁶.

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²⁸, which in the frequency domain at frequency ω/2π reads

Here, −e is the electron charge and the sample impedance Z(ω) is given by Z₀/[1 + n_sub + Z₀dσ(ω)] with the free-space impedance Z₀, the substrate refractive index n_sub ≈ 2 (ref. ⁵⁶) 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⁻⁷ Pa. To avoid oxidation, SiO₂ (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. ⁹).

Estimate of electronic temperatures

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

Here, T₀ = 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⁻³ K⁻² for W, 320 J m⁻³ K⁻² for Ni, 330 J m⁻³ K⁻² for Ti and 90 J m⁻³ K⁻² for Pt (ref. ⁵⁷).

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. ⁵⁸). Consequently, the absorbed fluence is determined by

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_kzt), 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

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

where z > 0, and

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

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³³. Assuming a spherical Fermi surface and isotropic ΔL_k0, we have ΔL_k0 ∝ δ(E_k − E_F) and v_kz = v_Fcosθ, 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

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. ²³). With a localized accumulation μ_L(z, t) ∝ δ(z) at time t ≈ 0, the accumulation disperses according to the well-known solution

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²³, j_L = −D∂μ_L/∂z to equation (11) and obtain

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², which implements the full-potential linearly augmented plane-wave (FLAPW) method⁵⁹. The exchange–correlation effect is included in the scheme of the generalized gradient approximation by using the Perdew–Burke–Ernzerhof functional⁶⁰. The lattice parameter of the cubic unit cell is set to 5.96a₀, where a₀ is the Bohr radius. For the muffin-tin potential, we set R_MT = 2.5a₀ 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⁶¹. 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

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 E_nk 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

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

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

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

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

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.