# Detecting Atmospheric Molecules of Temperate Terrestrial Exoplanets using
High-Resolution Spectroscopy in the Mid Infrared Domain

Yuka Fujii National Astronomical Observatory of Japan  Earth-Life Science
Institute, Tokyo Institute of Technology  Taro Matsuo  Nagoya University

(Accepted January 28, 2021)

###### Abstract

Motivated by the development of high-dispersion spectrographs in the mid-
infrared (MIR) range, we study their application to the atmospheric
characterization of nearby non-transiting temperate terrestrial planets around
M-type stars. We examine the detectability of CO2, H2O, N2O, and O3 in high-
resolution planetary thermal emission spectra at 12-18 $\mu$m assuming an
Earth-like profile and a simplified thermal structure. The molecular line
width of such planets can be comparable to or broader than the Doppler shift
due to the planetary orbital motion. Given the likely difficulty in knowing
the high-resolution MIR spectrum of the host star with sufficient accuracy, we
propose to observe the target system at two quadrature phases and extract the
differential spectra as the planetary signal. In this case, the signals can be
substantially suppressed compared with the case where the host star spectrum
is perfectly known, as some parts of the spectral features do not remain in
the differential spectra. Despite this self-subtraction, the CO2 and H2O
features of nearby ($\lesssim$ 5 pc) systems with mid-/late-M host stars would
be practical with a 6.5-meter-class cryogenic telescope, and orbital
inclination could also be constrained for some of them. For CO2 and N2O in a 1
bar Earth-like atmosphere, this method would be sensitive when the mixing
ratio is 1-103 ppm. The detectability of molecules except O3 is not
significantly improved when the spectral resolution is higher than
$\mathcal{R}\gtrsim 10,000$, although the constraint on the orbital
inclination is improved. This study provides some benchmark cases useful for
assessing the value of MIR high-resolution spectroscopy in terms of
characterization of potentially habitable planets.

planets and satellites: atmospheres — planets and satellites: terrestrial
planets — astrobiology

##  1 Introduction

Spectroscopy of Earth-sized planets around the so-called habitable zones (HZs;
e.g., Kasting et al., 1993) is one of the key near-future targets in
astronomical observations. Several approaches have been put forward so far.
The first major approach is direct (high-contrast) imaging. Since the late
20th century, space-based direct imaging possibilities—with coronagraphs in
the visible and near-infrared domains, or using nulling interferometers in the
mid-infrared (MIR) domain—have been examined and were encapsulated in the
proposed missions like TPF-C, TPF-I, and Darwin (e.g., Levine et al., 2009;
Beichman et al., 1999; Léger et al., 1996). These studies are inherited by the
ongoing discussions for the future space missions for direct imaging of
scattered light of Earth-sized planets (e.g., Starshade with Roman Space
Telescope, LUVOIR, HabEx, LIFE).

In the meanwhile, the successful atmospheric characterization of Hot Jupiters
through transmission spectroscopy, eclipse spectroscopy, and detection of
phase variations has encouraged its application to smaller Earth-sized planets
around late-type stars—this is the second appraoch (see a review by Deming &
Seager, 2017, and references therein). It has actually been applied to the
Hubble Space Telescope (HST) observations of TRAPPIST-1 planets (Gillon et
al., 2016, 2017), the Earth-sized transiting planets around the habitable zone
of an M8-type star, while only the upper limit of atmospheric features have
been obtained (de Wit et al., 2018; Zhang et al., 2018). The upcoming James
Webb Space Telescope (JWST) will be used for more sensitive transit
spectroscopy of smaller planets in a wider wavelength range (Beichman et al.,
2014; Greene et al., 2016). The high-dispersion spectrographs on the next-
generation ground-based telescopes (i.e., extremely large telescopes) may also
be used for transit transmission spectroscopy to identify the atmospheric
features of Earth-sized exoplanets buried in a number of telluric absorption
lines (e.g., Snellen et al., 2013).

The third approach to study potentially habitable exoplanets is to combine
high-contrast imaging and high-resolution spectroscopy with these next-
generation ground-based telescopes (e.g., Snellen et al., 2015). Such
observations using ground-based telescopes are planned mainly in the near-
infrared wavelengths, so as to avoid thermal background that is inevitable for
ground-based observations at longer wavelengths, while achieving the spatial
resolution necessary for high-contrast imaging.

Lastly, the atmospheric features of potentially habitable planets may be
identified in the combined spectrum of the star and the planet in the MIR
range, without a planetary transit or a starlight-suppressing instrument
(e.g., Kreidberg & Loeb, 2016; Snellen et al., 2017). Indeed, the planet-to-
star flux ratio of potentially habitable planets is significantly improved in
the MIR domain ($\sim$a few tens of parts-per-million or larger) compared with
that in the shorter visible domain ($\sim 10^{-7}$ or less) as shown in the
upper panel of Figure 1, and it is not unrealistic to detect planetary signal
in the combined spectra. Approximately estimating the wavelength-dependence of
the signal-to-noise ratio for the planet signal per wavelength resolution
element by
$\mathcal{N}_{p}\Delta\lambda/\sqrt{\mathcal{N_{\star}}\Delta\lambda}=\sqrt{\mathcal{N}_{p}C\Delta\lambda}=\sqrt{\mathcal{N}_{p}C\lambda/\mathcal{R}}$
($\mathcal{N}_{p}$ and $\mathcal{N}_{\star}$ are the photon spectrum of the
planet and the star, respectively, $\Delta\lambda$ is the wavelength width of
the resolution element, $C$ is the planet-to-star flux ratio, i.e.,
$C\equiv\mathcal{N}_{p}/\mathcal{N}_{\star}$, and $\mathcal{R}$ is the fixed
spectral resolution), the MIR domain, specifically around 13-100 $\mu$m, is
most useful (the lower panel of Figure 1). The unique advantages of this MIR
approach are that (1) unlike high-contrast technique it does not require
specialized instruments to occult stellar flux, and that (2) unlike
transmission or eclipse spectroscopy it can be applied to both transiting and
non-transiting planets. The latter is critical to increase the number of
targets, as the transit probability of habitable-zone planets is up to about
5%.

Figure 1: Upper panel: A schematic figure of the photon count per wavelength
resolution element assuming a black body planet spectrum with 288 K (i.e.
$\mathcal{N}_{\rm p}\Delta\lambda$; red dashed) and the planet-to-star flux
ratio assuming a black body spectrum mimicking a mid-M host star (i.e. $C$;
blue solid). The wavelength dependence of the planet-to-star flux ratio does
not depend on the spectral types of the star as long as the stellar spectrum
can be approximated by Rayleigh-Jeans law. Lower panel: Signal-to-noise ratio
per wavelength resolution element as a function of wavelength, estimated by
the product of the planet photon count and the planet-to-star flux ratio.

Along this line, Kreidberg & Loeb (2016) proposed low-resolution spectroscopy
of Proxima Centauri systems to try to detect 9.6 $\mu$m O3 features originated
from the atmosphere of Proxima Centauri b, the nearest non-transiting
potentially habitable planets (Anglada-Escudé et al., 2016).

Snellen et al. (2017) proposed that medium-resolution spectroscopy (MRS) in
the MIR range can be used to identify high-frequency features due to planetary
atmospheric molecules in the combined spectra, and they estimated that CO2
features of Proxima Centauri b would be detected after 5 days of observations
with MRS mode of JWST/MIRI. Although these are encouraging possibilities, such
observations rely on the precise knowledge of the stellar spectrum as well as
the sensitivity of the detector elements.

While the above studies have in mind the upcoming JWST that have only low-
resolution and medium-resolution capabilities in the MIR, technologies for
high-dispersion spectrograph in the MIR have been recently developed for
future cryogenic space telescopes, including SPICA (e.g., Sarugaku et al.,
2012). The technologies are expected to largely reduce the size of the MIR
spectrograph and enhance its throughput. In the further future, Origins Space
Telescope (OST) is also contemplating high-resolution spectrograph with a
larger aperture (Sakon et al., 2018). With high-resolution spectroscopy
($\mathcal{R}\sim 30,000$), not only the fine structures of absorption bands
of the planetary spectrum can be identified, but also the Doppler shift of a
HZ planet can be resolved (unless its orbit is close to face-on). Thanks to
the Doppler-shift, it is possible to identify planetary features even without
the assumption for the stellar spectrum, and in theory the orbital inclination
can also be constrained. Although these new opportunities will be useful for
characterizing atmospheres of potentially habitable planets, its potential
remains largely unexplored yet.

In this paper, we examine the utility of High-Resolution Spectroscopy (HRS) in
the MIR for characterizing atmospheres of nearby potentially habitable
planets. In Section 2, we begin by demonstrating the MIR spectral features of
potentially habitable planets using a simplified 1-dimensional atmospheric
model. We then move to detectability assessment in Section 3. We propose to
observe the target system at two quadrature phases, and discuss how this
method could constrain the presence of molecular features and the orbital
inclination under the realistic assumptions on the observational noise,
highlighting the difference between the two analysis strategies that reflect
the possible uncertainty in the prior knowledge of the stellar spectrum. The
dependence on the spectral resolution, bandpass, and the abundance of
molecules are examined in Section 4. Section 5 discusses the range of the
targets and the synergies with other techniques to characterize potentially
habitable planets (Section 5.1), the effects of other noise (Section 5.2), the
effect of different thermal profiles of the planet (§5.3), and the comparison
with the previously studied medium-resolution spectroscopy (§5.4). Finally,
Section 6 summarizes our results.

##  2 Atmospheric features in MIR high-resolution spectra

In this section, we model high-resolution planetary spectra of an Earth-like
planet and discuss the characteristics of the spectral features of molecules
of interest. These modeled spectra will be used as input for the detectability
analyses in Section 3.

###  2.1 Model Atmosphere

Figure 2: The assumed vertical profiles of the mixing ratios of the molecules
considered (left) and temperature (right). The mixing ratio of molecules are
based on “US standard” model. The temperature profile is determined by the dry
adiabatic lapse rate (9.8 K/km) in the lower atmospheres below 0.1 bar, above
which an isothermal profile is assumed.

