# Purely antiferromagnetic frustrated Heisenberg model in spin ladder compound
BaFe2Se3

A. Roll  Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides,
91405, Orsay, 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
A. Forget  Université Paris-Saclay, CEA, CNRS, SPEC, 91191, Gif-sur-Yvette,
France.  D. Colson  Université Paris-Saclay, CEA, CNRS, SPEC, 91191, Gif-sur-
Yvette, France.  A. Banerjee  Neutron Scattering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA  Department of Physics and
Astronomy, Purdue University, West Lafayette, IN 47906, USA  P. Foury-
Leleykian  Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides,
91405, Orsay, France.  V. Balédent  Université Paris-Saclay, CNRS, Laboratoire
de Physique des Solides, 91405, Orsay, France.  victor.baledent@universite-
paris-saclay.fr

###### Abstract

The spin dynamics in the block magnetic phase of the iron-based ladder
compound BaFe2Se3 has been studied by means of single crystal inelastic
neutron scattering. Using linear spin wave theory and Monte-Carlo simulations,
our analysis points to a magnetic Heisenberg model with effective frustrated
antiferromagnetic couplings only, able to describe both the exotic block order
and its dynamics. This new and purely antiferromagnetic picture offers a
fruitful perspective to describe multiferroic properties but also understand
the origin of the stripe-like magnetic instability observed under pressure as
well as in other parent compounds with similar crystalline structure.

Magnetism plays a prominent role in many exotic properties of condensed
matter, from unconventional superconductivity to multiferroicity and spin
liquid phases. Copper has long been the focus of attention for studying these
exotic phases, thanks to its half-integer spin, which promotes the emergence
of these quantum states. On the other hand, research on multiferroicity has
focused on manganese-based compounds, whose large magnetic moment has made it
possible to obtain sufficiently strong magnetoelectric couplings to allow
applications. More recently, iron-based compounds have found a place at the
interface between these two paradigms. Indeed, unconventional
superconductivity has been discovered in pnictides, and the large iron
magnetic moment combined with magnetic interactions allows to obtain
multiferroic materials at quite high temperatures. BaFe2Se3 embodies this new
path, by showing multiferroic properties below 250K, and superconductivity
above 10 GPa (see phase diagram Fig. 1).

At ambient pressure, BaFe2Se3 crystallizes in the polar Pm space group Zheng
_et al._ (2020); Weseloh  _et al._ (2022), with a weak distorsion from the
average Pnma space goup. For sake of clarity, we shall use the orthorhombic
average space group, in which there are two ladders per unit cell, each formed
by two adjacent legs running along the $b$ axis. Below $T_{N}=$ 250 K, an
exotic magnetic order develops, called ”block order”, consisting in an
antiferromagnetic (AFM) arrangement of square blocks of four parallel iron
spins, pointing mainly along the $a$ direction, i.e. perpendicular to the
ladders. Along the legs, the structure thus shows an up-up-down-down pattern
Caron  _et al._ (2012). This model has been refined using single crystal
neutron diffraction, revealing a tilt of the moments with respect to the $a$
axis, and resulting in an umbrella-like magnetic order Zheng  _et al._
(2022a). Under pressure, this exotic magnetic state gives way to a more
conventional order consisting in stripes of an up-down-up-down pattern along
the legs Zheng  _et al._ (2022b), and analogous to the one observed in parent
compounds with similar crystalline structure. It is worth noting that this
classical up-down-up-down pattern is also stable close to the superconducting
dome (see dark-green region in Fig. 1).

Figure 1: (Color online) Schematic pressure temperature phase diagram of
BaFe2Se3. The block order (up-up-down-down) occupies the yellow region (black
squares), while the stripe order (up-down-up-down) is observed in the green
region (black triangles). Superconductivity appears above 10 GPa (black
cricles from refs. Ying  _et al._ (2017); Zheng  _et al._ (2022b)). The
structural transition from Pm to Cmcm space group is observed at 4 GPa, and
progressive metallization occurs between 4 and 10 GPa.

