# Gravitational-wave glitches in chaotic extreme-mass-ratio inspirals

Kyriakos Destounis<EMAIL_ADDRESS>Theoretical
Astrophysics, IAAT, University of Tübingen, Germany  Arthur G. Suvorov
Theoretical Astrophysics, IAAT, University of Tübingen, Germany  Kostas D.
Kokkotas  Theoretical Astrophysics, IAAT, University of Tübingen, Germany

###### Abstract

The Kerr geometry admits the Carter symmetry, which ensures that the geodesic
equations are integrable. It is shown that gravitational waveforms associated
with extreme-mass-ratio inspirals involving a non-integrable compact object
display ‘glitch’ phenomena, where the frequencies of gravitational waves
increase abruptly, when the orbit crosses certain spacetime regions known as
Birkhoff islands. The presence or absence of these features in data from
upcoming space-borne detectors will therefore allow not only for tests of
general relativity, but also of fundamental spacetime symmetries.

_Introduction._ One of the main targets for space-borne gravitational-wave
(GW) detectors, such as the Laser Interferometer Space Antenna (LISA) and
Taiji lisafreq ; taiji , are extreme-mass-ratio inspirals (EMRIs). These
binaries involve a (super-)massive black hole (BH) and a companion whose mass
is negligible (mass ratio $\lesssim 10^{-4}$) relative to the primary. This
latter aspect, which implies that the inspiral dynamics are well approximated
as the geodesic motion of the companion through the gravitational field
generated by the primary, even close to the final plunge ori00  (though see
also Ref. berry16 ), means that EMRIs provide unparalleled information about
the governing theory of gravity ryan95 ; glamp06 ; bara20 . In particular,
general relativity (GR) predicts that the Kerr metric uniquely describes the
gravitational field of astrophysically stable BHs rob75 . Therefore, if EMRIs
with significantly non-Kerr trajectories present themselves in the LISA or
Taiji data, this would represent a ‘smoking gun’ for the breakdown of GR
bert05 ; Barack:2006pq ; schutz09 .

One interesting possibility of a non-GR signature relates to the integrability
of the orbit apo09 ; luk ; cont11 ; luk14 ; glamp19 ; dest20 . The Kerr metric
is stationary and axisymmetric and therefore possesses two Killing vectors,
which lead to the conservation of energy and angular momentum for relativistic
particles traversing the gravitational field. Additionally, the Kerr metric is
sufficiently special, algebraically speaking, so as to admit a non-trivial
Killing tensor carter68 , which provides an additional constant of motion and
leads to the Liouville-integrability of the geodesic equations; Kerr geodesics
do not display chaotic phenomena wilkins72 . Distortions in the Hamiltonian
describing geodesic motion in the Kerr spacetime may either preserve the
integrability or not depending on their exact character, in accord with the
Kolmogorov-Arnold-Moser (KAM) theorem moser62 ; arnold63 .

Even when a non-integrable perturbation is introduced into the Kerr
Hamiltonian, some orbits necessarily remain periodic poin12 ; birk13 . From a
phase-space perspective however, when the perturbation is non-integrable,
small ‘islands’ of stability come to surround the surviving periodic orbits
arnold89 . These islands have the property that, as a particle passes through
them, the ratio of longitudinal and transverse frequencies associated with the
orbit remain constant Brink:2013nna ; brink13 . In this way, the orbital
dynamics are expected to display a transient plateau feature cont02  as
gradual radiation reaction can cause orbits to weave in and out of islands.
While this leads to clear signatures from a dynamical systems perspective
apo09 ; luk ; cont11 ; luk14 ; dest20 , we show here that occupancy in an
island leads to a noticeable ‘glitch’ (i.e., a rapid though short-lived
increase, akin to that which is observed in pulsars) in the GW frequency.
However, as observed in the LISA pathfinder mission path1 , instrumental noise
artifacts can also lead to similar phenomena in the data stream (see Fig. 2 in
Ref. ed20 ). The main contribution of this Letter is to demonstrate that
abrupt frequency increases during an EMRI may have a genuine astrophysical
origin related to the fundamental spacetime symmetries associated with
supermassive compact objects, and discarding them _a priori_ as instrumental
artifacts may miss crucial physics.

_Spacetime geometry._ Theoretical deviations from a Kerr description for dark,
compact objects can be broadly classified into two categories: those which
introduce generic, though integrable, deformations into the geodesic
Hamiltonian, and those which introduce non-integrable deformations vig11 ;
joh13 ; krz16 . To give a concrete example of how these features appear at the
level of the spacetime metric, we consider the geometry introduced in Ref.
dest20 . In Boyer-Lindquist coordinates $(t,r,\theta,\phi)$, the metric we
operate with reads

$\displaystyle ds^{2}=$
$\displaystyle-\frac{\Sigma[\left(\alpha_{Q}/r\right)M^{3}+\Delta-a^{2}A(r)^{2}\sin^{2}\theta]}{[(r^{2}+a^{2})-a^{2}A(r)\sin^{2}\theta]^{2}}dt^{2}$
(1)
$\displaystyle-\frac{2a[(r^{2}+a^{2})A(r)-\Delta]\Sigma\sin^{2}\theta}{[(r^{2}+a^{2})-a^{2}A(r)\sin^{2}\theta]^{2}}dtd\phi$
$\displaystyle+\frac{\left(\alpha_{Q}/r\right)M^{3}+\Sigma}{\Delta}dr^{2}+\Sigma
d\theta^{2}$
$\displaystyle+\frac{\Sigma\sin^{2}\theta\left[(r^{2}+a^{2})^{2}-a^{2}\Delta\sin^{2}\theta\right]}{[(r^{2}+a^{2})-a^{2}A(r)\sin^{2}\theta]^{2}}d\phi^{2},$

where $\Sigma=r^{2}+a^{2}\cos^{2}\theta$, $\Delta=r^{2}-2Mr+a^{2}$, and
$A(r)=1+r^{-2}\alpha_{22}M^{2}$ for mass $M$ and spin $a$. The line element
(1) resembles the Kerr geometry, though contains two additional parameters.
One of these, $\alpha_{22}$, was first considered by Johannsen joh13 , and
introduces a generic, though integrable, deformation from Kerr (henceforth
referred to as a ‘deformed-Kerr’ metric). The other, $\alpha_{Q}$, breaks the
Carter symmetry and leads to a non-integrable perturbation (‘non-Kerr’ metric)
dest20 . In the limit where these two latter parameters vanish, the Kerr
metric is recovered and the spacetime is an exact solution in GR and various
other theories of gravity psal08 . For a wide range of values of $\alpha_{22}$
and $\alpha_{Q}$, the metric described by (1) exists as an exact, vacuum
solution in some particular family of mixed scalar-$f(R)$ theories suv20 .

_Geodesic dynamics and non-integrability._ Noting that the metric components
$g_{\mu\nu}$ within (1) are functions of $r$ and $\theta$ only, so that the
spacetime is stationary and axisymmetric, Hamilton’s equations for a particle
of mass $\mu$ with momenta $p^{\mu}=\dot{x}^{\mu}$ imply that
$\dot{p}_{t}=0=\dot{p}_{\phi}$, where the overhead dot denotes differentiation
with respect to proper time. In general, there exist three constants of motion
associated with geodesic dynamics: the orbital energy, $E$, angular momentum,
$L_{z}$, and the particle mass itself, $\mu$. The Kerr spacetime
$(\alpha_{22}=\alpha_{Q}=0$), in addition to admitting the three linear
symmetries detailed above, also possesses a quadratic symmetry due to the
existence of a non-trivial rank-2 Killing tensor carter68 . This additional
integral of motion, giving rise to the Carter constant (which generally exists
even for $\alpha_{Q}=0,\,\alpha_{22}\neq 0$), implies the integrability of the
geodesic equations as a whole cont02 .

