11institutetext:  Institut für Astrophysik und Geophysik, Georg-August-
Universität, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany
email<EMAIL_ADDRESS>Department of Astronomy,
University of Science and Technology of China, Hefei 230026, China
33institutetext: Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany
44institutetext: Exzellenzcluster Origins, Boltzmannstraße 2, 85748 Garching,
Germany 55institutetext: Department of Physics and Astronomy, Uppsala
University, Box 516, 75120 Uppsala, Sweden 66institutetext: Thüringer
Landessternwarte Tautenburg, Sternwarte 5, 07778 Tautenburg, Germany
77institutetext: Institut de Recherche en Astrophysique et Planétologie,
Université de Toulouse, CNRS, IRAP/UMR 5277, 14 avenue Edouard Belin, F-31400,
Toulouse, France 88institutetext: Hamburger Sternwarte, Universität Hamburg,
Gojenbergsweg 112, 21029 Hamburg, Germany 99institutetext: Max-Planck-Institut
für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany
1010institutetext: European Southern Observatory, Alonso de Cordova 3107,
Vitacura, Santiago, Chile 1111institutetext: European Southern Observatory,
Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany
1212institutetext: Instituto de Astrofísica de Andalucía - CSIC, Glorieta de
la Astronomía s/n, 18008 Granada, Spain

# Retrieval of the dayside atmosphere of WASP-43b with CRIRES+

F. Lesjak 11 L. Nortmann 11 F. Yan 22 D. Cont 334411 A. Reiners 11 N. Piskunov
55 A. Hatzes 66 L. Boldt-Christmas 55 S. Czesla 66 U. Heiter 55 O. Kochukhov
55 A. Lavail 7755 E. Nagel 118866 A. D. Rains 55 M. Rengel 99 F. Rodler 1010
U. Seemann 1111 1 1 D. Shulyak 1212

(Received June 11, 2023; accepted July 17, 2023)

Accurately estimating the C/O ratio of hot Jupiter atmospheres is a promising
pathway towards understanding planet formation and migration, as well as the
formation of clouds and the overall atmospheric composition. The atmosphere of
the hot Jupiter WASP-43b has been extensively analysed using low-resolution
observations with HST and Spitzer, but these previous observations did not
cover the K band, which hosts prominent spectral features of major carbon-
bearing species such as CO and CH4. As a result, the ability to establish
precise constraints on the C/O ratio was limited. Moreover, the planet has not
been studied at high spectral resolution, which can provide insights into the
atmospheric dynamics.

In this study, we present the first high-resolution dayside spectra of
WASP-43b with the new CRIRES+ spectrograph. By observing the planet in the K
band, we successfully detected the presence of CO and provide evidence for the
existence of H2O using the cross-correlation method. This discovery represents
the first direct detection of CO in the atmosphere of WASP-43b. Furthermore,
we retrieved the temperature-pressure profile, abundances of CO and H2O, and a
super-solar C/O ratio of 0.78 by applying a Bayesian retrieval framework to
the data. Our findings also shed light on the atmospheric characteristics of
WASP-43b. We found no evidence for a cloud deck on the dayside, and recovered
a line broadening indicative of an equatorial super-rotation corresponding to
a jet with a wind speed of $\sim$ 5 km s-1, matching the results of previous
forward models and low-resolution atmospheric retrievals for this planet.

###### Key Words.:

Planets and satellites: atmospheres - techniques: spectroscopic - planets and
satellites: individuals: WASP-43b

##  1 Introduction

The study of hot Jupiters has opened up exciting opportunities and posed
significant new challenges for the analysis of exoplanetary atmospheres. These
gas giants provide an excellent opportunity for detailed observations of their
atmospheres with current telescopes and instrumentation. Due to their large
size and close proximity to the host star, they exhibit a relatively high
planet-to-star flux ratio at near-infrared wavelengths, while the short
orbital period allows for observations of full phase curves. Additionally, the
high day-night temperature difference can cause fast winds on a global scale
that redistribute heat and atmospheric constituents.

WASP-43b is a 2.0 $M_{\mathrm{Jup}}$ hot Jupiter that was discovered by
Hellier et al. (2011). The planet orbits a young K7 dwarf with a
gyrochronological age estimate of 400${}^{+200}_{-100}$ Myr (Hellier et al.
2011), which shows chromospheric and coronal emission consistent with a
correspondingly high level of stellar activity (Czesla et al. 2013; Staab et
al. 2017). Due to its short orbital period of only 0.81 d and the small
orbital distance of 0.015 AU (see Table 1), WASP-43b is thought to be tidally
locked, making it a suitable target to analyse variabilities across its
atmosphere based on phase curve observations. Previous observational data
include a low-resolution transmission spectrum and emission phase curves from
the Hubble Wide Field Camera 3 (Bean 2013), as well as two datasets in the
infrared region measured with the Spitzer InfraRed Array Camera (Blecic et al.
2014). Kreidberg et al. (2014) first found evidence for H2O absorption and
emission, while Chubb et al. (2020) recently reanalysed the data and reported
evidence for the presence of AlO being favoured over all other investigated
molecules. This points towards the existence of disequilibrium effects in the
planet’s atmosphere causing the heavy AlO molecules to be detectable in the
upper regions. While the presence of other molecules, such as CO and CH4, has
been suggested in previous studies (Kreidberg et al. 2014; Feng et al. 2016;
Stevenson et al. 2017), the results are not conclusive.

Despite WASP-43b’s ultra-short period, its equilibrium temperature of
$T_{\mathrm{eq}}\approx 1400\,$K is considerably lower than that of many other
comparable planets due to the low temperature of its host star
($T_{\mathrm{eff}}\approx 4500\,$K). This opens up the possibility of cloud
formation even on the dayside of WASP-43b, but the existence of such cloud
layers is still under investigation. While Helling et al. (2020) and Venot et
al. (2020) predict the presence of silicate and metal-oxide clouds based on 3D
atmospheric modelling, Scandariato et al. (2022) analysed the geometric albedo
and found no evidence for condensation of reflective clouds on the planetary
dayside.

Exploring the C/O ratio of exoplanet atmospheres has been identified as a
promising avenue towards unravelling the mysteries of planet formation and
migration. Öberg et al. (2011) suggested that the C/O ratio could offer clues
to the planet’s origin within the protoplanetary disk, and several subsequent
studies have explored this hypothesis (e.g. Madhusudhan et al. 2016; Mordasini
et al. 2016; Booth & Ilee 2019; Cridland et al. 2019). The composition of
solid and gaseous material varies with distance from the host star, owing to
differences in their condensation temperatures, and as the planet forms and
undergoes migration its final composition and atmospheric makeup are shaped by
the materials it accumulates along the way. Using simulations, Turrini et al.
(2021) find a significant difference in C/O based on whether the migration is
gas dominated or solid enriched, while the dependence on the formation
distance is less pronounced.

Whether the atmosphere is above or below the tipping point of C/O = 1 has
major implications for the entire composition, especially for hot and ultra-
hot Jupiters. If there is more carbon than oxygen present, the entire oxygen
supply is bound in CO and is not available to form other oxides. For more
moderate atmospheric temperatures, however, the dominant carbon-bearing
species is CH4, and thus the oxygen is free to form H2O (e.g. Giacobbe et al.
2021).

Additionally, the observable C/O is linked to the formation of clouds, as the
condensation of certain elements removes them from the gaseous phase that is
accessible with transmission and emission spectroscopy. According to Helling
et al. (2021), cloud particles such as MgO and Al2O3 can remove significant
amounts of oxygen from the local atmosphere, leading to an increase in C/O.

The previous observations of WASP-43b with the Hubble Space Telescope (HST)
and Spitzer are suitable to measure the H2O content in hot Jupiter
atmospheres, but due to the limited wavelength coverage, the determination of
other molecules with the majority of their strong emission lines in the
infrared region is more difficult. This is the case for CO and CH4, which emit
most strongly at wavelengths longer than 2$\,\mu$m. However, an estimation of
the C/O ratio necessarily relies on robust estimates of the abundances of
these main carbon-bearing molecules. Studies of one and the same planet that
are based on different observations, reduction methods, or model assumptions
can lead to vastly different C/O measurements, as was the case for example for
WASP-127b (e.g. Spake et al. 2021; Boucher et al. 2023) or the ultra-hot
Jupiter WASP-12b (e.g. Madhusudhan et al. 2011; Crossfield et al. 2012;
Stevenson et al. 2014).

In the case of WASP-43b, the results are less controversial as they generally
agree on a C/O ¡ 1. Kreidberg et al. (2014) retrieved the abundances of the
main carbon- and oxygen-bearing species using HST observations, but only the
H2O content was well constrained. Based on this result, they conclude that the
C/O ratio is consistent with the solar value. Benneke (2015) found an upper
limit of C/O ¡ 0.87 in a reanalysis of the data, while Irwin et al. (2020)
included the Spitzer/IRAC observations to retrieve C/O $\sim$ 0.91. However,
they caution that the abundances of CO and CH4 effectively only depend on the
two Spitzer datasets. The data reduction and assumed thermal profile thus have
a strong influence on the retrieved C/O ratio, such that Changeat et al.
(2021) recovered C/O = 0.68 with slightly different methods and assumptions.

Using the large simultaneous wavelength coverage and high resolution (R $\sim$
100 000) of ground-based spectrographs such as CRIRES+ mounted on the ESO
VLT/UT3 telescope, it becomes possible to circumvent dependencies on a limited
number of measurements and the ensuing ambiguities. Observing the spectral
bands of CO and CH4 in the mid-infrared promises the tightest constrains on
their abundances due to the large number and intensity of spectral lines,
consequently allowing for a more precise determination of the C/O ratio.

Table 1: Stellar and planetary parameters of the WASP-43 system.

Parameter | Symbol | Value
---|---|---
Planet |  |
Radiusa |  $R_{p}$ ($R_{\mathrm{Jup}}$) | 1.006 $\pm$ 0.017
Massa |  $M_{p}$ ($M_{\mathrm{Jup}}$) | 1.998 $\pm$ 0.079
Orbital perioda |  $P_{\mathrm{orb}}$ (days) | 0.813473978
Orbital inclinationa |  $i$ (∘) | 82.109 $\pm$ 0.088
Orbital eccentricitya | $e$ | 0
Semi-major axisa |  $a$ (AU) | 0.01504 $\pm$ 0.00029
Time of mid-transitb |  $T_{0}$ (BJD) | 2458555.80567
RV semi-amplituded |  $K_{p}$ (km s-1) | 202 $\pm$ 4
Surface gravityd |  $\log{g}$ (cgs) | 3.77
Equil. temperaturea |  $T_{\mathrm{eq}}$ (K) | 1426.7 $\pm$ 8.5
Star |  |
Radiusa |  $R_{*}$ ($R_{\odot}$) | 0.6506 $\pm$ 0.0054
Effective temperaturec |  $T_{\mathrm{eff}}$ (K) | 4400 $\pm$ 200
Systemic velocitya |  $v_{\mathrm{sys}}$ (km s-1) | -3.5955 $\pm$ 0.0043

(a)Esposito et al. (2017), (b)Patel & Espinoza (2022), (c)Bonomo et al.
(2017), (d) Calculated from the other orbital parameters.

In this work, we present high-resolution dayside observations of WASP-43b by
employing the CRIRES+ (CRyogenic InfraRed Echelle Spectrograph +, Dorn et al.
2014, 2023) instrument. This spectrograph is the upgraded version of the
original CRIRES (Kaeufl et al. 2004), offering a significantly larger
simultaneous wavelength coverage and new infrared detectors. It is mounted on
an 8 m-class telescope, which allows for high S/N observations with sufficient
time resolution. In contrast to previous observations using low-resolution
instruments, the high resolution of CRIRES+ offers the opportunity to directly
measure individual absorption and emission lines, while the capability to
observe in the infrared region enables the study of species such as CO and
CH4, which do not produce significant spectral features in the visible
wavelength range.

In Sect. 2 we describe the observations and our data reduction procedure, and
in Sect. 3 we go into the details of the model generation, cross-correlation
and retrieval analysis. We present our results and compare them to previous
observations and simulations in Sect. 4, and conclude in Sect. 5.

##  2 Observation and data reduction

Table 2: Observation log of the two nights. The variation of observational conditions over time is shown in Fig. A. Night | 1 | 2
---|---|---
Date | 2022-03-11 | 2023-02-20
Phase coverage | 0.52 - 0.66 | 0.28 - 0.49
Exposure time per frame | 300 s | 300 s
Nspectra | 34 | 48
Spectral resolution (R) |  $\sim$ 120 000 |  $\sim$ 97 000
Mean S/N per exposure | 48 | 34
Airmass | 1.43 - 1.04 | 1.03 - 1.75
| Median FWHM of the PSF
---
(in spatial direction)
0.16 $\arcsec$ | 0.33 $\arcsec$

