S Omid Sayedaghaee

Charles Paillard

Sergey Prosandeev

Sergei Prokhorenko

Yousra Nahas

Bin Xu

L Bellaiche

Electric-field-control of electromagnons' frequency in multiferroics

Electromagnons, which are coupled polar and magnetic excitations in magnetoelectric materials, are of large interest for electronic and computing technological devices. Using Molecular Dynamics simulations based on an ab-initio effective Hamiltonian, we predict that the frequency of several electromagnons can be tuned by the application of electric fields in the model multiferroic BiFeO 3 , with this frequency either increasing or decreasing depending on the selected electromagnon. In particular, we show that the frequency of electromagnons localized at ferroelectric domain walls can be tuned over a 200 GHz range by realistic dc electric fields. We interpret the realized frequency increase (respectively, decrease) by local hardening (respectively, softening) of the associated polar phonons which couples to the applied electric field. The increase versus decrease of the electromagnons' frequency is further found to be correlated with the real-space localization of such phonons.

I. INTRODUCTION

Electromagnons, a coupled oscillation wave of electrical and magnetic dipoles in magnetoelectric materials, have bolstered large interest since they were first discussed in 1970s [START_REF] Baryakhtar | Quantum Theory of Oscillations In A Ferroelectric Ferromagnet[END_REF]2].

They have remained elusive for a long time, except in a few materials such as rare-earth (R) manganites RMnO 3 [START_REF] Pimenov | Coupling of phonons and electromagnons in GdMnO3[END_REF][START_REF] Pimenov | Possible evidence for electromagnons in multiferroic manganites[END_REF][START_REF] Rovillain | Magnetic Field Induced Dehybridization of the Electromagnons in Multiferroic TbMnO3[END_REF], RMn 2 O 5 [START_REF] Sushkov | Electromagnons in Multiferroic YMn2O5 and TbMn2O5[END_REF][START_REF] Sushkov | Electromagnons in multiferroic RMn2O5 compounds and their microscopic origin[END_REF] or ferrite perovskite oxides [START_REF] Stanislavchuk | electromagnon excitations in multiferroic DyFeO3[END_REF][START_REF] Cazayous | Possible Observation of Cycloidal Electromagnons in BiFeO3[END_REF][START_REF] Rovillain | Polar phonons and spin excitations coupling in multiferroic BiFeO3 crystals[END_REF][START_REF] Skiadopoulou | Spin and lattice excitations of a BiFeO3 thin film and ceramics[END_REF]. Among them, BiFeO 3 , a room temperature multiferroic [START_REF] Lebeugle | Room-temperature coexistence of large electric polarization and magnetic order in BiFeO3 single crystals[END_REF], has been one of the most intensely considered magnetoelectrics for electromagnon detection, characterization and technological applications. Among the various experimental and theoretical studies of electromagnons, most have focused on so-called electro-active magnons, i.e., control of spin waves by electric fields [START_REF] Sousa | Electrical control of magnon propagation in multiferroic BiFeO3 films[END_REF][START_REF] Chang | Prediction of electromagnons in antiferromagnetic ferroelectrics from first-principles: The case of BiFeO3[END_REF][START_REF] Seki | Electromagnons in the Spin Collinear State of a Triangular Lattice Antiferromagnet[END_REF][START_REF] Chen | Electric-field control of electromagnon propagation and spinwave injection in a spiral multiferroic/ferromagnet composite[END_REF][START_REF] Davydova | Spin dynamics, antiferrodistortion and magnetoelectric interaction in multiferroics. The case of BiFeO3[END_REF]. On the other hand, much less work has been devoted to the study of magnetic control of polarization waves, with a few theoretical works realized in BiFeO 3 [START_REF] Sayedaghaee | Domain-wall-induced electromagnons in multiferroics[END_REF][START_REF] Sayedaghaee | Novel dynamical magnetoelectric effects in multiferroic BiFeO3[END_REF][START_REF] Sayedaghaee | Strain-induced resonances in the dynamical quadratic magnetoelectric response of multiferroics[END_REF] and manganites [START_REF] Stenberg | Model for twin electromagnons and magnetically induced oscillatory polarization in multiferroic RM nO 3[END_REF]. The present work goes one step further by bridging the two approaches, as we intend to demonstrate the possibility of resonantly exciting electromagnons (induced by ferroelectric domain walls) using ac magnetic fields while concurrently manipulating their frequency with dc electric fields.

Note that ferroelectric domain walls can now routinely be written, erased and reconfigured using, for instance, PiezoForce Microscopy [START_REF] Lachheb | Surface and bulk ferroelectric phase transition in super-tetragonal BiFeO3 thin films[END_REF][START_REF] Cherifi | Imaging ferroelectric domains in multiferroics using a low-energy electron microscope in the mirror operation mode[END_REF][START_REF] Gruverman | Piezoresponse force microscopy and nanoferroic phenomena[END_REF][START_REF] Kholkin | Review of Ferroelectric Domain Imaging by Piezoresponse Force Microscopy[END_REF]. The ability to generate magnetically polarization waves localized at domain walls, as proposed in Ref. [START_REF] Sayedaghaee | Domain-wall-induced electromagnons in multiferroics[END_REF], already promises the tantalizing possibility of reconfigurable nanometer size electrical circuits, whose power is switched on and sustained remotely by ac magnetic fields. If now one is able to act on the electrical polarized waves localized at the ferroelectric domain walls, for instance using local electric fields, one could dream of achieving reconfigurable logical elements for computing or detection.

