Introduction

Owing to recent advances in laser technology it is now possible to access nonlinear phonon dynamics, which can be used to control physical properties of solids. For example, an intense ultrashort optical pulse can generate coherent lattice vibrations with large displacements via the impulse Raman-scattering process1,2,3,4. By using the phonon anharmonicity under photoexcitation, several methods for the control of material properties have been explored. In particular, owing to wavelength conversion techniques5, infrared-active phonons can be driven efficiently6,7,8. Such strongly driven phonons can cause an atomic rearrangement in the crystal and can change the material functionality, e.g. they can reveal an otherwise hidden ferroelectricity in quantum paraelectric materials9. The driven phonons can also induce exotic electron phenomena in strongly correlated systems, which is important for the field of nonlinear phononics10,11,12,13. These phenomena can be enhanced by an electric-field enhancement with metallic sub-wavelength periodic structures14.

Here, we propose another method to control materials by phonon excitation: we consider the possibility of a martensitic transformation that is triggered by the nonlinear response of phonons under excitation with terahertz (THz) pulses. The term martensitic transformation describes a particular crystal-structure transformation that occurs via shear deformation. The obtained crystal structure is exceptionally hard and the transformation is usually induced by rapid thermal quenching of suitable alloys. This transformation has been intensively investigated in steel alloys by applying mechanical stress15 or magnetic fields16, and has been widely used for shape memory alloys17. Since softening of transverse acoustic (TA) phonons with large wavenumbers has been reported in the premartensitic phase17,18, TA phonons are considered to play an important role in the transformation. Acoustic phonons at the zone boundary are accessible via the Klemens process from the \(\Gamma\)-point transverse optical (TO) phonon19,20. Note that the Klemens process has also been identified as a decay channel of coherent phonons21,22. Thus, we expect that large-amplitude TA phonons generated from resonant TO phonon excitation via the Klemens process can cause a martensitic transformation in certain alloys. To demonstrate this, we use tetragonal zirconia (t-ZrO2). To populate the TA phonon mode that triggers the martensitic transformation, we use a THz free electron laser (THz-FEL). Recently, crystalline-amorphous phase switching of Ge2Sb2Te5 has been demonstrated using this THz source (in Ge2Sb2Te5, photoinduced phase transitions are also possible)23. Additionally, a high-efficiency elastic wave generation at the surface of water has been reported using the same THz source24,25, suggesting that generation of a mechanical force in t-ZrO2 is possible.

In this work, we demonstrate the martensitic transformation in t-ZrO2 induced by the THz excitation resonant for the lowest infrared-active phonons. The morphology of the irradiation scar, Raman scattering spectra, the photographs of transmission electron microscopy (TEM), and the excitation frequency dependence of the irradiation scars provide clear evidence for a THz-induced martensitic transformation, which is not accessible by thermal effects.

Results

Verification of suitable phonon modes

Zirconia is an ionic material with a large phonon anharmonicity and it has been shown that zirconia doped with about 3 mol% Y2O3 consists mainly of a meta-stable tetragonal phase at room temperature26,27. In this work, we used commercial partially stabilized zirconia ceramic plates (3.5 mol% Y2O3) with a high volume-fraction of t-ZrO2. t-ZrO2 becomes cubic (c-ZrO2) at higher temperatures and can become monoclinic (m-ZrO2) due to mechanical strain26,27,28. Figure 1a shows the crystal structures of ZrO2 (different projections are shown in Supplementary Fig. 1). The skeleton of the m-ZrO2 lattice exhibits a shear distortion of 9° parallel to the basal plane of the unit cell of t-ZrO226. Since the tetragonal to monoclinic transformation is accompanied by a volume expansion, the microcracks that are normally caused by the deformation of a polycrystalline material, are suppressed via this stress-induced transformation of t-ZrO226,27,28,29. Figure 1b shows the infrared reflection and Raman spectra of one the zirconia ceramic plates. Note that t-ZrO2 has infrared phonon modes at 4.2 THz (Tu), 10.2 THz (A2u) and near 13 THz (Tu)30. In this work, either the infrared phonon mode at 4.2 THz or that at 10.2 THz was excited resonantly by appropriately tuning the THz-FEL parameters (see Methods). The X-ray diffraction (XRD) and scanning electron microscopy (SEM) data of our sample are provided in Supplementary Fig. 2.

