Electronic properties and materials – Page 2

Observation of 2D-magnesium-intercalated gallium nitride superlattices

[ad_1]

Thermal annealing of metallic Mg on GaN

This intercalation occurs at low cost and with accessible conditions, specifically the annealing of a bare metallic Mg film on GaN in either a nitrogen or argon atmosphere at atmospheric pressure. This process is characterized by a rapid diffusion that requires only a few minutes at temperatures in the range 550–800 °C, which is fully compatible with existing semiconductor fabrication processes.

We first deposited metallic Mg thin film, without a capping layer, onto single-crystalline wurtzite GaN (hexagonal P6₃mc space group) using conventional physical vapour-deposition methods, such as electron-beam evaporation or sputtering. Subsequently, these samples were then annealed at elevated temperatures in atmospheric pressure, typically exceeding 550 °C for 10 min. Notably, higher annealing temperatures corresponded to reduced times required for the intercalation process. When heated above the melting point of Mg (650 °C), molten Mg also proved effective in producing this effect. The rapid interstitial diffusion of Mg within GaN resulted in its segregation into single monolayer atomic sheets. These sheets then expanded in size, resembling climbing motion of edge dislocations, and aligned vertically (along the c axis) with each other in an even and orderly manner. However, given the spontaneous and diffusion-driven nature of this reaction, the nucleation of Mg intercalants within the GaN matrix tended to be inhomogeneous, leading to a non-uniform distribution of the resulting MiGs domains (Extended Data Fig. 1).

Atmospheric pressure is strongly recommended to prevent the early sublimation of Mg before it reacts with GaN. Because Mg evaporates easily at low pressure, for example in the high-vacuum condition in a typical molecular-beam epitaxy (MBE) growth environment (10⁻⁸ to 10⁻¹² torr, which is roughly 10⁻⁶ to 10⁻¹⁰ Pa), the boiling point of Mg is below 200 °C⁵¹, whereas under 1 atmosphere, the boiling point of Mg increases to more than 1,000 °C, which is much higher than the melting point of 650 °C⁵². Therefore, if a 50-nm-thick Mg film is annealed in a rapid thermal annealing furnace at 700 °C under a vacuum of 10 Pa, given that the boiling point of Mg is approximately 500 °C, the Mg should have sublimated during the temperature ramp-up phase.

Before Mg deposition, if the GaN samples contain a top p-type GaN epitaxial layer that was initially doped with Mg during epitaxy, the samples were first subjected to a conventional annealing process (750 °C for 10 min) for Mg activation and hydrogen depassivation⁵³.

We then deposited an amorphous, pure, bare metallic Mg thin film, approximately 50 nm thick, onto the GaN samples using an ACS-4000 (ULVAC) sputtering system with a 3 N (99.9% purity) Mg target.

After deposition, the handling of metallic Mg film requires special attention to prevent oxidation and dissolution during the lift-off process in liquids. Generally, the oxidation resistance of Mg is considered good in dry air up to approximately 400 °C, and up to about 350 °C in moist air, with humidity having an important role in Mg corrosion as it forms Mg(OH)₂ (ref. ⁵⁴). Based on our experience, at room temperature the oxidation resistance of Mg is robust; the as-deposited Mg thin film (50 nm) retains its metallic lustre when kept in air for days. We did not take special measures to prevent its oxidation, and the thermal annealing of the Mg film was typically conducted within one day after deposition. However, Mg is prone to oxidation in dry air above 500 °C, so pure nitrogen or argon is needed during annealing. Nonetheless, we acknowledge that the Mg film surface was probably oxidized and covered by a thin oxide layer as a result of its exposure to air. This might also explain the extra oxygen that was incorporated into the MiGs structure in the secondary ion mass spectrometry (SIMS) profiles (Supplementary Fig. 3).

For photolithography and patterning of Mg, the lift-off process of photoresist and Mg film was done in liquid solution. We found that Mg nanofilm tended to dissolve in deionized water, losing its metallic colour after just a few minutes, so we avoided using deionized water. Instead, we opted for pure isopropyl alcohol (IPA), finding that Mg remains stable in IPA, which is sufficient for the lift-off process. IPA is routinely used instead of deionized water in semiconductor processing practices. We also used N-methyl-2-pyrrolidone (or acetone) followed by methanol for the lift-off process before using IPA.

