#  Quantum-informed plasmonics for strong coupling: the role of electron
spill-out

Ida Juliane Bundgaard  POLIMA—Center for Polariton-driven Light-Matter
Interactions, University of Southern Denmark, Campusvej 55, DK-5230 Odense M,
Denmark  Christian Nicolaisen Hansen  POLIMA—Center for Polariton-driven
Light-Matter Interactions, University of Southern Denmark, Campusvej 55,
DK-5230 Odense M, Denmark  P. Elli Stamatopoulou  0000-0001-9121-911X
POLIMA—Center for Polariton-driven Light-Matter Interactions, University of
Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark  Christos Tserkezis
0000-0002-2075-9036<EMAIL_ADDRESS>POLIMA—Center for Polariton-driven Light-
Matter Interactions, University of Southern Denmark, Campusvej 55, DK-5230
Odense M, Denmark

###### Abstract

The effect of nonlocality on the optical response of metals lies at the
forefront of research in nanoscale physics and, in particular, quantum
plasmonics. In alkali metals, nonlocality manifests predominantly as electron
density spill-out at the metal boundary, and as surface-enabled Landau
damping. For an accurate description of plasmonic modes, these effects need be
taken into account in the theoretical modelling of the material. The resulting
modal frequency shifts and broadening become particularly relevant when
dealing with the strong interaction between plasmons and excitons, where
hybrid modes emerge and the way they are affected can reflect modifications of
the coupling strength. Both nonlocal phenomena can be incorporated in the
classical local theory by applying a surface-response formalism embodied by
the Feibelman parameters. Here, we implement surface-response corrections in
Mie theory to study the optical response of spherical plasmonic–excitonic
composites in core–shell configurations. We investigate sodium, a jellium
metal dominated by spill-out, for which it has been anticipated that nonlocal
corrections should lead to an observable change in the coupling strength,
appearing as a modification of the width of the mode splitting. We show that,
contrary to expectations, the influence of nonlocality on the anticrossing is
minimal, thus validating the accuracy of the local response approximation in
strong-coupling photonics.

##  I Introduction

Plasmonics deals with the collective oscillations of free electrons in the
bulk or on the surface of metallic structures due to an external
electromagnetic (EM) excitation [1, 2, 3, 4]. The quasiparticles of plasma
oscillations on the surface of metallic nanoparticles (NPs), called localized
surface plasmons (LSPs), are known for their ability to strongly enhance and
confine the incident light in subwavelength volumes [5]. Owing to these
properties, plasmonic NPs constitute ideal platforms for various applications,
such as in optical communications [6], photovoltaics [7] and sensors [8]. To
increase the practical benefit from such systems, it has long been a trend to
minimize the dimensions of the NPs involved, in order to obtain a more
dramatic confinement of light. This goal is achievable as modern technologies
allow the fabrication and characterisation of metallic structures of ever
decreasing nanometric dimensions. At the extreme nanoscale of sizes $\lesssim
10\ \mathrm{nm}$, the local response approximation (LRA), based on the Drude
model [9, 10], or on the experimentally measured permittivity of the bulk
metal, can no longer accurately predict the optical response of the system.

Quantum-informed plasmonics [11] aims to incorporate nonlocal effects of
quantum-mechanical origin into the description of the system [12], including
screening, Landau damping, and electron spill-out [13]. Screening, in the case
of a metallic NP, refers to the free electrons experiencing a smaller
restoring force, due to charge shielding from the bound electrons in low-lying
bands, leading to the electron density spilling inwards at the assumed metal
boundary [14]. Landau damping describes the decay of the plasma oscillations
into electron–hole pairs, and is enhanced in small NPs due to collisions of
conduction electrons with the metal boundary [15, 16]. Electron spill-out
describes the fact that the electron density does not terminate abruptly, but
rather has a smoother profile, that spills out of the formal boundary of the
metal [17]. One of the most widely used nonlocal models that takes screening
into account is the Hydrodynamic Drude model (HDM) [18, 19, 20]; by extending
HDM to the generalized nonlocal optical response (GNOR) model [21], we can
additionally reproduce Landau damping. However, both HDM and GNOR typically
rely on a hard-wall boundary condition, which assumes that the electron
density is zero beyond the formal boundary of the NP [22, 23], meaning that
electron spill-out is overlooked; self-consistent variants of HDM that can
take this effect into account have been proposed, but they remain
computationally heavy [24, 25, 26].

Whether screening or spill-out is the dominant mechanism determining the
optical response of the metal depends on its work function. In alkali metals,
that are characterized by low work functions, spill-out is most significant
[27, 28, 29]. An elegant way to incorporate this effect in the description of
the system is offered by the Feibelman formalism, in which nonlocality is
encompassed by means of surface-response functions [30]. This framework
invokes new boundary conditions at the air–metal interface using the surface-
response parameters $d_{\perp}$ and $d_{\parallel}$, which correspond to an
induced surface charge and current density, respectively [31]. It thus allows
the bulk response of the plasmonic material to still be locally modelled,
while introducing a discontinuity due to the surface polarization [30]. In
general, the $d$-parameters are complex numbers with dimensions of length; the
real parts of the $d$-parameters express the position of the charge/current
centroid, whereas the imaginary parts account for surface-enabled Landau
damping [32].

Due to the strong field enhancement and confinement that plasmonic NPs
achieve, they are often used as (open) nanocavities [33], in order to enhance
the excitation and emission processes of nearby emitters [34, 35, 36, 37, 38].
When a cavity is coupled strongly to an emitter, modal hybridization and
reversible energy exchange can occur between the two components [39, 40, 41,
42]. One advantage of plasmonic nanocavities is the large control over the
plasmonic mode by the modulation of their size, which is enabled by
nanofabrication techniques [43]. To achieve high coupling strengths, plasmons
are most commonly combined with excitonic materials that sustain collective
optical excitations, such as low-dimensional semiconductors or organic dye
molecules [44, 45, 46, 47]. By tuning the plasmonic mode and reducing the size
of the metallic NPs—and thus the mode volume, plasmon–exciton geometries show
strong potential for generating wide Rabi splitting, and thus enabling the
study of quantum effects at room temperature [48, 49].

The effect of screening and Landau damping in strongly-coupled, composite
core–shell NPs consisting of a noble metal and an excitonic material has been
studied earlier, within the HDM and GNOR models. In particular, it was found
that the width of the Rabi splitting is not affected by such corrections. On
the one hand, both hybrid modes shift in energy by the same amount; on the
other hand, for such small NPs the increase of the effective mode volume due
to Landau damping is typically rather small, and its influence on the
anticrossing is counterbalanced by the nonlocal blueshift of the modes [50].
Nevertheless, it was speculated that the situation should be different in
cases where spill-out is the dominant quantum effect. In this work, the
surface-response formalism (SRF) framework is applied to determine whether
electron spill-out indeed affects the coupling strength in NPs consisting of
sodium and excitonic components. Specifically, we explore whether it
influences the width of the mode anticrossing, when compared to the LRA
prediction.

##  II Results and discussion

In what follows, we study the optical response of core–shell NPs in two
different configurations, focusing on corrections pertaining to electron
spill-out and Landau damping. Sodium, in which spill-out is non-negligible
[51, 52], is the plasmonic material considered throughout the paper, since its
simple atomic structure has rendered it the workhorse of theoretical quantum
plasmonics [53, 54, 55]. Within LRA, sodium is modelled with the simple Drude
model, expressing the relative permittivity $\varepsilon$ as a function of
angular frequency $\omega$ as

$\displaystyle\varepsilon=1-\frac{\omega_{\mathrm{p}}^{2}}{\omega\left(\omega+\mathrm{i}\gamma_{\mathrm{p}}\right)},$
(1)

with plasma frequency $\omega_{\mathrm{p}}$ and damping rate
$\gamma_{\mathrm{p}}$, given as $\hbar\omega_{\mathrm{p}}=5.89132$ eV and
$\hbar\gamma_{\mathrm{p}}=0.1$ eV. From now on, all frequencies will appear
normalized to $\omega_{\mathrm{p}}$.

Figure 1: (a) Real (blue line) and imaginary (red line) part of the
perpendicular Feibelman parameter $d_{\perp}$ for a flat sodium–air interface
[31]. (b), (c) Contour maps of the extinction cross section (in logarithmic
scale), normalized to the geometric cross section $\pi R_{1}^{2}$, for varying
radius of a spherical sodium NP, obtained with (b) LRA and (c) SRF, sharing a
common color scale. The frequency of the incident light $\omega$ is normalized
to the plasma frequency $\omega_{\mathrm{p}}$ of sodium in all panels. The
white dashed line in (b) marks the quasistatic frequency
$\omega_{\mathrm{p}}/\sqrt{3}$; in (c), the white dashed line gives the
corresponding frequency according to Eq. (4).

The optical response of the NPs is investigated using Mie theory [56]; the
incident, scattered, and internal fields are expanded in a sum of spherical
waves of multipolar order $\ell$, and the connection between their amplitudes
is given by the scattering Mie coefficients $T_{E\ell}$ and $T_{H\ell}$ (where
$E,H$ stands for transverse-electric and transverse-magnetic polarization,
respectively). The Mie coefficients are straightforward to calculate through
the usual pillbox boundary conditions at the particle (medium 1)–environment
(medium 2) interface, which, in the absence of surface charges and currents,
translate into continuity of the tangential components of the electric and
magnetic fields, $\mathbf{E}_{2\parallel}-\mathbf{E}_{1\parallel}=0$ and
$\mathbf{H}_{2\parallel}-\mathbf{H}_{1\parallel}=0$ [9].

In SRF, the Feibelman $d$-parameters are introduced in a modified Mie
description. These parameters, typically calculated for an infinite flat
metal–dielectric interface, are defined as [30]

$\displaystyle
d_{\perp}(\omega)=\frac{\int\mathrm{d}x\;x\rho(x,\omega)}{\int\mathrm{d}x\;\rho(x,\omega)}~{},\quad
d_{\parallel}(\omega)=\frac{\int\mathrm{d}x\;x\partial_{x}J_{y}(x,\omega)}{\int\mathrm{d}x\;\partial_{x}J_{y}(x,\omega)},$
(2)

where $x$ is the direction normal to the interface, $\rho$ is the induced
charge density, and $\mathbf{J}$ the corresponding current density. They can
be incorporated into the boundary conditions of the problem [31],

