# Influence of the electron spill-out and nonlocality on gap-plasmons in the
limit of vanishing gaps

M. Khalid<EMAIL_ADDRESS>Istituto Italiano di Tecnologia, Center for
Biomolecular Nanotechnologies, Via Barsanti 14, 73010 Arnesano, Italy.  O.
Morandi  Dipartimento di Matematica e Informatica U. Dini, University of
Florence, via S. Marta 3, 50149, Florence, Italy  E. Mallet  Clermont
Université, Institut Pascal (IP), BP 10448, F-63000 Clermont-Ferrand, France
CNRS, UMR 6602, IP, F-63171 Aubière, France  P. A. Hervieux  Université de
Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg,
UMR 7504, F-67000 Strasbourg, France  G. Manfredi  Université de Strasbourg,
CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504,
F-67000 Strasbourg, France  A. Moreau  Clermont Université, Institut Pascal
(IP), BP 10448, F-63000 Clermont-Ferrand, France  CNRS, UMR 6602, IP, F-63171
Aubière, France  C. Ciracì<EMAIL_ADDRESS>Istituto Italiano di
Tecnologia, Center for Biomolecular Nanotechnologies, Via Barsanti 14, 73010
Arnesano, Italy.

(September 2, 2024)

###### Abstract

We study the effect of electron spill-out and of nonlocality on the
propagation of light inside a gap between two semi-infinite metallic regions.
We first present a simplified physical model for the spill-out phenomenon, an
approach sufficient to show that the propagation of the gap-plasmon becomes
impossible in the tunneling regime. However, in the limit of very small gaps,
only a Quantum Hydrodynamic Theory (QHT) approach, taking into account both
the electron spill-out and nonlocality, is able to accurately model the gap-
plasmon characteristics and to correctly retrieve the refractive index of the
bulk metal as the limit of the effective index of the gap-plasmon for
vanishing gaps. Finally, we analyze the relation between different models and
show that up to a certain size it is possible to predict the correct gap-
plasmon effective index by considering a properly resized effective gap.

##  I Introduction

The possibility to design extremely miniaturized devices, even optical
circuits, has fostered research in plasmonics for decades. Surface plasmons
can in fact propagate with very large wavevectors and effective indexes at
frequencies right below the plasma frequency. Indeed, the typical size of a
resonator is roughly half the effective wavelength – the wavelength in vacuum
divided by the effective index – of the guided mode responsible for the
resonance. Ohmic losses in metals however hinder the propagation of the guided
mode so that, in practice, surface plasmons can only reach very modest
effective index.

It has been realized quite early [1, 2, 3, 4] that reaching larger effective
index requires for light to propagate inside metal as much as possible,
because of the plasmonic drag metals provide [5]. Especially, plasmonic
structures exhibiting narrow dielectric gaps, known in the literature as
metal-insulator-metal (MIM) structures, support a fundamental mode with a
particularly large effective index called the gap-plasmon [6, 7]. The
excitation of such modes in finite systems thus leads to resonances in
exceptionally small volumes, while the quality factor of the resonance [8] and
its cross-section [9] are essentially preserved. Gap-plasmon resonators have
had an enormous impact on the field of plasmonics [10]. For instance, since
reduced mode volumes lead to a large Purcell effect, gap-plasmon resonators
have allowed to reach record Purcell effects [11, 12], making spontaneous
emission faster than stimulated emission.

The dielectric-filled gaps can be so small (typically smaller than 1 nm [13])
that many questions arise concerning the accuracy of an electromagnetic
description at such a scales. Indeed, at these scales the classical local
response approximation (LRA), based on the Drude model, reaches its limits
[14], as spatial dispersion due to the repulsion between electrons inside the
metal has a direct impact on the effective index of the gap-plasmon [15, 16].
Electron-electron interactions can be taken into account in a simple way in
the framework of the hydrodynamic theory (HT) with hard-wall boundary
conditions (i.e., no electron spill-out is allowed), by adding the Thomas-
Fermi (TF) electron pressure term [17, 18]. In the LRA the effective index of
the gap-plasmon diverges in the limit of vanishing gaps. The TFHT, while
rectifying the divergence, nevertheless predicts a monotonically increasing
behavior of effective index as gap closes. This is not quite satisfactory
since in this limit, one would expect to retrieve the bulk characteristics,
i.e. no propagating mode at all.

When dealing with the near-atomic sub-nanometer length scales, it becomes
inevitable to consider the quantum nature of electrons, and the quantum
mechanical effects such as electron spill-out and quantum tunnelling become
important.

In the case of the gap plasmon, a Local Response Approximation [19] taking
into account the spill out is possible, but is characterized by unexpected
features in the electric field [20, 21]. The TFHT cannot address these
phenomena properly and, therefore, more sophisticated theories are required
for an accurate description of electron dynamics [22, 23, 24]. Ab-initio
quantum mechanical approaches provide an appealing solution to such problems
although they are often computationally complex and demanding, especially for
large plasmonic structures [25, 26, 27]. However, in the past few years, a
hybrid approach called quantum hydrodynamic theory (QHT)—which adds a $\nabla
n$-dependent correction ($n$ being the charge density) to the TF kinetic
energy—has been developed to study plasmonic systems more accurately [28, 29,
30, 31, 32].

Here we study the impact of the spill-out of the electron gas outside of metal
and of nonlocality induced by the repulsion between electrons on the
propagation of the gap-plasmon in the vanishing gap limit. More precisely, we
show that the spill-out is the phenomenon hindering the effective index from
diverging but that the most surprising features associated to a purely local
approach tend to disappear when nonlocality is taken into account. The correct
limit, as the gap closes, is thus retrieved only when both the spill-out and
nonlocality are both taken into account.

In section II, we detail the theoretical background of our study, by deriving
simplified equations describing the propagation of the gap-plasmon in the
presence of spill-out, for the sake of the physical analysis, and then
detailing the QHT framework. In section III, we discuss our numerical results
by comparing the behaviour of the effective index for the different approaches
of interest (LRA, TFHT, simplified spill-out model and QHT). Finally, in
section IV, we discuss the implications of our work regarding the definition
of the actual optical edges of a metal.

##  II Theoretical framework

The propagation of the electric and magnetic fields, $\bf E$ and $\bf H$, of
the gap-plasmon is described by Maxwell’s equations:

$\displaystyle\begin{aligned} &\nabla\times{\bf E}=-\partial_{t}\mu_{0}{\bf
H},\\\ &\nabla\times{\bf
H}=\varepsilon_{0}\varepsilon_{\infty}\partial_{t}{\bf E}+{\bf
J},\end{aligned}$ (1)

where $\mu_{0}$ and $\varepsilon_{0}$ are the vacuum permeability and
permittivity, respectively; $\varepsilon_{\infty}$ is the relative
permittivity accounting for any dielectric local response. $\bf J$ is the
current inside the electron gas. This current is given by ${\bf J}=-en{\bf u}$
where $\bf u$ is the velocity of the electron gas, considered as a fluid, $e$
is the electron charge in absolute value, and $n$ is the electron density. In
this framework the velocity satisfies the following equations [33, 34]:

$\displaystyle\begin{aligned} &\partial_{t}n+{\nabla}\cdot(n{\bf u})=0\\\
&\partial_{t}{\bf u}+\gamma{\bf u}+{\bf u}\cdot{\bf\nabla}{\bf u}=\\\
&\quad=-\dfrac{e}{m}({\bf E}+{\bf u}\times\mu_{0}{\bf
H})-\frac{1}{m}\nabla\frac{\delta G[n]}{\delta n},\end{aligned}$ (2)

where $m$ is the electron mass and the energy functional $G[n]$ contains all
the internal energy of the electronic system (here the operator
$\frac{\delta}{\delta n}$ indicates the functional derivative with respect to
the density $n$). The last term in the second equation gives rise to a
nonlocal contribution to the current $\bf J$. Using the previous equations, we
obtain the following evolution equation for the current:

$\displaystyle\begin{aligned} \partial_{t}{\bf J}+\gamma{\bf
J}&-\frac{1}{e}\left(\frac{\bf J}{n}{\bf\nabla}\cdot{\bf J}+{\bf
J}\cdot{\bf\nabla}\frac{\bf J}{n}\right)=\\\ &=\frac{ne^{2}}{m}{\bf
E}-\frac{\mu_{0}e}{m}{\bf J}\times{\bf H}+\frac{en}{m}\nabla\frac{\delta
G[n]}{\delta n},\end{aligned}$ (3)

Neglecting higher order terms and taking $n({\bf r},t)=n_{0}({\bf
r})+n_{1}({\bf r},t)$, and ${\bf E}({\bf r},t)={\bf E}_{0}({\bf r})+{\bf
E}_{1}({\bf r},t)$, with $n_{0}$ being the equilibrium charge density and
${\bf E}_{0}$ the electrostatic field, we first obtain an equation for the
dynamic part,

$\displaystyle\begin{aligned} &(\partial_{t}+\gamma){\bf
J}=\frac{e^{2}}{m}\left(n_{0}{\bf E}_{1}+n_{1}{\bf E}_{0}\right)+\\\
&\quad\quad+\frac{e}{m}\left[n_{0}\nabla\left(\frac{\delta G[n]}{\delta
n}\right)_{1}+n_{1}\nabla\left(\frac{\delta G[n]}{\delta
n}\right)_{0}\right],\end{aligned}$ (4)

where $\left(\frac{\delta G[n]}{\delta n}\right)_{0}$ and $\left(\frac{\delta
G[n]}{\delta n}\right)_{1}$ are the zero-th and first order terms of the
potential, respectively. Second, we get an equation linking the static
electron density to the static electric field,

$\nabla\left(\frac{\delta G[n]}{\delta n}\right)_{0}+e{\bf E}_{0}=0,$ (5)

where ${\bf E}_{0}$ and $n_{0}$ can be found using Gauss’s law:

$\nabla\cdot\left(\varepsilon_{\infty}{\bf
E}_{0}\right)=-\frac{e}{\varepsilon_{0}}(n_{0}-n^{+}),$ (6)

$n^{+}$ being the positive background density.

Finally, combining Eqs. (1) and assuming a harmonic field oscillating at
$\omega$ we obtain:

$\displaystyle\nabla\times\nabla\times{\bf
E}_{1}-\varepsilon_{\infty}k_{0}^{2}{\bf E}_{1}=i\omega\mu_{0}{\bf J}$ (7)
$\displaystyle(\gamma-i\omega){\bf J}=\frac{e^{2}}{m}\left(n_{0}{\bf
E}_{1}+n_{1}{\bf E}_{0}\right)+$
$\displaystyle\quad\quad+\frac{e}{m}\left[n_{0}\nabla\left(\frac{\delta
G[n]}{\delta n}\right)_{1}+n_{1}\nabla\left(\frac{\delta G[n]}{\delta
n}\right)_{0}\right],$ (8)

where $k_{0}=\omega/c$.

Combining Eqs. (5) and (6), we obtain:

$\displaystyle\begin{aligned}
&\nabla\left[\varepsilon_{\infty}\cdot\nabla\left(\frac{\delta G}{\delta
n}\right)_{0}\right]-\frac{e^{2}}{\varepsilon_{0}}\left(n_{0}-n^{+}\right)=0,\end{aligned}$
(9)

Note also that using Eq. (5) we can simplify Eq. (8) as:

$\displaystyle\begin{aligned} &(\gamma-i\omega){\bf
J}=\frac{e^{2}}{m}n_{0}{\bf E}_{1}+\frac{en_{0}}{m}\nabla\left(\frac{\delta
G[n]}{\delta n}\right)_{1}.\end{aligned}$ (10)

###  II.1 Simplified local equations