Understanding the magnetic couplings, at the heart of the different properties
of this compound, requires the study of its excitations. On the basis of a
first inelastic neutron scattering measurement performed on powder samples,
the authors of Ref.Mourigal  _et al._ (2015) have suggested an alternation of
ferromagnetic and antiferromagnetic couplings along the legs. Although this
model clearly stabilizes the correct magnetic structure due to the absence of
magnetic frustration, it raises several questions. While the exchange paths
are admittedly different between up-up and up-down iron bonds along a given
leg Weseloh  _et al._ (2022), they hardly justify the change in both sign and
amplitude, hence the nature of the exchange interactions, as also underlined
in theoretical literature on spin dynamics Herbrych  _et al._ (2018).
Furthermore, the ferromagnetic nature of some interactions seems incompatible
with the Goodenough-Kanamori orbital rules : as the Fe-Se-Fe angle (denoted by
$\psi$) strongly departs from 90 degrees (see Tab. 1), the antiferromagnetic
super-exchange is expected to prevail over the ferromagnetic direct coupling.
Another physical argument in favour of purely antiferromagnetic bonds along
the ladder is the vicinity of the classical antiferromagnetic up-down-up-down
order observed under moderate pressure (4 GPa) and in all other members of the
family presenting a similar crystalline structure, such as KFe2Se3, CsFe2Se3
and BaFe2S3 Caron  _et al._ (2012); Hawai  _et al._ (2015); Takahashi  _et
al._ (2015).

$\psi_{Fe-Se-Fe}$ | 66(7) | 71(6) | 65.8(19) | 101(3) | 105(3)
---|---|---|---|---|---
Involved coupling | J1 | J${}^{{}^{\prime}}_{1}$ | JR | J2 | J${}^{{}^{\prime}}_{2}$
Table 1: Angles (in degrees) between the Fe-Se-Fe bonds for the different
possible paths within the ladder, according to the structure published in ref.
Zheng  _et al._ (2020). The couplings correspond to the exchange interactions
affected by these angles, as represented in Fig. 2g.

In this letter we present new inelastic neutron scattering (INS) measurements
performed on a single crystal of BaFe2Se3 providing unique details on the spin
wave dispersion. The modeling of the data leads to a purely antiferromagnetic
Heisenberg Hamiltonian, allowing one to describe both the magnetic order and
the dynamics. This picture also provides a natural understanding of the
evolution of the magnetism under pressure as well as the magnetic ground state
of parent compounds.

Figure 2:  a) Color map showing the dispersion along $(1/2,K,1/2)$ at 20 K
($E_{i}$ = 125 meV). Neutron intensity is in logscale. b) Experimental energy
cut at $Q=(1/2,1/2,1/2)$ taken at 20 K ($E_{i}$ = 25 meV). The elastic line
was removed (see . ) c) Experimental energy cut at $Q=(1/2,7/2,1/2)$ taken at
20 K ($E_{i}$ = 125 meV). The integration range in Q-space was adapted to
improve visibility ($\Delta$h=0.5, $\Delta$k=0.5, $\Delta$l=10). d) Simulated
color map of $S_{SWT}(Q,\omega)$ along $(1/2,K,1/2)$. Intensity is in
logscale. e) and f) Simulation of the energy cut at $Q=(1/2,1/2,1/2)$. g)
Schematic sketch of the magnetic structure and the different exchange
interactions used in the simulation. Left : projection in the (b,c) plane.
Right : projection in the (a,c) plane. Crystallographic axis correspond to the
$Pnma$ space group setting.

