# Bridging the gap between luminous red novae and common envelope evolution:
the role of recombination energy and radiation force

Zhuo Chen (陈卓) Institute of Advanced Study, Tsinghua University, Beijing
100084, China  Department of Astronomy, Tsinghua University, Beijing 100084,
China  Natalia Ivanova Department of Physics, University of Alberta, Edmonton,
AB T6G 2E1, Canada

###### Abstract

Luminous red novae (LRNe) and their connection to common envelope evolution
(CEE) remain elusive in astrophysics. Here, we present a radiation
hydrodynamic model capable of simulating the light curves of material ejected
during a CEE. For the first time, the radiation hydrodynamic model
incorporates complete recombination physics for hydrogen and helium. The
radiation hydrodynamic equations are solved with Guangqi. With time-
independent ejecta simulations, we show that the peaks in the light curves are
attributed to radiation-dominated ejecta, while the extended plateaus are
produced by matter-dominated ejecta. To showcase our model’s capability, we
fit the light curve of AT2019zhd. The central mass object of $6M_{\odot}$ is
assumed based on observations and scaling relations. Our model demonstrates
that the ejecta mass of AT2019zhd falls within the range of $0.04M_{\odot}$ to
$0.1M_{\odot}$. Additionally, we demonstrate that recombination energy and
radiation force acceleration significantly impact the light curves, whereas
dust formation has a limited effect during the peak and plateau phases.

Transient sources(1851) — Common envelope evolution(2154)

††software: Guangqi 光启 (Chen & Bai, in prep), Matplotlib (Hunter, 2007), Petsc
(Balay et al., 1997, 2019)

##  1 Introduction

Since the discovery of V1309 Sco (Tylenda et al., 2011), a luminous red nova
(LRN), also a confirmed binary merger, many more new LRNe have been discovered
(Kurtenkov et al., 2015; Blagorodnova et al., 2017; Cai et al., 2019;
Pastorello et al., 2019; Blagorodnova et al., 2021; Pastorello et al., 2021a,
b; Cai et al., 2022a, b), to name a few. It is argued that, at least for some
of the LRNe, the progenitors are binary stars that have undergone common
envelope evolution (CEE) (Ivanova et al., 2013a) \- one of the most mysterious
events in binary evolution. Bridging the gap between the observables of LRNe
and theoretical models of CEE could be crucial to getting a better
understanding of binary evolution (Chen et al., 2024), including cataclysmic
variables (Warner, 2003), X-ray binaries (Reig, 2011), gravitational wave
sources (Renzo et al., 2021), and type Ia supernovae (Wang & Han, 2012; Liu et
al., 2023).

Observationally, LRNe exhibit three prominent characteristic features.
Firstly, they all have extended plateaus with decreasing effective
temperatures in the light curves. This phenomenon suggests a substantial mass
ejection from the central object, which cools as it expands outward. During
this phase, hydrogen and helium may recombine and release their latent heat,
contributing to the thermal energy of the plateau (Ivanova et al., 2013b,
2015; Matsumoto & Metzger, 2022). Secondly, a notable number of LRNe present a
peak preceding the plateau in their light curves. Modeling this peak is
challenging due to its exponential increasing and decreasing phases. Pejcha et
al. (2017) employed a hydrodynamic model with post-processing to fit the
slowly increasing luminosity before the peak, attributing the peak to an
eruption-like event. Although Matsumoto & Metzger (2022) successfully fitted
the rapidly decreasing luminosity after the peak with their hot ejecta model
in 1D, the rapidly increasing phase remains a puzzle. Thirdly, LRNe frequently
exhibit strong H$\alpha$ emission (Munari et al., 2002; Stritzinger et al.,
2020; Blagorodnova et al., 2021; Pastorello et al., 2021a), indicating a
possible collision between the ejecta and circumstellar matter (Blagorodnova
et al., 2020).

A prevailing theory on the origin of LRNe implicates the progenitor binary
undergoes CEE and experiences a plunge-in phase accompanied with ejection of
some envelope material. Numerous 3D simulations have been conducted to study
the CEE dynamics (Nandez et al., 2014, 2015; Ivanova & Nandez, 2016; Nandez &
Ivanova, 2016; Ohlmann et al., 2016; Chamandy et al., 2018; Iaconi et al.,
2019; Prust & Chang, 2019; Reichardt et al., 2019; Sand et al., 2020;
González-Bolívar et al., 2023; Röpke & De Marco, 2023; Chamandy et al., 2024).
However, due to the complexity of physical processes and the multi-scale
nature of the CEE problem, there is still no consensus on how much mass is
ejected during a CEE, even if 3D simulations start with the same initial
conditions. Furthermore, predicting the observational appearance of these
simulations is inherently difficult (Hatfull et al., 2021). Notably, radiation
hydrodynamics, a crucial factor in CEE, is often absent from these 3D models
due to its difficulty and is easier to consider in 1D models (Soker et al.,
2018; Bronner et al., 2023; O’Connor et al., 2023). Recent work by Matsumoto &
Metzger (2022) demonstrates the effectiveness of a cooling shell model, devoid
of radiation transport and radiation hydrodynamics, in producing LRNe-like
light curves.

However, there is still a lack of first principle models that can relate CEE
to LRNe. In CEE, the plunge-in phase marks the rapid conversion of the
gravitational potential energy into the kinetic energy. Some kinetic energy
would be converted to thermal energy via a shock between the plunge-in star
and the envelope (MacLeod et al., 2017). We would expect the shocked gas to be
close to adiabatic because of the high optical depth inside a common envelope
(CE). The shocked gas may be accelerated in the radial direction due to high
radiation pressure and radiation force acceleration – this has been overlooked
previously – and become an eruption-like event, eventually appearing as an
LRN. CEE is a 3D problem by nature. However, solving 3D radiation hydrodynamic
equations with complex physics is computationally demanding. Therefore, as the
first step, and hoping to resolve the microphysics better, we start by
approximating the radiation hydrodynamic problem with 1D spherical symmetry.
In a follow-up work, we will revisit this problem with 2D axisymmetric models.

This Letter introduces the first radiation hydrodynamic model that
comprehensively incorporates recombination energies of H and He, radiation
transport, and radiation force acceleration. The model can successfully
produce both the peak and plateau phases of the light curve, allowing for the
estimation of ejecta mass through curve fitting. The organization of the
Letter is as follows: Section 2 presents the physical model, governing
equations, initial and boundary conditions, and numerical setups, including
Adaptive Mesh Refinement (AMR) criteria. Simulation results are presented in
Section 3. The Letter concludes in Section 4.

##  2 Methods

###  2.1 Physical model

We adopt the following radiation hydrodynamic equations to model the evolution
of the LRNe in a 1D spherical coordinate,

$\displaystyle\partialderivative{\rho}{t}+\frac{1}{r^{2}}\partialderivative{r}(r^{2}\rho
v)$ $\displaystyle=$ $\displaystyle 0,$ (1)
$\displaystyle\partialderivative{\rho
v}{t}+\frac{1}{r^{2}}\partialderivative{r}(r^{2}\rho v^{2})$ $\displaystyle=$
$\displaystyle-\partialderivative{p}{r}+\rho(a_{\rm{rad}}-g),$ (2)
$\displaystyle\partialderivative{E}{t}+\frac{1}{r^{2}}\partialderivative{r}[r^{2}(E+p)v]$
$\displaystyle=$ $\displaystyle\rho v(a_{\rm{rad}}-g)+\mathbb{G},$ (3)
$\displaystyle\partialderivative{E_{r}}{t}+\frac{1}{r^{2}}\partialderivative{r}[r^{2}({F}_{r}+E_{r}v)]$
$\displaystyle=$ $\displaystyle-\mathbb{G},$ (4)

where $\rho,v,p$, and $E$ are the density, radial velocity, pressure, and
total energy of the gas. In addition, $a_{\rm{rad}}$ (to be explained later)
and $g=GM_{\star}/r^{2}$ are the radiation and gravitational force
acceleration, respectively, where $G$ is the gravitational constant and
$M_{\star}$ is the mass of the central object. The time and coordinate are
denoted by $t$ and $r$. The radiation-related variables $E_{r}$, ${F}_{r}$,
and $\mathbb{G}$ are explained later in this section.

We adopt a simple hydrogen and helium mixture equation of state (EoS, see
Appendix A),

$\displaystyle E$ $\displaystyle=$
$\displaystyle\rho\frac{v^{2}}{2}+e_{\text{g}}(\rho,T_{\text{g}})$ (5)
$\displaystyle p$ $\displaystyle=$
$\displaystyle\sum_{i}n_{i}k_{b}T_{\text{g}}$ (6)

where $i$ is the species index, $n_{i}$ is the number density of species $i$,
$e_{\text{g}}$ and $T_{\text{g}}$ are the internal energy and temperature of
the gas, $k_{b}$ is the Boltzmann constant. Throughout this Letter, we assume
the hydrogen mass ratio $X=0.74$ and the helium mass ratio $Y=0.26$ for
simplicity, because metal does not contribute much to gas thermodynamics. The
hydrodynamics, together with the complex EoS, is solved by the approximate
HLLC Riemann solver in Chen et al. (2019).

To solve the radiation transport problem, we use the flux-limited diffusion
(FLD) approximation, which relates the radiation flux ${F}_{r}$ to the
gradient of radiation energy $\partial E_{r}/\partial r$ with a flux limiter
$\lambda$ (Levermore & Pomraning, 1981), i.e.,

$\displaystyle{F}_{r}$ $\displaystyle=$
$\displaystyle-\frac{c\lambda(R)}{\kappa_{\rm{R}}\rho}\partialderivative{E_{r}}{r},$
(7) $\displaystyle\lambda(R)$ $\displaystyle=$
$\displaystyle\frac{2+R}{6+3R+R^{2}},$ (8) $\displaystyle R$ $\displaystyle=$
$\displaystyle\frac{|\partial E_{r}/\partial r|}{\kappa_{\rm{R}}\rho E_{r}},$
(9)

where $\kappa_{\rm{R}}$ is the Rosseland mean opacity (see Appendix B). The
flux limiter has the property that

$\frac{c\lambda}{\kappa_{\rm{R}}\rho}\partialderivative{E_{r}}{r}\rightarrow\begin{cases}\frac{c}{3\kappa_{\rm{R}}\rho}\partialderivative{E_{r}}{r}&\quad\quad\rm{optically\
thick,}\\\ cE_{r}&\quad\quad\rm{optically\ thin.}\\\ \end{cases}$ (10)

The radiation and gas energy coupling is modeled by solving the following
equations implicitly with sub-timesteps,

$\displaystyle\partialderivative{e_{\text{g}}}{t}$ $\displaystyle=$
$\displaystyle\mathbb{G},$ (11)
$\displaystyle\partialderivative{E_{r}}{t}+\frac{1}{r^{2}}\partialderivative{r}[r^{2}({F}_{r}+E_{r}v)]$
$\displaystyle=$ $\displaystyle-\mathbb{G},$ (12)

where $\mathbb{G}=\kappa_{\rm{P}}\rho c(E_{r}-a_{\rm{r}}T_{\text{g}}^{4})$ is
the energy coupling strength, $\kappa_{\rm{P}}$ is the Planck mean opacity,
and $a_{\rm{r}}$ is the radiation constant. Meanwhile, the radiation and
gravitational acceleration are integrated explicitly through,

$\partialderivative{\rho v}{t}=\rho(a_{\rm{rad}}-g),$ (13)

where $a_{\rm{rad}}=\kappa_{\rm{R}}{F}_{r}/c$ is the radiation force
acceleration and $c$ is the speed of the light.

We use Guangqi(Chen & Bai, in prep) to solve Equation 1-4. Guangqi is a
second-order in time and space accurate and finite volume radiation
hydrodynamic code. It has HLLC Riemann solvers, realistic EoS (Chen et al.,
2019), and adaptive mesh refinement (AMR). Guangqi solves the radiation
transport problem with FLD approximation implicitly (similar to Kolb et al.
(2013)), using iterative solvers from Petsc (Balay et al., 1997, 2019).
Currently, Guangqi has spherical and Cartesian geometry in 1D and 2D.

###  2.2 Initial and boundary conditions

We set our initial condition to be an outflow at escape velocity with a
constant mass loss rate,

$\displaystyle\rho_{\rm{init}}(r)$ $\displaystyle=$
$\displaystyle\rho_{0}(r_{\rm{in}}/r)^{3/2},$ (14) $\displaystyle
v_{\rm{init}}(r)$ $\displaystyle=$ $\displaystyle\sqrt{2GM_{\star}/r}$ (15)
$\displaystyle T_{\rm{init}}(r)$ $\displaystyle=$ $\displaystyle
T_{0}(r_{\rm{in}}/r)$ (16)

where $\rho_{0}=10^{-13}$g$\cdot$cm-3, $T_{0}=1000$K, and
$r_{\rm{in}}=10\text{R}_{\odot}$ is the inner radius of the computational
domain. The constant mass loss rate is $\dot{M}=2.30\times
10^{-7}$M${}_{\odot}\cdot$yr-1, which is low compared to the mass loss rate of
the ejecta of the CEE. Radiation’s initial condition is assumed to be in local
thermal equilibrium (LTE) with gas. The initial mass and energy in the
computational domain are $8.52\times 10^{-7}\text{M}_{\odot}$ and $3.68\times
10^{35}$erg, respectively. They are significantly smaller than the ejecta’s
mass and energy, and our light curve results are insensitive to the initial
condition.

The outer boundary is free, i.e., the gas and radiation can leave the
computational domain freely. In our simulations, we set the out boundary at
$r_{\rm{out}}=$4000R⊙. The free boundary for the radiation is,

$\displaystyle\partialderivative{(r^{2}E_{r})}{r}=0.$ (17)