In this work, we investigate how dc electric fields affect the domain-wall-induced electromagnons evidenced in Ref. [START_REF] Sayedaghaee | Domain-wall-induced electromagnons in multiferroics[END_REF] in multiferroic BiFeO 3 . Using Molecular Dynamics based on an ab-initio effective Hamiltonian, we reveal that the frequency of these electromagnons is rather sensitive to these dc fields, either increasing or decreasing with them depending on the chosen electromagnon. This latter different behavior (namely, increase versus decrease) is found to be correlated with the real-space localization of the optical phonon associated with these electromagnons. The energetics and properties of the studied system are modeled via the use the effective Hamiltonian framework detailed in Ref. [START_REF] Prosandeev | Novel nanoscale twinned phases in perovskite oxides[END_REF] within a molecular dynamics (MD) approach [START_REF] Prosandeev | Displacement current in domain walls of bismuth ferrite[END_REF], in order to obtain dynamical properties. Technically, Newtonian equations of motion are used to investigate the dynamics of the ionic degrees of freedom of this effective Hamiltonian -that are local modes which are directly proportional to local electric dipoles inside each 5-atom cell; pseudo-vectors representing the antiferrodistortive motions inside such 5atom cells; and both homogeneous and inhomogeneous strain components. Regarding the dynamics of the magnetic moments, the approach of Ref. [START_REF] Wang | Atomistic molecular dynamic simulations of multiferroics[END_REF] is followed, for which such dynamics are treated through the Landau-Lifshitz-Gilbert (LLG) equation [START_REF] García-Palacios | Langevin-dynamics study of the dynamical properties of small magnetic particles[END_REF].

II. METHODS

Here and as shown in

DC electric fields with magnitude ranging between 1.0 × 10 7 V/m and 1.0 × 10 8 V/m are applied along the [ 111] direction to such multidomain, that is currently mimicked by using a 24 × 24 × 6 supercell. It should be noted that an electric field of the magnitude of 2.0×10 8 V/m is strong enough to reorient all dipole moments along its direction and convert the configuration to monodomain. The MD simulations are conducted at the temperature of 10 K. To ensure that the magnetic moments only slightly fluctuate about the same direction during the course of simulations, a dc magnetic field of the magnitude of 245 T is applied along the [211] direction. A lower magnetic field (which is accessible in practical applications) could be applied. In fact, we also performed simulations with a dc magnetic field having a magnitude of 2.45 T along the [211] direction, and obtained similar results. Such a large magnetic field along with low temperature were chosen here to avoid large fluctuations of magnetic properties and see the results more clearly. Each MD simulation was performed for 1, 000, 000 steps of 0.5 f s each, while both homogeneous and inhomogeneous strains were relaxed.

III. RESULTS AND DISCUSSION

It is important to realize that, under zero dc electric field, there is no macroscopic polarization, as a result of the cancellations between the polarizations of the D1 and D2 domains.

On the other hand, we numerically found that the application of our considered dc electric fields results in the increase of the magnitude of electric dipoles in the D1 regions and its reduction in the D2 areas, therefore yielding now a finite overall polarization along the direction of the applied field which is [ 111]. It is worth mentioning that we observed the change in the structure from a 180-degree multidomain to a monodomain when the applied electric field exceeds a threshold (E threshold = 2 × 10 8 V /m) which is rather abrupt -at least with the resolution of the applied electric field change here, being E = 1 × 10 7 V /m. One may thus need a much smaller step in field's magnitude to observe motion of the domain walls.

One can also increase the temperature to see such motion. For instance, and as shown in the Supplementary Material, motion of the domain walls does occur at T = 300 K when the applied electric field varies from E = 3 × 10 7 V /m to E = 5 × 10 7 V /m.

We then perform Fourier analysis of the temporal evolution of this resulting electrical polarization along the [ 111] direction, and that of the magnetization along the [211] direction.

Figure 2 shows the resulting Fast Fourier Transformations (FFT) of both the electrical polarization and magnetization for a dc electric field of 1.0 × 10 7 V/m, which reveals the existence of frequency modes presenting phonons and magnons. Here, in particular, some peaks are observed at the same frequency in the FFT spectrum of both polarization and magnetization, which is a signature of electromagnons. Within the electromagnonic modes that are observed here, four modes are of particular interest due to their noticeable frequency shift as a response to the application of electric fields of various magnitudes. These modes are denoted by Mode 1, Mode 2, Mode 3, and Mode 4 in Figure 2 (Note that these four modes have very similar frequencies than those one can guess to be around 90 cm -1 (≃ 2700 GHz) and 140 cm -1 (≃ 4200 GHz) in the Supplemental materials of Ref. [START_REF] Hlinka | Terahertz-range polar modes in domainengineered BiFeO 3[END_REF] for precisely 180 • domains of BFO). The electric-field-induced frequency shift of both the phonon and magnon associated with the electromagnonic Mode 1 is demonstrated in Figs. 3a and3b, respectively. For this mode, as indicated by purple solid lines in Figs. 3a and3b, when the applied electric field is 1.0 × 10 7 V/m the phonon peak is at 2490 GHz (Figure 3a) and the magnon peak is exactly at the same frequency location (Figure 3b). When the applied electric field increases to higher magnitudes, the frequency of the phonon and magnon decrease concurrently which results in a softening for Mode 1. It should be noted that, in multiferroics, strain can contribute to the emergence of so-called electroacoustic magnons as a mixture of acoustic phonons, optical phonons, and magnons [START_REF] Sayedaghaee | Novel dynamical magnetoelectric effects in multiferroic BiFeO3[END_REF][START_REF] Sayedaghaee | Strain-induced resonances in the dynamical quadratic magnetoelectric response of multiferroics[END_REF]. Here, we numerically found (not shown here) that the four aforementioned modes are present in the frequency spectrum of both the polarization and magnetization even if the homogeneous strain is clamped during the course of MD simulations, which clarifies that these modes are electromagnons consisting of optical phonons and magnons. It should be noted that we do not believe that these computed intensities of Figure 3 For each of these four considered modes, the evolution of their frequency as a function of the magnitude of the dc electric field is shown in Figure 4. One can see that Mode 1 experiences a rather strong decrease of its frequency when the electric field varies from 1.0 × 10 7 V/m to 1.0 × 10 8 V/m, namely by about 200 GHz from 2490 GHz to 2290 GHz.