We first begin with a simplified model, which is not fully accurate when
compared to QHT, but it allows to better grasp the physics at play since the
dispersion relation of the gap-plasmon is essentially analytical. We begin by
neglecting the pressure term, i.e. $G[n]\simeq 0$, such that our description
is purely local. The impact of nonlocality in the case of vanishing gaps (and
hard-wall boundaries) has already been investigated in Ref. [16], where it has
been shown that the effective index of the gap-plasmon increases although
staying finite (contrarily to the diverging behaviour of the LRA case [3,
35]). Yet, it does not tend to the expected value of the bulk metal refractive
index.

Starting from Eq. (10) we get

$\displaystyle\begin{aligned} &{\bf
J}=\frac{1}{\gamma-i\omega}\dfrac{e^{2}}{m}n_{0}{\bf E}_{1}.\end{aligned}$
(11)

The current can always be incorporated into Eqs. (1) as an effective
polarization $\textrm{P}_{f}$, by simply writing
$\partial_{t}\textrm{P}_{f}=\textrm{J}$. This finally allows to introduce a
local susceptibility by writing, just as for a Drude model, that

$\mathbf{P}_{f}=-\varepsilon_{0}\frac{\frac{e^{2}n_{0}(\mathbf{r})}{\varepsilon_{0}m}}{\omega^{2}-i\gamma\omega}{\mathbf{E}}_{1}=\varepsilon_{0}\chi_{f}{\mathbf{E}}_{1},$
(12)

except that the susceptibility is proportional to the electron density and
thus depends on the position here.

Now if we consider the propagation of a guided mode in the $z$ direction, and
assume that all fields are invariant in the $x$ direction, Maxwell’s equations
in the $p$-polarization reduce to:

$\displaystyle\partial_{y}E_{z}-\partial_{z}E_{y}=i\omega\mu_{0}\mu H_{x},$
(13)
$\displaystyle\partial_{z}H_{x}=-i\omega\varepsilon_{0}\varepsilon(y)E_{y},$
(14) $\displaystyle-$
$\displaystyle\partial_{y}H_{x}=-i\omega\varepsilon_{0}\varepsilon(y)E_{z},$
(15)

where $\varepsilon(y)=\varepsilon_{\infty}+\chi_{f}(y)$ is the local
permittivity.

Combining the previous equations to get rid of $H_{x}$ we obtain,

$\displaystyle\partial_{y}\partial_{z}E_{z}-\partial^{2}_{z}E_{y}$
$\displaystyle=i\omega\mu_{0}\partial_{z}H_{x}=k_{0}^{2}\varepsilon(y)E_{y}$
(16) $\displaystyle\partial_{y}^{2}E_{z}-\partial_{y}\partial_{z}E_{y}$
$\displaystyle=i\omega\mu_{0}\partial_{y}H_{x}=-k_{0}^{2}\varepsilon(y)E_{z}$
(17)

Since we are looking for a guided mode we take $\partial_{z}=ik_{z}$ so that
we obtain

$(k_{z}^{2}-\varepsilon(y)k_{0}^{2})E_{y}=-ik_{z}\partial_{y}E_{z},$ (18)

and finally manage to keep only the $E_{z}$ field to write the wave equation
under the form

$\partial_{y}^{2}E_{z}+k_{0}^{2}\,\varepsilon(y)\,E_{z}+\partial_{y}\left(\frac{k_{z}^{2}}{\varepsilon(y)\,k_{0}^{2}-k_{z}^{2}}\partial_{y}E_{z}\right)=0.$
(19)

We underline that such an approach is not self-consistent, as it requires to
choose arbitrarily the electron density. The choice of electron density put
aside, our simplified model is essentially identical to the model used in [20]
and yields almost identical results. It is close to a well known approach
[19].

We have then converted this equation into a numerical problem using a finite
difference scheme, in order to look for values of $k_{z}$ in the complex plan
which allows to make the determinant of the finite difference matrix vanish.
As such a problem is extremely unstable numerically, obviously because of the
different scales involved (the gap being much smaller than the skin depth), we
have used the exponentially decaying analytic solution inside the metal. This
solution is known when $\varepsilon(y)$ is uniform and corresponds to the
classical solution of Maxwell’s equations considering Drude’s model for the
metal. Connecting the finite difference problem with the analytic solution
allows to reduce the size of the problem and to improve the stability, finally
allowing to compute the wave vector of the gap-plasmon, even for vanishing
gaps. We compared it with a finite element based approach, and both are in
almost perfect agreement.

###  II.2 Self-consistent quantum hydrodynamic theory

Equation (9) with Eqs. (7) and (10) constitute the basis of the QHT [30, 28]
linear response and allow to calculate the detailed electron gas dynamic at
the metal surface. In this work, we employ the following approximation:

$\displaystyle G[n]=T_{\textrm{TF}}[n]+\lambda_{\rm
vW}T_{\textrm{vW}}[n,\nabla n]+E_{\textrm{XC}}[n]$ (20)

where $T_{\rm TF}$ is the TF kinetic energy, $T_{\rm vW}$ is the gradient-
dependent correction, namely the von Weizsäcker term, to the TF functional,
which allows to take into account effects depending on the gradient of the
electron density $n$ and thus should not be neglected when spatially dependent
densities are considered. $E_{\rm XC}$ is the exchange-correlation energy
functional. The parameter $\lambda_{\rm vW}$, weighting the von Weizsäcker
contribution, is usually taken in the literature as, $1/9\leq\lambda_{\rm
vW}\leq 1$. In simple words, this parameter defines the distance an electron
can travel away from the metal surface during the time of an optical cycle.
Thus, choosing $\lambda_{\rm vW}=1$ signifies a higher electron spill-out and
slower decay of $n_{0}$ from a metal surface, whereas $\lambda_{\rm vW}=1/9$
gives a lower spill-out and a faster decay of the charge density. The explicit
expressions of $G[n]$ in Eq. (20) and its functional derivative, i.e.,
$\left(\frac{\delta G}{\delta n}\right)_{1}$, can be found in Ref. 30.
Considering a spatially-constant $n_{0}=n_{b}$, $n_{b}$ being the charge
density of the bulk metal, and $G[n]=0$ leads to the usual Drude-type relation
(LRA) for the polarization, and using $T_{\textrm{vW}}=E_{\textrm{xc}}=0$
returns $G[n]=T_{\rm TF}$, the standard hydrodynamic theory in the Thomas-
Fermi approximation (TFHT).

The QHT is highly dependent on the spatially varying equilibrium charge
density $n_{0}$ and, consequently, approximation of the optical response of a
nanoplasmonic system strongly relies on the description of $n_{0}$. Over the
past few years, QHT has been applied to probe linear plasmonic properties of a
variety of individual nanoparticles as well as nearly-touching nanostructures
(i.e., in the tunnelling regime) [36, 37, 38] and a fairly good agreement
between the QHT and time-dependent density functional theory calculations has
been reported [36]. Very recently, electron spill-out effects in singular
metasurfaces have been analyzed using QHT, showing that the spill-out
effectively blunts the sharp singularities [39]. QHT has also been applied to
study nonlinear optical properties of thick metal films and it has been shown
that there exist spill-out induced resonances which can be exploited to
achieve very large second-harmonic generation efficiencies [40]. These
resonances cannot be excited in the LRA and TFHT, as electron spill-out is
overlooked in these approximations. Therefore, it becomes important to take
into account the nonlocal and quantum effects when dealing with plasmonic
structure with sub-nanometer geometrical features. Despite the fact that QHT
gives a very good prediction of the plasmon resonances by adding a $\nabla
n-$dependent correction to the kinetic energy, however, it also suffers some
limitations. For example, it predicts unexpected modes between the surface and
bulk plasmon frequencies lying in the far tail region of the exponentially
decaying charge density. Very recently, a development in the QHT theory based
on the Laplacian-level kinetic energy functionals has been presented [41]. The
introduction of Laplacian-dependent charge density results in more robust
numerical solutions but its implementation becomes more complex. Since in the
present study we only consider very small gaps, the tail of the charge density
does never extends to the problematic region. Therefore, we expect that the
QHT approach in the limit of the von Weizsäcker approximation used herein
gives us results with the same level of accuracy and with the advantage of
avoiding additional numerical complexity, arising from the Laplacian-level
functionals. We adopt the jellium approximation [42] in implementation of the
QHT, which assumes that the electrons in a metal are confined by a constant
positive background charge $n^{+}=(\frac{4}{3}\pi r_{s}^{3})^{-1}$, where
$r_{s}$ indicates the Wigner-Seitz radius ($r_{s}=4$ for Na and $r_{s}=3$ for
Ag).

In order to compute the gap-plasmon mode and its dispersion relation we first
compute the space-dependent equilibrium charge density using Eq. (9) and then
solve Eqs. (7) and (10) assuming a solution of the type
$\mathbf{E}(\mathbf{r})=\tilde{\mathbf{E}}(y)\,e^{ik_{\rm gp}z}$ (with similar
expressions for $\mathbf{J}$), where $k_{\rm gp}$ is the propagation constant
of the gap plasmon mode. For a given angular frequency $\omega$, we solve this
system of equations using the finite-element method to compute the associated
propagation constant $k_{\rm gp}(\omega)$, as well as the mode
$\tilde{\mathbf{E}}(y)$.

##  III Results

Let us consider a MIM configuration as schematically shown in Fig. 1, that is,
a dielectric gap $g$ sandwiched between two semi-infinite metallic regions.
For simplicity we use air as the dielectric. We examine the impact of
nonlocality and the electron spill-out on the gap-plasmon, the unique mode
guided in $p$-polarization in the structure and propagating along the $z$
direction. More precisely, we compute the effective index $n_{\rm
eff}=\frac{k_{z}}{k_{0}}$ on the gap-plasmon according to the different models
described above and compare the results, in particular in the gap size range
($<$ 1 nm) where the charges densities from the two surfaces significantly
overlap, as schematically shown in Fig. 1. We begin with the simplified model
of spill-out for simple Drude metals, i.e., sodium (Na), which we compare with
the LRA and TFHT approximations as well as the QHT. Finally, we consider the
case of silver (Ag) in the QHT framework and compare to the commonly used LRA
and TFHT without spill-out.

Figure 1: Visual representation of overlapping of the equilibrium charge
densities $n_{0}(y)$ due to electron spill-out in the dielectric gap $g$
between two semi-infinite bulk metals. Deep inside the bulk metal, $n_{0}$ is
equal to the constant density of the bulk metal $n_{b}$. The magnetic mode
profile $H_{x}(y)$ of a $p-$polarized light propagating along the
$z-$direction is also shown. The structure is invariant along the $x-$ and
$z-$directions.

###  III.1 Simplified model

We begin with the simplified case presented above in which the equilibrium
electron density $n_{0}(\mathbf{r})$ can be described by an exponentially
decaying profile with the following analytical function,

$n_{0}(y)=n_{b}\left[\frac{1}{1+e^{\kappa\left(y+\frac{g}{2}\right)}}+\frac{1}{1+e^{-\kappa\left(y-\frac{g}{2}\right)}}\right]$
(21)

where the parameter $\kappa$ defines the spill-out or the decay of the charge
density from the metal surface. In what follows, we set
$\kappa\simeq(0.7a_{0})^{-1}$, with $a_{0}$ being the Bohr radius, which we
get by fitting of the analytical function with the decay of the self-
consistent charge density (see Eq. (9)), considering $\lambda_{\rm vW}=1/9$.
We take $r_{s}=4$, $\varepsilon_{\infty}=1$ and $\hbar\gamma=0.16$ eV, values
that correspond to Na. Although it is impossible to work with Na in practical
applications, due to its high reactivity, it presents an ideal Drude response
as the impact of interband transitions can be omitted, thus offering a
convenient theoretical study platform.