We observed WASP-43b on 11 March 2022 and 20 February 2023 with the CRIRES+
instrument as part of the guaranteed time observations of the CRIRES+
consortium. The two observing runs were conducted shortly before and shortly
after the secondary eclipse, respectively, to observe the dayside of WASP-43b
(see Fig. 1 for the phase coverage). We used the nodding mode technique to
observe the target alternating between two different positions (A and B) along
the slit to allow for the removal of detector artefacts and the sky background
during the extraction. In our further analysis, we combined the A and B
spectra to form a single time series for each night. We observed in the K2192
wavelength setting covering parts of the wavelength region from 1972 nm to
2486 nm with an exposure time of 300 s. The CRIRES+ science detector array
consists of 3 chips, and thus each spectral order is separated into three
segments. In the analysis, we treated each segment individually, resulting in
20 wavelength segments in the chosen wavelength setting.

We employed metrology to improve the wavelength calibration by fine-tuning the
position of a defined set of emission lines from Kr- and Ne-lamps, and the
0.2$\arcsec$ slit was chosen to ensure a high spectral resolution. The first
observing run was conducted under good weather with stable seeing and relative
humidity ($\sim$30%), and adaptive optics (AO) was employed. The second night
was hampered by a varying coverage of thick and thin clouds, especially during
the first half of the observation, which prohibited the use of AO and caused a
substantial loss of flux. Although the seeing was similar in both nights, this
led to a vastly increased point spread function (PSF) during the second night.
The observational conditions that could influence the data quality are shown
in Fig. A for the two nights.

The raw spectra were reduced with the ESO CRIRES+ data reduction pipeline
(version 1.2.3) using the
EsoRex111https://www.eso.org/sci/software/cpl/esorex.html tool. This pipeline
includes the standard reduction steps for echelle time series observations,
such as dark, bias and flat field correction, removal of bad pixels, and a
calculation of the wavelength solution. In addition, we attempted to refine
the wavelength solution by using Molecfit (Smette et al. 2015) to fit the
telluric lines in the spectra, but we did not find any significant offset to
the solution provided by the pipeline. In addition, by cross-correlating the
telluric lines between individual spectra we measured the drift of the
gratings, which manifested themselves in a slight shift of the spectra in
time. The calculated drifts were consistently smaller than 0.1 pixels
(corresponding to a velocity shift of $\sim$$0.1\,$km s-1), and therefore we
did not apply any further correction.

In the case of good observing conditions that allow an excellent AO correction
and coherence time, the PSF can be smaller than the minimum slit width of 0.2
$\arcsec$, causing an increase in the spectral resolution (Dorn et al. 2023).
The median FWHM of the PSF during night 1 was 2.8 pixel (corresponding to 0.16
$\arcsec$), while in night 2 it was 5.8 pixel (0.33 $\arcsec$). To determine
the resolution, we assumed that the width of the PSF in dispersion direction
is similar to the one in spatial direction. Then we converted the PSF from
pixel space into wavelength space, which allows for the calculation of the
resolution for each segment, and found an average resolution of R $\sim$ 120
000 for night 1. The PSF in night 2 was significantly larger, and in this case
the resolution was instead constrained by the slit width to R $\sim$ 97 000.
This shows the effect of the cloud coverage and consequent lack of AO in the
second observation, resulting in a significantly broader PSF. Details of the
two observations are summarised in Table 2.

Figure 1: Orbital phase intervals of the two observations. The grey shaded
area indicates the region in which the primary and secondary eclipses occur.
The filled circles represent the planet WASP-43b, and the empty circle
represents its host star. The semi-major axis and stellar radius are to scale.

###  2.1 Normalisation

In the following, we describe the normalisation of the spectra, which consists
of removing outliers, masking any sky emission lines, and removing the
continuum variations. For the outlier removal, any data points deviating by
more than 5$\sigma$ in the time series of each pixel were removed and linearly
interpolated. If more than three values of a pixel’s time series would have to
be replaced, the entire pixel was masked for all spectra of the time series.

To mask the sky emission lines, a rough continuum approximation was computed
and all regions with a flux above 103% of the continuum were masked. To
compute such a continuum for each segment, a 90% percentile filter with a
width of 100 pixel was applied to the spectra, before smoothing the result
with a 50 pixel wide Gaussian filter.

To remove flux variations in the continuum, we determined some points
following the shape of the continuum, to which we can fit a polynomial. To
this end, we computed a master spectrum for each order as the time average of
all exposures. We divided the master spectrum into 100 wavelength bins of
equal size, and assigned to each bin the 90% percentile value of all flux
values in the bin. Subsequently, the 100 points found in this way were binned
again into only ten bins, and from each of these bins we selected the second
largest value. The ten resulting points follow the overall shape of the blaze
function while not being skewed by improperly removed emission lines. A third
degree polynomial was fitted to these ten data points, and all spectra and
errors were divided by the fit. Subsequently, each individual spectrum of each
segment was corrected to the same continuum level by division of a linear fit
to the entire spectrum. As the last normalisation step, deep telluric lines
that drop below 40% of the continuum flux were masked.

###  2.2 SYSREM

We used SYSREM (Tamuz et al. 2005; Birkby et al. 2013, see also Nugroho et al.
2017; Alonso-Floriano et al. 2019; Cabot et al. 2019; Gibson et al. 2020) for
the removal of the telluric and stellar lines from the spectra. This principal
component analysis algorithm takes into account the uncertainty of each data
point to construct a model of the data consisting of linear components in
wavelength and time. The stellar, telluric and planetary components of the
observed spectra all have different Doppler shifts. While the stellar and
telluric lines can be considered nearly static in velocity space over the
course of the observation, the planetary radial velocity changes significantly
of the order of several tens of km s-1. The SYSREM algorithm iteratively
subtracts linear trends in time from each spectral pixel, which primarily
serves to remove the (quasi-)static components of the spectra. This leaves
residuals consisting only of noise and the planetary signal, which shifts in
wavelength over time due to the planetary motion and is therefore not modelled
by SYSREM. During each iteration, the linear model created by SYSREM is
continuously refined until a convergence criterion is met, after which the
model is subtracted from the data and the computation of the next iteration
starts. We assumed the model of each iteration to be converged when the
average relative change is below 0.01.

The number of SYSREM iterations applied to the data can have a large influence
on the recovered signal strength, and we explain our process for choosing this
number in Sect. 3.3. The effect of the normalisation and application of SYSREM
on the spectra is illustrated in Fig. 2.

Figure 2: Data reduction steps of a representative wavelength range. Top
panel: Unprocessed spectra as they are produced by the CRIRES+ pipeline.
Middle panel: After continuum normalisation and masking of deep telluric
lines. Bottom panel: After the removal of stellar and telluric lines with
SYSREM.

##  3 Methods

Our analysis consisted of generating synthetic planet spectra with key
molecular species of interest, cross-correlating these with our cleaned
CRIRES+ spectra and subsequently applying a retrieval framework to constrain
the atmospheric properties.

###  3.1 Generating planetary model spectra

The synthetic model spectra for the cross-correlation analysis were computed
with the radiative transfer code petitRADTRANS (Mollière et al. 2019), using
the CO, H2O and CO2 line lists of Rothman et al. (2010), the CH4 line list of
Hargreaves et al. (2020), the HCN line list of Harris et al. (2006) and Barber
et al. (2014), and the NH3 line list of Yurchenko et al. (2011). We assumed a
non-inverted temperature-pressure ($T$-$p$) profile characterised by the
parametric model from Guillot (2010). This model is determined by the
irradiation temperature $T_{\mathrm{irr}}$, the internal temperature
$T_{\mathrm{int}}$, the infrared opacity $\kappa_{\mathrm{IR}}$ and the ratio
of visible to infrared opacity $\gamma$ (see his Eq. 29). Based on the assumed
volume mixing ratios (VMRs), we calculated the mean molecular weight and
subsequently the mass fractions of these molecules, assuming a solar-like VMR
ratio between H2 and He. We included collision-induced absorption (CIA) of
H2-H2 and H2-He (Borysow 2002; Richard et al. 2012, and the references
therein), which we found to have a substantial effect on the line strength for
this planet. For analysing the influence of clouds, we additionally added a
grey cloud deck at varying pressures, which blocks the contribution of deeper
layers. The resulting spectra were converted to a planet-to-star flux ratio.
In addition we accounted for rotational broadening of the lines, which has a
non-negligible effect for high-resolution spectroscopy. Following Eq. (3) in
Díaz et al. (2011), a rotational profile was computed based on the equatorial
rotation velocity $v_{\mathrm{eq}}$. Similar to Yan et al. (2023), we assumed
a linear limb darkening law with a coefficient $\epsilon=1$ and an inclination
angle of $\sin{i}\approx 1$ as the planet is expected to be tidally locked and
thus the equatorial plane is aligned with the orbital plane. The model was
then convolved with this profile in logarithmic wavelength space.
Subsequently, the spectra were convolved to the previously determined
resolution of CRIRES+ (Sect. 2) using a Gaussian instrumental profile, and
continuum normalised.

We computed individual synthetic spectra for the species analysed with cross-
correlation, and incorporated the model calculation in our retrieval
framework. Figure 3 shows the CO, H2O and CH4 models used for cross-
correlation, before rotational broadening is applied.

Figure 3: Synthetic model spectra of CO (top panel), H2O (middle panel) and
CH4 (bottom panel) in the wavelength range of the observations, calculated
using petitRADTRANS. The shown models are convolved to a resolution of R = 100
000 and are not rotationally broadened. The grey regions indicate the
wavelength ranges covered by the K2192 setting.

###  3.2 Cross-correlation

We calculated the weighted cross-correlation function (CCF) from the residual
spectra and the model following the work from Cont et al. (2022):

$\displaystyle\mathrm{CCF}(v,t)=\sum_{i=0}^{N}\frac{R_{i}(t)\cdot
M_{i}(v)}{\sigma_{i}(t)^{2}}\,,$ (1)

where $t$ is the time index, $i$ is the pixel index, $R$ are the residual
spectra, $M$ is the model spectrum and $\sigma$ is the uncertainty of $R$. The
model spectrum was Doppler-shifted with velocities $v$ ranging from -500 km
s-1 to +500 km s-1. This created a two-dimensional CCF map for each spectral
segment. Combining the information of all segments into a single map was done
by computing the mean CCF map. In this step we excluded some segments that
were substantially contaminated by tellurics and in which no planetary signal
could be recovered. We selected these segments by first visually inspecting
the spectra, and subsequently calculating the signal strength of an injected
model of CO and H2O as outlined in Sect. 3.3 to determine whether the signal
improves when a certain segment is excluded. Based on these criteria, we
excluded segments 1 and 20 in night 1 and segments 1, 2 and 20 in night 2. An
example of the CCF map for night 1 is shown in Fig. 4.

Following this, a $K_{p}$-$v_{\mathrm{sys}}$ map was calculated from the CCFs.
To this end, we considered a range of orbital semi-amplitude values ($K_{p}$),
and for each value we shifted the CCFs into the corresponding planetary frame
according to the Doppler velocity:

$\displaystyle
v_{p}=K_{p}\sin\,(2\pi\,\phi)+v_{\mathrm{sys}}-v_{\mathrm{bary}}+v_{\mathrm{max}}\,,$
(2)

where $\phi$ is the orbital phase, $v_{\mathrm{sys}}$ and $v_{\mathrm{bary}}$
are the systemic and barycentric velocities and $v_{\mathrm{max}}$ is a
potential deviation from the planetary rest frame. Subsequently, the shifted
CCFs were collapsed into a one-dimensional vector by averaging along the time
axis. Stacking these vectors for every trial $K_{p}$ value resulted in a two-
dimensional $K_{p}$-$v_{\mathrm{sys}}$ map, which was converted into units of
S/N by dividing by the standard deviation of the map in regions far away from
the expected signal position ($|v|>50$ km s-1).

Figure 4: Cross-correlation function of CO for night 1 after five sysrem
iterations. The dashed lines indicate the expected planetary trail after the
eclipse (which is located between the two lines). The noise structure at
$\varv=0$ km$s^{-1}$ is caused by tellurics that could not be completely
removed.

###  3.3 Selecting the number of SYSREM iterations

The removal of stellar and telluric lines with SYSREM is an iterative process,
where the choice of the number of iterations influences the final result. Too
few iterations lead to strong residuals in the $K_{p}$-$v_{\mathrm{sys}}$ map
which may hinder the detection of the planetary spectral signature. With too
many iterations, SYSREM begins to affect and remove the planetary trail as
well. In order to identify the optimal number of iterations in an objective
way, we followed the method proposed by Cheverall et al. (2023, see their
Sect. 3.8). We first applied SYSREM to the normalised data and calculated the
cross-correlation function CCFobs. Then we repeated these analysis steps on
the same data, but with a planetary model spectrum injected at the expected
Doppler velocity. This yielded the signal-injected cross-correlation function
CCFinj. Subsequently, we derived the differential cross-correlation function
$\Delta$CCF = CCFinj \- CCFobs, and calculated the $K_{p}$-$v_{\mathrm{sys}}$
map and the S/N for each iteration. Finally, we chose the SYSREM iteration
that maximises the S/N of $\Delta$CCF, and applied the same number of
iterations to the non-injected data. The results of this selection process are
summarised in Table 3 and Fig. 5.