Note that it is known that the effective Hamiltonians of BFO overestimate the magnitude of electric fields by a factor of 23 [START_REF] Xu | Designing lead-free antiferroelectrics for energy storage[END_REF]. Hence, our maximum applied electric field would correspond experimentally to approximately 43 kV/cm, which is easily sustained in BiFeO 3 thin films [START_REF] Boyn | Learning through ferroelectric domain dynamics in solid-state synapses[END_REF]. Note also that the decrease of the Mode 1 frequency appears to be quadratic in nature. Mode 2 also adopts a reduction of its frequency, but at a smaller extent (that is by about 55 GHz from 2715 GHz to 2660 GHz) and in a linear fashion, with a rate of 0.06 GHz/(kV/cm). Interestingly, such linear variation has indeed been observed in BiFeO 3 for some electromagnons. For instance, the frequency of the so-called extra-cyclon mode ψ 2 may increase or decrease (depending on whether the electric field increases or decreases the polarization) with a rate of approximately 1.3 GHz/(kV/cm) [START_REF] Rovillain | Electric-field control of spin waves at room temperature in multiferroic BiFeO3[END_REF]. In the same reference, the cyclon mode ϕ 2 shows an opposite behavior with an electrical rate of control of the magnon frequency of ≈ 0.24 GHz/(kV/cm). Note that we focus here on domain-wallinduced electromagnons with an antiferromagnetic structure, while the cyclon and extracyclon modes considered in Ref [START_REF] Rovillain | Electric-field control of spin waves at room temperature in multiferroic BiFeO3[END_REF] are modes in a monodomain single crystal having a magnetic cycloidal state. In addition, the commonly known overestimation of electric fields in effective Hamiltonian models may also contribute to the difference between computational and experimental rates for the dc-field-induced change in frequency. As a matter of fact, rescaling our theoretical fields by dividing them by 23 (as indicated in Ref. [START_REF] Xu | Designing lead-free antiferroelectrics for energy storage[END_REF]) results in a rate for our Mode 2 that goes from 0.06 GHz/(kV/cm) to 1.38 GHz/(kV/cm), which is similar to the observed magnitude of such rate for the ψ 2 mode in Ref [START_REF] Rovillain | Electric-field control of spin waves at room temperature in multiferroic BiFeO3[END_REF]. Similarly, Mode 3 has the same kind of qualitative behavior than Mode 2 but with about twice the slopethat is, a linear decrease of its frequency by ≃ 117 GHz when the dc electric field strengthens from 1.0 × 10 7 V/m to 1.0 × 10 8 V/m. Strikingly, Mode 4, whose frequency is basically that of Mode 3 for an interpolated zero field, adopts a mirror behavior with respect to Mode 3, in the sense that its frequency concomitantly linearly increases by a similar amount of ≃ 117 GHz. Note that we also performed simulations with opposite dc electric fields (i.e., along [1 11 ]), which allows us to further determine that the frequencies of Modes 1 and 2 depend on the magnitude of the electric field along [ 111] or [1 11 ] while those of Modes 3 and 4 linearly depend on the projection of the electric field along [ 111] (i.e., on the magnitude but also sign of this projection).

In order to understand all these behaviors and demonstrate their relationship with realspace localization, a layer-by-layer analysis is performed at the different planes that are parallel to the domain wall. More precisely, the Fourier transform of the average of po-larization of each of these planes is computed for the frequencies associated with the four aforementioned modes for a dc field of 1.0×10 7 V/m, and is shown in Figure 5. For instance, Panel (a) of Figure 5 tells us that Mode 1 is a mode that is strongly localized at the domain walls. Furthermore, Panel (b) of Figure 5 reveals that Mode 2 also localizes near the domain walls but to a smaller extent, as seen by comparing its vertical scale with that of Figure 5 Interestingly, Figures 5(c) and 5(d) tell us that Modes 3 and 4 are rather different from Modes 1 and 2, in the sense that they prefer to localize inside the domains rather that at the domain walls. More precisely, Mode 3 reaches its maximum of the Fourier transform of the polarization in the D2 region inside which the polarization is antiparallel to the applied electric field. Consequently, applying such field will decrease the magnitude of the polarization in the D2 area, which corresponds to a local softening of the optical phonon, hence explaining the decrease of the frequency seen in Fig. 4c for Mode 3. In contrast, Mode 4 preferentially localizes in the D1 area for which the polarization is parallel to the dc electric field, and, as a result, such polarization increases in magnitude when the field strengthens. Such increase leads to a local hardening of Mode 4, therefore to a frequency that now increases with the field. Note that Modes 3 and 4 (whose frequency are about 4210 GHz and 4245 GHz for a field of 1.0 × 10 7 V/m, which correspond to 140 cm -1 and 142 cm -1 , respectively) can be thought as both originating from the known zone-center optical phonon of BiFeO 3 monodomain, that then splits in two under dc electric fields because of the existence of domain walls and two different types of domains in our studied system. Based on previous works on BiFeO 3 monodomain single crystals and the fact that we numerically further found (not show here) that Modes 3 and 4 have rather small FFT of the component of the polarization along the [110]-direction (which is perpendicular to the polarizations of both D1 and D2), one can suggest that Modes 3 and 4 originate from the A 1 (LO) mode, rather than the E(T O) mode, of BiFeO 3 monodomain [START_REF] Wang | Atomistic molecular dynamic simulations of multiferroics[END_REF][START_REF] Cazayous | Electric field effect on Bi Fe O 3 single crystal investigated by Raman spectroscopy[END_REF][START_REF] Lobo | Infrared phonon dynamics of a multiferroic BiFeO3 single crystal[END_REF][START_REF] Haumont | Phonon anomalies and the ferroelectric phase transition in multiferroic BiFeO3[END_REF][START_REF] Kamba | Infrared and terahertz studies of polar phonons and magnetodielectric effect in multiferroic BiFeO3 ceramics[END_REF].

IV. SUMMARY

In summary, we used an atomistic effective Hamiltonian to reveal that the frequency of electromagnons can be significantly tuned by applying dc electric fields in the prototypical BiFeO 3 multiferroic adopting ferroelectric domains. Such finding is promising towards the design of novel devices taking advantage of the dual electric and magnetic natures of electromagnons, with the additional conveniences demonstrated here that (1) it should thus be possible to select the desired operating frequency by choosing the right combination of ac magnetic field's frequency and magnitude of the dc electric field (in order to activate such electromagnons at this desired frequency); and (2) some of these electromagnons are localized at the ferroelectric domain walls, therefore rendering feasible the application of local electric fields for realizing reconfigurable logical elements for computing or detection. These domain-wall-induced electromagnons are further found to either increase or decrease their frequencies under the dc electric fields, depending on the real-space localization of their associated phonons-that is at the ferroelectric domain walls or at the "up" versus "down" domains. We therefore hope that the present study deepens the fascinating fields of electromagnons, ferroelectric domains and magnonics.

