# Self-consistent treatment of thermal effects in neutron-star post-mergers:
observational implications for third-generation gravitational-wave detectors

Verónica Villa-Ortega  Instituto Galego de Física de Altas Enerxías,
Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galicia,
Spain  Ana Lorenzo-Medina  Instituto Galego de Física de Altas Enerxías,
Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galicia,
Spain  Juan Calderón Bustillo ∗ Instituto Galego de Física de Altas Enerxías,
Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galicia,
Spain  Department of Physics, The Chinese University of Hong Kong, Shatin,
N.T., Hong Kong  Milton Ruiz† Departament d’Astronomia i Astrofísica,
Universitat de València, C/ Dr Moliner 50, 46100, Burjassot (València), Spain
Davide Guerra  Departament d’Astronomia i Astrofísica, Universitat de
València, C/ Dr Moliner 50, 46100, Burjassot (València), Spain  Pablo Cerdá-
Durán  Departament d’Astronomia i Astrofísica, Universitat de València, C/ Dr
Moliner 50, 46100, Burjassot (València), Spain  José A. Font  Departament
d’Astronomia i Astrofísica, Universitat de València, C/ Dr Moliner 50, 46100,
Burjassot (València), Spain  Observatori Astronòmic, Universitat de València,
C/ Catedrático José Beltrán 2, 46980, Paterna (València), Spain

###### Abstract

We assess the impact of accurate, self-consistent modelling of thermal effects
in neutron-star merger remnants in the context of third-generation
gravitational-wave detectors. This is done through the usage, in Bayesian
model selection experiments, of numerical-relativity simulations of binary
neutron star (BNS) mergers modelled through: a) nuclear, finite-temperature
(or “tabulated”) equations of state (EoSs), and b) their simplifed piecewise
(or “hybrid”) representation. These cover four different EoSs, namely SLy4,
DD2, HShen and LS220. Our analyses make direct use of the Newman-Penrose
scalar $\psi_{4}$ outputted by numerical simulations. Considering a detector
network formed by three Cosmic Explorers, we show that differences in the
gravitational-wave emission predicted by the two models are detectable with a
natural logarithmic Bayes Factor $\log{\cal{B}}\geq 5$ at average distances of
$d_{L}\simeq 50$Mpc, reaching $d_{L}\simeq 100$Mpc for source inclinations
$\iota\leq 0.8$, regardless of the EoS. This impact is most pronounced for the
HShen EoS. For low inclinations, only the DD2 EoS prevents the detectability
of such modelling differences at $d_{L}\simeq 150$Mpc. Our results suggest
that the usage a self-consistent treatment of thermal effects is crucial for
third-generation gravitational wave detectors.

Neutron-star mergers – Gravitational waves – Numerical relativity – Bayesian
inference

##  1 Introduction

Mergers of stellar compact binaries, having at least a neutron star as one of
the components, lie at the intersection of several very active fields of
research in relativistic astrophysics. They are prime targets for multi-
messenger observations, being visible in both gravitational wave (GW) and
electromagnetic (EM) radiation. The LIGO-Virgo-KAGRA (LVK) Collaboration has
reported observations from both types of systems (see e.g. Abbott et al.
(2021)). In particular, the detection of the first binary neutron star (BNS)
merger, GW170817, along with its post-merger EM emission (Abbott et al.,
2017a, b, c), has been instrumental to begin addressing some long-standing
questions. It has placed tight constraints on: i) the equation of state (EoS)
at supranuclear densities (Shibata et al., 2017; Margalit & Metzger, 2017;
Rezzolla et al., 2018; Ruiz et al., 2018); ii) the radius and tidal
deformability of a spherical neutron star (Bauswein et al., 2017; Most et al.,
2018); and iii) the speed of GW, which in turn disfavours a large class of
scalar-tensor theories and other theories predicting varying GW speed
(Ezquiaga & Zumalacárregui, 2017; Pardo et al., 2018). Besides, it also
provided an independent measure for the expansion of the Universe (Abbott et
al., 2017d; Dietrich et al., 2020) and the most direct evidence that BNS
mergers are progenitors of the central engines of short gamma-ray bursts
(sGRBs), followed by a longer optical transient afterglow known as a kilonova,
powered by the radioactive decay of heavy r-process nuclei (Metzger et al.,
2010; Cowperthwaite et al., 2017; Kasen et al., 2017; Pian et al., 2017).

Many of those results are drawn from the analysis of the information encoded
in the GW signal at the late BNS inspiral stage. While searches have been
conducted by the LVK Collaboration (Abbott et al., 2017e) no observational
evidence of the post-merger GW signal is yet available, as the sensitivity of
current detectors at the kHz frequency range is severely hampered by photon
shot noise. Estimates on the detectability of the post-merger signal with
third-generation detectors have been reported by Clark et al. (2016) and
Miravet-Tenés et al. (2023a).

The BNS inspiral is accurately modelled assuming a cold (zero temperature)
EoS. However, during and after merger, shocks heat up the binary remnant to
temperatures $\gtrsim 10~{}\rm MeV$ adding a thermal pressure support that may
change the internal structure of the remnant and its subsequent evolution. It
is expected that these thermal effects are encoded in the emitted GWs (Raithel
et al., 2021). Therefore, observations during the post-merger stage may shed
light on the microphysical EoS of NSs at finite temperature. These
observations may become possible with third-generation GW observatories
(Chatziioannou et al., 2017).

Advances in our knowledge of the dynamics of BNS mergers and the subsequent
post-merger evolution rely on numerical relativity simulations (for recently
reviews see e.g. Baiotti & Rezzolla (2017); Shibata & Hotokezaka (2019); Sarin
& Lasky (2021); Bauswein et al. (2020); Ruiz et al. (2021a)). In those
simulations NS matter is described using mainly two approaches. The first one
is a “hybrid” approach that assumes that pressure and internal energy can be
divided into two contributions (Janka et al., 1993; Dimmelmeier et al., 2002):
i) a cold part, computed through a zero-temperature EoS, $P_{\rm
cold}=\kappa_{i}\,\rho{{}_{0}}^{\Gamma_{i}}$, in a set of intervals in rest-
mass density $\rho_{0}$ (typically referred as a piecewise-polytropic
representation of a nuclear EoS (Read et al., 2009)); and ii) a thermal
(ideal-gas-like) part to account for shock heating, $P_{\rm th}=\epsilon_{\rm
th}\,(\Gamma_{\rm th}-1)$, which is computed through a thermal index
$\Gamma_{\rm th}$. Here, $\kappa_{i}$ and $\Gamma_{i}$ are the polytropic
constant and index, respectively, and $P_{\rm th}$ and $\epsilon_{\rm th}$ the
thermal pressure and thermal energy density. Moreover, $\Gamma_{\rm th}$ is a
constant ranging between $1$ and $2$ (Constantinou et al., 2015). Above half
nuclear saturation density $\Gamma_{\rm th}$ strongly depends on the nucleon
effective mass (Lim & Holt, 2019). Hence a constant value of $\Gamma_{\rm th}$
may overestimate the thermal pressure by a few orders of magnitude (Raithel et
al., 2021) which, in turn, may induce significant changes in the GW
frequencies (Bauswein et al., 2010; Figura et al., 2021). The second approach
is based on the use of microphysical, finite-temperature EoS tables
constructed using “tabulated” data from observations and nuclear physics
experiments (Oechslin et al., 2007; Bauswein et al., 2010; Sekiguchi et al.,
2011; Fields et al., 2023; Espino et al., 2022; Werneck et al., 2023). This
approach provides a self-consistent method for probing the impact of thermal
effects on the fate of the binary remnant.

In this work, we assess the importance of a self-consistent (or “tabulated”)
treatment of thermal effects in BNS post-merger remnants within a fully
Bayesian framework. We use the numerical-relativity waveforms recently
obtained by Guerra et al. (2023) both as “reference” signals (or injections)
and as recovery templates, considering sources where thermal effects are
modelled both through a simplified “hybrid” approach and a more realistic
“tabulated” one. Our sample of EoSs includes the SLy4, DD2, HShen and LS220
cases. Given reference data consisting on a synthetic “tabulated” signal, we
can obtain Bayesian evidences for both the “tabulated” and “hybrid” waveform
models. With this, we evaluate the level to which the difference between the
two treatments of thermal effects leads to observable effects in the GW
emission. We do this through a novel technique that makes direct use of the
Newman-Penrose $\psi_{4}$ scalar directly outputted by numerical-relativity
simulations (Bustillo et al., 2022a), as opposed to traditional analyses
making use of the GW strain $h(t)$. This prevents the occurrence of systematic
numerical errors that might arise when computing the latter from the former
through a double time integration (Reisswig & Pollney, 2010). We show that,
regardless of the EoS, differences between “tabulated” and “hybrid” treatments
of thermal effects lead to differences in the post-merger GW that are
observable by third-generation detectors at source distances $d_{L}\leq 50$
Mpc (averaged over the source sky-location). Such distances reach $150$ Mpc
for orbital inclinations $\iota<0.8$, with the exception of BNSs modelled with
the DD2 EoS.

On a similar note, Fields et al. (2023) also studied the impact of thermal
effects on the GWs from BNS mergers by changing the effective masses of
neutrons and protons, and thus the specific heat capacity of the binary. They
found that an increased heat capacity results in denser, cooler remnants,
which leaves imprints on the GWs. However, it is not clear if these imprints
are actually due to thermal effects or to changes in the nuclei properties
triggered by changes in the effective mass. In addition, Miravet-Tenés et al.
(2023b) have recently used the same simulations of Guerra et al. (2023) to
study the detectability of thermal effects in post-merger signals using
Bayeswave (Cornish & Littenberg, 2015), a Bayesian data-analysis algorithm
that reconstructs signals injected into noise through a morphology-independent
approach. Miravet-Tenés et al. (2023b) find that differences in the
distribution of the main frequency peaks in the spectra in hybrid and
tabulated models can be resolved in third-generation detectors up to distances
similar to those reported in this work.

##  2 Numerical Setup

Our study is based on the recent numerical-relativity simulations performed by
Guerra et al. (2023). We briefly describe them here, addressing interested
readers to Guerra et al. (2023) for full details. The binaries consist of two
equal-mass, irrotational NSs modeled by finite-temperature, microphysical EoSs
(see Table 1). The initial temperature is fixed to $T=0.01\,\rm MeV$. We use
the tables by Schneider et al. (2017), freely available at
stellarcollapse.org. We choose four different EoS that span a suitable range
of NS central densities, radii and maximum gravitational masses for
irrotational NSs. We also consider a piecewise-polytropic representation of
the cold part of these EoSs using a piecewise regression as in Pilgrim (2021)
with seven polytropic pieces (Read et al., 2009). To this cold part we add an
ideal-gas part with a constant adiabatic index $\Gamma_{\rm th}=1.8$ to
account for thermal effects. Initial data are obtained using LORENE
(Gourgoulhon et al., 2001; Taniguchi & Gourgoulhon, 2002).