$\displaystyle\mathbf{E}_{2\parallel}-\mathbf{E}_{1\parallel}$
$\displaystyle=-d_{\perp}\nabla_{\parallel}\left(E_{2\perp}-E_{1\perp}\right),$
(3a) $\displaystyle\mathbf{H}_{2\parallel}-\mathbf{H}_{1\parallel}$
$\displaystyle=\mathrm{i}\omega
d_{\parallel}\left(\mathbf{D}_{2\parallel}-\mathbf{D}_{1\parallel}\right)\times\mathbf{\hat{n}},$
(3b)

where $\mathbf{D}$ is the displacement field, $\nabla_{\parallel}$ is the
surface nabla operator, $\perp$ denotes the field component normal to the
interface, and $\mathbf{\hat{n}}$ the unit vector normal to the interface. It
should be noted that the use of Feibelman parameters for curved interfaces
imposes requirements on the dimensions of the system. For a spherical particle
of radius $R$, we require $R\gg\ell d_{\perp,\parallel}$ [31]; the
applicability of the formalism can further be extended by considering nonlocal
Feibelman parameter [57, 58]. The modified Mie coefficients used throughout
this paper can be found in the Appendix.

Figure 2: (a, b) Extinction cross sections, normalized to the geometric cross
section $\pi R^{2}$, for a sodium core of radius (a) $2$ nm and (b) $7.3$ nm,
both with an excitonic shell of $2$ nm thickness obtained from SRF (red solid
line), compared to LRA (black solid line), showing the two hybrid modes
$\omega_{\pm}$. The dashed lines correspond to the sodium core alone. (c, d)
Hybrid mode resonance energies within SRF (red dots) and within LRA (black
triangles) as a function of (c) the sodium core radius and (d) the sodium core
resonance frequency $\omega_{\mathrm{r}}$, normalized to
$\omega_{\mathrm{p}}$. The red lines show the trend of the resonance energies
in SRF for both hybrid modes compared to LRA (black lines). does not signify
an increase in coupling strength but the shift of the uncoupled resonance
frequency $\omega_{\mathrm{r}}$. (e, d) Logarithmically scaled extinction
contour maps as a function of resonance frequency $\omega_{\mathrm{r}}$ for
(e) LRA and (f) SRF, sharing a common color scale.

###  II.1 Surface Response Formalism for Sodium Nanoparticles

To demonstrate the effect of spill-out on the optical response of plasmonic
NPs, we first consider a sodium sphere of radius $R_{1}$ embedded in air.
Following Ref. [31], the $d_{\perp}$ parameter is given by Lorentzian fitting
to time-dependent density-functional theory calculations, where sodium was
treated as a jellium of Wigner-Seitz radius $r_{\mathrm{s}}=4$ [59]. A plot of
the real and imaginary part of $d_{\perp}$ is presented in Fig. 1a. We set
$d_{\parallel}=0$, which is appropriate for charge-neutral surfaces [30, 32].

Figs. 1b and c show the extinction cross sections of a sodium NP surrounded by
air, for varying radii $R_{1}$, calculated within LRA and SRF, respectively.
Up to a radius of around $R_{1}=25$ nm, the only mode contributing to the
extinction spectra for both LRA and SRF is the electric dipolar (ED) LSP. As
the radius is decreased, the ED resonance energy in LRA converges to the
quasistatic limit
$\omega_{1}^{\mathrm{(cl)}}/\omega_{\mathrm{p}}=1/\sqrt{3}\simeq 0.58$. This
value is indicated with the dashed white line in Fig. 1b, and captures the
absence of retardation of the free electrons in the NP, which experience the
incident field as spatially constant. In SRF, in which the electron density
varies over the NP volume, and thus the field that the NP experiences cannot
be assumed constant, the ED resonance energy deviates from the local
quasistatic limit and follows the relation

$\displaystyle\omega_{1}^{\mathrm{(q)}}=\omega_{1}^{\mathrm{(cl)}}\sqrt{1-2d_{\perp}/R_{1}},$
(4)

as indicated by the dashed white line in Fig. 1c. As the NP dimensions
increase, and retardation becomes important, higher multipole contributions
emerge. The electric quadrupolar (EQ) mode is seen in Figs. 1b and c at
approximately $0.61\omega_{\mathrm{p}}$. The effect of spill-out and Landau
damping is most prominent for NP radii below $5$ nm, manifested in SRF as a
red-shift and a broadening of the plasmon resonance in the extinction
spectrum, compared to the local model.

###  II.2 Sodium Core with an Excitonic Shell

Having now established the corrections to the plasmonic core alone, we
introduce an excitonic material to obtain strong coupling. Here, we
investigate a spherical sodium core of varying radius $R_{1}$, encapsulated in
an excitonic shell with a fixed thickness $W=2$ nm, embedded in air. The
excitonic shell is modelled with the Drude-Lorentz model,

$\displaystyle\varepsilon_{\mathrm{exc}}(\omega)=1-\frac{f\omega_{\mathrm{exc}}}{\omega^{2}-\omega_{\mathrm{exc}}^{2}+\mathrm{i}\omega\gamma_{\mathrm{exc}}},$
(5)

with excitonic energy $\hbar\omega_{\mathrm{exc}}=3.38$ eV, damping rate
$\hbar\gamma_{\mathrm{exc}}=0.06$ eV, and oscillator strength $f=0.02$ [60].
The scattering Mie coefficients core–shell NP are given by (11) and (12) in
the Appendix.

The spectra of the extinction cross section for NPs with $R_{1}=2$ nm and
$R_{1}=7.3$ nm sodium core are shown in Fig. 2a and Fig. 2b, respectively. The
SRF (LRA) extinction spectrum is given by the red (black) line, with the
dashed lines denoting the extinction spectra for the sodium core in the
absence of the excitonic shell. All spectra are normalized to the geometric
cross section $\pi R^{2}$, where $R=R_{1}+W$ for the core–shell NP, and
$R=R_{1}$ for the sodium core alone. The interaction of the excitonic material
with the LSP forms two hybrid modes, as expected from literature [61], with a
mixed plasmonic-excitonic chracter. Their energies are denoted as
$\omega_{\pm}$ for the lower ($-$) and upper ($+$) hybrid mode, respectively.
Between the two hybrid modes there is a weak peak at
$\omega/\omega_{\mathrm{p}}=0.57$. This is a geometric resonance of the charge
oscillations that matches the excitonic energy of the system, where the
plasmonic core and excitonic shell act as a uniform material [61], and it does
not participate in the coupling.

For both core radii, the introduction of nonlocal effects through SRF produces
a redshift of both hybrid modes compared to the local model, albeit a more
significant shift for the upper hybrid mode than for the lower one.
Additionally, both modes in the nonlocal model experience broadening compared
to the local model, again more significant for the upper hybrid mode. The
redshift and damping of the hybrid modes is larger for the smaller core radius
of $R_{1}=2$ nm in Fig. 2a, compared to $R_{1}=7.3$ nm in Fig. 2b. This can be
explained by the electron spill-out effectively increasing the core radius,
allowing for LSPs of longer wavelengths. This, in turn, introduces a detuning
of the plasmon energy for SRF compared to LRA, red-shifting the hybrid modes.

To better visualize the trends of the hybrid modes, a comparison of SRF and
LRA is provided in Figs. 2c and d, showing the resonance energies of the
hybrid modes as functions of select sodium-core radii and of select sodium-
core resonance frequencies $\omega_{\mathrm{r}}$, respectively. The resonance
frequency $\omega_{\mathrm{r}}$ is given by the LSP of the sodium core in the
absence of the excitonic shell for varying $R_{1}$. The radius of the sodium
core in Fig. 2c varies from $1$ to $35$ nm, while the width of the excitonic
shell is fixed at $2$ nm. We observe that for both hybrid modes, the
differences in energy between LRA and SRF hybrid modes are most significant
for a decreasing core radius $R_{1}$ smaller than about $R_{1}\approx 3$ nm.
At a $1$ nm core radius, which lies at the limits of validity of SRF, the
energy difference between LRA and SRF is $0.0035\omega_{\mathrm{p}}$ for the
upper mode, and $0.0037\omega_{\mathrm{p}}$ for the lower mode. For LRA, a
decreasing radius corresponds to an increase in $\omega_{\mathrm{r}}$. But for
SRF, a decrease in core radius only leads to an increase in
$\omega_{\mathrm{r}}$ up to about $0.574\omega_{\mathrm{p}}$. At around this
frequency, referencing Fig. 1b, further decrease in $R_{1}$ now leads to a
decrease in $\omega_{\mathrm{r}}$. As both peaks of the SRF hybrid modes are
redshifted in energy compared to LRA, the energy difference, and thus the
coupling strength between the hybrid modes remains practically unaffected in
SRF compared to LRA at each core radius.

Figure 3: (a) Extinction cross section of a spherical excitonic core with
radius $R_{1}=5$ nm embedded in air. (b, c) Extinction contour maps as a
function of shell width for (b) a vacuum core with a radius of $5$ nm with a
sodium shell, and (c) the combined excitonic core–sodium shell. The color bar
is shared for (b) and (c). All results are obtained within LRA.

When considering mode splitting, we commonly address the spectrum where zero
detuning between the two modes occurs, namely, where the resonance
$\omega_{\mathrm{r}}$ of the sodium core matches the resonance of the
excitonic shell. Zero detuning in the local model occurs for a core radius of
$R_{1}=7.3$ nm, with the energy difference between the hybrid mode peaks found
to be $0.048\omega_{\mathrm{p}}$. In the nonlocal model, zero detuning occurs
at a lower core radius than for LRA, at $R_{1}=6.0$ nm, with the energy
difference between hybrid modes being $0.051\omega_{\mathrm{p}}$ ($\simeq 300$
meV). Comparing SRF to LRA, this is a difference of
$0.0027\omega_{\mathrm{p}}$, with the larger anticrossing occurring in the
nonlocal model. With practically no change in the coupling strength from LRA
to SRF at each core radius, this change in anticrossing is attributed to the
fact that zero detuning occurs for two different radii.

Lastly, to examine the effect of resonance broadening due to Landau damping,
we consult the contour plots of Figs. 2e and f, where we plot the extinction
cross section of the coupled system as a function of the resonance energy
$\omega_{\mathrm{r}}$, in LRA and SRF, respectively. Clearly, no significant
broadening can be observed, and the two plots are nearly identical, since SRF
reproduces a similar coupling strength. For high frequencies of the incident
light and low resonance frequencies $\omega_{\mathrm{r}}$, we see also the EQ
($\ell=2$) resonance for the sodium core, which, due to its large detuning to
the excitonic energy, does not participate in the coupling.

