Terahertz electric-field-driven dynamical multiferroicity in SrTiO3

Sample details

The sample considered was a 10 mm × 10 mm, 500-µm-thick SrTiO₃ crystal substrate from MTI Corporation, with the [001] crystallographic direction normal to the cut direction. Both sides were polished. The sample was mounted tilted at 45° with respect to the free space propagation direction to measure the geometry through reflection. However, the large dielectric constant \((| \widetilde{\varepsilon }| \approx 100)\) of STO at the resonant frequency of the soft phonon mode causes such a large refraction of the pump terahertz beam that the propagation within the crystal is always orthogonal with respect to the sample surface.

Experimental methods

Broadband single-cycle terahertz radiation was generated by optical rectification in a DSTMS (4-N,N-dimethylamino-4′-N′-methyl-stilbazolium 2,4,6-trimethylbenzenesulfonate) crystal¹⁵ of a 40-fs-long, 800 μJ near-infrared laser pulse centred at a wavelength of 1,300 nm. This pulse was generated by optical parametric amplification from a 40-fs-long, 6.3 mJ pulse centred at 800 nm wavelength, which was produced by a 1 kHz regenerative amplifier. The terahertz pulses were focused onto the sample by three parabolic mirrors to a rounded beam of approximately 0.5 mm in diameter. Knowing the exact size of the beam is not crucial when estimating the fluence as we characterized the electric field of the radiation. A pair of wire-grid polarizers was used to tune the field amplitude without affecting the pulse shape. Narrowband terahertz radiation was obtained by filtering the broadband field with a 3 THz band-pass filter resulting in a peak frequency of 3 THz and a full-width at half-maximum of 0.5 THz. The probe beam was a 40-fs-long pulse at 800 nm wavelength, produced by the same 1 kHz regenerative amplifier used to generate the pump radiation. A β-barium borate crystal was used to convert the probe wavelength to 400 nm. The probe polarization was set by a nanoparticle linear film polarizer. The probe size at the sample was approximately 100 μm in diameter, substantially smaller than the terahertz pump. To record the change in the polarization state of the probe beam, a Wollaston prism was used to implement a balanced detection scheme with two photodiodes. A half-wave plate was used before the Wollaston prism to detect the Kerr rotation. The signals from the photodiodes were fed to a lock-in amplifier, whose reference frequency (500 Hz) came from a mechanical chopper mounted along the pump path.

Characterizing the terahertz electric field

The electric field component of the terahertz pulse at the sample location was characterized by electro-optical sampling⁵¹ in a 50-μm-thick (110)-cut GaP crystal. In particular, the field strength was calculated using a standard procedure, namely by measuring the time-resolved birefringence at a wavelength of 800 nm, which was caused by the terahertz electric field in the GaP crystal. To do that, we used a balanced detection scheme with two photodiodes measuring the I₁ and I₂ signals produced by a quarter-wave plate and a Wollaston prism placed after the sample. Once the time trace was retrieved, we moved the delay stage to the time delay where the maximum of the terahertz electric field was found. Following ref. ⁵², we computed \({E}_{{\rm{THz}}}={\lambda }_{0}\,{\sin }^{-1}(({I}_{1}-{I}_{2})/({I}_{1}+{I}_{2}))/(2\pi {n}_{0}^{3}{r}_{41}{t}_{{\rm{GaP}}}L)\), where n₀ = 3.193 is the refractive index of GaP at 800 nm, L = 50 μm is the thickness of the crystal, λ₀ = 800 nm is the probe wavelength, r₄₁ = 0.88 pm V⁻¹ is the GaP electro-optic coefficient⁵³ and t_GaP = 0.4769 is the Fresnel coefficient for reflective loss at the GaP crystal. For the DSTMS-generated broadband pulse, the maximum measured terahertz peak electric field was approximately 1.15 MV cm⁻¹, and the peak frequency of the terahertz pump pulse was at 2.7 THz with measurable components extending up to approximately 5 THz. After filtering the field with the 3 THz band-pass filter, a typical measured terahertz peak electric field was around 200–300 kV cm⁻¹. The sampled terahertz pump traces are reported in Extended Data Fig. 2, which shows the narrowband data representing the field used to measure all the data in the main text.

