DESY-23-020
IFT-UAM/CSIC-23-14
RESCEU-2/23
TUM-HEP-1450/23

Gravitational waves from defect–driven

phase transitions: domain walls
Simone Blasia, Ryusuke Jinnob,c, Thomas Konstandind,
Henrique Rubirae, Isak Stombergd

a Theoretische Natuurkunde and IIHE/ELEM, Vrije Universiteit Brussel, & The
International Solvay Institutes, Pleinlaan 2, B-1050 Brussels, Belgium
b Instituto de Física Teórica IFT-UAM/CSIC
C/ Nicolás Cabrera 13-15, Campus de Cantoblanco, 28049, Madrid, Spain c
Research Center for the Early Universe (RESCEU), Graduate School of Science
The University of Tokyo, Tokyo 113-0033, Japan d Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany e Physik Department
T31, Technische Universität München
James-Franck-Straße 1, D-85748 Garching, Germany

We discuss the gravitational wave spectrum produced by first-order phase
transitions seeded by domain wall networks. This setup is important for many
two-step phase transitions as seen for example in the singlet extension of the
standard model. Whenever the correlation length of the domain wall network is
larger than the typical bubble size, this setup leads to a gravitational wave
signal that is shifted to lower frequencies and with an enhanced amplitude
compared to homogeneous phase transitions without domain walls. We discuss our
results in light of the recent PTA hints for gravitational waves.

##  1 Introduction

Cosmological first-order phase transitions (FOPTs) in the early universe can
source a stochastic gravitational wave background (SGWB) which current and
near-future observational searches could detect. The remarkable successes of
LIGO-Virgo [1, 2, 3] to measure GWs from compact binaries in the kHz frequency
band marks the beginning of an era of gravitational wave (GW) probes of the
universe. While LIGO-Virgo and other interferometer-based GW searches continue
to be developed around the world, the pulsar timing arrays and the associated
experiments, including NANOGrav, EPTA, PPTA, and IPTA, are simultaneously
making strides towards detecting a stochastic background and have recently
reported evidence for a stochastic common-spectrum process among the pulsars
[4, 5, 6, 7]. While further evidence is necessary to conclusively confirm or
reject the SGWB hypothesis, it has nevertheless been shown that current
observations could possibly be explained by GWs from a FOPT [8, 9].

Next in line we find LISA, the Laser Interferometer Space Antenna, designed to
measure GWs in the mHz frequency band [10, 11, 12]. While the strongest
candidate to source GWs in this band are supermassive black-hole binaries,
LISA offers to probe GWs from FOPTs occurring around the electroweak scale
(TeV). The possibility to measure the GW spectrum over various frequency bands
calls for accurate predictions of GW spectra that could discriminate various
cosmological scenarios.

In this paper, we investigate the influence of defects on the dynamics of
FOPTs and their GW production (see [13, 14, 15, 16, 17] for early works on
impurities, [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] for the case of black
holes and local over–densities, and [29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41] for strings, monopoles, and defects in general). One interesting
possibility is that a FOPT responsible for the stochastic GW production is
locally catalyzed, or _seeded_ , by a domain wall (DW) network [42]. Depending
on the underlying particle physics model, the nucleation probability can be
significantly larger on the surfaces of the DWs than in the bulk of the
Universe.

The presence of a DW network effectively introduces a new length scale,
$\xi_{\rm DW}$, denoting the mean separation of the DWs. Whenever $\xi_{\rm
DW}$ exceeds the mean separation of adjacent bubbles, $\langle r\rangle$, (see
Appendix A for details)

$\xi_{\rm DW}\gg 2\langle r\rangle\simeq 6v_{w}/\beta\,,$ (1)

we expect that the network imprints a signature in the GW spectrum. In
particular, we expect that the peak frequency is shifted to lower frequencies
and that the peak is enhanced, akin to GW spectra computed from FOPTs modified
by macroscopic thermal fluctuations as discussed in Ref. [43]. We furthermore
expect that the spectral shape is modified due to the spatially inhomogeneous
bubble distribution as induced by the DW network; rather than colliding
bubbles, the DW network results in the propagation of sheets moving away from
the walls which collide after the mean time $\xi_{\rm DW}/v_{w}$.

Notice that in the extreme limit (1), the bubble nucleations can be regarded
as essentially simultaneous everywhere on the DW network compared to the
timescale for bubbles to travel between the walls, and $\beta$ becomes an
irrelevant parameter. Moreover, the dynamics of the DWs can be neglected in
this limit. Since the resulting GW spectrum will still be well-defined in this
limit, all relevant features in the nucleation history and the GW spectrum
cannot depend on $\beta$ for large DW networks. This point will be discussed
in detail in Section 2. By virtue of the arguments above, we only need a
snapshot of the DW network with an appropriate correlation length to infer a
possible _bubble nucleation history_ , i.e. the bubble nucleation sites and
times, and the corresponding GW spectrum. Instead of proper simulations of DW
networks, we will mimic them with a Monte Carlo implementation of the Ising
model. Using Monte Carlo Methods to study DWs was already put forward in the
classic review on DWs [44] and more generally in the analysis of percolation
[45]. The above procedure should at least qualitatively yield a reasonable
distribution of bubble nucleation sites. This topic will be discussed further
in Section 2.

Once a bubble nucleation history is obtained, we use hydrodynamic simulations
to compute the GW signal following [46]. The naive expectation is, again, that
once the DW network becomes relevant in the limit (1), the GW spectrum peak
frequency is reduced and the amplitude enhanced compared to a conventional
FOPT without a DW network. The last point fits well with the observation that
the GW production is proportional to the correlation length of the energy-
momentum tensor. In this study, we aim to explore these dependencies and how
the overall shape of GW spectra is affected by a DW network. We also discuss
our results in light of the recent PTA hints for gravitational waves.

The outline of the manuscript is as follows: In Section 2 we discuss the
relevant properties of the DW networks. We also provide some results for the
Ising model (in the low-temperature regime) and how to make contact with DW
networks using the correlation length of the system. In Section 3 we discuss
the gravitational wave spectra produced by our hydrodynamic simulations and
discuss the results in Section 4.

##  2 Domain walls

In this section, we start by describing the DW formation mechanism and discuss
the regimes in which DWs can have relevant phenomenological consequences.
Later, we argue why the Ising model resembles a realistic setup to mimic
configurations of DW networks from which bubble nucleation histories of a FOPT
seeded by such a network can be deduced.

###  2.1 Formation and dynamics of domain walls

DWs arise whenever a discrete symmetry is broken during the course of the
Universe [47]. The dynamics of such DWs has been studied numerically since
several decades (see for example [48, 49, 50, 51, 52]). The consensus is that
sometime after the FOPT, the DW network enters a scaling regime where the
correlation length of the system approaches the Hubble parameter and the
energy density roughly scales as $\rho\sim H\sim 1/t$ 111There is an argument
in the literature whether there are logarithmic corrections to this scaling,
see [53, 49].. In the long term, this is disastrous, since the DW network
would begin to dominate the energy density. Moreover, the DWs would induce
fluctuations in the plasma. Both these effects are in conflict with CMB
observations. In order to alleviate these issues, the DWs have to disappear
before they reach a sizable fraction of the energy density of the Universe.
Several solutions are conceivable:

  * •

i) The discrete symmetry is explicitly slightly broken [54]. This leads to
different vacua that are not entirely degenerate and the global minimum will
dominate after a while.

  * •

ii) There is a bias in the initial conditions [55]. Also in this scenario,
depending on the initial conditions, the DWs will decay before they become
harmful.

  * •

iii) The exact discrete symmetry is initially spontaneously broken. Later, the
symmetry is again restored whereby the DWs cease to exist.

In all these scenarios, the DWs will produce gravitational waves [56, 50, 51,
52, 57]. The setup we have in mind is of type iii), more concretely we
envision a two-step phase transition where the second step is of first order.
Consider the potential

$V=\frac{\lambda}{4}(\phi^{2}-v^{2})^{2}+\frac{\lambda_{m}}{4}s^{2}\phi^{2}+\frac{\lambda_{m}}{4}s^{4}-\frac{1}{2}\mu_{S}^{2}s^{2}\,,$
(2)

where $\phi$ is a Higgs field and $s$ a gauge singlet endowed with a $Z_{2}$
symmetry. For some regions in parameter space, this system will have a two-
step phase transition: first the system transitions to a phase where $s$
obtains a VEV. This breaks the $Z_{2}$ symmetry spontaneously and produces
DWs. At somewhat lower temperatures, the system attains the ground state of
the system with a finite $\phi$-VEV but unbroken $Z_{2}$ symmetry
($\left<s\right>=0$). The DWs cease to exist once the system completely
transitions into the low-temperature phase. This avoids the over-closure issue
and excessive fluctuations in the plasma [58, 59].

In this setup, the DWs exist only for a relatively short period of time.
Ultimately, the gravitational waves are not predominantly produced by the DWs
or the scalar field that drives the phase transition, but by sound waves that
are generated during the phase transition. To quantify the spectrum of GWs, we
use hydrodynamic simulations. Compared to conventional FOPTs, the nucleation
history of bubbles is modified since bubbles preferentially nucleate on the
DWs [42]. This turns out to be a generic feature in all the two-step parameter
space, where the temperature for the seeded FOPT is always higher than the
would-be homogeneous nucleation temperature [42]. Intuitively, a faster
nucleation rate around the DWs can be understood by noticing that the energy
stored in the DW tension can be used to tunnel through the potential barrier.

A crucial parameter controlling the phenomenology of the FOPT is the number of
DWs per Hubble patch at the time of seeded nucleation, or equivalently the
correlation length of the system which describes the average curvature radius
as well as the mean distance among the walls. In the picture we are
considering, DWs are formed during the first step of the FOPT when the $Z_{2}$
symmetry is spontaneously broken. If this transition is of the second-order
type, the typical size of domains at formation will be given by the
microscopic scale associated with the singlet mass $\xi_{\text{DW}}\sim
1/m_{s}$.

