All-optical control and visualization of ultrafast two-dimensional atomic motions in a single crystal of bismuth

In a bulk solid, optical control of atomic motion provides a better understanding of its physical properties and functionalities. Such studies would benefit from active control and visualization of atomic motions in arbitrary directions, yet, so far, mostly only one-dimensional control has been shown. Here we demonstrate a novel method to optically control and visualize two-dimensional atomic motions in a bulk solid. We use a femtosecond laser pulse to coherently superpose two orthogonal atomic motions in crystalline bismuth. The relative amplitudes of those two motions are manipulated by modulating the intensity profile of the laser pulse, and these controlled motions are quantitatively visualized by density functional theory calculations. Our control-visualization scheme is based on the simple, robust and universal concept that in any physical system, two-dimensional particle motion is decomposed into two orthogonal one-dimensional motions, and thus it is applicable to a variety of condensed matter systems.

A tomic motions in the solid state can be triggered with optical pulses. If the pulse duration is shorter than the periods of lattice vibrations, the atoms start oscillating collectively. This collective oscillation is referred to as coherent phonons 1,2 . Coherent phonons have been observed in a wide variety of materials [3][4][5][6][7][8][9][10][11][12][13][14][15][16] and have been controlled onedimensionally with a train of multiple laser pulses [3][4][5][6][7][8][9] . One of these control experiments was demonstrated with two different phonon modes. The directions of atomic motions are, however, parallel between those two different modes, so that the control was one-dimensional 8 . Two-dimensional (2D) control has been performed with a pair of identical laser pulses with different polarizations applied to energetically degenerate orthogonal atomic motions 10 . Another approach is necessary to control orthogonal non-degenerate atomic motions, which is more common in most materials. Such 2D control should be important in photo-induced processes 17,18 and could lead to further development of optically engineered functionality of bulk solids 19,20 .
Here we demonstrate 2D control of atomic motions in crystalline bismuth. It is known that superposing two chirped femtosecond laser pulses causes intensity modulation whose characteristic period lies within the terahertz (THz) region 21,22 . This method is applied to atomic motions in the current study. Two orthogonal phonon modes of A 1g and E g symmetries 23 are superposed with the modulated pulse and the transient relative reflectivity change of the crystal surface is measured with another femtosecond probe pulse. The reflectivity changes show beat structures whose temporal evolution changes drastically as we change the delay between the chirped femtosecond pulses on the attosecond timescale. We map these beat structures into a 2D space of atomic displacements with our density functional theory (DFT) calculations, demonstrating direct 2D control of atomic motion in each unit cell and its visualization.