EoS | $M$ | $\mathcal{C}$ | $k_{2}$ | $\Lambda$ | $M_{\rm ADM}$ | $J_{\rm ADM}$ | $\Omega$
---|---|---|---|---|---|---|---
SLy4 | 1.28 | 0.13 | 0.086 | 536.09 | 2.54 | 6.63 | 1.77
DD2 | 1.29 | 0.11 | 0.105 | 1100.92 | 2.56 | 6.73 | 1.77
HShen | 1.30 | 0.10 | 0.109 | 1805.63 | 2.58 | 6.82 | 1.78
LS220 | 1.29 | 0.12 | 0.106 | 915.00 | 2.55 | 6.68 | 1.77
Table 1: Summary of the initial properties of the tabulated BNS
configurations: We list the EoS, the gravitational mass $M\,[M_{\odot}]$, the
compactness $\mathcal{C}\equiv M/R_{\rm eq}$, the second Love number $k_{2}$,
and the tidal deformability $\Lambda=(2/3)\kappa_{2}\,\mathcal{C}^{-5}$ for an
isolated NS with the same $M$. Here $R_{\rm eq}$ is the equatorial coordinate
stellar radius. The last three columns report the ADM mass $M_{\rm
ADM}\,[M_{\odot}]$, the ADM angular momentum $J_{\rm ADM}\,[M_{\odot}^{2}]$
and the angular velocity $\Omega\,[\rm krad/s]$. In all cases the NSs have a
rest-mass $M_{0}=1.4M_{\odot}$, and an initial temperature of $0.01\rm MeV$.
The initial binary coordinate separation is $\sim 44.3\,\rm km$. The hybrid
BNS models have similar initial properties (see Guerra et al. (2023)).  Figure
1: Each panel displays the rest-mass density (left) and thermal specific
internal energy (right) of the binary remnant on the equatorial plane at
selected times. Those are shown for the tabulated (top) and hybrid (bottom)
cases of the DD2 EoS (top row) and HShen (bottom row) EoS. Isocontours
correspond to rest-mass densities $\rho_{0}=\\{10^{11},10^{12},10^{14}\\}\,\rm
g\,cm^{-3}$. The boundary of the bulk of the star is displayed as dashed
(contour) lines defined as regions where $\rho_{0}=10^{-3}\,{\rho}_{\rm
0,max}$, where ${\rho}_{\rm 0,max}$ is the initial maximum value of the rest-
mass density.

Evolutions were performed using the publicly available IllinoisGRMHD code
(Werneck et al., 2023; Etienne et al., 2015) embedded on the Einstein Toolkit
infrastructure (Loffler et al., 2012). The code evolves the
Baumgarte–Shapiro–Shibata–Nakamura equations for the spacetime fields
(Baumgarte & Shapiro, 1998; Shibata & Nakamura, 1995) coupled to the puncture
gauge conditions setting the damping coefficient appearing in the shift
condition to $1/M_{\rm ADM}$. The illinoisGRMHD code adopts the Valencia
formalism for the general relativistic hydrodynamics equations (Banyuls et
al., 1997) which are integrated numerically with a state-of-the-art finite-
volume algorithm. The simulations used fourth-order spatial stencils and a
fourth-order Runge-Kutta scheme for time integration with a Courant-
Friedrichs-Lewy factor of 0.5.

##  3 Analysis setup

To assess the observational importance of an accurate implementation of
thermal effects in BNS mergers we perform parameter inference and model
selection on numerically simulated GW signals $h(\theta_{\rm true})$, with
source parameters $\theta_{\rm true}$. These signals correspond to the very
late inspiral, merger and post-merger emission of the BNS models of Table 1.
We inject them in an idealised three-detector network composed of three Cosmic
Explorers (Reitze et al., 2019; Abbott et al., 2017f) placed at the locations
of LIGO Hanford, LIGO Livingston and Virgo. We perform Bayesian parameter
inference on these signals using numerically simulated templates
$h_{\cal{M}}(\theta)$ according to two different emission models ${\cal{M}}$,
respectively including the two alternative implementations of thermal effects
(i.e. hybrid and tabulated).

The posterior Bayesian probability for source parameters $\theta$ according to
a waveform model ${\cal{M}}$ is given by

$p_{\cal{M}}(\theta|\theta_{\rm
true})=\frac{\pi(\theta){\cal{L}}(\theta|h_{{\cal{M}}}(\theta_{\rm
true}))}{{\cal{Z}}_{{\cal{M}}}(\theta|\theta_{\rm true})}.$ (1)

Here, ${\cal{L}}(\theta|h(\theta_{\rm true}))$ denotes the likelihood for the
parameters $\theta$ while $\pi(\theta)$ denotes their prior probability.
Finally, the term ${\cal{Z}}_{\cal{M}}$ denotes the Bayesian evidence for the
model ${\cal{M}}$, given by

${\cal{Z}}_{\cal{M}}=\int\pi(\theta){\cal{L}}(\theta|h_{{\cal{M}}}(\theta_{\rm
true}))d\theta.$ (2)

For two competing waveform models, their relative Bayes Factor is simply given
by

${\cal{B}}^{{\cal{M}}_{1}}_{{\cal{M}}_{2}}=\frac{{\cal{Z}}_{{\cal{M}}_{1}}}{{\cal{Z}}_{{\cal{M}}_{2}}}.$
(3)

It is commonly considered that ${\cal{M}}_{1}$ is strongly preferred with
respect to ${{\cal{M}}_{2}}$ when
$\log{\cal{B}}^{{\cal{M}}_{1}}_{{\cal{M}}_{2}}>5$.

Figure 2: Quadrupolar GW signals of two BNS mergers of our sample, HShen EoS
(left) and DD2 EoS (right). Yellow and red waveforms refer to hybrid and
tabulated models, respectively.

We use the standard frequency-domain likelihood for GW transients (Finn, 1992;
Romano & Cornish, 2017)

$\log{\cal{L}}(\theta|h(\theta_{\rm true}))\propto-\sum_{N}(h(\theta_{\rm
true})-h(\theta)|h(\theta_{\rm true})-h(\theta)),$ (4)

where $N$ runs over the different detectors of our network. As usual, $(a|b)$
represents the inner product (Cutler & Flanagan, 1994)

$(a|b)=4\Re\int_{f_{min}}^{f_{max}}\frac{\tilde{a}(f)\tilde{b}^{*}(f)}{S_{n}(f)}df,$
(5)

where $\tilde{a}(f)$ denotes the Fourier transform of $a(t)$ and ∗ the complex
conjugate. The factor $S_{n}(f)$ is the one-sided power spectral density of
the detector. We use the predicted Cosmic Explorer sensitivity (Abbott et al.,
2017g) with a lower frequency cutoff of $f_{\rm min}=600$ Hz and a sampling
frequency of $16$ kHz so that $f_{\rm max}=8$ kHz.

Our analysis only includes the merger and post-merger emission, together with
the very late inspiral of the process, which is where thermal effects will be
most impactful. We note that this part of the signal provides little
information about the masses and spins of the binary. In addition, given our
usage of numerical relativity simulations, we cannot sample over the
individual masses and spins, but only over the extrinsic source parameters,
namely the source sky-location, orientation, luminosity distance, polarization
angle and time-of-arrival. Nevertheless, it is sensible to assume that masses
and spins would be accurately measured from the minutes-long inspiral signal
(Smith et al., 2021; Branchesi et al., 2023). While the same is true for the
extrinsic parameters, we choose to actually sample over these, so that we
obtain “conservatively pessimistic” results that actually underestimate the
importance of a “tabulated” treatment of thermal effects. We assume standard
prior probabilities for the sky-location, source orientation and polarisation
angle, together with a distance prior uniform in co-moving distance and a
uniform prior on the time-of-arrival, with a width of $0.2\rm\ s$, centered on
the true value. We sample the likelihood across the parameter space using a
custom version of the publicly available software Parallel Bilby (Ashton et
al., 2019; Smith et al., 2020), sampling the parameters with the algorithm
Dynesty (Speagle, 2020).

Finally, we note we do not perform our analysis under the strain formalism
using $h(t)$. Instead, we perform injection and recovery of the Newman-Penrose
scalar $\psi_{4}={d^{2}h(t)}/{dt^{2}}$ directly outputted by numerical
simulations, using the formalism described in Bustillo et al. (2022a). This is
done to avoid potential systematic errors arising during the computation of
the GW strain $h$ from $\psi_{4}$ (Reisswig & Pollney, 2010).

##  4 Results

###  4.1 Merger dynamics and waveforms

The initial data for the BNS reported in Table 1 were evolved in Guerra et al.
(2023) for over $t-t_{\rm mer}\sim 140\,\rm ms$ after merger. Each panel in
Fig. 1 displays two snapshots in the evolution of the rest-mass density (left
half) and the thermal specific energy (right half) for both, tabulated (top
half) and hybrid cases (bottom half). The top two panels correspond to the DD2
EoS and the bottom ones to the HShen EoS. As discussed in Guerra et al.
(2023), we use the thermal specific energy $\epsilon_{\rm th}$ as a proxy for
the temperature. Only at low densities/high temperatures, the definition of
temperature usually employed in BNS merger evolutions based on piecewise
polytropes, $T=(\Gamma_{\rm th}-1)\,\epsilon_{\rm th}$ (see e.g. De Pietri et
al. (2020)) is equivalent to that of microphysical, finite-temperature EOS
evolutions. Fig. 1 shows that for tabulated evolutions the BNS remnant is
hotter near its surface, defined by the density isocontour where
$\rho=10^{-3}\rho_{\rm 0,max}$, with $\rho_{\rm 0,max}$ being the initial
maximum value of the rest-mass density. This induces a slightly more extended
remnant than that of the hybrid evolutions (Guerra et al., 2023).

Fig. 2 displays the dominant $\ell=m=2$ GW quadrupole mode of the late
inspiral, merger and post-merger emission of the simulations corresponding to
the same two EoSs of Fig. 1. The yellow (red) curves indicate hybrid
(tabulated) evolutions. Differences between the two EoSs are noticeable. The
effect of different thermal treatments are clearly more visible in the post-
merger signal emitted by the HShen case (left panel) for which the amplitude,
associated with the oscillations of the remnant, decays significantly faster.
Moreover, small differences between hybrid and tabulated evolutions are
noticeable in the last cycles of the inspiral signals also in the case of the
HShen EoS. This suggests that for this EoS thermal effects might alter the
tidal deformability in a different way for both treatments, possibly leading
to some observational difference.