In Fig. 2 we compare the results of our simplified LRA with spill-out against
the commonly used LRA without any spill-out. The input equilibrium charge
density profiles used in the both approaches are plotted in Fig. 2(a) for
$g=1$ nm. In Fig. 2(b) and (c), we compare the behaviour of the real and
imaginary parts of the effective mode index. As expected, for large gaps
($g>1\leavevmode\nobreak\ $nm) spill-out can be neglected, the LRA with and
without electron spill-out agree very well. The spill-out has a noticeable
impact on the real part of the effective index for gaps that are typically
around 1 nm or less: it increases the effective index, signaling a plasmonic
slow down of the mode which is larger than without any form of spill-out. This
is understandable, as the Poynting vector inside a plasma is opposite to the
propagation direction, leading to a plasmonic drag [5], to a lower group
velocity and finally to a larger wave vector.

Figure 2: A comparison between the local response approximation (LRA) with and
without electron spill-out from the metal surfaces. (a) Equilibrium electron
density $n_{0}$ normalized with the bulk density $n_{\mathrm{b}}$ for
$g=1\leavevmode\nobreak\ $nm. The real (b) and (c) imaginary parts of the
effective refractive index $n_{\mathrm{eff}}$ plotted at
$E=1\leavevmode\nobreak\ $eV. The inset in (c) depicts a zoom-in of the region
where the spill-out has a significant impact on the mode index. The horizontal
dashed-grey line marks the refractive index of the bulk metal, $n_{\rm Bulk}$.
The dotted-red lines represent the region where it becomes hard to find a
proper solution due to numerical artifacts. (d)-(f) present a comparison
between the electric and magnetic field components considering
$g=1\leavevmode\nobreak\ $nm and $E=1\leavevmode\nobreak\ $nm. The plots in
the right column are normalized with the field value at the center of the gap,
except in (e) where field is zero at the center of the gap and the maximum
value of the field is considered for normalization.

When the gap is as small as 0.5 nm it begins to “close”, electromagnetically
speaking: the electron density inside the gap becomes so large that the
propagation of the gap-plasmon is hindered. The real part of the effective
index decreases quickly, while its imaginary part is much higher when the
spill-out is taken into account. It is important to note that in the limit of
vanishing gaps, where the spill-out is expected to virtually close the gap,
both real and imaginary parts of $n_{\rm eff}$ should ideally converge to the
index of the bulk metal. In Fig. 2(c), we notice that when
$g<0.35\leavevmode\nobreak\ $ nm, i.e., the gap distances where the gap is
expected to vanish completely, the imaginary parts of $n_{\rm eff}$ tends to
bend back. However, we have found to be very hard to obtain a proper solution
in this region. Indeed the dotted line as shown in Figs. 2(b) and (c)
represents a region where we could not find a solution due to numerical
artifacts. When $g<0.28\leavevmode\nobreak\ $nm, we find $n_{\rm eff}=n_{\rm
Bulk}$, where $n_{\rm Bulk}$ being the refractive index of the bulk metal.
These results are in agreement with the ones reported in [20].

We attribute the large imaginary part of the effective index to peaks in the
electric field that can be seen when the electron density reaches a value so
that $\chi_{f}\simeq 0$, meaning a vanishing local permittivity
$\varepsilon(y)$. The components of the electric field in the transverse and
longitudinal directions with respect to mode propagation are shown in Figs.
2(d) and (e), respectively. The extremely narrow peaks the electric field
exhibits have already been seen in a similar context [20] with a different
electron gas density. This means that they are fundamentally linked to the
assumption that a local permittivity is able to describe the spill-out
accurately. A clue that nonlocality may play an important role here is that
$\varepsilon=0$ is the condition for which volume plasmons are expected to be
supported inside the electron gas.

The peaks in the electric field drastically increase the losses experienced by
the gap plasmon, since the local absorption by the electron gas is
proportional to $\mathbf{E}_{1}\cdot\mathbf{J}^{*}$, where ∗ denotes the
complex conjugate [43]. This explains the larger imaginary part. We underline
that, on the contrary, the magnetic field remains smooth (see Fig. 2(f)).

###  III.2 QHT for simple metals

In the framework of QHT, we compute the self-consistent ground-state density
of a symmetric MIM waveguide from Eq. (9). First we consider that the guiding
structure is composed of an air gap sandwiched between two semi-infinite Na
metals, to better grasp the difference with the simplified spill-out model.

The self-consistent ground-state density $n_{0}(\mathbf{r})$ for a Na-air-Na
configuration is shown in Fig. 3 (a) considering different gaps. The metal-
dielectric interfaces for each gap size are indicated by the dashed lines. It
can be seen that when the inter-distance of the metal surfaces is sufficiently
large, there is no overlap of the charge densities. In this scenario, the
system can be well described using a classical formulation. However, when the
gap spacing between the metal boundaries decreases, the charge density
profiles start to overlap and the electrons from the two surfaces start to
tunnel across the gap. For sufficiently small gaps the charge transfer between
the metals is large enough to close the gap optically speaking, despite the
fact that the metal ion surfaces do not touch physically. Figure 3 (b)
presents a comparison of $n_{0}$ between the LRA, TFHT and QHT for gap size
$g=1$ nm. In the LRA and TFHT approximations, we assume a constant
$n_{0}=n_{b}$ in the metal and zero outside. In the QHT case, the space-
dependent $n_{0}(\mathbf{r})$ is plotted for both $\lambda_{\rm vW}=1/9$ and
$1$. The parameter $\lambda_{\rm vW}=1/9$ represents a lower spill-out and a
faster decay of the charge density, whereas $\lambda_{\rm vW}=1$ gives a
higher spill-out and a slower decay, as shown in Fig. 3 (c).

Figure 3: Equilibrium charge density $n_{0}$ calculated self-consistently for
an air gap sandwiched between the Na metals. (a) For different gap size ($g$).
The vertical dashed lines represent the metal boundaries for the corresponding
gaps. (b) A comparison between the LRA, TFHT and QHT for the spill-out
parameter $\lambda_{\rm vW}=1/9$ and $\lambda_{\rm vW}=1$, considering $g=1$
nm and (c) the decay of the charge density at $\lambda_{\rm vW}=1/9$ and $1$.
The results are normalized with the bulk density $n_{b}$ and the shaded grey
area represents the metal region.

In Fig. 4, we plot the effective refractive index as a function of gap
distance $g$. The results computed within the QHT, considering the self-
consistent space-dependent charge density, are compared against the LRA and
TFHT without spill-out. It is interesting to note that for extremely small gap
separation, LRA predicts very large values of both real and imaginary parts of
$n_{\rm eff}$, whereas the TFHT gives relatively smaller values. Nonetheless,
the trend is quite similar, i.e. $n_{\rm eff}$ increases monotonically as the
gap decreases. In contrast, in the QHT approximation, the $n_{\rm eff}$ first
increases with the decreasing gap size and then rapidly decrease after a
certain gap distance, getting a value of the refractive index of the bulk
metal. This is due to fact that at a certain gap size, a reasonable charge
transfer through the gap occurs between the metal surfaces due to the
overlapping of the charge densities and a virtual fade-out of the physical gap
appears when the spacing between the metal boundaries is further decreased.
For sufficiently small gaps, the mode virtually vanishes due to the electron
tunnelling and both the real and imaginary parts of the $n_{\rm eff}$
converges to the real and imaginary parts of the refractive index of the bulk
metal, as shown in the inset of Fig. 4(a) and (b), respectively. The expected
limit for the effective index of the gap-plasmon is in that case actually
retrieved.

Figure 4: Effective refractive index $n_{\rm eff}$ of Na-air-Na configuration
as a function of gap distance $g$. A comparison between the LRA, TFHT and QHT
is given. (a) Real and (b) imaginary parts of $n_{\rm eff}$, plotted at $E=1$
eV. The horizontal dashed-grey lines represent the refractive index of the
bulk metal.

In the QHT case, results for both $\lambda_{\rm vW}=1/9$ and $\lambda_{\rm
vW}=1$ are presented. Since $\lambda_{\rm vW}=1/9$ entails the lower spill-
out, the gap virtually vanishes when $g<0.3$ nm, however, the gap disappears
even at higher separation ($g<0.45$ nm) for $\lambda_{\rm vW}=1$, as it
implies higher electron spill-out.

Figure 5: (a) Normalized electric $|E_{y}|$ and (b) magnetic field $|H_{x}|$
field components plotted along the y-direction, i.e., in the normal direction
to the interfaces, plotted uder different approaches. The inset of (b) is a
close-up view of the $|H_{x}|$ at the interfaces. The plots are normalized by
the field value at the center of the gap. (c) and (d) show decay of the fields
at longer distances inside the metal. In these plots, we have considered $g=1$
nm, $E=1$ eV and $\lambda_{\rm vW}=1/9$ in the QHT case.

Now, to show the influence of the electron spill-out on the mode profile, we
plot here the electric $|E_{y}|$ and magnetic $|H_{x}|$ fields in the
orthogonal direction to the metal-air interfaces under different
approximations for $g=1$ nm and $E=1\leavevmode\nobreak\ $eV, as shown in Fig.
5. Both in the LRA and TFHT, $|E_{y}|$ is constant in the center of the gap
and then exponentially decays in the metal, as can be seen in Figs. 5 (a) and
(c). QHT, on the other hand, predicts that $|E_{y}|$ has a maximum value at
the center of the gap and it decreases as moving away from the center. The
effect of the electron pressure in the metal in the TFHT can also be observed
on the field profile. The $|H_{x}|$, as plotted in Fig. 5 (b), remains quite
similar in the gap in all approximations expect near the interfaces where QHT
predict a bit lower value. The decay of fields is given in Fig. 5 (c) and (d),
which show that the fields decay faster in QHT and slower in TFHT with respect
to the LRA fields, essentially because the wavevector along $z$ and thus the
decay along the $y$ axis is slightly different.

Finally, it is interesting to compare the QHT and the simplified LRA with
spill-out, as shown in Fig. 6. The real part of $n_{\rm{eff}}$ in both cases
accord well until the gap closes (Fig. 6 (a)) whereas, as described
previously, the imaginary part is much higher in the simplified spill-out
model (Fig. 6 (b)). This really means that nonlocality avoids the existence of
the peaks in the electric field when the electron density reaches a particular
value, making the imaginary part much closer to the prediction of the LRA
without spill-out.

Figure 6: (a) Real and (b) imaginary parts of $n_{\mathrm{eff}}$ within the
LRA with electron spill-out and full self-consistent QHT for $\lambda=1/9$.
The inset shows a zoom-in of the effective mode index in the tunneling regime.
The dotted-red lines represent the gap sizes where it becomes hard to find a
proper solution due to numerical artifacts, as was discussed in the main text
of Figs. 2(b) and (c).

###  III.3 QHT for noble metals

Noble metals, such as Ag, show in general a more complex behavior compared to
Drude-type metals due the contribution of interband transitions and core
electrons. In what follows, we use $r_{s}=3$ a.u. and $\hbar\gamma=0.03$ eV
for Ag, and account for interband effects with a constant background
dielectric constant $\varepsilon_{\infty}=5.9$ [40], which adds a local
contribution to the free-electron metal permittivity. Effective refractive
index for the Ag-air-Ag configuration as a function of the gap size computed
using different models are plotted in Fig. 7. The real and imaginary parts of
$n_{\rm eff}$ are plotted in Fig. 7 (a) and (b), respectively, showing quite
similar behavior as seen for the Na-Air-Na case in the previous section. LRA
and TFHT clearly show unphysical increasing character of effective index when
sub-nanometer gaps are in question. QHT, on the other hand, predicts that in
the tunneling regime, for sufficiently small metal separation, the gap
virtually fades away completely and the effective index converges to the
refractive index of the bulk Ag. Note that in Fig. 7, the QHT gives different
threshold, for the gap sizes for virtually absence of the any gap, depending
on the amount of electron spill-out considered from the metal surface, i.e.,
$\lambda_{\rm vW}=1/9$ or $\lambda_{\rm vW}=1$.