The subsequent evolution depends strongly on whether the DW motion is damped
by friction or if instead the walls can freely oscillate at (mildly)
relativistic velocities. In the latter case the DW network is expected to
increase its correlation length to saturate the causality bound from the
Hubble expansion, yielding $\xi_{\text{DW}}\sim t$ [48]. This corresponds to
the standard scaling regime where one expects $\mathcal{O}(1)$ DWs per Hubble
volume. On the other hand, when friction dominates the correlation length of
the system, it is supposed to increase at a slower pace. Standard arguments
based on a balance between the tension force,
$\sigma_{\text{DW}}/\xi_{\text{DW}}$, and the friction force that scales with
$T^{4}v$ (where $v$ is the DW velocity), yielding to
$\xi_{\text{DW}}\sim(G\sigma_{\text{DW}})^{\frac{1}{2}}t^{\frac{3}{2}}$, where
$G$ is the Newton’s constant [60]. Following this estimate, the correlation
length of the system will still be $\xi_{\text{DW}}\ll H^{-1}$ at the time of
seeded nucleation, as the two steps of the phase transition are usually
expected to be separated at most by a factor of a few in temperature [42]. In
this case, the phenomenology of the FOPT will resemble the standard
homogeneous case in which the duration is controlled by the time dependence of
the nucleation rate.

In our analysis we will not engage in a model–dependent determination of the
number of DWs at the seeded nucleation temperature, but rather use the
correlation length of the network as an input to determine the
phenomenological implications for the GW spectrum. In doing so we will
consider an intermediate case in which the DW size has grown to macroscopic
scales but the standard scaling regime has not yet been reached, so that
$\xi_{\text{DW}}<H^{-1}$. As the correlation length of the network always
stays smaller than the Hubble size, the expansion of the Universe can be
neglected in the hydrodynamic simulations of the GW production. Ultimately,
the expansion of the Universe might only lead to a moderate suppression of the
GW signal, but neglecting the expansion simplifies the numerical
implementation of the simulation.

From simulations of DWs, it is known that the main properties of DW networks
are encoded in a single quantity: the correlation length. Instead of
simulating a full DW network, it is tempting to mimic it with a simpler system
with similar properties. The prototypical example of this is the Ising model.
Close to the critical temperature, the system displays universality and is
applicable to many systems with phase transitions. Here we rely on the fact
that the outcome of the Ising model resembles very well how DWs appear in
genuine simulations of DW networks. There are several arguments why a snapshot
of the Ising model should suffice for our purpose. Whether a single snapshot
of the DW network provides enough information to simulate the bubble
nucleation history, or if alternatively the dynamics of the DW network is
relevant. Simulations in the scaling regime show that DWs can oscillate with
relativistic velocity. From the first nucleation of bubbles until the end of
the FOPT, the bubble wall of the new phase has to travel the correlation
length of the DW network, $\Delta t\sim\xi_{\rm DW}/v_{w}$. However, the time
window during which actual bubble nucleations happen is rather small, $\Delta
t\sim 10/\beta$. Hence, as long as one considers the limit (1), the motion of
the DWs is irrelevant.

In principle, there could be an issue with DWs propagating faster than the
velocity of the bubble wall interfaces, $v_{w}$, changing the bubble
nucleation history in the process. However, in the broken (new) phase, the DWs
cease to exist, since the $Z_{2}$ symmetry is restored. Overall, a snapshot of
the DW network should suffice to deduce a realistic bubble nucleation history.

###  2.2 Modeling the domain wall network

Even if the typical correlation length of the DW network coincides with the
one from the Ising model, one might wonder to what extent the substructure of
the two systems agree. For example, small islands could be present in the DW
network whose size distribution is not captured by the scale of the large
walls. These are, for instance, generated when two big walls collide and a
closed wall is produced in the aftermath. These walls would further collapse,
producing small structures (‘spheres’ or ‘bags’), similar to loops in networks
of cosmological strings.

A large number of closed islands would have a strong impact on the
gravitational wave production, as in this case the duration of the FOPT would
be controlled by the average distance between the islands rather than the
distance between the ”infinite” walls. However, those substructures are not
observed in sizable numbers in simulations (see [48] and [49] for visual
representations). In fact, differently to the case of cosmic strings, the
energy losses allowing the DW network to reach the scaling regime are most
likely due to particle radiation (see also Ref. [61] where it was found that
the ‘chopping’ parameter for a DW network is indeed compatible with zero).

In order to generate the actual bubble nucleation history from the Ising
model, one has to map this Ising model data to a realistic DW network. As
mentioned before, the details of this procedure, as we will see, are not
essential, since in the limit (1) (when we expect the GW signal to be
affected), the bubble nucleation happens simultaneously everywhere on the
network. This is because the time scale of the FOPT is much shorter than the
separation of the DWs in this limit, in which case the broken phase will grow
uniformly out of the DW network with a velocity given by the usual propagation
speed of the bubble walls $v_{w}$. In practice, the correlation length of the
DWs will only be one order larger than the mean bubble size and we will
explore two different possibilities for this mapping (see Appendix B.1).

Figure 1: Some representative DW configurations for cells with 0 to 4 corners
in the opposite phases. The corresponding surface weights of the DW network
are $\\{0,\sqrt{3}/8,1/\sqrt{2},1/\sqrt{2}+\sqrt{11}/8,1\\}$. The
representative of the empty cell is omitted.

Ideally, one would like to assign to every grid cell in the simulation a
weight that represents the surface area of the DW network in the cell. We
implemented two different methods to assign weights to the grid densities,
both leading to identical results. We now introduce one of those methods. The
presentation and comparison with the second (and to some extent simpler)
method is delivered in Appendix B.2.

One, and perhaps the most straightforward, way to assign weights to cells is
by observing for every cell the state of the eight corners and then
constructing a surface that separates the two phases. Such a surface runs
through the midpoints of the edges between corners of opposite states (see
Fig. 1).

While, in principle, every corner configuration is possible, only those which
are maximally connected are likely to occur, the others being suppressed by
virtue of the Ising model. A slight simplification to the method of
constructing surfaces through midpoints is, therefore, to assume that corners
of equal states are most compactly connected. We hence chose a configuration
that is representative of all cells of the same corner counts. Identification
of the weight can then proceed by means of simply counting the number of
corners in each phase since in this case, the mapping between corner counts
and weight is one-to-one. Consider $k,l$ corners in each phase. Then, the
corresponding weights $w_{kl}$ are $w_{08}=0$, $w_{17}=\sqrt{3}/8$,
$w_{26}=1/\sqrt{2}$, $w_{35}=1/\sqrt{2}+\sqrt{11}/8$, and $w_{44}=1$, which
are proportional to the enclosed DW area. While this simplification induces a
systematic error in the estimation of the DW area, we show in Appendix B.2
that the final results are insensitive to this error. The probability to
nucleate a bubble in cell $ijk$ at time $t$ is then proportional to
$w_{ijk}\exp(t\beta)$, where $w_{ijk}$ is the weight of cell $ijk$, whereby a
bubble nucleation history can be constructed as in the standard picture
without DWs (see Sec 2.1 of [46]). Such bubble nucleation histories will serve
as inputs to the hydrodynamic simulations, which will be discussed in the next
section. Notice that, since $w_{ijk}\neq 0$ only for cells containing a DW
surface, bubbles nucleate along the DWs only.

As emphasized before, the DW network is characterized by the mean separation
between the DWs, $\xi_{\rm DW}$. One way to estimate this measure is to
consider the total surface of the DW network, $S_{\rm DW}$. The mean
separation is then approximately

$\xi_{\rm DW}\sim\frac{V}{S_{\rm DW}}$ (3)

where $V=L^{3}$ denotes the volume of the simulation. For example, for the DW
network in Fig. 2, we obtain $\xi_{\rm DW}\sim 0.1\times L$, where $L$ denotes
the box size of the simulation. Starting from a random initial state and using
the Metropolis algorithm, we perform $4\times 10^{8}$ phase alterations on a
grid of size $128^{3}$ with the parameter $J/T=2$ to yield this correlation
length (see App. B).

In this estimate, the total surface area of the DW network was extracted using
the weights discussed above. Alternatively, one can extract the correlation
length using the correlation function

$\xi(\vec{x})=\left<\phi(\vec{y})\phi(\vec{y}+\vec{x})\right>\,.$ (4)

The field $\phi$ denotes the order parameter of the DW network that attains
$\pm v$ in the two degenerate minima. Moreover, one can study correlations in
Fourier space leading to the power spectrum of the order parameter

$P_{\phi}(\vec{k})=\int
d^{3}x\,\xi(\vec{x})\,\exp(-i\,\vec{x}\cdot\vec{k})\,.$ (5)

For large volumes the system should be rotationally symmetric such that
$P_{\phi}(k)$ and $\xi(x)$ only depend on absolute values and are related by

$\xi(x)=\frac{1}{2\pi^{2}}\int dk\,k^{2}P_{\phi}(k)\frac{\sin(kx)}{kx}\,.$ (6)

In order to build intuition and make the connection with the correlation
length $\xi_{\rm DW}$, we discuss some limiting cases of these quantities. The
contributions to the correlation function $\xi(x)$ are $\pm v^{2}$, depending
on whether the order parameter $\phi$ is in the same vacuum at $\vec{y}$ and
$\vec{x}+\vec{y}$ or not. Hence, as long as the thickness of the DWs is
neglected, one obtains $\xi(0)=v^{2}$.

For increasing $x$ the correlation function decreases. As long as
$x\ll\xi_{\rm DW}$, the curvature of the DW can be neglected and the
contribution to the correlation function depends on the probability that the
locations $\vec{y}$ and $\vec{x}+\vec{y}$ are on different sides of the DW. A
simple estimate gives

$\xi(x)\simeq v^{2}\left(1-\frac{|x|S_{\rm
DW}}{V}\right)=v^{2}\left(1-\frac{|x|}{\xi_{\rm DW}}\right)\,,$ (7)

in the limit $x\ll\xi_{\rm DW}$, while in the large scale limit $x\gg\xi_{\rm
DW}$, $\xi(x)=0$. On the other hand, the power spectrum $P_{\phi}(k)$ is
expected to go to a constant for $k\to 0$ and vanish for $k\to\infty$ (because
the DW domains have no substructures). Notice that the relation (6) seems at
first sight at odds with the fact the $\xi(x)$ has a linear term in (7). The
correlation function in (6) is even, $\xi(-x)=\xi(x)$, and expanding the sine
function seems to imply that the correlation function can only have even
powers of $x$. There is however a caveat in case the power spectrum does not
falls off faster than $k^{-4}$ in the UV. In this case the correlation
function is still even, $\xi(-x)=\xi(x)$, but the expansion of the $\sin$
function is not valid. In particular, the correlation function will develop a
kink at $x=0$.