EoS |  $\log{\cal B}^{\text{T}}_{\text{H}}$ ($d_{L}$ [Mpc], $\rho_{\rm opt}^{\rm net}$)
---|---
| $\iota=0.3$ | $\iota=0.8$ | $\iota=\pi/2$
SLy | 10.1 (100, 15.3) | 22.1 (50, 22.3) | 38.6 (10, 27.9)
LS220 | 9.7 (100, 14.3) | 22.1 (50, 20.9) | 33.7 (10, 26.1)
HShen | 16.2 (100, 15.1 ) | 37.0 (50, 22.1 ) | 65.2 (10, 27.4 )
DD2 | 5.6 (100, 15.3) | 11.7 (50, 22.4) | 14.3 (10, 27.5)
Table 2: Impact of thermal-effects implementation in gravitational-wave
observations We show the natural log Bayes Factors between tabulated and
hybrid models, obtained when these are used to recover a signal from a BNS
post-merger using the tabulated EoS. We show results for source inclinations
of $\iota=\\{0.3,0.8,\pi/2\\}$ and respective distances of
$d_{L}=\\{100,50,10\\}$Mpc. We assume a three-detector network composed by
three Cosmic Explorers.

###  4.2 Target binary neutron star sources

We choose target sources corresponding to four equal-mass BNSs with total mass
$M\approx 2.56M_{\odot}$ (see Table 1) characterized by four different EoSs.
For each source, we consider three orbital inclinations, namely $\iota=0.3$,
$\iota=0.8$ and $\iota=\pi/2$. Varying the orbital inclination changes the
contribution of higher-order harmonics to the observed waveform, which may
influence the impact of the thermal effects. While higher-order harmonics are
highly suppressed for equal-mass systems, it has been shown that these can be
triggered by the so-called “one-arm” spiral instability (East et al., 2016;
Radice et al., 2016; Lehner et al., 2016), facilitating the estimation of the
source orientation (Bustillo et al., 2021b). In particular, our waveforms
include the dominant quadrupole modes, $(\ell,m)=(2,\pm 2)$, together with the
$(\ell,m)=(2,\pm 1),(3,\pm 3),(3,\pm 2)$ modes.

Ideally, we would perform a wide injection campaign considering sources placed
at a large variety of sky-locations, distances and source inclinations.
However, the high computational cost of our parameter estimation runs, which
make use of numerical relativity, prevents this. For this reason, we only use
one fiducial sky-location for our runs. Moreover, the luminosity distance and
sky-location do mainly control the signal loudness, without altering the mode
content of the signal. Given this, it is expected that the ratio of the
maximum signal-to-noise ratios (SNRs) $\rho^{\rm net}$ across the detector
network (and therefore the ratio of the respective maximum likelihoods)
obtained by the tabulated and hybrid models, should show a weak dependence on
these parameters. Since the likelihood, which is what controls the Bayes
Factor ${\cal{B}}$, goes as $(\rho^{\rm net})^{2}/2$, we expect the relative
Bayes Factor for the hybrid and tabulated models to go as $(\rho_{\rm
opt})^{2}/2$. Next, the Bayes Factor for each analysis roughly goes as
$\log{\cal{B}}\simeq\log{\cal{L}}_{\rm max}-{\cal{C}}$, where ${\cal{C}}$
accounts for the Occam penalty paid by the model. Finally, $\log{\cal{L}}_{\rm
max}\simeq{\cal{M}}\times(\rho^{\rm net}_{\rm opt})^{2}/2$, with ${\cal{M}}$
denoting the match between our injection and the best-fitting template, and
$\rho^{\rm net}_{\rm opt}$ indicating the SNR of the injection, equal to the
maximal (or optimal) SNR that any analysis can recover (Lange et al., 2017).

In this situation, for each inclination we perform parameter inference for
three fiducial combinations of sky-location and distances yielding reasonably
different $\rho_{\rm opt}$. With this, we perform a linear fit
$\log{\cal{B}}=\alpha\rho_{\rm max}^{2}+\beta$. Next, for each of our target
sources, we compute the average SNR over sky-locations at a fiducial distance
$d_{L}^{\rm ref}$. Finally, using that the SNR is inversely proportional to
the distance, we compute the distance $d_{L}^{\rm det}$ at which the averaged
relative Bayes Factor between the tabulated and hybrid models satisfies
$\log{\cal{B}}^{\rm T}_{\rm H}=5$, which we will call “detection” distance.

###  4.3 Model selection

Table 2 shows $\log{\cal B}^{\rm T}_{\rm H}$ for selected fiducial runs
corresponding to each of our selected source inclinations. The parentheses
show the luminosity distance chosen for the injection and the corresponding
optimal SNR across the detector network. We find that the HShen EoS is most
impacted by the choice of implementation of thermal effects, for all source
inclinations. This is not a trivial result, as the impact of thermal-effects
treatment in the waveform for a given EoS may strongly depend on the source
inclination. This may cause some EoS to produce more observable differences
for different inclinations. This is observed for the two next cases that are
most impacted by the treatment of thermal-effects, namely the SLy and LS220
EoSs. Both cases return similar values of $\log{\cal B}^{\rm T}_{\rm H}$ for
the two lowest inclinations. However, the SLy EoS is most affected for edge-on
systems. Finally, the emission from BNS mergers with DD2 EoS is the least
influenced.

Figure 3: Impact of thermal effects implementation averaged over source sky-
location and azimuthal angle. Average source distance for which the difference
between the log-evidences obtained by the “tabulated” and “hybrid” models
reaches $\log{\cal{B}}^{\rm T}_{\rm H}=5$, when recoverying a signal modeled
through the “tabulated” implementation of thermal effects (solid lines).
Averaging is performed over source sky-locations and observer’s azimuthal
angle. Distances at which $\log{\cal{B}}^{\rm T}_{\rm H}=8$ are shown in
dashed lines. We show results for and different EoSs and varying source
inclinations.

Figure 3 shows the average “detection” distance $d_{L}^{\rm det}$, as a
function of EoS and source inclination, obtained after performing the fit
described in the previous section. We have checked empirically that typical
deviations from such linear fit are of order $\Delta\log{\cal{B}}\simeq 1$. To
add conservative results we also show lines corresponding to $\log{\cal
B}^{\rm T}_{\rm H}=8$. We find that for all EoS, the signal features coming
from the “tabulated” implementation of thermal effects can be detected at
distances $d_{L}\leq 50$ Mpc regardless of the source inclination, and at
distances $d_{L}=100$ Mpc for inclinations $\iota<0.8$. In other words, the
correct analysis of a source like GW170817 would require of the more realistic
tabulated modeling of its EoS, were the source consistent with one of the EoS
of our sample and assuming that other effects our simulations do not yet
incorporate, like magnetic fields, bulk and shear viscosity, and neutrinos, do
not play a major role111It should be noted that the accurate numerical
modelling of the actual GW spectrum of BNS post-merger remnants is still a
matter of debate. In particular, the development of MHD-driven turbulence in
the remnant, triggered by the magneto-rotational instability, is likely to
affect the evolution of the system, influencing the emission of high-frequency
GWs (see e.g. Shibata & Kiuchi (2017); Ruiz et al. (2021b); Chabanov &
Rezzolla (2023) for recent studies).. Moreover, the same is true even for weak
edge-on sources, for averaged distances of order $50$ Mpc.

To close this section, we note that the use of a simplified hybrid thermal
treatment leads, in certain cases, to significant parameter biases in the
sampled parameters, specially the luminosity distance. The biases are due to
the different signal amplitudes predicted by tabulated and hybrid models, as
it is clear in Fig. 2, specially in the case of the HShen EoS. Nevertheless,
we stress again that these parameters are expected to be measured with great
accuracy from the minute-duration inspiral signal.

##  5 Conclusions

Using a combination of a Bayesian framework and numerical-relativity
simulations, we have shown the dramatic importance of an accurate treatment of
thermal effects in the post-merger GW emission of BNS mergers. We have found
that, for all EoSs and source inclinations considered, self-consistent
“tabulated” implementations lead to modifications of the post-merger GW
emission with respect to simplified “hybrid” approaches that are observable
with a detector network formed by three Cosmic Explorers, for sources at
distances $d_{L}\leq 50$ Mpc. In particular, for inclinations $\iota\leq 0.8$,
consistent with GW170817, the differences are visible for distances $d_{L}\leq
150$ Mpc, with the exception of the DD2 EoS. Out of the four EoSs considered
and within the limitations of our simulations, the post-merger GW signal of
the HShen EoS is the most influenced by a “hybrid” implementation of thermal
effects.

We note that our results rather underestimate the importance of a self-
consistent implementation of thermal effects, as we perform inference on
parameters that are expected to be accurately constrained from the long
inspiral signal, like sky-location or source orientation. This allows hybrid
models to exploit parameter degeneracies to achieve higher Bayesian evidences
than they would do otherwise, if such parameters were accurately constrained
through the inspiral signal.

Our work can be regarded as a “proof of principle” of the application of
parameter inference to BNS merger remnants using numerical-relativity
waveforms, which is obviously limited by the fact that we do not sample over
the actual EoS. A natural extension of our work would be to compare injections
including thermal effects for a given EoS with a large set (ideally
continuous) of templates covering a wide range of EoSs (as somewhat done in
Bustillo et al. (2021a, 2022b) for the case of collisions of a different
matter source, Proca-stars222Note that, in the Proca-star case, the
simulations do not span a continuous range of EoSs, as for neutron stars.
Instead, they span a very dense grid in the individual frequencies of the
Proca fields. ). However, the high cost of our numerical simulations, together
with the intrinsic discreteness and low number of simulations including
thermal effects prevents such extension at present. Nevertheless, our study
proves that waveforms sourced by the EoSs of our sample are sufficiently
different so as to allow the identification of the underlying EoS for the
range of signal loudness we consider.

##  6 Acknowledgements