A generic orbit for a stationary spacetime can be characterized by two
libration-like frequencies, $\omega_{r}$ and $\omega_{\theta}$, which describe
the transition rate between periastron and apastron and longitudinal
oscillations around a given plane, respectively Brink:2013nna ; brink13 .
Those orbits for which $\omega_{r}$ is an integer multiple of
$\omega_{\theta}$, or vice versa, are called resonant. At the level of orbital
dynamics, the KAM and Poincaré-Birkhoff theorems moser62 ; arnold63  together
imply that half of the resonant orbits remain stable while half become
unstable when introducing a non-integrable ‘perturbation’ into the Kerr (or
any other regular) Hamiltonian. In the phase space, small islands of stability
(‘Birkhoff islands’) come to surround each of the stable orbits, while the
unstable orbits form chaotic layers which surround the islands. A key property
of these islands is that the ratio $\omega_{r}/\omega_{\theta}$, often called
the rotation number, remains constant there and forms a dynamical ‘plateau’
cont02 .

The quadratic symmetry described above is preserved for any value of
$\alpha_{22}$, however, and thus the equations of motion remain integrable. On
the other hand, for non-zero values of the parameter $\alpha_{Q}$ and spin,
the spacetime (1) does not admit the Carter symmetry. Therefore, for
astrophysical objects described by the line element (1) with $\alpha_{Q}\neq
0$, Birkhoff islands form within the phase space dest20 . Depending on the
exact value of the parameters $M,a,\alpha_{22}$, and $\alpha_{Q}$, the
location and width of the islands vary luk ; dest20 . More details on the
properties of the metric (1) and its relationship with existing astrophysical
constraints can be found in Ref. dest20 .

It is the purpose of this work to explore the theoretical manifestation of a
plateau at the level of the gravitational waveform. In particular, a spacetime
for which the geodesic equations are non-integrable will possess a scattered
series of islands. When a companion object passes through one of these islands
it undergoes an abrupt orbital evolution, which leads to a sudden jump in the
GW frequencies. The frequency jumps we observe are similar to that which is
sometimes seen in pulsar timing experiments, namely ‘glitches’. For better or
worse, jumps of this sort are also known to occur due to effects of
instrumental origin path1 ; ed20 . It is important therefore that the
phenomenology of non-integrability is better understood, as meaningful physics
may be discarded if it is _a prori_ assumed that frequency jumps are not of
astrophysical origin.

_Orbital evolution and the kludge scheme._ The orbital dynamics governing an
EMRI are not purely geodesic since GW emission leads to a dissipation of the
particle 4-momenta. In general, the equations describing the motion of a
particle in a spacetime are the Einstein equations (or some appropriate
generalization) together with the conversation laws
$T^{\alpha\beta}_{\phantom{\alpha\beta};\beta}=0$. Treating backreaction at
the linear level, the conservation laws lead to the MiSaTaQuWa equations misa1
; misa2 , which are known to be equivalent to the geodesic equations on a
_modified_ spacetime $\tilde{g}_{\mu\nu}=g_{\mu\nu}+h_{\mu\nu}^{R}$, where the
superscript $R$ stands for the regularized metric perturbation (see Sec. 19 of
Ref. pois  for a discussion). This latter term, although necessarily small and
satisfying $||h_{\mu\nu}^{R}||\ll||g_{\mu\nu}||$, is time-dependent and
therefore allows for the particle to potentially enter and leave Birkhoff
islands during its orbital lifetime.

In any case, the above shows that self-force sourced by backreaction can be
modeled by introducing time-dependent shifts into the momenta $p^{\alpha}$ of
the particle, which can be related to the specific energy and angular momentum
through $E=-p^{t}\partial_{t}$ and $L_{z}=p^{\phi}\partial_{\phi}$. This
observation forms the basis for the _adiabatic_ approximation introduced by
Mino mino  (see also Ref. gal ), valid when the change in any quantity that
characterizes the orbit (such as energy or angular momentum) is sufficiently
small over a single orbit sago . In this sense, taking an average of the
MiSaTaQuWa equations essentially leads to some evolution equations for $E$ and
$L_{z}$ and the remaining momenta. These equations are still, however,
relatively difficult to work with. A hybrid _kludge_ scheme kludge1 ; glamp07
can instead be used where the aforementioned equations are expanded up to some
desired post-Newtonian and multipolar orders.

In this Letter, we adopt the _kludge_ scheme described above to perform the
orbital evolutions. Up to second post-Newtonian order the relevant equations
for the evolution of an orbit are lengthy, though are given explicitly by
equations (37)–(39) in Ref. gair06 . These are the equations used here, though
with one important modification.

In particular, at second post-Newtonian order, the _kludge_ scheme involves
the mass quadrupole moment $M_{2}$ of the system. To construct the non-Kerr
inspiral we apply a modification to the mass quadrupole moment, similar to
that in Refs. Gair:2007kr ; luk , where
$M_{2}=-Ma^{2}-\frac{1}{3}M^{3}\alpha_{Q}$ for the metric (1). The parameter
$\alpha_{22}$ introduces current multiple moments that appear at higher-PN
order, and therefore does not explicitly enter into the _kludge_ equations of
motion. While the exact definition for the multipole moments can depend on the
theory under consideration papp15 , we have checked that our results are
qualitatively unchanged for different choices of $M_{2}$ (e.g.,
$M_{2}=M_{2}^{\text{Kerr}}$). We further linearize the adiabatic energy and
angular momentum evolution as in Ref. Canizares , viz.
$E_{1}=E_{0}+dE/dt|_{0}N_{r}\,T_{r}$ and
$L_{z,1}=L_{z,0}+dL_{z}/dt|_{0}N_{r}\,T_{r}$, where $E_{0},\,L_{z,0}$ are the
initial energy and angular momentum, respectively,
$dE/dt|_{0},\,dL_{z}/dt|_{0}$ are the radiation loss rates calculated using
the _kludge_ expressions at the beginning of the inspiral, respectively, and
$T_{r}$ is the time that the orbit takes to travel from the periastron to
apoastron and back. The aforementioned equations are updated every $N_{r}$
cycles for the whole EMRI evolution.

_Gravitational waves and frequency evolution._ We model the gravitational
waveform using the Einstein-quadrupole approximation to understand the main
phenomenology of the transient features acquired through passage of an island.
A more sophisticated approach would involve solving the (appropriately
generalized, see, e.g., Refs. suv19a ; suv19b ) Teukolsky equations directly
to deduce the GW characteristics through the Weyl scalars. However, such an
analysis is considerably more complicated and does not help to elucidate the
main new features presented in this work: the existence of ‘glitches’ from
chaotic inspirals.

In general, an incoming GW can be projected onto the mutually orthogonal $+$
and $\times$ polarization states by introducing two vectors
$\boldsymbol{p}=\boldsymbol{n}\times{\boldsymbol{Z}}/|\boldsymbol{n}\times{\boldsymbol{Z}}|$
and $\boldsymbol{q}=\boldsymbol{p}\times\boldsymbol{n}$, which are defined in
terms of a unit vector $\boldsymbol{n}$ which points from the source to the
detector for particle position $Z^{i}(t)$. Note that some modified theories of
gravity predict the existence of up to six polarization states, a detection of
which would also conclusively signal the breakdown of GR liu20 . These extra
polarization states are not, however, tied to the symmetries of the spacetime,
since even a Kerr black hole can emit scalar-like waves in an $f(R)$ theory,
for instance suv19a , and we assume they are negligible. At a luminosity
distance $d$ from the source and in the quadrupole approximation, the GW
amplitudes read
$h_{+,\times}(t)=2\mu\epsilon^{+,\times}_{ij}\left[a^{i}(t)Z^{j}(t)+v^{i}(t)v^{j}(t)\right]/d$
for velocity, $v^{i}(t)=dZ^{i}/dt$, and acceleration
$a^{i}(t)=d^{2}Z^{i}/dt^{2}$ Canizares , where the $+$ and $\times$
polarization tensors have components
$\epsilon_{+}^{ij}=p^{i}p^{j}-q^{i}q^{j}$, and
$\epsilon_{\times}^{ij}=p^{i}q^{j}+p^{j}q^{i}$, respectively.