Figure 7: Effective refractive index for Ag-air-Ag configuration as a function
of gap distance. A comparison between the LRA, TFHT and QHT is given. Upper
and lower panels shows the real and imaginary parts of the effective index,
respectively. The results are plotted at $E=1$ eV. The dashed grey line
represents the refractive index of the bulk metal.

So far, we have analyzed the effective mode index of the MIM waveguide as
function of air gap between the metals. It is also important to examine the
behavior of $n_{\rm eff}$ as a function of energy. Figure 8 (a) and (b)
present the real and imaginary parts of $n_{\rm eff}$, respectively, as a
function of energy ($E$ in eV) for a Ag-air-Ag waveguide, considering a gap
size $g=1$ nm – a width which is experimentally rather easy to obtain [14].

This underlines the differences between the predictions of the different
approaches regarding the effective index of the mode. The LRA predicts a bend
back for the real part and a much larger imaginary part in this region. The
TFHT on the contrary predicts a much lower imaginary part at any frequency.
The QHT approach predicts a more reasonable behaviour of the gap-plasmon with
no bend back in the effective index and values of its imaginary part that can
be much lower than for the LRA at higher frequencies.

The lower the imaginary part of the effective index, the larger the typical
propagation length of the gap-plasmon, and finally the larger the quality
factor of an eventual cavity for the gap-plasmon [44]. This underlines the
importance of an accurate prediction at in this regime.

Figure 8: Effective refractive index as a function of energy (in eV) at
$g=1\leavevmode\nobreak\ $nm. A comparison between the LRA and QHT is given.
Upper and lower panels shows the real and imaginary parts of the effective
index. The results are normalized with the screened plasma frequency,
$\omega_{p}^{\prime}=\omega_{p}/\sqrt{\varepsilon_{\infty}}$.

##  IV Effective gap

To summarize qualitatively the above results, the spill-out seems to be able
to close the gap electromagnetically speaking before it is actually spatially
or mechanically closed. In addition, we have seen that when the spill-out is
taken into account, whether in the simplified model or within the QHT
framework, the effective index of the gap-plasmon is slightly larger than when
the spill-out is neglected, however it presents a very similar trend. All
these results suggest that everything occurs as if the gap were narrower when
the spill-out is taken into account, which seems physically sound. The fact
that, for the simplified model, peaks seem to appear when the electron density
reaches a certain value suggests that it should be possible to define the
electromagnetic edges of the metal, as well as an effective value for the gap
width. We define this value as the width that should be considered in the LRA
model or in the TFHT approach for the gap-plasmon to fit the effective index
of the gap-plasmon provided by the QHT.

The effective gap is thus given by

$g_{e}=g-\delta_{g}$

where $g$ is the actual gap between the edges of the metal. Said otherwise, in
order to retrieve the results of the QHT, one can use the LRA (or the TFHT)
with hard wall boundaries, but simply use a diminished size of the gap to take
the spill-out into account.

For silver, shifting the LRA curve of $\delta_{g}=0.15$ nm or the TFHT curve
of $\delta_{g}=0.22$ nm allows to retrieve the effective index of the gap-
plasmon as computed with the QHT extremely accurately, as long as the gap does
not start closing (which occurs around $g\simeq 0.6$ nm). The results for the
QHT reported in Fig. 9 are for $\lambda=1/9$, however, the QHT result for
$\lambda=1$ can also be nicely reproduced from the LRA and TF approximations
considering $\delta_{g}\simeq 0.25$ nm and $\delta_{g}\simeq 0.32$ nm,
respectively. The idea that the gap-plasmon is sensitive to the effective gap
between the metals is thus well founded, as it yields quantitative results.
This can be understood because the spill out makes the response of the
electron gas outside the metal similar to the response inside the metal, as
soon as the electron density is high enough. However, we point out that we
have not found another way to define the effective gap based on the profile of
the electron density or of the fields. The values we give are simply the ones
that yield the best fit and there does not seem to be any simple way to
predict such values.

This discussion finally underlines that at such scales, for the gap-plasmon,
the tiny difference (a rule of thumb would be, according to our results, to
consider a 0.1 nm) which exists between the actual edge of the metal and the
”electromagnetic edge” has its importance and must be considered in the limit
of vanishing gaps.

Figure 9: The effective mode index for Ag-air-Ag waveguide (a) real and (b)
imaginary part as a function of gap distance $g$ between the metal edges
plotted at $E=1\leavevmode\nobreak\ $eV. The LRA and TFHT curves are shifted
by $\delta_{g}$ to cope with the QHT results, which allows to introduce the
concept of effective gap that can be defined as $g_{e}=g-\delta_{g}$. The
$g_{e}$ seems to appear smaller than the mechanical gap $g$ by an amount of
$\delta_{g}$. In particular, using $\delta_{g}=0.15$ nm in the LRA and
$\delta_{g}=0.22$ nm in the TF approximation gives good fit with QHT. The
insets plot $n_{\rm{eff}}$ as a function of actual gap.

##  V Conclusions

In this article, we analyzed the impact of both electron spill-out and
nonlocality on the propagation of gap-plasmons in MIM waveguides using the
full self-consistent QHT, which can efficiently take into account the spatial
variation of the electron density outside of the metal. We also compare our
results with the commonly used approaches (LRA and TFHT) without any spill-out
and with a LRA model for the spill out.

Peaks in the electric field profile and a large imaginary part for the
effective index of the gap plasmon are predicted by the latter approach. Such
features are however absent in the framework of the QHT, showing that
nonlocality is key to an accurate description of the gap plasmon in the
tunneling regime. We underline that, in this framework, the effective mode
index converges to the refractive index of the bulk metal, as expected from
physical considerations. These points reinforce the idea that QHT is the right
tool for plasmonics at the subnanometer level.

Our study thus shows that the spill-out indeed imposes an upper limit on the
effective index of the gap-plasmon, confirming previous studies [20]. A gap-
plasmon is very unlikely to present an effective index larger than 15 for
noble metals in the visible.

Importantly, whatever model is chosen to describe the spill-out, it seems
possible to introduce an effective gap width between metals, smaller than the
gap size defined as the distance between the atoms of the metal across the
gap. This effective gap width allows to take simply the effect of the spill-
out into account, as long as the gap can be considered open to the propagation
of the mode. In many experiments, the definition of the gap between metals is
of crucial importance [14]. This reasoning suggests using optical
characterization methods (like ellipsometry typically) to determine the
thickness of a material deposited at the surface of a metal is particularly
relevant in the context of gap-plasmon resonators, as the optical thickness
which is retrieved probably takes partially the spill-out into account. The
gap width deduced from such measurements is in that way probably closer to the
effective gap width we have defined here. This may thus explain why hard-wall
approaches have been so far surprisingly successful at predicting the
experimental behaviour of gap-plasmons in the limit of vanishing gaps.

###### Acknowledgements.

Antoine Moreau is an Academy CAP 20-25 chair holder. He acknowledges the
support received from the Agence Nationale de la Recherche of the French
government through the program “Investissements d’Avenir” (16-IDEX-0001 CAP
20-25).