It is a pleasure to thank Zachariah Etienne and Leonardo Werneck for useful
discussions. This work was supported by a fellowship from “la Caixa”
Foundation (ID100010434), the European Union’s Horizon2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
847648 (fellowship code LCF/BQ/PI20/11760016), the Generalitat Valenciana
through the grants CIDEGENT/2021/046 and Prometeo CIPROM/2022/49, and the
Spanish Agencia Estatal de Investigación through the grants PRE2019-087617,
PID2020-118635GB-I00 and PID2021-125485NB-C21 funded by
MCIN/AEI/10.13039/501100011033 and ERDF A way of making Europe. Further
support has been provided by the EU’s Horizon 2020 Research and Innovation
(RISE) programme H2020-MSCA-RISE-2017 (FunFiCO-777740) and by the EU Staff
Exchange (SE) programme HORIZON-MSCA-2021-SE-01 (NewFunFiCO-101086251). VVO
acknowledges support from the Ministry of Science, Innovation and Universities
(MICIU) of Spain (FPI grant PRE2018-085436 and ) The authors acknowledge
computational resources provided by the LIGO Laboratory and supported by
National Science Foundation Grants PHY-0757058 and PHY0823459, the support of
the NSF CIT cluster for the provision of computational resources for our
parameter inference runs, and the computational resources and technical
support of the Spanish Supercomputing Network through the use of MareNostrum
at the Barcelona Supercomputing Center (AECT-2023-1-0006) where the BNS merger
simulations were performed. This material is based upon work supported by
NSF’s LIGO Laboratory which is a major facility fully funded by the National
Science Foundation. This document has LIGO DCC number LIGO-P2300344.

## References

  * Abbott et al. (2017a) Abbott, B. P., et al. 2017a, Astrophys. J., 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2017b) —. 2017b, Astrophys. J. Lett., 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2017c) —. 2017c, Astrophys. J. Lett., 850, L39, doi: 10.3847/2041-8213/aa9478
  * Abbott et al. (2017d) —. 2017d, Nature, 551, 85, doi: 10.1038/nature24471
  * Abbott et al. (2017e) —. 2017e, Astrophys. J. Lett., 851, L16, doi: 10.3847/2041-8213/aa9a35
  * Abbott et al. (2017f) —. 2017f, Classical and Quantum Gravity, 34, 044001, doi: 10.1088/1361-6382/aa51f4
  * Abbott et al. (2017g) —. 2017g, Classical and Quantum Gravity, 34, 044001, doi: 10.1088/1361-6382/aa51f4
  * Abbott et al. (2021) Abbott, R., et al. 2021.  https://arxiv.org/abs/2111.03606
  * Ashton et al. (2019) Ashton, G., et al. 2019, Astrophys. J. Suppl., 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Baiotti & Rezzolla (2017) Baiotti, L., & Rezzolla, L. 2017, Rept. Prog. Phys., 80, 096901, doi: 10.1088/1361-6633/aa67bb
  * Banyuls et al. (1997) Banyuls, F., Font, J. A., Ibáñez, J. M., Martí, J. M., & Miralles, J. A. 1997, ApJ, 476, 221, doi: 10.1086/303604
  * Baumgarte & Shapiro (1998) Baumgarte, T. W., & Shapiro, S. L. 1998, Phys. Rev. D, 59, 024007, doi: 10.1103/PhysRevD.59.024007
  * Bauswein et al. (2010) Bauswein, A., Janka, H. T., & Oechslin, R. 2010, Phys. Rev. D, 82, 084043, doi: 10.1103/PhysRevD.82.084043
  * Bauswein et al. (2017) Bauswein, A., Just, O., Janka, H.-T., & Stergioulas, N. 2017, Astrophys. J. Lett., 850, L34, doi: 10.3847/2041-8213/aa9994
  * Bauswein et al. (2020) Bauswein, A., Blacker, S., Vijayan, V., et al. 2020, Phys. Rev. Lett., 125, 141103, doi: 10.1103/PhysRevLett.125.141103
  * Branchesi et al. (2023) Branchesi, M., et al. 2023, JCAP, 07, 068, doi: 10.1088/1475-7516/2023/07/068
  * Bustillo et al. (2021a) Bustillo, J. C., Lasky, P. D., & Thrane, E. 2021a, Physical Review D, 103, doi: 10.1103/physrevd.103.024041
  * Bustillo et al. (2021b) Bustillo, J. C., Leong, S. H. W., Dietrich, T., & Lasky, P. D. 2021b, The Astrophysical Journal Letters, 912, L10, doi: 10.3847/2041-8213/abf502
  * Bustillo et al. (2022a) Bustillo, J. C., Wong, I. C. F., Sanchis-Gual, N., et al. 2022a, doi: https://doi.org/10.48550/arXiv.2205.15029
  * Bustillo et al. (2022b) Bustillo, J. C., Sanchis-Gual, N., Leong, S. H. W., et al. 2022b, doi: https://doi.org/10.48550/arXiv.2206.02551
  * Chabanov & Rezzolla (2023) Chabanov, M., & Rezzolla, L. 2023.  https://arxiv.org/abs/2307.10464
  * Chatziioannou et al. (2017) Chatziioannou, K., Clark, J. A., Bauswein, A., et al. 2017, Phys. Rev. D, 96, 124035, doi: 10.1103/PhysRevD.96.124035
  * Clark et al. (2016) Clark, J. A., Bauswein, A., Stergioulas, N., & Shoemaker, D. 2016, Classical and Quantum Gravity, 33, 085003, doi: 10.1088/0264-9381/33/8/085003
  * Constantinou et al. (2015) Constantinou, C., Muccioli, B., Prakash, M., & Lattimer, J. M. 2015, Phys. Rev. C, 92, 025801, doi: 10.1103/PhysRevC.92.025801
  * Cornish & Littenberg (2015) Cornish, N. J., & Littenberg, T. B. 2015, Classical and Quantum Gravity, 32, 135012, doi: 10.1088/0264-9381/32/13/135012
  * Cowperthwaite et al. (2017) Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJ, 848, L17, doi: 10.3847/2041-8213/aa8fc7
  * Cutler & Flanagan (1994) Cutler, C., & Flanagan, É. E. 1994, Physical Review D, 49, 2658, doi: 10.1103/physrevd.49.2658
  * De Pietri et al. (2020) De Pietri, R., Feo, A., Font, J. A., et al. 2020, Phys. Rev. D, 101, 064052, doi: 10.1103/PhysRevD.101.064052
  * Dietrich et al. (2020) Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Science, 370, 1450, doi: 10.1126/science.abb4317
  * Dimmelmeier et al. (2002) Dimmelmeier, H., Font, J. A., & Müller, E. 2002, A&A, 388, 917, doi: 10.1051/0004-6361:20020563
  * East et al. (2016) East, W. E., Paschalidis, V., & Pretorius, F. 2016, Classical and Quantum Gravity, 33, 244004, doi: 10.1088/0264-9381/33/24/244004
  * Espino et al. (2022) Espino, P. L., Bozzola, G., & Paschalidis, V. 2022.  https://arxiv.org/abs/2210.13481
  * Etienne et al. (2015) Etienne, Z. B., Paschalidis, V., Haas, R., Mösta, P., & Shapiro, S. L. 2015, Class. Quant. Grav., 32, 175009, doi: 10.1088/0264-9381/32/17/175009
  * Ezquiaga & Zumalacárregui (2017) Ezquiaga, J. M., & Zumalacárregui, M. 2017, Phys. Rev. Lett., 119, 251304, doi: 10.1103/PhysRevLett.119.251304
  * Fields et al. (2023) Fields, J., Prakash, A., Breschi, M., et al. 2023, Astrophys. J. Lett., 952, L36, doi: 10.3847/2041-8213/ace5b2
  * Figura et al. (2021) Figura, A., Li, F., Lu, J.-J., et al. 2021, Phys. Rev. D, 103, 083012, doi: 10.1103/PhysRevD.103.083012
  * Finn (1992) Finn, L. S. 1992, Physical Review D, 46, 5236, doi: 10.1103/physrevd.46.5236
  * Gourgoulhon et al. (2001) Gourgoulhon, E., Grandclement, P., Taniguchi, K., Marck, J.-A., & Bonazzola, S. 2001, Phys. Rev. D, 63, 064029, doi: 10.1103/PhysRevD.63.064029
  * Guerra et al. (2023) Guerra, D., et al. 2023, in preparation
  * Janka et al. (1993) Janka, H. T., Zwerger, T., & Moenchmeyer, R. 1993, A&A, 268, 360
  * Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Nature, 551, 80, doi: 10.1038/nature24453
  * Lange et al. (2017) Lange, J., et al. 2017, Phys. Rev., D96, 104041, doi: 10.1103/PhysRevD.96.104041
  * Lehner et al. (2016) Lehner, L., Liebling, S. L., Palenzuela, C., & Motl, P. M. 2016, Physical Review D, 94, doi: 10.1103/physrevd.94.043003
  * Lim & Holt (2019) Lim, Y., & Holt, J. W. 2019.  https://arxiv.org/abs/1909.09089
  * Loffler et al. (2012) Loffler, F., et al. 2012, Class. Quant. Grav., 29, 115001, doi: 10.1088/0264-9381/29/11/115001
  * Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, Astrophys. J. Lett., 850, L19, doi: 10.3847/2041-8213/aa991c
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650, doi: 10.1111/j.1365-2966.2010.16864.x
  * Miravet-Tenés et al. (2023a) Miravet-Tenés, M., Castillo, F. L., De Pietri, R., Cerdá-Durán, P., & Font, J. A. 2023a, Phys. Rev. D, 107, 103053, doi: 10.1103/PhysRevD.107.103053
  * Miravet-Tenés et al. (2023b) Miravet-Tenés, M., Ruiz, M., Guerra, D., Cerdá-Durán, P., & Font, J. A. 2023b, in preparation
  * Most et al. (2018) Most, E. R., Weih, L. R., Rezzolla, L., & Schaffner-Bielich, J. 2018, Phys. Rev. Lett., 120, 261103, doi: 10.1103/PhysRevLett.120.261103
  * Oechslin et al. (2007) Oechslin, R., Janka, H. T., & Marek, A. 2007, Astron. Astrophys., 467, 395, doi: 10.1051/0004-6361:20066682
  * Pardo et al. (2018) Pardo, K., Fishbach, M., Holz, D. E., & Spergel, D. N. 2018, JCAP, 07, 048, doi: 10.1088/1475-7516/2018/07/048
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67, doi: 10.1038/nature24298
  * Pilgrim (2021) Pilgrim, C. 2021, The Journal of Open Source Software, 6, 3859, doi: 10.21105/joss.03859
  * Radice et al. (2016) Radice, D., Bernuzzi, S., & Ott, C. D. 2016, Physical Review D, 94, doi: 10.1103/physrevd.94.064011
  * Raithel et al. (2021) Raithel, C., Paschalidis, V., & Özel, F. 2021, Phys. Rev. D, 104, 063016, doi: 10.1103/PhysRevD.104.063016
  * Read et al. (2009) Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Phys. Rev., D79, 124032, doi: 10.1103/PhysRevD.79.124032
  * Reisswig & Pollney (2010) Reisswig, C., & Pollney, D. 2010, doi: 10.1088/0264-9381/28/19/195015
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019
  * Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, Astrophys. J. Lett., 852, L25, doi: 10.3847/2041-8213/aaa401
  * Romano & Cornish (2017) Romano, J. D., & Cornish, N. J. 2017, Living Reviews in Relativity, 20, doi: 10.1007/s41114-017-0004-1
  * Ruiz et al. (2018) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, Phys. Rev. D, 97, 021501, doi: 10.1103/PhysRevD.97.021501
  * Ruiz et al. (2021a) —. 2021a, Front. Astron. Space Sci., 8, 39, doi: 10.3389/fspas.2021.656907
  * Ruiz et al. (2021b) Ruiz, M., Tsokaros, A., & Shapiro, S. L. 2021b, Phys. Rev. D, 104, 124049, doi: 10.1103/PhysRevD.104.124049
  * Sarin & Lasky (2021) Sarin, N., & Lasky, P. D. 2021, Gen. Rel. Grav., 53, 59, doi: 10.1007/s10714-021-02831-1
  * Schneider et al. (2017) Schneider, A. S., Roberts, L. F., & Ott, C. D. 2017, Phys. Rev. C, 96, 065802, doi: 10.1103/PhysRevC.96.065802
  * Sekiguchi et al. (2011) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2011, Phys. Rev. Lett., 107, 051102, doi: 10.1103/PhysRevLett.107.051102
  * Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., et al. 2017, Phys. Rev. D, 96, 123012, doi: 10.1103/PhysRevD.96.123012
  * Shibata & Hotokezaka (2019) Shibata, M., & Hotokezaka, K. 2019, Ann. Rev. Nucl. Part. Sci., 69, 41, doi: 10.1146/annurev-nucl-101918-023625
  * Shibata & Kiuchi (2017) Shibata, M., & Kiuchi, K. 2017, Phys. Rev. D, 95, 123003, doi: 10.1103/PhysRevD.95.123003
  * Shibata & Nakamura (1995) Shibata, M., & Nakamura, T. 1995, Phys. Rev. D, 52, 5428, doi: 10.1103/PhysRevD.52.5428
  * Smith et al. (2021) Smith, R., Borhanian, S., Sathyaprakash, B., et al. 2021, Physical Review Letters, 127, doi: 10.1103/physrevlett.127.081102
  * Smith et al. (2020) Smith, R. J. E., Ashton, G., Vajpeyi, A., & Talbot, C. 2020, Monthly Notices of the Royal Astronomical Society, 498, 4492, doi: 10.1093/mnras/staa2483
  * Speagle (2020) Speagle, J. S. 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132, doi: 10.1093/mnras/staa278
  * Taniguchi & Gourgoulhon (2002) Taniguchi, K., & Gourgoulhon, E. 2002, Phys. Rev. D, 66, 104019, doi: 10.1103/PhysRevD.66.104019
  * Werneck et al. (2023) Werneck, L. R., et al. 2023, Phys. Rev. D, 107, 044037, doi: 10.1103/PhysRevD.107.044037