The GaN samples, after undergoing blank deposition or pattering of Mg film, were annealed at various temperatures: 500 °C, 550 °C, 600 °C, 700 °C and 800 °C, all for 10 min. Unless specified otherwise in comparative experiments for annealing temperature, we set the default annealing temperature at 800 °C. We found that the formation of superlattice structure was characterized by a rapid diffusion process that required only 10 min or less at temperatures ranging from 550 °C to 800 °C. Temperatures between 800 °C and 1,000 °C are feasible but made no substantial difference compared with temperatures of around 800 °C. However, annealing at temperatures higher than 1,000 °C requires a masking layer to prevent incongruent decomposition of GaN. Furthermore, annealing at temperatures of more than 1,000 °C favours substitutional diffusion of Mg⁵⁵, which could lead to thermal instability of the superlattice structure and decomposition of the interstitial Mg sheets. Temperatures lower than 550 °C might be possible but could require much longer annealing times⁴⁵.

After cooling, the samples were cleaned with either hydrochloric acid or aqua regia. This cleaning step was essential for removing any Mg residues and compounds that were soluble in acid, which might have formed during exposure to moisture, hot air or hot nitrogen.

Preparation of GaN samples

We used a diverse range of GaN samples to demonstrate the accessibility and versatility of the Mg intercalation process. Our findings indicate that, regarding the formation and morphology of the superlattice structure, there is no great difference between n-type GaN and p-type GaN, except for the initial Mg content in these samples. Specifically, n-type GaN does not contain Mg initially, whereas p-type GaN is already enriched with Mg. These pre-existing Mg atoms could contribute to the formation of Mg sheets. Consequently, the presence of Mg in p-type GaN may slightly lower the diffusion temperature or the time required to form the MiGs structure. Thus, including studies on both n-type and p-type GaN samples in our research may highlight the general applicability of this method, regardless of their conductivity type and film quality. As a result, our GaN samples varied in initial doping concentration (heavily n-doped, which is a commercial n-type GaN substrate without re-epitaxy, and lightly n-doped, which is unintentionally doped GaN, lightly p-doped and heavily p-doped) and substrate type (commercial GaN templates, consisting of GaN buffer layers on sapphire substrates, and commercial GaN substrates).

1.

For samples used for the HAADF-STEM and iDPC-STEM observations, except for the STEM observation with I–V characteristics in Fig. 4, the unintentionally doped GaN 2 µm thick was grown by metal–organic vapour phase epitaxy on a commercial (0001) heavily doped n-type GaN substrate 400 μm thick and 2 inches (51 mm) in diameter with a free electron concentration of around 10¹⁸ cm⁻³.
2.

For samples used for the STEM-cathodoluminescence (CL) characterizations, these GaN substrates with epitaxy-ready smoothness were used.
3.

For GaN samples used for plan-view CL and cross-sectional polarized CL measurements, we used epitaxy layers consisting of 200-nm-thick p-type GaN with [Mg] = 6 × 10¹⁸ cm⁻³ followed by 1-μm-thick unintentionally doped GaN on the above-mentioned GaN substrate.
4.

For GaN samples used for terahertz time-domain ellipsometry, the epitaxy layers consisted of p-type GaN around 450 nm thick, one with [Mg] = 6 × 10¹⁸ cm⁻³ and another with [Mg] = 2 × 10¹⁹ cm⁻³, followed by 3.9-μm-thick unintentionally doped GaN ([C] = 2 × 10¹⁶ cm⁻³, [O] = 1 × 10¹⁶ cm⁻³ and [Si] below the detection limit) on commercial sapphire substrates 2 inches in diameter. We also prepared another GaN sample with the same 3.9-μm-thick unintentionally doped GaN (without p-type GaN) on the above-mentioned sapphire substrate and a bare sapphire substrate.
5.

For GaN samples used for C–V profiling, the epitaxy layers consisted of unintentionally doped GaN roughly 2 μm thick with a net doping concentration of 6 × 10¹⁶ cm⁻³ on the above-mentioned GaN substrate.
6.

For GaN samples used for backside secondary-ion mass spectrometry (SIMS), the epi-layers consisted of 400-nm-thick p-type GaN ([Mg] = 7 × 10¹⁸ cm⁻³) and 2.5-μm-thick unintentionally doped GaN with [Si] = 7 × 10¹⁶ cm⁻³ on the above-mentioned GaN substrate.
7.