Circular polarization of the terahertz beam

After filtering the broadband terahertz radiation with a terahertz band-pass filter, the linear polarization can be converted into circular polarization with a terahertz quarter-wave plate oriented at ±45° to obtain the opposite helicities. For this purpose, we chose a Tydex quarter-wave plate made of x-cut terahertz-grade crystal quartz, whose thickness was adjusted to provide a π/2 phase shift at 3 THz. However, for wave plates, the phase shift is very sensitive to the radiation frequency. Moreover, we were not working with a monochromatic beam. Nonetheless, in our case the use of a quarter-wave plate was justified because the bandwidth of the filtered pulse was narrow enough to allow the wave plate to operate according to its design. To characterize the polarization state of the circular pump beam, we measured both the E_x and E_y electric field components with electro-optical sampling by rotating the GaP crystal by 90° around the light propagation axis. This allowed us to get the sensitivity to two orthogonal components of the terahertz pump, as shown in ref. ⁵¹. The time-domain traces obtained are shown in Extended Data Fig. 2 for the two helicities, labelled as LCP and RCP. From these traces, the polarization state can be unambiguously identified by calculating the Stokes parameters in the frequency domain, as shown in Extended Data Fig. 3. The S₃ parameter is associated with circular polarization, and a change of sign represents opposite helicities. The following inequality holds for a broadband pulse⁵⁴: \({({S}_{1}^{* }/{S}_{0}^{* })}^{2}+{({S}_{2}^{* }/{S}_{0}^{* })}^{2}+{({S}_{3}^{* }/{S}_{0}^{* })}^{2}\le 1\), where \({S}_{0}^{* }={\sum }_{i}{S}_{0,i}\), \({S}_{1}^{* }={\sum }_{i}{S}_{1,i}\), \({S}_{2}^{* }={\sum }_{i}{S}_{2,i}\), \({S}_{3}^{* }={\sum }_{i}{S}_{3,i}\) and i represents the ith frequency. The \({({S}_{3}^{* }/{S}_{0}^{* })}^{2}\) quantity gives an indication of the average amount of circular polarization in the terahertz pump pulse. Considering only the 0.5 THz full-width at half-maximum region of the peak, we found that the beam was 85–90% circularly polarized, as summarized in Extended Data Fig. 3.

Evaluating the complex refractive index

The complex refractive index \(\widetilde{n}=n+{\rm{i}}k\) of SrTiO₃ was derived from a combination of previous ellipsometry measurements on STO thin films⁵⁵ and hyper-Raman scattering in bulk STO¹⁸. In particular,

where \(\widetilde{\varepsilon }={\varepsilon }_{1}+{\rm{i}}{\varepsilon }_{2}\) is the complex permittivity. To estimate the permittivity in the experimental temperature range 160 K < T < 375 K for our specific sample, we first used the experimental data of ref. ⁵⁵ at T = 300 K (Extended Data Fig. 4), which contains a broadband response that can be fitted with a Lorentz oscillator. Then, to adjust it to our case, we rigidly shifted the curve, moving the peak from 3 to 2.7 THz, in accordance with our own data and ref. ¹⁷ on bulk samples. All centre frequency and linewidth values at the different temperatures are listed in Extended Data Table 1. All values for n and k are reported in Extended Data Fig. 4.

Modelling the terahertz reflectance, transmittance and absorptance

The electric field reflection, absorption and transmission properties were calculated for an air/STO/air stack using the analytical formulas for optical trilayers at normal incidence⁵⁶:

where \(\mathop{n}\limits^{ \sim }\) is the refractive index of STO, d is the thickness of the STO sample, λ is the wavelength and the refractive index of air is considered to be 1. The reflectance, transmittance and absorptance are given, respectively, by \(R={|r|}^{2}\), \(T={|t|}^{2}\) and A = 1 − R − T. Considering d = 500 μm, λ = 100 μm (3 THz) and \(\mathop{n}\limits^{ \sim }\) = 3.8 + i6.4 from ref. ¹⁷ at a temperature of 300 K, we get R ≈ 0.76, T ≈ 0 and A ≈ 0.24. As T ≈ 0, it is also interesting to estimate the decay length l_decay of the electric field inside STO, which indicates how much of the pump radiation penetrates into the sample:

The estimated penetration depth in the experimental temperature range 160 K < T < 375 K is listed in Extended Data Fig. 4.

Estimating the polarization rotation and magnetic field

Measuring the probe polarization rotation allowed us to calculate the magnetic field induced in STO. According to theory, the Faraday rotation ϑ_F and magnetic field are connected through the equation²⁰

as l_decay ≪ d, where d is the STO thickness. The parameter B represents the amplitude of the magnetic field at the surface, l_decay is the decay length of the pump field and V is the Verdet constant. The factor of 2 in the exponential function appears because the induced magnetic field is proportional to the square of the pump electric field. To extract the magnetic moment generating ϑ_F, we exploited the relation B = μ₀M, where M is the magnetization induced by the pump. Considering the STO lattice parameter a = 3.9 Å, the magnetic moment per unit cell μ is given by

Even if the measurements reported in the main text are performed in reflection (Kerr rotation ϑ_K), we evaluated the magnetic field considering a transmission measurement (Faraday geometry). Those were the only reliable values for the Verdet constant that we could find, which we were able to validate ourselves, as shown in Extended Data Fig. 5. To confirm that the reflection and transmission geometries give comparable responses, we measured the Faraday rotation during the same set of experiments described in the main text. We found that the absolute measured signal is within a factor of 2 compared to the Kerr rotation. Moreover, the pump penetration depth was still the limiting factor for the decay length l_decay to be considered in the above equation, as even in reflection, the STO thickness probed by the probe pulse is more than the pump penetration depth. The thickness contributing to the probe signal in reflection for an ultrafast pulse can be estimated through the distance travelled in the material during the pulse duration. For a 400 nm probe pulse (n_STO = 2.6 from ref. ⁵⁷) of 50 fs duration, the distance travelled in the STO during that interval corresponds to approximately 5.8 µm, which is longer than all the pump field penetration depths listed in Extended Data Fig. 4. According to ref. ⁵⁸, we have V ≈ 250 rad m⁻¹ T⁻¹, so that for ϑ_K = 10 μrad and l_decay = 2.49 μm, the magnetic field at the surface B ≈ 0.032 T. The average energy ϵ stored per unit surface in such a magnetic field is given by:

where the first factor of 1/2 in the definition of ϵ is due to the time average of the square of a sine wave, as we approximate the slow oscillation in Fig. 2a with a sinusoidal function. The integral takes into account that the induced magnetic field does not fill the whole sample volume but has a finite penetration depth. The energy ϵ is delivered by the pump pulse, and its fluence can be calculated by integrating the square of the trace shown in Extended Data Fig. 2 to give approximately 60 μJ cm⁻², which is much higher than the energy per unit surface in the generated magnetic field.

The pump fluence can be used to compute the absorbed energy density and give an estimate of the related temperature variation. As stated above, at 300 K, the decay length of the terahertz pump electric field is 2.49 μm and the absorptance is 0.24, which leads to an estimate of the average energy density absorbed by the sample of 115.7 mJ cm⁻³ = 0.043 meV per unit cell. The temperature increase for such an energy density can be obtained from the heat capacity and density of STO. Considering a density of 5.18 g cm⁻³ (MTI Corporation), a heat capacity of 100 J K⁻¹mol⁻¹ (ref. ⁵⁹) and a molar mass of 183.5 g mol⁻¹, the temperature increase is expected to be approximately 0.04 K, which could be neglected during the temperature-dependent measurements presented in the main text.