## References

  * Bozhevolnyi [2006] S. I. Bozhevolnyi, Effective-index modeling of channel plasmon polaritons, Opt. Express 14, 9467 (2006).
  * Bozhevolnyi and Søndergaard [2007] S. I. Bozhevolnyi and T. Søndergaard, General properties of slow-plasmon resonant nanostructures: nano-antennas and resonators, Optics express 15, 10869 (2007).
  * Bozhevolnyi and Jung [2008] S. I. Bozhevolnyi and J. Jung, Scaling for gap plasmon based waveguides, Opt. Express 16, 2676 (2008).
  * Oulton  _et al._ [2008] R. F. Oulton, V. J. Sorger, D. A. Genov, D. F. P. Pile, and X. Zhang, A hybrid plasmonic waveguide for subwavelength confinement and long-range propagation, Nature Photonics 2, 496 (2008).
  * Ajib  _et al._ [2019] R. Ajib, A. Pitelet, R. Pollès, E. Centeno, Z. Ajaltouni, and A. Moreau, The energy point of view in plasmonics, JOSA B 36, 1150 (2019).
  * Wang and Kempa [2007] X. Wang and K. Kempa, Plasmon polaritons in slot waveguides: Simple model calculations and a full nonlocal quantum mechanical treatment, Phys. Rev. B 75, 245426 (2007).
  * Jung  _et al._ [2009] J. Jung, T. Søndergaard, and S. I. Bozhevolnyi, Gap plasmon-polariton nanoresonators: Scattering enhancement and launching of surface plasmon polaritons, Phys. Rev. B 79, 035401 (2009).
  * Yang  _et al._ [2012] J. Yang, C. Sauvan, A. Jouanin, S. Collin, J.-L. Pelouard, and P. Lalanne, Ultrasmall metal-insulator-metal nanoresonators: impact of slow-wave effects on the quality factor, Optics Express 20, 16880 (2012).
  * Moreau  _et al._ [2012] A. Moreau, C. Ciracì, J. J. Mock, R. T. Hill, Q. Wang, B. J. Wiley, A. Chilkoti, and D. R. Smith, Controlled-reflectance surfaces with film-coupled colloidal nanoantennas, Nature 492, 86 (2012).
  * Baumberg  _et al._ [2019] J. J. Baumberg, J. Aizpurua, M. H. Mikkelsen, and D. R. Smith, Extreme nanophotonics from ultrathin metallic gaps, Nature materials 18, 668 (2019).
  * Akselrod  _et al._ [2014] G. M. Akselrod, C. Argyropoulos, T. B. Hoang, C. Ciracì, C. Fang, J. Huang, D. R. Smith, and M. H. Mikkelsen, Probing the mechanisms of large purcell enhancement in plasmonic nanoantennas, Nature Photonics 8, 835 (2014).
  * Hoang  _et al._ [2015] T. B. Hoang, G. M. Akselrod, C. Argyropoulos, J. Huang, D. R. Smith, and M. H. Mikkelsen, Ultrafast spontaneous emission source using plasmonic nanoantennas, Nature communications 6, 1 (2015).
  * Søndergaard  _et al._ [2012] T. Søndergaard, S. M. Novikov, T. Holmgaard, R. L. Eriksen, J. Beermann, Z. Han, K. Pedersen, and S. I. Bozhevolnyi, Plasmonic black gold by adiabatic nanofocusing and absorption of light in ultra-sharp convex grooves, Nature communications 3, 1 (2012).
  * Ciraci  _et al._ [2012] C. Ciraci, R. T. Hill, J. J. Mock, Y. Urzhumov, A. I. Fernandez-Dominguez, S. A. Maier, J. B. Pendry, A. Chilkoti, and D. R. Smith, Probing the Ultimate Limits of Plasmonic Enhancement, Science 337, 1072 (2012).
  * Moreau  _et al._ [2013] A. Moreau, C. Ciracì, and D. R. Smith, Impact of nonlocal response on metallodielectric multilayers and optical patch antennas, Phys. Rev. B 87, 045401 (2013).
  * Raza  _et al._ [2013] S. Raza, T. Christensen, M. Wubs, S. I. Bozhevolnyi, and N. A. Mortensen, Nonlocal response in thin-film waveguides: Loss versus nonlocality and breaking of complementarity, Phys. Rev. B 88, 115401 (2013).
  * 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).
  * Ciracì  _et al._ [2013] C. Ciracì, J. B. Pendry, and D. R. Smith, Hydrodynamic Model for Plasmonics: A Macroscopic Approach to a Microscopic Problem, ChemPhysChem 14, 1109 (2013).
  * Esteban  _et al._ [2012] R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, Bridging quantum and classical plasmonics with a quantum-corrected model, Nat Commun 3, 825 (2012).
  * Skjølstrup  _et al._ [2018] E. J. H. Skjølstrup, T. Søndergaard, and T. G. Pedersen, Quantum spill-out in few-nanometer metal gaps: Effect on gap plasmons and reflectance from ultrasharp groove arrays, Phys. Rev. B 97, 115429 (2018).
  * Skjølstrup  _et al._ [2019] E. J. H. Skjølstrup, T. Søndergaard, and T. G. Pedersen, Quantum spill-out in nanometer-thin gold slabs: Effect on the plasmon mode index and the plasmonic absorption, Phys. Rev. B 99, 155427 (2019).
  * Teperik  _et al._ [2013] 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 (2013).
  * 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, The Journal of Physical Chemistry C 117, 8941 (2013), https://doi.org/10.1021/jp401887y .
  * Zhu  _et al._ [2016] W. Zhu, R. Esteban, A. G. Borisov, J. J. Baumberg, P. Nordlander, H. J. Lezec, J. Aizpurua, and K. B. Crozier, Quantum mechanical effects in plasmonic structures with subnanometre gaps, Nature Communications 7, 11495 (2016).
  * Zuloaga  _et al._ [2009] J. Zuloaga, E. Prodan, and P. Nordlander, Quantum description of the plasmon resonances of a nanoparticle dimer, Nano Letters 9, 887 (2009), pMID: 19159319, https://doi.org/10.1021/nl803811g .
  * Ullrich [2011] C. A. Ullrich,  _Time-Dependent Density-Functional Theory: Concepts and Applications_ (Oxford University Press, Oxford, 2011).
  * Aguirregabiria  _et al._ [2018] G. Aguirregabiria, D. C. Marinica, R. Esteban, A. K. Kazansky, J. Aizpurua, and A. G. Borisov, Role of electron tunneling in the nonlinear response of plasmonic nanogaps, Phys. Rev. B 97, 115430 (2018).
  * Toscano  _et al._ [2015] G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl, F. Evers, H. Xu, N. Asger Mortensen, and M. Wubs, Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics, Nat Commun 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).
  * Ciracì [2017] C. Ciracì, Current-dependent potential for nonlocal absorption in quantum hydrodynamic theory, Phys. Rev. B 95, 245434 (2017).
  * Ding and Chan [2017] K. Ding and C. T. Chan, Plasmonic modes of polygonal rods calculated using a quantum hydrodynamics method, Phys. Rev. B 96, 125134 (2017).
  * Manfredi and Haas [2001] G. Manfredi and F. Haas, Self-consistent fluid model for a quantum electron gas, Physical Review B 64, 075316 (2001).
  * Crouseilles  _et al._ [2008] N. Crouseilles, P.-A. Hervieux, and G. Manfredi, Quantum hydrodynamic model for the nonlinear electron dynamics in thin metal films, Physical Review B 78, 155412 (2008).
  * Smith  _et al._ [2015] C. L. C. Smith, N. Stenger, A. Kristensen, N. A. Mortensen, and S. I. Bozhevolnyi, Gap and channeled plasmons in tapered grooves: a review, Nanoscale 7, 9355 (2015).
  * Khalid  _et al._ [2018] M. Khalid, F. Della Sala, and C. Ciracì, Optical properties of plasmonic core-shell nanomatryoshkas: A quantum hydrodynamic analysis, Opt. Express 26, 17322 (2018).
  * Khalid and Ciracì [2019] M. Khalid and C. Ciracì, Numerical Analysis of Nonlocal Optical Response of Metallic Nanoshells, Photonics 6, 39 (2019).
  * Ciracì  _et al._ [2019] C. Ciracì, R. Jurga, M. Khalid, and F. Della Sala, Plasmonic quantum effects on single-emitter strong coupling, Nanophotonics 8, 1821 (01 Oct. 2019).
  * Yang and Ding [2021] F. Yang and K. Ding, Electron spill-out effect in singular metasurfaces, Photonics 8, 10.3390/photonics8050154 (2021).
  * Khalid and Ciracì [2020] M. Khalid and C. Ciracì, Enhancing second-harmonic generation with electron spill-out at metallic surfaces, Communications Physics 3, 1 (2020).
  * Baghramyan  _et al._ [2021] H. M. Baghramyan, F. Della Sala, and C. Ciracì, Laplacian-level quantum hydrodynamic theory for plasmonics, Phys. Rev. X 11, 011049 (2021).
  * Brack [1993] M. Brack, The physics of simple metal clusters: Self-consistent jellium model and semiclassical approaches, Rev. Mod. Phys. 65, 677 (1993).
  * Sakat  _et al._ [2021] E. Sakat, A. Moreau, and J.-P. Hugonin, Generalized electromagnetic theorems for nonlocal plasmonics, Phys. Rev. B 103, 235422 (2021).
  * Lemaître  _et al._ [2017] C. Lemaître, E. Centeno, and A. Moreau, Interferometric control of the absorption in optical patch antennas, Scientific Reports 7, 2941 (2017).


# Evolving Self-Assembling Neural Networks: From Spontaneous Activity to
Experience-Dependent Learning

Erwan Plantec, Joachim W. Pedersen, Milton L. Montero, Eleni Nisioti
Sebastian Risi