Results
Coherent phonons in crystalline bismuth. The crystal structure of bismuth is rhombohedral 24 and its primitive unit cell consists of two atoms as shown in Fig. 1a. There exist three optical phonon modes with A 1g and doubly degenerate E g symmetries at the G-point in the reciprocal lattice space 23 . According to our definition of the axes as given in Fig. 1a, the A 1g mode corresponds to the longitudinal motion and the E g mode corresponds to the lateral motion of the Bi atoms.
Experimental setup. Figure 1b schematically illustrates our pump-probe experiments performed on the (0001) surface of a single-crystal Bi. The crystal was cooled in a liquid helium cryostat to 5.2 K to reduce damping and dephasing of the coherent phonons. The output of a Ti:sapphire laser (centre wavelength: B802 nm, bandwidth: B34 nm) was split with a partial beam splitter into two beams with 9:1 ratio, which were used as pump and probe pulses, respectively. The pump pulse was input to our homemade highly stabilized Michelson-type interferometer 25,26 to produce a phased pair of laser pulses. The intensity of each pulse was B6 mW (fluence: B2 mJ cm À 2 on the sample). These pulses (pulse width: B30 fs in their Fourier transform (FT) limits) were chirped with transmitting optics inserted on their optical path and are hereafter referred to as chirped subpulses (CSPs) 27 . The group delay dispersion (GDD) parameter of each CSP is hereafter referred to as f 00 . Those CSPs were temporally superposed with their delay t mod to produce a modulated pump pulse whose temporal intensity envelope I(t) shows periodical undulations. The FT of I(t) shows a single peak in the THz region where its centre frequency and peak width are functions of t mod and f 00 . A combination of these two degrees of freedom in the centre frequency and peak width allows for arbitrary tuning of the relative intensities of any two different frequency components in the THz region. Superposing CSPs thus serves as a practical technique for THz intensity modulation with near-infrared laser pulses. The train of femtosecond pulses has often been used in the previous experiments to control coherent phonons one-dimensionally [3][4][5][6][7][8][9] . The difference between such a pulse train and the current CSPs will be described later from the viewpoint of 2D control. The polarization of the modulated pump pulse was rotated by 90°, so that the pump and probe pulses were focused on the Bi sample with their polarizations perpendicular to each other. The reflectivity change DR of the probe pulse induced by the modulated pump pulse was monitored with a pair of balanced photodiodes (PD1 and PD2). The differential signal from the PDs was amplified with a current amplifier (SRS SR570) and averaged in a digital oscilloscope. The temporal evolution of DR was measured by scanning the delay t probe of the probe pulse from the first CSP repetitively at 20 Hz (APE Scan Delay) with a fixed t mod . A bandpass filter from 3 to 300 kHz was used in all of the measurements reported in this study, so that the nonoscillatory background in the temporal evolution, which corresponds to the featureless electronic response, was filtered out. The intensity of the modulated pump pulse was monitored with a photodiode (PD3) at each t mod . Other details of signal processing are given in the Methods section.
The temporal full width at half maximum (FWHM) of the probe pulse at the sample was measured to be B50 fs, giving its GDD parameter f 00 B430 fs 2 . The FWHM and GDD parameter of the CSPs were estimated from the cross-correlation between one of the CSPs and the probe pulse to be B260 fs and f 00 B2,600 fs 2 , respectively.
Control of 2D atomic motions. Figure 1c shows an example of the temporal evolution DR(t probe )/R with t probe scanned and t mod fixed to B0 fs. The beat structure with a period of B330 fs arises from coherent phonon motion and is described with a linear combination of two damped harmonic oscillators corresponding to A 1g (longitudinal) and E g (transverse) phonons 23 : where A 0 is the baseline, and DR i /R, G i , n i and d i with {i} ¼ {a, e} are the phonon amplitudes, damping factors, oscillation frequencies and initial phases of A 1g and E g phonons, respectively.
Assuming each CSP has a Gaussian spectral profile with the GDD parameter f 00 , the spectrum of the electric field is described as 28 where o L is the central angular frequency of the pulse and Do is its spectral FWHM. The modulated pump pulse is the superposition of two CSPs with their interpulse delay t mod , so that its intensity envelope I(t) is ARTICLE where E þ CSP ðtÞ is given by The envelope I(t) is Fourier-transformed to give its THz spectrum I(o) asĨ IðtÞe iot dt: ð5Þ Figure 2a shows the intensity envelopes I(t) of the modulated pump pulse simulated with three different CSP delays t mod ¼ 49.4, 92.2 and 93.6 fs combined with the CSP parameters o L , Do and f 00 , taken from the current experimental conditions. The dotted curve represents the intensity envelope of a CSP pair with the delay t mod ¼ 0 fs. Figure 2b shows the THz spectraĨ(o) Fourier-transformed from those simulated envelopes I(t), demonstrating that the spectrumĨ(o) is modulated by changing the CSP delay t mod , so that we can control the relative intensities of two different THz components resonant with A 1g and E g phonons, respectively. Figure 2c shows the observed temporal evolutions DR(t probe )/R with three different CSP delays t mod similar to those for the simulations shown in Fig. 2a,b. It is clearly seen that the beat structure changes drastically as we change the CSP delay. The FTs of the traces in Fig. 2c are plotted in Fig. 2d in which two peaks at B3.0 and B2.1 THz correspond to the A 1g and E g phonons, respectively, demonstrating that their relative intensities are clearly controlled.
Visualization of 2D atomic motions. The real-space atomic motions have been reconstructed in Fig. 2e-g from the measured beats shown in Fig. 2c, based on the correspondence between the nuclear displacement 29-31 and the reflectivity change DR. This correspondence has been obtained by linear response DFT calculations 32 with the full-potential linearized augmented plane wave programme WIEN2k 33,34 , which is among the most accurate DFT codes presently available. The two quantities that we have computed are the potential energy surface as a function of the absorbed laser energy and optical properties in dependence of the Bi atomic coordinates. Details of our total-energy computations have been published in ref. 35. For the optical calculations, we have used the module described in refs 36 and 37, which is based on the Lindhard equation. With this package, we computed the imaginary part of the dielectric tensor for photon energies up to 27 eV. We obtained the real part using the Kramers-Kronig relation, where we added an imaginary part of 0.1 eV to the energy to account for the finite lifetime of excited states. It is important to realize that both our total-energy computations and the optical calculations rely on the interpretation of the Kohn-Sham energies as single-electron excitation energies 34,36,37 , which, although not strictly valid within DFT, is widely used and appears to work well for many materials including bismuth 34,37 . As the optical calculations are very sensitive to the quality of the k mesh, in this part of our theoretical work we have included 40 Â 40 Â 40 k points.
To check the validity of our approach, in Fig. 3 we have plotted a number of calculated dependencies and we have compared our predictions to published experimental data. In particular, the quasiequilibrium coordinate z eq , which is associated with the A 1g phonons, changes with increasing pump-laser fluence. Simultaneously, the A 1g phonon frequency softens. In Fig. 3a, we have plotted the functional dependence of both effects, which can also be experimentally assessed through time-resolved X-ray spectroscopy 29,30 . Clearly, the agreement with these experiments is very good. In Fig. 3b, we show the computed optical conductivity for the electronic ground state. It qualitatively reproduces the essential features of the measured data of ref. 38, in particular it has a maximum around 0.95 eV, a shoulder at 1.5 eV, a dip at 2.1 eV, followed by a maximum and a decay. At the photon energy of 1.55 eV (800 nm), to which we shall restrict ourselves in the following, the computed optical conductivity is in good quantitative agreement with the data of ref. 38. As a final check, in Fig. 3c, we show how the average reflectivity (sometimes called background reflectivity, because it ignores the superimposed oscillatory effect due to excited coherent phonons) and the A 1g phonon frequency change simultaneously. Our prediction of Fig. 3c is compared with optical measurements of refs 4,39,40. The good agreement between our calculations and the experimental data confirms that our approach reliably reproduces the experimental results.
We now turn to a description of our computed results for the optical reflectivity R of Bi at a wavelength of 800 nm. For the ground state (z eq ¼ 0.2344c), where c is the lattice parameter, we have obtained R ¼ (R xx þ R yy þ R zz )/3 ¼ 74.4%, in perfect agreement with the experimental value of 74.3% (ref. 39). In Fig. 4a, we show the relative change in the reflectivity R xx for small displacements of the atoms in the z direction. For comparison, a value of Dz ¼ 18.4 pm would undo the Peierls distortion in Bi. The calculated reflectivity change shows a clear dependence on the displacement in the z direction with a fitted slope of q(DR/R)/ qz ¼ 0.0164 per pm. Moving the Bi nuclei in the x direction changes the reflectivity tensor anisotropically with a slope of ± 0.0046 per pm, as shown in Fig. 4b.
Utilizing these parameters, we have converted our measured beats DR(t probe )/R shown in Fig. 2c to the real-time atomic motions in the z and x directions. Details of the conversion process are described in the Methods section. Figure 2e-g shows those atomic motions in the Bi unit cell converted from the observed beats shown in Fig. 2c. In this conversion, the observed beats between t probe ¼ 0.82 and 10.48 ps are fitted with the oscillatory function in equation (1), giving the amplitude parameters, DR a /R and DR e /R, transformed to the absolute displacement of the Bi atom. Details of this fitting procedure and the relevant parameters are given in the Methods section. It is clearly seen from Fig. 2e-g that the ultrafast 2D atomic motions in crystalline bismuth have been actively controlled by shaping the pump pulse.
Comparison between CSP and pulse-train techniques. The train of femtosecond pulses has often been used in the previous experiments to control coherent phonons [3][4][5][6]8 . It should be interesting to understand the difference between those pulsetrain techniques and the current CSP technique from the viewpoint of THz intensity modulation. Figure 5a shows the intensity envelopes I(t) of the two superposed CSPs and the train of FT-limited pulses with their similar peak-to-peak intervals. Figure 5b shows the THz spectraĨ(o) Fourier-transformed from the envelopes I(t) shown in Fig. 5a. The clear difference is seen in the high-frequency region above 10 THz, where the second and third harmonics of the lowest frequency component around 6 THz are filtered out in the CSP case. The CSP technique thus allows for more flexible tuning of the relative intensities of different frequency components in the THz region than the conventional pulse-train approach, making CSPs ideal for controlling ultrafast 2D atomic motions in bulk solids. Spatial light modulators or acousto-optic modulators could be other candidates to produce temporal profiles similar to the one of the CSPs.

