Nonlinear phononics in 2D SnTe: a ferroelectric material with phonon dynamical amplification of electric polarization

Ultrafast optical control of ferroelectricity using intense terahertz fields has attracted significant interest. Here we show that the nonlinear interactions between two optical phonons in SnTe, a two-dimensional in-plane ferroelectric material, enables a dynamical amplification of the electric polarization within subpicoseconds time domain. Our first principles time dependent simulations show that the infrared-active out-of-plane phonon mode, pumped to nonlinear regimes, spontaneously generates in-plane motions, leading to rectified oscillations in the in-plane electric polarization. We suggest that this dynamical control of ferroelectric material, by nonlinear phonon excitation, can be utilized to achieve ultrafast control of the photovoltaic or other non-linear optical responses.


INTRODUCTION
Light-matter interaction has become the main gateway through which the microscopic quantum materials properties are interacted, and hence engineered, with tools of laboratories in the macroscopic world. [1][2][3][4] Optical responses of non-centrosymmetric insulators are paradigmatic examples in this context. 5 They can host various photo-electric dynamical phenomena, such as shift current, second harmonic generation and circular photo-galvanic effect. [6][7][8][9] Moreover, various driven states induced by intense low-frequency external fields attracted substantial interest within various different context in recent years. 10,11 For example, the nonlinear Hall current in inversion-asymmetric materials shows a second order responses to low-frequency driving fields, controlled by asymmetric distribution of Berry curvature (more specifically, the Berry curvature dipole). 10 Defining characteristics of ferroelectricity is the existence of switchable bi-stable configurations, which share the same space group but with opposite electric polarization. 12,13 However, recent studies using optical pump, in terahertz or low infra-red band, have expanded the scope of ferroelectricity to the non-equilibrium regime of driven states. 3,4,14,15 In this field, low-dimensional materials or nano-structured systems are more attractive in view of few-atom scale local switching that can lead to even enhanced non-volatility of state. 3,4,[11][12][13][14][15][16][17] An exemplary system in this area is the monolayer SnTe: the in-plane ferroelectricity of the atom-thick layer has Curie temperature of 270K, three-times higher than that of the corresponding three-dimensional ferroelectric bulk structure (100K). 13 In this study, we investigate the interacting phonon dynamics in this ferroelectric SnTe monolayer. Through real-time ab initio molecular dynamics simulations, we find that the coherently pumped out-of-plane A u phonon mode, up to nonlinear regimes, spontaneously entails the in-plane E u,x vibrations. Remarkably, the entangled phonon dynamics induces the in-plane polarization to increase over time in a rectified manner. This non-linear phonon interaction is explained in terms of the modification of the potential energy surface by the phonon modes, and the microscopic contributions to the polarization changes are analyzed on the standpoint of the electronic part of the Born effective charge. We discuss that this dynamically amplified ferroelectricity makes the variations of the Berry curvature field more steep, which results in enhanced Berry curvature dipole and also the increased stability of the ferroelectric phase.

