Strain-induced resonances in the dynamical quadratic magnetoelectric response of multiferroics

For the last few years, the research interest in magnetoelectric (ME) effect, which is the cross-coupling between ferroelectric and magnetic ordering in multiferroic materials, has experienced a significant revival. The extensive recent studies are not only conducted towards the design of sensors, actuators, transducers, and memory devices by taking advantage of the cross-control of polarization (or magnetization) by magnetic (or electric) fields, but also aim to create a clearer picture in understanding the sources of ME responses and the novel effects associated with them. Here we derive analytical models allowing to understand the striking and novel dynamics of ME effects in multiferroics and further confirm it with atomistic simulations. Specifically, the role of strain is revealed to lead to the existence of electroacoustic magnons, a new quasiparticle that mixes acoustic and optical phonons with magnons, which results in resonances and thus a dramatic enhancement of magnetoelectric responses. Moreover, a unique aspect of the dynamical quadratic ME response under a magnetic field with varying frequencies, which is the second harmonic generation (SHG), has not been discussed prior to the present work. These SHGs put emphasis on the fact that nonlinearities should be considered while dealing with such systems.


INTRODUCTION
Magnetoelectric (ME) materials couple electrical and magnetic degrees of freedom. Discovered in 1960 1,2 , researchers' interest in them faded after a decade due to difficult integration into actual technological devices and lack of microscopic understanding. Yet, a new and vigorous tide of research 3-6 is rising to design sensors, actuators, transducers, and memory devices by taking advantage of the cross-control of polarization (magnetization) by magnetic (electric) fields [7][8][9][10] . Such tide also aims at creating a deeper understanding of ME effects [11][12][13] to properly engineer those aforementioned devices. In particular, dynamical ME effects, which are critical for fast computing operations, have not received the attention they deserve. In this work, an analytical model is derived to shed light on the dynamical strain-induced aspects of the ME effect in multiferroic materials. The model is then verified by firstprinciples-based simulations within an atomic effective Hamiltonian scheme.
In ME devices, one is often interested in the response of the electric polarization P (or, conversely, the magnetization M) to an applied magnetic field H (or, conversely, an electric field E). For the longest time, the direct control of polarization by magnetic fields, hence the direct coupling of spins with dipole moments, has been the main focus of ME research in single-phase materials. It has however been reckoned that the use of heterostructures and laminar piezoelectric-magnetostrictive composites 10,[15][16][17] could also lead to control of the polarization, by using strain as a mediator between electric dipoles and spins 18,19 . Those separate mechanisms lead, in the former case of a "direct" coupling, to the formation of socalled electromagnon quasiparticles 12,20,21 (mix of optical phonons and magnons). On the other hand, the latter "indirect" coupling was recently proposed to generate a new quasiparticle, the electroacoustic magnon 14 , which essentially mixes acoustic and optical phonons with magnons. We note that most studies of the strainmediated dynamical ME effect focused on (i) the linear ME effect and (ii) heterostructures stacking piezoelectric and ferromagnetic films [22][23][24][25] . Our work, in comparison, focuses on strain-induced resonances in the dynamical quadratic effect in a single-phase material. This is potentially important as (i) the linear ME effect can vanish or be negligible even in ME materials 26 and (ii) working with a single-phase material may lead to easier processing and lower manufacturing costs.