###  II.3 Excitonic Core with a Sodium Shell

We now investigate the inverted geometry, i.e. a spherical excitonic core with
a fixed radius $R_{1}=5$ nm, encapsulated in a sodium shell of varying width
$W$, and subsequently varying outer radius $R=R_{1}+W$. The excitonic core is
now modelled with an excitonic energy of $\hbar\omega_{\mathrm{exc}}=5.2$ eV
($\omega_{\mathrm{exc}}=0.88\omega_{\mathrm{p}}$), while all other parameters
remain the same as in section II.2. Because the new geometry has two plasmonic
surfaces, Feibelman corrections must apply on both. The nonlocal scattering
Mie coefficients are now given by (13) and (14), while for the inner boundary
we use (9) and (10) in the Appendix. When considering the inner boundary of
the sodium shell, it should be noted that the spill-out occurs radially
towards the centre of the coordinate system. The Feibelman parameters,
expressing lengths evaluated with respect to the metal–environment interface,
are directionally dependent, with the positive direction previously assumed to
be radially outwards. Therefore, the sign of the Feibelman parameter shown in
Fig. 1 must be inverted when dealing with the inner boundary of the sodium
shell.

Fig. 3a shows the extinction cross section of the excitonic core with a radius
of $5\ \mathrm{nm}$, while Fig. 3b shows a contour map of the LRA extinction
cross section of a sodium shell thickness $W$ ranging between $1$ and $15$ nm,
in the absence of the excitonic core. The spectrum of the excitonic core is a
singular peak at $\omega_{\mathrm{exc}}/\omega_{\mathrm{p}}\simeq 0.88$. The
sodium shell has two distinct LSP modes; the higher energy mode is a cavity
mode, originating from LSPs at the inner surface of the shell, while the lower
mode is the NP mode, originating from the outer surface of the shell [62]. A
contour map of the LRA extinction cross section for the sodium shell in the
presence of the excitonic core is presented in Fig. 3c, showing an
anticrossing emerging between the cavity mode and the excitonic transition, at
the excitonic energy, indicating strong coupling. As the cavity mode is the
one participating in the coupling, in what follows we only consider this mode
when discussing the effect of the Feibelman corrections.

Figure 4: (a, b) Extinction cross sections, normalized to the geometric cross
section $\pi R^{2}$, for a sodium shell of width (a) $3.14$ nm and (b) $13.0$
nm around an excitonic core with a $5$ nm radius. The hybrid modes
$\omega_{\pm}$ for SRF (LRA) are shown as the solid red (black) lines. The
dashed red (black) lines are the extinction cross section of the sodium shell
within SRF (LRA) with an air core instead. (c, d) Logarithmically scaled
extinction contour maps as a function of resonance frequency
$\omega_{\mathrm{r}}$ for (c) LRA and (d) SRF, sharing a common color scale.
(e) Hybrid mode resonance energies, shown as red dots (black triangles) within
SRF (LRA), as a function of the sodium shell resonance energy. The red lines
show the trend of the resonance frequencies of SRF for both hybrid modes.

Figs. 4a and b show the extinction spectrum for the sodium shell with an air
or excitonic core, for shell widths $W=3.14$ nm and $W=13.0$ nm, respectively.
When considering nonlocal effects, the uncoupled cavity mode of the plasmonic
shell experiences damping and blueshift compared to the LRA spectrum, though
it can be seen that the blueshift is greater for the smaller shell width; the
upper hybrid modes on both graphs follow the same trend. Looking at the lower
hybrid modes in the two spectra, they both experience damping in SRF compared
to LRA, but in Fig. 4b we see a redshift for the lower hybrid mode. The
opposite shifts are the result of the nonlocality affecting both the coupling
of the cavity mode to the NP mode and the cavity mode alone. For this
geometry, the electron spill-out effectively increases the shell width
compared to LRA, leading to a spectral redshift of the cavity mode compared to
LRA. At the same time, at the inner shell boundary, the electron spill-out
effectively decreases the inner radius $R_{1}$ of the shell, meaning that the
LSPs in the cavity for the SRF model will have shorter wavelengths than those
in the LRA model, leading to a blueshift of the SRF cavity mode spectrum
compared to LRA. In the quasistatic limit, we find the correction to the
electric dipole resonance frequency of the cavity mode to be

$\displaystyle\omega_{1\mathrm{c}}^{(\mathrm{q})}=\omega_{1\mathrm{c}}^{(\mathrm{cl})}\sqrt{1-d_{\perp}/R_{1}},$
(6)

where the quasistatic LRA cavity-mode resonance frequency is
$\omega_{1\mathrm{c}}^{\mathrm{(cl)}}/\omega_{\mathrm{p}}=\sqrt{2/3}\simeq
0.82$. These two nonlocal contributions result in an overall blueshift of the
lower mode at most shell widths, such as the one shown in Fig. 4a, indicating
that the contribution of the correction to the coupling is more significant.
However, at the larger shell widths, such as the one shown in Fig. 4b with an
overall redshift of the lower mode, the correction to the NP mode becomes the
more significant contribution.

To compare the mode splitting quantatively, we consider both models at zero
detuning, i.e. when the cavity mode of the sodium shell in the absence of the
excitonic core has resonance frequency $\omega_{\mathrm{r}}$ that equals the
excitonic frequency $\omega_{\mathrm{exc}}$. The tuning of the energy is done
by varying the width of the sodium shell while keeping the inner radius
constant. Zero detuning occurs for the LRA shell, seen in Fig. 4a at width
$W=3.14$ nm, where the difference in energy between the resonance peaks is
$312$ meV. For SRF, zero detuning occurs at a slightly larger width $W=3.26$
nm with a difference in energy of $314$ meV. Comparing the two values, it is
apparent that the energy difference between the corresponding coupling
strengths is negligible.

Focusing on the coupled cavity mode across the entire resonance frequency
spectrum, Figs. 4c and d show the contour maps of the extinction cross section
for LRA and SRF, respectively, as a function of the plasmonic resonance
frequency $\omega_{\mathrm{r}}/\omega_{\mathrm{p}}$ of the uncoupled cavity
mode. The plasmonic resonance energy is calculated as the resonance of the
cavity mode of the sodium shell in the absence of the excitonic core for a
varying shell width $W$ and a fixed inner radius $R_{1}=5$ nm, and the shell
width $W$ varies from $1$ to $15$ nm. On both panels, it can be seen that the
splitting is centred around the excitonic energy
$\omega_{\mathrm{exc}}/\omega_{\mathrm{p}}\simeq 0.88$, as expected. Comparing
Figs. 4c and d, the most notable difference is the broadening of the hybrid
modes in SRF compared to LRA due to Landau damping, especially for the lower
mode. Otherwise, the hybrid modes are still discernible in the strong coupling
regime.

To better visualize the splitting in the SRF and LRA models, Fig. 4e shows the
resonance energies of the hybrid modes of select plasmonic resonance energies
$\omega_{\mathrm{r}}$. We see that, compared to the energy difference between
the hybrid modes, the prior described shifts of the SRF modes from the LRA
modes become negligible. This means that the energy split between the hybrid
modes is practically not affected by nonlocality for all the shown plasmonic
resonance energies.

##  III Conclusion

We examined the optical response of different core–shell configurations of
plasmonic–excitonic composites and compared the results obtained with standard
LRA to those produced with a nonlocal surface-response approach. Specifically,
we chose to work within the Feibelman framework, which is ideal for describing
the spill-out in alkali metals, together with surface-enabled Landau damping.
Even though energy shifts and resonance broadenings of the hybrid modes are
evident in the spectra, both models generate very similar anticrossings in all
configurations, since both hybrid modes shift within SRF by the same amount.
We thus found that LRA provides a reliable estimation of the coupling
strength, even when the dimensions of metallic–excitonic systems are driven
down to the limits of nanofabrication capabilities.

## Acknowledgments

We acknowledge funding from VILLUM Fonden (Villum Investigator, N. A.
Mortensen, grant No. 16498). The Center for Polariton-driven Light–Matter
Interactions (POLIMA) is sponsored by the Danish National Research Foundation
(Project No. DNRF165).

## Appendix

Mie theory is used throughout this paper to calculate extinction cross-
sections and predict the optical response of the spherical NPs presented. The
coefficients and equations for the cross sections that we have used are
presented in this appendix. For all materials modelled, the relative
permeability is always assumed $\mu=1$.

For a spherical NP with permittivity $\varepsilon_{1}$ and permeability
$\mu_{1}$ embedded in a host medium with permittivity $\varepsilon_{2}$ and
permeability $\mu_{2}$, the scattering Mie coefficients using an LRA framework
are [63, 64]

$\displaystyle T_{E\ell}$
$\displaystyle=\frac{j_{\ell}(x_{11})\Psi^{\prime}_{\ell}(x_{21})\varepsilon_{1}-j_{\ell}(x_{21})\Psi^{\prime}_{\ell}(x_{11})\varepsilon_{2}}{h^{(1)}_{\ell}(x_{21})\Psi^{\prime}_{\ell}(x_{11})\varepsilon_{2}-j_{\ell}(x_{11})\xi^{\prime}_{\ell}(x_{21})\varepsilon_{1}}$
(7) $\displaystyle T_{H\ell}$
$\displaystyle=\frac{j_{\ell}(x_{11})\Psi^{\prime}(x_{21})\mu_{1}-j_{\ell}(x_{21})\Psi_{\ell}^{\prime}(x_{11})\mu_{2}}{h^{(1)}_{\ell}(x_{21})\Psi^{\prime}(x_{11})\mu_{2}-j_{\ell}(x_{11})\xi_{\ell}^{\prime}(x_{21})\mu_{1}},$
(8)

where $x_{i1}=k_{i}R_{1}$ for wavenumbers $k_{1}$ and $k_{2}$ of the core and
host material, respectively. The functions $\Psi_{\ell}^{\prime}(z)$ and
$\xi_{\ell}^{\prime}(z)$ are the derivatives with respect to their argument
$z$ of the Ricatti-Bessel functions $\Psi_{\ell}(z)=zj_{\ell}(z)$ and
$\xi_{\ell}(z)=zh^{(1)}_{\ell}(z)$, where $j_{\ell}(z)$ and
$h^{(1)}_{\ell}(z)=j_{\ell}(z)+\mathrm{i}y_{\ell}(z)$ are the spherical Bessel
functions and spherical Hankel functions of the first kind, respectively (with
$y_{\ell}(z)$ being the spherical Bessel function of the second kind, a.k.a
Neumann function).