Single crystals were grown following the method described in Ref. Zheng  _et
al._ (2020), and fully characterized . . The INS experiments were carried out
on the wide angular-range time-of-flight chopper spectrometer (ARCS) at the
Spallation Neutron Source (SNS), ORNL, US. The measurements were conducted at
20 K with incident beam energies $E_{i}$ = 25 and 125 meV yielding an energy
resolution of about 1 and 6 meV respectively. Slices were taken from the full
$(Q,\omega)$ dataset to produce maps as a function of energy transfer $\omega$
and wavevector $Q$ along the 3 main crystallographic directions $a,b$ and $c$.
Only dispersion along the $b$ direction, corresponding to the ladder, was
visible, as expected for a quasi uni-dimensional system . . This allowed to
integrate along the other directions to improve statistics. The measured
dispersion obtained with $E_{i}$=125 meV at 20 K along $(1/2,K,1/2)$ is shown
in Fig. 2a. An acoustic mode is clearly visible, dispersing up to 60 meV. The
data obtained with $E_{i}$=25 meV allows a closer look at the low energy part.
By subtracting the elastic part . , a gap is evidenced in the energy cut at
$Q=(1/2,1/2,1/2)$ as shown in Fig. 2b. The gap value around 6 meV is in good
agreement with the one inferred from powder sample data Mourigal  _et al._
(2015). At higher energy, three dispersion-less optical modes are observed. In
order to make them more visible, an integration along H ($\Delta$H=0.5), K
($\Delta$K=0.5) and L ($\Delta$L=10) was performed, yielding the energy
dependence displayed in Fig. 2c. They appear at 92(2), 101(2) and 110(2) meV.
Only two of these modes were resolved in the powder experiment Mourigal  _et
al._ (2015), providing a strong additional constraint on the subsequent
modeling.

To model the spin dynamics, an exhaustive survey of the structures stabilized
by different sets of couplings is necessary. To this end, classical Monte-
Carlo was used to determine the phase diagram of a single ladder (see details
in . ). To describe the system, we used the following Heisenberg model:

$\centering
H_{mag}=\sum_{i<j}J_{ij}\textbf{S}_{i}.\textbf{S}_{j}+\frac{3}{2}D\sum_{i}(\textbf{S}_{i}.\textbf{n}_{i})^{2}\@add@centering$
(1)

the first sum runs over iron magnetic ions, $\textbf{S}_{i}$ denote the spin
at site $i$ and $J_{ij}$ are exchange integrals between spins at different
sites. As sketched in Fig. 2g, $J_{1}$ and $J^{{}^{\prime}}_{1}$ correspond to
nearest neighbor interactions, and $J_{3}$ to next nearest neighbor
interactions within a given leg of the ladder. These couplings may lead to
frustration in a purely antiferromagnetic model. $J_{R}$, $J_{2}$ and
$J^{{}^{\prime}}_{2}$ couple nearest and next nearest neighbor spins
respectively in adjacent legs. The second sum describes the magnetic
anisotropy : $3D/2$ is the single ion anisotropy and $\textbf{n}_{i}$ defines
the local anisotropy axis. The results of our calculations are presented in
Fig. 3, with maps of the propagation vector value as a function of the various
couplings in units of $J_{3}$. The unit of the propagation vector is defined
by a lattice consisting in a single magnetic site per unit cell. This means
that the purely antiferromagnetic order corresponds to $k_{y}$=0.5 (doubling
of the unit cell), and the block order, with 4 spins per unit cell,
corresponds to $k_{y}$=0.25. As can be intuitively expected, provided
$J_{2}=J^{\prime}_{2}=J_{R}$=0, a fully ferromagnetic order with $k_{y}$=0 is
stabilized if $J_{1}$ and $J^{\prime}_{1}$ are ferromagnetic and tend to
infinity ($J_{1},J^{\prime}_{1}\rightarrow-\infty$). In contrast, a classical
antiferromagnetic order with $k_{y}$=0.5 is stabilized if these parameters are
antiferromagnetic ($J_{1},J^{\prime}_{1}\rightarrow+\infty$), as presented in
Fig. 3a. The large green surfaces correspond to the $k_{y}$=0.25 propagation
vector, hence to the up-up-down-down pattern typical of the block order. As
expected, this peculiar ordering becomes the ground state as soon as $J_{1}$
and $J^{{}^{\prime}}_{1}$ depart from each other. The symmetry with respect to
the diagonal $J_{1}=J^{\prime}_{1}$ is due to the symmetric role of $J_{1}$
and $J^{\prime}_{1}$. The block order is eventually obtained taking into
account $J^{{}^{\prime}}_{2}$ and $J_{2}$, which couple the two legs of the
ladder. The map presented in Fig. 3b shows the propagation vector along the
ladder for a set of $J^{{}^{\prime}}_{1}$ and $J_{1}$ which stabilizes the up-
up-down-down sequence. The values used for $J^{{}^{\prime}}_{1}$ $J_{1}$ and
$J_{R}$ are given in Tab. 2. Interestingly, the block order remains stable
provided that $J^{{}^{\prime}}_{2}$ or $J_{2}$ is antiferromagnetic. This
indeed contributes to the stabilization of identical up-up-down-down patterns
in adjacent legs.

