# Actinide-boosting $r$ Process in Black Hole–Neutron Star Merger Ejecta

Shinya Wanajo<EMAIL_ADDRESS>Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), Am Mühlenberg 1, D-14476
Potsdam-Golm, Germany  Sho Fujibayashi  Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), Am Mühlenberg 1, D-14476
Potsdam-Golm, Germany  Kota Hayashi  Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto, 606-8502, Japan  Kenta Kiuchi  Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), Am Mühlenberg 1, D-14476
Potsdam-Golm, Germany  Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University,
Kyoto, 606-8502, Japan  Yuichiro Sekiguchi  Center for Gravitational Physics
and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto, 606-8502, Japan  Department of Physics, Toho University,
Funabashi, Chiba 274-8510, Japan  Masaru Shibata  Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), Am Mühlenberg 1, D-14476
Potsdam-Golm, Germany  Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University,
Kyoto, 606-8502, Japan

###### Abstract

We examine nucleosynthesis in the ejecta of black hole-neutron star mergers
based on the results of self-consistent, long-term neutrino-radiation-
magnetohydrodynamics simulations for the first time. We find that the
combination of dynamical and post-merger ejecta reproduces a solar-like
$r$-process pattern. Moreover, the enhancement level of actinides is highly
sensitive to the distribution of both electron fraction and the velocity of
the dynamical ejecta. Our result implies that the mean electron fraction of
dynamical ejecta should be $\sim 0.05$–0.08 in order to reconcile the
nucleosynthetic abundances with those in $r$-process-enhanced, actinide-
boosted stars. This result provides an important constraint for nuclear
equations of state.

Introduction.—One of the long-standing issues in nuclear astrophysics is to
determine the mechanism that ensures a robust solar-like $r$-process pattern
as can be found in $r$-process-enhanced metal-poor stars (e.g., Ref. Cowan
_et al._ (2021)). Recent work Fujibayashi  _et al._ (2022) based on long-term
neutrino-radiation-hydrodynamics simulations of neutron star merger remnants
suggests that a combination of dynamical and post-merger ejecta leads to a
solar-like $r$-process pattern (see also Just  _et al._ (2015); Martin  _et
al._ (2015); Kullmann  _et al._ (2022)). However, such nucleosynthesis studies
have not yet been carried out for black hole–neutron star mergers.

Another salient feature observed in $r$-process-enhanced stars is the presence
of stars with an excess of the radioactive species Th with respect to a stable
element, e.g., Eu. This excess is known as the “actinide boost”. The Th/Eu
ratio ranges over a factor of five, $0.14\leq\mathrm{Th/Eu}\leq 0.76$ Ji and
Frebel (2018); Holmbeck  _et al._ (2018). Among those, about one-third of
$r$-process-enhanced stars are classified as actinide-boost stars with
$\mathrm{Th/Eu}>0.5$ (the production ratio of $\mathrm{Th/Eu}>0.9$ assuming a
stellar age of 13 Gyr). A nucleosynthesis study Holmbeck  _et al._ (2019)
using a single thermodynamic trajectory has shown that actinides are
overproduced if the initial neutron-richness is sufficient for fission
recycling. However, it is currently an open question if the actinide boost can
be also realized when adopting a large number of thermodynamic trajectories
based on the hydrodynamics simulations of binary neutron (or black
hole–neutron) star mergers.

In this Letter, we present the first nucleosynthesis study based on self-
consistent, long-term neutrino-radiation-magnetohydrodynamics simulations in
Refs. Hayashi  _et al._ (2022a, b). The aim of this work is two-fold. One is
to examine if black hole-neutron star mergers can reproduce solar-like
$r$-process patterns as can be seen in $r$-process-enhanced stars. Second we
attempt to explore the conditions in which actinide boost is realized. Unlike
binary neutron star mergers (Fujibayashi  _et al._ , 2022), thermal effects
(before $r$-processing) are unimportant in the cold tidal ejecta (i.e.,
without weak interaction) of black hole-neutron star mergers. This fact
drastically simplifies the physical conditions relevant to nucleosynthesis,
and also enables us to testify the neutron-richness in the inner crust of a
neutron star.

Models.—We adopt the results of neutrino-radiation-magnetohydrodynamics
simulations of black hole–neutron star mergers in Refs. Hayashi  _et al._
(2022a, b). In their models, the temperature-dependent, tabulated nuclear
equations of state DD2 Banik  _et al._ (2014) or SFHo Steiner  _et al._ (2013)
were adopted, and early dynamical and late post-merger mass ejections were
self-consistently calculated in a single three-dimensional computational
domain.

In this study, we take their models of Q4B5H with DD2 and SFHo (hereafter
referred to as Q4B5H-DD2 and Q4B5H-SFHo) and Q6B5L with DD2 (Q6B5L-DD2). Here,
Q4 and Q6, B5, and H and L stand for the ratios of black-hole masses
($5.4\,M_{\odot}$ and $8.1\,M_{\odot}$) with respect to that of a neutron star
($1.35\,M_{\odot}$), the initial maximum magnetic-field strength ($5\times
10^{16}$ G), and the grid spacing for the finest refinement level (270 m and
400 m), respectively. The ejecta masses were approximately saturated for all
models at the end of simulations ($\gtrsim 1$–2 s). The total ejecta masses
are $M_{\mathrm{ej}}/M_{\odot}=0.070$, 0.031, and 0.067 for Q4B5H-DD2,
Q4B5H-SFHo, and Q6B5L-DD2, respectively.

Figure 1: Mass histograms with respect to the total ejecta masses for
$Y_{\mathrm{e}}$ (top) with an interval of $\Delta Y_{\mathrm{e}}=0.01$, and
$v/c$ (bottom) with an interval of $\Delta(v/c)=0.01$.

For nucleosynthesis calculations during post-processing, about 1400 (400 for
the dynamical component) tracer particles are generated for each model as in
Ref. Fujibayashi  _et al._ (2022). Fig. 1 displays the mass histograms with
respect to the total ejecta mass for electron fraction (proton number per
nucleon) $Y_{\mathrm{e}}$ (top) and velocity $v$ with respect to the speed of
light $c$ (bottom). We find both narrow and broad distributions for each model
at $Y_{\mathrm{e}}\sim 0.05$–0.06 and 0.15–0.4 originating from dynamical and
post-merger ejecta, respectively. The masses of the dynamical ejecta (defined
as those with $Y_{\mathrm{e}}<0.08$) are $M_{\mathrm{ej}}/M_{\odot}=0.039$
(56% of the total ejecta mass), 0.012 (39%), and 0.046 (69%) for Q4B5H-DD2,
Q4B5H-SFHo, and Q6B5L-DD2, respectively. We also find narrow and broad
distributions at $v/c\sim 0.03$–0.07 and 0.1–0.35, which come from post-merger
and dynamical ejecta, respectively.