Figure 5: S/N detection strength as a function of SYSREM iterations for the
data of night 1. The measured signal strengths are a result of using a pure CO
model (blue colour), a pure H2O model (orange colour), and a pure CH4 model
(black colour). The top panel shows the signal strength for the real data,
while the bottom panel shows the strength of an injected signal using the
differential $\Delta$CCF according to the method described in Sect. 3.3. The
diamond shape indicates the SYSREM iteration with the strongest detection of
the injected signal.

###  3.4 Retrieval framework

The properties of the atmosphere of WASP-43b were analysed by performing an
atmospheric retrieval, in which the observed data is compared to a large
number of different models. Using a sampling algorithm to determine the
posterior distributions of the model parameters allows for the statistical
measurement of atmospheric characteristics. In this way, any properties that
have an effect on the observed emission spectra can be analysed. Both the
$T$-$p$ profile and the abundances influence the line strengths, while the
position of the lines over the course of the observing run depends on the
$K_{p}$ and systemic velocity. Additionally, the equatorial velocity
$v_{\mathrm{eq}}$ can be determined as a rotation increases the line width by
means of rotational broadening. The existence of clouds at a certain pressure
Pcloud truncates the spectral lines across the entire wavelength range as
contributions from regions below the cloud deck are blocked.

For the retrieval, we assumed that the abundances are isobaric and that the
temperature is described along the parameterisation of Guillot (2010). While
this profile does not have enough degrees of freedom to capture the
temperature variations across the entire height of the atmosphere, we expect
it to be sufficient to model the region with the steepest temperature
gradient, from which the majority of the signal originates. We did not include
Rayleigh scattering as this process did not have any appreciable effect on the
models.

We followed the method of Yan et al. (2023), which is based on the work of Yan
et al. (2020) and inspired by Brogi & Line (2019), Shulyak et al. (2019) and
Gibson et al. (2020). The forward models were computed as described in Sect.
3.1. To be able to compare these models to the data, it is important to
consider any factors that impact the planetary signal. Therefore we modelled
the effect of smeared spectral lines and applied corrections to account for
distortions introduced by the application of SYSREM.

####  3.4.1 Accounting for smearing of the signal

Due to its short orbital period, the radial velocity of WASP-43b changes
significantly during a single exposure of 300 s. For night 1, the velocity
difference between start and end of an exposure was between 2 - 5 km s-1,
causing a substantial smearing of the spectral lines. To account for this
effect, the model was interpolated to the same wavelength grid as the 2D data
(wavelength $\times$ time), resulting in a model matrix $M$. The model could
then be shifted according to the planetary Doppler velocity (see Eq. 2). For
each exposure, we calculated these radial velocities at the beginning and end
time, shifted the model to ten evenly spaced planetary rest-frames between
these two velocities, and computed the smeared model as the mean over the ten
sub-exposures. The effect of smearing on the shape of the spectral lines is
shown in Fig. 6. For a model without rotational broadening, the smearing
changes the line shape significantly, while in the case of a strongly
rotationally broadened model, the smearing is barely noticeable.

Figure 6: Profile of a single CO line for models with different combinations
of smearing and rotational broadening. All models have the instrumental
broadening applied, and were additionally smeared with a velocity change of 5
km s-1 (orange), rotationally broadened with $v_{\mathrm{eq}}=11$ km s-1
(purple), or had both effects combined (brown).

####  3.4.2 Filtering the models

Applying SYSREM not only removes static signals in the data, but also distorts
the remaining planetary signal. Therefore, it is important to account for this
distortion effect when comparing the planetary signal with synthetic model
spectra via Bayesian inference algorithms to retrieve precise atmospheric
properties. One could apply the same detrending procedure to both data and
model, but this is computationally expensive. Instead, Gibson et al. (2022)
proposed a fast alternative filtering technique using just the output of
SYSREM, which is briefly summarised in the following:

During each iteration, SYSREM calculates column and row vectors u and w as a
best approximation of the data. Combining all column vectors into a matrix U
(with the dimensions $N_{\mathrm{exposures}}$ $\times$
$N_{\mathrm{iterations}}$), the correction for a model
$M_{\mathrm{unfiltered}}$ can be computed:

$\displaystyle M^{\prime}=U(\Lambda U)^{\dagger}(\Lambda
M_{\mathrm{unfiltered}})\,,$ (3)

where $A^{\dagger}=(A^{T}A)^{-1}A^{T}$ is the Moore-Penrose inverse of a
matrix $A$. Each entry in the diagonal matrix $\Lambda$ represents the
reciprocal mean uncertainty in wavelength direction of the corresponding
exposure. We precomputed this filter and applied it to every model during the
retrieval, yielding the filtered model $M_{\mathrm{unfiltered}}-M^{\prime}$
which was then compared to the data.

####  3.4.3 Likelihood function

To retrieve the atmospheric properties, we applied a Markov chain Monte Carlo
(MCMC) algorithm using the EMCEE package (Foreman-Mackey et al. 2013) to the
logarithm of the likelihood function (Hogg et al. 2010; Yan et al. 2020):

$\displaystyle\ln(L)=-\frac{1}{2}\sum_{i,j}\left(\frac{(R_{i,j}-M_{i,j})^{2}}{(\beta\sigma_{i,j})^{2}}+\ln(2\pi(\beta\sigma_{i,j})^{2})\right)\,,$
(4)

where $R_{i,j}$ is the residual spectrum at pixel index $i$ and time $j$,
$M_{i,j}$ is the 2D matrix of a filtered model spectrum, $\sigma_{i,j}$ are
the uncertainties of the residuals and $\beta$ is a scaling factor to
$\sigma_{i,j}$. For the MCMC, we used 32 walkers with 15 000 steps each. We
chose uniform priors for all input parameters, which are summarised in Table
4.

##  4 Results

###  4.1 Cross-correlation results

We calculated $K_{p}$-$v_{\mathrm{sys}}$ maps for CO, H2O and CH4 as detailed
in Sect. 3.2. The models used for the cross-correlation were computed assuming
abundances of $\log_{10}$ CO = $-4$, $\log_{10}$ H2O = $-4$ and $\log_{10}$
CH4 = $-4$, respectively. We calculated the temperature-pressure profile using
the parameterisation of Guillot (2010) with $T_{\mathrm{irr}}=1400\,$K,
$T_{\mathrm{int}}=800\,$K, log${}_{10}(\kappa_{\mathrm{ir}})$= $-2$ and
log${}_{10}(\gamma)=-0.4$, and calculated individual model spectra for these
species. We only included instrumental broadening and did not account for any
other broadening effects for these models. The cross-correlation analysis was
conducted separately for each night. Table 3 summarises the resulting S/N of
the signal in the $K_{p}$-$v_{\mathrm{sys}}$ maps, as well as the number of
SYSREM iterations, determined as outlined in Sect. 3.3. Figures 7 and 12 show
the $K_{p}$-$v_{\mathrm{sys}}$ maps for CO and H2O of night 1 and night 2,
respectively, while the results for CH4 are shown in Fig. 13.

Table 3: S/N of the signal in the $K_{p}$-$v_{\mathrm{sys}}$ map. | Night | CO | H2O | CH4
---|---|---|---|---
S/N | 1 | 5.02 | 4.03 | 2.14
| 2 | 0.82 | 2.47 | 2.46
# iterations | 1 | 5 | 6 | 10
| 2 | 8 | 8 | 8

We detected the presence of CO and evidence for H2O in the first night.
However, we were not able to find either of the species in night 2. This does
not come as a surprise, because the data quality of night 2 is considerably
worse (see Sect. 2). We assessed the quality of night 2 by injecting a signal
with the expected strength, and were not able to recover it. Therefore the
non-detection in night 2 is most likely due to the poor data quality. Pino et
al. (2022) showed that principal component analysis algorithms remove much of
the planetary signal at phases close to quadrature, as there is only a small
shift in the line positions. This could affect especially the signal in night
2, as the observation started at a phase of 0.28. Using an injection-recovery
test, we determined that SYSREM removes a substantial amount of the planet
signal at phases up to 0.4. We ran the analysis again and excluded any
exposures prior to this point in orbit, but we were still not able to detect
CO or H2O in night 2.

We found no evidence for CH4 in either of the nights. In addition we searched
for other molecules that are expected to be present based on equilibrium
chemistry calculations (CO2, HCN and NH3), but we found no evidence for any of
these species in either night. We note that there are spurious peaks which
have a S/N of 3-4 visible in some of the $K_{p}$-$v_{\mathrm{sys}}$ maps.
However, these are far away from the expected position of a real signal, and
are generally not stable against small changes in the model. We tested
different models with variations in the $T$-$p$ profile or the inclusion of
rotational broadening, and found that the strength and position of these
spurious peaks is not consistent. On the other hand, the CO and H2O signal in
night 1 remain at the same position and at a stable S/N for all of the tested
models. This shows that one has to be careful in the interpretation of signals
with S/N $\sim$ 4, as these can in principle be created by spurious noise
structures, and the stability of such signals against changes in the model
should be investigated.

The results of night 1 agree with earlier studies of WASP-43b using HST and
Spitzer observations. Since the first detection of H2O by Kreidberg et al.
(2014), water features have been continuously reconfirmed in subsequent
analyses (Stevenson et al. 2017; Chubb et al. 2020). While these studies found
no conclusive evidence for CO absorption or emission lines, this is not
surprising as the K band, which hosts the most prominent spectral bands of CO,
was not covered by these previous observations.

Figure 7: $K_{p}$-$v_{\mathrm{sys}}$ map for the CO signal after five SYSREM
iterations (top panel) and H2O signal after six iterations (bottom panel) in
night 1. The white dotted lines indicate the expected position of a signal at
the literature $K_{p}$ value in the planetary rest-frame. The CO signal has a
S/N of 5.02, and the H2O signal a S/N of 4.03.  Figure 8: Posterior
distributions of the atmospheric and orbital parameters for WASP-43b using the
data of night 1. The top right panel shows the median $T$-$p$ profile, which
was computed from 10 000 random samples of the MCMC walkers. The blue shaded
region indicates the 1$\sigma$ uncertainty interval, while the vertical dotted
line shows the equilibrium temperature of WASP-43b
($T_{\mathrm{eq}}=1426.7\,$K). The mean contribution function for the
wavelength range of our observation is shown in grey.

###  4.2 Retrieval results

Table 4: Priors and retrieved values of the retrieval using only the data from night 1. Parameter | Prior | Retrieved value | Unit
---|---|---|---
$T_{\mathrm{int}}$ | [100, 3000] | 1100${}^{+600}_{-700}$ | K
$T_{\mathrm{irr}}$ | [300, 3000] | 1100 $\pm$ 500 | K
$\log_{10}$ $\kappa_{\mathrm{IR}}$ | [-15, 4] | $-1.7^{+1.4}_{-0.6}$ |  $\log_{10}$ cm2 s-1
$\log_{10}$ $\gamma$ | [-2, 2] | $-0.6^{+1.4}_{-1.0}$ | …
$\log_{10}$ CO | [-10, 0] | $-2.8^{+1.4}_{-0.8}$ | …
$\log_{10}$ H2O | [-10, 0] | $-4.3^{+1.0}_{-0.4}$ | …
$\log_{10}$ Pcloud | [-8, 2] | $0.5^{+1.0}_{-1.4}$ |  $\log_{10}$ bar
$K_{p}$ | [170, 230] | $188^{+5}_{-6}$ | km s-1
$v_{\mathrm{max}}$ | [-20, 20] | $-2\pm 3$ | km s-1
$v_{\mathrm{eq}}$ | [0.1, 20] | $11.7^{+2.8}_{-2.1}$ | km s-1
$\beta$ | [0.2, 5] | $0.6807\pm 0.0005$ | …

While the cross-correlation method can confirm the presence of a species, it
does not provide robust estimates of the abundances. We therefore applied the
retrieval framework detailed in Sect. 3.4 to constrain the atmospheric
properties. As we did not find evidence for a signal of CO or H2O (or any of
the other species) in night 2, we applied the retrieval only to the data of
night 1.