The first of the deformation parameters introduced in (1), $\alpha_{22}$, has
the strongest effect on observables when the object is near the extremal limit
$a\lesssim M$ joh13 . Since accretion torques are expected to spin-up super-
massive BHs that exist within galactic centers thorne74 , even a small value
of $\alpha_{22}$ can lead to significant modifications in the orbital dynamics
dest20 . We thus consider a near-extremal object with $a=0.99M$ throughout.

Figure 1: Power spectral density of a Kerr (green), a non-Kerr (red) and a
deformed-Kerr (blue) EMRI consisting of a light companion with $\mu=M_{\odot}$
and a supermassive compact object with $M=10^{6}M_{\odot}$ ($\mu/M=10^{-6}$)
and $a=0.99M$ at a luminosity distance $d=50\text{ Mpc}$, over-plotted on the
LISA power spectral density curve with $\text{SNR}=1$ (dashed curve). All
small mass companions begin their inspiral with $E_{0}=0.95\mu$,
$L_{z,0}=3M\mu$ and initial conditions
$r(0)=5.70903M,\,\dot{r}(0)=0.1,\,\theta(0)=\pi/2$, with the remaining
$\dot{\theta}(0)$ being defined by the constraint equation (7) in dest20 . The
most prominent Fourier peaks are presented in the zoomed inlays, while the
remaining peaks correspond to higher harmonics. The evolution time of the EMRI
is $t=2\times 10^{6}M$, which corresponds, roughly, to $3$ months of
detection. The deformation parameters are $\alpha_{Q}=-1.8$ for the non-Kerr
and $\alpha_{22}=-2.2$ for the deformed-Kerr inspiral, while for the Kerr
inspiral $\alpha_{Q}=\alpha_{22}=0$. Figure 2: Periodogram of a deformed-Kerr
(left) and a non-Kerr EMRI (right) GW, plotted below the most prominent GW
peak, as depicted in the right inlay of Fig. 1.

Although the LISA data stream consists of two linearly independent channels,
with each one being better suited to the detection of particular types of
signals cutler98 , we will work, for simplicity, within the single channel
approximation and neglect any noise in the data stream to demonstrate the main
features of detection. In what follows, we compare frequency evolution and
power spectral densities (PSDs) of two EMRIs, consisting of a small-mass
companion (a stellar mass BH or a neutron star) with $\mu=M_{\odot}$ and a
supermassive compact object with $M=10^{6}M_{\odot}$ ($\mu/M=10^{-6}$); the
first evolves on a non-Kerr spacetime and crosses a Birkhoff island of
$2/3$-resonance, while the second evolves on a deformed-Kerr spacetime and
passes through a $2/3$-resonance, though no Birkhoff islands are present
anywhere in the phase space. The detectability of the PSD of such EMRIs is
demonstrated in Fig. 1, where we assume a source at $d=50\text{ Mpc}$, where
the LISA PSD is also plotted with signal-to-noise ratio (SNR) equal to unity.
For completeness, we further include the PSD of a Kerr EMRI. The most
prominent peaks, including various higher harmonics, appear in the range of
maximal sensitivity of LISA, between $10^{-3}-10^{-2}$ Hz. Although the PSD
peaks of Kerr and deformed-Kerr EMRIs, which share the same fundamental
spacetime symmetries, have smooth profiles (see inlays in Fig. 1), non-Kerr
EMRIs, which cross Birkhoff islands, will typically undergo a frequency
modulation during the crossing, where the PSD peaks (and the amplitude of the
Fourier spectrum of the GW) abruptly decrease by up to $2$ orders of magnitude
(see right inlay in Fig. 1). This substantial drop in amplitude should lead to
a significant effect in the frequency evolution of a non-Kerr EMRI.

In Fig. 2 we show the frequency evolution of the most prominent spectral peaks
from a deformed-Kerr and a non-Kerr EMRI with the same initial conditions (we
have chosen the deformation parameters for each EMRI so that the initial
rotation numbers of the orbits agree to $\sim 0.1\%$). When the underlying
spacetime symmetry is similar to that of Kerr, the frequency evolution follows
a linear growth as the object inspirals towards its supermassive companion.
Although the deformed-Kerr EMRI passes through a $2/3$-resonance, no obvious
effect is present. On the other hand, when a non-Kerr EMRI crosses through a
Birkhoff island of $2/3$-resonance a ‘glitch’ interrupts the linear evolution
of the GW frequencies. Although we only demonstrate the effect of the most
prominent resonant island (corresponding to a width of $\sim 0.05M$ and
average crossing time $\sim 4\times 10^{4}M\sim 3$ days for the parameters
used), other islands surrounding less prominent resonant stable orbits still
exhibit a similar, but less profound, glitch.

_Discussion._ Future space-borne GW detectors, such as LISA and Taiji lisafreq
; taiji , will unlock the detection realm to a wider range of GW sources. The
detection of GWs from EMRIs, which consist of small mass companions orbiting
around supermassive compact objects, in much wider orbits than those already
observed by ground-based detectors, will provide significant information on
the validity of GR in the strong-field region. If the astrophysical
environment around EMRIs ‘modifies’ the underlying theory of gravity in such a
way that the integrability of the equations of motion is broken, then the
orbital phase space of the small companion will contain a series of Birkhoff
islands, by virtue of the chaotic dynamics arnold89 .

By employing the hybrid _kludge_ scheme kludge1 ; glamp07 ; gair06  to evolve
EMRIs in a deformed and a non-Kerr spacetime, introduced in dest20 , and the
Einstein-quadrupole approximation thorne80  to model GW emission, we have
explored the detectability of transient phenomena in the gravitational
waveform, which designate a crossing through a Birkhoff island, and thus probe
spacetime symmetry and constrain a potential departure from a Kerr
description.

Our results indicate that the continuous evolution of GW frequencies is
abruptly, though consistently, broken during a Birkhoff island crossing, which
leads to a glitch, similar to that seen in pulsars and in the LISA pathfinder
data stream, which until now is assumed to be instrumental noise. The glitches
displayed here are present in the majority of the periodograms of the
prominent GW frequencies and have a clear astrophysical origin.

Taking into account that a plethora of initial orbital parameters can
eventually lead an inspiral through an island of strong resonance apo09 ; luk
, and that an object can potentially occupy an island for up to a week Lukes-
Gerakopoulos:2014dpa , we speculate that space-borne detectors should be able
to unveil or constrain the existence of chaotic phenomena in EMRIs, which may
be associated with deviations from GR (though cf. Refs. bara07 ; amaro11 ;
kiu04 ; spin1 ). Furthermore, we argue that such abrupt frequency jumps may
have a true astrophysical origin, thus discarding them _a priori_ from the
data stream as instrumental artifacts may miss potential ‘smoking gun’
physics.

###### Acknowledgements.

_Acknowledgments._ The authors would like to thank Georgios Lukes-Gerakopoulos
for helpful discussions. AGS is supported by the Alexander von Humboldt
Foundation.