We now introduce the Feibelman correction to a spherical NP of plasmonic
material. The scattering Mie coefficients in this SRF framework become [31]

$\displaystyle
T_{E\ell}=\frac{C^{\mathrm{n}}_{E\ell}+\bar{\varepsilon}[\bar{d}_{\perp}j_{\ell}(x_{21})j_{\ell}(x_{11})+\bar{d}_{\parallel}\Psi_{\ell}^{\prime}(x_{21})\Psi_{\ell}^{\prime}(x_{11})]}{C^{\mathrm{d}}_{E\ell}-\bar{\varepsilon}[\bar{d}_{\perp}h^{(1)}_{\ell}(x_{21})j_{\ell}(x_{11})+\bar{d}_{\parallel}\xi_{\ell}^{\prime}(x_{21})\Psi_{\ell}^{\prime}(x_{11})]}$
(9) $\displaystyle
T_{H\ell}=\frac{C^{\mathrm{n}}_{H\ell}+[\mu_{2}x_{11}^{2}-\mu_{1}x_{21}^{2}]\bar{d}_{\parallel}j_{\ell}(x_{21})j_{\ell}(x_{11})}{C^{\mathrm{d}}_{H\ell}-[\mu_{2}x_{11}^{2}-\mu_{1}x_{21}^{2}]\bar{d}_{\parallel}h^{(1)}_{\ell}(x_{21})j_{\ell}(x_{11})},$
(10)

where $\bar{\varepsilon}=\varepsilon_{1}-\varepsilon_{2}$,
$\bar{d}_{\perp}={\ell(\ell+1)}d_{\perp}/{R_{1}}$ and
$\bar{d}_{\parallel}=d_{\parallel}/{R_{1}}$. The coefficients
$C^{\mathrm{n}}_{P\ell}$ and $C^{\mathrm{d}}_{P\ell}$ are the numerator and
denominator from Eqs. (7) and (8),

$\displaystyle C^{\mathrm{n}}_{P\ell}=$ $\displaystyle
j_{\ell}(x_{11})\Psi^{\prime}_{\ell}(x_{21})p_{1}-j_{\ell}(x_{21})\Psi^{\prime}_{\ell}(x_{11})p_{2}$
$\displaystyle C^{\mathrm{d}}_{P\ell}=$ $\displaystyle
h^{(1)}_{\ell}(x_{21})\Psi^{\prime}_{\ell}(x_{11})p_{2}-j_{\ell}(x_{11})\xi^{\prime}_{\ell}(x_{21})p_{1},$

where $p_{m}=\varepsilon_{m}$ if $P=E$ and $p_{m}=\mu_{m}$ if $P=H$. Note that
if $d_{\perp,\parallel}=0$, the equations will become identical to Eqs. (7)
and (8).

We further add an excitonic shell to our plasmonic core NP; for a
configuration with a non-plasmonic shell the scattering Mie-coefficients
become

$\displaystyle T_{E\ell}$
$\displaystyle=\frac{[\Psi_{\ell}^{\prime}(x_{32})h^{(1)}_{\ell}(x_{22})\varepsilon_{2}-j_{\ell}(x_{32})\xi_{\ell}^{\prime}(x_{22})\varepsilon_{3}]T_{E\ell}^{\prime}+S^{\mathrm{n}}_{E\ell}}{[h^{(1)}_{\ell}(x_{32})\xi_{\ell}^{\prime}(x_{22})\varepsilon_{3}-\xi_{\ell}^{\prime}(x_{32})h_{\ell}(x_{22})\varepsilon_{2}]T_{E\ell}^{\prime}+S^{\mathrm{d}}_{E\ell}}$
(11) $\displaystyle T_{H\ell}$
$\displaystyle=\frac{[\Psi_{\ell}^{\prime}(x_{32})h^{(1)}_{\ell}(x_{22})\mu_{2}-j_{\ell}(x_{32})\xi_{\ell}^{\prime}(x_{22})\mu_{3}]T_{H\ell}^{\prime}+S^{\mathrm{n}}_{H\ell}}{[h^{(1)}_{\ell}(x_{32})\xi_{\ell}^{\prime}(x_{22})\mu_{3}-\xi_{\ell}^{\prime}(x_{32})h_{\ell}(x_{22})\mu_{2}]T_{H\ell}^{\prime}+S^{\mathrm{d}}_{H\ell}},$
(12)

where the arguments $x_{i2}=k_{i}R$, are in terms of the outer radius
$R=R_{1}+W$ and the wavenumbers $k_{2}$ and $k_{3}$ in the excitonic shell and
host medium, air, respectively. $T^{\prime}_{P\ell}$ for polarization $P=E,H$
are the Mie coefficients for the core, using (7) and (8) for LRA and (9) and
(10) for an SRF core. The coefficients $S^{\mathrm{n}}_{P\ell}$ and
$S^{\mathrm{d}}_{P\ell}$ are given by

$\displaystyle S^{\mathrm{n}}_{P\ell}$
$\displaystyle=\Psi_{\ell}^{\prime}(x_{32})j_{\ell}(x_{22})p_{2}-j_{\ell}(x_{32})\Psi_{\ell}^{\prime}(x_{22})p_{3}$
$\displaystyle S^{\mathrm{d}}_{P\ell}$
$\displaystyle=h^{(1)}_{\ell}(x_{32})\Psi_{\ell}^{\prime}(x_{22})p_{3}-\xi_{\ell}^{\prime}(x_{32})j_{\ell}(x_{22})p_{2},$

where $p_{m}=\varepsilon_{m}$ if $P=E$ and $p_{m}=\mu_{m}$ if $P=H$. These
coefficients are the LRA Mie coefficients evaluated at the shell-host
boundary.

Inverting the geometry such that we have a non-plasmonic core and a plasmonic
shell, the scattering Mie coefficients become

$\displaystyle T_{E\ell}$
$\displaystyle=\frac{K^{\mathrm{n}}_{E\ell}+S^{\mathrm{n}}_{E\ell}+(\varepsilon_{2}-\varepsilon_{3})(D^{\mathrm{n}}_{E\ell}+N^{\mathrm{n}}_{E\ell})}{K^{\mathrm{d}}_{E\ell}+S^{\mathrm{d}}_{E\ell}-(\varepsilon_{2}-\varepsilon_{3})(D^{\mathrm{d}}_{E\ell}+N^{\mathrm{d}}_{E\ell})}$
(13) $\displaystyle T_{H\ell}$
$\displaystyle=\frac{K^{\mathrm{n}}_{H\ell}+S^{\mathrm{n}}_{H\ell}+(x_{22}^{2}\mu_{3}-x_{32}^{2}\mu_{2})(D^{\mathrm{n}}_{H\ell}+N^{\mathrm{n}}_{H\ell})}{K^{\mathrm{d}}_{H\ell}+S^{\mathrm{d}}_{H\ell}-(x_{22}^{2}\mu_{3}-x_{32}^{2}\mu_{2})(D^{\mathrm{d}}_{H\ell}+N^{\mathrm{d}}_{H\ell})},$
(14)

where the coefficients $N^{n,d}_{P\ell}$ are

$\displaystyle N^{\mathrm{n}}_{E_{\ell}}$
$\displaystyle=[\bar{d}_{\perp}h_{\ell}^{(1)}(x_{22})h_{\ell}^{(1)}(x_{32})+\bar{d}_{\parallel}\xi^{\prime}_{\ell}(x_{22})\xi^{\prime}_{\ell}(x_{32})]T_{E\ell}^{\prime}$
$\displaystyle N^{\mathrm{d}}_{E\ell}$
$\displaystyle=[\bar{d}_{\perp}h_{\ell}^{(1)}(x_{22})h_{\ell}^{(1)}(x_{32})+\bar{d}_{\parallel}\xi^{\prime}_{\ell}(x_{22})\xi^{\prime}_{\ell}(x_{32})]T_{E\ell}^{\prime}$
$\displaystyle N^{\mathrm{n}}_{H\ell}$
$\displaystyle=\bar{d}_{\parallel}h_{\ell}^{(1)}(x_{22})j_{\ell}(x_{32})T_{H\ell}^{\prime}$
$\displaystyle N^{\mathrm{d}}_{H\ell}$
$\displaystyle=\bar{d}_{\parallel}h_{\ell}^{(1)}(x_{22})h^{(1)}_{\ell}(x_{32})T_{H\ell}^{\prime}$

We further introduce the coefficients $D^{n,d}_{P\ell}$ as

$\displaystyle D^{\mathrm{n}}_{E\ell}$
$\displaystyle=[\bar{d}_{\perp}j_{\ell}(x_{22})j_{\ell}(x_{32})+\bar{d}_{\parallel}\Psi^{\prime}_{\ell}(x_{22})\Psi^{\prime}_{\ell}(x_{32})]$
$\displaystyle D^{\mathrm{d}}_{E\ell}$
$\displaystyle=[\bar{d}_{\perp}j_{\ell}(x_{22})h_{\ell}^{(1)}(x_{32})+\bar{d}_{\parallel}\Psi^{\prime}_{\ell}(x_{22})\xi^{\prime}_{\ell}(x_{32})]$
$\displaystyle D^{\mathrm{n}}_{H\ell}$
$\displaystyle=\bar{d}_{\parallel}j_{\ell}(x_{22})j_{\ell}(x_{32})$
$\displaystyle D^{\mathrm{d}}_{H\ell}$
$\displaystyle=\bar{d}_{\parallel}j_{\ell}(x_{22})h^{(1)}_{\ell}(x_{32}),$

which are terms also shown in (9) and (10), they represent a similar
correction to the shell–host boundary as previously considered for the
core–host boundary. Lastly, we introduce the $K^{n,d}_{P\ell}$ coefficients