Fig. 1: Phase transformation of t-ZrO2.
figure 1

a The crystal structures of t-, c-, and m-ZrO2. The green and red spheres correspond to Zr and O atoms, respectively. b Infrared reflection spectrum (red curve) and Raman spectrum (blue curve) for the sintered 3.5-mol% Y2O3-stabilized zirconia plate, which mainly shows characteristics of the tetragonal phase. The imaginary part of the dielectric constant was determined from the reflection spectrum with the assumption of four Lorentz oscillators.

Full size image

Irradiation with a single THz macropulse that is resonant with the 4.2-THz phonon

For this experiment, a partially stabilized zirconia ceramic plate was irradiated with a single macropulse as shown in Fig. 2a. The excitation pulse had a spectrum with a peak at 4.2 THz, which is resonant with the lowest optical phonon. Each photograph in Fig. 2b–e shows a scar on the surface of the t-ZrO2 plate generated by a macropulse with a different energy. The photograph of the whole plate is shown in Supplementary Fig. 3. Irradiation scars were already clearly visible above macropulse energies of 2 mJ (3 J cm−2). We confirmed that field electron emission does not occur under these experimental conditions by verifying the absence of electron-emission-induced luminescence from the air around the sample31. Figure 2g shows that, when a macropulse energy of 29 mJ is used, the height of the center of the irradiation spot lies below the initial level while the outer regions lie above. These features indicate the occurrence of melting; we interpret that the t-ZrO2 was melted in the center by thermal heating and a strong elastic wave propagated towards the outer region, and the rebound of the elastic wave results in the higher region around the depressed center area. This suggests that the maximum temperature reached for the 7.5-mJ macropulse is much lower than the melting point. Nevertheless, the height of the center of the irradiation spot of t-ZrO2 rises as shown in Fig. 2f. There are several possible explanations for this latter result, including a phase transformation of zirconia and crack formation due to rapid volumetric expansion and contraction. In other words, changes in the surface morphology do not necessarily imply a phase transformation, and therefore this data needs to be evaluated in conjunction with another method such as Raman spectroscopy.

Fig. 2: Scars on the surface of the t-ZrO2 sintered plate after irradiation with 4.2-THz pulses.
figure 2

a Schematic of our experimental setup. b–e Photographs of scars generated by single 4.2 THz pulses with different excitation pulse energies. f, g Two photographs of scars and the corresponding height profiles of the irradiation scar at the position indicated by the red line. Red and blue hatches represent rise and depressions in the irradiation spot, respectively.

Full size image

In order to investigate the THz-induced crystal-structure change, we performed Raman spectroscopy using a confocal microscope. Figure 3a shows the Raman spectra obtained at the centers of irradiation scars generated by different macropulse energies. For the region without irradiation, Raman peaks that are characteristic of t-ZrO2 appear30,32. On the other hand, for the irradiation scar generated by a macropulse with an energy of 7.5 mJ (Fig. 3a; red curve), new Raman peaks appear at 170, 184, and 369 cm−1. These peaks are characteristic of m-ZrO230,32, suggesting that a phase transformation contributed to the overall height increase in Fig. 2f. Above excitation pulse energies of 20 mJ (Fig. 3a; blue curve), the Raman peaks that are characteristic of m-ZrO2, are less pronounced and all Raman peaks become broad due to crystal melting.

Fig. 3: Characterization of the atomic structure at different irradiation spots.
figure 3

a Raman spectra at the center positions of irradiation scars generated by excitation pulse energies of 29 (blue curve) and 7.5 mJ (red curve). The reference Raman spectrum without the irradiation (W/O) is shown by the black curve. The lower panel shows the Raman-active modes of t- and m-ZrO230. b TEM image of the t-ZrO2 sample before irradiation. c TEM image of a sample region that was irradiated with a weak 4.2-THz FEL pulse.

