Simulation of Hot-Carrier Dynamics and Terahertz Emission in Laser-Excited Metallic Bilayers

TL;DR

This paper introduces a multiscale simulation model for understanding and optimizing terahertz emission from laser-excited metallic bilayers, specifically Fe/Pt heterostructures, by modeling spin-current dynamics and electromagnetic wave propagation.

Contribution

It presents a novel multiscale approach combining microscopic charge carrier transport with electromagnetic modeling to predict THz emission and guide design improvements.

Findings

The model accurately reproduces experimental THz emission characteristics.

Stacked layer configurations can significantly enhance THz generation efficiency.

Optimization strategies include using anti-reflection coatings and multilayer stacking.

Abstract

Metallic bilayer structures have been shown to emit strong terahertz (THz) pulses. We present a predictive multiscale mode that simulates optically induced spin-currents in a Fe/Pt-heterostructure and the emitted electric field. Electronic effects are treated on the nano scale using the Boltzmann transport equation for the dynamics of out-of-equilibrium charge carriers, numerically solved with a particle-in-cell code. The optical effects are simulated with a formalism that bridges the nanometer scale of the structure to the micrometer scale of the emitted waves. The approach helps to understand recent experimental findings on the basis of microscopic scattering effects and transport phenomena. Our theory's versatility allows it to be readily adapted to a wide spectrum of spintronic THz emitter designs. As an example, we show how the THz generation efficiency, defined as output to input power ratio, can be improved and optimized using serially stacked layers in conjunction with THz anti-reflection coatings.

Equations

$$A_m(z)=-\frac{d{P}_m(z)}{dz},$$

$$P_m(z)=\mathrm{Re}[{\mathcal{E}}_m(z){\mathcal{H}}_m^{\ast }(z)]$$

$$A=\sum _m{\int }_0^{d_m}{A}_m(z)dz,$$

$$\begin{split}\bigg{[}\frac{\partial}{\partial t}+v_{\sigma}(E)\cos(\theta)\frac{\partial}{\partial z}\bigg{]}&g_{\sigma}(z,E,\theta,t)=S_{\sigma}(z,E,t)-\frac{g_{\sigma}(z,E,\theta,t)}{\tau^{\mathrm{eff}}_{\sigma}(z,E)}\\ +&\sum_{\sigma^{\prime}}\int\frac{d\Omega^{\prime}}{4\pi}\int dE^{\prime}\,w(z,\sigma^{\prime},E^{\prime};\sigma,E)g_{\sigma^{\prime}}(z,E^{\prime},\theta^{\prime},t)\rho_{\sigma^{\prime}}(z,E^{\prime})\,,\end{split}$$

$${\tau }_{\sigma }^{\mathrm{eff}}(E)=1/({\tau }_{\mathrm{el}}^{-1}+{\tau }_{\sigma }^{-1}(E)),$$

$$\begin{array}{rl}{j}_s(z,t)=\int \frac{d\Omega }{4\pi }& \int dEv(E)\cos (\theta )\\ \times & \left[g_{\uparrow }(z,E,\theta ,t)-g_{\downarrow }(z,E,\theta ,t)\right].\end{array}$$

$${\mathbf{j}}_c=\gamma {\mathbf{j}}_s\times \frac{\mathbf{M}}{|\mathbf{M}|},$$

$$\left[\frac{d^2}{d{z}^2}+\frac{{\omega }^2}{c_0^2}{\tilde{n}}_{\mathrm{THz}}(z{)}^2\right]\mathcal{E}(z,\omega )=-i{\mu }_0\omega j_c(z,\omega ),$$

$$\mathcal{E}(z_{-},\omega )=\mathcal{E}(z_{+},\omega )=-\frac{1}{2}{\mu }_0{c}_0{\mathcal{J}}_0(\omega ),$$

Full text

I Introduction