IT University of Copenhagen, Denmark
{erpl<EMAIL_ADDRESS>
###### Abstract

Biological neural networks are characterized by their high degree of
plasticity, a core property that enables the remarkable adaptability of
natural organisms. Importantly, this ability affects both the synaptic
strength and the topology of the nervous systems. Artificial neural networks,
on the other hand, have been mainly designed as static, fully connected
structures that can be notoriously brittle in the face of changing
environments and novel inputs. Building on previous works on Neural
Developmental Programs (NDPs), we propose a class of self-organizing neural
networks capable of synaptic and structural plasticity in an activity and
reward-dependent manner which we call _Lifelong Neural Developmental Program_
(LNDP). We present an instance of such a network built on the graph
transformer architecture and propose a mechanism for pre-experience plasticity
based on the spontaneous activity of sensory neurons. We demonstrate the
model’s ability to learn from experiences in different control tasks starting
from randomly connected or empty networks. We further show that structural
plasticity is advantageous in environments necessitating fast adaptation or
with non-stationary rewards.

## Introduction

Millions of years of open-ended evolutionary processes have led to the
emergence of complex cognitive abilities, which are primarily supported by
brains composed of continually adapting neurons and synapses. These highly
plastic biological neural networks enable organisms to quickly adapt in
ontogenetic timescales through learning, a major evolutionary transition in
the evolution of cognition (Ginsburg and Jablonka,, 2021). Furthermore, by
growing through a temporally extended developmental phase from a DNA
representation into an extremely complex system, brains also exhibit a high
degree of adaptability on a phylogenetic timescale.

On the other hand, artificial neural networks (ANNs) have been mainly built as
fixed fully connected structures that remain frozen after training, unable to
adapt to any unexpected change. Indeed, even though much research has been
devoted to the topic of online and meta-learning Hochreiter et al., (2001);
Andrychowicz et al., (2016); Wang et al., (2018), it is still the case that
the most modern neural network systems use offline learning as it is much
simpler to use in combination with backpropagation.

Fundamentally, biological and artificial neural networks seem to follow
diametrically opposed design philosophies: self-organization and engineering
(Ha and Tang,, 2022). In the past decade, self-organization has gained a lot
of interest in the AI community as a path towards closer to natural artificial
intelligence (Risi,, 2021). Several works have attempted to design ANNs
endowed with more of the characteristics found in living organisms in order to
replicate properties such as fast adaptability and lifelong learning, or
generalization capabilities (Kudithipudi et al.,, 2022; Parisi et al.,, 2019).
For instance, Evolved Plastic Artificial Neural Networks (EPANNs) have been
proposed as a research framework aiming at building such neural networks
through neuroevolution (Soltoggio et al.,, 2018), while work on indirect
encodings of artificial neural networks has proposed several models of self-
assembling neural networks (Najarro et al.,, 2022) (see Stanley and
Miikkulainen, (2003) for a taxonomy of developmental encodings for artificial
neural networks).

Recently, Najarro et al., (2023) proposed a model capable of assembling
functional artificial neural networks represented as graphs, which they call
Neural Developmental Programs (NDPs). However, these models are temporally
limited to a pre-environmental phase and do not account for lifetime, and even
less lifelong learning.

In this work, we propose a way to address those limitations by extending the
NDP framework. Specifically, we incorporate a mechanism that enables
plasticity and structural changes throughout the lifetime of an agent. It
achieves this by performing local computations which depend on both the local
activity of each neuron in the ANN and the global reward function from the
environment. The resulting system thus defines a family of plastic neural
networks bridging the gap between indirect developmental encodings and meta-
learned plasticity rules which we call _Lifelong Neural Developmental
programs_ (LNDP).

LNDPs are defined by a set of possibly parametric components defining neural
and synaptic dynamics. Motivated by works on architectural priors (Bhattasali
et al.,, 2022) and the role of structural plasticity in learning (Caroni et
al.,, 2012; May,, 2011), LNDPs are also endowed with structural plasticity
(i.e. synapses can be dynamically added or removed). We further propose an
instance of an LNDP based on a Graph Transformer layer (Dwivedi and Bresson,,
2021), modeling neuronal communication through synaptic and extra-synaptic
channels, which enables neurons to self-organize and differentiate. Neural and
synaptic dynamics are then modeled with Gated Recurrent Units (GRUs) (Chung et
al.,, 2014). These components together define a large family of learnable
LNDPs which we can optimize to solve different reinforcement learning tasks.

Motivated by works exploring the role of spontaneous activity (SA) in
developing biological (Luhmann et al.,, 2016; Leighton and Lohmann,, 2016;
Kirkby et al.,, 2013) and artificial (Raghavan and Thomson,, 2019; Bednar and
Miikkulainen,, 1998; Valsalam et al.,, 2007) neural networks, we further
extend the system with a mechanism that enables pre-experience development
based on SA. Here we model SA with a simple learnable stochastic process of
sensory neurons. Using SA enables the reuse of the same components during both
pre-experience development and experience-dependent learning.

In this work, we demonstrate the effectiveness of the LNDP in driving
artificial neural networks to self-organize in an activity and experience-
dependent manner into functional networks able to solve control tasks,
starting from randomly connected or empty networks. We further present
evidence that structural plasticity improves results in environments requiring
fast adaptation or with non-stationary dynamics requiring continuous
adaptation. Additionally, the effectiveness of pre-environmental spontaneous
activity-driven developmental phases to drive the network’s self-organization
into functional units is demonstrated.

## Background

This section reviews relevant related works, such as meta-learning synaptic
plasticity, developmental encodings, and approaches that investigate the
intersection between structural and synaptic plasticity.

### Meta-learning Synaptic Plasticity

Multiple works have tried to meta-learn synaptic plasticity rules to discover
new learning rules for ANNs potentially enabling faster and more sample
efficient adaptation of artificial agents. For instance, several works have
meta-trained gradientless synaptic plasticity rules enabling faster and more
sample-efficient learning than traditional stochastic gradient descent
(Sandler et al.,, 2021; Kirsch and Schmidhuber,, 2022). Recently, Pedersen et
al., have proposed a model called Structurally Flexible Neural Networks
(SFNNs), which are able to learn from rewards in different environments with
different input and output dimensions (Pedersen et al.,, 2024). More
biologically inspired self-organized networks, based on Spike Timing-Dependent
Plasticity-like rules have also been developed for reservoir computing (Lazar
et al.,, 2009). Aswolinskiy and Pipa, also proposed a reward-modulated version
of the same networks for reinforcement learning (Aswolinskiy and Pipa,, 2015).
Najarro and Risi, (2020) evolved Hebbian synaptic plasticity rules allowing
agents to perform well in locomotion tasks from initially random networks and
even adapt to damages unseen during training.

### Developmental Encodings

Developmental encodings are a subset of indirect encodings characterized by a
process unrolling in time while indirect encodings can be direct “one-shot”
mappings (like hypernetworks (Ha et al.,, 2022)). For instance, L-systems are
able to grow artificial neural networks through a process of
development/rewriting (de Campos et al.,, 2011). In Lowell and Pollack,
(2016), the authors used a model of gene regulatory networks for this purpose
and showed it led to the emergence of modular networks. Neural Cellular
Automata (NCA) have also been demonstrated to allow the control of neural
network development (Najarro et al.,, 2022). More recently, Najarro et al.,
(2023) proposed NDPs which are instances of developmental encodings based on
graph neural networks.

## Lifelong Neural Developmental Programs

On a high level, the Neural Developmental Program (NDP) introduced in Najarro
et al., (2023), defines a set of neural components that can be optimized to
direct the growth of an artificial neural network into functional networks in
a fully decentralized manner (i.e. neurons perceive the state of the neural
network through a graph neural network layer and can independently of each
other decide to divide). However, the approach was so far restricted to
growing neural structures that were fixed during the lifetime of the agent.
This section presents the Lifelong Learning NDP (LNDP) approach, which enables
learning to modify the weights and the connectivity of the neural network
while the agent interacts with the environment. This ability is especially
important for environments with non-stationary rewards that require lifetime
learning. A similar objective has been addressed in (Ferigo and Iacca,, 2023).
However, in their approach, the number of trainable parameters of their model
is tied to the size of the produced networks while this is not the case for
NDPs.

In more detail, we represent a neural network at a given time $t$ as a
directed graph $\mathcal{G}^{t}$ with $N$ nodes (neurons) which is fully
described by the tuple $<A^{t},I,O,h^{t},e^{t},v^{t},w^{t}>$. $A^{t}$ is the
$N\times N$ binary adjacency matrix indicating the presence of an edge in
between two nodes. $I$ and $O$ are the sets of input and output nodes
respectively. $h^{t}\in H^{N}$ and $e^{t}\in E^{N^{2}}$ are the nodes and
edges states respectively with $H\equiv\mathbb{R}^{dh}$ and
$E\equiv\mathbb{R}^{de}$. Edge states are masked by the adjacency matrix, i.e.
set to zero wherever the adjacency matrix is 0. $v^{t}\in\mathbb{R}^{N}$ are
the nodes’ activations and $w^{t}$ is the real-valued weight matrix. We refer
as $h_{i}$ to the state of node $i$ and $e_{ij}$ to the state of the edge
between nodes $i$ and $j$.

In the following, we describe the components that together make a LNDP. These
components can be parameterised learnable functions. We note the full set of
parameters as $\theta$.

#### Initialization

We note $\mathcal{G}^{0}\equiv<A^{0},I,O,h^{0},e^{0},v^{0},w^{0}>$ as the
initial network, which is sampled from a distribution
$\Psi=\mathbb{P}(\mathcal{G}^{0}=g)$. Input and output neurons stay fixed
during the agent’s lifetime.

#### Nodes

$f_{\theta}^{h}:\mathcal{G}\to H^{N}$ the node function, which updates the
nodes’ states based on the complete graph state, is defined as:

$h^{t+1}=f_{\theta}^{h}(\mathcal{G}^{t}).$ (1)

Using the complete graph state to update the nodes enables taking information
from the network structure through interactions between nodes (e.g with graph
neural networks). While many different types of information can be used in the
node update, taking activity states $v^{t}$ into account allows for activity-
dependent mechanisms. Intuitively, one can see the nodes’ states as the
molecular or membrane states of neurons which can, in turn, affect plasticity
and neuronal dynamics.

#### Edges

Edges’ states are updated through a function $f_{\theta}^{e}:H\times
H\times\mathbb{R}\to E$, which updates their states based on the states of
connected nodes and a reward signal $r^{t}$ received from the environment. We
refer to this function as the edge or synapse model:

$e_{ij}^{t+1}=f_{\theta}^{e}(e_{i}j^{t},h^{t+1}_{i},h^{t+1}_{j},r^{t}).$ (2)

Following, a weight model function $f_{\theta}^{w}:E\to\mathbb{R}$ associates
a scalar weight to each edge based on its state. Note that if nodes are
updated using activities, then some edge update functions could be able to
learn activity-dependent rules such as Hebbian learning:

$w_{ij}^{t}=f^{w}_{\theta}(e^{t}_{ij}).$ (3)

#### Structural plasticity

Structural plasticity in the network is achieved by allowing new edges to be
added (i.e. synaptogenesis) and existing edges to be removed (i.e. pruning).
$f_{\theta}^{+}:H^{2}\to[0,1]$ is the synaptogenesis function which defines
the probability of adding an edge (synapse) between two nodes at each time
step:

$\mathbb{P}(A_{ij}\leftarrow 1)=f_{\theta}^{+}(h_{i}^{t},h_{j}^{t}).$ (4)

Similarly, $f_{\theta}^{-}:E\to[0,1]$ is the pruning function and defines the
probability of removing an edge in the graph (i.e. setting its entry in the
adjacency matrix to $0$) as a function of its state:

$\mathbb{P}(A_{ij}\leftarrow 0)=f_{\theta}^{-}(e_{ij}^{t}).$ (5)

#### Network dynamics

We define the network dynamics through a dynamical model $\phi$ updating the
nodes’ activations based on past activations, the current weight matrix
$w^{t}$, nodes’ states $h^{t}$, and observation $o^{t}$:

$v^{t+1}=\phi(v^{t},w^{t},h^{t},o^{t})$ (6)

Node states can be used to define neuron parameters such as biases.

### Implementation

(a) Node model

(b) Edge model

(c) Network topology

Figure 1: LNDP components. (a) Node model Node features $h^{t}$, activations
$v^{t}$ as well as additional graph structural features are fed through a
single Graph Transformer layer whose output is fed as input to a GRU to obtain
new nodes states $h^{t+1}$. (b) Edge model Edges are also modelled with GRUs
and take as input pre and post-synaptics nodes’ states and the last reward
received $r^{t}$. (c) Network topology Networks are divided in input (blue),
hidden (black) and output (red) neurons. Connections can only exist from input
to hidden, hidden to hidden, and hidden to output. Some nodes can have no
connections at all and the total number of hidden nodes is constant.
Hyperparameters $\mu^{conn}$ and $\sigma^{conn}$ define the distribution
(truncated normal) of the initial network density.

We describe here an instance of LNDPs based on Graph Transformer layers (GT)
(Dwivedi and Bresson,, 2021) for the node update. More precisely, nodes are
updated using Gated Recurrent Units (GRUs) (Chung et al.,, 2014) and take as
input the output of a GT layer as depicted in Figure 1a, which enables nodes
to “perceive” the network and self-organize. The GT layer takes as input the
concatenation of nodes activation $v^{t}$, nodes’ states $h^{t}$, and
structural graph features associated with the nodes. Structural features are
the in-degree, out-degree, and total degree, as well as a one-hot encoding
indicating if the node is an input, hidden, or output node. Moreover, we
introduce edge features in the attention layer as proposed in Dwivedi and
Bresson, (2021). We also augment edge features with structural features which
are 2 bits indicating if there is a forward or backward connection between
nodes and a bit indicating if the edge is a self-loop. Intuitively, the graph
transformer layer is used to model node-node interactions allowing them to
self-organize and differentiate. We chose the GT rather than another graph
neural network architecture as it seems easier to train. Importantly, it
allowed to keep nodes’ states from collapsing (i.e. all becoming the same), a
recurrent problem with NDPs (Pedersen et al.,, 2024). Node embeddings
outputted by the GT are then fed into the GRU unit to update nodes’ states
($h_{i}^{t+1}=GRU(h^{t}_{i},x_{i}^{t})$ where $x_{i}^{t}$ is the output of the
GT layer). Importantly, all nodes share the same GRU parameters. Synapses are
also modeled as GRUs with shared parameters taking as input the concatenation
of the pre and post-synaptic neurons’ states as well as their activity and the
reward received at the last time step $r$ as depicted in Figure 1b.

Pruning and synaptogenesis functions ($f_{\theta}^{-}$ and $f_{\theta}^{+}$
respectively) are modeled with multi-layer perceptrons (MLPs) with 1 hidden
layer and ReLU activation function followed by a sigmoid to output the
probability of respective actions. In this work, the dynamical model is
defined as follows:

$\phi(v^{t},w^{t},o^{t})=tanh(\hat{v}^{t}\cdot w^{t}),$ (7)

where $\hat{v}^{t}$ is the the activation vector with input nodes’ activations
clamped to the values of the observation $o^{t}$. The weight matrix $w^{t}$ is
obtained by taking the first element of the edges states
($w^{t}_{ij}=e^{t}_{ij,0}$).

The initial network topology (i.e. adjacency matrix, see Figure 1c) is
randomly initialized according to connection parameters $\mu_{conn}$ and
standard deviation $\sigma_{conn}$, which together define a truncated normal
distribution from which the connection probability is sampled:
$\mathbb{P}(A^{0}_{ij}=1)\sim\mathcal{N}^{[0,1]}(\mu^{conn},\sigma^{conn})$.
The connection probability defines the probability that an edge exists from
input to hidden, hidden to hidden, and hidden to output neurons at time $t=0$.
There is a fixed number of nodes $N$ but some nodes can have no connections at
all. Nodes’ activations are initialized as zeros. Nodes’ and edges’ initial
features ($h^{0}$ and $e^{0}$ respectively) are sampled from a uniform
distribution between $-1$ and $1$.

While the model can differentiate between input, hidden, and output neurons as
this information is given to the GT layer, it cannot know a priori which of
its input (output) nodes corresponds to which information (action). This is
because neurons are initialized following the same distribution and have the
same connection probabilities and thus the same expected connectivity degree.
Thus, for the model to solve a task, it has to differentiate its inputs and
output neurons in order to make the right connections. This can be achieved
for input neurons by looking at specific activation patterns, either from the
environment or from spontaneous activity (see next paragraph). For output
neurons, this can be achieved by finding specific dynamical relations between
output neurons, rewards and/or input neuron activations to infer the effect of
each of the output neurons.

### Spontaneous activity

We also study a mechanism for pre-experience (i.e. before receiving feedback
from the environment) development of networks based on spontaneous activity,
which has been shown to help the self-organization and learning of
convolutional neural networks (Raghavan and Thomson,, 2019). In our case, we
implement spontaneous activity through a learnable generative model generating
spontaneous activity of the input nodes replacing observations from the
environment, which we denote as $o_{SA}^{t}$. While many different models
could be used, we chose to train a Ornstein-Uhlenbeck stochastic process
defined as:

$\begin{cases}o_{SA}^{t+1}=o_{SA}^{t}+\alpha(\mu-o_{SA}^{t})+W^{t},\\\
W^{t}\sim\mathcal{N}(0,\Sigma).\end{cases}$ (8)

Mean $\mu\in\mathbb{R}^{do}$, covariance $\Sigma\in\mathbb{R}^{do\times do}$,
and time constants $\alpha\in\mathbb{R}^{do}$ are learnable parameters of the
model where $do$ is the observation space dimensionality. This phase takes
place before the environmental phase and lasts for $T^{SA}$ steps. Using
spontaneous activity for pre-experience development enables the use of the
same synaptic plasticity mechanisms as in the environmental phase.

Figure 2: Training curves of the LNDP with varying initialization
distributions (where $\mu^{conn}$ is the mean connection probability and
$\sigma^{conn}$ its variance) and structural plasticity (SP) enabled (red) and
disabled (blue). For all conditions, structurally plastic LNDPs outperform non
structurally plastic ones in Cartpole and Foraging. Models without PS are not
evaluated on empty networks ($\mu^{conn}=0$) as they would obviously fail (the
network would remain empty).

## Methods

We optimized LNDPs on a set of reinforcement learning tasks using the
Covariance Matrix Adaptation Evolutionary strategy (CMA-ES). The following
sections give more details about the environments in use and the training
procedure. LNDP hyperparameters used in this work are listed in Table 1. The
code used to produce the results is available at:
https://github.com/erwanplantec/LNDP.

$N$ | # nodes | $\in\\{32,64\\}$
---|---|---
$dh$ | node features | $8$
$de$ | edge features | $4$
$T^{SA}$ | developmental steps | $\in\\{0,100\\}$
$\mu^{conn}$ | initial density average | $\in\\{0,0.5,1.\\}$
$\sigma^{conn}$ | initial density variance | $\in\\{0,0.1\\}$
(a) LNDP hyperparameters

Population size | $128$
---|---
Initial CMA-ES variance | $0.065$
MC trials | $3$
(b) CMA-ES hyperparameters

Table 1: Model and training hyperparameters.

### Environments

We use three classical control tasks (Cartpole, Acrobot and Pendulum) and one
foraging task with non-stationary dynamics requiring lifetime adaptivity from
the agent. For all tasks, the fitness of an agent is defined as the average
return of the agent over three different trials (i.e. different random seeds).

  * •

Cartpole: discrete control task in which the agent controls a cart with a pole
on top. The agent has to balance the pole to keep it straight. This task
features a 4-dimensional input space and a 2-dimensional discrete action
space. The action is determined by the highest activated output nodes. The
optimal reward is $500$. Importantly, it is difficult for an agent to find a
solution directly solving the environment when starting from a random network
(as is the case with LNDPs) as the agent needs to go from a random to a
functional network in a short amount of time. Thus, we repeat the episodes
without resetting the network states to give time for the model to adapt.

  * •

Acrobot: discrete control task where the agent has to swing up a double
pendulum above a certain height by applying torque to the joint, with the
pendulum initially hanging downward. This environment features a 6-dimensional
input space and a 3-dimensional discrete output space.

  * •

Pendulum: continuous control task where the agent has to balance a pendulum in
an upright position starting from a downward position. The observation space
is 3-dimensional and there is a single continuous output.

  * •

Foraging: a simple 1D grid with 5 cells and discrete 3-dimensional action
space corresponding to moving right, left, and staying still. The agent starts
in the center cell and can only observe which cell it is currently in. A food
source is randomly placed at one end of the grid. When the agent reaches the
food, it receives a reward of 10, and its position is reset. If the agent has
not found the food after 10 steps, the environment is also reset. Each time
the environment is reset there is a probability $p_{switch}$ (set to $0.5$ if
not stated otherwise) that the food location is changed, making this
environment non-stationary.

### Training

We optimize the model parameters on the different tasks using evosax (Lange,
2022a, ) implementation of the Covariance Matrix Adaptation Evolutionary
Strategy (CMA-ES) (Hansen,, 2023) with a population size of $128$ for all
tasks except for the foraging environment, for which we use a population size
of $256$. Fitness is obtained by averaging over three independent Monte-Carlo
trials. Models are optimized for $10,000$ generations, with all presented
results averaged over 5 different runs with different random seeds. Except for
the foraging task, we use the gymnax library implementation of the
environments (Lange, 2022b, ). We compare the performance of the same model
under different hyperparameters sets (see Table 1a) and with/without
structural plasticity. CMA-ES hyperparamters can be found in Table 1b.

Figure 3: Structural Plasticity. Changes during the agent lifetime of output
nodes’ in-degrees in the foraging task. Red and blue lines correspond to left
and right options respectively while shaded regions indicate the current
location of the reward (red for left and blue for right). We can see that new
connections are created towards the current best option while synapses towards
the other action are pruned.

## Results

Figure 4: LNDPs in CartPole. (a) Training curves showing average reward of the
population for each episode separately in CartPole. (b) Distribution of
rewards at each episode (from top to bottom) obtained by one evolved LNDP
trained on 5 episodes, episodes 5 to 9 (shown in red) are thus out of the
training distribution. (c) Evolution of a network through its SA phase. Node
colors code for the states $h^{t}$ and edge colors for weights $w^{t}$. Inputs
are on the left and outputs on the right. (d) 1,000 developmental trajectories
of a single evolved LNDP in the weights distribution space formed by the first
two statistical moments (mean and variance of the weights). Trajectories
correspond to different random initializations and colors code for the agent
return (bluer lines correspond to failures and greener to successes), with the
$\star$ indicating the starting point. (e) Trajectory of the network
activations $v^{t}$ through time (dimensions are reduced through PCA).

Results from the evolutionary runs are shown in Figure 2. Solutions were
discovered in all of the tested environments. We found differences between SP
models and non-SP ones in the Cartpole environment, which represent a
particularly difficult challenge: discovering sample efficient learning
algorithms from very few samples. In our setting, the number of steps the
agent experiences – and thus the number of update steps of the model –
directly depends on its skills as the episode immediately ends after the pole
has lost its balance. Consequently, a bad agent will see much fewer steps than
a good agent. This creates a hard optimization problem at the evolutionary
scale where an agent needs high sample efficiency for fast adaptation and this
sample efficiency must itself be discovered from very few samples.
Interestingly, while models without SP completely fail here, model endowed
with SP are able to find well-performing solutions (it should be noted that
while the average return over the population does not reach the highest
reward, i.e. 500, the maximum fitness per generation does).

When looking at the rewards obtained at each episode we observe that SP models
starting from empty networks reach higher rewards at the first episode and
higher overall fitness showing fast adaptation capabilities as networks are
able to self-organize into functional networks before the pole has lost its
balance (for indication, a random policy would lose the balance in around $20$
steps). While no evolved LNDP is able to always find a functional network in
the first episode, they are able to reach such a solution using lifetime
learning. It should be noted that models with SP have more parameters to
optimize ($1,854$ vs $1,468$ for non-SP models) and also have a much larger
output space as they describe dynamics over weights and topologies.

In the foraging environment with non-stationary rewards, we found that SP
models consistently reached higher average population fitness than non-SP
models while both reached similar maximum fitness. It should be noted that
looking at average fitness is more informative as maximum fitness can easily
be reached by chance. These results indicate that SP promotes better
adaptivity in non-stationary settings. Most of the found solutions go on the
side where the reward has been last obtained and go back to the other side if
it has not been found, which is an optimal strategy when the switching
probability is inferior or equal to $50\%$. Interestingly, these changes can
be implemented as weight changes or structural changes where excitatory
connections are created between hidden nodes and the output node corresponding
to the new best option as shown in Figure 3.

No significant differences have been found in Acrobot except in the case of
fully connected networks ($\mu^{conn}=1$) where models without SP reached
higher average fitness than models with SP. Finally, the pendulum environment,
which is the only continuous control problem, appeared as the most challenging
and the only one for which solutions where found only in a few runs and
uniquely for SP models starting from empty networks (i.e. $\mu^{conn}=0$) with
initialization variance $\sigma^{conn}=0.1$.

In the CartPole environment, it is particularly difficult for models without
SA to reach good performance in the first episode as the domain requires fast
adaptation for the pole to not directly fall over. Models with SA, on the
other hand, display innate skills in solving the task at the first episode,
demonstrating the ability of the model to reach target functional networks in
a reward-independent and self-organized manner (Figure 4). We also observe
that lifetime learning appears in later generations during training for models
with SA, the first $\sim 2,000$ generations being characterized by an increase
in innate skills only (Figure 4a). This could be explained by reduced pressure
on learning abilities. When training models without SA for much longer (i.e.
$50\times 10^{3}$ generations), models were found that can adapt quickly
enough to succeed in the first episode, showing increased sample efficiency.
Furthermore, we observed that while LNDPs without a SA phase tend to not
improve further after the second episode, models trained with SA display
stronger differences in between all episodes (Figure 4a). The latter
observation could indicate synergies between lifetime learning and SA-driven
development. We also compared results with an NDP, i.e. an LNDP not able to
modify its state during its lifetime. Results shown in Figure 5 demonstrate
that the LNDP outperforms the NDP in Foraging and CartPole. In the Foraging
environment, the NDP is not able to adapt at all which it could only achieve
by keeping track of environmental states in the activation patterns.

While we observed a high degree of variability across runs, in terms of
dynamics and operating mechanisms, some features were commonly observed.
First, networks tend to display short critical periods during which most of
the changes happen. Interestingly, for some models, we observed that failures
in the task were associated with abnormal bifurcations in the developmental
trajectory happening during this critical period as shown in Figure 4(d).
Having critical periods after which almost no change happens could represent a
trade-off between adaptivity and stability. In particular, we observed that
LNDPs evolved in the CartPole environment both with and without SA were very
consistent through time, even in episodes unseen during training (see Figure
4b where the agent is trained on 5 episodes and evaluated on 10) were some
solutions keep improving.

Drawing general conclusions about the discovered learning mechanisms is
difficult because of the high variability between seeds. While we have not
observed the emergence of noticeable structures in the grown networks, most of
them display less than fully-connected structures (Average network densities
at the end of training: Cartpole: $0.79\pm 0.23$, Foraging: $0.72\pm 0.16$,
Acrobot: $0.69\pm 0.22$, Pendulum: $0.24\pm 0.36$, densities are corrected to
account for impossible connections). Interestingly, we found that failed
training runs, mostly observed in the Pendulum environment, are characterized
by a fast convergence of the population towards producing empty networks,
which explains the lower average density found for Pendulum).