Full size image

We also used TEM to evaluate the changes of the atomic arrangement in the irradiation scars. Figure 3b is a TEM image of the original partially stabilized zirconia, indicating a tetragonal atomic arrangement. Figure 3c shows a TEM image of a sample region that was irradiated with a THz pulse using a pulse energy of 8 mJ. In Fig. 3c, there are many areas with an apparently tilted atomic arrangement, which is characteristic of monoclinic crystals. We also performed XRD measurements as shown in Supplementary Figs. 4 and 5. Because the X-ray beam spot was much larger than the irradiation scar, the XRD measurements only provided averaged information. In addition, since partially stabilized zirconia contains many twins, it is not possible to extract the monoclinic features by conventional Rietveld analysis33,34. However, by using the two peak intensities that characterize the monoclinic crystal structure as indices, results consistent with the TEM images were obtained.

We also investigated the irradiation scar on an 8-mol% Y2O3-stabilized zirconia sintered plate, which is intrinsically cubic and does not exhibit a phase transformation under pressure or exposure to high temperatures27. The photographs, height profiles, and Raman spectra of the irradiation scars for different macropulse energies are shown in Supplementary Fig. 6. The Raman spectroscopy data clarifies that a phase transformation has not occurred in this sample. The presence of innumerable peeling structures on the surface suggests a thermal volume expansion and its rebound without phase transformation.

Temporal evolution of the surface roughness