For GaN samples used for I–V characteristics and the combined STEM observations, the samples were the same with those in item (6). The conventional top degenerately doped p-type contact layer (the p⁺ contact layer) intended for ohmic contact purposes was not grown intentionally.

STEM-based characterizations

STEM samples were prepared using a Thermo Fisher Helios 660 focused ion beam (FIB) STEM. Before FIB sample preparation, a protective carbon layer was deposited on the sample surface to reduce Ga FIB-induced damage. Samples were thinned with starting and final energies of 30 kV and 1 kV, respectively. The resulting lamellae were less than 50 nm thick. HAADF-STEM imaging was performed on a probe corrected Hitachi HD 2700 STEM operated at 200 kV with a convergence semi-angle of 23.5 mrad.

The iDPC-STEM imaging was performed on a probe corrected Thermo Fisher Themis 80-300 kV STEM operated at 200 kV with a convergence semi-angle of 25.3 mrad using a four-quadrant DF4 detector⁵⁶.

For the elemental EDS maps of Ga, N and Mg, the measurements were taken using a Thermo Fisher Super-X system in the same STEM set-up as that used for the iDPC data. The Super-X system consists of four silicon drift detectors for improved X-ray collection efficiency and offers good sensitivity and high spatial resolution capability. A radial Wiener filter was used for the elemental EDS maps of Ga, N and Mg.

STEM-CL characterization involved sample preparation by FIB using a Thermo Fisher Scientific Helios 660 with lamella specimens approximately 500 nm thick. The STEM-CL measurements were done on a JEOL JEM-2100 F analytical electron microscope at an acceleration voltage of 120 kV and a beam diameter of 1.0 nm. High-resolution STEM observations associated with the STEM-CL measurements were performed on a JEOL JEM-ARM200F analytical electron microscope at an acceleration voltage of 200 kV and a beam diameter of 0.15 nm.

Strain determination

For the in-plane and out-of-plane strains, they are defined by the following formulas, respectively:

$${\varepsilon }_{a}( \% )=\frac{{a}_{{\rm{s}}}-{a}_{{\rm{r}}}}{{a}_{{\rm{r}}}}\times 100$$

(2)

$${\varepsilon }_{c}( \% )=\frac{{c}_{{\rm{s}}}-{c}_{{\rm{r}}}}{{c}_{{\rm{r}}}}\times 100$$

(3)

where the subscript s refers to GaN in the strain-mapping region (GaN between periodic 2D-Mg_i sheets excluding the tip region) and the subscript r denotes GaN in the reference region (outside the intercalation region). The terms a and c represent the in-plane and out-of-plane lattice constants of hexagonal GaN, respectively. These constants are obtained on the basis of the measured average distance from atom column to atom column. The average in-plane distance, d_in-plane, and out-of-plane distance, d_out-of-plane, have the following relationship with the lattice constants in a hexagonal lattice:

$${d}_{{\rm{out}} \mbox{-} {\rm{of}} \mbox{-} {\rm{plane}}}=\frac{1}{2}c$$

(4)

$${d}_{{\rm{in}} \mbox{-} {\rm{plane}}}=\frac{\sqrt{3}}{2}a$$

(5)

The direction of the in-plane strain is parallel to the m axis (the $ < 1\bar{1}00 > $ direction) and also equivalently to the a axis (the $ < 11\bar{2}0 > $ direction) in this work, owing to the symmetry of strain produced by the Mg sheets parallel to both the m axis and the a axis. The direction of out-of-plane strain is parallel to the c axis (the <0001> direction). The directions of the m, a and c axes are indicated in the lower left corner of Fig. 3b. When GaN is uniaxially strained along the c axis, the Poisson ratio, ν, is defined as:

$$\nu =-\frac{{\varepsilon }_{\perp }}{{\varepsilon }_{\parallel }}=-\frac{{\varepsilon }_{a}}{{\varepsilon }_{c}}$$

(6)

For atom-column by atom-column distance measurements, images were acquired using a frame-averaging approach to limit the effects of sample drift. Each frame-averaged dataset consisted of a minimum of 12 individual images that were then position-corrected using cross-correlation and averaged for the final image. To further limit the effects of drift, frame-averaged image series were acquired with the fast scan direction parallel to the measurement direction. The atom-column by atom-column measurements were performed by fitting each atom column position to a two-dimensional Gaussian distribution and measuring the distance between columns directly using the Atomap software package⁵⁷.