# Bessel-beam direct-write of the etch-mask in a nano-film of alumina for
high-efficiency Si solar cells

Tomas Katkus1, Soon Hock Ng1, Haoran Mu1, Nguyen Hoai An Le1, Dominyka
Stonytė2, Zahra Khajehsaeidimahabadi1, Gediminas Seniutinas1, Justas
Baltrukonis3, Orestas Ulčinas3, Mindaugas Mikutis3, Vytautas Sabonis3,
Yoshiaki Nishijima4, Michael Rienäcker5, Jan Krügener6, Robby Peibst5,6,
Sajeev John7, Saulius Juodkazis1,2,8 1 Optical Sciences Centre, ARC Training
Centre in Surface Engineering for Advanced Materials (SEAM), Swinburne
University of Technology, Hawthorn, Victoria 3122, Australia  2 Laser Research
Center, Physics Faculty, Vilnius University, Saulėtekio Ave. 10, 10223
Vilnius, Lithuania  3 Altechna R&D, Mokslininku̧ str. 6A, 08412 Vilnius,
Lithuania  4Department of Electrical and Computer Engineering, Graduate School
of Engineering, Yokohama National University, 79-5 Tokiwadai, Hodogaya-ku,
Yokohama, 240-8501, Japan  5 Institut für Solarenergieforschung Hameln (ISFH),
Am Ohrberg 1, 31860 Emmerthal, Germany  6 Institute of Electronic Materials
and Devices, Leibniz Universität Hannover, Schneiderberg 32, 30167 Hannover,
Germany  7Department of Physics, University of Toronto, 60 St. George Street,
Toronto, ON, M5S 1A7, Canada  8 WRH Program International Research Frontiers
Initiative (IRFI) Tokyo Institute of Technology, Nagatsuta-cho, Midori-ku,
Yokohama, Kanagawa 226-8503 Japan

###### Abstract

Large surface area applications such as high-efficiency $>26\%$ solar cells
require surface patterning with 1-10 $\mu$m periodic patterns at high fidelity
over 1-10 cm2 areas (before up scaling to 1 m2) to perform at, or exceed, the
Lambertian (ray optics) limit of light trapping. Here we show a pathway to
high-resolution sub-1 $\mu$m etch mask patterning by ablation using direct
femtosecond laser writing performed at room conditions (without the need for a
vacuum-based lithography approach). A Bessel beam was used to alleviate the
required high surface tracking tolerance for ablation of 0.3-0.8 $\mu$m
diameter holes in $\sim 40$ nm alumina Al2O3-mask at high writing speed, 7.5
cm/s; a patterning rate 1 cm2 per 20 min. Plasma etching protocol was
optimised for a zero-mesa formation of photonic crystal (PhC) trapping
structures and smooth surfaces at the nanoscale level. Scaling up in area and
throughput of the demonstrated approach is outlined.

Si solar cell, Lambertian limit, high-efficiency solar-to-electrical energy
conversion, Bessel beam

##  I Introduction

Si solar cells are the most developed practical solution for harvesting solar
energy and are widely available at a $>23\%$ efficiency of solar-to-electrical
power conversion at the commercial photo-voltaic (PV) module level for a
household and industrial usage as install-ready $\sim 2$ m2 panels [1, 2, 3].
The technology of Si solar cells is well-matured and produces energy at a cost
lower than that traditionally available from national power grids distributing
electricity from power stations (coal, gas, nuclear) [4]. Recently, the world
record performance has been improved from the previously achieved efficiency
of Si solar cells 26.81% [5] and stands now at 27.1% measured by the National
Renewable Energy Laboratory (NREL) [6]. An improvement by 1% in the efficiency
took $\sim 25$ years from M. Green’s established 25% benchmark for the
monocrystalline Si demonstrated in 1999 [7]. None of the modern solar cells
are performing at the Lambertian (ray optics) limit of light trapping
$4n^{2}$, where $n$ is the refractive index of solar cell material [8]. The
practical limit is predicted at 29.1%, which is lower than the theoretical
Shockley–Queisser 32% boundary for Si at one Sun [7].

Random micro-textured surfaces are used to reduce reflectivity and enhance
light trapping, however, they fundamentally cannot perform above the
Lambertian limit, which requires harnessing the interference, slow light, and
normal-to-parallel redirection of incident light. By harnessing those wave
properties provide an enhanced absorbance at strongly localized resonances,
especially at the energies close to the bandgap [9]. The first experimental
demonstration of light absorption above the Lambertian limit in Si slab was
made using PhC structures with etch masks patterned by electron beam
lithography (EBL) with subsequent wet etching in KOH and plasma [10]. The
possibility to define etch mask by stepper projection lithography and direct
laser write (DLW), which both have the capability of larger area patterning,
was demonstrated for Si wafer [11]. The same PhC pattern, which performed the
above Lambertian light trapping, was made on actual Si solar cells and
compared to a neighboring central cell of the same geometry and structure on
the same wafer but without PhC patterns (7 solar cells in total per wafer)
[12]. This first demonstration of PhC light trapping on a Si solar cell showed
23.1% efficiency on the PhC textured region (independently verified by ISFH
CalTeC) as compared to 21.1% for the neighboring central cell on the same
wafer. Another active trend in solar cells is making them thin, flexible and
low-weight [13] including use of 2D materials [14]. PhC light trapping can be
used to achieve 95% absorption of light within 10 $\mu$m thick Si solar cell
to expand their application potential [15]. Flexible Si solar cells of 54
$\mu$m thickness showed the highest power-to-weight ratio of 1.9 W/g for $\sim
274$ cm2 area [16]; the light-to-electrical power converison was 26.06%.