We also checked that the measured magneto-optical effect was not affected by the probe wavelength being too close to the bandgap of the material. In Extended Data Fig. 6, we present the total measured Faraday effect at both 400 and 800 nm probe wavelengths. Apart from an overall scaling factor consistent with the different values of the Verdet constant, the scaled response is identical for the two wavelengths, excluding wavelength-dependent artefacts.

Modelling the total Kerr effect

In ref. ²³ it was shown that the EKE response is given by

where E_x and E_y are the components of the pump pulse along generic x and y orthogonal directions, ϑ is the angle that x and y form with respect to the main crystallographic axes (i and j) and \({\Delta \chi \equiv \chi }_{{iiii}}^{(3)}-3{\chi }_{{iijj}}^{(3)}\), as \({\chi }_{{iijj}}^{(3)}=0.47{\chi }_{{iiii}}^{(3)}\) are the only two independent tensor components of the χ⁽³⁾ tensor in cubic STO from ref. ⁵⁸. If ϑ = 45°, then \(\Delta {\varGamma }^{{\rm{e}}}\propto 2{E}_{x}{E}_{y}\left[{\chi }_{iijj}^{(3)}+\frac{1}{2}\Delta \chi \right]\) and the signal is proportional to the product of the terahertz pump field components along perpendicular directions. For circularly polarized light of opposite helicities (left and right), the signal difference \(\Delta {\varGamma }^{{\rm{e}}}\left(\,{\rm{LCP}}\right)-\Delta {\varGamma }^{{\rm{e}}}\,({\rm{RCP}})\) is still proportional to E_xE_y, as only one of the two pump components changes sign.

Besides the EKE, it has been shown that an additional contribution, IKE, associated with the nonlinear excitation of the infrared-active soft phonon mode, is present²³. The IKE response \(\Delta {\varGamma }^{{\rm{ph}}}\) can be effectively modelled by replacing the E_x and E_y components in ∆Γ^e with a convolution between the pump and the single- or two-phonon propagators to account for the intermediate second-order excitation of the soft mode. Moreover, the χ⁽³⁾ tensor should be replaced with an effective nonlinear coupling between the pump and probe pulses and the infrared-active phonon.

An ab initio estimation of the effective nonlinear coupling is needed to estimate the relative weights of the IKE and the EKE in a rigorous way. This would require a state-of-the-art extension of the available density functional theory (DFT) codes, which has been investigated only recently⁶⁰ and goes far beyond the scope of this work. For this reason, in the main text we decided to model the full Kerr response with only the electronic contribution by assuming that the ionic contribution has a similar spectral content, so as not to introduce any free adjustable parameter into our simulations. For completeness, the full Kerr effect, including both the EKE and the IKE, is shown in Fig. 4a,b. The relative weight between the electronic and ionic contributions has been fixed to better reproduce the experimental time traces. This was done using the Kerr response measured with linearly polarized terahertz pulses as a reference.

Extended Data Fig. 7 compares the experimental and calculated responses of the material to linearly and circularly polarized terahertz pump fields. These measurements allowed us to isolate the dependence of the response on the polarization that is beyond the one captured by the third-order susceptibility, in particular for the EKE description. In Extended Data Fig. 7a,b, we use the linearly polarized pump data to match the experimental and calculated amplitudes. With the same scaling factor applied to all the data, in Extended Data Fig. 7c, we plot the difference between the experimental and simulated Kerr rotation for the circularly polarized pump case. The part of the signal before the zero-crossing point at approximately 2.4 ps can be explained in terms of the IKE effect (which we can also model as discussed above but was left out for simplicity of reasoning), whereas the negative dip after the zero-crossing point is the signature of dynamical multiferroicity.

Ab initio calculations

