Proposal for midinfrared light--induced ferroelectricity in oxide paraelectrics
aa r X i v : . [ c ond - m a t . s t r- e l ] D ec Proposal for midinfrared light–induced ferroelectricity in oxide paraelectrics
Alaska Subedi
Centre de Physique Theorique, Ecole Polytechnique, CNRS,Universit´e Paris-Saclay, F-91128 Palaiseau, France andColl`ege de France, 11 place Marcelin Berthelot, 75005 Paris, France (Dated: November 12, 2018)I show that a nonequilibrium paraelectric to ferroelectric transition can be induced using mid-infrared pulses. This relies on a quartic lQ z Q x coupling between the lowest ( Q l z ) and highest( Q h x ) frequency infrared-active phonon modes of a paraelectric material. Density functional cal-culations show that the coupling constant l is negative, which causes a softening of the Q l z modewhen the Q h x mode is externally pumped. A rectification along the Q l z coordinate that stabilizesthe nonequilibrium ferroelectric state occurs only above a critical threshold for the electric field ofthe pump pulse, demonstrating that this is a nonperturbative phenomenon. A first principles calcu-lation of the coupling between light and the Q h x mode shows that ferroelectricity can be induced inthe representative case of strained KTaO by a midinfrared pulse with a peak electric field of 17 MVcm − and duration of 2 ps. Furthermore, other odd-order nonlinear couplings make it possible toarbitrarily switch off the light-induced ferroelectric state, making this technique feasible for all-opticdevices. PACS numbers: 77.80.Fm,78.20.Bh,63.20.Ry,78.47.J-
I. INTRODUCTION
Living organisms have used light to observe the prop-erties of materials since the evolutionary development ofcomplex eyes. However, ultrafast light-control of mate-rials properties has only become feasible after the con-struction of high-powered lasers in the previous century,and this field has flourished because light-induced pro-cesses have the potential to lead to new devices andphysical phenomena. Many examples of light-inducedphase transitions using near-visible sources have been ob-served, including discontinuous volume changes in poly-mer gels, quasiionic-to-quasineutral transition in organicmolecular compounds, low-spin to high-spin transitionin metal organic frameworks, insulator-to-metal transi-tion in perovskite manganites, and, remarkably, a tran-sition to a hidden metastable state in 1 T -TaS . All ofthese examples involve transition to higher-temperatureor metastable phases, and the quest to stabilize a phasewith less symmetry or more order using light remainselusive.More recently, intense midinfrared pulses have beenused to directly control the dynamical degrees of freedomof the crystal lattice. Such mode-selective vibrationalexcitations have been used to induce insulator-to-metaltransitions and melting of orbital and magnetic or-ders. Midinfrared excitations have so far not caused tran-sitions to more-ordered phases, although light-inducedsuperconductivity has been claimed in several cupratecompounds and K C . However, these claims relyon interpreting the two-dimensional response functionΣ( ω, τ ) measured in the pump-probe experiments asthe optical conductivity σ ( ω ) that is measured in time-domain spectroscopy. It is unclear whether such an inter-pretation is justified, especially at low frequencies, whenthe light-induced state is short-lived, as is the case in these experiments. In any case, a light-induced tran-sition to a lower-symmetry phase is not observed in anyof these experiments. Nonetheless, midinfrared excita-tion should be an effective tool to drive materials tobroken-symmetry phases because selective and coherentexcitation of the low-energy structural degrees of freedomshould cause minimal dissipation as heat.The mechanism for mode-selective light-control of ma-terials was proposed by F¨orst et al.
This involves ex-citing an infrared-active phonon mode Q IR of a materialusing an intense light pulse which then causes the lat-tice to displace along a fully symmetric A g Raman modecoordinate Q R due to a nonlinear coupling Q R Q be-tween the two modes. A quantitative microscopic theoryof this phenomenon was developed in Ref. 19, and calcu-lations based on this theory in combination with a time-resolved x-ray diffraction experiment was used to resolvethe midinfrared light–induced changes in the structure ofYBa Cu O . . In addition to the historically discussedcubic order Q R Q IR Q IR coupling, Subedi et al. haveshown that a sizeable quartic order Q Q coupling canoccur and studied the dynamics due to such a coupling. They found that such a quartic coupling exhibits variousdistinct regimes of dynamics, including transient modesoftening and dynamic stabilization in a rectified state.In contrast to the case of the cubic coupling, the dis-placement along the Q R coordinate occurs only above acritical pump amplitude threshold for the quartic Q Q coupling. A more recent work has reproduced several as-pects of the dynamics of this coupling. Unlike the cu-bic Q R Q IR Q IR coupling, the Q Q coupling cancause a rectification along a symmetry breaking mode, but such a light-induced symmetry breaking of a crystalstructure has so far not been observed.It has recently been predicted that ferroelectric polar-ization can be switched using midinfrared pulses. Inthis paper, I extend that work to the paraelectric phaseand show that ferroelectricity can also be induced in tran-sition metal oxide paraelectrics. This relies on a quartic lQ z Q x coupling and is a nonperturbative effect that oc-curs only above a critical pump amplitude. Here, Q h x is ahigh-frequency infrared-active phonon mode that shouldbe externally pumped and Q l z is the lowest frequencyinfrared-active mode that is transverse to the pumpedmode. I find that the sign of the coupling constant l isnegative in several transition metal oxides, which causesthe Q l z mode to soften as the Q h x mode is externallypumped. But other quartic order couplings t Q x Q h x , t Q x Q x , and t Q l x Q x between Q h x and the lowestfrequency mode Q l x that is longitudinal to the pumpedmode are larger in magnitude. In the cubic materials,the couplings between Q l x and Q h x modes are such thatthe Q l z mode may not develop a light-induced dynami-cal instability. However, I find that the couplings in thelongitudinal direction can be effectively reduced by ap-plying strain so that a light-induced ferroelectric state isstabilized by rectification along the Q l z coordinate.I illustrate this theory for the representative case ofKTaO to show that light-induced ferroelectricity can begenerated in the strained version of this material when apump pulse with an electric field of ∼
17 MV cm − andpulse duration of 2 ps is used. Interestingly, this valueis noticeably smaller than what is expected for the crit-ical pump amplitude due to a Q z Q x coupling. I findthat this reduction is due to the presence of substan-tial sixth order Q z Q x and Q z Q x couplings. Further-more, I show that the light-induced rectification can bearbitrarily suppressed by pumping the highest frequencyinfrared-active mode Q h x that is longitudinal to Q l x withanother weak pulse. Such a control is necessary for ap-plications in devices. In addition to KTaO , I find simi-lar nonlinear couplings in SrTO , and LaAlO , and thistechnique could be generally applied to many transitionmetal oxide paralectrics. II. APPROACHA. Computational details
The phonon frequencies and eigenvectors, nonlinearcouplings between different normal mode coordinates,and the coupling between light and pumped infraredmode were all obtained from first principles using densityfunctional calculations as implemented in the vasp soft-ware package. I used the projector augmented wave pseu-dopotentials provided with the package with the elec-tronic configurations 3 s p s (K), 5 p s d (Ta), and2 s p (O, normal cut-off). A plane-wave cut-off of 550eV for basis-set expansion, an 8 × × k -point grid forBrillouin zone sampling, and the PBEsol version of thegeneralized gradient approximation was used. The calculations were done using the relaxed lattice pa-rameters for the cubic structure. For the strained struc-ture, the c lattice parameter that minimized the total en- ergy for the given strain was used. A very small energyconvergence criteria of 10 − eV was used in the calcula-tions to ensure high numerical accuracy. After relaxingthe lattice parameters, I calculated the phonon frequen-cies and eigenvectors using the frozen phonon method asimplemented in the phonopy software package. Af-ter the normal mode coordinates were identified, totalenergy calculations were performed as a function of the Q l z , Q l x , and Q h x coordinates for values ranging between − √ amu with a step of 0.1 ˚A √ amu to obtainthe energy surfaces V ( Q l z , Q l x , Q h x ). These were thenfitted to polynomials given in Eq. A1 to obtain normalmode anharmonicities and nonlinear couplings betweenthe three coordinates. For the materials that I explored,polynomials with anharmonicities up to twentieth orderand nonlinearities up to eighth order were needed to en-sure accurate fit to the calculated energy surfaces. Sincethe polynomial fits the calculated energy surfaces almostexactly, there are no approximations in the calculationsof the nonlinear couplings, beyond that for the exchange-correlation functional.The Born effective charges were calculated using den-sity functional perturbation theory, and a larger 16 × × k -point grid was used in these calculations. Thecalculated Born effective charges and phonon mode eigen-vectors were used to calculate the mode effective charge Z ∗ m that determines the strength of the coupling of lightto the pumped phonon mode from first principles. Thecoupled equations of motion for the three coordinateswere numerically solved using the lsode subroutine ofthe octave software package. B. Identifying light-induced ferroelectricity
Phase transitions cannot occur at short timescales innonequilibrium conditions, and any light-induced ferro-electricity will disappear once the external light sourcevanishes. Therefore, it is necessary to establish an un-equivocal protocol for identifying light-induced ferroelec-tricity. Examining the intensity and phase of the secondharmonic generation of the transmitted probe pulse is aconvenient way to study ferroelectricity in pump-probeexperiments, and it will be necessary to distinguishbetween light-induced ferroelectricity and a long–time-period excitation that both generate second harmonicsif the probe pulse is shorter than the period of the low-frequency mode. For the purpose of this discussion, alight-induced ferroelectric state is deemed to have oc-curred both if the phase of the second harmonics doesnot change and the intensity of the second harmonicsshows at least two peaks over the full width at half max-imum (FWHM) duration of the pump pulse. Therefore,the pump pulse duration should in general be larger thanthe period of the equilibrium-condition lowest frequencymode to establish light-induced ferroelectricity. However,this is not a strict condition, and other well-defined crite-ria could also be specified. In particular, the lowest fre-
FIG. 1. (Color online) Displacement patterns of the (a) low-est frequency Q l z and (b) highest frequency Q h z modes ofthe cubic phase of KTaO . The x and y components of thesetriply degenerate modes can be obtained by appropriate rota-tion. The big, medium, and small spheres denote K, Ta, andO, respectively. quency oscillations could (and indeed does) occur with alarger frequency in the rectified state, and any method(such as time resolved x-ray diffraction) that can distin-guish oscillations about a displaced position can establishlight-induced ferroelectricity. III. RESULTS AND DISCUSSIONSA. Cubic KTaO The paraelectric phase of several AB O perovskite ox-ides occurs in the cubic structure. So it is natural toask if ferroelectricity can be induced in these cubic para-electrics by a midinfrared excitation of their infrared-active phonon modes. These materials have five atomsper unit cell, and they thus have four triply degenerateoptical phonon modes at the zone center. Factor groupanalysis shows that three of these modes have the irre-ducible representation T u , and these are infrared active.The remaining one has the irreducible representation T u and is optically inactive. Ferroelectricity is generally as-cribed to a dynamical instability of an infrared-activetransverse optic phonon mode. Indeed, most ferroelec-tric materials show a characteristic softening of an in-frared transverse optic mode as the transition tempera-ture is approached. Here I investigate if a similar soft-ening and instability of the lowest frequency T u modecan be achieved by an intense laser-induced excitationof the highest frequency T u mode in the representativecase of cubic KTaO .The calculated phonon frequencies of cubic KTaO using the relaxed PBEsol lattice parameter of 3.99˚A are Ω l = 85 cm − and Ω h = 533 cm − for thelowest and highest frequency T u modes, respectively.These are in good agreement with previously calcu-lated values. They also agree well with the frequenciesobtained from hyper-Raman scattering experiments at V ( , Q l x , Q h x ) ( e V ) Q l x (Å ‘‘‘‘(cid:214) amu) Q h x = 0.0Q h x = 0.5Q h x =−0.5 FIG. 2. (Color online) Total energy as a function of the longi-tudinal Q l x coordinate for several values of the Q h x coordinatefor cubic KTaO . room temperature. The atomic displacement patternsdue to these two modes are shown in Fig. 1. Withoutloss of generality, I consider the case where the x compo-nent of the highest frequency T u mode Q h x is pumpedby an intense light source and study how such an exci-tation changes the dynamics of the lowest frequency T u mode along the longitudinal Q l x and transverse Q l z co-ordinates. I ignore the dynamics along the second trans-verse coordinate Q l y as its dynamics will be qualitativelysimilar to that of the Q l z coordinate.
1. Dynamics of the lowest frequency longitudinal component
Fig. 2 shows several total energy V ( Q l z , Q l x , Q h x )curves along the projection Q l z = 0. The curves arenot symmetric upon reflection at Q l x = 0 and the + Q h x and − Q h x curves do not overlap. This indicates the pres-ence of coupling terms that have odd orders of Q l x and Q h x . A polynomial fit of the energy surface shows thatthe coupling terms t Q x Q h x , t Q x Q x and t Q l x Q x areall large relative to the harmonic term Ω of the low-est frequency mode (see Table I). The presence of thesecouplings is consistent with the symmetry requirements.Since the equilibrium structure has inversion symme-try and we are considering two odd modes along thesame direction, any term Q m l x Q n h x is allowed as long as m + n = even . The next allowed order of coupling is Q m l x Q n h x with m + n = 6. These are an order of magni-tude smaller than the m + n = 4 terms (see Table II inthe Appendix), but they are comparable in magnitude tothe harmonic term Ω .The nonlinear couplings between the Q l x and Q h x modes impart a force equal to − ∂V /∂Q l x along the Q l x coordinate. This force is finite and large whenthe Q h x mode is externally excited by an intense lightsource. The lowest order nonlinear terms of this force TABLE I. The coefficients of the harmonic and nonlinear cou-pling terms of cubic and strained KTaO . The units of a Q m Q n term are meV ˚A − ( m + n ) amu − ( m + n )2 . The sign of thecoupling is relevant only when the coordinates come with evenpowers.coefficient order cubic strainedΩ z Q z .
06 1 . x Q x .
06 55 . Q x .
77 1136 . t Q x Q h x − .
35 97 . t Q x Q x .
00 208 . t Q l x Q x − .
58 195 . l Q z Q x − . − . m Q z Q x − . − . m Q z Q x − . − . are − ∂V /∂Q l x = − t Q x Q h x − t Q l x Q x − t Q x . The − t Q x term acts as a nonresonant oscillating force tothe Q l x mode. The effect of the − t Q x Q h x term wouldaverage over the slow oscillation of the Q l x mode rela-tive to that of the Q h x mode. The − t Q l x Q x termaffects a time-dependent modulation of the frequency ofthe Q l x mode, and it does not cancel over the slow oscil-lation of the Q l x mode because Q x has a nonzero timeaverage. Unfortunately, the sign of t is positive, sothe frequency of the Q l x mode increases as the Q h x modeis pumped. A similar analysis of the next order Q m l x Q n h x terms with m + n = 6 also shows that the Q l x mode doesnot soften due to the effects of nonlinear coupling terms.
2. Dynamics of the lowest frequency transverse component
What about the dynamics of the transverse compo-nent Q l z of the lowest frequency mode? Fig. 3 showsseveral total energy V ( Q l z , Q l x , Q h x ) curves along theprojection Q l x = 0. One immediately notices that thecurves are symmetric upon reflection at Q l z = 0 andthat the − Q h x and Q h x curves overlap, showing thatonly even powers of both Q l z and Q h x appear in thenonlinear coupling terms. This is again consistent withthe symmetry requirements, which does not allow prod-ucts with odd powers of mutually perpendicular compo-nents Q l z and Q h x . The coefficients of the lowest or-der nonlinear terms lQ z Q x , m Q z Q x , and m Q z Q x are given in Table I. They are all at least twenty timessmaller than the magnitude of the quartic order cou-plings between the Q l x and Q h x coordinates. Neverthe-less, the sign of the coupling coefficients between the Q l z and Q h x modes are such that these terms softenthe frequency of the Q l z mode, as one sees by analyz-ing the forcing terms due to these nonlinear couplings − ∂V /∂Q l z = − lQ l z Q x − m Q z Q x − m Q l z Q x .Each term in the previous expression has even pow-ers of the Q h x coordinate, which ensures that their ef-fects are not averaged over the slow oscillation of the V ( Q l z , , Q h x ) ( e V ) Q l z (Å ‘‘‘‘(cid:214) amu)Q h x = 0.0Q h x = 1.0 Q h x = 2.5Q h x =−2.5 FIG. 3. (Color online) Total energy as a function of the trans-verse Q l z coordinate for several values of the Q h x coordinatefor cubic KTaO . Q l z mode. Furthermore, all these terms are propor-tional to odd powers of the Q l z coordinate, which causesthe frequency of the Q l z mode to change as Ω z → Ω z [1 + (2 lQ x ( t ) + 4 m Q z ( t ) Q x ( t ) + 2 m Q x ( t )) / Ω z ].Since the coupling constants are negative, this shouldlead to a softening of the transverse Q l z mode when the Q h x mode is externally pumped.The above discussion is not sufficient to convincinglyargue that the transverse Q l z mode will become dynami-cally unstable when the Q h x mode is externally pumped.There are two counteracting processes that may precludethis from happening. First, the coupling between the Q l x and Q h x modes is at least twenty times larger. As a re-sult, the Q l z component may receive a much smaller pro-portion of the external force due to nonlinear couplingwith the Q h x mode that is not sufficient to make thismode dynamically unstable when the latter is pumped.Moreover, I find that the coupling term pQ z Q x to bepositive and larger than the term lQ z Q x (see Table IIin the Appendix). When the Q l x component oscillateswith a large amplitude, this will provide an additive fac-tor that increases the frequency of the Q l z component.To settle this issue, I numerically solved the coupledequations of motion of the three coordinates Q l z , Q l x ,and Q h x , which are¨ Q h x + γ h ˙ Q h x + Ω Q h x = − ∂V nh ( Q l z , Q l x , Q h x ) ∂Q h x + F ( t )¨ Q l x + γ l ˙ Q l x + Ω Q l x = − ∂V nh ( Q l z , Q l x , Q h x ) ∂Q l x ¨ Q l z + γ l ˙ Q l z + Ω Q l z = − ∂V nh ( Q l z , Q l x , Q h x ) ∂Q l z . (1)Here, V nh ( Q l z , Q l x , Q h x ) is the nonharmonic part of thepolynomial that fits the calculated energy surface, andit includes both the anharmonicities of each coordinatesas well as the nonlinear couplings between these coor- (a) (b) − − − − − − − Q l z ( Å √ a m u ) time delay (ps) Q l z − − − − − − − − Q h x ( Å √ a m u ) E ( M V c m − ) time delay (ps) EQ h x FIG. 4. (Color online) Dynamics of the (a) Q l z and (b) Q h x coordinates of cubic KTaO after the Q h x coordinateis pumped by a pump pulse E with FWHM of 2 ps as shownin (b). dinates. The full expression for V nh ( Q l z , Q l x , Q h x ) isgiven in Eq. A1 in the Appendix. In addition to thenumerically large nonlinear couplings discussed above, itincludes anharmonicities up to the sixteenth order andnonlinear couplings up to the eighth order. γ h and γ l are the damping coefficients of the highest and lowestfrequency T u modes, respectively. They are taken tobe ten percent of the respective harmonic terms. F ( t ) = Z ∗ h x E sin(Ω t ) e − t / σ/ √ is the external force expe-rienced by the Q h x coordinate due to a light pulse of peakelectric field E . The calculated mode effective charge ofthe Q h x mode of cubic KTaO is Z ∗ h x = − . e amu − .A pump with a frequency of Ω = 1 . h and FWHM of σ = 2 . / Ω l .The result of the numerical integration of Eq. 1 in ahighly nonlinear regime is shown in Fig. 4. This wasobtained with a large peak electric field of E = 30 MVcm − that caused the pumped Q h x mode to oscillate witha maximum amplitude of 1.1 ˚A √ amu [Fig. 4(b)]. In thisregime, the Q l x mode oscillates about the equilibriumposition with a maximum amplitude of 0.25 ˚A √ amu (notshown). The frequency of its transverse counterpart Q l z does soften by around ∼ Q l z component, and its oscillationsabout the equilibrium position are damped even duringthe duration of the pump pulse [Fig. 4(a)]. I performedsimilar calculations for pump fields up to 100 MV cm − but was not able to find any instances where the Q l z mode becomes dynamically unstable.These calculations show that a dynamical instability ofthe lowest frequency infrared mode of cubic KTaO can-not be achieved by a midinfrared excitation of its highestfrequency infrared mode. However, this does not allow usto infer that light-induced dynamical instability cannotoccur in any cubic paraelectric. Indeed, if I artificiallyincrease the coefficient of the lQ z Q x term by six times,I am able to obtain a solution where the Q l z coordinateoscillates about a displaced position during the durationof the pump pulse. Such a large coupling between thepumped high-frequency mode and the transverse compo-nent of the low-frequency mode may exist in some mate- rials. B. Strained KTaO If the coupling between the externally pumped high-est frequency T u mode and the component of the lowestfrequency T u mode longitudinal to the pumped modecould be weakened in cubic KTaO , the transverse com-ponent of the lowest frequency mode would develop alight-induced dynamical instability. An effective way ofachieving this is by raising the frequency of the longitu-dinal component relative to that of the transverse com-ponent. This can be accomplished by applying a biaxialstrain on KTaO via an epitaxial growth on an appropri-ate substrate.I performed calculations on KTaO with 0.6% com-pressive biaxial strain. This can be achieved, for example,by growing KTaO on a GdScO substrate. The calcu-lated PBEsol lattice parameters of thus strained KTaO are a = b = 3 .
965 and c = 4 . T u mode of the cubic phasesplits into a nondegenerate A u mode and a doubly de-generate E u mode. The A u phonons involve atomic mo-tions along the z axis while the atoms move in the xy plane for the E u phonons. The calculated values for thelowest frequency A u and E u modes are Ω l z = 20 andΩ l x = Ω l y = 122 cm − , respectively. The highest fre-quency E u phonon that should be externally pumped hasa frequency of Ω h = 556 cm − .To find out whether the lowest frequency A u mode Q l z of strained KTaO develops a dynamical instabil-ity when the x component of the highest frequency E u mode Q h x is intensely excited by a light source, I againstarted my investigation by calculating the total energysurface V ( Q l z , Q l x , Q h x ) as a function of the three coordi-nates using density functional calculations. The nonlin-ear couplings between the Q l z , Q l x , and Q h x coordinatesof strained KTaO have the same symmetry requirementsas discussed for the cubic case, and a fit of a general poly-nomial to the calculated first-principles V ( Q l z , Q l x , Q h x )shows that same orders of nonlinearities are present inboth the cases. As a comparison of the numbers pre-sented in Table I shows (see also Table II in the Ap-pendix), the nonlinear couplings in the two cases do notdiffer by a large amount. The crucial difference betweenthe two cases is that the frequencies of the Q l z and Q l x coordinates are different in the strained case (Ω l z = 20and Ω l x = 122 cm − ), whereas they are equal in the cu-bic case (Ω l z = Ω l x = 85 cm − ). This has a profoundeffect in the dynamics of the Q l z coordinate because theforces experienced by a coordinate due to the nonlinearcouplings are weighted by the square of the frequency ofthe coordinate. We can see that ∂V∂Q l z is likely to bemuch larger than ∂V∂Q l x or ∂V∂Q l x . In simple words,the Q l z coordinate of strained KTaO gets much largerproportion of the force than the Q l z coordinate of cubicKTaO because the frequency of the Q l z mode is much − − − − Q l z ( Å √ a m u ) time delay (ps) Q l z − − − − Q l z ( Å √ a m u ) time delay (ps) Q l z − − − − − − Q l z ( Å √ a m u ) time delay (ps) Q l z (a) (b)(c) (d) E = 17E = 60 E = 100E = 10 − − − − − − − − Q h x ( Å √ a m u ) E ( M V c m − ) time delay (ps) EQ l z FIG. 5. (Color online) Dynamics of the Q l z coordinate ofstrained KTaO after the Q h x coordinate is pumped by apump pulse E with FWHM of 2 ps. The dynamics for fourdifferent values of the peak electric field E (MV cm − ) of thepump pulse are shown. smaller in the strained structure compared to the cubicstructure. Is this change big enough to result in a light-induced dynamical instability of the Q l z mode in strainedKTaO ?I again solved the coupled equations of motion of thethree coordinates Q h x , Q l x , and Q l z as given by Eq. 1 forthe case of strained KTaO . This time I used the poten-tial V nh ( Q l z , Q l x , Q h x ) obtained for strained KTaO fromfirst principles. The polynomial expression used in thecalculations and the numerical values of the coefficientsfor all the terms in the polynomial that fit the calculatedenergy surface are given in the Appendix. A pump pulsewith an FWHM of 2 ps ( > / Ω l z ) and frequency 1 . h is again used to excite the Q h x mode. The mode effec-tive charge of the Q h x mode in the strained structure is Z ∗ h x = − . e amu − .Fig. 5 shows the results of the numerical integration ofthese equations for three different regimes of dynamicsof the Q l z coordinate. At relatively small peak electricfields of the pump ( E < − ), the Q l z mode os-cillates about the equilibrium position with its harmonicfrequency (not shown). As the peak electric field is in-creased, the frequency of the Q l z mode decreases dur-ing the duration that the Q h x mode is being pumped[Fig. 5(a)]. As discussed above, this is due to the negativevalues of the coefficients of the coupling terms lQ z Q x , m Q z Q x , and m Q z Q x that cause a light-inducedsoftening of the Q l z coordinate. Since the duration ofthe pump pulse is finite, I naturally do not observe theperiod of the Q l z mode diverge. Instead, beyond a criticalvalue of the peak electric field of the pump ( E c ≃
17 MVcm − for σ = 2 ps), the Q l z coordinate oscillates abouta displaced position and has a non-zero value while the Q h x mode is being pumped [Figs. 5(b, c)]. In this recti-fied regime, the average potential felt by the Q h z mode has a double-well structure, and this mode is oscillatingabout one of the minima. The displacement along the Q l z coordinate is also amplified strongly in this regime.Since the Q l z mode is infrared active, this implies thatthe material is in a broken symmetry state with a finitedipole moment while the Q h x mode is externally pumped.The frequency of the Q l z oscillations in the rectifiedstate increases as the peak electric field of the pump isincreased beyond the critical threshold. This is evidentfrom a comparison of Figs. 5(b) and (c), which showsthat the frequency of the Q l z mode doubles as the peakelectric field E is increased from 17 to 60 MV cm − .This increase occurs because the double-well potentialfor the Q l z coordinate becomes deeper as the amplitudeof the Q h x oscillations increases.When the peak electric field is increased further ( E >
75 MV cm − ), the Q l z mode oscillates with a large ampli-tude and high frequency about the equilibrium position[Fig. 5(d)]. In this regime, the kinetic energy impartedto the Q l z mode is larger than the depth of the doublewells. As a result, the oscillation of the Q l z mode stopsbeing confined to one of the double wells, and the rectifiedbehavior along the Q l z coordinate is no longer observed.Even though the light-induced broken-symmetry phase isstabilized only for a range of values of the peak electricfield of the pump, this range 17 < E <
75 MV cm − isboth wide and approachable enough to make the light-induced ferroelectric state experimentally accessible.The existence of a critical threshold above which the Q l z coordinate is rectified and the presence of three dif-ferent regimes for the dynamics of this coordinate is con-sistent with the analysis of a Q Q nonlinear couplingbetween two different normal mode coordinates as pre-sented in Ref. 19. These features should be present in theexperiments to confirm the predictions made in this work.The critical pump amplitude depends on the frequenciesof the Q l z and Q h x modes and the coupling coefficient,as well as the pump pulse length and the initial condition(i.e. the r.m.s. displacement of the Q l z mode at a partic-ular temperature). For a pump pulse with FWHM of 2ps, I find that the Q l z mode starts to get rectified whenthe peak electric field is E = 17 MV cm − . With thispump pulse, the Q h x mode is oscillating with an ampli-tude of 0.9 ˚A √ amu, which corresponds to a maximumchange in the Ta-apical O bond length of 0.2 ˚A (thatis, 10%). The Q h x mode oscillates with an amplitude of1.3 ˚A √ amu when the peak electric field is E = 75 MVcm − . This is a modest increase in the energy of the Q h x mode due to the pump, and it indicates that a signifi-cant fraction of the pumped energy goes to maintainingthe rectified state along the Q l z coordinate. However, thelight-induced ferroelectric displacement along the Q l z co-ordinate is quite small because of the small magnitude ofthe couplings between Q l z and Q h x modes. The aver-age displacement along Q l z is ∼ ∼ √ amu for E = 17 and 75 MV cm − , respectively, which results inthe change of Ta-apical O distance by 0.015–0.030 ˚A.Curiously, the critical pump threshold obtained for (a) (b) − − − − − − − − Q h x ( Å √ a m u ) E ( M V c m − ) time delay (ps) E x E z Q l z − − − − − − − − Q h x ( Å √ a m u ) E ( M V c m − ) time delay (ps) E x E z Q l z FIG. 6. (Color online) Dynamics of the Q l z coordinate ofstrained KTaO for the delays of (a) 0.5 and (b) 0.0 ps be-tween the E h x and E h z pulses that pump Q h x and Q h z modes,respectively. strained KTaO is noticeably smaller than what is ex-pected for a Q z Q x coupling. In the total energy calcu-lations, the Q l z mode starts developing instability when Q h x is above 0.7 ˚A √ amu. So the critical Q h x amplitudeshould be 0 . √ . √ amu. Instead, I find that the Q l z mode becomes unstable when the Q h x amplitude is0 . √ amu. This reduction in the critical threshold isdue to the presence of a large and negative sixth ordercoupling terms m Q z Q x and m Q z Q x . Both theseterms give a subtractive contribution to the effective,light-induced frequency of the Q l z mode, which hastensits instability as a function of the pump intensity. C. Abruptly halting light-induced ferroelectricity
For light-induced ferroelectricity to be useful in appli-cations, it is necessary to be able to control the light-induced phase at will in an all-optical setup. In thiscontext, this means having the capability to switch offthe rectification of the Q l z mode while the Q h x modeis being pumped. The quartic order odd Q z Q h z and Q l z Q z couplings in the longitudinal direction can beused to our advantage for this purpose. To investigatethis possibility, I consider an experiment where an over-lapping pulse polarized along Q h z comes at an arbitrarydelay with respect to the rectification-causing pulse thatpumps the Q h x mode. I study the resulting dynamicsalong the Q l z coordinate by solving the coupled equa-tions of motion for the four coordinates ( Q l z , Q l x , Q h x , and Q h z ). The equations of motions are obtained fromthe potential V nh ( Q l z , Q l x , Q h x ) + V nh ( Q l z , Q h z ). Forcomputational efficiency, I do not consider the full po-tential V nh ( Q l z , Q l x , Q h x , Q h x ) spanned by the four co-ordinates.The results for the delays of 0.5 and 0.0 ps betweenthe pump pulses E h x and E h z that excite the Q h x and Q h z coordinates, respectively, are shown in Fig. 6. Thepeak electric fields of E h x and E h z are 60 and 3 MVcm − , respectively, and their FWHM is 2 ps. The pumpfrequencies are 1.01 times the respective phonon frequen-cies and the mode effective charge of Q h z is Z ∗ h z = − . e amu − . Note that excitation by only E h x causes rectifi-cation of the Q l z mode during the FWHM of the pulse [Fig. 5(c)]. However, an overlapping excitation by an-other weak pulse E h z immediately suppresses the light-induced rectification of the Q l z mode. Even a weak lon-gitudinal pump is efficient in halting the rectification be-cause the quartic order odd couplings between Q l z and Q h z modes are much larger than the couplings in thetransverse direction. IV. SUMMARY AND CONCLUSIONS
In summary, I have shown that midinfrared pulsescan be used to stabilize nonequilibrium ferroelectricityin strained KTaO , which is paraelectric at equilibriumconditions. This phenomenon relies on a quartic lQ z Q x coupling between the highest frequency infrared-activephonon mode Q h x and the lowest frequency infared-active mode Q l z that is transverse to Q h x . Density func-tional calculations show that the coupling constant l isnegative, which causes the Q l z mode to soften when the Q h x mode is externally pumped. The rectification alongthe Q l z coordinate occurs only above a critical electricfield of the pump pulse, demonstrating that this light-induced symmetry breaking is a unique nonperturbativeeffect. Such a threshold behavior should be observed inexperiments to corroborate the predictions made in thispaper. Additionally, the Q z Q x and Q z Q x couplingsare large, and this makes the rectified regime more ac-cessible. A first principles calculation of the couplingbetween light and the Q h x mode shows that ferroelectric-ity can be induced in strained KTaO by a midinfraredpulse with a peak electric field of 17 MV cm − and aduration of 2 ps. Furthermore, large odd quartic cou-plings Q z Q h z and Q l z Q z between Q l z and the highestfrequency infrared-active mode Q h z longitudinal to Q l z makes it possible to arbitrarily switch off the inducedferroelectricity by pumping the Q h z mode with anotherweak pulse. I find that similar nonlinear interactions ex-ist in SrTiO and LaAlO , and this technique could begenerally applied to other transition metal oxide para-electrics.At a more basic level, I have shown that materials canexhibit various nonlinear interactions between differentdynamical degrees of freedom that have hitherto beenoverlooked. These interactions enable us to induce andcontrol broken-symmetry phases using light, whose oscil-lating electric and magnetic fields average to zero by def-inition. Furthermore, I have demonstrated that the non-linear interactions can be effectively modified by applyingstrain. This motivates experiments that combine the dis-parate fields of nonlinear optics and heterostructuring toachieve materials control in an interesting manner. Ina broader perspective, these nonlinear interactions mayalso be present in other classes of systems, and they mightallow us to influence the dynamics of these systems in anunusual way. ACKNOWLEDGMENTS
I am grateful to Indranil Paul for helpful discus-sions. This work was supported by the European Re-search Council grants ERC-319286 QMAC and ERC-61719 CORRELMAT and the Swiss National Supercom-puting Center (CSCS) under project s575.
Appendix A: Expressions for total energy surfaces
Only low order nonlinear couplings that are relativelylarge were discussed in the main text. However, if onlylow order couplings and anharmonicities are considered,the fit to the calculated total-energy surfaces are not sat-isfactory. The dynamics of the coordinates with andwithout using the high order couplings also show largedifferences, especially at the nonlinear regime. Since theuse of the full polynomial expression in the solutions ofthe equations of motion are not computationally demand-ing, all the numerical results discussed in this paper wereobtained using the full expression given below.For cubic KTaO , the following polynomial V ( Q l z , Q l x , Q h x ) accurately fits the calculated totalenergy surface spanned by the three coordinates forvalues between − . . √ amu. V = 12 Ω Q z + a Q z + a Q z + a Q z + a Q z + a Q z + a Q z + a Q z + 12 Ω Q x + a Q x + a Q x + a Q x + a Q x + a Q x + a Q x + a Q x + 12 Ω Q x + c Q x + c Q x + c Q x + c Q x + c Q x + lQ z Q x + m Q z Q x + m Q z Q x + n Q z Q x + n Q z Q x + n Q z Q x + t Q x Q h x + t Q x Q x + t Q l x Q x + u Q x Q h x + u Q x Q x + u Q x Q x + u Q x Q x + u Q l x Q x + pQ z Q x + q Q z Q x + q Q z Q x + r Q z Q x + r Q z Q x + r Q z Q x + dQ z Q l x Q h x + e Q z Q x Q h x + e Q z Q x Q x + e Q z Q l x Q x + f Q z Q x Q h x + f Q z Q x Q x + f Q z Q x Q x + f Q z Q x Q x + f Q z Q l x Q x + gQ z Q l x Q h x + h Q z Q x Q h x + h Q z Q x Q x + h Q z Q l x Q x . (A1)For strained KTaO , the potential has additional a Q z and a Q z terms. Also, the coefficients of the Q n l z and Q n l x terms are different in the strained case. The coeffi-cients of the Q n l x terms for the strained case are denotedby b n in Table II.The polynomial V ( Q l z , Q h z ) that fits the energy sur-face spanned by the Q l z and Q h z coordinates is given by V = 12 Ω z Q z + a Q z + a Q z + a Q z + a Q z + a Q z + a Q z + a Q z + a Q z + a Q z + 12 Ω z Q z + d Q z + d Q z + d Q z + d Q z + d Q z + v Q z Q h z + v Q z Q z + v Q l z Q z + w Q z Q h z + w Q z Q z + w Q z Q z + w Q z Q z + w Q l z Q z . (A2)The nonharmonic potential V nh defined in the maintext is V without the harmonic Ω Q terms. The valuesof all the coefficients in Eqs. A1 and A2 obtained from afit to the calculated energy surfaces of cubic and strainedKTaO are given in Table II. I note that values lowerthan the magnitude of 10 − are below the accuracy ofthe density functional calculations. They are kept so thatthat the highest order anharmonicity has a positive sign,which keeps the numerical solution of the equation ofmotions stable. Appendix B: Mode effective charges
The mode effective charge vector Z ∗ m,α = ∂F m,α /∂E α relates the force F m,α experienced by the normal modecoordinate Q m due to an electric field E α along the direc-tion α . It is related to the Born effective charges Z ∗ κ,αβ of atoms κ in the unit cell of a material by Z ∗ m,α = X κ,β Z ∗ κ,αβ U m ( κ, β ) , where U m ( κ, β ) is the q = 0 eigendisplacement vectornormalized as X κ,β M κ [ U m ( κ, β )] ∗ U n ( κ, β ) = δ mn . Here M κ is the mass of the atom κ . The eigendisplace-ment vector is related to the eigenvector w m ( κ, β ) of thedynamical matrix by U m ( κ, β ) = w m ( κ, β ) √ M κ . Note that this definition of the mode effective charge isslightly different from the one used in Ref. 31. Here, Z ∗ m,α is related to the change in the value of the normalmode coordinate rather than the change in the atomicdisplacements due to a motion along the normal modecoordinate. This gives a different normalization factorfor Z ∗ m,α , and this quantity is expressed in the units of e amu − . Its sign is arbitrary because the eigenvectorof the dynamical matrix is defined up to a multiplicativeconstant. TABLE II. The coefficients of the harmonic, anharmonic and nonlinear coupling terms of cubic and strained KTaO . The unitsof a Q m Q n Q p term are meV ˚A − ( m + n + p ) amu − ( m + n + p )2 . The sign of the coupling is relevant only when the coordinates comewith even powers.coefficient order cubic strained coefficient order strainedΩ z Q z .
06 1 .
39 Ω z Q z . x Q x .
06 55 . b Q x . Q x .
77 1136 . b Q x − . a Q z .
55 51 . b Q x . a Q z − . − . b Q x − . × − a Q z .
47 2 . b Q x . × − a Q z − . × − − . × − b Q x − . × − a Q z . × − . × − b Q x . × − a Q z − . × − − . × − d Q z . a Q z . × − . × − d Q z − . × − a Q z − . × − d Q z . × − a Q z . × − d Q z − . × − c Q x .
17 78 . d Q z . × − c Q x − . × − − . v Q z Q h z . c Q x . × − . v Q z Q z . c Q x − . × − − . × − v Q l z Q z . c Q x . × − . × − w Q z Q h z . l Q z Q x − . − . w Q z Q z . m Q z Q x − . − . w Q z Q z . m Q z Q x − . − . w Q z Q z . n Q z Q x . × − . × − w Q l z Q z . n Q z Q x . × − . n Q z Q x − . × − − . × − t Q x Q h x − .
35 97 . t Q x Q x .
00 208 . t Q l x Q x − .
58 195 . u Q x Q h x − .
72 1 . u Q x Q x .
64 6 . u Q x Q x − .
81 18 . u Q x Q x .
38 24 . u Q l x Q x − .
70 15 . p Q z Q x .
29 6 . q Q z Q x − . − . q Q z Q x − . − . r Q z Q x . × − . × − r Q z Q x . × − . × − r Q z Q x . × − . × − d Q z Q l x Q h x − .
09 15 . e Q z Q x Q h x . − . e Q z Q x Q x − . − . e Q z Q l x Q x . − . f Q z Q x Q h x . × − − . × − f Q z Q x Q x − . × − − . × − f Q z Q x Q x . − . f Q z Q x Q x − . − . f Q z Q l x Q x . × − − . × − g Q z Q l x Q h x . − . × − h Q z Q x Q h x − . × − . × − h Q z Q x Q x . × − . × − h Q z Q l x Q x − . × − . × − ǫ αβ (Ω) = ǫ ∞ αβ + 4 πV X m Z ∗ m,α Z ∗ m,β Ω m − Ω , where V is the unit cell volume and Ω m is the frequencyof the mode m . This expression shows that the oscillatorstrength measured in optical spectroscopy is the squareof the mode effective charge. A. Suzuki and T. Tanaka, Nature , 345 (1990). S. Koshihara, Y. Tokura, T. Mitani, G. Saito, and T. Koda,Phys. Rev. B , 6853(R) (1990). S. Decurtins, P. Gutlich, C. P. Kohler, H. Spiering, and A.Hauser, Chem. Phys. Lett. , 1 (1984). K. Miyano, T. Tanaka, Y. Tomioka, and Y. Tokura, Phys.Rev. Lett. , 4257 (1997). L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D.Svetin, S. Brazovskii, and D. Mihailovic, Science , 177(2014). M. Rini, R. Tobey, N. Dean, J. Itatani, Y. Tomioka, Y.Tokura, R. W. Schoenlein, and A. Cavalleri, Nature ,72 (2007). A. D. Caviglia, R. Scherwitzl, P. Popovich, W. Hu, H.Bromberger, R. Singla, M. Mitrano, M. C. Hoffmann, S.Kaiser, P. Zubko et al. , Phys. Rev. Lett. , 136801(2012). R. I. Tobey, D. Prabhakaran, A. T. Boothroyd, and A.Cavalleri, Phys. Rev. Lett. , 197404 (2008). M. F¨orst, R. I. Tobey, S. Wall, H. Broberger, V. Khanna,A. L. Cavalieri, Y.-D. Chuang, W. S. Lee, R. Moore, W.F. Schlotter et al. , Phys. Rev. B , 241104(R) (2011). M. F¨orst, A. D. Caviglia, R. Scherwitzl, R. Mankowsky, P.Zubko, V. Khanna, H. Bromberger, S. B. Wilkins, Y.-D.Chuang, W. S. Lee et al.
Nat. Mater. , 883 (2015). D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M.C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A.Cavalleri, Science , 189 (2011). S. Kaiser, C. R. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Y.Liu, M. Le Tacon, T. Loew, D. Haug, B. Keimer, and A.Cavalleri, Phys. Rev. B , 184516 (2014). W. Hu, S. Kaiser, D. Nicoletti, C. R. Hunt, I. Gierz, M.C. Hoffmann, M. Le Tacon, T. Loew, B. Keimer, and A.Cavalleri, Nat. Mater. , 705 (2014). M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A.Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Ricc, S.R. Clark, D. Jaksch, and A. Cavalleri, Nature , 461(2016). J. T. Kindt and C. A. Schmuttenmaer, J. Chem. Phys. , 8589 (1999). J. Orenstein and J. S. Dodge, Phys. Rev. B , 134507(2015). M. F¨orst, C. Manzoni, S. Kaiser, Y. Tomioka, Y. Tokura,R. Merlin, and A. Cavalleri, Nat. Phys. , 854 (2011). M. F¨orst, R. Mankowsky, H. Bromberger, D. M. Fritz, H. Lemke, D. Zhu, M. Chollet, Y. Tomioka, Y. Tokura, R.Merlin, J. P. Hill, S. L. Johnson, and A. Cavalleri, SolidState Commun. , 24 (2013). A. Subedi, A. Cavalleri, and A. Georges, Phys. Rev. B ,220301(R) (2014). R. Mankowsky, A. Subedi, M. F¨orst, S. O. Mariager, M.Chollet, H. T. Lemke, J. S. Robinson, J. M. Glownia, M.P. Minitti, A. Frano et al. , Nature , 71 (2014). R. F. Wallis and A. A. Maradudin, Phys. Rev. B , 2063(1971). T. P. Martin and L. Genzel, Phys. Status Solidi B , 493(1974). M. Fechner and N. A. Spaldin, Phys. Rev. B , 134307(2016). D. M. Juraschek, M. Fechner, and N. A. Spaldin, cond-mat, arXiv:1607.01653 (2016). Ref. 24 claims that a Q R Q IR Q IR (with 1 = 2) couplingcauses a rectification along Q R . However, such a couplingactually leads to an oscillation of the Q R mode with anenvelop frequency Ω IR − Ω IR . A. Subedi, Phys. Rev. B , 214303 (2015). J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov,G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke,Phys. Rev. Lett. , 136406 (2008). K. Parlinski, Z-.Q. Li, and Y. Kawazoe, Phys. Rev. Lett. , 4063 (1997). A. Togo, F. Oba, and I. Tanaka, Phys. Rev. B , 134106(2008). S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi,Rev. Mod. Phys. , 515 (2001). X. Gonze and C. Lee, Phys. Rev. B , 10355 (1997). J. W. Eaton, D. Bateman, S. Hauberg, and R. We-hbring,
GNU Octave version 4.0.0 manual: a high-level interactive language for numerical computations K. Takahashi, N. Kida, and M. Tonouchi, Phys. Rev. Lett. , 117402 (2006). D. Talbayev, S. Lee, S.-W. Cheong, and A. J. Taylor, Appl.Phys. Lett. , 212906 (2008). J. F. Scott, Rev. Mod. Phys. , 83 (1974). D. J. Singh, Phys. Rev. B , 176 (1996). H. Vogt and H. Uwe, Phys. Rev. B29