Atomic geometry and electronic structure of monolayer SnTe
Monolayer SnTe consists of alternating Sn and Te atoms in a puckered rectangular lattice, as shown in Fig. 1a. Previous experimental studies revealed that the in-plane ferroelectricity of monolayer SnTe originates from the staggered atomic displacement between Sn and Te, which otherwise comprises the cubic rock-salt structure. 13 In the ferroelectric configuration, the inversion symmetry is lifted along with the removal of the mirror plane perpendicular to the x-axis while the mirror plane along the y-axis (M y ) is well preserved, as depicted by the dashed line in Fig. 1a. This equilibrium structure (one of the energy minima shown in Fig. 1b) is to be characterized by the space group of P mn2 1 . 13,16,18 To visualize the double-well potential energy with respect to the in-plane polarization, as depicted in Fig. 1b, the atomic structures are linearly interpolated (or extrapolated) from the two equilibrium geometries at each equilibrium valley. The optimized geometry of monolayer SnTe, obtained by minimizing the total energies of density functional theory (DFT), is distorted by 70mÅ along the x-axis as shown in Fig. 1a, producing the in-plane polarization of P x = ±13.9 × 10 −12 Cm −1 . 16 The monolayer SnTe is a seminconductor with 0.6 eV indirect band gap, which is underestimated value by DFT calculation comparing with experimental gap (1.6eV), 13,19 between the valence band maximum at a point on the ΓX line, and conduction band minimum at a point on ΓY, as indicated by upward and downward arrows in Fig. 1c, respectively. This geometric nature of the monolayer SnTe is well reflected in the distribution of the Berry curvature which is asymmetric with respect to the same mirror plane (M y ), as indicated by the dashed line in Fig. 1d. Previous studies have suggested that this asymmetric Berry curvature can be attributed to the orbital Rashba effect. 16 Upon hole doping, this asymmetry in the Berry curvature grants sizable Berry curvature dipole (Λ y ) in the in the y-direction which can be manifested helicity-dependent second-order electro-optical responses. 10,11 Non-linear lattice dynamics and electric polarization We now show how the in-plane polarization is delicately intertwined with selected phonon dynamics. We first briefly describe the phonon structures of the monolayer SnTe, as depicted in Fig. 2a. 20 The optical phonons can be classified into three groups: the in-plane oscillations at around 1.2 to 1.6 THz, the out-of-plane modes at higher frequency range from 4.0 to 4.8 THz, and the intermediate bands around 3.0 THz that comes form the mixed motions. In the present work, we focus on the phonons that can be directly coupled with the electronic dipole oscillation. In this monolayer SnTe, the in-plane motions (E u,x and E u,y ) and the out-of-plane A u mode are active for infra-red field (IR-active). The A u phonon mode mainly consists of the alternating atomic motion between the Sn and Te in the out-of-plane direction with the period of 216 fs period as shown in Fig. 2a (upper). The E u,x phonon is characterized by the alternating motion between Sn and Te along the x-direction, has a period of 742 fs, as shown in Fig. 2a Figure 2b shows that the obtained z-direction motions in the BO dynamics indeed prove the A u mode vibration with the frequency of 216 fs. On the different instantaneous atomic configuration along this lattice dynamics, we calculated the electric polarization, of which the in-plane component is shown in Fig. 2c. Unexpectedly, instead of vibrating back and forth around the equilibrium point, the induced in-plane polarization oscillates in a rectified manner leading to increased polarizations at all time. Furthermore, this polarization oscillation complies with the frequency of the in-plane E u,x phonon mode (742 fs) rather than that of the driving of A u phonon mode (216 fs). This dynamical polarization is mainly contributed by the electronic part, while the change in the ionic contribution is negligible. Similarly, the atomic displacements in the x-direction, recorded in Fig. 2d, mainly follows the pattern of E u,x phonon: the major oscillation reveal the period of 742 fs, while the minor secondary oscillations with higher frequencies are deemed to reflect the beating between the two phonon modes. However, it is worth noting that this induced in-plane vibration deviates from the static equilibrium point. The oscillation centers for two Sn atoms shifted positively (indicated by the upward arrow in Fig. 2d), whereas that of two Te atoms are relocated negatively (indicated by the negative arrow in To access the non-linear interaction between the E u,x and A u phonon modes, we calculated the total energy of the lattice with the additional displacement along the E u,x mode on top of a given distortion along the A u mode: the potential energy surface along the additional E u,x coordinate hereafter is denoted as ∆E Q Au [Q Eu,x ], which is plotted in Fig. 3b. Notice that the potential energy shown in the E u,x coordinate is plotted with respect to the minimum energy of configurations with given Q Au . Irrespective of whether the movement in the A u mode is positive or negative (Q Au = +80mÅ and Q Au = −80mÅ in Fig. 3b), the energy minimum shifts positively in the E u,x coordinate (as indicated by the positive arrow in Fig. 3b). This can be explained by considering the Coulomb interaction between the asymmetric ionic chain in the primitive cell. At the equilibrium, the Sn and Te atoms alternates in the x-axis with small in-plane displacement, as indicated by x in the cross-sectional view, shown in  Fig. 3c), which provide repulsive and attractive forces between Sn and Te atoms, respectively. The repulsive restoring force on the compressed bond is depicted in Fig. 3c, and the attractive restoring force is expected from the elongated part. The stronger repulsive force produces the in-plane translations of Sn and Te atom from their equilibrium position, resulting in an increased in-plane displacement (Q Eu,x ). As depicted in the middle and bottom parts of We now quantitatively estimate the effect of this dynamical polarization enhancement on the stability of the ferroelectric phase. The Curie temperature for ferroelectric transitions in two-and three-dimensional structures, such as monolayer group-IV monochalcogenides and perovskite transition metal oxides, have been very accurately predicted from the double well potential energy surface by using forth-order Landau theory. 12,21,22 The free energy (F ) can be written in terms of spontaneous electric polarization (P ) with positive constants α and β : At the equilibrium polarization, the free energy minimization condition (dF/dP = 0) requires the Curie temperature to be T c = 2β α P 2 . By fitting the free energy form to the double well potential given in Fig. 1b, we obtained the constants 2β α = 0.967eV /(10 −10 Cm −1 ) 2 , and thus we have T c = 287K, which is well matched with experimental observation (∼ 270K). 13 Notice that, while thermal average phonons make the lattice vibrates with respect to the equilibrium, the coherently pumped A u phonons shifts the center of vibration, granting an Assuming an impulsive kick, the lattice was initially distorted into an A u mode eigenvector (Q init Au = 40mÅ), and the electronic states were prepared in its ground state at t = 0. The corresponding Fourier components are presented in Figs. 4c and 4d. On this initial impulsive kick into the A u phonon eigenvector, the out-of-plane oscillation (J z (t)) is dominated by the ω Au component with a sizable band width (Fig. 4c). On the other hand, the J x (ω) exhibits an apparent ω Eu,x peak, together with 2ω Au signal, as shown in Fig. 4d.
The apparent E u,x excitation directly evidences the nonlinear coupling of the E u,x mode to the A u mode, as explained above. The mechanism that underlies the frequency doubling (2ω Au ) can be ascribed to the parabolic dependence of the polarization on the A u mode displacement, as shown in the bottom panel of Fig. 3a: on the harmonic vibration in the A u mode (Q Au (t) ∼ sin(ω Au t)), the parabolic polarization should convey the second harmonic responses in the current from J x (t) = dP elec x /dt ∼ dQ 2 Au /dt ∼ sin(2ω Au t). This nonlinear phononic interaction can be detected from the radiation field by tracing the in-plane component of the polarity in-plane polarity, as schematically sketched in left panel of Fig. 4e.
In realistic experiment, the perfect impulsive kick is hard to realize, but a pulse with a substantial band width centered at ω Au can be applied. If the pulse is sufficiently short, which band width covers the range between ω Eu,x and ω Au as depicted in the Fig. 4e, the proposed phenomena can be experimentally observed.
The dynamical amplification of the polarization also has a substantial effect on the photovoltaic optical responses. Among various characteristics of non-centrosymmetric materials, such as second harmonic generation and circular photogalvanic current, which sharply depends on the polarization, here we examine the effect of the nonlinear phononics on the shift current. [5][6][7][8][9] The response function for shift current can be evaluated as follows; where f nm , R a nm , r b nm and ω nm are occupation difference, shift vector, Berry connection, and energy difference between n-th and m-th states, respectively. 6 The shift vector, which describes the charge center difference between occupied and unoccupied bands, is directly proportional to the internal polarization, 6,9,25 and thus will be most directly affected by this nonlinear phononics of the A u phonon. As an example, the shift current response function σ xxx (ω), the x-direction photoconductivity in response to the incoming light with x-direction polarization, is summarized in Fig. 4f with respect to various distortions into the A u mode.
The responses to a particular frequency of incoming light (ω = 0.6eV ) and the maximum of σ xxx over frequencies are presented in Fig. 4f. This enhanced shift current is attributed to the increase in the in-plane polarization, which is induced by the nonlinear phononics and thus irrespective of whether the A u phonon is positively or negatively distorted (see Fig. 2c).
Whether the bulk photovoltaic effect of solid can be altered by phonons has been questioned previously, 26 and our finding of the rectified oscillation of the polarization and the amplified photovoltaic responses can be considered as a new example in this light of search.