RESULTS AND DISCUSSION
There are various ways to show that long-wavelength acoustic phonons (i.e., strain) can mediate ME effects. Here we start from a model Hamiltonian that bears resemblance with a Landau expansion of the free energy in terms of powers of the polarization P, magnetization M, and strain η: with m P and m η representing a mass for the polarization and strain degrees of freedom, whereas k and b are constants. ϵ is the coupling constant of the direct bi-quadratic ME effect, Q is an electrostrictive coupling constant, λ is a magnetostrictive coupling constant, and c represents the elastic constant of the system. E represents the external electric field, which we assume to vanish here. We use a kinetic energy that is quadratic in the time derivative of strain, because it mimics what is in the effective Hamiltonian on which our MD equations are based. Such form comes from the integration of Newton's equations. Deriving the Newtonian dynamical equations for polarization and strain (and including terms of viscosity Àγ P _ P and Àγ η _ η), one gets: We then use the Fourier transformsP,η, andM of those quantities and obtain an equation for the polarization as a function of magnetization (see Supplementary Information for details): We now derive Eq. (4) with respect to a magnetic field H(ω 0 ) having angular frequency ω 0 and then derive with respect to a second magnetic field H ω 0 0 À Á having another (or the same) angular frequency ω 0 0 , and consequently apply the following small field-limit condition where β ω 0 0 ; ω 0 À Á is the dynamical (general) quadratic magnetic susceptibility. The full expression and the explicit derivations of the quadratic ME constants are detailed in the Supplementary Information. To simplify the analysis, we will consider, as in BiFeO 3 (BFO), that the linear ME constant is negligible as compared with quadratic ME terms 26 . We can then express the dynamical quadratic ME coefficient as We note that in Eq. 6, the term proportional to ξ PM , representing the direct coupling between polarization and magnetization, is exactly the term derived by Laguta et al. 27 . In a case where the applied field is a combination of an ac field (with frequency ω) superimposed on a dc bias field, the quadratic ME effect will generate periodic changes of the polarization with frequency ω (associated with the case ω 0 = ω and ω 0 0 ¼ 0, i.e., a coupling of the ac and dc field) and 2ω (corresponding to the situation which are proportional to the following ME constants The first component, β(0, ω), is the quadratic ME coefficient, which couples the dc and ac fields, and generates a frequency at ω for the polarization. This component can be considered as an effective linear ME term (as the linear term also generates a signal at the frequency ω) and experiences resonances under the following circumstances, according to Eq. (7): (i) χ M (ω) has a resonance, i.e., ω ≈ ω magnon , (ii) the frequency of the applied field corresponds to an optical phonon frequency (ω ¼ω P ), or (iii) the frequency of applied field corresponds to a strain resonance mode (ω = ω η ), the latter being shown in Fig. 1a.  Fig. 1 Mixing of long-wavelength acoustic phonons, optical phonons, and magnons. Polarization P is usually manipulated by electric fields E. Yet, direct coupling between optical polar phonons and magnetic moments and their quasiparticle-like excitation (called magnon) allows to control polarization with magnetic field H(ω). Such coupled electric and magnetic excitations are called electromagnons. Another pathway to control the polarization with magnetic fields is to use the mediation of acoustic phonon through magnetic moments-strain couplings (magnetostriction and piezomagneticity) and strain-polarization couplings (electrostriction and piezoelectricity). This strain-mediated magnetoelectric coupling involves a combined magnetic, mechanic, and electric excitation called an electroacoustic magnon 14 . In the case of a, the combined application of a static and dynamical magnetic field with frequency ω, the output change in polarization has frequency ω, and shows strain-induced resonance when the frequency of the dynamical field matches the strain resonance frequency ω η . When a dynamical field H(ω) is applied, a signal at 2ω can also be observed (b) and resonance occurs when the condiction 2ω = ω η is satisfied. At last, in the general case (c), two dynamical fields with frequencies ω 0 and ω 0 0 can be concurrently applied to yield a polarization oscillation with frequency ω 0 þ ω 0 0 ; in this latter case, resonance can be induced when The second component, β(ω, ω), is the quadratic ME coefficient that couples the ac field with itself and produces a frequency at 2ω for the polarization (second harmonic generation). This component also experiences resonances under three circumstances, according to Eq. (8): (i) χ M (ω) has a resonance, i.e., ω ≈ ω magnon -exactly as for the first case resonance of β(0, ω); (ii) the frequency of the applied magnetic field corresponds to half of the phonon frequency (ω ¼ω P 2 ); or (iii) the frequency of the applied field is half of the strain frequency (ω ¼ω η 2 ). It is noteworthy that, contrary to the case of β(0, ω) for which the phonon and strain resonances happen atω P andω η , in the case of β(ω, ω) the phonon-and strain-related resonances take place at ω ¼ω P 2 and ω 0 ¼ω η 2 , respectively, which is depicted in Fig.  1b for the latter. It is also worth noting that neither β(ω, ω) nor the more general β ω 0 0 ; ω 0 À Á of Eq. (6) were discussed in ref. 14 , despite their potential importance for nonlinear physics.
Our focus here is primarily on the quadratic response, β. Note that, due to the combined presence of a DC and an AC magnetic field, the quadratic term β(0, ω) is actually an effective linear ME term from a dynamical perspective (i.e., because of the existence of the dc field which results in the breaking of the magnetic symmetry, we have a ME effect that is linear in the ac field). As shown in the Supplementary Information, we also checked that the linear coefficient, α, does not represent a significant part of the effective linear term β(0, ω).
The terms involving strain-induced resonances are proportional toQ Pλ in Eqs. (7) and (8), which we recognize as the product of electro-mechanical and magneto-mechanical couplings. It then becomes clear that the electroacoustic magnon, discussed in ref. 14 (but for which the origin was not clearly explained), is presently revealed as a mixing of long-wavelength acoustic phonons, optical phonons, and magnons, as it results from the combined operation of piezoelectric and magnetostrictive couplings (which we depict in Fig. 1). This new additional handle provides a different path, as depicted in Fig. 1, to enhance the ME effect at variable frequencies as compared with "pure" electromagnons (understood here as resulting from the direct coupling of electrical and magnetic excitations). In practice, one can tailor the mechanical resonance frequency associated with the electroacoustic magnon by adapting the shape of the material, as commonly done in piezoelectric-induced dielectric resonances 28 .
To verify the validity of our model, we performed molecular dynamics (MD) simulations using an atomistic effective Hamiltonian (see "Methods" section for details) on a prototypical multiferroic material, i.e., BFO [29][30][31][32][33] . We neglect the inhomogeneous electric field arising from the time-dependent magnetic field for two main reasons: (1) to avoid the polarization response generated by the dielectric response of BFO to the electric field, i.e., our results truly reflect the ME nature of the response (in other words, we neglect this E-field, as we are solely focusing on the ME response); (2) we use periodic supercells in all three Cartesian directions that practically prevents us from modeling such inhomogeneous electric field. In addition, we assume an arbitrary (large but not infinite) value for the mass of strains and this mass is naturally linked with the frequencies of the electroacoustic magnons. The mass of the strains can be affected by varying the size and shape of the sample to induce a piezoelectric resonance and, therefore, adjust the electroacoustic magnon frequencies. We ran two types of simulations as follows: (1) unclamped simulations for which the homogeneous strain is allowed to change and (2) clamped simulations, with frozen supercell lattice vectors. The temporal behaviors of polarization, magnetization, and the components of strain under the influence of magnetic fields were then studied by Fourier analysis. The results of such analyses are shown in Figs. 2 and 3, when the frequency of the applied ac field is 160 GHz. More precisely, Fig. 2a shows the Fourier transform of the polarization when the system is clamped and unclamped. In both cases, one observes two peaks at 160 GHz and 320 GHz, which are the first and second harmonic responses to the applied magnetic field, respectively, as consistent with the fact that BFO is an ME material. Other peaks also occur, in both of these cases, at 4320 and 7008 GHz, which are optical phonon excitations impacting the polarization evolution. Interestingly, those latter peaks can also be seen in the Fourier transform of magnetization (see Fig. 2b for clamped and unclamped systems) and are thus in fact electromagnons 34 . Interestingly, when comparing the clamped and unclamped response of the polarization (Fig. 2a), two additional Fourier components can be noticed at 90 and 267 GHz in the second case only. The Fourier analysis also reveals that the peaks at 267 and 90 GHz are respectively associated with the vibrations of diagonal and shear components of homogeneous strain (see Fig. 3). As reckoned in ref. 14 , those peaks are thus linked to acoustic phonon excitation (i.e., strain resonances). There is also a peak at 100 GHz in Fig. 3, which is a mode that is observed at any applied frequency of magnetic field and even under no applied magnetic field. It is intrinsic to the strain, but does not couple with magnetization, polarization, and consequently with the ME response, β. The two additional frequencies at 267 and 90 GHz are also present in the Fourier decomposition of the magnetization response in Fig. 2b (a) (b) Fig. 2 Fourier analyses of polarization and magnetization. a Fourier analyses of polarization when homogeneous strain is clamped and unclamped. b Fourier analyses of magnetization when homogeneous strain is clamped and unclamped (for all cases, the frequency of the applied ac magnetic field is 160 GHz).
for the unclamped case, therefore showing that the overall response of the system at this frequency is a mixing of electric, magnetic, and strain response, hence the name electroacoustic magnon given in ref. 14 .
Furthermore, Fig. 4a and Fig. 4b show the results of our MD calculations for the value of the first and second components of ME coefficient β(0,ω) and β(ω,ω). In agreement with our aforementioned analytical model, we do observe resonances in β(0,ω) when the strain is unclamped in the simulations and when the frequency of the ac field corresponds to a resonance of the strain (i.e., 90 or 267 GHz). It is noteworthy that the numerically calculated values obtained in our simulations at small frequencies are in agreement with previously reported measurements for β(0, ω) 35 , which attests of the accuracy of our effective Hamiltonian for BFO.
Moreover, as it was discussed and predicted earlier in our phenomenological model for the case of β(ω,ω) and as numerically confirmed by Fig. 4 when the strain is unclamped, two strong strain-mediated resonances now take place at 45 and 133 GHz (i.e., precisely half frequencies with respect to the straininduced resonances of β(0,ω) in he second component of ME response. We also note that resonance frequencies are visible at the strain resonance frequencies of β(0,ω) (i.e., 90 GHz and 267 GHz) in β(ω,ω) too, albeit with smaller intensity. The detailed model in Eq. 40 of the Supplementary Information predicts that such latter resonances must happen in β(ω,ω) and are proportional to the linear ME coupling. It is thus natural to also observe a smaller intensity compared with strain half-frequency resonances in β(ω,ω), as such linear ME coupling is rather weak in BFO 26 .
Strain-mediated dynamical resonances in both components of the quadratic ME coefficients β(0, ω) and β(ω, ω), which we described at first using a simple analytical model, and then verified using more complex and complete numerical simulations, are thus a promising route to enhance ME couplings. Indeed, as the homogenous strain can vary depending on the size and shape of a sample, the resonances in ME responses can change accordingly 36 , and thus be tailored to operate at a specific speed. The dependence of the resonance frequencies on strain creates a potential to engineer samples with enhanced ME responses to be used towards fast (by tuning the resonance frequency) and efficient ME devices 37 . They also give rise to fascinating new physics, as the indirect couplings between strain and electrical and magnetic degrees of freedom allows the emergence of a new type of quasiparticle, i.e., the electroacoustic magnon.

METHODS
MD simulations within the framework of an effective Hamiltonian, which was developed for BFO 38 , were performed to study the quadratic ME coefficient, β, as well as the temporal behavior of electrical polarization and magnetization under the influence of combined dc and ac magnetic fields with various frequencies. The magnetic arrangement considered is Gantiferromagnetic; in other words, we assume the cycloid not to be present, which is the case at large magnetic fields 39 or under strain 40 . A 12 × 12 × 12 supercell with perovskite unit cells (8640 atoms) was subjected to MD simulations with periodic boundary conditions. The degrees of freedom that are considered in the total energy of such effective Hamiltonian scheme include the following:    Fig. 4 Real and imaginary parts of the first and second components of the quadratic magnetoelectric coefficien. a Real part of the first and second components of the quadratic magnetoelectric coefficient, β(0,ω) and β(ω,ω), when the homogeneous strain is clamped and unclamped. a Imaginary part of the first and second components of the quadratic magnetoelectric coefficient, β(0,ω) and β(ω,ω), when the homogeneous strain is clamped and unclamped.