# The PEPSI-LBT Exoplanet Transit Survey (PETS).111Based on data acquired with
the Potsdam Echelle Polarimetric and Spectroscopic Instrument (PEPSI) using
the Large Binocular Telescope (LBT) in Arizona.
II. A Deep Search for Thermal Inversion Agents in KELT-20 b/MASCARA-2 b with
Emission and Transmission Spectroscopy

Marshall C. Johnson Department of Astronomy, The Ohio State University, 4055
McPherson Laboratory, 140 West 18th Ave., Columbus, OH 43210 USA  Ji Wang (王吉)
Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory,
140 West 18th Ave., Columbus, OH 43210 USA  Anusha Pai Asnodkar Department of
Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 West 18th
Ave., Columbus, OH 43210 USA  Aldo S. Bonomo  INAF – Osservatorio Astrofisico
di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, Italy  B. Scott
Gaudi Department of Astronomy, The Ohio State University, 4055 McPherson
Laboratory, 140 West 18th Ave., Columbus, OH 43210 USA  Thomas Henning  Max-
Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
Ilya Ilyin  Leibniz-Institute for Astrophysics Potsdam (AIP), An der
Sternwarte 16, D-14482 Potsdam, Germany  Engin Keles  Leibniz-Institute for
Astrophysics Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany
Luca Malavolta  INAF – Osservatorio Astronomico di Padova, Vicolo
dell’Osservatorio 5, I-35122 Padova, Italy  Dipartimento di Fisica e
Astronomia “Galileo Galilei”, Universita degli Studi di Padova, I-35122
Padova, Italy  Matthias Mallonn Leibniz-Institute for Astrophysics Potsdam
(AIP), An der Sternwarte 16, D-14482 Potsdam, Germany  Karan Molaverdikhani
Universitäts-Sternwarte, Ludwig-Maximilians-Universität München,
Scheinerstrasse 1, D-81679 München, Germany  Exzellenzcluster Origins,
Boltzmannstraße 2, 85748 Garching, Germany  Max-Planck-Institut für
Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany  Valerio Nascimbeni
INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122
Padova, Italy  Jennifer Patience  School of Earth and Space Exploration,
Arizona State University, 660 S. Mill Ave., Tempe, Arizona 85281, USA  Katja
Poppenhaeger  Leibniz-Institute for Astrophysics Potsdam (AIP), An der
Sternwarte 16, D-14482 Potsdam, Germany  Institute of Physics & Astronomy,
University of Potsdam, Karl-Liebknecht-Str. 24/25, D-14476 Potsdam, Germany
Gaetano Scandariato INAF – Osservatorio Astrofisico di Catania, via S. Sofia
78, I-95123 Catania, Italy  Everett Schlawin  Steward Observatory, University
of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA  Evgenya Shkolnik
School of Earth and Space Exploration, Arizona State University, 660 S. Mill
Ave., Tempe, Arizona 85281, USA  Daniela Sicilia  INAF – Osservatorio
Astrofisico di Catania, via S. Sofia 78, I-95123 Catania, Italy  Alessandro
Sozzetti  INAF – Osservatorio Astrofisico di Torino, Strada Osservatorio 20,
I-10025 Pino Torinese, Italy  Klaus G. Strassmeier  Leibniz-Institute for
Astrophysics Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany
Institute of Physics & Astronomy, University of Potsdam, Karl-Liebknecht-Str.
24/25, D-14476 Potsdam, Germany  Christian Veillet  Large Binocular Telescope
Observatory, 933 N. Cherry Ave., Tucson, AZ 85721, USA  Fei Yan Institute of
Astrophysics, University of Göttingen, Friedrich-Hund-Platz 1 , D-37077
Göttingen, Germany  Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117
Heidelberg, Germany

###### Abstract

Recent observations have shown that the atmospheres of ultra hot Jupiters
(UHJs) commonly possess temperature inversions, where the temperature
increases with increasing altitude. Nonetheless, which opacity sources are
responsible for the presence of these inversions remains largely
observationally unconstrained. We used LBT/PEPSI to observe the atmosphere of
the UHJ KELT-20 b in both transmission and emission in order to search for
molecular agents which could be responsible for the temperature inversion. We
validate our methodology by confirming a previous detection of Fe i in
emission at $15.1\sigma$; however, we are unable to reproduce published
detections of Fe ii, Cr i, or Si i. We attribute the non-detection of Si i to
the lack of lines in our bandpass, but the non-detections of Fe ii and Cr i
are puzzling due to our much higher signal-to-noise ratio than previous works.
Our search for the inversion agents TiO, VO, FeH, and CaH results in non-
detections. Using injection-recovery testing we set $4\sigma$ upper limits
upon the volume mixing ratios for these constituents as low as $\sim 1\times
10^{-10}$ for TiO. For TiO, VO, and CaH, our limits are much lower than
expectations from an equilibrium chemical model, while FeH is lower than the
expectations only from a super-Solar metallicity model. We thus rule out TiO,
VO, and CaH as the source of the temperature inversion in KELT-20 b, while FeH
is disfavored only if KELT-20 b possesses a high-metallicity atmosphere.

††software: astropy (Astropy Collaboration et al., 2013, 2018), jupyter
(Kluyver et al., 2016), matplotlib (Hunter, 2007), molecfit (Smette et al.,
2015; Kausch et al., 2015), numpy (van der Walt et al., 2011), petitRADTRANS
(Mollière et al., 2019), PyFastChem (Stock et al., 2018), uncertainties
(Lebigot, 2017)

##  1 Introduction

Hot Jupiters (HJs) offer some of our best chances at present to characterize
exoplanetary atmospheres in detail. These planets’ large radii, hot
temperatures, short orbital periods, and often bright host stars mean that we
can obtain more frequent observations and overall higher signal-to-noise
ratios than are possible for smaller, cooler, or longer-period planets.

Several ultra hot Jupiters (UHJs; planets with equilibrium temperatures in
excess of 2200 K) have recently been shown to have atmospheric temperature
inversions through the presence of emission lines of atomic species in their
spectra; if there was no inversion, the spectrum would present absorption
lines instead. A temperature inversion occurs when some absorber exists within
the atmosphere, heating the atmosphere and causing the temperature to
increase, rather than decrease, with increasing altitude. UHJs had previously
been proposed to harbor thermal inversions, with the transition temperature in
this regime occurring anywhere from 1700 to 2500 K (e.g., Baxter et al., 2020;
Lothringer et al., 2018). Individual planets with detections of line emission
due to an inversion include WASP-33 b (Nugroho et al., 2017), WASP-121 b
(Evans et al., 2017), WASP-18 b (Sheppard et al., 2017), WASP-189 b (Yan et
al., 2020, 2022a), KELT-9 b (Kasper et al., 2021), and KELT-20 b (Cont et al.,
2021a). The number of possible inversion agents in these planets is limited,
as only a few molecular species can survive at these high temperatures.

Nonetheless, direct evidence for which species cause these inversions has been
lacking. Nugroho et al. (2017) and Cont et al. (2021b) detected TiO in WASP-33
b, but Herman et al. (2020) and Serindag et al. (2021) were unable to
confidently detect TiO with their own data or methods. Prinoth et al. (2022)
detected TiO in absorption in WASP-189 b. VO has been detected at low
resolution in WASP-121 b, but high-resolution observations have been unable to
confirm this detection, likely due to inaccuracies with the VO line list
(Hoeijmakers et al., 2020a; Merritt et al., 2021; de Regt et al., 2022). No
molecular inversion agent candidate has been detected in either KELT-9 b or
KELT-20 bdespite the presence of atomic emission lines, motivating further
searches; KELT-9 b, however, may be too hot for any molecular agent to exist.

A number of molecular species, including TiO, VO (Fortney et al., 2008), SiO,
metal hydrides (Lothringer et al., 2018), AlO, CaO, NaH, and MgH Gandhi &
Madhusudhan (2019), have been proposed to cause inversions. Alternately, in
the hottest HJs, atomic species such as Fe or Mg, or continuum H- opacity,
could give rise to inversions (Lothringer et al., 2018).

KELT-20 b/MASCARA-2 b (hereafter KELT-20 b for brevity) is an UHJ discovered
independently by both the Kilodegree Extremely Little Telescope (KELT) and
Multi-site All-Sky CAmeRA (MASCARA) ground-based transit surveys (Lund et al.,
2017; Talens et al., 2018). With a $V$-band magnitude of 7.58, the host star
HD 185603 is one of the five brightest stars to host a known transiting giant
planet. It has thus been a target for atmospheric observations using all
available methods (e.g., Wong et al., 2021; Casasayas-Barris et al., 2018;
Cont et al., 2021a).

Many observations of the atmosphere of KELT-20 b have already been obtained
via high-resolution transmission spectroscopy. These have allowed the
detection of several atomic species in the planet’s atmosphere: H i, Na i,
(Casasayas-Barris et al., 2018), Ca ii, Fe ii (Casasayas-Barris et al., 2019;
Nugroho et al., 2020; Bello-Arufe et al., 2022), Fe i (Stangret et al., 2020),
Mg i, and Cr ii (Hoeijmakers et al., 2020b). Critically, however, Nugroho et
al. (2020) reported non-detections of the potential thermal inversion agents
NaH, MgH, AlO, SH, CaO, VO, FeH, and TiO, while Kesseli et al. (2020) found a
low-confidence signal for FeH.

An alternative method is to use emission spectroscopy. Although the emission
signals are overall smaller than the transmission signals, emission
spectroscopy directly probes the altitudes around the temperature inversion,
since this is where the emission lines must be formed.

Several species have already been detected in emission for KELT-20 b at high
resolution: Fe i (Yan et al., 2022b; Borsa et al., 2022), Si i (Cont et al.,
2021a), Fe ii, and Cr i (Borsa et al., 2022). Fu et al. (2022) also detected
H2O and CO emission at low resolution. These detections of emission lines not
only add to the detected atmospheric species but also demonstrate the presence
of an inverted temperature profile. A monotonically-decreasing pressure-
temperature (P-T) profile would result in absorption rather than emission
lines.

Here, we aim to detect or constrain which species may be responsible for this
inversion. Yan et al. (2022b) already obtained non-detections for TiO, VO, and
FeH in their CARMENES data, while Borsa et al. (2022) did the same for AlO,
TiO, and VO with HARPS-N; however, both of those spectrographs are fed by
4-m-class telescopes, motivating an independent search with 8-m-class
telescopes and including additional potential inversion agents.

The PEPSI-LBT Exoplanet Transit Survey (PETS) is a large ($\sim$400 hour)
survey of transiting exoplanets using the PEPSI high-resolution spectrograph
on the Large Binocular Telescope (LBT). The first paper resulting from the
survey set limits on the atmospheric presence of atomic species that could be
vaporized from the crust of the super-Earth 55 Cnc e (Keles et al., 2022). In
this second paper in the series, we analyze our emission and transmission
spectra of KELT-20 b.

##  2 Observations

We obtained our data on KELT-20 b with the Potsdam Échelle Polarimetric and
Specrographic Instrument (PEPSI; Strassmeier et al., 2015) on the $2\times
8.4$ m Large Binocular Telescope (LBT) located on Mt. Graham, Arizona, USA. We
observed using the 1.5” fiber, giving a resolving power of $R=130,000$. We
used Cross Dispersers (CD) III and V, covering 4800-5441 Å and 6278-7419 Å in
the blue and red arms of the spectrograph, respectively. The 2019 transit
observations were unbinned and used 600 second exposures, while the 2021
emission observations used $2\times 1$ binning on the chip and 300 second
exposures.

We observed one transit of KELT-20 b and two $\sim 4.5$-hour segments near
secondary eclipse. We list the key parameters of our datasets in Table 1, and
show the phases covered by our observations graphically in Fig. 1.

We reduced the data using the Spectroscopic Data Systems (SDS) pipeline, as
described in Ilyin (2000) and Strassmeier et al. (2018). The pipeline performs
bias subtraction and flat field correction, order tracing and extraction,
cosmic ray correction, wavelength calibration, and normalization and
combination of all orders in each arm (red versus blue). The pipeline also
estimates the variance in each pixel, which we pipe through our analysis
procedure described below.

Table 1: Log of Observations