Figure 3: Phase diagram obtained based on the propagation vector determined by
Metropolis-Hasting simulations of a single ladder as a function of different
couplings. The colorscale indicates the value of $k_{y}$ along the ladder. In
a, the calculations are performed for $J_{2}$=$J^{\prime}_{2}$=$J_{R}$=0,
while in b, $J_{1}$, $J^{\prime}_{1}$ and $J_{R}$ values are given in Tab. 2.

Having reduced the possible range for relative exchange interaction strengths
stabilizing the ground state, the spin dynamics was simulated using linear
spin wave theory Petit (2011). Each coupling constant $J_{i}$ was varied
independently from the others within the appropriate region of the Monte-Carlo
phase diagram to disentangle their role in shaping the dispersion. First,
$J_{1}$, $J_{R}$ and $J_{2}$ mainly affect the optical modes bandwidth and
energies. $J_{3}$, $J^{{}^{\prime}}_{1}$ and $J^{{}^{\prime}}_{2}$ affect both
acoustic and optical modes. Finally, $J_{4}$ , $J_{5}$ and the local
anisotropy have both an effect on the low energy gap of the acoustic mode. The
local anisotropy axis orientation $\textbf{n}_{i}$ has been chosen to
reproduce the umbrella-like magnetic ground state Zheng  _et al._ (2022a),
with especially a non-zero component along both the $b$ and $c$ directions.
The exchange couplings that best fit the experiment are listed in Table 2. The
simulated dispersion is represented Fig. 2d. As can be seen, the acoustic
branches are well reproduced and the 3 optical modes are present. Using the
same values, we simulated the energy cut at $Q=(1/2,1/2,1/2)$, displayed in
Fig. 2e, and showing a good agreement with the experimental cut displayed in
Fig. 2b. A powder average has been performed for the high energy branches to
compare with the large q-integration of the experimental cut shown in Fig. 2c.
As can be seen in Fig. 2f, the three experimental modes are reproduced.

Table 2: Values of exchange interactions (in meV) determined using a purely antiferromagnetic effective Heisenberg model $J_{1}$ | 5.5 | $J_{R}$ | 1.3 | $J_{4}$ | 0.15
---|---|---|---|---|---
$J^{\prime}_{1}$ | 19.2 | $J^{\prime}_{2}$ | 18.4 | $J_{5}$ | 2.0
$J_{3}$ | 15.3 | $J_{2}$ | 6.3 | $D$ | -3.9

In the present model, the block order is stabilized by strong
$J^{{}^{\prime}}_{1}$, $J^{{}^{\prime}}_{2}$ and $J_{3}$. The other couplings
are weaker and frustrated. This is at variance with the model proposed in
Mourigal  _et al._ (2015), which assumes the existence of alternating
ferromagnetic and antiferromagnetic couplings along each leg of the ladder.
Note that single ion anisotropy is also relatively strong in our model. This
is at the origin of the strong gap at the magnetic zone center
$q=(1/2,1/2,1/2)$. Based on our previous analysis of the influence of each
coupling, these behaviors lead us to believe that optical modes are related to
intra-blocks and frustrated couplings, while the acoustic mode is associated
with inter-block couplings and strong anisotropy. It is interesting to note
that the energy of this gap nearly coincides with the low energy phonons
observed in infrared spectroscopy around 60 cm-1 Weseloh  _et al._ (2022).
This particular polar phonon mode seems to be affected by the magnetic
ordering below $T_{N}$, suggesting a coupling with the magnetic excitations.
Such coupling between phonon and magnon could indicate the presence of an
electromagnon, the famous excitation bearing electroactive and magnetic
character Petit  _et al._ (2013); Vaunat  _et al._ (2021).