In order to extract the correlation length, we use for the power spectrum

$P(k)\propto\frac{k_{0}^{4}}{k_{0}^{4}+\kappa k^{2}+k^{4}}\,.$ (8)

In the case $\kappa=2k_{0}$, one finds $\xi(x)=v^{2}\exp(-x/\xi_{\rm DW})$,
while in the general case, this yields a correlation function of the form

$\xi(x)=v^{2}\exp(-x/\xi_{\rm DW})\frac{\sin(x/x_{o})}{x/x_{o}}\,.$ (9)

In Fig. 3 we show the power spectrum and correlation function of the DW
network. We also show the fits according to the shapes we expect. The
extracted correlation length for the fit in (9) is $\xi_{\rm DW}=0.0499L$
compared to the value $\xi_{\rm DW}=0.05L$ that is obtained via relation (3).
The scale related to the oscillations turns out to be $x_{o}=0.035L<\xi_{\rm
DW}$ which leads to a slightly negative parameter $\kappa$ in the power
spectrum. Notice that we used a smaller correlation length in these
simulations to allow for more domains and better statistics in our
measurement. In the end, the power spectrum analysis shows that the
correlation length of the Ising model encodes all relevant information and
resembles the DW system.

Figure 2: The figure shows a slice of the simulation of the Ising model at
various time steps. The left panel shows the interface between the two domains
and nucleated bubbles while the right panel displays the kinetic energy in the
fluid. The used parameters are $v_{w}=0.8$, $L=160v_{w}/\beta$, $\xi_{\rm
DW}=0.1L$.

Figure 3: The power spectrum, $P(k)$, and correlation function $\xi(x)$ of the
Ising model in arbitrary units. The fit can be used to extract the correlation
length, $\xi_{\rm DW}=0.0499L$, very close to the value inferred from the
surface of the DW network, $\xi_{\rm DW}=0.05L$.

##  3 Simulations of first-order phase transitions

To perform hydrodynamic simulations of FOPTs, we use the _Higgsless_ approach,
as recently presented in [46]. In this approach, the scalar field is not
dynamical but acts as a background that changes the equation of state as a
bubble wall passes over a certain location during the FOPT. The fluid responds
to this change in the equation of state whereby latent energy is transferred
from the vacuum to the fluid and self-similar bubble profiles develop. The
bubble locations and times are predetermined by the bubble nucleation history
according to an exponential-in-time nucleation probability and the bubble wall
velocity is a fixed external parameter. The advantage of this approach is that
the bubble wall does not need to be resolved in the simulation as long as the
numerical algorithm remains stable and converges despite the presence of
shocks in the plasma. Overall, this provides a very efficient framework to
predict the GW spectra of FOPTs.

We calculate the GW spectrum in the simulations using the dimensionless GW
power [62]

$Q^{\prime}(q)=\frac{q^{3}\beta}{w^{2}VT}\int\frac{d\Omega_{k}}{4\pi}\left[\Lambda_{ij,kl}T_{ij}(q,\vec{k})T_{kl}^{*}(q,\vec{k})\right]_{q=k}\,,$
(10)

with $w$ representing the enthalpy in the false vacuum. The frequency and
momentum of the GW waves are represented respectively by $q$ and $\vec{k}$,
and $T$ is the simulation time. The tensor $\Lambda$ projects the energy-
momentum tensor $T_{ij}$ on its transverse traceless part. The observed GW
power spectrum is then given by

$\frac{\Omega_{\rm GW}}{Q^{\prime}}=\frac{4H\tau_{\rm
sw}}{3\pi^{2}}\frac{H}{\beta}\,,$ (11)

which is directly proportional to the lifetime of the sound waves $\tau_{\rm
sw}$ (determined either by the onset of turbulence or Hubble damping).

In the presence of DWs, we use relatively large box sizes, $L=160v_{w}/\beta$,
to make sure that we are in the regime where the mean separation between the
DWs is larger than the mean bubble size in the FOPT without the DWs 222 In the
simulations, it is most efficient to scale the box size in the different runs
with the wall velocity. This way, the number of nucleated bubbles (which is
important to limit boundary effects) is kept constant and even bubble
nucleation histories can be reused across simulations with different wall
velocities [62, 46]. Likewise, we scale the correlation length $\xi_{\rm DW}$
with the box size and hence also the wall velocity. .

Figure 4: Final spectra of the gravitational waves with (left) and without
(right) a DW network. The strength of the FOPT is $\alpha=0.05$, and the
velocities of the bubble walls are $v_{w}=0.4$, $0.55$ and $0.8$. The green
points indicate the part of the spectrum that is used in the fit (shown in
red).

As mentioned before, the FOPT then happens quickly enough that the separation
of the DWs will be imprinted in the GW spectrum. Figure 4 shows the final
spectra of the simulations. The simulations with (without) bubbles have been
run with box sizes $80v_{w}/\beta$ ($160v_{w}/\beta$) and simulation durations
$64/\beta$ ($128/\beta$) in order to accommodate the proper size of the
nucleated bubbles at percolation.

As expected, the GW signal is shifted to smaller frequencies and enhanced by a
similar factor. Moreover, the spectrum seems somewhat steeper in the IR. In
order to extract the features of the spectra, we use a similar family of
functions as in [46].

$S_{f}(q)=S_{0}\times\frac{(q/q_{0})^{n}}{1+(q/q_{0})^{(n-1)}[1+(q/q_{1})^{4}]}\times
e^{-(q/q_{e})^{2}}\,.$ (12)

Compared to [46], the Ansatz is slightly generalized to allow for an IR
behavior different than $q^{3}$. Since the FOPT is relatively strong, there is
no clear separation of the scales related to the bubble size and the shell
thickness, so $q_{0}\sim q_{1}$ in all fits and $q_{0}$ is basically the peak
of the spectrum.

Table 1 summarizes our quantitative findings. $\Omega$ denotes the frequency-
integrated GW power, while $q_{0}$ and $n$ are fitting parameters according to
(12). Quantities with a bar are computed from simulations without DWs.

$v_{w}$ | $q_{0}/\beta$ | $\xi_{\rm DW}/\beta^{-1}$ | $q_{0}\cdot\xi_{\rm DW}$ | $q_{0}/\bar{q}_{0}$ | $\Omega/{\bar{\Omega}}$ | $n$ | $\bar{n}$
---|---|---|---|---|---|---|---
0.4 | 0.50 | 6.4 | 3.19 | 3.69 | 2.37 | 5.51 | 4.87
0.55 | 0.35 | 8.8 | 3.07 | 3.95 | 2.75 | 5.03 | 4.63
0.8 | 0.22 | 12.8 | 2.81 | 3.40 | 3.44 | 4.78 | 4.14
Table 1: Parameters extracted from the spectra in Fig. 4. Quantities with a
bar are computed without DWs. The frequency $q_{0}$ is obtained fitting Eq.
(12) to simulation data. The domain wall correlation length $\xi_{\rm DW}$ is
obtained from the Ising model snapshots (notice that it scales with $v_{w}$).

The peak of the spectrum closely correlates with the correlation length of the
DW network, $q_{0}\simeq 3/\xi_{\rm DW}$. As expected, in presence of DWs, the
peak is shifted to smaller frequencies while the amplitude is enhanced by a
similar factor. In the fits, the IR tail is somewhat steeper but all
simulation boxes are in principle too small to reliably extract the slope of
IR behavior $n$.

##  4 Discussion and conclusions

We studied the gravitational wave spectrum from a first-order phase transition
seeded by domain walls. When bubbles nucleate preferentially on the domain
walls, the characteristic length scale will be the correlation length of the
domain wall network if it is larger than the typical bubble size in a
homogeneous first-order phase transition. In fact, the inverse duration of the
first-order phase transition, $\beta$, will cease to be relevant once the
correlation length of the domain walls network is large enough.

Accordingly, the correlation length of the sound waves in the plasma will be
increased. This moves the peak frequency to smaller values and enhances the
peak amplitude by a similar factor. Overall, a reasonable estimate for the
gravitational wave spectrum is given by using the established formulae for the
spectrum from a first-order phase transition and making the replacement

$\beta\to 3/\xi_{\rm DW}\simeq 4\,{\rm max}(v_{w},c_{s})/\xi_{\rm DW}\,,$ (13)

by just matching the mean bubble size of the homogeneous first-order phase
transition to the correlation length of the domain walls (see Fig. 4 and Table
1). The factor on the right side takes into account that for detonations the
bubble size relates to the wall velocity, while for deflagrations the shock
always propagates with the speed of sound. Moreover, the IR tail of the
spectrum appears a bit steeper in our simulations than in conventional first-
order phase transition. However, the IR tale of the spectrum suffers from
quite large statistical uncertainties and the effect is still tentative right
now. Also notice that the temperature of the seeded phase transition is
somewhat higher since the seeded nucleation probability dominates over the
homogeneous one [42].

Finally, let us comment on a possible application of this effect. Recently,
NanoGrav and other pulsar timing array experiments hinted at a signal of a
stochastic gravitational background [4, 5, 6, 7]. One possible explanation for
such a signal could be a first-order phase transition around a few MeV [8, 9].
However, a cosmological first-order phase transition at such low temperatures
tends to produce tension with the number of effective degrees of freedom
measured in big bang nucleosynthesis (see e.g. [63]). In many models, the
strength and duration of the first-order phase transition are correlated. If
the first-order phase transition is seeded by domain walls and the
gravitational wave spectrum peaks at the correlation length of the domain wall
network, this correlation can be broken and the first-order phase transition
temperature can be accordingly somewhat larger. This can allow explaining the
pulsar timing array excess without being in tension with big bang
nucleosynthesis constraints.

## Acknowledgment