Date (UT) | Type | $N_{\mathrm{spec}}$ | Exp. Time (s) | Airmass Range | Phases Covered | SNRblue | SNRred
---|---|---|---|---|---|---|---
2019 May 4 | Transmission | 23 | 600 | 1.00-2.01 | -0.023-0.034 | 288 | 308
2021 May 1 | Emission | 47 | 300 | 1.01-2.03 | 0.529-0.582 | 301 | 340
2021 May 18 | Emission | 45 | 300 | 1.00-1.48 | 0.422-0.473 | 347 | 397

Note. — $N_{\mathrm{spec}}$ is the number of spectra obtained on that night.
Exp. time is the exposure time in seconds. SNRblue and SNRred are the nightly
average of the $95^{\mathrm{th}}$ quantile per-pixel signal-to-noise ratios in
the blue and red arms, respectively.

Figure 1: Graphic illustrating the phases at which our observations of KELT-20
b took place. Different nights are marked by different colors. The planetary
orbit and the stellar surface are the large and small circles, respectively.
The observer is viewing from the bottom, and the illustration is looking down
on the north pole of the star, such that the planet moves in a counter-
clockwise direction; transit occurs at the bottom, and secondary eclipse at
the top. The dotted lines denote the edges of the secondary eclipse and the
transit.

##  3 Analysis

###  3.1 Overall Procedure

After receiving the reduced and normalized PEPSI spectra, more steps are
necessary in order to detect the atmospheric signal. Our overall procedure is
inspired by that of Nugroho et al. (2017), although with significant
differences. We first run the molecfit package (Smette et al., 2015; Kausch et
al., 2015) on the red arm of each PEPSI spectrum to model out the telluric
lines, as is described in more detail in §3.3; no such treatment is necessary
for the largely telluric-free blue arm. We median-stack the telluric-corrected
spectra and subtract this median spectrum from each of the time-series spectra
in order to remove the stellar lines and the time-invariant component of the
telluric lines. We then run the SYSREM algorithm (Tamuz et al., 2005) on these
residuals, as described in §3.3. Next, we cross-correlate a model spectrum
generated as described in §3.2 with each spectrum. The cross-correlation
algorithm is adapted from that in the BANZAI-NRES
package222https://github.com/LCOGT/banzai-nres. These procedures, and the
tests that we conducted to show that this is the optimal method, are detailed
in §3.3. We show the spectra and the results of the molecfit in Fig. 2.

Figure 2: Example PEPSI data displaying our procedure to telluric-correct the
spectra. The figure shows one PEPSI spectrum from 2021-05-01 in black, with
the blue and red arms in the top and bottom panels, respectively. The
molecfit-corrected spectrum is shown in red for the red arm. For that arm we
also highlight the wavelength regions fit by molecfit in gray, and those which
we mask out when performing our analyses due to poor telluric correction in
red.

Finally, for a grid of possible radial velocity (RV) semi-amplitude $K_{P}$
values for the planet, we shift each of the time-series cross-correlation
functions to remove the expected Doppler shift from that $K_{P}$ at the
observation time; that is, we shift the spectra into the planetary rest frame.
We assumed the ephemeris from Lund et al. (2017). We then stack the resulting
shifted CCFs. For the transmission data we only combine the CCFs derived from
the in-transit spectra. In this step we combine the CCFs derived from the blue
and red arm data, and data taken on different nights. The exceptions are for
FeH, where there are no lines present in the blue arm so we use only the red
arm, and CaH, where we use only the blue arm due to the lack of red lines. We
combine the CCFs using a weighted sum, where the weights are the 95th
percentile SNR from that arm multiplied by the total equivalent width in the
model emission or transmission lines within that bandpass.

We search for a CCF peak near the expected values of $K_{P}$,
$v_{\mathrm{sys}}$ in this space, and measure the SNR of the peak (if any)
present at the location expected from the known parameters of the planet. We
estimate the significance of any peaks using the standard deviation of the
points in the shifted and combined CCFs with $|v|>100$ km s-1 from the
expected location of the planetary signal.

###  3.2 Model Spectra

We generate model spectra using the petitRADTRANS package (Mollière et al.,
2019). We assumed the stellar and planetary parameters from Lund et al.
(2017), which are listed in Table 2. We assumed a P-T profile of the form
given by Eqn. 29 of Guillot (2010) as implemented in petitRADTRANS, with
values of $\kappa_{\mathrm{IR}}=0.01$ and $\gamma=5$ chosen to approximate the
P-T profiles retrieved for the dayside atmosphere of KELT-20 b by Yan et al.
(2022b), Fu et al. (2022), and Borsa et al. (2022). We show this profile in
Fig. 3. For simplicity we use the same P-T profile for both the emission and
transmission spectra; although in reality the P-T profile is likely to be
different on the dayside versus the terminator, the dependence of the
transmission spectrum on the details of the P-T profile is weak.

Figure 3: Pressure-temperature profile used to generate our model spectra,
approximating that retrieved by Yan et al. (2022b), Fu et al. (2022), and
Borsa et al. (2022), showing a strong temperature inversion. Table 2: System
and Atmospheric Model Parameters

Parameter | Value | Source
---|---|---
$R_{P}$ ($R_{\oplus}$) | 19.51 | Lund et al. (2017)
$M_{P}$ ($M_{\oplus}$) | 1072 | Lund et al. (2017)
$T_{\mathrm{eq}}$ (K) | 2262 | Lund et al. (2017)
$P$ (d) | 3.4741085 | Lund et al. (2017)
$T_{0}$ (BJD_TDB) | 2457485.74965 | Lund et al. (2017)
$R_{\star}$ ($R_{\odot}$) | 1.565 | Lund et al. (2017)
$M_{\star}$ ($M_{\odot}$) | 1.76 | Lund et al. (2017)
$T_{\mathrm{eff}}$ (K) | 8720 | Lund et al. (2017)
$[$Fe/H$]$ | -0.29 | Lund et al. (2017)
$\kappa_{\mathrm{IR}}$ | 0.01 | …
$\gamma$ | 50 | …
$P_{0}$ (bar) | 1 | …
$X_{\mathrm{H}_{2}}$ | 0.7496 | …
$X_{\mathrm{He}}$ | 0.2504 | …
VMR(H-) | $1\times 10^{-9}$ | …

Note. —

We generate spectra over the wavelength range 3850-7500 Å, covering the full
PEPSI CD III+V bandpass, with an even spacing of 0.01 Å. We set up a model
atmosphere with 130 grid points spanning pressures from 102 to 10-10 bar.

We include continuum opacity due to $H^{-}$ bound-free absorption, H2-H2 and
H2-He collisions, and Rayleigh scattering due to H2 and He. We assume a cloud-
free atmosphere. We conservatively assumed an H- volume mixing ratio (VMR) of
$1\times 10^{-9}$, similar to that found by Arcangeli et al. (2018) for
WASP-18 b, a planet several hundred K hotter than KELT-20 b. This is the major
source of continuum opacity in the optical in our models, and is degenerate
with the VMR ratio limits that we set using our transmission spectra (§3.4;
§4.2). Including this level of H- absorption worsens the transmission VMR
limits by approximately two orders of magnitude over no H- absorption; this
has no significant effect on the emission results. We also note that we assume
a reference pressure level $P_{0}$ of 1 bar; changing $P_{0}$ to 1 mbar has a
similar effect upon the VMR limits to including the noted H- absorption. We
conservatively assume the high H- abundance in order to account for the
unknown $P_{0}$ level and H- abundance. The single line absorber in question
is the only line opacity source we include in our models.

The treatment of the emission versus transmission spectra is slightly
different. For the emission spectra, after calculating the planetary flux
using petitRADTRANS, we generate a blackbody spectrum corresponding to the
stellar $T_{\mathrm{eff}}$ and use this to estimate the planet-to-star flux
ratio $F_{P}/F_{\star}$. It is this quantity that we use as as our model
spectrum for the cross-correlation.

For the transmission spectra, petitRADTRANS calculates the planetary radius as
a function of wavelength $R_{P}(\lambda)$. As our spectra have been normalized
by the data reduction pipeline and are not flux-calibrated, we do not have any
information in our data on the planetary radius, only the differential change
in radius as a function of wavelength. We therefore estimate the transmission
spectrum as:

$F(\lambda)=1-\frac{R_{P}(\lambda)^{2}-R_{P0}^{2}}{R_{\star}^{2}}$ (1)

where $R_{P0}$ is the broadband planetary radius, which we adopt for Lund et
al. (2017).

We list the species considered and the sources of the opacity data used in
Table 3. We chose to search for the “classical” inversion agents TiO and VO
(Fortney et al., 2008), and the metal hydrides FeH and CaH. petitRADTRANS does
not currently have high-resolution opacity tables for several other species of
interest, such as AlO, CaO, NaH, or MgH (Gandhi & Madhusudhan, 2019). In the
interests of simplicity, for TiO we only use the single isotope 48TiO, as it
makes up the majority of Ti atoms, and ignore the other stable isotopes.
Likewise, for FeH the opacity table only considers the most common
isotopologue, but unlike for TiO there are no petitRADTRANS opacity tables
available for the other isotopologues.

Table 3: Opacity Data Sources

Species | Reference 1 | Reference 2
---|---|---
TiO | McKemmish et al. (2019) | B. Plez (priv. comm.)
VO | McKemmish et al. (2016) | B. Plez (priv. comm.)
FeH | Wende et al. (2010) | …
CaH | Li et al. (2012) | …
Fe i | Kurucz | …
Fe ii | Kurucz | …
Cr i | Kurucz | …
Si i | Kurucz | …

The reference Kurucz refers to data from http://kurucz.harvard.edu/.

Note. — Reference 1 is the opacity table used to generate the CCF template,
and Reference 2 is that used to generate the spectrum injected for injection-
recovery testing, if different. The opacity data for all species were obtained
from the petitRADTRANS
website333https://petitradtrans.readthedocs.io/en/latest/content/available_opacities.html

Figure 4: Model spectra generated with petitRADTRANS. The wavelength ranges
covered by PEPSI CD III and V are shown in blue and red, respectively. The
left column shows the emission spectra of, from top to bottom, Fe i, Fe ii, Si
i, and Cr i. The middle and right columns show the emission and transmission
spectra, respectively, of (top to bottom) TiO, VO, FeH, and CaH. The spectra
shown for the molecular constituents are generated with the minimum VMR which
would give a 4$\sigma$ detection in our injection-recovery testing; see the
main text for more details.

###  3.3 Comparison of Telluric and Systematics Correction Methods

In order to detect the minute planetary signal, it is necessary to correct for
systematic sources of error, most notably the presence of telluric absorption
lines in the spectra, but also systematics due to detector imperfections, etc.
Most recent studies of exoplanetary atmospheres at high resolution use one of
two methods to correct systematics: SYSREM (e.g., Nugroho et al., 2020;
Merritt et al., 2021; Cont et al., 2021a) or molecfit (e.g., Stangret et al.,
2021; Cabot et al., 2021; Kesseli et al., 2022). The SYSREM algorithm (Tamuz
et al., 2005) is a PCA-based algorithm to remove common linear or time-varying
systematics. The molecfit package (Smette et al., 2015; Kausch et al., 2015)
fits and removes the telluric lines in the spectrum. After implementing both
methods, we performed tests to evaluate which method works best for our PEPSI
data.

In order to implement SYSREM, we modified the source code from the PySysRem
package444https://github.com/stephtdouglas/PySysRem to handle the PEPSI data.
We added code to formally propagate uncertainties through the SYSREM
algorithm. For all runs, we used the airmass of the exposure as the initial
guess for the starting values of the first systematic removed. For the red
arm, we split the spectra into regions with telluric contamination versus no
telluric contamination, with a threshold of $>1\%$ absorption over multiple
nearby lines, varying the number of systematics removed for the former regions
and removing one systematic from the latter. Our tests indicated that the
recovered SNR did not change significantly if we removed two versus one
systematics from the telluric-free regions.

Although molecfit was designed to work on data from ESO facilities, it has
also been used for data from other spectrographs (e.g., Cabot et al., 2021).
We only used molecfit to correct the PEPSI red arm data, as there are no
significant telluric lines in the blue arm. We fit three narrow regions
spanning the red arm (6290-6320 Å, 6470-6520 Å, and 7340-7410 Å) containing
moderately-strong telluric lines, and included H2O and O2 in our models. While
the telluric model removes weak lines adequately, it does not fit the strong
lines to any good precision as they are saturated, and we consequently mask
out regions of the spectra containing poorly-corrected lines in our further
analysis. The telluric correction is sufficiently good for weak lines that any
remaining residuals can be removed well by SYSREM (see below).

