\section{Introduction} \label{sec:intro}

Star formation occurs both within low-mass dark clouds and inside large molecular
cloud complexes, but the role radiation plays in the two types of regions
differs markedly \citep[e.g.][for a review]{tan14}.  Unlike their more quiescent counterparts,
giant molecular clouds typically form O and B stars, and ultraviolet
radiation from these massive stars injects enough energy into the surrounding
medium to create large H~II regions. The resulting volumes of expanding ionized gas remove
material that might otherwise accrete onto a protostar, while radiation photoevaporates
the outer portions of disks and may alter the spectrum of
stellar masses created in the region \citep{winter18,krumholz11}.  On the other hand,
density waves driven into molecular clouds by radiation from massive stars
could also potentially trigger new stars to form, either directly as the fronts
compress existing clumps into gravitationally
unstable cores, or indirectly if the fronts increase the overall density in a cloud to the
point where subsequent collisions between cloud fragments induce collapse \citep[e.g.][]{haworth11}.
Radiative feedback and winds from newly-formed stars into the nascent cloud material also
drive structures on scales of tens or even hundreds of pc, and create
sequential episodes of star formation \citep[e.g.][]{venuti18} and clear out chimneys and bubbles
perpendicular to the galactic plane \citep{rob18}. On smaller scales,
how massive stars form remains a topic of intense theoretical work, as
models must include the effects of strong radiation pressure and photoevaporation
from massive protostars that both limit the available timescale for accretion
and may widen existing outflow cavities \citep{tanaka17,rosen16}.

To understand quantitatively how intense radiation fields from young massive stars affect star formation
in molecular clouds, we must assess parsec scales, where filaments respond to the
influence of large-scale magnetic fields, smaller scales of order 0.1~pc relevant to
core formation and fragmentation, and protostellar disk scales of several hundred AU \citep[e.g.][]{li18}.
The most intense radiation environments occur in regions that produce early O-stars, and
even the closest of these are over 1.5~kpc away, so to resolve structures on $\sim$ 100~AU
scales requires spatial resolutions better than 0.1\arcsec . This increase in resolution
of an order of magnitude over typical non-adaptive-optics images holds the key to bridging
interstellar scales to Oort cloud scales.

Recent advances in instrumentation, especially with adaptive optics (AO), have made it
possible to acquire nearly diffraction-limited images in the near-infrared over
a field of view greater than an arcminute in size. Such observations are ideal
for tracing large-scale outflows and shock
waves in star-forming regions where subarcsecond resolution uncovers morpologies
crucial to interpreting the physics at work in these systems
\citep[e.g.][]{bally15}. Wide-field near-infrared AO images presage what
the James Webb Space Telescope (JWST) should deliver when it becomes operational,
and complement the highest-resolution ALMA maps of molecular cloud filaments and cores,
which also have $\lesssim$ 1\arcsec\ resolution \citep[e.g.][]{cheng18}.

In this paper we use the GSAOI adaptive optics
imager on Gemini to study the Western Wall (G287.38-0.62) in Carina,
one of the nearest \citep[2.3 kpc; ][]{lim19}
and most strongly-irradiated interfaces known to occur in a region of massive star
formation. This interface appears as an unassuming dark cloud at optical wavelengths,
but it stands out as the brightest feature in the entire Carina complex when imaged
in the 2.12 $\mu$m H$_2$ and Br-gamma emission lines \citep[][Fig.~\ref{fig:brg}]{hartigan15}.
The H$_2$ 2.12 $\mu$m line is excited by UV fluorescence, and provides
a means to trace irradiated interfaces in regions like Carina where UV emission lines
are absorbed by dust along the line of sight.  The observed spatial offsets
between H$_2$ and Br-$\gamma$ emission, the latter
arising in a photoevaporative flow as UV radiation
ionizes hydrogen, generally agree well with theoretical model predictions \citep{ch18}.
Over a dozen O-stars in the nearby open cluster Tr~14 and a similar number in the cluster Tr~16,
including an O2 star (HD 93129A) and an O3.5 star (HD~93129B) in Tr~14 and an O2.5 star (HD~93162)
and O3.5 star (ALS~15210) in Tr~16, contribute to the FUV and EUV flux at the
Western Wall (Fig.~\ref{fig:brg}), as do other O and B stars scattered throughout the region \citep{gagne11,wu18}.
The Wall has a particularly favorable geometry in that
it is a convex structure, which makes it easier to observe interface shapes than in
bowl-like concave interfaces like the Orion Bar \citep{goi19,goi16,pelligrini09}.

In what follows, we describe our image processing steps and
combine the Br-$\gamma$, H$_2$, and K-cont AO images into a color composite in Sec.~2.
The images reveal a new H$_2$ outflow and several new features along the interface, including long
ridges, sinusoidal-like waves, and fragments that we discuss in Sec.~3,
where we consider what the new observations tell us about instabilities, triggering, and
fragmentation processes.  A summary of the work is given in Sec.~4.

\section{Image Processing} \label{sec:data}

We acquired images of Carina's Western Wall with the Adaptive Optics Imager
(GSAOI) on Gemini South between Jan 15, 2018 and Jan 19, 2018 for a program with ten hours of queue time.
The images used three narrowband filters,
H$_2$ 1-0 S(1) (G1121; $\lambda_\circ$ = 2.122$\mu$m; $\Delta\lambda$ = 0.032$\mu$m),
Br-$\gamma$ (G1122; $\lambda_\circ$ = 2.166$\mu$m; $\Delta\lambda$ = 0.032$\mu$m), and
K-Continuum (G1112; $\lambda_\circ$ = 2.270$\mu$m; $\Delta\lambda$ = 0.034$\mu$m), and
employed two pointings offset in declination from one another by 85\arcsec\
to cover the region of interest. The center of the composite image is
$\alpha$ = 10:43:30.6, $\delta$ = $-$59:35:20 (epoch 2000).  Within each pointing
we used multiple dithers of up to $\pm$ 8\arcsec\
to remove hot pixels and to fill in the $\sim$ 3\arcsec\ gaps
between the four Rockwell 2048 $\times$ 2048 arrays. Exposure times per dither
for the H$_2$, Br-$\gamma$, and Cont-K filters were 120, 120, and 100 seconds, respectively,
with 8 non-destructive reads for each exposure. The total exposure times in
the northern position were 1320, 1560, and 1400 seconds, respectively, for
H$_2$, Br-$\gamma$, and Cont-K, while the corresponding exposure times for the southern
position were 1920, 1680, and 1200 seconds. We took sets of dithered sky frames without AO-correction
and with the same exposure times as the object in
a relatively blank field located about 20 arcminutes away from the object.
In all, the data set includes 56 sky frames, enabling us to correct for
changing atmospheric backgrounds in each filter as described below.

We combined dome and sky flats with the IRAF packages {\it gemcombine} and {\it gemarith}. For each
object frame we used the IRAF command {\it gemexpr} to flatfield the data, and
subtracted a sky composite acquired on the same night in that filter.
After using the `disco-stu' routines developed by the Gemini staff to correct for distortion,
we combined the dithers for each filter together with a median filter and
corrected for zero-point offsets in the background by
devising our own bad pixel masks for the data to eliminate residual
edge effects and bad columns. The two pointings then combine into a single
image by merging their data in the region of overlap. The resulting images have a plate
scale of 0.0195\arcsec\ per pixel, and cover a region 95\arcsec\ in RA by 173\arcsec\ in DEC.
The FWHM of the point-spread-function of the final images varies between 0.065\arcsec\ near
the center of each of the pointings, to 0.110\arcsec\ near the edges, where stellar
images become slightly distorted. These numbers compare well with the theoretical Rayleigh limit
of 0.069\arcsec\ for the Gemini 8-m telescope at 2.2$\mu$m. For comparison, JWST should
achieve a resolution of 0.085\arcsec\ at this wavelength.

\begin{figure}
\plotone{fig1.pdf}
\caption{Br-gamma image of the central portion of the Carina star formation region
from \citep{hartigan15}. Locations of the Trumpler 14, 15, and 16 clusters are circled. The boxed
area shows the region imaged in Fig.~\ref{fig:overview}.}
\label{fig:brg}
\end{figure}

\begin{figure}
\plotone{fig2.pdf}
\caption{Images of Carina's Western Wall. Left: A non-AO image of the region
from \citep{hartigan15} taken in 0.8\arcsec\ seeing, with H$_2$ in red, Br-$\gamma$ in green and
[O~III] $\lambda$5007 in blue. White boxes identify the regions expanded in Fig.~\ref{fig:zoom}. The scale bar is 0.1~pc
(8.97\arcsec).
Right: The new adaptive optics images acquired with GSAOI, where H$_2$ is in red, Br-$\gamma$ in green and K-cont
in blue. The FWHM of the stellar images ranges from 0.065\arcsec\ to 0.11\arcsec.}
\label{fig:overview}
\end{figure}

To improve the signal-to-noise we rebinned the AO images by 3$\times$3 pixels, so the images
we present here have square pixels of size 0.0585\arcsec. There are 80 stars in
the GAIA-DR2 catalog that have unsaturated profiles in our composite image, and the scatter in
their positions was rather large, $\pm$ 0.3\arcsec. Using the IRAF routines {\it geomap} and {\it geotran}
we redid the distortion corrections and were able to reduce the scatter in the
GAIA-DR2 stellar coordinates by a factor of 15, to $\pm$ 0.02\arcsec\ in each of the three filters.
Fig.~\ref{fig:overview} compares our final AO composite with a ground-based non-AO image of the same
region.

\section{Discussion}

\subsection{Structure and Outflows Within The Western Wall}

The new AO images resolve previously unseen structure in the Western Wall
on scales between the seeing limit of the ground-based images (0.8\arcsec; 1840~AU) and
the diffraction limit of Gemini at 2.2$\mu$m (0.065\arcsec; 150~AU). The composite image in
Fig.~\ref{fig:overview} illustrates that this size range is particularly rich in detail
across the field of view, especially in H$_2$, which traces
the photodissociation fronts.  In the H$_2$ image, the apex area of the Western Wall
appears as a series of long ridges that run along the interface, with projected
spacings of $\sim$ 2000~AU (Fig.~\ref{fig:zoom}).
A large triangular-shaped dark clump marked `Nose'
protrudes from the wall near the apex.
The northern boundary of the Wall appears more chaotic than the apex does, and
is home to a series of swept-back fragments of width $\sim$ 1800~AU. To the south,
the Wall exhibits a remarkable set of wave-like interfaces, including
one that appears almost perfectly sinusoidal with a half-wavelength of 2500~AU.

From the color composite in Fig.~\ref{fig:overview} one might easily get the
impression that the Wall is transparent where H$_2$ emission is absent,
because Br-$\gamma$ is present between the filamentary H$_2$ emission.
However, in optical images the entire region where H$_2$ radiates is opaque,
so it is likely that the Br-$\gamma$ located within the Wall
comes from photoevaporation from its near-side in the direction of
the observer. If this interpretation is correct, this Br-$\gamma$ emission
lies superposed upon the H$_2$, and should be blueshifted relative to the
Br-$\gamma$ emission to the east of the Wall by $\sim$ 10 km$\,$s$^{-1}$,
the sound speed of the photoevaporative flow.