$\displaystyle K^{\mathrm{n}}_{P\ell}$
$\displaystyle=[\Psi_{\ell}^{\prime}(x_{32})h_{\ell}^{(1)}(x_{22})p_{2}-j_{\ell}(x_{32})\xi_{\ell}^{\prime}(x_{22})p_{3}]T_{P\ell}^{\prime}$
$\displaystyle K^{\mathrm{d}}_{P\ell}$
$\displaystyle=[h_{\ell}(x_{32})\xi_{\ell}^{\prime}(x_{22})p_{3}-\xi_{\ell}^{\prime}(x_{32})h_{\ell}^{(1)}(x_{22})p_{2}]T_{P\ell}^{\prime},$

which are recognisable from (11) and (12) for the core–shell NP with non-
plasmonic shell. Note that if $d_{\parallel,\perp}=0$ the equations will
become identical to (11) and (12).

The scattering $\sigma_{\mathrm{sc}}$, absorption $\sigma_{\mathrm{abs}}$ and
extinction $\sigma_{\mathrm{ext}}$ cross-sections for spherical NPs are
calculated using Mie coefficients according to the following equations

$\displaystyle\sigma_{\mathrm{sc}}$
$\displaystyle=\frac{2}{(kR)^{2}}\sum_{\ell=1}^{\infty}(2\ell+1)(|T_{E\ell}|^{2}+|T_{H\ell}|^{2})$
(15) $\displaystyle\sigma_{\mathrm{ext}}$
$\displaystyle=-\frac{2}{(kR)^{2}}\sum_{\ell=1}^{\infty}(2\ell+1)\mathrm{Re}(T_{E\ell}+T_{H\ell})$
(16) $\displaystyle\sigma_{\mathrm{abs}}$
$\displaystyle=\sigma_{\mathrm{ext}}-\sigma_{\mathrm{sc}},$ (17)

where $R$ is the outer radius, $k$ the wave vector of the host medium and
$T_{E\ell}$ and $T_{H\ell}$ are the scattering Mie coefficients for the
electric and magnetic multipoles, respectively.

## References

  * Pines and Bohm [1952] D. Pines and D. Bohm, A collective description of electron interactions: II. Collective $\mathrm{vs}$ individual particle aspects of the interactions, Phys. Rev. 85, 338 (1952).
  * Ritchie [1957] R. H. Ritchie, Plasma losses by fast electrons in thin films, Phys. Rev. 106, 874 (1957).
  * Economou [1969] E. N. Economou, Surface plasmons in thin films, Phys. Rev. 182, 539 (1969).
  * Teng and Stern [1967] Y.-Y. Teng and E. A. Stern, Plasma radiation from metal grating surfaces, Phys. Rev. Lett. 19, 511 (1967).
  * Oulton  _et al._ [2008] R. F. Oulton, G. Bartal, D. F. P. Pile, and X. Zhang, Confinement and propagation characteristics of subwavelength plasmonic modes, New Journal of Physics 10, 105018 (2008).
  * de Waele  _et al._ [2007] R. de Waele, A. F. Koenderink, and A. Polman, Tunable nanoscale localization of energy on plasmon particle arrays, Nano Letters 7, 2004 (2007).
  * Pillai  _et al._ [2007] S. Pillai, K. R. Catchpole, T. Trupke, and M. A. Green, Surface plasmon enhanced silicon solar cells, Journal of Applied Physics 101, 093105 (2007).
  * Willets and Van Duyne [2007] K. A. Willets and R. P. Van Duyne, Localized surface plasmon resonance spectroscopy and sensing, Annu. Rev. Phys. Chem. 58, 267 (2007).
  * Jackson [1998] J. D. Jackson,  _Classical Electrodynamics_ , 3rd ed. (Wiley, 1998).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Harcourt, Orlando, 1976).
  * Tserkezis [2021] C. Tserkezis, Towards quantum nanophotonics: from quantum-informed plasmonics to strong coupling (2021), dr. Techn. thesis (University of Southern Denmark).
  * Tame  _et al._ [2013] M. S. Tame, K. R. McEnery, Ş. K. Özdemir, J. Lee, S. A. Maier, and M. S. Kim, Quantum plasmonics, Nature Physics 9, 329 (2013).
  * Mortensen [2021] N. A. Mortensen, Mesoscopic electrodynamics at metal surfaces, Nanophotonics 10, 2563 (2021).
  * Hohenester [2020] U. Hohenester,  _Nano and Quantum Optics: An Introduction to Basic Principles and Theory_ (Springer, 2020).
  * Li  _et al._ [2013] X. Li, D. Xiao, and Z. Zhang, Landau damping of quantum plasmons in metal nanostructures, New Journal of Physics 15, 023011 (2013).
  * Stamatopoulou and Tserkezis [2022] P. E. Stamatopoulou and C. Tserkezis, Finite-size and quantum effects in plasmonics: manifestations and theoretical modelling, Opt. Mater. Express 12, 1869 (2022).
  * Kreibig and Vollmer [1995] U. Kreibig and M. Vollmer,  _Optical Properties of Metal Clusters_ (Springer, 1995).
  * Ruppin [1973] R. Ruppin, Optical properties of a plasma sphere, Phys. Rev. Lett. 31, 1434 (1973).
  * Barton [1979] G. Barton, Some surface effects in the hydrodynamic model of metals, Rep. Prog. Phys. 42, 963 (1979).
  * García de Abajo [2008] F. J. García de Abajo, Nonlocal effects in the plasmons of strongly interacting nanoparticles, dimers, and waveguides, J. Phys. Chem. C 112, 17983 (2008).
  * Mortensen  _et al._ [2014] N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, and S. I. Bozhevolnyi, A generalized non-local optical response theory for plasmonic nanostructures, Nature Communications 5, 3809 (2014).
  * Melnyk and Harrison [1970] A. R. Melnyk and M. J. Harrison, Theory of optical excitation of plasmons in metals, Phys. Rev. B 2, 835 (1970).
  * Raza  _et al._ [2011] S. Raza, G. Toscano, A.-P. Jauho, M. Wubs, and N. A. Mortensen, Unusual resonances in nanoplasmonic structures due to nonlocal response, Phys. Rev. B 84, 121412 (2011).
  * Toscano  _et al._ [2015] G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl, F. Evers, H. Xu, N. A. Mortensen, and M. Wubs, Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics, Nature Communications 6, 7132 (2015).
  * Yan [2015] W. Yan, Hydrodynamic theory for quantum plasmonics: linear-response dynamics of the inhomogeneous electron gas, Phys. Rev. B 91, 115416 (2015).
  * Ciracì and Della Sala [2016] C. Ciracì and F. Della Sala, Quantum hydrodynamic theory for plasmonics: Impact of the electron density tail, Phys. Rev. B 93, 205405 (2016).
  * Weick  _et al._ [2006] G. Weick, G.-L. Ingold, R. A. Jalabert, and D. Weinmann, Surface plasmon in metallic nanoparticles: Renormalization effects due to electron-hole excitations, Phys. Rev. B 74, 165421 (2006).
  * Stella  _et al._ [2013] L. Stella, P. Zhang, F. J. García-Vidal, A. Rubio, and P. García-González, Performance of nonlocal optics when applied to plasmonic nanostructures, J. Phys. Chem. C 117, 8941 (2013).
  * Teperik  _et al._ [2013a] T. V. Teperik, P. Nordlander, J. Aizpurua, and A. G. Borisov, Quantum effects and nonlocality in strongly coupled plasmonic nanowire dimers, Opt. Express 21, 27306 (2013a).
  * Feibelman [1982] P. J. Feibelman, Surface electromagnetic fields, Progress in Surface Science 12, 287 (1982).
  * Gonçalves  _et al._ [2020] P. A. D. Gonçalves, T. Christensen, N. Rivera, A.-P. Jauho, N. A. Mortensen, and M. Soljačić, Plasmon–emitter interactions at the nanoscale, Nature Communications 11, 366 (2020).
  * Christensen  _et al._ [2017] T. Christensen, W. Yan, A.-P. Jauho, M. Soljačić, and N. A. Mortensen, Quantum corrections in nanoplasmonics: Shape, scale, and material, Phys. Rev. Lett. 118, 157402 (2017).
  * Hugall  _et al._ [2018] J. T. Hugall, A. Singh, and N. F. van Hulst, Plasmonic cavity coupling, ACS Photonics 5, 43 (2018).
  * Bitton  _et al._ [2019] O. Bitton, S. N. Gupta, and G. Haran, Quantum dot plasmonics: from weak to strong coupling, Nanophotonics 8, 559 (2019).
  * Balci  _et al._ [2014] S. Balci, E. Karademir, C. Kocabas, and A. Aydinli, Absorption enhancement of molecules in the weak plasmon-exciton coupling regime, Opt. Lett. 39, 4994 (2014).
  * Lee  _et al._ [2004] J. Lee, A. O. Govorov, J. Dulka, and N. A. Kotov, Bioconjugates of CdTe nanowires and Au nanoparticles: Plasmon-exciton interactions, luminescence enhancement, and collective effects, Nano Letters 4, 2323 (2004).
  * Kumar  _et al._ [2012] A. Kumar, R. Srivastava, P. Tyagi, D. Mehta, and M. Kamalasanan, Efficiency enhancement of organic light emitting diode via surface energy transfer between exciton and surface plasmon, Organic Electronics 13, 159 (2012).
  * Koenderink [2010] A. F. Koenderink, On the use of Purcell factors for plasmon antennas, Opt. Lett. 35, 4208 (2010).
  * Tserkezis  _et al._ [2020] C. Tserkezis, A. I. Fernández-Domínguez, P. A. D. Gonçalves, F. Todisco, J. D. Cox, K. Busch, N. Stenger, S. I. Bozhevolnyi, N. A. Mortensen, and C. Wolff, On the applicability of quantum-optical concepts in strong-coupling nanophotonics, Rep. Prog. Phys. 83, 082401 (2020).
  * Törmä and Barnes [2014] P. Törmä and W. L. Barnes, Strong coupling between surface plasmon polaritons and emitters: a review, Rep. Prog. Phys. 78, 013901 (2014).
  * Zhao  _et al._ [2022] Q. Zhao, W.-J. Zhou, Y.-H. Deng, Y.-Q. Zheng, Z.-H. Shi, L. K. Ang, Z.-K. Zhou, and L. Wu, Plexcitonic strong coupling: unique features, applications, and challenges, J. Phys. D: Appl. Phys. 55, 203002 (2022).
  * Zengin  _et al._ [2015] G. Zengin, M. Wersäll, S. Nilsson, T. J. Antosiewicz, M. Käll, and T. Shegai, Realizing strong light-matter interactions between single-nanoparticle plasmons and molecular excitons at ambient conditions, Phys. Rev. Lett. 114, 157401 (2015).
  * Gates  _et al._ [2005] B. D. Gates, Q. Xu, M. Stewart, D. Ryan, C. G. Willson, and G. M. Whitesides, New approaches to nanofabrication: Molding, printing, and other techniques, Chemical Reviews 105, 1171 (2005).
  * Bellessa  _et al._ [2004] J. Bellessa, C. Bonnand, J. C. Plenet, and J. Mugnier, Strong coupling between surface plasmons and excitons in an organic semiconductor, Phys. Rev. Lett. 93, 036404 (2004).
  * Chikkaraddy  _et al._ [2016] R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess, and J. J. Baumberg, Single-molecule strong coupling at room temperature in plasmonic nanocavities, Nature 535, 127 (2016).
  * Stührenberg  _et al._ [2019] M. Stührenberg, S. Munkhbat, D. G. Baranov, J. Cuadra, A. B. Yankovich, T. J. Antosiewicz, E. Olsson, and T. Shegai, Strong light–matter coupling between plasmons in individual gold bi-pyramids and excitons in mono- and multilayer WSe2, Nano Lett. 18, 5938 (2019).
  * Geisler  _et al._ [2019] M. Geisler, X. Cui, J. Wang, T. Rindzevicius, L. Gammelgaard, B. S. Jessen, P. A. D. Gonçalves, F. Todisco, P. Bøggild, A. Boisen, M. Wubs, N. A. Mortensen, S. Xiao, and N. Stenger, Single-crystalline gold nanodisks on WS2 mono- and multilayers for strong coupling at room temperature, ACS Photonics 6, 994 (2019).
  * Hobson  _et al._ [2002] P. A. Hobson, W. L. Barnes, D. G. Lidzey, G. A. Gehring, D. M. Whittaker, M. S. Skolnick, and S. Walker, Strong exciton–photon coupling in a low-Q all-metal mirror microcavity, Applied Physics Letters 81, 3519 (2002).
  * Bellessa  _et al._ [2009] J. Bellessa, C. Symonds, K. Vynck, A. Lemaitre, A. Brioude, L. Beaur, J. C. Plenet, P. Viste, D. Felbacq, E. Cambril, and P. Valvin, Giant Rabi splitting between localized mixed plasmon-exciton states in a two-dimensional array of nanosize metallic disks in an organic semiconductor, Phys. Rev. B 80, 033303 (2009).
  * Tserkezis  _et al._ [2018] C. Tserkezis, M. Wubs, and N. A. Mortensen, Robustness of the Rabi splitting under nonlocal corrections in plexcitonics, ACS Photonics 5, 133 (2018).
  * Ekardt [1984] W. Ekardt, Dynamical polarizability of small metal particles: Self-consistent spherical jellium background model, Phys. Rev. Lett. 52, 1925 (1984).
  * Townsend and Bryant [2012] E. Townsend and G. W. Bryant, Plasmonic properties of metallic nanoparticles: The effects of size quantization, Nano Letters 12, 429 (2012).
  * Zuloaga  _et al._ [2010] J. Zuloaga, E. Prodan, and P. Nordlander, Quantum plasmonics: optical properties and tunability of metallic nanorods, ACS Nano 4, 5269 (2010).
  * Zuloaga  _et al._ [2009] J. Zuloaga, E. Prodan, and P. Nordlander, Quantum description of the plasmon resonances of a nanoparticle dimer, Nano Lett. 9, 887 (2009).
  * Teperik  _et al._ [2013b] T. V. Teperik, P. Nordlander, J. Aizpurua, and A. G. Borisov, Robust subnanometric plasmon ruler by rescaling of the nonlocal optical response, Phys. Rev. Lett. 110, 263901 (2013b).
  * Mie [1908] G. Mie, Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen, Ann. Phys. 330, 377 (1908).
  * Babaze  _et al._ [2022] A. Babaze, E. Ogando, P. E. Stamatopoulou, C. Tserkezis, N. A. Mortensen, J. Aizpurua, A. G. Borisov, and R. Esteban, Quantum surface effects in the electromagnetic coupling between a quantum emitter and a plasmonic nanoantenna: time-dependent density functional theory vs. semiclassical feibelman approach, Opt. Express 30, 21159 (2022).
  * Babaze  _et al._ [2023] A. Babaze, T. Neuman, R. Esteban, J. Aizpurua, and A. G. Borisov, Dispersive surface-response formalism to address nonlocality in extreme plasmonic field confinement, Nanophotonics 12, 3277 (2023).
  * Yan  _et al._ [2015] W. Yan, M. Wubs, and N. Asger Mortensen, Projected dipole model for quantum plasmonics, Phys. Rev. Lett. 115, 137403 (2015).
  * Fofang  _et al._ [2008] N. T. Fofang, T.-H. Park, O. Neumann, N. A. Mirin, P. Nordlander, and N. J. Halas, Plexcitonic nanoparticles: Plasmon-exciton coupling in nanoshell-J-aggregate complexes, Nano Letters 8, 3481 (2008).
  * Antosiewicz  _et al._ [2014] T. J. Antosiewicz, S. P. Apell, and T. Shegai, Plasmon–exciton interactions in a core–shell geometry: From enhanced absorption to strong coupling, ACS Photonics 1, 454 (2014).
  * Tserkezis  _et al._ [2008] C. Tserkezis, G. Gantzounis, and N. Stefanou, Collective plasmonic modes in ordered assemblies of metallic nanoshells, J. Phys. Condens. Matter 20, 075232 (2008).
  * Bohren and Huffman [1983] C. F. Bohren and D. R. Huffman,  _Absorption and Scattering of Light by Small Particles_ (Wiley, 1983).
  * Born and Wolf [1999] M. Born and E. Wolf,  _Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light_, 7th ed. (Cambridge University Press, 1999).