Figure 5: Performance of LNDP and NDP in Foraging and Cartpole environments.
The NDP is obtained by ablating the network updates during lifetime (i.e.
after the SA phase). Both models go through $100$ SA steps. While the
performance of both approaches is more similar in the CartPole domain, the
Foraging task requires an agent to adapt during its lifetime, which only the
LNDP is capable of.

## Discussion

Extending previous approaches on NDPs (Najarro et al.,, 2023) to lifetime
experience-dependant learning, we showed that LNDPs can be evolved to drive
the self-organization of ANNs into functional networks able to solve control
tasks with various properties. While structural plasticity is a poorly
addressed property in the context of artificial neural networks, we showed
evidence that it enables LNDPs to display better adaptivity in (1) an
environment with non-stationary dynamics requiring continuous adaptation, and
(2) in an environment requiring fast adaptation (with the additional
difficulty of having to learn this model from only a few experiences).

Spontaneous activity (SA) driven development, while being an important
mechanism in biological neural networks (Luhmann et al.,, 2016; Leighton and
Lohmann,, 2016; Kirkby et al.,, 2013), has not been well explored in ANNs (see
Raghavan and Thomson, (2019) for an example). We showed in this work that SA,
modeled as a simple stochastic process generating input neuron activity, can
be effectively used in concert with LNDPs to drive the network’s self-
organization in a pre-environmental phase allowing agents to display innate
skills and better adaptivity. Additionally, no extra component is needed
during the developmental phase; in other words, the same components drive the
network’s self-organization regardless of whether inputs come from SA or the
environment. This property enabled us to uncover synergies between the two
phases.

LNDPs represent a promising path toward self-assembling neural networks able
to adapt quickly. However, work remains to achieve performances comparable to
architectures traditionally used in the context of reinforcement learning.
While we have shown that LNDPs can be trained to solve and implement learning
algorithms in simple environments, future research must determine how they
scale to more complex ones with higher dimensions. As instances of complex
self-organized systems, their training represents a challenging optimization
problem characterized by many interactions over long periods of time.

