# Determination of energy-dependent neutron backgrounds using shadow bars

S. N. Paneru<EMAIL_ADDRESS>Facility for Rare Isotope Beams, East
Lansing, Michigan 48824, USA  K. W. Brown  Facility for Rare Isotope Beams,
East Lansing, Michigan 48824, USA  Department of Chemistry, Michigan State
University, East Lansing, MI 48824, USA  F.C.E Teh  Facility for Rare Isotope
Beams, East Lansing, Michigan 48824, USA  Department of Physics and Astronomy,
Michigan State University, East Lansing, MI 48824, USA  K. Zhu  National
Superconducting Cyclotron Laboratory, Michigan State University, East Lansing,
Michigan 48824, USA  Department of Physics and Astronomy, Michigan State
University, East Lansing, MI 48824, USA  M.B. Tsang  Facility for Rare Isotope
Beams, East Lansing, Michigan 48824, USA  Department of Physics and Astronomy,
Michigan State University, East Lansing, MI 48824, USA  D. Dell’Aquila
National Superconducting Cyclotron Laboratory, Michigan State University, East
Lansing, Michigan 48824, USA  Z. Chajecki  Department of Physics, Western
Michigan University, Kalamazoo, MI 49008, USA  W.G. Lynch  Facility for Rare
Isotope Beams, East Lansing, Michigan 48824, USA  Department of Physics and
Astronomy, Michigan State University, East Lansing, MI 48824, USA  S. Sweany
National Superconducting Cyclotron Laboratory, Michigan State University, East
Lansing, Michigan 48824, USA  Department of Physics and Astronomy, Michigan
State University, East Lansing, MI 48824, USA  C.Y. Tsang  National
Superconducting Cyclotron Laboratory, Michigan State University, East Lansing,
Michigan 48824, USA  Department of Physics and Astronomy, Michigan State
University, East Lansing, MI 48824, USA  A.K. Anthony  Facility for Rare
Isotope Beams, East Lansing, Michigan 48824, USA  Department of Physics and
Astronomy, Michigan State University, East Lansing, MI 48824, USA  J. Barney
National Superconducting Cyclotron Laboratory, Michigan State University, East
Lansing, Michigan 48824, USA  Department of Physics and Astronomy, Michigan
State University, East Lansing, MI 48824, USA  J. Estee  National
Superconducting Cyclotron Laboratory, Michigan State University, East Lansing,
Michigan 48824, USA  Department of Physics and Astronomy, Michigan State
University, East Lansing, MI 48824, USA  I. Gasparic  Division of Experimental
Physics, Rudjer Boskovic Institute, Zagreb, Croatia  G. Jhang  National
Superconducting Cyclotron Laboratory, Michigan State University, East Lansing,
Michigan 48824, USA  O.B. Khanal  Department of Physics, Western Michigan
University, Kalamazoo, MI 49008, USA  J. Mandredi  National Superconducting
Cyclotron Laboratory, Michigan State University, East Lansing, Michigan 48824,
USA  Department of Physics and Astronomy, Michigan State University, East
Lansing, MI 48824, USA  C.Y. Niu  National Superconducting Cyclotron
Laboratory, Michigan State University, East Lansing, Michigan 48824, USA  R.S.
Wang  National Superconducting Cyclotron Laboratory, Michigan State
University, East Lansing, Michigan 48824, USA  J.C. Zamora  Facility for Rare
Isotope Beams, East Lansing, Michigan 48824, USA

###### Abstract

Understanding the neutron background is essential for determining the neutron
yield from nuclear reactions. In the analysis presented here, the shadow bars
are placed in front of neutron detectors to determine the energy dependent
neutron background fractions. The measurement of neutron spectra with and
without shadow bars is important to determine the neutron background more
accurately. The neutron background, along with its sources and systematic
uncertainties, are explored with a focus on the impact of background models
and their dependence on neutron energy.

###### pacs:

##  I Introduction

Neutron detection provides essential insights into the nuclear dynamics of a
nuclear reaction and nuclear structure. As the availability of access to
radioactive beams has increased, allowing access to more neutron rich nuclei,
the problem of properly measuring the neutron yields has become increasingly
important.Thus more dedicated neutron arrays Luther  _et al._ (2003);
Valiente-Dobón  _et al._ (2019); Peters  _et al._ (2016); Boretzky  _et al._
(2021); Nakamura and Kondo (2016) have become available in recent years. The
problem of measuring neutron yields is complicated for many reasons such as
neutron background, low neutron detection efficiency, varying behaviour of
neutrons with energy, etc. A critical component in the analysis of neutron
yields is the ability to accurately model and subtract the neutron background
from the experimentally measured yield.

Usually, a bar or block, called a “shadow bar” hereafter, is placed between
the target and the neutron detector to prevent neutrons from taking a direct
path to reach the detector. Comparing the number of neutrons detected in the
blocked and unblocked areas provides an estimate of the neutron background.
This method is limited by the assumption that the estimated background is
independent of neutron energy and the neutron distribution is isotropic. The
previous works on neutron background determination in Refs. Peters  _et al._
(2016); Perdikakis  _et al._ (2012); Heilbronn  _et al._ (2007); Kei  do not
explore the consequences of these assumptions. In this paper, we determine the
energy dependent neutron background using the shadow bar for the Large Area
Neutron Array (LANA) P.D Zecher and A Galonsky and J.J Kruse and S.J Gaff and
J Ottarson and J Wang and F Deák and Á Horváth and Á Kiss and Z Seres and K
Ieki and Y Iwata and H Schelin (1997) constructed at the National
Superconducting Cyclotron Laboratory. This technique can be applied to other
neutron detectors.

##  II EXPERIMENT

The Large Area Neutron Array (LANA) is composed of two large neutron walls,
each possessing 25 independent detection cells that cover a total area of 2 x
2 $\text{m}^{2}$. Each detector cell is a 76′′-long Pyrex glass container,
with a rectangular cross section of height $h$ = 7.62 cm and depth $w$= 6.35
cm, filled with liquid organic scintillator NE-213. To maximize detection
efficiency, LANA features large active detection area, small dead space
between scintillators and minimum amount of inactive mass through which the
neutrons must pass. For heavy ion collision experiments, it is designed to
detect neutron from about 10 MeV to about 300 MeV when it is placed at 4.5 m
from the target. LANA has been used in many previous nuclear experiments, for
example in Ref Coupland (2013); Famiano  _et al._ (2006). The analysis
reported here will focus on a single neutron wall.

Figure 1: Experimental Setup showing the shadow bar. The red rectangular box
indicates the position of shadow bars.

In the experiment, LANA was used to measure neutrons with beams of 40,48Ca at
56,140 MeV/u on targets of 58,64Ni and 112,124Sn. The calibration of LANA was
performed using cosmic rays and an AmBe neutron source as explained in Ref.
Zhu  _et al._ (2020). The neutrons were distinguished from charged particles
using the charged particle veto wall consisting of thin bars of plastic
scintillator that were placed in front of the LANA. The neutron-gamma
separation was performed with value-assigned pulse shape discrimination method
as explained in Ref. Teh  _et al._ (2021). The energy of neutrons is
determined using the time of flight method.

To determine the neutron background, we block some areas of LANA and then
compare the number of neutrons detected in blocked and not blocked areas of
LANA. Four 30-cm-long brass shadow bars are used to block about 2$\%$ of front
LANA detection area. They are placed about 1.7 m from the target, outside the
chamber as shown in Fig. 1. The shadow bars are thick enough to stop neutrons
of energies $\leq$300 MeV that take a direct path to the neutron wall. Figure
2 shows one of the four shadow bars used in the experiment. The bottom bar of
LANA is not used in our experiment because of the shadowing by stand holding
the shadow bars. Figure 3 shows the neutron hit distribution obtained for 48Ca
beam impinging on 124 Sn at E/A = 140 MeV in lab $\theta$-$\phi$ coordinates
on LANA, showing the decreased amount of neutrons in the shadow regions. The
positions A, B, C, and D represent the shadow bar positions as shown in Fig.
1. The neutron bars are counted from bottom to the top. The shadow bars at
positions C and D fully cover the vertical cross section of bar 8 and
partially cover bars 7 and 9, while the shadow bars at positions A and B fully
cover bar 16 and partially cover bars 15 and 17. The shadow bars in position A
and C cover larger angles in $\theta$-$\phi$ distribution while the shadow
bars in positions B and D cover more forward angles. The neutron spectrum is
mostly unaffected by the thin aluminum mount holding the shadow bar as their
interaction probability in this material is small and the energy loss of
neutrons is also negligible. Example spectra for the detected neutron yield as
a function of angle for one of the shadowed detector bars are shown in Fig. 8,
both with and without the shadow bars in place.

Figure 2: Shadow bar used in experiment. Figure 3: The neutron hit pattern
with 4 shadow bar installed in $\theta$-$\phi$ lab coordinates. The 4 square
areas with less yield correspond to the area shadowed by the shadow bars. The
neutron bars are counted from bottom to the top. The shadow bars at the bottom
fully cover bar 8 and partially cover bars 7 and 9, while the shadow bars at
the top fully cover bar 16 and partially cover bars 15 and 17. The reaction
system for this plot is 48Ca +124 Sn at E/A = 140 MeV .