For all our tests, we inject a model emission spectrum of 48TiO or FeH into
the PEPSI data. The large-scale line pattern of VO is similar to that of TiO,
and so we do not test it separately. We constructed the model spectra as
described in §3.2, using a VMR of $3.2\times 10^{-9}$ for TiO and $1.0\times
10^{-7}$ for FeH, chosen to give a moderate but not overwhelming detection
significance. For the transmission spectra we used VMRs of $3.2\times
10^{-11}$ and $5.6\times 10^{-8}$ for TiO and FeH, respectively. We passed
these spectra through SYSREM and then attempted to recover the signal using
the cross-correlation methodology described in §3.1. We evaluated different
numbers of SYSREM systematics removed in order to find the optimal number.

We show the results of these tests in Fig. 5. We find that, for TiO, using
only molecfit or only SYSREM to remove the tellurics offer similar
performance. These two methods result in detections with peak SNR of 7.6 and
7.4, respectively, for the 2021-05-01 data, and 7.4 and 7.6 for the 2021-05-18
data. SYSREM requires 10 and 7 systematics removed from the telluric-affected
regions in order to achieve the peak SNR.

We find, however, that we are able to achieve even better performance by using
a hybrid of the two techniques–i.e., first correcting the spectra with
molecfit and then correcting the remaining systematics with SYSREM. For
instance, the two datasets mentioned above, we can achieve peak SNRs of 8.7
and 7.8 after 2 and 3 systematics removed, respectively. We therefore proceed
with this hybrid approach for the TiO and VO for remainder of the paper. We
summarize the optimal number of systematics removed that we use for the
remainder of the paper in Table 4.

In contrast, for FeH we find the best performance when using only molecfit
(right panel of Fig. 5). Furthermore, we could not recover the FeH signal at
$>3\sigma$ using only SYSREM even after removing 20 systematics. We therefore
use only molecfit, without any SYSREM iterations, for FeH for the remainder of
this paper.

For the blue arm, we are unable to detect the injected TiO signal in the
2021-05-18 data. For the 2021-05-01 dataset, we reach a peak SNR of 5.8 with 3
systematics removed. We remove 3 systematics from all of the blue arm datasets
for the remainder of this paper. We did not inject the FeH signal in the blue
arm because it does not have any significant lines in that region of the
spectrum.

We performed a similar analysis for the transmission data; here we assumed a
VMR of$3.2\times 10^{-11}$ for TiO $5.6\times 10^{-8}$ for FeH. We show the
results in the middle panel of Fig. 5 for TiO and the bottom panel for FeH.
For both molecules we were unable to detect the red arm signal using molecfit
alone, likely due to higher telluric contamination than the other datasets.
Here we found overall similar results to the emission data, with the optimal
number of SYSREM iterations varying by no more than 2 between different dates
for the same species. We summarize the results and the number of iterations we
use going forward in Table 4.

Overall, the optimal strategy that we have used as is follows. For TiO and VO,
we use both molecfit and SYSREM. For FeH, we use only molecfit, while for CaH,
we use only SYSREM, since its lines are only present in the blue arm.

Figure 5: Results of the tests of the systematics-correction methods. Lines of
the same color show the same method and PEPSI arm: blue arm SYSREM (blue), red
arm SYSREM (red), red arm molecfit (dark red), and red arm molecfit \+ SYSREM
(orange). In the emission plots (top and bottom), circles and stars show the
2021-05-01 dataset, and squares and diamonds the 2021-05-18 dataset. The
x-axis shows the number of systematics that have been removed with SYSREM. The
top panel shows the results for the TiO injection in the emission data, the
middle panel the TiO injection in the transmission data, and the right for FeH
in the emission data. We could not detect the injected FeH emission signal or
the TiO transmission signal with SYSTEM even after removing 20 systematics,
and so the SYSREM-only lines are not shown on these plots. Table 4: SYSREM
Iterations Used

Date | Species | $N_{\mathrm{SYSREM}}$ | $N_{\mathrm{SYSREM}}$
---|---|---|---
|  | (blue arm) | (red arm)
2019-05-04 | TiO, VO, CaH | 4 | 1
2019-05-04 | FeH | … | 1
2021-05-01 | TiO, VO, CaH | 3 | 2
2021-05-01 | FeH | … | 0
2021-05-18 | TiO, VO, CaH | 3 | 3
2021-05-18 | FeH | … | 0

.

Note. —

###  3.4 Injection-Recovery Testing

In cases where we obtain a non-detection, it is useful to attempt to quantify
these non-detections. Do the data constrain the non-detection to an
interesting level, or are the data simply not high enough quality to place
useful limits? We attempt to address this issue through injecting a model
spectrum into our data before treating the data with SYSREM, and then
recovering it with our cross-correlation process described in §3.1. For each
species, we inject models with VMRs spaced 0.25 dex apart, i.e., VMR=[1.0,
1.8, 3.2, 5.6]$\times 10^{x}$, and find the smallest VMR for which we can
recover the signal at $>4\sigma$ (cf. Kesseli et al., 2020). We note that
these results are contingent upon the choice of P-T profile, as the strength
of the inversion affects the strength of the emission lines, and thus the
inferred VMR limit. However, Yan et al. (2022b), Fu et al. (2022), and Borsa
et al. (2022) all independently retrieved broadly similar P-T profiles using
independent methodology and data, giving us some confidence in the broad
accuracy of our P-T profile. Additionally, the presence of additional or
stronger continuum opacity sources than we have assumed could affect the
limits. These results should therefore only be taken as approximate under the
assumption of the P-T profile and other model atmospheric parameters.

One other potential issue is with the inaccuracy of line lists, which can
cause systematic reductions in signals strong enough to result in a non-
detection when a detection would otherwise be expected (e.g., de Regt et al.,
2022). In order to attempt to account for this issue, for species with
multiple opacity tables derived from different line lists available (i.e., TiO
and VO), we generate separate model spectra using the two line lists, inject
one into the data, and attempt to recover it using the other line list. We
have found that this worsens the detectable VMR by approximately 1 order of
magnitude.

##  4 Results

###  4.1 Recovery of Previous Detections

In order to perform a final validation of our methodology and data, we attempt
to reproduce the detections of emission from Fe i from Yan et al. (2022b) and
Borsa et al. (2022), Si i from Cont et al. (2021a), and Fe ii and Cr i from
Borsa et al. (2022). We estimated VMRs for these models using a FastChem
(Stock et al., 2018) equilibrium chemical model as implemented in
PyFastChem555https://github.com/exoclime/FastChem. We assumed a solar
abundance mixture, and the same P-T profile and pressure grid that we used to
construct our model spectra. We assumed the maximum VMRs present in the
FastChem models for these species, which range from $5.4\times 10^{-5}$ for Fe
i to $7.4\times 10^{-7}$ for Cr i. We show the model spectra in the left
panels of Fig. 4, and the CCFs shifted and stacked according to the expected
planetary orbit for all of these species in Fig. 6.

We were able to recover a signal from Fe i near the expected $K_{P}$, $v$
values at a significance of $15.1\sigma$. Yan et al. (2022b) obtained a
$7.7\sigma$ detection, and Borsa et al. (2022) at $7.1\sigma$. It is
unsurprising that we obtain a stronger detection, due to the much larger
aperture of LBT ($2\times 8$ m) versus the smaller telescopes used by those
works (3.5 m at Calar Alto Observatory, and 3.58 m at the Telescopio Nazionale
Galileo), although our bandpass is smaller than that of these instruments,
partially offsetting this advantage.

We were unable, however, to detect the other atomic species. Our Fe i
detection is so much stronger than that of Borsa et al. (2022) that we might
have been expected to easily detect Fe ii and Cr i as well. Indeed, injection-
recovery testing with VMRs of these species set to their FastChem values at 1
bar, similar to our procedure for the molecular species, suggest that we
should have detected both Fe ii and Cr i at high confidence, while Si i should
also have been detectable. Borsa et al. (2022) only detected Fe ii and Cr i in
their post-occultation dataset, not their pre-occultation dataset, which they
ascribed to a better viewing angle of a hot spot in the planetary atmosphere.
We verified that we did not detect these species when considering only our
post-occultation dataset (2021-05-01). There could be a number of reasons for
the non-detection, including our different bandpass and different choice of
P-T profile. We are preparing a companion paper on the atomic emission signals
in KELT-20 b and defer further discussion to that work.

Figure 6: Shifted and combined CCFs for emission from Fe i (top left), Si i
(top right), Fe ii (bottom left) and Cr i. We detect the Fe i emission at a
significance of $14.4\sigma$, reproducing the results of Yan et al. (2022b)
and Borsa et al. (2022), but are unable to reproduce the Si i detection from
Cont et al. (2021a), nor the detections of Fe ii and Cr i from Borsa et al.
(2022).

###  4.2 Limits on Inversion Agents

We searched for TiO, VO, FeH, and CaH emission and transmission in our spectra
of KELT-20 b, and did not make any detections with a signficance of
$>4\sigma$. Although the TiO transmission data do contain a bump near the
correct ($v$, $K_{P}$) peaking at $2.94\sigma$ (Fig. 8), due to the presence
of other similarly-strong peaks elsewhere in the ($v$, $K_{P}$) space we do
not consider this to be a robust detection. We show the data and minimum
recovered peaks in the Appendix.

Instead, we set limits upon the VMRs of these species as described in §3.4.
These limits range from $3.2\times 10^{-8}$ for FeH to $1.0\times 10^{-10}$
for TiO for the emission data, and from $1.8\times 10^{-8}$ for FeH to
$1.8\times 10^{-12}$ for TiO for the transmission data. We list the full
results in Table 5, and show the results graphically in Fig. 7.

We derive the pressure levels for the VMR limits shown in Fig. 7 from the
contribution functions calculated by petitRADTRANS; we calculate the pressure
of peak contribution and show as the error bars the range of pressures over
which $68\%$ of the contribution takes place. Counterintuitively, in our
models the transmission spectra probe deeper into the atmosphere than the
emission spectra. This is because, in our models, the lack of other optical
opacity sources besides the absorber in question allows the transmission rays
to penetrate deep into the atmosphere, below the temperature inversion. The
emission lines, on the other hand, by definition are formed only in the
inversion. Although the nominal transmission and emission models are
consistent with each other, they are unlikely to be completely realistic. We
also calculated spectra with an artificially increased VMR of H- to $1\times
10^{-9}$, similar to the maximum VMR found by Arcangeli et al. (2018) for
WASP-18 b, a planet $\sim 600$ K hotter than KELT-20 b. Including this
artificially increased H- absorption does not significantly affect the
emission results, but for the transmission results the limits worsen by 1-2
orders of magnitude, and the contribution moves by $\sim 1$ dex lower in
pressure. These points are shown in Fig. 7; they do not qualitatively affect
our overall conclusions.

Table 5: Summary of Results

Species | Emission VMR limit | Recovered SNR | Transmission VMR limit | Recovered SNR | Expected VMR
---|---|---|---|---|---
TiO | $1\times 10^{-10}$ | 4.38 | $1.8\times 10^{-10}$ | 5.01 | $1.4\times 10^{-7}$
VO | $1\times 10^{-9}$ | 5.48 | $1.0\times 10^{-9}$ | 5.06 | $8.0\times 10^{-9}$
FeH | $3.2\times 10^{-8}$ | 5.38 | $3.2\times 10^{-7}$ | 5.43 | $9.3\times 10^{-9}$
CaH | $1.8\times 10^{-9}$ | 5.51 | $5.6\times 10^{-10}$ | 4.64 | $3.4\times 10^{-7}$

.

Note. — The “recovered SNR” is the SNR of the CCF peak recovered for the
injected model spectrum at the quoted VMR.

