# Entanglement detection in coupled particle plasmons

Javier del Pino  Departamento de Física Teórica de la Materia Condensada and
Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
Madrid E-28049, Spain  Instituto de Física Fundamental, IFF-CSIC, Calle
Serrano 113b, Madrid E-28006, Spain  Johannes Feist  Departamento de Física
Teórica de la Materia Condensada and Condensed Matter Physics Center (IFIMAC),
Universidad Autónoma de Madrid, Madrid E-28049, Spain  F.J. García-Vidal
Departamento de Física Teórica de la Materia Condensada and Condensed Matter
Physics Center (IFIMAC), Universidad Autónoma de Madrid, Madrid E-28049, Spain
Juan Jose García-Ripoll  Instituto de Física Fundamental, IFF-CSIC, Calle
Serrano 113b, Madrid E-28006, Spain

###### Abstract

When in close contact, plasmonic resonances interact and become strongly
correlated. In this work we develop a quantum mechanical model for an array of
coupled particle plasmons. This model predicts that when the coupling strength
between plasmons approaches or surpasses the local dissipation, a sizable
amount of entanglement is stored in the collective modes of the array. We also
prove that entanglement manifests itself in far-field images of the plasmonic
modes, through the statistics of the quadratures of the field, in what
constitutes a novel family of entanglement witnesses. Finally, we estimate the
amount of entanglement, the coupling strength and the correlation properties
for a system that consists of two or more coupled nanospheres of silver,
showing evidence that our predictions could be tested using present-day state-
of-the-art technology.

###### pacs:

42.50.Dv, 73.20.Mf, 03.67.Mn

Surface plasmons are hybrid light-matter excitations confined at the interface
between a metal and a dielectric. Due to their small mode volume and strong
electromagnetic (EM) fields, surface plasmons interact very strongly with
quantum optical emitters Chang  _et al._ (2006); Dzsotjan  _et al._ (2010);
Andersen  _et al._ (2011); Gonzalez-Tudela  _et al._ (2011), such as quantum
dots Akimov  _et al._ (2007), NV-centers Kolesov  _et al._ (2009) or inorganic
Gomez  _et al._ (2010) and organic molecules Bellessa  _et al._ (2004);
Schwartz  _et al._ (2011). This, together with their broadband nature, small
size and their inherent quantum properties make them a promising platform for
future integrated quantum information technologies Tame  _et al._ (2013).
However, a very important problem lies in the characterization and control of
those quantum properties. So far, several experiments have demonstrated that
coupling photons in and out of plasmonic resonances preserves quantum features
such as single-photon excitations and anti-bunching Akimov  _et al._ (2007),
photon-photon entanglement Altewischer  _et al._ (2002), energy-time
entanglement Fasel  _et al._ (2005) and squeezing Huck  _et al._ (2009). In
this work we focus on the quantum properties of the surface plasmon themselves
and in particular in how many-body entanglement can be engineered using arrays
of coupled plasmonic modes.

In this Letter, we present a plasmonic setup that intrinsically exhibits many-
body entanglement and provide a recipe for characterizing it experimentally.
Our results build on a quantum mechanical model for a 1D or a 2D array of
coupled nanoparticles Maier  _et al._ (2002a, b); Weick  _et al._ (2013) that
includes the dipole-dipole interaction between particle plasmons, the losses
in each nanoparticle and the possibility of injecting energy via coherent or
incoherent light. Using this model we can not only study the transport of
excitations through the plasmonic band, but we also demonstrate the emergence
of stationary entanglement in the array at room temperature. Moreover, we
argue that this entanglement can be detected by measuring fluctuations in the
far-field from the light that is emitted from the plasmonic array.

We introduce three important theoretical ideas. The first one is a quantum
mechanical model for the nanoparticle array that consists of an array of
coupled oscillating dipoles with nearest-neighbor interaction and a local
dissipation that accounts for the losses. This model results in a master
equation for the density matrix associated with the plasmonic array. The
second important idea is that, under very general circumstances, this density
matrix will be Gaussian Weedbrook  _et al._ (2012) and all properties of the
array can be deduced from expectation values or “moments” of a finite set of
operators. In practice this implies a single set of exactly solvable ordinary
differential equations that fully describes the evolution of the quantum
surface plasmons. This technique allows us to make predictions not only on the
dynamics of the dipoles (i.e., absorption and transport of energy) but also
about their correlations and the resulting entanglement.

The final idea in this work is a formal study of the experimental observables
that can detect the presence of entanglement in the plasmonic array, the so-
called entanglement witnesses Gühne and Tóth (2009); Horodecki  _et al._
(1996); Lewenstein  _et al._ (2000); Terhal (2000); Bruß  _et al._ (2002). To
this end, we study the plasmonic band and compute the fluctuations of the EM
field in momentum space. We formally prove that the presence of squeezing in
the light with opposite momenta is a signature of entanglement. From an
experimental point of view, this implies that by refocusing the far-field
light emitted from the structure and studying its quantum fluctuations [cf.
Fig. 1], the amount of entanglement that is present in the plasmonic array can
be quantified. This general result is valid even when the Gaussian assumption
or our underlying quantum model breaks down.

Figure 1: An array of interacting nanoparticles gives rise to a set of coupled
plasmonic modes. The far field emission of these modes is collected by a lens.
By correlating the properties of the light at different points in the focal
plane, we get information about the multipartite entanglement.