Following a 15 000 step MCMC sampling, we discarded the initial 5000 steps as
burn-in, and subsequently calculated the posterior distributions of the
atmospheric and orbital parameters, which are illustrated in Fig. 8. The
retrieved parameters and their uncertainties are given in Table 4. We
retrieved a non-inverted $T$-$p$ profile (see Fig. 8), which agrees in the
upper atmospheric layers with the equilibrium temperature of
$T_{\mathrm{eq}}=1426.7\pm 8.5\,$K. We found a planet velocity semi-amplitude
$K_{p}=188^{+5}_{-6}\,$km s-1 and an offset in velocity from the planetary
rest frame of $v_{\mathrm{max}}=-2\pm 3\,$km s-1. The retrieved $K_{p}$ value
is slightly below the literature value of $K_{p,\mathrm{lit}}=202\pm 4\,$km
s-1, which can be explained by the elongated diagonal structure of the H2O
detection peak in the $K_{p}$-$v_{\mathrm{sys}}$ map and the dependency of
$K_{p,\mathrm{lit}}$ on the semi-major axis, which is not precisely known.

Figure 9: Posterior distributions of the CO abundance (top panels) and H2O
abundance (bottom panels). The left panels show the correlation between the
VMRs and the cloud deck pressure, and the right panels show the total
distribution of the VMRs divided into two subsets: with $\log_{10}$ Pcloud ¿ 0
(blue and red) and with $\log_{10}$ Pcloud ¡ 0 (green and purple). The
threshold (indicated by the dashed line) is caused by collision-induced
absorption, which blocks any signals originating in the deeper layers. Figure
10: Line depth of the deepest line in a model with CO and H2O as a function of
cloud deck pressure, for a model with collision-induced absorption (orange)
and without this effect (blue).

The inclusion of a cloud deck leads to a muting of the signal, as
contributions of layers below the cloud deck pressure $P_{\mathrm{cloud}}$
become inaccessible to observations. This introduces a degeneracy between the
cloud height and abundances, as a muting of the signal due to clouds can be
counterbalanced by increased volume mixing ratios. Figure 9 shows the
correlation between $P_{\mathrm{cloud}}$ and the VMRs of CO and H2O, and the
strong degeneracy for cloud deck pressures below one bar is clearly visible.
Above this threshold, the VMRs remain constant with decreasing cloud height,
because collision-induced absorption (CIA) of H2-H2 and H2-He acts as an
additional continuum opacity source, preventing us from probing deeper into
the atmosphere. This effect is illustrated in Fig. 10, showing the relation
between line strength and cloud deck pressure for models with and without CIA.
At pressures of one bar, the line strength of the model with CIA levels off,
and a further increase of $P_{\mathrm{cloud}}$ does not affect the model.

Previous retrievals of low-resolution observations found no evidence for the
existence of clouds on the dayside of WASP-43b. Fraine et al. (2021) used HST
WFC3/UVIS to study the planet’s reflected light and were not able to detect
the eclipse, leading them to conclude that a significant cloud coverage is
unlikely to exist above the probed pressures of around one bar. Phase curve
observations with CHEOPS, TESS and HST constrained the maximum cloud deck
pressure to a similar value (Scandariato et al. 2022), and the retrieval of
Chubb & Min (2022) using Spitzer data also found a good fit using cloud-free
models. Murphy et al. (2023) found that, when comparing their models to the
Spitzer phase curve, either thin or very thick clouds were favoured. Although
they are unable to reject either of the two scenarios due to a limited model
resolution and degeneracies, they suggest a relatively cloud-free dayside to
be most likely.

Contrary to these results, pure general circulation and cloud formation models
of WASP-43b’s atmosphere by Venot et al. (2020) and Helling et al. (2020)
predict several species of clouds at pressures above the detection limit of
the aforementioned observational studies, although these results could
strongly depend on the initial assumptions and fidelity of the simulations.

Our retrieval results favour an atmosphere free of clouds down to pressures of
at least one bar. In addition, due to the unphysically large CO volume mixing
ratios when assuming larger cloud deck pressures, and the absence of
observational evidence for clouds in the low-resolution data, we deem it most
likely that clouds are either not present or located below this threshold
pressure on the dayside. Constraining the posterior distributions to include
only MCMC steps with $\log_{10}P_{\mathrm{cloud}}<0$, we retrieve a super-
solar carbon monoxide abundance of $\log_{10}$(CO) = $-3.0^{+0.8}_{-0.6}$, in
agreement with the findings of Irwin et al. (2020) and Chubb & Min (2022). If
the H2O signal is real and not created by spurious noise, the detected feature
corresponds to a water abundance of $\log_{10}$(H2O) = $-4.4\pm$ 0.3. This H2O
abundance is consistent within $\sim 1\sigma$ with the results of Stevenson et
al. (2017) and Welbanks & Madhusudhan (2019, 2022), but lower than the
detections of Kreidberg et al. (2014) and Irwin et al. (2020).

We found a broadening of the spectral lines corresponding to an equatorial
velocity $v_{\mathrm{eq}}=11.7^{+2.8}_{-2.1}\,$km s-1, which is larger than
the rotational velocity $v_{\mathrm{rot}}=2\pi R_{p}/P=6.29\pm 0.11\,$km s-1
assuming the planet is tidally locked. The difference in velocity can be
explained by a high velocity jet causing equatorial super-rotation with a wind
speed of $v_{\mathrm{jet}}=5.4^{+2.8}_{-2.1}\,$km s-1. For the rotational
profile, we assumed a limb darkening coefficient of $\epsilon=1$, which is
justified by the existence of a hot spot creating a strong contrast between
the centre and limb region. Setting this parameter to significantly smaller
values does not influence the results by more than 1-2 km s-1.

The existence of such a jet is supported both by previous observations and
simulations. Multiple studies of Spitzer observations found an eastward shift
of the hot spot by $\sim$ 10∘ (e.g. Stevenson et al. 2017; Chubb & Min 2022;
Murphy et al. 2023), while Scandariato et al. (2022) analysed TESS phase
curves and found marginal evidence for a $\left(50^{+30}_{-20}\right)^{\circ}$
eastward phase offset. The formation of prograde winds at the equator is
generally predicted by 3D general circulation models of tidally locked hot
Jupiters (e.g. Showman et al. 2009; Rauscher & Menou 2010; Heng et al. 2011;
Mayne et al. 2014). In simulations specifically of WASP-43b, Kataria et al.
(2015) find that the jet can in principle reach velocities upwards of 4 km
s-1, while Zhang et al. (2017) report a wind speed of 6.3 km s-1 at pressure
levels associated with CO lines. There has been the notion that both prograde
and retrograde streams exist in different regions and heights, and that the
retrograde equatorial flow could dominate the overall dynamics (Carone et al.
2020). However, the follow-up study by Schneider et al. (2022) suggests that
the retrograde flow may be an artefact of imposed temperature forcing. In
their refined model, the prograde jet emerges as the dominant structure with a
velocity of $\sim$ 5 km s-1. Our retrieved equatorial wind speed thus agrees
well with previous observations and various simulations of the atmospheric
dynamics.

We investigated the influence of smeared spectral lines (see Sect. 3.4.1) on
the retrieved $v_{\mathrm{eq}}$, but we found no significant difference
between retrievals with and without smearing. As shown in Fig. 6, the addition
of smearing is relevant in the absence of rotational broadening, but only has
a marginal effect on the shape of the lines when a strong rotational
broadening is present.

###  4.3 Determining the C/O ratio

We conducted a second retrieval, in which the C/O ratio and metallicity [M/H]
were set as free parameters instead of the abundances of CO and H2O. We used
the poor_mans_nonequ_chem subpackage from petitRADTRANS to calculate chemical
equilibrium abundances from C/O ratios and metallicities. Because we concluded
previously that there is no evidence for clouds in the probed region of the
atmosphere, we did not include a cloud deck in this second retrieval.

We found C/O = 0.78 $\pm$ 0.09 and [M/H] = 0.2${}^{+0.7}_{-0.6}$, while the
other parameters ($T$-$p$ profile, $K_{p}$, $v_{\mathrm{max}}$,
$v_{\mathrm{eq}}$, $\beta$) were consistent with our initial retrieval with
free abundances (see Fig. C). Calculating the abundances of CO and H2O based
on this C/O ratio and metallicity also gave VMRs close to the results of our
initial retrieval. According to these calculations, there could be a
significant amount of CH4 present ($\log_{10}$(CH${}_{4})\approx-3.2$) given
the retrieved atmospheric conditions, which we were not able to detect with
our data. We conducted an injection-recovery test on the data of night 1 with
a pure CH4 model of that abundance, and were not able to recover any signal.
The line strength of the model had to be increased by a factor of 3 to be able
to tentatively detect the injected signal in the data. Therefore CH4 could be
present in the atmosphere with an abundance consistent with our retrieved C/O
ratio, without being detectable in the data.

Our result is in agreement with the general consensus that the gaseous C/O
ratio is significantly smaller than unity (Benneke 2015; Changeat et al. 2021;
Chubb & Min 2022). The retrieved super-solar metallicity matches with findings
of Kreidberg et al. (2014, [M/H] = 0.4 - 3.5 $\times$ solar) and Stevenson et
al. (2017, 0.4 - 1.7 $\times$ solar), while it is lower than the estimates of
Kataria et al. (2015, 5 $\times$ solar) and Changeat et al. (2021, [M/H] =
1.81).

Many previous studies used the Spitzer datasets to constrain the C/O ratio,
and these results show a strong dependency on the reduction method and assumed
thermal profiles (Changeat et al. 2021). The observations presented in this
work cover the CO band head in the K band with the high resolution of CRIRES+.
Therefore the determination of the CO abundance does not rely on individual
data points, as it is informed by a large number of CO absorption lines.
Nevertheless, similar to previous works the retrieved C/O ratio is only based
on two molecular species, as we were not able to detect other carbon- and
oxygen-bearing molecules. This could be either due to small overall
atmospheric abundances, or because the observations were not sensitive to the
relevant wavelength or pressure regions for these species.

WASP-43b is a target of the JWST Cycle 1 ERS Program (DD-ERS 136, PIs: N.
Batalha, J. Bean, K. Stevenson; Stevenson et al. 2016; Bean et al. 2018),
which observed a full phase curve of the planet with JWST/MIRI in the infrared
5-12$\,\mu$m range (Venot et al. 2020). This observation will be used to study
the planetary night side, and to constrain abundances and properties of
clouds. A joint retrieval using these not yet published observations with JWST
will offer the best bet for an accurate determination of the abundances and
the C/O ratio.

##  5 Conclusions

We analysed the thermal emission spectra of WASP-43b observed with the high-
resolution spectrograph CRIRES+ by implementing the cross-correlation method.
We found a robust CO absorption signal, the first direct detection of this
species in WASP-43b’s atmosphere, as well as evidence for H2O absorption
lines. For this, we removed stellar and telluric contributions using SYSREM,
and optimised the number of iterations based on an injection-recovery test.
Subsequently, forward-model retrievals were conducted to constrain the
atmospheric and orbital parameters of WASP-43b. Our retrieval accounts for the
smearing of spectral lines due to the planetary motion during exposures,
applies a filter to the model to account for distortion effects introduced by
SYSREM, and includes rotational broadening due to the tidally locked rotation
and equatorial winds.

We recovered a non-inverted $T$-$p$ profile that agrees well with the
equilibrium temperature in the upper regions of the atmosphere. The retrieved
CO and H2O abundances are generally consistent with results of previous
studies, although there are some studies that report a notably greater H2O
VMR. Assuming equilibrium chemistry in the atmosphere, we found a super-solar
metallicity and a C/O ratio of 0.78 $\pm$ 0.09. This result is based on the
detected signal of CO and H2O only, and does not account for any other carbon-
and oxygen-bearing species that were not detectable with our data. In
agreement with simulations of tidally locked hot Jupiters, we found a
rotational broadening consistent with a fast super-rotating equatorial jet
with a speed of $v_{\mathrm{jet}}=5.4^{+2.8}_{-2.1}$ km s-1. Furthermore,
there was no evidence for a significant cloud deck coverage on the dayside
above pressures of one bar based on our data. A combination of these results
with the JWST observations of WASP-43b will lead to a larger wavelength
coverage and more robust abundance measurements not only of CO and H2O, but
potentially also of other carbon- and oxygen-bearing species.

###### Acknowledgements.

CRIRES+ is an ESO upgrade project carried out by Thüringer Landessternwarte
Tautenburg, Georg-August Universität Göttingen, and Uppsala University. The
project is funded by the Federal Ministry of Education and Research (Germany)
through Grants 05A11MG3, 05A14MG4, 05A17MG2 and the Knut and Alice Wallenberg
Foundation. Based on observations collected at the European Organisation for
Astronomical Research in the Southern Hemisphere under ESO programmes
108.22PH.004 and 110.249Y.004. F.L. acknowledges the support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – Project number
314665159. D.C. is supported by the LMU-Munich Fraunhofer-Schwarzschild
Fellowship and by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany´s Excellence Strategy – EXC 2094 – 390783311. M.R.
and S.C. acknowledge the support by the DFG priority program SPP 1992
“Exploring the Diversity of Extrasolar Planets” (DFG PR 36 24602/41 and CZ
222/5-1, respectively). D.S. acknowledges funding from project
PID2021-126365NB-C21(MCI/AEI/FEDER, UE) and financial support from the grant
CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033.