Figure 1 ,

 1 we simulate a multidomain system of BFO with 180 • domain configuration inside which two types of domains alternate along the [110] pseudocubic direction. The first type of domain, denoted as D1, is shown in red in Fig. 1 and possesses electric dipoles aligned along the [ 111] direction. The second kind of domain, coined D2, is displayed in blue in this figure and exhibits electric dipoles lying along the opposite [1 11 ] direction. In both domains the magnetic moments and antiferromagnetic (AFM) moments are along the perpendicular [211] and [0 11] directions, respectively.

  have a real physical meaning in our simulations, since such intensities are not monotonic with the fields and can vary when slightly changing some technical details of the Fourier Transform (especially, considering the small values of the vertical scale in the Figure). On the other hand, the frequency position of these peaks is insensitive to such details and does carry physical significance.

  (a). Consequently, by looking at Figs 4(a), 4(b), 5(a) and 5(b), one can conclude that modes localizing near the domain walls soften under a dc electric field, that is they have their frequency decreasing when the field increases -and such decrease is larger when the localization near the walls is stronger.

FIG. 1 : 3 FIG. 2 : 4 FIG. 3 :

 13243 FIG. 1: (Color online) Schematic representation of the zero-field electric dipole moments' pattern for our studied 24 × 24 × 6 supercell of BiFeO 3 . The blue and red vectors are used to represent electric dipole moments along the [1 11 ] and [ 111] directions, respectively.

FIG. 4 :

 4 FIG.4: (Color online) Frequency shift of four electromagnons. In each panel, the frequency shift of the corresponding mode versus the applied electric fields is shown. A polynomial fit of second order for Mode 1, and of first order for Modes 2, 3, and 4 is applied. |∆ω| is the magnitude of the difference between the highest and the lowest frequencies of the fitted lines for our chosen range of applied electric fields.

4 FIG. 5 :

 45 FIG. 5: (Color online) The degree of localization of the electromagnons possessing significant frequency shift under the application of electric fields. In each panel, the degree of localization of the corresponding mode at the planes parallel to the domain walls within the investigated supercell is shown. The shaded area corresponds to the domains possessing electric dipoles pointing along the [1 11 ] direction, while the white area corresponds to the domains possessing electric dipoles pointing along the opposite [ 111] direction. The arrow in each domain indicates the direction of the z-component of electric dipole moments present in the domain. The applied electric field for all cases is E = 1.0 × 10 7 V /m along the [ 111] direction which is parallel to the direction of the electric dipoles in D1 and antiparallel to the direction of the electric dipoles in D2.

ACKNOWLEDGEMENTS

The authors thank the Vannevar Bush Faculty Fellowship (VBFF) Grant No. N00014-20-1-2834 from the Department of Defense and the ARO Grant No. W911NF-21-1-0113. C.P. thanks the ANR Grant No. ANR-21-CE24-0032 SUPERSPIN, and a public grant overseen by the ANR as part of the Investissements d'Avenir program (Reference: ANR-10-LABX-0035, Labex-NanoSaclay). S.P. also acknowledges ONR Grant No. N00014-21-1-2086. B.X. further acknowledges the financial support from National Natural Science Foundation of China under Grant No. 12074277, the startup fund from Soochow University, and the support from Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions. S.P. and Y.N. thank DARPA Grant No. HR0011727183-D18AP00010 (under the TEE Program). The authors also acknowledge the High Performance Computing Center


Carmel Dansou

Charles Paillard

Laurent Bellaiche

Controlling properties of PbTiO 3 /SrTiO 3 superlattices by photo-excited carriers

The effect of light on properties of PbTiO3/SrTiO3 superlattices is studied using first-principles calculations in their paraelectric phase p as well as in the, c, r and aa phases that are defined by polarization lying along [001], [11v] and [110] directions, respectively. Whilst the r phase is the ground state in dark conditions, illumination triggers a rotation of the polarization towards the outof-plane direction, which results in the r phase transforming into a c phase at large concentration of photoexcited carriers. We further reveal that a space charge region of free carriers builds up at the interface under illumination, implying that these superlattices can be thought as acting as p-i-n junctions under light. These free charges screen the polarization discontinuity at the interface, effectively reducing the electrostatic depolarizing field.

I. INTRODUCTION

(PbTiO 3 /SrTiO 3 ) m superlattices ( (PTO/STO) m SLs) have gained a lot of attention [START_REF] Dawber | Unusual behavior of the Ferroelectric polarization in PbTiO3/SrTiO3 superlattices[END_REF][START_REF] Dawber | Tailoring the properties of artificially layered ferroelectric superlattices[END_REF][START_REF] Aguado-Puente | Interplay of couplings between antiferrodistortive, ferroelectric, and strain degrees of freedom in monodomain PbTiO3/SrTiO3 superlattices[END_REF][START_REF] Bousquet | Improper ferroelectricity in perovskite oxide artificial superlattices[END_REF][START_REF] Misirlioglu | Space charge contribution to the apparent enhancement of polarization in ferroelectric bilayers and multilayers[END_REF][START_REF] Cooper | Polarization enhancement in short period superlattices via interfacial intermixing[END_REF][START_REF] Das | Observation of room-temperature polar skyrmions[END_REF][START_REF] Yadav | Observation of polar vortices in oxide superlattices[END_REF][START_REF] Abid | Creating polar antivortex in PbTiO3/SrTiO3 superlattice[END_REF][START_REF] Aguado-Puente | Structural and energetic properties of domains in PbTiO3/SrTiO3 superlattices from first principles[END_REF][START_REF] Zubko | Ferroelectric Domains in PbTiO3/SrTiO3 superlattices[END_REF]. It was found that, by turning the epitaxial strain and the PTO volume fraction, a large variety of fascinating properties can be engineered in these SLs [START_REF] Dawber | Unusual behavior of the Ferroelectric polarization in PbTiO3/SrTiO3 superlattices[END_REF][START_REF] Dawber | Tailoring the properties of artificially layered ferroelectric superlattices[END_REF][START_REF] Bousquet | Improper ferroelectricity in perovskite oxide artificial superlattices[END_REF]. They have, e.g., been found to host non trivial topological structures [START_REF] Das | Observation of room-temperature polar skyrmions[END_REF][START_REF] Yadav | Observation of polar vortices in oxide superlattices[END_REF][START_REF] Abid | Creating polar antivortex in PbTiO3/SrTiO3 superlattice[END_REF][START_REF] Aguado-Puente | Structural and energetic properties of domains in PbTiO3/SrTiO3 superlattices from first principles[END_REF] and negative capacitance [START_REF] Zubko | Negative capacitance in multidomain ferroelectric superlattices[END_REF], which are of fundamental as well as practical interest [START_REF] Das | A new era in ferroelectrics[END_REF].