This outer boundary condition means that the radiation flux is optically thin
and outward. We can confirm that the outer region of our computational domain
is indeed optically thin (see Section 3.2) and always outward. We calculated
the luminosity at the outer boundary by,

$L=4\pi r_{\rm{out}}^{2}{F}_{r}(r_{\rm{out}})$ (18)

The inner boundary is time-dependent. We specify the time-dependent density,
velocity, and temperature of the ejecta at the inner boundary, i.e.,
$[\rho(t),v_{\rm{ej}}(t),T_{\text{g}}(t)]$.

$\displaystyle v_{\rm{ej}}$ $\displaystyle=$ $\displaystyle f_{\rm
ej}(t)\sqrt{2GM_{\star}/r_{\rm{in}}},$ (19) $\displaystyle\rho$
$\displaystyle=$ $\displaystyle\dot{M}/(4\pi r_{\rm{in}}^{2}v_{\rm{ej}}),$
(20) $\displaystyle T_{\text{g}}$ $\displaystyle=$ $\displaystyle\alpha
m_{\ce{H}}v_{\rm{ej}}^{2}/2k_{b},$ (21)

where $f_{\rm ej}(t)$ can be a time-dependent factor, $\dot{M}$ is the mass
loss rate, and we introduce $\alpha$ to characterize the temperature of the
ejecta. When the ejection of the CEE stops, the inner boundary of the gas is
changed to free to let the fallback gas pass through the inner boundary,
therefore, we do not model any fallback shocks and fallback accretion disks.
Meanwhile, radiation transport is turned off at the inner boundary.
Mathematically, we use the zero-gradient ($\partial E_{r}/\partial r=0$)
radiation inner boundary condition.

We assume the gas and radiation are in LTE at the inner boundary because the
density and opacity are very high inside the CE, and the thermal timescale is
short.

###  2.3 Numerical setups

We adopt a uniform base grid with spherical geometry. The computational domain
is $r\in[10,4000]$R⊙, and the base resolution is $N=1536$. We add 5 levels of
static mesh refinement where $r\in[10,15]$R⊙to resolve the strong gradient of
the gravitational potential. Each level of mesh refinement doubles the
resolution. We also adaptively refine zones with temperature gradients up to 5
levels. The mesh refine and derefine criterion are as follows,

$\displaystyle\text{refine}\quad\frac{|T_{g,i+1}-T_{g,i-1}|}{2T_{g,i}}$
$\displaystyle>$ $\displaystyle 0.03,$ (22)
$\displaystyle\text{derefine}\quad\frac{|T_{g,i+1}-T_{g,i-1}|}{2T_{g,i}}$
$\displaystyle<$ $\displaystyle 0.001,$ (23)

where $i\pm 1$ in the subscript represents the cell index (not to be confused
with species index). We can capture the shock, radiative layer, and radiation-
dominated zones with photon trapping with AMR. The finest cell has a length of
$8.12\times 10^{-2}$R⊙. The Courant–Friedrichs–Lewy number is 0.95 in our
simulations. The simulation time is $4.5\times 10^{6}$s=52.08 days.

##  3 Results

We first show some simple simulations and get a sense of the correspondence
between the properties of the ejecta and the light curves. After that, we fit
the light curve of AT2019zhd with a more complex ejecta.

###  3.1 Simple ejecta

For simplicity, in this subsection, we set $M_{\star}=6$M⊙,
$r_{\rm{in}}=10$R⊙. We consider seven models and list the physical properties
of the seven models in Table 1. In what follows, we refer to a model as
radiation-dominated if $E_{r}/e_{\text{g}}\gg 1$, and as matter-dominated if
$E_{r}/e_{\text{g}}<1$. In particular, we calculate the total energy
$E_{\rm{total}}$ of the ejecta by,

$E_{\rm{total}}=4\pi
r_{\rm{in}}^{2}(E_{r}+e_{\text{g}})f_{\rm{ej}}\sqrt{2GM_{\star}/r_{\rm{in}}}t_{\rm{ej}},$
(24)

where $t_{\rm{ej}}$ is the duration of the ejection.

model | $f_{\rm{ej}}$ | $t_{\rm{ej}}$ | $\dot{M}$ | $\Delta M$ | $\rho$ | $\alpha$ | $T_{\text{g}}$ | $E_{r}/e_{\text{g}}$ | $E_{\rm{total}}$
---|---|---|---|---|---|---|---|---|---
|  | [days] | [$\text{M}_{\odot}\cdot\text{yr}^{-1}$] | [M⊙] | [g$\cdot$cm-3] |  | [K] |  | [erg]
m01a18 | 1 | 1 | 0.1 | $2.74\times 10^{-4}$ | $2.17\times 10^{-8}$ | 0.018 | $2.50\times 10^{5}$ | 19.9 | $7.76\times 10^{44}$
m02a045v2 | 2 | 0.5 | 0.2 | $2.74\times 10^{-4}$ | $2.17\times 10^{-8}$ | 0.0045 | $2.50\times 10^{5}$ | 19.9 | $7.76\times 10^{44}$
m01a09 | 1 | 1 | 0.1 | $2.74\times 10^{-4}$ | $2.17\times 10^{-8}$ | 0.009 | $1.25\times 10^{5}$ | 2.01 | $6.91\times 10^{43}$
m05a18 | 1 | 1 | 0.5 | $1.37\times 10^{-3}$ | $1.08\times 10^{-7}$ | 0.018 | $2.50\times 10^{5}$ | 3.98 | $9.24\times 10^{44}$
m05a09 | 1 | 1 | 0.5 | $1.37\times 10^{-3}$ | $1.08\times 10^{-7}$ | 0.009 | $1.25\times 10^{5}$ | 0.40 | $1.61\times 10^{43}$
m25a18 | 1 | 1 | 2.5 | $6.85\times 10^{-3}$ | $5.42\times 10^{-7}$ | 0.018 | $2.50\times 10^{5}$ | 0.80 | $1.66\times 10^{45}$
m25a09 | 1 | 1 | 2.5 | $6.85\times 10^{-3}$ | $5.42\times 10^{-7}$ | 0.009 | $1.25\times 10^{5}$ | 0.08 | $6.20\times 10^{44}$
Table 1: From left to right, model name, ejecta’s velocity factor, ejecta
duration $t_{\rm{ej}}$, mass loss rate $\dot{M}$, cumulative mass loss $\Delta
M$, ejecta’s density, $\alpha$, the temperature of the ejecta, the ratio of
$E_{r}$ and $e_{\text{g}}$, and the total energy of the ejecta of the seven
simple models.

Figure 1 shows the light curves of the simple models. We can see that the
light curves of m02a045v2 and m01a18 have pronounced peaks, i.e., rapid,
exponential increases and decreases, while m25a18 and m25a09 have long
plateaus. The light curves of other models are in between. Because the
physical system is highly nonlinear, we provide a qualitative discussion here.
The optical depth of an ejecta decreases as it expands, and the ejecta cools
faster as the optical depth decreases. The more mass in an ejecta, the higher
the optical depth; the faster the ejecta, the faster the optical depth
decreases. Therefore, m25a18 and m25a09 have long plateaus because their
ejecta are the most massive ones among the seven simple models; m25a18 has a
longer plateau than m25a09 because m25a18 has a higher energy budget and thus
a longer cooling timescale. On the other hand, m02a045v2 has a sharper peak
than m01a18 because its expansion speed is faster, resulting in a more rapid
decrease in the optical depth and a shorter cooling timescale. Because
m02a045v2 and m01a18 have the same amount of energy budget, m02a045v2 has a
higher luminosity peak.

Overall, we can also relate that we observe the luminosity peak formation
(exponential increase and decrease) in the case of radiation-dominating
ejecta. In contrast, the plateau is the feature observed in the case of
matter-dominated ejecta.

Figure 1: The light curves of simple simulations with different mass and
temperature of the ejecta. The physical parameters of the ejecta are listed in
Table1.

###  3.2 Fitting AT2019zhd

We choose AT2019zhd as a fitting example because it is a recently well-
observed LRN and its light curve resembles V1309 Sco, whose progenitor is a
confirmed binary (Tylenda et al., 2011; Pastorello et al., 2021a). The mass of
the progenitor of AT2019zhd is unclear (Pastorello et al., 2021a). We take the
following steps to estimate the total mass of the merger product.

  1. 1.

The total mass of V1309 Sco is estimated to be $1-2M_{\odot}$ (Tylenda et al.,
2011).

  2. 2.

The peak luminosity of AT2019zhd is roughly $L_{\rm{AT2019zhd}}=2.08\times
10^{39}$erg$\cdot$s-1, and the peak captured luminosity of V1309 Sco is
$L_{\rm{V1309Sco}}=1.26\times 10^{38}$erg$\cdot$s-1. Therefore,
$L_{\rm{AT2019zhd}}/L_{\rm{V1309Sco}}=16.5$.

  3. 3.

The total radiation energy released during the peak and the plateau of
AT2019zhd is roughly 10 times more than V1309 Sco.

If we assume that the luminosity energy mainly comes from the release of the
gravitational potential energy, it may scale as,

$\displaystyle L\sim M^{2}.$ (25)

This scaling relation motivates us to set the mass of the central object to be
the same as was adopted for the simple models, $6M_{\odot}$, and we set
$\alpha=0.018$ in this subsection.

Unlike the case of simple models, the inner boundary conditions here are time-
dependent, see Figure 2. As a result, the time-dependent ejecta is initially
radiation-dominated (denoted by the orange color in Figure 2) and then
transits to matter-dominated. We anticipate that the initial high-temperature
ejecta is produced by a shock between the rapidly plunge-in companion star and
the envelope (MacLeod et al., 2017). The plunge-in speed $v_{p}$ is comparable
to the Keplerian speed at the plunge-in radius $r_{p}$, but much larger than
the envelope’s speed $v_{\rm{env}}$. On the other hand, the kinetic energy
required to produce $E_{r}$ and $e_{\text{g}}$ in the ejecta, converted to
speed, can be calculated by,

$v_{\rm{source}}=\sqrt{\frac{2(E_{r}+e_{\text{g}})}{\rho}},$ (26)

where $\rho$ is the ejecta’s density. Due to energy conversion,

$v_{p}^{2}\approx\frac{GM_{\star}}{r_{p}}>v_{\rm{env}}^{2}+v_{\rm{source}}^{2}.$
(27)

The maximum $v_{\rm{source}}$ in Figure 2 is 448km$\cdot$s-1, which means that
$r_{p}<5.69$R⊙. More sophisticated stellar and binary evolution analysis
should be considered to further narrow down $r_{p}$ and $M_{\star}$ (Ge et
al., 2010, 2015, 2020).

We carry out simulations of two sub-models: the shock model and the shock-free
model. The shock model has an ejecta with a slightly increasing speed during
the late stage and the shock-free model has an ejecta with a slightly
decreasing speed during the late stage. Figure 2 shows the time-dependent
variables of the two sub-models. The functions that generate these variables
can be found in Appendix C.

Figure 2: The solid and dashed lines are the time-dependent inner boundary
conditions of the shock and shock-free models, respectively. The red line
indicates the escape velocity at the inner boundary. The orange region denotes
the radiation-dominated ejecta. On the left panel, from top to bottom, each
plot shows the velocity of the ejecta, $v_{\rm{source}}$ calculated from
Equation 26, and the temperature of the ejecta. On the right panel, from the
top to the bottom, they are the mass loss rate, cumulative mass loss, and
$E_{r}/e_{\text{g}}$, respectively. The shock model and shock-free model also
differ in the duration of the ejecta.

Figure 3 shows the light curves of the shock and shock-free models. We have
shifted the time axis to let the peak be located at $t=0$. We can see that the
light curves of the shock and shock-free models both resemble AT2019zhd.
Consequently, the time evolution of the shock-free model is similar to the
shock model, therefore, we just show the time evolution of the shock model in
Figure 4 to save space. In particular, we calculate the optical depth of the
radiation flux by,

$\tau_{R}(r)=\int_{r}^{r_{\rm{out}}}\rho\kappa_{\rm{R}}dr.$ (28)

The time axis is adjusted to be consistent with the light curves. The time
evolution of H and He species can be found in Appendix D, since they are
determined by $\rho$ and $T_{\text{g}}$.

Figure 3: The solid line with black dots is the observed light curve of AT
2019zhd (Pastorello et al., 2021a), and the red and blue solid lines are the
light curves of the shock and the shock-free models, respectively. Figure 4:
The time evolution of the shock model. In each panel, the x-axis is time in
days, and the y-axis is radius in 1000R⊙. From panel 1 to 6, they are
$\log_{10}\rho$ with contours of $[-10,-11,-12]$, $\log_{10}T_{\text{g}}$ with
contours of $[3.1,3.7]$, and $\log_{10}(a_{\rm{rad}}/g)$ with contours of
$[-2,0]$, the optical depth $\log_{10}\tau_{R}$ with contours of $[0,1,3]$,
$\log_{10}(E_{r}/e_{\text{g}})$ with contours of $[-1,0,1]$, and
$\log_{10}\kappa_{\rm{R}}$ with contours of $[-2,0]$.

There are several prominent results in the time evolution plots, see Figure 5.

  1. 1.

In panel 2, the temperature profile has a spike near the peak of the
luminosity. This is when the radiation-dominated ejecta becomes optically thin
and the radiation heats the surrounding gas.

  2. 2.

The observationally important ejecta material is the one that has a
temperature of about 5000K. The ejecta’s matter that is colder than 5000K is
closer to the observer but is mostly transparent, revealing behind it the
ejecta’s matter with a temperature of about 5000K. The time-evolution of the
surface $\log_{10}T_{\text{g}}=3.7$ (see panel 2) resembles the shape of the
light curve in Figure 3.

  3. 3.