Interaction strength and the effect of hole doping
The larger initial pumping along the A u mode produces the stronger entailing in-plane motion in the E u,x mode. Here, we quantify the non-linear coupling strength between the two optical phonons and analyze the electronic origin of the coupling in terms of mode effective charge. The interaction between these two optical phonons can be described through a nonlinear harmonic oscillator coupling model, as treated in previous literature: 12,27 Here, M and g are effective mass and the coupling constant, respectively. The equations of motion for Q Eu,x (t) and Q Au (t) can be derived from the Hamiltonian dynamics, and their time series can be integrated through the Verlet algorithm. As we have assumed an impulsive kick in the ab initio dynamics study above, the time trajectories were evaluated from various initial displacements (Q init Au ) with zero initial velocity. By fitting the obtained Q Eu,x (t) to the BO molecular dynamics results, shown in Fig. 2d, we determined the coupling coefficients : g 0 = 1.7Å −1 and g 1 = 220Å −3 , when M = 1.12 × 10 6 m e . For example, the real-time profile obtained from the displacement of Q init Au = 40mÅ is presented in the inset of Fig. 5a. The amplitude for the induced vibration, denoted as Q ind Eu,x , is deduced from the time series Q Eu,x (t) and shown with respect to the various initial displacement of Q init Au , as shown in Fig. 5a. The lattice dynamics, and the ensuing polarization variation, can be described with various choices of coordinates, such as bond angles as used in the previous literature. 12 In the present work, in order to efficiently manifest the interaction between phonon modes, the displacements along the phonon modes are selected. The efficiency of this selection of coordinates can be observed from the molecular dynamics calculation results (Fig. 5a); the dynamical path obtained from the molecular dynamics is quite nicely described with the model Hamiltonian written in terms of these phonon coordinates. A comparison with Fig. 5a and the Fig. 3a explains that the nonlinear coupling is negligible for small Q Au , but increases almost parabolically as Q Au ≥ 20mÅ.
We now show that the coupling strength between the two phonons can be adjusted by an electrostatic gating: the hole doping. We evaluate the same Q ind Eu,x , as introduced above, by Q init Au = 80mÅ with the change of Fermi level, as shown in Fig. 5b. Overall, the increase in the hole concentration enhances the nonlinear phononic coupling strength. The electronic origin for this variation with hole-doping can be analyzed with the electric part of Born effective charge: Z β α = V cell dP elec α dR I,β . 28 Assuming the adiabatic evolution of the electronic states along a specific phonon mode Q, one can obtain the k-resolved Q mode effective charge ) 29 for the n-th band using a Kubo formula: 30-32  Figure 5c shows the x-component of the A u mode effective charge, which is sharply asymmetric against the sign change in Q Au . As a result, we can deduce the fact that, for a light hole doping (in the range −0.17 ≤ F ≤ 0 in Fig. 5c) Besides optical responses of this inversion-asymmetric insulators, as mentioned above, the carrier dynamics in the hole-doped material have been intriguing in terms of nonlinear Hall effect and second-order responses. 10,11,16 The Berry curvature dipole have attracted intense interests for last couple of years. 10,11 The coupling between an external electric field and the Berry curvature dipole of the system can result in a second-order response current of J = e 3 τ 2(1+iωτ )ẑ × E * (Λ · E), where Λ is Berry curvature dipole (Λ a = k f occ (k)∂ a Ω(k)), z is the normal unit vector perpendicular to the plane made up of D and E, and f occ , τ are occupation of Bloch state and constant. 10 As discussed previously, the Berry curvature dipole directly scales with the polarization that measures the amount of inversion-symmetry breaking of the lattice. 10,16 Here, we focus on the fact that the lattice distortion along the A u mode enhances the net in-plane polarization, and leads to an increase in the Berry curvature dipole. The static DFT calculation of the Berry curvature dipoles are presented in Fig. 5d, which confirms an almost monotonic increase of the Berry curvature dipole with the lattice displacement in A u mode eigenvector. With the enhanced coupling strength, as shown in Fig. 5b, the nonlinear phononic interaction is to be more apparently revealed by the nonlinear Hall current and second harmonic generations of the Fermi level carriers of the hole-doped system, which is mediated by the Berry curvature dipole. These results indicate that second-order optical responses, such as nonlinear Hall effect and shift current, can be enhanced by nonlinear phonon interactions.