There have also been significant developments in the interaction of light with ferroelectric (FE) materials during the last decade [START_REF] Kundys | Photostrictive materials[END_REF][START_REF] Paillard | Photovoltaics with Ferroelectrics: Current Status and Beyond[END_REF], with light having the capability to induce structural change (photostriction) and even phase transition in classical FE materials [START_REF] Kundys | Photostrictive materials[END_REF][START_REF] Paillard | Photostriction in ferroelectrics from density functional theory[END_REF][START_REF] Chen | Deterministic control of ferroelectric polarization by ultrafast laser pulses[END_REF][START_REF] Haleoot | Photostrictive Two-Dimensional Materials in the Monochalcogenide Family[END_REF][START_REF] Paillard | Ab initio approach to photostriction in classical ferroelectric materials[END_REF]. More recently it was also reported that one can control the topology of phonons in FE materials with light [START_REF] Peng | Topological phonons in oxide perovskites controlled by light[END_REF].

It is now natural and timely to ask whether light can modulate physical properties of FE SLs. For example, it was reported that under illumination, PTO/STO SLs transform from a polydomain structure to a uniform domain configuration [START_REF] Ahn | Photoinduced domain pattern transformation in ferroelectric-dielectric superlattices[END_REF]. This transition was accompanied by the expansion of the out-of-plane lattice of the SLs. A light-induced polarization rotation was also evidenced at the nanoscale in multidomains [START_REF] Lee | Structural evidence for ultrafast polarization rotation in ferroelectric/dielectric superlattice nanodomains[END_REF], as well as an optically-induced phase transition from a vortex and domain state to a so-called super crystal phase in PTO/STO SLs [START_REF] Stoica | Optical creation of a supercrystal with three-dimensional nanoscale periodicity[END_REF]. However, light interaction with FE SLs, has been mostly experimentally studied so far. To the best of our knowledge, theoretical investigations were limited to phenomenologically varying the dielectric constant of the dielectric layer to reproduce the super crystal phase [START_REF] Stoica | Optical creation of a supercrystal with three-dimensional nanoscale periodicity[END_REF]. Microscopic insights are thus lacking. Many questions remain unresolved, such as what are the microscopic mechanisms favoring photo-induced transitions in these SLs or change in their physical properties under light.

Here, we explore light interaction with monodomain (PbTiO 3 /SrTiO 3 ) m SLs. Technically, we apply our recently developed scheme to mimic the effect of photo-induced thermalized carriers [START_REF] Paillard | Photoinduced phase transitions in ferroelectrics[END_REF] on the relative energy and phase stability of different phases having different electrical polarization orientation. We observe photoinduced phase transitions and change in properties, and further reveal their microscopic origins via the development of an analytical model demonstrating that these SLs can be thought as p-i-n junctions when under illumination.

II. METHODS

A. Constraint density functional theory

The constraint density functional theory (c-DFT) is a purpose built method to simulate photo-excitation in pump-probe optical experiment within the known DFT framework. Photo-excited thermalized electrons and holes are described using Fermi-Dirac distribution with two chemical potentials, µ e for electrons and µ h for holes. During the self consistent field DFT cycles, the electronic density is converged with the constraint of having n e = n ph electrons per formula unit in the conduction band and n h = n ph holes in the valence band. Explicitly one solves numerically for the couple of chemical potentials (µ e , µ h ) using the following equations at each SCF cycle:

n h = i≤Nv k,σ w k [1 -f ( ik,σ , µ h )] (1
)

n e = i>Nv k,σ w k f ( ik,σ , µ e ) ( 2
)
w k is the weight of the k-point k, f ( ik,σ ) is the occupation number, σ stands for spin and N v is the index of highest occupied band [START_REF] Paillard | Photoinduced phase transitions in ferroelectrics[END_REF]. This scheme recently implemented in the Abinit package has proven to successfully mimic the effects of photo-excited carriers in ferroelectrics [START_REF] Paillard | Photoinduced phase transitions in ferroelectrics[END_REF].

B. Computational details

We performed density functional theory (DFT) calculations as implemented in the Abinit package [START_REF] Romero | Abinit: Overview and focus on selected capabilities[END_REF] with the projector augmented wave method [START_REF] Torrent | Implementation of the projector augmented-wave method in the ABINIT code: Application to the study of iron under pressure[END_REF]. The generalized gradient approximation with the Perdew-Burke-Ernzerhof parameterization revised for solid was used as the exchange-correlation functional [START_REF] Perdew | Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces[END_REF]. A 20 Ha plane wave cut-off was employed in our calculations. A 6×6×1 k-meshes equivalent to a 6 × 6 × 6 mesh in a five-atom unit cell allows to reach convergency for both forces and energies. Density was considered converged when the difference of forces on atoms between two self-consistent field iterations is less that 1 × 10 -7 Ha.Bohr -1 . Both cell shape and ionic positions were further optimized until the maximum forces on the ions were less than 2 × 10 -5 Ha.Bohr -1 .