The changes in the surface condition of the irradiation spot were investigated in the time domain by detecting the intensity of a visible laser reflected from the irradiation spot. Figure 4a shows the experimental setup. In this experiment, we reduced the macropulse width to 1 μs to improve the temporal resolution. Fifty consecutive shots of a weak macropulse (repetition rate 5 Hz) were used and we measured the temporal evolution of the change in the reflected intensity, \(\varDelta I\), for each macropulse. Figure 4b shows the photographs of the three irradiation scars generated by the irradiation with fifty macropulses at three different pulse energies (0.99, 1.7, and 2.4 mJ). Note that irradiation scars appeared even below the threshold macropulse energy of 2 mJ for the single-shot excitation, implying a partial phase transformation even by low-energy macropulses. Figure 4c shows the three time profiles of the changes in \(\varDelta I/{I}_{0}\) observed after the first macropulse for 0.99, 1.7, and 2.4 mJ, respectively (\({I}_{0}\) is the initial reflected-beam intensity. Here, \(\varDelta I/{I}_{0}\) decreases sharply within the macropulse width, indicating that the THz pulse irradiation causes a rapid volume expansion. The reflected intensity decreases, because the surface roughness increases as shown in Fig. 2f. In all three curves, this reduction of \(\varDelta I/{I}_{0}\) continues up to 10 μs after excitation. This is the time required by all modes, including high-frequency phonons, to reach thermal equilibrium. The inset plots the macropulse energy dependence of \(\varDelta I/{I}_{0}\) measured at a time of 10 μs after the excitation.

Fig. 4: THz-induced temporal evolution of the surface roughness.
figure 4

a The experimental setup. The green arrows indicate the reflection measurement using visible light. b The photographs of the three spots after irradiation with the fifty THz macropulses for the three pulse energies used (0.99, 1.7, and 2.4 mJ). The contrast of these are adjusted to show the changes in the irradiation spot clearly. The length of the black scale bar corresponds to 100 μm. c The time dependences of the changes in \(\varDelta I/{I}_{0}\) observed after the first macropulse for 0.99, 1.7, and 2.4 mJ. The inset shows the macropulse energy dependence of \(\varDelta I/{I}_{0}\) at a time of 10 μs after excitation. d Dependences of the irreversible component in \(\varDelta I/{I}_{0}\) on the number of macropulse shots for 0.99, 1.7, and 2.4 mJ. For this, we plotted the \(\varDelta I/{I}_{0}\) observed immediately before the following macropulse as shown by the arrow in the inset. The inset shows the \(\varDelta I/{I}_{0}\) profile in Fig. 4c for 0.99 mJ.

Full size image

These observed changes of \(\varDelta I/{I}_{0}\) in the early time region were mainly reversible and they decayed within several milliseconds. We extracted the magnitude of the THz-induced irreversible change of each macropulse by measuring the \(\varDelta I/{I}_{0}\) that occurs immediately before irradiation with the next macropulse about 200 ms later. Figure 4d plots the dependences of the irreversible \(\varDelta I/{I}_{0}\) component on the number of macropulse shots. For pulse energies larger than 2 mJ, the irreversible \(\varDelta I/{I}_{0}\) component is positive after the first pulse and it becomes negative in subsequent shots. For 1.7 mJ, a gradual increase of the irreversible \(\varDelta I/{I}_{0}\) component with the shot number can be confirmed. A plausible interpretation is that a partial phase transformation compensated the voids in the ceramic27.

Irradiation with laser pulses that are not resonant with the lowest phonon

Figure 5a, b shows the irradiation scars generated by irradiation with ultraviolet (UV) and near-infrared (near-IR)-femtosecond pulses, respectively. Clear traces of ablation due to electron heating and evaporation are observed. The black spots are characteristic of oxygen-deficient zirconia35. These results reproduce the results of a previous report, in which laser ablation of zirconia was investigated from the aspect of laser machining and laser pulse deposition techniques36. In spite of the large bandgap energy (>5 eV)35, laser-induced plasma is generated via multi-photon absorption, which accelerates the desorption of oxygen atoms. Figure 5c shows the irradiation scar generated by irradiation with mid-infrared (mid-IR) pulses. The large cracks on the surface are due to thermal stress. We found that similar cracks are also apparent in the irradiation scar generated by THz-FEL pulses with a center frequency of 10.2 THz (Fig. 5d). The peeling microstructure in Fig. 5d is similar to that of the irradiation scar on the 8-mol% Y2O3-stabilized zirconia sintered plate (Supplementary Fig. 6). Therefore, we consider that the microcracks generated by the 10.2-THz-FEL pulses were caused by rapid thermal volume expansion and not primarily due to a phase transformation.

Fig. 5: Characterization of the irradiation scars generated by four different lasers.
figure 5

The photographs of scars generated by UV (a) 4 mJ (28 J cm−2), near-IR-femtosecond (b) 0.2 mJ (0.2 J cm–2), mid-IR (c) 0.6 mJ (3 J cm–2) and 10.2-THz-FEL (d) 0.2 mJ (8 J cm–2) irradiation, respectively. The data below the photographs show the height profiles of the irradiation scars at the positions indicated by the red line. Red and blue hatches represent rise and depressions in the irradiation spot, respectively. e Raman spectra at the center positions of the irradiation scars generated by UV, near-IR-femtosecond, mid-IR and 10.2-THz-FEL pulses. The red shaded region shows the Raman spectra at the center position of the irradiation scar generated by the 4.2-THz FEL pulse with a pulse energy of 7.5 mJ (12 J cm–2) (Fig. 3). The reference Raman spectrum without the irradiation (W/O) is shown by the black curve. Black dots and red crosses indicate the wavenumbers of the Raman-active modes of t- and m-ZrO228.

Full size image

Figure 5e shows the Raman spectra at the center positions of the irradiation scars generated by the UV, near-IR-femtosecond, mid-IR, and 10.2-THz FEL pulses. Tiny Raman peaks that are characteristic of the monoclinic phase can be confirmed, which indicates a phase transformation due to stress induced by rapid heating and/or cooling. However, these peak intensities are much weaker than those observed in the spectrum for the 4.2-THz-FEL irradiation (the red shaded area in Fig. 5e). Therefore, we conclude that the martensitic transformation to the monoclinic phase occurs after the irradiation with the 4.2-THz FEL pulse because the required optical phonon mode is directly excited, while this transformation cannot be efficiently induced by using the UV, near-IR-femtosecond, mid-IR, and 10.2-THz FEL pulses.

Discussion

Our experimental results suggest that the irradiation with the 4.2-THz FEL pulse induced a martensitic transformation to the monoclinic phase. This transformation observed for relatively weak pulse energies is different from the thermal stress-induced martensitic transformation. Note that the irradiation with the 10.2-THz FEL pulse, which enables resonant excitation of a different phonon mode in the crystal, did not induce a phase transformation with similar efficiency, but formed microcracks due to rapid photothermal expansion and contraction. This is direct evidence of the unique nature of the excitation using the THz pulse resonant with the 4.2-THz optical phonons.

To explain this phenomenon, we focus on TA phonons with large wavenumbers, and not on optical phonons, because the former phonons are also important for martensitic transformation in metal alloys18. A prominent difference between the t- and m-ZrO2 crystal structures in Fig. 1a is the tilt in the lattice skeleton, which corresponds to a displacement described by TA phonons at the zone edge. Two mechanisms have been proposed to explain the martensitic phase transformation in zirconia involving phonons at zone edge. One is the condensation of two phonons at the M-point37. The other mechanism is attributed to a phonon instability at the Z-point38. Numerical simulations have suggested that the frequency of the phonon at the Z-point can become imaginary in a certain volume depending on the pressure39. Therefore, selective relaxation from the optical phonon to either a phonon at the M-point or to a phonon at the Z-point is required to interpret our experimental results. In general, phonon relaxation is governed by phonon–phonon scattering and phonon–dislocation scattering40. The excitation frequency dependence of the phase transformation shown in our work implies that higher-order phonon scattering and phonon–dislocation scattering are not significant. Only the lower-order phonon–phonon scattering at low frequencies restricts the relaxation channel (according to the wavenumber conservation law). The lowest-order process is the Klemens process. Figure 6a shows the schematic of the Klemens process, which has been intensively investigated by ultrafast spectroscopy21,22. Since its efficiency is determined by the phonon density of states, relaxation to zone-edge phonons becomes stimulated.

Fig. 6: The phonon modes of t-ZrO2 that are related to the martensitic transformation.
figure 6

a The schematic of the Klemens process. b The phonon branches are shown along high-symmetry directions in the Brillouin zone. The open red circles show the phonon modes at half of the frequency of the infrared-active TO phonon mode shown by the closed blue circle in (b). In the right panel, the total phonon density of states in t-ZrO2 is indicated by the black line. The partial phonon density of states (DOS) of the Zr and O atoms are indicted by the green and red lines, respectively. Ion displacements induced by infrared-active TO phonon modes at the Γ-point (c); 4.73 THz and (d); 10.5 THz), and those induced by the lowest phonon mode at the Z-point (e), 1.99 THz). These three phonons are shown by the closed blue circle, the closed blue square, and the closed red circle in (b), respectively.