## References

  * (1) P. Amaro-Seoane et al., [LISA Collaboration] arXiv:1702.00786.
  * (2) W. R. Hu and Y. L. Wu, Natl. Sci. Rev. 4, 685 (2017).
  * (3) A. Ori and K. S. Thorne, Phys. Rev. D 62.12 (2000): 124022.
  * (4) C. P. L. Berry, et al., Phys. Rev. D 94, 124042 (2016).
  * (5) F. D. Ryan, Phys. Rev. D 52, 5707 (1995).
  * (6) K. Glampedakis and S. Babak, Classical and Quantum Gravity 23.12 (2006): 4167.
  * (7) E. Barausse, et al., Gen. Rel. Grav. 52, 1 (2020).
  * (8) D. C. Robinson, Phys. Rev. Lett. 34.14 (1975): 905.
  * (9) E. Berti, A. Buonanno, and C. M. Will, Classical and Quantum Gravity 22.18 (2005): S943.
  * (10) L. Barack and C. Cutler, Phys. Rev. D 75, 042003 (2007).
  * (11) B. F. Schutz, Classical and Quantum Gravity 26.9 (2009): 094020.
  * (12) T. Apostolatos, G. Lukes-Gerakopoulos, and G. Contopoulos, Phys. Rev. Lett. 103, 111101 (2009).
  * (13) G. Lukes-Gerakopoulos, T. A. Apostolatos, and G. Contopoulos, Phys. Rev. D 81, 124005 (2010).
  * (14) G. Contopoulos, G. Lukes-Gerakopoulos, and T. A. Apostolatos, International Journal of Bifurcation and Chaos 21, 2261 (2011).
  * (15) A. Cárdenas-Avenda$\tilde{\text{n}}$o et al., Class. Quant. Grav. 35, 165010 (2018).
  * (16) K. Glampedakis and G. Pappas, Phys. Rev. D 99, 124041 (2019).
  * (17) K. Destounis, A. G. Suvorov, and K. D. Kokkotas, Phys. Rev. D 102, 064041 (2020).
  * (18) B. Carter, Phys. Rev. 174, 1559 (1968).
  * (19) D. C. Wilkins, Phys. Rev. D 5, 814 (1972).
  * (20) J. Moser, Nachr. Akad. Wiss. Göttingen Math. Phys. Kl. II, 1 (1962).
  * (21) V. I. Arnold, Russian Math. Survey 18, 9 (1963).
  * (22) H. Poincaré, Rend. Circ. Mat. Palermo 33, 375 (1912).
  * (23) G. D. Birkhoff, Trans. Am. Math. Soc. 14, 14 (1913).
  * (24) V. I. Arnold, _Mathematical Methods of Classical Mechanics_ , (Springer, New York, 1989), 2nd ed.
  * (25) J. Brink, M. Geyer and T. Hinderer, Phys. Rev. Lett. 114, no.8, 081102 (2015).
  * (26) J. Brink, M. Geyer, and T. Hinderer, Phys. Rev. D 91, 083001 (2015).
  * (27) G. Contopoulos, _Order and Chaos in Dynamical Astronomy_ (Springer, New York, 2002).
  * (28) M. Armano, et al., Phys. Rev. Lett. 120, 061101 (2018).
  * (29) M. C. Edwards, et al., Phys. Rev. D 102, 084062 (2020).
  * (30) T. Johannsen, Phys. Rev. D 88, 044002 (2013).
  * (31) S. Vigeland, N. Yunes, and L. C. Stein, Phys. Rev. D 83, 104027 (2011).
  * (32) R. Konoplya, L. Rezzolla, and A. Zhidenko, Phys. Rev. D 93, 064015 (2016).
  * (33) D. Psaltis, et al., Phys. Rev. Lett. 100, 091101 (2008).
  * (34) A. G. Suvorov, Gen. Rel. Grav. 53, 1 (2021).
  * (35) Mino, M. Sasaki, and T. Tanaka, Phys. Rev. D 55, 3457 (1997).
  * (36) T. C. Quinn and R. M. Wald, Phys. Rev. D 56, 3381 (1997).
  * (37) E. Poisson, A. Pound, and I. Vega, Living Reviews in Relativity 14, 7 (2011).
  * (38) Y. Mino, Phys. Rev. D 67, 084027 (2003).
  * (39) D. V. Gal’tsov, Journal of Physics A: Mathematical and General 15, 3737 (1982).
  * (40) N. Sago et al., Progress of Theoretical Physics 114, 509 (2005).
  * (41) L. Barack and C. Cutler, Phys. Rev. D 69, 082005 (2004).
  * (42) S. Babak et al., Phys. Rev. D 75, 024005 (2007).
  * (43) J. R. Gair and K. Glampedakis, Phys. Rev. D 73, 064037 (2006).
  * (44) J. R. Gair, C. Li and I. Mandel, Phys. Rev. D 77, 024035 (2008).
  * (45) G. Pappas and T. P. Sotiriou, Phys. Rev. D 91, 044011 (2015).
  * (46) P. Canizares, J. R. Gair, and C. F. Sopuerta, Phys. Rev. D 86, 044010 (2012).
  * (47) A. G. Suvorov, Phys. Rev. D 99, 124026 (2019).
  * (48) A. G. Suvorov, Acta Phys. Pol. B Proc. Suppl. 13, 283 (2020).
  * (49) C. Liu, W. H. Ruan, and Z. K. Guo, Phys. Rev. D 102, 124050 (2020).
  * (50) K. S. Thorne, Astrophys. J. 191, 507 (1974).
  * (51) C. Cutler,.Phys. Rev. D 57, 7089 (1998).
  * (52) K. S. Thorne, Reviews of Modern Physics 52, 299 (1980).
  * (53) G. Lukes-Gerakopoulos, G. Contopoulos and T. A. Apostolatos, Springer Proc. Phys. 157, 129-136 (2014).
  * (54) E. Barausse et al., Phys. Rev. D 75.6 (2007): 064026.
  * (55) P. Amaro-Seoane et al., The Astrophysical Journal Letters 744, L20 (2011).
  * (56) K. Kiuchi and K. Maeda, Phys. Rev. D 70, 064036 (2004).
  * (57) G. Lukes-Gerakopoulos, J. Seyrich, and D. Kunst, Phys. Rev. D 90, 104019 (2014).


# Electronic ground-state hysteresis under magnetic field in GdMn2O5

V. Balédent  Université Paris-Saclay, CNRS, Laboratoire de Physique des
Solides, 91405, Orsay, France<EMAIL_ADDRESS>A.
Vaunat  Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides,
91405, Orsay, France.  Synchrotron SOLEIL, L’Orme des Merisiers, Saint Aubin
BP 48, 91192, Gif-sur-Yvette, France  Laboratoire Léon Brillouin, CEA, CNRS,
Université Paris-Saclay, 91191, Gif sur Yvette, France  S. Petit  Laboratoire
Léon Brillouin, CEA, CNRS, Université Paris-Saclay, 91191, Gif sur Yvette,
France  L. Nataf  Synchrotron SOLEIL, L’Orme des Merisiers, Saint Aubin BP 48,
91192, Gif-sur-Yvette, France  S. Chattopadhyay  Dresden High Magnetic Field
Laboratory (HLD-EMFL), Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden,
Germany  UGC-DAE Consortium for Scientific Research Mumbai Centre, 246-C CFB,
BARC Campus, Mumbai 400085, India  S. Raymond  Université. Grenoble Alpes,
CEA, IRIG, MEM, MDN, 38000 Grenoble, France  P. Foury-Leylekian  Université
Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405, Orsay, France.