Our model atmosphere is based on Earth’s atmosphere except for the temperature
profile. We consider four molecules that are present in Earth’s atmosphere and
are radiatively active in the mid infrared range: CO2, H2O, N2O, and O3. The
vertical profiles of the mixing ratios of these molecules are taken from the
“US standard” model and are shown in the left panel of Figure 2. The surface
temperature is set at 288 K, again referring to the “US standard” model. The
vertical temperature profile is, however, replaced by a simplified one
comprised of a troposphere with a constant lapse rate ($\Gamma=g/C_{p}$ where
$g$ is the gravity and $C_{p}$ is the specific heat capacity) and an
isothermal stratosphere. The removal of stratospheric thermal inversion is
motivated by the fact that O3, even if it exists, does not lead to a strong
thermal inversion under the irradiation of M-type stars due to the reduced
near-UV flux. The effects of the atmospheric profile on the detectability of
molecules are discussed in Section 5.3.1. For simplicity, the effects of
clouds are ignored.

The surface pressure is fixed at 1 bar, and the tropopause is set at 0.1 bar.
This tropopause pressure is consistent with the observations of Solar system
planets (e.g., Robinson & Catling, 2014) and 3D climate simulations for
habitable planets around M-type stars (e.g., Fujii et al., 2017).

description | symbol | value
---|---|---
surface pressure | $P_{\rm surf}$ | 1 bar
surface temperature | $T_{\rm surf}$ | 288 K
tropopause | $P_{\rm tp}$ | 0.1 bar
surface gravity | $g$ | 9.8 m/s2
tropospheric lapse rate | $\Gamma(=g/C_{p})$ | 9.8 K/km
Table 1: Assumptions for atmospheric profiles.

Given a vertical profile, the top-of-atmosphere outgoing radiance at
wavelength $\lambda$ and the cosine of the zenith angle $\mu$,
$L(\lambda,\mu)$, is computed by

$\displaystyle L(\lambda,\mu)$ $\displaystyle=$ $\displaystyle\int
B(T;\lambda)\exp(-\tau/\mu)d\tau/\mu$ (1) $\displaystyle\tau$
$\displaystyle\equiv$ $\displaystyle\int_{z}^{\infty}k[T(z),P(z)]x(z)n(z)dz$
(2) $\displaystyle=$
$\displaystyle\int_{0}^{p}k[T(P),P]\,x(P)\frac{dP}{\mu_{\rm atm}g}$ (3)

assuming a no-scattering atmosphere.

The cross sections of molecules are based on HITRAN2016 (Gordon et al., 2017)
and the lines are broadened by the Voigt functions using the algorithm of
Zaghloul & Ali (2011) with the partition functions adopted from HAPI program
(Kochanov et al., 2016). We impose the cut-off of the line wings at 100 cm-1
apart from the line centers. We do not include the continuum absorption (e.g.,
H2O continuum), as it does not significantly affect the detectability of high-
resolution features.

For the efficient evaluation of equation (1) at a large number of the
wavelength grid points, we employed GPU computation through PyCUDA. The
integral in equation (3) is performed with 50 equi-distributed points in $\log
P$ space and that in equation (1) with 50 equi-distributed points in
$\log\tau$ space (from $10^{-2}$ to $10^{2}$), using the trapezoidal rule.

Thermal emission spectrum is calculated with the opacities of all four
molecules included. In addition, we also calculate the spectra with the
opacity of only one molecule turned on, in order to discuss the spectral
characteristics and detectability of individual molecules in isolation.

To obtain the total thermal emission from the planet, equation (1) needs to be
integrated over the planetary disk. With several trials, we find that the
disk-averaged radiance that integrates equation (1) over different $\mu$ with
the weight of the projected area is close to the radiance at $\mu=0.6$. Thus,
we represent the disk-averaged radiance by the radiance with $\mu=0.6$.
Doppler broadening due to planet rotation is not included; HZ planets around
M-type stars are likely to be tidally locked, which means the rotational
velocity of $\sim 10$–100 m/s, and the corresponding to the wavelength shift,
$\sim 10^{-6}\;\mu$m, is sufficiently smaller than the wavelength resolution
elements considered in this paper.

###  2.2 Model Spectra and Their Characteristics

Figure 3: High-resolution ($\mathcal{R}=30,000$) features of molecules in
thermal emission spectra assuming the vertical temperature profile and
molecular mixing ratio shown in Figure 2. In the uppermost panel the opacities
of all molecules are included, while in the lower panels the opacity of one
molecule is assumed at a time. The left panels show the broadband spectra
while the right panels show the zoom-in where the horizontal bars
approximately indicate the maximum range of the Doppler shift (corresponding
to $\sim$100 km/s).

Figure 3 presents the major features in high-resolution spectra
($\mathcal{R}=30,000$) assuming the temperature profile shown in Figure 2. The
uppermost panel includes the opacities of all four molecules, while the lower
panels turn off the opacities of molecules other than the indicated one,
showing more clearly the spectral features of individual molecules.

CO2 has the prominent rotational-vibrational features around 15 $\mu$m
associated with the fundamental bending mode ($\nu_{2}$). The strongest lines
on the P- and R-branches are separated by $\sim 0.04~{}\mu$m, corresponding to
the rotational energy $\sim 1.5$ cm-1, well resolved with
$\mathcal{R}=30,000$. On the other hand, H2O has rotational lines that broadly
spread in the MIR range. The 17 $\mu$m N2O band is mainly due to the bending
mode ($\nu_{2}$) and has a structure similar to the 15 $\mu$m CO2 band. Among
these, CO2 largely contribute to shaping the overall shape of the all-included
spectrum, although the minor features due to H2O and N2O can also be seen.

Compared to these bands, the famous 9.7 $\mu$m O3 band with the overlapping
fundamental vibration modes (symmetric $\nu_{1}$ and asymmetric $\nu_{3}$) is
densely populated with lines. A weaker O3 band exists around 14.5 $\mu$m,
corresponding to the bending mode, which also features clustered lines in a
narrow bandpass. The latter is largely masked by CO2 in Earth’s thermal
emission spectra. We find that this band start to kick in when CO2 is smaller
than $\sim 1$ ppm. Do such CO2-poor atmosphere exist? It is possible, as CO2
abundance in general depends on the carbon cycles, as well as how much the
planet acquires and retains carbon. Planets which develops higher weathering
rate may turn into CO2-poor worlds (e.g., Nakayama et al., 2019). In the
absence of CO2 features, 14.5 $\mu$m O3 band may be more useful than the 9.7
$\mu$m band, due to the better intrinsic detectability (Fig. 1); see more
discussions in Section 4.2.

Here we highlight a characteristics of these high-resolution features—the
relative width of the lines. As indicated in the right panel of Figure 3, the
broad lines extend beyond the wavelength resolution (5$\times 10^{-4}$ $\mu$m
at 15 $\mu$m) and the width of the Doppler shift ($\sim$0.005 $\mu$m at 15
$\mu$m; horizontal bars). This is partly due to the intrinsic broadening of
the lines. The half width of the Lorentz profile (collisional broadening),
which is important at pressures higher than $\sim 0.01$ bar, is approximately
constant in wavenumber for a given set of pressure and temperature, i.e., the
width relative to the wavelength is wider at longer wavelengths. At 1 bar and
288 K, it is approximately 0.1 cm-1, corresponding to 0.002 $\mu$m. The line
shape also depends on the abundance of the molecule and the vertical
temperature gradient. Lines become broad if the opacity is so large that the
atmosphere becomes optically thick at the far wings; this happens for CO2 and
H2O with our model atmosphere. The broad lines have notable influence on the
analysis, as we will see in the next section.

##  3 Detectability

In this section, we examine the detectability of MIR molecular features of
temperate Earth-sized planets presented in Section 2, assuming a high-
dispersion spectrograph mounted on a cryogenic telescope. We try to identify
molecular features in the combined spectrum of the host star and the planet.
The assumptions for mock observations and noise estimate are given in Section
3.1, which is followed by the description of our analysis procedures in
Section 3.2. The resultant constraints on the contrast and the orbital
inclination are discussed in Section 3.3.

###  3.1 Mock observation

####  3.1.1 Targets

Table 2: Assumptions for the host star and the planetary orbit. description | M5 star | M8 star
---|---|---
star radius ($R_{\star}$) | 0.14 $R_{\odot}$ | 0.10$R_{\odot}$
star temperature ($T_{\star}$) | 3000 K | 2500 K
planet/star flux ratioaaEvaluated at 15$\mu$m, assuming a black body spectrum with 288 K for the planet spectrum.  |  $\sim 70$ ppm |  $\sim 200$ ppm
planet orbital radius | 0.0485 au | 0.0146 au
planet orbital period | 11.26 days | 2.27 days
planet orbital velocity | 46.83 km/s | 69.70 km/s

The prime target of this study is temperate rocky planets around M-type stars,
because those around earlier-type stars have too large planet-to-star flux
ratio for planetary signals to be detectable in a reasonable amount of time.
We initially considered three types of host star: early-, mid- and late-M
stars. However, the typical planet-to-star flux ratio with early-M stars is
smaller than 10 ppm, making it extremely challenging to detect planetary
features. Therefore, we focus on mid-M and late-M stars.

Our fiducial model assumes an M5 or M8 star at 5 parsecs. The assumed host
star and planetary parameters are summarized in Table 2. Additionally, the
case study for Proxima Centauri b is also presented, where the host star
properties are the same as those of M5 star in Table 2, the distance is set to
1.3 parsecs, and the planetary radius is assumed to be 1.1 Earth radius (e.g.,
Snellen et al., 2017).

The inclination of the planetary orbit, which controls the amplitude of the
planetary line-of-sight velocity, is fixed at $60^{\circ}$.

The host star spectra are taken from the BT-Settl model (Allard et al., 2012)
with the corresponding effective temperatures, while assuming $\log g=5.0$ and
[Fe/H]$=0.0$. This is an update from the previous studies on the detectability
of spectral features of Proxima Centauri b where the black-body spectrum is
assumed (Snellen et al., 2017; Kreidberg & Loeb, 2016).

####  3.1.2 Configuration

Figure 4 illustrates the observational configuration. We assume that the mock
observations are carried out when the planet is near $\phi=90^{\circ}$ and
near $\phi=270^{\circ}$, where $\phi$ is the orbital longitude (Figure 4).
This is because the orbital inclination would be constrained best when the
data cover the orbital phases where the radial velocity changes the most.