Electronic band-structure calculation

The electronic band structures of GaN were calculated according to density functional theory (DFT) using the Vienna ab initio Simulation Package. The nuclei and core electrons were simulated by pseudopotentials generated by the projector-augmented-wave method. The exchange-correlation energy was treated in the generalized gradient approximation using the Perdew–Burke–Ernzerhof functional. In the electronic-structure calculations with the Heyd–Scuseria–Ernzerhof hybrid functional, the exact exchange amount was set to 26%, yielding a bandgap of 3.4 eV for unstrained wurtzite GaN, which aligns well with the experimental value for bulk GaN.

DFT simulations were not performed on the broader GaN and Mg superlattice for two main reasons. First, the MiGs structure, like the Guinier–Preston zone, is categorized as a structure rather than a compound material, the latter being exemplified by MAX phases. The non-stoichiometry of a structure greatly increases the computational expense of the relevant calculations. Merely constructing an Mg sheet, as depicted in Fig. 1e, would suggest that Mg sheets extend throughout the entire layer, which would result in substitutional rather than interstitial intercalation. Consequently, to honour the characteristic of interstitial intercalation, the Mg sheets must terminate within the GaN matrix. This necessitates a supercell comprising hundreds of atoms at least. Second, the gradual transition of polarity in the GaN lattice, as illustrated in Fig. 2f, substantially complicates the construction of a supercell model that properly reflects the unique features of such a structure. As such, DFT simulations on the GaN and Mg superlattice may require much more effort in future studies.

Scanning electron microscopy (SEM)-CL mapping and polarization-resolved CL spectroscopy

SEM-CL was performed on an Allalin spectroscopic (analytical) platform from Attolight. A focused electron beam scanned the sample with normal incidence and the optical emission was collected. Angles rotationally symmetric between 18° and 46° from the incident beam were captured. The optical beam was then analysed by using a spectrometer to obtain a luminescence spectrum for each pixel of the image. The optical analysis apparatus consisted of a monochromator fitted with a 1,024×256 charge-coupled device high-speed camera adapted for ultraviolet-visible detection (wavelengths of 200–1,100 nm). The camera enabled almost-instant acquisition of the entire emission spectrum. The samples were measured in both plan and cross-sectional views (for polarization-resolved CL only) at room temperature. The voltage of the incident electron beam was set to 10 kV. The scans were square and their size varied from 5 μm to 50 μm, depending on the width of the measured line on the sample, and the number of data points was kept at 128×128. The integration time at each pixel was set to 10 ms, resulting in a total acquisition time for one image of 165 s. Acquired luminescence spectra resulted from electronic transitions around the bandgap of GaN. After acquisition, the spectrum at each pixel was approximatively fitted with a Gaussian curve. Peak energy was determined from the Gaussian centre.

THz-TDE

Free carrier absorption exists in the terahertz region, which allows the characterization of electrical conductivity⁵⁸. This method is preferred to the conventional van der Pauw method because the latter only measures sheet resistivity (the inverse of in-plane conductivity, σ_⊥), whereas THz-TDE encompasses both in-plane (σ_⊥) and out-of-plane (σ_∥) components⁵⁹.

Our THz-TDE was performed using the Tera Evaluator, a commercial terahertz time-domain ellipsometer from Nippo Precision, with a reliable spectral range of 1.0–2.5 THz. The angle of incidence was fixed at 70° (ref. ⁵⁸). To characterize the change in polarization of the terahertz waves after reflection on the sample, the ellipsometric parameters, namely amplitude ratio (tan Ψ) and phase difference (Δ) between the p and s polarization components, were measured⁵⁹.

These ellipsometric parameters were initially measured from a bare sapphire substrate. This was followed by measurements on a 3.99-μm-thick unintentionally doped (lightly n-doped) GaN on the aforementioned sapphire substrate, and measurements of three 430-nm-thick p-type/3.99-μm-thick unintentionally doped (lightly n-doped) GaN/sapphire samples (one of them has incorporated MiGs nanostructures, formed by annealing of 50 nm metallic Mg film onto the entire GaN at 800 °C and 10 min). Three p-type GaN samples formed two control groups: the first consisted of p-type GaN samples with varying initial Mg doping concentrations (7 × 10¹⁸ cm⁻³ and 2 × 10¹⁹ cm⁻³); the second comprised the lower-Mg-doped p-type GaN samples, with and without MiGs phases, named samples 1 and 2, respectively. The higher-Mg-doped GaN sample without MiGs phases was known as sample 3 (Extended Data Fig. 6a).