Proofs - Causal Function sec.dynamic.causal

We firstly provide Lemma lem: interpolation, which is fundamental to our
subsequent proofs. Lemma lem: interpolation is similar to the universal
approximation theorem (UAT)
haykin1998neural,hassoun1995fundamentals,pinkus1999approximation. We do not
directly use the UAT as it is for neural networks but our results are not
bounded to any specific learning paradigm or algorithm. But due to the UAT,
our results are applicable to learning in neural networks. The UAT guarantees
that $\sup_{x\in K}||f(x)-g(x)||<\eta$, where $g(x)$ is the target function
and $K$ is compact. Applying the UAT in our setting of Lemma lem:
interpolation, we have $\sup_{1\leq i\leq n}||f(x_{i})-y_{i}||<\eta$,
blackwhere $n$ is the number of training samples. The UAT has a condition that
$K$ is compact. In Lemma lem: interpolation, since we have a stronger
condition that $K=\\{(x_{i},y_{i})|\,1\leq i\leq n\\}$ is finite, we provide a
stronger result that $\sup_{1\leq i\leq n}||f(x_{i})-y_{i}||=0$.

[A Simple Interpolating Function] Let $(X,d_{X})$ be a metric space. Let
$(Y,||\cdot||)$ be a Banach space. For any $D\subseteq\\{(x,y)|\,x\in X,y\in
Y\\}$, if $x\neq x^{\prime},\,\forall\,(x,y),(x^{\prime},y^{\prime})\in D$,
then there exists a continuous approximation function $f:X\rightarrow Y$, s.t.
$f(x)=y,\,\forall\,(x,y)\in D$. [Proof of Lemma 1] Let

$\epsilon=\min_{(x,y),(x^{\prime},y^{\prime})\in D}d_{X}(x,x^{\prime}).$

Define

$K_{\epsilon}(x,x^{\prime})=\frac{\epsilon}{d_{X}(x,x^{\prime})}.$

Denote $D=\\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\\}$, let

$f(x)=\frac{\sum_{i=1}^{n}y_{i}K_{\epsilon}(x,x_{i})}{\sum_{i=1}^{n}K_{\epsilon}(x,x_{i})},\,x\in
X\backslash\\{x_{1},\dots,x_{n}\\}.$ (1)

For $\forall\,1\leq i_{0}\leq n$ ($i_{0}$ is any integer), since
$\lim_{x\rightarrow x^{\prime}}K_{\epsilon}(x,x^{\prime})=+\infty$,
$\sup_{d_{X}(x,x^{\prime})>\epsilon}K_{\epsilon}(x,x^{\prime})<1$ and

$f(x)=\frac{y_{i_{0}}K_{\epsilon}(x,x_{i_{0}})}{\sum_{i=1}^{n}K_{\epsilon}(x,x_{i})}+\frac{\sum_{i\neq
i_{0}}y_{i}K_{\epsilon}(x,x_{i})}{\sum_{i=1}^{n}K_{\epsilon}(x,x_{i})},$

it’s obvious that $\forall\,\eta>0$, $\exists\,\delta>0$, s.t.
$\forall\,d_{X}(x,x_{i_{0}})<\delta$, $||f(x)-y_{i_{0}}||<\eta$, which is
$\lim_{x\rightarrow x_{i_{0}}}f(x)=y_{i_{0}}$. Therefore, we can define $f$ on
$X$ and it’s obvious that $f$ is a continuous approximation function.
[Proof of Theorem LABEL:thm:_full_dynamic_function] In the theorem presented
in Section LABEL:sec:dp in the main text, we wrote $D=Y\times X$ for
simplicity. In detail, we have $D=\\{((y,x),f(y,x))|\,(y,x)\in Y\times X\\}$.
By Lemma 1, let $Y\times X$ be $X$ and $Y$ be $Y$, there exists $\Hat{f}$ s.t.
$\Hat{f}(y,x)=f(y,x),\,\forall\,(y,x)\in D$. Since $D=Y\times X$, the proof is
done.
[Proof of Corollary LABEL:thm:_not_full_dynamic_function] In the corollary
presented in Section LABEL:sec:dp in the main text, we wrote $D\neq Y\times X$
for simplicity. In detail, we have $D\neq\\{((y,x),f(y,x))|\,(y,x)\in Y\times
X\\}$. For simplicity, denote $A=\\{((y,x),f(y,x))|\,(y,x)\in Y\times X\\}$.
Since $D\neq A$, let $((y_{0},x_{0}),f(y_{0},x_{0}))\in A\backslash D$. Since
$|Y|>1$, let $y^{err}\in Y$ s.t. $y^{err}\neq f(y_{0},x_{0})$. Let