We simulated monodomains in (PTO) 5 /(STO) 5 SLs that are built by stacking along the [001] direction five unit cells of PTO and five unit cells of STO. The choice of a 5/5 PTO/STO SLs as a model for this study is inspired by previous studies and our computational resources. When simulating light interaction with the superlattices, calculations become very heavy for larger supercells. As such we had to choose the size of the superlattices to be thick enough but in accordance with the resource's constraints. Also, a 5/5 PTO/STO film has been experimentally grown [START_REF] Zubko | Ferroelectric Domains in PbTiO3/SrTiO3 superlattices[END_REF], which demonstrates the experimental relevance of our choice. The in-plane lattice constant was kept fixed at the theoretical optimized lattice constant of STO, in order to mimic a growth on a STO substrate. First, the SLs were optimized to find a reference paraelectric structure with mirror symmetry called the p phase [START_REF] Pertsev | Effect of mechanical boundary conditions on phase diagrams of epitaxial ferroelectric thin films[END_REF]. Three other phases are then considered: the c phase with a polarization along the out-of-plane [001] direction; the aa phase with a polarization lying along the in-plane [110] direction; and the r phase with a polarization being along [uuv] directions [START_REF] Pertsev | Effect of mechanical boundary conditions on phase diagrams of epitaxial ferroelectric thin films[END_REF]. Starting from the optimized structure under dark conditions, light-interaction with these phases was then simulated using the c-DFT. When under light, both cell shape and ions positions of the structures were again relaxed with the aforementioned parameters.

III. RESULTS

A. Energy and stability

On Figure 1 , we show the relative energies of the different phases (with respect to the p reference phase) for various excited charged carriers concentration (n ph ). In dark condition (n ph = 0 electron per formula unit (e/f.u)), the most stable phase is the r phase where the polarization points along a specific [11v] direction, such that it makes an angle θ ≈ 72 • with the z-axis. As we increase the light intensity (n ph ), the relative energy of the c phase gets closer to that of the r phase until the two basically coincide for n ph larger than 0.2 e/f.u. As we will see later, such behaviors originate from the facts that, under light, the r phase first rotates its polarization towards the [001] direction and then transforms into the c phase that becomes the ground state. Meanwhile, the energy of the aa phase increases away from its ground state's energy, as n ph increases.

B. Light-induced polarization and lattice dynamism in the c phase

To further determine effects of light on the different phases, we also compute the polarization in layer-by-layer fashion using the Born effective charges of the atoms in each local unit cell via:

P i = e Ω i j∈i w j δ j Z * j (3)
where Ω i , e, and δ j are the volume of the unit cell i, the elementary charge and the atomic displacement from the reference paraelectric p phase, respectively. Z * j is the Born effective charge of the atom j, with j running through all the atoms in the cell i and w j being a weight corresponding to the number of five-atoms unit cell that share the atom j. The Born charges Z * j are taken from Refs. [START_REF] Zhong | Giant LO-TO splittings in perovskite ferroelectrics[END_REF].

The data are shown in Figure 2 (a) for the c phase. In dark condition, the electrical polarization throughout the SLs is nearly constant, as consistent with previous reports [START_REF] Neaton | Theory of polarization enhancement in epitaxial BaTiO3/SrTiO3 superlattices[END_REF][START_REF] Zhu | Firstprinciple study of ferroelectricity in PbTiO3/SrTiO3 superlattices[END_REF] and which arises from the minimization of the electrostatic cost associated with bound charges at the interfaces. Under illumination, the electrical polarization in PTO layers increases while the polarization in STO layers decreases. The polarization in the interfacial regions (unit cells 5 and 10 on Figure 2 and (d) and as also found in optical excitation experiment on PTO/STO SLs having multidomains [START_REF] Ahn | Photoinduced domain pattern transformation in ferroelectric-dielectric superlattices[END_REF]. The local out-of-plane distortion is quantified as (z i A2 -z i A1 )/a where z i A1 , z i A2 are the zpositions of the A j=1,2 = Pb or Sr atoms at the corners of the cell i in the SLs with z i A1 < z i A2 ; and a is the in-plane lattice constant of the SLs.

C. Light-induced polarization dynamism in the aa phase

We also looked at the aa phase and the effect of light on its polarization. On Figures 3 (a&b), we show the total and local in-plane polarizations of this aa phase, as function of light intensity. The total electrical polarization of the aa phase decreases with n ph , indicating that light penalizes in-plane polar displacements. This finding is also consistent with the steady increase in the energy of this aa phase seen in Fig. 1. Furthermore, inspection of the layer-by-layer polarization in the aa phase, as shown on Figure 3 (b), indicates that the decrease of its polarization mostly comes from the STO layers that tend to become paraelectric under light (in dark condition, the STO layers are ferroelectric as a consequence of their proximity with PTO layers within our considered supercell).

D. Light-induced polarization rotation in the r phase

The observations from the c and aa phases imply that in the r phase, the in-plane component of the electrical polarization should be reduced while its out-of plane component should be increased by illumination. This should thus lead to the polarization in the r phase rotating towards the out-of-plane direction upon illumination, and to the transition from the r to the c phase observed in the energy diagram on Figure 1. We map out the angle between the z-axis and the electrical polarization within the r phase as function of n ph . The results shown on Figs. 3 (c &d) are well consistent with our above explanations. We therefore have a light-induced polarization rotation within the r phase. Similar light-induced polarization rotation has been seen at domain boundaries in multidomain PTO/STO SLs where the polarization rotates from an r-like orientation to a c-like orientation upon optical excitation [START_REF] Lee | Structural evidence for ultrafast polarization rotation in ferroelectric/dielectric superlattice nanodomains[END_REF], in-line with our predictions. In Fig. 3 (d), we show the variation of the angle of the polarization with the z axis across the SLs. Under illumination, the polarization in the STO layer rotates faster towards the z axis and this angle is almost zero for all value of n ph > 0.06 e/f.u which indicates that, beyond a certain light intensity, the polarization in the STO layer stays practically along the z axis. Any further rotation of the polarization in the SLs beyond n ph > 0.06 e/f.u, therefore, emanates mostly from the rotation of the local polarization in the PTO layers.

E. Light-induced depolarization field screening

Next, we look at the macro-averaged [START_REF] Junquera | Nanoscale smoothing and the analysis of interfacial charge and dipolar densities[END_REF][START_REF] Baldereschi | Band Offsets in Lattice-Matched Heterojunctions: A Model and First-Principles Calculations for GaAs/AlAs[END_REF] (also detailed in Ref.

[34], Sec. 1) electrostatic potential and charges density inside the c phase as function of n ph . From the macro-averaged potential, we deduce the macro-averaged internal electric field by numerical finite difference differentiation of the potential data with respect to z [START_REF] Kong | Python Programming and Numerical Methods: A Guide for Engineers and Scientists[END_REF]. In dark condition and as shown in Figure 4(a) via the blue curve, there is a finite electric field in both PTO and STO layers, consistent with previous reports on FE-dielectrics SLs [START_REF] Neaton | Theory of polarization enhancement in epitaxial BaTiO3/SrTiO3 superlattices[END_REF][START_REF] Sepliarsky | Long-ranged ferroelectric interactions in perovskite superlattices[END_REF], with short-circuit boundary conditions requiring the fields in the PTO and STO layers to be opposite in sign. It is the polarization mismatch between PTO and STO layers that induces this electric field in the STO layer [START_REF] Lim | Electrostatic coupling and interface intermixing in ferroelectric superlattices[END_REF][START_REF] Roytburd | Dielectric anomaly due to electrostatic coupling in ferroelectricparaelectric bilayers and multilayers[END_REF]. These internal fields act as depolarizing fields limiting the polarization in the SLs. As we photo-excite charged carriers, they The different potentials for n ph = 0.4 e/f.u and n ph = 0 e/f.u; b) the difference of the macroscopic average excited charged carriers density between the case of n ph = 0.4 e/f.u and n ph = 0 e/f.u; c) E-fields for n ph = 0.4 e/f.u and n ph = 0 e/f.u; d) total E-field in different part of the Sls as a function of light intensity. mainly populate the interface regions of the PTO/STO SLs as evidenced by the difference, shown in Fig. 4 (b), in the macroscopic average in the charged carriers density between one illumination (n ph = 0.4 e/f.u) and dark conditions (the same trend is obtained for other values of n ph ). Interfaces are populated by opposite but equal amount of excited charged carriers. Our results agree with some existing experimental reports indicating that excited charges carriers migrate to the interfaces to screen the depolarizing field [START_REF] Ahn | Photoinduced domain pattern transformation in ferroelectric-dielectric superlattices[END_REF][START_REF] Lee | Structural evidence for ultrafast polarization rotation in ferroelectric/dielectric superlattice nanodomains[END_REF]. However, for the considered SLs thickness, the excited charged carriers are rather significantly localized at the interfaces while they tail into PTO and STO layers. This feature has been observed in previous studies on electron gas at LaAlO 3 /SrTiO 3 where for thin film, the free gas was observed to also tail into the STO layers [START_REF] Delugas | Spontaneous 2-Dimensional Carrier Confinement at the n-Type SrTiO3/LaAlO3 interface[END_REF][START_REF] Stengel | First-Principles Modeling of Electrostatically Doped Perovskite Systems[END_REF]. Consequently, the internal electric fields in PTO and STO layers reduce their value under illumination (see the red curve versus the blue curve in Figure 4 (c)) and even quench for n ph equal or larger than 0.2 e/f.u -as revealed by Fig. 4 (d). For each considered finite n ph , one can also define a lightinduced potential that consists of the difference between the total electrostatic potential for such finite n ph and the one for n ph = 0, as shown via black solid lines in Figure 4 (a). Such light-induced potential can be understood as a light-induced electric field (see Figure 4 (c) for the case of n ph = 0.4 e/f.u. ). In both PTO and STO layers, light produces fields opposite to the intrinsic internal field leading to a net zero field in these layers (Figure 4 (c)). We therefore have a light-induced depolarizing field screening. In contrast, light has very negligible effect on the fields at the two interfaces, as demonstrated in both Figures 4 (c) and 4 (d).