Neutrons produced in the reaction target can also be scattered on the room
walls (or other volume), and later detected in LANA. Any neutrons other than
the ones directly hitting the LANA are considered as background neutrons. In
order to investigate the various sources of neutron multiple scattering, a
Monte-Carlo simulation was performed with the code PHITS Sato  _et al._
(2018). A realistic geometry of the experimental hall, with the most
representative scattering elements, was included in the simulation. Primary
neutrons were generated in the target location by using the experimental
energy distribution. The neutron angles were uniformly distributed in the
4$\pi$ solid angle. Secondary neutron reactions were accounted in the
simulation with the respective nucleon-nucleus models Boudard  _et al._
(2013); Iida  _et al._ (2007). Neutrons were tagged on event-by-event basis to
separate the multiple-scattering component from the time of flight (ToF)
spectra. The simulated ToF spectra for various neutron sources are shown in
Fig. 4. The sharp edges in the total ToF spectrum (around 30 and 90 ns)
correspond to the limits of the energy distribution used to generate the
neutrons (between 20 to 200 MeV). The neutrons scattered off the shadow bars
and hitting LANA could be contributing to the background. However, the
simulation results show that the fraction of neutrons detected from this
process contribute a very small component ($\sim$1/1000) to the total neutron
spectrum (see dotted blue histogram in Fig. 4). Another major source of
background is the scattering of fast neutrons on the concrete volumes of the
experimental hall such as floor, ceiling, walls, etc. The simulation estimated
this type of component to only contribute at small neutron energies (or
equivalently large time of flight) as shown by the dashed red histogram in
Fig. 4. Neutrons produced in the beam dump have a longer ToF and do not
contribute significantly to the neutron background, since they are located at
the low energy region of the distribution. The dominant background background
contribution thus arise from secondary neutrons resulting from the charged
particle induced reactions on various components of the experimental setup.
The present simulation has a limitation to quantify realistically this
contribution given that the number of reaction channels that would produce
neutrons is relatively large.

Figure 4: The time of flight spectra deduced from the Monte Carlo simulation
of neutron scattering. The blue dotted histogram shows the contribution of
neutron scattering of the shadow bar. The red dashed histogram shows the total
neutron scattering contribution from roof, walls and beam dump. The solid
black histogram shows the neutron spectra observed from simulation.

##  III Analysis

Figure 5: Energy independent shadow bar fitting. Top: The total neutron
spectra is fitted using Eq. (1) for bar 16 for 48Ca + 124Sn system at beam
energy of 140 MeV/u. Bottom: The fit spectra is scaled by the shape of the
spectra around the shadowed regions to determine the neutron background. The
dashed and dotted lines are explained in the text.

Ideally, the spectrum of a shadowed neutron wall element would include a
square gap with a flat bottom, where the bottom of the gap shows the
background, and outside the gap is the result of both background and signal.
In reality, the edges of the gap are smeared due to the position resolution of
the detector as evident from Fig. 8. This is represented by Gaussian
convoluted square gap characterized by the error function, $\mathrm{erf}$. So,
we assume the fit function is characterized by

$\begin{split}F(x)=\mathcal{A}*\left(1.0+\mathrm{erf}\left(-\frac{x-x_{0}}{\sqrt{2}\sigma_{0}}\right)+\mathrm{erf}\left(-\frac{x-x_{1}}{\sqrt{2}\sigma_{1}}\right)\right)\\\
+(b+c*x),\end{split}$ (1)

where
$\mathcal{A}*\left(\mathrm{erf}\left(-\frac{x-x_{0}}{\sigma_{0}}\right)+\mathrm{erf}\left(-\frac{x-x_{1}}{\sigma_{1}}\right)\right)$
is the component representing the area between $x_{0}$ and $x_{1}$ Coupland
(2013). $x_{0}$ and $x_{1}$ are the position of edges of the shadow bar as
shown in Fig. 5, with associated standard deviation $\sigma_{0}$ and
$\sigma_{1}$ coming from the LANA position resolution. The $\sigma_{0}$ and
$\sigma_{1}$ values are varied between 0.31 and 0.45 corresponding to the
position resolution for bars 8 and 16 in LANA as reported in Ref. Zhu (2020).
The $\mathcal{A}+b+c*x$ component, referred as “kin” hereafter, represents the
expected shape of the neutron spectrum in the absence of the shadow bar
Coupland (2013). The neutron spectra for a completely shadowed bar is fitted
using Eq. (1), where $\chi^{2}$ is minimized. The result of the fit for bar 16
for 48Ca + 124Sn system at beam energy of 140 MeV/u is shown in Fig. 5. The
best fit line is represented by the black dashed line in the figure. To get
the background fraction Eq. (1) is scaled by $\mathcal{A}+b+c*x$ as shown in
bottom panel of Fig. 5. The background fraction is then determined at the
center between the $x_{0}$ and $x_{1}$ points. The uncertainty in the fitting
parameters is propagated to determine the upper limit and lower limits of the
best fit line represented by the red dashed and red dotted lines,
respectively, in upper panel of Fig. 5, and their corresponding “kin”
components represented by dashed pink and dotted cyan lines, respectively. The
error bar in neutron background is then determined as the difference between
maximum and minimum values inferred as shown in bottom panel of Fig. 6. This
method is referred as “Method 1” hereafter in the text.

The neutron background is dependent on neutron energy or equivalently the time
of flight. The background fraction in the past has been determined without
taking into consideration the energy dependence of neutrons as shown in Fig.
5, this is equivalent to method 1. In this work, we extend method 1 by
dividing the neutron spectra is divided into equal energy bins and determining
the background fraction for each bin. For example, Fig. 6 shows the background
scattering fraction determined for neutrons with energies $E_{n}$=30 $\pm$ 10
MeV. The energy bin width of 20 MeV was chosen such that we have sufficient
statistics for the fitting process. The figures shown throughout this paper is
for 48Ca + 124Sn system at beam energy of 140 MeV/u for bar 16 unless
otherwise stated.

Figure 6: Energy dependent shadow bar fitting. The background scattering
fraction determined for neutron energies, $E_{n}$=30 $\pm$ 10 MeV for bar 16
for 48Ca + 124Sn system at beam energy of 140 MeV/u. The lines are as
described in Fig. 5.

##  IV Systematics

The error in background fraction determined from the algorithm discussed in
Sec. III is the statistical error. The dip in the neutron shadowed region
might not be characterized by the flat bottom as shown in Fig. 6. So, another
way of calculating the errors in the background fraction is to calculate the
standard deviation of the background fraction within the dip region. The
error-bars obtained from this approach is much smaller than the statistical
error-bars except in the case where we are limited by the statistics.

###  IV.1 Model dependent systematic uncertainty

In another approach to determine the background fraction, one could envision
scaling the neutron spectra with shadow bars by the neutron spectra without
shadow bars. Fitting the region around the gap with kin produces the expected
shape of the spectrum in the absence of the shadow bar. The spectrum is then
scaled by this expected shape of the spectrum in absence of shadow bar. The
scaled spectrum is then fitted using the fitting function given by,

$\begin{split}F(x)=S*\left(1.0+\mathrm{erf}\left(-\frac{x-x_{0}}{\sqrt{2}\sigma_{0}}\right)+\mathrm{erf}\left(-\frac{x-x_{1}}{\sqrt{2}\sigma_{1}}\right)\right)\\\
+Bg,\end{split}$ (2)

where $S$ is the signal fraction and $Bg$ is the background fraction. This
model assumes a constant neutron background. The uncertainty in the background
fraction is the statistical uncertainty in $Bg$.

The background fractions determined using the above described two methods
produces identical results, which is not surprising as both of the methods are
based on the assumption that the angular distribution of neutrons is explained
by a straight line for a small range of angles. However, the validity of this
assumption needs to be verified. Since the two methods discussed so far, we
used only the data points from Method 1 to characterize the energy dependent
neutron background fraction.

Instead of scaling the neutron spectra around the shadowed regions with a
straight line, we could use the neutron spectra from the experiment where the
shadow bars are not used. This is termed as the $``$spectral ratio” method
hereafter, where the ratio of the spectra with and without shadow bar is taken
and the scaled spectra is fitted with the fitting function characterized by
Eq. (2) to get the background fraction, as shown in Fig. 8 and Fig. 8. Figure
8 shows the observed neutron spectra with and without the shadow bars,
represented by the blue and red histograms respectively. The blue histogram is
then divided by the red histogram to get the scaled neutron spectra,
represented by the histogram in Fig. 8, which is then fitted using Eq. (2) to
determine the background fraction. The background fractions extracted from
this method are shown in Fig. 9 for the shadowed positions along with the
background fractions determined from Method 1 for comparison. The background
fractions determined by the spectra ratio method are systematically lower than
the ones determined using Method 1. This highlights the fact that
understanding the angular distribution of neutrons is important to determine
the background fraction. However, running experiments with and without shadow
bars might not be practically feasible. In the experiment discussed in
Section. II, there is not enough data to follow the spectral ratio method for
all the systems.

[]

Figure 7: The neutron spectrum, for $E_{n}$=30 $\pm$ 10 MeV for bar 16, with
and without shadow bars are represented by blue and red histograms
respectively. The red histogram is scaled by a factor such that the counts in
$37\leq\theta\leq 42$ are same for both histograms.

[]