$D^{err}=D\cup\\{((y_{0},x_{0}),y^{err})\\}.$

By Lemma 1, let $Y\times X$ be $X$ and $Y$ be $Y$, there exists $\Hat{f}$ s.t.
$\Hat{f}(y,x)=y^{\prime},\,\forall\,((y,x),y^{\prime})\in D^{err}$. Since
$(y_{0},x_{0}),y^{err})\in D^{err}$, $\Hat{f}$ makes a wrong prediction as

$\Hat{f}(y_{0},x_{0})=y^{err}\neq f(y_{0},x_{0}).$

Since $D\subseteq D^{err}$, $\Hat{f}$ is correct on training dataset as

$\Hat{f}(y,x)=f(y,x),\,\forall\,((y,x),f(y,x))\in D.$

[Proof of Corollary LABEL:thm:_not_full_dynamic_function_inf_x] Denote
$D\subseteq\\{((y,x),f(y,x))|\,(y,x)\in Y\times X\\}$. For simplicity, denote
$A=\\{((y,x),f(y,x))|\,(y,x)\in Y\times X\\}$. Since $|X|=\infty$,
$|A|=\infty$. Since $|D|<\infty$, we know $D\neq A$ and $|A\backslash
D|=\infty$. Let $((y_{i},x_{i}),f(y_{i},x_{i}))\in A\backslash
D,\,i=0,1\dots,m$. Since $|Y|>1$, let $y_{i}^{err}\in Y$ s.t. $y^{err}_{i}\neq
f(y_{i},x_{i})$. Let

$D^{err}=D\cup\\{((y_{i},x_{i}),y_{i}^{err}),\,i=0,\dots,m)\\}.$

By Lemma 1, let $Y\times X$ be $X$ and $Y$ be $Y$, there exists $\Hat{f}$ s.t.
$\Hat{f}(y,x)=y^{\prime},\,\forall\,((y,x),y^{\prime})\in D^{err}$. Since
$((y_{i},x_{i}),y_{i}^{err})\in D^{err}$, $\Hat{f}$ makes $m+1$ mistakes as

$\Hat{f}(y_{i},x_{i})=y_{i}^{err}\neq f(y_{i},x_{i}),\,i=0,1,\dots,m.$

Since $D\subseteq D^{err}$, $\Hat{f}$ is correct on training dataset as

$\Hat{f}(y,x)=f(y,x),\,\forall\,((y,x),f(y,x))\in D.$

[Proof of Theorem LABEL:thm:_full_causal_function] By Lemma 1, let
$X^{\sup|s(v)|}$ be $X$ and $X$ be $Y$, there exists $\Hat{f}$ s.t.
$\Hat{f}(s(v))=v^{\prime},\,\forall\,(v^{\prime},s(v))\in D$. Since
$D=\\{(f(s(v)),s(v))|\,v\in V\\}$, the proof is done.
[Proof of Corollary LABEL:thm:_not_full_causal_function] In the corollary
presented in Section LABEL:sec:_dag in the main text, we wrote $D\neq
X^{\sup|s(v)|}$ for simplicity. In detail, we have
$D\neq\\{(f(s(v)),s(v))|\,v\in V\\}$. For simplicity, denote
$A=\\{(f(s(v)),s(v))|\,v\in V\\}$. Since $D\neq A$, let
$(f(s(v_{0})),s(v_{0}))\in A\backslash D$. Since $|\\{f(s(v)|\,s(v)\in
X^{\sup|s(v)|})\\}|>1$, let $v^{err}\in\\{f(s(v)|\,s(v)\in X^{\sup|s(v)|})\\}$
s.t. $v^{err}\neq f(s(v_{0}))$. Let

$D^{err}=D\cup\\{(v^{err},s(v_{0}))\\}.$

By Lemma 1, let $X^{\sup|s(v)|}$ be $X$ and $X$ be $Y$, there exists $\Hat{f}$
s.t. $\Hat{f}(s(v))=v^{\prime},\,\forall\,(v^{\prime},s(v))\in D^{err}$. Since
$(v^{err},s(v_{0}))\in D^{err}$, $\Hat{f}$ makes a wrong prediction as

$\Hat{f}(s(v_{0}))=v^{err}\neq f(s(v_{0})).$

Since $D\subseteq D^{err}$, $\Hat{f}$ is correct on training dataset as

$\Hat{f}(s(v))=f(s(v)),\,\forall\,(f(s(v)),s(v))\in D.$

[Proof of Corollary LABEL:thm:_not_full_causal_function_inf_x] Denote
$D\subseteq\\{(f(s(v)),s(v))|\,v\in V\\}$. For simplicity, denote
$A=\\{(f(s(v)),s(v))|\,v\in V\\}$. Since $\max(|X|,\sup|s(v)|)=\infty$,
$|X^{\sup|s(v)|}|=\infty$, $|A|=\infty$. Since $|D|<\infty$, we know $D\neq A$
and $|A\backslash D|=\infty$. Let $(f(s(v_{i}),v_{i}))\in A\backslash
D,\,i=0,1\dots,m$. Since $|\\{f(s(v)|\,s(v)\in X^{\sup|s(v)|})\\}|>1$, let
$v_{i}^{err}\in\\{f(s(v)|\,s(v)\in X^{\sup|s(v)|})\\}$ s.t. $v^{err}_{i}\neq
f(s(v_{i}))$. Let

$D^{err}=D\cup\\{(f(s(v_{i}),v_{i})),\,i=0,\dots,m)\\}.$

By Lemma 1, let $X^{\sup|s(v)|}$ be $X$ and $X$ be $Y$, there exists $\Hat{f}$
s.t. $\Hat{f}(s(v))=v^{\prime},\,\forall\,(v^{\prime},s(v))\in D^{err}$. Since
$(v_{i}^{err},s(v_{i}))\in D^{err}$, $\Hat{f}$ makes $m+1$ mistakes as

$\Hat{f}(s(v_{i}))=v_{i}^{err}\neq f(s(v_{i})),\,i=0,1,\dots,m.$

Since $D\subseteq D^{err}$, $\Hat{f}$ is correct on training dataset as

$\Hat{f}(s(v))=f(s(v)),\,\forall\,(f(s(v)),s(v))\in D.$

##  Proofs - Recursive Formula

[Proof of Theorem LABEL:thm:_recursive_markov_dp_solved] By Theorem
LABEL:thm:_full_dynamic_function, there exists $\Hat{f}$ s.t.
$\Hat{f}(y,x)=f(y,x),\,\forall\,(y,x)\in Y\times X$. Since $\tau$ is a Markov
dynamic process, by the definition in Eq. eq: def_dp, we have