Full size image

To verify this model, we determined the phonon dispersion curves of t-ZrO2 by using density functional theory calculations within the strongly constrained and appropriately normed (SCAN) meta-generalized gradient approximation (SCAN meta-GGA)41,42. Figure 6b shows the calculated phonon dispersion curves along the major high-symmetry directions and the corresponding the phonon density of states. These results almost reproduce the previously reported data30. Figure 6c, d visualizes the atomic displacements due to the infrared-active TO phonons at 4.73 (Tu) and 10.5 THz (A2u), which are indicated by the closed blue circle and the closed blue square in Fig. 6b, respectively. These two phonon modes cannot induce shear deformation. Therefore, we consider the relaxation channels for these optical phonons. There are many relaxation channels for the 10.5-THz phonon, while those for the 4.73-THz phonon are considerably fewer. The phonons accessible via the Klemens process are shown by the open circles in Fig. 6b. The relaxation to other phonons near the Γ point (Ridley process43) is also possible, but the contribution of this channel is small due to the small phonon density of states. Among those phonons accessible via the Klemens process, we can confirm that there are several low-frequency TA phonons at the zone boundary near the Z-point (along Z–A and Z–R), but not near the M-point. The mode number at the final state of the Klemens process is large owing to \(\partial \omega /\partial k=0\) at the zone edge. Moreover, because there is almost no dissipation by elastic wave propagation during pulse-train excitation, an accumulation of phonons will occur. Figure 6e shows the atomic displacements induced by the low-frequency TA phonons near the Z-point. This phonon mode can also cause a tilting motion of ions that is matched to the crystal skeleton of t-ZrO2. However, since the phonons at the zone boundary have relatively short wavelengths, an additional mechanism is required for the macroscopic phase transformation. Note that softening of TA phonons with large wavenumbers has been reported in premartensitic phase in metal alloys17,18. Therefore, we consider that a similar softening of TA phonons also occurred in our demonstration.