Nanoscale spintronic structures excited by ultrashort optical pulses provide a novel emitter technology for broadband terahertz (THz) radiation Kampfrath et al. (2013); Seifert et al. (2016). Emission strength and efficiency has been intensively studied in structures of various materials and the generated pulses exceed the power of conventional THz emitters. As they are comparatively cheap to fabricate Seifert et al. (2016), they open up new possibilities to study intricate material parameters Cramer et al. (2018). In metallic heterostructures, the emission can be traced back to a charge current perpendicular to the layer stacking in the non-magnetic part of the slab Seifert et al. (2016). This charge current is generated by an optically induced spin-current through the inverse spin-Hall effect. The structure then acts as an antenna for these oscillating charge-currents and emits a short pulse with broad spectrum in the THz range. Subsequent studies found a strong dependence of the emitted pulse strength on the spin-Hall angle of the non-magnetic material Seifert et al. (2016). Both the thickness dependence and the laser excitation wavelength in Fe/Pt structures have been studied as well Torosyan et al. (2018); Papaioannou et al. (2018); Herapath et al. (2018).

Currently, there is no complete model to describe all the underlying processes from excitation to emission. Experimental studies concentrate on extracting the correct shape of the spin-current from the measured field pulse to compare with simulations Seifert et al. (2016). Other studies have used simplified approaches to extract the spin-diffusion length that describes the decay of the optically induced spin-currents Torosyan et al. (2018).

The purpose of this paper is to present a model that simulates the carrier dynamics after ultrashort laser excitation using the Boltzmann transport equation (BTE) in Fe/Pt multilayers Nenno et al. (2016). The BTE has been proven to be an adequate tool to simulate excited carrier dynamics in metallic structures on the nanoscale Hurst et al. (2018); Manfredi and Hervieux (2005); Nenno et al. (2016, 2018). To solve for the carrier-mediated spin-current, we use the particle-in-cell approach that introduces numerical superparticles to represent the hot-carrier distribution. Calculations of field propagation are both used as an input to estimate the excitation strength as well as to calculate spectrum and shape of the emitted pulse. Thus, we account for both laser absorption within the layers and the subsequent dynamics as well as for the propagation of the emitted fields through the structure.

The paper is organized as follows. We describe our approach to the laser absorption in the individual metal layers. We then introduce the guiding equation for the carrier dynamics, sketch its numerical solution and introduce the inhomogeneous Helmholtz equation for the emitted optical fields. We then present results for the induced spin-currents in a Fe/Pt heterostructure and the structure of the emitted field on the example of varying laser pulse duration. We show that the emitted field shape is directly proportional to the shape of the current. In the last part of the paper, we propose a serial stacking of spintronic emitters involving anti-reflective coatings to increase mono-energetic output in a feasible way.

II Model

We organize the different theoretical approaches presented in this section by the temporal order of the physical processes involved in the terahertz generation. All three parts, excitation, electron dynamics and emission as well as the structure are sketched in Fig. 1. First, an optical laser pulse irradiates the structure.

Hence, in a first step we calculate the number of excited carriers over the depth of the metallic bilayer. These excited carriers then traverse the structure, scatter and induce a spin current into the non-magnetic metal. The dynamics are described in a second step using a Boltzmann transport approach. Due to the inverse spin-Hall angle, the spin current induced in the platinum layer is converted into a perpendicular charge current. In the final step, we will use the calculated charge current as an input to solve the inhomogeneous Helmholtz equation, which we solve on the millimeter scale to obtain the emitted field and its spectrum.

II.1 Laser Absorption

The absorbed laser power is calculated for an excitation at 800 nm wavelength, which corresponds to laser wavelengths typically used in experiments Herapath et al. (2018); Papaioannou et al. (2018); Torosyan et al. (2018).

We use a standard finite-difference approach to solve the Helmholtz equation for the complex field amplitudes $\mathcal{E}$ and $\mathcal{H}$, denoting electric and magnetic field respectively. The boundary condition is an incoming wave at $\lambda_{0}=$ 800 nm. From the resulting fields, we calculate the negative derivative of the absorbed power in one layer $m$ Prentice (2000),

[TABLE]