We model our coupled particle plasmons as a set of $N$ oscillating dipoles
forming a linear 1D array, which interact through nearest-neighbor dipole
coupling and may be subject to external driving. The Hamiltonian reads
($\hbar=1$)

$H=\sum_{n=1}^{N}\frac{\omega}{2}(p_{n}^{2}+x_{n}^{2})+\sum_{\langle
n,m\rangle}gx_{n}x_{m}+\sum_{n=1}^{N}f_{n}(t)x_{n},$ (1)

where $f_{n}(t)$ is a driving force, $x_{n}$ is the dipole moment of the
particle plasmon and $p_{n}$ its associated canonical momentum. $g$ is the
coupling strength between neighboring sites, $\left\langle{n,m}\right\rangle$,
which are separated by a distance $\Lambda$.

We introduce local dissipation by means of a master equation to describe the
evolution of the quantum state or density matrix, $\rho$. This equation groups
all plasmonic losses in a single parameter, $\gamma$, and reads

$\displaystyle\partial_{t}\rho=-\frac{i}{\hbar}[H,\rho]+\sum_{n=1}^{N}\frac{\gamma}{2}(2a_{n}\rho
a_{n}^{\dagger}-a_{n}^{\dagger}a_{n}\rho-\rho a_{n}^{\dagger}a_{n}),$ (2)

where $a_{n}=\frac{1}{\sqrt{2}}(x_{n}+ip_{n})$ are the Fock operators that
diagonalize each individual harmonic oscillator.

Due to the quadratic nature of the problem, we can assume that the ground
state of the array is Gaussian Weedbrook  _et al._ (2012), as is usually done
in linear optics. This implies that the density matrix $\rho$ can be
reconstructed from the expectation values,
$\left\langle{O}\right\rangle:=\mathrm{tr}(O\rho)$, of the operators
$O\in\\{x_{n},p_{n},x_{n}x_{m},p_{n}p_{m},x_{n}p_{m}\\}$. Moreover, the
evolution equations for these “moments” form a closed set of first order
different equations that can be exactly solved, as described in detail in
section I of the Supplemental Material. Let us start with the first moment
equations, which describe the dynamics of the effective dipoles
$d_{n}=\left\langle{x_{n}}\right\rangle$. It is straightforward to find a set
of coupled driven classical harmonic oscillators subject to friction

$\ddot{d}_{n}=-\left(\omega^{2}+\frac{\gamma^{2}}{4}\right)d_{n}-2\omega
g\sum_{l}d_{l}-\gamma\dot{d}_{n}+f_{n},$ (3)

where the sum over $l$ is over nearest neighbors of $n$. This is a classical
model that has already been used to describe a particle plasmon array
Brongersma  _et al._ (2000) and shows the compatibility of our master equation
with earlier theoretical studies. In particular, our equations must describe
the transport of excitations and absorption of energy by the plasmonic array.
In fact, we can use the available experimental results to extract quantitative
information about the three parameters $g,\omega$ and $\gamma$, which
characterize our modeling.

Regarding transport, let us assume a coherent driving on the first site,
$f_{1}(t)\sim\sin(\nu t)$, and study the asymptotic state of the dipoles as a
function of the distance. From this calculation we can extract a propagation
length, $\xi$, defined as

$\xi=\frac{\sum_{n=1}^{N}n\Lambda\left|\left\langle{x_{n}}\right\rangle\right|}{\sum_{n=1}^{N}\left|\left\langle{x_{n}}\right\rangle\right|}.$
(4)

For the case of a very long chain, this propagation length would determine the
exponential decay of the plasmon population,
$\left\langle{x_{n}}\right\rangle\sim e^{-n\Lambda/\xi}$. In Fig. 2a we show
the propagation length in units of particle spacing, $\xi/\Lambda$, obtained
numerically for a chain of $N=20$ oscillators, as a function of the coupling
strength $g$ and plasmonic loss $\gamma$, under quasi-resonant driving
($\nu=0.99\omega$). Dissipation leads to a finite propagation length, which
grows with $g$ and diverges at the critical point $g/\omega=1/2,\gamma=0$,
where the current model becomes unphysical.

Figure 2: (a) Average propagation length (in units of $\Lambda$) in the 1D
chain of $N=20$ nanoparticles versus coupling strength, $g$, and local
dissipation, $\gamma$. (b) Entanglement in the chain measured by the
logarithmic negativity. (c) Entanglement witness in momentum space.

While the first order moments reproduce predictions of the classical theory,
the second order moments contain information about the non-classicality of the
many-body particle plasmon state. In particular, the matrix of second order
correlations, or covariance matrix, can also be exactly computed (see section
I of the Supplemental Material) and used to quantify the amount of
entanglement present in the plasmonic array. For this purpose let us eliminate
the driving $f_{n}(t)$, whose role is merely to displace the different
oscillator modes, without adding entanglement. In the absence of this driving,
we focus on the second order moment for the covariance matrix

$\sigma_{i,j}=\frac{1}{2}\left\\{\langle R_{i}R_{j}\rangle-\langle
R_{i}\rangle\langle R_{j}\rangle\right\\},$ (5)