###### Abstract

In this paper, we investigate the physical properties of the type II
multiferroic GdMn2O5 material by means of neutrons scattering, electric
polarization and magnetization measurements. A complex (T,H) phase diagram
shows up, with especially a field induced magnetic transition around 12 T at
low temperature. The high field phase is accompanied by an additional electric
polarization along both the $a$ and $b$ directions, as authorized by symmetry,
but never observed experimentally up to now. While the magnetic properties
recover their initial states after driving the field back to zero, the
polarization along $a$ shows a significant increase. This behavior is observed
for all directions of the magnetic field. It constitutes a novel and striking
manifestation of the magneto-electric coupling, resulting in the establishment
of a new ground state at zero magnetic field.

##  I Introduction

The search for new materials with remarkable properties is a major concern for
many condensed matter physicists. In this quest, a simple idea consists in
combining, within the same material, two different properties, even those
which are a priori mutually exclusive. Magneto-electrical multiferroics are
one of the outstanding examples. These materials allow for the manipulation of
the magnetic state using an electric field via the coupling between
ferroelectricity and magnetism, which represents a great potential in the
field of spintronics or information storage. Different routes have been
explored to obtain such properties. Artificial materials, such as
heterostructures Garcia  _et al._ (2015) for instance, consist in alternating
layers with different properties. The coupling is then induced by proximity
effect. Another extensively studied route focuses on materials whose magnetic
and ferroelectric properties coexist in the bulk. Even if quite difficult to
find, their list is getting longer and longer. If ferroelectricity and
magnetism have a distinct microscopic origin, these materials are called type
I multiferroics. One of the main interests of this family is to offer
multiferroic properties at room temperature, as in BiFeO3, but with the
disadvantage of a weak coupling inherent to the distinct origin of the two
orders. In type II multiferroics, ferroelectricity is induced by magnetism,
the two orders are then intrinsically strongly coupled, as in RMn2O5
manganites for instance. This last variety of compounds has received much
attention because of the fundamental problem posed by the origin of this
intrinsic coupling, which results in a complex ground state where lattice,
electronic and magnetic degrees of freedom are entangled.

In this study, we report the possibility of manipulating the electronic ground
state of a multiferroic material in a non-reversible manner using a magnetic
field. We show that a first order magnetic transition occurs at 12 T,
involving both charge and spin degrees of freedom. While the electronic
properties exhibit significant hysteresis, the magnetism exhibits only a weak
1 T wide loop. This opens up a new avenue for creating a new functional ground
state in multiferroic materials.

The type II RMn2O5 multiferroic family attracted attention for multiple
reasons. First, while the average space group is $Pbam$ with a=7.2931 Å, b=
8.5025 Å, c=5.6743 Å, this family crystallizes in a non-centrosymmetric
structure at room temperature Balédent  _et al._ (2015). The actual space
group $Pm$, with the mirror perpendicular to the $c$ axis of $Pbam$
effectively allows for electrical polarization in the $(a,b)$ plane. While
such room temperature polarization has not been confirmed experimentally, it
has been calculated using DFT Dai  _et al._ (2020). This polarization shows
two components, electronic and ionic, whose amplitudes depend on the path
followed in the temperature-electric field phase diagram. Second, one of the
magnetic transitions, below TN= 33K, is accompanied by the rise of an
additional contribution along the $b$ direction Lee  _et al._ (2013), that
adds to the initial polarization. It has been proven experimentally that the
exchange-striction mechanism is at play in this family Yahia  _et al._ (2018).
This additional spin-induced contribution to the polarization is maximized in
GdMn2O5 , where it culminates at 360 nC.cm-2. Third, magnetic field has proven
very efficient in modifying the electronic ground state in the whole RMn2O5
family. Indeed, a modulation of polarization has been observed below 2 T in
GdMn2O5 Lee  _et al._ (2013) and the emergence of a polarization together with
ferromagnetism has been reported in PrMn2O5 above 15 T Chattopadhyay  _et al._
(2020). All these reasons motivate the present high magnetic field study of
GdMn2O5 .

The possibility to modify electronic properties with magnetic field along $a$
has been reported recently Ponet  _et al._ (2022). The authors propose that
the magnetic field can induce a transition between four different topological
magnetoelectric states, formed by the two zero field states and the two high
field ones. The orientation of the magnetic field is also important. For a
field purely along the $a$ direction, the final state is identical to the
initial state, the polarization along $b$ is unchanged and the magnetic state
either. For a field at 10 degrees from the $a$-axis in the $(a,b)$ plane, the
initial and final states however are different, changing the polarization
along the $b$-axis. The magnetic structure is also modified, but is
experimentally indistinguishable from the initial structure: there are two
identical degenerate states respecting the very same symmetries and magnetic
space group. The existence of these two equivalent magnetic structures at zero
field giving rise to different polarizations had already been pointed out by
DFT calculationsDai  _et al._ (2020). The possibility to switch from one to
the other using a magnetic field along the $a$-axis to modify the polarization
along the $b$-axis seems to be linked to the topological properties of these
states. We present here a combination of magnetization, electric polarization
and neutron diffraction to investigate further this mechanism.

##  II Results

Figure 1: Magnetization along the three different crystallographic axes as a
function of magnetic field at T=2 K.

###  II.1 Magnetization

We performed magnetization measurements with the magnetic field applied along
the three crystallographic axes of $GdMn_{2}O_{5}$. These field-sweep
measurements were carried out up to $\sim$57 T at 2 K using pulsed magnetic
fields available at the Dresden High Magnetic Field Laboratory (HLD-EMFL). The
duration of each pulse was $\sim$20 ms. Reproducibility of the data was
ensured by repeating the measurements for each direction. As can be observed
in Fig. 1, the magnetization along the $b$ direction shows a clear hysteresis
loop around 10 T. Such hysteresis is also visible for the field along $a$ and
$c$ around 5 and 11 T respectively. This corroborates the results obtained for
the $a$ direction in ref Lee  _et al._ (2013).

###  II.2 Electric polarization

Figure 2: Variation of the electric polarization along $a$, $b$ and $c$ with
respect to zero field initial state as a function of increasing and decreasing
magnetic field at 2 K. $\Delta P_{a}$ at zero field after field cycling is
indicated for each direction of the magnetic field. Its average among the 3
orientations is around 10 nC.cm-1

Motivated by the strong magneto-electric coupling in this material, we turned
to the evolution of the electric polarization $P$ under magnetic field.
Similar to the magnetization, the measurements were also performed up to
$\sim$57 T using a pyroelectric technique Chattopadhyay  _et al._ (2020) at
the HLD with a pulse duration of $\sim$20 ms. Both sides of the $a$, $b$, and
$c$ thinned samples were covered using silver paste with gold contact wires
attached to each of the sides. For each direction of the magnetic field we
measured the components of the electric polarization along the 3
crystallographic directions : in total we performed 9 configurations at 2 K.
The field variations of the pyroelectric current $I$ were recorded by
measuring the voltage variation across a 1 M$\Omega$ shunt resistor connected
in series with the measurement circuit by a digital oscilloscope (Yokogawa
DL750). To calculate the field dependent electric polarization for each of the
nine configurations, we integrated the $I$($H$) data. All the measurements
were repeated to verify the reproducibility of the data.