HR acknowledges Oriol Pujolàs and Diego Blas for helpful discussions. SB is
thankful to Prateek Agrawal, Alberto Mariotti and Michael Nee for useful
discussions. SB is supported by FWO-Vlaanderen through grant number 12B2323N.
HR is supported by the Excellence Cluster ORIGINS, which is funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311. TK and IS are supported
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. We
acknowledge the support of the European Consortium for Astroparticle Theory in
the form of an Exchange Travel Grant.

##  Appendix A Statistical quantities on the bubble distribution

In this appendix, we estimate several statistical quantities related to the
bubble distribution.

The mean bubble separation can be calculated from the expected number of
bubbles $N_{b}$ nucleating in a box of volume $V$. The differential
probability $dP$ for a bubble to nucleate in the infinitesimal volume $d^{3}x$
around $\vec{x}$ in the time interval $[t,t+dt]$ is

$\displaystyle dP$
$\displaystyle={\rm(prob.~{}for~{}}\vec{x}{\rm~{}to~{}be~{}in~{}the~{}false~{}vacuum)}$
$\displaystyle\quad\times{\rm(prob.~{}for~{}one~{}bubble~{}to~{}nucleate~{}in~{}}[t,t+dt]{\rm)}$
$\displaystyle=\exp\left[-\int_{-\infty}^{t}dt^{\prime}~{}\frac{4\pi}{3}v_{w}^{3}(t-t^{\prime})^{3}\Gamma(t^{\prime})\right]\times\Gamma(t)\,dt\,d^{3}x.$
(14)

Using $\Gamma(t)=\Gamma_{*}e^{\beta(t-t_{*})}$, one obtains

$\displaystyle N_{b}$ $\displaystyle=\int dP=\frac{\beta^{3}V}{8\pi
v_{w}^{3}}.$ (15)

Thus the mean separation is

$\displaystyle R_{*}$
$\displaystyle=\left(\frac{V}{N_{b}}\right)^{\frac{1}{3}}=(8\pi)^{\frac{1}{3}}\frac{v_{w}}{\beta}\simeq\frac{3v_{w}}{\beta}.$
(16)

The average bubble radius $\langle r\rangle$ can be calculated from the
distribution of bubble radii when an arbitrary spatial point is passed by a
wall for the first time. To derive it, consider the differential probability
$dP$ for an arbitrary spatial position $\vec{x}$ to be passed for the first
time by a bubble with radius $r$ in the time interval $[t,t+dt]$. It is
written as

$\displaystyle dP$
$\displaystyle={\rm(prob.~{}for~{}}\vec{x}{\rm~{}to~{}be~{}in~{}the~{}false~{}vacuum)}$
$\displaystyle\quad\times{\rm(prob.~{}for~{}one~{}bubble~{}to~{}nucleate~{}in~{}}[t-r,t-r+dt]{\rm~{}on~{}the~{}past~{}cone)}$
$\displaystyle=\exp\left[-\int_{-\infty}^{t}dt^{\prime}~{}\frac{4\pi}{3}v_{w}^{3}(t-t^{\prime})^{3}\Gamma(t^{\prime})\right]\times
4\pi r^{2}\Gamma(t-r/v_{w})\,dt.$ (17)

Again using $\Gamma(t)=\Gamma_{*}e^{\beta(t-t_{*})}$, one obtains the radius
distribution

$\displaystyle P(r)$ $\displaystyle=\int
dP=\frac{\beta^{3}}{2v_{w}^{3}}r^{2}e^{-\beta r/v_{w}}.$ (18)

Thus the average bubble radius is, therefore, estimated as

$\displaystyle\langle r\rangle$ $\displaystyle=\int
dr~{}rP(r)=\frac{3v_{w}}{\beta}.$ (19)

##  Appendix B Variants of the Ising model

In this section, we study the effect of parameter choices for the Ising model
and DW weight assignments to assess to what extent our choices in the main
text influence our findings.

###  B.1 Ising Model parameters

In Fig. 5 we show the dependence of the GW spectrum on the value of $J/T$. In
Fig. 5, we show the dependence of the GW spectrum on the value of the Ising
model parameters, the interaction $J$ and the temperature $T$ (or more
concretely on the ratio $J/T$). The Hamiltonian of the system is given by

$H=\sum_{i,j}J\,\sigma_{i}\,\sigma_{j}\,,$ (20)

where $\sigma_{i}=\pm 1$ denotes the state on site $i$.

The orange data corresponds to $J/T=12$ while the blue data uses $J/T=2$. For
each choice, we average over three different (randomized) initial conditions
of the DW network. The error shows the variance between these three
realizations with the same parameters. The correlation length of the DW
network was chosen to be relatively large, $\xi_{\rm DW}=16v_{w}/\beta$, and
$v_{w}=0.8$, such that it dominates the bubble distribution according to (1).

Figure 5: Variants of the Ising model. The orange data corresponds to $J/T=12$
while the blue data uses $J/T=2$. For each choice, we average over three
different (randomized) initial conditions of the DW network. The error shows
the variance between these three realizations with the same parameters.

###  B.2 Cell weights assignment

Another possible modification of our methods concerns the mapping from Ising
model to a DW network and corresponding cell weight assignment. In the main
text (Sec. 2.2), we argued for a relatively simple procedure using just the
count of the states for the eight corners in each cell. An even simpler model
consists of just weighting the DW network proportionally to the number of
corners of opposite phases in the individual cells of the Ising model
simulation. For example, in the case that of the eight corners 0 to 4 corners
are of opposite states, one could simply assign the values 0 to 4 to the
weights. Various GW spectra comparing the two different methods for the
weights are shown in 6, for various parameter values. Differences are always
minute.

Figure 6: Different methods to assign DW weights to the grid cells. The method
from the main text is shown in orange, and the simplistic method described in
B.2 in blue.