For a planet around a M5 (M8) star, mock observations are continued for 1
(0.25) day centered at $\phi=90^{\circ}$ and for another 1 (0.25) day centered
at $\phi=270^{\circ}$, which cover approximately $32^{\circ}$ ($40^{\circ}$)
on each side. These two cadences are repeated until the parameters are
constrained, and the total integration time is recorded.

Figure 4: Observational configuration showing the geometrical parameters. The
data near $\phi=90^{\circ}$ and $\phi=270^{\circ}$ are used in this study.
Figure 5: Schematic figures showing the signals to be detected. Analysis (A)
(left) utilizes the high-frequency features of the spectra extracted by
subtracting the moving average from the spectrum, while Analysis (B) (right)
utilizes the difference of the Doppler shifted spectra from the average
spectra.

####  3.1.3 Instruments and observational configurations

Mock observations are carried out with a high-resolution spectrograph at 12-18
$\mu$m with the resolving power of $\mathcal{R}=30,000$. This bandpass and
resolving power are motivated by the latest specification of SPICA’s SMI/HR
($\mathcal{R}\sim 33,000$ in 12-18 $\mu$m) (e.g., Kaneda et al., 2018) and by
the expected resolution of OST/MISC (e.g., Sakon et al., 2018). The
dependences on the spectral resolution and the bandpass are examined in
Section 4.

The total throughput including the quantum efficiency of the detector is
assumed to be 0.2. We note that the throughput of the optical system of
existing high-resolution spectrographs installed on ground-based telescopes
have reached approximately 60% (e.g., Ikeda et al., 2016, 2018), and the
throughput of future instruments could be higher than what we assume here. It
is trivial to scale our results by throughput; see Section 3.3.1.

We collect data every 1800 sec (i.e., exposure time is assumed to be 1800
sec). Due to the change of the planetary radial velocity, the planetary
spectrum is Doppler shifted relative to the host star spectrum on the detector
plane. Over the course of the planetary orbital motion, the planetary spectrum
moves beyond the resolution elements unless the orbital inclination is close
to zero (i.e., face-on orbit), as the radial velocity amplitudes shown in
Table 2 (and the assumed orbital inclination of $60^{\circ}$) are larger than
the velocity corresponding to the resolution element, $c/\mathcal{R}\sim 10$
km/s. )

####  3.1.4 Signal and noise

The total photoelectron count that the detector receives is the summation of
the planetary light ($\mathcal{N}_{{\rm p}}$), stellar light
($\mathcal{N}_{\star}$), zodiacal light ($\mathcal{N}_{{\rm zodi}}$), thermal
background of the telescope ($\mathcal{N}_{{\rm tele}}$), and the dark current
($\mathcal{N}_{{\rm dark}}$). We assume that the contribution from the zodical
light and the dark current are perfectly subtracted through a post-processing.
This leaves $\mathcal{N_{\rm total}}=\mathcal{N}_{{\star}}+\mathcal{N}_{{\rm
p}}$ alone as a signal.

description | symbol | value
---|---|---
spectral resolution | $\mathcal{R}$ | 30,000
telescope diameter | $D$ | 6.5 m
total throughput | $\xi$ | 0.15
distance to target | $d$ | 5 pc
planetary radius | $R_{\rm p}$ |  $R_{\oplus}(=6.371\times 10^{6}$ m)
exposure time | $\tau_{\rm exp}$ | 1800 sec
Table 3: Fiducial values for observational parameters. The scaling of the
required integration time by these parameters are given in equation (13).

The shot noise from all of these factors contributes to the observational
noise, although the shot noise due to the planetary flux is negligible
compared to that from the stellar flux. An additional factor that is taken
into account is the read noise. The systematic noises such as a fringe of the
detector are assumed to be perfectly removed; the effect of systematic noise
is discussed in Section 5.2.2. These assumptions imply the Gaussian random
noise with the following standard deviation for $j$-th wavelength element at
wavelength $\lambda_{j}$:

$\displaystyle\sigma^{2}_{{\rm photon},j}=$ $\displaystyle\mathcal{N}_{{\rm
star},\,j}(1+\xi b_{j}(T_{\star}))+\mathcal{N}_{{\rm zodi},\,j}$ (4)
$\displaystyle+\mathcal{N}_{{\rm tele},\,j}+\mathcal{N}_{{\rm
dark},\,j}+\sigma_{\rm read}^{2}$

where $\sigma_{\rm read}^{2}$ represents the read noise, $b(T)$ is the Bose
factor, and $\xi$ is the total throughput.

The Bose factor, $b(T)$, is to take account of the sub-Poissonian nature of
the photon count statistics at $h\nu\ll kT$ (e.g., Boyd, 1982):

$b_{j}(T)=\frac{1}{\exp\left(\frac{hc}{\lambda_{j}kT}\right)-1}$ (5)

For MIR observations ($\lambda\sim 15\,\mu$m) of M-type stars ($\sim 3000$ K),
$b\sim 2.7$. Multiplied by the throughput ($\xi=0.2$), the Bose factor
increases the variance of the photon count by about 50%. The Bose factor for
the zodiacal light is not included due to the low temperature.

In the fiducial cases studied in this paper, the dominant noise source is the
stellar flux. The contributions of these factors in general cases are
presented in Figure 6, as a function of the spectral type of the star and the
distance to the target.

Figure 6: Photon count of the zodiacal light (black) and the starlight
(colors) at wavelength 10 $\mu$m (solid line) and at 20 $\mu$m (dashed line),
per exposure time (1800 sec) per wavelength resolution element
($\lambda/\mathcal{R}$). These numbers can also be compared with dark current
($360e^{-1}$/exposure) read noise, $14\times 4e^{-}$/exposure.

Each term in equation (4) is expanded below.

##### Planet/star spectra

The photoelectron counts of the planet and the host star per exposure are
simply:

$\displaystyle\mathcal{N}_{{\rm p}}(\lambda)$ $\displaystyle=$
$\displaystyle\mathcal{F}L_{\rm p}(\lambda)\cdot\pi R_{\rm
p}^{2}\left(\frac{1}{d}\right)^{2}$ (6)
$\displaystyle\mathcal{N}_{\star}(\lambda)$ $\displaystyle=$
$\displaystyle\mathcal{F}L_{\star}(\lambda)\cdot\pi
R_{\star}^{2}\left(\frac{1}{d}\right)^{2}$ (7) $\displaystyle\mathcal{F}$
$\displaystyle\equiv$
$\displaystyle\pi\left(\frac{D}{2}\right)^{2}\xi\tau_{\rm
exp}\frac{\lambda}{\mathcal{R}}$ (8)

where $L_{\rm p}$ and $L_{\star}$ represent the radiance of the planet and
star, respectively, while $R_{\rm p}$ and $R_{\star}$ represent the radius of
the planet and the star, respectively. The meanings of other parameters are
summarized in Table 3 together with our fiducial values. Here $\mathcal{F}$ is
used to denote the common factor related to the observational configuration.

##### Zodiacal light

The zodiacal flux of the Solar System is set to 15 MJy/sr (Glasse et al.,
2015) and is assumed to be constant over the observing wavelength range. The
zodiacal light in the target system is not included in our simulation because
it is much fainter than the other sources. Thus,

$\mathcal{N}_{{\rm zodi}}=\mathcal{F}f_{\rm zodi}\cdot\pi(\theta_{\rm
aperture})^{2}$ (9)

where $f_{\rm zodi}$ is the zodiacal flux per wavelength per steradian (rather
than per frequency as represented in Jansky per steradian) and $\theta_{\rm
aperture}$ denotes the aperture radius. The aperture radius should be
determined by the balance between the shot noise of the thermal background and
the systematic noise (Matsuo et al., 2018). Based on the analytical
formulation on the relation between the aperture radius and systematic noise
(Itoh et al., 2017), the aperture radius is set to 1.85 arcsecond,
corresponding to 4 times the diffraction limit at 15 $\micron$ such that the
systematic noise is reduced down to 100 ppm under a condition that the
pointing jitter of the Origins Space Telescope is 22 milli-arcsecond (RMS),
corresponding to approximately 0.05 $\lambda/D$ at 15 µm (Leisawitz et al.,
2018).

##### Telescope background

The thermal light and stray light from a telescope assembly also contributes
to the background light. For the JWST MIRI, the telescope background can be
approximately fitted by a combination of several blackbody radiations with
temperatures ranging from 50 to 70 K (Glasse et al., 2015). When the telescope
is cooled down to below 10 K, as planned for SPICA and OST, and the telescope
background becomes negligible at $<$ 100 $\mu$m.

##### Dark current

The dark current is assumed to be 0.2 $e^{-}$ s-1 based on the performance of
the Si:As detector for the JWST MIRI (Rieke et al., 2015), but this never
becomes of relative importance in our study.

##### Read noise

The read noise is assumed to be 14 $e^{-}$/read assuming the Fowler-eight
sampling (Rieke et al., 2015). The loss time due to the readout is not
considered in the paper. Each resolution element is sampled by 4 pixels. With
these assumptions, read noise never becomes dominant for the exposure time of
1800 sec. We note that 1800 sec is longer than the typical values for JWST.
The reduction in the exposure time can lead to a significant contribution from
read noise.

###  3.2 Analysis

Because the planetary signal is a tiny portion in the total spectrum, it is
critical to subtract the stellar spectra precisely from the total spectrum.
Previous studies discussing the detectability of atmospheric features of non-
transiting habitable planets assumed that the host star spectrum is determined
precisely (e.g., Kreidberg & Loeb, 2016; Snellen et al., 2017). The analysis
based on this optimistic assumption is denoted by Analysis (A) and is
described in Section 3.2.1.

However, M-type stars are rich in spectral features as shown in Figure 7 and
its MIR stellar spectrum may not be determined to a satisfactory accuracy, due
to the inaccurate line list and to the uncertainties in the stellar
atmospheric structure. Indeed, the models of star spectra do not fit the near-
infrared observed spectrum of M-type stars to the precision of observed data
(e.g., Zhang et al., 2018; Wakeford et al., 2019). In the case where we do not
have any prior knowledge about the stellar spectrum, we have to focus on the
Doppler-shifted components in the combined spectrum of the star and the planet
to extract planetary one. Such an analysis is denoted by Analysis (B) and is
described in Section 3.2.2.