Results are shown in Fig. 2. First, regardless of the direction and intensity
of the magnetic field, no polarization has been measured along the $c$
direction ($P_{c}\equiv 0$, blue curves). The observed small variations are
most likely due to imperfect alignment of the crystal and the electrical
contacts. As expected owing to the $Pm$ space group with the mirror
perpendicular to $c$ axis, the electric polarization lies within the $(a,b)$
plane ($P_{a}$ and $P_{b}\neq 0$). $P_{b}$ (green curve), however shows a
step-like increase, above 12, 10 and 15 T for increasing field along $a$, $b$
and $c$ respectively. Decreasing the field from 57 T, the polarization goes
back to its initial value with a hysteresis loop.

For H along $b$, this behavior is similar to the one observed in magnetization
at the same critical field around 10 T. Hence, the high field magnetic phase
induces an additional ferroelectric component $P_{b}=$4 nC.cm-2, on top of the
one initially present at zero field (360 nC.cm-2). Interestingly, such
increase in $P_{b}$ is also visible for magnetic fields along $a$ and $c$.
It’s worth noting that only $P_{b}$ polarization has been experimentally
measured at low temperature in the $CM1$ phase. Although the $Pm$ space group
also allows non-zero $P_{a}$ polarization, this has never been demonstrated
experimentally. Applying a magnetic field changes the game : as can be seen in
Fig. 2, $P_{a}$ increases with increasing the field, similarly to $P_{b}$.
Even more astonishing, with decreasing the field back to zero, $P_{a}$ does
not go back to zero, but shows a residual polarization $\Delta P_{a}$ between
6 and 15 nC.cm-2 at 0 T depending on the field orientation and thus around 10
nC.cm-1 in average.

###  II.3 Neutron scattering

In order to investigate further the magnetic properties under magnetic field
of GdMn2O5, the temperature-magnetic field phase diagram was investigated by
means of single crystal neutron diffraction. The experiment was carried out on
the IN12 cold triple axis spectrometer, a CEA-Juelich CRG installed at ILL
(Grenoble, France). Since the magnetic propagation vector is $(0.5,0,0)$ or
$(0.5,0,\delta)$ depending on the temperature in all members of the RMn2O5
family, the natural scattering plane is the $(a,c)$ plane, allowing access to
magnetic peaks of the form $(H,0,L)$. The magnetic field was applied along the
vertical direction, hence the $b$ direction. The evolution of selected
reciprocal space regions was then recorded as a function of magnetic field and
temperature, leading to the schematic phase diagram shown in Fig. 3. Several
structures appear, characterized by different Bragg reflections and that we
shall put in different categories:

  * •

$q_{ICM1}=(3/2,0,\epsilon\approx 0.19)$ (gray), $q_{CM1}=(3/2,0,0)$ (red) and
$q_{ICM2}=(3/2,0,\delta\approx 0.43)$ (blue) are characteristic of the zero
field, high, intermediate and low temperature structures respectively, already
reported in literature Yahia  _et al._ (2018). The incommensurability along
the $c$ axis indicates that the magnetic moments wrap around the $c$ axis,
forming a helicoidal type structure, and reflecting an exchange frustration
along this direction. It is worth mentioning that the $q_{CM1}$ phase is
fragile: strikingly, the spin wave dispersion, measured on the very same
crystal in the 5 to 30 K temperature range Vaunat  _et al._ (2021), does not
go soft at $q_{CM1}$ but at $q_{ICM2}$, indicating that the incommensurate
phase is the incipient ground state.

  * •

$q_{FM}=(1,0,0)$ (green) reflects a ferromagnetic component, imposed by the
field, with magnetic moments also along the applied magnetic field.

  * •

$q_{CM2}=(3/2,0,1/4)$ (magenta) and $q_{ICM3}=(3/2,0,0.2-0.3)$ (yellow) are
characteristic of field induced high and low temperature structures
respectively. The magnetic structure inferred from $q_{CM2}$ is a quadrupling
of the $q_{CM1}$ magnetic unit cell along the $c$-axis, and can also be seen
as a commensurate lock-in of the periodicity along the same direction with
respect to the $q_{ICM2}$ phase. Interestingly, this $q_{CM2}$ magnetic
propagation vector corresponds to the ground state of the low temperature
magnetic phase observed in numerous other members of the family (R=Tb, Dy, Ho,
Er, Tm and Yb).

Overall, and apart from the appearance of the ferromagnetic component, the
effect of the field is to restore the $q_{CM2}=(3/2,0,1/4)$ magnetic phase,
common to the other members of the RMn2O5 family. Here, such cycloidal
structure wraps around the $c$ axis, allowing to align 2 spins out of 4 in a
direction perpendicular to the field, which limits the Zeeman energy and
satisfies the antiferromagnetic exchanges. The obtained phase diagram is also
similar to the one derived from magnetization and dielectric constant
measurements on this compound Bukhari  _et al._ (2016), with comparable number
of phases and consistent phase boundaries in the field-temperature diagram.
Details of the measured temperature and field evolution of those different
Bragg reflections are gathered in Fig. 4 and 6. Fig. 4 especially shows
Q-scans concatenated to produce maps as a function of wavevector and either
temperature or field. Fig 6 shows the magnetic field dependence of the
intensity and/or q-position determined from a fit of those scans.

Figure 3: Schematic temperature-magnetic field phase diagram deduced from the
neutron diffraction measurements. Points represents the region in the phase
diagram where a magnetic signal has been observed at the different wavevectors
indicated in the labels. The presence of multiple colored stripes denotes the
coexistence of several magnetic peaks with different propagation vectors.
Figure 4: Color map of the magnetic intensity on the (1,0,L) and (1.5,0,L)
Bragg reflections as a function of temperature and magnetic field. Color
scales are not comparable to each other.

Evolution at 2K under magnetic field.  At 2 K, a contribution appears at
$q_{FM}$=(1,0,0), where nuclear contribution is forbidden by the average
structure $Pbam$ (see Fig 4a). As represented in Fig. 6c, the intensity of
this peak increases and reaches a maximum for H=4 T and decreases
progressively until 12 T, where it stabilizes at a minimum value up to 15 T.
Fig. 4e-h shows that the transition temperature increases from 15 K at 1 T to
more than 45 K at 15 T. As it preserves the lattice translation symmetry, it
can be interpreted as a ferromagnetic contribution emerging progressively
under magnetic field. Owing to the absence of magnetic anisotropy at the Gd
site, contrary to the Mn site, and since the coupling between the Gd and Mn
moments is the weakest in this compound Vaunat  _et al._ (2021), it is very
likely that Gd is the main contributor to this signal. To confirm this
interpretation, we performed X-ray Magnetic Circular Dichroism on a GdMn2O5
powder sample at the ODE beamline, synchrotron SOLEIL, at both Mn K edge and
Gd L3 edge Baudelet  _et al._ (2011, 2016). As reported in Fig. 5, the
measurement at 1.3 T and 6 K exhibits a clear ferromagnetic contribution from
the Gd but no sizeable contribution from the Mn moments with experimental
sensitivity. By comparing with the literature Mito  _et al._ (2009); Sorg  _et
al._ (2007); Ramos  _et al._ (2013), we were able to estimate the magnitude of
the moment to be about 2.5 $\mu_{B}$ for Gd and below 0.05 $\mu_{B}$ for Mn.
These results confirm the rise of a ferromagnetic component under magnetic
field and further attest that this contribution is due to the Gd moments. This
can easily be explained with the same mechanism reported for PrMn2O5
Chattopadhyay  _et al._ (2020) replacing the exchange interaction values by
the one of GdMn2O5 Vaunat  _et al._ (2021).

Figure 5: X-ray Absorption Spectra and associated dichroic signal at Gd (left)
and Mn (right) edge, performed at 1.3 T and 6 K. Dichroic signals are
multiplied by a factor 100 and 1000 respectively for a sake of visibility.