## References

  * Alonso-Floriano et al. (2019) Alonso-Floriano, F. J., Sánchez-López, A., Snellen, I. A. G., et al. 2019, Astronomy & Astrophysics, 621, A74
  * Barber et al. (2014) Barber, R. J., Strange, J. K., Hill, C., et al. 2014, Monthly Notices of the Royal Astronomical Society, 437, 1828
  * Bean (2013) Bean, J. 2013, 13467
  * Bean et al. (2018) Bean, J. L., Stevenson, K. B., Batalha, N. M., et al. 2018, Publications of the Astronomical Society of the Pacific, 130, 114402
  * Benneke (2015) Benneke, B. 2015, arXiv e-prints, arXiv:1504.07655
  * Birkby et al. (2013) Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, Monthly Notices of the Royal Astronomical Society: Letters, 436, L35
  * Blecic et al. (2014) Blecic, J., Harrington, J., Madhusudhan, N., et al. 2014, The Astrophysical Journal, 781, 116
  * Bonomo et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, Astronomy & Astrophysics, 602, A107
  * Booth & Ilee (2019) Booth, R. A. & Ilee, J. D. 2019, Monthly Notices of the Royal Astronomical Society, 487, 3998
  * Borysow (2002) Borysow, A. 2002, Astronomy and Astrophysics, 390, 779
  * Boucher et al. (2023) Boucher, A., Lafreniére, D., Pelletier, S., et al. 2023, MNRAS, 522, 5062
  * Brogi & Line (2019) Brogi, M. & Line, M. R. 2019, The Astronomical Journal, 157, 114
  * Cabot et al. (2019) Cabot, S. H. C., Madhusudhan, N., Hawker, G. A., & Gandhi, S. 2019, Monthly Notices of the Royal Astronomical Society, 482, 4422
  * Carone et al. (2020) Carone, L., Baeyens, R., Mollière, P., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 3582
  * Changeat et al. (2021) Changeat, Q., Al-Refaie, A. F., Edwards, B., Waldmann, I. P., & Tinetti, G. 2021, The Astrophysical Journal, 913, 73
  * Cheverall et al. (2023) Cheverall, C. J., Madhusudhan, N., & Holmberg, M. 2023, Monthly Notices of the Royal Astronomical Society, 522, 661
  * Chubb & Min (2022) Chubb, K. L. & Min, M. 2022, Astronomy & Astrophysics, 665, A2
  * Chubb et al. (2020) Chubb, K. L., Min, M., Kawashima, Y., Helling, C., & Waldmann, I. 2020, Astronomy & Astrophysics, 639, A3
  * Cont et al. (2022) Cont, D., Yan, F., Reiners, A., et al. 2022, A&A, 668, A53
  * Cridland et al. (2019) Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2019, Astronomy & Astrophysics, 632, A63
  * Crossfield et al. (2012) Crossfield, I. J. M., Barman, T., Hansen, B. M. S., Tanaka, I., & Kodama, T. 2012, The Astrophysical Journal, 760, 140
  * Czesla et al. (2013) Czesla, S., Salz, M., Schneider, P. C., & Schmitt, J. H. M. M. 2013, Astronomy and Astrophysics, 560, A17
  * Díaz et al. (2011) Díaz, C. G., González, J. F., Levato, H., & Grosso, M. 2011, Astronomy & Astrophysics, 531, A143
  * Dorn et al. (2014) Dorn, R. J., Anglada-Escude, G., Baade, D., et al. 2014, The Messenger, 156, 7
  * Dorn et al. (2023) Dorn, R. J., Bristow, P., Smoker, J. V., et al. 2023, Astronomy & Astrophysics, 671, A24
  * Esposito et al. (2017) Esposito, M., Covino, E., Desidera, S., et al. 2017, Astronomy & Astrophysics, 601, A53
  * Feng et al. (2016) Feng, Y. K., Line, M. R., Fortney, J. J., et al. 2016, The Astrophysical Journal, 829, 52
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Fraine et al. (2021) Fraine, J., Mayorga, L. C., Stevenson, K. B., et al. 2021, The Astronomical Journal, 161, 269
  * Giacobbe et al. (2021) Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205
  * Gibson et al. (2020) Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 2215
  * Gibson et al. (2022) Gibson, N. P., Nugroho, S. K., Lothringer, J., Maguire, C., & Sing, D. K. 2022, MNRAS, 512, 4618
  * Guillot (2010) Guillot, T. 2010, Astronomy and Astrophysics, Volume 520, id.A27, 13 pp., 520, A27
  * Hargreaves et al. (2020) Hargreaves, R. J., Gordon, I. E., Rey, M., et al. 2020, The Astrophysical Journal Supplement Series, 247, 55
  * Harris et al. (2006) Harris, G. J., Tennyson, J., Kaminsky, B. M., Pavlenko, Ya. V., & Jones, H. R. A. 2006, Monthly Notices of the Royal Astronomical Society, 367, 400
  * Hellier et al. (2011) Hellier, C., Anderson, D. R., Cameron, A. C., et al. 2011, Astronomy & Astrophysics, 535, L7
  * Helling et al. (2020) Helling, C., Kawashima, Y., Graham, V., et al. 2020, Astronomy & Astrophysics, 641, A178
  * Helling et al. (2021) Helling, C., Lewis, D., Samra, D., et al. 2021, Astronomy & Astrophysics, 649, A44
  * Heng et al. (2011) Heng, K., Menou, K., & Phillipps, P. J. 2011, Monthly Notices of the Royal Astronomical Society, 413, 2380
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv e-prints, arXiv:1008.4686
  * Irwin et al. (2020) Irwin, P. G. J., Parmentier, V., Taylor, J., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 106
  * Kaeufl et al. (2004) Kaeufl, H.-U., Ballester, P., Biereichel, P., et al. 2004, 5492, 1218
  * Kataria et al. (2015) Kataria, T., Showman, A. P., Fortney, J. J., et al. 2015, The Astrophysical Journal, 801, 86
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, The Astrophysical Journal, 793, L27
  * Madhusudhan et al. (2016) Madhusudhan, N., Agúndez, M., Moses, J. I., & Hu, Y. 2016, Space Science Reviews, 205, 285
  * Madhusudhan et al. (2011) Madhusudhan, N., Harrington, J., Stevenson, K. B., et al. 2011, Nature, 469, 64
  * Mayne et al. (2014) Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014, Astronomy and Astrophysics, 561, A1
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, Astronomy & Astrophysics, 627, A67
  * Mordasini et al. (2016) Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, The Astrophysical Journal, 832, 41
  * Murphy et al. (2023) Murphy, M. M., Beatty, T. G., Roman, M. T., et al. 2023, The Astronomical Journal, 165, 107
  * Nugroho et al. (2017) Nugroho, S. K., Kawahara, H., Masuda, K., et al. 2017, The Astronomical Journal, 154, 221
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, The Astrophysical Journal, 743, L16
  * Patel & Espinoza (2022) Patel, J. A. & Espinoza, N. 2022, The Astronomical Journal, 163, 228
  * Pino et al. (2022) Pino, L., Brogi, M., Désert, J. M., et al. 2022, A&A, 668, A176
  * Rauscher & Menou (2010) Rauscher, E. & Menou, K. 2010, The Astrophysical Journal, 714, 1334
  * Richard et al. (2012) Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, Journal of Quantitative Spectroscopy and Radiative Transfer, 113, 1276
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, Journal of Quantitative Spectroscopy and Radiative Transfer, 111, 2139
  * Scandariato et al. (2022) Scandariato, G., Singh, V., Kitzmann, D., et al. 2022, A&A, 668, A17
  * Schneider et al. (2022) Schneider, A. D., Carone, L., Decin, L., et al. 2022, A&A, 664, A56
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, The Astrophysical Journal, 699, 564
  * Shulyak et al. (2019) Shulyak, D., Rengel, M., Reiners, A., Seemann, U., & Yan, F. 2019, Astronomy & Astrophysics, 629, A109
  * Smette et al. (2015) Smette, A., Sana, H., Noll, S., et al. 2015, Astronomy & Astrophysics, 576, A77
  * Spake et al. (2021) Spake, J. J., Sing, D. K., Wakeford, H. R., et al. 2021, MNRAS, 500, 4042
  * Staab et al. (2017) Staab, D., Haswell, C. A., Smith, G. D., et al. 2017, Monthly Notices of the Royal Astronomical Society, 466, 738
  * Stevenson et al. (2014) Stevenson, K. B., Bean, J. L., Madhusudhan, N., & Harrington, J. 2014, The Astrophysical Journal, 791, 36
  * Stevenson et al. (2016) Stevenson, K. B., Lewis, N. K., Bean, J. L., et al. 2016, Publications of the Astronomical Society of the Pacific, 128, 094401
  * Stevenson et al. (2017) Stevenson, K. B., Line, M. R., Bean, J. L., et al. 2017, The Astronomical Journal, 153, 68
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, Monthly Notices of the Royal Astronomical Society, 356, 1466
  * Turrini et al. (2021) Turrini, D., Schisano, E., Fonte, S., et al. 2021, The Astrophysical Journal, 909, 40
  * Venot et al. (2020) Venot, O., Parmentier, V., Blecic, J., et al. 2020, The Astrophysical Journal, 890, 176
  * Welbanks & Madhusudhan (2019) Welbanks, L. & Madhusudhan, N. 2019, The Astronomical Journal, 157, 206
  * Welbanks & Madhusudhan (2022) Welbanks, L. & Madhusudhan, N. 2022, The Astrophysical Journal, 933, 79
  * Yan et al. (2023) Yan, F., Nortmann, L., Reiners, A., et al. 2023, A&A, 672, A107
  * Yan et al. (2020) Yan, F., Pallé, E., Reiners, A., et al. 2020, Astronomy & Astrophysics, 640, L5
  * Yurchenko et al. (2011) Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1828
  * Zhang et al. (2017) Zhang, J., Kempton, E. M.-R., & Rauscher, E. 2017, The Astrophysical Journal, 851, 84

##  Appendix A Observational conditions

Figure 11: Observational conditions for night 1 (blue) and night 2 (orange) as
a function of phase. Shown are the seeing (first panel), relative humidity
(second panel), airmass (third panel), total flux per exposure (fourth panel)
and median FWHM of the PSF in spatial direction (fifth panel). The worse data
quality of night 2 can be contributed to clouds (as seen by the drop in total
flux) and the lack of AO causing a significantly larger PSF.

##  Appendix B Additional $K_{p}$-$v_{\mathrm{sys}}$ maps

Figure 12: $K_{p}$-$v_{\mathrm{sys}}$ map for the CO signal (top panel) and
H2O signal (bottom panel) in night 2 after eight SYSREM iterations. The white
dotted lines indicate the expected position of a signal at the literature
$K_{p}$ value in the planetary rest-frame. The absence of a significant signal
in either of these species can be attributed to the poor data quality of night
2. Figure 13: $K_{p}$-$v_{\mathrm{sys}}$ map for the CH4 signal in night 1
after ten SYSREM iterations (top panel) and night 2 after eight
iterations(bottom panel). The white dotted lines indicate the expected
position of a signal at the literature $K_{p}$ value in the planetary rest-
frame.

##  Appendix C Posterior distribution of C/O retrieval

Figure 14: Posterior distributions of the retrieval using C/O and [M/H] as
free parameters instead of the CO and H2O abundances.


# Long-living excited states of a 2D diamagnetic exciton

R. E. Putnam, Jr. and M. E. Raikh  Department of Physics and Astronomy,
University of Utah, Salt Lake City, UT 84112

###### Abstract

Hydrogenic excited states of a 2D exciton are degenerate. In the presence of a
weak magnetic field, the $S$-states with a zero momentum of the center of mass
get coupled to the $P$-states with finite momentum of the center of mass. This
field-induced coupling leads to a strong modification of the dispersion
branches of the exciton spectrum. Namely, the lower branch acquires a shape of
a “mexican hat” with a minimum at a finite momentum. At certain magnetic
field, exciton branches exhibit a linear crossing, similarly to the spectrum
of a 2D electron in the presence of spin-orbit coupling. While spin is not
involved, degenerate $S$ and $P$ states play the role of the spin projections.
Lifting of degeneracy due to diamagnetic shifts and deviation of electron-hole
attraction from purely Coulomb suppresses the linear crossing.

##  I Introduction

Diamagnetic exciton, a bound state of electron and hole in a magnetic field,
was studied theoretically and experimentallySeisyan  for many years.
Originally, the attention was focused on the bulk semiconductors.Gorkov
Later, the interest has shifted to the two-dimensional systems.Lozovik ;
Kallin ; Macdonald1986 ; Butov ; Zimmermann  It was established that, for
interband excitonsLozovik  as well as for inter-Landau-level excitons,Kallin
that the exciton dispersion law has a local minimum at momenta of the order of
the inverse magnetic length.