## References

  * [1] LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * [2] LIGO Scientific, Virgo collaboration, _GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence_ , _Phys. Rev. Lett._ 116 (2016) 241103 [1606.04855].
  * [3] LIGO Scientific, Virgo collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101 [1710.05832].
  * [4] NANOGrav collaboration, _The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background_ , _Astrophys. J. Lett._ 905 (2020) L34 [2009.04496].
  * [5] S. Chen et al., _Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search_ , _Mon. Not. Roy. Astron. Soc._ 508 (2021) 4970 [2110.13184].
  * [6] B. Goncharov et al., _On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 917 (2021) L19 [2107.12112].
  * [7] J. Antoniadis et al., _The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background_ , _Mon. Not. Roy. Astron. Soc._ 510 (2022) 4873 [2201.03980].
  * [8] NANOGrav collaboration, _Searching for Gravitational Waves from Cosmological Phase Transitions with the NANOGrav 12.5-Year Dataset_ , _Phys. Rev. Lett._ 127 (2021) 251302 [2104.13930].
  * [9] X. Xue et al., _Constraining Cosmological Phase Transitions with the Parkes Pulsar Timing Array_ , _Phys. Rev. Lett._ 127 (2021) 251303 [2110.03096].
  * [10] eLISA collaboration, _The Gravitational Universe_ , 1305.5720.
  * [11] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry et al., _Science with the space-based interferometer LISA. V: Extreme mass-ratio inspirals_ , _Phys. Rev. D_ 95 (2017) 103012 [1703.09722].
  * [12] LISA Cosmology Working Group collaboration, _Cosmology with the Laser Interferometer Space Antenna_ , 2204.05434.
  * [13] P. J. Steinhardt, _Monopole and Vortex Dissociation and Decay of the False Vacuum_ , _Nucl. Phys. B_ 190 (1981) 583.
  * [14] P. J. Steinhardt, _Monopole Dissociation in the Early Universe_ , _Phys. Rev. D_ 24 (1981) 842.
  * [15] L. G. Jensen and P. J. Steinhardt, _DISSOCIATION OF ABRIKOSOV-NIELSEN-OLESEN VORTICES_ , _Phys. Rev. B_ 27 (1983) 5549.
  * [16] Y. Hosotani, _Impurities in the Early Universe_ , _Phys. Rev. D_ 27 (1983) 789.
  * [17] E. Witten, _Cosmic separation of phases_ , _Phys. Rev. D_ 30 (1984) 272.
  * [18] W. A. Hiscock, _CAN BLACK HOLES NUCLEATE VACUUM PHASE TRANSITIONS?_ , _Phys. Rev. D_ 35 (1987) 1161.
  * [19] D. R. Green, E. Silverstein and D. Starr, _Attractor explosions and catalyzed vacuum decay_ , _Phys. Rev. D_ 74 (2006) 024004 [hep-th/0605047].
  * [20] R. Gregory, I. G. Moss and B. Withers, _Black holes as bubble nucleation sites_ , _JHEP_ 03 (2014) 081 [1401.0017].
  * [21] P. Burda, R. Gregory and I. Moss, _Vacuum metastability with black holes_ , _JHEP_ 08 (2015) 114 [1503.07331].
  * [22] K. Mukaida and M. Yamada, _False Vacuum Decay Catalyzed by Black Holes_ , _Phys. Rev. D_ 96 (2017) 103514 [1706.04523].
  * [23] D. Canko, I. Gialamas, G. Jelic-Cizmek, A. Riotto and N. Tetradis, _On the Catalysis of the Electroweak Vacuum Decay by Black Holes at High Temperature_ , _Eur. Phys. J. C_ 78 (2018) 328 [1706.01364].
  * [24] D.-C. Dai, R. Gregory and D. Stojkovic, _Connecting the Higgs Potential and Primordial Black Holes_ , _Phys. Rev. D_ 101 (2020) 125012 [1909.00773].
  * [25] B. K. El-Menoufi, S. J. Huber and J. P. Manuel, _Black holes seeding cosmological phase transitions_ , 2006.16275.
  * [26] N. Oshita, M. Yamada and M. Yamaguchi, _Compact objects as the catalysts for vacuum decays_ , _Phys. Lett. B_ 791 (2019) 149 [1808.01382].
  * [27] R. Balkin, J. Serra, K. Springmann, S. Stelzl and A. Weiler, _Density Induced Vacuum Instability_ , 2105.13354.
  * [28] D.-C. Dai, D. Minic and D. Stojkovic, _Interaction of cosmological domain walls with large classical objects, like planets and satellites, and the flyby anomaly_ , _JHEP_ 03 (2022) 207 [2105.01894].
  * [29] J. Preskill and A. Vilenkin, _Decay of metastable topological defects_ , _Phys. Rev. D_ 47 (1993) 2324 [hep-ph/9209210].
  * [30] U. A. Yajnik, _PHASE TRANSITION INDUCED BY COSMIC STRINGS_ , _Phys. Rev. D_ 34 (1986) 1237.
  * [31] U. A. Yajnik and T. Padmanabhan, _ANALYTICAL APPROACH TO STRING INDUCED PHASE TRANSITION_ , _Phys. Rev. D_ 35 (1987) 3100.
  * [32] I. Dasgupta, _Vacuum tunneling by cosmic strings_ , _Nucl. Phys. B_ 506 (1997) 421 [hep-th/9702041].
  * [33] B. Kumar and U. A. Yajnik, _On stability of false vacuum in supersymmetric theories with cosmic strings_ , _Phys. Rev. D_ 79 (2009) 065001 [0807.3254].
  * [34] B.-H. Lee, W. Lee, R. MacKenzie, M. B. Paranjape, U. A. Yajnik and D.-h. Yeom, _Battle of the bulge: Decay of the thin, false cosmic string_ , _Phys. Rev. D_ 88 (2013) 105008 [1310.3005].
  * [35] B.-H. Lee, W. Lee, R. MacKenzie, M. B. Paranjape, U. A. Yajnik and D.-h. Yeom, _Tunneling decay of false vortices_ , _Phys. Rev. D_ 88 (2013) 085031 [1308.3501].
  * [36] I. Koga, S. Kuroyanagi and Y. Ookouchi, _Instability of Higgs Vacuum via String Cloud_ , _Phys. Lett. B_ 800 (2020) 135093 [1910.02435].
  * [37] B. Kumar and U. Yajnik, _Graceful exit via monopoles in a theory with O’Raifeartaigh type supersymmetry breaking_ , _Nucl. Phys. B_ 831 (2010) 162 [0908.3949].
  * [38] B. Kumar, M. B. Paranjape and U. A. Yajnik, _Fate of the false monopoles: Induced vacuum decay_ , _Phys. Rev. D_ 82 (2010) 025022 [1006.0693].
  * [39] P. Agrawal and M. Nee, _The Boring Monopole_ , _SciPost Phys._ 13 (2022) 049 [2202.11102].
  * [40] D. I. Dunsky, A. Ghoshal, H. Murayama, Y. Sakakihara and G. White, _Gravitational Wave Gastronomy_ , 2111.08750.
  * [41] A. Strumia, _Black holes don’t source fast Higgs vacuum decay_ , 2209.05504.
  * [42] S. Blasi and A. Mariotti, _Domain Walls Seeding the Electroweak Phase Transition_ , _Phys. Rev. Lett._ 129 (2022) 261303 [2203.16450].
  * [43] R. Jinno, T. Konstandin, H. Rubira and J. van de Vis, _Effect of density fluctuations on gravitational wave production in first-order phase transitions_ , _JCAP_ 12 (2021) 019 [2108.11947].
  * [44] A. Vilenkin, _Cosmic Strings and Domain Walls_ , _Phys. Rept._ 121 (1985) 263.
  * [45] D. Stauffer, _Scaling theory of percolation clusters_ , _Phys. Rept._ 54 (1979) 1.
  * [46] R. Jinno, T. Konstandin, H. Rubira and I. Stomberg, _Higgsless simulations of cosmological phase transitions and gravitational waves_ , 2209.04369.
  * [47] T. W. B. Kibble, _Some Implications of a Cosmological Phase Transition_ , _Phys. Rept._ 67 (1980) 183.
  * [48] W. H. Press, B. S. Ryden and D. N. Spergel, _Dynamical Evolution of Domain Walls in an Expanding Universe_ , _Astrophys. J._ 347 (1989) 590.
  * [49] T. Garagounis and M. Hindmarsh, _Scaling in numerical simulations of domain walls_ , _Phys. Rev. D_ 68 (2003) 103506 [hep-ph/0212359].
  * [50] T. Hiramatsu, M. Kawasaki and K. Saikawa, _Gravitational Waves from Collapsing Domain Walls_ , _JCAP_ 05 (2010) 032 [1002.1555].
  * [51] M. Kawasaki and K. Saikawa, _Study of gravitational radiation from cosmic domain walls_ , _JCAP_ 09 (2011) 008 [1102.5628].
  * [52] T. Hiramatsu, M. Kawasaki and K. Saikawa, _On the estimation of gravitational wave spectrum from cosmic domain walls_ , _JCAP_ 02 (2014) 031 [1309.5001].
  * [53] M. Hindmarsh, _Analytic scaling solutions for cosmic domain walls_ , _Phys. Rev. Lett._ 77 (1996) 4495 [hep-ph/9605332].
  * [54] P. Sikivie, _Of Axions, Domain Walls and the Early Universe_ , _Phys. Rev. Lett._ 48 (1982) 1156.
  * [55] S. E. Larsson, S. Sarkar and P. L. White, _Evading the cosmological domain wall problem_ , _Phys. Rev. D_ 55 (1997) 5129 [hep-ph/9608319].
  * [56] M. Gleiser and R. Roberts, _Gravitational waves from collapsing vacuum domains_ , _Phys. Rev. Lett._ 81 (1998) 5497 [astro-ph/9807260].
  * [57] K. Saikawa, _A review of gravitational waves from cosmic domain walls_ , _Universe_ 3 (2017) 40 [1703.02576].
  * [58] J. R. Espinosa, T. Konstandin and F. Riva, _Strong Electroweak Phase Transitions in the Standard Model with a Singlet_ , _Nucl. Phys. B_ 854 (2012) 592 [1107.5441].
  * [59] H. H. Patel and M. J. Ramsey-Musolf, _Stepping Into Electroweak Symmetry Breaking: Phase Transitions and Higgs Phenomenology_ , _Phys. Rev. D_ 88 (2013) 035013 [1212.5652].
  * [60] A. Vilenkin and E. P. S. Shellard, _Cosmic Strings and Other Topological Defects_. Cambridge University Press, 7, 2000.
  * [61] C. J. A. P. Martins, I. Y. Rybak, A. Avgoustidis and E. P. S. Shellard, _Extending the velocity-dependent one-scale model for domain walls_ , _Phys. Rev. D_ 93 (2016) 043534 [1602.01322].
  * [62] R. Jinno, T. Konstandin and H. Rubira, _A hybrid simulation of gravitational wave production in first-order phase transitions_ , _JCAP_ 04 (2021) 014 [2010.00971].
  * [63] N. Ramberg, W. Ratzinger and P. Schwaller, _One $\mu$ to rule them all: CMB spectral distortions can probe domain walls, cosmic strings and low scale phase transitions_, 2209.14313.


# Radiative Heat Transfer
in Free-Standing Silicon Nitride Membranes

Chang Zhang, Mathieu Giroux, Thea Abdul Nour, Raphael St-Gelais

###### Abstract

Free-standing silicon nitride (SiN) mechanical resonators are of central
interests in applications such as temperature and mass sensing, and for
fundamental optomechanical research. Understanding thermal coupling between a
membrane resonator and its environment is required for predicting thermal
noise, frequency noise, as well as sensors responses to temperature changes.
In this work, we provide closed-form derivations of intrinsic thermal coupling
quantities in free-standing thin films–namely total thermal conductance with
the surroundings, thermal response time, and the relative contribution of
thermal radiation. Our model is valid for any free-standing thin film anchored
on all sides, although we particularly emphasize the specific case of SiN for
which spectral emissivity is thoroughly investigated as a function of
thickness and temperature. We find that radiative heat exchanges can play a
non-negligible role, and even dominate thermal coupling for membranes of sizes
commonly employed in optomechanics experiments. We experimentally confirm the
validity of our model by measuring radiative thermal coupling between a SiN
mechanical resonator and a ceramic heater in high vacuum.

##  1 Introduction

Thin-film silicon nitride (SiN) membranes are heavily used as mechanical and
optical resonators in both fundamental opto-mechanical studies [1, 2, 3] and
many state-of-the-art sensing technologies, including ultra-sensitive mass
sensing [4], gas detection [5], nanoparticulate mass detection [6], thermal
radiation sensing [7, 8] and pressure sensing [9]. Some of these works involve
cooling [1, 3], while in other cases, the temperature dependence of membrane
stress is used as a sensing mechanism [8, 7], or an active technique for
controlling resonator frequency [10].

Despite, the strong influence of temperature in SiN membranes, there is still
no closed-form expression describing thermal coupling with their environment,
as well as their characteristic thermal response time ($\tau$). As such, our
goal is to provide expression for $\tau$; for the heat conduction between the
membrane and the environment ($G$, in $\mathrm{W/K}$); and for the fraction
($x_{rad}$) of this conduction that occurs via thermal radiation. We provide
these expressions for the specific case SiN films of square and circular shape
in vacuum (i.e., in the absence of convection heat transfer).

Analytical expression for thermal coupling of SiN membranes are notably needed
for understanding their fundamental noise mechanisms. Values of $G$ and $\tau$
are central in the calculation of noise processes such as temperature
fluctuation noise and temperature-induced frequency fluctuation [11, 12, 13].

Likewise, in the context of nanomechanical radiation sensors [7],
understanding the ratio of heat transfer occurring via radiation ($x_{rad}$)
is key for determining ultimate sensor detectivity and noise equivalent power
[14]. The fundamental performance limit of a thermal-based radiation sensor is
notably reached when thermal coupling between the sensor and the surroundings
is dominated by radiation (i.e., $x_{rad}\approx 1$) [7, 10, 15].

Efforts have been devoted to investigating the thermal properties [16, 17, 18]
of SiN membranes such as thermal conductivity, heat capacity, diffusivity and
thermal expansion coefficient. Emissivity of metal coating on SiN membranes
have also been investigated in the context of electron microscopy [19].
Meanwhile, others have provided finite element simulation [20] of temperature
profiles SiN membranes. Nevertheless, to the best of our knowledge, there is
still no rigorous closed-form analytical expression describing heat transfer
in such membranes.

##  2 Heat transfer in free-standing thin films