Figure 8: The ratio of the neutron spectra, for $E_{n}$=30 $\pm$ 10 MeV for
bar 16, with and without shadow bars for 48Ca + 124Sn system at beam energy of
140 MeV/u fitted with the functional form characterized by Eq. (2)

The discussion so far about the systematic uncertainties on neutron background
have been focused on the method of its determination. The additional
contribution comes from the position where the background is being measured.
The background fraction is position dependent as evident in Fig. 9 which is
not surprising due to the fact that the neutron angular distribution is not
flat. The background fraction increases with decreasing neutron energy. The
energy dependence of neutron background is also verified by the Monte Carlo
Simulation. We used the exponential decay function
[$F(E)=\text{Amp}*\text{exp}(-E/\tau)$] to characterize this behaviour. The
neutron background measured at each position A, B, C and D are fitted with the
exponential decay functions are represented in Fig. 9. In order to
characterize the background as a function of energy for a given system, we
used the combined fit to the background fractions determined via method 1
shown in Fig. 9 and come up with a single set of parameters (Amp,$\tau$) with
relaxed error bars to describe the positional dependence. It was found that
the 3-$\sigma$ error band in the combined fit will explain most of the data
points as shown in Fig. 10. The relative error in the determined background is
calculated. The systematic uncertainties coming from fitting the backward and
forward angles separately is negligible compared to the systematic uncertainty
due to the choice of models. The fits to all data points from the spectral
ratio method using the exponential decay function will yield a systematically
lower background fraction functional curve as represented by the black dashed
lines in Fig. 10 with similar 3-$\sigma$ bands. The systematic uncertainty due
to different models is a function of neutron energy which increases with
increasing neutron energy. The systematic uncertainty at neutron energies of
125 MeV due to the choice of model is a maximum of 15$\%$.

Figure 9: Positional dependence of neutron background. The neutron background
determined using both Method 1 and spectra ratio method, represented by
circles and squares respectively, for 48Ca + 124Sn system at beam energy of
140 MeV/u. A, B, C, D represents the position of shadow bars. The upper row is
for bar 16 and bottom row is for bar 8, respectively. The parameters shown in
each figure are obtained from individual fits (dashed red line) of neutron
background to data from Method 1 at each position respectively.  Figure 10:
Combined fit to neutron background data using the exponential decay function.
The shaded region shows the 3-$\sigma$ error band. The cyan dashed line
represents the best fit line to all data points obtained from method 1 using
the exponential decay function. The black dotted line represents the best fit
line to all data points obtained from the spectral ratio method using the
exponential decay function. The blue solid line represents the fit to all data
points obtained from method 1 using the first order polynomial fitting
function which yields a reduced $\chi^{2}$ value of 1.06 compared to a reduced
$\chi^{2}$ value of 1.30 for exponential decay function represented by dashed
cyan line.

The choice of exponential decay function is motivated by the fact that the
background fraction should not be negative at large neutron energies. If we
choose the first order polynomial for the combined fit, the background
fraction at large neutron energies is negative which is not physically
possible, but could be interpreted as no neutron background at those higher
energies. The first order polynomial gives the best fit for all the systems
compared to the exponential decay as shown in Fig. 10. The reduced $\chi^{2}$
for exponential fitting function is 1.30 compared to 1.06 that of first order
polynomial. However, the systematic uncertainty due to choice in exponential
decay function versus the first order polynomial is less than 5$\%$ where the
background fraction data points are available and is not a dominant
contributor to the overall systematic uncertainty.

Although the systematic uncertainty at large neutron energies is higher, but
the fact that the background fraction at these energies is much smaller, the
overall uncertainty to the neutron spectra at large energies is less than
6$\%$ for most of the systems we analyzed.

##  V Conclusion

In summary, we describe the method of determining the background fraction
using the shadow bars in the heavy-ion collision experiment. The systematic
uncertainty of the method is described and found that understanding the
neutron angular distribution is important to determine the background
fraction. The position or angle dependence of the neutron background is
observed. We would recommend measuring the neutrons with and without the
shadow bars to determine the neutron background fraction more accurately with
spectral ratio method. The non-linear angular distributions of neutrons
impacts the neutrons background determined using method 1.

##  VI ACKNOWLEDGEMENTS

This work was supported in part by the U.S. National Science Foundation under
Grants PHY-1565546, PHYS-1712832, and PHY-2110218; in part by the U.S.
Department of Energy (the Office of Science) under Grant DE-NA0003908; and in
part by the National Research Foundation of Korea (NRF) under Grant
2016K1A3A7A09005578, Grant 2018R1A5A1025563, and Grant 2013M7A1A1075764.

## References

  * Luther  _et al._ (2003) B. Luther, T. Baumann, M. Thoennessen, J. Brown, P. DeYoung, J. Finck, J. Hinnefeld, R. Howes, K. Kemper, P. Pancella, G. Peaslee, W. Rogers, and S. Tabor, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 505, 33 (2003), proceedings of the tenth Symposium on Radiation Measurements and Applications.
  * Valiente-Dobón  _et al._ (2019) J. Valiente-Dobón, G. Jaworski, A. Goasduff, F. Egea, V. Modamio, T. Hüyük, A. Triossi, M. Jastrząb, P. Söderström, A. Di Nitto, G. de Angelis, G. de France, N. Erduran, A. Gadea, M. Moszyński, J. Nyberg, M. Palacz, R. Wadsworth, R. Aliaga, C. Aufranc, M. Bézard, G. Baulieu, E. Bissiato, A. Boujrad, I. Burrows, S. Carturan, P. Cocconi, G. Colucci, D. Conventi, M. Cordwell, S. Coudert, J. Deltoro, L. Ducroux, T. Dupasquier, S. Ertürk, X. Fabian, V. González, A. Grant, K. Hadyńska-Klęk, A. Illana, M. Jurado-Gomez, M. Kogimtzis, I. Lazarus, L. Legeard, J. Ljungvall, G. Pasqualato, R. Pérez-Vidal, A. Raggio, D. Ralet, N. Redon, F. Saillant, B. Sayğı, E. Sanchis, M. Scarcioffolo, M. Siciliano, D. Testov, O. Stezowski, M. Tripon, and I. Zanon, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 927, 81 (2019).
  * Peters  _et al._ (2016) W. Peters, S. Ilyushkin, M. Madurga, C. Matei, S. Paulauskas, R. Grzywacz, D. Bardayan, C. Brune, J. Allen, J. Allen, Z. Bergstrom, J. Blackmon, N. Brewer, J. Cizewski, P. Copp, M. Howard, R. Ikeyama, R. Kozub, B. Manning, T. Massey, M. Matos, E. Merino, P. O’Malley, F. Raiola, C. Reingold, F. Sarazin, I. Spassova, S. Taylor, and D. Walter, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 836, 122 (2016).
  * Boretzky  _et al._ (2021) K. Boretzky, I. Gašparić, M. Heil, J. Mayer, A. Heinz, C. Caesar, D. Kresan, H. Simon, H. Törnqvist, D. Körper, G. Alkhazov, L. Atar, T. Aumann, D. Bemmerer, S. Bondarev, L. Bott, S. Chakraborty, M. Cherciu, L. Chulkov, M. Ciobanu, U. Datta, E. De Filippo, C. Douma, J. Dreyer, Z. Elekes, J. Enders, D. Galaviz, E. Geraci, B. Gnoffo, K. Göbel, V. Golovtsov, D. Gonzalez Diaz, N. Gruzinsky, T. Heftrich, H. Heggen, J. Hehner, T. Hensel, E. Hoemann, M. Holl, A. Horvat, . Horváth, G. Ickert, D. Jelavić Malenica, H. Johansson, B. Jonson, J. Kahlbow, N. Kalantar-Nayestanaki, A. Kelić-Heil, M. Kempe, K. Koch, N. Kozlenko, A. Krivshich, N. Kurz, V. Kuznetsov, C. Langer, Y. Leifels, I. Lihtar, B. Löher, J. Machado, N. Martorana, K. Miki, T. Nilsson, E. Orischin, E. Pagano, S. Pirrone, G. Politi, P.-M. Potlog, A. Rahaman, R. Reifarth, C. Rigollet, M. Röder, D. Rossi, P. Russotto, D. Savran, H. Scheit, F. Schindler, D. Stach, E. Stan, J. Stomvall Gill, P. Teubig, M. Trimarchi, L. Uvarov, M. Volknandt, S. Volkov, A. Wagner, V. Wagner, S. Wranne, D. Yakorev, L. Zanetti, A. Zilges, and K. Zuber, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1014, 165701 (2021).
  * Nakamura and Kondo (2016) T. Nakamura and Y. Kondo, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 376, 156 (2016), proceedings of the XVIIth International Conference on Electromagnetic Isotope Separators and Related Topics (EMIS2015), Grand Rapids, MI, U.S.A., 11-15 May 2015.
  * Perdikakis  _et al._ (2012) G. Perdikakis, M. Sasano, S. M. Austin, D. Bazin, C. Caesar, S. Cannon, J. Deaven, H. Doster, C. Guess, G. Hitt, J. Marks, R. Meharchand, D. Nguyen, D. Peterman, A. Prinke, M. Scott, Y. Shimbara, K. Thorne, L. Valdez, and R. Zegers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 686, 117 (2012).
  * Heilbronn  _et al._ (2007) L. Heilbronn, C. J. Zeitlin, Y. Iwata, T. Murakami, H. Iwase, T. Nakamura, T. Nunomiya, H. Sato, H. Yashima, R. M. Ronningen, and K. Ieki, Nuclear Science and Engineering 157, 142 (2007), https://doi.org/10.13182/NSE07-A2719 .
  * (8) “Symposium on Nuclear Data 2020,” https://indico2.riken.jp/event/3359/contributions/15449/attachments/9807/14096/38_SND2020_FFAG_Com.pdf, accessed: 2021-10-06.
  * P.D Zecher and A Galonsky and J.J Kruse and S.J Gaff and J Ottarson and J Wang and F Deák and Á Horváth and Á Kiss and Z Seres and K Ieki and Y Iwata and H Schelin (1997) P.D Zecher and A Galonsky and J.J Kruse and S.J Gaff and J Ottarson and J Wang and F Deák and Á Horváth and Á Kiss and Z Seres and K Ieki and Y Iwata and H Schelin, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 401, 329 (1997).
  * Coupland (2013) D. D. S. Coupland,  _Probing the nuclear symmetry with heavy ion collisions_ , Ph.D. thesis, Michigan State University (2013).
  * Famiano  _et al._ (2006) M. A. Famiano, T. Liu, W. G. Lynch, M. Mocko, A. M. Rogers, M. B. Tsang, M. S. Wallace, R. J. Charity, S. Komarov, D. G. Sarantites, L. G. Sobotka, and G. Verde, Phys. Rev. Lett. 97, 052701 (2006).
  * Zhu  _et al._ (2020) K. Zhu, M. Tsang, D. Dell’Aquila, K. Brown, Z. Chajecki, W. Lynch, S. Sweany, F. Teh, C. Tsang, C. Anderson, A. Anthony, J. Barney, J. Crosby, J. Estee, I. Gasparic, G. Jhang, O. Khanal, S. Kodali, J. Manfredi, C. Niu, and R. Wang, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 967, 163826 (2020).
  * Teh  _et al._ (2021) F. C. E. Teh, J.-W. Lee, K. Zhu, K. W. Brown, Z. Chajecki, W. G. Lynch, M. B. Tsang, A. Anthony, J. Barney, D. Dell’Aquila, J. Estee, B. Hong, G. Jhang, O. B. Khanal, Y. J. Kim, H. S. Lee, J. W. Lee, J. Manfredi, S. H. Nam, C. Y. Niu, J. H. Park, S. Sweany, C. Y. Tsang, R. Wang, and H. Wu, IEEE Transactions on Nuclear Science 68, 2294 (2021).
  * Sato  _et al._ (2018) T. Sato, Y. Iwamoto, S. Hashimoto, T. Ogawa, T. Furuta, S. ichiro Abe, T. Kai, P.-E. Tsai, N. Matsuda, H. Iwase, N. Shigyo, L. Sihver, and K. Niita, Journal of Nuclear Science and Technology 55, 684 (2018), https://doi.org/10.1080/00223131.2017.1419890 .
  * Boudard  _et al._ (2013) A. Boudard, J. Cugnon, J.-C. David, S. Leray, and D. Mancusi, Phys. Rev. C 87, 014606 (2013).
  * Iida  _et al._ (2007) K. Iida, A. Kohama, and K. Oyamatsu, Journal of the Physical Society of Japan 76, 044201 (2007), https://doi.org/10.1143/JPSJ.76.044201 .
  * Zhu (2020) K. Zhu,  _Measuring neutrons in heavy ion collisions_ , Ph.D. thesis, Michigan State University (2020).