$\tau_{f}(x_{[0:\infty]})=(y_{1},y_{2},\dots),\ \text{where}\
\left\\{\aligned&y_{1}=f(y_{0},x_{0}),\\\ \dots\\\ y_{n+1}=f(y_{n},x_{n}).\\\
\right.$

By induction, we have

$\aligned y_{1}=f(y_{0},x_{0})&=\Hat{f}(y_{0},x_{0})=\Hat{y}_{1},\\\ \dots\\\
y_{n+1}=f(y_{n},x_{n})=\Hat{f}(\Hat{y}_{n},x_{n})=\Hat{y}_{n+1}.\\\ $

[Proof of Theorem LABEL:thm:_resursive_dag_solved] By Theorem
LABEL:thm:_full_causal_function, there exists $\Hat{f}$, s.t.
$\Hat{f}(s(v))=f(s(v)),\,\forall\,(f(s(v)),s(v))\in\\{(f(s(v)),s(v))|\,v\in
V\\}$. By the definition in Eq. eq: def_dag in Section LABEL:sec:_dag, for an
arbitrary topological ordering of the graph $G$, we have

$G_{f}(\\{v_{i},i\leq|G||\,d(v_{i})=0\\})=(v_{1},\dots,v_{|G|}),\text{where}\
\left\\{\aligned&v_{1}=f(s(v_{1})),\\\ \dots\\\ v_{n}=f(s(v_{n})).\right.$

Given the graph $G$ and the topological ordering, by induction, we have

$\aligned v_{1}=f(s(v_{1}))&=\Hat{f}(s(v_{1}))=\Hat{v}_{1},\\\ \dots\\\
v_{n}=f(s(v_{n}))=\Hat{f}(s(\Hat{v}_{n}))=\Hat{v}_{n}.$

##  Proof - Maximal Input Element Distance of a Reasoning step

[9]r2.8in  Two examples of the ko problem.  [width=0.8]fig/example_ko.pdf
Before presenting the proof for Theorem
LABEL:thm:_find_input_vertices_is_possible, we first provide an example to
show that a $4R$-length window is not enough to decide which elements should
be used in the next reasoning or calculation step. Let’s consider a simple
one-dimensional ko problem. Let $s=\dots 001001111\dots$ be a sequence of
$0$’s and $1$’s. We say $s_{i}$ is captured if $s_{i}\neq s_{i-1},\,s_{i}\neq
s_{i+1}$. We say $s_{i}$ is not a ko if both $s_{i-1}$ and $s_{i+1}$ are not
captured. We say $s_{i}$ is a ko if $s_{i}$ is captured and at least one of
$s_{i-1}$ and $s_{i+1}$ is captured. When some $s_{i}$ is captured and is not
a ko, we set $s_{i}=1-s_{i}$. We continue the process until $s$ won’t change,
which is the final settlement of the problem. In this problem, if $s_{i}$ is
captured and is not a ko, it should be acted upon next (changing $s_{i}$’s
value). The neighbors $s_{i-1}$ and $s_{i+1}$ are enough to decide the value
of $s_{i}$, i.e. $\max_{v^{\prime}\in(s\circ
t)(s_{i})}d(v^{\prime},s_{i})=\max(d(s_{i-1},s_{i}),d(s_{i+1},s_{i}))=1$.
Therefore, $R=1$ in this problem. However, a $4R$-length window is not enough
to decide whether $s_{i}$ should be reasoned or acted upon next. For instance,
let $s=11010$. In this problem, $s_{[2]}=0$ is captured and is a ko, so it
cannot be changed. Therefore, when we only consider $4$ elements
$s_{[0:4]}=1101$, we will receive the ground truth label $y_{[0:4]}=0000$,
where $0$ indicates the element shouldn’t be changed and $1$ indicates the
element should be changed. However, when we consider a problem
$s^{\prime}=11011$, $s_{[2]}^{\prime}=0$ is captured and is not a ko, so it
should be changed. Therefore, when we only consider $4$ elements
$s_{[0:4]}^{\prime}=1101$, the ground truth label is
$y_{[0:4]}^{\prime}=0010$. In these two problems, they share an identical
$4R$-length window $1101$, but the labels of this window are different, which
is obviously ambiguous. [Proof of Theorem
LABEL:thm:_find_input_vertices_is_possible] Denote the training dataset as

$D\subseteq\\{(\Tilde{g}(s),s)|s=s_{0}s_{1}\dots s_{4R},s=g^{-1}(G)\\},$

where $\Tilde{g}(s)=1_{0}(s_{0})1_{1}(s_{1})\dots 1_{4R}(s_{4R})$ is a
$(4R+1)$-length binary sequence that $1_{i}(s_{i})$ indicates whether $s_{i}$
should be reasoned next. By Lemma 1, let $X^{4R+1}$ to be $X$ and $2^{4R+1}$
to be $Y$, there exists $\Hat{f}$ s.t. $\Hat{f}(s)=I,\,\forall\,(I,s(v))\in
D$. Since $D=\\{(\Tilde{g}(s),s)|s=s_{0}s_{1}\dots s_{4R},s=g^{-1}(G)\\}$, the
proof is done.

Table :  Experimental data examples

[arctan] 1.1  Input[0]: a, b Output[0]: arctan(a/b)

[arithmetic in prime field $F_{7}$] 1.1  Input[0]: (0+4-(2-3*6))*(4+0)
Output[0]: ( 4 -(2- 4 ))* 4 Input[1]: (4-(2-4))*4 Output[1]: (4- 5 )*4
Input[2]: (4-5)*4 Output[2]: 6 *4 Input[3]: 6*4 Output[3]: 3

1.1 [1-line addition]  Input[0]: 285+9805=? Output[0]: 285+9805=c0 Input[1]:
285+9805=c0 Output[1]: 285+9805=?90 Input[2]: 285+9805=?90 Output[2]:
285+9805=c090 Input[3]: 285+9805=c090 Output[3]: 285+9805=10090 [1-line
multiplication] 1.1  Input[0]: 1*3=? Output[0]: 1*3=1+? Input[1]: 1*3=1+?
Output[1]: 1*3=1+1+? Input[2]: 1*3=1+1+? Output[2]: 1*3=1+1+1 Input[3]:
1*3=1+1+1 Output[3]: 1*3=2+1 Input[4]: 1*3=2+1 Output[4]: 1*3=3

[3-line addition] 1.1  Input[0]: 89283 3360 ? Output[0]: ?3 Input[1]: 89283
3360 ?3 Output[1]: c43 Input[2]: 89283 3360 c43 Output[2]: ?643 Input[3]:
89283 3360 ?643 Output[3]: c2643 Input[4]: 89283 3360 c2643 Output[4]: 92643

##  Experimental Data and Implementation Details

The arithmetic in $F_{7}$ problem is defined in the prime field $F_{7}$. The
1-line addition problem and the 3-line addition problem use ‘?’ to represent
$0$ being carried from the right and ‘c’ to represent $1$ being carried from
the right. In all problems, * is equivalent to $\times$. We use * instead of
$\times$ for ease of aligning chars. All the problems are trained with $50k$
batches. Each batch contains $256$ CoT steps. For each problem, we first
randomly generate a question and then its detailed CoT steps. Each CoT step is
a pair (Input[i], Output[i]), as shown in Table 1. We put the CoT steps into a
batch until the batch size reaches $256$. When testing the performance after
training, we test $6$ different datasets for each problem, which is shown in
Table LABEL:tab:experiment_setting. Each testing dataset is independently
generated with $1k$ questions. We solve each question by CoT using the trained
model. For the arctan problem, since it only has one step, we only inference
in one step. For the arithmetic problem, we stop the CoT output generation in
a step when (i) the output of the step has only one token/element, or (ii) the
output of the step is identical to the input of the step. For 1-line addition
and 3-line addition, we stop the CoT output generation if (i) ‘?’ and ‘c’
aren’t shown in the output, or (ii) the number of CoT steps is greater than
the number of digits. For the 1-line multiplication problem, we stop the CoT
output generation if (i) ‘+’ is not shown in the output, or (ii) the number of
CoT steps is greater than the multipliers. For the arctan problem, the
predicted value is considered correct if the absolute error is smaller than
$0.01$. For the other problems, the final output is considered correct only
when it is identical to the ground truth. For the dataset of each problem, the
accuracy is the number of correctly answered questions divided by the total
number of questions. The model of the arctan problem has 3 fully connected
layers. The model of the arithmetic problem, the 1-line addition problem, the
3-lines addition problem and the 1-line multiplication problem has 3
Transformer encoders with relative position embedding. The optimizer is Adam
and the learning rate is $0.0001$. The training data for each task contains
$12.8$M CoT steps and was trained for $1$ epoch.

##  Related Empirical Work

Evaluations and limitations of LLM reasoning. Continuing with the discussion
about evaluations of the reasoning capabilities of LLMs in Section LABEL:sec-
intro, we present a more extensive literature survey here. Note that Sections
LABEL:sec-intro and LABEL:sec.overview have discussed empirical works about
length generalization [anil2022exploring, dziri2023faith, zhang2022unveiling].
In these papers, the authors also tried to mitigate the problem through
improved training [anil2022exploring], improved promoting and fine-tuning of
LLMs [zhang2022unveiling], and curriculum learning [abbe2023generalization].
However, none of them studied the length generalization problem theoretically
as we do. Below, we focus on surveying other empirical works. Many of them
identified limitations of LLMs in solving different reasoning problems, but
few have characterized the limitations in a formal manner to facilitate
theoretical investigation. [meadows2023symbolic] created a dataset
specifically for mathematical reasoning that can be perturbed. They showed
that perturbations of the tasks heavily affect the results, reducing F1 score
from 97% to 17%, which suggests that inference is likely to be dominated by
surface-level patterns unrelated to the deeper understanding of the
mathematical operators. However, this evaluation was done using only BERT
[devlin2018bert] based models, but not on more recent LLMs like ChatGPT and
GPT4. [wu2023reasoning] used “counterfactual” tasks that deviate from the
standard reasoning tasks to evaluate LLMs. It was found that the performance
degrades substantially compared to the default conditions, which again
suggests that while LLMs can perform reasoning to some extend, they often rely
on narrow, non-transferable procedures or surface patterns for task-solving. A
counterfactual based evaluation was also done in li2023counterfactual, which
reached the same conclusion. [liu2023evaluating] evaluated ChatGPT and GPT-4
on logical reasoning. The results showed that they do relatively well on well-
known public domain datasets, but their performances drop substantially when
newly released and out-of-distribution datasets are used. [xu2023large] too
evaluated LLMs using logical reasoning (deductive, inductive, abductive and
mixed-form reasoning) and gave pros and cons of LLMs. [she2023scone] created a
dataset for reasoning involving negations and evaluated LLMs and showed poor
results. [ando2023evaluating] created a dataset, originally designed for
psychological experiments to assess human logical abilities in syllogistic
reasoning. The authors examined three types of biases observed in human
syllogistic reasoning: belief biases, conversion errors, and atmosphere
effects. The evaluation on LLMs showed that they struggle with problems
involving these biases too. [tan2023towards] created a dataset to evaluate
LLMs on temporal reasoning and showed some weaknesses of LLMs. They then
proposed an approach to improve the results. Chain of thoughts (CoT) and
variants. Earlier prompting for solving reasoning problems using LLMs only
states the question and the answer. They found that these two pieces of
information are insufficient for LLMs to learn to perform effective reasoning,
which we have proved. Then chain of thought (CoT) prompting wei2022chain was
proposed to improve the situation. CoT basically contains the detailed
intermediate reasoning steps between the question and the answer for fine-
tuning the LLMs, which significantly enhance LLMs’ reasoning capabilities
[chung2022scaling, hsieh2023distilling, mukherjee2023orca, fu2023chain].
[saparov2022language] created a synthetic dataset generated based on first-
order logic. They then parsed the generated CoT into symbolic proofs for
formal analysis. It was shown that LLMs are capable of reasoning. The success
of CoT has encouraged researchers to refine the technique and also propose
variations of the technique. For example, [chen2023measuring] proposed a
metric to measure the effectiveness of CoT and a technique to improve CoT for
vision-language models. [wang2023making] studied using multiple reasoning
paths and positive and negative answers to improve CoT reasoning.
[zhang2023cumulative] proposed cumulative reasoning, which employs LLMs in a
cumulative and iterative manner to emulate the human thought process.
[qi2023art] proposed a divide-and-conquer algorithm that simulates the self-
questioning and recursive thinking process of humans to improve CoT.
[wang2023learning] investigated how to incorporate into relatively small LMs
the capabilities of multi-step reasoning and CoT. [wang2022towards] found that
even logically invalid CoT also helps reasoning. This was confirmed in
schaeffer2023invalid. To deal with unsound inferences, [poesia2023certified]
introduced a class of external tools for LLMs called guides that use states
and incremental constraints to guide the generation in reasoning. A related
work on using external tools was done in [xu2023rewoo]. [wang2022self]
improved CoT using multiple-paths and consistency check. [ling2023deductive]
studied the verification of CoT. [stolfo2023understanding] identified part of
an LLM responsible for reasoning. In a different direction, yang2022chain
argued that the prevailing approach to CoT prompt selection through trial and
error is unsatisfactory. They then proposed a principled approach for multi-
domain LLM CoT prompt selection. Several researchers also broadened the CoT
method and proposed the neural symbolic code prompting hu2023code, program of
thoughts chen2022program,cheng2023binding, tree-of-thoughts
yao2023tree,long2023large, tree-of-mixed-thoughts hu2023tree, tree of
uncertain thoughts mo2023tree, hypergraph-of-thoughts yao2023thinking,
recursion of thoughts lee2023recursion, chain of knowledge wang2023boosting,
chain of simultaneous thoughts shao2022chaining, graph-of-thoughts
yao2023beyond and faithful chain of thoughts lyu2023faithful. Further,
[bi2023program] proposed a complexity measure and chose the optimal complexity
to improve the program of thoughts chen2022program. [wang2023exploring]
proposed a method to improve the generation of equations from natural language
questions as the intermediate step to answer the original question.
[gao2023pal] combined CoT and Program-Aided Language Models (PAL) for improved
reasoning.