Figure 1:  (a) Square-shape SiN membrane and silicon frame. (b) Surface
temperature profile of the SiN film with uniform internal generation.
$q_{rad}$ and $q_{cond}$ denote the heat flux transferred via radiation and
conduction, respectively.

In order to calculate thermal coupling between a freestanding film and its
environment, we consider volumetric heat generation
($\dot{q},\textrm{in}\mathrm{W/m^{3}}$) occurring uniformly within the
membrane. From a heat transfer standpoint, such internal generation is
mathematically analogous to uniform absorption of radiation from an external
radiation source (e.g., in the context of radiation sensors). As a result of
this internal generation, membrane temperature ($T_{m}$) rises relative to
ambient ($T_{\infty}$), and heat leaves the membrane by conduction to the
silicon frame ($q_{cond}$, in W), and re-emission of radiation to the
environment ($q_{rad}$, in W), as schematized in Fig. 1(b). For example, an
extreme case in which the membrane would be perfectly coupled to its
environment via radiation would yield $q_{cond}=0$, and $q_{rad}=q$, where
$q=\dot{q}V$ and $V$ is the membrane volume. We neglect the contribution of
convective heat transfer, an assumption representative of a system in vacuum
as in most optomechanics and thermal radiation sensing experiments. In steady
state, heat flux out of the system must equal heat generation inside the
system, yielding the general heat equation:

$-k\cdot\nabla^{2}T_{m}+\dot{q}_{rad}=\dot{q},$ (1)

where $k$ is the membrane conductivity, and $T_{m}$ is the position-dependent
membrane temperature. In Eq. (1), $\dot{q}_{rad}$ is the radiative exchange,
per unit volume, between membrane and the environment at $T_{\infty}$; is it
given by:

$\dot{q}_{rad}=\frac{2\sigma\epsilon(T_{m}^{4}-T_{\infty}^{4})}{d},$ (2)

Where $d$ is the membrane thickness, $\sigma$ is Stefan-Boltzmann constant,
and $\epsilon$ is the total hemispherical emissivity of SiN. The factor 2
accounts for emission on both faces of the membrane. We consider Dirichlet
boundary conditions, setting the temperature of the SiN membrane edges to be
the same as the ambient temperature ($T_{\infty}$).

The non-linear nature of radiative heat exchange ($\dot{q}_{rad}\propto
T_{m}^{4}$) prevents direct derivation of a closed-form solution of the heat
equation. We consequently linearize Eq. (2) by considering a small temperature
difference between the membrane and the environment, yielding:

$\dot{q}_{rad}\approx\frac{8\sigma\epsilon
T_{\infty}^{3}}{d}(T_{m}-T_{\infty}).$ (3)

Due to the minute scale thickness of the SiN film, we consider uniform
temperature along the direction normal to the surface, thus reducing Eq. (1)
to a two-dimensional problem.

For a circular membrane of radius $r_{0}$, the solution to the linearized heat
equation in cylindrical coordinates is conveniently simple:

$T_{m}(r)=[1-\frac{I_{0}(\beta\cdot r)}{I_{0}(\beta\cdot
r_{0})}]\cdot\frac{\dot{q}}{k\cdot\beta^{2}}+T_{\infty},$ (4)

where

$\beta=\sqrt{\frac{8\sigma\epsilon T_{\infty}^{3}}{kd}},$ (5)

and $I_{N}$ is the $N^{th}$ order modified Bessel function of the first kind.
From this temperature profile, we calculate heat transfer by conduction at the
boundaries, using Fourier law of conduction:

$q_{cond}=-2k\pi r_{0}d\cdot\frac{\partial T_{m}}{\partial
r}\Bigr{|}_{r_{0}},$ (6)

which yields:

$q_{cond}=\frac{2q}{\beta\cdot r_{0}}\cdot\frac{I_{1}(\beta\cdot
r_{0})}{I_{0}(\beta\cdot r_{0})}.$ (7)

From Eq. (7), we can finally calculate the fraction of heat that leaves the
membrane radiation ($x_{rad}$):

$x_{rad}=\frac{q_{rad}}{q}=1-\frac{q_{cond}}{q}=1-\frac{2}{\beta\cdot
r_{0}}\cdot\frac{I_{1}(\beta\cdot r_{0})}{I_{0}(\beta\cdot r_{0})}.$ (8)

In Eq. (8), we note that the right-hand side depends only on intrinsic
membrane properties (i.e., independent of $\dot{q}$). As such, solving for
$x_{rad}>0.5$ yields the properties required for a membrane to be thermally
coupled to its environment more strongly via radiation than via solid-state
conduction.

From the temperature profile of the membrane [see Eq. (4)], we can also
express the thermal time constant of the membrane ($\tau$, in s) and its
overall thermal conductance with the environment ($G$,
$\textrm{in}\mathrm{W/K}$), both of which are of particular importance for use
in sensors [14], and for predicting noise profiles in micro resonators [11,
12, 13]. Note that $G=G_{cond}+G_{rad}$ includes heat transfer both by
conduction in the supporting frame ($G_{cond}$) and by radiation ($G_{rad}$).
We obtain $G$ and $\tau$ by evaluating the average temperature
($\overline{T}_{m}$) of the membrane [i.e., by integrating Eq. (4)] and using:

$\displaystyle G=\frac{q}{\overline{T}_{m}-T_{\infty}}=\frac{k\pi(\beta\cdot
r_{0})^{2}d}{1-\frac{2}{\beta\cdot r_{0}}\cdot\frac{I_{1}(\beta\cdot
r_{0})}{I_{0}(\beta\cdot r_{0})}},$ (9) $\displaystyle\tau=\frac{c_{p}\rho
V}{G}=\frac{c_{p}\rho[1-\frac{2}{\beta\cdot r_{0}}\cdot\frac{I_{1}(\beta\cdot
r_{0})}{I_{0}(\beta\cdot r_{0})}]}{k\beta^{2}},$ (10)

where $c_{p}$ and $\rho$ are respectively the specific heat capacity and the
material density. We note that, when heat transfer becomes entirely dominated
by radiation (i.e. for large areas), Bessel terms become negligible and both
expressions simplify to radius-independent quantities:

$\displaystyle G_{rad}=8A\sigma\epsilon T_{\infty}^{3},$ (11)
$\displaystyle\tau_{rad}=\frac{\rho c_{p}d}{8\sigma\epsilon T_{\infty}^{3}}.$
(12)

For a square membrane, a geometry much more frequently encountered in
practice, an analytical solution to the linearized heat equation also exists,
but comprises an infinite amount of Fourier terms in order to respect the
boundary conditions ($T=T_{\infty}$ at the membrane edges). Rather than using
this complex solution, we numerically solve for the temperature profile of a
square membrane using finite element analysis. From this solution (see
supplemental Fig. S1), we determine that the radiative thermal coupling ratio
($x_{rad}=q_{rad}/q$) of a square membrane of side length $L$ matches that of
a circular membrane if we consider an effective radius:

$r_{eff}=1.252\frac{L}{2}.$ (13)

Logically, this effective radius falls between half of a square membrane side
length ($L/2$), and half of its diagonal ($L/\sqrt{2}$), as shown
schematically in Fig. 1(a). By replacing Eq. (13) in Eq. (8), the fraction of
heat transfer occurring by radiation in a square membrane ($x_{rad}$) is given
by:

$x_{rad}=1-\frac{2}{\beta\cdot r_{eff}}\cdot\frac{I_{1}(\beta\cdot
r_{eff})}{I_{0}(\beta\cdot r_{eff})}.$ (14)

The total thermal conductance of a square membrane is subsequently:

$G=\frac{G_{rad}}{x_{rad}}=\frac{8L^{2}\sigma\epsilon
T_{\infty}^{3}}{x_{rad}},$ (15)

while $\tau$ is now given by:

$\tau=\frac{c_{p}\rho d}{8\sigma\epsilon T_{\infty}^{3}}x_{rad}.$ (16)

Eq. (16) yields an aberrant $\tau\approx 0$ when $x_{rad}\approx 0$, but this
result occurs for $L\approx 0$, in which case our assumption of a 2-D problem
(i.e., $L>>d$) does not hold.

##  3 Emissivity of SiN films

From Eq. (8) to Eq. (16), it is obvious that the total hemispherical
emissivity ($\epsilon$) is a key parameter governing radiative heat transfer
in free-standing membranes. We calculate $\epsilon$ for the specific case of
SiN using Kirchhoff law–i.e., directional spectral emissivity is equal to
absorption:

$\epsilon_{\lambda,\theta}(\lambda,\theta)=\alpha_{\lambda,\theta}(\lambda,\theta).$
(17)

We calculate
$\alpha_{\lambda,\theta}(\lambda,\theta)=1-R_{\lambda,\theta}(\lambda,\theta)-T_{\lambda,\theta}(\lambda,\theta)$
using conventional optical multi-layer calculation [21] with the complex
permittivity of SiN taken from [22]. $R$ and $T$ denote the optical power
reflection and transmission coefficients, respectively. We find that these
coefficients depend non-negligibly on the angle, such that we integrate the
emissivity according to Lambert’s cosine law to obtain [23]:

$\epsilon_{\lambda}(\lambda)=2\int_{0}^{\pi/2}\epsilon_{\lambda,\theta}(\lambda,\theta)\cdot
cos(\theta)\cdot sin(\theta)d\theta,$ (18)

Figure 2:  (a) Comparison between normal and hemispherical emissivity for
free-standing SiN membranes. (b) Normalized emission spectrum of a 100 nm
thick SiN membrane at different temperatures. Dashed lines represent blackbody
spectra. As temperature decreases, the peak of the emission spectrum moves
towards longer wavelengths. (c) Hemispherical total emissivity of free-
standing SiN membranes at different temperatures. (d) Hemispherical total
emissivity normalized by film thicknesses for different temperatures. This
quantity scales with of the fraction of heat transfer occurring by radiation
($x_{rad}$) and is maximized for lower membrane thicknesses.