# First-principle study of multiple metastable charge ordering states in
La1/3Sr2/3FeO3

Nam Nguyen1, Alex Taekyung Lee2,3, Vijay Singh4, Anh T. Ngo2,3, and Hyowon
Park1,2 1Department of Physics, University of Illinois at Chicago, Chicago, IL
60607, USA
2Materials Science Division, Argonne National Laboratory, Argonne, IL, 60439,
USA
3Department of Chemical engineering, University of Illinois at Chicago,
Chicago, IL 60608, USA
4GITAM School of Science, Bangalore 561203, India

###### Abstract

La doped SrFeO3, La1/3Sr2/3FeO3, exhibits a metal-to-insulator transition
accompanied by both antiferromagnetic and charge ordering states along with
the Fe-O bond disproportionation below a critical temperature near 200K.
Unconventionally slow charge dynamics measured in this material near the
critical temperature [Nature Communications, 9 1799 (2018)] shows that its
excited charge ordering states can exhibit novel electronic structures with
nontrivial energy profiles. Here, we reveal possible metastable states of
charge ordering structures in La1/3Sr2/3FeO3 using the first-principle and
climbing image nudged elastic band methods. In the strong correlation regime,
La1/3Sr2/3FeO3 is an antiferromagnetic insulator with a charge ordering state
of the big-small-big pattern, consistent with the experimental measurement of
this material at the low temperature. As the correlation effect becomes weak,
we find at least two possible metastable charge ordering states with the
distinct Fe-O bond disproportionation. Remarkably, a ferroelectric metallic
state emerges with the small energy barrier of $\sim$7meV, driven by a
metastable CO state of the small-medium-big pattern. The electronic structures
of these metastable charge ordering states are noticeably different from those
of the ground-state. Our results can provide an insightful explanation to
multiple metastable charge ordering states and the slow charge dynamics of
this and related oxide materials.

##  I Introduction

Charge ordering (CO) or charge density wave (CDW) is an intriguing material
property driven by a spontaneous symmetry breaking of the periodicity in
crystals. In strongly correlated materials, the charge degree of freedom is
typically coupled to other degrees of freedom including spin, orbital, or
lattice. While the origin of CDW can be purely electronic and the electronic
correlation plays an important role, it is often accompanied by structural
distortions such as the bond-order or the Peierls transition, possibly leading
to ferroelectricity. Indeed, the combination of CDW, spin density wave (SDW),
and the bond order has been proposed as the mechanism of ferroelectricity van
den Brink and Khomskii (2008); Efremov  _et al._ (2005); Wan  _et al._ (2016).

La1/3Sr2/3FeO3 (LSFO) is a transition metal oxide with a perovskite structure
undergoing a weakly first-order transition at a temperature, $T$=200K, from a
paramagnetic metallic state with the average valence state of
Fe${}^{3.67+}(d^{4.3})$ at a high temperature to an antiferromagnetic (AFM)
insulating state with a CO sate of Fe${}^{3+}(d^{5})$: Fe${}^{5+}(d^{3})$=2:1
at a low temperature Park  _et al._ (1999). Structural properties of LSFO with
or without CO phases have been characterized experimentally using X-ray
diffraction, neutron diffraction, and electron microscopy. The studies of
X-ray and neutron diffraction Battle  _et al._ (1990); Sabyasachi  _et al._
(2012); Dann  _et al._ (1994); Yang  _et al._ (2003); Takano  _et al._ (1981,
1981); Okamoto  _et al._ (2010); Li  _et al._ (2018) showed that bulk LSFO
forms a rhombohedral structure in the space group of $R\bar{3}c$ (see Fig. 1)
with the lattice constants $a=5.47$ Å and $c=13.35$ Å . A sign of the CDW
spanning the periodicity of three Fe ions accompanied by SDW with a
periodicity of six Fe ions was measured along the pseudocubic $[111]$
direction, but there was no clear evidence of structural distortions. Later,
the electron microscopy study by Li et al. Li  _et al._ (1997) revealed a
structural distortions along the pseudocubic [111] direction in the real space
upon the CDW transition. Finally, the neutron diffraction studies by
Sabyasachi et al. Sabyasachi  _et al._ (2012) and Yang et al.Yang  _et al._
(2003) also showed a possibility of the meta-stable CO state due to multiple
neutron peaks below the critical temperature.

Electronic properties of LSFO at the low-temperature CO phase haven been
characterized by various experiments. The study of optical spectroscopy by
Ishikawa et al. Ishikawa  _et al._ (1998) showed the optical gap of LSFO was
about 0.13 eV at low temperature. The studies of Mössbauer spectrocopy
Sabyasachi  _et al._ (2012); Dann  _et al._ (1994); Yang  _et al._ (2003);
Takano  _et al._ (1981); Chernyshev  _et al._ (2014); Battle  _et al._ (1990);
Li  _et al._ (2018); Shimony and Knudsen (1966) captured two kinds of Fe ions
with different hyperfine fields, confirming the charge disproportionation
below the critical temperature. Recent ultrafast X-ray measurement in LSFO by
Zhu et al has shown that the noticeable slowdown occurs during the relaxation
of CO near the critical temperature Zhu  _et al._ (2018). They argued that the
photoexcitation due to an ultrafast pump can drive a ground state of
La1/3Sr2/3FeO3 into metastable states with different spin/charge orderings,
which can be the origin of slowdown in the relaxation process. According to
Yamamoto et al., Yamamoto  _et al._ (2022) these metastable or transient
states are the CO in sequence of Fe4+Fe3+Fe4+. However, the magnetic moments
as well as the spin states, i.e. high spin (HS) or low spin (LS), of these
Fe4+ and Fe3+ ions were unknown. In general, the slow dynamics of CO can be
originated from the multiple meta-stable CO states accessible during the
relaxation process.