where $\mathbf{R}^{T}=(x_{1},\ldots,x_{L},p_{1},\ldots,p_{L})$ is a vector
that groups all positions and momenta.

Let us consider a bipartition of the plasmonic array into two subarrays, A and
B. It is clear that the covariance matrix can be split into boxes that group
the operators of one or the other array,

$\sigma=\left(\begin{matrix}\sigma_{AA}&\sigma_{AB}\\\
\sigma_{BA}&\sigma_{BB}\end{matrix}\right),$ (6)

together with some off-diagonal terms, $\\{\sigma_{AB},\sigma_{BA}\\}$ that
imply some correlation (quantum or classical) between the two arrays. In order
to quantify purely quantum correlations, we compute the so called negativity
Weedbrook  _et al._ (2012), $E_{N}[\sigma;A,B]$. A value of
$E_{N}[\sigma;A,B]$ above zero means that the plasmonic array is entangled at
least with respect to this bipartition. Subsequent application of this
criterion to different partitions of the array can be used to ensure true
multipartite entanglement.

The results of this calculation are shown in Fig. 2b for a 1D array of 20
nanoparticles divided into two blocks of 10 consecutive particles. We plot the
negativity as a function of the coupling strength $g$ and the plasmonic loss
$\gamma$. As expected, entanglement grows with $g$ and becomes maximum at the
critical point $g=\omega/2,\gamma=0$, where the propagation length diverges.
The effect of dissipation is to decrease the entanglement, which remains
sizable for moderate coupling strengths, $g\simeq\gamma$.

Unfortunately, the negativity is not an observable. It may be estimated from
the full covariance matrix if a sufficiently accurate reconstruction of this
matrix is available, but this is an experimentally daunting task. It would
therefore be interesting to have an experimental criterion that allows the
detection of entanglement in the plasmonic chain with the least number of
measurements, while being robust to noise and imperfections.

For this task we suggest what is called an entanglement witness Gühne and Tóth
(2009); Horodecki  _et al._ (1996); Lewenstein  _et al._ (2000); Terhal
(2000); Bruß  _et al._ (2002). A witness is an observable $W$ such that when
its expectation value $\langle W\rangle=\mathrm{Tr}(W\rho)$ becomes negative,
we can positively assure that the state $\rho$ is not separable. There are
several such entanglement criteria in the literature of quantum optics. One of
them is the so-called Duan criterion for detecting two-mode squeezing Duan
_et al._ (2003), which was later extended by Hyllus and Eisert Hyllus and
Eisert (2006) to include multipartite entanglement. In this work we develop a
very general but simpler version of this last protocol.

Theorem: Let us take two vectors $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ which
satisfy the following conditions: (i) they are normalized,
$\|{\mathbf{u}_{i}}\|=1$, (ii) have the same modulus element-wise
($|u_{1,i}|=|u_{2,i}|$) and (iii) define two pairs of canonical variables,

$X_{k}=\sum_{j=1}^{L}u_{k,j}x_{j},\ \mathrm{and}\
P_{k}=\sum_{j=1}^{L}u_{k,j}p_{j}.$ (7)

If the two opposite quadratures are squeezed

$\langle\Delta{X}_{1}^{2}\rangle+\langle\Delta{P}_{2}^{2}\rangle<1,$ (8)

then the state is entangled. The demonstration of this theorem is presented in
the Supplemental Material, section II.

While the conditions (i)-(ii) might seem rather artificial, they can be
satisfied by the normal modes of the plasmonic array. The undriven part of
Hamiltonian (1) can be diagonalized using normal modes $\\{X_{k},P_{k}\\}$
(see details in section III of the Supplemental Material)

$H_{0}=\sum_{k}\frac{\omega}{2}(P_{k}^{2}+\lambda_{k}X_{k}^{2}),$ (9)

where $k$ represents the quantized momentum, $k=\pi j/[(N+1)\Lambda]$ with $j$
running from $1$ to $N$. The magnitude $\lambda_{k}=1+2(g/\omega)\cos
k\Lambda$ determines the plasmonic dispersion band,
$\omega_{k}=\omega\sqrt{\lambda_{k}}$.

Therefore, in the case of a 1D linear chain (corresponding to open boundary
conditions) and for a very large number of nanoparticles, $\mathbf{u}_{1}$ and
$\mathbf{u}_{2}$ of the theorem could be the wavefunctions associated to two
eigenmodes with opposite momenta $(k,k^{\prime})=(k,\pi/\Lambda-k)$, which are
equal in modulus and only differ in the fact that one has alternating signs
and the other does not, $u_{1,j}=(-1)^{j}u_{2,j}$. From a practical point of
view, this means that we can detect entanglement by looking for squeezing
among states with momenta $k$ and $(\pi/\Lambda-k)$. In other words, we can
define our entanglement witness

$W_{k}:=\mathrm{min}\\{0,\langle\Delta{X_{k}}^{2}\rangle+\langle\Delta{P}_{\pi/\Lambda-k}^{2}\rangle-1\\},$
(10)

so that $W_{k}<0$ implies entanglement. For the particular case $k=0$, i.e.,
the extrema of the dispersion band, we can find an analytical expression for
the entanglement witness (see details in section III of the Supplemental
Material)