which is plotted in Fig. 2(a) for the specific example of a 100 nm thick
membrane. In comparison, the normal directional spectral emissivity
$\epsilon_{\lambda,\theta}(\lambda,0)$ is roughly $25\%$ higher than the
integrated value $\epsilon_{\lambda}(\lambda)$ . In other words, emission and
absorption are slightly stronger at normal incidence. This may be beneficial
for radiation sensing applications [14] where a sensor should be more strongly
coupled at normal incidence (i.e. pointing at the object to be detected). We
finally obtain the total hemispherical emissivity $\epsilon$ by weighting
$\epsilon_{\lambda}(\lambda)$ with the blackbody emission spectrum
$E_{\lambda,b}(\lambda,T)$ at temperature $T$:

$\epsilon(T)=\frac{\int_{0}^{\infty}\epsilon_{\lambda}(\lambda)\cdot
E_{\lambda,b}(\lambda,T)d\lambda}{E_{b}(T)},$ (19)

where $E_{b}(T)=\sigma T^{4}$. For concision, we use the notation
$\epsilon(T)=\epsilon$ in this work. This weighting is shown, for various
membrane temperatures, in Fig. 2(b) while $\epsilon$ as a function of
thicknesses, is presented in Fig. 2(c). We find that both the emission
distribution [Fig. 1(b)] and $\epsilon$ [Fig. 1(c)] weakly depend on
temperature unless cryogenic membrane temperatures (e.g., 100 K) are
considered.We willingly omit calculations at very low temperatures
($<1\,\mathrm{K}$) as this would presumably require a different model for
material properties of SiN. We also note that the hemispherical emissivity for
thin films is a strong function of thickness. This emphasizes that the common
assumption of $\epsilon\approx 0.6$ [7] for SiN appears appropriate only for
bulk SiN materials and not for thin films.

We note that thicker films lead to higher emissivity [see Fig. 2(c)]–however
as shown in Eq. (5), the relevant quantity for estimating the contribution of
radiation to the total heat transfer is the ratio $\epsilon/d$. This is
plotted in Fig. 2(d), from which we conclude that, for a given surface area,
thinner membranes are more efficiently coupled via radiation. The thickness
dependence is the strongest in the 100 nm–1 $\mathrm{\mu m}$ range, while the
relation reaches a relative plateau for thickness commonly employed in opto-
mechanics experiments ($d<100\,\mathrm{nm}$) [24].

Figure 3:  (a) Fraction of total heat transfer occurring via radiation
($x_{rad}$) in free-standing membranes as a function of $\beta\cdot r_{eff}$.
The model applies both to circular ($r_{eff}=r_{0}$) and square
($r_{eff}=1.252L/2$) membranes. For $x_{rad}$ larger than 0.5, membrane
thermal coupling is radiation-dominated. (b) $\beta$ as a function of membrane
thickness for the specific case of SiN

Having calculated the $\epsilon/d$ ratio, we can determine the critical
membrane length for which heat transfer is dominated by radiation (i.e.,
$x_{rad}>0.5$). Using the closed-form relation developed in Eq. (14), we plot
$x_{rad}$ as a function of $\beta\cdot r_{eff}$ in Fig. 3(a). The $\beta$
value for SiN, for given $d$ and $T_{m}$, is given in Fig. 3(b). We note that
Fig. 3(a) is a universal relation for any freestanding thin film membrane
anchored on all sides, while Fig. 3(b) accounts for the specific case of SiN.
In Fig. 3(a), we find that thermal coupling of free-standing thin films is
radiation dominated for $\beta\cdot r_{eff}\approx 3.33$. From this, we
finally derive a simple expression for the threshold of a radiation-coupled
membrane:

$r_{eff}>\frac{3.33}{\beta}=3.33\sqrt{\frac{kd}{8\sigma\epsilon
T_{\infty}^{3}}},$ (20)

for which the $\epsilon/d$ ratio can be obtained graphically from Fig. 2(d),
in the specific case of SiN. We note that this equation is the same for a
circular membrane, using $r_{eff}=r_{0}$ in lieu of $r_{eff}=1.252L/2$ for a
square membrane.

##  4 Experimental results

We verify our model by conducting experiments using a commercially available
($3\times 3$ mm side length, 200 nm thickness) low stress SiN membrane. We
correlate the membrane temperature to its mechanical resonance frequency
($f_{r}$) using the relation given in [10]. The change in temperature of the
membrane ($\Delta T$) is measured by comparing its instantaneous resonance
($f_{r}$) with its initial room-temperature resonance ($f_{0}$):

$\frac{f_{r}}{f_{0}}=\sqrt{1-\frac{E\alpha}{\sigma_{0}}\Delta T}$ (21)

where SiN Young’s modulus is $E=300\,\textrm{GPa}$, the thermal expansion
coefficient is $\alpha=3.27\times 10^{-6}\,\mathrm{K^{-1}}$ [18] and the
tensile stress is $\sigma_{0}\approx 100\,\mathrm{Mpa}$. We consider a $25\%$
uncertainty for this relation given the variability of material constants and
on the membrane dimensions [see error bars in Fig. 4(d)]. The experiment is
conducted in high vacuum ($1.5\times 10^{-6}\,\mathrm{torr}$) to eliminate
convection heat transfer and viscous damping by air.

Figure 4:  (a) Schematic of temperature-dependent frequency shift mechanism
for SiN membranes. (b) Photograph of the experimental setup inside the vacuum
chamber. For the actual experiment, the heater is placed closer (5 mm) from
the membrane but is displaced here for a clearer picture (c) Resonance
frequency of mechanical mode 2, 2 as a function of heater temperature. (d)
Membrane temperature as a function of heater temperature. Experimental results
agree with our model within a 25% error on the membrane material constants.
Predicted theoretical results for $x_{rad}=0.5$ and 1 are also shown for
illustrative purpose, for the same $F=0.55$ view factor. Shaded areas are
bounded by the range $\epsilon_{heater}=0.6$ and 0.8.

We vary the temperature of the SiN membrane by exposing it to a rough-surfaced
aluminum oxide heater [see Fig 4.(b)] placed within a short distance (5 mm)
from the membrane. We infer the heater temperature by correlating it to its
electrical resistance using a temperature coefficient of resistance of
$4.7\times 10^{-3}\,\mathrm{K^{-1}}$. We measure this value in a separate
experiment by placing the heater on a hot plate and by measuring its
resistance as a function of temperature. Due to the variability in documented
value for the total hemispherical emissivity of rough-surfaced aluminum oxide
[23], we assume an upper bound value of $\epsilon_{heater}=0.8$ and a lower
bound of 0.6. A reflective aluminum shield is placed between the membrane and
the heater to prevent the silicon frame of the membrane from heating up upon
absorption of radiation. A piezoelectric actuator is attached to the same
glass slide as the membrane to excite its mechanical resonance. Membrane
displacement is measured using an optical interferometer setup [25]. The
instantaneous resonance ($f_{r}$) of the SiN membrane shifts by over 1 KHz
when subject to a 70 K increment in heater temperature, as presented in Fig.
4(c).

By neglecting photons having more than one interaction with the heater and
membrane (due to the relatively small size of the membrane, the diffusive
surface of the heater, and its high emissivity/absorptivity) the temperature
of the membrane is correlated to the heater temperature simply by:

$\Delta T_{m}\approx\frac{\Delta T_{heater}}{2}\cdot x_{rad}\cdot
F\cdot\epsilon_{heater}$ (22)

where, $F=0.55$ is the geometrical view factor [23]. The detailed thermal
equivalent circuit from which Eq. (22) derived is given in supplementary Fig.
S2. From Eq. (22), correlating the membrane temperature with the heater
temperature allows measurement of $x_{rad}$ and validation of our model. This
is presented in Fig. 4(d), where the experimental points fall within the
theoretically expected values from Eq. (14). For illustrative purpose, we also
plot the projected values for a membrane that would be on the threshold
($x_{rad}=0.5$) or completely radiation-dominated ($x_{rad}=1$), for the same
geometrical view factor ($F=0.55$).

##  5 conclusion

We expect that our work will be interest for achieving high performance
radiation sensors exploiting the high temperature sensitivity of SiN
mechanical resonance [10]. In such sensors, one would ideally want to achieve
$x_{rad}\approx 1$ to reach the highest possible detectivity [14]. Our work
shows that achieving such $x_{rad}$ value is feasible using realistic membrane
dimensions. We readily achieve $x_{rad}=0.18$ using commercially available SiN
membranes of non-optimized dimensions. We also expect that the provided
closed-form expressions for $\tau$ and $G$ will be of great use for predicting
frequency noise in high Q-factor SiN [11, 12, 13]. Given the high temperature
sensitivity of resonance frequency in SiN resonators, we expect temperature
fluctuation noise–which is directly linked to $\tau$ and $G$ [13]–to have a
non-negligible contribution to frequency noise. Finally, we expect that
outlining the non-negligible contribution of radiation heat transfer may be
useful in experiments involving cooling of SiN membranes. As a striking
example, a membrane with $x_{rad}>0.5$ could be more efficiently cooled by a
cold object facing it, than by direct contact cooling of its supporting
silicon frame.