Unlike those various experimental characterizations, theoretical studies of
LSFO have been rather limited. The Hartree-Fock study by Matsuno et al Matsuno
_et al._ (2002) captured an energy gap of $0.14$ eV , which was in a good
agreement with the experimental gap at low temperature. The first-principle
study of density functional theory plus the Hubbard $U$ (DFT$+U$) by Zhu et
al. Zhu  _et al._ (2018) and Saha-Dasgupta et al. Saha-Dasgupta  _et al._
(2005) verified the presence of structural modulation or oxygen breathing
distortions accompanied by CO of Fe ions in a sequence of Fe3+Fe5+Fe3+. They
also found that another sequence of CO is possible, namely Fe4+Fe3+Fe4+. These
CO states are strongly coupled to the spin states as the Fe ion with a larger
charge state shows the high-spin state with the Fe-O bond elongation. Finally,
the possibility of ferroelectricity in La1/3Sr2/3FeO3 was pointed out by Park
et al. by rearranging the La/Sr layers Park  _et al._ (2019). Nevertheless,
the effect of electronic correlations on the stability of CO states and the
emergence of novel metastable states, such as ferroelectricity, which can be
accessible from the photo-excitation experiment, have not been studied from
first-principle.

In this work, we study the effect of electron correlations on structural and
electronic properties of LSFO having the strong charge-spin-lattice coupling
by adopting the first-principle DFT+$U$ method. In particular, we explore
possible meta-stable CO phases driven by a new pattern of structural
distortions by adopting the climbing image nudged elastic band (CINEB) method
along with DFT+$U$. Remarkably, we find a new electronic phase in LSFO
exhibiting the ferroelectricity driven by a small-medium-big CO pattern and a
distinct Fe-O bond disproportionation with the small-medium-big magnetic
moments. This new meta-stable phase has almost the degenerate energy compared
to previously known CO phases with a small energy barrier of $5\sim 7$ meV,
implying the promising tunability of this material as a future electronic
device.

##  II METHODS

###  II.1 First-principle calculation

To perform the structural relaxation and the band structure calculations of
LSFO, we adopt DFT+$U$ Liechtenstein  _et al._ (1995) based on the projected-
augmented wave (PAW) method Blöchl (1994) as implemented in the Vienna ab
initio simulation package (VASP) Kresse and Furthmüller (1996); Kresse and
Joubert (1999). The exchange-correlation energy functional was treated using
generalized gradient approximation (GGA) by adopting the Perdew-Burke-
Ernzerhof (PBE) functional Perdew  _et al._ (1997). The cutoff energy for the
plane-wave basis was used as 600 eV, and the Gamma-centered
8$\times$8$\times$2 $k-$point mesh Monkhorst and Pack (1976) was used for all
calculations. For structural relaxations, the Hellmann-Feynman force Feynman
(1939) on each atom was set to be smaller than 0.01 eV/Å. To treat the
correlation effect of Fe $d$ orbitals, we impose the Hubbard $U$ and the
Hund’s coupling $J$ within DFT+U.

As noted from the previous study of Ref. Zhu  _et al._ , 2018, two distinct CO
structures (CO1 and CO3, see Fig. 2) can be obtained in LSFO by relaxing the
crystal structure imposing different $U$ values on the Fe ions. For the
ground-state CO1 structure, we used $U$=5eV and $J$=1eV, while a distinct CO3
structure is obtained using $U$=3eV and $J$=0.6eV. Both the crystal shape and
ionic positions were relaxed during the structural relaxation, while the
crystal volume of LSFO was fixed to 329.24 Å3. To obtain a meta-stable CO
phase (CO2, see Fig. 2), we adopt the CINEB method along with the DFT+$U$
using $U$=3.62eV (see Sec. III.2). We also explore the effect of $U$ values
($J=0.2U$) on the stability of different CO phases (see Sec. III.3).

###  II.2 Energy calculation along a structural path

To obtain the minimum energy curve along a structural path and explore
possible metastable structures, we adopt the CINEB method along with DFT$+U$.
The nudged elastic band (NEB) method is an efficient tool for finding the
minimum energy path between two stable structures, i.e. a given initial
(reactant) and final (product) state Mills and Jónsson (1994); Mills  _et al._
(1995). The CINEB method is a small modification of the NEB method without
adding any significant computational method Uberuaga and Jonsson (2000). The
CINEB method yields a rigorous convergence of a saddle point, which has a
maximum energy along the band but minimum in the other directions. The other
images in the band serve for the purpose of defining one degree of freedom for
which the energy of the climbing image is maximized along the band. In this
work, we adopt the CINEB method to explore metastable CO states with distinct
structural distortions following a computed structural path and compute the
energy barrier along the path. We obtain the structural path by defining two
stable CO structures relaxed with different initial conditions and
constructing an energy path between two structures using the CINEB method.

###  II.3 Order Parameter

While ferroelectricity is a phenomenon driven by the spontaneous polarization
of materials, the polarization calculation in a periodic system requires a
careful treatment of the formula Vanderbilt and King-Smith (1993). At the same
time, the inversion symmetry breaking of a structure is a clear indication of
the spontaneous polarization. While we are not interested in obtaining the
quantitative value of the polarization in this work, we will investigate the
displacements of Fe and O planes in the rhombohedral unit cell along the
$[111]_{c}$ direction (Figure 1), where the Fe plane distortion occurs below
the critical temperature.

The displacements of Fe and O planes were investigated in the following way.
First, we confirm that the CO1 and CO3 structures are centrosymmetric and
define the central plane $C$ as the midway between Fe1 and Fe6 planes. Next,
we generated the other dashed-line planes which are equidistant and correspond
to the Fe and O planes of the undistorted high-temperature structure (see
Figure 2). Then, we can quantify how much Fe and O planes are displaced from
the dashed lines. We define the total displacements per unit cell
($\Delta_{tot}$) for these Fe and O planes (see Table 1). For the CO2
structure, the $\Delta_{tot}$ is finite due to the inversion symmetry
breaking, also implying the emergence of ferroelectricity.

##  III Results and Discussions

###  III.1 Structural Relaxations

The study of the neutron diffraction measurement by Battle et al Battle  _et
al._ (1990) at the room temperature showed that bulk LSFO forms a rhombohedral
structure of the space group $R\bar{3}c$ (Figure 1) with lattice constants of
$a=5.47$ Å and $c=13.35$ Å . The rhombohedral unit cell of LSFO has 30 atoms
including two La, four Sr, six Fe, and eighteen O ions. The c-axis of the
rhomboheral unit cell is equivalent to the $[111]_{c}$ direction of the
pseudocubic one, which is conventionally adopted in literatures. Thus, the
cubic $[111]_{c}$ direction will be also adopted in this paper. Above the CO
critical temperature, all Fe ions are equivalent and have the same Fe-O bond
lengths. As the temperature is lowered below TCO, both SDW and CDW orders
develop along the [111]c direction. While the CDW order spans the periodicity
of three Fe ions, the antiferromagnetic SDW repeats in the unit-cell of six Fe
ions, which is commensurate with the crystal lattice periodicity. As a result,
the space group of the crystal structure is lowered to $P\bar{3}m1$ (trigonal,
No. 164) with the point group symmetry of $D_{3d}$ while the crystal remains
centrosymmetric. Here, we find that three distinct CO phases can be stable in
LSFO with the same commensurate modulations of the SDW and CDW, and the
stabilities of these CO structures are dependent on the electronic correlation
effect (the Hubbard $U$ values).

Figure 1: The crystal structure of La1/3Sr2/3FeO3 in a rhombohedral unit cell. The rhombohedral c-axis is equivalent to the $[111]_{c}$ direction of the cubic unit cell.  Figure 2: The schemetics of Fe magnetic moments and the displacements of Fe/O planes for different CO1, CO2, and CO3 phases along the $[111]_{c}$ direction. The displacements are the changes of atomic positions from their undistorted structures (grey dash lines). The central plane C (grey solid line) is midway between Fe1 and Fe6 planes.  Table 1: The relaxed cell parameters ($a$ and $c$), the space group, the Fe plane displacements, and the total displacement ($\Delta_{tot}$) of LSFO at each CO phase. Both CO1 and CO3 structures (space group: $P\bar{3}m1$) have the inversion symmetry, while CO2 and CO$\bar{2}$ structures (space group: $P3m1$) do not. CO phase | c (Å) | a (Å) |  $\Delta_{Fe1}/\Delta_{Fe3}$/$\Delta_{Fe2}$ [Å] |  $\Delta_{tot}$ [Å]
---|---|---|---|---
CO1 | 13.12 | 5.38 | -0.01/0.01/0.00 | 0.00
CO3 | 13.22 | 5.36 | 0.03/-0.03/0.00 | 0.00
CO2 | 13.19 | 5.37 | 0.01/-0.04/-0.03 | -0.22
CO$\bar{2}$ | 13.19 | 5.37 | 0.04/-0.01/0.03 | 0.22