Discussion
The intensity of the E g mode in Fig. 2d should be correlated with the THz intensity at its resonance frequency, because the E g mode is excited directly by the pump pulse 1,2,15,16 . The coupling between the A 1g and E g phonons, which has been discussed in ref. 34, is of third order in the atomic displacements and can be safely neglected at the low fluence in the present experiments.
The A 1g phonons originate from a laser-induced shift of the potential energy surface minimum in the z direction in Fig. 1a (displacive excitation of coherent phonons (DECP) model) 11 , whereas the E g phonons originate from impulsive stimulated Raman scattering (ISRS model) 41 . When irradiated with a FTlimited laser pulse, the DECP mechanism launches coherent phonons that are cosine-like with an initial phase of zero degree, whereas the ISRS launches sine-like ones that are phase-shifted by 90°(ref. 11). The present initial phases of the A 1g and E g phonons are determined by a combination of their generation mechanisms (DECP and ISRS, respectively) and the intensity envelope I(t) of our modulated pump pulse. The effect of an intensity envelope on the initial phase of coherent phonons has been recently discussed by Shimada et al. 9 . The trajectory starting with these initial phases will eventually reach any point in the 2D configuration space accessible with a given ratio of A 1g and E g amplitudes as long as the phonon lifetime is infinite. However, the lifetime is finite, so that there could be points that are not reached by this trajectory. It is important to note that this trajectory can be modified to reach an arbitrary point in that accessible configuration space if the initial phases of those two phonon modes are actively controlled independently. This control could be implemented by tuning the chirp rates of the two CSPs independently with a dispersive optic inserted in each of the two optical paths of the interferometer. This additional control should be useful when the phonon lifetime is short, and the trajectory needs to be driven to a particular atomic configuration that is effective in inducing a specific dynamical process such as a phase transition. It should be also noted, however, that the relative phase of A 1g and E g modes recurs to its initial value every 1.1 ps, which is shorter than the lifetimes of those phonons by a factor of B10 in the present case. Each of the controlled trajectories shown in Fig. 2e-g thus covers almost all possible trajectories within the possible 2D space given by each relative amplitude of those two phonon modes. The current amplitude control, therefore, allows for almost full control of 2D atomic motions.
The present combination of the ultrafast optical measurements and the DFT calculation of optical response as a function of atomic displacements offers a new method to visualize ultrafast structural changes in solids, all-optical alternative to timeresolved X-ray and electron diffraction techniques 30,42 .
Our control-visualization scheme is based on the simple, robust and universal concept that in any physical system, 2D motion of a particle is determined by the projections of its pathway onto two orthogonal axes. The scheme is thus applicable to a variety of condensed matter systems.