Panel 3 shows that the $a_{\rm{rad}}$ exceeds the $g$ at small radii during
the first 30 days, and at large radii during the late stage. This phenomenon
is consistent with the evolution of $\kappa_{\rm{R}}$. When $10^{4}\lesssim
T_{\text{g}}\lesssim 10^{5}$K, most of the gas is ionized and
$\kappa_{\rm{R}}$ is large, radiation flux can be significant. When
$T_{\text{g}}<1400$K, dust may form and provide the opacity, and the
$a_{\rm{rad}}$ becomes significant again. At this stage, the LRNe is similar
to the asymptotic-giant-branch (AGB) stars (Höfner & Olofsson, 2018).

  4. 4.

In panel 4, the radius of $\tau_{R}=1$ expands rapidly when the inner region
of the ejecta cools off because dust can form at low temperatures and obscure
the object.

  5. 5.

The inner region ($r<2000$R⊙) transits from radiation-dominated to matter-
dominated as the ejecta expands into the ambient and cools off. It can be seen
in panel 5 of Figure 4.

Our model is comprehensive in the sense that we incorporate many physical
processes. To address the impact of each specific physics on the light curve,
we run some companion simulations by turning off the corresponding physics in
Appendix E, F, and G. We show that our calculations are converged in Appendix
H.

##  4 Conclusions

We present a 1D radiation hydrodynamic model that can obtain light curves of
LRNe events. For the first time, our model incorporates recombination physics,
radiation transport, and radiation force acceleration of a CEE into
consideration. There are several key new physics that we determined to be
intrinsic to CEE and LRNe.

First, we find that the light curve’s peak and plateau are formed differently.
Specifically, the peak is formed primarily by radiation-dominated
($E_{r}/e_{\text{g}}\gg 1$), and the plateau is formed by matter-dominated
ejecta ($E_{r}/e_{\text{g}}<1$).

Secondly, radiation force plays a significant role in driving the ejecta
outwards (panel 3 of Figure 4). This occurs where opacity is high,
specifically when both H and He are still ionized. Radiation force may help to
eject matter that initially has sub-escape velocity.

At last, rapid dust formation, assuming it occurs, absorbs radiation and
accelerates during the late LRNe stages, resembling environments near AGB
stars. While the radiation force helps to accelerate the dusty matter, this
has minimal impact on the peak and first plateau of the light curve.

There are several limitations in our model right now. For example, the
temperature of the ejecta is described by Equation 21. Thus, we can not
strictly distinguish the effect of changing $M_{\star}/r_{\rm{in}}$ and
$\alpha$. However, one can constrain $M_{\star}$ by a scaling relation, for
example, Equation 25, and can constrain $r_{\rm{in}}$ (which is affected by
$r_{s}$) with long-term binary evolution models. Moreover, one can further
constrain the thermalization efficiency $\alpha$ with realistic 3D radiation
hydrodynamic calculations.

The effect of changing initial $e_{\text{g}}$ and $E_{r}$ in the ejecta is
also not fully independent, but radiation transport can, to some extent, help
to distinguish the radiation-dominated cases from the matter-dominated cases.

We also find that shock and shock-free models can produce similar results. In
our models, shock and shock-free models refer to whether the speed of the
ejecta was increasing or decreasing with time. Indeed, a shock can generate
heat and prolong the plateau phase, but adding mass and internal energy to the
ejecta can also result in a prolonged plateau in the light curve. We may need
other observables, such as H$\alpha$ emission, to distinguish these two
scenarios. We can also explore the effect of pre-CEE asymmetric matter
distribution around the central binary object; this will be the subject of the
follow-up 2D work.

Even with a 1D model, we can fit the AT2019zhd light curve reasonably well.
Motivated by observations and scaling relations (Section 3.2), we adopted the
mass of the central object to be $M_{\star}=6M_{\odot}$, and show that our
choice of $M_{\star}$ and $r_{\rm{in}}$ are reasonable. With this choice, we
estimate that the mass of the ejecta could be $\Delta
M\in[0.04,0.1]M_{\odot}$111A different choice of the central object mass will
lead to different inferred CE ejecta properties if such a fit is found. While
our model fits the observations considering reasonable assumptions about
ejecta’s density, internal and kinetic energies, we do not exclude a
possibility of other fits with CE ejecta properties being far outside of
considered in this paper.. As a result, the 1D model that we proposed in this
Letter can provide a crucial bridge between CEE theory and the observables of
LRNe.

##

We thank the referee for the careful review and the constructive suggestions
that improve the quality of this Letter. The authors thank Andrea Pastorello
for sharing AT2019zhd’s light curve data with us. Z.C. enjoyed discussing with
Xuening Bai (白雪宁), Shude Mao (毛淑德), and Wei Zhu (祝伟). Z.C. is grateful for the
support from the National Natural Science Foundation of China (NSFC
No.12103028, No.12342501), Tsinghua University Dushi Program, and Shuimu
Tsinghua Scholar Program. N.I. acknowledges funding from NSERC Discovery under
Grant No. NSERC RGPIN-2019-04277.

## References

  * Balay et al. (1997) Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. 1997, in Modern Software Tools in Scientific Computing, ed. E. Arge, A. M. Bruaset, & H. P. Langtangen (Birkhäuser Press), 163–202
  * Balay et al. (2019) Balay, S., Abhyankar, S., Adams, M. F., et al. 2019, PETSc Users Manual, Tech. Rep. ANL-95/11 - Revision 3.11, Argonne National Laboratory
  * Banerjee et al. (2015) Banerjee, D. P. K., Nuth, Joseph A., I., Misselt, K. A., et al. 2015, ApJ, 814, 109, doi: 10.1088/0004-637X/814/2/109
  * Bermúdez-Bustamante et al. (2024) Bermúdez-Bustamante, L. C., De Marco, O., Siess, L., et al. 2024, arXiv e-prints, arXiv:2401.03644, doi: 10.48550/arXiv.2401.03644
  * Blagorodnova et al. (2017) Blagorodnova, N., Kotak, R., Polshaw, J., et al. 2017, ApJ, 834, 107, doi: 10.3847/1538-4357/834/2/107
  * Blagorodnova et al. (2020) Blagorodnova, N., Karambelkar, V., Adams, S. M., et al. 2020, MNRAS, 496, 5503, doi: 10.1093/mnras/staa1872
  * Blagorodnova et al. (2021) Blagorodnova, N., Klencki, J., Pejcha, O., et al. 2021, A&A, 653, A134, doi: 10.1051/0004-6361/202140525
  * Bronner et al. (2023) Bronner, V. A., Schneider, F. R. N., Podsiadlowski, P., & Roepke, F. K. 2023, arXiv e-prints, arXiv:2311.06332, doi: 10.48550/arXiv.2311.06332
  * Cai et al. (2022a) Cai, Y., Reguitti, A., Valerin, G., & Wang, X. 2022a, Universe, 8, 493, doi: 10.3390/universe8100493
  * Cai et al. (2019) Cai, Y. Z., Pastorello, A., Fraser, M., et al. 2019, A&A, 632, L6, doi: 10.1051/0004-6361/201936749
  * Cai et al. (2022b) —. 2022b, A&A, 667, A4, doi: 10.1051/0004-6361/202244393
  * Chamandy et al. (2024) Chamandy, L., Carroll-Nellenback, J., Blackman, E. G., et al. 2024, MNRAS, 528, 234, doi: 10.1093/mnras/stae036
  * Chamandy et al. (2018) Chamandy, L., Frank, A., Blackman, E. G., et al. 2018, MNRAS, 480, 1898, doi: 10.1093/mnras/sty1950
  * Chen et al. (2024) Chen, X., Liu, Z., & Han, Z. 2024, Progress in Particle and Nuclear Physics, 134, 104083, doi: 10.1016/j.ppnp.2023.104083
  * Chen & Bai (in prep) Chen, Z., & Bai, X. in prep
  * Chen et al. (2019) Chen, Z., Coleman, M. S. B., Blackman, E. G., & Frank, A. 2019, Journal of Computational Physics, 388, 490, doi: 10.1016/j.jcp.2019.03.016
  * Ge et al. (2010) Ge, H., Hjellming, M. S., Webbink, R. F., Chen, X., & Han, Z. 2010, ApJ, 717, 724, doi: 10.1088/0004-637X/717/2/724
  * Ge et al. (2015) Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2015, ApJ, 812, 40, doi: 10.1088/0004-637X/812/1/40
  * Ge et al. (2020) —. 2020, ApJ, 899, 132, doi: 10.3847/1538-4357/aba7b7
  * González-Bolívar et al. (2023) González-Bolívar, M., De Marco, O., Bermúdez-Bustamante, L. C., Siess, L., & Price, D. J. 2023, MNRAS, doi: 10.1093/mnras/stad3748
  * Hatfull et al. (2021) Hatfull, R. W. M., Ivanova, N., & Lombardi, J. C. 2021, MNRAS, 507, 385, doi: 10.1093/mnras/stab2140
  * Höfner & Olofsson (2018) Höfner, S., & Olofsson, H. 2018, A&A Rev., 26, 1, doi: 10.1007/s00159-017-0106-5
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iaconi et al. (2019) Iaconi, R., Maeda, K., De Marco, O., Nozawa, T., & Reichardt, T. 2019, MNRAS, 489, 3334, doi: 10.1093/mnras/stz2312
  * Iaconi et al. (2020) Iaconi, R., Maeda, K., Nozawa, T., De Marco, O., & Reichardt, T. 2020, MNRAS, 497, 3166, doi: 10.1093/mnras/staa2169
  * Ivanova et al. (2013a) Ivanova, N., Justham, S., Avendano Nandez, J. L., & Lombardi, J. C. 2013a, Science, 339, 433, doi: 10.1126/science.1225540
  * Ivanova et al. (2015) Ivanova, N., Justham, S., & Podsiadlowski, P. 2015, MNRAS, 447, 2181, doi: 10.1093/mnras/stu2582
  * Ivanova & Nandez (2016) Ivanova, N., & Nandez, J. L. A. 2016, MNRAS, 462, 362, doi: 10.1093/mnras/stw1676
  * Ivanova et al. (2013b) Ivanova, N., Justham, S., Chen, X., et al. 2013b, A&A Rev., 21, 59, doi: 10.1007/s00159-013-0059-2
  * Kamiński et al. (2015) Kamiński, T., Mason, E., Tylenda, R., & Schmidt, M. R. 2015, A&A, 580, A34, doi: 10.1051/0004-6361/201526212
  * Kolb et al. (2013) Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80, doi: 10.1051/0004-6361/201321499
  * Kurtenkov et al. (2015) Kurtenkov, A. A., Pessev, P., Tomov, T., et al. 2015, A&A, 578, L10, doi: 10.1051/0004-6361/201526564
  * Levermore & Pomraning (1981) Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321, doi: 10.1086/159157
  * Liu et al. (2023) Liu, Z.-W., Röpke, F. K., & Han, Z. 2023, Research in Astronomy and Astrophysics, 23, 082001, doi: 10.1088/1674-4527/acd89e
  * MacLeod et al. (2017) MacLeod, M., Antoni, A., Murguia-Berthier, A., Macias, P., & Ramirez-Ruiz, E. 2017, ApJ, 838, 56, doi: 10.3847/1538-4357/aa6117
  * MacLeod et al. (2022) MacLeod, M., De, K., & Loeb, A. 2022, ApJ, 937, 96, doi: 10.3847/1538-4357/ac8c31
  * Malygin et al. (2014) Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91, doi: 10.1051/0004-6361/201423768
  * Matsumoto & Metzger (2022) Matsumoto, T., & Metzger, B. D. 2022, ApJ, 938, 5, doi: 10.3847/1538-4357/ac6269
  * Metzger & Pejcha (2017) Metzger, B. D., & Pejcha, O. 2017, MNRAS, 471, 3200, doi: 10.1093/mnras/stx1768
  * Munari et al. (2002) Munari, U., Henden, A., Kiyota, S., et al. 2002, A&A, 389, L51, doi: 10.1051/0004-6361:20020715
  * Nandez & Ivanova (2016) Nandez, J. L. A., & Ivanova, N. 2016, MNRAS, 460, 3992, doi: 10.1093/mnras/stw1266
  * Nandez et al. (2014) Nandez, J. L. A., Ivanova, N., & Lombardi, J. C., J. 2014, ApJ, 786, 39, doi: 10.1088/0004-637X/786/1/39
  * Nandez et al. (2015) Nandez, J. L. A., Ivanova, N., & Lombardi, J. C. J. 2015, MNRAS, 450, L39, doi: 10.1093/mnrasl/slv043
  * Nicholls et al. (2013) Nicholls, C. P., Melis, C., Soszynski, I., et al. 2013, MNRAS, 431, L33, doi: 10.1093/mnrasl/slt003
  * O’Connor et al. (2023) O’Connor, C. E., Bildsten, L., Cantiello, M., & Lai, D. 2023, ApJ, 950, 128, doi: 10.3847/1538-4357/acd2d4
  * Ohlmann et al. (2016) Ohlmann, S. T., Röpke, F. K., Pakmor, R., & Springel, V. 2016, ApJ, 816, L9, doi: 10.3847/2041-8205/816/1/L9
  * Omukai et al. (2005) Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627, doi: 10.1086/429955
  * Pastorello et al. (2019) Pastorello, A., Mason, E., Taubenberger, S., et al. 2019, A&A, 630, A75, doi: 10.1051/0004-6361/201935999
  * Pastorello et al. (2021a) Pastorello, A., Fraser, M., Valerin, G., et al. 2021a, A&A, 646, A119, doi: 10.1051/0004-6361/202039952
  * Pastorello et al. (2021b) Pastorello, A., Valerin, G., Fraser, M., et al. 2021b, A&A, 647, A93, doi: 10.1051/0004-6361/202039953
  * Pejcha et al. (2017) Pejcha, O., Metzger, B. D., Tyles, J. G., & Tomida, K. 2017, ApJ, 850, 59, doi: 10.3847/1538-4357/aa95b9
  * Prust & Chang (2019) Prust, L. J., & Chang, P. 2019, MNRAS, 486, 5809, doi: 10.1093/mnras/stz1219
  * Reichardt et al. (2019) Reichardt, T. A., De Marco, O., Iaconi, R., Tout, C. A., & Price, D. J. 2019, MNRAS, 484, 631, doi: 10.1093/mnras/sty3485
  * Reig (2011) Reig, P. 2011, Ap&SS, 332, 1, doi: 10.1007/s10509-010-0575-8
  * Renzo et al. (2021) Renzo, M., Callister, T., Chatziioannou, K., et al. 2021, ApJ, 919, 128, doi: 10.3847/1538-4357/ac1110
  * Röpke & De Marco (2023) Röpke, F. K., & De Marco, O. 2023, Living Reviews in Computational Astrophysics, 9, 2, doi: 10.1007/s41115-023-00017-x
  * Sand et al. (2020) Sand, C., Ohlmann, S. T., Schneider, F. R. N., Pakmor, R., & Röpke, F. K. 2020, A&A, 644, A60, doi: 10.1051/0004-6361/202038992
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611, doi: 10.1051/0004-6361:20031279
  * Soker et al. (2018) Soker, N., Grichener, A., & Sabach, E. 2018, ApJ, 863, L14, doi: 10.3847/2041-8213/aad736
  * Stritzinger et al. (2020) Stritzinger, M. D., Taddia, F., Fraser, M., et al. 2020, A&A, 639, A104, doi: 10.1051/0004-6361/202038019
  * Tomida et al. (2013) Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6, doi: 10.1088/0004-637X/763/1/6
  * Tylenda et al. (2011) Tylenda, R., Hajduk, M., Kamiński, T., et al. 2011, A&A, 528, A114, doi: 10.1051/0004-6361/201016221
  * Wakelam et al. (2017) Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Molecular Astrophysics, 9, 1, doi: 10.1016/j.molap.2017.11.001
  * Wang & Han (2012) Wang, B., & Han, Z. 2012, New A Rev., 56, 122, doi: 10.1016/j.newar.2012.04.001
  * Warner (2003) Warner, B. 2003, Cataclysmic variable stars, Vol. 28 (Cambridge University Press)