The new set of exchange parameters determined in the present study presents
numerous advantages with respect to the previous proposed model Mourigal  _et
al._ (2015). First, all interactions being antiferromagnetic, the model is
compatible with the Goodenough-Kanamori rule and the bonds angle listed in
Tab. 1. There is also no need to resort to ferromagnetic interaction, an
artificial solution that remains rather difficult to justify physically.
Second, the new set of parameters describes the 3 modes observed at high
energy, and that were not resolved in the previous experiment performed on
powder Mourigal  _et al._ (2015). The strong amplitude of the nearest-neighbor
coupling $J_{3}$ could be partially explained by the chemistry of the
chalcogenide. Selenium orbitals are known to be particularly delocalized in
space, likely responsible for this large second-neighbor interaction. In the
same vein, Sulfur orbitals are expected to be less delocalized, resulting in a
weaker second-neighbor interaction in BaFe2S3. This explains why the nearest-
neighbor frustrated coupling is lifted in this compound, leading to an
antiferromagnetic stripe-like order. Third, the magnetic transition observed
under pressure in Zheng  _et al._ (2022b) finds a natural explanation in the
present model. Indeed, as shown by Monte Carlo calculations, the magnetic
propagation wavevector is quite sensitive to the $J_{1}$/$J_{3}$ ratio, with
the stripe-like phase being found when $J_{3}$ decreases. The intuitive
decrease of $J_{3}$ upon pressure is further supported by first-principle
calculations performed on the parent compounds of BaFe2Se3, namely FeSe and
FeTe Glasbrenner  _et al._ (2015). The authors claim that under pressure, the
ground state of these systems transitions to an antiferromagnetic stripe
order, as is the case for BaFe2Se3 at 4 GPa.

However, detailed calculations of the exchange couplings would be helpful :
their evolution as a function of pressure is related to the structure in a
non-trivial way. Fourth, the model detailed in this work explains the stripe
antiferromagnetic order observed in all other members of the series with
similar structure Caron  _et al._ (2012); Hawai  _et al._ (2015); Takahashi
_et al._ (2015). Finally, the frustration inherent to this model provides a
key ingredient to elucidate the magneto-electric coupling in this system.
Indeed, it has been shown that this compound is multiferroic Zheng  _et al._
(2020); Du  _et al._ (2020); Weseloh  _et al._ (2022), and exhibits a
significant structural change across the Néel temperature, with an increased
deviation from the non-polar average space group Zheng  _et al._ (2020). While
the magnetic structure chirality reported recently suggests that the inverse
Dzyaloshinski-Moriya interaction could be a relevant ingredient at play Zheng
_et al._ (2022a), the presence of frustration suggests that an exchange-
striction mechanism could also be at work. This reinforces the picture of a
strong magneto-elastic coupling in this systems as suggested in ref Weseloh
_et al._ (2022).

In conclusion, on the basis of single crystal inelastic neutron scattering
experiments, is infered an antiferromagnetic Heinsenberg model, providing a
natural explanation for numerous properties observed in this quasi
unidimensional iron based multiferroic compound. This new model can be easily
extended to other members of the family by tuning the parameters to describe
the underlying magnetic orders. It would be interesting to revisit the already
published spin dynamics data Wang  _et al._ (2016, 2017) in the light of this
new model and use it also for future studies of magnetic excitations. This
result finally provides a new starting point to study the remarkable
properties of these systems, among which multiferroicity and
superconductivity.

Aknowledgment : The research performed at ARCS used resources at Spallation
Neutron Source and was supported by the U.S. Department of Energy, Office of
Science, User Facilities Division, and the National Quantum Information
Science Research Centers, Quantum Science Center, operated by Oak Ridge
National Laboratory. This work was financially supported by ANR COCOM
20-CE30-0029 and FRAGMENT 19-CE30-0040. We acknowledge the MORPHEUS platform
at the Laboratoire de Physique des Solides for sample alignment.