During measurements, the equipment chamber was purged with dry air to minimize terahertz absorption by water vapour, and the ambient temperature was maintained at 22 °C. The terahertz spot size on the sample surface is 10 mm for the major diameter and 3 mm for the minor diameter. For the analysis, a three-layer optical model, approximating the sample structure, was constructed. This model was then used to calculate the predicted response on the basis of Fresnel equations, which describe the light reflection at the interface between each layer, taking into account thickness and optical constants. The calculated values were subsequently compared with the experimental data, and the solution was optimized using regression analysis.

In the MiGs-incorporated p-type GaN epilayer, the measured optical and transport properties represent a composite of both p-type GaN and the MiGs phase. Fitting the conductivity spectra to the Drude model reveals the inverse of scattering time at the crossing point of the σ′ and σ″ curves (Extended Data Fig. 6e). Sample 3 had a shorter scattering time, τ, and higher conductivity, σ, consistent with its higher ionized impurity concentration. However, the scattering time for sample 2 was similar to that of sample 1. The sheet resistivity and hole concentration for these samples, also measured by the van der Pauw method, are summarized in Extended Data Fig. 6i for comparison. We observed that the resistivity measured by THz-TDE, denoted as ρ_DC, of sample 1 and sample 3 closely aligns with their sheet resistivity measured using the van der Pauw method, validating the optical model used in our THz-TDE analysis. By contrast, sample 2 exhibits a substantial decrease in ρ_DC compared with its sheet resistivity. The six-fold conductivity increase suggests there are important factors beyond the increased hole concentration, minor changes in the dielectric constant and surface roughness. Given the consistent scattering time in samples 1 and 2, it implies that sample 2 has greatly enhanced out-of-plane mobility, potentially resulting from a reduced effective mass from the predominance of split-off holes along the c axis.

Based on the Drude model, in which conductivity (σ) is proportional to carrier concentration and scattering time (τ) and inversely proportional to effective mass, our preliminary estimates (using data from Extended Data Fig. 6i) indicate that the average effective mass of holes in sample 2 is reduced to approximately 30% of that in sample 1.

For more specific details, in THz-TDE, p– and s-polarized terahertz waves were irradiated onto a sample at an oblique angle, and the physical properties of the sample are derived from the change in the polarization state of terahertz waves on reflection. The measured values (tan Ψ and Δ) are defined from the ratio of the reflection coefficients for p and s polarizations, ${\widetilde{r}}_{p}$ and ${\widetilde{r}}_{s}$:

$$\mathop{\rho }\limits^{ \sim }=\,\tan \,\varPsi \exp i\varDelta =\frac{{\mathop{r}\limits^{ \sim }}_{p}}{{\mathop{r}\limits^{ \sim }}_{s}}=\left(\frac{{\mathop{E}\limits^{ \sim }}_{r,p}}{{\mathop{E}\limits^{ \sim }}_{i,p}}\right)/\left(\frac{{\mathop{E}\limits^{ \sim }}_{r,s}}{{\mathop{E}\limits^{ \sim }}_{i,s}}\right)$$

(7)

where ${\widetilde{E}}_{r}$ and ${\widetilde{E}}_{i}$ represent the reflected and incident electric fields, respectively. In the measurement, the incident field is oriented such that ${\widetilde{E}}_{i,p}={\widetilde{E}}_{i,s}$. Hence, equation (7) can be simplified to $\mathop{\rho }\limits^{ \sim }=\tan \,\varPsi \exp i\varDelta ={\mathop{E}\limits^{ \sim }}_{r,p}/{\mathop{E}\limits^{ \sim }}_{r,s}$, and it follows that

$$\tan \,\varPsi =\frac{|{\mathop{E}\limits^{ \sim }}_{r,p}|}{|{\mathop{E}\limits^{ \sim }}_{r,s}|}\hspace{2.22144pt}{\rm{a}}{\rm{n}}{\rm{d}}\,\varDelta ={\delta }_{r,p}-{\delta }_{r,s}$$

(8)