Figure 1:  SEM images of 38-nm-thick ($n+i\kappa=1.595+i\times 0.01054$) Al2O3
mask ablated by single pulses 515 nm/170 fs (Carbide, Light Conversion) using
an axicon lens, which formed the Bessel beam. Pulse energies were
$E_{p}\approx 32,44,56$ nJ (marked on images). Photo image of a Si wafer with
$\sim$40-nm-thick mask ablated in one and two pulses per site modes. The
circle marker of 629 nm corresponds to a barely recognizable modification at
the smallest pulse energy. This diameter is smaller than the
$1/e^{2}$-intensity cross-section of the central core of the Bessel beam which
was $\sim 880~{}$nm (the first minimum of the Bessel intensity envelope).
Figure 2:  Optical microscopy observation of etch time progression of Photonic
Crystal (PhC) light trapping surface on p-type Si (783 wafer); microscope
Nikon Optiphot-pol. with an objective lens of numerical aperture $NA=0.9$; see
top-markers in (a-f) for the 5-18 min span). Reactive ion etching (RIE; Samco
RIE-101iPH) at high etch rate of isotropic mode at 200 W of Inductively
Coupled Plasma (ICP) power and at 0 W Bias (directional etch). Chemistry: 50/5
cm3/min (or sccm) of SF6/O2; process pressure 2 Pa, He cooling pressure 1 Pa.
The inset in (c) is a Scanning Electron Microscopy (SEM) image of etched PhC
at 50/0 sccm of SF6/O2 after 15 min. It illustrates the mesa-structures and
Al2O3 mask adhered at the corner positions of the square pits (a rainbow-
marker). When those corner-structures become smaller than $\sim 1~{}\mu$m in
size (close to the diffraction limit of optical observation $0.61\lambda/NA$
at central visible wavelength $\lambda=530$ nm-green), they appear colorful
(see inset in (d)). This condition indicates an optimal time for the smallest
mesa segments. Longer etch cause failure of the alumina mask shown in (e).
Removal of 38-nm-thick alumina mask was made in ultra-sonic bath in millipore
water. Mask writing was made with the Bessel beam.  Figure 3:  SEM images of
etch cross sections at different Bias powers: 0 W (a), 2 W (b), 10 W (c) and
50 W (top-view in (d)). Etch time is shown at the bottom of panels; when
repeated etching was applied, times were shown with a summation marker “+”
(b). RIE was carried out at 2 Pa process pressure, 200 W ICP power at 50/0
sccm SF6/O2 gas mixture. For comparison of evolution of etched pattern oval
and right-triangular markers of the same size are shown in (a-c). Thickness of
Al2O3 mask was 28.4 nm; mask writing was with Bessel beam. Residue of oxidised
debris increased with larger bias powers and after repeated etching runs.
Figure 4:  Evolution of under-etching at decreasing ICP power from 100 W to 30
W (the final optimised condition) at 0 W Bias; SF6/O2 50/0 sccm, pressure of 2
Pa. The under-etching removes the post regions. The minimum size of the posts
also corresponded at the largest depth of the etched pattern (appeared
colorful under optical observation; Fig. A6). The timing of etching is the
main parameter for the optimal conditions of PhC formation. Thickness of Al2O3
mask was 38 nm; mask writing was with Bessel beam. The angle of $54.7^{\circ}$
between top surface (100) plane and (111) plane is shown by dashed-inclined
lines.  Figure 5:  PhC light trapping for Polysilicon on Oxide Interdigitated
Back Contacted - POLO2-IBC cell. The thickness of Si solar cell 190 $\mu$m.
(a) The mask (38-nm-thick Al2O3) was patterned by fs-laser direct write with
Bessel beam 515 nm/167 fs, pulse energy $E_{p}=43.8$ nJ over $2.2\times 2.2$
cm2 areas on the front-side above of an integrated Si cells with back-side
contacts. Plasma etch and PhC pattern evolution were observed at different
paused times during etching. The optimised protocol was used: 1 Pa process
pressure, 30 W ICP power, SF6/O2 50/0 ccm flow rate. (b) Center-image is the
final PhC trapping surface after mask removal in an IPA ultrasonic bath.
Insets show SEM images of the final PhC light trapping surface at different
locations. There was a slight difference in terms of the under-etch conditions
across the wafer recognizable in optical (a) and SEM (b) images. The right-
bottom inset image shows the back-side contacts.  Figure 6:  (a) Concept for a
better mechanical support of the nano-film mask over large areas for a deeper
etch of PhC patterns: laser inscribed mask in $\sim 30$-nm-thick Al2O3 with
periodic micro-posts (not ablated sites). (b) Finite difference time domain
(FDTD) calculations of absorption and scattering cross sections for the
pattern in (a) and same geometry only with four pyramidal pits. Periodic
boundary conditions are applied for the xy-plane (lateral dimension) and
perfectly matching layers in the longitudinal direction (along light
propagation). The thickness of the cell is 3 $\mu$m; Si permittivity
$(n+i\kappa)^{2}$ is from the Lumerical, Ansys database. Geometrical cross-
section (footprint) of the pyramid (or inverted pyramid) is $9.61\times
10^{6}$ nm2.  Figure 7:  Light intensity enhancement at two cross-sections
(monitors) at difference wavelengths simulated by FDTD with Periodic boundary
conditions applied for the lateral dimension and perfectly matching layers in
the longitudinal direction (along light propagation); same geometry as in Fig.
6. The thickness of the cell is 3 $\mu$m; Si permittivity $(n+i\kappa)^{2}$ is
from the Lumerical, Ansys database. Incident light intensity $|E|^{2}=1$.

Here we demonstrate the fabrication of PhC pattern for light trapping on the
face surface of actual Si solar cells. Dielectric Al2O3 mask of tens-of-nm
thickness was laser ablated by ultra-short laser pulses (515 nm/167 fs) using
an axicon lens for the formation of the Bessel beam [17]. Ablated holes in the
mask had a diameter smaller than the diffraction limit at the focusing
conditions used due to the threshold effect, i.e., the ablation threshold of
the mask was exceeded only in the very center of the focal spot. The dry
plasma etching protocol was optimized for a nano-smooth surface of PhC pattern
with minimal ridges between inverted pyramidal pits of PhC, without apparent
debris and oxidised nanoscale regions. Optimization of plasma etching over a
wide parameter space of chemistry used, process pressure, and flow rates of
gasses is presented. The PhC light trapping structures were fabricated on
Polysilicon on Oxide Interdigitated Back Contacted - POLO2-IBC cell; such
cells have the potential to perform at $>26\%$ efficiency of solar-to-
electrical power conversion [18].

##  II Methods and materials

The thickness of the electron (e)-beam (JKLesker AXXIS) coated Al2O3 was
determined by ellipsometry (Woolam, VASE). Masks of 30 nm and 40 nm were
investigated for optimisation of the laser ablation and plasma etch protocols.
The refractive index $n+i\kappa\approx 1.595+i0.011$ of Al2O3 was determined
by ellipsometry. Thinner 20-nm-masks showed up to a few percent reduction in
$\kappa$ while $n$ was proportionally larger.

Si wafers (783) were p-type (Boron doped) with resistivity
$1-10~{}\Omega\cdot$cm and were close to that of the actual POLO2-IBC p-type
solar cells. The wafer’s cut edge (see inset in Fig. 1) was used to align the
PhC mask pattern. The alignment is critical for plasma etching while low ICP
power is used and the etching result resembles structures, similar to those
obtained by wet etching in KOH with a distinct pattern of inverted pyramids.
Imperfect mask alignment causes defects in mesa ridges and corner posts, which
should be minimized for the best-performing light trapping.

The axicon used in this study was designed for $\lambda=515$ nm wavelength and
it was fs-laser inscribed inside the bulk of silica glass. The form-
birefringence of nano-gratings was used to engineer the refractive index of
concentric grating for the required phase retardance. The cone angle of the
flat axicon corresponded to 178∘. A $f=400$ mm lens combined with $\times 20$
magnification, $NA=0.4$ objective lens was used to scale down the Bessel beam
using an optical 4F relay scheme. The first minimum of the Bessel intensity
profile is at $d_{min}^{B}=0.383\lambda/\sin\gamma$, where $\lambda=515$ nm is
the laser wavelength, $\gamma=18.46^{\circ}$ is the cone angle with respect to
the optical axis along beam propagation;
$d_{min}^{B}=\frac{4.816}{k_{\perp}}$, where $k_{\perp}=k\sin\gamma$ is the
wavevector component perpendicular to the propagation ($k=2\pi/\lambda$). The
final cone angle on the sample was $\gamma=18.46^{\circ}$, the diameter of
$d_{min}^{B}=1.24~{}\mu$m at the first minimum or 0.88 $\mu$m at
$1/e^{2}$-intensity level. The length of the non-diffractig Bessel beam zone
was $\sim 400~{}\mu$m.

The central spot as well as every ring of the Bessel beam in the focal region
carries the same amount of energy. The average intensity and fluence per pulse
of the optical $E_{p}=44$ nJ energy used in mask patterning is estimated
considering the central spot of the diameter equal to the first minimum
$d_{min}^{B}$: $I_{p}\equiv F_{p}/t_{p}=21.76$ TW/cm2 and
$F_{p}=E_{p}/(\pi[d_{min}^{B}/2]^{2})=3.63$ J/cm2, respectively, for the
$t_{p}=167$ fs pulse ($ct_{p}=50~{}\mu$m axial extent). This calculation is
carried out the same way as for the Gaussian pulse and can be used for
protocol development. For the actual comparison of fluences between Bessel and
Gaussian pulses, the axial extent of the beams can be compared since all the
depth-of-focus (DoF) of a Bessel beam carries laterally distributed rings that
converge onto the optical axis. The same focal spot for a Gaussian beam can be
achieved with focusing using $NA_{G}$ lens defined by
$d^{B}_{min}=1.22\lambda/NA_{G}$. Then, the DoF or double the Rayleigh length
$2z_{R}=2\pi[d^{B}_{min}]^{2}/\lambda=18.8~{}\mu$m. This is $\sim 21$ times
shorter than the DoF of the Bessel beam. Hence, the corresponding fluence of
the Bessel beam reduced by this factor is 0.17 J/cm2 and is close to the
ablation threshold of Si. For the Gaussian pulse, the ablation threshold of Si
is 0.2 J/cm2 [19].

##  III Results

For this study, we chose Al2O3 mask due to its mechanical strength even at
nano-film conditions and being compatible with anti-reflection coating of Si
solar cells. Metal masks such as Cr or Ti make good adhesion to Si and can
withstand plasma etch used in this study [11], however, can deteriorate the
performance of Si in terms of increased surface recombination velocity.

###  III.1 Bessel beam patterning of Al2O3 mask

The Direct Laser Write (DLW) by 515 nm/167 fs laser pulses with Bessel beam
focusing was carried out on blank Si wafers coated by $30\pm 2$ nm and $40\pm
2$ nm Al2O3 films deposited by electron (e)-beam evaporation. Single and
double pulse per ablation site was used to write $1\times 1$ cm2 masks at 7.5
cm/s scan speed; the laser repetition rate was 0.6 MHz. It takes $\sim 20$ min
to pattern a $1\times 1$ cm2 area. These scan and repetition rate conditions
were chosen for the reliable positioning of ablation sites within $\pm 20$ nm
tolerance from the designated positions on a square grid with $\Lambda=3.1$
$\mathrm{\SIUnitSymbolMicro m}$ period as required for the PhC light trapping.
The writing speed can be further increased since $\sim 30$ nm positioning
errors were observed for the maximum tested writing speed of 10 cm/s; further
increae of a writing speed will be studied separately. The positioning errors
translate into mesa errors, especially at the corners of the etched pit. The
orientation of the mask was aligned to the base cut of Si wafer to have the
best etching conditions and smallest mesa regions. Irradiation of two pulses
per ablation site caused a larger debris field as well as holes in the mask
clogged by droplets of quenched melt (similar as in region-3 of Fig. 1). It is
noteworthy that debris on a mask and small changes in the shape of the hole-
in-mask was still acceptable for plasma etch. Thinner 30 nm Al2O3 masks showed
less debris and smaller holes, however, this thickness value was not selected
for the final fabrication on Si solar cells due to being prone to fracture
when undercut-etched. Thicker $\sim 50$ nm masks showed the presence of
brittle fracture at the rim of the opening and larger debris fields.