In order to interpret the results, we used a FastChem (Stock et al., 2018)
equilibrium chemical model as implemented in
PyFastChem666https://github.com/exoclime/FastChem. We constructed two models,
one assuming a solar abundance mixture and the other $10\times$ solar
metallicity, with the same P-T profile and pressure grid that we used to
construct our model spectra. We assumed that quenching occurs at a pressure of
1 bar, such that the VMRs of our species of interest are constant at $P<1$ bar
(cf. Zahnle & Marley, 2014; Marley & Robinson, 2015). FastChem does not
include photoionization or photodissociation, which may be important processes
in the upper atmosphere of KELT-20 b. Photodissociation should serve to
depress the concentrations of molecules including inversion agents, while
photoionization should increase the concentrations of ionized species and H-.
These effects are counteracted by the replenishing process of vertical mixing
when we assume a quench pressure at 1 bar. Although these assumptions
regarding quenching and the model’s disregard of irradiation effects are
likely to be only approximately correct, this nonetheless serves as a useful
tool to help interpret our limits.

We show these predicted VMRs from FastChem as vertical lines in Fig. 7. Our
VMR limits for TiO, VO, and CaH lie several orders of magnitude below the
expectations for those species at both metallicities. Our limits for FeH lie
above the expected concentrations for solar metallicity, and below for
$10\times$ solar. Most of the strong FeH lines lie redward of the PEPSI CD V
bandpass; new observations with PEPSI CD VI (which covers 7419-9067 Å;
Strassmeier et al., 2015) would be required to better constrain FeH with
PEPSI. We again note that our VMR limits are only approximate, and in detail
subject to our assumptions regarding the P-T profile, continuum opacity from
H-, and reference pressure level. Nonetheless, the large mismatch between the
expected VMRs and our limits for TiO, VO, and CaH suggest that these species
are depleted in the atmosphere of KELT-20 b and cannot be responsible for the
temperature inversion, while FeH may also be depleted if KELT-20 b’s
atmosphere has a super-solar metallicity.

Figure 7: Graphical summary of our limits on the abundances of inversion
agents. For each species, the left-pointing arrow shows our upper limit set
upon the VMR, while the error bar shows the range in altitudes which
contribute 68% of the emission or absorption from that species within the
PEPSI bandpass at the VMR limit. The filled point for each species is from the
emission data, and the empty points from the transmission data. The full-color
transmission points are for the nominal models, while the transparent points
show the effects of artificially increasing the H- VMR to $1\times 10^{-9}$
(see §3.2). The vertical lines show the VMRs expected from a chemical
equilibrium calculation at 1 bar, assuming that this is the quenching pressure
and that the VMR is constant above this point. Solid lines assume solar
metallicity, and dotted lines $10\times$ solar. The colors of the lines
correspond to those of the data points for each species. If the data point is
to the left of the colored bar, the observation is in conflict with the
expectation.

##  5 Conclusions

We presented observations of the atmosphere of KELT-20 b with LBT/PEPSI in
both emission and transmission. We recovered the previously-detected Fe i
emission from the planetary atmosphere, but were unable to confirm the
previous detections of Si i, Cr i, or Fe ii. Given abundances from FastChem
and injection-recovery testing we should have been able to detect these
species, but the somewhat different P-T profile that we assumed as compared to
previous works could have altered the line strengths enough to result in a
non-detection.

We set strict limits upon the presence of the possible inversion agents TiO,
VO, FeH, and CaH in the atmosphere of KELT-20 b (Fig. 7, Table 5), although
these limits are subject to our assumptions regarding the P-T profile and
continuum opacity. Nonetheless, together with expectations from simple
chemical models, this suggests that TiO, VO, and CaH are depleted in the
atmosphere of KELT-20 b and cannot be responsible for the temperature
inversion. If the atmosphere of KELT-20 b has a significantly super-solar
metallicity, then FeH is also depleted. Either a different molecular agent is
responsible for the inversion in KELT-20 b, or, as suggested by Nugroho et al.
(2020), atomic species such as Fe i could be responsible. Alternately, if
there is indeed significant H- opacity in the atmosphere of KELT-20 b it could
be the opacity responsible for the inversion (Lothringer et al., 2018).

In this paper we have presented our results on KELT-20 b. As part of the PETS
survey we have already gathered similar observations of other UHJs, with the
aim of conducting a population study of inversion agents, which we will
present in a future paper along with improved analysis techniques. These data
will also be more broadly useful for characterizing these planets. Given the
strength of our Fe i detection, we may be able to use the time-resolved
emission to constrain the night-side brightness and hotspot offset of KELT-20
b much like Herman et al. (2022) were able to do for WASP-33 b, but doing so
is beyond the scope of the present paper and will be presented in a future
work along with further analysis of the other atomic species.

Thanks to Wilson Cauley for technical assistance, and Francesco Borsa and
Lorenzo Pino for useful discussions.

B.S.G. was partially supported by a Thomas Jefferson Chair for Space
Exploration endowment at the Ohio State University. This research was
supported by the Excellence Cluster ORIGINS which is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311.

The LBT is an international collaboration among institutions in the United
States, Italy and Germany. LBT Corporation partners are: The Ohio State
University, representing OSU, University of Notre Dame, University of
Minnesota and University of Virginia; The University of Arizona on behalf of
the Arizona Board of Regents; Istituto Nazionale di Astrofisica, Italy; and
LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society,
The Leibniz Institute for Astrophysics Potsdam, and Heidelberg University.

LBT (PEPSI)

## References

  * Arcangeli et al. (2018) Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30, doi: 10.3847/2041-8213/aab272
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., SipH ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Baxter et al. (2020) Baxter, C., Désert, J.-M., Parmentier, V., et al. 2020, A&A, 639, A36, doi: 10.1051/0004-6361/201937394
  * Bello-Arufe et al. (2022) Bello-Arufe, A., Buchhave, L. A., Mendonça, J. M., et al. 2022, arXiv e-prints, arXiv:2203.04969.  https://arxiv.org/abs/2203.04969
  * Borsa et al. (2022) Borsa, F., Giacobbe, P., Bonomo, A. S., et al. 2022, arXiv e-prints, arXiv:2204.04948.  https://arxiv.org/abs/2204.04948
  * Cabot et al. (2021) Cabot, S. H. C., Bello-Arufe, A., Mendonça, J. M., et al. 2021, AJ, 162, 218, doi: 10.3847/1538-3881/ac1ba3
  * Casasayas-Barris et al. (2018) Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2018, A&A, 616, A151, doi: 10.1051/0004-6361/201832963
  * Casasayas-Barris et al. (2019) —. 2019, A&A, 628, A9, doi: 10.1051/0004-6361/201935623
  * Cont et al. (2021a) Cont, D., Yan, F., Reiners, A., et al. 2021a, arXiv e-prints, arXiv:2112.10461.  https://arxiv.org/abs/2112.10461
  * Cont et al. (2021b) —. 2021b, A&A, 651, A33, doi: 10.1051/0004-6361/202140732
  * de Regt et al. (2022) de Regt, S., Kesseli, A. Y., Snellen, I. A. G., Merritt, S. R., & Chubb, K. L. 2022, arXiv e-prints, arXiv:2203.03585.  https://arxiv.org/abs/2203.03585
  * Evans et al. (2017) Evans, T. M., Sing, D. K., Kataria, T., et al. 2017, Nature, 548, 58, doi: 10.1038/nature23266
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419, doi: 10.1086/528370
  * Fu et al. (2022) Fu, G., Sing, D. K., Lothringer, J. D., et al. 2022, arXiv e-prints, arXiv:2201.02261.  https://arxiv.org/abs/2201.02261
  * Gandhi & Madhusudhan (2019) Gandhi, S., & Madhusudhan, N. 2019, MNRAS, 485, 5817, doi: 10.1093/mnras/stz751
  * Guillot (2010) Guillot, T. 2010, A&A, 520, A27, doi: 10.1051/0004-6361/200913396
  * Herman et al. (2020) Herman, M. K., de Mooij, E. J. W., Jayawardhana, R., & Brogi, M. 2020, AJ, 160, 93, doi: 10.3847/1538-3881/ab9e77
  * Herman et al. (2022) Herman, M. K., de Mooij, E. J. W., Nugroho, S. K., Gibson, N. P., & Jayawardhana, R. 2022, arXiv e-prints, arXiv:2203.11188.  https://arxiv.org/abs/2203.11188
  * Hoeijmakers et al. (2020a) Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020a, A&A, 641, A123, doi: 10.1051/0004-6361/202038365
  * Hoeijmakers et al. (2020b) Hoeijmakers, H. J., Cabot, S. H. C., Zhao, L., et al. 2020b, A&A, 641, A120, doi: 10.1051/0004-6361/202037437
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ilyin (2000) Ilyin, I. V. 2000, PhD thesis, Astronomy Division Department of Physical Sciences P.O.Box 3000 FIN-90014 University of Oulu Finland
  * Kasper et al. (2021) Kasper, D., Bean, J. L., Line, M. R., et al. 2021, ApJ, 921, L18, doi: 10.3847/2041-8213/ac30e1
  * Kausch et al. (2015) Kausch, W., Noll, S., Smette, A., et al. 2015, A&A, 576, A78, doi: 10.1051/0004-6361/201423909
  * Keles et al. (2022) Keles, E., Mallonn, M., Kitzmann, D., et al. 2022, MNRAS, doi: 10.1093/mnras/stac810
  * Kesseli et al. (2020) Kesseli, A. Y., Snellen, I. A. G., Alonso-Floriano, F. J., Mollière, P., & Serindag, D. B. 2020, AJ, 160, 228, doi: 10.3847/1538-3881/abb59c
  * Kesseli et al. (2022) Kesseli, A. Y., Snellen, I. A. G., Casasayas-Barris, N., Mollière, P., & Sánchez-López, A. 2022, AJ, 163, 107, doi: 10.3847/1538-3881/ac4336
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Lebigot (2017) Lebigot, E. O. 2017, Uncertainties: a Python package for calculations with uncertainties.  http://pythonhosted.org/uncertainties/
  * Li et al. (2012) Li, G., Harrison, J. J., Ram, R. S., Western, C. M., & Bernath, P. F. 2012, J. Quant. Spec. Radiat. Transf., 113, 67, doi: 10.1016/j.jqsrt.2011.09.010
  * Lothringer et al. (2018) Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27, doi: 10.3847/1538-4357/aadd9e
  * Lund et al. (2017) Lund, M. B., Rodriguez, J. E., Zhou, G., et al. 2017, AJ, 154, 194, doi: 10.3847/1538-3881/aa8f95
  * Marley & Robinson (2015) Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279, doi: 10.1146/annurev-astro-082214-122522
  * McKemmish et al. (2019) McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019, MNRAS, 488, 2836, doi: 10.1093/mnras/stz1818
  * McKemmish et al. (2016) McKemmish, L. K., Yurchenko, S. N., & Tennyson, J. 2016, MNRAS, 463, 771, doi: 10.1093/mnras/stw1969
  * Merritt et al. (2021) Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2021, MNRAS, 506, 3853, doi: 10.1093/mnras/stab1878
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67, doi: 10.1051/0004-6361/201935470
  * Nugroho et al. (2020) Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, MNRAS, 496, 504, doi: 10.1093/mnras/staa1459
  * Nugroho et al. (2017) Nugroho, S. K., Kawahara, H., Masuda, K., et al. 2017, AJ, 154, 221, doi: 10.3847/1538-3881/aa9433
  * Prinoth et al. (2022) Prinoth, B., Hoeijmakers, H. J., Kitzmann, D., et al. 2022, Nature Astronomy, 6, 449, doi: 10.1038/s41550-021-01581-z
  * Serindag et al. (2021) Serindag, D. B., Nugroho, S. K., Mollière, P., et al. 2021, A&A, 645, A90, doi: 10.1051/0004-6361/202039135
  * Sheppard et al. (2017) Sheppard, K. B., Mandell, A. M., Tamburo, P., et al. 2017, ApJ, 850, L32, doi: 10.3847/2041-8213/aa9ae9
  * Smette et al. (2015) Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77, doi: 10.1051/0004-6361/201423932
  * Stangret et al. (2021) Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2021, arXiv e-prints, arXiv:2111.13064.  https://arxiv.org/abs/2111.13064
  * Stangret et al. (2020) —. 2020, A&A, 638, A26, doi: 10.1051/0004-6361/202037541
  * Stock et al. (2018) Stock, J. W., Kitzmann, D., Patzer, A. B. C., & Sedlmayr, E. 2018, MNRAS, 479, 865, doi: 10.1093/mnras/sty1531
  * Strassmeier et al. (2018) Strassmeier, K. G., Ilyin, I., & Steffen, M. 2018, A&A, 612, A44, doi: 10.1051/0004-6361/201731631
  * Strassmeier et al. (2015) Strassmeier, K. G., Ilyin, I., Järvinen, A., et al. 2015, Astronomische Nachrichten, 336, 324, doi: 10.1002/asna.201512172
  * Talens et al. (2018) Talens, G. J. J., Justesen, A. B., Albrecht, S., et al. 2018, A&A, 612, A57, doi: 10.1051/0004-6361/201731512
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466, doi: 10.1111/j.1365-2966.2004.08585.x
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Wende et al. (2010) Wende, S., Reiners, A., Seifahrt, A., & Bernath, P. F. 2010, A&A, 523, A58, doi: 10.1051/0004-6361/201015220
  * Wong et al. (2021) Wong, I., Kitzmann, D., Shporer, A., et al. 2021, AJ, 162, 127, doi: 10.3847/1538-3881/ac0c7d
  * Yan et al. (2020) Yan, F., Pallé, E., Reiners, A., et al. 2020, A&A, 640, L5, doi: 10.1051/0004-6361/202038294
  * Yan et al. (2022a) —. 2022a, A&A, 661, L6, doi: 10.1051/0004-6361/202243503
  * Yan et al. (2022b) Yan, F., Reiners, A., Pallé, E., et al. 2022b, A&A, 659, A7, doi: 10.1051/0004-6361/202142395
  * Zahnle & Marley (2014) Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41, doi: 10.1088/0004-637X/797/1/41