Figure 1: (Color online) Three branches of $n=2$ exciton spectrum are plotted
using Eq. (43). Different colors correspond to different solutions of the
cubic equation Eq. (II). Energy is measured in the units of $\varepsilon_{0}$
defined by Eq. (22), while the momentum $Q$ is measured in the units of
$\left(\frac{2M\varepsilon_{0}}{\hbar^{2}}\right)^{1/2}$. Plot a) corresponds
to the resonance $E_{SP}=\hbar\Omega_{-}=\frac{\varepsilon_{0}}{10}$ for which
the dispersion is linear as predicted by Eq. (45). Plot b) corresponds to
$E_{SP}=\frac{\varepsilon_{0}}{20}$ and
$\hbar\Omega_{-}=\frac{\varepsilon_{0}}{5}$. The structure of the spectrum
corresponds to the limiting case Eq. (27) with a well pronounced minimum in
the lower branch.

Recently, Zeeman ; CrookerNano ; Veligzhanin ; Crooker ; Luminescence ;
Luminescence1 ; Potemski1 ; Potemski2 ; report  the exciton spectroscopy in a
perpendicular magnetic field was applied to the novel van der Waals
monolayers. These materials host a series of exciton Rydberg states
corresponding to the principal quantum number $n=1,2,...$. This property is a
consequence of strong Coulomb interaction resulting from the reduced
dielectric screening. Experimental spectra in Refs. Zeeman, ; CrookerNano, ;
Veligzhanin, ; Crooker, ; Luminescence, ; Luminescence1, ; Potemski1, ;
Potemski2, ; report,  were interpreted as diamagnetic shifts of the $S$-states
of the exciton. The growth of $\Big{[}E_{n}(B)-E_{n}(0)\Big{]}$ with magnetic
field, ${\bf B}$, gets faster with increasing $n$. Note, that the $P$-states
of the exciton do not show up in the absorption experiments since the matrix
element between the vacuum and the $P$-states is zero. However, $P$-states can
manifest themselves in luminescence.

In the present paper we consider excited states of a diamagnetic exciton in a
weak magnetic field. Our main finding is that a non-quantizing magnetic field
still affects strongly the exciton dispersion. The underlying reason is that,
due to accidental degeneracy of the excited states, the motion of the center
of mass in a finite field couples different states of the internal motion. It
is this coupling, together with lifting of the accidental degeneracy due to
deviation from Coulomb attraction at short distances,Rytova ; Keldysh  that
modifies the dispersion law even in a weak magnetic field. We demonstrate that
such a modification can give rise to the loop of extrema in the dispersion
law. In turn, this loop of extrema leads to anomalous broadening of the
exciton absorption line. Another consequence is that, with a minimum in the
dispersion law, an exciton is trapped by arbitrarily weak impurity. Since
recombination of a trapped exciton with rapidly precessing center of mass
requires a big momentum transfer, this state is long-lived.

Figure 2: (Color online) Three branches of $n=2$ exciton spectrum are plotted
using Eq. (43). Energy and momentum are measured in the same units as in Fig.
1. Both panels correspond to $\hbar\Omega_{-}=\frac{\varepsilon_{0}}{30}$,
when two of the branches are almost degenerate at $Q=0$. According to Eq.
(23), all three effective masses are positive for $E_{SP}>\varepsilon_{0}$,
while for $E_{SP}<\varepsilon_{0}$ the effective mass of the lower branch is
negative. This is illustrated by plots $E_{SP}=\frac{6\varepsilon_{0}}{5}$,
Panel (a) and $E_{SP}=\frac{\varepsilon_{0}}{10}$, Panel (b).

##  II Dispersion law of a 2D diamagnetic exciton with $n=2$

We start with a standard Hamiltonian of diamagnetic excitonGorkov ; Lozovik ;
Butov ; Zimmermann

$\displaystyle\hat{H}=\frac{1}{2m_{e}}\Big{[}-i\hbar\nabla_{e}+\frac{e}{c}{\bf
A}({\bm{r}}_{e})\Big{]}^{2}$
$\displaystyle+\frac{1}{2m_{h}}\Big{[}-i\hbar\nabla_{h}-\frac{e}{c}{\bf
A}({\bm{r}}_{h})\Big{]}^{2}-\frac{e^{2}}{\kappa\big{|}{\bf r_{e}}-{\bf
r_{h}}\big{|}}.$ (1)

Here $m_{e}$ and $m_{h}$ are the masses of electron and hole, $\kappa$ is the
dielectric constant. Vector potential ${\bm{A}}({\bm{r}})$ is defined as
${\bm{A}}({\bm{r}})=\frac{1}{2}{\bm{B}}\times{\bm{r}}$.

Upon introducing the center of mass and the relative coordinates

${\bf R}=\frac{m_{e}{\bf r_{e}}+m_{h}{\bf r_{h}}}{m_{e}+m_{h}},~{}~{}{\bf
r}={\bf r}_{e}-{\bf r}_{h},$ (2)

the Hamiltonian Eq. (II) acquires the form

$\hat{H}=\hat{H}_{cm}({\bm{R}})+\hat{H}_{rel}({\bm{r}})+\hat{H}_{c}({\bm{r}},{\bm{R}}),$
(3)

where $\hat{H}_{cm}({\bm{R}})$ describes the motion of the center of mass,
$\hat{H}_{rel}({\bm{r}})$ describes the relative motion, and
$\hat{H}_{c}({\bm{r}},{\bm{R}})$ describes their coupling. Analytical
expressions for the three Hamiltonians are the following

$\displaystyle\hat{H}_{cm}=-\frac{\hbar^{2}}{2M}\Delta_{\bf R},$ (4)
$\displaystyle\hat{H}_{rel}=-\frac{\hbar^{2}}{2\mu}\Delta_{\bf
r}-\frac{e^{2}}{\kappa r}+\frac{\left(eBr\right)^{2}}{8\mu c^{2}}$
$\displaystyle+i\hbar\Bigl{(}\frac{eB}{2c}\Bigr{)}\Bigl{(}\frac{1}{m_{e}}-\frac{1}{m_{h}}\Bigr{)}\left(y\partial_{x}-x\partial_{y}\right),$
(5)
$\displaystyle\hat{H}_{c}=i\hbar\frac{eB}{Mc}\left(y\partial_{X}-x\partial_{Y}\right),$
(6)

where $M=m_{e}+m_{h}$ is the net mass and $\mu=m_{e}m_{h}/(m_{e}+m_{h})$ is
the reduced mass of the exciton.

In the absence of magnetic field, the states $n=2$ of the exciton are triply
degenerate. At finite field, the degeneracy is lifted by the third and the
fourth terms in $\hat{H}_{rel}$. Due to the third term, diamagnetic shifts are
different for the $S$ and $P$-states. Fourth term results in the repulsion of
the two $P$-states. In addition, the $S$ and $P$-states are coupled to each
other via the center-of-mass motion. This coupling is promulgated by the
Hamiltonian $\hat{H}_{c}$. Our main point is that the interplay of the three
effects leads to a nontrivial modification of the dispersion law of the
exciton. In contrast to Refs. Lozovik, ; Kallin, , this nontrivial
modification takes place at low fields, so that all three effects can be taken
into account perturbatively.

Explicit form of the $n=2$ wavefunctions is the following

$\displaystyle\psi_{0}({\bm{r}})=C_{0}\Big{(}1-\frac{r}{l}\Big{)}\exp\left[-\frac{r}{2l}\right],$
(7)
$\displaystyle\psi_{x}({\bm{r}})=C_{x}\left(\frac{x}{l}\right)\exp\left[-\frac{r}{2l}\right],~{}\psi_{y}({\bm{r}})=C_{y}\left(\frac{y}{l}\right)\exp\left[-\frac{r}{2l}\right],$
(8)

where $l$ is expressed via the Bohr radius $a_{B}=\frac{\hbar^{2}\kappa}{\mu
e^{2}}$ as $l=\frac{3}{4}a_{B}$. The wavefunctions Eq. (7) correspond to the
binding energy $\frac{4}{9}E_{B}$, where $E_{B}=\frac{\mu
e^{4}}{2\hbar^{2}\kappa^{2}}$ is the Bohr energy. Normalization constants are
the same for all three functions

$C_{0}=C_{x}=C_{y}=\frac{1}{\left(6\pi\right)^{1/2}l}.$ (9)

With the help of the wavefunctions Eq. (7) we calculate the diamagnetic shifts

$\displaystyle E_{S}=\frac{e^{2}B^{2}}{8\mu
c^{2}}\Big{<}r^{2}\Big{>}_{S},~{}~{}E_{P}=\frac{e^{2}B^{2}}{8\mu
c^{2}}\Big{<}r^{2}\Big{>}_{P}.$ (10)

Elementary integration yields

$\Big{<}r^{2}\Big{>}_{S}=26l^{2},~{}~{}~{}\Big{<}r^{2}\Big{>}_{P}=20l^{2}.$
(11)

Thus, the shifts of $2S$ and $2P$ states are related as $E_{S}=1.3E_{P}$, with
$2S$ state being higher.

Next we evaluate the coupling coefficient between the states $\psi_{0}$ and
$\psi_{x},\psi_{y}$. As follows from the form of the Hamiltonian
$\hat{H}_{c}$, this matrix element contains

$x_{SP}=\int d{\bf r}\psi_{0}({\bf r})x\psi_{x}({\bf r}).$ (12)

Using the wavefunctions Eq. (7), we find

$x_{SP}=C_{0}^{2}\int\limits_{0}^{\infty}drr\int_{0}^{2\pi}d\phi~{}\cos^{2}\phi~{}\frac{r^{2}}{l}\left(1-\frac{r}{l}\right)\exp\left[-\frac{r}{l}\right].$
(13)

Angular integration yields $\pi$. Performing the radial integration, we find

$x_{SP}=\pi
C_{0}^{2}\int\limits_{0}^{\infty}dr~{}\frac{r^{3}}{l}\left(1-\frac{r}{l}\right)\exp\left[-\frac{r}{l}\right]=-3l.$
(14)

Since the Hamiltonian $\hat{H}_{c}$ contains $\partial_{X}$ and
$\partial_{Y}$, the coupling between the $S$ and $P$-states depends on the
motion of the center of mass. If the center of mass moves with momentum, $Q$,
the general form of the exciton wavefunction can be presented as a linear
combination

$\displaystyle\Psi({\bf r},{\bf R})=\exp[i{\bf
QR}]\Big{[}A_{0}\psi_{0}({\bm{r}})+A_{x}\psi_{x}({\bm{r}})+A_{y}\psi_{y}({\bm{r}})\Big{]}.$
(15)

Substituting the form Eq. (15) into the Schrödinger equation,
${\hat{H}}\Psi=E\Psi$, yields the system of equations for the coefficients
$A_{0},A_{x}$, and $A_{y}$

$\displaystyle\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{S}-E\Bigg{)}A_{0}-\hbar\Omega_{+}x_{SP}(Q_{y}A_{x}-Q_{x}A_{y})=0,$
(16)
$\displaystyle\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{P}-E\Bigg{)}A_{x}-i\hbar\Omega_{-}A_{y}-\hbar\Omega_{+}Q_{y}x_{SP}A_{0}=0,$
(17)
$\displaystyle\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{P}-E\Bigg{)}A_{y}+i\hbar\Omega_{-}A_{x}-\hbar\Omega_{+}Q_{x}x_{SP}A_{0}=0.$
(18)

Here we have introduced two magnetic-field-induced energy scales

$\hbar\Omega_{-}=\hbar\frac{eB}{2c}\Big{(}\frac{1}{m_{e}}-\frac{1}{m_{h}}\Big{)},~{}\hbar\Omega_{+}=\hbar\frac{eB}{Mc}.$
(19)

Expressing $A_{x}$, $A_{y}$ from Eqs. (17), (18) and substituting them into
Eq. (16), we arrive at the closed equation for the dispersion $E(Q)$

$\displaystyle\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{S}-E\Bigg{)}\Bigg{[}\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{P}-E\Bigg{)}^{2}-(\hbar\Omega_{-})^{2}\Bigg{]}$
$\displaystyle=\Big{(}\hbar\Omega_{+}x_{SP}\Big{)}^{2}Q^{2}\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{P}-E\Bigg{)}.$
(20)

Three solutions of Eq. (II) determine three branches of the exciton
dispersion. We analyze this dispersion in the next Section.

##  III Limiting Cases

###  III.1 $m_{e}=m_{h}$