As already noted from Ref. Zhu  _et al._ , 2018, two distinct centrosymmetric
CO structures (CO1 and CO3, see Fig. 2) can be obtained by relaxing them with
different Hubbard $U$ values in DFT$+U$. To explore other metastable CO phases
and the energy barriers, these two CO1 and CO3 structures will be used for two
reference structures in the CINEB method. The first stable structure CO1
(charge ordering 1) was obtained in a strongly correlated regime with $U$=5 eV
and $J$=1 eV, which have been used for LSFO in literatures Saha-Dasgupta  _et
al._ (2005); Zhu  _et al._ (2018); Park  _et al._ (2019). Another stable
structure CO3 was obtained in a weakly correlated regime Saha-Dasgupta  _et
al._ (2005) with $U$=3 eV and $J$=0.6 eV, which was also used by Zhu et al Zhu
_et al._ (2018). Both CO1 and CO3 structures exhibit a sixfold (six Fe ions)
spin density wave (SDW) along the cubic [111]c direction, such that
Fe1$(\uparrow)$Fe2$(\uparrow)$Fe3$(\uparrow)$Fe4$(\downarrow)$Fe5$(\downarrow)$Fe6$(\downarrow)$
(Figure 2).

We find that the $<$Fe-O$>$ mean bond-lengths are closely related to the
magnetic moments. In the CO1 structure, the magnetic moments of Fe1(Fe4) and
Fe3(Fe6) are larger than the one of Fe2(Fe5), so the charge states of Fe1 and
Fe3 should be larger than the one of Fe2 (see Section III.4). As a result, the
high-spin (HS) Fe-O bond expands and the bond-disproportionation occurs.
Particularly, in the CO1 structure the $<$Fe-O$>$ mean bond-lengths of
Fe1-Fe2-Fe3 (similar to Fe4-Fe5-Fe6) show the big-small-big pattern, coupled
to the big-small-big magnetic moments of Fe1, Fe2, and Fe3 ions respectively
(Fig. 2). The $<$Fe-O$>$ mean bond-lengths of CO1 are 1.92 Å for Fe1 and 1.86
Å for Fe2, respectively.

In the case of CO3, the magnetic moments of Fe1 (Fe4) and Fe3 (Fe6) ions
become smaller than Fe2 (Fe5) ion and the bond-length changes to small-big-
small. The bond-lengths are 1.88 Å for Fe1 and 1.94 Å for Fe2. As a result of
the Fe-O bond-length disproportionation, the displacements of Fe and O planes
are also non-uniform as shown in Fig. 2.

In Table 1, we list the relaxed unit-cell parameters, the space group, the
displacements of Fe planes, and the total displacement of Fe and O planes
($\Delta_{tot}$) along the [111]c. Both of CO1 and CO3 have the space group of
$P\bar{3}m1$, implying CO1 and CO3 are centrosymmetric. Also, based on the
displacements of Fe and O planes of CO1 and CO3, the Figure 2 shows that the
reflection of the supercells, including Fe1O6-Fe6O6 cells, of CO1 and CO3
around the central plane C yields the same supercells respectively. Finally,
the total displacements of Fe and O planes in CO1 and CO3 are also zeros
(Table 1), meaning polarization is not induced in CO1 and CO3.

###  III.2 Energy Path along Multiple Charge Orderings

Figure 3: (a) The relative energies per formula unit of CO1, CO2, and CO3
phases as a function of the Hubbard $U$. (b) Comparison of CINEB and Linear
Interpolation energies vs Image structures calculated with DFT+U at $U$=3.62
eV and $J$=0.724 eV.

In this section, we compute the energy path between two energetically
degenerate CO phases (CO1 and CO3) to explore the possible metal-stable CO
states along the structural path. We first tuned $U$ ($J$=0.2$U$) for CO1 and
CO3 phases while relaxing crystal structures at a fixed volume to investigate
their stability and plot the relative energy
$\Delta$E${}_{CO1-CO3}(=E[\mathrm{CO1}]-E[\mathrm{CO3}])$ between them as a
function of $U$ in Fig. 3(a). Here, we find that the low-temperature
experimental ground-state structure (CO1) is stable when the $U$ value becomes
larger than 3.7 eV (see Fig. 5(a)). While DFT$+U$ is a zero-temperature
theory, we find that the energetics of different CO phases can be tuned by
changing the onsite Coulomb repulsion $U$, which can mimic the effect of
temperature, pressure, or photoinduced excitation. In principle, laser
excitation in experiment can modulate the electronic structure from the
ground-state, by affecting the exchange interaction Mentink (2017), and may
eventually trigger a phenomenon called “photoinduced structural phase
transition” Nasu (2004).

Fig. 3(a) shows that the energy difference between CO1 and CO3 structures can
be almost zero at $U_{c}=3.62$ eV ($J=0.724$ eV). This means that other meta-
stable CO structures could be found in the CINEB calculation near $U=3.62$eV.
At $U$=3.62 eV, CO1(CO3) still has the big-small-big (small-big-small) bond-
order (see Fig. 5). Here, we perform a CINEB calculation at $U=3.62$ eV using
both CO1 and CO3 as two reference structures. Remarkably, Fig. 3(b) shows that
the CINEB energy curve calculated with $U_{c}$=3.62 eV can capture a meta-
stable structure of CO2, whose energy is only 3meV above the CO1 or CO3
structure with the energy barrier of $\sim$7meV. This CO2 structure is
obtained by the spontaneous displacement of the Fe plane and it can not be
captured by the linear interpolation method where the image structures along
the path are obtained by linearly interpolating atomic positions between CO1
and CO3.

The obtained CO2 structure has the small-medium-big FeO6 octahedra (Fe-O bond
order), (see Fig. 5(b)), coupled to the magnetic ordering of 2.2 (small), 3.0
(medium), and 3.4 $\mu_{B}$ (big). Unlike CO1 and CO3, CO2 has the space group
of $P3m1$ (trigonal, No. 156) with the point group symmetry of $C_{3v}$ ,
which belongs to a polar point group Shi  _et al._ (2016). The reflection of
CO2 supercell around the C plane does not yield the same supercell (see Fig.
2), implying the broken inversion symmetry in CO2. Also, the total
displacement ($\Delta_{tot}$) of Fe and O planes in CO2 is not zero (see Table
1), resulting the spontaneous polarization. Remarkably, the CO2 phase is
metallic, which is not common for the ferroelectric materials due to the
screening of the polarizationAnderson and Blount (1965); Zhou and Ariando
(2020). We also find that the other meta-stable CO$\bar{2}$ structure can be
obtained by applying the inversion operation to the CO2 structure about the
central plane in Fig. 2 and shows the opposite polarization compared to the
CO2 case (see Table 1). The polarizations between CO2 and CO$\bar{2}$ are
easily switchable due to the very small energy barrier of $\sim$7meV between
them. The existence of ferroelectricity in LSFO was first pointed out by Park
et al. Park  _et al._ (2019), but its structure was obtained manually by
exchanging the La and Sr layers.

Figure 4:  Comparison between the CINEB and the Linear Interpolation paths:
(a) The displacement of Fe2 plane and the total displacement of Fe and O
planes. (b) The bond angle $\angle$O1-Fe2-O2.

To address the difference between the CINEB and Linear Interpolation results,
we compare the displacement ($\Delta_{Fe2}$) of Fe2 plane, the total
displacement ($\Delta_{tot}$) of Fe and O planes, and the bond angle along
O1-Fe2-O2 ($\angle$O1-Fe2-O2). Figure 4(a) shows that along the Linear
Interpolation path the displacement of Fe2 plane $\Delta_{Fe2}^{LI}$ and the
total displacement of Fe and O planes $\Delta_{tot}^{LI}$ are kept zero, while
along the CINEB path an abrupt change of the displacement of Fe2 plane
$\Delta_{Fe2}^{CINEB}$ and the total displacement $\Delta_{tot}^{CINEB}$
occurs at image number 5 and reach a minimum at image 9 where CO2 was
captured. This change of the Fe2 displacement along the CINEB path is also
accompanied by a sudden change of the bond angle ($\angle$O1-Fe2-O2)CINEB,
while the one along the Linear Interpolation path ($\angle$O1-Fe2-O2)LI
remains $180^{0}$ (Figure 4(b)). The existences of this CO2 phase might be
captured by Sabyasachi et al. and Yang et al., where the neutron diffraction
shows the multiple $Q$ plane magnetic reflections with equivalent intensities
Sabyasachi  _et al._ (2012); Yang  _et al._ (2003).

###  III.3 The dependence of structural parameters on $U$

Figure 5: The $<$Fe-O$>$ mean bond lengths vs the Hubbard $U$ ($J=0.2U$)
obtained for (a) CO1, (b) CO2, and (c) CO3 phases. The vertical dash lines
represent the $U$ values used for stabilizing each CO phase, namely CO1
($U=5.0$ eV), CO2 ($U=3.62$ eV), and CO3 ($U=3.0$ eV) as shown in Table 1. The
white (shaded) regions represent metallic (insulating) phases.

Our structural relaxation results show that the stability of the different CO
phases depends sensitively on the Hubbard $U$ values ($J=0.2U$). In general,
the correlation effect of the Hubbard $U$ is important to stabilize the
bond/charge disproportionation in many oxides including nickelates Park  _et
al._ (2012); Liao  _et al._ (2021), cobaltates Park  _et al._ (2020), ferrites
Greenberg  _et al._ (2018), and manganites Valli  _et al._ (2015). This is
because only parts of the transition metal $M$ sites undergo the spin-state
transition to the HS state with the $M-O$ bond elongation and more HS sites
are populated with the stronger $U$ values. Our calculation confirms that the
DFT-relaxed structure of LSFO shows no Fe-O bond disproportionation,
consistently as the experimental high-temperature structure, and the increase
of $U$ energetically favors the structures with more HS states in a non-
trivial way.