##  Appendix A Cross-Correlation Functions for Inversion Agents

In this Appendix we show the shifted and combined cross-correlation functions
for the inversion agents TiO (Fig. 8), VO (Fig. 9), FeH (Fig. 10), and CaH
(Fig. 11).

Figure 8: Shifted and combined CCFs for TiO, in emission (left panels) and
transmission (right panels). The top panels show the data, and the bottom
panels the data with a signal from the minimum VMR that can be recovered at
$>4\sigma$ injected.

Figure 9: Same as Fig. 8 but for VO.

Figure 10: Same as Fig. 8 but for FeH.

Figure 11: Same as Fig. 8 but for CaH.


# Strain-induced phase diagram of the $S=\frac{3}{2}$ Kitaev material
$\rm{CrSiTe_{3}}$

Zongsheng Zhou  School of Physical Science and Technology $\&$ Key Laboratory
for Magnetism and Magnetic Materials of the MoE, Lanzhou University, Lanzhou
730000, China  Lanzhou Center for Theoretical Physics and Key Laboratory of
Theoretical Physics of Gansu Province, Lanzhou University, Lanzhou 730000,
China.  Ken Chen  School of Physical Science and Technology $\&$ Key
Laboratory for Magnetism and Magnetic Materials of the MoE, Lanzhou
University, Lanzhou 730000, China  Lanzhou Center for Theoretical Physics and
Key Laboratory of Theoretical Physics of Gansu Province, Lanzhou University,
Lanzhou 730000, China.  Qiang Luo  Department of Physics, University of
Toronto, Toronto, Ontario M5S 1A7, Canada  Hong-Gang Luo<EMAIL_ADDRESS>School of Physical Science and Technology $\&$ Key Laboratory for Magnetism
and Magnetic Materials of the MoE, Lanzhou University, Lanzhou 730000, China
Lanzhou Center for Theoretical Physics and Key Laboratory of Theoretical
Physics of Gansu Province, Lanzhou University, Lanzhou 730000, China.  Beijing
Computational Science Research Center, Beijing 100084, China  Jize Zhao
<EMAIL_ADDRESS>School of Physical Science and Technology $\&$ Key
Laboratory for Magnetism and Magnetic Materials of the MoE, Lanzhou
University, Lanzhou 730000, China  Lanzhou Center for Theoretical Physics and
Key Laboratory of Theoretical Physics of Gansu Province, Lanzhou University,
Lanzhou 730000, China.

###### Abstract

The interplay among anisotropic magnetic terms, such as the bond-dependent
Kitaev interactions and single-ion anisotropy, plays a key role in stabilizing
the finite-temperature ferromagnetism in the two-dimensional compound
$\rm{CrSiTe_{3}}$. While the Heisenberg interaction is predominant in this
material, a recent work shows that it is rather sensitive to the compressive
strain, leading to a variety of phases, possibly including a sought-after
Kitaev quantum spin liquid [C. Xu, et. al., Phys. Rev. Lett. 124, 087205
(2020)]. To further understand these states, we establish the quantum phase
diagram of a related bond-directional spin-$3/2$ model by the density-matrix
renormalization group method. As the Heisenberg coupling varies from
ferromagnetic to antiferromagnetic, three magnetically ordered phases, i.e., a
ferromagnetic phase, a $120^{\circ}$ phase and an antiferromagnetic phase,
appear consecutively. All the phases are separated by first-order phase
transitions, as revealed by the kinks in the ground-state energy and the jumps
in the magnetic order parameters. However, no positive evidence of the quantum
spin liquid state is found and possible reasons are discussed briefly.

##  I Introduction

The bond-directional interactions have emerged as a focus of research in
hunting for nontrivial quantum states ever since Kitaev proposed his famous
spin-$1/2$ model on the honeycomb lattice in 2006, later dubbed Kitaev
honeycomb model Kitaev2006 . As a rare exactly solvable model which supports
non-Abelian statistics and Majorana fermion excitations, the Kitaev model is
believed to be a milestone in the area of quantum spin liquids (QSLs)
Anderson1973 ; Balents2010 ; WenNPJ2019 . However, experimental realization of
such interactions is highly nontrivial due to their bond-dependent Ising
nature. Later on, Jackeli and Khaliullin provided a feasible routine to
realize these interactions by the collaboration of spin-orbit coupling and
electron correlations Jackeli2009 . Of particular interest are the $5d^{5}$
iridates $A\rm{{}_{2}IrO_{3}}$$(A\rm{=Na,Li)}$ Chaloupka2010  and
$\alpha$-$\rm{RuCl_{3}}$ PlumbKW2014 ; KimHeungSik2015 ; Banerjee2016 , which
are recognized as effective $j=1/2$ Mott insulators with essential Kitaev
interactions. These works soon ignited prosperous experimental and theoretical
studies on such Kitaev materials (see WangDongetal2017 ; WinterNcom2018 ;
Maksimov2020 ; LiZhangWang2021 ; HermannsM2018 ; Takagi2019  and references
therein).

On the other hand, the non-Kitaev interactions are usually nonnegligible in
candidate materials. For example, the Heisenberg ($J$) interaction is
assumable to be large in $\rm{Na_{2}IrO_{3}}$ Chaloupka2010 , while the
symmetric off-diagonal $\Gamma$ term Rau2014a , together with the distortion-
induced $\Gamma^{\prime}$ term Rau2014b , is essential in
$\alpha$-$\rm{RuCl_{3}}$. This leads to a prototypical $\rm
JK\Gamma\Gamma^{\prime}$ model, which serves as a versatile playground to
explore exotic phases and collective phenomena. These include a multi-node QSL
in the vicinity of the Kitaev limit WangBL2019 , a gapless QSL proposed in the
honeycomb $\Gamma$ model LuoZhaoetal2021 ; CatunYWetal2018 ; GohlkeWYetal2018
, and a chiral-spin ordering favored by a tiny $\Gamma^{\prime}$ interaction
in the dominated $\Gamma$ region LuoCSL2020 .

While the spin-$1/2$ iridates and $\alpha$-$\rm{RuCl_{3}}$ have been
extensively studied, recent efforts are also extended to $3d^{7}$
SanoRyoya2018 ; LiuHuimei2018 ; Vivanco2020 ; Songvilay2020  and $f^{1}$
JangSeongHoon2019 ; JangSeongHoon2020 ; Motome2020B  electron states. In these
materials holding large effective spins, bond-directional Ising interactions
have also been proposed Stavropoulos2019 ; Stavropoulos2021 . Among them
$X\rm{{}_{3}Ni_{2}SbO_{6}}$$\left(X=\rm{Li,Na}\right)$ Zvereva2015 ;
Kurbakov2017  are candidates of $S=1$ Kitaev materials, while
$\rm{CrXTe_{3}\left(X=Si,Ge\right)}$ and $\rm{CrI_{3}}$ are $S=3/2$ ones
XuChangsong2018 ; XuChangsong2020 ; Lee2020 ; ChenLebing2020 ; Bacaksiz2021 .
In these high-spin materials, a single-ion anisotropy (SIA) coming from the
partially unquenched orbital moment of the magnetic ion should also be taken
into account XuChangsong2018 ; XuChangsong2020 ; Bacaksiz2021 ; Lado2017 .

Recently, density functional theory calculations on $\rm{CrSiTe_{3}}$ show
that the Heisenberg interaction can be tuned significantly while the other
magnetic terms are nearly unchanged under compressive strain XuChangsong2018 ;
XuChangsong2020 ; Bacaksiz2021 . Interestingly, a possible Kitaev QSL state,
which is inferred from the double-peak specific heat anomaly and the
associated thermal entropy plateau, is proposed to sit between a ferromagnetic
(FM) phase and an antiferromagnetic (AFM) phase XuChangsong2020 . Their
argument in mind is that the Kitaev QSL is known to display two peaks in the
specific heat NasuUM2015 ; KogaTN2018 , which are the reminiscence of the
itinerant Majorana fermions and vison excitations KnolleKCM2014 .
Nevertheless, the proposed QSL state is far from mature for the following
reasons. In the Kitaev model, intensity of the low-temperature peak in the
specific heat is usually comparable to or is even larger than the high-
temperature peak NasuUM2015 ; KogaTN2018 ; LuoHuKee2021 . However, it is
reported in Ref. [XuChangsong2020, ] that the low-temperature peak is rather
weak in $\rm{CrSiTe_{3}}$. In addition, occurrence of the double-peak specific
heat does not necessarily imply a QSL ground state, and several magnetically
ordered states are demonstrated to possess such a specific heat anomaly in the
extended Kitaev-Heisenberg model YamajiSY2016 . These observations motivate us
to revisit the QSL state in the relevant parameter region.

Figure 1: (a) A sketch diagram of the top view of the ${\rm CrSiTe_{3}}$
monolayer, a unit cell is marked by the dashed parallelogram. ${\rm Cr}$ atoms
locate at the centers of octahedra formed by six ${\rm Te}$ atoms, and these
${\rm Cr}$ atoms form a two-dimensional honeycomb network. In the
$\left\\{{XYZ}\right\\}$ coordinate system, the three directions are marked by
cyan, blue, and magenta arrows, respectively. The exchange matrix on each
color bond represented by $X,Y$ and $Z$ is given in the main text. The
orthogonal $\left\\{abc\right\\}$ coordinate system is given as follows. ${a}$
and ${b}$ are perpendicular and parallel to the Z bonds, respectively, ${c}$
is in the $\left[111\right]$ direction that is perpendicular to the honeycomb
network. (b) The 24-site $C_{6}$ cluster ($C_{6}$ is the symmetry of the
honeycomb lattice, not the model). Black dots represent $S=\frac{3}{2}$ spins,
three different Ising bonds are shown in different colors as in panel (a).
Similarly, panel (c) plots 48-site cluster, $L_{x}$ is perpendicular to the Z
bonds while $L_{y}$ is parallel to them. We adopt periodic boundary condition
along $L_{y}$ direction to roll the lattice into a cylinder, which is refered
as the YC cylinder. In (b) and (c), the dashed line at the edges means the
periodic boundary condition while the solid line means open boundary
condition. The snake-like light gray line illustrates the one-dimensional
mapping in our DMRG calculations.

In this work, we investigate the model proposed to describe $\rm{CrSiTe_{3}}$
by the density-matrix renormalization group (DMRG) method White1992 ;
Peschel1999 ; Schollwock2005 . As we tune the Heisenberg interaction from FM
to AFM as a mimicry of the strain process, the ground state changes from the
FM phase to the AFM phase, through a wide intermediate region of $120^{\circ}$
phase rather than a QSL proposed in the previous work on a 12-site cluster
XuChangsong2020 . To illuminate the discrepancy, we perform DMRG calculations
on four different 12-site clusters at the proposed parameter. Based on the
results, we conjecture that the evidence of the possible QSL state may
originate from the artificial frustration induced by the mismatch between the
$120^{\circ}$ order on small cluster size.