The planetary signals to be fitted in the case of Analysis (A) and (B) are
illustrated in Figure 5. Importantly, as described in more detail below, the
signal in Analysis (B) is smaller than that in Analysis (A) because some
fraction of the spectral features are cancelled out.

Figure 7: MIR spectra based on BT-Settl model (Allard et al., 2012).

####  3.2.1 Analysis (A): with a well-modeled stellar spectrum

In the first analysis, we follow the procedure of Snellen et al. (2017),
assuming that the fine structure of the stellar spectrum can be modeled
precisely. Even with this assumption, we cannot separate the host star and the
planet spectra from observations alone, so we take the following steps:

  1. 1.

The model stellar spectrum is fitted to the data (that is the combined
spectrum of the star and the planet), by varying the scaling of the spectrum
as well as the stellar parameters.

  2. 2.

The best-fit stellar spectrum is subtracted from the data.

  3. 3.

The residual spectrum is corrected by subtracting its moving-average.

  4. 4.

The model planet spectrum, which is also corrected through the moving-average
subtraction (the bottom panel of Figure 5), is fitted to the corrected
residual spectra.

In Step 1, we use the input BT-Settl model and fit for the absolute scale. The
best-fit scale is slightly larger than the input stellar model, due to the
contribution from the planet spectrum. In other words, some fraction of the
planetary spectrum is subtracted in Step 2. This results in a trend in the
residual spectrum, which is corrected in Step 3.

Finally in Step 4, the high-frequency features are fitted by the model planet
spectrum. For simplicity, our fitting model is the same as the input model,
similar to Snellen et al. (2017), and we consider two fitting parameters: the
planet-to-star contrast $C$ and the orbital inclination $I$. The parameter
estimate is based on the posterior probability,
$\mathcal{P}(c,i|\\{\mathcal{N}_{i}(t_{k})\\})$:

$\mathcal{P}(C,I|\\{\mathcal{N}_{{\rm
res},jk}\\})=\mathcal{L}(\\{\mathcal{N}_{{\rm res},jk}\\}|C,I)\Pi(C)\Pi(I)$
(10)

The $\mathcal{N}_{{\rm res},jk}$ represents the corrected residual spectra at
$j$-th wavelength elements and $k$-th observation epoch. The likelihood is
simply:

$\displaystyle\mathcal{L}(\\{\mathcal{N}_{{\rm
res},jk}|c,i)\propto\exp\left(-\frac{\chi^{2}}{2}\right)$ (11)
$\displaystyle\chi^{2}\equiv\sum_{\\{j,k\\}}\left(\frac{\mathcal{N}_{{\rm
res},jk}^{\rm obs}-C\mathcal{N}_{{\rm res},jk}^{\rm
theory}(I)}{\sigma_{j,k}}\right)^{2}.$ (12)

We assume a flat prior for contrast $C$ between 0 and 10, and a flat prior for
the orbital inclination $I$ between 0 to 90 degree; $\Pi(C)=const.$ and
$\Pi(I)=const.$ This means that the posterior probability is simply equivalent
to the likelihood function. The posterior probability is normalized so that
the total is unity. The $1\sigma$, $2\sigma$, and $3\sigma$ contours of the
posterior probability are close to those of $\Delta\chi^{2}$ measure (e.g.,
Snellen et al., 2017), with a slight difference due to the non-linear
$I$-dependence of the model.

####  3.2.2 Analysis (B): with unknown stellar spectrum

In the second analysis, we do not assume that the stellar spectrum is known a
priori. Instead, we consider the time average of the combined spectra of the
star and the planet and assume that the stellar spectrum is fully included in
this averaged spectrum, which is valid if the stellar spectrum is stable in
time. On the other hand, the spectral features of the planet should remain in
the average-subtracted residual spectrum, due to the Doppler shift caused by
its orbital motion.

The analysis procedure is as follows:

  1. 1.

The average of the spectra at different planetary orbital phases is obtained.

  2. 2.

The average spectrum is subtracted from the spectrum at different orbital
phases.

  3. 3.

The model planet spectrum, which is also Doppler-shifted and average-
subtracted (the bottom panel of Figure 5), is fitted to the set of residual
spectra.

The parameter estimate in Step 3 is performed in the same way as Analysis (A)
(equations (10)-(12)).

Figure 5 illustrates the signal that can be used in this analysis. The average
of the spectrum at $\phi=90^{\circ}$ and that at $\phi=270^{\circ}$ is shown
in the black line in the upper panel, and the average-subtracted spectra are
shown in the lower panel. During this subtraction process, some portions of
the planetary spectral features are cancelled out, reducing the signal level.
This effect is substantial in the case of MIR observations of potentially
habitable planets whose strongest lines are broader than the Doppler shift
(Figure 3). The situation is different from the typical high-resolution
spectroscopy of hot Jupiters at the shorter wavelengths, where the lines are
narrower and the Doppler-shift is larger. In such a case, the Doppler shift
moves the line beyond their intrinsic width and the average-subtracted spectra
are more similar to the original spectrum.

####  3.2.3 Notes on cross-correlation analysis

High-resolution spectroscopy of hot-Jupiter systems have been routinely
analyzed through cross-correlation function (Snellen et al., 2010; Brogi et
al., 2013). While this technique is useful for detecting high-frequency
features, its statistical treatment and the uncertainties in the estimated
model parameters are not straightforward (e.g., Brogi et al., 2017). On the
other hand, fitting observed data with theoretical models by minimizing the
sum of squared residuals offers a more plain interpretation (Snellen et al.,
2017). Thus, in this paper, we employ the latter method.

###  3.3 Results

Figure 8: The 1$\sigma$ (solid lines), 2$\sigma$ (dashed lines), and 3$\sigma$
(dotted lines) confidence intervals for planets around a mid-M star (upper
panels) and a late-M star (lower panels), based on analysis (A), i.e., when
the stellar spectrum is precisely determined. Both the mock data and the
fitting model are based on the model thermal emission spectra presented in
Figure 3. The constraints from two integration times are presented for some
cases, in order to show how the constraints are developed as the integration
time becomes longer. The blue colors imply that the presence of the molecule
is detected (i.e., the $C=0$ is rejected by 3$\sigma$) and the inclination is
not well constrained, while the red colors imply that both the contrast and
inclination are constrained. The corresponding signal-to-noise ratios of the
host star spectrum per wavelength resolution element (SN) are also reported.
Figure 9: Same as Figure 8, but the mock data are now analyzed by analysis
(B).

Figures 8 and 9 summarize the results of Analysis (A) and Analysis (B),
respectively, performed on the thermal emission spectra with all molecules
included and with individual molecules in isolation (Figure 3). The solid,
dashed, and dotted lines present the 1$\sigma$, 2$\sigma$, and 3$\sigma$
constraints in the contrast-inclination plane (i.e., 68.27%, 95.45%, and
99.73% probability), after the indicated integration time. The corresponding
average signal-to-noise ratio (SN) of the host star spectra per wavelength
resolution element is also presented.

As expected, detecting any molecule is easier for planets around late-M stars
than around mid-M stars due to the larger planet-to-star flux ratio and, for
analysis (B), to the larger radial velocity amplitude. In the model considered
here, the total observation time required for a M5-star system is by a few
times larger than a M8-star system, or 3 times smaller signal-to-noise ratio
of the host star.

When the host star spectrum is accurately known and corrected (Analysis (A);
Figure 8), the spectral features of an Earth-like atmosphere at 5 prsecs away
can be detected within a few days of total observation time, even around mid-M
systems (the leftmost column). Compared to this optimistic scenario, Analysis
(B) requires approximately 4 times longer observations to detect the spectral
features of the same target, or doubled signal-to-noise ratio of the host star
flux. This is because of the self-subtraction discussed in Section 3.2.2.

For comparison, we also perform the same analysis assuming only one of the
molecules as opacity source (the lower five panels of Figure 3); the results
of which are shown in the right four columns of Figures 8 and 9. Comparison of
the two figures clearly shows that the effect of the self-subtraction is
substantial except for O3. In addition, for Analysis (B), the contrast and the
inclination angle tend to degenerate for pure CO2 or pure H2O atmospheres.
This is because when a small orbital inclination (i.e., close to face-on) is
assumed, the Doppler shift of the planetary spectrum is small and the
amplitude of the differential spectrum (the lower panel of Figure 5) is also
small, which is then compensated by a large contrast. Note that the orbital
inclination is poorly constrained even with Analysis (A) for pure CO2 or pure
H2O atmospheres, due to the broad nature of the spectral lines.

Compared to these two molecules, the self-subtraction is not substantial for
O3 because the O3 bands are densely populated by narrow lines, and the
detectability through Analysis (B) is similar to that of Analysis (A).
Furthermore, the orbital inclination is well constrained as soon as the
contrast is constrained to non-zero. N2O features are also relatively good at
constraining the inclination, despite its relatively weak features.

In the spectrum with all the molecules included, the features of each molecule
tend to be more muted than the features of individual molecules in isolation,
due to the line overlaps. However, the spectrum become more rich in features,
and the constraints on the contrast is as good as for the pure-CO2 or pure-H2O
cases. Furthermore, the overlaps of the spectral features of individual
molecules break the degeneracy between the contrast and inclination for
Analysis (B), significantly improving the constraints on the orbital
inclination compared to the pure-CO2 or pure-H2O cases. As a result, both the
contrast and the orbital inclination can be reasonably constrained within
$\sim 1$ week of integration for M8 star systems within 5 parsecs even with
Analysis (B). For earlier type systems, constraining the orbital inclination
would become more challenging.

####  3.3.1 Scaling and the application to the known targets

While the estimated integration time would be too long to be practical except
for the case of late-M stars, it is possible to approximately scale these
numbers for varying observational configurations. When the stellar light is
the dominant noise source, the integration time that achieves certain noise
level, $\tau_{\rm esp,\,0}$, is proportional to:

$\tau_{\rm
esp,\,0}\propto\left(\frac{R_{p}}{R_{\oplus}}\right)^{-2}\left(\frac{d}{\mbox{5
pc}}\right)^{2}\left(\frac{D}{\mbox{6.5
m}}\right)^{-2}\left(\frac{\xi}{0.2}\right)^{-1}$ (13)

In particular, the observations of the nearest possible target, Proxima
Centauri b, at the distance of 1.3 pc and with the estimated radius of $\sim
1.1R_{\oplus}$, would only require $\sim$1 days of observation even without an
assumption on the stellar spectrum.

A more SPICA-like specification with $D=2.5$ m and a lowered total throughput
($\xi=0.1$) yield $\sim$14 days of total integration time.

##  4 Dependence of Detectability on Spectral Resolution, Molecular Abundance,
and Bandpass

In the assessment in Section 3, we have made several assumptions for the
specifications of the spectrographs and the atmospheric profile. In this
section, we discuss the effects of varying these assumptions on the
detectability of the molecules.

###  4.1 Effects of Spectral Resolution

First, we discuss the dependence on the spectral resolution of the
spectrograph. Let us represent the signal by the total area of the shaded
region shown in Figure 5, namely the integration of the absolute value of the
differential spectrum both in Analysis (A) and (B); For Analysis (A), the
differential spectrum corresponds to the difference between the original
spectrum and the moving average, while for Analysis (B) it is the difference
between the spectrum at a certain orbital phase and the time-averaged
spectrum.

Figure 10: The total areas of the filled region in Figure 5, one of the
diagnostics of the signal, as a function of the spectral resolution, using
Analysis (A) (left) or Analysis (B) (right). The fiducial atmospheric profile
is assumed.  Figure 11: 3$\sigma$ constraints in the contrast-inclination
plane based on the observations of Proxima Centauri b with CO2 features using
varying spectral resolution, after 2 days (left) or 7 days (right) of
integration.

Figure 10 shows the signals defined above for Analysis (A) and (B) as a
function of the spectral resolution. Black line shows the case of the spectrum
with all four molecules included (the bottom panel of Figure 3), while other
lines show the cases for the spectra of individual molecules. As expected, the
signal level monotonically increases as the spectral resolution increases.
However, it is mostly saturated above $\mathcal{R}\gtrsim 10,000$ except O3,
in both analysis (A) and (B). This is partly because the interval of the
strongest lines are $\sim 0.1-0.2\,\mu$m and these lines start to be resolved
with $\mathcal{R}\gtrsim 2000$. In addition, the intrinsic line broadening
starts to dominate over the broadening determined by the spectral resolution
for $\mathcal{R}\gtrsim 10,000$. The required resolution, $\mathcal{R}\gtrsim
10,000$, is about one order of magnitude smaller than that in the visible or
near-infrared domain ($\lambda\sim 1\,\mu$m), $\mathcal{R}\gtrsim 100,000$,
because the photon energy is about one order of magnitude lower while the
energy corresponding to the line broadening and the interval of lines do not
depend on the wavelengths (See the discussion in Section 2.2).

The saturation is less prominent for O3, where the narrow lines are densely
populated. Thus, detection O3 benefits from the increased spectral resolution
beyond $\mathcal{R}\sim 10,000$.

Figure 11 shows the constraints for the contrast and the orbital inclination
obtained through Analysis (B) of mock observations of Proxima Centauri b for 2
days (left) and 7 days (right) with varying resolution. As expected from the
above assessment, the detectability, namely the lower constraints on the
contrast, improves as the spectral resolution increases up to $\mathcal{R}\sim
10,000$, above which no significant improvement is seen. However, the
constraints on the orbital inclination is improved beyond $\mathcal{R}\sim
10,000$.

###  4.2 Choice of the bandapss

Our fiducial bandpass is set to 12-18 $\mu$m, based on the planned
specification of SPICA/SMI and OST/MISC. In this range of wavelengths, several
other molecules have prominent absorption bands, in addition to the molecules
studied in this manuscript. They includes NH3, SO2, and NO2 and they may also
be useful for characterizing the surface environment of potentially habitable
planets.

Although there is no planed high-resolution spectrograph that goes beyond this
bandpass, it would be worth considering the impact of the alternative choice
of bandpass for future instrumental designs. In the following, we briefly
discuss how the extension of the bandpass would affect the prospects.

If the bandpass is extended toward shorter wavelengths, a feature of
particular interest would be O3 9.7 $\mu$m band (the bottom panel of Figure
3). The peak line strength of this band is approximately one order of
magnitude larger than the 14.5 $\mu$m band and, combined with the fact that
the shorter wavelength is in the Wien’s regime of the Planck function, the
feature relative to the continuum is deeper. However, the reduced planet-to-
star contrast (Figure 1) does not make this band easy to detect. We find that
detection of O3 9.7 $\mu$m O3 band would require approximately 1.7 times
longer integration time (or 75% photon noise) than the detection of 14 $\mu$m
O3 band, in the absence of CO2 (not shown). In the presence of CO2 larger than
1 ppm, 9.7 $\mu$m band would be the only detectable band. The detectability of
O3 9.7 $\mu$m would be improved for planets with higher surface temperature.
The features at even shorter wavelengths is even less likely to be detectable
with this method due to the small planet-to-star contrast.

At the wavelengths longer than 18 $\mu$m, there are few vibrational modes of
simple molecules. However, the rotational lines of H2O extend beyond 18 $\mu$m
and adding longer wavelengths improves the detectability of H2O. Using 12-24
$\mu$m bandpass, the integration time required for the detection of H2O from
the H2O-only spectrum (the second panel of Figure 3) would be reduced to
approximately 60%.

###  4.3 Molecular abundances that can be most sensitively probed

The detectability also depends on the actual abundance of the molecules. While
in transmission spectroscopy the detectability of certain molecules generally
increases as the abundance increases (ignoring the effects of clouds), this is
not the case for the high-resolution method considered in this paper. It is
true that the signal would not be detected if the abundance is too small.
However, larger abundances does not always lead to the improved detectability,
because some of the features start to saturate. For example, with the assumed
Earth-like abundance of CO2 ($\sim$ 330 ppm) the strongest lines of CO2 around
15 $\mu$m are saturated (i.e., the spectrum is flattened) and do not
significantly contribute to the signal. With larger abundance, the flattened
region only increases, although the features associated with smaller opacities
kick in. Thus, there is the optimal range of molecular abundance.

Figure 12: Same as Figure 10 but with varying mixing ratio of CO2.

To see such a dependence of detectability of CO2 and N2O, we created model
spectrum with varying mixing ratio assuming a 1-bar atmosphere with vertically
uniform mixing ratio of each molecule. For each model spectrum, we repeated
the analysis of Section 4.1 to see the dependence of the “signal” as a
function of the spectral resolution. The result is shown in Figure 12. The
signal level is maximized when the mixing ratio is 1-103 ppm. A similar trend
was found for N2O, with the signal level being maximized for 1-102 ppm.
Therefore, this technique is sensitive to relatively small amount of
molecules. This would be rather complementary to scattered-light observations
where larger abundance of these molecules would be probed; see Section 5.1.

##  5 Discussion

###  5.1 Potential Targets and Synergy with other techniques

The prime targets of the method studied in this manuscript would be
potentially habitable planets around mid- to late-M stars within 5 parsecs.
The fact that this method does not require planetary transits is a significant
advantage over transmission or eclipse spectroscopy, given the low transit
probability ($\sim$5%). Within approximately 5 parsecs, 6 mid- to late-M stars
have already been found to harbor Earth-sized planets around habitable zones
as of the submission of this manuscript (Table 4), none of which are found to
transit in front of the host star. The high-resolution technique would provide
valuable opportunities to investigate the atmospheric compositions of these
known planets in thermal emission. The number of the targets may increase
thanks to the current vigorous radial velocity measurements in near-infrared
all over the world (There are 7 late-M stars (M6 or later) and 10 mid-M stars
(M3-M5) within 5 parsecs, according to Gaidos et al. (2014)). Note that we do
not necessarily limit our scope to 5 parsecs. Some of the features that are
easier to detect may extend its scope beyond 5 parsecs.

These nearby M-star systems are also considered to be the prime targets for
the scattered-light observations with the next-generation extremely-large
ground-based telescopes (Fujii et al., 2018). The combined efforts would allow
for the comprehensive characterization of the target systems, as the
information obtained from scattered-light spectra (e.g., surface liquid water,
clouds, certain molecules) and that obtained from thermal emission spectra
(temperature, radius, certain molecules) are complementary. For example,
thermal emission spectra could in principle constrain the planetary radius,
which critically helps the interpretation of the scattered-light spectrum.
Furthermore, for a given molecule, the abundance detectable in the scattered
light at $\sim$1 $\mu$m would be larger than the abundance probed by the high-
resolution technique (Section 4.3), thus the combination would constrain a
wide range of atmospheric properties.

The targets of MIR high-resolution technique are also complementary to the
prime targets for transmission and eclipse spectroscopy. The main scope of OST
is currently to characterize potentially habitable planets through transit
spectroscopy by achieving ultrahigh stability (see the technical report of
OST). There are pros and cons for transit spectroscopy and for the technique
studied in this paper. While the former may require less total observation
time for certain biosignature molecules (provided that the noise floor could
be down to a few ppm) than the latter, it has to visit the target 50 times for
TRAPPIST-1 planets (Tremblay et al., 2020) and more for less favorable
targets, likely spanning 1 year or longer. The number of visits of the high-
resolution method studied in this manuscript can be much smaller because of
the longer observation time per visit. In addition, the transmission
spectroscopy could be hampered by high-altitude thin clouds, which does not
significantly affect the high-resolution technique. Thus, combining both
techniques would expand our capability of investigating potentially habitable
worlds.

Table 4: Nearby systems that are found to harbor Earth-sized planets around habitable zones. Stellar parameters are based on Gaidos et al. (2014).  Star name | $d$ | Spectral type | $T_{\star}$
---|---|---|---
Proxima Cen | 1.3 pc | M7 | 2883 K
Ross 128 | 3.4 pc | M5 | 3145 K
Gliese 1061 | 3.7 pc | M6? | 3000 K?
Luyten’s | 3.8 pc | M4 | 3317 K
Teegarden’s | 3.9 pc | M7 | 2700 K
GJ 682 | 5.1 pc | M5 | 3190 K

###  5.2 Effects of other noise

####  5.2.1 Time variability of stellar spectra

The major concern of this method is the variability of high-resolution stellar
spectra. Aside from sporadic flares, stars generally exhibit photometric and
spectroscopic variabilities mainly due to starspots. Given that the signal
from the planet is a tiny fraction of the stellar spectrum, such variability
potentially cases serious problem in the analysis.

There is a growing amount of modeling efforts for the temperatures and
covering fractions of star spots/faculae of mid- to late-M stars based on
observations. One piece of evidence comes from Doppler imaging. Doppler
imaging of mid- to late-M stars (Morin et al., 2008; Barnes et al., 2015,
2017) estimates that around a few % of the surface is covered with spots that
are assumed to be cooler than other area by 200-400 K. Recently,
characterization of TRAPPIST-1, an M8 star harboring transiting temperate
Earth-sized planets, has attracted attention. The light curve of TRAPPIST-1
observed by Kepler/K2 (0.43-0.89 $\mu$m) shows $\sim 1$% peak-to-peak
variation amplitude, which may be explained by rotating dark spots or faculae.
In contrast, the variability of Spitzer at 4.5 $\mu$m do not show a clear
variability and is at least smaller than the photometric precision, which is
limited to the order of 0.1% owing to the shot noise and instrumental
systematic noise. Morris et al. (2018) find that the combination of these
light curves at different wavelengths can be explained by the existence of
several very small faculae with $>5300$ K rather than dark spots, and show
that the theoretical photometric variability at 4.5 $\mu$m with this model
seems to be random with an amplitude of about 100 ppm. This prediction is
encouraging because this is compatible or smaller to the contrast of the
planet to stellar signals. The impact of variability on this method will be
further mitigated by averaging the data over various stellar rotational
phases.

The covering fraction of the spots and faculae are also estimated based on the
multi-component fit to the stellar spectrum using the theoretical stellar
spectra with various temperatures (e.g., Zhang et al., 2018; Wakeford et al.,
2019). They typically result in relatively large covering fractions for darker
spots and faculae (about 30-50%). However, this method is sensitive to the
small-scale features that spread across the photosphere as well as the large-
scale asymmetric patterns, in contrast to the Doppler Imaging and
reconstruction from photometric variabilities that are only sensitive to the
latter. The small-scale features that distribute rather uniformly do not
affect the observations studied in this manuscript.

We note that the recent observations of Proxima Centauri by Suárez Mascareño
et al. (2020) show V-band photometric variability of $>$10% (peak-to-peak),
larger than the Kepler/K2 variability of TRAPPIST-1. In reality, stars exhibit
the diversity in their surface magnetic activities and the resultant
spots/faculae, and therefore it is critical to characterize individual stars
before characterizing the planets. With the launch of JWST that achieves
higher precision than Spitzer thanks to its larger aperture and better
attitude control, our understanding of variabilities of the stellar MIR
spectra will be developed and the techniques to cope with these variability
will be advanced. For example, simultaneous light-curve observations at
different wavelengths will allow us to better model the spots, as demonstrated
by Morris et al. (2018). Such practice will further enhance the feasibility of
the method presented in this manuscript.

####  5.2.2 Systematic Noise

We did not explicitly deal with the impact of the unknown instrumental
systematic errors. Because this technique utilizes high-frequency features of
the spectrum, the uncertainties in the absolute flux or the low-frequency
modulation of the continuum, even if it is substantial, will not affect our
results after applying some corrections (e.g., high-pass filter; Snellen et
al., 2017).

In contrast, residual flat-fielding, spectral fringe patterns, and telescope
pointing jitters form artificial absorption and emission lines (i.e., high-
frequency systematic noise) in the observed spectrum. Such artificial high-
frequency patterns will be problematic for both the analyses (A) and (B). The
high-frequency noise will contribute to the total noise floor because of its
random distribution along the wavelength. Since the high-frequency systematic
noise is independent of the shot noise, the total noise may be expressed by
$\sqrt{\sigma^{2}_{{\rm photon},j}+\sigma^{2}_{{\rm sys},j}}$, where
$\sigma_{{\rm photon},j}$ is defined by equation (4) and $\sigma_{{\rm
sys},j}$ denotes the high-frequency systematic noise at the $j$-th spectral
element. Given that the noise required for detection is order of 100 ppm, it
is necessary that the instrument systematic noise is less than this level. If
the systematic noise is $\sim 30$ ppm or less, the systematic noise will not
affect our estimate much. Conversely, if it is comparable to the photon noise,
a substantially larger integration time will be required in order to reduce
the photon noise further. We note that the high-resolution spectroscopy
observes the shift of the planetary Doppler signal on the detector, which
could mitigate an impact of the patterns fixed to the detector (owing to
relative sensitivity variations between adjacent spectral channels).

###  5.3 Variety of Planet Thermal Emission Spectrum

####  5.3.1 Stratospheric Thermal Inversion

Figure 13: Left: high-resolution ($\mathcal{R}=30,000$) thermal emission
spectra assuming the vertical profiles of temperature and CO2 that are same as
Earth’s “US standard” profile. Right: Parameter constraints based on the left
panel, assuming a late-M host star at 5 parsec.

This section discusses the variety of high-resolution thermal emission
spectra, and how they affect the parameter estimate.

Our fiducial model atmosphere has a relatively cool stratosphere,
approximately 150 K. A higher stratospheric temperature results in shallower
spectral features, increasing the required integration time for detection.

If the planet has thermal inversion in the upper atmosphere (like the ozone
layer of Earth), thermal spectra may exhibit emission lines. These lines can
be sharp, because the pressure in the upper atmosphere (where the emission
line are formed) is low and the line width is relatively narrow. This can
relax the degeneracy between the contrast and the orbital inclination which
would otherwise be present, like our fiducial cases for CO2 and H2O.

As a demonstration, Figure 13 is the modeled thermal emission spectrum that
assumes an Earth-like temperature profile with a stratospheric thermal
inversion as well as an Earth-like CO2 vertical profile. Note the narrow
emission lines within 14-16 $\mu$m. The analysis (B) performed on this
spectrum yields the constraints shown in the right panel of Figure 13,
assuming a late-M host star. While the total integration time required for
detection is longer than our fiducial case due to the reduced line depths
(corresponding to the warmer upper atmosphere), the degeneracy between the
contrast and the inclination is less severe.

####  5.3.2 Horizontal temperature gradient

We did not include the horizontal temperature gradient of the planet. In
reality, it can be substantial near the surface, depending on the heat re-
distribution efficiency of the atmosphere and/or ocean. The horizontal
temperature gradient would result in the time variation of the disk-averaged
thermal emission as the planet rotates (e.g., Knutson et al., 2007). While
this itself can be used to characterize planet atmosphere, such variations may
complicate the analysis of high-resolution spectroscopy. Fortunately, the
effect of the horizontal temperature gradient will be mitigated if the data at
the same phase angle are stacked and analyzed. Our analysis assumes
observations near 90 degrees and 270 degrees (i.e., the same phase angle) so
the impact will be limited unless there is an extreme asymmetry between the
eastern and western hemispheres.

####  5.3.3 Full Exploration of Parameter Space

Although in this paper we attempted to constrain the contrast and inclination
by fitting the mock data by specific model spectra, in reality, the spectral
shape of the planet thermal emission also has to be constrained. In other
words, the distribution of the molecules as well as the temperature structure
should be fitting parameters as well, and the dimension of the parameter space
is large. It is also possible that the constraints on certain parameters can
be affected by the prior probability when exploiting the Bayesian framework.
Full explorations of the parameter space and development of sophisticated
retrieval techniques will need further study. Also, our analysis implicitly
assumes that the parameter space explored includes the ground truth. If this
is not the case, the obtained constrains should be interpreted with a special
care, which is to be examined in future studies.

###  5.4 Significance of mid-infrared high-resolution spectroscopy for
characterization of potentially habitable planets

Because our analysis (A) follows the procedure of Snellen et al. (2017), we
compare our results with their estimate. Their modeled CO2 features in the
medium-resolution spectra is detected (by nearly 4 $\sigma$) after the photon
noise becomes $\sim$50 ppm. Using the same bandpass (13.2-15.8 $\mu$m) and a
similar temperature profile (left panel of Figure 13), our analysis (A) of
high-resolution spectra was able to detect the CO2 features when the photon
noise is $\sim$100 ppm. This is likely a consequence of about two times
sharper features at high-resolution spectra, as well as the minor features
that are smoothed in the medium resolution spectra.

Note that it is not trivial to compare our estimate for the integration time
to the integration time estimated in Snellen et al. (2017) for JWST/MIRI. This
is because the noise of JWST/MIRI, especially beyond 12 $\mu$m, is not
dominated by the photon noise of the star. Other noise, including the thermal
background of the telescope, zodiacal light (with the larger aperture size),
and/or the readout noise (due to the short exposure time assumed) suppresses
S/N per exposure by 2-3 times compared to the idealized case. This leads to
the substantial difference in the integration time required to achieve a
certain precision. Thus, it is critical to suppress the background noise in
order to be able to detect atmospheric molecules of nearby mid-M and late-M
stars in a reasonable amount of time.

##  6 Summary

In this paper, we study the use of an MIR high-resolution spectrograph mounted
on a cryogenic telescope for characterizing non-transiting temperate
terrestrial planets orbiting M-type stars. We modeled high-resolution thermal
emission spectra of an Earth-like atmosphere with CO2, H2O, N2O, and O3
assuming a simplified temperature profile composed of an dry adiabatic
troposphere and an isothermal stratosphere. We show that the MIR spectral
features of Earth-like atmospheres can be broader than the width of the
Doppler shift, depending on the atmospheric structure.

In order to reasonably identify the tiny planetary features in the combined
spectra of the star and the planet, it is critical to subtract the stellar
contribution precisely, as low-mass stars are particularly rich in spectral
features in the MIR. For non-transiting planets, the stellar spectrum cannot
be determined observationally as we always observe the star and planet
together. Considering this possible difficulty, we proposed to observe the
target at around $\phi=90^{\circ}$ and at $\phi=270^{\circ}$ and that the
differential spectra are fitted by the model. This process can reduce the
signal because some fraction of broad lines are self-subtracted. This effect
is substantial except for O3, increasing the total integration required for
detection by several times.

Nevertheless, the spectral features ($\mathcal{R}=30,000$) of an Earth-like
planet between 12-18 $\mu$m would allow us to constrain the contrast and the
orbital inclination within $\sim$1 days of total integration time, for an
M8-star system within 5 parsecs. For earlier-type stars, constraining the
orbital inclination may be challenging. Scaling these estimates through to
Proxima Centauri b, an Earth-like abundance of CO2 would be inferred in
$\sim$1 day of integration. In addition to Proxima Centauri b, several
potentially habitable Earth-sized planets have been already discovered within
5 parsecs, and they are the prime targets for this technique.

We find that the constraints on the contrast improves with higher spectral
resolution, but $\mathcal{R}\gtrsim 10,000$ do not result in significant
improvement except for O3 whose absorption bands are densely populated with
narrow lines. However, the constraints on the inclination would benefit from
the higher spectral resolution. We also find that this method is most
sensitive to relatively small amount of CO2 and N2O (1-103 ppm for a 1-bar
Earth-like atmosphere), where the higher abundance does not lead to better
detectability.

In this study, we did not include the stellar variability and systematic noise
in our simulation. Because the photon noise required for the detection of
these features is 100 ppm or larger, the stellar variability and systematic
noise must be suppressed to these levels in order to detect the planetary
features. The expected stellar variability in the mid-infrared wavelength
range based on the previous photometric observations of TRAPPIST-1 with
Spitzer is comparable or less than 100 ppm, which will not affect this method.
Suppression of the noise sources other than the shot noise from the stellar
flux (e.g., thermal background) is critical in reducing the required
integration time and making the observations realistic.

We thank Teruyuki Hirano for helpful discussions on high-resolution
spectroscopy and Klaus Pontoppidan for the information about the instrumental
noise of JWST. YF is supported by Grand-in-Aid from MEXT of Japan, No.
18K13601. TM is supported by Grand-in-Aid from MEXT of Japan, No. 19H00700.

## References

  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765
  * Anglada-Escudé et al. (2016) Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437
  * Barnes et al. (2017) Barnes, J. R., Jeffers, S. V., Haswell, C. A., et al. 2017, MNRAS, 471, 811
  * Barnes et al. (2015) Barnes, J. R., Jeffers, S. V., Jones, H. R. A., et al. 2015, ApJ, 812, 42
  * Beichman et al. (2014) Beichman, C., Benneke, B., Knutson, H., et al. 2014, PASP, 126, 1134
  * Beichman et al. (1999) Beichman, C. A., Woolf, N. J., & Lindensmith, C. A. 1999, The Terrestrial Planet Finder (TPF) : a NASA Origins Program to search for habitable planets
  * Boyd (1982) Boyd, R. W. 1982, Infrared Physics, 22, 157
  * Brogi et al. (2017) Brogi, M., Line, M., Bean, J., Désert, J. M., & Schwarz, H. 2017, ApJ, 839, L2
  * Brogi et al. (2013) Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2013, ApJ, 767, 27
  * de Wit et al. (2018) de Wit, J., Wakeford, H. R., Lewis, N. K., et al. 2018, Nature Astronomy, 2, 214
  * Deming & Seager (2017) Deming, L. D., & Seager, S. 2017, Journal of Geophysical Research (Planets), 122, 53
  * Fujii et al. (2017) Fujii, Y., Del Genio, A. D., & Amundsen, D. S. 2017, ApJ, 848, 100
  * Fujii et al. (2018) Fujii, Y., Angerhausen, D., Deitrick, R., et al. 2018, Astrobiology, 18, 739
  * Gaidos et al. (2014) Gaidos, E., Mann, A. W., Lépine, S., et al. 2014, MNRAS, 443, 2561
  * Gillon et al. (2016) Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456
  * Glasse et al. (2015) Glasse, A., Rieke, G. H., Bauwens, E., et al. 2015, PASP, 127, 686
  * Gordon et al. (2017) Gordon, I. E., Rothman, L. S., Hill, C., et al. 2017, J. Quant. Spec. Radiat. Transf., 203, 3
  * Greene et al. (2016) Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17
  * Ikeda et al. (2016) Ikeda, Y., Kobayashi, N., Kondo, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, 99085Z
  * Ikeda et al. (2018) Ikeda, Y., Kobayashi, N., Kondo, S., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, 107025U
  * Itoh et al. (2017) Itoh, S., Matsuo, T., Goda, S., Shibai, H., & Sumi, T. 2017, The Astronomical Journal, 154, 97
  * Kaneda et al. (2018) Kaneda, H., Ishihara, D., Oyabu, S., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Proc. SPIE, 106980C
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
  * Knutson et al. (2007) Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183
  * Kochanov et al. (2016) Kochanov, R. V., Gordon, I. E., Rothman, L. S., et al. 2016, J. Quant. Spec. Radiat. Transf., 177, 15
  * Kreidberg & Loeb (2016) Kreidberg, L., & Loeb, A. 2016, ApJ, 832, L12
  * Léger et al. (1996) Léger, A., Mariotti, J. M., Mennesson, B., et al. 1996, Ap&SS, 241, 135
  * Leisawitz et al. (2018) Leisawitz, D., Amatucci, E., Carter, R., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup, H. A. MacEwen, G. G. Fazio, N. Batalha, N. Siegler, & E. C. Tong, 1069815
  * Levine et al. (2009) Levine, M., Lisman, D., Shaklan, S., et al. 2009, arXiv e-prints, arXiv:0911.3200
  * Matsuo et al. (2018) Matsuo, T., Greene, T., Roellig, T. L., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Proc. SPIE, 1069844
  * Morin et al. (2008) Morin, J., Donati, J. F., Forveille, T., et al. 2008, MNRAS, 384, 77
  * Morris et al. (2018) Morris, B. M., Agol, E., Davenport, J. R. A., & Hawley, S. L. 2018, ApJ, 857, 39
  * Nakayama et al. (2019) Nakayama, A., Kodama, T., Ikoma, M., & Abe, Y. 2019, MNRAS, 488, 1580
  * Rieke et al. (2015) Rieke, G. H., Ressler, M. E., Morrison, J. E., et al. 2015, PASP, 127, 665
  * Robinson & Catling (2014) Robinson, T. D., & Catling, D. C. 2014, Nature Geoscience, 7, 12
  * Sakon et al. (2018) Sakon, I., Roellig, T. L., Ennico-Smith, K., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Proc. SPIE, 1069817
  * Sarugaku et al. (2012) Sarugaku, Y., Ikeda, Y., Kobayashi, N., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, 844257
  * Snellen et al. (2015) Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
  * Snellen et al. (2013) Snellen, I. A. G., de Kok, R. J., le Poole, R., Brogi, M., & Birkby, J. 2013, ApJ, 764, 182
  * Snellen et al. (2017) Snellen, I. A. G., Désert, J. M., Waters, L. B. F. M., et al. 2017, AJ, 154, 77
  * Suárez Mascareño et al. (2020) Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, A&A, 639, A77
  * Tremblay et al. (2020) Tremblay, L., Line, M. R., Stevenson, K., et al. 2020, AJ, 159, 117
  * Wakeford et al. (2019) Wakeford, H. R., Lewis, N. K., Fowler, J., et al. 2019, AJ, 157, 11
  * Zaghloul & Ali (2011) Zaghloul, M. R., & Ali, A. N. 2011, arXiv e-prints, arXiv:1106.0151
  * Zhang et al. (2018) Zhang, Z., Zhou, Y., Rackham, B. V., & Apai, D. 2018, AJ, 156, 178


# Spin dynamics simulation of the $Z_{2}$-vortex fluctuations

Yo P. Mizuta1 Kazuki Aoyama2 Keisuke Tomiyasu3 Masato Matsuura4 and Hikaru
Kawamura5<EMAIL_ADDRESS>School of Engineering
science1Graduate School of Engineering science Osaka University Osaka
University Toyonaka Toyonaka Osaka 560-8531 Osaka 560-8531 Japan
2Department of Earth and Space Science Japan
2Department of Earth and Space Science Graduate School of Science Graduate
School of Science Osaka University Osaka University Toyonaka Toyonaka Osaka
560-0043 Osaka 560-0043 Japan
3Nissan ARC Ltd. Japan
3Nissan ARC Ltd. Natsushima-cho 1 Natsushima-cho 1 Yokosuka 237-0061 Yokosuka
237-0061 Japan
4CROSS Neutron Science and Technology Center Japan
4CROSS Neutron Science and Technology Center IQBRC Bldg IQBRC Bldg Tokai Tokai
Ibaraki 319-1106 Ibaraki 319-1106 Japan
5Molecular Photoscience Research Center Japan
5Molecular Photoscience Research Center Kobe University Kobe University Kobe
Kobe Hyogo 657-8501 Hyogo 657-8501 Japan Japan

Motivated by the recent quasi-elastic neutron scattering experiment, we extend
the spin-dynamics simulation on the triangular-lattice Heisenberg
antiferromagnet, to observe a sharp central peak of its energy width $\sim
0.001J$ ($J$ the exchange coupling) of the $Z_{2}$-vortex origin, consistently
with the experiment.

Some time ago, Kawamura and Miyashita pointed out that the frustrated
isotropic Heisenberg magnets in two dimensions could possess a novel vortex
characterized by the parity-like two-valued topological number corresponding
only to its presence/absence, a $Z_{2}$ vortex, which drove a topological
transition at a finite temperature $T=T_{V}$ associated with its binding-
unbinding [1]. Interestingly, the spin correlation length stays finite even at
and below $T_{V}$, and the low-$T$ phase is the spin paramagnetic state with
topologically broken ergodicity, called the spin-gel state [2, 3].

While the direct experimental observation of the $Z_{2}$ vortex has remained
elusive for years, a very recent quasi-elastic neutron scattering (QENS)
experiment performed on the powder sample of the quasi-two-dimensional (2D)
triangular-lattice Heisenberg antiferromagnet NaCrO2 has succeeded in directly
probing the signature of the free Z2 vortex via the observation of a sharp
quasi-elastic (QE) scattering of its energy width as narrow as $\sim 10\mu$eV
corresponding to $\sim 0.001J$ ($J$ the exchange coupling) in the finite-$T$
range of 30-40 K [4]. As this material has been considered to be a promising
candidate of the $Z_{2}$-vortex-bearing system [5], it seems fully consistent
to identify the origin of the observed sharp QE scattering as the dynamics of
the free $Z_{2}$ vortex.

The theoretical proposal that the signature of the free $Z_{2}$ vortex might
be detectable via the appearance of a sharp central peak (QE scattering) in
the dynamical spin structure factor at $T\gtrsim T_{V}$ was made earlier [6].
Indeed, by the spin-dynamics simulation on the triangular-lattice Heisenberg
model with the nearest-neighbor (nn) antiferromagnetic (AF) coupling $J$, the
appearance of such $Z_{2}$-vortex-induced central peak of the width $\sim
0.01J$ at $T>T_{V}\simeq 0.285J$ was numerically obtained [6]. The central-
peak width in the recent experiment on NaCrO2, however, was even narrower than
that in the model simulation [6] by an order of magnitude. One might wonder if
it might be possible to numerically reproduce an order-of-magnitide sharper
central peak as observed in the recent experiment. [4]

In order to examine this issue, we extend the earlier simulation of
Ref.[OkuboKawamura] to longer simulation time and to larger lattices. Longer
simulation time is implemented to improve the $\omega$-resolution required to
detect the sharp structure as a function of $\omega$, and we implement here
twenty times longer simulation time than that of Ref.[OkuboKawamura]. Larger
lattice sizes enables us to deal with more isolated, long-lived free $Z_{2}$
vortex, and we treat here the lattices of its linear size twice as large as
that treated in Ref.[OkuboKawamura].

The model considered is the same as that in Ref.[OkuboKawamura], i.e., the AF
classical Heisenberg model on the 2D triangular lattice, whose Hamiltoanin is
given by ${\cal H}=J\sum_{i,j}\bm{S}_{i}\cdot\bm{S}_{j}$ ($J>0$), where the
sum is taken over all nn pairs on the $L\times L$ triangular lattice under
periodic boundary conditions. We set here $L=768$ and 1536.

Following Ref.[OkuboKawamura], the spin dynamics simulation is performed
according to the classical analogue of the Bloch equation. The temperature
effect is taken into account via the initial spin configurations generated by
the equilibrium Monte Carlo (MC) simulation at temperature $T$ based on the
combined heat-bath and over-relaxation methods. The time evolution of the
spins is solved by the fourth-order Runge-Kutta method, where the time mesh
$\Delta t$ is taken to be 0.01 (in the $\hbar=1$ unit), commonly with
Ref.[OkuboKawamura]. The maximum simulation time $t_{max}$ is taken to be
16,000, twenty times longer than $t_{max}$ of Ref.[OkuboKawamura]. Thermal
average is taken by averaging over 200-1000 independent runs with different
spin initial conditions.

We compute the dynamical spin structure factor,

$\displaystyle S(\bm{q},\omega)$ $\displaystyle=$
$\displaystyle\langle|\bm{S}_{\bm{q}}(\omega)|^{2}\rangle,$ (1)
$\displaystyle\bm{S}_{\bm{q}}(\omega)$ $\displaystyle=$ $\displaystyle\int{\rm
d}t\sum_{i}\bm{S}_{i}(t)\exp[-i(\bm{q}\cdot\bm{r}_{i}+\omega t)],$ (2)

where $\bm{q}$ is the wavevector, $\omega$ the angular frequency, and
$\langle\cdots\rangle$ denotes the thermal average.

Typical $\omega$-dependence of the computed $S(\bm{q},\omega)$ close to the
$K$ point is shown in Fig. 1(a), where the temperature $T/J$=0.295 is slightly
above $T_{V}/J$=0.285 and the wave vector $\bm{q}$ is slightly away from the
$K$ point in the direction of the $M$ point with
$|\bm{q}-\bm{q}_{K}|=\frac{2\pi}{192}$ (in the unit of $\frac{1}{a}$, $a$
being the lattice constant). As can be seen from Fig. 1(a), a central peak is
observed in addition to the side peak originating from the damped spin waves.
The data are collected for various temperatures and wavevectors lying on the
line connecting the $K$ and $M$ points, and are fitted by the form,

$C_{L}\frac{1}{\omega^{2}+\Gamma_{L}^{2}}+C_{{\rm
DHO}}\frac{\omega_{0}^{2}\gamma}{(\omega^{2}-\omega_{0}^{2})^{2}+\Gamma^{2}\omega^{2}}+C_{0},$
(3)

where the first Lorentzian term represents the QE central peak of the
amplitude $C_{L}$ and the width $\Gamma_{L}$, the second term the contribution
of the damped harmonic oscillator (DHO) of the frequency $\omega_{0}$, the
width $\Gamma$ and the amplitude $C_{{\rm DHO}}$, the last term being the
constant background. We show in Fig. 1(a) the fitting results of the raw
$S(\bm{q},\omega)$ data where the best fitted curve is shown together with the
contribution of the three terms in Eq. (3). Reasonably good fit is obtained.

Figure 1: (Color online) (a) The $\omega$ dependence of $S(\bm{q},\omega)$
computed at a temperature $T/J=0.295$ just above $T_{V}/J=0.285$, and at the
wavevector $\bm{q}=\bm{q}_{\rm K}+\frac{1}{64}(\bm{q}_{\rm K}-\bm{q}_{\rm M})$
close to the $K$-point for the size $L=768$. The best fitted curve based on
Eq. (3) is also shown, together with the contribution of each term, i.e., the
QE central peak, the DHO term, and the constant (const) term. (b) The energy
width and (c) the integrated intensity of the fitted QE central peak plotted
versus the distance from the $K$ point along the $KM$ line,
$|\bm{q}-\bm{q}_{\rm K}|$ in the unit of $\frac{1}{a}$, illustrated in the
hexagon representing the 1st Brillouin zone shown at the right bottom.

Figures 1(b) and 1(c) exhibit the width $\Gamma_{L}$ (b) and the integrated
intensity (c) of the QE central peak plotted versus the wavevector
$|\bm{q}-\bm{q}_{K}|$ for several temeratures around $T_{V}$. As can be seen
from Fig. 1(b), the width $\Gamma_{L}$ gets narrower as $T$ approaches $T_{V}$
from above, and as $\bm{q}$ approaches $\bm{q}_{K}$. In particular,
$\Gamma_{L}$ gets as small as $\sim 0.001J$ just above $T_{V}$ in the close
vicinity of the $K$ point. The earlier estimate $\sim 0.01J$ was obtained a
bit away from the $K$ point with the lower $\omega$-resolution [6], but is
numerically consistent with our present result. One can also see from Fig.
1(c) that the QE central-peak intensity tends to decrease as $T$ approaches
$T_{V}$, eventually vanishing at $T<T_{V}$, while it gets larger as $\bm{q}$
approaches $\bm{q}_{K}$. It was already pointed out in Ref.[OkuboKawamura]
that the width of the QE central peak is quite narrow in the vicinity of the
$K$ point and is much boarder far away from the $K$ point, where the former
was associated with the free $Z_{2}$ vortices, while the latter associated
with the bound $Z_{2}$-vortex pairs. Our present observation, which reveals
the systematic variation of the width against $\bm{q}$, is fully consistent
with such a previously proposed picture.

On the basis of our present observation, we wish to discuss experimental
implications. As the recent QENS measurements on NaCrO2 was performed on the
powder sample, the measured intensity is not the contribution from a single
particular $\bm{q}$, but rather the average over the contributions from
various $\bm{q}$’s. Our result suggests that, through the powder averaging
over various $\bm{q}$’s, the resulting QE scattering intensity would become a
superposition of many Lorentzians with continuously distributed widths. Since
the free $Z_{2}$ vortex near $T_{V}$ yields the central peak of the narrowest
width of $\sim 0.001J$ with a significant intensity among nearby $\bm{q}$’s,
it would be visible even after the powder average. Indeed, the recent QENS
experiment succeeded in clearly catching the corresponding sharp QE scattering
of its width 0.01meV $\sim$ 0.001$J$ [4].

We note that, in the analysis of Ref.[Tomiyasu], the superpositions of these
continuously distributed Lorentzians are approximated by the three
Lorentzians, i.e., very sharp $L_{1}$-component of $\sim$0.01meV width
describing the free $Z_{2}$ vortex, $L_{2}$-component of $\sim$0.1meV width
describing the typical $Z_{2}$-vortex pairs, and $L_{3}$-component of
$\sim$1meV width describing the damped spin waves. The analysis, though an
approximation, is expected to capture the essential part of relevant
fluctuations. Our present analysis then suggests that the $L_{1}$-component is
borne by the contributions close to the $K$-point, while the $L_{2}$-component
borne by those away from the $K$-point. Meanwhile, if the QENS experiment
could be performed on the single crystal, more detailed information including
the $\bm{q}$ dependence might well be obtained, which would make further
detailed comparison with the theory possible. Thus, single-crystal QENS
measurements on the $Z_{2}$-vortex-bearing magnet remain to be an interesting
future task.

We are thankful to ISSP, the University of Tokyo, for providing us with CPU
time. This study was financially supported by JSPS KAKENHI (JP17H06137,
JP18K03503).

## References

  * [1] H. Kawamura and S. Miyashita, J. Phys. Soc. Jpn. 53, 4138 (1984).
  * [2] H. Kawamura, A.Yamamoto, and T. Okubo, J. Phys. Soc. Jpn. 79, 023701 (2010).
  * [3] H. Kawamura, J. Phys.: Conf. Ser. 320, 012002 (2011).
  * [4] K. Tomiyasu, Y.P. Mizuta, M. Matsuura, K. Aoyaman, and H. Kawamura, arXiv:2110.15068.
  * [5] A. Olariu, P. Mendelse, F. Bert, B.G. Ueland, P. Schiffer, R.F. Berger, and R.J. Cava, Phys. Rev. Lett. 97, 167203 (2006).
  * [6] T. Okubo and H. Kawamura, J. Phys. Soc. Jpn. 79, 084706 (2010).