Nucleosynthesis.—Nucleosynthetic abundances in each tracer particle are
computed in post-processing using the nuclear reaction network code rNET
described in Ref. Wanajo  _et al._ (2018). Fission fragment distributions as
well as the number of prompt neutrons per event are computed with the code of
the GEF model Schmidt and Jurado (2010) (version 2021/1.1;
http://www.khschmidts-nuclear-web.eu/GEF-2021-1-1.html) for spontaneous and
neutron-induced channels, with the latter assuming an incident neutron energy
of 0.1 MeV. For the $\beta$-delayed channel, the fragment distributions for
the neutron-induced fission are adopted.

Thermodynamic histories are determined based on matter density in tracer
particles. The time evolution of the temperature is computed with the chemical
composition, density, and entropy by using the equations of state for a gas
consisting of non-degenerate ions, arbitrarily degenerate, relativistic
electrons and positrons, and photons Timmes and Swesty (2000). The energy
feedback from $\beta$ decay, $\alpha$ decay, and fission is reflected in the
increase of entropy at each time step, where 40% of the energy per $\beta$
decay is assumed to be lost as free-streaming neutrinos (Hotokezaka  _et al._
, 2016).

Figure 2: Left: total isobaric abundances (gray) as well as those from
dynamical (blue) and post-merger (red) ejecta components at the end of
simulation (1 yr; all trans-Pb nuclei except for Th and U are assumed to have
decayed) for Q4B5H-DD2. The black circles denote the solar $r$-residuals
Prantzos  _et al._ (2020), which are vertically shifted to match the
calculated total abundance of 153Eu. Right: same as left panel but taken from
dynamical ejecta with different fission-fragment distributions of GEF Schmidt
and Jurado (2010) (blue), KT75 Kodama and Takahashi (1975) (red), and
KT75-symmetric (olive).

The initial compositions for dynamical and post-merger ejecta are treated
separately as follows. For the dynamical (or cold) ejecta, the initial
condition is set at the neutron-drip density $\rho_{\mathrm{drip}}=4\times
10^{11}$ g cm-3. The initial entropy is taken at this density in a given
tracer particle (typically less than 1 $k_{\mathrm{B}}$/nucleon, where
$k_{\mathrm{B}}$ is Boltzmann constant). The matter is assumed to consist of
free neutrons and single representative heavy nuclei of $(Z,A)$ with the
numbers per nucleon of $Y_{\mathrm{n}}$ and $Y(Z,A)$, respectively. These
values satisfy mass conservation, $1=Y_{\mathrm{n}}+AY(Z,A)$, as well as
charge neutrality, $Y_{\mathrm{e}}=ZY(Z,A)$. We further assume that the free
energy of cold ($\ll 5$ GK) tidal ejecta is given solely by the mass of
species per nucleon, i.e.,

$\displaystyle f(Z,A)/c^{2}$ $\displaystyle=$ $\displaystyle
m_{\mathrm{n}}Y_{\mathrm{n}}+m(Z,A)Y(Z,A)$ $\displaystyle=$ $\displaystyle
m_{\mathrm{n}}(1-AY_{\mathrm{e}}/Z)+m(Z,A)Y_{\mathrm{e}}/Z,$

where $m_{\mathrm{n}}$ and $m(Z,A)$ are the masses of a neutron and a nucleus
$(Z,A)$, respectively. For a given $Y_{\mathrm{e}}$, minimization of $f(Z,A)$
gives the representative $(Z,A)$. We find that the nucleus is uniquely
determined to be $(Z,A)=(36,118)$ at the shell closure of $N=82$, i.e., 118Kr,
for $Y_{\mathrm{e}}<0.3$ with the nuclear masses of HFB-21. Note that the
equations of state adopted in this study also give the same representative
nucleus (see also the same result in Ref. Baym  _et al._ (1971)).

For post-merger (or hot) ejecta, each calculation is initiated from the
density at which the temperature decreases to 10 GK in a given tracer
particle. The initial compositions are set to $1-Y_{\mathrm{e}}$ and
$Y_{\mathrm{e}}$ for free neutrons and protons, respectively, where
$Y_{\mathrm{e}}$ here is taken from the final value in a given tracer
particle. At such high temperature, the compositions immediately relax to
those in nuclear statistical equilibrium.

$r$-Process Abundance Patterns.—The nucleosynthesis result for model Q4B5H-DD2
is presented in Fig. 2 (left). We find that the dynamical and post-merger
ejecta are responsible for the heavy ($A>130$) and light ($A<130$) components
of $r$-process abundances, respectively. The ensemble of both components is in
good agreement with the pattern of $r$-process residuals to the solar
abundances ($r$-residuals hereafter) Prantzos  _et al._ (2020) for $A\approx
90$–210.

The right panel of Fig. 2 compares the abundances in dynamical ejecta for
Q4B5H-DD2 with different fission-fragment distributions of GEF, KT75
(phenomenological double-Gaussian distributions in Ref. Kodama and Takahashi
(1975)), and KT75-symmetric (assuming symmetric fission with single-Gaussian
distributions). We find that the total abundances for lanthanides with
$A<170$, including the rare-earth peak, are entirely determined by the
fission-fragment distributions adopted. While the GEF prediction leads to an
excellent match of the lanthanide distribution to that of the solar
$r$-residuals, no agreement can be seen when the KT75 or KT75-symmetric
distributions are adopted. It is important to note that the fission properties
in the neutron-rich region including fragment distributions are currently very
uncertain Goriely  _et al._ (2015); Vassh  _et al._ (2019); Lemaître  _et al._
(2021). However, we regard the adopted GEF model as a reasonable choice for
the purposes in this study, in which robustly reproducing a solar-like
$r$-process pattern for lanthanides, including Eu, is particularly important.

Fig. 3 compares the elemental abundances normalized by those of Eu for all
models. It is noteworthy that the abundance patterns of all models are similar
for $Z=40$–81 ($A=90$–205), which are in good agreement with that of the solar
$r$-residuals. Normalized abundances of $r$-process-enhanced stars J0954+5246
Holmbeck  _et al._ (2018) (open circles), CS 31082-001 Siqueira Mello  _et
al._ (2013) (crosses), and DES J033523-540407 Ji and Frebel (2018) (filled
circles) are also plotted. These stars exhibit the highest, modest, and lowest
values of measured Th/Eu ratio, respectively. The abundances of Pb, Bi, Th,
and U at 13 Gyr are also shown by dotted lines, given that these stars were
born several 100 Myr after the big bang (e.g., Refs. Wanajo  _et al._ (2021);
Hirai  _et al._ (2022)). We find that the abundance patterns are in reasonable
agreement with those of $r$-process-enhanced stars. In particular, the
normalized abundances of Pb, Bi, Th, and U in model Q4B5H-DD2 are in good
agreement with those in CS 31082-001, one of the actinide-boost stars.

Actinide Boost.—Here, we examine the dependencies of nucleosynthetic outcomes
in dynamical ejecta (which are responsible for the production of lanthanides
and heavier elements) on $Y_{\mathrm{e}}$ and outflow velocity $v$. The left
panel of Fig. 4 shows the abundances of Eu and Th (at 1 yr; all trans-Pb
nuclei except for Th and U are assumed to have decayed) as well as their ratio
Th/Eu for model Q4B5H-DD2 as functions of $Y_{\mathrm{e}}$. Here,
$Y_{\mathrm{e}}$ is taken to be a free parameter (taking into account possible
variation due to the nuclear EOS adopted), replacing the original values by a
single $Y_{\mathrm{e}}$ (with an interval of 0.005). The Eu abundance becomes
minimal at $Y_{\mathrm{e}}=0.155$ as the nuclear flow proceeds toward heavier
nuclei. With a reduction of $Y_{\mathrm{e}}$, Eu becomes abundant owing to
fission recycling. The Th abundance reaches a maximum at
$Y_{\mathrm{e}}=0.125$. For $Y_{\mathrm{e}}<0.125$, we find the effects of
“fission waves” for Eu and Th in response to the second and third fission
recycling. However, the $Y_{\mathrm{e}}$-offsets of local maxima for Th and Eu
lead to a convergence of the Th/Eu ratio towards lower $Y_{\mathrm{e}}$ with a
local maximum at $Y_{\mathrm{e}}=0.07$.

Figure 3: Elemental abundances normalized by the abundances of Eu for all
models. The normalized solar $r$-residuals Prantzos  _et al._ (2020) are also
displayed by a gray curve. The dotted lines show the abundances of Pb, Bi, Th,
and U at 13 Gyr. The normalized stellar abundances (with error bars) of
J0954+5246 Holmbeck  _et al._ (2018) (open circles), CS 31082-001 Siqueira
Mello  _et al._ (2013) (crosses), and DES J033523-540407 Ji and Frebel (2018)
(filled circles) are also shown.

Figure 4: Abundances of Eu and Th (at 1 yr; all trans-Pb nuclei except for Th
and U are assumed to have decayed) as well as their ratio Th/Eu for Q4B5H-DD2
(dynamical ejecta) as a function of $Y_{\mathrm{e}}$. In the left panel, the
factor $10^{\sigma}$ is also shown, which indicates the deviation from the
solar $r$-residual distribution. In the right panel, Th/Eu ratios calculated
using tracer particles with selected outflow velocities are also shown. The
symbols indicate the original Th/Eu ratios for Q4B5H-DD2 (circle), Q4B5H-SFHo
(diamond), and Q6B5L-DD2 (square) at each mass-averaged $Y_{\mathrm{e}}$. The
horizontal lines (with errors shown by shaded areas) indicate the
observational ratios of J0954+5246 Holmbeck  _et al._ (2018) (magenta), CS
31082-001 Siqueira Mello  _et al._ (2013) (cyan), and DES J033523-540407 Ji
and Frebel (2018) (purple), where the abundances of Th are corrected to those
of 13 Gyr ago.

Assuming that black hole-neutron star mergers are, in part, the sources of
$r$-process elements, which have enriched $r$-process-enhanced stars with
solar-like abundance patterns, we limit the range of $Y_{\mathrm{e}}$ for
subsequent discussion. The agreement level of the computed abundances with the
solar $r$-residuals can be evaluated by the deviation factor $10^{\sigma}$
with $\sigma$ defined by
$\sigma^{2}=\sum_{A}\left[\log_{10}N(A)-\log_{10}N_{\odot}(A)\right]^{2}/n_{\mathrm{tot}}$
Fujibayashi  _et al._ (2022), where $N(A)$ and $N_{\odot}(A)$ are the
normalized isobaric abundances for the dynamical ejecta and the solar
$r$-residuals, respectively. The sum runs over the $n_{\mathrm{tot}}=67$
isobars from $A=139$ to 205, which are predominantly produced in the dynamical
ejecta. For Q4B5H-DD2, Q4B5H-SFHo, and Q6B5L-DD2 with original
$Y_{\mathrm{e}}$ values, $10^{\sigma}=2.2$, 2.1, and 2.0, respectively,
indicating agreements with the solar $r$-residual pattern within about a
factor of 2. Here, we regard the abundance patterns with $10^{\sigma}\lesssim
2$ as appropriate models for $r$-process sources (see also Refs. Goriely
(1999); Fujibayashi  _et al._ (2022)). Thus, hereafter we restrict the range
of our attention to $Y_{\mathrm{e}}<0.10$.

As can be seen in the right panel of Fig. 4, the Th/Eu ratio reaches the local
maximum at $Y_{\mathrm{e}}\approx 0.07$ and decreases for
$Y_{\mathrm{e}}\lesssim 0.07$. This can be understood as a result of the
fission waves in Eu and Th abundances described above. The symbols indicate
the Th/Eu ratios (at mass-averaged $Y_{\mathrm{e}}$) for Q4B5H-DD2 (circle),
Q4B5H-SFHo (diamond), and Q6B5L-DD2 (square). The plot also displays the Th/Eu
ratios for tracer particles with selected outflow velocities. We find
substantial variation of the Th/Eu ratio at a given $Y_{\mathrm{e}}$, with a
tendency for this ratio to take on a higher value for a slower ejecta
velocity. This is due to a higher temperature during an $r$-process in the
slowly (and nearly adiabatically) expanding material. The energy feedback
(chiefly from $\beta$ decay) increases temperature to $\sim 0.7$–2 GK, at
which an equilibrium between radiative neutron capture and its inverse is
established. Under this condition, a higher temperature shifts the $r$-process
path closer to $\beta$ stability, at which $\beta$-decay timescales are
longer. As a result, for a given $Y_{\mathrm{e}}$, more abundances tend to
populate in the range with $A=232$–252 (few fissile nuclei) than $A\geq 254$
(predominantly fissile nuclei) at the end of an $r$ process. The former and
the latter include, respectively, the progenitors of Th (by $\alpha$ decay)
and Eu (by fission). Therefore, the resulting Th/Eu ratio becomes higher for a
smaller $v$.

The Th/Eu ratio also tends to be higher for a higher $Y_{\mathrm{e}}$ (with a
fixed $v$) because of an initially larger amount of heavy nuclei (that
generate radioactive energies), although the effect of fission waves appears
to dominate for $Y_{\mathrm{e}}>0.05$. It is important to note that the Th/Eu
ratio depends on theoretical $\beta$-decay predictions Holmbeck  _et al._
(2019) (slower rates are favored for actinide boost). The adopted
$\beta$-decay rates (GT2 Tachibana  _et al._ (1990) based on the HFB-21 mass
prediction Goriely  _et al._ (2010)) appear to result in a relatively high
Th/Eu ratio according to Ref. Kullmann  _et al._ (2022).

The stellar Th/Eu ratios of J0954+5246 Holmbeck  _et al._ (2018), CS 31082-001
Siqueira Mello  _et al._ (2013), and DES J033523-540407 Ji and Frebel (2018)
are also presented by horizontal lines (and shaded areas) in the right panel
of Fig. 4, where the observational values of Th are corrected to those 13 Gyr
ago (a few Gyr of difference does not substantially change our conclusion).
Note that the errors for Th ($\pm 0.10$ dex, $\pm 0.10$ dex, and $\pm 0.23$
dex, respectively; equal or larger than those for Eu) are adopted from the
literature, although the errors for Th/Eu ratios are unknown. We find that the
models with DD2 result in actinide boost, i.e., Th/Eu $>0.9$, which reconcile
with the Th/Eu ratio of an actinide-boost star CS 31082-001 but not the
highest measured Th/Eu ratio of J0954+5246. However, if we allow an error of
0.2 dex for J0954+5246, the Th/Eu ratios for these models become consistent
with the stellar value. The model with SFHo does not result in actinide boost.

According to the results with single $Y_{\mathrm{e}}$ values for Q4B5H-DD2
(red curve), only a range of $Y_{\mathrm{e}}\sim 0.05$–0.08 meets the level of
actinide boost, Th/Eu $>0.9$, or $Y_{\mathrm{e}}\sim 0.06$–0.08 for the most
actinide-boosted star J0954+5246.

Conclusion.—We have conducted the first exploration of nucleosynthesis based
on self-consistent magnetohydrodynamics simulations of black hole–neutron star
mergers Hayashi  _et al._ (2022a, b). Lighter ($A<130$) and heavier ($A>130$)
$r$-process nuclei were synthesized in the early dynamical and late-time post-
merger components, respectively, the ensemble of these reproducing solar-like
$r$-process patterns. This indicates that in addition to binary neutron star
mergers, black hole-neutron star mergers can also be galactic $r$-process
sites.

Moreover, our result has demonstrated that the presence of actinide-boost
stars, including the star with the highest Th/Eu ratio, can be explained if
the range of bulk $Y_{\mathrm{e}}$ in the dynamical ejecta is $\sim
0.05$–0.08. This can be an important constraint on nuclear equations of state,
such as symmetry energy, under the condition that the high-density matter
(above the saturation density) is tidally ejected.

Acknowledgement.—We thank L. Held for his careful reading of the manuscript.
This work was in part supported by Grant-in-Aid for Scientific Research (Grant
No. JP20H00158) of Japanese MEXT/JSPS. Numerical computations were performed
on Cobra, Raven and Sakura clusters at Max Planck Computing and Data Facility
and on Yukawa21 at Yukawa Institute for Theoretical Physics, Kyoto University.

## References

  * Cowan  _et al._ (2021) J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langanke, G. Martínez-Pinedo, and F.-K. Thielemann, Reviews of Modern Physics 93, 015002 (2021).
  * Fujibayashi  _et al._ (2022) S. Fujibayashi, K. Kiuchi, S. Wanajo, K. Kyutoku, Y. Sekiguchi, and M. Shibata, arXiv e-prints , arXiv:2205.05557 (2022), arXiv:2205.05557 [astro-ph.HE] .
  * Just  _et al._ (2015) O. Just, A. Bauswein, R. Ardevol Pulpillo, S. Goriely, and H. T. Janka, Mon. Not. Roy. Astron. Soc. 448, 541 (2015), arXiv:1406.2687 [astro-ph.SR] .
  * Martin  _et al._ (2015) D. Martin, A. Perego, A. Arcones, F. K. Thielemann, O. Korobkin, and S. Rosswog, Astrophys. J. 813, 2 (2015), arXiv:1506.05048 [astro-ph.SR] .
  * Kullmann  _et al._ (2022) I. Kullmann, S. Goriely, O. Just, A. Bauswein, and H. T. Janka, arXiv e-prints , arXiv:2207.07421 (2022), arXiv:2207.07421 [astro-ph.HE] .
  * Ji and Frebel (2018) A. P. Ji and A. Frebel, Astrophys. J. 856, 138 (2018), arXiv:1802.07272 [astro-ph.SR] .
  * Holmbeck  _et al._ (2018) E. M. Holmbeck, T. C. Beers, I. U. Roederer, V. M. Placco, T. T. Hansen, C. M. Sakari, C. Sneden, C. Liu, Y. S. Lee, J. J. Cowan, and A. Frebel, Astrophys. J. Lett. 859, L24 (2018), arXiv:1805.11925 [astro-ph.SR] .
  * Holmbeck  _et al._ (2019) E. M. Holmbeck, T. M. Sprouse, M. R. Mumpower, N. Vassh, R. Surman, T. C. Beers, and T. Kawano, Astrophys. J. 870, 23 (2019), arXiv:1807.06662 [astro-ph.SR] .
  * Hayashi  _et al._ (2022a) K. Hayashi, S. Fujibayashi, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Phys. Rev. D 106, 023008 (2022a).
  * Hayashi  _et al._ (2022b) K. Hayashi, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, arXiv e-prints , arXiv:2211.07158 (2022b), arXiv:2211.07158 [astro-ph.HE] .
  * Banik  _et al._ (2014) S. Banik, M. Hempel, and D. Bandyopadhyay, Astrophys. J. Suppl. 214, 22 (2014), arXiv:1404.6173 [astro-ph.HE] .
  * Steiner  _et al._ (2013) A. W. Steiner, M. Hempel, and T. Fischer, Astrophys. J. 774, 17 (2013), arXiv:1207.2184 [astro-ph.SR] .
  * Wanajo  _et al._ (2018) S. Wanajo, B. Müller, H.-T. Janka, and A. Heger, Astrophys. J. 852, 40 (2018), arXiv:1701.06786 [astro-ph.SR] .
  * Schmidt and Jurado (2010) K.-H. Schmidt and B. Jurado, Phys. Rev. Lett. 104, 212501 (2010), arXiv:0912.3651 [nucl-th] .
  * Timmes and Swesty (2000) F. X. Timmes and F. D. Swesty, Astrophys. J. Suppl. 126, 501 (2000).
  * Hotokezaka  _et al._ (2016) K. Hotokezaka, S. Wanajo, M. Tanaka, A. Bamba, Y. Terada, and T. Piran, Mon. Not. Roy. Astron. Soc. 459, 35 (2016), arXiv:1511.05580 [astro-ph.HE] .
  * Prantzos  _et al._ (2020) N. Prantzos, C. Abia, S. Cristallo, M. Limongi, and A. Chieffi, Mon. Not. Roy. Astron. Soc. 491, 1832 (2020), arXiv:1911.02545 [astro-ph.GA] .
  * Kodama and Takahashi (1975) T. Kodama and K. Takahashi, Nucl. Phys. A 239, 489 (1975).
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Goriely  _et al._ (2015) S. Goriely, A. Bauswein, O. Just, E. Pllumbi, and H. T. Janka, Mon. Not. Roy. Astron. Soc. 452, 3894 (2015), arXiv:1504.04377 [astro-ph.SR] .
  * Vassh  _et al._ (2019) N. Vassh, R. Vogt, R. Surman, J. Randrup, T. M. Sprouse, M. R. Mumpower, P. Jaffke, D. Shaw, E. M. Holmbeck, Y. Zhu, and G. C. McLaughlin, Journal of Physics G Nuclear Physics 46, 065202 (2019), arXiv:1810.08133 [nucl-th] .
  * Lemaître  _et al._ (2021) J. F. Lemaître, S. Goriely, A. Bauswein, and H. T. Janka, Phys. Rev. C 103, 025806 (2021), arXiv:2102.03686 [nucl-th] .
  * Siqueira Mello  _et al._ (2013) C. Siqueira Mello, M. Spite, B. Barbuy, F. Spite, E. Caffau, V. Hill, S. Wanajo, F. Primas, B. Plez, R. Cayrel, J. Andersen, B. Nordström, C. Sneden, T. C. Beers, P. Bonifacio, P. François, and P. Molaro, Astron. Astrophys. 550, A122 (2013), arXiv:1212.0211 [astro-ph.SR] .
  * Wanajo  _et al._ (2021) S. Wanajo, Y. Hirai, and N. Prantzos, Mon. Not. Roy. Astron. Soc. 505, 5862 (2021), arXiv:2106.03707 [astro-ph.GA] .
  * Hirai  _et al._ (2022) Y. Hirai, T. C. Beers, M. Chiba, W. Aoki, D. Shank, T. R. Saitoh, T. Okamoto, and J. Makino, arXiv e-prints , arXiv:2206.04060 (2022), arXiv:2206.04060 [astro-ph.GA] .
  * Goriely (1999) S. Goriely, Astron. Astrophys. 342, 881 (1999).
  * Tachibana  _et al._ (1990) T. Tachibana, M. Yamada, and Y. Yoshida, Progress of Theoretical Physics 84, 641 (1990).
  * Goriely  _et al._ (2010) S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. C 82, 035804 (2010), arXiv:1009.3840 [nucl-th] .


spacing=nonfrench

# Two regimes of confinement in photonic nanocavities: bulk confinement versus
lightning rods

Marcus Albrechtsen<EMAIL_ADDRESS>Department of Photonics Engineering,
DTU Fotonik, Technical University of Denmark, Building 343, DK-2800 Kgs.
Lyngby, Denmark  Babak Vosoughi Lahijani  Department of Photonics Engineering,
DTU Fotonik, Technical University of Denmark, Building 343, DK-2800 Kgs.
Lyngby, Denmark  NanoPhoton - Center for Nanophotonics, Technical University
of Denmark, Ørsteds Plads 345A, DK-2800 Kgs. Lyngby, Denmark  Søren Stobbe
Department of Photonics Engineering, DTU Fotonik, Technical University of
Denmark, Building 343, DK-2800 Kgs. Lyngby, Denmark  NanoPhoton - Center for
Nanophotonics, Technical University of Denmark, Ørsteds Plads 345A, DK-2800
Kgs. Lyngby, Denmark

(August 27, 2024)

###### Abstract

We present a theoretical study of dielectric bowtie cavities and show that
they are governed by two essentially different confinement regimes. The first
is confinement inside the bulk dielectric and the second is a local lightning-
rod regime where the field is locally enhanced at sharp corners and may yield
a vanishing mode volume without necessarily enhancing the mode inside the bulk
dielectric. We show that while the bulk regime is reminiscent of the
confinement in conventional nanocavities, the most commonly used definition of
the mode volume gauges in fact the lightning-rod effect when applied to ultra-
compact cavities, such as bowties. Distinguishing between these two regimes
will be crucial for future research on nanocavities, and our insights show how
to obtain strongly enhanced light-matter interaction over large bandwidths.

##  I Introduction

The confinement of light in optical nanocavities [1, 2, 3] enhances the
inherently weak interaction between light and matter [4]. This is important
for a wide span of applications ranging from optomechanics over lasers to
biosensing and emerging quantum technologies [5, 4, 6, 7]. Each application
has unique figures of merit [3, 8], and here we consider enhancement of the
decay rate in a nanostructure compared to a homogeneous medium, i.e., the
ratio of the local density of optical states (LDOS) to the density of optical
states, which is commonly denoted the Purcell factor [9],

$P=\frac{3}{4\pi^{2}}\left(\frac{\lambda}{n}\right)^{3}\frac{Q}{V}.$ (1)

Here $\lambda$ is the resonant wavelength of the nanocavity, $n$ is the
refractive index, and $V$ is effective mode volume. Notably, both $n$ and $V$
depend on position and must be evaluated at the position of the emitter,
$\mathbf{r}_{0}$. The quality factor, $Q$, describes the localization in
frequency space of the cavity mode due to resonant features. Most previous
research on dielectric cavities focused on enhancing the Purcell factor
through high quality factors [2, 10] since it was believed that the mode
volume was bounded at the diffraction limit [11, 12, 13],
$V_{\text{DL}}=[\lambda/(2n)]^{3}$. This has provided important advances in
light-matter interaction during the past decades, but high-$Q$ cavities are
unsuited for broadband interactions with pulses required for, e.g.,
electronics-photonics integration [14]. Additionally, it has become
increasingly difficult to further increase $Q$ in nanostructures due to
structural disorder and other ubiquitous surface effects [15, 16]. Plasmonic
cavities can realize deep subwavelength $V$, but are intrinsically limited by
ohmic losses to $Q<100$ [17]. These considerations call for research into
decreasing the mode volume of dielectric nanocavities.

Dielectric cavities with mode volumes below the diffraction limit can be
realized using bowtie structures. These were first identified by Gondarenko et
al. [18] who used an inverse-design algorithm to produce a structure from
which the authors intuited resonant features that define a single cavity mode
with a bowtie in the center and mode volumes below the diffraction limit. In
later work, Gondarenko et al. combined this concept with larger, high-$Q$
cavities to further enhance the $Q/V$-ratio [19]. Later works on bowtie
cavities [20, 21, 8, 22, 23, 24, 25, 26] expanded on these concepts. Our work
is motivated by this line of research where in several cases mode volumes
below the diffraction limit are claimed even for modes extending over regions
much larger than the $\simeq($200\text{\,}\mathrm{nm}$)^{3}$ volume
corresponding to the diffraction limit at telecom wavelengths in silicon. This
apparent yet so far overlooked inconsistency calls for a detailed study of the
mode volume and the mechanisms for photon confinement in nanocavities. This is
what we report here and we find that two related yet very different effects
govern the confinement: Lightning-rod effects, which appear due to localized
surface fields at sharp features, and bulk confinement, which confines light
inside the dielectric material. We present detailed numerical investigations
of the two effects, which point to several new insights that were so far
overlooked. For example, we observe that lightning-rod effects are highly
localized and barely change the global structure of the mode, which explains
how subdiffraction mode volumes can appear in spatially extended modes. We
further show that the radius of curvature plays a crucial role, which, if
neglected, can result in arbitrarily small but unphysical mode volumes
regardless of whether or not the mode is globally confined.

##  II Theory

Figure 1:  Working principles of a dielectric bowtie. (a), Boundary effects on
the electric field (red arrow), E, and the displacement field (blue arrow),
$\mathbf{D}=\epsilon\mathbf{E}$, across dielectric interfaces. Larger arrows
imply enhancement of the associated field and the magenta arrows illustrate
the continuity requirement of the potential around corners, causing lightning-
rod effects and field singularities if the corners are sharp. (b), Two
oppositely facing lightning rods form a bowtie (top) where the field is
enhanced between them. In contrast, the field decays rapidly away from a
single lightning rod (bottom). (c), Bowties defined by a gap,
$g=$20\text{\,}\mathrm{n}\mathrm{m}$$, constructed of two facing cones of
angle $\theta=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$, controlled by a
radius of curvature of the tip, $R=$20\text{\,}\mathrm{n}\mathrm{m}$$, which
enhances the field everywhere between the tips in silicon (air) in the top
(bottom) row. The central (last) column shows the electric field amplitude
(energy density) calculated as a 2D electrostatic problem using a finite-
element method. (d-e), Electric energy of a single cone and a bowtie,
respectively, each with different $R$. (f), Normalized energy density, $u$,
along the center of (d-e). All colormaps (c)-(e) are normalized to their
maximum.

A resonant mode has a certain distribution of the electric field,
$\mathbf{E}$, with an associated energy density,
$u_{0}=\mathbf{D}\cdot\mathbf{E}=\epsilon_{0}\epsilon_{r}\mathbf{E}\cdot\mathbf{E}$,
and the total energy in the mode is $\int u_{0}\text{d}\mathbf{r}^{\prime}$.
This defines a normalized energy density with dimensionality of inverse
volume,

$u(\mathbf{r})=\frac{\epsilon_{r}(\mathbf{r})\mathbf{E}(\mathbf{r})\cdot\mathbf{E}(\mathbf{r})}{\int_{\mathbf{r^{\prime}}}\epsilon_{r}(\mathbf{r}^{\prime})\mathbf{E}(\mathbf{r}^{\prime})\cdot\mathbf{E}(\mathbf{r}^{\prime})\text{d}\mathbf{r^{\prime}}},$
(2)

where $\epsilon_{r}=n^{2}$ is the dielectric constant, and here we consider
silicon with $n=3.48$. In general, $u$ is complex but for high-$Q$ cavities it
is an excellent approximation to assume that $u$ is purely real [27]. To
obtain the strongest possible enhancement of spontaneous emission from a
dipole emitter in a cavity, the dipole should be placed at the most intense
point of the mode in the cavity, which defines the effective mode volume as
the inverse of the normalized energy density at that position [27, 4, 3],
$V(\mathbf{r}_{0})=1/\text{Re}\\{u(\mathbf{r}_{0})\\}$. More generally, the
emitter should be placed in the position with the highest LDOS, which is a
function of frequency, polarization, and position.

We note that besides enhancement of spontaneous emission, which is quantified
by the Purcell factor, cf. Eq. (1), other applications of light-matter
interactions such as two-photon absorption and Kerr nonlinearities are
governed by different mode volumes with, e.g., different powers of the field
or the volume integral only over regions with nonlinear materials [3]. For the
remainder of this work, we limit the discussion to the Purcell factor and how
its mode volume can be reduced using the boundary conditions of Maxwell’s
equations.

##  III Electrostatic analysis

###  III.1 Bowties versus lightning rods

Before proceeding with full three-dimensional (3D) electrodynamic
calculations, we first consider a two-dimensional (2D) finite-element
electrostatic model to gain insight into the scaling of the energy density,
and thus mode volume, with cavity dimensions and geometry. The electrostatic
model is a valid approximation to the full electrodynamic problem because we
consider deep-subwavelength features [13, 8]. The relevant figure of merit in
this model is the mode area, $A(\mathbf{r}_{0})$, which appears when the
integral over $\text{d}\mathbf{r}^{\prime}$ is evaluated in 2D [10].

Figure 1(a) shows the boundary conditions across a dielectric material
interface, which demand that the tangential component of the electric field,
E, and the normal component of the displacement field,
$\textbf{D}=\epsilon\textbf{E}$, are continuous. The opposite configuration
leads to discontinuous fields that differ by a factor of $\epsilon_{r}$. The
magenta arrows indicate geometrically sharp vertices at which the electric
field must diverge locally. Figure 1(b) shows two facing tips that are close
to each other and form a bowtie [18, 19, 20, 21, 8, 22, 23, 24, 25, 26], thus
enhancing the field between them. A single tip results in a locally enhanced
field that falls off rapidly away from it. The field divergences at sharp tips
imply lightning-rod singularities to ensure continuity of the electric
potential, which is required by the electromagnetic continuum model of
materials [28, 29, 30, 31, 32]. The potential remains continuous so the energy
is well defined. However, the field, and thus the energy density, locally
diverges at such points, which implies that local evaluations of the mode
volume based on the point of maximum $u$, is prone to gauge these divergent
fields. Since perfectly sharp corners are unphysical, so is the associated
mode volume.

Figure 2:  Scaling of the mode area with the radius of curvature, $R$, gap,
$g$, and bowtie angle, $\theta$. (a-b), Scaling of $A$ (in
${\mathrm{\SIUnitSymbolMicro m}}^{2}$) evaluated at the maximum (at material
boundary), $A_{\text{max}}$, and at the geometric center, $A_{\text{cen}}$,
respectively. (c), Linear ratio between (a) and (b),
$A_{\text{max}}/A_{\text{cen}}$. (d), Mode area $A_{\text{max}}$ (black line)
and $A_{\text{cen}}$ (magenta line) for
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$ bowtie angles as a function of $R$,
i.e., horizontal cuts through the center of (a-b). Wider angles (blue lines)
correspond to the geometry approaching a slot, which does not have a
lightning-rod regime but also cannot reduce $A_{\text{cen}}$ as much as a
sharper bowtie. More narrow angles (red lines) exhibit more pronounced
differences between the two confinement regimes and has minimal
$A_{\text{cen}}$ for $R>g$ depending on the angle. The smallest mode area is
obtained for $\theta=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$ with $R=g$.

In finite-element methods and other meshed numerical models, the calculated
mode volume will depend on the size of the mesh and will not converge.
Instead, to ensure converging results that describe the system instead of the
mesh, a finite sharpness must be imposed on these lightning-rods, i.e., using
a finite radius of curvature, $R$, as shown in Fig. 1(c). The mesh must
obviously be much smaller than the radius of curvature it aims to resolve. The
angle, $\theta$, is $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ everywhere in
this work unless specified. Here the field is enhanced between the tips and in
the limit where $R\rightarrow\infty$, the bowtie reduces to a slot [33, 34,
35], and when $R\rightarrow 0$, the bowtie becomes two facing lightning rods
with diverging fields. Figure 1(c) further compares bowties with solid (void)
centers in the top (bottom) row, where the central column shows $|\mathbf{E}|$
and the last column shows the energy density, $u$. Notably, while the electric
field appears more localized in the air center compared to the dielectric
center, the energy distributions of both bowtie types are identical. For this
reason, we will denote the central part between the bowtie tips as a gap,
regardless of whether the refractive index of the gap is smaller or larger
than that of the tips. For example, our considerations apply to both
dielectric gaps (i.e., bridges) surrounded by air tips and an air gap between
dielectric tips; only it should be noted that the polarization of the electric
field is rotated by $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ between the two
cases due to the different boundary conditions [8]. These insights show that
the radius of curvature plays a decisive role for nanocavities based on
boundary effects, which we explore in further detail below.

Figure 1(d-e) show $u$ for a lightning-rod and a bowtie respectively, with a
gap, $g=$10\text{\,}\mathrm{n}\mathrm{m}$$, and a radius of curvature of the
tip in the top (bottom), $R=$10\text{\,}\mathrm{n}\mathrm{m}$$
($R=$1\text{\,}\mathrm{n}\mathrm{m}$$). Figure 1(f) shows the normalized
energy density, $u$, evaluated at different positions along the cavity center
for the four structures shown in Figs 1(d-e). This figure illustrates a main
insight of our work: While the field can be enhanced locally by both lightning
rods and bowties, the field maximum is larger but falls off much more rapidly
for a lightning rod compared to a bowtie with the same radius of curvature.
This means that although a lightning rod creates a region of high field
intensity, the field is mostly enhanced at the surface. In addition, the
bowtie effect is stronger because the long tail of the lightning rod that
would have appeared if the bowtie had had only one boundary is suppressed by
the presence of the oppositely facing boundary, and the energy density is
normalized to the total energy in the system, cf. Eq. 2. Another important
insight from this figure is that while reducing the radius of curvature for a
bowtie enhances the maximum field, which occurs at the material interface, the
energy density $u$ in the bulk can be reduced when $R$ decreases. The blue-
shaded region highlights the range where this reduction in the light-matter
interaction occurs.

The sharpest possible tip induces the smallest mode volume at the interface,
which is relevant for light-matter interaction at the very surface. However,
this is seldomly desirable due to, e.g., nonradiative surface recombination.
In practice, such surface effects are extremely sensitive to, or even
indistinguishable from, structural disorder, and besides the extreme care
needed in design, nanofabrication, and numerical calculations to approach such
effects, they would also raise more fundamental questions about the validity
of the electromagnetic continuum model of materials, local field effects, and
other complex surface phenomena such as oxidation and condensation [16, 36].
Approaching this in a consistent way would be a daunting research challenge
that has not been attempted so far. To summarize the findings of the preceding
paragraphs: It has no physical meaning to discuss bowtie or lightning-rod
cavities without specifying the radius/radii of curvature. This is in part
because a radius of curvature of zero represents an unphysical extrapolation
of continuum theory and in part because the mode volume goes to zero when the
radius of curvature goes to zero. In addition, the smallest possible mode
volume in the center of a bowtie cavity, which is typically the relevant
position for experiments and technology, is not obtained for the smallest
possible radius of curvature but rather at some optimum value as discussed in
further detail below.

###  III.2 Scaling of bowties structures

The findings of the preceding section call for a systematic study of how the
mode area depends on the gap and the radius of curvature. Figure 2(a-b) shows
the mode area (in ${\mathrm{\SIUnitSymbolMicro m}}^{2}$) evaluated at the
point of maximum energy density (i.e., at the material interface),
$A_{\text{max}}$, and at the center of the bowtie, $A_{\text{cen}}$, as a
function of $R$ and $g$. We observe that reducing the gap always reduces both
mode areas, implying an enhanced light-matter interaction. In contrast,
reducing the radius of curvature always reduces $A_{\text{max}}$ but not
always $A_{\text{cen}}$. This is shown more directly in Fig. 2(c) where the
ratio $A_{\text{max}}/A_{\text{cen}}$ is shown. The diagonal contours along
$g=R$ show that the two evaluations can give significantly different results
for some parameters while they are equivalent for large $R$ and small $g$. We
note that these observations are consistent with Fig. 1(f). Figure 2(d) shows
the two mode areas as a function of $R$ for a gap of
$g=$1\text{\,}\mathrm{n}\mathrm{m}$$, i.e., the figure shows a horizontal scan
through Fig. 2(a-b). Additionally, $A_{\text{cen}}$ is shown for multiple
other bowtie angles with wider angles corresponding to a more slot-like system
and where narrow angles results in lightning rods for larger $R$ compared to
the angle, $\theta=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$. Notably,
$A_{\text{cen}}$ exhibits a minimum at $R_{\text{min}}=g$ for
$\theta=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$, indicating optimum
parameters for bowties confining light inside the material under the
electrostatic approximation. We further observe that $A_{\text{cen}}$ rapidly
increases for $R<R_{\text{min}}$ for narrow-angled bowties, while
$A_{\text{cen}}$ is substantially unaffected for the wide-angled bowties and
instead plateaus at a larger mode area compared to the optimal case of
$\theta=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$. This is consistent with
the scaling of a slot, which does not need a finite radius of curvature but
also cannot achieve as small mode areas as bowtie structures.

We denote the region $R<R_{\text{min}}$ the lightning-rod regime of bowtie
cavities, with $R_{\text{min}}=g$ in the optimal case of
$\theta=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$: Cavities aiming for
interactions at the surface could aim to make the sharpest possible tips to
enhance the interaction, while cavities relying on emitters inside the
material should aim for $R\sim R_{\text{min}}$. These effects can also be
relevant for conventional nanocavities, such as H0 cavities, but it is not
pronounced in cavities with a global antinode in the very center, i.e., inside
the bulk, such as the fundamental mode of an L3 cavity. Notably, it can be
important for higher-order modes, which may seem to have a smaller volume than
the fundamental mode due to field discontinuities at interfaces [37].

Cavities approaching the lightning-rod regime must carefully consider the
position at which their figures of merit are computed, such as the mode
volume. We observe that lightning rods cause local spikes in the field that do
not substantially change the volume integral, and, thus, the overall energy
distribution within the mode. In turn, they depend on the local polarization
of the field of the dominating cavity mode. In other words, lightning-rod
effects are local effects that may imply locally enhanced fields but they do
not imply that the mode is globally confined to a small volume. For example, a
sharp dielectric defect on the surface of a waveguide, or even just any corner
of a rectangular waveguide, would imply a mode volume that can be arbitrarily
small (given by the radius of curvature).

We denote the region $R\geq R_{\text{min}}$ the bulk confinement regime of
bowtie cavities in which the field falls off slowly within the gap, which
causes an enhanced field throughout the bulk material and not just at the
boundaries. This substantially affects $u$, since the volume integral
normalizes it, as well as the other cavity figures of merit, i.e., evaluation
in the center of the bowtie yields a value corresponding to confinement in the
bulk. To summarize these findings: The mode volume evaluated at the center of
a bowtie cavity describes confinement in the bulk and resembles the concept of
the mode volume commonly discussed in the literature for other types of
nanocavities. In contrast, the mode volume evaluated a the field maximum for
bowtie cavities describes a local lightning-rod effect at material boundaries.
This can be equivalent to $V$ inside the structure if $R\geq R_{\text{min}}$
but otherwise the two definitions are not comparable.

##  IV Electrodynamic analysis

###  IV.1 Lightning-rod effects in conventional cavities

Figure 3:  Mode volumes below the diffraction limit in an L3 cavity due to
lightning-rod effects at a cone far away from the cavity center. (a),
Illustration of an L3 cavity where 3 holes are removed along a line in the
photonic crystal. The blue square shows the position of a defect
($\mathbf{r}_{d}=(y=a\sqrt{3}-R,x=z=0)$), enlarged in (b). (c), The cavity
mode calculated without the defect, and (d) with the defect limited by a
radius of curvature, $R=$10\text{\,}\mathrm{p}\mathrm{m}$$. The inset in (d)
shows the mode around the tip of the defect. Note that the mode profile is
identical to (c) everywhere except around the defect but appears different as
the colormap is normalized to the field at $\mathbf{r}_{d}$. (e), The resonant
wavelength and quality factor for the structure without a defect (solid line)
and with a defect (circles) for a number defect radius of curvatures, $R$.
Evidently, the defect has no significant effect. (f), The mode volume
evaluated in the center, $V(\mathbf{r}_{c})$, and at the defect tip,
$V(\mathbf{r}_{d})$. The red line is a power-law fit $V(\textbf{r}_{d},R)$
showing an almost linear dependence on $R$ of the defect in this particular
case. (g), The norm of the electric field, $|\mathbf{E}(y,x=z=0)|$. The black
outline is for a defect with $R=$10\text{\,}\mathrm{p}\mathrm{m}$$ and the
blue shaded area is the mode without a defect. They are identical everywhere
except at $\mathbf{r}_{d}$ where the field has a local spike.

While 2D electrostatics can provide insights into the scaling of bowties, we
will henceforth consider 3D electrodynamic theory to illustrate the impact of
lightning-rod effects in cavities studied in the literature. We first
investigate an L3 photonic-crystal cavity [2] (with no shifted holes) shown in
Fig. 3(a). We consider a pitch of $a=$420\text{\,}\mathrm{n}\mathrm{m}$$, a
hole radius of $b=$120\text{\,}\mathrm{n}\mathrm{m}$$, and a membrane
thickness of $t=$250\text{\,}\mathrm{n}\mathrm{m}$$. The photonic crystal has
9 periods along $x$ and 7 periods along $y$. In addition to the conventional
design, we add a nanometer-scale defect to the simulation domain, acting as a
lightning-rod. Figure 3(b) shows the defect introduced in the side of the
cavity shaped as a cone of silicon with a height and base radius
$h=$1\text{\,}\mathrm{n}\mathrm{m}$$. The very tip of the defect is rounded to
a sphere of radius $R=$10\text{\,}\mathrm{p}\mathrm{m}$$. Figure 3(c-d) shows
the eigenmode of the structure without and with the defect, both calculated in
3D using a finite-element method. The mesh at the apex of the defect is much
smaller than $R$ such that the model converges with finer mesh, and instead
the divergence is controlled directly by $R$. The effect of the tip is
directly observed from the normalized field plots (Fig. 3(c) and (d)).

To show explicitly that the perturbation by small defects does not affect the
mode at the center of the L3 cavity, we first ensured convergence of the
unperturbed structure and then kept this mesh constant in the structure
$5\text{\,}\mathrm{n}\mathrm{m}$ away from the defect to avoid numerical
fluctuations and only show the effects of the defect. The small defect acts as
a lightning rod, enhancing the field locally at the tip in a vanishing volume,
controlled by the finite radius of curvature [29, 30]. This has a negligible
effect on the cavity figures of merit as shown in Fig. 3(e-f), i.e., the
resonant wavelength, the quality factor, and the mode volume at the center of
the L3 cavity. This is because the divergent field happens in a vanishing
volume, and therefore, the contribution from such a sharp lightning rod to the
volume integral fades with mesh refinement. Defects may in general reduce the
quality factor but this is not observed due to the modest quality factors of
the L3 cavity studied here. Figure 3(g) shows a line-scan through the cavity
where the defect is visible as a local enhancement of the electric field on
top of the global field intensity of the resonant structure. Notice that the
smallest mode volume is found in air, not silicon, since the defect is a
protrusion. It is possible to make a similar defect, which achieves a deep-
subwavelength $V$ in the high-index material instead by simply observing the
boundary conditions and the local polarization of the field of the L3 cavity
mode [29, 30] to make a reentrant defect instead. For example, a reentrant
defect at the first hole along $(x,y=z=0)$ results in $V<(\lambda/(2n))^{3}$
for a radius of curvature $R\leq$1\text{\,}\mathrm{p}\mathrm{m}$$. This study
verifies that the mode volume calculated by the maximum evaluation can go to
zero when the domain contains sharp features without affecting the global
structure of the cavity mode. In practice, however, such features are
indiscernible from structural disorder and numerically the divergence will be
resolved by the finite mesh size. Importantly, the addition of such features
does not change the light-matter interaction strength elsewhere in the cavity
including the center.

###  IV.2 Bulk confinement versus lightning rods in bowtie cavities

Figure 4:  Model of a nanocavity with an array of V-groove bowties that have a
gap, $g=$4\text{\,}\mathrm{n}\mathrm{m}$$, where lightning-rod effects imply
mode volumes deep below the diffraction limit at the V-groove tip and bulk
confinement slightly below the diffraction limit. (a), The cavity is a
suspended photonic-crystal nanobeam-cavity defined in silicon. The inset shows
a V-groove penetrating half way into the bowtie that creates a lightning-rod
effect. (b-c), Electric field norm, $|\mathbf{E}|$, and effective mode volume,
$V$, respectively, in the ($xy$) cross-section plane, shown at the center of
the structure, close to the interface, and at the interface. (d-e), Side-view
of $V$ in $yz$\- and $xz$-planes with insets showing the apex of the
lightning-rod. (f-g), Line-scans through the structure showing $1/V$
normalized to the inverse of the diffraction limit; the gray-shaded area
highlights the regions with a mode volume below the diffraction limit. The
cavity is solved as an eigenfrequency problem in 3D with a finite element
method using $\sim$900\text{\,}\mathrm{G}\mathrm{B}$$ RAM to resolve the many
narrow features with high aspect ratios with a mesh size well below $g$.
Figure 5:  Model of a nanocavity with an array of V-groove bowties that have a
gap, $g=$12\text{\,}\mathrm{n}\mathrm{m}$$, where lightning-rod effects imply
mode volumes below the diffraction limit at the V-groove tip and bulk
confinement above the diffraction limit. (a), The cavity is a suspended
photonic-crystal nanobeam-cavity defined in silicon. The inset shows a
V-groove penetrating half way into the bowtie that creates a lightning-rod
effect. (b-c), Electric field norm, $|\mathbf{E}|$, and effective mode volume,
$V$, respectively, in the ($xy$) cross-section plane, shown at the center of
the structure, close to the interface, and at the interface. (d-e), Side-view
of $V$ in $yz$\- and $xz$-planes with insets showing the apex of the
lightning-rod. (f-g), Line-scans through the structure showing $1/V$
normalized to the inverse of the diffraction limit; the gray-shaded area
highlights the regions with a mode volume below the diffraction limit.

The example in Fig. 3 is hypothetical, constructed to illustrate the
physically different regimes associated with evaluating mode volumes in the
geometric center of nanocavities compared to at the points with the highest
field intensity. However, our insights have significant impact also on current
experimental research. We now consider the recent experiments on dielectric
nanocavities by Hu et al. [23], where a photonic-crystal nanobeam cavity with
an array of V-groove bowtie unit cells was considered. As we will show, the
interplay between the bowtie effect and the lightning rod due to the V-groove
has decisive impact on the numerical challenges as well as on the conclusions
that can be drawn from the experiment. Due to significant differences between
the designed and the actually fabricated cavities in this experiment, we
consider both independently. Starting with the designed structure, Fig. 4(a)
shows the nanobeam cavity design consisting of a central bowtie, 20 tapered
unit cells on both sides, and 10 additional bowties of the largest tapered
bowtie forming a mirror to obtain a high quality factor; the bowtie has a
designed gap, $g=$4\text{\,}\mathrm{n}\mathrm{m}$$. The nanobeam is
$220\text{\,}\mathrm{n}\mathrm{m}$ thick and
$700\text{\,}\mathrm{n}\mathrm{m}$ wide. The inset shows the central unit cell
where a V-groove is introduced to penetrate half-way
($110\text{\,}\mathrm{n}\mathrm{m}$) through the membrane. This causes a shift
of confinement regime from bulk to lightning-rod as discussed further below.

We estimate an in-plane radius of curvature,
$R_{xy}=$30\text{\,}\mathrm{n}\mathrm{m}$$, from the geometry used in the
theoretical model by Hu et al. We find the fundamental mode of the cavity
using a finite-element method solving a 3D electrodynamic problem. Figure 4(b)
shows the norm of the electric field, $|\mathbf{E}|$, mid-way between the tip
of the V-groove and the bottom of the structure, i.e., at
$z=$55\text{\,}\mathrm{n}\mathrm{m}$$ for the five central unit cells. The
mode exhibits an array of hotspots at each bowtie, modulated by an envelope
along the tapered unit cells. Since $R_{xy}>g$ and
$\theta=$120\text{\,}\mathrm{\SIUnitSymbolDegree}$$ for the in-plane geometry,
the confinement would have been in the bulk regime, cf. Figs. 2(c) and (d), in
the absence of the V-groove. Figure 4(c) shows the normalized energy density,
proportional to the inverse mode volume, at different distances from the
V-groove tip of the central bowtie in Fig. 4(b), with
$z=$55\text{\,}\mathrm{n}\mathrm{m}$$ corresponding to the middle of the
structure between the V-groove and the bottom of the cavity shown in Fig. 4(b)
and $z=$110\text{\,}\mathrm{n}\mathrm{m}$$ corresponding to the tip of the
V-groove. These colormaps are normalized to the minimum mode volume of the
structure, which occurs at the very tip of the V-groove at
$z=$110\text{\,}\mathrm{n}\mathrm{m}$$, and show that the mode is visibly much
less enhanced in the bulk only $10\text{\,}\mathrm{n}\mathrm{m}$ away from the
tip.

Figure 4(d-e) shows $1/V$ at two planes cut through the center of the
structure in the $yz$\- and $xz$-planes, respectively, with insets showing the
lightning-rod tip caused by the small radii of curvature estimated from the
theoretical geometry, $R_{yz}=$2\text{\,}\mathrm{n}\mathrm{m}$$ and
$R_{xz}=g/2=$2\text{\,}\mathrm{n}\mathrm{m}$$. The mode is enhanced in the
dielectric material just at the interface and drops rapidly away from the
boundary, similarly to the effect observed in Fig. 1(d). Figure 4(f) shows
$1/V$ across the center at $z=$55\text{\,}\mathrm{n}\mathrm{m}$$, i.e., a
horizontal line-scan through the center of Fig 4(b). This confirms that the
symmetric bowtie effect in the $xy$-plane leads to a slowly varying field
inside the material and thus a bulk confinement regime in the absence of the
V-groove as already predicted from the geometric parameters. In contrast,
Figure 4(g) shows $1/V$ along $z$, i.e., vertically through the centers of
Figure 4(d) and (e), where the V-groove acts as a lightning rod enhancing the
mode volume at the surface, because both $R_{xz}$ and $R_{yz}$ are much
smaller than the height of the structure, as shown conceptually in Figs.
1(d)-(f). The shaded areas in Figs. 4(f) and (g) show regions where the mode
volume is below the diffraction limit, and it can be readily seen from Fig.
4(g) that the mode volume increases from $<(1/10)[\lambda/(2n)]^{3}$ to
$\sim(1/2)[\lambda/(2n)]^{3}$ from the lightning-rod tip and into the bulk,
i.e., almost an order of magnitude difference. Even though the mode volume
stems from the lightning-rod effect, it remains below the diffraction limit
inside the structure due to the small designed gap,
$g=$4\text{\,}\mathrm{n}\mathrm{m}$$. More specifically, we calculate the mode
volume of this theoretical structure at the tip of the V-groove,
$\mathbf{r}_{\text{min}}=(x=0,y=0,z=110)\ \text{nm}$, to be
$V_{\text{tip}}^{\text{t}}=0.073[\lambda/(2n)]^{3}$. This is in good agreement
with the claimed value, $V_{\text{min}}^{\text{ref}}=0.067[\lambda/(2n)]^{3}$.
However, the mode volume in the center, i.e.,
$\mathbf{r}_{\text{c}}=(x=0,y=0,z=55)\ \text{nm}$, is
$V_{\text{cen}}^{\text{t}}=0.51[\lambda/(2n)]^{3}$ which is much larger but
still below the diffraction limit.

We now consider the same cavity as in Fig. 4 but for the dimensions which Hu
et al. measured on the fabricated structure. The result is shown in Fig. 5.
All components are identical to Fig. 4, except the bowtie gap, which is
changed to the value measured by scanning electron microscopy,
$g=$12\text{\,}\mathrm{n}\mathrm{m}$$, and we estimate the following radii of
curvature from the scanning electron micrographs provided by Hu et al.:
$R_{xy}=$30\text{\,}\mathrm{n}\mathrm{m}$$,
$R_{yz}=$15\text{\,}\mathrm{n}\mathrm{m}$$, and
$R_{xz}=g/2=$6\text{\,}\mathrm{n}\mathrm{m}$$, which we use to build the
numerical model shown in Fig. 5(a). The resulting mode is shown in Figs.
5(b-g) and it is qualitatively similar to that found in Fig. 4 although with
important differences. In particular, the minimum mode volume of this
structure, $V_{\text{min}}^{\text{SEM}}=0.36[\lambda/(2n)]^{3}$, which occurs
at the tip of the fabricated V-groove as shown in Figs. 5(d)-(g), is $\sim 6$
times larger than for the designed structure due to the larger gap and,
importantly, the larger radii of curvature. We find the mode volume at the
center of the structure, i.e., the mode volume associated with bulk
confinement, $V_{\text{cen}}^{\text{ref}}=1.2[\lambda/(2n)]^{3}$, which is
above the diffraction limit and comparable to conventional nanocavities such
as H1 photonic-crystal cavities [37] with $V\sim 2.6[\lambda/(2n)]^{3}$. The
larger radii of curvature also imply that the field falls off slower compared
to the designed structure, which for example can be seen in Figs. 4(d), (e),
and (g). Again, the shaded areas highlights the regions where the mode volume
is below the diffraction limit. Our estimates of the radii of curvature are
associated with substantial uncertainties but we believe they are conservative
and therefore provide lower bounds to the calculated mode volumes.

##  V Discussion and conclusion

In summary, we have shown that the ratio of the gap to the radius of curvature
divides bowtie cavities into two confinement regimes: bulk confinement versus
lightning rods at surfaces. An important conclusion is that while smaller
gaps, $g$, always lead to smaller mode volumes, there is an optimum value of
the radius of curvature, which for the optimal bowtie angle of
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$ is $R_{\text{min}}=g$, and it is not
advantageous to decrease the radius of curvature below this threshold if the
goal is enhancing the light-matter interaction inside the material. For a
radius of curvature below $R_{\text{min}}$, the bowtie operates instead as two
opposing lightning rods. This distinction is crucial for research on
semiconductor devices such as nanolasers and nonlinear optical devices as well
as on electronic-photonic integration, which rely on field enhancements inside
the dielectric material. Lightning rods may have interesting applications for
enhancing surface interactions, but they do not necessarily imply a strong
light-matter interaction inside the material. Since defects or sharp corners
on the surface of any dielectric structure can lead to lightning-rod effects,
one might wonder why they are not dominating experiments on transmission in
waveguides, near-field spectroscopy, and essentially all optical experiments.
The answer may be that these surface effects are so localized and broadband
that they average out in actual experiments. In this context, it is
interesting to note that strongly localized fields also appear at sharp
metallic tips, but more advanced hydrodynamic models of the carriers have
shown that such surface effects arise from extrapolation of continuum theory
beyond its validity [38].

We presented detailed numerical models of nanocavities considered in recent
experiments to demonstrate the relevance of our new insights. Since the mode
volume can only be measured indirectly through measurements of the LDOS [39,
4], it may be more useful to infer the mode volume from the fabricated
geometry. This means in turn that it becomes critical to accurately estimate
the radius of curvature of the fabricated structures, such that theoretical
predictions of experiments remain accurate and optimized designs perform as
intended. The need for careful measurements of structural parameters is of
course not new in nanocavity research but there is a main difference between
conventional high-$Q$ nanocavities and cavities with mode volumes below the
diffraction limit: In the former case, the most important geometry parameters
govern structural disorder, which has a decisive impact on $Q$ while leaving
the mode volume approximately invariant. For bowtie cavities, however, the gap
and the radius of curvature are also essential fabrication constraints and
they have a decisive impact on $V$. The interplay between disorder, gap, and
radius of curvature and their impact on $Q$ has not been studied and would be
an important line of future research. Ideally, the fabrication constraints
should be included already at the design stage [22], which would require an
unprecedented integration of nanofabrication, design, and optical experiments.

## Acknowledgments

We gratefully acknowledge financial support from the Villum Foundation Young
Investigator Program (Grant No. 13170), the Danish National Research
Foundation (Grant No. DNRF147 - NanoPhoton), and Innovation Fund Denmark
(Grant No. 0175-00022 - NEXUS).

## References

  * Painter  _et al._ [1999] O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, Science 284, 1819 (1999).
  * Akahane  _et al._ [2003] Y. Akahane, T. Asano, B.-S. Song, and S. Noda, Nature 425, 944 (2003).
  * Notomi [2010] M. Notomi, Rep. Prog. Phys. 73, 096501 (2010).
  * Lodahl  _et al._ [2015] P. Lodahl, S. Mahmoodian, and S. Stobbe, Rev. Mod. Phys. 87, 347 (2015).
  * Vahala [2003] K. J. Vahala, Nature 424, 839 (2003).
  * Koenderink  _et al._ [2015] A. F. Koenderink, A. Alu, and A. Polman, Science 348, 516 (2015).
  * Ma  _et al._ [2020] Y. Ma, B. Dong, and C. Lee, Nano Converg. 7, 12 (2020).
  * Choi  _et al._ [2017] H. Choi, M. Heuck, and D. Englund, Phys. Rev. Lett. 118, 223605 (2017).
  * Purcell [1946] E. M. Purcell, Phys. Rev. 69, 681 (1946).
  * Kippenberg  _et al._ [2004] T. J. Kippenberg, S. M. Spillane, and K. J. Vahala, Appl. Phys. Lett. 85, 6113 (2004).
  * Coccioli  _et al._ [1998] R. Coccioli, M. Boroditsky, K. Kim, Y. Rahmat-Samii, and E. Yablonovitch, IEE Proc., Optoelectron. 145, 391 (1998).
  * Deotare  _et al._ [2009] P. B. Deotare, M. W. McCutcheon, I. W. Frank, M. Khan, and M. Lončar, Appl. Phys. Lett. 94, 121106 (2009).
  * Khurgin [2015] J. B. Khurgin, Nat. Nanotechnol. 10, 2 (2015).
  * Mork and Yvind [2020] J. Mork and K. Yvind, Optica 7, 1641 (2020).
  * Minkov and Savona [2014] M. Minkov and V. Savona, Sci. Rep. 4, 5124 (2014).
  * Sekoguchi  _et al._ [2014] H. Sekoguchi, Y. Takahashi, T. Asano, and S. Noda, Opt. Express 22, 916 (2014).
  * Wang and Shen [2006] F. Wang and Y. R. Shen, Phys. Rev. Lett. 97, 206806 (2006).
  * Gondarenko  _et al._ [2006] A. Gondarenko, S. Preble, J. Robinson, L. Chen, H. Lipson, and M. Lipson, Phys. Rev. Lett. 96, 143904 (2006).
  * Gondarenko and Lipson [2008] A. Gondarenko and M. Lipson, Opt. Express 16, 17689 (2008).
  * Liang and Johnson [2013] X. Liang and S. G. Johnson, Opt. Express 21, 30812 (2013).
  * Hu and Weiss [2016] S. Hu and S. M. Weiss, ACS Photonics 3, 1647 (2016).
  * Wang  _et al._ [2018] F. Wang, R. E. Christiansen, Y. Yu, J. Mørk, and O. Sigmund, Appl. Phys. Lett. 113, 241101 (2018).
  * Hu  _et al._ [2018] S. Hu, M. Khater, R. Salas-Montiel, E. Kratschmer, S. Engelmann, W. M. J. Green, and S. M. Weiss, Sci. Adv. 4, eaat2355 (2018).
  * Mignuzzi  _et al._ [2019] S. Mignuzzi, S. Vezzoli, S. A. R. Horsley, W. L. Barnes, S. A. Maier, and R. Sapienza, Nano Lett. 19, 1613 (2019).
  * Zhou  _et al._ [2019] J. Zhou, J. Zheng, Z. Fang, P. Xu, and A. Majumdar, Opt. Express 27, 30692 (2019).
  * Zhao  _et al._ [2020] Q. Zhao, L. Zhang, and O. D. Miller, arxiv:2008.13241 (2020).
  * Kristensen  _et al._ [2012] P. T. Kristensen, C. Van Vlack, and S. Hughes, Opt. Lett. 37, 1649 (2012).
  * Sommerfeld [2004] A. Sommerfeld,  _Mathematical Theory of Diffraction_, 1st ed. (Birkhäuser Boston, Boston, MA, 2004) translated by R. J. Nagem, M. Zampolli, G. Sandri. (Original work published Math. Ann. 47, 317–374 (1896).).
  * Andersen and Solodukhov [1978] J. Andersen and V. Solodukhov, IEEE Trans. Antennas Propag. 26, 598 (1978).
  * Van Bladel [1985] J. Van Bladel, IEEE Trans. Antennas Propag. 33, 893 (1985).
  * Landau  _et al._ [1984] L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii,  _Electrodynamics of Continuous Media_ , 2nd ed. (Pergamon, Oxford, 1984) translated by J. B. Sykes, J. S. Bell, M. J. Kearsley.
  * Jackson [1999] J. D. Jackson,  _Classical Electrodynamics_ , 3rd ed. (Wiley, N.J., 1999).
  * Almeida  _et al._ [2004] V. R. Almeida, Q. Xu, C. A. Barrios, and M. Lipson, Opt. Lett. 29, 1209 (2004).
  * Robinson  _et al._ [2005] J. T. Robinson, C. Manolatou, L. Chen, and M. Lipson, Phys. Rev. Lett. 95, 143901 (2005).
  * Schneider and Seidler [2016] K. Schneider and P. Seidler, Opt. Express 24, 13850 (2016).
  * Asano  _et al._ [2017] T. Asano, Y. Ochi, Y. Takahashi, K. Kishimoto, and S. Noda, Opt. Express 25, 1769 (2017).
  * Nakayama  _et al._ [2011] S. Nakayama, S. Ishida, S. Iwamoto, and Y. Arakawa, Appl. Phys. Lett. 98, 171102 (2011).
  * Mortensen  _et al._ [2014] N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, and S. I. Bozhevolnyi, Nat. Commun. 5, 3809 (2014).
  * García de Abajo [2010] F. J. García de Abajo, Rev. Mod. Phys. 82, 209 (2010).