This paper is organized as follows. In Sec. II, we briefly introduce the model
and present some crucial details of the numerical calculations. In Sec. III
and Sec. IV, we present the results of DMRG calculation and classical Monte
Carlo simulation, respectively. In Sec. V we show the DMRG results on four
different 12-site clusters for comparison. Finally, we summarize our work in
Sec.VI.

##  II Model and methods

The crystal structure of $\rm{CrSiTe_{3}}$ is illustrated in Fig. 1(a). The
$\rm{Cr^{3+}}$ ions form a two-dimensional honeycomb network, and the
effective low-energy state of $\rm{Cr^{3+}}$ ion is described by an $S=3/2$
spin. Anisotropic interactions resulting from the strong spin-orbit coupling
of the heavy ligand of $\rm{Cr}$, such as the Kitaev interaction, the off-
diagonal $\rm{\Gamma}$ and $\rm{\Gamma^{\prime}}$ interactions as well as the
SIA, are essential to stabilize the ferromagnetism at finite temperature. A
possible minimal model Hamiltonian up to nearest-neighbor interactions is
given by XuChangsong2018 ; XuChangsong2020 ; Stavropoulos2021 ; Bacaksiz2021

$\mathcal{H}=\sum_{\langle
ij\rangle}\mathbf{S}_{i}^{\rm{T}}\mathcal{J}_{ij}^{\alpha}\mathbf{S}_{j}+\sum_{i}\mathbf{S}_{i}^{\rm{T}}\mathcal{A}\mathbf{S}_{i}.$
(1)

The first term in Eq. (1) is known as the $\rm{JK\Gamma\Gamma^{\prime}}$
model. Here, $\langle{ij}\rangle$ are summed over nearest neighbors, and
$\alpha\left({=X,Y,Z}\right)$ represents the bond shown in Fig. 1(a). In the
second term the sum runs over all sites.
$\mathbf{S}_{i}=\left(S_{i}^{\rm{X}},S_{i}^{\rm{Y}},S_{i}^{\rm{Z}}\right)^{\rm{T}}$
represent the spin operators at the site $i$. The bond-dependent interactions
and SIA matrices in the $\\{XYZ\\}$ coordinate system have the following form
XuChangsong2020

$\mathcal{J}^{\rm{X}}=\left(\begin{array}[]{ccc}J+K&\Gamma^{\prime}&\Gamma^{\prime}\\\
\Gamma^{\prime}&J&\Gamma\\\
\Gamma^{\prime}&\Gamma&J\end{array}\right),\qquad\mathcal{J}^{\rm{Y}}=\left(\begin{array}[]{ccc}J&\Gamma^{\prime}&\Gamma\\\
\Gamma^{\prime}&J+K&\Gamma^{\prime}\\\
\Gamma&\Gamma^{\prime}&J\end{array}\right)$
$\mathcal{J}^{\rm{Z}}=\left(\begin{array}[]{ccc}J&\Gamma&\Gamma^{\prime}\\\
\Gamma&J&\Gamma^{\prime}\\\
\Gamma^{\prime}&\Gamma^{\prime}&J+K\end{array}\right),\qquad\mathcal{A}=\frac{A_{cc}}{3}\left(\begin{array}[]{ccc}1&1&1\\\
1&1&1\\\ 1&1&1\end{array}\right).$

Of important note is that in the crystallographic ${abc}$ coordinate system,
the SIA term could be rewritten elegantly as $A_{cc}(S_{i}^{c})^{2}$. For
$S=3/2$ Kitaev materials, the mechanism for anisotropic interactions is
different from the $S=1/2$ case Jackeli2009 ; XuChangsong2018 ;
Stavropoulos2021 . Therefore, the Kitaev interaction can be AFM and the off-
diagonal $\Gamma$ and $\Gamma^{\prime}$ terms could be comparable Maksimov2020
; XuChangsong2020 . We focus on the ground phase diagram of $\rm{CrSiTe_{3}}$
induced by compressive strain, and the parameters in (1) are referred from the
supplemental material of Ref. XuChangsong2020 . Since $K$, $\Gamma$,
$\Gamma^{\prime}$ and $A_{cc}$ are almost independent of the strength of the
strain, we take them as constants by setting $K=0.275$ meV, $\Gamma=0.020$
meV, $\Gamma^{\prime}=-0.083$ meV, $A_{cc}=0.222$ meV. Thus, $J$, in unit of
meV, is the only variable in the phase diagram.

Figure 2: (a) Ground-state energy per site $E_{g}/N$ is shown as functions of
$J$ for two different clusters (see Fig. 1(b) and (c)). In both cases, we can
find three regions separated by kinks in the ground-state energy. The
positions of these kinks for the two different clusters are almost the same,
which are, from left to right, $J=-0.282(3)/-0.274(2)$ and $0.022(2)$,
respectively. The transition point separating the FM phase and the
$120^{\circ}$ phase has a minor difference in these two clusters. This
distinction is marked out in light gray. (b) The typical SSSF in each phase on
$C_{6}$ cluster is shown. (c) The order parameters $M$ (see the text) are
shown as functions of $J$ for the YC cluster and ${C_{6}}$ cluster,
respectively. The characteristic momenta to extract the order parameters from
the SSSF are $\mathcal{Q}_{\Gamma}$, $\mathcal{Q}_{{\rm{K}}}$, and
$\mathcal{Q}_{\Gamma^{\prime}}$, which correspond to the FM, $120^{\circ}$ and
AFM orders, respectively. The positions of the discontinuities in the dominant
order parameters and the kinks in the ground-state energy agree well.

The phase diagram of the model (1) induced by the strain is identified through
the ground-state energy, order parameters as well as static spin structure
factors (SSSFs), which are readily available in our DMRG calculation. In the
simulations, we use different clusters to confirm that our numerical results
are reliable. One should be careful whether the clusters match the magnetic
orders. As shown in Fig. 1(b) and Fig. 1(c), 24-site $C_{6}$ cluster with
periodic boundary condition and 48-site YC cylinder ($6\times 4\times 2$) are
mainly considered in our simulations. In the calculations, the number of
states we keep ranges from 1200 to 1600, and at least 12 sweeps are
implemented in each calculation to ensure the convergence of DMRG results.
After the calculations are converged, the utmost magnitude of truncation error
in each sweep is $10^{-5}$.

##  III Quantum Phase Diagram

Our DMRG results are summarized in Fig. 2. In panel (a), we plot the ground-
state energy per site $E_{g}/N$ as a function of $J$ on the 24-site ${C_{6}}$
cluster (Fig. 1(b)) and the 48-site YC cylinder (Fig. 1(c)). Three phases are
clearly detected which are separated by kinks in the energy curve. A kink in
the energy curve indicates the level crossing and thus represents a first-
order phase transition. Moreover, the positions of these kinks agree well on
these two clusters. It is interesting to notice that the exact diagonalization
(ED) on a small cluster with 12 sites XuChangsong2020  already provides an
excellent estimate of the phase boundaries. Such a weak finite-size effect is
the typical feature of first-order phase transitions. These three phases turn
out to have FM, $120^{\circ}$ and AFM orders. One may notice that the energy
of the $C_{6}$ cluster is lower than that of the YC cluster. This is because
in the $C_{6}$ cluster PBC is used and more bonds are connected which
contribute negative energy to $E_{g}$.

To gain more information on these phases, we calculate the SSSF which is
defined as

$S^{\alpha\alpha}\left(\mathcal{Q}\right)=\frac{1}{N^{2}}\sum_{ij}\langle
S_{i}^{\alpha}S_{j}^{\alpha}\rangle
e^{i\mathcal{Q}\cdot\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)}$

where $N$ is the number of sites and $i$, $j$ run over all the sites.
$\langle\cdots\rangle$ is the expectation in the ground state. The order
parameter at the momenta $\mathcal{Q}$ is then given by
$M=\sqrt{\sum_{\alpha}{{S}^{\alpha\alpha}(\mathcal{Q})}}$. In the magnetic
ordered phase, it is finite at the characteristic momentum in the
thermodynamic limit. Therefore, we can tell the magnetic orders at the
characteristic momenta determined by the peaks in the SSSF. In panel (b), we
show the typical contour plots of the SSSF
$\sum_{\alpha}S^{\alpha\alpha}(\mathcal{Q})$ in the three different phases. To
keep the symmetry of the Brillouin zone, we use the data from the $C_{6}$
cluster. The data from the cylinder give the same conclusion except that the
symmetry of the pattern is lower. Hence, the characteristic momenta in these
three phases are $\mathcal{Q}_{\Gamma}=\left(0,0\right)$,
$\mathcal{Q}_{{\rm{K}}}=\left({2\pi}/{3},{2\pi}/(3\sqrt{3})\right)$, and
$\mathcal{Q}_{\Gamma^{\prime}}=\left({2\pi}/{3},{2\sqrt{3}\pi}/{3}\right)$,
respectively, as well as others associated by the symmetry transformation of
the Brilouin zone. In panel (c), we show the order parameters at the
characteristic momenta as a function of $J$. The order parameters at the
characteristic momenta in other phases are also shown for comparison. This
figure tells us that the three phases from left to right are an FM phase, a
$120^{\circ}$ phase and an AFM phase. At all the phase transition points, the
discontinuity in the order parameters is clearly visible, indicating a first-
order transition. This is consistent with the conclusions drawn from the
ground-state energy.

Figure 3: Extrapolation of the order parameter at $\mathcal{Q}_{\rm{K}}$ in
$120^{\circ}$ phase. The geometries used are rhombus clusters of 18 and 24
sites, and XC cylinder of $6\times 3\times 2$ and YC cylinder of $6\times
4\times 2$. The dashed line is the linear fitting and the solid line is the
quadratic fitting, and the magnetization is estimated as $M=0.67(2)$.

To check the tendency of the $120^{\circ}$ magnetic ordering against the
system size, we calculate the order parameter of the $120^{\circ}$ phase on
four clusters. Here, the Heisenberg interaction $J$ is taken as $-0.10$, which
is close to the strength of strain $-2.34\%$ XuChangsong2020 . As can be seen
from Fig. 3, the order parameter is rather robust and the magnetization is
around one-half of the spin value $3/2$, albeit with a decreasing trend as the
system size increases. The linear and quadratic fittings suggest that the
magnetization is approximately 0.67(2) in the thermodynamic limit. It is
worthing to note that the magnetization is obvious larger than the well-
recognized $120^{\circ}$ order in the spin-$1/2$ triangular-lattice Heisenberg
antiferromagnet QianLi2020 , demonstrating the robustness of the $120^{\circ}$
order in the present honeycomb model.

##  IV Classical phase diagrams

Large spins usually indicate weak quantum fluctuation. An $S=3/2$ spin is
considered to be near the boundary between the classical and the quantum world
Songvilay2021 . On the other hand, quantum fluctuation can be enhanced by
anisotropic magnetic terms, So to what extent does an $S=3/2$ spin model
behave classically in frustrated systems? We will try to answer this question
by a direct comparison between the quantum phase diagram and its classical
counterpart. For this purpose, in this section we will establish the classical
phase diagram by the parallel-tempering Monte Carlo method Hukushima1996 ;
Mitsutake2000 ; JanssenAV2016  . The simulations are performed on a
12$\times$12$\times$2 rhombic cluster. The spin configurations are updated by
the heat bath algorithm.

Figure 4: (a) The blue open circles are the ground energy obtained by the
Monte Carlo. There are four phases in the phase diagram, i.e., an FM phase, a
zigzag phase, a $120^{\circ}$ phase and an AFM phase from left to right. The
transition points are $J\approx{-0.32063},-0.26074,0.020$. The red lines are
analytical results of the ground-state energy. (b) The SSSFs in the
corresponding phases.

In the classical limit, each spin is regarded as an $O\left(3\right)$ vector
in real space. We use the Monte Carlo method to find the ground state. The
corresponding ground-state energy and spin configurations are then readily
available. In Fig. 4 (a), we plot the classical ground-state energy per spin
$E_{g}/\left(NS^{2}\right)$ as a function of $J$. Interestingly, it is easy to
see that the ground-state energy is a piece-wise function of $J$ and four
phases can be determined in the phase diagram. The spin configurations in each
phase directly tell us the characteristic features, which are an FM phase, a
zigzag (ZZ) phase, a $120^{\circ}$ phase and an AFM phase, respectively, from
left to right. Moreover, the expressions of the classical energy can be given
analytically