## References

  * Zheng  _et al._ (2020) W. Zheng, V. Balédent, M. Lepetit, P. Retailleau, E. Elslande, C. Pasquier, P. Auban-Senzier, A. Forget, D. Colson, and P. Foury-Leylekian, Phys. Rev. B 101, 020101 (2020).
  * Weseloh  _et al._ (2022) M. J. Weseloh, V. Balédent, W. Zheng, M. Verseils, P. Roy, J. B. Brubach, D. Colson, A. Forget, P. Foury-Leylekian, and M.-B. Lepetit, Journal of Physics: Condensed Matter 34, 255402 (2022).
  * Caron  _et al._ (2012) J. M. Caron, J. R. Neilson, D. C. Miller, K. Arpino, A. Llobet, and T. M. McQueen, Phys. Rev. B 85, 180405 (2012).
  * Zheng  _et al._ (2022a) W. G. Zheng, V. Balédent, E. Ressouche, V. Petricek, D. Bounoua, P. Bourges, Y. Sidis, A. Forget, D. Colson, and P. Foury-Leylekian, Phys. Rev. B 106, 134429 (2022a).
  * Zheng  _et al._ (2022b) W.-G. Zheng, V. Balédent, C. V. Colin, F. Damay, J.-P. Rueff, A. Forget, D. Colson, and P. Foury-Leylekian, Communications Physics 5, 1 (2022b).
  * Ying  _et al._ (2017) J. Ying, H. Lei, C. Petrovic, Y. Xiao, and V. V. Struzhkin, Phys. Rev. B 95, 241109 (2017).
  * Mourigal  _et al._ (2015) M. Mourigal, S. Wu, M. B. Stone, J. R. Neilson, J. M. Caron, T. M. McQueen, and C. L. Broholm, Phys. Rev. Lett. 115, 047401 (2015).
  * Herbrych  _et al._ (2018) J. Herbrych, N. Kaushal, A. Nocera, G. Alvarez, A. Moreo, and E. Dagotto, Nat. Commun. 9, 3736 (2018).
  * Hawai  _et al._ (2015) T. Hawai, Y. Nambu, K. Ohgushi, F. Du, Y. Hirata, M. Avdeev, Y. Uwatoko, Y. Sekine, H. Fukazawa, J. Ma, _et al._ , Physical Review B 91, 184416 (2015).
  * Takahashi  _et al._ (2015) H. Takahashi, A. Sugimoto, Y. Nambu, T. Yamauchi, Y. Hirata, T. Kawakami, M. Avdeev, K. Matsubayashi, F. Du, C. Kawashima, H. Soeda, S. Nakano, Y. Uwatoko, Y. Ueda, T. J. Sato, and K. Ohgushi, Nat. Mater. 14, 1008 (2015).
  * (11) ., See Supplemental Material at [URL will be inserted by publisher] for sample characterization, data analysis and Monte Carlo calculation details. ., . (.).
  * Petit (2011) S. Petit, École thématique de la Société Française de la Neutronique 12, 105 (2011).
  * Petit  _et al._ (2013) S. Petit, V. Balédent, C. Doubrovsky, M. B. Lepetit, M. Greenblatt, B. Wanklyn, and P. Foury-Leylekian, Phys. Rev. B 87, 140301 (2013).
  * 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).
  * Glasbrenner  _et al._ (2015) J. K. Glasbrenner, I. I. Mazin, H. O. Jeschke, P. J. Hirschfeld, R. M. Fernandes, and R. Valentí, Nature Physics 11, 953 (2015).
  * Du  _et al._ (2020) K. Du, L. Guo, J. Peng, X. Chen, Z.-N. Zhou, Y. Zhang, T. Zheng, Y.-P. Liang, J.-P. Lu, Z.-H. Ni, S.-S. Wang, G. Van Tendeloo, Z. Zhang, S. Dong, and H. Tian, npj Quantum Materials 5, 49 (2020).
  * Wang  _et al._ (2016) M. Wang, M. Yi, S. Jin, H. Jiang, Y. Song, H. Luo, A. D. Christianson, C. de la Cruz, E. Bourret-Courchesne, D.-X. Yao, D. H. Lee, and R. J. Birgeneau, Phys. Rev. B 94, 041111 (2016).
  * Wang  _et al._ (2017) M. Wang, S. J. Jin, M. Yi, Y. Song, H. C. Jiang, W. L. Zhang, H. L. Sun, H. Q. Luo, A. D. Christianson, E. Bourret-Courchesne, D. H. Lee, D. X. Yao, and R. J. Birgeneau, Physical Review B 95, 1 (2017).