where the amplitude of the Poynting vector is given by

[TABLE]

and is normalized to the amplitude of the incident field. The field amplitudes for $\mathcal{E}$ and $\mathcal{H}$ are given by $\mathcal{E}^{+}_{m}(z)+\mathcal{E}^{-}_{m}(z)$ and $\tilde{n}_{m}(\mathcal{E}^{+}_{m}(z)-\mathcal{E}^{-}_{m}(z))$ respectively. The right- ($\mathcal{E}^{+}$) and left-propagating ($\mathcal{E}^{-}_{m}$) contributions to the field can be easily extracted from the calculation by comparing field value and spatial derivative. The complex index of refraction $\tilde{n}_{m}$ at the laser wavelength of 800 nm is taken from Ref. Werner et al., 2009 for iron and platinum. The total absorbance of the structure is given by the sum over the contributions from individual layers Prentice (2000),

[TABLE]

where the integral runs over the thickness $d_{m}$ of each layer. For structures that are irradiated through a material with a refractive index that lies between that of Fe/Pt and vacuum, e.g., the frequently used substrate MgO, the obtained values for the total absorbed energy can be higher Torosyan et al. (2018).

Integrating over the total length of a Fe(10 nm)Pt(5 nm) slab, 41% of the laser pulse energy is absorbed when irradiated from the iron interface and 37% from the platinum interface. The depth-dependent absorption relative to the incident laser power is shown in Fig. 2(a).

Through the absorbed energy over the duration of the pulse, we can estimate the number of excited electrons as a function of position, assuming a linear relation for the generation rate Koopmans et al. (2000). The position-dependent hot-carrier number is then used as an input for the simulation of the electron dynamics.

II.2 Electron Dynamics

While the electromagnetic field propagates through the structure, electrons are excited to unoccupied states above the Fermi level. These hot carriers then move through the structure and relax back to equilibrium. This electronic dynamics is often called superdiffusive transport and we describe it with a Boltzmann transport equation Battiato et al. (2010); Nenno et al. (2016). On the femto- to picosecond timescale, we consider two types of scattering processes. Scattering events with impurities and phonons are considered elastic. Secondary-carrier generation due to hot-electron scattering with equilibrium carriers is included using inelastic carrier scattering times Nenno et al. (2016).

The formal treatment of electronic dynamics in the framework of the BTE starts with the distribution of hot carriers, $g_{\sigma}(z,E,\theta,t)$, which depends on spin $\sigma$, position $z$ in the slab, the particle energy $E$, its propagation angle with respect to the $z$-axis, $\theta$, and time $t$. We assume a universal polarization axis along which the spin aligns parallel or antiparallel to, such that $\sigma=\uparrow,\downarrow$. The evolution of this distribution function is given by Nenno et al. (2018)

[TABLE]

where we denote the carrier velocity by $v$ and the source term for excitation processes by $S$. The last two terms in Eq. (4) describe out- and in-scattering processes due to interactions with equilibrium carriers and many-particle excitations. The effective, spin-dependent lifetime is given by

[TABLE]

combining the elastic lifetime $\tau_{\mathrm{el}}$ and the inelastic scattering lifetime $\tau_{\sigma}(E)$. The integration over spin $\sigma$, energy $E$ and solid angle $\Omega$ in Eq. (4) distributes electrons in momentum space and includes secondary carrier generation. All scattering events are assumed to be local. We use material parameters from ab initio calculations presented in  Zhukov et al., 2006 and Ref. Kaltenborn and Schneider, 2014 and shown in Fig. 3. The numerical scheme to solve Eq. (4) is presented in Ref. Nenno et al., 2018. From the hot-electron distribution, we calculate the hot-carrier spin-current density

[TABLE]

The charge current density resulting from the inverse spin-Hall effect is calculated by

[TABLE]