$\displaystyle\begin{split}&E_{\rm{FM}}\\!=\\!\frac{S^{2}}{2}\\!\left[3J\\!+\\!K\\!-\\!\Gamma\\!-\\!2\Gamma^{\prime}\\!+\\!P(\Gamma,\\!\Gamma^{\prime},\\!A_{cc})\Theta\\!\left(P(\Gamma,\\!\Gamma^{\prime},\\!A_{cc})\right)\right],\\\
&E_{\rm{ZZ}}=\frac{S^{2}}{4}\left(2J\\!-\\!\Gamma\\!+\\!2\Gamma^{\prime}\\!+\\!2A_{cc}\right)+S^{2}Q\left(K,\Gamma,\Gamma^{\prime},A_{cc}\right),\\\
&E_{120^{\circ}}=\frac{S^{2}}{2}\left(-K-2\Gamma+2\Gamma^{\prime}\right),\\\
&E_{\rm{AFM}}=\frac{S^{2}}{2}\left(-3J-K+\Gamma+2\Gamma^{\prime}\right),\end{split}$

where $\Theta\left(x\right)$ is a step function and
$P\left(\Gamma,\Gamma^{\prime},A_{cc}\right)$ and
$Q\left(K,\Gamma,\Gamma^{\prime},A_{cc}\right)$ have the following forms

$\displaystyle P=$ $\displaystyle 3\Gamma+6\Gamma^{\prime}+2A_{cc},$
$\displaystyle Q=$
$\displaystyle-\frac{1}{12}\sqrt{32\left(K\\!-\\!\Gamma\\!+\\!\Gamma^{\prime}\right)^{2}\\!+\\!\left(6A_{cc}\\!+\\!2K\\!+\\!7\Gamma\\!+\\!2\Gamma^{\prime}\right)^{2}}.$

Here, the occurrence of the step function $\Theta(x)$ in the FM phase is
because that its classical moment direction could either lie in the $ab$ plane
(when $x>0$) or point along the $c$ direction (when $x<0$). In Fig. 4(a), we
also plot these analytical results for comparison, which are marked by red
lines. The overlap of these data demonstrates the reliability of our numerical
results. Just as in the quantum case, the phases can also be identified by the
SSSF, which are plotted in Fig. 4(b). The SSSF in the FM phase, the
$120^{\circ}$ phase and the AFM phase show a close resemblance to the quantum
cases in Fig. 2(b). In particular, the heights of characteristic peaks are
close in their value, indicating that the Eq. (1) in the quantum case is close
to its classical limit.

Figure 5: Classical ground-state energy is plotted as a function of $J$ for
$A_{cc}=0.000,0.100,0.150,0.180$ and $0.222$. We can see that when
$A_{cc}\lesssim 0.150$, zigzag phase disappears.

From the Fig. 2 and Fig. 4, we can see that the transition point from the
$120^{\circ}$ phase to the AFM phase agrees well in both the quantum and
classical cases, which are $J=0.022(2)$ and $J=0.020$, respectively. However,
an unexpected zigzag phase emerges in the classical phase diagram, which sits
between the FM phase and $120^{\circ}$ phase. We have checked that when
$A_{cc}$ is smaller than $0.150$, the zigzag phase disappears as shown in Fig.
5. This is not surprising because the classical spin corresponds to
$S\rightarrow\infty$ and in the $S=3/2$ case the quantum fluctuation can not
be completely neglected and thus is expected to shift the phase boundary.

##  V Finite-size effect in 12-site clusters

Generally speaking, the QSL is known as a nonmagnetic state with
fractionalized excitations, which could be manifested by topological
entanglement entropy KtvPre2006 ; LvnWen2006 , continuum excitation spectrum,
and thermodynamic signatures DoPark2017 . In the Kitaev QSL, the spins
fractionalize into localized and itinerant Majorana fermions, resulting in a
double-peak specific heat anomaly NasuUM2015 ; KogaTN2018 . Nevertheless, the
double-peak structure of specific heat found in the previous work on a 12-site
cluster can not simply be interpreted as a signature of spin fractionalization
XuChangsong2020 . Instead, it may originate from the artificial frustration
due to the incompatible cluster with respect to the underlying magnetic
ordering.

We have illustrated in Sec. III that the intermediate region is a
$120^{\circ}$ magnetically state rather than the QSL. For comparison, we fix
$J$ = $-0.10$ again where the ground state locates deep in the $120^{\circ}$
phase. To be compatible with the six-sublattice $120^{\circ}$ ordering, the
number of lattice sites in a cluster should be a multiple of six. Even though
such condition is satisfied, the cluster shape should be chosen carefully so
as to match the magnetic order. As stated in Ref. [XuChangsong2020, ], the 12A
cluster in Fig. 6(a) is likely used in the specific heat calculation.
Unfortunately, this geometry does not coordinate the $120^{\circ}$ order in
the vertical boundary, which inevitably introduces the artificial frustration
and ruins the magnetic ordering significantly in the small cluster. To confirm
this, we have performed the DMRG calculations on the 12A cluster. The total
ground-state energy is -4.835103988, and the corresponding SSSF is shown in
the first panel of Fig. 6(b). There is not a well-defined peak in the
Brillouin zone, in accord with the proposed QSL inferred from the specific
heat.

By contrast, a slightly modification of the geometry will remove the mismatch
of the spins at boundaries, see 12B cluster in Fig. 6(a). Its SSSF in the
second panel of Fig. 6(b) exhibits a well-defined peak at
$\mathcal{Q}_{{\rm{K}}}$ point, in line with the desirable $120^{\circ}$
order. We have also considered two other 12-site clusters (termed 12C and 12D)
which match the $120^{\circ}$ ordering, and their SSSFs with peaks at
$\mathcal{Q}_{{\rm{K}}}$ point are shown in the last row in Fig. 6(b).
Notably, All the clusters that coordinate the $120^{\circ}$ spin textures
(i.e., 12B, 12C, and 12D) share a same total energy of -5.535684265, which is
considerably lower than that of the 12A cluster. This demonstrates again that
the QSL based on the 12A cluster is not the true ground state. We want to
stress that even in the 12B, 12C and 12D clusters which have the same ground-
state energy the difference in their SSSFs is clearly visible. In particular,
the SSSFs in 12C and 12D are rather diffusive in the momentum space in
comparison with that of 12A case, demonstrating that the results on large
clusters are necessary to reveal the underlying phase.

Figure 6: (a) Four different 12-site clusters with full periodic boundary
conditions. The background shows the classical spin configuration of the six-
sublattice $120^{\circ}$ ordering. The 12A cluster have a energy of
-4.835103988, while the last three clusters (12B, 12C, and 12D) have a same
but lower total energy of -5.535684265. (b) SSSF on the four clusters, which
are 12A (top left), 12B (top right), 12C (bottom left) and 12D (bottom right).

##  VI summary and outlook

In this work, we study an $S=3/2$ $\rm JK\Gamma\Gamma^{\prime}$+SIA model
proposed to describe $\rm{CrSiTe_{3}}$ by using the DMRG method. Upon
increasing the strain, the Heisenberg ($J$) interaction varies significantly
while the remaining ones almost keep unchanged. we thus map out the
$J$-dependent phase diagram as a mimicry of the strain process. When $J$
varies from negative to positive, three magnetically ordered phases are found,
including an FM phase, a $120^{\circ}$ phase and an AFM phase. All the phase
transitions are of the first order, characterized by the kinks in the ground-
state energy and discontinuities in the order parameters. Experientially, the
quantum fluctuation is expected to be strongly suppressed in the models with a
large spin of $S=3/2$, and the physics should bear resemblance to its
classical counterpartLee2020 . We have strengthened this conclusion by
performing the classical Monte Carlo simulations on the same model, with the
exception that a zigzag phase is found between the FM phase and the
$120^{\circ}$ phase. The zigzag phase is sensitive to the SIA, and its
territory shrinks gradually and it disappears eventually as SIA decreases.

It was argued in a previous work XuChangsong2020  that there is a possibly
QSL, which is inferred from the double-peak structure of specific heat and the
related plateau of thermal entropy on a 12-site cluster. Curiously, the
proposed QSL region happens to coincide with the $120^{\circ}$ phase
identified in our work. In fact, the double-peak specific heat anomaly solely
could not legitimate the Kitaev QSL. To begin with, the low-temperature peak
found in Ref. [XuChangsong2020, ] is rather weak when compared to the high-
temperature analogy, which is in contrast to the current wisdom that the
magnitudes of the two peaks are comparable in the Kitaev QSL NasuUM2015 ;
KogaTN2018 . In addition, it is unclear if such a low-temperature peak will
still persist with the increasing of the system size. Even though the double-
peak structure indeed exists in the ground state, it is still not enough to
authenticate a QSL as magnetically ordered states could also exhibit such an
abnormal phenomenon YamajiSY2016 . In this regard, it is imperative to check
if the ground state is disordered or not. Our large-scale DMRG calculation,
however, suggests that the ground state is a magnetically ordered
$120^{\circ}$ state with a large order parameter.

Although our work implies that QSL is unfavorable in $\rm{CrSiTe_{3}}$, three
magnetic ordered phases are found under the strain. These ordered states may
persist to finite temperature because the Mermin-Wagner theorem does not
forbid such orders in the model (1). In this sense, this material offers a
playground to study these magnetic orders and phase transitions in low
dimensionsLiXuetal2021 . Recent experiment and ab initio calculations support
$\rm{Cr}$-based ferromagnetic van der Waals materials are $S=3/2$ Kitaev
materials XuChangsong2018 ; XuChangsong2020 ; Lee2020 ; Bacaksiz2021 , the
model here just takes the nearest-neighbor interactions into account, the
other possible interactions are not considered. Due to the complexity of real
material, the neutron spectra of $\rm{CrI_{3}}$ can be fitted by different
models ChenLebing2020 . In $\rm{CrSiTe_{3}}$, we can not exclude other
possible terms such as the next-nearest neighbors in the Hamiltonian.
Moreover, the model Hamiltonian (1) and corresponding parameters remain to be
verified by experiments. We expect our results can be verified in future
experiments and be helpful to understand the magnetic behavior in two-
dimensional ferromagnetism in high-spin Kitaev materials.

###### Acknowledgements.

We thank Y. Qi, F. Zheng, C. Xu and L. Bellaiche for helpful discussion. This
work is supported by the National Natural Science Foundation of China (Grant
No. 11874188, 12047501, 11774300). The computational resource is partly
supported by Supercomputing Center of Lanzhou University.