##  Appendix A Equation of state

Realistic EoS can be highly nonlinear. Here, we present a simplified EoS model
of H and He mixture, i.e., the sum of the mass fraction of hydrogen element
$X$ and helium element $Y$ equals 1. This model considers the latent heat of
hydrogen dissociation and all the ionization transitions of the H and He. To
simplify the calculation, we assume that H, He, and e- have the same
temperature, but H and He have independent thermodynamic subsystems. This
thermodynamic system is not as complete as in (Tomida et al., 2013) but is
much easier to solve.

For a thermodynamic system with multi-species, its entropy $s$, pressure $p$,
internal energy $e_{\text{g}}$, and heat capacity $C_{V}$ are,

$\displaystyle s$ $\displaystyle=$ $\displaystyle
k_{b}\sum_{i}n_{i}\bigg{[}1+\frac{d\ln z_{i}}{d\ln
T_{\text{g}}}-\ln\bigg(\frac{n_{i}}{z_{i}}\bigg{missing})\bigg{]},$ (A1)
$\displaystyle p$ $\displaystyle=$
$\displaystyle\sum_{i}n_{i}k_{b}T_{\text{g}},$ (A2) $\displaystyle
e_{\text{g}}$ $\displaystyle=$ $\displaystyle\sum_{i}n_{i}\xi_{i},$ (A3)
$\displaystyle C_{V}$ $\displaystyle=$
$\displaystyle\sum_{i}\partialderivative{n_{i}}{T_{\text{g}}}\xi_{i}+\sum_{i}n_{i}\partialderivative{\xi_{i}}{T_{\text{g}}},$
(A4)

where $i$ is the summation index of all the species, $n_{i}$, $\xi_{i}$, and
$z_{i}$ are the number density, particle energy, and partition function of
species $i$. The particle energy of hydrogen and helium species, in the regime
of the ideal gas law, are described by,

$\xi_{i}=3k_{b}T_{\text{g}}/2+\phi_{i},$ (A5)

where $\phi_{i}$ is each particle’s dissociation/ionization binding energy,
which is also the source of the latent heat. For example, the ground state of
He is chosen to be $\phi_{\ce{He}}=0$, and the binding energy of He^2+ is
$\phi_{\ce{He^{2}+}}=24.59+54.42=80.01$eV. The partition functions, not
considering the atomic structure, are described by,

$z_{i}=\bigg{(}\frac{2\pi
m_{i}k_{b}T_{\text{g}}}{h^{2}}\bigg{)}^{3/2}\exp(\frac{-\phi_{i}}{k_{b}T_{\text{g}}}),$
(A6)

where $m_{i}$ is the mass of each particle species and $h$ is the Planck
constant.

We also need the adiabatic sound speed $c_{s}$ and adiabatic index $\gamma$
for the Riemann solver (Chen et al., 2019), which is calculated by,

$\displaystyle c_{s}$ $\displaystyle=$
$\displaystyle\sqrt{\bigg{(}\partialderivative{p}{\rho}\bigg{)}_{T}+\bigg{(}\partialderivative{p}{T}\bigg{)}_{\rho}\bigg{(}\partialderivative{T}{\rho}\bigg{)}_{s}}$
(A7) $\displaystyle=$
$\displaystyle\sqrt{\bigg{(}\partialderivative{p}{\rho}\bigg{)}_{T}-\bigg{(}\partialderivative{p}{T}\bigg{)}_{\rho}\frac{(\partial
s/\partial\rho)_{T}}{(\partial s/\partial T)_{\rho}}},$ $\displaystyle\gamma$
$\displaystyle=$ $\displaystyle\rho c_{s}^{2}/p$ (A8)

Thus, to calculate $c_{s}$, we need to know $(\partial s/\partial\rho)_{T}$
and $(\partial s/\partial\rho)_{\rho}$, which means $(\partial
n_{i}/\partial\rho)_{T}$ and $(\partial n_{i}/\partial\rho)_{\rho}$. They can
be calculated by thermodynamic equilibrium equations, i.e., the Saha
equations. To further simplify the EoS model, we divide the $\rho-T$ diagram
into pure and mixture states (Chen et al., 2019). We define pure states if the
mass fraction of a single species (not including electron) exceeds $0.99999$
and all other cases are mixture states. Figure 5 shows the division of H and
He. The blue color represents the pure states of H and He. In doing so, we can
solve fewer Saha equations.

Figure 5: The $\rho$-$T$ diagram of H and He are divided into 6 zones to
improve the computational efficiency. The pure states are colored with blue,
and the rest are mixture states that require solving Saha equations.

As a first example, let us consider hydrogen’s ionization and dissociation
reactions (the red color zone in the first panel of Figure 5).

$\displaystyle\ce{H}$ $\displaystyle\longleftrightarrow$
$\displaystyle\ce{H+}+\ce{e-},$ (A9) $\displaystyle\ce{H2}$
$\displaystyle\longleftrightarrow$ $\displaystyle 2\ce{H}.$ (A10)

The corresponding Saha equations are

$\displaystyle\frac{n_{\ce{H}}^{2}}{n_{\ce{H2}}}$ $\displaystyle=$
$\displaystyle\frac{z_{\ce{H}}^{2}}{z_{\ce{H2}}}=q_{1},$ (A11)
$\displaystyle\frac{n_{\ce{H+}}^{2}}{n_{\ce{H}}}$ $\displaystyle=$
$\displaystyle\frac{z_{\ce{H+}}z_{\ce{e-}}}{z_{\ce{H}}}=q_{2},$ (A12)

where we have used the charge neutrality equation

$\displaystyle n_{\ce{H+}}=n_{\ce{e-}}$ (A13)

We also have the number conservation equation,

$\displaystyle 2n_{\ce{H2}}+n_{\ce{H}}+n_{\ce{H+}}=X\rho/m_{\ce{H}}=\bar{X},$
(A14)

where $X$ is the mass fraction of $\ce{H}$ of the gas. Substitute the Saha
equations into the number conservation equation, and we get a quadratic
equation

$\displaystyle\frac{2n_{\ce{H+}}^{4}}{q_{1}^{2}q_{2}}+\frac{n_{\ce{H+}}^{2}}{q_{2}}+n_{\ce{H+}}=\bar{X}.$
(A15)

The solution of Equation A15 is $n_{\ce{H+}}$, and $n_{\ce{H}}$ and
$n_{\ce{H2}}$ can be calculated from the Saha equations.

We use the differential form of Equation A11, A12, and A14 to obtain
$(\partial n_{i}/\partial\rho)_{T}$ and $(\partial n_{i}/\partial\rho)_{\rho}$
of H.

$\displaystyle\begin{pmatrix}-1&2&0\\\ 0&-1&2\\\
2n_{\ce{H2}}&n_{\ce{H}}&n_{\ce{H+}}\end{pmatrix}\begin{pmatrix}d\ln
n_{\ce{H}}\\\ d\ln n_{\ce{H+}}\\\ d\ln
n_{\ce{H2}}\end{pmatrix}=\begin{pmatrix}d\ln q_{1}\\\ d\ln q_{2}\\\
\bar{X}d\ln\bar{X}\end{pmatrix}$ (A16)

If only one reaction is present, the thermodynamic system is trivial to solve.

For He, let us consider the two reactions (the red color zone in the second
panel of Figure 5)

$\displaystyle\ce{He}$ $\displaystyle\longleftrightarrow$
$\displaystyle\ce{He+}+\ce{e-},$ (A17) $\displaystyle\ce{He+}$
$\displaystyle\longleftrightarrow$ $\displaystyle\ce{He^{2+}}+\ce{e-}$ (A18)

The corresponding Saha equations are

$\displaystyle\frac{n_{\ce{He+}}n_{\ce{e-}}}{n_{\ce{He}}}$ $\displaystyle=$
$\displaystyle\frac{z_{\ce{He+}}z_{\ce{e-}}}{z_{\ce{He}}}=q_{3}$ (A19)
$\displaystyle\frac{n_{\ce{He^{2+}}}n_{\ce{e-}}}{n_{\ce{He+}}}$
$\displaystyle=$
$\displaystyle\frac{z_{\ce{He^{2+}}}z_{\ce{e-}}}{z_{\ce{He+}}}=q_{4}$ (A20)

The number conservation equation and charge neutrality equations are,

$\displaystyle n_{\ce{He}}+n_{\ce{He+}}+n_{\ce{He^{2+}}}$ $\displaystyle=$
$\displaystyle Y\rho/m_{\ce{He}}=\bar{Y}.$ (A21) $\displaystyle
n_{\ce{He+}}+2n_{\ce{He^{2+}}}$ $\displaystyle=$ $\displaystyle n_{\ce{e-}}$
(A22)

Substitute Equation A19, A20 and A22 into A21, we get a cubic equation,

$\displaystyle
n_{\ce{He^{2+}}}^{3}(q_{3}-4q_{4})+n_{\ce{He^{2+}}}^{2}q_{1}\bar{Y}$
$\displaystyle-$ $\displaystyle n_{\ce{He^{2+}}}q_{3}q_{4}(q_{4}+2\bar{Y})$
(A23) $\displaystyle+$ $\displaystyle q_{3}q_{4}^{2}\bar{Y}=0$

Solve Equation A23 gives us $n_{\ce{He^{2+}}}$, then $n_{\ce{He+}}$ and
$n_{\ce{e-}}$ can be solved by Equation A20 and A22, and $n_{\ce{He}}$ is
calculated from A21. Meanwhile, we write down the differential form of
Equation A19, A20, A21, and A22 and calculate $(\partial
n_{i}/\partial\rho)_{T}$ and $(\partial n_{i}/\partial\rho)_{\rho}$ of helium
species.

$\displaystyle\begin{pmatrix}-1&1&0&1\\\ 0&-1&1&1\\\
n_{\ce{He}}&n_{\ce{He+}}&n_{\ce{He^{2+}}}&0\\\
0&n_{\ce{He+}}&n_{\ce{He^{2+}}}&-n_{\ce{e-}}\end{pmatrix}\begin{pmatrix}d\ln
n_{\ce{He}}\\\ d\ln n_{\ce{He+}}\\\ d\ln n_{\ce{He2+}}\\\ d\ln
n_{\ce{e-}}\end{pmatrix}=$ $\displaystyle\begin{pmatrix}d\ln q_{3}\\\ d\ln
q_{4}\\\ \bar{Y}d\ln\bar{Y}\\\ 0\end{pmatrix}$ (A24)

If only one reaction is present, the thermodynamic system is trivial to solve.

Figure 6: The first panel shows the combined specific heat capacity of
$X=0.74$ and $Y=0.26$, and the bottom panel shows the adiabatic index.

Finally, we can calculate all the thermodynamic variables. Figure 6 shows the
specific heat capacity and adiabatic index we use in this Letter.

##  Appendix B Opacity

The opacity table combines the gas opacity table of Malygin et al. (2014) that
dominates over $1500$K, and the dust opacity table (1% of dust to gas ratio)
of Semenov et al. (2003) that dominates below 1100-1200K, depending on the
density. For temperatures in between, the maximum of the value of the two
tables is taken. Figure 7 shows $\kappa_{\rm{R}}$ at different densities and
$\kappa_{\rm{P}}$ is similar.

Figure 7: The Rosseland mean opacity that includes dust opacity at different
densities. Solid and dashed lines denote the gas and dust opacities,
respectively. The red zone indicates the dust formation temperature, the gray
zone indicates the temperature of hydrogen dissociation and H- formation, and
the orange zone indicates the ionization temperature of H and He.

##  Appendix C time-dependent boundary conditions

