Dynamic multiferroicity of a ferroelectric quantum critical point
K. Dunnett, J.-X. Zhu, N. A. Spaldin, V. Juricic, A. V. Balatsky
DDynamic multiferroicity of a ferroelectric quantum critical point
K. Dunnett, J.-X. Zhu, N. A. Spaldin, V. Juriˇci´c, and A. V. Balatsky
1, 4 Nordita, KTH Royal Institute of Technology and Stockholm University,Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden T-4 and CINT, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA Materials Theory, ETH Zurich, Wolfgang-Pauli-Strasse 27, CH-8093 Z¨urich, Switzerland Department of Physics, University of Connecticut, Storrs, CT 06269, USA (Dated: February 14, 2019)Quantum matter hosts a large variety of phases, some coexisting, some competing; when twoor more orders occur together, they are often entangled and cannot be separated. Dynamicalmultiferroicity, where fluctuations of electric dipoles lead to magnetisation, is an example wherethe two orders are impossible to disentangle. Here we demonstrate elevated magnetic response of aferroelectric near the ferroelectric quantum critical point (FE QCP) since magnetic fluctuations areentangled with ferroelectric fluctuations. We thus suggest that any ferroelectric quantum criticalpoint is an inherent multiferroic quantum critical point. We calculate the magnetic susceptibilitynear the FE QCP and find a region with enhanced magnetic signatures near the FE QCP, andcontrolled by the tuning parameter of the ferroelectric phase. The effect is small but observable -we propose quantum paraelectric strontium titanate as a candidate material where the magnitudeof the induced magnetic moments can be ∼ × − µ B per unit cell near the FE QCP. Quantum matter exhibits a plethora of novel phasesand effects upon driving [1], one of which is the strongconnection between the quantum critical point (QCP) ofone order parameter and the presence of another phase.The discussion has often focussed on the relation betweensuperconductivity and one or more magnetic phases [2–4]. However, other fluctuation-driven phase transitions,for example nematic phases in iron-based superconduc-tors [3, 5], have also received significant attention. We fo-cus here on the ferroelectric (FE) QCP which is a key partof the discussion of FE behaviour, particularly in dis-placive quantum paraelectrics [6, 7]. The behaviours thatmay occur near or as a result of such an FE QCP havebeen explored in various contexts [6–13], and the list ofsystems where the effects of quantum fluctuations can beobserved is expanding, with temperatures up to ∼ m can be induced by time-dependent oscillationsof electric dipole moments p : m = λ p × ∂ t p = C n × ∂ t n . (1)For magnetic moments to be induced, p has to exhibittransverse fluctuations; we therefore focus on rotationaldegrees of freedom of electric dipole moments [22]. Theunit direction vector of the constant amplitude electricdipole moment is n ≡ n ( r , t ), with time derivative ∂ t n , and C = λ | p | in terms of the electric dipole moments p (we use estimates from uniform polarisation P = | p | V with volume V in FE phases), and coupling λ = π/e ,with e the electric charge. Generally, we expect that or-ders entangled with the underlying static order can beexcited dynamically. One possibility is to use externaldriving mechanisms such as light, magnetic field or lat-tice strain to induce transient excitations of the entangledorders [2]. The present work addresses the complemen-tary case where inherent FE quantum fluctuations induceentangled ferromagnetic order fluctuations without anyexternal drive.In this Letter, we demonstrate that i) the fluctuatingdipoles can induce magnetic fluctuations that surroundthe FE QCP, as shown in Fig. 1a. The mechanism forthis effect is the induction of magnetic moments by fluc-tuating electric dipoles, described by Eq. (1), near theFE QCP and therefore describes inherent dynamic mul-tiferroicity. We support this scenario by calculating themagnetic susceptibility that, as we show, diverges in theparaelectric phase (Figs. 2 and 3), indicating a transitionto a new regime, labelled ‘Multiferroic PE’ in Fig. 1a. Wethus surmise that any FE QCP is a multicritical multifer-roic (MF) QCP with elevated magnetic fluctuations. Westress that the proposed effect is not due to permanentintrinsic magnetic moments, for example from unpairedelectrons on ions, but arises solely due to dynamics of theferroelectric order. While the proposed effect is general,we will consider the specific implications for magnetismin strontium titanate (STO) and provide estimates rel-evant to STO. ii) Within the approximations used, theapplication of a magnetic field B does not introduce astatic, B -dependent mass term to the effective action for p , and the position of the FE QCP is therefore indepen-dent of B . The Zeeman splitting of the FE active phonon a r X i v : . [ c ond - m a t . s t r- e l ] F e b FIG. 1. (a) Phase diagram near a ferroelectric quantum criti-cal point (QCP, where x = x cr ) with the magnetic susceptibil-ity (red line, dashed in FE phase) at ω = 0 . ω which divergesat the vertical dashed line, leading to a new ‘Multiferroic PE’phase. The ferroelectric quantum critical region (pale green)is now a dynamical multiferroic quantum critical region. Inboth the PE (white background) and FE (yellow background)phases, qualitatively similar behaviours of χ m are expected,despite the different underlying orders. The blue shading indi-cates the main regions where induced magnetic signatures areexpected: a dome around the FE QCP due to quantum fluc-tuations; a narrow band around the finite temperature phasetransition line due to thermally induced fluctuations. (b) Asimple experiment using a SQUID could detect magnetic sig-natures resulting from rotating electric dipoles in a systemtowards its ferroelectric phase transition. Here, the electricdipoles are constrained to the horizontal plane and lead to anout-of-plane magnetic moment m and susceptibility. modes [15, 22] meanwhile does affect the magnetic sus-ceptibility χ m and, in higher order approximations, isexpected to lead to a B term in the free energy, shift-ing the FE QCP. iii) We estimate the typical inducedmagnetic moment from a single rotating electric dipoleto be | m | ≈ × − µ B , where µ B is the Bohr magneton.This is for coupling λ = π/e and a dipole with charge 4 e and length 1 × − ˚A, rotating with frequency 0 . § I [30]). The overall contribution of thefluctuating FE order is diamagnetism near the FE QCP.
Model:
The system considered consists of fluctuatingelectric dipoles close to the PE-FE phase transition, in-ducing a magnetic moment via Eq. (1). In the absenceof external fields, the generic description of the system of rotating electric dipoles consists of the paraelectricphase: L P E = (cid:0) ω − ω q (cid:1) p ω,q p − ω, − q , where ω q is thedispersion of the phonon mode relevant for ferroelectric-ity, p ω,q is the rotating electric dipole moment writtenin Fourier (energy ω - momentum q ) space. The para-electric phase has negligible intrinsic magnetic contribu-tion and we therefore ignore intrinsic magnetisation alto-gether. However, the dynamic induction of m , Eq. (1),will lead to magnetic susceptibility of the paraelectricnear the FE QCP.The interaction between induced magnetic momentscan be neglected in the PE phase since the lowest ordercontribution | m | ∝ | p | . We assume optical phonons,relevant for the PE-FE transition in STO [27], with dis-persion ω q given by: ω q = ω (cid:18) − xx cr (cid:19) + bq = ω δ x + bq , (2)where δ x = x/x cr describes the distance to the ferroelec-tric QCP at x cr ; ω is energy at the zero momentum ofthe soft mode when x = 0, i.e. with no driving of thesystem towards the FE QCP, and q is the momentum. x is a tuning parameter that controls the paraelectric-ferroelectric phase transition at zero temperature, suchas doping. If the system is very close to the FE QCP, themomentum dependence is negligible and a flat dispersionwith b = 0 can be used. The system is paraelectric for δ x > δ x < p are present near the FE QCP, we willignore the amplitude fluctuations, so the time depen-dence of p is contained entirely in the unit direction vec-tor n . In this model, at the boundary between the PEand FE phases instead of | p | →
0, the dipoles rotate.That is: in the PE phase, finite-sized electric dipoles arepresent, but not aligned so the net polarisation is zero,and in the FE phase the dipoles align. n is linearisedas: n = n + ˜n ( t ) with ∂ t n = 0 and (cid:104) ˜ n (cid:105) = 0. Thezero-temperature Green’s function of the n field in ω − q space reads (cid:104) ˜ n jω ˜ n m − ω (cid:105) = A j δ jm G ( iω, q ) , (3)with A j as a constant factor. To find dynamic sus-ceptibilities, we use the retarded Green’s function G R ,obtained by analytical continuation to real frequencies( iω → ω + iη , SM, § II [30]): G R ( ω, q ) = Re (cid:18) ω q − ω (cid:19) + iπ ω q [ δ ( ω q − ω ) − δ ( ω q + ω )] . (4)We now calculate the magnetic susceptibility χ m in thePE phase: χ m = (cid:104) m ( r , t ) m ( r , t ) (cid:105) ≡ χ (1) + χ (2) (5)where m is given by Eq. (1). The two contributions to χ m are χ (1) ∝ (cid:104) ˜ n k ˜ n n (cid:105) and χ (2) ∝ (cid:104) ˜ n j ˜ n k ˜ n m ˜ n n (cid:105) .The quadratic contribution in ω − q space is: χ (1) il = C n j n m A k (cid:15) ijk (cid:15) lmk ω G R ( ω ) , (6)with G R ( ω ) given by Eq. (4), the factor C = λ V P gives the size of the magnetic susceptibility in terms ofthe coupling λ for the induced magnetic moments, andthe polarisation P of a sample of volume V . n i are thecomponents of n around which the fluctuations are ex-panded. The factor ω comes from the Fourier transformof (cid:104) ∂ t ˜ n k ∂ t ˜ n n (cid:105) .The quartic contribution to the magnetic susceptibilitycorresponds to a one-loop diagram as discussed in theSM, § III [30], with the real part given byRe[ χ (2) ii ] = − C δ il A j A k Λ πω x f ( ω ) (7)where f ( ω ), given in full in the SM, § III [30], contains δ -functions at 2 ω √ δ x ± ω and ω with weights ω or ω √ δ x ,and Λ is a momentum cut off. The imaginary part is:Im[ χ (2) ii ] = C δ il A j A k Λ π ( ω − ω x ) (cid:26) ω − ω x ω − ω x (cid:27) . (8)If the energy ω is written in terms of the q = 0 phononenergy ω , the size of the χ (2) contribution is determinedby Λ /ω . In STO, areas of coherent fluctuations are lim-ited to tetragonal domains, ∼ µ m [31], in which caseΛ /ω ∼ × , for ω = 0 . Results:
The total magnetic susceptibility χ m from Eq.(5) is plotted in Figs. 2 and 3 with the overall scale givenby the shared prefactor C = λ V P set to unity in allplots. In STO samples, the value of C can be estimatedfrom experimental data of samples tuned through the FEphase transition by applied strain or O isotope substi-tution, which indicates the possible size of the dipole mo-ments in the PE phase: C ∼ × − C m , for bulkSTO crystals, and C ∼ × − C m for 10 µ m tetrag-onal domains [23–25].Considering a sample with a singleinduced magnetic moment of 5 × − µ B per unit celland a sample volume of 1 µm , smaller than the tetrago-nal domains in STO, gives a sample magnetic moment of8 × µ B , well within spin sensitivities of 200 µ B / √ Hzof current SQUIDs [34].We consider tuning towards the FE QCP at a con-stant energy (fixed ω/ω ) first. In Fig. 2a, far from theFE QCP, the system is dielectric with Re[ χ m ] > χ m divergesand changes sign at δ x = ( ω/ω ) ; this indicates a phasetransition into a region where magnetic signatures canbe expected. As the energy is decreased, the divergencemoves towards the FE QCP and the magnetic featuresare confined into a narrower range of the tuning param-eter. -100 0 100 200 (a) R e [ χ m ] ( λ V P ) ω/ω I m [ χ m ] ( x λ V P ) δ x FIG. 2. The total magnetic susceptibility in units of the com-mon prefactor C = λ V P , and with Λ /ω = 1 × , asa function of δ x at several energies. (a) the real part; (b) theimaginary part. The behaviour in the FE phase δ x < χ (2) prefactorΛ /ω and the individual contributions of χ (1) and χ (2) arediscussed in the SM, § IV [30].
There are two contributions to the peaks in the realpart of the susceptibility: one is from the poles in Re[ χ (1) ]resulting in the large derivative feature at δ x = ( ω/ω ) ;the other comes from the δ -functions in Im[ G R ] that leadto peaks in Re[ χ (2) ] at δ x = ( ω/ ω ) . After the initialdivergence, Re[ χ m ] is negative apart from the sharp peakat δ x = ( ω/ ω ) , below which it quickly reaches a con-stant negative value independent of ω .The imaginary part of χ m , plotted in Fig. 2b, alsodiverges as expected at the border of the magnetic re-gion. This is followed by a divergence at δ x = ( ω/ ω ) corresponding to the peaks originating from χ (2) in thereal part. In the limit of δ x → Im [ χ m ] reaches apositive value that depends on the energy ω considered.It is important to note that the details of both the realand imaginary parts once the magnetic phase transitionhas been passed, that is for δ < ( ω/ω ) , are determinedsolely by the higher order χ (2) contributions.Changing energy while at a fixed distance from the FEQCP is considered in Fig. 3. The peaks and divergencesat finite ω are exactly those seen in Fig. 2, with an extra,artificial, divergence of both the real and imaginary partsat ω = 0, originating from calculating χ (2) in the contin- -200 0 200 400 (a) R e [ χ m ] ( λ V P ) δ x I m [ χ m ] ( x10 λ V P ) ω/ω FIG. 3. Total magnetic susceptibility, in units of C = λ V P , and with Λ /ω = 1 × , as a function of ω/ω atseveral distances from the FE QCP. (a) real part; (b) imagi-nary part. The peaks in the imaginary part ( δ − functions plot-ted as Lorentzian functions) corresponding to the divergenceand sign change of the real part are weak for Λ /ω = 1 × used; their presence is more easily seen at smaller values ofΛ /ω (SM, § IV [30]). uum limit. At the lowest energies, Re [ χ m ] <
0, it thenincreases and diverges at ω/ω = √ δ x , thus signallingthe phase transition with magnetic signatures expectedabove a critical energy scale. We note that upon increas-ing δ x to move away from the QCP, the onset of the tran-sition moves to higher energy. The energy dependence ofthe size of the imaginary part is seen particularly clearlyin Fig. 3b. In the FE phase, we expect qualitativelysimilar features, despite the underlying FE order, due tofluctuations of the ordered dipoles.A magnetic field B , applied perpendicular to the planeof the rotating dipoles, will have two effects. Firstly, thephonon Zeeman effect splits the phonon modes with alinear dependence on B [15] and moves the divergence of χ m [which occurs at δ = ( ω/ω ) for B = 0]. Secondly,an additional term in the Lagrangian for the interactionof the induced magnetic moments with the B , B · m = λ B · ( p × ∂ t p ) [22], can be treated as a perturbation tothe paraelectric system. Calculating the correspondingsecond order diagram (SM, § V [30]) does not introduce astatic, B -dependent mass term, but may do so at higherorders. Experimental proposal:
Strontium titanate (STO) maybe a suitable candidate material for the observation ofmagnetic signatures on tuning towards the FE QCP be- cause of its incipient ferroelectric nature below ∼
35K and its quantum paraelectric nature below 4 K [35]where the zero-point motion of the soft transverse op-tical phonon mode is high enough to prevent ferroelec-tricity even at zero temperature [36]. In O substitutedSTO, ω q =0 ( T ) becomes constant below 4 −
10K depend-ing on the distance from the FE QCP [37–41]. Thus,rotating electric dipoles could be present over an appre-ciable temperature range. Additional flexibility existsbecause several methods are available for tuning STOtowards the FE QCP, such as Ca doping [42], O sub-stitution [24, 25, 43], strain or applied pressure [23, 28].A simple experimental set up, consisting of a super-conducting quantum interference device (SQUID) abovean STO sample, that may permit the observation of theregion of pronounced magnetic fluctuations is sketchedin Fig. 1b. Strain is a particularly flexible meansof tuning STO samples towards the FE QCP, and bi-axial strain in STO thin films can confine polarisa-tion to the plane perpendicular to the tetragonal c-axis, but does not unambiguously determine the polari-sation direction [44–47], a favourable condition for theobservation of the magnetic signatures proposed here.Although strained STO is considered here, other FEQCPs and tuning mechanisms could be studied, e.g.:Ca − x Pb x TiO [48], strained KTaO [49]. The quan-tum dipole phase of the triangular lattice Mott insulator κ − (BEDT − TTF) Hg(SCN) Br [50] may also exhibitmagnetic signatures of inherent dynamical multiferroic-ity. The crucial ingredient for inherent dynamical mul-tiferroicity is incipient ferroelectricity (or quantum para-electricity) and only weak anisotropy between at leasttwo in-plane polarisation directions, to allow the fluctu-ations to well-defined circulating ions.
Discussion:
Including the long range interactions be-tween electric dipoles, such as those resulting from twinboundaries between tetragonal domains with differentlyoriented c-axes [51, 52], would introduce off-diagonalterms to the Green’s function [47]. The immediate ef-fect is a non-zero average magnetisation (cid:104) M (cid:105) ∝ (cid:104) n × ∂ t n (cid:105) .Alongside this, the off-diagonal components of the dielec-tric susceptibility χ ije = (cid:104) p i p j (cid:105) ∝ (cid:104) ˜ n i ˜ n j (cid:105) would also benon-zero at the twin boundaries, leading to a finite Kerreffect [53]. Further, the motion of twin boundaries maybe a means to induce relevant fluctuations of the elec-tric dipoles [54]. Scanning SQUID measurements able toresolve the individual tetragonal domains would be re-quired to investigate the effects of domain structures onthe magnetic signals. Again, STO is a potential candi-date material since tetragonal domains form naturally oncooling through the antiferrodistortive structural phasetransition at 105 K and their distribution can be con-trolled by applied pressure [32].The situation examined here is distinct from that re-cently considered in the context of multiferroic criti-cality [12] and other systems where the quantum criti-cal points of two or more types of order can be tunedby the same or different parameters leading to a fanwhere the quantum fluctuations of both orders are im-portant [12, 55]. In our model, the magnetic order doesnot exist independently of the ferroelectric order, leadingto an FE quantum critical region that is surrounded by aregion of strong magnetic fluctuations. While distinctfrom the nematic phase transitions seen in iron pnic-tides [5, 56], the multiferroic paraelectric region is an-other realization of competing orders near a QCP. Theinteraction between the induced magnetic moments andan external magnetic field is expected to mostly affectthe nature of the FE phase transition, as discussed formagnetic phase transitions [57–59]. Conclusions:
We have expanded the framework of dy-namic multiferroicity [15], and predict strongly enhancedferromagnetic (FM) susceptibility in a paraelectric ma-terial near its FE QCP. The induced magnetic suscep-tibility diverges at a finite distance from the FE QCP.The predicted effect points to another way for entangledquantum orders to appear. On the approach to the FEQCP, the fluctuations of the entangled (FM) order areenhanced as the static FE order develops quantum fluctu-ations. We thus suggest that any FE QCP may be an in-herent multiferroic QCP with entangled ferroelectric and(much weaker but present) ferromagnetic fluctuations.We expect magnetic signatures of fluctuating dipoles tobe observable experimentally, e.g. in SQUID measure-ments and could lead to additional signatures in opticalKerr and Faraday effects. Our results are applicable toany ferroelectric-paraelectric transition including classi-cal transitions at finite temperatures, where the fluctu-ations will be confined to a narrow Ginzburg-Levanyukregion near the transition. The effect will become pro-nounced near the T = 0 QCP. Finally, to illustrate thisscenario, we have considered STO as a system that canbe tuned towards its FE QCP.We are grateful to G. Aeppli, J. Lashley, and I.Sochnikov for useful discussions. The work was sup-ported by the US DOE BES E3B7, by VILLUMFONDEN via the Centre of Excellence for Dirac Ma-terials (Grant No. 11744) and by The Knut and AliceWallenberg Foundation (2013.0096). [1] D. N. Basov, D. A. Averitt, and D. Hsieh, Nature Ma-terials , 1077 (2017).[2] P. Gegenwart, Q. Si, and F. Steglich, Nature Physics ,186 (2008).[3] T. Shibauchi, A. Carrington, and Y. Matsuda, AnnualReview of Condensed Matter Physics , 113 (2014).[4] C. M. Varma, Reports on Progress in Physics , 082501(2016).[5] R. M. Fernandes, A. V. Chubukov, and J. Schmalian,Nature Physics , 97 (2014). [6] S. E. Rowley, L. J. Spalek, R. P. Smith, M. P. M. Dean,M. Itoh, J. F. Scott, G. G. Lonzarich, and S. S. Saxena,Nature Physics , 367 (2014).[7] P. Chandra, G. G. Lonzarich, S. E. Rowley, and J. F.Scott, Reports on Progress in Physics , 112502 (2017).[8] L. P´alov´a, P. Chandra, and P. Coleman, Phys. Rev. B , 075101 (2009).[9] J. M. Edge, Y. Kedem, U. Aschauer, N. A. Spaldin, andA. V. Balatsky, Phys. Rev. Lett. , 247002 (2015).[10] S. Rowley, M. Hadjimichael, M. N. Ali, Y. C. Durmaz,J. C. Lashley, R. J. Cava, and J. F. Scott, Journal ofPhysics: Condensed Matter , 395901 (2015).[11] C. W. Rischau, X. Lin, C. P. Grams, D. Finck, S. Hams,J. Engelmayer, T. Lorenz, Y. Gallais, B. Fauqu´e, J. Hem-berger, and K. Behnia, Nature Physics , 643 (2017).[12] A. Narayan, A. Cano, A. V. Balatsky, and N. A. Spaldin,Nature materials , 1 (2018).[13] J. R. Arce-Gamboa and G. G. Guzm´an-Verri, Phys. Rev.Materials , 104804 (2018).[14] S. Horiuchi, K. Kobayashi, R. Kumai, N. Minami, F. Ka-gawa, and Y. Tokura, Nature Communications , 7469(2015).[15] D. M. Juraschek, M. Fechner, A. V. Balatsky, and N. A.Spaldin, Phys. Rev. Materials , 014401 (2017).[16] I. Dzyaloshinsky, Journal of Physics and Chemistry ofSolids , 241 (1958).[17] T. Moriya, Phys. Rev. Lett. , 228 (1960).[18] T. Moriya, Phys. Rev. , 91 (1960).[19] H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev.Lett. , 057205 (2005).[20] I. A. Sergienko and E. Dagotto, Phys. Rev. B , 094434(2006).[21] S.-W. Cheong and M. Mostovoy, Nature Materials , 13(2007).[22] I. E. Dzyaloshinskii and D. L. Mills, Philosophical Mag-azine , 2079 (2009).[23] W. J. Burke and R. J. Pressley, Solid State Communica-tions , 191 (1971).[24] M. Itoh, R. Wang, Y. Inaguma, T. Yamaguchi, Y.-J.Shan, and T. Nakamura, Phys. Rev. Lett. , 3540(1999).[25] R. Wang and M. Itoh, Phys. Rev. B , 174104 (2001).[26] W. A. Atkinson, P. Lafleur, and A. Raslan, Phys. Rev.B , 054107 (2017).[27] Y. Yamada and G. Shirane, Journal of the Physical So-ciety of Japan , 396 (1969).[28] H. Uwe and T. Sakudo, Phys. Rev. B , 271 (1976).[29] A. Yamanaka, M. Kataoka, Y. Inaba, K. Inoue,B. Hehlen, and E. Courtens, EPL (Europhysics Letters) , 688 (2000).[30] See Supplemental Material for further details. Containsadditional Refs. [60–62].[31] H. Noad, E. M. Spanton, K. C. Nowack, H. Inoue,M. Kim, T. A. Merz, C. Bell, Y. Hikita, R. Xu, W. Liu,A. Vailionis, H. Y. Hwang, and K. A. Moler, Phys. Rev.B , 174516 (2016).[32] J. Hemberger, M. Nicklas, R. Viana, P. Lunkenheimer,A. Loidl, and R. B¨ohmer, Journal of Physics: CondensedMatter , 4673 (1996).[33] M. Honig, J. A. Sulpizio, J. Drori, A. Joshua, E. Zeldov,and S. Ilani, Nature materials , 1112 (2013).[34] M. E. Huber, N. C. Koshnick, H. Bluhm, L. J.Archuleta, T. Azua, P. G. Bjrnsson, B. W. Gard-ner, S. T. Halloran, E. A. Lucero, and K. A. Moler, Review of Scientific Instruments , 053704 (2008),https://doi.org/10.1063/1.2932341.[35] K. A. M¨uller and H. Burkard, Phys. Rev. B , 3593(1979).[36] E. Courtens, B. Hehlen, G. Coddens, and B. Hennion,Physica B: Condensed Matter
219 & 220 , 577 (1996).[37] H. Vogt, Phys. Rev. B , 8046 (1995).[38] Y. Yamada, N. Todoroki, and S. Miyashita, Phys. Rev.B , 024103 (2004).[39] M. Takesada, M. Itoh, and T. Yagi, Phys. Rev. Lett. ,227602 (2006).[40] H. Taniguchi, M. Itoh, and T. Yagi, Phys. Rev. Lett. ,017602 (2007).[41] S. E. M. Tchouobiap, Physica Scripta , 025702 (2014).[42] J. G. Bednorz and K. A. M¨uller, Phys. Rev. Lett. ,2289 (1984).[43] M. Itoh and R. Wang, Applied Physics Letters , 221(2000).[44] N. A. Pertsev, A. K. Tagantsev, and N. Setter, Phys.Rev. B , R825 (2000).[45] L. Ni, Y. Liu, C. Song, W. Wang, G. Han, and Y. Ge,Physica B: Condensed Matter , 4145 (2011).[46] F. Sun, H. Khassaf, and S. P. Alpay, Journal of MaterialsScience , 5978 (2014).[47] R. Roussev and A. J. Millis, Phys. Rev. B , 014105(2003).[48] V. V. Lemanov, A. V. Sotnikov, E. P. Smirnova, andM. Weihnacht, Applied Physics Letters , 886 (2002).[49] M. Tyunina, J. Narkilahti, M. Plekh, R. Oja, R. M.Nieminen, A. Dejneka, and V. Trepakov, Phys. Rev.Lett. , 227601 (2010).[50] N. Hassan, S. Cunningham, M. Mourigal, E. I. Zhilyaeva,S. A. Torunova, R. N. Lyubovskaya, J. A. Schlueter, andN. Drichko, Science , 1101 (2018).[51] H. J. H. Ma, S. Scharinger, S. W. Zeng, D. Kohlberger,M. Lange, A. St¨ohr, X. R. Wang, T. Venkatesan,R. Kleiner, J. F. Scott, J. M. D. Coey, D. Koelle, andAriando, Phys. Rev. Lett. , 257601 (2016).[52] Y. Frenkel, N. Haham, Y. Shperber, C. Bell, Y. Xie,Z. Chen, Y. Hikita, H. Y. Hwang, E. K. H. Salje, andB. Kalisky, Nature Materials , 1203 (2017).[53] B. Gu, S. Takahashi, and S. Maekawa, Phys. Rev. B ,214423 (2017).[54] R. T. Brierley and P. B. Littlewood, Phys. Rev. B ,184104 (2014).[55] C. Morice, P. Chandra, S. E. Rowley, G. Lonzarich, andS. S. Saxena, Phys. Rev. B , 245104 (2017).[56] R. M. Fernandes, A. V. Chubukov, J. Knolle, I. Eremin,and J. Schmalian, Phys. Rev. B , 024534 (2012).[57] J.-H. She, J. Zaanen, A. R. Bishop, and A. V. Balatsky,Phys. Rev. B , 165128 (2010).[58] U. Karahasanovic, F. Kr¨uger, and A. G. Green, Phys.Rev. B , 165111 (2012).[59] B. Liu and J. Hu, Phys. Rev. A , 013606 (2015).[60] D. J. Griffiths, Introduction to electrodynamics , 3rd ed.(Pearson Benjamin Cummings, San Francisco, 2008).[61] G. D. Mahan,
Many-particle Physics , 2nd ed., Physics ofsolids and liquids (Plenum Publishers, New York, 1990).[62] R. H. White,
Quantum Theory of Magnetism , Springerseries in solid-state sciences, Vol. 32 (Springer, BerlinHeidelberg, 2007).
SUPPLEMENTAL MATERIALI: Coupling strength
To calculate the strength of the coupling in the dy-namical multiferroic set up, we consider a current I flow-ing round a loop of radius r = a B (the Bohr radius)and a period τ B such that the energy (cid:126) /τ B is one Ryd-berg R y = (cid:126) / m e a B . The magnetic dipole moment isperpendicular to the plane of the loop and has magni-tude [60] m B = Iπr = eπa B τ B = e (cid:126) π m e = πµ B , (9)where µ B = e (cid:126) / m e is the Bohr magneton.The coupling strength λ is obtained by equating themagnitude of the magnetic dipole moment of the currentloop ( m B , above) with that of the dynamical multifer-roic set up [Eq. (1) of the main text] for electric dipolemoments of charge e and length a B rotating with period τ B : m dyn = λe a B τ B . (10)Requiring that these magnetic moments are equal, m B = m dyn , gives the coupling strength: λ = πµ B τ B e a B = πea B τ B τ B e a B = πe . (11)The ratio of the induced magnetic moment of any rotat-ing electric dipole to the Bohr model [Eq. (9)] is: m F E m B = ( n q ) ( n d ) τ B τ F (12)where τ F is the rotation period of the electric dipole(s), n q and n d are the charge and size of the dipoles in unitsof the electron charge and Bohr radius respectively. II: Model
The analytical continuation from the Matsubara fre-quencies to real frequencies consists of replacing iω inthe Green’s function by ω + iη [61]: G ( iω, q ) = 1 ω q − ( iω ) → G R ( ω, q ) G R ( ω, q ) = 1 ω q − ( ω + iη ) = 12 ω q (cid:20) ω q + ω + iη + 1 ω q − ω − iη (cid:21) (13) This is then evaluated using the principal value integralsto give: G R ( ω, q ) = 12 ω q (cid:40) ω q + ω + 1 ω q − ω + iπ [ δ ( ω q − ω ) − δ ( ω q + ω )] (cid:41) = 1 ω q − ω + iπ ω q (cid:2) δ ( ω q − ω ) − δ ( ω q + ω ) (cid:3) , (14)where ω q is the dispersion of the ferroelectric (FE)phonons: ω q = (cid:112) ω δ x + bq . Near the ferroelectricquantum critical point (FE QCP), the momentum de-pendence can be neglected [47]. III: Calculation of χ (2) The second contribution to the magnetic susceptibilityis: χ (2) m,il = C (cid:15) ijk (cid:15) lmn (cid:104) ˜ n j ( t ) ∂ t ˜ n k ( t )˜ n n ( t ) ∂ t ˜ n m ( t ) (cid:105) (15)where the temporal (and spatial) arguments have beenincluded explicitly (( t ) ⇒ ( r , t )), and C = λ V P .The easiest way to evaluate this is to recognise that theangular bracket corresponds to the loop of the diagram: A B ω, Q ω, Q Ω , qω − Ω , Q − qi jk lnm in which the wiggly lines represent magnetic propagatorswhile the plain curves with direction arrows are ferroelec-tric propagators. The equivalence of the internal linesmeans δ jm and δ kn and the (cid:15) ijk (cid:15) lmn prefactor becomes (cid:15) ijk (cid:15) ljk = + δ il d !. Integration over the internal degreesof freedom is the integral (cid:82) d Ω d d q/ (2 π ) d +1 which is tobe computed. The coupling constant λ has already beenfactored out so the factor at A is i Ω, and that at B is i ( ω − Ω) which we assume are independent of momen-tum. The internal lines give contributions of the (FEphonon) Green’s functions: A k G R (Ω) for the lower and A j G R ( ω − Ω) for the upper parts of the loop respectively.Evaluating the diagram corresponds to calculating theintegral: D jk = − A j A k (cid:90) d d q (2 π ) d d Ω2 π G R (Ω) G R ( ω − Ω)Ω( ω − Ω) , (16)with χ (2) il = C d ! δ il D jk .The Green’s function given by Eq. (14) has both real and imaginary parts, so G R (Ω) G R ( ω − Ω) leads to severalterms: G R (Ω) G R ( ω − Ω) = 1 ω q − Ω (cid:16) ω Q − q − ( ω − Ω) (cid:17) − π ω q ω Q − q (cid:104) δ ( ω q − Ω) δ ( ω Q − q − ω + Ω) − δ ( ω q + Ω) δ ( ω Q − q − ω + Ω) − δ ( ω q − Ω) δ ( ω Q − q + ω − Ω) + δ ( ω q + Ω) δ ( ω Q − q + ω − Ω) (cid:105) + iπ (cid:40) δ ( ω Q − q − ω + Ω) − δ ( ω Q − q + ω − Ω) ω Q − q ( ω q − Ω ) + δ ( ω q − Ω) − δ ( ω q + Ω) ω q ( ω Q − q − ( ω − Ω) ) (cid:41) . The integral over the internal energy is calculated first,then one considers that the momentum dependence ofthe phonon spectrum is irrelevant near the FE QCP so ω Q − q = ω q = ω √ δ x . Assuming spherical symmetry, theintegral over d d q becomes (cid:82) ∞ dqq in d = 3; evaluatingup to some cut-off value Λ introduces a Λ weight to χ (2) .The real part comes from the middle term of G R (Ω) G R ( ω − Ω):Re[ χ (2) ] = − C δ il A j A k Λ πω x f ( ω ) (17)with f ( ω ) = ω (cid:104) δ ( ω ) + δ (cid:16) ω (cid:112) δ x + ω (cid:17) − δ (cid:16) ω (cid:112) δ x − ω (cid:17)(cid:105) + ω (cid:112) δ x (cid:104) δ (cid:16) ω (cid:112) δ x − ω (cid:17) + δ (cid:16) ω (cid:112) δ x + ω (cid:17)(cid:105) . (18)The imaginary part, from the first (via Cauchy’s residuetheorem) and third ( δ − functions) terms of G R (Ω) G R ( ω − Ω), is given in the main text [Eq. (8)]. Dimensionalanalysis gives the dimensions of the A k factors in χ (2) ,and the n factor in χ (1) contains an implicit integral over dωd d q to ensure dimensional consistency. For simplicity, | A k | = 1 and | n | = 1 are used.In the limit of static dipoles ( ω = 0), both the realand imaginary parts of χ (1) are zero due to the ω fac-tor. Meanwhile, Re[ χ (2) ] diverges with an overall nega-tive factor, thus determining the static behaviour. Theimaginary part of χ (2) also diverges. At the FE QCPwhere δ x = 0, the real part of χ (1) is a negative constantand the imaginary part is zero. The real part of χ (2) goes to −∞ exactly at the FE QCP, but is zero on theapproach, and the contribution to Im[ χ m ] is a positive,energy dependent constant.In the opposite limit, ω → ∞ , the real part of χ (1) isa negative constant, while Re[ χ (2) ] diverges to −∞ . Theimaginary parts are Im[ χ (1) ] = 0 unless δ x = ∞ too andIm[ χ (2) ] = 0 as 1 /ω + 1 /ω . The system is well behaved in that the rate of energy absorption, as quantified byIm[ χ m ] is finite, even in the limit of infinite energies [62]. IV: χ (1) and χ (2) contributions The contribution from χ (2) depends on a momentum-dependent factor Λ /ω . As seen in Fig. 4, the regionsof positive χ m after the initial divergence at δ x = ω /ω are suppressed for small Λ /ω , corresponding to largedistances (or sample size) for a given phonon frequency ω . The observation of these features close to the FEQCP will depend strongly on the distances and energiesconsidered.The contributions of χ (1) , given by the Green’s func-tion, Eq. (14), and χ (2) , Eq. (17) here and Eq. (8) of themain text, are plotted as functions of δ x in Fig. 5 (a) and(b), and as functions of energy in Fig. 5 (c) and (d). Thishighlights the origin of the features seen in Figs. 2 and3 of the main text. In all plots, both here and the maintext, the δ -functions from Im[ G R ( ω )] have been replacedby Lorentzian functions. V: Behaviour in a magnetic field
The additional term in the Lagrangian describing theinteraction between a magnetic moment and an externalfield, L B = B · m = λ B · ( p × ∂ t p ), can be treated asa perturbative term in the full Lagrangian, the secondorder expansion of which gives the diagram A B ω, Q ω, Q Ω , qω − Ω , Q − q j ik lm n -200-100 0 100 (a) ℜ [ χ m ] ( λ V P ) Λ /ω I m [ χ m ] ( λ V P ) δ x Λ /ω χ FIG. 4. (a) Real part of χ m at ω/ω = 0 . /ω =1 × , × , × . Reducing the significance of the χ (2) contribution narrows the width of the peak in Re[ χ m ] nearthe FE QCP and ensures diamagetic behaviour (Re[ χ m ] > χ m at ω/ω = 0 . /ω = 1 × , × , × ; for the last 2 × − Im[ χ m ] is plotted:the position of the peak at δ x = (cid:112) ω/ω does not move, butit becomes less noticeable as Λ /ω is increased. Using the standard diagrammatic rules, this correspondsto the integral: I B = Lω (cid:90) d d q (2 π ) d d Ω2 π Ω G R (Ω , q ) G B ( ω − Ω , Q − q ) , (19)where the overall positive sign is from ( i ) (cid:15) ijk (cid:15) ikn , and L ≡ C A k A n d ! δ jn . G R is the usual Green’s functionfor the ferroelectric propagator, Eq. (14), and G B is themagnetic propagator: G B ( ω − ω (cid:48) , q − q (cid:48) ) ≡ (cid:104) B i ( − ω, − q ) | B l ( ω (cid:48) , q (cid:48) ) (cid:105) = B δ il δ ( ω − ω (cid:48) ) δ ( q − q (cid:48) ) . (20)In the diagram, have ω − Ω, and Q − q . Thus the calcu-lation to be performed is: I jnB = LωB (cid:90) d d q (2 π ) d d Ω2 π δ ( ω − Ω) δ ( Q − q )Ω ×× (cid:18) ω q − Ω + iπ ω q [ δ ( ω q − Ω) − δ ( ω q + Ω)] (cid:19) (21)The δ -functions from G B lead to evaluation of the inte-grals at Ω = ω and q = Q . The integral over energy isstraightforward, and after assuming spherical symmetry, d d q = dqq d − S d ; with S = 4 π , one finds I jnB = Lω Q B π G R ( ω, Q ) . (22)A magnetic field therefore affects the energy of the FEphonons, but does not, at this level of approximation,move the FE QCP.0 -20-10 0 10 20 30 40 (a) R e [ χ i ] χ (1) χ (2) -50050100150 01(b) I m [ χ i ] δ x χ (1) χ (2) (c) χ (1) χ (2) ω/ω χ (1) χ (2) FIG. 5. Real and imaginary parts of χ (1) (solid purple lines) and χ (2) (dashed green lines) contributions to χ m . (a), (b) as afunction of distance from the FE QCP ( δ x = 0) at fixed ω = ω ; scale of χ (2) is 5000. (c), (d) as a function of energy at fixed δ x = 0 .
4; scale of χ (2) is 1000, see Fig. 4 for the effect of changing weight of χ (2) on the total susceptibility. In all cases, thescale is in terms of the common size λ V P , and the finite width and height of the peaks in Re[ χ (2) ] and Im[ χ (1) ] are theresult of replacing the δ −−