## References

  * (1) A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006).
  * (2) P. W. Anderson, Resonating valence bonds: A new kind of insulator?, Mater. Res. Bull. 8, 153 (1973).
  * (3) L. Balents, Spin liquids in frustrated magnets, Nature (London) 464, 199 (2010).
  * (4) J. Wen, S.-L. Yu, S. Li, W. Yu, and Jian-Xin Li, Experimental identification of quantum spin liquids, npj Quantum Mater. 4, 12 (2019).
  * (5) G. Jackeli and G. Khaliullin, Mott insulators in the strong spin-orbit coupling limit: From Heisenberg to a quantum compassand Kitaev models, Phys. Rev. Lett. 102, 017205 (2009).
  * (6) J. Chaloupka, G. Jackeli, and G. Khaliullin, Kitaev-Heisenberg model on a honeycomb lattice: Possible exotic phases in iridium oxides ${\rm A_{2}IrO_{3}}$, Phys. Rev. Lett. 105, 027204 (2010).
  * (7) K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, $\alpha$–${\rm RuCl}_{3}$: A spin-orbit assisted mott insulator on a honeycomb lattice, Phys. Rev. B 90, 041112 (2014).
  * (8) H.-S. Kim, V. S. V., A. Catuneanu, and H.-Y. Kee, Kitaev magnetism in honeycomb $\alpha$–$\mathrm{RuCl_{3}}$ with intermediate spin-orbit coupling, Phys. Rev. B 91, 241110 (2015).
  * (9) A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li, M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Proximate Kitaev quantum spin liquid behaviour in a honeycomb magnet, Nat. Mater. 15, 733 (2016).
  * (10) Wei Wang, Zhao-Yang Dong, Shun-Li Yu, and Jian-Xin Li, Theoretical investigation of magnetic dynamics in $\rm\alpha-RuCl_{3}$, Phys. Rev. B 96, 115103 (2017).
  * (11) S. M. Winter, K. Riedl, P. A. Maksimov, A. L. Chernyshev, A. Honecker, and R. Valenti, Breakdown of magnons in a strongly spin-orbital coupled magnet, Nat. Commun. 8, 1152 (2017).
  * (12) P. A. Maksimov and A. L. Chernyshev, Rethinking $\alpha$–$\mathrm{RuCl_{3}}$, Phys. Rev. Research 2, 033011 (2020).
  * (13) H. Li, H.-K. Zhang, J. Wang, H.-Q. Wu, Y. Gao, D.-W. Qu, Z.-X. Liu, S.-S. Gong, and W. Li, Identification of Magnetic Interactions and High-field Quantum Spin Liquid in $\alpha$–$\mathrm{RuCl_{3}}$, Nat. Commun. 12, 4007 (2021).
  * (14) M. Hermanns, I. Kimchi, and J. Knolle, Physics of the Kitaev model: Fractionalization, dynamic correlations, and material connections, Annu. Rev. Condens. Matter Phys. 9, 17 (2018).
  * (15) H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and realization of Kitaev quantum spin liquids, Nat. Rev. Phys. 1, 264 (2019).
  * (16) J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Generic spin model for the honeycomb iridates beyond the Kitaev limit, Phys. Rev. Lett. 112, 077204 (2014).
  * (17) J. G. Rau and H.-Y. Kee, Trigonal distortion in the honeycomb iridates: Proximity of zigzag and spiral phases in ${\rm Na_{2}IrO_{3}}$, arXiv:1408.4811 (2014).
  * (18) J. Wang, B. Normand, and Z.-X. Liu, One Proximate Kitaev Spin Liquid in the $K$-$J$-$\Gamma$ Model on the Honeycomb Lattice, Phys. Rev. Lett. 123, 197201 (2019).
  * (19) Q. Luo, J. Zhao, H.-Y. Kee and X. Wang, Gapless quantum spin liquid in a honeycomb $\Gamma$ magnet, npj Quantum Mater. 6, 57 (2021).
  * (20) A. Catuneanu, Y. Yamaji, G. Wachtel, Y.-B. Kim, and H.-Y. Kee, Path to stable quantum spin liquids in spin-orbit coupled correlated materials, npj Quantum Mater. 3, 23 (2018).
  * (21) M. Gohlke, G. Wachtel, Y. Yamaji, F. Pollmann, and Y. B. Kim, Quantum spin liquid signatures in Kitaev-like frustrated magnets, Phys. Rev. B 97, 075126 (2018).
  * (22) Q. Luo, P. P. Stavropoulos, and H.-Y. Kee, Spontaneous Chiral-Spin Ordering in Spin-Orbit Coupled Honeycomb Magnets, arXiv:2010.11233 (2020).
  * (23) R. Sano, Y. Kato, and Y. Motome, Kitaev-Heisenberg hamiltonian for high-spin $d^{7}$ mott insulators, Phys. Rev. B 97, 014408 (2018).
  * (24) H. Liu and G. Khaliullin, Pseudospin exchange interactions in $d^{7}$ cobalt compounds: Possible realization of the Kitaev model, Phys. Rev. B 97, 014407 (2018).
  * (25) H. K. Vivanco, B. A. Trump, C. M. Brown, and T. M. McQueen, Competing antiferromagnetic-ferromagnetic states in a $d^{7}$ Kitaev honeycomb magnet, Phys. Rev. B 102, 224411 (2020).
  * (26) M. Songvilay, J. Robert, S. Petit, J. A. Rodriguez-Rivera, W. D. Ratcliff, F. Damay, V. Balédent, M. Jiménez-Ruiz, P. Lejay, E. Pachoud, A. Hadj-Azzem, V. Simonet, and C. Stock, Kitaev interactions in the Co honeycomb antiferromagnets $\mathrm{Na_{3}Co_{2}SbO_{6}}$ and $\mathrm{Na_{2}Co_{2}TeO_{6}}$, Phys. Rev. B 102, 224429 (2020).
  * (27) S.-H. Jang, R. Sano, Y. Kato, and Y. Motome, Antiferromagnetic Kitaev interaction in $f$-electron based honeycomb magnets, Phys. Rev. B 99, 241106 (2019).
  * (28) S.-H. Jang, R. Sano, Y. Kato, and Y. Motome, Computational design of $f$-electron Kitaev magnets: Honeycomb and hyperhoneycomb compounds $A\mathrm{{}_{2}PrO_{3}}$ ($A$ = alkali metals), Phys. Rev. Materials 4, 104420 (2020).
  * (29) Y. Motome, R. Sano, S. Jang, Y. Sugita, and Y. Kato, Materials design of Kitaev spin liquids beyond the Jackeli-Khaliullin mechanism, J. Phys. Condens. Matter 32, 404001 (2020).
  * (30) P. P. Stavropoulos, D. Pereira, and H.-Y. Kee, Microscopic mechanism for a higher-spin Kitaev model, Phys. Rev. Lett. 123, 037203 (2019).
  * (31) P. P. Stavropoulos, X. Liu, and H.-Y. Kee, Magnetic anisotropy in spin-3/2 with heavy ligand in honeycomb mott insulators: Application to $\mathrm{CrI_{3}}$, Phys. Rev. Research 3, 013216 (2021).
  * (32) E. A. Zvereva, M. I. Stratan, Y. A. Ovchenkov, V. B. Nalbandyan, J.-Y. Lin, E. L. Vavilova, M. F. Iakovleva, M. Abdel-Hafiez, A. V. Silhanek, X.-J. Chen, A. Stroppa, S. Picozzi, H. O. Jeschke, R. Valentí, and A. N. Vasiliev, Zigzag antiferromagnetic quantum ground state in monoclinic honeycomb lattice antimonates $A\mathrm{{}_{3}Ni_{2}SbO_{6}(}A\mathrm{=Li,Na)}$, Phys. Rev. B 92, 144401 (2015).
  * (33) A. I. Kurbakov, A. N. Korshunov, S. Y. Podchezertsev, A. L. Malyshev, M. A. Evstigneeva, F. Damay, J. Park, C. Koo, R. Klingeler, E. A. Zvereva, and V. B. Nalbandyan, Zigzag spin structure in layered honeycomb $\mathrm{Li_{3}Ni_{2}SbO_{6}}$: A combined diffraction and antiferromagnetic resonance study, Phys. Rev. B 96, 024417 (2017).
  * (34) C. Xu, J. Feng, H. Xiang, and L. Bellaiche, Interplay between Kitaev interaction and single ion anisotropy in ferromagnetic $\mathrm{CrI_{3}}$ and $\mathrm{CrGeTe_{3}}$ monolayers, npj Comput. Mater. 4, 57 (2018).
  * (35) C. Xu, J. Feng, M. Kawamura, Y. Yamaji, Y. Nahas, S. Prokhorenko, Y. Qi, H. Xiang, and L. Bellaiche, Possible Kitaev quantum spin liquid state in 2D materials with $S$ = 3/2, Phys. Rev. Lett. 124, 087205 (2020).
  * (36) I. Lee, F. G. Utermohlen, D. Weber, K. Hwang, C. Zhang, J. van Tol, J. E. Goldberger, N. Trivedi, and P. C. Hammel, Fundamental spin interactions underlying the magnetic anisotropy in the kitaev ferromagnet $\mathrm{CrI_{3}}$, Phys. Rev. Lett. 124, 017201 (2020).
  * (37) L. Chen, J.-H. Chung, T. Chen, C. Duan, A. Schneidewind, I. Radelytskyi, D. J. Voneshen, R. A. Ewings, M. B. Stone, A. I. Kolesnikov, B. Winn, S. Chi, R. A. Mole, D. H. Yu, B. Gao, and P. Dai, Magnetic anisotropy in ferromagnetic $\mathrm{CrI_{3}}$, Phys. Rev. B 101, 134418 (2020).
  * (38) C. Bacaksiz, D. Šabani, R. M. Menezes, and M. V. Milošević, Distinctive magnetic properties of ${\rm CrI_{3}}$ and ${\rm{CrBr_{3}}}$ monolayers caused by spin-orbit coupling, Phys. Rev. B 103, 125418 (2021).
  * (39) J. L. Lado and J. Fernández-Rossier, On the origin of magnetic anisotropy in two dimensional $\mathrm{CrI_{3}}$, 2D Mater. 4, 035002 (2017).
  * (40) J. Nasu, M. Udagawa, and Y. Motome, Thermal fractionalization of quantum spins in a Kitaev model: Temperature-linear specific heat and coherent transport of Majorana fermions, Phys. Rev. B 92, 115122 (2015).
  * (41) A. Koga, H. Tomishige, and J. Nasu, Ground-state and Thermodynamic Properties of an $S=1$ Kitaev Model, J. Phys. Soc. Jpn. 87, 063703 (2018).
  * (42) J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moessner, Dynamics of a Two-Dimensional Quantum Spin Liquid: Signatures of Emergent Majorana Fermions and Fluxes, Phys. Rev. Lett. 112, 207203 (2014).
  * (43) Q. Luo, S. Hu, and H.-Y. Kee, Unusual excitations and double-peak specific heat in a bond-alternating spin-1 $K$-$\Gamma$ chain, arXiv:2104.10725 (2021) [Phys. Rev. Research (to be published)].
  * (44) Y. Yamaji, T. Suzuki, T. Yamada, S.-i. Suga, N. Kawashima, and M. Imada, Clues and criteria for designing Kitaev spin liquid revealed by thermal and spin excitations of honeycomb iridates Na2IrO3, Phys. Rev. B 93, 174425 (2016).
  * (45) S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * (46) I. Peschel, X. Q. Wang, M. Kaulke, and K. Hallberg, Density-matrix renormalization (Springer, Berlin, 1999).
  * (47) U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * (48) Qian Li, Hong Li, Jize Zhao, Hong-Gang Luo, Z. Y. Xie, Magnetization of the spin-1/2 Heisenberg antiferromagnet on the triangular lattice, arXiv:2009.03765 (2020).
  * (49) M. Songvilay, S. Petit, F. Damay, G. Roux, N. Qureshi, H. C. Walker, J. A. Rodriguez-Rivera, B. Gao, S.-W. Cheong, and C. Stock, From one to two-magnon excitations in the $S$ = 3/2 magnet $\beta$–$\mathrm{CaCr_{2}O_{4}}$, Phys. Rev. Lett. 126, 017201 (2021).
  * (50) K. Hukushima and K. Nemoto, Exchange Monte Carlo method and application to spin glass simulations, J. Phys. Soc. Jpn. 65, 1604 (1996).
  * (51) A. Mitsutake and Y. Okamoto, Replica-exchange simulated tempering method for simulations of frustrated systems, Chem. Phys. Lett. 332, 131 (2000).
  * (52) L. Janssen, E. C. Andrade, and M. Vojta, Honeycomb-Lattice Heisenberg-Kitaev Model in a Magnetic Field: Spin Canting, Metamagnetism, and Vortex Crystals. Phys. Rev. Lett. 117, 277202 (2016).
  * (53) A. Kitaev and J. Preskill, Topological Entanglement Entropy, Phys. Rev. Lett. 96, 110404 (2006).
  * (54) Michael Levin and Xiao-Gang Wen, Detecting Topological Order in a Ground State Wave Function, Phys. Rev. Lett. 96, 110405 (2006).
  * (55) Seung-Hwan Do, Sang-Youn Park, Junki Yoshitake, Joji Nasu, Yukitoshi Motome, Yong Seung Kwon, D. T. Adroja, D. J. Voneshen, Kyoo Kim, T.-H. Jang, J.-H. Park, Kwang-Yong Choi, and Sungdae Ji, Majorana fermions in the Kitaev quantum spin system $\alpha$-RuCl3, Nat. Phys. 13, 1079 (2017).
  * (56) Zefang Li, Dong-Hong Xu, Xue Li, Hai-Jun Liao, Xuekui Xi, Yi-Cong Yu, Wenhong Wang, Abnormal critical fluctuations revealed by magnetic resonance in the two-dimensional ferromagnetic insulators, arXiv:2101.02440.