When $t<t_{1}$, the ejecta is radiation-dominated,

$\displaystyle\dot{M}=\dot{M}_{1}+\frac{(\dot{M}_{2}-\dot{M}_{1})\exp{x}}{\exp{x}+1}+\dot{M}_{3}\bigg{(}\frac{t}{t_{1}}\bigg{)}^{2}$
(C1) $\displaystyle
f_{\rm{ej}}=\begin{cases}1+\frac{(v_{1}-1)t}{x_{2}t_{1}},\quad\quad\quad\quad\quad\quad\quad\quad
t\leq x_{2}t_{1}\\\
1+(v_{3}-1)\exp{-\big{(}\frac{x_{1}(t-x_{2}t_{1})}{t_{1}}\big{)}^{2}},\text{else}\end{cases}$
(C2)

When $t_{1}\leq t\leq t_{1}+t_{2}$, the ejecta becomes matter-dominated,

$\displaystyle\dot{M}$ $\displaystyle=$
$\displaystyle\dot{M}_{2}-\frac{(\dot{M}_{2}-\dot{M}_{1})\exp{-x}}{\exp{-x}+1}-\delta\dot{M}\frac{t-t_{1}}{t_{2}}$
(C3) $\displaystyle f_{\rm{ej}}$ $\displaystyle=$ $\displaystyle
v_{2}+\big{(}1+(v_{3}-1)\exp{-(x_{1}(1-x_{2}))^{2}}-v_{2}\big{)}$ (C4)
$\displaystyle\exp{\frac{t-t_{1}}{dt_{2}}}+\delta
v\big{(}\frac{t-t_{1}}{t_{2}}\big{)},$

where

$\displaystyle v_{\rm{esc}}$ $\displaystyle=$
$\displaystyle\sqrt{2GM_{\star}/r_{\rm{in}}},$ (C5)
$\displaystyle\delta\dot{M}$ $\displaystyle=$
$\displaystyle(\dot{M}_{2}-\dot{M}_{1}),$ (C6) $\displaystyle x$
$\displaystyle=$ $\displaystyle 20(t-t_{1})/t_{1}.$ (C7)

There are 13 parameters in this model, we list all the parameters just to show
the functional form of the ejecta. In terms of the shock and shock-free
models, they have many parameters in common. The common parameters are:
$t_{1}=1.3$ days, $dt_{2}=1.6$ days, $v_{1}=1.06$,
$\dot{M}_{1}=0.16$M${}_{\odot}\cdot$yr-1,
$\dot{M}_{3}=0.2$M${}_{\odot}\cdot$yr-1, $x_{1}=6$, and $x_{2}=0.55$. We list
the different parameters in Table2.

| $t2$ | $v2$ | $v3$ | $\delta v$ | $\dot{M}_{2}$ | $\delta\dot{M}$
---|---|---|---|---|---|---
| [day] |  |  |  | [M${}_{\odot}\cdot$yr-1] | [M${}_{\odot}\cdot$yr-1]
s | 8 | 0.87 | 0.94 | 0.01 | 1.3 | 0
sf | 20 | 0.875 | 0.955 | -0.04 | 1.8 | -0.5
Table 2: The parameters that are different in shock (s) and shock-free (sf)
models. In particular, $\delta v>0$ cases generally belong to shock models and
$\delta v<0$ cases belong to shock-free models.

We anticipate that there could be other fitting functions with more/less
parameters. One should explore different functional forms and a thorough
parameter space study is deemed appropriate in the future.

##  Appendix D The time evolution of hydrogen and helium species

Figure 8 shows the time evolution of the mass fraction of H and He species of
the shock model. For example, the mass fraction of H2 is defined as,

$\chi_{\ce{H2}}=\frac{2n_{\ce{H2}}}{2n_{\ce{H2}}+n_{\ce{H}}+n_{\ce{H+}}}$ (D1)

Figure 8: From the left to the right, the first row shows the time evolution
of the mass fraction of He, He+, and He^2+, the second row shows the time
evolution of the mass fraction of H2, H, and H+.

Although the mass fraction of each species can be derived from $\rho$ and
$T_{\text{g}}$ already, there are two features we would like to clarify.

  1. 1.

H2 forms in the late stage. There are two channels for H2 formations. At high
density, H can form H2 directly (Omukai et al., 2005). When dust is present,
H2 can form on dust (Wakelam et al., 2017). In the time evolution, we find
that the H2 formation region largely coincide with the dust formation region
(can be seen from the $\log_{10}\kappa_{\rm{R}}$ panel in Figure 4) or with
the high-density region (can be seen from the $\log_{10}\rho$ panel in Figure
4). Therefore, Saha’s EoS model is suitable for this problem.

  2. 2.

There is a line showing the partial ionization of H+ in the H+ panel. This is
due to Zeldovich’s spikes in the radiation hydrodynamic calculation after a
shock. The shock is caused by the fast ejecta colliding with the ambient.
However, we tend not to ascribe the H$\alpha$ emission to this line because
the decreasing width of H$\alpha$ in observations (Pastorello et al., 2021a)
indicates that the shock decelerates. We believe that the H$\alpha$ emission
is caused by a marginally escaping ejecta colliding with the pre-existing
decreation disk (Pejcha et al., 2017; Metzger & Pejcha, 2017). However, the 1D
model is only shell-like and cannot model disks.

##  Appendix E The impact of realistic EoS

Figure 9: The upper and lower panels show the light curves of H and He mixture
EoS and perfect gas EoS with various $\gamma$ of the shock (s) and shock-free
(sf) models, respectively.

EoS can play an important role in radiation hydrodynamics. We compare our
simulations with realistic EoS to the widely used in hydrodynamic codes
$\gamma$-law EoS, described by

$\displaystyle e_{\text{g}}$ $\displaystyle=$ $\displaystyle\rho
k_{b}T_{\text{g}}/(\gamma-1),$ (E1) $\displaystyle p$ $\displaystyle=$
$\displaystyle\rho k_{b}T_{\text{g}}/(\mu m_{H}),$ (E2)

where $\mu=1$ is the mean molecular weight.

With $\gamma$-law EoS, gas energy does not contain latent heat. As a result,
the cooling timescale is shorter, and all the light curves ran with
$\gamma$-law EoS last for a shorter time. We can also see that the peak is
less affected by the change in the EoS because the radiation-dominated ejecta
governs the peak. The internal energy $e_{\text{g}}$ plays a minor role in
this stage.

##  Appendix F The impact of radiation force

Figure 10: The black line with dots shows the light curve of AT2019zhd. The
red and blue colors indicate the shock (s) and shock-free (sf) models. The
solid and dashed lines indicate our original models and their counterparts
when we turn off the radiation force.

Radiation force was largely ignored in previous hydrodynamical simulations of
CEE. However, we find that the $a_{\rm{rad}}$ is larger than $g$ at high
($>10000$K) and low temperatures ($<1400$K); see panel 3 in Figure 4. To
demonstrate the importance of radiation force, we compare our original
simulations’ light curves to those without radiation force ($a_{\rm{rad}}=0$).
The results are presented in Figure 10.

In the case of no radiation force in the momentum equation, the luminosity
drops off faster. This faster drop is due to more material in the ejecta
falling back to the central object under the influence of the gravitational
force. The drop happens so early that the gas is still fully ionized and
should be pushed away by the radiation force. In our model, the fallback
material will pass the inner boundary and no longer contribute to the
luminosity.

##  Appendix G The impact of dust formation

Figure 11: Upper panel: light curves of the last 10 days of shock (s) and
shock-free (sf) models and their corresponding dust deficit models. Bottom
panel: the time evolution of $\log_{10}(\kappa_{\rm{R}}{F}_{r}/(cg))$ and
$\log_{10}\kappa_{\rm{R}}$ of the shock dust deficit model.

It is argued that molecules and dust form during the late stage of CEE
(Nicholls et al., 2013; Kamiński et al., 2015; Banerjee et al., 2015;
Blagorodnova et al., 2020; Iaconi et al., 2020; MacLeod et al., 2022).
Numerical models also show that dust may provide an additional push away from
the binary merger due to radiation acceleration (González-Bolívar et al.,
2023). To test the impact of dust formation, we modify the opacity table
(Figure 7) by decreasing $\kappa_{\rm{R}}$ and $\kappa_{\rm{P}}$ of dust by
100 times and name it the dust deficit model. We find that modification of
dust opacities does now affect light curves during the peak and the plateau.
This is unsurprising, as the temperature is too high during the early stage of
LRNe for dust to form. The minor difference is that dust-deficit models have
slightly higher luminosity during the very late stage (see Figure 11), mainly
because the optical depth of the dust region is smaller.

Figure 11 also shows the time evolution of $\log_{10}(a_{\rm{rad}}/g)$ and
$\log_{10}\kappa_{\rm{R}}$ of the dust-deficit model. Comparing Figure 11 to
Figure 4, we notice that the difference is mainly in the upper right corner of
the figure, caused by the lowered dust opacity. The time and radial location
of the dust formation is not changed. It indicates that the dust formation may
have a limited impact on the radiation hydrodynamics of LRNe events during the
early stage. In the late stage, due to radiation force acceleration, dust
formation may play a role in shaping the ejecta, but still not a major
physical mechanism of unbinding the ejecta (Bermúdez-Bustamante et al., 2024).

##  Appendix H convergence

Convergence is very important for simulations that have iterative solvers and
mesh refinement. It is also well-known that radiation hydrodynamic equations
may be sensitive to the choice of timesteps. We show that our simulations
approach convergence by comparing the light curves of our original resolution
to a lower base resolution ($N=768$), which also doubles the timesteps. Figure
12 has the light curves of the high and low resolutions of the shock and
shock-free models. We can see that they are similar, indicating that our
results are close to convergence.

Figure 12: The black line with dots shows the light curve of AT2019zhd. The
red and blue colors indicate the shock (s) and shock-free (sf) models. The
solid and dashed lines indicate the original and low-resolution models,
respectively.


# Fabrication of high-quality PMMA/SiOx spaced planar microcavities for strong
coupling of light with monolayer WS2 excitons

Tinghe Yun  ARC Centre of Excellence in Future Low-Energy Electronics
Technologies and Department of Quantum Science and Technology, Research School
of Physics, The Australian National University, Canberra, ACT 2601, Australia
Songshan Lake Materials Laboratory, Dongguan 523808, Guangdong, China
Institute of Physics, Chinese Academy of Science, Beijing, 100190, China
Eliezer Estrecho  ARC Centre of Excellence in Future Low-Energy Electronics
Technologies and Department of Quantum Science and Technology, Research School
of Physics, The Australian National University, Canberra, ACT 2601, Australia
Andrew G. Truscott  Department of Quantum Science and Technology, Research
School of Physics, The Australian National University, Canberra, ACT 2601,
Australia  Elena A. Ostrovskaya<EMAIL_ADDRESS>ARC Centre of
Excellence in Future Low-Energy Electronics Technologies and Department of
Quantum Science and Technology, Research School of Physics, The Australian
National University, Canberra, ACT 2601, Australia  Matthias J. Wurdack
<EMAIL_ADDRESS>ARC Centre of Excellence in Future Low-Energy
Electronics Technologies and Department of Quantum Science and Technology,
Research School of Physics, The Australian National University, Canberra, ACT
2601, Australia

###### Abstract

Exciton polaritons in atomically-thin transition metal dichalcogenide crystals
(monolayer TMDCs) have emerged as a promising candidate to enable topological
transport, ultra-efficient laser technologies, and collective quantum
phenomena such as polariton condensation and superfluidity at room
temperature. However, integrating monolayer TMDCs into high-quality planar
microcavities to achieve the required strong coupling between the cavity
photons and the TMDC excitons (bound electron-hole pairs) has proven
challenging. Previous approaches to integration had to compromise between
various adverse effects on the strength of light-matter interactions in the
monolayer, the cavity photon lifetime, and the lateral size of the
microcavity. Here, we demonstrate a scalable approach to fabricating high-
quality planar microcavities with an integrated monolayer WS2 layer-by-layer
by using polymethyl methacrylate/silicon oxide (PMMA/SiOx) as a cavity spacer.
Because the exciton oscillator strength is well protected against the required
processing steps by the PMMA layer, the microcavities investigated in this
work, which have quality factors of above $10^{3}$, can operate in the strong
light-matter coupling regime at room temperature. This is an important step
towards fabricating wafer-scale and patterned microcavities for engineering
the exciton-polariton potential landscape, which is essential for enabling
many proposed technologies.

Planar microcavities with embedded exciton-hosting materials based on the
design of vertical-cavity surface-emitting lasers (VCSELs) have enabled the
regime of strong light-matter coupling between excitons and photons [1]. The
formation of exciton polaritons (polaritons herein), part-light part-matter
bosonic quasiparticles, in these semiconductor microcavities has triggered
significant research efforts, which led to the creation of Bose-Einstein
condensates (BEC) [2] and superfluids [3], the demonstration of non-Hermitian
effects [4], the development of electrically-driven ultra-efficient lasers [5]
and the realisation of topological edge transport [6] on microchips. By using
semiconductors with large exciton binding energies, some of these observations
were successfully realised at room temperature [7, 8, 9, 10, 11, 12, 13, 14].