\begin{figure}
\plotone{fig3.pdf}
\caption{Expanded views of sections of the H$_2$ adaptive optics image of Carina's Western Wall. Scale
bars are 5000~AU (2.17\arcsec), assuming a distance of 2.3~kpc.  Left: The apex of the Wall consists of
3 $-$ 4 long parallel ridges and a large triangular-shaped protuberance labeled `Nose'.
Top-right: The northern edge of the Wall reveals several swept-back fragments along
the boundary of the object. Top-middle: An outflow is defined by a
series of H$_2$ bow shocks.  Bottom-right panels: The southern edge of the Wall
has a complex, and at times sinusoidal-like boundary and several overlapping arcs.}
\label{fig:zoom}
\end{figure}

In addition to tracing fluorescence along irradiated interfaces, H$_2$ also emits
in shock waves along bipolar flows from stellar jets.
Our images reveal that an H$_2$ feature noted previously
as a possible shocked object \citep{tapia06}, resolves into
four nested bow shocks, labeled B, C, D, and E in Figure 3,
and a short linear feature, denoted A, located upstream (east)
of the bow shocks.  These objects differ from the irradiated interfaces in that they have no
associated Br-$\gamma$ or K-cont emission (pure red color in Fig~\ref{fig:overview}).
The region also stands out from the other H$_2$ emission in the non-AO ground-based images
\citep[][ left panel of Fig.~\ref{fig:overview}]{tapia06,hartigan15}, but with much less clarity.

The driving source of the jet is unknown. There is a narrow dark gap between knot
A and a faint nebula $\sim$ 0.6\arcsec\ to the east that could represent an
opaque disk, and another possibility for a driving source is the star labeled 'neb'
in Fig~\ref{fig:zoom} located 11.5\arcsec\ at PA 66.5 degrees from knot A. This star
has a conical-shaped reflection nebula aligned in the direction of the H$_2$ bow shocks,
but no nebulae connect the star with the H$_2$ shocks. The Vela-Carina
legacy survey with Spitzer did not detect any mid-infrared point
sources at the position of knot A, though the nebulous star is visible faintly
at 3.6 $\mu$m, 4.5 $\mu$m, and 5.8 $\mu$m.
MHO~1630 is not obviously associated with a nearby filament of
faint optical emission lines \citep{tapia06}, or with any
sources visible in Spitzer or Herschel images \citep[e.g.][]{tapia15}.
Future ALMA mapping of the continuum and molecular line emission
in the vicinity of MHO~1630 would clarify whether or not protostars exist there.
Table~1 presents coordinates for several of the features shown in Fig.~\ref{fig:zoom}.

\begin{deluxetable}{lcc}
\tablewidth{0pt}
\tablenum{1}
\tablecolumns{5}
\tablecaption{Coordinates for Objects in Fig.~2}
\tablehead{}
\startdata
\noalign{\smallskip}
Name&$\alpha$ (2000)&$\delta$ (2000)\\
neb&10:43:29.548&$-$59:34:10.695\\
MHO-1630 A&10:43:28.179&$-$59:34:15.154\\
MHO-1630 B&10:43:27.860&$-$59:34:15.570\\
MHO-1630 C&10:43:27.367&$-$59:34:16.132\\
MHO-1630 D&10:43:27.294&$-$59:34:16.805\\
MHO-1630 E&10:43:27.195&$-$59:34:16.578\\
Nose $^a$&10:43:35.813&$-$59:34:48.752\\
Wave $^{b, c}$&10:43:31.684&$-$59:34:33.942\\
Bay $^c$&10:43:31.772&$-$59:36:02.958\\
\noalign{\smallskip}
\enddata
\tablenotetext{a}{Apex location}
\tablenotetext{b}{Southern of the two waves marked in Fig.~\ref{fig:zoom}}
\tablenotetext{c}{Center of feature}
\end{deluxetable}

\subsection{Theoretical Implications}

The AO images in Figs.~\ref{fig:overview} and~\ref{fig:zoom}
allow us to observe the irradiated interfaces in the Western
Wall with unprecedented resolution, with implications for models of
triggered star formation.  In a simple one-dimensional model, we
expect the D-front associated with an irradiated interface to increase
the density in the cloud as the ionization front propagates, leading generally to
conditions that are more favorable to gravitational collapse. In agreement
with this general picture, ALMA observations
along the southern portion of the Western Wall show that the densest portions of
the molecular cloud lie along the dissociation front \citep{hartigan20b}.
In three dimensions, compression behind a D-front
becomes more complex.  On a timescale of $\sim$ 10$^5$ yrs, as a
front passes through clumps it will compress them,
and shadowing leaves behind a pillar with a dense knot at the apex which might then collapse
to form a new star \citep[e.g. Fig.~2 of][]{mackey10}.
Compression is highest when the front wraps around the clumps \citep{tremblin12},
and the external magnetic field can also influence the overall shape of the pillar \citep{henney09}.

If stars are indeed forming all along the Wall in response to the D-front, we would expect
pillars to trail behind the newly-formed stars away from the direction of
radiation.  However, in contrast to the pillar-rich southern region of the Carina Nebula \citep{smith00},
the AO images do not show any indication that the Western Wall breaks into pillars
on size scales down to $\sim$ 150~AU.  Instead, we observe more of a smooth, slightly
wavy and sometimes ridged surface along much of the Wall down to the resolution of the images.
The `Nose' in Fig.~\ref{fig:zoom} is the only feature where one can make a
case that radiation may be wrapping around a clump.  While this feature is striking and
may form one or more pillars at some point, it is currently a rather broad structure,
and it is unclear whether or not protostars or multiple cores exist there and, if
so, whether those clumps predated the passage of the front.
The available ALMA data on the Western Wall cover only its southern parts, and
the Nose has not been observed yet. Future high-resolution millimeter observations
might reveal the presence of protostars in this region.

The regularly-spaced ridges that run parallel to the photodissociation front
are a striking feature of the Western Wall.  Their geometry implies a preferred
direction for the ridges, characteristic of a dynamically significant large-scale
magnetic field.  Radiation-MHD simulations of H~II regions show that there should
be a large-scale magnetic field in the neutral gas that aligns along the
direction of the D-front \citep{arthur11}.
Comparing with the results of \citet{mackeylim11}, we see that their model R8
most closely resembles the ridges we see here. These simulations imply there should be
a large-scale dynamically significant magnetic field oriented northeast to southwest
in the Western Wall, parallel to the ridges. This prediction can be tested by
measuring the dust polarization along the Wall.

The flocculent nature of the northern edge of the Wall differs markedly
from the generally smoother, ridged morphology present near the apex.
If the northern edge were strongly sculpted by the incident radiation
from Tr 14, then the dense structures here should form pillars that
point back towards the radiation sources to the east.  Instead, the northern
interface resembles one subject to Kelvin-Helmholtz (KH)
instabilities \citep{frank96,jones12}. If the
orientation of the global magnetic field is approximately constant, and
aligned as described above, then the field will be nearly perpendicular to the northern
edge and thus will not tend to stabilize the system against the KH instability.
The remaining ingredient for the KH instability is, of course, a shear flow.
Radiation from nearby O stars in Tr~14 and Tr~16 is capable of causing just such a
shear flow along this edge and, indeed, the radiation itself may destabilize the
system further \citep[e.g.][]{shadmehri12}.

On a final note, it is particularly intriguing that the projected width of the
ridges referred to above, the scale of the density structures along the northern
edge of the Wall, and the wavelength of the sinusoidal structure referred to in
Sect 3.1 are all approximately 2000 AU in size. This is well within the
scales we would expect to be influenced by ambipolar diffusion in the undisturbed
molecular cloud, and therefore below the lengthscales expected to be produced by
turbulence \citep[e.g.][]{vanloo08,dos11,downes12,xu19}. Ambipolar
diffusion is usually significant only in weakly ionised regions and it is
therefore tempting to suppose that structures of this scale are created by flows
at the surface of the molecular cloud, where the ionisation fraction is high,
rather than being ``uncovered'' as photoevaporation of the cloud proceeds.

\section{Summary}

The H$_2$, Br-$\gamma$ and K-cont AO images Carina's Western Wall
presented in this paper reveal several unexpected structures within
this classic irradiated molecular cloud that are difficult or impossible
to discern from non-AO ground-based images. The Wall has a convex geometry so
the interfaces appear as sharp tangent lines across the cloud, and are
ideal for resolving interface shapes that provide clues to MHD instabilities.
The apex of the Wall, which points in the general direction of the irradiating
sources, consists of series of ridges that follow along the front. Wave-like
structures also occur along the front. The most chaotic morphologies occur along
the northern and southern edges of the Wall, where pillar-like fragments are
present. These fragments do not point towards the irradiating sources, but
are swept back as expected in a shear flow.
The new images also resolve the H$_2$ object MHO~1630 into an outflow
that consists of a series of nested bow shocks without a clear driving source.

The observed ridging in the Wall argues for the presence of a large-scale
dynamically-important magnetic field, while the flocculent appearance along
the northern and southern boundaries is a characteristic of
Kelvin-Helmholtz instabilities. We find no direct evidence that the ionization
front triggers new stars to form, in that there are no obvious cases where the
front wraps around a dense clump to form a protostar. Nonetheless, the
ionization front does encompass the entire Wall, and increases the
density along the surface where radiation is absorbed.

These observations make it clear that size scales between $\sim$ 200~AU and
2000~AU are rich with detail for an irradiated dark cloud, and in many ways
hold the key to understanding the complex dynamical processes that sculpt these
systems. The ridges, fragments, and waves we found all have sizes in this range,
and require subarcsecond resolution to resolve even for the nearest
regions of massive star formation. This capability now currently exists with
AO imaging, but should expand greatly once JWST becomes operational.
The interplay between large-scale, global physics such as magnetic field
strengths and direction, rotation, and Jeans lengths, with the small scale three-dimensional
morphologies and complex motions induced by radiative feedback and turbulence
continues to make star formation a challenging and engaging field of study.

\acknowledgments
We thank Rodrigo Carrasco of the Gemini staff for his assistance with observation
scheduling and for routines we used for initial distortion corrections, and
D.~Froebrich for the official new MHO number.
Based on observations obtained for program GS-2018A-Q-123
at the international Gemini Observatory, a program of NSF’s NOIRLab,
which is managed by the Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory
partnership: the National Science Foundation (United States), National Research Council (Canada),
Agencia Nacional de Investigaci\'on y Desarrollo (Chile), Ministerio de Ciencia, Tecnologia e Innovacion
(Argentina), Minist\'erio da Ciência, Tecnologia, Inovacoes e Comunicacoes (Brazil), and Korea
Astronomy and Space Science Institute (Republic of Korea).

\vspace{5mm}
\facilities{Gemini (GSAOI)}
\clearpage


\section{Introduction}

Over the past decade, quantum plasmonics \cite{maier-np13} underwent a rapid development fueled by a host of recently discovered phenomena such as strong exciton-plasmon coupling effects \cite{bellessa-prl04,sugawara-prl06,wurtz-nl07,fofang-nl08,hakala-prl09,gomez-nl10}, plasmon-assisted hot carrier generation \cite{park-nl11,melosh-nl11,halas-science11}, plasmonic laser (spaser) \cite{bergman-prl03,stockman-natphot08,noginov-nature09,zhang-nature09,gwo-science12}, plasmon tunneling \cite{baumberg-nature12,aizpurua-nl12,dionne-nl13,tan-science14} and more, along with a growing number of applications. Surface plasmons are collective electron excitations living at the metal-dielectric interfaces which can interact strongly with light and localized electron excitations such as excitons in molecules or semiconductors \cite{stockman-review}. Although classical description of many experiments in terms of local field enhancement has  largely been successful, a growing number of topics and applications require a more rigorous quantum approach \cite{hohenester-prb08,nordlander-nl09,stockman-jo10,waks-pra10,garcia-prl13,garcia-optica17,aizpurua-nanophot20}.

In nanoscale systems, the local fields can change strongly over the length scale well below the diffraction limit, and  so the plasmon interactions with the electromagnetic (EM) field and  excitons depend sensitively on the system  parameters such as the geometry of a metal-dielectric structure and the exciton position relative to it. While  the coupling parameters, characterizing these interactions, have been suggested in several forms by using analogy with the cavity modes \cite{khitrova-nphys06,andreani-prb12,shahbazyan-nl19,mortensen-rpp20}, these parameters have yet to emerge within a consistent quantization approach for interacting  plasmons.  Yet another longstanding  challenge for quantum plasmonics is to  account properly  for strong optical dispersion and losses in metals that give rise to non-Markovian effects \cite{nori-prb09,tejedor-prb10,thanopulos-prb17,moradi-sr18,molmer-acsph19}.

Within the canonical quantization scheme, localized plasmon modes with discrete frequency  spectrum $\omega_{m}$   are described by the Hamiltonian
\begin{equation}
\label{H-plasmon}
\hat{H}_{\rm pl}=\sum_{m}\hbar\omega_{m}\hat{a}^{\dagger}_{m}\hat{a}_{m},
\end{equation}
where $\hat{a}^{\dagger}_{m}$ and $\hat{a}_{m}$ are, respectively, the plasmon creation and annihilation operators obeying the canonical commutation relations $[\hat{a}_{m},\hat{a}^{\dagger}_{n}]=\delta_{mn}$. Plasmon interactions with quantum emitters (QEs), modeled hereafter by two-level systems, are usually described by the Jaynes-Cummings  interaction Hamiltonian
\begin{equation}
\label{H-pl-qe}
\hat{H}_{\rm pl-qe}=\sum_{im}\hbar g_{im}(\hat{\sigma}^{\dagger}_{i}\hat{a}_{m}+\hat{a}^{\dagger}_{m}\hat{\sigma}_{i}),
\end{equation}
where  $\hat{\sigma}^{\dagger}_{i} $ and $\hat{\sigma}_{i} $ are, respectively, the raising and lowering operators for the $i$th QE (related to the Pauli matrices for Fermionic systems), while $g_{im}$ is the QE-plasmon coupling, which, within the canonical approach, is an \textit{ad hoc}  parameter. Although widely employed, the canonical scheme has significant limitations when used in metal-dielectric structures characterized by a complex dielectric function $\varepsilon (\omega,\bm{r})=\varepsilon' (\omega,\bm{r})+i\varepsilon'' (\omega,\bm{r})$ since it ignores the medium optical dispersion and, hence, is unsuitable for describing non-Markovian effects in plasmonics.

On the other hand, the medium optical dispersion is accurately accounted for in the macroscopic electrodynamics approach based upon the fluctuation-dissipation (FD) theorem \cite{welsch-pra98,welsch-p00,philbin-njp10}. In this framework, the EM fields are quantized in terms of reservoir noise operators $\hat{\bm{f}}(\omega,\bm{r})$  driven by the Hamiltonian
\begin{equation}
\label{H-noise}
\hat{H}_{N}= \!\int_{0}^{\infty}\! d\omega \!\int\! dV \,\hbar \omega \hat{\bm{f}}^{\dagger}(\omega,\bm{r})\!\cdot\!\hat{\bm{f}}(\omega,\bm{r})
\end{equation}
and obeying the commutation relations
\begin{equation}
[\hat{\bm{f}}(\omega,\bm{r}),\hat{\bm{f}}^{\dagger}(\omega',\bm{r}')]=\bm{I}\delta(\omega-\omega')\delta(\bm{r}-\bm{r}'),
\end{equation}
where $\bm{I}$ is the unit tensor. Interactions with QEs are described by the Hamiltonian term
$\hat{H}_{\rm int}=-\sum_{i}\hat{\bm{p}}_{i} \cdot \hat{\bm{E}}(\bm{r}_{i})$,
where $\hat{\bm{p}}_{i}$ and $\hat{\bm{E}}(\bm{r}_{i})$ are, respectively,  the QE dipole moment and  electric field operators.  The latter is given by
\begin{equation}
\label{field-dissip}
\hat{\bm{E}}(\bm{r})=\!\int_{0}^{\infty}\! d\omega \! \!\int\! dV'\! \bm{D}(\omega;\bm{r},\bm{r}')\hat{\bm{P}}_{N}(\omega,\bm{r}') + \text{H.c.},
\end{equation}
where $\hat{\bm{P}}_{N}(\omega,\bm{r}) = (i/2\pi)\sqrt{\hbar\varepsilon'' (\omega,\bm{r})}\hat{\bm{f}}(\omega,\bm{r})$ is the noise polarization vector operator and  $\bm{D}(\omega;\bm{r},\bm{r}')$ is the EM dyadic Green function defined as
\begin{align}
\label{green-em}
\bm{\nabla}\!\times\! \bm{\nabla}\!\times\!\bm{D}(\omega;\bm{r},\bm{r}')-\frac{\omega^{2}}{c^{2}}\varepsilon(\omega,\bm{r}) \bm{D}(\omega;\bm{r},\bm{r}')~~~~~~
\nonumber
\\
 =\frac{4\pi\omega^{2}}{c^{2}} \bm{I}\delta(\bm{r}-\bm{r}').
\end{align}
The macroscopic FD approach has been extensively used to model spontaneous emission,  strong coupling effects and non-Markovian dynamics in metal-dielectric structures \cite{welsch-pra00,dzsotjan-prb10,hughes-prb12,andreani-prb12,hughes-prb13,zubairy-prb14,garcia-prl14,rousseaux-prb18,sivan-prb19,nori-prb09,tejedor-prb10,thanopulos-prb17,moradi-sr18,molmer-acsph19}.  Its major drawback in relation to plasmonics is that, even though surface plasmons reside primarily at the metal-dielectric interfaces, the eigenstates of $\hat{H}_{N}$ extend over the entire system reservoir, i.e., the Hilbert space,  spanned by the operators $\hat{\bm{f}}(\omega,\bm{r})$, is excessively large. Furthermore, the plasmons only  appear  as resonances in the classical EM Green function $\bm{D}$, so that, in practical terms, this approach is limited to systems that allow  $\bm{D}$ to be evaluated analytically or numerically.

In principle, the  Hamiltonians (\ref{H-plasmon}) and (\ref{H-pl-qe}), along with the canonical commutations relations and QE-plasmon coupling, should be obtained within the macroscopic FD framework starting with a suitable mode expansion for the EM Green function that would define the basis set for normal mode expansion of the electric field operator (\ref{field-dissip}) \cite{varguet-ol16,dzsotjan-pra16,hughes-prl19}. However, for systems of general shape, a straightforward expansion of $\bm{D}$ over a  discrete  set of EM modes that include radiation, such as quasinormal modes \cite{hughes-njp14,hughes-acsphot14,lalanne-pra14,hughes-pra15,muljarov-prb16,lalanne-prb18,lalanne-lpr18}, gives rise to  dissipation-induced coupling between the modes which compicates the   commutations relations \cite{hughes-prl19}. On the other hand, on the length scale well below the diffraction limit, the plasmon interactions with QEs take place primarily via the near field coupling rather than  photon exchange, while the rate of nonradiative losses, determined by $\varepsilon'' (\omega,\bm{r})$, by far exceeds the plasmon radiative decay rate. This implies that in nanoscale systems, the plasmons should be treated as primarily electronic excitations  interacting with the EM field and QEs, and so the quantization of plasmons should be carried on its own, i.e., separate from the radiation field. Furthermore, since the electron motion in metals is unretarded, the plasmon modes in such systems can be described within quasistatic approximation \cite{stockman-review}, which allows, as we show below,  one to define the orthogonal basis set  that is free of  dissipation-induced coupling.

In this paper, we develop a consistent quantum approach to localized plasmons  interacting with quantum emitters and the electromagnetic field  which accurately accounts for   medium optical dispersion and losses.  First,  starting within the macroscopic FD framework, we employ the near-field plasmon Green function \cite{shahbazyan-prl16,shahbazyan-prb18} do define the normal mode expansion for the plasmon field operators free of dissipation-induced coupling. Using this basis set, we explicitly obtain the plasmon Hamiltonian (\ref{H-plasmon}), along with the equal-time canonical commutations relations, and show that the canonical quantization scheme is valid  only in the Markov approximation, i.e., if the dielectric function $\varepsilon (\omega,\bm{r})$ is replaced by its value at the plasmon frequency (i.e., $\omega=\omega_{m}$). In this way, we also obtain the microscopic coupling parameters that define the plasmon interactions with the EM field and QEs in terms of the plasmon local fields, system geometry and QE positions relative to the plasmonic structure. Second,  moving beyond the Markov approximation, we present a new approach that bridges the macroscopic and canonical quantization schemes while accounting accurately for the medium optical dispersion and losses. In this approach, quantum plasmons are described in terms of a discrete set of bosonic modes with linear dispersion which reside at the metal-dielelectric interfaces, and whose operators span a reduced Hilbert space obtained by projecting the full reservoir states upon localized plasmon modes. The interactions of such projected reservoir modes  with the EM field and QEs are mediated by \textit{classical} plasmons, so that the classical plasmon enhancement effects are encoded in the Hamiltonian coupling parameters.

\section{From macroscopic to canonical quantization of surface plasmons}

In this section, starting within the macroscopic quantization scheme \cite{welsch-pra98,welsch-p00,philbin-njp10}, we  obtain the canonical Hamiltonian (\ref{H-plasmon}) and the  equal-time commutation relations for localized surface plasmon in metal-dielectric structures of arbitrary shape.

\subsection{Quasistatic modes and plasmon Green function}
We consider a metal-dielectric structure characterized by dielectric function  of the form $\varepsilon (\omega,\bm{r})=\sum_{i}\theta_{i}(\bm{r})\varepsilon_{i}(\omega)$, where $\theta_{i}(\bm{r})$ is  unit step function that vanishes outside the connected region, metal or dielectric, of volume $V_{i}$ with uniform dielectric function $\varepsilon_{i}(\omega)$.  For unretarded electron motion in metals, the potentials $\Phi_{m}(\bm{r})$ and frequencies  $\omega_{m}$  of plasmon  modes   are determined by the quasistatic Gauss law as \cite{stockman-review}
\begin{equation}
\label{gauss}
\bm{\nabla}\cdot\left [\varepsilon' (\omega_{m},\bm{r})\bm{\nabla} \Phi_{m}(\bm{r})\right ]=0.
\end{equation}
Accordingly, the plasmon mode fields,  which we choose to be real,  are defined as $\bm{E}_{m}(\bm{r})=-\bm{\nabla} \Phi_{m}(\bm{r})$. Importantly, the plasmon mode fields   are orthogonal in \textit{each}  region $V_{i}$,
\begin{equation}
\label{superorthogonality}
\int\! dV_{i} \bm{E}_{m}(\bm{r})\!\cdot\!\bm{E}_{n}(\bm{r})=\delta_{mn}\int\! dV_{i} \bm{E}_{m}^{2}(\bm{r}),
\end{equation}
so that $\int \! dV\varepsilon''(\omega,\bm{r})\bm{E}_{m}(\bm{r})\bm{E}_{n}(\bm{r})=0$ for $m\neq n$, implying no dissipation-induced coupling (see Appendix).

The near-field Green function that defines the  field operator  (\ref{field-dissip}) can be split into free-space and plasmon parts as \cite{shahbazyan-prl16,shahbazyan-prb18} $\bm{D}=\bm{D}_{0}+\bm{D}_{\rm pl}$. When inserted into  Eq.~(\ref{field-dissip}), the first term yields the electric field due to noise fluctuations, while the second term defines  the normal mode expansion of the plasmon field operator. In this paper, we  focus only on the plasmonic sector of the Hilbert space. In the absence of dissipation-induced coupling, the plasmon Green function  can be derived \textit{exactly} in the following form (see Appendix):
\begin{equation}
\label{gauss-green-plasmon-full}
\bm{D}_{\rm pl}(\omega;\bm{r},\bm{r}')=\sum_{m}D_{m}(\omega)\bm{E}_{m}(\bm{r}) \bm{E}_{m}  (\bm{r}')
\end{equation}
where
\begin{align}
\label{gauss-green-plasmon}
D_{m}(\omega) =
\dfrac{4\pi}{\int\! dV \bm{E}_{m}^{2}(\bm{r})}
-  \dfrac{4\pi}{\int\! dV \varepsilon (\omega,\bm{r})\bm{E}_{m}^{2}(\bm{r})}.
\end{align}
The first  term ensures that $\bm{D}_{\rm pl}=0$ for  $\varepsilon=1$ (or, in the limit $\omega\rightarrow\infty$). The plasmon Green functions exhibits poles in the complex frequency plane defined by the condition $\int\! dV \varepsilon(\omega,\bm{r})\bm{E}_{m}^{2}(\bm{r})=0$. In the following, we restrict ourselves to the plasmonics frequency domain  $\varepsilon''(\omega)/\varepsilon'(\omega)\ll 1$. Since  $\int\! dV \varepsilon'(\omega_{m},\bm{r})\bm{E}_{m}^{2}(\bm{r})=0$ due to the Gauss law, we can expand  $\varepsilon' (\omega,\bm{r})$ in Eq.~(\ref{gauss-green-plasmon}) near $\omega_{m}$ to present the plasmon Green function as a sum over the plasmon poles \cite{shahbazyan-prb18} (see Appendix)
\begin{align}
\label{dyadic-plasmon}
\bm{D}_{\rm pl}(\omega;\bm{r},\bm{r}')
=
\sum_{m}
\frac{\omega_{m}}{4 U_{m}}\frac{\bm{E}_{m}(\bm{r}) \bm{E}_{m}  (\bm{r}')}{\omega_{m}-\omega -\frac{i}{2}\gamma_{m}(\omega)},
\end{align}
where $U_{m}$ is the plasmon mode energy \cite{landau},
\begin{align}
\label{energy-mode}
U_{m}
= \frac{1}{16\pi}
\!\int \!  dV
\dfrac{\partial[\omega_{m}\varepsilon'(\omega_{m},\bm{r})]}{\partial \omega_{m}}
\, \bm{E}_{m}^{2}(\bm{r}),
\end{align}
and $\gamma_{m}(\omega)$ is the frequency-dependent decay rate
\begin{equation}
\label{mode-decay}
\gamma_{m}(\omega)=\dfrac{2\!\int \!  dV
\varepsilon''(\omega,\bm{r})\bm{E}_{m}^{2}(\bm{r}) }{\!\int \!  dV
[\partial\varepsilon'(\omega_{m},\bm{r})/\partial \omega_{m}]
 \bm{E}_{m}^{2}(\bm{r}) }.
\end{equation}
In structures with a single metallic component,

the decay rate takes the form \cite{stockman-review} $\gamma_{m}(\omega)=2\varepsilon''(\omega)/[\partial\varepsilon'(\omega_{m})/\partial \omega_{m}]$. Finally, with help of Eqs.~(\ref{superorthogonality}) and (\ref{mode-decay}), we obtain the following relation
\begin{align}
\label{FDT}
\int \! dV  \varepsilon''(\omega,\bm{r}) \bm{D}_{\rm pl}^{*}(\omega;\bm{r},\bm{r}')\bm{D}_{\rm pl}(\omega;\bm{r},\bm{r}'')
~~~~~~~~~~
\nonumber\\
=4\pi \text{Im}\bm{D}_{\rm pl}(\omega;\bm{r}',\bm{r}''),
\end{align}
which ensures consistency with the FD theorem \cite{welsch-pra98,welsch-p00,philbin-njp10}.

\subsection{Plasmon Hamiltonian and canonical commutation relations}

The normal mode expansion of the plasmon field operator is obtained upon inserting the plasmon Green function (\ref{dyadic-plasmon}) into Eq.~(\ref{field-dissip}):
$\hat{\bm{E}}_{\rm pl}(\bm{r})=\sum_{m}\hat{\bm{E}}_{m}(\bm{r})$, where
\begin{align}
\label{field-plasmon-normal}
 \hat{\bm{E}}_{m}(\bm{r})
=
\sqrt{\frac{\hbar\omega_{m}}{4 U_{m}}}\,\bm{E}_{m}(\bm{r})   (\hat{a}_{m}+\hat{a}^{\dagger}_{m}),
\end{align}
is  the field operator for an individual mode. Here, $\hat{a}_{m}$ is the plasmon annihilation operator defined as
\begin{equation}
\label{a-plasmon}
\hat{a}_{m}=i\!\int_{0}^{\infty}\! \frac{d\omega}{\sqrt{2\pi}}\frac{\hat{f}_{m}(\omega)}{\omega-\omega_{m} +\frac{i}{2}\gamma_{m}(\omega)},
\end{equation}
where $\hat{f}_{m}$ is  noise operator \textit{projected} on plasmon mode:
\begin{align}
\label{f-proj}
\hat{f}_{m}(\omega)
&=i\sqrt{\frac{\pi\omega_{m}}{2\hbar U_{m}}}\!\int\!dV \bm{E}_{m}(\bm{r}) \! \cdot \! \hat{\bm{P}}_{N}(\omega,\bm{r})
\\
&=-\sqrt{\frac{\omega_{m}}{8\pi U_{m}}}\!\int\! dV\!\sqrt{\varepsilon'' (\omega,\bm{r})}
\,\bm{E}_{m}(\bm{r}) \! \cdot \! \hat{\bm{f}}(\omega,\bm{r}).
\nonumber
\end{align}
Commutation relations for  $\hat{f}_{m}$ follow from those for $\hat{\bm{f}}$ and from Eqs.~(\ref{superorthogonality}) and (\ref{mode-decay}),
\begin{equation}
\label{f-comm}
[\hat{f}_{m}(\omega),\hat{f}_{n}^{\dagger}(\omega')]=\delta_{mn}\delta(\omega-\omega') \gamma_{m}(\omega),
\end{equation}
where commutativity of operators at $m\neq n$ is insured by the absence of dissipation-induced coupling. Now, using Eqs.~(\ref{a-plasmon}) and (\ref{f-comm}), we obtain the commutation relations for plasmon operators as
\begin{equation}
\label{a-comm}
[\hat{a}_{m}, \hat{a}_{n}^{\dagger}]=\delta_{mn} \!\int_{0}^{\infty}\!  \frac{d\omega }{2\pi}
\frac{ \gamma_{m}(\omega)}{(\omega_{m}-\omega)^{2} +\gamma_{m}^{2}(\omega)/4}.
\end{equation}
In the Markov approximation, by replacing $\gamma_{m}(\omega)$ with $\gamma_{m}\equiv \gamma_{m}(\omega_{m})$  and extending the integral to negative frequencies, we obtain  the  canonical commutation relations  $[\hat{a}_{m}, \hat{a}_{n}^{\dagger}]=\delta_{mn}$. The  plasmon Hamiltonian (\ref{H-plasmon}) follows from  the normal mode expansion (\ref{field-plasmon-normal})   by verifying that the normal-ordered Hamiltonian for individual  modes has the form
\begin{equation}
\label{H-plasmon-mode}
\hat{H}_{m}=\frac{1}{8\pi}
\!\int \!  dV
\dfrac{\partial
(\omega_{m}\varepsilon')
}{\partial \omega_{m}}
 \hat{\bm{E}}_{m}^{2}
 =\hbar\omega_{m} \hat{a}_{m}^{\dagger} \hat{a}_{m},
\end{equation}
where we dropped the terms $\hat{a}_{m} \hat{a}_{m}$ and $\hat{a}_{m}^{\dagger} \hat{a}^{\dagger}_{m}$ and disregarded the zero-point energy. The factor $1/2$ difference between Eqs.~(\ref{energy-mode}) and (\ref{H-plasmon-mode}) reflects the presence of both positive and negative frequency terms in $\hat{\bm{E}}_{m}(\bm{r})$ \cite{landau}. We stress that, by using the  plasmon Green function (\ref{dyadic-plasmon}), \textit{both}  the Hamiltonian (\ref{H-plasmon})  and canonical commutation relations have been explicitly obtained.

Turning to the plasmon dynamics, the time-evolution of projected noise operators (\ref{f-proj}) is determined by the Heisenberg equations,
\begin{equation}
\label{noise-dynamics}
\dot{\hat{f}}_{m}(\omega) =-(i/\hbar)[\hat{f}_{m}(\omega), \hat{H}_{N}]=-i\omega \hat{f}_{m}(\omega),
\end{equation}
where the dot stands for   time derivative. From this relation and Eq.~(\ref{a-plasmon}), the Heisenberg  equations  for   plasmon operators  readily follow (in the Markov approximation),
\begin{equation}
\label{langevin}
\dot{\hat{a}}_{m}(t)  =-(\gamma_{m}/2+i\omega_{m}) \hat{a}_{m}(t) +\hat{f}_{m}(t),
\end{equation}
where
$\hat{f}_{m}(t)\! = \! (2\pi)^{-1/2}\!\int_{0}^{\infty}\!  d\omega  \hat{f}_{m}(\omega)e^{-i\omega t}$ is time-domain projected noise operator. The commutation relations for $\hat{f}_{m}(t)$ are obtained from Eq.~(\ref{f-comm}) as
\begin{equation}
\label{f-comm-time}
[\hat{f}_{m}(t),\hat{f}_{n}^{\dagger}(t')]= \delta_{mn}\gamma_{m}\delta(t-t'),
\end{equation}
where the Markov approximation was used again. Thus, the Markovian dynamics of plasmon operators $\hat{a}_{m}(t)$ is described by quantum Langevin equation (\ref{langevin}) with  white-noise source $\hat{f}_{m}(t)$,  which  guarantees \cite{scully-book}  the equal-time  commutation relations: $[\hat{a}_{m}(t), \hat{a}_{n}^{\dagger}(t)]=\delta_{mn}$.

\section{Plasmon interactions with  quantum emitters and the electromagnetic field}

Consider now the plasmon coupling to QEs and the EM field.  In contrast to cavity modes, the plasmons are  localized at the scale well below the  diffraction limit and, therefore,  interact with the EM field $\bm{\mathcal{E}}(t)$ in a way similar to point dipoles. The interaction Hamiltonian has the form  $H_{\rm pl-em}=-\sum_{m}\hat{\bm{p}}_{m}\cdot \bm{\mathcal{E}}(t)$, where   $\hat{\bm{p}}_{m}=\!\int\! dV\hat{\bm{P}}_{m}(\bm{r})$ is the plasmon dipole moment  and $\hat{\bm{P}}_{m}(\bm{r})$ is the polarization vector operators.  To determine $\hat{\bm{P}}_{m}(\bm{r})$, we present the Gauss law (\ref{gauss}) as $\bm{\nabla}\!\cdot\!\bm{E}_{m}(\bm{r})+4\pi \bm{\nabla}\!\cdot\!\bm{P}_{m}(\bm{r})=0$, where $\bm{P}_{m}(\bm{r})=\chi' (\omega_{m},\bm{r})\bm{E}_{m}(\bm{r})$ is the  mode   polarization vector and $\chi(\omega,\bm{r})$ is the  system  susceptibility. In the Markov approximation, converting this relation to  operator form as $\hat{\bm{P}}_{m}(\bm{r})= \chi' (\omega_{m},\bm{r})\hat{\bm{E}}_{m}(\bm{r})$ and using Eq.~(\ref{field-plasmon-normal}), for the EM field of the form $\bm{\mathcal{E}}(t)=\bm{\mathcal{E}}e^{-i\omega_{L} t}+\bm{\mathcal{E}}^{*}e^{i\omega_{L} t} $ that is uniform on the system scale, we obtain in the rotating wave approximation (RWA)
\begin{equation}
\label{H-pl-em}
H_{\rm pl-em}=-\sum_{m}\left (\bm{\mu}_{m}\cdot\bm{\mathcal{E}}e^{-i\omega_{L} t} \, \hat{a}_{m}^{\dagger} + {\rm H.c.}\right ),
\end{equation}
where  $\bm{\mu}_{m}\equiv \bm{\mu}_{m}(\omega_{m})$, and we introduced the frequency-dependent transition matrix element,
\begin{equation}
\label{mode-matrix-element}
\bm{\mu}_{m}(\omega)=\frac{1}{2}\sqrt{\frac{\hbar\omega_{m}}{U_{m}}}
\!\int\! dV\chi' (\omega,\bm{r})\bm{E}_{m}(\bm{r}).
\end{equation}
to be used in the following section. The scaling factor $\sqrt{\hbar\omega_{m}/U_{m}}$ in Eq.~(\ref{mode-matrix-element})   converts the plasmon  energy $U_{m}$  to $\hbar\omega_{m}$ in order to match the energy of the EM field. With matrix element (\ref{mode-matrix-element}), the plasmon radiative decay rate is given by the standard expression
\begin{equation}
\label{mode-decay-rad}
\gamma_{m}^{rad}= W_{m}^{rad}/U_{m}=(4\omega_{m}^{3} \mu_{m}^{2})/(3\hbar c^{3}),
\end{equation}
where $W_{m}^{rad}=p_{m}^{2}\omega_{m}^{4}/3c^{3}$ is the power radiated by dipole $\bm{p}_{m}$ \cite{scully-book}.

We now turn to  the interactions  between plasmons and QEs   situated at positions $\bm{r}_{i}$  with excitation frequency $\omega_{e}$ and dipole moments $\hat{\bm{p}}_{i}= \bm{\mu}_{i}(\hat{\sigma}^{\dagger}_{i}+\hat{\sigma}_{i})$, where $\bm{\mu}_{i}=\mu_{e} \bm{n}_{i}$   is the transition matrix element ($\bm{n}_{i}$ is the dipole orientation).  Employing the mode expansion (\ref{field-plasmon-normal}) in the interaction Hamiltonian $\hat{H}_{\rm pl-qe}=-\sum_{i}\hat{\bm{p}}_{i}  \cdot \hat{\bm{E}}_{\rm pl}(\bm{r}_{i})$, we obtain, in RWA, the coupling Hamiltonian   (\ref{H-pl-qe})  with   $g_{im}$ given by
\begin{equation}
\label{coupling}
\hbar g_{im}
=
-\sqrt{\frac{\hbar\omega_{m}}{4U_{m}}}\,\bm{\mu}_{i}\!\cdot\!\bm{E}_{m}(\bm{r}_{i}).
\end{equation}
Using Eq.~(\ref{energy-mode}), the QE-plasmon coupling $g_{im}$ can be recast in the form similar to cavities
\begin{equation}
\label{coupling-mode-volume}
g^{2}_{im}
=\frac{2\pi \mu_{e}^{2}\omega_{m}}{\hbar{\cal V}^{(i)}_{m}},
~~
\frac{1}{{\cal V}^{(i)}_{m}}
= \frac{2[\bm{n}_{i}\!\cdot\!\bm{E}_{m}(\bm{r}_{i})]^{2}}{\int \! dV [\partial (\omega_{m}\varepsilon')/\partial \omega_{m}]\bm{E}_{m}^{2}},
\end{equation}
where ${\cal V}^{(i)}_{m}$ is the projected  plasmon mode volume  \cite{shahbazyan-acsphot17,shahbazyan-prb18}, which characterizes  the plasmon field confinement at a point $\bm{r}_{i}$ in the direction $\bm{n}_{i}$. Since the Gauss equation (\ref{gauss}) is scale-invariant \cite{stockman-review},  the coupling parameters  (\ref{mode-matrix-element}) and  (\ref{coupling}) are independent of the overall field normalization. By rescaling the fields as $\tilde{\bm{E}}_{m}(\bm{r})=\sqrt{\hbar\omega_{m}/4U_{m}} \bm{E}_{m}(\bm{r})$, these parameters are brought to a more familiar form
\begin{equation}
\label{coupling-normalized}
g_{im}=-\bm{\mu}_{i}\!\cdot\!\tilde{\bm{E}}_{m}(\bm{r}_{i})/\hbar,
~~~
\bm{\mu}_{m}=\!\int\! dV\chi' (\omega_{m},\bm{r})\tilde{\bm{E}}_{m}(\bm{r}).
\end{equation}

To illustrate the sensitivity of QE-plasmon coupling to the system parameters, we  present the results of numerical calculations of the coupling parameter $g_{em}$, given by Eq.~(\ref{coupling-mode-volume}), for a single QE situated at a distance $d$ from the tip of an Au nanorod in water (see Fig.~\ref{fig1}).  The nanorod   was modeled by a prolate spheroid with semi-major and semi-minor axes $a$ and $b$, respectively, at aspect ratio $a/b=3.0$, the standard spherical harmonics were used for calculating the longitudinal plasmon fields, the QE dipole orientation was chosen along the nanorod symmetry axis,  and the Au experimental dielectric function was used in all calculations. The calculated QE-plasmon coupling is normalized by full plasmon decay rate $\gamma_{m}$, which includes the nonradiative decay rate (\ref{mode-decay}) and radiative decay rate (\ref{mode-decay-rad}) taken at plasmon frequency $\omega=\omega_{m}$. The coupling parameter  behaves as $g_{em}\propto 1/\sqrt{{\cal V}_{m}}$, where, for small nanostructures, the plasmon mode volume ${\cal V}_{m}$ scales as the metal volume \cite{shahbazyan-prb18,shahbazyan-nl19}, resulting in the enhanced coupling for smaller nanorods. The coupling sharply increases as  QE approaches hot spot at the nanorod tip characterized by very large plasmon field confinement (small mode volume).

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.95\columnwidth]{fig1}
\caption{\label{fig1} The QE-plasmon coupling  $g_{em}$ for a  QE near the tip of an Au nanorod in water is plotted against the distance to the tip for several values of the nanorod overall size. Inset: Schematics of a QE near  the Au nanorod tip.
 }
\end{center}
\vspace{-3mm}
\end{figure}

Finally, the full canonical Hamiltonian for plasmons interacting with the EM field and QEs has the form $H=H_{\rm pl}+H_{\rm pl-qe}+H_{\rm pl-em}+H_{\rm qe}+H_{\rm qe-em}$, where we  included the standard QE Hamiltonian $H_{\rm qe}=\hbar\omega_{e}\!\sum_{i}\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i}$  and the coupling term $H_{\rm qe-em}=-\sum_{i}(\bm{\mu}_{i}\!\cdot\!\bm{\mathcal{E}}e^{-i\omega_{L} t} \hat{\sigma}_{i}^{\dagger} + {\rm H.c.})$ describing  the QEs' interactions with the EM field.   We stress that for plasmons, the canonical quantization scheme is valid only in the Markov approximation that ignores the dielectric function dispersion.  However, in metal-dielectric structures, the role of materials optical dispersion can be very significant, and so a quantum description that includes such effects is necessary, as we discuss in the following section.

\section{Beyond Markov approximation: Intermediate quantization scheme}

In this section, we present yet another plasmon quantization approach that bridges the macroscopic and canonical schemes while combining the advantages of both. In this "intermediate" scheme, the Hilbert space is restricted to plasmonic degrees of freedom, in contrast to macroscopic approach, while the medium optical dispersion  is still accurately accounted for,  in contrast to the canonical scheme.

\subsection{Recasting quantum plasmons via projected reservoir modes}

The Langevin equation (\ref{langevin}) with the source (\ref{f-proj}) imply that each plasmon mode is driven  by a fraction of the  reservoir that   overlaps with the plasmon  electric field.  These projected reservoir modes (PRM) form a discrete subspace of the full reservoir Hilbert space  spanned by the operators $\hat{b}_{m}(\omega)=\hat{f}_{m}(\omega)/\sqrt{\gamma_{m}(\omega)}$ or,  using Eqs.~(\ref{mode-decay})  and  (\ref{f-proj}),
\begin{equation}
\label{b-definition}
\hat{b}_{m}(\omega)
=-\frac{\int\! dV\sqrt{\varepsilon'' (\omega,\bm{r})}
\,\bm{E}_{m}(\bm{r}) \! \cdot \! \hat{\bm{f}}(\omega,\bm{r})}
{\left [\int\! dV\varepsilon'' (\omega,\bm{r})
\,\bm{E}_{m}^{2}(\bm{r})\right ]^{1/2}},
\end{equation}
satisfying the commutation relations
\begin{equation}
\label{b-comm}
[\hat{b}_{m}(\omega),\hat{b}_{n}^{\dagger}(\omega')]=\delta_{mn}\delta(\omega-\omega').
\end{equation}
Time-evolution in the reduced Hilbert space  is driven by the PRM Hamiltonian  with linear dispersion
\begin{equation}
\label{H-reduced-noise}
\hat{H}_{\rm b}=\sum_{m}  \!\int_{0}^{\infty}\!   d\omega \,\hbar\omega \,\hat{b}^{\dagger}_{m}(\omega)\hat{b}_{m}(\omega),
 \end{equation}
which leads to accurate Heisenberg equations  of the form $\dot{\hat{b}}_{m}(\omega)=-(i/\hbar)[\hat{b}_{m}(\omega), \hat{H}_{\rm b}]=-i\omega\hat{b}_{m}(\omega)$ [compare to  Eq.~(\ref{noise-dynamics})]. The PRMs and plasmons can be set as independent variables by adding the  coupling Hamiltonian term:
$\hat{H}_{\rm pl-b}=i\hbar\sum_{m} \! \int_{0}^{\infty}\! \! d\omega \kappa_{m}(\omega) [\hat{a}^{\dagger}_{m}\hat{b}_{m}(\omega) -\hat{b}_{m}^{\dagger}(\omega)\hat{a}_{m}]$, where  $\kappa_{m}(\omega)=\sqrt{\gamma_{m}(\omega)/2\pi}$ is  the plasmon-PRM coupling. Then, upon tracing the PRMs out, one would arrive, in the standard way, at the master equation for density matrix \cite{scully-book}. Here, we chose a different approach and instead describe the system \textit{directly} in terms of PRMs.

The interaction Hamiltonian between PRMs and QEs is obtained from the QE-plasmon coupling term (\ref{H-pl-qe}) by using the relation  (\ref{a-plasmon}) between  the plasmon and PRM operators, with $\hat{f}_{m}(\omega)=\hat{b}_{m}(\omega)\sqrt{\gamma_{m}(\omega)}$. We then obtain
\begin{equation}
\label{H-b-qe}
\hat{H}_{\rm b-qe}=\sum_{im} \int_{0}^{\infty}\!  d\omega\left [\hbar q_{im}(\omega)\,\hat{\sigma}^{\dagger}_{i}\,\hat{b}_{m}(\omega) +\text{H.c.}
\right ],
\end{equation}
where $q_{im}(\omega)$ is the QE-PRM coupling,
\begin{align}
\label{coupling-prm-qe}
q_{im}(\omega)
= \sqrt{\frac{\gamma_{m}(\omega)}{2\pi}}\frac{ i g_{im} }{\omega -\omega_{m}+\frac{i}{2}\gamma_{m}(\omega)},
\end{align}
with $g_{im}$ given by Eq.~(\ref{coupling}).

The PRM coupling to the EM field $\bm{\mathcal{E}}(t)$ that is uniform on the system scale is described by the Hamiltonian $H_{\rm int}=-\text{Re} \int dV \hat{\bm{E}}_{\rm pl}(\bm{r})  \cdot  \bm{P}(t,\bm{r})$, where $\bm{P} =\hat{\chi} \bm{\mathcal{E}}$ is the induced polarization vector. For a monochromatic field,  using Eqs.~(\ref{field-plasmon-normal}) and (\ref{a-plasmon}), we obtain
\begin{equation}
\label{H-b-em}
\hat{H}_{\rm b-em}= -\!\!\sum_{m}\! \int_{0}^{\infty}\!\! \! d\omega \! \left [\bm{d}_{m}^{*}(\omega_{L},\omega)\!\cdot \!\bm{\mathcal{E}}e^{-i\omega_{L}t}\, \hat{b}_{m}^{\dagger}(\omega) +\text{H.c.}
\right ],
\end{equation}
where $\bm{d}_{m}(\omega_{L},\omega)$ is the optical transition matrix element for PRMs [compare to Eq. (\ref{coupling-prm-qe})],
\begin{align}
\label{coupling-prm-em}
\bm{d}_{m}(\omega_{L},\omega)
= \sqrt{\frac{\gamma_{m}(\omega)}{2\pi}}\frac{i\bm{\mu}_{m}(\omega_{L})}{\omega-\omega_{m}+\frac{i}{2}\gamma_{m}(\omega)},
\end{align}
and the plasmon transition matrix element $\bm{\mu}_{m}(\omega)$ is given by Eq. (\ref{mode-matrix-element}).

\subsection{Classical plasmonic enhancement effects and first-order transition probability rates}

To elucidate the mechanism behind the QE-PRM interaction, let us evaluate the first-order transition probability rate for PRM excitation by a QE. For an excited QE with energy $\hbar \omega$, the transition rate for an individual PRM has the form
\begin{equation}
\label{rate-abs}
\Gamma_{im}(\omega)=\frac{2\pi}{\hbar}
\int_{0}^{\infty}\!  d\omega'
\left | \hbar q_{im}(\omega')\right |^{2}\delta(\hbar\omega'-\hbar\omega),
\end{equation}
where the  integration is taken over the PRM's final states. Evaluating the  integral, we obtain
\begin{equation}
\label{rate-et-mode}
\Gamma_{im} (\omega)
= 2\pi \left | q_{im}(\omega) \right |^{2}=\frac{g_{im}^{2} \gamma_{m}(\omega)}{(\omega_{m}-\omega)^{2} +\frac{1}{4}\gamma_{m}^{2}(\omega)},
\end{equation}
where we used Eq.~(\ref{coupling-prm-qe}). The full transition rate is obtained by summing Eq.~(\ref{rate-et-mode}) over all PRMs: $\Gamma_{i} (\omega)=\sum_{m}\Gamma_{im} (\omega)$. In fact, the above expression represents the rate of energy transfer (ET) from a QE to plasmons evaluated, in a standard way,   using the classical plasmon Green function (\ref{dyadic-plasmon}) with help of Eq.~(\ref{coupling})  \cite{shahbazyan-prl16},
\begin{equation}
\label{rate-et}
\Gamma_{i} (\omega)
=\dfrac{2}{\hbar}\, \text{Im} \left [\bm{\mu}_{i} \bm{D}_{\rm pl}(\omega;\bm{r}_{i},\bm{r}_{i})\bm{\mu}_{i} \right ] =\sum_{m}\Gamma_{im} (\omega),
\end{equation}
which indicates that QE-PRM interactions are \textit{mediated by classical plasmons} absorbing the QE energy. Thus, the classical effect of resonance ET from a QE to plasmons emerges from the Hamiltonian (\ref{H-b-qe}) in the lowest order.

Turning to PRM interactions with the EM field,  the transition probability rate for excitation of a PRM  by the incident monochromatic  light of frequency $\omega$ is
\begin{equation}
\label{rate-abs-mode}
\Gamma_{m}(\omega)=\frac{2\pi}{\hbar}
\!\int_{0}^{\infty}\!  d\omega'
\left |  \bm{d}_{m}(\omega,\omega')\!\cdot \!\bm{\mathcal{E}}\right |^{2}\delta(\hbar\omega'-\hbar\omega),
\end{equation}
which, after evaluating the integral, can be presented as [compare to Eqs.~(\ref{rate-et}) and (\ref{rate-et-mode})]
\begin{equation}
\label{rate-abs-pol}
\Gamma_{m} (\omega)
=\dfrac{2\pi}{\hbar^{2}}\left |  \bm{d}_{m}(\omega,\omega)\!\cdot \!\bm{\mathcal{E}}\right |^{2}
=\dfrac{2}{\hbar}\, \text{Im} \left [\bm{\mathcal{E}}^{*}\bm{\alpha}_{m}(\omega)\bm{\mathcal{E}} \right ].
\end{equation}
Here, $\bm{\alpha}_{m} (\omega)$ is the optical polarizability tensor of a plasmon mode that defines its response to an external field \cite{shahbazyan-prb18} (see Appendix):
\begin{equation}
\label{polar-mode}
\bm{\alpha}_{m} (\omega)=\frac{1}{\hbar}\frac{\bm{\mu}_{m}(\omega)\bm{\mu}_{m}(\omega)}{\omega_{m}-\omega -\frac{i}{2}\gamma_{m}(\omega)}.
\end{equation}
The full absorption rate is obtained by summing over all PRM modes  as $\Gamma_{\rm pl} (\omega)
=(2/\hbar)\text{Im} \left [\bm{\mathcal{E}}^{*}\bm{\alpha}_{\rm pl}(\omega)\bm{\mathcal{E}} \right ]$, where $\bm{\alpha}_{\rm pl}(\omega)=\sum_{m}\bm{\alpha}_{m}(\omega)$ is the full optical polarizability of plasmonic structure. Within RWA, the absorbed power is given by $P(\omega)=\hbar\omega \Gamma_{\rm pl}(\omega)/4=(\omega/2)\text{Im} \left [\bm{\mathcal{E}}^{*}\bm{\alpha}_{\rm pl}(\omega)\bm{\mathcal{E}} \right ]$. Upon normalizing it by the incident flux $S_{0}=(c/8\pi) \mathcal{E}^{2}$, we obtain the absorption cross-section as $\sigma_{\rm pl}^{\rm abs}=\sum_{m}\sigma_{m}^{\rm abs}$, where $\sigma_{m}^{\rm abs}(\omega)$ is the absorption cross-section for an individual mode,
\begin{equation}
\label{cross-abs-mode}
\sigma_{m}^{\rm abs}
=\dfrac{4 \pi \omega}{c} \text{Im} \left [\bm{\varepsilon}\bm{\alpha}_{m}(\omega)\bm{\varepsilon} \right ]
=\dfrac{2 \pi \omega}{\hbar c}\frac{\bm{|\varepsilon}\!\cdot\!\bm{\mu}_{m}(\omega)|^{2}\gamma_{m}(\omega)}{(\omega-\omega_{m})^{2}+\frac{1}{4}\gamma_{m}^{2}(\omega)},
\end{equation}
and $\bm{\varepsilon}$ is the incident light polarization.Thus, the PRM-EM interaction is mediated by  plasmon  absorption of the incident light, while the classical absorption cross-section is obtained from the interaction Hamiltonian (\ref{coupling-prm-em}) in the lowest order. Obviously, if the incident light frequency $\omega$ is close to a plasmon mode frequency $\omega_{m}$, this interaction is enhanced due to resonant plasmon absorption.

To illustrate the role of optical dispersion in the above first-order effects, in Fig.~\ref{fig2} we compare the normalized absorption cross-section for an Au nanorod in water, calculated using Eq.~(\ref{cross-abs-mode}), with the result of Markov approximation obtained by setting $\omega=\omega_{m}$ in the coupling parameters  $\bm{\mu}_{m}(\omega)$ and $\gamma_{m}(\omega)$. With optical dispersion included, the plasmon resonance spectral shape deviates slightly from the Lorentzian  by exhibiting a red shift of the peak position along with overall enhancement of the lower frequency range. Although, in this case, the effect appear to be relatively small, much stronger non-Markovian effects are expected in higher orders, including in strongly-coupled QE-plasmon systems to be discussed elsewhere.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.95\columnwidth]{fig2}
\caption{\label{fig2} The calculated absorption cross-section of an Au nanorod in water, given by the full expression (\ref{cross-abs-mode}), is compared to the Markov approximation result obtained from the same expression but with  $\omega=\omega_{m}$ in the Au dielectric function. Inset: Schematics of the Au nanorod.
 }
\end{center}
\vspace{-6mm}
\end{figure}

The Hamiltonian $\hat{H}=\hat{H}_{\rm b}+\hat{H}_{\rm b-qe}+\hat{H}_{\rm b-em}+\hat{H}_{\rm qe}+\hat{H}_{\rm qe-em}$ provides a starting point for studying quantum correlations and non-Markovian dynamics in hybrid plasmonic systems involving localized plasmons interacting with QEs and the EM field. Within this framework, classical plasmons are encoded  in the coupling parameters (\ref{coupling-prm-em}) and (\ref{coupling-prm-qe}), respectively,  to mediate resonant coupling between the system components.  Namely, the classical  effects of  resonant plasmon excitation by  the EM field and resonance ET between the QEs and plasmons, which underpin most of the plasmon-enhanced spectroscopy phenomena,  now emerge in the lowest order of quantum perturbation theory. In higher orders, these classical effects will modulate quantum correlations and non-Markovian dynamics in hybrid plasmonic systems.

\section{Conclusions}

In summary, we developed a quantization approach for surface plasmons in metal-dielectric structures that  accounts for optical dispersion of the medium complex dielectric function and, hence, is suitable for describing non-Markovian effects in quantum plasmonics. Starting within the macroscopic quantization framework based on the dissipation-fluctuation theorem, we derived, in the Markov approximation, the canonical quantization scheme for interacting surface plasmons  and obtained the relevant coupling parameters in terms of local fields and system geometry.  Beyond Markov approximation, we developed a quantum approach in terms of discrete set of bosonic modes with linear dispersion whose interactions with  the electromagnetic field and quantum emitters are mediated by the classical plasmons, providing resonant coupling between the system components. We have shown that, within this approach, the classical plasmonic enhancement effects, such as resonant plasmon excitation by incident light and resonance energy transfer from a quantum emitter to  plasmonic structure, are obtained in the first-order of perturbation theory.  The quantum dynamics of bosonic modes is restricted to the reduced Hilbert space of reservoir states which overlaps with the electric fields of plasmon modes.

\acknowledgments
This work was supported in part by the National Science Foundation  Grants  Nos. DMR-2000170, DMR-1856515,  DMR-1826886 and HRD-1547754.


\section{Introduction}
\label{intro}
Understanding the nature of dark matter is a major goal of modern physics. There is a distinct possibility that dark matter is an as-yet undiscovered fundamental particle; there are several avenues being explored that target the direct detection of particle dark matter \cite{Klasen_2015}. Ultra-low background, underground direct detection searches are a prominent example. Of these, most experiments search for an excess of interaction events beyond known backgrounds in an active detector volume, that may be due to galactic dark matter.

To unambiguously attribute any measured excess to dark matter, an astrophysical signature will be required. For non-directional terrestrial experiments, one such signature is the annual modulation of the rate of dark-matter induced events in the detector, owing to the orbit of the earth around the sun \cite{Freese_2013}. However, any unexplained signal excess that modulates annually could also be due to a poorly understood seasonal systematic effect. Indeed, the DAMA experiment has observed a modulation which could be attributed to dark matter \cite{Bernabei2018}, but their result remains contested due to tension with other experiments and suggested seasonal systematics \cite{Davis_2014, mckinsey2018dama}. There are a number of NaI:Tl based detectors currently operating, including COSINE and ANAIS, that have recently reported their first annual modulation searches~\cite{Adhikari_2019_modulation, Amare_2019_modulation}; in both cases they are not yet sensitive enough to confirm or refute the DAMA result. In addition, both of these detectors are located in the Northern Hemisphere, whereas a measurement in both hemispheres can control for seasonal effects to verify the astrophysical origin of any modulation in a DAMA-like experiment. This is a major goal for the present SABRE experiment, as discussed in Section~\ref{sabre}.

For detectors with a directional capability, there is the possibility to measure an additional diurnal modulation due to the rotation of the Earth, which modulates the detector orientation with respect to the average incident dark matter velocity \cite{Ohare2015}. The magnitude of this modulation depends upon the latitude of the detector. A southern-hemisphere experiment with a latitude equal to Stawell, Victoria has a larger diurnal modulation than any existing underground laboratory~\cite{urquijo2016southern}. Since it does not possess a directional capability we do not expect to observe the diurnal modulation with the SABRE detector \cite{BernabeiConf2018}.

\section{The Stawell Underground Physics Laboratory and the Australian Research Council (ARC) Centre of Excellence for Dark Matter Particle Physics}\label{supl}
The Stawell Underground Physics Laboratory (SUPL) is a new, general-purpose, underground laboratory currently under construction, and is the first such laboratory in the Southern Hemisphere. SUPL is situated 1025 m underground in the Stawell Gold Mine in regional Victoria, Australia at a latitude and longitude of (-36.060, 142.801). The laboratory has similar overburden (2900 m water equivalent) and radioactive background to the Gran Sasso National Laboratory (LNGS). The laboratory design includes a radon suppression system and a clean area. Details on the SUPL design are given elsewhere \cite{urquijo2016southern}.

The creation of SUPL adds a major new capability to Australia's research infrastructure. Furthermore, it will be exploited through the ARC Centre of Excellence for Dark Matter Particle Physics that will support a broad experimental program to leverage SUPL and other Australian areas of strength, that are relevant to dark matter particle physics, in a coordinated way. The Centre will be led by the University of Melbourne, who will work with the Australian National University, University of Adelaide, University of Sydney, University of Western Australia and the Swinburne University of Technology, as well as a large number of Australian and international partner organisations. The Centre will prosecute research across four main themes: direct detection searches for dark matter with underground experiments at SUPL and above-ground experiments targeting axion-like particles; precision metrology to facilitate the ultra-pure, ultra-sensitive experimental measurement program; searches for dark matter production at the Large Hadron Collider; and a theory program focused on dark matter phenomenology. Details of the new Centre are given on its website~\cite{CoEDMPP}.

\section{SABRE}\label{sabre}
An important question facing experimental searches for dark matter is the origin of the signal observed by the DAMA experiment and claimed by them to be direct detection of dark matter. DAMA has measured an annual modulation in the rate of 2-6 keV events in their NaI:Tl detector for many years \cite{Bernabei2018}, which may be interpreted as being due to WIMP dark matter. Given the null results from other experiments (see summaries, for example, in Refs.~\cite{Klasen_2015,PDG2019}), DAMA's measured modulation is strongly inconsistent with most WIMP interaction models. However, the DAMA result has yet to be verified or refuted in a model-independent way with an identical target material.

Given the potential significance of the DAMA measurement, it is essential to verify the observed signal in a separate experiment.
SABRE is a NaI:Tl-based experiment that aims to conclusively determine whether the DAMA observation originates from dark matter or an experimental systematic effect.

SABRE will consist of two detectors, one located in Gran Sasso National Laboratory (LNGS), Italy and the other in the Stawell Underground Physics Laboratory (SUPL), Australia. Having detectors in the Northern and Southern Hemispheres permits SABRE to control for seasonal systematic effects and provides a powerful means to verify whether the observed modulation has a galactic origin.

Each SABRE detector will consist of 50 kg of ultrapure NaI:Tl contained within a liquid scintillator veto and low-radioactivity passive shielding (Fig.~\ref{fig:SSDrawing}). The full SABRE design and expected sensitivity is outlined in Refs.~\cite{Antonello_2019_MC, Antonello_2019_PoP}, and is expected to confirm or refute the DAMA result at 5$\sigma$ significance within 3 years of operation

\begin{figure}[!tb]
\centering
\includegraphics[width=5cm]{Vessel_Render7_feb18_cropped.eps}
\caption{A rendering of the SABRE South detector. The central NaI:Tl detectors are surrounded by a vessel filled with liquid scintillator, while the vessel itself will be surrounded by low-radioactivity passive shielding (not shown).}
\label{fig:SSDrawing}
\end{figure}

\subsection{Quenching Factor Measurements}\label{QF}
Many WIMP models predict that the most prevalent interaction with galactic dark matter is via coherent elastic nuclear scattering. Thus, it is essential to understand a dark matter detector's nuclear recoil energy scale to interpret a possible signal \cite{Collar2010}. For scintillation detectors, the nuclear recoil light yield is reduced with respect to an electron recoil of the same energy, due to losses to lattice excitations and differences in the electronic stopping power \cite{BIRKS1964}. This can be quantified via the quenching factor; the ratio of the light yields from a nuclear recoil compared to an electron recoil with an identical energy.

While the NaI:Tl quenching factor has been measured previously \cite{Collar2013, Xu2015, stiegler2017QF, joo2018quenching}, there is tension between these measurements that exceeds the experimental uncertainties. The origin of the difference in the reported values is unclear, as different methods have been used and possible causes of systematic effects
were not always reported. The DAMA experiment has assumed an energy-independent Na quenching factor of 0.3 for their analysis, which is also at odds with the most recent quenching factor measurements.

In preparation for investigating the quenching factor of the SABRE crystals, we have performed measurements using an existing NaI:Tl crystal originally produced in the USSR circa 1988. Its dimensions were 40 mm high by 40 mm diameter and the Tl doping concentration was not specified on the datasheet. Two high quantum efficiency 1"x1" photomultipliers (Hamamatsu H11934-200, peak Q.E. of 43\%) were coupled to the single quartz window for readout. The low-energy, electron-recoil calibration used $^{241}$Am, $^{137}$Cs, and $^{133}$Ba sources. The calibration was corrected for the known nonlinear NaI:Tl electron recoil response~\cite{ValentineAndRooney1998, Khodyuk_2010}.

Nuclear recoils were generated by elastic neutron scattering at the ANU 14UD Heavy Ion Accelerator. A pulsed proton beam (<2~ns on, 749~ns off) with energies of 3 MeV, 5.2 MeV, and 6 MeV was used to irradiate a 500~$\mu$g/cm$^{2}$ LiF target backed with sufficient Ta to stop the beam. Quasi-monoenergetic neutrons produced by the $^{7}$Li(p,n) reaction were collimated onto the NaI:Tl crystal using a polyethylene collimator. Six EJ-309 liquid scintillator detectors were located at 12$^{\circ}$, 22.5$^{\circ}$, 40$^{\circ}$, 67.5$^{\circ}$, 112.5$^{\circ}$, and 135$^{\circ}$ with respect to the beam axis. These were placed at various distances ranging from 15 to 150~cm away from the NaI:Tl detector, with the closest having an angular acceptance of 0.2 Sr. Waveforms in the NaI:Tl and liquid scintillator detectors were acquired using a PIXIE16 digitizer from XIA~\cite{XIA}, sampling at 500 MSPS with 12 bit resolution. The waveform capture was gated on coincidences between the NaI:Tl and any of the liquid scintillator detectors. The beam timing signal from the accelerator was also recorded to facilitate time-of-flight calculations.

The digitized waveforms were processed to estimate the pulse arrival time, charge, and (for the liquid scintillation detectors) a particle identification metric. Baseline noise was removed using a Wiener filter \cite{smith2003DSP} built using the averaged single photoelectron waveform as the signal template. The filtering greatly improved the pulse arrival time estimate, which was determined using a threshold. The pulse charge was calculated by integrating the unfiltered waveforms. However, to avoid integrating noise on the NaI:Tl waveforms in particular, the charge was evaluated for regions that exceeded a threshold, including 4~ns on either side of these regions. For the liquid scintillation detectors, the charge was evaluated in two fixed integration windows relative to the inferred pulse arrival time. The ratio of the charges in these windows was adopted as the particle identification metric. The window lengths were chosen to maximise a figure of merit given by $\frac{\mu_{n} - \mu_{\gamma}}{\sigma_{n} + \sigma_{\gamma}}$, where $\mu$ and $\sigma$ are the mean and standard deviation of peaks in the particle identification metric corresponding to the neutron- and gamma-induced events.

\begin{figure}
    \centering
    \includegraphics[width=8.1cm]{RecoilSpec_2.eps}
    \caption{Measured nuclear recoil spectra, after applying the analysis cuts, for detectors at various angles and with a 3 MeV proton beam.}
    \label{fig:RecoilSpec}
\end{figure}

Analysis cuts were applied so as to select nuclear recoil events in the NaI:Tl detector that are a result of elastic neutron scattering.  The times-of-flight between the pulsed beam arrival at the target and signals occurring in both the NaI:Tl and the liquid scintillator, were required to be consistent with a neutron, while the particle-identification metric measured by the liquid scintillator detector was gated to select neutron-like events.

Example recoil distributions after applying cuts are shown in Fig.~\ref{fig:RecoilSpec}. There are peaks in the spectra that can be attributed to nuclear recoils from Na and I and have the expected dependence upon recoil angle. The peak at $\sim$60~keV that is present in all spectra is from the inelastic excitation of $^{127}$I. The recoil peaks were fit by a Gaussian with a linear plus constant background. The incident neutron energy was calculated using PINO \cite{Reifarth_2009}. The recoil-energy calculation used the mean angle of each liquid scintillation detector with respect to the beam. A full Monte Carlo treatment of the expected distribution of nuclear recoil events in the NaI:Tl is underway, and will be reported in a future publication.

\begin{figure}
    \centering
    \includegraphics[width=8.1cm]{QFplot.eps}
    \caption{The quenching factors for Na and I recoils in NaI:Tl measured both in this work and in  previous measurements \cite{Xu2015, stiegler2017QF} taken with a similar methodology to that used here. The black and red dashed lines show the constant quenching factor values adopted by DAMA for Na and I, respectively.}
    \label{fig:QFs}
\end{figure}

The quenching factors determined from these measurements are given in Fig.~\ref{fig:QFs}, along with values from two other publications~\cite{Xu2015, stiegler2017QF}. The present results are in agreement with the high-energy extrapolation of both of these prior datasets, but do not yet distiniguish between them at the lower energies where the DAMA signal has been observed. The present 1980s detector has provided an effective proof-of-principle and it is anticipated that future measurements with a SABRE crystal (cut to a 1-inch cube to specifically match the high QE PMTs) will result in a low-energy threshold that probes quenching factors at the low energies required.

\subsection{Veto Scintillator Characterisation}
\label{sec:LS}
The SABRE veto functions both as high-purity active shielding and as a means of eliminating decay events from intrinsic radioimpurities in the inner detector volume that have correlated gamma-ray emissions. SABRE North will use a pseudocumene-based liquid scintillator from the Borexino experiment \cite{Benziger_2008}, which is situated adjacent to SABRE at LNGS, and has well-understood properties. Due to operational constraints associated with handling flammable liquids in a mine environment, SABRE South will not use pseudocumene, and will instead use a linear alkylbenzene (LAB) based scintillator. While the final fluorophore concentration is being optimised for the experiment, 2,6-diphenyl oxazole (PPO) will be used as the primary fluorophore, and 1,4-Bis~(2-methylstyryl)~benzene (bis-MSB) being evaluated as a secondary fluorophore. We do not expect that the use of different veto materials across SABRE North and South will significantly affect our sensitivity to annual modulations in two hemispheres.

To facilitate characterisation studies of LAB scintillator, we have purchased a bulk quantity (1 m$^{3}$) of industrial grade LAB (Jintung Petrochemical) and developed a simple purification process using vacuum distillation. Acceptable distillation results were achieved with a temperature of 120 $^\circ$C and a pressure of 0.05 mBar. The important material parameters from a particle detection point of view are the scintillation light yield and the optical attenuation length; measurements of these are presented below.

\begin{figure}
    \centering
    \includegraphics[width=8cm]{UVVISimpurities.eps}
    \caption{The optical absorbance of linear alkylbenzene at various levels of purification via vacuum distillation. ``Boiling Flask'' refers to the remnant material that remains after distillation.}
    \label{fig:UVVISimpurities}
\end{figure}

\subsubsection{Optical Absorbance}
Measurements of the optical absorbance are achieved using a dual-beam UV-VIS spectrophotometer (Shimadzu UV-2600), and are performed relative to deionised water. Several unknown optical impurities are visible in the absorbance spectrum of unpurified LAB, and measurements of samples taken from the distillate and boiling flask confirm that these are less volatile than the LAB (see Fig.~\ref{fig:UVVISimpurities}).

\subsubsection{Light Yield}
Light yield measurements are performed with the dedicated experimental setup shown in Fig.~\ref{fig:LYsetup} housed in a custom-made dark box. Liquid scintillator samples of constant volume in a standard vial are placed in a fixed sample holder above a high quantum efficiency (peak quantum efficiency 35\%) 1.5" diameter photomultiplier (Hamamatsu R9420-100). A 10 $\mu$Ci $^{137}$Cs source is used to irradiate the scintillator from a fixed location, producing predominantly Compton electrons that deposit energy in the scintillator volume. Photomultiplier signals are collected and amplified before collecting a spectrum with a multi-channel analyser (FastComTec). Since the distribution of energy deposits between samples should be identical, for a fixed scintillation light yield, the measured spectra should also be identical. This can be quantified by comparing spectra measured under the same conditions using an analysis inspired by Ref.~\cite{Bignell_2015} and illustrated in Fig.~\ref{fig:LYdemo}. One spectrum is taken as a fixed reference, while the other is scaled by factors to emulate the spectrum that would be collected were the light yield higher or lower. The scaling factor which minimises the $\chi^2$ between the two spectra, calculated over a pre-defined spectral range, is taken as the `best fit' relative light yield. Repeated measurements of a single scintillator sample confirm that we can reproducibly measure light yields to 1\%.

\begin{figure}[h]

\begin{subfigure}[b]{2cm}
\includegraphics[width=2cm]{LYholder.eps}
\caption{}
\label{fig:LYsetup}
\end{subfigure}
\begin{subfigure}[b]{6.4cm}
\includegraphics[width=6.4cm]{LYdemo.eps}
\caption{}
\label{fig:LYdemo}
\end{subfigure}
\caption{(a) Rendering of the light yield test setup that is contained in a dark box. The lid of the vial is shown, while the $^{137}$Cs source is placed in the holder (small slit) atop the mounting post. (b) A reference spectrum (blue) is compared with a test spectrum (orange) over a limited range with various scaling factors. The green trace represents the rescaled test spectrum with minimum $\chi^2$ that best reproduces the blue spectrum.}
\label{fig:LYsubfigs}
\end{figure}

\subsubsection{Compatibility Tests}
While a number of characterisation studies relating to the SABRE veto scintillator are ongoing, we report here on material compatibility tests that have been underway for over 100 days at the time of writing. The SABRE veto design includes numerous materials that will be in contact with the liquid scintillator over the multi-year operating life of the experiment. To ensure chemical compatibility, tests are underway using material samples submerged in a LAB-based test scintillator. The test scintillator was composed of purified LAB, 3 g/L PPO, and 15 mg/L bis-MSB (see Section~\ref{sec:LS}).

A summary of the tested materials is given in Table~\ref{tab:Compatibility}. One important question considered in these tests was whether the steel vessel would be fully compatible with the LAB scintillator. The steel vessel represents the largest surface area in contact with the LAB out of any material, so to ensure its compatibility we considered the possibility of coating with a chemically inert fluoropolymer. Many of the remaining materials in Table~\ref{tab:Compatibility} relate to the potted bases of the veto PMTs that will be required to allow the high voltage and detector signal cables to be used when the PMTs are immersed in LAB.

\begin{table}[!t]
\centering
\caption{The materials used in the veto scintillator compatibility tests. Abbreviations: polytetrafluoroethylene (PTFE), fluoroelastopolymer (FEP), polyolefin (PO), photomultiplier tube (PMT).}
\label{tab:Compatibility}
\begin{tabular}{lll}
\hline
Label & Material & Relevance  \\\hline
SS & Stainless steel & Veto vessel \\
SS-PTFE & PTFE-coated steel & Vessel coating \\
SS-FEP & FEP-coated steel & Vessel coating \\
Stud & Steel stud bolt & PMT mount point \\
Viton & Viton o-ring & Vessel flange seals \\
BlackCable & Coaxial cable & veto PMT \\
Epoxy & Epoxy sealant & veto PMT base \\
Potting & Potting compound & veto PMT base \\
FEPtube & FEP heat-shrink & veto PMT base \\
POtube & PO heat-shrink & veto PMT base \\
\hline
\end{tabular}
\end{table}

Two samples of each material were placed in contact with the scintillator inside glass bottles capped with a PTFE-lined lid and left in an air-conditioned room. Several control samples with no material in contact with the scintillator were also prepared. For the first 10 weeks of the compatibility tests, weekly measurements of the light yield relative to the control and the optical attenuation length relative to the control, were performed on one of these samples. The second `reserve' sample of each material was not measured in order to control for contamination during the measurement process. After 10 weeks, a comparison was made between the two samples, and the measurement cadence shifted to monthly.

Two of the materials (potting compound and polyolefin heat shrink) were observed to affect the light yield (see Fig.~\ref{fig:LYcompatibility}), with the reserve samples confirming the results within uncertainty. Only one material (polyolefin heat shrink) was observed to significantly affect the optical absorbance (see Fig.~\ref{fig:UVVIScompatibility}). The reserve polyolefin tube sample also exhibited degradation at a somewhat lower level, although it is not clear whether the excess in the measured sample was due to contamination or normal variation between the polyolefin samples.

\begin{figure}
    \centering
    \includegraphics[width=8cm]{LY_latest.eps}
    \caption{The light yield of all liquid scintillator samples placed in contact with materials to be submerged in the liquid scintillator veto, relative to a control with no material present.}
    \label{fig:LYcompatibility}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=8cm]{UVVIS_POtube.eps}
    \caption{The optical absorbance of the liquid scintillator containing submerged polyolefin heat-shrink tube, relative to a control with no material. The legend labels denote the measurement date.}
    \label{fig:UVVIScompatibility}
\end{figure}

The polyolefin tube will not be used in the potted PMT base design for SABRE. The potting compound will be used, but is sealed inside an acrylic box with the tested epoxy compound and thus should not be directly exposed to the scintillator. A key outcome from these tests was the confirmation that the veto vessel need not be coated to ensure compatibility; if any degradation scales as the surface area to volume ratio and scales linearly with time in contact with the material, then the lack of observable change (at the 1\% level) after 100 days in our test geometry, translates to a stability at the level of 0.6\% over the 3-year life of the experiment.

\section{Future Work}
We have given an overview of detector characterisation studies that have been done to date for SABRE at the ANU. A more complete analysis of the quenching factors with Monte Carlo calculations of the expected distribution of recoil energies will be included in a future publication.
The quenching factor measurements also provide a labelled dataset of nuclear recoils and electron recoils at low energy in NaI:Tl. Studies that use this dataset are underway to investigate particle identification approaches in NaI:Tl. Effective particle identification can help to further suppress backgrounds for SABRE.
The LAB compatibility tests are ongoing, and a new measurement campaign with further materials to be submerged in the veto liquid will commence shortly.

\begin{acknowledgement}
The authors thank A. Duffy and S. Krishnan for use of the NaI:Tl crystal from Swinburne University of Technology, and the ANU technical staff for their support preparing the beamline for the quenching factor measurements. I.M. and F.N acknowledge travel support from AINSE for their attendance at the ANU quenching factor measurements. This research was supported by the Australian Research Council, grant numbers DP170101675, LE160100080, LE170100162 and LE190100196.
\end{acknowledgement}