It is seen from Eq. (19) that the scales $\Omega_{-}$ and $\Omega_{+}$ are
strongly different under the condition $|m_{e}-m_{h}|\ll\mu$. When this
condition is satisfied, we can neglect $\Omega_{-}$ in Eq. (II). After that,
three branches of the exciton spectrum can be easily found. While one branch
is purely parabolic: $E=E_{P}+\frac{\hbar Q^{2}}{2M}$, two other branches
satisfy the quadratic equation

$\displaystyle\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{S}-E\Bigg{)}\Bigg{(}\frac{\hbar^{2}Q^{2}}{2M}+E_{P}-E\Bigg{)}=\varepsilon_{0}\frac{\hbar^{2}Q^{2}}{2M},$
(21)

where we have introduced the energy

$\varepsilon_{0}=2M\Omega_{+}^{2}x_{SP}^{2}.$ (22)

Consider the expression for the lower branch

$\displaystyle
E-\frac{\hbar^{2}Q^{2}}{2M}=\frac{E_{S}+E_{P}}{2}-\Bigg{[}\frac{(E_{S}-E_{P})^{2}}{4}+\frac{\hbar^{2}Q^{2}}{2M}\varepsilon_{0}\Bigg{]}^{1/2}.$
(23)

Our main point is that the spectrum, $E(Q)$, has a minimum for large enough
$\varepsilon_{0}$. From Eq. (23) we find the position of minimum

$Q_{\text{min}}=\Big{(}\frac{M}{2\varepsilon_{0}\hbar^{2}}\Big{)}^{1/2}\Big{[}\varepsilon_{0}^{2}-\left(E_{S}-E_{P}\right)^{2}\Big{]}^{1/2}.$
(24)

Substituting $Q_{\text{min}}$ back into Eq. (23), we find the depth of the
minimum

$E(Q_{\text{min}})-E(0)=-\frac{1}{4\varepsilon_{0}}\Big{[}E_{S}-E_{P}-\varepsilon_{0}\Big{]}^{2}.$
(25)

We see that the minimum emerges when
$\varepsilon_{0}>\left(E_{S}-E_{P}\right)$. Note that both $\varepsilon_{0}$
and $\left(E_{S}-E_{P}\right)$ depend on magnetic field as $B^{2}$. Then the
condition of minimum takes the form

$x_{SP}^{2}>\frac{1}{16}\left(\frac{M}{\mu}\right)\Bigg{[}\Big{<}r^{2}\Big{>}_{S}-\Big{<}r^{2}\Big{>}_{P}\Bigg{]}.$
(26)

With the help of Eqs. (11), (14) the above condition reduces to $M<24\mu$. We
thus conclude that the minimum in the exciton spectrum is quite generic.

###  III.2 $\hbar\Omega_{-}\gg\left(E_{S}-E_{P}\right)$

In this limit we can neglect the difference, $\left(E_{S}-E_{P}\right)$, in
Eq. (II). Then the first branch is still purely parabolic,
$E=E_{P}+\frac{\hbar^{2}Q^{2}}{2M}$, while two other branches are given by

$E-\frac{\hbar^{2}Q^{2}}{2M}=E_{P}\pm\Big{[}\left(\hbar\Omega_{-}\right)^{2}+\frac{\hbar^{2}Q^{2}}{2M}\varepsilon_{0}\Big{]}^{1/2}.$
(27)

Similarly to Eq. (23) the lower branch develops a minimum when the condition
$\varepsilon_{0}>2\hbar\Omega_{-}$ is met. Note, that $\varepsilon_{0}$ grows
with $B$ quadratically, while $\Omega_{-}$ grows with $B$ linearly. Thus, the
minimum can be enforced upon increasing magnetic field.

Figure 3: (Color online) The shifts, $\delta_{S}$ (blue curve) and
$\delta_{P}$ (red curve), of the exciton levels are plotted versus the
dimensionless screening parameter, $r_{0}/l$, from Eqs. (35) and (37),
respectively. At small $r_{0}/l$ the shift $\delta_{S}$ grows linearly with
$r_{0}$, while $\delta_{P}$ grows quadratically.

##  IV Modification of $E_{S}$ and $E_{P}$ by the Keldysh potential

At short distances between electron and hole their attraction deviates from
purely Coulomb as a result of potential created by polarization charges
created in the 2D plane, see e.g. Refs. Rytova, ; Keldysh, ; Cudazzo, ;
Lewenkopf, ; Dery, . The modified interaction can be parametrized by a single
parameter, $r_{0}$, proportional to polarizability, and has a form

$\displaystyle V_{\scriptscriptstyle eff}(\boldsymbol{r})=\frac{\pi
e^{2}}{2r_{0}}\Bigg{[}H_{0}\Big{(}\frac{{r}}{r_{0}}\Big{)}-Y_{0}\Big{(}\frac{r}{r_{0}}\Big{)}\Bigg{]},$
(28)

where $H_{0}$ and $Y_{0}$ are the Struve and second-kind Bessel functions,
respectively. With modified interaction, the accidental degeneracy of the
$n=2$ exciton level is lifted even in the absence of magnetic field. Since the
shifts of the $S$ and $P$-levels are relatively small, the deviation of
$V_{\scriptscriptstyle eff}(\boldsymbol{r})$ from the Coulomb attraction can
be taken into account perturbatively. Corrections to the energy of $S$ and
$P$-states read

$\displaystyle\delta_{S}=\int d{\bf r}\Bigl{(}V_{\scriptscriptstyle
eff}(r)-\frac{e^{2}}{\kappa r}\Bigr{)}\psi_{0}({\bm{r}})^{2},$
$\displaystyle\delta_{P}=\int d{\bf r}\Bigl{(}V_{\scriptscriptstyle
eff}(r)-\frac{e^{2}}{\kappa r}\Bigr{)}\psi_{x}({\bm{r}})^{2}.$ (29)

It is convenient to evaluate $\delta_{S}$ and $\delta_{P}$ in the momentum
space. For this purpose, we present the Fourier transform of
$V_{\scriptscriptstyle eff}(\boldsymbol{r})$ in the form

$\Phi(q)=-\frac{2\pi e^{2}}{\kappa q}\Big{(}1-\frac{qr_{0}}{1+qr_{0}}\Big{)}$
(30)

and treat the second term as a perturbation. Then the expressions for
$\delta_{S}$, $\delta_{P}$ take the form

$\displaystyle\delta_{S}=\frac{2\pi e^{2}r_{0}}{\kappa}\int\frac{d{\bf
q}~{}q}{1+qr_{0}}U_{0}(q),~{}\delta_{P}=\frac{2\pi
e^{2}r_{0}}{\kappa}\int\frac{d{\bf q}~{}q}{1+qr_{0}}U_{x}({\bf q}),$ (31)

where $U_{0}(q)$ and $U_{x}({\bf q})$ stand for the Fourier transforms

$\displaystyle U_{0}(q)=\int d{\bf r}\psi_{0}({\bm{r}})^{2}\exp\left(i{\bf
qr}\right),$ $\displaystyle U_{x}({\bf q})=\int d{\bf
r}\psi_{x}({\bm{r}})^{2}\exp\left(i{\bf qr}\right).$ (32)

In the first transform, angular integration of $\exp\left(i{\bf qr}\right)$
yields the zero-order Bessel function, $J_{0}(qr)$. Then the radial
integration is performed with the help of the identity

$\displaystyle\int\limits_{0}^{\infty}d{t}J_{0}(t)\exp\left({-\gamma
t}\right)=\frac{1}{(\gamma^{2}+1)^{1/2}}.$ (33)

The result reads

$U_{0}(q)=\frac{1-3q^{2}l^{2}+q^{4}l^{4}}{\left(1+q^{2}l^{2}\right)^{7/2}}.$
(34)

Substituting this result into Eq. (IV), we rewrite $\delta_{S}$ in the
dimensionless form

$\displaystyle\delta_{S}=\frac{e^{2}}{\kappa
l}F_{S}\left(\frac{r_{0}}{l}\right),~{}F_{S}=\frac{r_{0}}{l}\int\limits_{0}^{\infty}\frac{dQQ}{1+\frac{r_{0}}{l}Q}\frac{1-3Q^{2}+Q^{4}}{\left(1+Q^{2}\right)^{7/2}}.$
(35)

We see that the integral is a function of the dimensionless ratio $r_{0}/l$.

Angular integration in $U_{x}({\bf q})$ is less trivial. This is because both
the exponent $\exp\left(i{\bf qr}\right)$ and the pre-exponential factor
contain the angle $\phi_{\bf r}$. One has

$\displaystyle U_{x}({\bf
q})=C^{2}_{0}\\!\\!\int\limits_{0}^{2\pi}\int\limits_{0}^{\infty}d{\phi_{\bf
r}}drr\Big{(}\frac{r\cos\phi_{\boldsymbol{r}}}{l}\Big{)}^{2}\\!\exp\Big{(}-\frac{r}{l}\Big{)}$
$\displaystyle\times\exp\Big{(}iqr\cos(\phi_{{\bf
q}}-\phi_{\boldsymbol{r}})\Big{)}.$ (36)

The result of integration contains a constant part and two other parts
proportional to $\cos(2\phi_{\bf q})$ and $\sin(2\phi_{\bf q})$. Only a
constant part contributes to $\delta_{P}$. To calculate this constant part it
is sufficient to replace $\cos\phi_{\boldsymbol{r}}^{2}$ by $1/2$. Subsequent
steps are similar to those in calculation of $\delta_{S}$ and leads to the
result

$\delta_{P}=\frac{e^{2}}{\kappa
l}F_{P}\left(\frac{r_{0}}{l}\right),~{}F_{P}=\frac{r_{0}}{l}\int\limits_{0}^{\infty}\frac{dQQ}{1+\frac{r_{0}}{l}Q}\frac{1-\frac{3}{2}Q^{2}}{\left(1+Q^{2}\right)^{7/2}}.$
(37)

At small $r_{0}\ll l$ the integral is proportional to $r_{0}$. This is a
consequence of the fact that the wavefunction of the $P$-state turns to zero
at the origin. The functions $F_{S}\left(\frac{r_{0}}{l}\right)$ and
$F_{P}\left(\frac{r_{0}}{l}\right)$ are plotted in Fig. 3. Realistic value of
$\frac{r_{0}}{l}$ for transition-metal dichalcogenides can be estimated e.g.
using the data of Ref. Crooker, . From the spectroscopic measurements the
value $r_{0}\approx 3.5$nm was inferred, while the radius of the ground-state
wavefunction $a_{B}\approx 1.5$nm. Thus, for the first excited state, the
ratio $\frac{r_{0}}{l}$ is $\sim 1$.

##  V General case

Corrections $\delta_{S}$, $\delta_{P}$ are independent of magnetic field,
while $E_{S}$ and $E_{P}$ are proportional to $B^{2}$. Thus, the difference

$E_{SP}=\Big{(}E_{S}+\delta_{S}\Big{)}-\Big{(}E_{P}+\delta_{P}\Big{)}$ (38)

is a linear function of magnetic field. Two other parameters,
$\hbar\Omega_{-}$ and $\varepsilon_{0}$, which enter into the equation Eq.
(II), are proportional to $B$ and to $B^{2}$, respectively. Overall, the
evolution of the exciton branches with magnetic field is quite nontrivial. To
examine this evolution, we turn to the analytical solutions of the cubic
equation. The shortest way to arrive to these solutions is performing the
following substitution in Eq. (II)

$E=\frac{\hbar^{2}Q^{2}}{2M}+\frac{E_{SP}}{3}-\Bigg{[}\frac{E_{SP}^{2}}{3}+\left(\hbar\Omega_{-}\right)^{2}+\varepsilon_{0}\frac{\hbar^{2}Q^{2}}{2M}\Bigg{]}^{1/2}\eta,$
(39)

Then the cubic equation for $\eta$ assumes the form

$\eta^{3}-\eta+f=0.$ (40)

Here the dimensionless parameter $f$ is the following combination of $E_{SP}$,
$\hbar\Omega_{-}$, and $\varepsilon_{0}$

$f(Q)=\frac{E_{SP}\Big{[}\frac{2}{9}E_{SP}^{2}-2\left(\hbar\Omega_{-}\right)^{2}+\varepsilon_{0}\frac{\hbar^{2}Q^{2}}{2M}\Big{]}}{3\Big{[}\frac{1}{3}E_{SP}^{2}+\left(\hbar\Omega_{-}\right)^{2}+\varepsilon_{0}\frac{\hbar^{2}Q^{2}}{2M}\Big{]}^{3/2}}.$
(41)

Three solutions of Eq. (40) can be easily expressed via the phase $\varphi$
defined as

$\varphi(Q)=\arctan\left\\{\frac{1}{f(Q)}\left(\frac{4}{27}-f(Q)^{2}\right)^{1/2}\right\\}.$
(42)

Analytical form of these solutions is the following