A promising family of inorganic semiconductors for future polaritonic devices
are atomically-thin TMDCs [15, 16, 17, 18, 19]. However, integrating these
monolayers into planar microcavities to create polaritons has turned out to be
challenging because they are extremely fragile and easily damaged by various
fabrication techniques [20]. Many attempts at TMDC integration by using
metallic top mirrors [21, 19] or depositing the top structures directly on top
of the monolayer with aggressive material deposition techniques [16, 22], have
led to dramatic degradation of the photon lifetimes, exciton quantum yields,
or exciton-photon coupling strengths. These properties are significantly
improved by recent techniques for integrating TMDCs into all-dielectric
microcavities. However, these approaches are based on micrometer-sized
distributed Bragg reflector (DBR) chips [23, 24, 25, 26, 27] or mechanically
exfoliated hexagonal Boron Nitride (hBN) crystals [28, 29], which massively
limit its technological scalability. While all of these approaches have
enabled the demonstration of many striking effects of TMDC polaritons [16,
18], which include signatures of bosonic condensation [30, 26], the spin-
valley Hall effect [23], and ballistic propagation at room temperature [25],
fabrication recipes which allow for both scalability and high performance of
the final device to enable future TMDC polariton based technologies are still
to be developed.

In this work, we introduce a method for integrating monolayer WS2 into high-
quality planar microcavities with the lateral size mainly limited by the size
of the monolayer. In particular, we layer PMMA/SiOx on top of the monolayer
before depositing the top DBR via plasma enhanced chemical vapour deposition
(PECVD). This is in contrast to using metal as the top mirror of the cavity as
was done for the heavily investigated DBR/metal cavities [21, 31, 32, 19, 33,
34]. We study the effects of each processing step on the optical properties of
the WS2 monolayer showing that the exciton oscillator strength is well
preserved inside these structures, which enables strong light-matter coupling
at room temperature. Finally, we demonstrate that the quality factors
(Q-factors) of these microcavities can reach values on the order of $10^{3}$,
therefore exceeding the Q-factors of DBR/metal cavities by $1-2$ orders of
magnitude [21, 31, 32, 34] (see Supplementary Table S1). TMDC based
microcavities with comparable Q-factors were previously only realised on
microchips by using the technologies that inherently limit the lateral size of
the devices [28, 23, 24, 25] or dramatically strain the monolayer [35, 29].

Figure 1: Design and fabrication procedure for the microcavity. (a-d)
Schematics of the fabrication steps starting from (a) assembly of the base
structure consisting of a TiO2/SiO2 DBR substrate and a monolayer WS2,
followed by (b) the deposition of the PMMA and (c) SiOx layers, and finished
by (d) the deposition of the top SiNx/SiOx DBR via PECVD.

The fabrication procedure for the PMMA/SiOx spaced planar microcavities tested
and discussed in this work is schematically demonstrated in Fig. 1a-d. The
TiO2/SiO2 DBR substrate, which consists of 17.5 mirror pairs, is fabricated
via sputtering [36] and finishes with a SiO2 capping layer that has a
thickness corresponding to half of the cavity spacer. After transferring the
monolayer onto the substrate (see Fig. 1a), a PMMA layer is spin-coated on
top, which acts both as a protection layer and a cavity spacer (see Fig. 1b).
The remaining part of the cavity-spacer (SiOx) and the top DBR (SiNx/SiOx),
starting with the SiNx layer, are deposited via PECVD at 150 ∘C to finalise
the microcavity (see Fig. 1c-d and Supplementary Section S2 for details). The
SiOx spacer allows us to fine-tune the cavity energy and also prevents
cracking of the top DBR, which we found to occur when directly depositing it
onto PMMA. Most effective confinement of photons with the wavelength
$\lambda_{C}$ is achieved when the total cavity spacer and each DBR layer have
a thickness of $\lambda_{C}/2n$ and $\lambda_{C}/4n$, respectively, where $n$
is the respective refractive index of each layer [37]. In principle, the DBR
substrate can be fabricated with any deposition technique suitable for making
highly reflective DBRs with a smooth, high-quality layer finish.

Figure 2: Effects of PMMA on the optical properties of WS2. (a) Schematics of
the PMMA encapsulation process. (b-c) Comparison of (b) the PL spectra, (c)
the corresponding PL maps, and (d) the derivative of the reflectance contrast
spectra for the as-exfoliated monolayer, monolayer capped by PMMA, and the
PMMA-capped monolayer baked at 180 ∘C. The scale bar size in panel (c) is
$10~{}\mu\mathrm{m}$.

To understand the effect of PMMA deposition on the optical properties of
monolayer WS2 (see Fig. 1b), we mechanically exfoliated a monolayer (bulk
crystals sourced from HQ Graphene [38]) on commercial SiO2/Si substrates
(sourced from Nova materials [39]). A $80~{}\mathrm{nm}$ layer PMMA (molecular
weight $\sim$ 950,000 g/mol dissolved in anhydrous anisole with a 2.2%
concentration) was then spin-coated on top of the monolayer at 6000 rpm for 60
s (see Fig. 2a). To fully remove the anisole residues and ensure a high-
quality surface finish, the PMMA was hard baked at 180 ∘C for 90 s [21].

The tests of the optical response of the WS2 excitons were performed via real-
space photoluminescence (PL) and white-light reflectivity measurements under
ambient conditions. For the PL studies, we excited the samples with a large
Gaussian spot ($d\approx 25~{}\mu\mathrm{m}$) from a frequency doubled ND:YAG
continuous wave (cw) 532-nm laser source with the excitation energy
(E$\approx$ 2.33 eV [40]) above the A-valley exciton energy of monolayer WS2
[41]. The reflectivity spectra were measured under illumination of the sample
surface with a tungsten halogen white light source. The PL spectra were
obtained by averaging the signal from a small sample area, and the derivative
of the reflectance contrast spectra $d(\Delta R/R_{\rm ref})/dE$ was derived
by also measuring the white light reflection next to the monolayer as
reference $R_{\rm ref}$, with $\Delta R=R-R_{\rm ref}$. Since the product of
linewidth and amplitude of the reflectance contrast spectrum at the excitonic
resonance scales with the exciton oscillator strength, it quantifies the
light-matter interaction in the monolayer [42, 18] required for the strong
exciton-photon coupling regime.

The PL spectra in Fig. 2b show the PL quenching after the capping of monolayer
WS2 with PMMA, and further degradation after baking it. This change is also
seen in the corresponding PL maps (see Fig. 2c). The pronounced PL quenching
effect is possibly due to the chemical instability of this material. WS2
crystals are normally n-type semiconductors due to the presence of sulphur
vacancies [43, 44, 45, 46], and therefore, they can be strongly affected by
the surrounding molecules [47], such as organic and water molecules that can
act as dopants [48]. As a consequence, the PMMA capping process can modify the
doping level of the monolayer effectively causing a reduction of the exciton
quantum yield [47]. The additional decrease of the exciton PL after the baking
process is likely caused by the elevated temperature accelerating the aging
process (oxidation) of the monolayer, which further decreases the exciton
quantum yield [49]. However, the derivative of the reflectance contrast
spectra in Fig. 2d reveals that the exciton oscillator strength is barely
affected after the PMMA capping and baking, which indicates that PMMA
encapsulation will have a negligible effect on the exciton-photon
interactions.

Figure 3: Effects of PECVD deposition on the optical properties of WS2 with
and without PMMA passivation. (a-b) Comparison of (a) the PL and (b) the
derivative of the reflectance contrast spectra for monolayer WS2 capped with
PMMA, and further overgrown with SiOx of different thicknesses by PECVD. (c-d)
Comparison of (c) the PL spectra and (d) the derivative of the reflectance
contrast spectra of monolayer WS2 after deposition of (black) PMMA, (blue)
PMMA/SiOx and (red) SiOx.

Since the monolayers are fragile and easily damaged by many fabrication
techniques including PECVD [20], we further tested how well the PMMA layer
protects the monolayer against PECVD of SiOx. The influence of the deposition
on the exciton PL and absorption in the PMMA-capped monolayer WS2 is
demonstrated in Fig. 3a-b. Clearly, the exciton PL intensity significantly
decreases after the deposition process. The observed quenching is stronger
with a larger SiOx thickness, which scales with the duration of the
deposition. Between the layer thicknesses of 1.2 $\mu$m SiOx and 1.8 $\mu$m
SiOx the PL intensity remains constant and yields $\sim 25\%$ of the PL
intensity of the PMMA-capped WS2, which shows that PECVD induced PL quenching
saturates at this level. Similarly, the exciton absorption shown in Fig. 3b
also decreases after material deposition. This degradation of the PL and
oscillator strength could be attributed to the penetration of gaseous plasma
(leakage of conductive ions) through the microporous PMMA layer [50]. Once the
plasma reaches the monolayer surface, it can degrade the material by creating
additional defects, which explains the observed behaviour. Nevertheless, the
decrease of excitonic absorption after material deposition is much less
compared to that of the excitonic PL. Accordingly, the oscillator strength
remains relatively high after deposition of PMMA and SiOx in contrast with the
exciton quantum yield.

To further test the protective properties of the PMMA layer, we prepared an
uncapped monolayer sample, on which we directly deposited $40~{}\mathrm{nm}$
of SiOx via PECVD. This thickness was chosen to maintain some optical response
from the exciton in the WS2/SiOx structure [20]. The PL spectra in Fig. 3c
reveal that the PL emission of PMMA-capped monolayers with subsequent 1.2
$\mu$m SiOx deposition is almost twice as strong as that of the uncapped
monolayer after deposition of 40 nm SiOx. The reflectance contrast spectra
(see Fig. 3d) also show much stronger degradation of the exciton oscillator
strength for the uncapped monolayer compared to the capped one. Furthermore,
the uncapped monolayer shows a pronounced shoulder in the PL spectrum at
$E\approx 1.97~{}\mathrm{eV}$ likely stemming from trion emission [51, 52,
44], which indicates a significant defect-induced charge doping effect caused
by the deposition process [46, 20, 53]. These results clearly highlight that
the PMMA layer protects the monolayer against the dielectric deposition by
PECVD. Since the thickness of deposited SiOx is directly related to the
duration of exposure to plasma, the PMMA-capped monolayer was exposed to the
plasma for 30 times longer, but still maintained a larger excitonic response,
which confirms the protective function of the PMMA layer.

Figure 4: Design and characterisation of an empty microcavity. (a) The SEM
image of a SiOx/SiNx DBR fabricated by PECVD on top of PMMA as marked in
(inset) the schematics of the empty microcavity. The scale bar size is
$500~{}\mathrm{nm}$. (b) The (blue) refractive indices of the layers around
the centre of the microcavity and (red) the simulated E-field distribution of
the cavity mode. (c,d) Angle-resolved reflectivity spectra of microcavities
with the resonance energies (c) negatively and (d) positively detuned from the
(X) WS2 exciton energy at room temperature. (e) Reflectivity spectra at normal
incidence with (black) 8.5, (red) 10.5 and (blue) 14.5 SiNx/SiOx top DBR
mirror pairs. (f) Q-factor of the microcavity as a function of top DBR mirror
pairs.

In addition to testing how PMMA and dielectric material deposition affect the
optical response of the WS2 excitons at room temperature, we determined the
quality of the photon mode that can be achieved in our cavity without the
integration of a monolayer. To ensure that the PMMA/SiOx cavity provides a
suitable substrate for the top DBR, we perform scanning electron microscope
(SEM) imaging on the cross-section of a SiNx/SiOx DBR grown on top of a PMMA
layer via PECVD (see Fig. 4a). The SEM image reveals that the SiOx and SiNx
layers have a consistent thickness throughout the mirror and a smooth lateral
profile on the scale of the light wavelength that is crucial for high optical
performance.

Figure 4b presents the simulation of the optical E-field distribution inside
an empty cavity with a cavity-length of $\lambda_{C}/2n$ fulfilling the Bragg
condition for $\lambda_{C}=615~{}\mathrm{nm}$, performed with the transfer
matrix method [54]. As usual for $\lambda/2$-cavities [37], the strongly
confined resonant E-field has its maximum near the centre of the cavity, which
allows for maximum energy exchange between photons and excitons when the
monolayer is placed at this position. This is ensured by our sample design,
where the last layer of the DBR substrate corresponds to the first half of the
cavity spacer. By varying the thickness of the SiOx layer on top of the PMMA
layer, both having a similar refractive index, we can fine-tune the cavity
length $d$, and therefore, the cavity mode energy EC. This, in turn, enables
control of the exciton-photon detuning: $\Delta$ = EC-EX, where EX is the
exciton energy.

The angle-resolved reflectivity spectra for empty cavities with
$110~{}\mathrm{nm}$ and $95~{}\mathrm{nm}$ of PMMA/SiOx spacer are presented
in Fig. 4c-d, respectively, with the top DBRs consisting of 13.5 SiNx/SiOx
mirror pairs. The spectra clearly demonstrate the characteristic parabolic
dispersion of the cavity photons and show that the photon energy can be tuned
well across the exciton energy. The linewidth of the cavity photon ${\rm
FWHM}_{C}$ can be controlled with the number of mirror pairs of the top DBR
(see Fig. 4e), which allows us to increase the quality factor $Q=E_{C}/{\rm
FWHM}_{C}$ to $Q\approx 4\cdot 10^{3}$ reaching the resolution limit of our
optical setup (see Fig. 4f). These results demonstrate that this sample design
enables a good control over both the cavity photon energy and linewidth and
allows us to reach much higher Q-factors compared to those of the similarly
constructed DBR/metal cavities, e.g., [21, 19, 32] (see Supplementary Table
S1). However, to change the exciton-photon detuning, the sample needs to be
reproduced with a different thickness of the cavity spacer, and the achievable
Q-factors are smaller than the theoretically expected values (see Fig. 4f).
The growth of the Q-factor with the number of mirror pairs in the top DBR
decreases above 12 pairs, likely due to scattering and absorption processes in
the dielectric layers.