F. Electro-static model of light-induced dynamism

In the following, we propose a simple electrostatic model to account for the light-induced change in electrostatics in the SLs. For simplicity, we assume that the excited charged carriers are significantly localized and homogeneously populate the interface regions, and that PTO and STO are mainly recombination regions where there is no net excited charged carriers. In one unit cell of the SLs, the volume charge density of excited carriers can be defined as follows:

ρ l (z) =          0, z 0 ≤ z ≤ z 1 +ρ v , z 1 ≤ z ≤ z 2 0, z 2 ≤ z ≤ z 3 -ρ v , z 3 ≤ z ≤ z 4 , ρ v = z2 z1 ρ(z)dz z 2 -z 1 . (4)
where ρ(z) are the difference in the macro-averaged charged density (see, e.g., Fig. 4 (b) for the case of (n ph = 0.4 e/f.u); z 0 and z 1 the boundaries of the PTO layer; z 1 and z 2 the boundaries of PTO-STO interface; z 2 and z 3 the boundaries of the STO layer and z 3 and z 4 those of the STO-PTO interface. The values of z 1 , z 2 , z 3 , z 4 are extracted from our computational data and are (z 0 , z 1 , z 2 , z 3 , z 4 ) ≈ (0.32, 16.10, 20.06, 35.84, 39.63) Å.

The above-defined charge density is similar to charges density in a typical p-i-n junction in semiconductors [START_REF] Berz | A simplified theory of the p-i-n diode[END_REF][START_REF] Streetman | Solid state electronic devices[END_REF][START_REF] Baliga | Fundamentals of power semiconductor devices[END_REF]. Using basic electrostatic equations:

dE l (z) dz = - d 2 V l (z) dz 2 = ρ v . (5)
where is the static permittivity which value changes throughout the SLs, i.e = P T O in the PTO layer, = ST O in the STO layer, and the average of the dielectric permitivities of PTO and STO is taken as an approximation of at the interfaces. Note, however, that we did not use any previously computed values of but rather extracted the ratio ρ v / ≈ 0.02 V.

Å-2 from the slope of the computed light-induced field in the interface regions (Figure 4 (c)). The light-induced electrostatic field and potential are as follows:

E l (z) =              -E l 0 , z 0 ≤ z ≤ z 1 -E l 0 + z -z 1 ρ v , z 1 ≤ z ≤ z 2 +E l 0 , z 2 ≤ z ≤ z 3 E l 0 - z -z 3 ρ v , z 3 ≤ z ≤ z 4 ; V l (z) =                (z -z 0 )E l 0 + V 0 , z 0 ≤ z ≤ z 1 V l 1 + (z -z 1 )E l 0 - (z -z 1 ) 2 2 ρ v , z 1 ≤ z ≤ z 2 V l 2 -(z -z 2 )E l 0 , z 2 ≤ z ≤ z 3 V l 3 + (z 3 -z)E l 0 + (z -z 3 ) 2 2 ρ v , z 3 ≤ z ≤ z 4 (6)
with E l 0 ≈ 0.036 V. Å-1 being the strength of the lightinduced field at z 0 = 0.32 Å and (V l 0 , V l 1 , V l 2 , V l 3 ) ≈ (-0.26, 0.26, 0.26, -0.26) V being the electrostatic potentials induced by light at z 0 , z 1 , z 2 , z 3 respectively -as extracted from DFT data plotted on Figs. 4 (a&b). On Figure 5, we show the results of our analytical model and compared them with the data from our DFT computations for the highest studied light intensity (n ph = 0.4 e/f.u). The agreement is rather good for both the light-induced field and electrostatic potential. One can then conclude that, at the microscopic level, illuminated PTO/STO SLs act like p-i-n junctions with (1) the PTO-STO interface being the p-region; (2) the STO-PTO interface being the n-region; and (3) PTO and STO layers being the i-regions (intrinsic zone of neutral charges sandwiched between the p and n regions). The photo-induced electric fields produced by this (p-i-n)-like junction in PTO and STO layers screen out the depolarization field inside the PTO and STO layers. The supplemental material reports the results of another model that is similar in spirit with the one described by equation (4) (that is, excited charged carriers are strongly localized at the interfaces) but with different form for ρ(z) (See Ref.

[34], Sec. 2).

IV. DISCUSSION

Here we would like to point to some existing models proposed in literature to produce screening layer at PTO/STO interface by turning PTO films thickness [START_REF] Aguado-Puente | Model of twodimensional electron gas formation at ferroelectric interfaces[END_REF][START_REF] Yin | Twodimensional electron gas at the PbTiO3/SrTiO3 interface: An ab initio study[END_REF][46]. In those models, the free carriers can only be produced when the electrical polarization is normal to the interface. As such, one cannot use these models to study the impact of free carriers on the electrical polarization in the aa and the r phases considered in our study. Also, we did not consider oxygen octahedral rotation or tilting in this work that are are known to compete with polarization in some systems. This is mostly because we first wanted to know the effect of light on polarization alone in the studied superlattices but also because neither PTO nor STO possesses such rotation or tilting at room temperature (which is the temperature at which measurements are typically done). The effect of excited charged carriers on these antiferrodistortive modes and their coupling with polarization in these SLs, as well as the additional effect of varying the epitaxial strain, are interesting questions to explore for future studies.

Finally, note that in the present study, we limit ourselves to monodomains in these SLs in order to first reveal and understand the impact of photoexcitation on different polar phases before tackling the more complex impact of light on polydomains. The extension of this study to polydomain SLs is left for future exploration.

V. CONCLUSION

In summary, we report a first-principle investigation of the impact of light on energetic and physical properties of monodomain (PTO) 5 /(STO) 5 SLs. From energetics and electrostatic analysis, we now understand that the polarization in close to open-circuit conditions will grow under illumination as photo-excited carriers will screen the polarization bound charges. On the other hand, the polarization in directions close to short-circuit conditions reduces under illumination, as reminiscent of results in bulk ferroelectrics [START_REF] Paillard | Photostriction in ferroelectrics from density functional theory[END_REF]. In our specific case, this results in the polarization rotation of the r phase of monodomain PTO/STO SLs towards that of the c phase under illumination. The favoring of an out-of-plane polarization under illumination thus helps explaining the transition to a monodomain state observed in Ref. [START_REF] Ahn | Photoinduced domain pattern transformation in ferroelectric-dielectric superlattices[END_REF] or the domain boundary rotation discussed in Ref. [START_REF] Lee | Structural evidence for ultrafast polarization rotation in ferroelectric/dielectric superlattice nanodomains[END_REF]. We therefore hope that the present study deepens the active field of light-matter interaction and can be used to design novel efficient devices using light as a physical handle since, e.g., rotation of polarization is known to dramatically enhance many physical properties, including piezoelectricity, dielectric responses as well as elasto-optical and electro-optical conversions [START_REF] Bellaiche | Finitetemperature properties of PbZr (1-x) TixO3 alloys from first principles[END_REF][START_REF] Bellaiche | Low-temperature properties of PbZr (1-x) TixO3 solid solutions near the morphotropic phase boundary[END_REF][START_REF] Chen | Large elasto-optic effect in epitaxial PbTiO3 films[END_REF][START_REF] Paillard | Strain engineering of electro-optic constants in ferroelectric materials[END_REF].

FIG. 1 .

 1 FIG.1. Relative energies of the different phases. The zero of energy is chosen to be that of the p phase, for any illumination.

  (a)) remains