$W_{0}=1+\frac{\frac{2g}{\omega}(\frac{2g}{\omega}-1)}{\frac{\gamma^{2}}{\omega^{2}}+4(1-\frac{2g}{\omega})}-\frac{\frac{2g}{\omega}}{\frac{\gamma^{2}}{\omega^{2}}+4(1+\frac{2g}{\omega})}.$
(11)

Fig. 2c presents the numerical results corresponding to $W_{0}$. As it shown
in the plot, the growth of the witness follows the same trend as that of the
negativity, hence providing the same amount of information.

Figure 3: Absorption versus frequency for a single silver nanosphere (red
line) and a dimer (blue line). In these calculations the radii of the
nanoparticles is set to $R=25$ nm whereas the separation between nanoparticles
in the dimer case is $2$ nm. The dashed grey line represents a Lorentzian fit
to the absorption spectrum of the single nanosphere that is used to estimate
$\gamma$.

In what follows we describe how this entanglement could be measured using
present-day state-of-the-art technology. Squeezing in the plasmonic band is
related to entanglement, and the same applies to far-field images of the
lattice. The light emitted by the plasmons maps the quadratures in the
collective variables $\\{X_{k},P_{k}\\}$ onto the equivalent variables of the
field propagating along directions $\pm k$. This light can be collected by a
large aperture lens, so that each value of the momentum is mapped to a
different point on the focal plane of the lens, as sketched in Fig. 1.
Selecting the photons with the appropriate momenta, we can perform homodyne
detection Weedbrook  _et al._ (2012); Welsch  _et al._ (1999) to measure the
quadratures and recover the value of $W_{k}$ mentioned above. Moreover, two
important features make this a very useful protocol. The first one is that our
choice of witness (i.e., momentum pairs) is not relevant, as we get similar
results for other values of the momentum. This is a signature that the state
is indeed many-body entangled. The second one is that while we have estimated
$W_{k}$ using Gaussian states, the entanglement witness is valid for any
physical state. In other words, measuring $W_{k}$ detects entanglement
irrespectively of the underlying physical model.

The proposed measurements could be realized using different types of coupled
plasmonic modes. One interesting possibility is provided by already existing
setups with gold or silver nanoparticles Maier  _et al._ (2002a, b). Earlier
experiments with such nanoparticles revealed short propagation lengths,
discouraging the use of such arrays for the transport of quantum information.
However, in Fig. 2 it can be appreciated that, while the plasmon propagation
length is related to the coupling strength and local loss, there can be a non-
zero amount of entanglement even when the surface plasmons do not propagate
efficiently. As an example and to provide a quantitative and realistic
estimation, we have calculated the EM coupling between two silver nanospheres
of radii $R=25~{}$nm and separated by a distance of $2~{}$nm. As shown in Fig.
3, we obtain a coupling strength of around $g/\omega\approx 0.15$. By looking
at the absorption spectrum for a single nanoparticle we can also extract a
value for the loss coefficient, $\gamma/\omega\approx 0.08$. These two values
for $g$ and $\gamma$ are fully compatible with earlier works studying larger
arrays Harris  _et al._ (2009). For this coupling and the associated plasmonic
loss, we expect a measurable amount of squeezing, $12\%$ (see Fig. 2c), which
would be a conclusive evidence of many-body entanglement within the plasmonic
array.

Summing up, in this work we have studied a quantum model for an array of
particle plasmons. The model, which can be extended to any system of
interacting plasmonic resonances, not only describes the collective resonances
and the transport of excitations through the system, but it also predicts the
existence of many-body entanglement in the system. Using the formalism of
Gaussian states and entanglement witnesses we have provided an experimental
protocol to detect this entanglement and estimated the strength of the
measurement outcomes for realistic setups. The entanglement witness developed
in this work is quite general, as it detects entanglement in far-field images
even for states that are not Gaussian, including coupled surface plasmons that
do not fall within our model. Moreover, some of these ideas can be exported to
other fields, such as nanophotonics, matter waves and the study of coupled
resonators in superconducting circuits.

This work has been funded by the European Research Council (ERC-2011-AdG
Proposal No. 290981). We also acknowledge financial support from EU FP7
project PROMISCE, CAM Research Consortium QUITEMAD (S2009-ESP-1594) and
Spanish MINECO projects FIS2012-33022 and MAT2011-28581-C02-01.