While the architecture presented in this paper is intentionally general, using
versatile components such as GRUs and Graph Transformer layers, and
introducing very few inductive biases, future works should explore how to
efficiently constrain the LNDPs output space while keeping the possibility of
discovering unanticipated models. Inspiration could be drawn from
developmental neurosciences to introduce useful biases in learning and
developmental mechanisms. For instance, different synaptic plasticity rules
have been recently proposed leading to the emergence of properties usually
found in biological neural networks (Lynn et al.,, 2024; Achterberg et al.,,
2023). The general form of these rules could be used to guide the learning of
LNDPs.

Moreover, while we used classical evolutionary strategies together with
standard training methods, we believe that more focus should be given to this
aspect for more efficient training. Methods such as novelty search (Lehman and
Stanley,, 2011), quality diversity (Pugh et al.,, 2016) or intrinsically
motivated optimization (Laversanne-Finot et al.,, 2021) could represent
particularly relevant approaches for the problem at hand. Graph generation,
which is in our case an indirect task of the model, is a notoriously hard
problem necessitating large models and large amounts of data in a supervised
setting with much more informative feedback (Zhu et al.,, 2022). Moreover, we
found that training the wiring rules and the synaptic plasticity model made
the optimization particularly difficult because of the strong dependencies
between the two mechanisms.

Overall, we believe that by accounting for both pre-experience learning
through SA-driven development and experience-dependent adaptation, LNDPs
represent a promising path toward new types of self-organizing ANNs that can
narrow the gap between natural and artificial adaptivity in the future.

## Acknowledgements

This project was supported by a European Research Council (ERC) grant (GA no.
101045094, project ”GROW-AI”).

## References

  * Achterberg et al., (2023) Achterberg, J., Akarca, D., Strouse, D. J., Duncan, J., and Astle, D. E. (2023).  Spatially embedded recurrent neural networks reveal widespread links between structural and functional neuroscience findings.  Nature Machine Intelligence, 5(12):1369–1381.
  * Andrychowicz et al., (2016) Andrychowicz, M., Denil, M., Gómez, S., Hoffman, M. W., Pfau, D., Schaul, T., Shillingford, B., and de Freitas, N. (2016).  Learning to learn by gradient descent by gradient descent.  In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc.
  * Aswolinskiy and Pipa, (2015) Aswolinskiy, W. and Pipa, G. (2015).  RM-SORN: A reward-modulated self-organizing recurrent neural network.  Frontiers in Computational Neuroscience, 9.
  * Bednar and Miikkulainen, (1998) Bednar, J. A. and Miikkulainen, R. (1998).  Pattern-Generator-Driven Development in Self-Organizing Models.  In Bower, J. M., editor, Computational Neuroscience: Trends in Research, 1998, pages 317–323. Springer US, Boston, MA.
  * Bhattasali et al., (2022) Bhattasali, N., Zador, A. M., and Engel, T. (2022).  Neural Circuit Architectural Priors for Embodied Control.  Advances in Neural Information Processing Systems, 35:12744–12759.
  * Caroni et al., (2012) Caroni, P., Donato, F., and Muller, D. (2012).  Structural plasticity upon learning: Regulation and functions.  Nature Reviews Neuroscience, 13(7):478–490.
  * Chung et al., (2014) Chung, J., Gulcehre, C., Cho, K., and Bengio, Y. (2014).  Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling.
  * de Campos et al., (2011) de Campos, L. M. L., Roisenberg, M., and de Oliveira, R. C. L. (2011).  Automatic design of Neural Networks with L-Systems and genetic algorithms - A biologically inspired methodology.  In The 2011 International Joint Conference on Neural Networks, pages 1199–1206.
  * Dwivedi and Bresson, (2021) Dwivedi, V. P. and Bresson, X. (2021).  A Generalization of Transformer Networks to Graphs.
  * Ferigo and Iacca, (2023) Ferigo, A. and Iacca, G. (2023).  Self-building neural networks.  In Proceedings of the Companion Conference on Genetic and Evolutionary Computation, pages 643–646.
  * Ginsburg and Jablonka, (2021) Ginsburg, S. and Jablonka, E. (2021).  Evolutionary transitions in learning and cognition.  Philosophical Transactions of the Royal Society B: Biological Sciences, 376(1821):20190766.
  * Ha et al., (2022) Ha, D., Dai, A. M., and Le, Q. V. (2022).  HyperNetworks.  In International Conference on Learning Representations.
  * Ha and Tang, (2022) Ha, D. and Tang, Y. (2022).  Collective intelligence for deep learning: A survey of recent developments.  Collective Intelligence, 1(1):26339137221114874.
  * Hansen, (2023) Hansen, N. (2023).  The CMA Evolution Strategy: A Tutorial.
  * Hochreiter et al., (2001) Hochreiter, S., Younger, A. S., and Conwell, P. R. (2001).  Learning to Learn Using Gradient Descent.  In Dorffner, G., Bischof, H., and Hornik, K., editors, Artificial Neural Networks — ICANN 2001, pages 87–94, Berlin, Heidelberg. Springer.
  * Kirkby et al., (2013) Kirkby, L. A., Sack, G. S., Firl, A., and Feller, M. B. (2013).  A Role for Correlated Spontaneous Activity in the Assembly of Neural Circuits.  Neuron, 80(5):1129–1144.
  * Kirsch and Schmidhuber, (2022) Kirsch, L. and Schmidhuber, J. (2022).  Meta Learning Backpropagation And Improving It.
  * Kudithipudi et al., (2022) Kudithipudi, D., Aguilar-Simon, M., Babb, J., Bazhenov, M., Blackiston, D., Bongard, J., Brna, A. P., Chakravarthi Raja, S., Cheney, N., Clune, J., Daram, A., Fusi, S., Helfer, P., Kay, L., Ketz, N., Kira, Z., Kolouri, S., Krichmar, J. L., Kriegman, S., Levin, M., Madireddy, S., Manicka, S., Marjaninejad, A., McNaughton, B., Miikkulainen, R., Navratilova, Z., Pandit, T., Parker, A., Pilly, P. K., Risi, S., Sejnowski, T. J., Soltoggio, A., Soures, N., Tolias, A. S., Urbina-Meléndez, D., Valero-Cuevas, F. J., van de Ven, G. M., Vogelstein, J. T., Wang, F., Weiss, R., Yanguas-Gil, A., Zou, X., and Siegelmann, H. (2022).  Biological underpinnings for lifelong learning machines.  Nature Machine Intelligence, 4(3):196–210.
  * (19) Lange, R. T. (2022a).  Evosax: JAX-based Evolution Strategies.
  * (20) Lange, R. T. (2022b).  gymnax: A JAX-based reinforcement learning environment library.
  * Laversanne-Finot et al., (2021) Laversanne-Finot, A., Péré, A., and Oudeyer, P.-Y. (2021).  Intrinsically Motivated Exploration of Learned Goal Spaces.  Frontiers in Neurorobotics, 14.
  * Lazar et al., (2009) Lazar, A., Pipa, G., and Triesch, J. (2009).  SORN: A Self-Organizing Recurrent Neural Network.  Frontiers in Computational Neuroscience, 3.
  * Lehman and Stanley, (2011) Lehman, J. and Stanley, K. O. (2011).  Novelty Search and the Problem with Objectives.  In Riolo, R., Vladislavleva, E., and Moore, J. H., editors, Genetic Programming Theory and Practice IX, pages 37–56. Springer New York, New York, NY.
  * Leighton and Lohmann, (2016) Leighton, A. H. and Lohmann, C. (2016).  The Wiring of Developing Sensory Circuits—From Patterned Spontaneous Activity to Synaptic Plasticity Mechanisms.  Frontiers in Neural Circuits, 10.
  * Lowell and Pollack, (2016) Lowell, J. and Pollack, J. (2016).  Developmental Encodings Promote the Emergence of Hierarchical Modularity.  In ALIFE 2016, the Fifteenth International Conference on the Synthesis and Simulation of Living Systems, pages 344–351. MIT Press.
  * Luhmann et al., (2016) Luhmann, H. J., Sinning, A., Yang, J.-W., Reyes-Puerta, V., Stüttgen, M. C., Kirischuk, S., and Kilb, W. (2016).  Spontaneous Neuronal Activity in Developing Neocortical Networks: From Single Cells to Large-Scale Interactions.  Frontiers in Neural Circuits, 10.
  * Lynn et al., (2024) Lynn, C. W., Holmes, C. M., and Palmer, S. E. (2024).  Heavy-tailed neuronal connectivity arises from Hebbian self-organization.  Nature Physics, 20(3):484–491.
  * May, (2011) May, A. (2011).  Experience-dependent structural plasticity in the adult human brain.  Trends in Cognitive Sciences, 15(10):475–482.
  * Najarro and Risi, (2020) Najarro, E. and Risi, S. (2020).  Meta-Learning through Hebbian Plasticity in Random Networks.  In Advances in Neural Information Processing Systems, volume 33, pages 20719–20731. Curran Associates, Inc.
  * Najarro et al., (2022) Najarro, E., Sudhakaran, S., Glanois, C., and Risi, S. (2022).  HyperNCA: Growing Developmental Networks with Neural Cellular Automata.
  * Najarro et al., (2023) Najarro, E., Sudhakaran, S., and Risi, S. (2023).  Towards Self-Assembling Artificial Neural Networks through Neural Developmental Programs.  In ALIFE 2023: Ghost in the Machine: Proceedings of the 2023 Artificial Life Conference. MIT Press.
  * Parisi et al., (2019) Parisi, G. I., Kemker, R., Part, J. L., Kanan, C., and Wermter, S. (2019).  Continual lifelong learning with neural networks: A review.  Neural Networks, 113:54–71.
  * Pedersen et al., (2024) Pedersen, J. W., Plantec, E., Nisioti, E., Montero, M., and Risi, S. (2024).  Structurally Flexible Neural Networks: Evolving the Building Blocks for General Agents.
  * Pugh et al., (2016) Pugh, J. K., Soros, L. B., and Stanley, K. O. (2016).  Quality Diversity: A New Frontier for Evolutionary Computation.  Frontiers in Robotics and AI, 3.
  * Raghavan and Thomson, (2019) Raghavan, G. and Thomson, M. (2019).  Neural networks grown and self-organized by noise.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
  * Risi, (2021) Risi, S. (2021).  The future of artificial intelligence is self-organizing and self-assembling.  sebastianrisi.com.
  * Sandler et al., (2021) Sandler, M., Vladymyrov, M., Zhmoginov, A., Miller, N., Madams, T., Jackson, A., and Arcas, B. A. Y. (2021).  Meta-Learning Bidirectional Update Rules.  In Proceedings of the 38th International Conference on Machine Learning, pages 9288–9300. PMLR.
  * Soltoggio et al., (2018) Soltoggio, A., Stanley, K. O., and Risi, S. (2018).  Born to learn: The inspiration, progress, and future of evolved plastic artificial neural networks.  Neural Networks, 108:48–67.
  * Stanley and Miikkulainen, (2003) Stanley, K. O. and Miikkulainen, R. (2003).  A Taxonomy for Artificial Embryogeny.  Artificial Life, 9(2):93–130.
  * Valsalam et al., (2007) Valsalam, V. K., Bednar, J. A., and Miikkulainen, R. (2007).  Developing Complex Systems Using Evolved Pattern Generators.  IEEE Transactions on Evolutionary Computation, 11(2):181–198.
  * Wang et al., (2018) Wang, J. X., Kurth-Nelson, Z., Kumaran, D., Tirumala, D., Soyer, H., Leibo, J. Z., Hassabis, D., and Botvinick, M. (2018).  Prefrontal cortex as a meta-reinforcement learning system.  Nature Neuroscience, 21(6):860–868.
  * Zhu et al., (2022) Zhu, Y., Du, Y., Wang, Y., Xu, Y., Zhang, J., Liu, Q., and Wu, S. (2022).  A Survey on Deep Graph Generation: Methods and Applications.  In Proceedings of the First Learning on Graphs Conference, pages 47:1–47:21. PMLR.