First-principles phonon calculations of cubic SrTiO₃ were performed within DFT using the Vienna Ab initio Simulation Package⁶¹, which implements the projector augmented-wave method⁶². The adopted projector augmented-wave potentials treat Sr 4s² 4p⁶ 5s², Ti 3s² 3p⁶ 4s², and O 2s² 2p² as the valence states. An energy cutoff of 550 eV was used, and the Brillouin-zone integration was performed with a 12 × 12 × 12 gamma-centred k-point mesh. The Heyd–Scuseria–Ernzerhof hybrid functional (HSE06; ref. ⁶³) was adopted to give an accurate description of the phonon potential energy surface. The lattice constant was optimized within HSE06; the optimized value of 3.900 Å agrees well with the experimental value, 3.905 Å.

The effective phonon frequencies and eigenvectors at room temperature were calculated based on SCP theory, as implemented in the software ALAMODE (ref. ⁶⁴). The harmonic and fourth-order interatomic force constants (IFCs), which are necessary as inputs to the SCP calculation, were calculated with the real-space supercell approach using a 2 × 2 × 2 supercell. The harmonic IFCs were estimated by systematically displacing each atom in the supercell from its equilibrium site by 0.01 Å, calculating forces by DFT and fitting the harmonic potential to the displacement–force datasets. The fourth-order IFCs were estimated using the compressive sensing method, for which 40 training structures were generated by combining DFT and molecular dynamics with random displacements following ref. ⁶⁵.

After obtaining the harmonic and anharmonic IFCs, the effective phonon frequency ω with branch index ν at wavevector q and the corresponding eigenvector were obtained by solving the SCP equation:

where \({\varLambda }_{{\bf{q}}}^{\left({\rm{HA}}\right)}={\rm{diag}}\left({\widetilde{\omega }}_{{\bf{q}},1}^{2},\ldots ,{\widetilde{\omega }}_{{\bf{q}},\nu }^{2}\right)\) with harmonic frequencies \(\widetilde{{\omega }}\). C is a unitary transformation matrix that modifies the polarization vector at finite temperature, Φ^SCP is the fourth-order anharmonic force constant and \({n}_{q,\nu }{(\omega }_{q,\nu })\) is the Bose–Einstein distribution. The equation was solved numerically by iteratively updating the effective frequency \({\omega }_{{\bf{q}}\nu }(T)\) and the unitary matrix C_q for the phonon modes at the gamma-centred 2 × 2 × 2 q points.

The summation over the \({{\bf{q}}}^{{\prime} }\) points was conducted with the denser 10 × 10 × 10 \({{\bf{q}}}^{{\prime} }\) points, which was sufficient to achieve convergence. The quartic coupling coefficient \({\varPhi }^{{\rm{SCP}}}(-{\bf{q}}\nu \,;{\bf{q}}\nu \,;-{{\bf{q}}}^{{\prime} }{\nu }^{{\prime} }\,;{{\bf{q}}}^{{\prime} }{\nu }^{{\prime} })\) was obtained from the fourth-order IFCs with the Fourier interpolation. The splitting of longitudinal optical and transverse optical modes was considered in the SCP calculation. The obtained SCP frequencies agree well with the inelastic neutron scattering data, as shown in Extended Data Fig. 8.