## References

  * Chang  _et al._ (2006) D. E. Chang, A. S. Sørensen, P. R. Hemmer, and M. D. Lukin, Phys. Rev. Lett. 97, 053002 (2006).
  * Dzsotjan  _et al._ (2010) D. Dzsotjan, A. S. Sørensen, and M. Fleischhauer, Phys. Rev. B 82, 075427 (2010).
  * Andersen  _et al._ (2011) M. L. Andersen, S. Stobbe, A. S. Sørensen, and P. Lodahl, Nat. Phys. 7, 215 (2011).
  * Gonzalez-Tudela  _et al._ (2011) A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Phys. Rev. Lett. 106, 020501 (2011).
  * Akimov  _et al._ (2007) A. V. Akimov, A. Mukherjee, C. L. Yu, D. E. Chang, A. S. Zibrov, P. R. Hemmer, H. Park, and M. D. Lukin, Nature 450, 402 (2007).
  * Kolesov  _et al._ (2009) R. Kolesov, B. Grotz, G. Balasubramanian, R. J. Stöhr, A. A. Nicolet, P. R. Hemmer, F. Jelezko, and J. Wrachtrup, Nat. Phys. 5, 470 (2009).
  * Gomez  _et al._ (2010) D. Gomez, K. Vernon, P. Mulvaney, and T. Davis, Nano Lett. 10, 274 (2010).
  * Bellessa  _et al._ (2004) J. Bellessa, C. Bonnand, J. C. Plenet, and J. Mugnier, Phys. Rev. Lett. 93, 036404 (2004).
  * Schwartz  _et al._ (2011) T. Schwartz, J. A. Hutchison, C. Genet, and T. W. Ebbesen, Phys. Rev. Lett. 106, 196405 (2011).
  * Tame  _et al._ (2013) M. Tame, K. McEnery, Ş. Özdemir, J. Lee, S. Maier, and M. Kim, Nat. Phys. 9, 329 (2013).
  * Altewischer  _et al._ (2002) E. Altewischer, M. Van Exter, and J. Woerdman, Nature 418, 304 (2002).
  * Fasel  _et al._ (2005) S. Fasel, F. Robin, E. Moreno, D. Erni, N. Gisin, and H. Zbinden, Phys. Rev. Lett. 94, 110501 (2005).
  * Huck  _et al._ (2009) A. Huck, S. Smolka, P. Lodahl, A. S. Sørensen, A. Boltasseva, J. Janousek, and U. L. Andersen, Phys. Rev. Lett. 102, 246802 (2009).
  * Maier  _et al._ (2002a) S. A. Maier, M. L. Brongersma, P. G. Kik, and H. A. Atwater, Phys. Rev. B 65, 193408 (2002a).
  * Maier  _et al._ (2002b) S. A. Maier, P. G. Kik, and H. A. Atwater, Appl. Phys. Lett. 81, 1714 (2002b).
  * Weick  _et al._ (2013) G. Weick, C. Woollacott, W. L. Barnes, O. Hess, and E. Mariani, Phys. Rev. Lett. 110, 106801 (2013).
  * Weedbrook  _et al._ (2012) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
  * Gühne and Tóth (2009) O. Gühne and G. Tóth, Phys. Rep. 474, 1 (2009).
  * Horodecki  _et al._ (1996) M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Lett. A 223, 1 (1996).
  * Lewenstein  _et al._ (2000) M. Lewenstein, B. Kraus, J. Cirac, and P. Horodecki, Phys. Rev. A 62, 052310 (2000).
  * Terhal (2000) B. M. Terhal, Phys. Lett. A 271, 319 (2000).
  * Bruß  _et al._ (2002) D. Bruß, J. I. Cirac, P. Horodecki, F. Hulpke, B. Kraus, M. Lewenstein, and A. Sanpera, J. Mod. Optics 49, 1399 (2002).
  * Brongersma  _et al._ (2000) M. L. Brongersma, J. W. Hartman, and H. A. Atwater, Phys. Rev. B 62, R16356 (2000).
  * Duan  _et al._ (2003) L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, in  _Quantum Information with Continuous Variables_ (Springer, 2003) pp. 145–153.
  * Hyllus and Eisert (2006) P. Hyllus and J. Eisert, New J. Phys. 8, 51 (2006).
  * Welsch  _et al._ (1999) D.-G. Welsch, W. Vogel, and T. Opatrný, in  _Progress in Optics_, Vol. 39, edited by E. Wolf (Elsevier, 1999) pp. 63 – 211.
  * Harris  _et al._ (2009) N. Harris, M. D. Arnold, M. G. Blaber, and M. J. Ford, J. Phys. Chem. C 113, 2784 (2009).

Supplemental Material

## I. Moment equations

We develop the general framework for studying the steady state of Hamiltonians
with quadratic bosonic operators. In the first place, it is convenient to
define $\boldsymbol{R}^{T}=\left(x_{1},\ldots,x_{N},p_{1},\ldots,p_{N}\right)$
and write $H$ as a quadratic form

$H=\frac{1}{2}\boldsymbol{R}^{T}\boldsymbol{B}\boldsymbol{R}+\boldsymbol{F}(t)^{T}\boldsymbol{R},$
(12)

where $\boldsymbol{B}$ is a real, symmetric matrix and
$\boldsymbol{F}(t)^{T}=\left(f_{1}(t),\ldots,f_{N}(t),0,\ldots,0\right)$
accounts for possible driving forces.

If the time evolution of a density matrix $\rho$ is given by

$\displaystyle\partial_{t}\rho=-i[H,\rho]+\sum_{n}\frac{\gamma}{2}(2a_{n}\rho
a_{n}^{\dagger}-a_{n}^{\dagger}a_{n}\rho-\rho a_{n}^{\dagger}a_{n}),$ (13)

where $a_{n}=\frac{1}{\sqrt{2}}(x_{n}+ip_{n})$, we can easily show that the
time evolution of the mean value of a time-independent operator $O$ is given
in a compact form as

$\partial_{t}\langle
O\rangle=-i\langle\left[O,H\right]\rangle+\sum_{n}\frac{\gamma}{2}\langle[a_{n}^{\dagger},O]a_{n}+a_{n}^{\dagger}\left[O,a_{n}\right]\rangle.$
(14)