Figure 5: Characterisation of the microcavity with an integrated monolayer
WS2. (a) PL spectrum of a monolayer WS2 on top of the DBR substrate fitted
with (black) a three peak Voigt function accounting for (XL) localised, (XC)
charged and (X0) neutral excitons. (b) Reflectance contrast spectrum of the
monolayer WS2 and its derivative. (c) Reflectivity spectrum of the cavity
around the photon energy at normal incidence next to the position of the
monolayer fitted with a (red) Lorentzian function; (blue) theoretical
reflectivity spectrum calculated via transfer matrix calculations. (inset)
Calculated photon energies at $k_{||}=0$ with and without the monolayer (ML).
The cavity energy is shifted by $\sim 11.4~{}\mathrm{meV}$ due to the optical
thickness of the monolayer ($d_{\rm WS_{2}}\approx 0.7~{}\mathrm{nm}$ [41],
$\varepsilon_{B}^{\rm WS_{2}}\approx 15.831$[55]). (d) Angle-resolved PL
spectrum of the full microcavity. The dotted black lines correspond to the
experimentally determined exciton and photon energies (see panels a-c), the
red dotted lines are the resulting lower and upper polariton branches Eq. (1)
with $V=(15.3\pm 0.3)~{}\mathrm{meV}$ and the red circles are the fitted peak
energies (see Supplementary Section S3). (Inset) (black) measured and (blue)
calculated PL peak of the lower polariton at $k_{||}=0$.

Finally, we fabricated the whole microcavity with an integrated monolayer WS2,
as schematically shown in Fig. 1. Figure 5a presents the PL spectrum of a
monolayer WS2 placed on the DBR substrate, after deposition of the PMMA spacer
(see Fig. 1b). When fitting the spectrum with a three peak Voigt-function
accounting for the localised, charged and neutral exciton emission [44] (see
Fig. 5a), we can extract the energy and the linewidth of the neutral exciton
peak: $E_{X}\approx 2.013~{}\mathrm{eV}$ and ${\rm FWHM}_{X}\approx
35~{}\mathrm{meV}$, with a contribution of homogeneous broadening ${\rm
FWHM}_{X}^{H}\approx 12~{}\mathrm{meV}$ and inhomogeneous broadening ${\rm
FWHM}_{X}^{IH}\approx 28~{}\mathrm{meV}$. The large spectral weight of the
neutral exciton is an evidence of the high quality of the DBR substrate
causing negligible doping [53] or strain effects [56, 57] in the monolayer,
which can reduce exciton-photon interactions and prevent strong light-matter
coupling.

The reflectance contrast spectrum and its derivative with respect to the
energy (see Fig. 5b) reveal the exciton energy $E_{X}\approx
2.018~{}\mathrm{eV}$ and linewidth $FWHM_{X}\approx 35~{}\mathrm{meV}$, which
is in agreement with the Stokes shifted exciton peak of the PL spectrum, and
show a significant exciton oscillator strength in the monolayer that ensures a
robust exciton-photon interaction strength.

After deposition of the SiOx spacer and the top DBR consisting of 13.5
SiNx/SiOx layer pairs, we characterise the cavity mode away from the monolayer
area. Figure 5c shows the measured and theoretical reflectivity spectra of the
microcavity with the cavity mode energy $E_{C}(k_{||}=0)\approx
2.0032~{}\mathrm{eV}$ and linewidths ${\rm FWHM}_{C}^{\rm Exp}=(565\pm
13)~{}\mu\mathrm{eV}$ and ${\rm FWHM}_{C}^{\rm Th}=257~{}\mu\mathrm{eV}$,
respectively. While the corresponding measured Q-factor $Q^{\rm Exp}\approx
3500$ is smaller than the theoretical Q-factor $Q^{\rm Th}\approx 7800$, it is
similar to the Q-factor of the empty microcavity (see Fig. 4f) and in general,
comparable to the largest values reached for microcavities with embedded TMDCs
[58, 24, 28, 29, 25] (see Supplementary Table S1).

Figure 5d presents the angle-resolved PL measurement of the whole structure,
which was obtained via excitation of the monolayer with a tightly focused
laser spot ($\lambda=532~{}\mathrm{nm}$) and collecting the PL in momentum
space with a NA=0.65 objective, projected into a 2D spectrometer. The measured
dispersion, which decreases in intensity at larger energies due to
thermalisation [21, 25], can be approximated with the model of two coupled
oscillators [59]:

$\bar{{E}}_{U/L}=0.5\left(\bar{E}_{X}+\bar{E}_{C}\pm\sqrt{(\Delta-i\delta)^{2}+4V^{2}}\right),$
(1)

where $\bar{E}=E-i\gamma$ are the complex energies of the lossy modes, $V$ is
the strength of coupling between the cavity photons and the excitons, and
$\delta=\gamma_{C}-\gamma_{X}$ (see Fig. 5d). Note, that $\gamma$ corresponds
to ${\rm HWHM}$ of the measured modes [59].

Using the energies and linewidths of the exciton and cavity resonances (see
Fig. 5a-c), we find an excellent agreement between Eq. (1) and the measured
dispersion with $V=(15.3\pm 0.3)~{}\mathrm{meV}$ (see Fig. 5d). The
corresponding Hopfield coefficient of the lower polariton branch

$|X|^{2}=0.5\left[1+\frac{\Delta}{\sqrt{\mathfrak{Re}\left(\sqrt{4V^{2}-\delta^{2}}\right)}+\Delta^{2}}\right]$
(2)

allows us to approximate the polariton PL spectrum with a Voigt profile, where
${\rm FWHM}_{L}^{IH/H}=|X|^{2}{\rm FWHM}_{X}^{IH/H}+\left(1-|X|^{2}\right){\rm
FWHM}_{C}^{IH/H}$, which reproduces the measured PL well (see Fig. 5d, inset).
Since $V>|\delta|/2\approx 8.6~{}\mathrm{meV}$ and
$V>0.5\sqrt{\gamma_{X}^{2}+\gamma_{C}^{2}}\approx 8.8~{}\mathrm{meV}$, the
sample operates in the strong exciton-photon coupling regime, exhibiting
distinct exciton-polariton modes [59]. Hence, the sample fabrication
technology presented here is suitable for making high-Q planar microcavities
with an integrated atomically-thin WS2 crystal operating in the strong light-
matter coupling regime.

In summary, we have developed a scalable approach for making high-quality
planar microcavities with integrated atomically-thin semiconductors suitable
for hosting exciton-polaritons at room temperature. This is achieved by
depositing the two DBRs, a monolayer WS2 and a PMMA/SiOx spacer layer-by-
layer. We demonstrate that the PMMA layer has only a minor effect on the WS2
exciton oscillator strength and effectively protects the monolayer against
further deposition of the top DBR, which allows for strong light-matter
interactions in the structure. Since the PMMA/SiOx spacer and the top DBR
cover the whole microchip homogeneously, the functional area of the loaded
microcavity is mainly limited by the size of the monolayer, which can be
scaled up to cm-scale with recently developed synthesis [60, 61] and
exfoliation [62, 63] techniques. Further, the PMMA can be easily patterned by
lithography methods, paving the way towards creating polariton lattices,
waveguides and other 2D potential landscapes on a microchip, as previously
demonstrated for microcavities with embedded perovskites [14, 64].

## Acknowledgements

###### Acknowledgements.

This work was supported by the Australian Research Council (ARC) grants
CE170100039 and DE220100712. We acknowledge the Australian National
Fabrication Facility (ANFF) OptoFab at its node in the Australian Capital
Territory (ACT) for fabricating the DBR substrate and technical support for
sample fabrication from the ANFF ACT node.