following Ref. Seifert et al., 2016. Here, the spin-Hall angle is denoted by $\gamma$. In the following we assume a geometry shown in Fig. 1. The magnetization points in positive $x$-direction, so that Eq. (7) is simplified to $j_{c}=\gamma j_{s}$, with $j_{c}$ denoting the magnitude of the charge current in $y$-direction. The spin current $j_{s}$ flows in $z$-direction. The current that is driven by the charge to spin conversion is then used to calculate the emitted field.

II.3 Terahertz Emission

We determine the electric field $\mathcal{E}(z)$ emitted by the time-varying charge current in a metal layer using the inhomogeneous Helmholtz equation (see, e.g., p. 6 of Ref. Meystre and Sargent, 1991)

[TABLE]

where $\omega$ denotes the THz photon frequency, $c_{0}$ the vacuum speed of light, $\mu_{0}$ the vacuum permeability, $\tilde{n}_{\mathrm{THz}}(z)$ the complex refractive index of the materials at THz frequencies, and $j_{c}$ is given by Eq. (7). We stress that this is a completely different calculation from that for the absorption of the optical field, as typical thicknesses of spintronic emitters are on the order of nanometers, whereas the emitted THz waves have wavelengths in the micro to sub-millimeter range. Thus, on the THz wavelength scale, the charge current can be treated in the delta-layer approximation, $j_{c}(z,\omega)=\mathcal{J}_{0}(\omega)\delta(z-z_{0})$. The emitted field adjacent to the delta-layer then obeys the condition

[TABLE]

where $z_{\pm}=z_{0}\pm\eta$ ($\eta\downarrow 0$) is on the vacuum side adjacent to each vacuum-metal interface. This shows that the spectrum $I(\omega)\propto|\mathcal{E}(\omega)|$ and the temporal profile of the emitted field are directly proportional to the spectrum and profile of the current density, respectively. We note that the delta-layer approximation is a convenient tool to bridge nanoscale to millimeter-scale problems, allowing one to simulate more complex structures as used in experiments, in which dispersive elements (e.g. MgO substrate and detector components) modify the spectrum. In this approach, conventional optical transfer matrix techniques can easily be combined with the transfer matrix for the delta-layer source.

III Results

III.1 General Results

We consider first a Fe(10 nm)Pt(5 nm) structure surrounded by vacuum and irradiated by an infrared laser pulse. We choose first a pulse length of 10 fs, centered at 25 fs, and a photon energy of 1.5 eV Seifert et al. (2016). Subsequently, the pulse length will be varied to study its effect on the emitted spectrum. The laser excitation conditions are the same as in Ref. Nenno et al., 2018.

Fig. 2(b)-(c) show the calculated temporal profile of the spin current in platinum, as well as its spectrum. The spin current shows a characteristic bipolar shape, with a maximum close to the maximum of the laser pulse intensity at $25$ fs. We explain the bipolar shape as follows Nenno et al. (2018). Carriers with spin polarization are excited in iron and start move into the platinum slab, inducing a positive spin current. However, after a very short time ($\approx 10$ fs), carriers traverse the whole layer and change their propagation direction, effectively canceling any current. It is the small difference in propagating electron distribution undergoing scattering that generates an effective current. This current is dominated by majority carriers that enter the slab early on. When these majority carriers decay, the signal changes sign. The less intense tail of the spin current is then due to remaining minority carriers. The induced spin-current signal duration of approx. 150 fs for a laser pulse of 10 fs duration. For direct emission from the non-magnetic platinum layer, the spectrum is directly proportional to the charge current and therefore by Eq. (7) the spin current.

III.2 Laser Pulse Duration

In Fig. 4(a) shows emitted THz spectra for different temporal widths of the exciting laser pulse. We vary the FWHM of the pulse between 10 fs and 100 fs, covering the range around relevant scattering processes (e.g. $\tau_{\mathrm{el}}=30$ fs). We find two notable effects. The spectral maximum gets smaller and moves to lower frequencies for longer pulse duration. Fig. 4(b) shows both quantities plotted for all simulated excitation pulses. Longer excitation pulses induce hot-carrier distributions that persist over a longer period of time and thus, the resulting transport and relaxation dynamics result in a broader spin current with smaller magnitude (not shown). Hence, the spectrum of the emitted pulses peak at lower frequencies and with a smaller amplitude.