Here, we used the cyclic invariance of the trace and
$\mathrm{Tr}\left\\{\dot{\rho}O\right\\}=\partial_{t}\langle O\rangle$. We
apply this idea to the first and second moments of the quadratures,
$O\in\\{x_{n},p_{n},x_{n}x_{m},p_{n}p_{m},x_{n}p_{m}\\}$. Writing Eq. 14 in
the quadrature basis, we get

$\partial_{t}\langle
O\rangle\\!=\\!-i\langle[O,H]\rangle\\!+\\!\sum_{nm}\frac{\Gamma_{nm}}{2}\langle[R_{m}^{\dagger},O]R_{n}\\!+\\!R_{m}^{\dagger}[O,R_{n}]\rangle,$
(15)

where
$\boldsymbol{\Gamma}=\bigoplus_{n=1}^{N}\frac{\gamma}{2}\left(\begin{array}[]{cc}1&-i\\\
i&1\end{array}\right)$ is a matrix that contains the effective dissipation
rates corresponding to operators $R_{n}$ and $R_{m}$ in this expression. The
$\bigoplus$ symbol denotes the direct sum of matrices, so for a set of
matrices $\\{A_{n}\\}$,
$\bigoplus_{n}A_{n}=\mathrm{diag}(A_{1},A_{2},\ldots,A_{n})$. If we now make
use of the commutation relations for quadrature operators, written in compact
form as

$[R_{n},R_{m}]=i\Omega_{nm}\quad\mathrm{with}\quad\boldsymbol{\Omega}=\left(\begin{array}[]{cc}0&\boldsymbol{1}_{N}\\\
-\boldsymbol{1}_{N}&0\end{array}\right),$ (16)

where $\boldsymbol{1}_{N}$ denotes the $N\mathrm{-dimensional}$ identity
matrix, it is straightforward to arrive to a closed set of $2N$ equations for
the first moments $\langle\boldsymbol{R}\rangle$, which we write in matrix
form as

$\partial_{t}\langle\boldsymbol{R}\rangle=\left(\boldsymbol{W}+\boldsymbol{\Omega
F}(t)\right)\langle\boldsymbol{R}\rangle$ (17)

where $\boldsymbol{W}=\boldsymbol{\Omega
B}+\frac{i}{2}(\boldsymbol{\Omega\Gamma}+(\boldsymbol{\Gamma\Omega})^{T})$.

As discussed in the main text, the second moments give information about the
nonclassical properties of the plasmonic array. In the same spirit, we can
write the following $2N\times 2N$ equations for the second moments
$\langle\boldsymbol{C}\rangle=\langle\boldsymbol{R}\boldsymbol{R}^{T}\rangle$

$\partial_{t}\langle\boldsymbol{C}\rangle=\boldsymbol{W}\langle\boldsymbol{C}\rangle+\langle\boldsymbol{C}\rangle\boldsymbol{W}^{T}-2\left(\boldsymbol{\Omega\Gamma\Omega}\right){}^{T}.$
(18)

Here, for simplicity, we have set $\boldsymbol{F}=0$ since its role is merely
to displace the first moments as we discuss in the main text and has no effect
on correlations.

Now, for instance, let us consider our particular case in which,
$\boldsymbol{B}=\boldsymbol{A}\oplus\boldsymbol{1}_{N}$, i.e., the Hamiltonian
may be written in the form

$H=\frac{\omega}{2}\mathbf{p}^{T}\mathbf{p}+\frac{\omega}{2}\mathbf{x}^{T}A\mathbf{x}+\mathbf{f}(t)^{T}\mathbf{x},$
(19)

where $\mathbf{x}^{T}=(x_{1},\ldots,x_{N})$,
$\mathbf{p}^{T}=(p_{1},\ldots,p_{N})$,
$\mathbf{f}(t)^{T}=(f_{1}(t),\ldots,f_{N}(t))$ and $\boldsymbol{A}$ is a
sparse matrix whose diagonal is unity, and the only other nonzero elements are
those connecting nearest neighbor sites, which are given by $2g/\omega$. In
this case we can write for the first moments

$\displaystyle\partial_{t}\langle x_{n}\rangle$ $\displaystyle=\omega\langle
p_{n}\rangle-\frac{\gamma}{2}\langle x_{n}\rangle,$ (20)
$\displaystyle\partial_{t}\langle p_{n}\rangle$ $\displaystyle=-\omega\langle
x_{n}\rangle-\frac{\gamma}{2}\langle p_{n}\rangle-2g\sum_{l}\langle
x_{l}\rangle+f_{n},$ (21)

where the sum over $l$ is over nearest neighbors of $n$. We then obtain the
equation for the effective dipole operator $d_{n}=\langle x_{n}\rangle$,

$\ddot{d}_{n}=-\left(\omega^{2}+\frac{\gamma^{2}}{4}\right)d_{n}-2\omega
g\sum_{l}d_{l}-\gamma\dot{d}_{n}+f_{n},$ (22)

which describe the dynamics of a set of coupled, driven harmonic oscillators
subject to friction. Additionally, the equations for the second moments read