## References

  * Weisbuch  _et al._ [1992] C. Weisbuch, M. Nishioka, A. Ishikawa, and Y. Arakawa, Observation of the coupled exciton-photon mode splitting in a semiconductor quantum microcavity, Phys. Rev. Lett. 69, 3314 (1992).
  * Kasprzak  _et al._ [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. Keeling, F. M. Marchetti, M. H. Szymánska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Bose-Einstein condensation of exciton polaritons, Nature 443, 409 (2006).
  * Amo  _et al._ [2009] A. Amo, J. Lefrère, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdré, E. Giacobino, and A. Bramati, Superfluidity of polaritons in semiconductor microcavities, Nature Physics 5, 805 (2009).
  * Gao  _et al._ [2015] T. Gao, E. Estrecho, K. Y. Bliokh, T. C. Liew, M. D. Fraser, S. Brodbeck, M. Kamp, C. Schneider, S. Höfling, Y. Yamamoto, F. Nori, Y. S. Kivshar, A. G. Truscott, R. G. Dall, and E. A. Ostrovskaya, Observation of non-Hermitian degeneracies in a chaotic exciton-polariton billiard, Nature 526, 554 (2015).
  * Schneider  _et al._ [2013] C. Schneider, A. Rahimi-Iman, N. Y. Kim, J. Fischer, I. G. Savenko, M. Amthor, M. Lermer, A. Wolf, L. Worschech, V. D. Kulakovskii, I. A. Shelykh, M. Kamp, S. Reitzenstein, A. Forchel, Y. Yamamoto, and S. Höfling, An electrically pumped polariton laser, Nature 497, 348 (2013).
  * Klembt  _et al._ [2018] S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, R. Ge, M. A. Bandres, M. Emmerling, L. Worschech, T. C. Liew, M. Segev, C. Schneider, and S. Höfling, Exciton-polariton topological insulator, Nature 562, 552 (2018).
  * Christopoulos  _et al._ [2007] S. Christopoulos, G. B. H. von Högersthal, A. J. D. Grundy, P. G. Lagoudakis, A. V. Kavokin, J. J. Baumberg, G. Christmann, R. Butté, E. Feltin, J.-F. Carlin, and N. Grandjean, Room-temperature polariton lasing in semiconductor microcavities, Phys. Rev. Lett. 98, 126405 (2007).
  * Guillet  _et al._ [2011] T. Guillet, M. Mexis, J. Levrat, G. Rossbach, C. Brimont, T. Bretagnon, B. Gil, R. Butté, N. Grandjean, L. Orosz, F. Réveret, J. Leymarie, J. Zúñiga-Pérez, M. Leroux, F. Semond, and S. Bouchoule, Polariton lasing in a hybrid bulk ZnO microcavity, Applied Physics Letters 99, 161104 (2011).
  * Su  _et al._ [2017] R. Su, C. Diederichs, J. Wang, T. C. H. Liew, J. Zhao, S. Liu, W. Xu, Z. Chen, and Q. Xiong, Room-temperature polariton lasing in all-inorganic perovskite nanoplatelets, Nano Letters 17, 3982 (2017).
  * Dietrich  _et al._ [2016] C. P. Dietrich, A. Steude, L. Tropf, M. Schubert, N. M. Kronenberg, K. Ostermann, S. Höfling, and M. C. Gather, An exciton-polariton laser based on biologically produced fluorescent protein, Science Advances 2, 8 (2016).
  * Dusel  _et al._ [2021] M. Dusel, S. Betzold, T. H. Harder, M. Emmerling, J. Beierlein, J. Ohmer, U. Fischer, R. Thomale, C. Schneider, S. Höfling, and S. Klembt, Room-temperature topological polariton laser in an organic lattice, Nano Letters 21, 6398 (2021).
  * Lerario  _et al._ [2017] G. Lerario, A. Fieramosca, F. Barachati, D. Ballarini, K. S. Daskalakis, L. Dominici, M. D. Giorgi, S. A. Maier, G. Gigli, S. Kéna-Cohen, and D. Sanvitto, Room-temperature superfluidity in a polariton condensate, Nature Physics 13, 837 (2017).
  * Su  _et al._ [2021] R. Su, E. Estrecho, D. Biegańska, Y. Huang, M. Wurdack, M. Pieczarka, A. G. Truscott, T. C. H. Liew, E. A. Ostrovskaya, and Q. Xiong, Direct measurement of a non-hermitian topological invariant in a hybrid light-matter system, Science Advances 7, 45 (2021).
  * Su  _et al._ [2020] R. Su, S. Ghosh, J. Wang, S. Liu, C. Diederichs, T. C. Liew, and Q. Xiong, Observation of exciton polariton condensation in a perovskite lattice at room temperature, Nature Physics 16, 301 (2020).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically thin ${\mathrm{MoS}}_{2}$: A new direct-gap semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * Liu  _et al._ [2014] X. Liu, T. Galfsky, Z. Sun, F. Xia, E. C. Lin, Y. H. Lee, S. Kéna-Cohen, and V. M. Menon, Strong light-matter coupling in two-dimensional atomic crystals, Nature Photonics 9, 30 (2014).
  * Sanvitto and Kéna-Cohen [2016] D. Sanvitto and S. Kéna-Cohen, The road towards polaritonic devices, Nature Materials 15, 1061 (2016).
  * Schneider  _et al._ [2018] C. Schneider, M. M. Glazov, T. Korn, S. Höfling, and B. Urbaszek, Two-dimensional semiconductors in the regime of strong light-matter coupling, Nature Communications 9, 2695 (2018).
  * Gu  _et al._ [2019] J. Gu, B. Chakraborty, M. Khatoniar, and V. M. Menon, A room-temperature polariton light-emitting diode based on monolayer WS2, Nature Nanotechnology 14, 1024 (2019).
  * Yun  _et al._ [2021] T. Yun, M. Wurdack, M. Pieczarka, S. Bhattacharyya, Q. Ou, C. Notthoff, C. K. Nguyen, T. Daeneke, P. Kluth, M. S. Fuhrer, A. G. Truscott, E. Estrecho, and E. A. Ostrovskaya, Influence of direct deposition of dielectric materials on the optical response of monolayer WS2, Applied Physics Letters 119, 133106 (2021).
  * Lundt  _et al._ [2016a] N. Lundt, S. Klembt, E. Cherotchenko, S. Betzold, O. Iff, A. V. Nalitov, M. Klaas, C. P. Dietrich, A. V. Kavokin, S. Höfling, and C. Schneider, Room-temperature Tamm-plasmon exciton-polaritons with a WSe2 monolayer, Nature Communications 7, 2695 (2016a).
  * Chen  _et al._ [2017] Y. J. Chen, J. D. Cain, T. K. Stanev, V. P. Dravid, and N. P. Stern, Valley-polarized exciton-polaritons in a monolayer semiconductor, Nature Photonics 11, 431 (2017).
  * Lundt  _et al._ [2019] N. Lundt, Łukasz Dusanowski, E. Sedov, P. Stepanov, M. M. Glazov, S. Klembt, M. Klaas, J. Beierlein, Y. Qin, S. Tongay, M. Richard, A. V. Kavokin, S. Höfling, and C. Schneider, Optical valley Hall effect for highly valley-coherent exciton-polaritons in an atomically thin semiconductor, Nature Nanotechnology 14, 770 (2019).
  * Rupprecht  _et al._ [2021] C. Rupprecht, N. Lundt, M. Wurdack, P. Stepanov, E. Estrecho, M. Richard, E. A. Ostrovskaya, S. Höfling, and C. Schneider, Micro-mechanical assembly and characterization of high-quality Fabry–Pérot microcavities for the integration of two-dimensional materials, Applied Physics Letters 118, 103103 (2021).
  * Wurdack  _et al._ [2021] M. Wurdack, E. Estrecho, S. Todd, T. Yun, M. Pieczarka, S. K. Earl, J. A. Davis, C. Schneider, A. G. Truscott, and E. A. Ostrovskaya, Motional narrowing, ballistic transport, and trapping of room-temperature exciton polaritons in an atomically-thin semiconductor, Nature Communications 12, 5366 (2021).
  * Anton-Solanas  _et al._ [2021] C. Anton-Solanas, M. Waldherr, M. Klaas, H. Suchomel, T. H. Harder, H. Cai, E. Sedov, S. Klembt, A. V. Kavokin, S. Tongay, K. Watanabe, T. Taniguchi, S. Höfling, and C. Schneider, Bosonic condensation of exciton–polaritons in an atomically thin crystal, Nature Materials 20, 1233 (2021).
  * Zhao  _et al._ [2022] J. Zhao, A. Fieramosca, R. Bao, W. Du, K. Dini, R. Su, J. Feng, Y. Luo, D. Sanvitto, T. C. Liew, and Q. Xiong, Nonlinear polariton parametric emission in an atomically thin semiconductor based microcavity, Nature Nanotechnology 17, 396 (2022).
  * Knopf  _et al._ [2019] H. Knopf, N. Lundt, T. Bucher, S. Höfling, S. Tongay, T. Taniguchi, K. Watanabe, I. Staude, U. Schulz, C. Schneider, and F. Eilenberger, Integration of atomically thin layers of transition metal dichalcogenides into high-Q, monolithic Bragg-cavities: an experimental platform for the enhancement of the optical interaction in 2D-materials, Optical Materials Express 9, 598 (2019).
  * Shan  _et al._ [2021] H. Shan, L. Lackner, B. Han, E. Sedov, C. Rupprecht, H. Knopf, F. Eilenberger, J. Beierlein, N. Kunte, M. Esmann, K. Yumigeta, K. Watanabe, T. Taniguchi, S. Klembt, S. Höfling, A. V. Kavokin, S. Tongay, C. Schneider, and C. Antón-Solanas, Spatial coherence of room-temperature monolayer WSe2 exciton-polaritons in a trap, Nature Communications 12, 6406 (2021).
  * Waldherr  _et al._ [2018] M. Waldherr, N. Lundt, M. Klaas, S. Betzold, M. Wurdack, V. Baumann, E. Estrecho, A. Nalitov, E. Cherotchenko, H. Cai, E. A. Ostrovskaya, A. V. Kavokin, S. Tongay, S. Klembt, S. Höfling, and C. Schneider, Observation of bosonic condensation in a hybrid monolayer MoSe2-GaAs microcavity, Nature Communications 9, 3286 (2018).
  * Lundt  _et al._ [2017] N. Lundt, P. Nagler, A. Nalitov, S. Klembt, M. Wurdack, S. Stoll, T. H. Harder, S. Betzold, V. Baumann, A. V. Kavokin, C. Schüller, T. Korn, S. Höfling, and C. Schneider, Valley polarized relaxation and upconversion luminescence from Tamm-plasmon trion-polaritons with a MoSe2 monolayer, 2D Materials 4, 025096 (2017).
  * Wurdack  _et al._ [2017] M. Wurdack, N. Lundt, M. Klaas, V. Baumann, A. V. Kavokin, S. Höfling, and C. Schneider, Observation of hybrid Tamm-plasmon exciton-polaritons with GaAs quantum wells and a MoSe2 monolayer, Nature Communications 8, 259 (2017).
  * Gu  _et al._ [2021] J. Gu, V. Walther, L. Waldecker, D. Rhodes, A. Raja, J. C. Hone, T. F. Heinz, S. Kéna-Cohen, T. Pohl, and V. M. Menon, Enhanced nonlinear interaction of polaritons via excitonic rydberg states in monolayer WSe2, Nature Communications 12, 2269 (2021).
  * Datta  _et al._ [2021] B. Datta, M. Khatoniar, P. Deshmukh, R. Bushati, S. De Liberato, S. K. Cohen, and V. M. Menon, Highly non-linear interlayer exciton-polaritons in bilayer MoS2, arXiv:2110.13326 (2021).
  * Zhao  _et al._ [2021] J. Zhao, R. Su, A. Fieramosca, W. Zhao, W. Du, X. Liu, C. Diederichs, D. Sanvitto, T. C. H. Liew, and Q. Xiong, Ultralow threshold polariton condensate in a monolayer semiconductor microcavity at room temperature, Nano Letters 21, 3331 (2021).
  * [36] The fabrication of the high-reflectivity substrate was peformed at the ACT OptoFab node of the Australian National Fabrication Facility - a company established under the National Collaborative Research Infrastructure Strategy to provide nano and microfabrication facilities for Australia’s researchers .
  * Kavokin  _et al._ [2017] A. Kavokin, J. Baumberg, G. Malpuech, and F. Laussy,  _Microcavities_ , Series on Semiconductor Science and Technology (OUP Oxford, 2017).
  * [38] http://www.hqgraphene.com.
  * [39] http://www.novawafers.com.
  * Zeng  _et al._ [2013] H. Zeng, G. B. Liu, J. Dai, Y. Yan, B. Zhu, R. He, L. Xie, S. Xu, X. Chen, W. Yao, and X. Cui, Optical signature of symmetry variations and spin-valley coupling in atomically thin tungsten dichalcogenides, Scientific Reports 3, 1608 (2013).
  * Mueller and Malic [2018] T. Mueller and E. Malic, Exciton physics and device application of two-dimensional transition metal dichalcogenide semiconductors, npj 2D Materials and Applications 2, 29 (2018).
  * Lundt  _et al._ [2016b] N. Lundt, A. Maryński, E. Cherotchenko, A. Pant, X. Fan, S. Tongay, G. Sek, A. V. Kavokin, S. Höfling, and C. Schneider, Monolayered MoSe2: a candidate for room temperature polaritonics, 2D Materials 4, 015006 (2016b).
  * Tongay  _et al._ [2013a] S. Tongay, J. Suh, C. Ataca, W. Fan, A. Luce, J. S. Kang, J. Liu, C. Ko, R. Raghunathanan, J. Zhou, F. Ogletree, J. Li, J. C. Grossman, and J. Wu, Defects activated photoluminescence in two-dimensional semiconductors: Interplay between bound, charged, and free excitons, Scientific Reports 3, 2657 (2013a).
  * Shang  _et al._ [2015] J. Shang, X. Shen, C. Cong, N. Peimyoo, B. Cao, M. Eginligil, and T. Yu, Observation of excitonic fine structure in a 2D transition-metal dichalcogenide semiconductor, ACS Nano 9, 647 (2015).
  * Carozo  _et al._ [2017] V. Carozo, Y. Wang, K. Fujisawa, B. R. Carvalho, A. McCreary, S. Feng, Z. Lin, C. Zhou, N. Perea-López, A. L. Elías, B. Kabius, V. H. Crespi, and M. Terrones, Optical identification of sulfur vacancies: Bound excitons at the edges of monolayer tungsten disulfide, Science Advances 3, 4 (2017).
  * Sebait  _et al._ [2021] R. Sebait, C. Biswas, B. Song, C. Seo, and Y. H. Lee, Identifying defect-induced trion in monolayer WS2 via carrier screening engineering, ACS Nano 15, 2849 (2021).
  * Tongay  _et al._ [2013b] S. Tongay, J. Zhou, C. Ataca, J. Liu, J. S. Kang, T. S. Matthews, L. You, J. Li, J. C. Grossman, and J. Wu, Broad-range modulation of light emission in two-dimensional semiconductors by molecular physisorption gating, Nano Letters 13, 2831 (2013b).
  * Zhang  _et al._ [2019] F. Zhang, Y. Lu, D. S. Schulman, T. Zhang, K. Fujisawa, Z. Lin, Y. Lei, A. L. Elias, S. Das, S. B. Sinnott, and M. Terrones, Carbon doping of WS2 monolayers: Bandgap reduction and p-type doping transport, Science Advances 5, 5 (2019).
  * Gao  _et al._ [2016] J. Gao, B. Li, J. Tan, P. Chow, T.-M. Lu, and N. Koratkar, Aging of transition metal dichalcogenide monolayers, ACS Nano 10, 2628 (2016).
  * Charlot  _et al._ [2011] B. Charlot, S. Gauthier, A. Garraud, P. Combette, and A. Giani, Pvdf/pmma blend pyroelectric thin films, Journal of Materials Science: Materials in Electronics 22, 1766 (2011).
  * Mak  _et al._ [2013] K. F. Mak, K. He, C. Lee, G. H. Lee, J. Hone, T. F. Heinz, and J. Shan, Tightly bound trions in monolayer MoS2, Nature Materials 12, 207 (2013).
  * Zhu  _et al._ [2015] B. Zhu, X. Chen, and X. Cui, Exciton binding energy of monolayer WS2, Scientific Reports 5, 9218 (2015).
  * Martín-Sánchez  _et al._ [2018] J. Martín-Sánchez, A. Mariscal, M. D. Luca, A. T. Martín-Luengo, G. Gramse, A. Halilovic, R. Serna, A. Bonanni, I. Zardo, R. Trotta, and A. Rastelli, Effects of dielectric stoichiometry on the photoluminescence properties of encapsulated WSe2 monolayers, Nano Research 11, 1399 (2018).
  * Born  _et al._ [2000] M. Born, E. Wolf, and A. Bhatia,  _Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light_ (Cambridge University Press, 2000).
  * Li  _et al._ [2014] Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E.-M. Shih, J. Hone, and T. F. Heinz, Measurement of the optical dielectric function of monolayer transition-metal dichalcogenides: ${\mathrm{MoS}}_{2}$, $\mathrm{Mo}\mathrm{S}{\mathrm{e}}_{2}$, ${\mathrm{WS}}_{2}$, and $\mathrm{WS}{\mathrm{e}}_{2}$, Phys. Rev. B 90, 205422 (2014).
  * Khatibi  _et al._ [2018] Z. Khatibi, M. Feierabend, M. Selig, S. Brem, C. Linderälv, P. Erhart, and E. Malic, Impact of strain on the excitonic linewidth in transition metal dichalcogenides, 2D Materials 6, 015015 (2018).
  * He  _et al._ [2013] K. He, C. Poole, K. F. Mak, and J. Shan, Experimental demonstration of continuous electronic structure tuning via strain in atomically thin MoS2, Nano Letters 13, 2931 (2013).
  * Sidler  _et al._ [2017] M. Sidler, P. Back, O. Cotlet, A. Srivastava, T. Fink, M. Kroner, E. Demler, and A. Imamoglu, Fermi polaron-polaritons in charge-tunable atomically thin semiconductors, Nature Physics 13, 255 (2017).
  * Savona  _et al._ [1995] V. Savona, L. Andreani, P. Schwendimann, and A. Quattropani, Quantum well excitons in semiconductor microcavities: Unified treatment of weak and strong coupling regimes, Solid State Communications 93, 733 (1995).
  * Lee  _et al._ [2017] J. Lee, S. Pak, P. Giraud, Y.-W. Lee, Y. Cho, J. Hong, A.-R. Jang, H.-S. Chung, W.-K. Hong, H. Y. Jeong, H. S. Shin, L. G. Occhipinti, S. M. Morris, S. Cha, J. I. Sohn, and J. M. Kim, Thermodynamically stable synthesis of large-scale and highly crystalline transition metal dichalcogenide monolayers and their unipolar n–n heterojunction devices, Advanced Materials 29, 1702206 (2017).
  * Choi  _et al._ [2022] S. H. Choi, S. J. Yun, Y. S. Won, C. S. Oh, S. M. Kim, K. K. Kim, and Y. H. Lee, Large-scale synthesis of graphene and other 2D materials towards industrialization, Nature Communications 13, 1484 (2022).
  * Desai  _et al._ [2016] S. B. Desai, S. R. Madhvapathy, M. Amani, D. Kiriya, M. Hettick, M. Tosun, Y. Zhou, M. Dubey, J. W. Ager III, D. Chrzan, and A. Javey, Gold-mediated exfoliation of ultralarge optoelectronically-perfect monolayers, Advanced Materials 28, 4053 (2016).
  * Liu  _et al._ [2020] F. Liu, W. Wu, Y. Bai, S. H. Chae, Q. Li, J. Wang, J. Hone, and X.-Y. Zhu, Disassembling 2d van der waals crystals into macroscopic monolayers and reassembling into artificial lattices, Science 367, 903 (2020).
  * Su  _et al._ [2018] R. Su, J. Wang, J. Zhao, J. Xing, W. Zhao, C. Diederichs, T. C. H. Liew, and Q. Xiong, Room temperature long-range coherent exciton polariton condensate flow in lead halide perovskites, Science Advances 4, 10 (2018).