The above discussion identified the importance of the low-frequency TA phonons near the Z-point, but these phonons might also induce a transformation to a different phase in which the tilt direction of each unit cell changes. The transformation from t- to m-ZrO2 requires a translational shift of the crystal cell as well as a crystal tilt44. The destructive spallation of the zirconia microcrystal seen in Fig. 2b, c shows the importance of the translational shift in the phase transformation. In Fig. 6b, we have shown that some TA and longitudinal acoustic (LA) phonons with small wavenumbers are also generated via the Klemens process. Their amplitudes are smaller than those of the low-frequency TA phonons near the Z-point, but the translational shift related to these modes is sufficient to assist the transformation to m-ZrO2. Thus, we conclude that the phonons generated via the Klemens process can trigger the phase transformation.

In conclusion, we have demonstrated the THz-induced martensitic transformation in a t-ZrO2 sintered plate. From the height profiles and the Raman spectroscopy data of the irradiation scar generated by the THz-FEL pulse resonant with the 4.2-THz optical phonon, we found that the drastic volume expansion due to the martensitic transformation causes a destructive spallation of the zirconia microcrystal. This martensitic transformation is caused by shear deformation and a translational shift due to an acoustic phonon that is generated from the optical phonon via the Klemens process. We calculated the phonon modes in t-ZrO2 and confirmed an effective channel from the infrared-active phonon mode to the mode required for the transformation. Since the THz pulse excitation can access specific channels of the transformation (which is in contrast to the conventional martensitic transformation via thermal expansion and/or mechanical stress), it can be used to elucidate complex martensitic transformation mechanisms in various materials with dynamical approaches. Our demonstration sheds light on the potential of nonlinear phononics in ferroelastic materials for crystallization control, material strengthening and strain-control of shape memory ceramics45.

Methods

Sample characterization

We used a commercial 3.5-mol% Y2O3-doped zirconia ceramic plates (MTI corporation) with dimensions of 10\(\times\)10\(\times\)0.5 mm and a polished surface. We characterized the samples with infrared reflection spectroscopy (Fig. 1b), Raman spectroscopy (Fig. 1b), and XRD measurements and SEM (Supplementary Fig. 2 in the Supplementary Information).

Terahertz source

For the irradiation of the sample with THz pulses, we used a THz-FEL based on the L-band electron linear accelerator at the Research Laboratory for Quantum Beam Science, Institute of Scientific and Industrial Research, Osaka University46. The THz-FEL generates a linearly polarised THz pulse train consisting of ca. 30–150 pulses (THz-FEL micropulses) with a micropulse interval of 37 ns. This forms the THz-FEL macropulse with a duration of 1–6 μs. The macropulse repetition rate of this system is 5 Hz. By electro-optic sampling using a synchronized Ti:sapphire laser with pulse duration of 100 fs, we determined that the micropulse duration can be controlled in the range from 2 to 6 ps31. We carefully confirmed that the observed phenomena are insensitive to the macropulse and micropulse widths. If a 15-MeV electron beam is used, the center frequency of the FEL pulse can be varied from 2.8 THz to 8.0 THz by changing the gap length of the wiggler magnet. In the present experiments, the peak frequency of the THz-FEL was tuned in the range from 4.2 to 5.2 THz, because the irradiation scars appeared above frequencies of 3.8 THz. As the THz-FEL characteristics depend on the optical cavity length, the sample was set at the position where the macropulse energy was the highest. We summarized the data of the THz-FEL pulses in Supplementary Fig. 7. The THz-FEL was focused by a 50-mm gold-coated 90° off-axis parabolic mirror. Since the beam diameter was 20 mm, the spot size at the focal point was about 0.2 mm. During the irradiation experiments, the sample was moved with a constant speed (5 mm s−1) along the y axis to obtain many individual scars on the surface, each generated by a single macropulse. In addition, by using an 18-MeV electron beam, we generated 10.2-THz FEL pulses for a reference experiment. We summarized the details of the THz-FEL pulses based on the 18-MeV electron beam in Supplementary Fig. 8.