$\displaystyle\partial_{t}\langle x_{n}x_{m}\rangle$
$\displaystyle=\omega\langle x_{n}p_{m}\\!+\\!p_{n}x_{m}\rangle-\gamma\langle
x_{n}x_{m}\rangle$ (23) $\displaystyle\partial_{t}\langle p_{n}p_{m}\rangle$
$\displaystyle=-\omega\langle x_{n}p_{m}+p_{n}x_{m}\rangle-\gamma\langle
p_{n}p_{m}\rangle-$ $\displaystyle-2g\sum_{l}(\langle
p_{n}x_{l}\rangle+\langle p_{m}x_{l}\rangle)$ (24)
$\displaystyle\partial_{t}\langle x_{n}p_{m}\rangle$
$\displaystyle=-\omega\langle x_{n}x_{m}\rangle+\omega\langle
p_{n}p_{m}\rangle-$ $\displaystyle-\frac{\gamma}{2}\langle
x_{n}p_{m}+p_{n}x_{m}\rangle-2g\sum_{l}\langle x_{n}x_{l}\rangle$ (25)
$\displaystyle\partial_{t}\langle x_{n}p_{m}\rangle$
$\displaystyle=\partial_{t}\langle p_{m}x_{n}\rangle.$ (26)

From their steady state solution in the absence of driving,
$\partial_{t}\langle\boldsymbol{R}\rangle=0=\partial_{t}\langle\boldsymbol{C}\rangle$,
we build the covariance matrix.

$\sigma_{i,j}=\frac{1}{2}\left\\{\langle R_{i}R_{j}\rangle-\langle
R_{i}\rangle\langle R_{j}\rangle\right\\}.$ (27)

Then, we consider a bipartition of the plasmonic array into two subarrays, $A$
and $B$, and compute the so-called logarithmic negativity between 2
bipartitions $E_{N}\left[\sigma;\,A,B\right]$. This quantity can be given in
terms of the absolute value of the eigenvalues of the matrix
$i\boldsymbol{\Omega\sigma}$, after performing a non-physical operation known
as partial transposition, as is discussed in detail in Weedbrook  _et al._
(2012). It can be shown that this non-physical operation is equivalent to
changing the sign of the $p_{i}$ components of one of the subsystems.

## II. Entanglement witness: Proof

Let us take two vectors $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ which satisfy
these conditions: (i) they are normalized, $\|{\mathbf{u}_{i}}\|=1$, (ii) have
the same modulus element-wise ($|u_{1,i}|=|u_{2,i}|$) and (iii) define two
pairs of canonical variables, $X_{j}=\sum_{i}u_{ji}x_{i},\ \mathrm{and}\
P_{j}=\sum_{i}u_{ji}p_{i}$ with $j=1,2$. If we now compute the fluctuations of
these operators assuming that the state is simply separable,
$\rho=\bigotimes_{i}\rho_{i}$, we have

$\displaystyle\langle{(\Delta{X_{1}})^{2}}\rangle+\langle{(\Delta{P_{2}})^{2}}\rangle=$
$\displaystyle\langle{X_{1}^{\dagger}X_{1}}\rangle-\langle{X^{\dagger}_{1}}\rangle\langle{X_{1}}\rangle+\langle{P_{2}^{\dagger}P_{2}}\rangle-\langle{P_{2}^{\dagger}}\rangle\langle{P_{2}}\rangle$
(28) $\displaystyle=$ $\displaystyle\sum_{i,j}\left[u_{i,1}^{*}u_{j,1}\langle
x_{i}x_{j}\rangle+u_{i,2}^{*}u_{j,2}\langle
p_{i}p_{j}\rangle\right]-\sum_{i,j}\left[u_{i,1}^{*}u_{j,1}\langle
x_{i}\rangle\langle{x_{j}}\rangle+u_{i,2}^{*}u_{j,2}\langle
p_{i}\rangle\langle{p_{j}}\rangle\right]$
$\displaystyle\stackrel{{\scriptstyle(1)}}{{=}}$
$\displaystyle\sum_{i}\left[|u_{i,1}|^{2}\langle\Delta{x}_{i}^{2}\rangle+|u_{i2}|^{2}\langle\Delta{p}_{i}^{2}\rangle\right]\stackrel{{\scriptstyle(2)}}{{=}}\sum_{i}|u_{i,1}|^{2}\left[\langle\Delta{x}_{i}^{2}\rangle+\langle\Delta{p}_{i}^{2}\rangle\right]\stackrel{{\scriptstyle(3)}}{{>}}=\sum_{i}|u_{i,1}|^{2}=1.$

Here we have used various key ideas: In $(1)$ we use the fact that the state
is separable and thus
$\langle{x_{i}x_{j}}\rangle=\langle{x_{i}}\rangle\langle{x_{j}}\rangle$
whenever $i\neq j$. In $(2)$ we use the fact that both vectors have the same
modulus element-wise , $|u_{i1}|=|u_{i2}|$. Finally, in $(3)$ we use the fact
that $\langle\Delta{A}^{2}\rangle+\langle\Delta{B}^{2}\rangle>=\|[A,B]\|$ and
the normalization of the vectors.

This proof can be extended to treat fully all possible cases of separable
states, $\rho=\sum p_{i}\otimes_{j}\rho_{j}^{i}$, which are convex linear
combinations of the previous situation we have shown. In this case the only
difference is that there appear additional cross-terms due to the linear
combinations, but these terms can be shown to be larger than zero, thus
increasing the fluctuations Duan  _et al._ (2003).