where δ represents phase. Therefore, tan Ψ and Δ represent the amplitude ratio and phase difference, respectively, between the p and s polarizations. These parameters are referred to as ellipsometric parameters.

The reflection coefficients for p and s polarizations are described by Fresnel equations:

$$\begin{array}{c}{\widetilde{r}}_{jk,p}=\frac{{\widetilde{n}}_{k}\cos {\theta }_{j}-{\widetilde{n}}_{j}\cos {\theta }_{k}}{{\widetilde{n}}_{k}\cos {\theta }_{j}+{\widetilde{n}}_{j}\cos {\theta }_{k}}\,\\ {\widetilde{r}}_{jk,s}=\frac{{\widetilde{n}}_{j}\cos {\theta }_{j}-{\widetilde{n}}_{k}\cos {\theta }_{k}}{{\widetilde{n}}_{j}\cos {\theta }_{j}+{\widetilde{n}}_{k}\cos {\theta }_{k}}\end{array}$$

(9)

where the subscripts j and k represent the media of light propagation, $\widetilde{n}$ is the complex refractive index expressed by $\widetilde{n}=n-i\kappa $, θ_j is the angle of incidence and θ_k is the transmission angle, which can be calculated using Snell’s law.

For a bulk sample that consists of only a single layer, the complex refractive index can be obtained directly from the measured ellipsometric parameters:

$${\widetilde{n}}_{1}={\widetilde{n}}_{0}\,\sin {\theta }_{0}\sqrt{{\left(\frac{1-\widetilde{\rho }}{1+\widetilde{\rho }}\right)}^{2}{\tan }^{2}{\theta }_{0}+1}$$

(10)

where ${\widetilde{n}}_{0}$ is the complex refractive index of the ambient (air) and θ₀ is the angle of incidence.

For a thin film/substrate structure, the complex refractive index of the film can be obtained numerically using the following expression:

$$\widetilde{\rho }=\frac{{\widetilde{r}}_{p}}{{\widetilde{r}}_{s}}=\left[\frac{{\widetilde{r}}_{01,p}+{\widetilde{r}}_{12,p}\exp (-i2\beta )}{1+{\widetilde{r}}_{01,p}{\widetilde{r}}_{12,p}\exp (-i2\beta )}\right]/\left[\frac{{\widetilde{r}}_{01,s}+{\widetilde{r}}_{12,s}\exp (-i2\beta )}{1+{\widetilde{r}}_{01,s}{\widetilde{r}}_{12,s}\exp (-i2\beta )}\right]$$

(11)

where β represents the film phase thickness given by $\beta =2{\rm{\pi }}d{\widetilde{n}}_{1}\cos {\theta }_{1}/\lambda $, λ is the wavelength and d is the film thickness.

For a three-layer sample structure, the ellipsometric parameter can be expressed as:

$$\begin{array}{c}\mathop{\rho }\limits^{ \sim }=\frac{{\mathop{r}\limits^{ \sim }}_{p}}{{\mathop{r}\limits^{ \sim }}_{s}}=\left[\frac{{\mathop{r}\limits^{ \sim }}_{01,p}+{\mathop{r}\limits^{ \sim }}_{123,p}\exp (-i2{\beta }_{1})}{1+{\mathop{r}\limits^{ \sim }}_{01,p}{\mathop{r}\limits^{ \sim }}_{123,p}\exp (-i2{\beta }_{1})}\right]\,/\,\left[\frac{{\mathop{r}\limits^{ \sim }}_{01,s}+{\mathop{r}\limits^{ \sim }}_{123,s}\exp (-i2{\beta }_{1})}{1+{\mathop{r}\limits^{ \sim }}_{01,s}{\mathop{r}\limits^{ \sim }}_{123,s}\exp (-i2{\beta }_{1})}\right]\\ \,=\,\tan \,\varPsi \exp i\varDelta \end{array}$$

(12)

where ${\widetilde{r}}_{01}$ represents the reflection coefficient for the ambient–first layer interface and ${\widetilde{r}}_{123}$ is the reflection coefficient for the second layer and substrate, and is in turn given by

$${\widetilde{r}}_{123}=\frac{{\widetilde{r}}_{12}+{\widetilde{r}}_{23}\exp (-i2{\beta }_{2})}{1+{\widetilde{r}}_{12}{\widetilde{r}}_{23}\exp (-i2{\beta }_{2})}$$