Scar characterization

The irradiation scars on the sample were characterized by using a digital microscope (Keyence, VHX6000) and a confocal microscope (OPTELICS, H1200). The structure of the samples was characterized by XRD (Rigaku, UltimaIV, CuKα). Since the X-ray irradiation spot had an extension of \(1\times 5\) mm, we prepared a surface with many irradiation scars as explained above. The photograph of the sample is shown in Supplementary Fig. 4a. XRD measurements were used to characterize two types of scars: The scars that include a melted region were created by irradiating a 40-mJ single macropulse with a spot size of 0.2 mm. The scars without melted regions were created by irradiating a 40-mJ single macropulse with a spot size of 0.5 mm. We also characterized the scars by cross-sectional high-resolution transmission electron microscopy (HR-TEM: JEOL, JEM-ARM200F). In this measurement, we used the plate shown in Supplementary Fig. 3. For the Raman spectroscopy, we used a commercial Raman spectroscopy system (Tokyo Instruments, Nanofinder30) with a 488-nm laser (COHERENT, Sapphire-488-20). In order to reliably measure the Raman spectrum at the irradiation scar generated by the 10.2-THz FEL pulse, we prepared alignment marks by photoresist techniques. Owing to the high chemical resistance of zirconia, no phase transformation occurred during this process. A photograph of the sample with the alignment marks is shown in Supplementary Fig. 9. In the measurement of the surface roughness change in the time domain, a 532-nm beam from a diode-pumped solid-state green laser (LSR532NL) was focused on the THz irradiation spot by a lens with a focus length of 300 mm and the reflected beam was detected by a photodiode (Thorlabs, DET10A) at a distance of 300 mm from the sample without any lens.

The laser sources for the reference experiments

For the UV irradiation, we used the 4th harmonic of a pulsed YAG laser (Spectra-Physics, Quanta-Ray) with a center wavelength of 266 nm, a repetition rate of 10 Hz, an output pulse energy of 8 mJ and a pulse duration of 10 ns. For the near-IR-femtosecond irradiation, we used an amplified Yb: YAG laser (HighQ laser, SC-1042-1000fs Reg Amp HE) with a center wavelength of 1.04 µm, a repetition rate of 1 kHz, an output pulse energy of 1.0 mJ and a pulse duration of 440 fs. For the mid-IR irradiation, we used a CO2 laser (Synrad 48-1) with a center wavelength of 10.6 µm, a repetition rate of 5 kHz, an output pulse energy of 0.6 mJ and a pulse duration of 200 μs. Also in this experiment, the sample was moved during the irradiation experiment by a motorized x–y stage.

Calculation of the phonon dispersion curves

The phonon dispersion of t-ZrO2 was derived by using density-functional theory within the SCAN meta-GGA41,42. The calculations were performed by applying finite displacements in the supercell using VASP47,48,49,50,51 and PHONOPY52. Convergence with respect to the k-point sampling and the energy cut-off values (Ew = 400 eV and Ec = 605.4 eV for the wave function and the augmentation charge, respectively) for the plane-wave expansion was confirmed, and a uniform k-point mesh of 8 × 8 × 8 in the 1st Brillouin zone was used. For the force field calculation, a 4 × 4 × 4 supercell was adopted.