Methods
Data processing. The reflectivity change obtained as the differential signal of the two balanced PDs was first amplified by a current amplifier and then transferred to a digital oscilloscope at each pump-probe delay t probe scanned with the fast scan mechanical stage. The displacement of this stage was proportional to the voltage applied, and the voltage was swept sinusoidally in each scan. This temporal evolution of the voltage was associated with electric noises and was fitted with a sine function to calibrate t probe . The data points thus sampled along the t probe axis were not equidistantly distributed. We interpolated these measured data points with cubic spline functions to give equidistant points to be fast-Fourier-transformed with the Hanning window function and zero padding. The temporal evolutions of the reflectivity change shown in this study are composed of these interpolated data.
The calibration of the reflectivity change signal was performed as follows. The phonon signal accumulated by the digital oscilloscope was transferred to the personal computer with a 16-bit integer format. To convert these integer-formatted data to the voltage unit, we independently performed a reference signal acquisition of a sine curve generated by a function generator as an input source. The range of the oscilloscope was set to the same scale as we have used in the phonon measurement, and data were saved in both the raw voltage format and 16-bit integer format. Each data was fitted with a sine function, and by comparing the fitted amplitude parameter of each data we obtain the proportionality parameter to transform the integer data to voltage unit. Using this parameter, the phonon signal in 16-bit integer format was rescaled to the voltage unit and then further transformed to the current unit according to the current amplifier setting. This value was divided by the recorded current value of each PD to give the reflectivity change data as plotted in the main text.
Conversion process from phonon signals to 2D atomic motions. The coherent phonon signals obtained with three different t mod timings as shown in Fig. 2c were fitted with the model function equation (1) to give the amplitude of each of the A 1g (longitudinal) and E g (transverse) phonon modes. Part of the phonon signal within the time window from t probe ¼ 0.82 to 10.48 ps was used for the fitting procedure.
The frequency and decay rate of each phonon mode was determined as follows. First, the green trace of Fig. 2c was fitted by equation (1), neglecting the contribution of E g mode. Next, using the parameters G a and n a fixed to the value determined by this process, we fitted the blue and red traces of Fig. 2c. The parameters G e and n e obtained from the blue and red traces were averaged to give the final value of G e and n e . Finally, we performed the fitting procedure again with these four parameters fixed to determine other parameters. Supplementary Table S1 summarizes the resulting fitting parameters for each phonon signal. The experimental phonon signals and the best-fitting results are shown in Supplementary Fig. S1. The amplitude parameters, DR a /R and DR e /R, were then converted to the displacements of the Bi atom in the z and x directions by using the coefficient q(DR/R)/qz ¼ 0.0164 per pm and q(DR/R)/qx ¼ 0.0046 per pm. These 2D displacements are plotted in Fig. 2e-g.