## III. Normal modes

First, in this section, we are going to show how to diagonalize the effective
model in the dissipation-free case. The diagonalization of the matrix
$\boldsymbol{A}$ through an orthogonal transformation,
$\boldsymbol{A}=\boldsymbol{U}^{T}\boldsymbol{D}\boldsymbol{U}$, allows us to
define new canonical variables

$X_{k}=\sum_{i}u_{k,i}x_{i},\;P_{k}=\sum_{i}u_{k,i}p_{i}.$ (29)

In these new quadratures, the Hamiltonian becomes

$H=\sum_{k}\frac{1}{2}\omega\left(P_{k}^{2}+\lambda_{k}X_{k}^{2}\right)+\sum_{i}u_{k,i}f_{i}(t),$
(30)

where $\lambda_{k}=D_{k,k}$ are the eigenvalues of the quadratic form and we
introduce the effective drivings in momentum space,
$\tilde{f}_{k}=\sum_{i}u_{k,i}f_{i}$. Note that in absence of driving, the
normal frequencies of the problem will be

$\omega_{k}=\omega\lambda_{k},$ (31)

and the new Fock operators will be related to the original ones by a
complicated squeezing transformation

$\tilde{a}_{k}=\sum_{i}\left(\lambda_{k}^{1/2}u_{k,i}x_{i}+i\lambda_{k}^{-1/2}u_{k,i}p_{i}\right),$
(32)

that is the source of the entanglement of this problem.

As a particular instance of the lattice of coupled plasmons we will consider
the case of a one-dimensional lattice of regularly spaced nanoparticles, with
period $\Lambda$, that corresponds to the 1D open-boundary condition case of
$\boldsymbol{A}$. This tridiagonal matrix is diagonalized with the orthogonal
transformation

$u_{i,k}=u_{k,i}=\sqrt{\frac{2}{N+1}}\sin(k_{j}\Lambda i),\ \ \;i,j=1\ldots N$
(33)

where $N$ is the total lattice size and the quasimomenta $k_{j}=\pi
j/(N+1)\Lambda$ determine the eigenfrequencies

$\lambda_{k}=1+2(g/\omega)\cos(k_{j}\Lambda).$ (34)

Notice how for small momenta, when we reach the critical value $g=\omega/2$,
we recover a linear dispersion relation of photon-like quasiparticles with
diverging correlations. In practice, however, $g$ is below this limit and we
obtain a band of massive excitations with a finite correlation length.

Now, as the form of (15) is preserved under the transformation (33) we apply
the ideas exposed previously to the first and second moments of the canonical
variables, obtaining

$\displaystyle\partial_{t}\langle X_{k}\rangle$ $\displaystyle=\omega\langle
P_{k}\rangle-\frac{\gamma}{2}\langle X_{k}\rangle$ (35)
$\displaystyle\partial_{t}\langle P_{k}\rangle$
$\displaystyle=-\omega\lambda_{k}\langle X_{k}\rangle-\frac{\gamma}{2}\langle
P_{k}\rangle-\tilde{f}_{k}(t).$

For the second moments, taking $f=0$ to simplify the expressions, we get the
following closed set

$\displaystyle\partial_{t}\langle X_{k}^{2}\rangle$
$\displaystyle\\!=\omega\langle X_{k}P_{k}+P_{k}X_{k}\rangle-\gamma\langle
X_{k}^{2}\rangle+\frac{\gamma}{2}$ (36) $\displaystyle\partial_{t}\langle
P_{k}^{2}\rangle$ $\displaystyle\\!=\\!-\omega\lambda_{k}\langle
X_{k}P_{k}\\!+\\!P_{k}X_{k}\rangle-\gamma\langle
P_{k}^{2}\rangle+\frac{\gamma}{2}$ (37) $\displaystyle\partial_{t}\langle
X_{k}P_{k}\rangle$ $\displaystyle\\!=\\!-\omega\lambda_{k}\langle
X_{k}^{2}\rangle+\omega\langle
P_{k}^{2}\rangle\\!-\\!\frac{\gamma}{2}\langle\\!X_{k}P_{k}\\!+\\!P_{k}X_{k}\rangle$
(38) $\displaystyle\partial_{t}\langle P_{k}X_{k}\rangle$
$\displaystyle\\!=\partial_{t}\langle X_{k}P_{k}\rangle.$ (39)

In this case, its steady-state solution allows us to compute the fluctuations
needed when computing the proposed entanglement witness

$W_{k}:=\mathrm{min}\\{0,\langle\Delta{X_{k}}^{2}\rangle+\langle\Delta{P}_{(\pi/\Lambda)-k}^{2}\rangle-1\\}$
(40)

given $(\Delta O)^{2}=\langle O^{\dagger}O\rangle-\langle
O^{\dagger}\rangle\langle O\rangle$. In the case of infinite chain length
$N\gg 1$, if we compute this quantity between the extrema of the band we
arrive to

$W_{0}=1+\frac{\frac{2g}{\omega}(\frac{2g}{\omega}-1)}{\frac{\gamma^{2}}{\omega^{2}}+4(1-\frac{2g}{\omega})}-\frac{\frac{2g}{\omega}}{\frac{\gamma^{2}}{\omega^{2}}+4(1+\frac{2g}{\omega})}.$
(41)