###  III.2 Plasma etch: optimization of conditions

Firstly, different chemical mixtures for dry plasma etching were investigated
with the aim of identifying those that yield the smoothest surfaces of Si
production. One of the simplest mixtures of SF6/O2 was the most promising
since the addition of CF4 and/or CHF3 resulted in more nano-textured surfaces
even if the etch rate was faster (for illustration see selected conditions
shown in Fig. A1).

The parameter space of optimal etch protocol was investigated with a change of
Inductively Coupled Plasma (ICP) power, which favors isotropic all-direction
material removal as well as Bias power, which defines anisotropic directional
(between electrodes) etching. Figure 2 shows optical images of pattern
evolution at large ICP power with bias at 0 W. When the corner posts of Si
below the mask become sub-1 $\mathrm{\SIUnitSymbolMicro m}$ in the lateral
cross-section, they appear colorful under optical inspection with an objective
with numerical aperture $NA=0.9$ (see inset in Fig. 2(d)). Different colors
are related to the small differences in shape of pillars which evolves from
small changes of ablated holes in the mask (see SEM inset image in (c)). This
stage of color appearance in pattern formation is indicative that the plasma
etch should be finished, since the Al2O3 mask breaks when the undercut etch
progresses and some of the posts will not support the mask (38-nm-thick in
Fig. 2). A mask is removed by immersion of the sample into an ultra-sonic bath
of isopropyl alcohol (IPA) for five minutes. If the mask is not fully removed
from the Si post, it appears colorful under optical inspection. Since alumina
is part of the anti-reflection coating in the final design of the solar
trapping surface, this is not considered a critical issue. For the final etch
protocol, the used slow etch rate resulted in posts of the same color (dark-
blue) under microscope observation, which indicates size and shape uniformity
(Fig. A6).

The PhC pattern was defined by etch mask on a $1\times 1$ cm2 area which had
unablated $4\times 4$ holes every 1 mm in x and y-directions of the pattern by
design. Those posts were added to support the Al2O3 membrane during dry plasma
etch (this action is illustrated in Fig. 2(e) where mask is fixed to the post
after it fails in the free-suspended region). Their contribution to an
increased reflection from Si solar cell with PhC light trapping patterns can
be considered small.

Even though ICP is known for favoring isotropic etch, at lower ICP powers the
etch becomes slower and the crystalline structure of Si becomes prominent on
the etched surface (Fig. A2). We chose low 30 W ICP power for the optimized
protocol and no oxygen was used to reduce surface oxidation, which caused non-
homogeneous islands/residue formation apparent under SEM inspection.

Next, the influence of the bias power was investigated. Figure 3 shows side
view cross sections aligned to the original surface of Si (a-c) etched under
increasing bias, i.e. directional plasma bombardment of the surface.
Elliptical markers are added as eye guides to facilitate a small increase of
depth at large powers. The effect of bias was not very prominent
understandably due to the small sub-1 $\mathrm{\SIUnitSymbolMicro m}$ opening
in the mask through which the etching of Si was made. An increased bias power
was also affecting the thickness of Al2O3 mask, which is undesirable. Zero
bias was used for the final protocol. During SEM observation of samples after
different times of plasma etch, it was observed that charging nano-debris are
increasingly present (Fig. 3(c,d)) after repeated exposures of the sample to
the air (in class 1000 cleanroom). This residue of etchants or/and oxidized Si
is minimized with reduced steps of air exposures as in the final fabrication.
With a lower content of oxygen in SF6/O2 mixture, the surface had less of the
oxidized residue and smoother surfaces. For the final protocol, no oxygen was
selected.

The process pressure in the chamber is another important parameter, especially
due to the ICP-dominated etch, which was selected for the PhC pattern. The He
gas flow below the sample was used for cooling the sample and was negligibly
contributing to the process pressure. The reactivity of plasma and the
isotropic nature of etching depended on the process pressure (Fig. A3). The
low pressure facilitated the smooth morphology of Si with the largest depth
(Fig. A3), which are both desirable virtues selected for the final protocol
(shown in Fig. A3(a)). The sensitivity of etched surface roughness and
presence of oxidized textures is illustrated in Fig. A4 where ambient pressure
changes are shown to strongly affect the morphology at larger ICP powers and
O2 flow even at 0 W Bias.

The key attribute of the proposed plasma etch protocol using a large area
laser patterning of 40-nm-thick Al2O3 mask is timing. Etching should be
stopped at the moment of the smallest pillar posts’ formation (Fig. 4). The
posts provide mechanical support for the mask. Undercut etching and narrow
mesa formation are already present even before posts are etched out. With the
posts of $\sim 150\pm 50$ nm in cross-section, the deepest PhC pattern is
formed with minimized mesa ridges. The final optimized recipe of PhC etch was
selected for the slow etch, which caused the formation of the smallest posts
(same deep blue color; Fig. A6), no oxygen in the etching gas mixture (only
SF6), 0 W Bias, 30 W ICP, and 1 Pa of process pressure. Such etching revealed
the (111) Si planes typically observed in wet KOH etching. Minimum oxidation
debris is formed on the surface, which is subsequently processed via wet etch
KOH passivation and anti-reflection coating steps for the final solar cell.

Figure 5(a) shows the appearance of etched PhC pattern using Al2O3 mask
patterned by fs-Bessel beam on the surface of actual solar cells. The total
etch time to achieve minimal mesa ridges (Fig. 5(b)) was approximately twice
as long for the entire Si wafer as compared with the test samples of $\sim
1\times 1$ cm2 patterned masks used for development of the etch protocol.
However, the quality of the etched surface, mechanical integrity of the mask,
and residue formation were the same as per the optimized recipe. One important
difference of the wafer scale etch was the slight radial variation of undercut
etching, which is an important and favorable attribute of the plasma
processing of PhC light trapping surface (Fig. 5). Dedicated wafer-size
optimization of the etch protocol is required. For that, we foresee the
requirement of a mechanically stronger mask for longer etching with an
undercut etch conditions as discussed in the next section.

##  IV Discussion

Figures 5 and A7 show currently achieved patterning of Polysilicon on Oxide
Interdigitated Back Contacted - POLO2-IBC cell. This POLO2-IBC architecture,
same materials, and fabrication protocols achieved 26.1% solar-to-electric
power conversion efficiency [18]. The mask writing throughput by the Bessel
beam was 1 cm2 in 20 min. A further increase is possible by combining xy-stage
scan with glavano-scanning [20]. For the Bessel beam with the axial extent of
the $\sim 400~{}\mu$m used in this study, there is a high tolerance for axial
position of the focal region on the Si surface with 40-nm-thick mask even for
an existing tilt of the sample and surface unevenness. Curved surfaces can be,
in principle, patterned with the Bessel beam. In this study, the pulse
duration corresponds to 50 $\mu$m in air as it travels through the non-
diffracting section of $400~{}\mu$m.

What was observed in the full wafer etch (Fig. 5), is a twice longer etch time
required for the same depth of the PhC pattern. This longer etch puts more
stringent requirements against mechanical failure of the patterned mask. A
higher ICP power would facilitate a faster etch rate, however, mask thinning
with more erosive plasma would take place [21]. A practical solution to make
better support for the etch mask would be to omit every second ablated hole in
every second line (see concept in Fig. 6(a)). Finite difference time domain
(FDTD) modeling (Lumerical, Ansys) of a toy-model of such PhC texture for
3-$\mu$m-thick Si is shown in (b) for the spectra of absorption and scattering
cross-sections. Such a pattern has very similar cross-sections as compared
with an inverted pyramid PhC. Figure 7 shows light intensity distributions at
very different wavelengths for a normally incident plane wave. The incident
intensity is normalized and the $|E|$ shows the enhancement of electrical
E-field due to interference inside Si especially at the longer wavelengths.
This toy-model simulation qualitatively captures the wave nature of PhC light
trapping which can perform above the Lambertian $4n^{2}$ ray optics trapping
used in solar cells [10].

The proposed model of the etch mask is amenable over large areas (10-inch
wafer) and larger. Due to a mechanically stronger mask, it can potentially be
thinner and smaller ablation holes can be printed if required for a smaller PC
period. The very same etch chemistry at a higher ICP and bias were also tested
for thinning of Si from 0.5 to 0.1 mm (Fig. A8). Such thinning opens the
possibility to use a compatible process to thin Si solar cells with backside
contacts to the required thickness before adding the PhC trapping surface,
which is the most efficient in thin tens-of-micrometers Si solar cells [15].

##  V Conclusions and outlook

Direct laser writing of a hard etch mask using ultra-short fs-laser pulses is
demonstrated for fabrication of PhC light trapping surface on actual high-
efficiency Polysilicon on Oxide Interdigitated Back Contacted - POLO2-IBC
cells. By using tens-of-nm Al2O3 hard mask ablated by fs-laser Bessel beam at
high $7.5$ cm/s scan speed at room conditions, makes this approach scalable
for solar cell applications. Periodic PhC patterns are required to break the
ray optics Lambertian limit of light trapping. A detailed account of explored
parameter space to optimize plasma etching is presented for different etch
chemistries, powers of ICP and Bias, and process pressure. Only an optical
inspection of the evolution of PhC formation is required for the determination
of etching time and degree of under-etching.

Scalability of such lithography-less and vacuum-less micro-patterning of PhC
textures over even larger areas is outlined and feasibility of light trapping
is qualitatively evaluated using a numerical toy model. Detailed performance
of these PhC-POLO2-IBC cells will be reported in a dedicated study. Here, the
experimental account for the selection of hard-mask thickness for direct laser
writing by Bessel beam and dry plasma etching were all synergetically
optimized and solved for an up-scalable application.

###### Acknowledgements.

We wish to acknowledge the support via ARC DP190103284 ”Photonic crystals: the
key to breaking the silicon-solar cell efficiency barrier” project (S.Ju.,
S.Jo.). This work was performed in part at the Melbourne Centre for
Nanofabrication (MCN) in the Victorian Node of the Australian National
Fabrication Facility (ANFF). Support by the German Federal Ministry for
Economic Affairs and Climate Protection, research project “27Plus6”
(FKZ03EE1056A) is acknowledged.