$\displaystyle\eta_{0}=-\frac{2}{\sqrt{3}}\cos\Big{(}\frac{\varphi}{3}\Big{)},~{}\eta_{\pm}=-\frac{2}{\sqrt{3}}\cos\Big{(}\frac{\varphi}{3}\pm\frac{2\pi}{3}\Big{)}.$
(43)

We see that the character of solutions changes as $f$ passes through the value
$\frac{2}{3^{3/2}}$, when the phase passes through zero. In particular, for
$Q=0$ the value $f=\frac{2}{3^{3/2}}$ is achieved under the condition
$E_{SP}=\hbar\Omega_{-}$. In fact, this condition corresponds to the linear
crossing of the two branches. Introducing the deviation

$\hbar\Omega_{-}-E_{SP}=\Delta$ (44)

and expanding Eq. (II) at small $Q$ and $\Delta$, we find behavior of two
close branches near the condition $E_{SP}=\hbar\Omega_{-}$

$E=\frac{\hbar^{2}Q^{2}}{2M}+\frac{\Delta}{2}\pm\Bigg{[}\left(\frac{\Delta}{2}\right)^{2}+\frac{\varepsilon_{0}}{2}\left(\frac{\hbar^{2}Q^{2}}{2M}\right)\Bigg{]}^{1/2}.$
(45)

We see that a gap of a width, $\Delta$, opens in the spectrum at finite
$\Delta$.

Finally, by expanding Eq. (II), we find the behavior of all three branches
near $Q=0$

$\displaystyle
E=E_{SP}+\frac{\hbar^{2}Q^{2}}{2M}\Bigg{[}1+\frac{\varepsilon_{0}E_{SP}}{E_{SP}^{2}-\left(\hbar\Omega_{-}\right)^{2}}\Bigg{]},$
(46) $\displaystyle
E=-\hbar\Omega_{-}+\frac{\hbar^{2}Q^{2}}{2M}\Bigg{[}1-\frac{\varepsilon_{0}}{2\left(E_{SP}+\hbar\Omega_{-}\right)}\Bigg{]},$
(47) $\displaystyle
E=\hbar\Omega_{-}+\frac{\hbar^{2}Q^{2}}{2M}\Bigg{[}1-\frac{\varepsilon_{0}}{2\left(E_{SP}-\hbar\Omega_{-}\right)}\Bigg{]}.$
(48)

We conclude that negative effective mass at $Q=0$ emerges at
$\left(E_{SP}\pm\hbar\Omega_{-}\right)<\frac{1}{2}\varepsilon_{0}$.
Approaching of denominators in Eqs. (47), (48) to zero signals the proximity
to the exciton resonance.

##  VI Discussion

Our most nontrivial finding is the exciton resonance originating from the
accidental degeneracy of the hydrogen-like excited levels. The resonant
condition reads $E_{SP}=\pm\hbar\Omega_{-}$. It corresponds to the magnetic
field at which the mismatch of $S$ and $P$ exciton levels due to diamagnetic
shift as well as due to field-independent Keldysh potential, is equal to the
field-induced splitting of the degenerate $P$-states. All three branches of
the bare exciton spectrum ($S$-branch and two $P$-branches) are involved into
the formation of the resonance. Under the resonant condition, two branches of
the modified spectrum cross linearly, as illustrated in Fig. 1a. Away from the
resonance, see Figs. 1b and 2b, the lowest branch of the spectrum evolves into
a “mexican hat” as predicted by Eqs. (27), (23) for the limiting cases
$E_{SP}\ll\hbar\Omega_{-}$ and $E_{SP}\gg\hbar\Omega_{-}$, respectively. The
main condition for the strong mixing of the $S$ and $P$-branches is
$\varepsilon_{0}>\left(E_{SP}+\hbar\Omega_{-}\right)$. We have demonstrated
that, neglecting the Keldysh shift, and for $m_{e}=m_{h}$ this condition is
satisfied. Physically, it requires that the matrix element, $x_{SP}$, defined
by Eq. (12) is big enough. When this condition is violated, the effective
masses for all three branches of the modified spectrum are positive, as
illustrated in Fig. 1a.

With regard to physical consequences of the spectrum modification into a
mexican hat, it is knownEfros ; Galstyan ; Mkhitaryan ; Voloshin  that the
density of states near the minimum behaves as
$\big{(}E-E_{min}\big{)}^{-1/2}$, i.e. it has a one-dimensional character. As
a result, even a weak attractive impurity can trap an exciton. The wave
function of the trapped exciton is centered around $Q_{min}$ in the momentum
space. With $Q_{min}$ exceeding the momentum of a photon required for
radiative recombination, these trapped states are long-lived.

Throughout the paper we considered the simplest model of diamagnetic exciton
without account for the valley effects which lead to the polarization
dependence of the optical absorption as well as the spin-orbit effects leading
to the dark-bright exciton mixing.Durnev  We have also assumed that the
coupling of the $S$ and $P$-exciton states is exclusively due to magnetic
field and not by the peculiarities of the bandstructure of dichalcogenide
monolayers.Glazov

While the optical absorption probes only the $S$-states of the exciton,
finite-momentum states can be probed in the microcavity setting.microcavity
Another consequence of the field-induced modification of the exciton spectrum
is anomalous broadening of $n=2$ absorption line. The underlying
mechanismsemimagnetic ; WithDisorder  is the elastic scattering of $Q=0$ state
into the states with finite center-of-mass momentum, $Q$. This scattering is
enabled by the mexican-hat shape of the spectrum.

##  VII Acknowledgements

The work was supported by the Department of Energy, Office of Basic Energy
Sciences, Grant No. DE- FG02- 06ER46313.

## References

  * (1) see e.g. the review R. P. Seisyan and B. P. Zakharchenya, “Interband magneto-optics of semiconductors as diamagnetic exciton spectroscopy,” in Modern Problems in Condensed Matter Sciences, Elsevier 27, 345 (1991).
  * (2) L. P. Gorkov and I. E. Dzyaloshinskii, “Contribution to the theory of the Mott exciton in a strong magnetic field,” Zh. Eksp. Teor. Fiz. 53, 717 (1967) [Sov. Phys. JETP 26, 449 (1968)].
  * (3) I. V. Lerner and Yu. E. Lozovik, “Mott exciton in a quasi-two-dimensional semiconductor in a strong magnetic field,” Zh. Eksp. Teor. Fiz. 78, 1167 (1980) [Sov. Phys. JETP 51, 588 1980].
  * (4) C. Kallin and B. I. Halperin, “Excitations from a filled Landau level in the two-dimensional electron gas,” Phys. Rev. B 30 5655, (1984).
  * (5) A. H. MacDonald and D. S. Ritchie, “Hydrogenic energy levels in two dimensions at arbitrary magnetic fields,” Phys. Rev. B 33 8336 (1986).
  * (6) Yu. E. Lozovik, I. V. Ovchinnikov, S. Yu. Volkov, L. V. Butov, and D. S. Chemla, “Quasi-two-dimensional excitons in finite magnetic fields,” Phys. Rev. B, 65, 235304 (2002).
  * (7) M. Grochol, F. Grosse, and R. Zimmermann, “Exciton wave function properties probed by diamagnetic shift in disordered quantum wells,” Phys. Rev. B 71, 125339 (2005).
  * (8) Y. Li, J. Ludwig, T. Low, A. Chernikov, X. Cui, G. Arefe, Y. D. Kim, A. M. van der Zande, A. Rigosi, H. M. Hill, S. H. Kim, J. Hone, Z. Li, D. Smirnov, and T. F. Heinz, “Valley Splitting and Polarization by the Zeeman Effect in Monolayer MoSe2,” Phys. Rev. Lett. 113, 266804 (2014).
  * (9) A. V. Stier, N. P. Wilson, G. Clark, X. Xu, and S. A. Crooker, “Probing the Influence of Dielectric Environment on Excitons in Monolayer WSe2: Insight from High Magnetic Fields,” Nano Lett. 16, 7054 (2016).
  * (10) A. V. Stier, N. P. Wilson, K. A. Velizhanin, J. Kono, X. Xu, and S. A. Crooker, “Magnetooptics of Exciton Rydberg States in a Monolayer Semiconductor,” Phys. Rev. Lett. 120, 057405 (2018).
  * (11) M. Goryca, J. Li, A. V. Stier, T. Taniguchi, K. Watanabe, E. Courtade, S. Shree, C. Robert, B. Urbaszek, X. Marie, and S. A. Crooker, “Revealing exciton masses and dielectric properties of monolayer semiconductors with high magnetic fields,” Nat. Commun. 10, 4172 (2019).
  * (12) E. Liu, J. van Baren, T. Taniguchi, K. Watanabe, Y.-C. Chang, and C. H. Lui, “Magnetophotoluminescence of exciton Rydberg states in monolayer WSe2,” Phys. Rev. B 99, 205420 (2019).
  * (13) S.-Y. Chen, Z. Lu, T. Goldstein, J. Tung, A. Chaves, J. Kunstmann, L. S. R. Cavalcante, T. Woznyak, G. Seifert, D. R. Reichman, T. Taniguchi, K. Watanabe, D. Smirnov, and J. Yan, “Luminescent emission of excited Rydberg excitons from monolayer WSe2,” Nano Lett. 19, 2464 (2019).
  * (14) M. R. Molas, A. O. Slobodeniuk, K. Nogajewski, M. Bartos, L. Bala, A. Babinski, K. Watanabe, T. Taniguchi, C. Faugeras, and M. Potemski, “Energy spectrum of two-dimensional excitons in a non-uniform dielectric medium,” Phys. Rev. Lett. 123, 136801 (2019).
  * (15) P. Kapuscinski, A. Delhomme, D. Vaclavkova, A. Slobodeniuk, M. Grzeszczyk, M. Bartos, K. Watanabe, T. Taniguchi, C. Faugeras, and M. Potemski, “Rydberg series of dark excitons and the conduction band spin-orbit splitting in monolayer WSe2” arXiv:2102.12286.
  * (16) J. Jadczak, J. Kutrowska-Girzycka, T. Smolenski, P. Kossacki, Y. S. Huang, and L. Bryja, “Exciton binding energy and hydrogenic Rydberg series in layered ReS2,” Sci. Rep. 9, 1578 (2019).
  * (17) N. S. Rytova, “Screened potential of a point charge in a thin film,” Moscow University Physics Bulletin 3, 30 (1967)
  * (18) L. V. Keldysh, “Coulomb interaction in thin semiconductor and semimetal films,” JETP Lett. 29, 658 (1979).
  * (19) P. Cudazzo, I. V. Tokatly, and A. Rubio, “Dielectric screening in two-dimensional insulators: Implications for excitonic and impurity states in graphane,” Phys. Rev. B 84, 085406 (2011).
  * (20) E. Ridolfi, C. H. Lewenkopf, and V. M. Pereira, “Excitonic structure of the optical conductivity in MoS2 monolayers, Phys. Rev. B 97, 205409 (2018).
  * (21) D. V. Tuan, M. Yang, and H. Dery, “The Coulomb interaction in monolayer transition-metal dichalcogenides,” Phys. Rev. B 98, 125308 (2018).
  * (22) M. E. Raikh and Al. L. Efros, “Broadening of an exciton line in solid solutions with a degenerate valence band,” JETP Lett. 45, 380 (1987).
  * (23) A. G. Galstyan and M. E. Raikh, “Disorder-Induced Broadening of the Density of States for 2D Electrons with Strong Spin-Orbit Coupling,” Phys. Rev. B 58, 6736 (1998).
  * (24) V. V. Mkhitaryan and M. E. Raikh, “Disorder-induced tail states in a gapped bilayer graphene,” Rev. B 78, 195409 (2008).
  * (25) B. Skinner, B. I. Shklovskii, and M. B. Voloshin, “Bound state energy of a Coulomb impurity in gapped bilayer graphene,” Phys. Rev. B 89, 0414059(R) (2014).
  * (26) M. V. Durnev and M. M. Glazov, “Excitons and trions in two-dimensional semiconductors based on transition metal dichalcogenides,” Phys. Usp. 61, 825 (2018).
  * (27) M. M. Glazov, L. E. Golub, G. Wang, X. Marie, T. Amand, and B. Urbaszek, “ Intrinsic exciton-state mixing and nonlinear optical properties in transition metal dichalcogenide monolayers,” Phys. Rev. B 95, 035311 (2017).
  * (28) V. Walther, R. Johne, and T. Pohl, “Giant optical nonlinearities from Rydberg excitons in semiconductor microcavities,” Nature communications, Nat. Commun. 9, 1309 (2018).
  * (29) M. E. Raikh and Al. L. Efros, “Exciton line broadening in a semimagnetic semiconductor,” Sov. Phys.-Solid State 30, 1708 (1988).
  * (30) M. Dwedari, S. Brem, M. Feierabend, and E. Malic, “Disorder-induced broadening of excitonic resonances in transition metal dichalcogenide,” Phys. Rev. Mater. 3, 074004 (2019).