Coming back to neutron scattering results at 2 K, Fig. 4b and Fig. 6e show
further that the (commensurate) magnetic signal at $q_{CM1}$=(3/2,0,0)
decreases rapidly up to 5 T while a very small intensity remains up to HC=11
T. This $q_{CM1}$ intensity coexists at low field with the magnetic phase
$q_{ICM2}$=($\frac{3}{2}$,0,$\delta$=0.425) below 4 K, as can be seen on the
maps Fig. 4b and 4i. The evolution of the intensity and position of $q_{ICM2}$
is reported in Fig. 6b and Fig. 6a respectively. From 0 to 10 T, the intensity
grows and the position shifts from $\delta$=0.425 to $\delta$=0.37. Between 0
and 1 T the position remains the same, suggesting a possibility that the 0.425
position may correspond to a commensurate order with q=3/7$\approx 0.428$. At
HC=11 T a magnetic transition occurs, the $q_{ICM2}$ signal disappears, and
another magnetic signal appears at $q_{CM2}$=($\frac{3}{4}$,0,$\frac{1}{4}$).
This critical value of magnetic field is similar to the observed transition in
both magnetization measurements and electric polarization for field along $b$
direction.

Evolution above 33 K under magnetic field.  The evolution of the high
temperature magnetic order as function of the field was also investigated. At
35 K, Fig. 4c shows that the magnetic peak $q_{ICM1}$=($\frac{3}{2}$,0,
$\epsilon$) shifts from $\epsilon$=0.17 to $\epsilon$=0.22 with a decreasing
intensity from 0 to 12 T. Interestingly, a commensurate peak appears above H=5
T at $q_{CM2}$=($\frac{3}{2}$,0,$\frac{1}{4}$), the very same position as the
high field low temperature one, showing that this new magnetic order
stabilizes at lower field at high temperature. Performing the same study at 40
K shows no sign of the commensurate phase at $q_{CM2}$, while $\epsilon$ moves
from 0.19 to 0.22 between 0 and 15 T (see Fig. 4d).

Temperature evolution of field induced magnetic phases. Despite a strong
reduction of the intensity of the ferromagnetic contribution at $q_{FM}$ above
4 T, the temperature range of existence of such contribution increases from
below 10 K at 1 T up to 40 K at 15 T (Fig. 4e-h). However, the transition
temperature of $q_{CM1}$ slowly decreases from 33 K at zero field down to 25 K
at 10 T (Fig. 4i-j). On the contrary, $q_{ICM2}$ shows an increase of its
transition temperature from 3.75 to 40 K between 1 and 15 T (Fig. 4j-m). Both
disappear above the transition field HC=11 T. At 10 T however, $q_{CM2}$
appears in the narrow temperature range between 25 and 40 K. Above 12 T, this
phase is present below 35 K down to 2 K (Fig. 4k). An additional subtlety
appears at 12 T around this $q_{CM2}$ magnetic peak. Indeed, two satellite
peaks appear at $q_{ICM3}$=(1.5,0,0.2) and (1.5,0,0.3) from 2 to 15 K with a
maximum at 10 K (Fig. 4l). These satellites are still present at 15 T but
disappear above 10 K. Further investigations are necessary to better identify
the corresponding magnetic structures but we anticipate that these satellites
may be the signatures of some kind of roughening of the cycloid around the
$c$-axis or discommensurations.

Hysteretic behaviour. In order to investigate the nature of the transition at
HC=11 T, we performed the same measurement while decreasing the field. Fig. 6
shows the field evolution of these peaks. A clear hysteresis is visible on the
intensity of $q_{ICM2}$ (6b) and $q_{CM2}$ (6d) with a width around 2 T. This
behavior, together with the coexistence of both orders between 10 and 12 T,
strongly suggests a first order transition.

Figure 6: Evolution as a function of the magnetic field along the $b$ axis and
T=2 K of the incommensurate position $\delta$ of the
$q_{ICM2}$=$(\frac{3}{2},0,\delta)$ phase (a), its associated intensity
$I_{ICM2}$ (b), the intensity of the ferromagnetic contribution $I_{FM}$ (c),
the intensity $I_{CM2}$ at $q_{CM2}$=$(\frac{3}{2},0,\frac{1}{4})$ (d) and the
$q_{CM1}$=$(\frac{3}{2},0,0)$ intensity $I_{CM1}$.

##  III Discussion

A similar hysteresis in electric polarization as a function of magnetic field
was reported on the same GdMn2O5 compounds Ponet  _et al._ (2022), with the
difference that the magnetic field was oriented 10 degrees off the $a$
direction in the $(a,b)$ plane, with a hysteresis effect on the polarization
along $b$. The authors interpreted this result arguing the presence of 4
topological states, that can be manipulated by the orientation of the applied
magnetic field. The authors further relate these 4 states to 4 different spin
configurations within the same zero-field unit cell, doubled along the $a$
direction, corresponding to the $q_{CM1}$ propagation wave vector. The
orientation of the magnetic field in the $(a,b)$ plane seems crucial and the
angle between the field and $a$-axis very specific to obtain such topological
transition. The proposed analysis is based on a simple model consisting in two
spin chains per unit cell, taking into account the anisotropy and only two
exchange interactions (intra-chain along $a$ and inter-chains along $b$). This
simple model is also based on the assumption that the intra-chain exchange
interaction is dominant in as observed in most RMn2O5.

Although this simplified model allows to introduce the interesting topological
aspects of these different ground states, the present results show that the
model lacks some ingredients to capture the physics of the high magnetic field
phase. First we show here that the topological transition can also be observed
for a totally different field direction (here along $b$), with an effect on
the polarization along $a$. Secondly, inelastic neutron scattering experiments
single out GdMn2O5 as an exception among the RMn2O5 family, with a weak intra-
chain coupling Vaunat  _et al._ (2021) thus not dominant. Finally, the present
results show the importance of the coupling along the $c$ axis and call for a
3D model. Indeed, i) the magnetic transition is accompanied by a quadrupling
of the unit cell in the $c$ direction (from $q_{CM1}$ to $q_{CM2}$) and ii)
the anomaly at Hc is systematically observed in susceptibility and
polarization measurements, whatever the direction of the magnetic field, which
indicates a change of magnetic order whatever the direction of the magnetic
field.

Figure 7: a) Schematic representation of the electric polarization
contribution from all degenerate states at room temperature (red arrows) taken
from ref Dai  _et al._ (2020). Blue arrow correspond to the additional spin
induced contribution at low temperature. Summing all contribution lead to a
b-component only polarization. b) Proposed remaining states after an
hysteretic cycle with magnetic field along b above Hc=11 T. While the $b$
component along $b$ remain unchanged, a finite component appears along $a$.