## References

  * Wilson et al. [2009] D. J. Wilson, C. A. Regal, S. B. Papp, and H. J. Kimble.  Cavity optomechanics with stoichiometric sin films.  _Phys. Rev. Lett._ , 103:207204, Nov 2009.  doi:10.1103/PhysRevLett.103.207204.
  * Regal and Lehnert [2011] C. A. Regal and K. W. Lehnert.  From cavity electromechanics to cavity optomechanics.  In _Journal of Physics: Conference Series_ , 2011.  doi:10.1088/1742-6596/264/1/012025.
  * Chakram et al. [2014] S. Chakram, Y. S. Patil, L. Chang, and M. Vengalattore.  Dissipation in ultrahigh quality factor sin membrane resonators.  _Phys. Rev. Lett._ , 112:127201, Mar 2014.  doi:10.1103/PhysRevLett.112.127201.
  * Hanay et al. [2012] M. S. Hanay, S. Kelber, A. K. Naik, D. Chi, S. Hentz, E. C. Bullard, E. Colinet, L. Duraffourg, and M. L. Roukes.  Single-protein nanomechanical mass spectrometry in real time.  _Nature Nanotechnology_ , 7(9):602–608, 2012\.  ISSN 17483395.  doi:10.1038/nnano.2012.119.
  * Brucoli et al. [2014-02-24] Giovanni Brucoli, Patrick Bouchon, Riad Haïdar, Mondher Besbes, Henri Benisty, and Jean-Jacques Greffet.  High efficiency quasi-monochromatic infrared emitter.  _Applied Physics Letters_ , 104(8), 2014-02-24.  ISSN 0003-6951.
  * Larsen et al. [2013] Tom Larsen, Silvan Schmid, Luis G. Villanueva, and Anja Boisen.  Photothermal analysis of individual nanoparticulate samples using micromechanical resonators.  _ACS Nano_ , 7(7):6188–6193, 2013.  ISSN 19360851.  doi:10.1021/nn402057f.
  * Zhang et al. [2013] X. C. Zhang, E. B. Myers, J. E. Sader, and M. L. Roukes.  Nanomechanical torsional resonators for frequency-shift infrared thermal sensing.  _Nano Letters_ , 13(4):1528–1534, 2013.  ISSN 15306984.  doi:10.1021/nl304687p.
  * Zhang et al. [2019] C Zhang, M Giroux, T A Nour, and R St-Gelais.  Thermal radiation sensing using high mechanical Q-factor silicon nitride membranes.  In _2019 IEEE SENSORS_ , pages 1–4, oct 2019.  doi:10.1109/SENSORS43011.2019.8956551.
  * Zhu et al. [2013] Shou En Zhu, Murali Krishna Ghatkesar, Chao Zhang, and G. C.A.M. Janssen.  Graphene based piezoresistive pressure sensor.  _Applied Physics Letters_ , 102(16):2011–2014, 2013.  ISSN 00036951.  doi:10.1063/1.4802799.
  * St-Gelais et al. [2019] Raphael St-Gelais, Simon Bernard, Christoph Reinhardt, and Jack C. Sankey.  Swept-Frequency Drumhead Optomechanical Resonators.  _ACS Photonics_ , 2019.  ISSN 23304022.  doi:10.1021/acsphotonics.8b01519.
  * Sansa et al. [2016] Marc Sansa, Eric Sage, Elizabeth C. Bullard, Marc Gély, Thomas Alava, Eric Colinet, Akshay K. Naik, Luis Guillermo Villanueva, Laurent Duraffourg, Michael L. Roukes, Guillaume Jourdan, and Sébastien Hentz.  Frequency fluctuations in silicon nanoresonators.  _Nature Nanotechnology_ , 2016.  ISSN 17483395.  doi:10.1038/nnano.2016.19.
  * Fong et al. [2012] King Y. Fong, Wolfram H. P. Pernice, and Hong X. Tang.  Frequency and phase noise of ultrahigh $q$ silicon nitride nanomechanical resonators.  _Phys. Rev. B_ , 85:161410, Apr 2012.  doi:10.1103/PhysRevB.85.161410.
  * Cleland and Roukes [2002] A. N. Cleland and M. L. Roukes.  Noise processes in nanomechanical resonators.  _Journal of Applied Physics_ , 2002.  ISSN 00218979.  doi:10.1063/1.1499745.
  * Rogalski [2003] Antoni Rogalski.  Infrared detectors: Status and trends.  _Progress in Quantum Electronics_ , 27(2-3):59–210, 2003.  ISSN 00796727.  doi:10.1016/S0079-6727(02)00024-1.
  * Blaikie et al. [2019] Andrew Blaikie, David Miller, and Benjamín J. Alemán.  A fast and sensitive room-temperature graphene nanomechanical bolometer.  _Nature Communications_ , 2019.  ISSN 20411723.  doi:10.1038/s41467-019-12562-2.
  * Sikora et al. [2012] A. Sikora, H. Ftouni, J. Richard, C. Hébert, D. Eon, F. Omnès, and O. Bourgeois.  Highly sensitive thermal conductivity measurements of suspended membranes (SiN and diamond) using a 3$\omega$-Völklein method, 2012.  ISSN 00346748.
  * Zhang and Grigoropoulos [1995] Xiang Zhang and Costas P. Grigoropoulos.  Thermal conductivity and diffusivity of free-standing silicon nitride thin films.  _Review of Scientific Instruments_ , 66(2):1115–1120, 1995.  ISSN 00346748.  doi:10.1063/1.1145989.
  * Tien and Lin [2012] Chuen-Lin Tien and Tsai-Wei Lin.  Thermal expansion coefficient and thermomechanical properties of SiN x thin films prepared by plasma-enhanced chemical vapor deposition .  _Applied Optics_ , 2012.  ISSN 1559-128X.  doi:10.1364/ao.51.007229.
  * Van Zwol et al. [2015] P. J. Van Zwol, D. F. Vles, W. P. Voorthuijzen, M. Péter, H. Vermeulen, W. J. Van Der Zande, J. M. Sturm, R. W.E. Van De Kruijs, and F. Bijkerk.  Emissivity of freestanding membranes with thin metal coatings.  _Journal of Applied Physics_ , 2015.  ISSN 10897550.  doi:10.1063/1.4936851.
  * Revaz et al. [2003] B. Revaz, B. L. Zink, D. O’Neil, L. Hull, and F. Hellman.  Numerical simulation of the heat transfer in amorphous silicon nitride membrane-based microcalorimeters.  _Review of Scientific Instruments_ , 2003.  ISSN 00346748.  doi:10.1063/1.1605498.
  * Edalatpour and Francoeur [2013-03] Sheila Edalatpour and Mathieu Francoeur.  Size effect on the emissivity of thin films.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ , 118:75,85, 2013-03.  ISSN 0022-4073.
  * Cataldo et al. [2012] Giuseppe Cataldo, James A Beall, Hsiao-mei Cho, Brendan Mcandrew, Michael D Niemack, and Edward J Wollack.  Infrared dielectric properties of low-stress silicon nitride.  37(20):4200–4202, 2012.
  * L.Bergman and S.Lavine [2017] Theodore L.Bergman and Adrienne S.Lavine.  _Fundamentals of Heat and Mass Transfer_.  Wiley, 8th editio edition, 2017.  ISBN 9781119330103.
  * Zwickl et al. [2008] B. M. Zwickl, W. E. Shanks, A. M. Jayich, C. Yang, A. C.Bleszynski Jayich, J. D. Thompson, and J. G.E. Harris.  High quality mechanical and optical properties of commercial silicon nitride membranes.  _Applied Physics Letters_ , 92(10):2006–2009, 2008.  ISSN 00036951.  doi:10.1063/1.2884191.
  * Rugar et al. [1989] D. Rugar, H. J. Mamin, and P. Guethner.  Improved fiber-optic interferometer for atomic force microscopy.  _Applied Physics Letters_ , 1989.  ISSN 00036951.  doi:10.1063/1.101987.

## Supplementary Information

###  S.1 Numerical simulation

The 2-D steady-state temperature profile of the squre-sized SiN membrane is
simulated using MATLAB PDE solver, such that the closed-form analytical
results can be compared with the numerical simulation. We first define the
computational geometry to be a square which has the same side length (3 mm) as
the sample SiN membrane. We then specify the coefficients for the PDE model by
rearranging the heat equation as:

Figure S1:  (a) 2-D temperature simulation for $3\times 3$mm, 200 nm SiN
membrane. (b) $x_{rad}$ as a function of side length which discrete points are
given by the numerical simulation, whereas the continuous lines are plotted by
closed-form expression.

$-\nabla\cdot(\nabla
T_{m})+{\beta^{2}}T_{m}=\frac{\dot{q}}{k}+\beta^{2}T_{\infty}$ (S1)

where $\sigma$ is the Stefan-Boltzmann constant, $\epsilon$ is the total
hemisphereical emissivity of the 200 nm ($d$) SiN membrane which is calculate
to be 0.11, $T_{\infty}$ is the ambient temperature (293 K), $k$ is the
conductivity of the SiN membane which is found to be 12
$\mathrm{W}/\mathrm{m}\cdot\mathrm{K}$. Note that $\dot{q}$ is the heating
power being absorbed by the SiN membrane per unit volume. Here we arbitrarily
set $\dot{q}=6.5\times 10^{7}\mathrm{W/m^{3}}$. The simulated temperature
profile is shown in Fig. S1(a).

With the numerical simulation, we can obtain the fraction of heat that leaves
the membrane by radiation ($x_{rad}$) using the simulated temperature profile.
By plotting this result and the analytical results for many side lengths
values, we find that the correction ratio between the effective radius and
half of the side length to be 1.252 which is presented in Eq. (13)

###  S.2 Thermal equivalent circuit

Figure S2:  Thermal equivalent circuit considering the heater, the membrane
and the vacuum environment. (a) Thermal circuit schematic. (b) Simplified
thermal circuit.

The analogy of the electric circuit is useful to understand the heat transfer
of the SiN membrane with its environment. Here we neglect multiple heat
exchanges between the SiN membrane and the heater by only considering single
photon interaction. In other words, photons leaving the membrane and reflected
by the surface of the heater will be directed to the environment, which
results in resistance $R_{\infty 2}$ in Fig. S2. This assumption is based on
the fact that the surface of the heater is highly diffusive, and the surface
area of the membrane is substantially smaller than the area of the heater.
Hence, photons emitted by the membrane and then hit the surface of the heater
have minimal chance of being reflected back to the membrane.

By rearranging the thermal circuit schematic in Fig. S2(a), we reach the
simplified thermal circuit shown in Fig. S2(b) for which the resistances are
expressed as:

$R_{h}=\frac{1}{\theta_{rad}\epsilon_{m}\epsilon_{h}F},$ (S2)
$R_{\infty,b}=\frac{1}{\theta_{rad}\epsilon_{m}},$ (S3)
$R_{\infty,1}=\frac{1}{\theta_{rad}\epsilon_{m}(1-F)},$ (S4)
$R_{\infty,2}=\frac{1}{\theta_{rad}\epsilon_{m}F(1-\epsilon_{h})},$ (S5)
$\theta_{rad}=4\sigma A_{m}T_{\infty}^{3}.$ (S6)

Here, $\epsilon_{m}$ is the total hemispherical emissivity of the SiN
membrane, $\epsilon_{h}$ is the total hemispherical emissivity of the ceramic
heater, $F$ is the geometrical view factor from the membrane to the heater and
$A_{m}$ is surface area of the SiN membrane. By following the thermal circuit
in Fig. S2(b), we obtain Eq. (22) in the main text.