DISCUSSION
In summary, we demonstrated that nonlinear couplings between optical phonons in monolayer SnTe, a two-dimensional in-plane ferroelectric system, produce a polarization oscillations of the in-plane polarization in response to the pumped out-of-plane phonon. We expect that this nonlinear phonon interaction can be generally observed from group-IV monochalco-

Computational details for first principle calculation
To investigate the electronic structure and effect of phonons of monolayer SnTe, we performed DFT calculation using the Quantum Espresso package. 34 To describe the electronelectron exchange correlation potential, Perdew-Burke-Ehrenof type generalized gradient approximation is employed. 35 The full-relativistic norm-conserving pseudopotential is considered to describe spin-orbit interaction. The plane wave basis set with 50 Ry energy cut-off is used to describe wavefunction under three-dimensional periodic boundary condition. The lattice parameters are employed from the experimental observation with vacuum slab up to 20Å. 13,16 The Bloch states are considered with k-point sampling up to 50 × 50 × 1, and convergence is checked up to 80 × 80 × 1 k-point mesh. The electric polarization of the 2D periodic system is evaluated following equation: is the electric polarization of the electrons, and Z i and R i are the ionic charge and the position of the i-th atom, respectively. 36 For the BO molecular dynamics simulation, Verlet algorithm is employed with time step dt = 0.48 fs. We also checked on that the total system energy is conserved within 0.1 meV during molecular dynamics time of 2.1 ps. Real-time propagation of Bloch state were performed by using our in-house computation packages. 23 Detailed techniques related to the time-integration of the time-dependent DFT and Kohn-Sham equations are well described in our previous works. 23,24 Data availability statement All relevant data is included in the main manuscript. All the data generated and analysed during this study are available from the corresponding authors upon reasonable request.