We propose a simpler reasoning based on DFT calculations of the electronic
ground state Dai  _et al._ (2020) at zero field. According to the symmetry of
the room temperature $Pm$ space group, the two components of polarization
along $a$ ($P_{a}$) and $b$ ($P_{b}$) are allowed. This results in 4
degenerate “high temperature” configurations depending on the relative sign of
these components : $(++),(+-),(-+)$ and $(--)$ as represented by the red
arrows in Fig. 7a. These four configurations are associated with two different
spin structures that are equivalent from a symmetry point of view and cannot
be discriminated. Interestingly, these two spin configurations correspond to
the two topological phases proposed at low field in Ref Ponet  _et al._
(2022). According to DFT calculations Dai  _et al._ (2020), a spin-induced
contribution adds up below the magnetic transition temperature (TN=33 K) :
$-\Delta P_{a}$ and $+\Delta P_{b}$ for the $(++)$ and $(+-)$ states, and
$+\Delta P_{a}$ and $+\Delta P_{b}$ for the $(-+)$ and $(--)$ states (see blue
arrows in Fig. 7a). The net resulting polarization below TN, summing over the
4 states, is thus along $b$ only, as measured experimentally. Above 12 T, the
$q_{CM2}$ magnetic state induces a new electronic ground state as revealed
experimentally by the changes in the $P_{a}$ and $P_{b}$ component of the
polarization. In an analogous way to the proposed switching model between
topological statesPonet  _et al._ (2022), one can expect the magnetic field to
select only certain states once it has returned to zero. In order to reproduce
the experimental observations, only two possible choices for the final states
seem possible : $(++)+(+-)$ or $(--)+(-+)$ (see Fig. 7b). From a magnetic
point of view, the initial state and the final state after field-sweep, are
equivalent. This is consistent with the present magnetization and neutron
scattering results. From an electronic point of view, however, the final state
can be discriminated from the initial state. Indeed, while $P_{b}$ remains
unchanged (see the $b$ component in Fig. 7.b), $P_{a}$ is no longer
compensated by the two other contributions (coming from $(++)$ and $(+-)$).
DFT calculations predicted a contribution of $\Delta P_{a}\approx$ 10 nC.cm-1
Dai  _et al._ (2020) which is in perfect agreement with our results in Fig. 2.
The detailed mechanism to describe how the field along $b$ selects only
certain electronic states remains to be unveiled.

##  IV Conclusion

In conclusion, the present study reveals a first order magnetic transition
under magnetic field in GdMn2O5 . This transition recalls metamagnetic
transitions characteristic of heavy fermions Stryjewski and Giordano (1977);
Aoki  _et al._ (2013). For H along $b$, the transition is characterized by a
new propagation vector $q_{CM2}$=$(\frac{1}{2},0,\frac{1}{4})$, which is
common among the RMn2O5 family. The new high field phase presents an
additional contribution to the electric polarization along both $a$ and $b$
directions. Surprisingly, the polarization along $a$ does not go back to zero
while decreasing the magnetic field, but remains finite. This indicates a
change of the electronic ground state after magnetic field cycling, while no
indication for such a change of magnetic ground state could be evidenced
neither by magnetic susceptibility nor neutron diffraction. The establishment
of a new electronic ground state after H cycling is further observed for all
directions of the magnetic field. These results challenge previous model of
magneto-electric switching between different topologically protected states.
They should motivate further experimental and theoretical work to investigate
the stability of these new states for other directions of the magnetic field,
and also study the consequence on the dynamical properties such as the
electromagnon. Our work confirms the importance of the path followed in the
three dimensions (T,E,H) phase diagram in the establishment of the ground
state in multiferroic materials as suggested in recent workDai  _et al._
(2020), and should be extended to other external parameters such as pressure.

##  V Acknowledgments

We thank W. Knafo for fruitfull discussions. This study was supported by
grants from LLB and SOLEIL synchrotron, and LabEx PALM through Contract No.
ANR-10-LABX-0039-PALM. Experiments at ILL were sponsored by the French Neutron
Federation (2FDN), with Data References 10.5291/ILL-DATA.CRG-2839 and
10.5291/ILL-DATA.4-01-1700. We acknowledge SOLEIL for provision of synchrotron
beamtime (proposal number 20220714) on ODE beamlines and the support of the
HLD at HZDR, member of the European Magnetic Field Laboratory (EMFL proposals
DMA14-219 and DMA15-219) for magnetization and polarization measurements, and
the MORPHEUS platform at the Laboratoire de Physique des Solides for sample
alignment.

## References

  * Garcia  _et al._ (2015) V. Garcia, M. Bibes, and A. Barthélémy, Comptes Rendus Physique 16, 168 (2015).
  * Balédent  _et al._ (2015) V. Balédent, S. Chattopadhyay, P. Fertey, M. B. Lepetit, M. Greenblatt, B. Wanklyn, F. O. Saouma, J. I. Jang, and P. Foury-Leylekian, Phys. Rev. Lett. 114, 117601 (2015).
  * Dai  _et al._ (2020) J.-Q. Dai, J. Yuan, and C. Ke, Journal of Magnetism and Magnetic Materials 516, 167373 (2020).
  * Lee  _et al._ (2013) N. Lee, C. Vecchini, Y. J. Choi, L. C. Chapon, A. Bombardi, P. G. Radaelli, and S.-W. Cheong, Phys. Rev. Lett. 110, 137203 (2013).
  * Yahia  _et al._ (2018) G. Yahia, F. Damay, S. Chattopadhyay, V. Balédent, W. Peng, S. W. Kim, M. Greenblatt, M. B. Lepetit, and P. Foury-Leylekian, Physical Review B 97, 5 (2018).
  * Chattopadhyay  _et al._ (2020) S. Chattopadhyay, V. Balédent, S. K. Panda, S. Yamamoto, F. Duc, T. Herrmannsdörfer, M. Uhlarz, T. Gottschall, O. Mathon, Z. Wang, C. Strohm, M. Greenblatt, P. Foury-Leylekian, and J. Wosnitza, Phys. Rev. B 102, 094408 (2020).
  * Ponet  _et al._ (2022) L. Ponet, S. Artyukhin, T. Kain, J. Wettstein, A. Pimenov, A. Shuvaev, X. Wang, S. W. Cheong, M. Mostovoy, and A. Pimenov, Nature 607, 81 (2022).
  * Vaunat  _et al._ (2021) A. Vaunat, V. Balédent, S. Petit, P. Roy, J. B. Brubach, G. Giri, E. Rebolini, P. Steffens, S. Raymond, Q. Berrod, M. B. Lepetit, and P. Foury-Leylekian, Phys. Rev. B 103, 174434 (2021).
  * Bukhari  _et al._ (2016) S. H. Bukhari, T. Kain, M. Schiebl, A. Shuvaev, A. Pimenov, A. M. Kuzmenko, X. Wang, S.-W. Cheong, J. Ahmad, and A. Pimenov, Phys. Rev. B 94, 174446 (2016).
  * Baudelet  _et al._ (2011) F. Baudelet, Q. Kong, L. Nataf, J. D. Cafun, A. Congeduti, A. Monza, S. Chagnot, and J. P. Itié, High Pressure Research 31, 136 (2011).
  * Baudelet  _et al._ (2016) F. Baudelet, L. Nataf, and R. Torchio, High Pressure Research 36, 429 (2016).
  * Mito  _et al._ (2009) M. Mito, K. Matsumoto, Y. Komorida, H. Deguchi, S. Takagi, T. Tajiri, T. Iwamoto, T. Kawae, M. Tokita, and K. Takeda, Journal of Physics and Chemistry of Solids 70, 1290 (2009).
  * Sorg  _et al._ (2007) C. Sorg, A. Scherz, K. Baberschke, H. Wende, F. Wilhelm, A. Rogalev, S. Chadov, J. Minár, and H. Ebert, Phys. Rev. B 75, 064428 (2007).
  * Ramos  _et al._ (2013) A. Y. Ramos, H. C. N. Tolentino, M. M. Soares, S. Grenier, O. Bunău, Y. Joly, F. m. c. Baudelet, F. Wilhelm, A. Rogalev, R. A. Souza, N. M. Souza-Neto, O. Proux, D. Testemale, and A. Caneiro, Phys. Rev. B 87, 220404 (2013).
  * Stryjewski and Giordano (1977) E. Stryjewski and N. Giordano, Advances in Physics 26, 487 (1977).
  * Aoki  _et al._ (2013) D. Aoki, W. Knafo, and I. Sheikin, Comptes Rendus Physique 14, 53 (2013), physics in High Magnetic Fields / Physique en champ magnétique intense.