The second consequence of increasing pulse length is that the spectrum of the emitted THz radiation becomes narrower. In a simplified picture, the duration of the excitation would trigger dynamics on the same timescale and thus the spectral width of the emitted field should be proportional to the inverse of that duration. To analyze this in greater depth, we show the spectral width at 50% of its maximum for all laser pulse durations in Fig. 4(c). The results are compared to the estimate $1/\mathrm{FWHM}$. We find that for short excitation lengths, the simple estimate deviates drastically from the simulated results, whereas they agree much better for longer laser pulses.

We note that the simulation of longer pulses needs to include thermalized spin currents which are not accessible in our current theoretical framework. The dynamics of shorter pulses can, in principle, be described, but such pulses typically have intensities that trigger demagnetization dynamics outside the scope of a linear charge accumulation.

III.3 Multilayer

The energy of a single laser pulse is not completely absorbed within one emitter. For a thin emitter, e.g. Fe(3 nm)Pt(2 nm), approx. 40% is absorbed for normal incidence from either side of the structure. Based on this observation, we propose an optimized structure of serial spintronic emitters (SE) coated with an anti-reflective coating (ARC) optimized to improve emission of a fixed THz frequency. The goal is to improve the power conversion rate and propose a structure suitable for on-chip terahertz opto-spintronic devices Wu et al. (2018). Our model allows one to easily calculate the absorbed energy in these more complex structures, as well as the propagation and enhancement of the emitted radiation. We optimize the output for 4 THz, close the maximum of the spectrum for an incident laser pulse of 70 fs. The ARC is assumed to be non-absorbing with $n=3$ (polycrystalline diamond, for example, has $\tilde{n}_{THz}\approx 2.3+i3\times 10^{-4}$, see Ref. Dore et al., 1998 or $n\approx 3.0-3.4$ for crystalline Sapphire and Silicon, see e.g. Ref. Grischkowsky et al., 1990). Its thickness is fixed to $\lambda/4$. In general, an optimization with respect to the specific material is possible. In this study, we fix the configuration of the spintronic emitter, the optical parameters of the ARC and vary the spacing between ARC–SE–ARC segments and their number. Fig. 5 shows a sketch of the proposed heterostructure and the ratio between single emitter peak intensity and output of structures with a different number of spintronic emitter layers for varying distance between the emitters. We see that a ratio of almost 2 can be obtained with 3 to 4 segments spaced at $\Delta\approx\lambda/2$.

IV Conclusion

We presented a multiscale theoretical framework for the processes that contribute to THz emission from spintronic heterostructures excited by optical pulses. We presented a microscopic analysis of optically induced spin- and charge currents that determine a spectrum in the THz range, which is directly proportional to the emitted pulse spectrum in the thin-film limit. Intrinsic material timescales in the range of THz oscillations are key for the physical processes that convert power from the exciting laser field to the THz emission. Our theory is based on a flexible solution of the BTE and does not contain any fit parameters and only a few well-established transport parameters from ab initio calculations. The model can be extended to different structures with more challenging interface effects Jungfleisch et al. (2018); Dewhurst et al. (2018) or ferrimagnetic emitter layers Schneider et al. (2018); Seifert et al. (2017) by enhancing or replacing our method for calculating the electron dynamics in future studies. We showed that up to a number of four layers, a serial structure can result in an increased input- to output power ratio. The study can easily be extended further to perform an automated optimization with respect to the used material and laser parameters and sample configurations.

Acknowledgements.

We are grateful for discussions with Garik Torosyan, René Beigang, Evangelos Papaiouannou, Sascha Keller and Laura Scheuer. This work was supported by the German Science Foundation in the framework of the SFB/TRR 173 Spin+X (Project B03). D. M. N. thankfully acknowledges financial support from the Graduate School of Excellence MAINZ (Excellence Initiative DFG/GSC 266).