Fig. 5(a) shows that the CO1 structure as shown in Table 1 can be stable only
when $U>3.7$ eV ($J$=0.74 eV) and the structural transition to CO3 occurs
along with the insulator-metal transition. The CO2 structure as shown in Table
1 is meta-stable in a narrow $U$ range of $3.55\leq U\leq 3.7$ eV and evolves
into a distorted CO3 phase as $U$ becomes lower than 3.55 eV. The CO3
structure can be stable in a wide-range of $U$ values although this phase is
energetically lower than CO1 or CO2 phases when $U\leq 3.62$ eV. Both of CO2
and CO3 structures converge to the high-temperature structure without the Fe-O
disproportionation as $U$ becomes smaller than 2eV. We find that the
insulating phase in LSFO occurs only in the CO1 structure with $U>$ 3.7eV.

###  III.4 Electronic Structure and magnetism in LSFO

Here, we investigate electronic structures of LSFO at different CO states
computed using DFT+U. Due to the AFM structure, the Fe1/Fe2/Fe3 density of
states (DOS) is equivalent to the Fe4/Fe5/Fe6 one in LSFO once their spins are
flipped. For CO1 and CO3, the crystal structures are centrosymmetric and we
show only Fe1 and Fe2 DOS since Fe1 (Fe4) and Fe3 (Fe6) are equivalent. To
distinguish the importance of electronic correlations from the structure
effect, we compare $U=3.62$ eV ($J=0.724$ eV) and $U=4$ eV ($J=0.8$ eV) DOS at
the fixed structure of each CO phase.

At $U$=4 eV, the CO1 phase is an insulating state with the spectral gap size
of $\sim$120 meV (see Fig. 6(a)), consistent with the optical gap measurement
in LSFO at a low temperature Ishikawa  _et al._ (1998). In the Fe1 ion, both
eg and t2g bands are half-filled with the gap size comparable to $U$ behaving
as a typical Mott insulator. However, only the t2g bands of Fe2 are half-
filled, while the eg bands are almost empty (see Fig. 6(a)). This is
consistent with the high-spin picture of the charge-ordering state between Fe1
($d^{5}$; t${}_{2g\uparrow}^{3}$e${}_{g\uparrow}^{2}$) and Fe2 ($d^{3}$;
t${}_{2g\uparrow}^{3}$e${}_{g\uparrow}^{0}$) ions. As the correlation becomes
weaker ($U$=3.62 eV), the DOS for CO1 becomes metallic as the Fe1 eg (Fe2 t2g)
state is less (more) occupied and the spectral gap at the Fermi energy is
closed.

In CO3 at $U$=4eV, the charge-ordering pattern changes for Fe1 ($d^{4}$;
t${}_{2g\uparrow}^{3}t_{2g\downarrow}^{1}e_{g\uparrow}^{0}$) and Fe2 ($d^{5}$;
t${}_{2g\uparrow}^{3}t_{2g\downarrow}^{1}e_{g\uparrow}^{1}$). The spin state
for Fe1 changes to the low-spin, while the Fe2 spin is close to the
intermediate one. This is because the crystal field splitting of the Fe1 ion
becomes larger due to the smaller octahedron size compared to Fe1 in the CO1
phase. As a result, both Fe1 t2g and Fe2 t2g states are partially filled and
the DOS becomes metallic (see Fig. 6(a)). As the correlation becomes weak
($U$=3.62 eV), the CO3 phase remains metallic.

Similar to CO3, CO2 is metallic at both $U$=4eV and 3.62eV. As the Fe1 $d$ DOS
of CO2 is similar to the Fe1 $d$ DOS of CO3 and the Fe1-O bond-lengths of CO2
and CO3 are similar to each other as well, we expect that the local electronic
configuration of Fe1 should be similarly given as the low-spin $d^{4}$
(t${}_{2g\uparrow}^{3}t_{2g\downarrow}^{1}e_{g\uparrow}^{0}$). Moreover, the
Fe-O bond-length of the Fe1 ion is the smallest, while those of Fe2 and Fe3
ions are close to each other implying the similar electronic structure between
Fe2 and Fe3. Nevertheless, the evidence of the CO can be found near
$E\approx-1eV$ where the occupied Fe3 eg states have slightly more DOS than
the Fe2 one, while their t2g DOS are similar. This implies that the local
electronic configurations of Fe2 and Fe3 ions should be Fe(3.5+δ)+
(t${}_{2g\uparrow}^{3}t_{2g\downarrow}^{1}e_{g\uparrow}^{0.5-\delta}$) and
Fe(3.5-δ)+
(t${}_{2g\uparrow}^{3}t_{2g\downarrow}^{1}e_{g\uparrow}^{0.5+\delta}$),
respectively.

Figure 6: The DOS plots of CO1, CO2, and CO3 phases calculated with DFT$+U$.
(a) $U=4.0$ eV and $J=0.8$ eV and (b) $U=3.62$ eV and $J=0.724$ eV. Schematic
energy diagrams of Fe t2g and eg orbitals are also shown in the insets.

The calculated magnetic moments of Fe ions ($m_{Fe1}$, $m_{Fe2}$, and
$m_{Fe3}$) are coupled to the above valence states and these values for CO1,
CO2, and CO3 are shown in Table 2. The magnetic moments in CO1 calculated with
DFT+U ($U$=4 eV, $J$=0.8 eV) are in a good agreement with the experimental
ones recently obtained by Li et al Li  _et al._ (2018) at low temperature. The
calculated value of $m_{Fe1}$ is rather screened from the electronic
configuration estimation based on the DOS since we expect $m_{Fe1}$
(t${}_{2g\uparrow}^{3}$e${}_{g\uparrow}^{2}$) = 5$\mu_{B}$, while the
$m_{Fe2}$ value is consistent (t${}_{2g\uparrow}^{3}$e${}_{g\uparrow}^{0}$ =
3$\mu_{B}$). The expected moments of the Fe1 and Fe2 ions in CO3 are $m_{Fe1}$
(t${}_{2g\uparrow}^{3}$t${}_{2g\downarrow}^{1}$e${}_{g\uparrow}^{0}$) =
2$\mu_{B}$ and $m_{Fe2}$
(t${}_{2g\uparrow}^{3}$t${}_{2g\downarrow}^{1}$e${}_{g\uparrow}^{1}$) =
3$\mu_{B}$, respectively. However, since CO3 is metallic and the magnetic
moments calculated with DFT+$U$ also depends on $U$ and $J$, our calculated
moments of 2.42$\mu_{B}$ and 3.52$\mu_{B}$ at $U$=4 eV are larger than these
expected values. We confirmed that the magnetic moments of Fe1 and Fe2 are
reduced at $U$=3 eV as 2.08 and 3.14 $\mu_{B}$, similar to the expected values
respectively.

Similarly, the magnetic moment of Fe1 at CO2 calculated with $U$=4 eV is 2.70
$\mu_{B}$, which is still large for a LS state of Fe4+ (2.0 $\mu_{B}$). We
find that this moment computed using $U$=3.62eV is more consistent as
2.18$\mu_{B}$. For CO2, the magnetic moments of $m_{Fe1}$, $m_{Fe2}$, and
$m_{Fe3}$ show the small-medium-big pattern, which is consistent with the
charge-ordering pattern of Fe4+-Fe(3.5+δ)+-Fe(3.5-δ)+.

Table 2: Magnetic moments of Fe1, Fe2, and Fe3 ions, namely $m_{Fe1}$,
$m_{Fe2}$, and $m_{Fe3}$ calculated with DFT$+U$ and compared with the
experimental ones obtained by Li et alLi  _et al._ (2018)

LSFO | $m_{Fe1}[\mu_{B}]$ | $m_{Fe2}[\mu_{B}]$ | $m_{Fe3}[\mu_{B}]$
---|---|---|---
CO1 (U=4.0 eV) | 3.74 (Fe3+) | 3.12 (Fe5+) | 3.74 (Fe3+)
Experiment.Li  _et al._ (2018) | 3.67 (Fe3+) | 3.26 (Fe5+) | 3.67 (Fe3+)
CO3 (U=4.0 eV) | 2.42 (Fe4+) | 3.52 (Fe3+) | 2.42 (Fe4+)
CO2 (U=4.0 eV) | 2.70 (Fe4+) | 3.26 (Fe(3.5+δ)+) | 3.60 (Fe(3.5-δ)+)
CO2 (U=3.62 eV) | 2.18 (Fe4+) | 2.99 (Fe(3.5+δ)+) | 3.36 (Fe(3.5-δ)+)

##  IV Conclusion

In conclusion, we studied the structural and electronic properties of charge-
ordered La doped SrFeO3, La1/3Sr2/3FeO3 (LSFO) systematically using DFT+U
along with the antiferromagnetic order. We find that metastable structures
with distinct CO phases in LSFO can be obtained by relaxing the structures
with the different $U$ values varying the correlation effect. The DFT+U
calculation of LSFO with $U$=5eV can capture the low temperature CO phase (CO1
in the main text) of the big-small-big pattern, where the enhanced charge
density is accompanied by the large magnetic moment with the Fe-O bond
elongation. The ground-state is insulating as the spectral function at the
Fermi energy opens a Mott gap driven by the high-spin states of Fe ions.