(13)

The measured ellipsometric parameters for the sapphire substrate and n-type GaN/sapphire sample are shown in Extended Data Fig. 6f. Using equations (10) and (11), the complex refractive index of sapphire and of the n-type GaN film were obtained, respectively, as shown in Extended Data Fig. 6g. The extinction coefficient of the n-type GaN film shows an increasing trend towards low frequencies, which can be attributed to free carrier absorption. Extended Data Fig. 6h shows the complex refractive-index spectra obtained for the p-type GaN layers in the p-GaN/n-GaN/sapphire samples using equation (12). The spectra are fitted to the Drude model. In these measurements, the optic axis of the samples is normal to the surface.

It is important to note that the optical constants were calculated under the assumption of an isotropic model, despite GaN being a uniaxial crystal. Thus, the solutions are only an approximation of the dielectric function and do not represent exact optical constants, although they can be more appropriately regarded as pseudo-dielectric functions. Regardless, they provide insights into the differences between the optical responses of samples with and without MiGs phases and offer semi-quantitative evidence of an enhancement in conductivity in Mg-intercalated GaN.

Van der Pauw method

After the THz-TDE measurements, Hall-effect measurements (the van der Pauw method) at room temperature were performed on the above-mentioned p-type GaN samples (ResiTest8400 system) using a magnetic field with an amplitude of 0.5 T. The samples were diced into small pieces and the van der Pauw pattern for each measurement consisted of 1.5 mm × 1.5 mm. Ni/Au electrodes were formed by electron-beam evaporation. The contact metal was annealed at 525 °C for 5 min in O₂ ambient to reduce the contact resistance.

SIMS

SIMS was performed at Eurofins EAG Laboratories. For backside SIMS measurement, samples were first polished and thinned to approximately 1 μm before being sputtered from the backside. Cs⁺ and O⁻ primary ion beams were used to detect negative secondary ions (from C, H, O and Si) and positive secondary ions (Mg), respectively.

Electrical characterizations

For C–V profiling, Ni (70 nm)/Au (70 nm) stacks were electron-beam evaporated to form circular Schottky-contact electrodes (around 220 µm in diameter) without sintering. A shadow mask was used for patterning, chosen over photolithography to prevent potential damage from developing agents to the sample surface. Ti (20 nm)/Au (100 nm) stacks were electron-beam evaporated on the backside of the n-type GaN substrates for ohmic contact without sintering. All C–V measurements were performed at a frequency of 100 kHz using an E4980A LCR meter by Agilent.

For measurement of I–V characteristics combined with STEM observations, the samples were prepared with half of the region blank deposited with Mg and the other half was patterned with Mg using maskless photolithography (Nanosystem Solutions DL-1000). After this, the samples underwent the previously mentioned annealing process and were subsequently cleaned with acid. In the next step, using maskless photolithography, Ni (20 nm)/Au (150 nm) was electron-beam evaporated onto regions approximately identical to those where Mg had previously been deposited. This was done with the aid of alignment photolithography, achieving a precision of around 2 μm, to ensure the absence of Mg in the regions between the Ni/Au electrodes. Finally, sintering was carried out at 525 °C for 5 min in an O₂ atmosphere to form the ohmic contact. After this, the samples were diced such that the half entirely covered with Mg-annealed surface was designated for STEM observations, whereas the other half, with selectively patterned Mg, was used for evaluating ohmic contact resistance through I–V characterization. An Agilent B1505A analyser was used to measure I–V characteristics in this work.

[ad_2]

Source link

Terahertz electric-field-driven dynamical multiferroicity in SrTiO3

[ad_1]

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,

$$n=\sqrt{\frac{\left|\widetilde{\varepsilon }\right|+{\varepsilon }_{1}}{2}},\qquad k=\sqrt{\frac{\left|\widetilde{\varepsilon }\right|-{\varepsilon }_{1}}{2}},$$

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⁵⁶:

$$\begin{array}{c}r=\frac{C-C\,\exp (2{\rm{i}}\delta )}{1-{C}^{2}\,\exp (2{\rm{i}}\delta )},\,t=\frac{(4\mathop{n}\limits^{ \sim }/{(1+\mathop{n}\limits^{ \sim })}^{2})\,\exp ({\rm{i}}\delta )}{1-{C}^{2}\,\exp (2{\rm{i}}\delta )},\\ C=\frac{1-\mathop{n}\limits^{ \sim }}{1+\mathop{n}\limits^{ \sim }},\,\delta =\frac{2\pi d}{\lambda }\mathop{n}\limits^{ \sim },\end{array}$$

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:

$${l}_{{\rm{d}}{\rm{e}}{\rm{c}}{\rm{a}}{\rm{y}}}(\mathop{n}\limits^{ \sim })=\frac{\lambda }{2\pi {\rm{\Im }}(\mathop{n}\limits^{ \sim })}\approx 2.49\,{\rm{\mu }}{\rm{m}}.$$

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²⁰

$${{\vartheta }}_{{\rm{F}}}=VB{\int }_{0}^{d}\exp \left(-2\frac{z}{{l}_{{\rm{decay}}}}\right)\,{\rm{d}}z=VB\frac{{l}_{{\rm{decay}}}}{2},$$

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

$$\mu =\frac{2{{\vartheta }}_{{\rm{F}}}{a}^{3}}{{\mu }_{0}V{l}_{{\rm{decay}}}}.$$

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:

$${\epsilon }=\frac{1}{2}\frac{1}{2{\mu }_{0}}{B}^{2}{\int }_{0}^{d}\exp \left(-4\frac{z}{{l}_{{\rm{decay}}}}\right)\,{\rm{d}}z=\frac{1}{16{\mu }_{0}}{B}^{2}{l}_{{\rm{decay}}}\approx 0.013\,{\rm{\mu }}{\rm{J}}\,{{\rm{cm}}}^{-2},$$

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

$$\Delta {\varGamma }^{{\rm{e}}}\propto \frac{1}{4}\left[{E}_{x}^{2}-{E}_{y}^{2}\right]\Delta \chi \sin \left(4{\vartheta }\right)+2{E}_{x}{E}_{y}\left[{\chi }_{iijj}^{(3)}+\frac{1}{2}\Delta {\chi \sin }^{2}(2{\vartheta })\right],$$

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:

$$\begin{array}{l}{\omega }_{{\bf{q}}\nu }^{2}(T)={({C}_{{\bf{q}}}^{\dagger }{\varLambda }_{{\bf{q}}}^{{\rm{HA}}}{C}_{{\bf{q}}})}_{\nu \nu }+\frac{1}{2}\sum _{{{\bf{q}}}^{{\prime} },{\nu }^{{\prime} }}{\varPhi }^{{\rm{SCP}}}(-{\bf{q}}\nu \,;{\bf{q}}\nu \,;-{{\bf{q}}}^{{\prime} }{\nu }^{{\prime} }\,;{{\bf{q}}}^{{\prime} }{\nu }^{{\prime} })\\ \,\,\,\,\times \frac{\hbar (1+2n({\omega }_{{{\bf{q}}}^{{\prime} }{\nu }^{{\prime} }}))}{2{\omega }_{{{\bf{q}}}^{{\prime} }{\nu }^{{\prime} }}},\end{array}$$

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

$${Z}_{\nu ,\alpha }^{* }=\sum _{\kappa \beta }{Z}_{\kappa ,\alpha \beta }^{* }\frac{{e}_{\nu }^{\beta }\left(\kappa \right)}{\sqrt{{m}_{\kappa }}},$$

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:

$$V\left({Q}_{1},{Q}_{2}\right)=\frac{1}{2}{\omega }^{2}{Q}_{1}^{2}+\frac{1}{2}{\omega }^{2}{Q}_{2}^{2}+\frac{1}{4}\text{k}{Q}_{1}^{4}+\frac{1}{4}\text{k}{Q}_{2}^{4}+\chi {Q}_{1}^{2}{Q}_{2}^{2}+\psi {Q}_{1}^{3}{Q}_{2}+\psi {Q}_{1}{Q}_{2}^{3},$$

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:

$${\ddot{Q}}_{i}+\frac{\partial V}{\partial {Q}_{i}}+\varGamma {\dot{Q}}_{i}={Z}^{* }{\widetilde{E}}_{i}^{{\rm{THz}}},\qquad i=1,2,$$

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:

$$\mu =\gamma Q\times \dot{Q}=\sum _{i}{\gamma }_{i}{Q}_{i}\times {\dot{Q}}_{i}=\sum _{i}{{\gamma }_{i}L}_{i},$$

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¹⁷.