## References

  * of Michigan [2024] U. of Michigan, Photovoltaic energy factsheet, pub. no. css07-08 2023 (accessed: 09.03.2024).
  * Saga [2010] T. Saga, Advances in crystalline silicon solar cell technology for industrial mass production, NPG Asia Materials 2, 96–102 (2010).
  * Battaglia  _et al._ [2016] C. Battaglia, A. Cuevas, and S. De Wolf, High-efficiency crystalline silicon solar cells: status and perspectives, Energy Environ. Sci. 9, 1552 (2016).
  * Ram  _et al._ [2018] M. Ram, M. Child, A. Aghahosseini, D. Bogdanov, A. Lohrmann, and C. Breyer, A comparative analysis of electricity generation costs from renewable, fossil fuel and nuclear sources in G20 countries for the period 2015-2030, Journal of Cleaner Production 199, 687 (2018).
  * Lin  _et al._ [2023] H. Lin, M. Yang, X. Ru, G. Wang, S. Yin, F. Peng, C. Hong, M. Qu, J. Lu, L. Fang, C. Han, P. Procel, O. Isabella, P. Gao, Z. Li, and X. Xu, Silicon heterojunction solar cells with up to 26.81% efficiency achieved by electrically optimized nanocrystalline-silicon hole contact layers, Nat Energy 8, 789 (2023).
  * National Renewable Energy Laboratory (2024) [NREL] National Renewable Energy Laboratory (NREL), Best research-cell efficiency chart (accessed: 20.03.2024).
  * Yoshikawa  _et al._ [2017] K. Yoshikawa, H. Kawasaki, W. Yoshida, T. Irie, K. Konishi, K. Nakano, T. Uto, D. Adachi, M. Kanematsu, H. Uzu, and K. Yamamoto, Silicon heterojunction solar cell with interdigitated back contacts for a photoconversion efficiency over 26%, Nature Energy 2, 17032 (2017).
  * Andreani  _et al._ [2019] L. Andreani, A. Bozzola, P. Kowalczewski, M. Liscidini, and L. Redorici, Silicon solar cells: toward the efficiency limits, Advances in Physics: X 4, 1548305 (2019).
  * Aydin  _et al._ [2011] K. Aydin, V. Ferry, and R. Briggs, Broadband polarization-independent resonant light absorption using ultrathin plasmonic super absorbers, Nat Commun 2, 517 (2011).
  * Hsieh  _et al._ [2020] M.-L. Hsieh, A. Kaiser, S. Bhattacharya, S. John, and S.-Y. Lin, Experimental demonstration of broadband solar absorption beyond the lambertian limit in certain thin silicon photonic crystals, Scientific Reports 10, 11857 (2020).
  * Maksimovic  _et al._ [2022] J. Maksimovic, J. Hu, S. Ng, T. Katkus, G. Seniutinas, T. Pinedo Rivera, M. Stuiber, Y. Nishijima, S. John, and S. Juodkazis, Beyond Lambertian light trapping for large-area silicon solar cells: fabrication methods, Opto-Electron Adv 5, 210086 (2022).
  * Rienäcker  _et al._ [2023] M. Rienäcker, J. Maksimovic, Y. Larionova, S.-H. Ng, T. Katkus, T. Pinedo Rivera, M. Stuiber, J. Krügener, F. Haase, R. Brendel, S. John, S. Juodkazis, and R. Peibst, Integration of photonic crystals in highly efficient POLO2-IBC cells - interplay between structural, optical, and electrical properties, in  _Program_ , 13th Int. Conf. on Crystalline Silicon Photovoltaics (Delft, Netherlands, 2023) presented at SiliconPV April 11-14 2023.
  * Yamamoto [2003] K. Yamamoto, Thin-film crystalline silicon solar cells, JSAP International 7, 12 (2003).
  * Bernardi  _et al._ [2013] M. Bernardi, M. Palummo, and J. C. Grossman, Extraordinary sunlight absorption and one nanometer thick photovoltaics using two-dimensional monolayer materials, Nano Letters 13, 3664 (2013).
  * Bhattacharya and John [2018] S. Bhattacharya and S. John, Designing high-efficiency thin silicon solar cells using parabolic-pore photonic crystals, Phys. Rev. Appl. 9, 044009 (2018).
  * Li  _et al._ [2024] Y. Li, X. Ru, M. Yang, Y. Zheng, S. Yin, C. Hong, F. Peng, M. Qu, C. Xue, J. Lu, L. Fang, C. Su, D. Chen, J. Xu, C. Yan, Z. Li, X. Xu, and Z. Shao, Flexible silicon solar cells with high power-to-weight ratios, Nature 626, 105–110 (2024).
  * McGloin and Dholakia [2005] D. McGloin and K. Dholakia, Bessel beams: Diffraction in a new light, Contemporary Physics 46, 15 (2005).
  * Haase  _et al._ [2018] F. Haase, C. Hollemann, S. Schäfer, A. Merkle, M. Rienäcker, J. Krügener, R. Brendel, and R. Peibst, Laser contact openings for local poly-si-metal contacts enabling 26.1%-efficient POLO-IBC solar cells, Solar Energy Materials and Solar Cells 186, 184 (2018).
  * Bonse  _et al._ [2002] J. Bonse, S. Baudach, J. Krüger, W. Kautek, and M. Lenzner, Femtosecond laser ablation of silicon–modification thresholds and morphology, Applied Physics A 74, 19 (2002).
  * Jonušauskas  _et al._ [2019] L. Jonušauskas, D. Gailevičius, S. Rekštytė, T. Baldacchini, S. Juodkazis, and M. Malinauskas, Mesoscale laser 3d printing, Opt. Express 27, 15205 (2019).
  * Maksimovic  _et al._ [2023] J. Maksimovic, H. Mu, D. Smith, T. Katkus, M. Vaičiulis, R. Aleksiejūnas, G. Seniutinas, S. H. Ng, and S. Juodkazis, Laser-patterned alumina mask and mask-less dry etch of si for light trapping with photonic crystal structures, Micromachines 14 (2023).

##  Appendix A Appendix

Here, a set of supplementary figures illustrates and summarises experimental
observations with a selection of qualitative illustrations of the main
consequence of particular etch parameter change. Figure captions are self-
explanatory with technical details and main conclusions.

Figure A1:  Typical PhC structures at different etch chemistries. (a) CF4
test: 2 Pa process pressure, 200 W ICP, 0 W Bias, 50/5 ccm of CF4/O2 for 60
min. Strong under-etch. Thickness of Al2O3 mask 37.6 nm, Bessel beam. (b) same
as (a) only 50 W ICP for 150 min. (c) Etch mixture with CHF3: 2.5 Pa process
pressure , 180 W ICP, 0 W Bias, 50/10/10 ccm of SF6/O2/CHF3. Thickness of
Al2O3 mask 28.4 nm; Bessel beam was used for ablation. The inset in (c) shows
optical images of etched pattern evolution, which turns into colored posts
after 60 min of etch.  Figure A2:  Dependence of etch on ICP power in SF6/O2
mixture (ratios are shown at the bottom of panels) at 0 W Bias. Representative
SEM images of top-view after different times of etching. Presence of oxygen
shows propensity of oxidation on the etched surface. At low ICP powers,
crystolographic structure typical for the Si $\left<001\right>$ becomes
discernable, see (a) vs. (d). The thickness of Al2O3 was 38 nm; mask writing
by Bessel beam.  Figure A3:  SEM images of etched PhCs on Si at different
process pressures: 1 Pa (a), 2 Pa (b), 4 Pa (c) and 2 Pa (d). Top-insets show
side-view cross-section. Etch time is shown at the bottom of panels; when
repeated etching was applied, times were shown with a summation marker “+” at
the bottom of panels. RIE was carried out at 0 W Bias, 30 W ICP power (a-c)
and 15 W (d), at 50/0 ccm SF6/O2 gas mixture. Thickness of Al2O3 mask was 28.4
nm; mask writing was with Bessel beam. Conditions in (a) were selected for
fabrication of PhC light trapping patterns on solar cells.  Figure A4:  SEM
images of etched PhCs on Si at different process pressures: 1 Pa (a), 2.5 Pa
(b), 5 Pa (c) and 10 Pa (d) and 0.5 Pa (e). Etch time is shown at the bottom
of panels. RIE was carried out at 0 W Bias in all cases, while ICP was 200 W
and 100 W (marked); Ar at 50 ccm was added for (c). The ratio of flow rate of
SF6/O2 gas mixture is shown below the panels. Thickness of Al2O3 mask was 28.4
nm; mask writing was with Bessel beam.  Figure A5:  Time, the ratio of flow
rate of SF6/O2 gas mixture, and pressure are show. Bias power was 0 W and that
of ICP 30 W. Single run of etching (a) and with three repeated inspections (b)
to achieve minimal mesa regions and under etching. Top-inset in (b) shows the
typical opening in the mask used for this etching with its projections on the
etched inverted pyramids as central circles. Thickness of Al2O3 mask was 38
nm; mask writing was with Bessel beam. The angle of $54.7^{\circ}$ between top
surface (100) plane and (111) plane is shown. Repeated exposure to room
conditions caused a larger amount of the oxidised residue shown in (b). Under
etching conditions (b) affected the depth and slope of side walls and was
stronger at larger ambient pressure (see Fig. A3).  Figure A6:  The optimal
etch conditions. Optical images of etching progression for a pressure of 1 Pa,
30 W ICP, 0 W Bias, and 50/0 ccm of SF6/O2 at the same position on the sample.
The smallest corner posts are of the same colour (dark blue) observed after 90
min etch; see different colors at conditions of faster etch in Fig. 2(d). At
this condition, etching is stopped and the mask is removed in ultra-sonic
bath. Thickness of Al2O3 mask was 38 nm; mask writing was with Bessel beam.
Figure A7:  (a) Back-side view of Polysilicon on Oxide Interdigitated Back
Contacted - POLO2-IBC cell with the top-inset image of patterned etch mask by
Bessel beam 515 nm/167 fs, pulse energy $E_{p}=43.8$ nJ. The thickness of Si
solar cell 190 $\mu$m. Insets show laser patterned mask and final PhC light
trapping surface and optimised plasma etch and mask removal. (b) Optical dark-
field image of laser inscribed mask in 38-nm-thick Al2O3.  Figure A8:
Thinning of $500~{}\mu$m wafers down to sub-100 $\mu$m thickness (Samco ICP-
RIE). Plasma thinning of p-type Si at different bias powers, when ICP power
was 700 W, process pressure 1 Pa, He back pressure 2 kPa, etching chemistry
SF6. Etching time was 90 min; error bars $15\%$. There was a strong dependence
on the size of Si substrate, i.e. a cm-sized sample was etched $\sim 50\%$
faster than a wafer. The most clean debris-free surface was obtained at the
highest bias power of 160 W, when the etch rate was $4.6~{}\mu$m/min and the
roughness was $\sim 25$ nm. Noteworthy, n-type Si had a strongly light
scattering appearance at the visible spectral range after etching; the
intrinsic and p-type Si was etched identically.