The anharmonic coupling coefficients of the triply-degenerate Γ₁₅ modes at room temperature were obtained by transforming the anharmonic IFCs into the normal coordinate basis. In this study, the anharmonic coupling terms up to the fourth-order were included in V(Q₁, Q₂). The normal coordinate Q_ν at finite temperature is given as \({Q}_{\nu }={\sum }_{\kappa }\sqrt{{m}_{\kappa }}{{\bf{e}}}_{\nu }\left(\kappa \right)\cdot {\bf{u}}(\kappa )\), where m_κ is the mass of atom κ and \({u}^{\alpha }(\kappa )\) is its displacement in the α direction. The polarization vector at room temperature e_ν(κ) was calculated as \({{\bf{e}}}_{\nu }\left(\kappa \right)={\sum }_{{\nu }^{{\prime} }}{\widetilde{{\bf{e}}}}_{{\nu }^{{\prime} }}(\kappa ){[{C}_{{\bf{q}}=0}]}_{{\nu }^{{\prime} }\nu }\), where \({\widetilde{{\bf{e}}}}_{\nu }\left(\kappa \right)\) is the harmonic polarization vector and C_q is the unitary matrix obtained as a solution to the SCP equation. As polarization mixing is significant in STO, the temperature dependence of the polarization vectors is noteworthy, as shown in Extended Data Table 2 for 300 K. As each atomic site of cubic STO is an inversion centre, all cubic coefficients became exactly zero. The effective charges of the Γ₁₅ modes were calculated as

with \({Z}_{\kappa ,\alpha \beta }^{* }\) being the Born effective charge of atom κ. For 300 K, we have that \({Z}_{{\rm{Si}},\alpha \beta }^{* }=2.553{\delta }_{\alpha \beta }\), \({Z}_{{\rm{Ti}},\alpha \beta }^{* }=6.704{\delta }_{\alpha \beta }\), \({Z}_{{\rm{O,\perp }}}^{* }=-\,1.941\) and \({Z}_{{\rm{O,\parallel }}}^{* }=-\,5.375\), in units of the electron charge. The effective charge of the oxygen atom is different when considering the direction perpendicular (⊥) or parallel (∥) to the nearest titanium atom, and δ_αβ is the Kronecker delta.

Phenomenological model for anharmonically coupled oscillators

To model the driven circular excitation of the ferroelectric mode along the two orthogonal directions, we derived the effective phonon potential for two of the threefold degenerate modes at q = 0, where the anharmonic coupling is included up to fourth order:

where Q is the normal coordinate in real space, the indices 1 and 2 refer to the two degenerate branches of the soft phonon along [100] and [010], k is the anharmonic contribution to the potential, and χ and ψ are the phonon–phonon coupling terms. As Q₁ and Q₂ are orthogonal to each other and the phonon potential spanned by them has a C₄ symmetry, ψ = 0. The resulting potential, with calculated parameters stated in Extended Data Table 1, represents two coupled anharmonic oscillators. The solution of this model is obtained by numerical integration of its equation of motion:

where Γ accounts for the lifetime of phonons and \({Z}^{\ast }{\mathop{E}\limits^{ \sim }}_{i}^{{\rm{T}}{\rm{H}}{\rm{z}}}\) is the oscillator coupling to the driving field through the mode effective charge \({Z}^{* }\). The effective field in the sample is expressed through the term \({\widetilde{E}}_{i}^{{\rm{THz}}}=\alpha {E}_{i}^{{\rm{THz}}}\) where α quantifies the amount of field actually experienced (not screened) by the sample. The value of \({E}_{i}^{{\rm{THz}}}\) was fixed from our experiment, whereas the values of Γ and ω at room temperature were taken from hyper-Raman measurements on bulk STO (ref. ¹⁷). Finally, the induced magnetic moment can be calculated via:

where i now represents the ith atom in the unit cell (Sr, Ti, O, O, O), \({\gamma }_{i}=e{Z}_{i}^{* }/2{m}_{i}\) is the gyromagnetic ratio and \({L}_{i}={Q}_{i}\times {\dot{Q}}_{i}\) is the angular momentum. The calculated magnetic moment per unit cell μ is shown in Fig. 4a in the time domain and in Fig. 4b in the frequency domain using an approach identical to that used to process the experimental data. All parameters used to solve the equation of motion were fixed, except for α, which was set to 0.7. The mode effective charge \({Z}^{* }\) and potential parameters k, χ and ψ were calculated from first principles (Extended Data Table 1), whereas the excitation field \({\widetilde{E}}_{i}^{{\rm{THz}}}\) and the phonon frequency ω and lifetime Γ are those obtained in experiments¹⁷.