As the correlation effect becomes weak by reducing the $U$ value in DFT+U, we
can capture other metastable CO phases with distinct Fe-O bond patterns. One
CO phase (CO3 in the main text) shows the crystal structure with the same
space group as the CO1 phase, while the CO pattern changes to small-big-small.
The other metastable CO phase (CO2 in the main text) can be obtained by
interpolating the structural path between CO1 and CO3 phases using the CINEB
calculation. Remarkably, the CO2 phase stabilize a lower symmetry crystal
structure along with the inversion symmetry breaking and it shows the
ferroelectric metallic state driven by the big-medium-small CO pattern. This
CO2 phase can not be captured by the linear interpolation method as it
requires the spontaneous displacement of Fe ions at the symmetric points. The
electronic structures of these metastable CO states are notably changed as
both CO2 and CO3 phases are metallic while the ground-state CO1 phase is
insulating. The energy barrier of this CO2 phase along the structural path is
only $\sim$7meV, which can be easily accessible from experiments by applying
the pressure or the optical excitation.

Our results suggest that the strong correlation effect plays an important role
to study and stabilize the multiple CO phases of transition metal oxides
accompanying the mixed valence and the metal-oxygen bond disproportionation.
The CINEB method combined with the energy and force calculations based on the
first-principles can capture such metastable CO phases along with the distinct
electronic structure from their ground state. While DFT+U is an efficient
static method to incorporate the correlation effect, it can generally suffer
from the convergence problem in systems with multiple correlated states Park
_et al._ (2020); Ratcliff  _et al._ (2022). More advanced first-principle
method such as dynamical mean field theory (DMFT) can be a promising way to
study metastable phases in strongly correlated materials driven by both the
structural distortion and the strong correlations especially when the CINEB
method is combined with the energy and force calculations within DMFT Leonov
_et al._ (2014); Haule and Pascut (2016).

## Acknowledgement

We thank Yue Cao for fruitful discussions. NN, AL, AN, and HP acknowledge
financial support from the US Department of Energy, Office of Science, Office
of Basic Energy Sciences, Materials Science and Engineering Division. VS was
supported from NSF SI2-SSE Grant 1740112. We gratefully acknowledge the
computing resources provided on Bebop, a high- performance computing cluster
operated by the Laboratory Computing Resource Center at Argonne National
Laboratory.

## References

  * van den Brink and Khomskii (2008) J. van den Brink and D. I. Khomskii, Journal of Physics: Condensed Matter 20, 434217 (2008).
  * Efremov  _et al._ (2005) D. Efremov, J. Brink, and D. Khomskii, Nature materials 3, 853 (2005).
  * Wan  _et al._ (2016) X. Wan, Q. He, S. Savrasov, and C.-G. Duan, Scientific Reports 6, 22743 (2016).
  * Park  _et al._ (1999) S. K. Park, T. Ishikawa, Y. Tokura, J. Q. Li, and Y. Matsui, Phys. Rev. B 60, 10788 (1999).
  * Battle  _et al._ (1990) P. Battle, T. Gibb, and P. Lightfoot, Journal of Solid State Chemistry 84, 271 (1990).
  * Sabyasachi  _et al._ (2012) S. Sabyasachi, M. Patra, S. Majumdar, S. Giri, S. Das, V. S. Amaral, O. Iglesias, W. Borghols, and T. Chatterji, Phys. Rev. B 86, 104416 (2012).
  * Dann  _et al._ (1994) S. Dann, D. Currie, M. Weller, M. Thomas, and A. Al-Rawwas, Journal of Solid State Chemistry 109, 134 (1994).
  * Yang  _et al._ (2003) J. Yang, X. Zhou, Z. Chu, W. Hikal, Q. Cai, J.-C. Ho, D. Kundaliya, W. Yelon, W. James, H. Anderson, H. Hamdeh, and S. Malik, Journal of Physics: Condensed Matter 15, 5093 (2003).
  * Takano  _et al._ (1981) M. Takano, J. Kawachi, N. Nakanishi, and Y. Takeda, Journal of Solid State Chemistry 39, 75 (1981).
  * Okamoto  _et al._ (2010) J. Okamoto, D. J. Huang, K. S. Chao, S. W. Huang, C.-H. Hsu, A. Fujimori, A. Masuno, T. Terashima, M. Takano, and C. T. Chen, Phys. Rev. B 82, 132402 (2010).
  * Li  _et al._ (2018) F. Li, V. Pomjakushin, T. Mazet, R. Sibille, B. Malaman, R. Yadav, L. Keller, M. Medarde, K. Conder, and E. Pomjakushina, Phys. Rev. B 97, 174417 (2018).
  * Li  _et al._ (1997) J. Q. Li, Y. Matsui, S. K. Park, and Y. Tokura, Phys. Rev. Lett. 79, 297 (1997).
  * Ishikawa  _et al._ (1998) T. Ishikawa, S. K. Park, T. Katsufuji, T. Arima, and Y. Tokura, Phys. Rev. B 58, R13326 (1998).
  * Chernyshev  _et al._ (2014) A. Chernyshev, S. Petrov, and N. Uvarov, Physics of the Solid State 56, 2519–2523 (2014).
  * Shimony and Knudsen (1966) U. Shimony and J. M. Knudsen, Phys. Rev. 144, 361 (1966).
  * Zhu  _et al._ (2018) Y. Zhu, J. Hoffman, C. E. Rowland, H. Park, D. A. Walko, J. W. Freeland, P. J. Ryan, R. D. Schaller, A. Bhattacharya, and H. Wen, Nature Communications 9, 1799 (2018).
  * Yamamoto  _et al._ (2022) K. Yamamoto, T. Tsuyama, S. Ito, K. Takubo, I. Matsuda, N. Pontius, C. Schüßler-Langeheine, M. Minohara, H. Kumigashira, Y. Yamasaki, H. Nakao, Y. Murakami, T. Katase, T. Kamiya, and H. Wadati, New Journal of Physics 24, 043012 (2022).
  * Matsuno  _et al._ (2002) J. Matsuno, T. Mizokawa, A. Fujimori, Y. Takeda, S. Kawasaki, and M. Takano, Phys. Rev. B 66, 193103 (2002).
  * Saha-Dasgupta  _et al._ (2005) T. Saha-Dasgupta, Z. Popovic, and s. Satpathy, Phys. Rev. B 72 (2005), 10.1103/PhysRevB.72.045143.
  * Park  _et al._ (2019) S. Y. Park, K. M. Rabe, and J. B. Neaton, Proceedings of the National Academy of Sciences 116, 23972 (2019).
  * Liechtenstein  _et al._ (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Feynman (1939) R. P. Feynman, Phys. Rev. 56, 340 (1939).
  * Mills and Jónsson (1994) G. Mills and H. Jónsson, Phys. Rev. Lett. 72, 1124 (1994).
  * Mills  _et al._ (1995) G. Mills, H. Jónsson, and G. K. Schenter, Surface Science 324, 305 (1995).
  * Uberuaga and Jonsson (2000) B. Uberuaga and H. Jonsson, J. Chem. Phys. 113, 9901 (2000).
  * Vanderbilt and King-Smith (1993) D. Vanderbilt and R. D. King-Smith, Phys. Rev. B 48, 4442 (1993).
  * Mentink (2017) J. H. Mentink, Journal of Physics: Condensed Matter 29, 453001 (2017).
  * Nasu (2004) K. Nasu,  _Photoinduced phase transitions_ (World Scientific, 2004).
  * Shi  _et al._ (2016) P.-P. Shi, Y.-Y. Tang, P.-F. Li, W.-Q. Liao, Z.-X. Wang, Q. Ye, and R.-G. Xiong, Chem. Soc. Rev. 45, 3811 (2016).
  * Anderson and Blount (1965) P. W. Anderson and E. I. Blount, Phys. Rev. Lett. 14, 217 (1965).
  * Zhou and Ariando (2020) W. X. Zhou and A. Ariando, Japanese Journal of Applied Physics 59, SI0802 (2020).
  * Park  _et al._ (2012) H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev. Lett. 109, 156402 (2012).
  * Liao  _et al._ (2021) X. Liao, V. Singh, and H. Park, Phys. Rev. B 103, 085110 (2021).
  * Park  _et al._ (2020) H. Park, R. Nanguneri, and A. T. Ngo, Phys. Rev. B 101, 195125 (2020).
  * Greenberg  _et al._ (2018) E. Greenberg, I. Leonov, S. Layek, Z. Konopkova, M. P. Pasternak, L. Dubrovinsky, R. Jeanloz, I. A. Abrikosov, and G. K. Rozenberg, Phys. Rev. X 8, 031059 (2018).
  * Valli  _et al._ (2015) A. Valli, H. Das, G. Sangiovanni, T. Saha-Dasgupta, and K. Held, Phys. Rev. B 92, 115143 (2015).
  * Ratcliff  _et al._ (2022) L. E. Ratcliff, L. Genovese, H. Park, P. B. Littlewood, and A. Lopez-Bezanilla, Journal of Physics: Condensed Matter 34, 094003 (2022).
  * Leonov  _et al._ (2014) I. Leonov, V. I. Anisimov, and D. Vollhardt, Phys. Rev. Lett. 112, 146401 (2014).
  * Haule and Pascut (2016) K. Haule and G. L. Pascut, Phys. Rev. B 94, 195146 (2016).