FIG. 2 .

 2 FIG. 2. Layer-by-layer polarization (a), local out-of-plane lattice distortion (b), total polarization (c) and relative lattice expansion in the c phase.

  )

FIG. 3 .

 3 FIG. 3. Total in plane polarization (a) and local polarization (b) in the aa phase; angle of the electrical polarization with the z-axis in the r phase: a) global angle, b) local angle.

2 E x c i t e d c a r r i e r s d e n s i t y ( 1

 21 FIG. 4.The different potentials for n ph = 0.4 e/f.u and n ph = 0 e/f.u; b) the difference of the macroscopic average excited charged carriers density between the case of n ph = 0.4 e/f.u and n ph = 0 e/f.u; c) E-fields for n ph = 0.4 e/f.u and n ph = 0 e/f.u; d) total E-field in different part of the Sls as a function of light intensity.

5 E x c i t e d c a r r i e r s d e n s i t y ( 1 2 E x c i t e d c a r r i e r s d e n s i t y ( 1 FIG. 5 .

 51215 FIG. 5. Comparison between the DFT data (first row) and the presently developed analytical model (second row) about light-induced electrostatics for the case of n ph = 0.4 e/f.u. a) and d) correspond to difference in the macro-averaged charged density, b) and e) represent the light-induced electric field, c) and f) correspond to the light-induced electrostatic potential.

ACKNOWLEDGMENTS

The work is supported by the ARO Grants No. W911NF-21-2-0162 (ETHOS) and W911NF-21-1-0113, and the Vannevar Bush Faculty Fellowship (VBFF) grant no. N00014-20-1-2834 from the Department of Defense. C.P. thanks the support from a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (Labex NanoSaclay, reference: ANR-10-LABX-0035).