We distinguish two cases and assume first that $N_0 < N^{5/6 + \delta}$ for some $\delta>0$. In this case applications of (<ref>) and Lemma <ref> show that the first term in (<ref>) is bounded by a constant times $L^{-2} N^{1/2+\delta}$.  Inspection of (<ref>) (recall that $\widetilde{N}_0 = \int |z|^2 \zeta(z) \de z$) shows $|\widetilde \mu | \lesssim 1/(\beta \widetilde{N}_0) + \widetilde{N}_0/(L^2N)$. We use this estimate, $N_0 \geq N^{2/3}$ (this implies, by Lemma <ref>, $\widetilde N_0 \gtrsim N^{2/3}$ ), and again Lemma <ref> to bound the second term in (<ref>) by a constant times $L^{-2} [ N^{\delta_{\mathrm{H}}} + N^{-1/6 + \delta_{\mathrm{H}} + \delta} + N^{1/3 + 2 \delta} ]$. If $N_0 \geq N^{5/6 + \delta}$ we apply part (a) of Lemma <ref> to bound the second term in (<ref>) from above by

\begin{align}
8 \pi \mathfrak{a}_N L^{-3} \widetilde{N}_0 (\widetilde{N}_0 - N_0) +  C \exp(-c N^{\delta}) \leq -8 \pi \mathfrak{a}_N L^{-3}\widetilde{N}_0 \sum_{p \in P_{\mathrm{B}}} \left( \gamma(p) - \gamma_0(p) \right) + C N^{\delta_{\mathrm{H}}}. \label{eq:Sec65}
\end{align}

To obtain the second bound, we also used $\sum_{p \in P_{\mathrm{L}}^{\mathrm{c}}} \gamma(p) \lesssim \exp(-c N^{2\delta_{\mathrm{L}}})$ for some $c>0$. We highlight that the first term on the r.h.s. of (<ref>) cancels the first term in (<ref>). We collect the above considerations, make the assumption $0 < \delta < 1/6$, and find

\begin{equation}
\frac{8\pi\frak{a}}{NL^3}\biggl[ \widetilde{N}_0 \sum_{p \in P_{\mathrm{B}}} \left( \gamma(p) - \gamma_0(p) \right) \biggr] + \widetilde{\mu}(\widetilde{N}_0 - N_0) \lesssim L^{-2} [ N^{\delta_\mathrm{H}} + N^{1/2 + \delta} ].
\label{eq:Sec66}
\end{equation}

It remains to collect our results.

In combination, (<ref>)–(<ref>), (<ref>) and (<ref>) imply the final upper bound

\begin{align}
\tr[\mathcal{H}_N \Gamma] - \frac{1}{\beta} S(\Gamma) &\leq \frac{1}{\beta} \sum_{p \in \Lambda_+} \ln\left( 1 - \exp(- \beta(p^2 - \mu_0)) \right) + \mu_0 (N - N_0) + 8 \pi \mathfrak{a}_N L^3 (\varrho^2 - \varrho_0^2) + F^{\mathrm{BEC}}(\beta,N_0,L,\mathfrak{a}_N) \nonumber \\
&- \frac{1}{2 \beta} \sum_{p \in \Lambda^*_+} \left[ \frac{16 \pi \mathfrak{a}_N \varrho_0(\beta,N,L)}{p^2} - \ln\left( 1 + \frac{16 \pi \mathfrak{a}_N \varrho_0(\beta,N,L)}{p^2} \right) \right] \nonumber \\
&+ C L^{-2} [ N^{1- \delta_{\mathrm{H}}} + N^{\delta_{\mathrm{H}} + 2 \delta_{\mathrm{B}}} + N^{-1/3 + \delta_{\mathrm{H}} + 2 \delta_{\mathrm{L}}} + N^{1/3 + \delta_{\mathrm{B}}} + N^{1/2 + \delta} + N^{2/3 - \delta_{\mathrm{B}}} ]. \label{eq:Sec67}
\end{align}

The parameters $\delta_{\mathrm{L}}, \delta$ need to be strictly positive but can otherwise be chosen as small as we wish. The requirements $\delta_{\mathrm{L}} \leq 1/6, \delta \leq 1/12$ assure that they play no role in the optimization. The optimal choice $\delta_{\mathrm{H}} = 1/2 - \delta_{\mathrm{B}}$ with error $N^{1/2+\delta_{\mathrm{B}}}$ follows by combining the first and the second term. Moreover, the optimal choice $\delta_{\mathrm{B}} = 1/12$ results if we combine $N^{1/2+\delta_{\mathrm{B}}}$ and the last term in the last line of (<ref>). This leads to an overall error term that is bounded by a constant times $L^{-2} N^{7/12} \ll L^{-2} N^{2/3}$. We recall that the above bound holds under the assumption $N_0 \geq N^{2/3}$.

We now prove a second bound with another trial state (see also Remark <ref>.(g)) that holds without a restriction on $N_0$ (it, however, captures the correct behavior of the free energy only if $N_0 \ll N^{5/6}$). As undressed trial state we choose the Gibbs state

\begin{equation}
G_0 = \frac{\exp(-\beta( \de \Gamma(-\Delta - \mu_0) ))}{\tr_{\mathscr{F}}[\exp(-\beta( \de \Gamma(-\Delta - \mu_0) ))]}.
\label{eq:Sec68}
\end{equation}

We define the dressed trial state $\widetilde{\Gamma}$ as in (<ref>) with $\Gamma_0$ replaced by $G_0$. To obtain an upper bound for the free energy of $\widetilde{\Gamma}$ we can use a simpler version of the above proof. This is related to the following facts: (a) A coherent state in the definition of our trial state is not needed and the pairing function of $G_0$ equals zero. (b) The eigenfunctions of $G_0$ are also eigenfunctions of $\de \Gamma(-\Delta)$. Accordingly, the special treatment of momentum modes in $P_{\mathrm{B}}$ at several places in the proof is not needed. (c) Since $[G_0,\mathcal{N}] = 0$, we have $\tr[\mathcal{N} \widetilde{\Gamma}] = \tr[\mathcal{N} G_0]$. We therefore simply state the result and leave further details to the reader:

\begin{align}
\tr[\mathcal{H}_N \widetilde{\Gamma}] - \frac{1}{\beta} S(\widetilde{\Gamma}) \leq F_0(\beta,N,L) + 8 \pi \mathfrak{a}_N L^{3} \varrho^2 + C L^{-2} N^{1/2}
\label{eq:Sec69}
\end{align}

with $F_0$ defined above (<ref>).

We are now prepared to provide the missing proofs in Section <ref>. Theorem <ref> follows from (<ref>), (<ref>), Proposition <ref> and fact that the absolute value of the term in the second line of (<ref>) is bounded by a constant times $N_0^2/(\beta N^2)$. The proof of Proposition <ref> is provided in Appendix <ref>, see Proposition <ref>. Finally, Corollary <ref> is a direct consequence of (<ref>), (<ref>) and Proposition <ref>.

A. D. gratefully acknowledges funding from the Swiss National Science Foundation (SNSF) through the Ambizione grant PZ00P2 185851. It is our pleasure to thank Marco Caporaletti, Phan Thành Nam, Marcin Napiórkowski and Robert Seiringer for inspiring discussions.

— Appendix —

§ THE SCATTERING EQUATION

In this appendix we collect some known properties of the finite volume scattering equation (<ref>). It is convenient to define $f(Nx)=f_N(x)$, where $f$ satisfies the eigenvalue equation

\begin{equation} \label{eq:NP1}
\bigg[-\Delta + \frac{v}{2}\bigg] f = \lambda_{\ell} f
\end{equation}

on the ball $|x| \leq N\ell$ with Neumann boundary conditions. It is normalized such that $f(x)=1$ holds for $|x|=N\ell$. By scaling, we have $N^2\lambda_\ell=\lambda_N$. In the next Lemma we collect the properties of $f_N,f$ and $\lambda_\ell$ that are useful for our analysis. The proof can be found in <cit.>.

Let $v \in L^3(\mathbb{R}^3)$ be nonnegative, compactly supported and spherically symmetric. Fix $0<\ell <L/2 $ and let $f$ denote the solution to (<ref>) and $f_N$ the solution to (<ref>). For $N \in \mathbb{N}$ large enough the following properties hold true.

* We have

\begin{equation}\label{eq:el}
\lambda_{\ell} = \frac{3\frak{a}}{(N\ell)^3}(1+\mathcal{O}\big(\frak{a}/\ell N)\big).
\end{equation}

* We have $0\leq f_{\ell}\leq 1$. Moreover there exists a constant $C > 0$ such that

\begin{equation} \label{eq:sl}
\left\lvert \int v(x)f(x)dx - 8\pi \frak{a} \right\lvert \leq \frac{C\frak{a}^2}{N\ell}.
\end{equation}

* There exists a constant $C > 0$ such that, for all $x \in \mathbb{R}^3$,

\begin{equation}\label{eq:Dw}
1-f(x) \leq \frac{C}{1 + |x|} \quad \text{and} \quad \left|\nabla f(x)\right| \leq \frac{C}{1 + x^2}.
\end{equation}

* There exists a constant $C > 0$ such that, for all $p \in \L^*_+$,

\begin{equation}\label{eq:wp}
\big|\widehat{(1-f_N)}(p)\big| \leq \frac{C}{N p^2}.
\end{equation}

By definition (<ref>), the function $\eta_p=-\widehat{(1-f_N)}(p)$ solves the equation

\begin{equation} \label{npscatt}
p^2\eta_p + \frac{\hat v_N(p)}{2} + \frac{1}{2L^3} \sum_{q \in \Lambda^*} \hat{v}_N(p-q)\eta_q =  \frac{\lambda_{N}}{L^3} \sum_{q \in \Lambda^*}\hat{f}_N(p-q)\hat{\mathds{1}}_{|x|\leq\ell}(q),
\end{equation}

where $\hat{\mathds{1}}_{|x|\leq\ell}(q)$ is the Fourier coefficient of the characteristic function of the ball with radius $\ell$. Note that we have reinstated units in (<ref>). Moreover, by (<ref>) and (<ref>) we have

\begin{equation} \label{eq:etabound}
|\eta_p| \lesssim \frac{L}{Np^2}. %\quad \text{and} \quad |\eta_0| \leq C\ell^2N^{-1}.
\end{equation}

Inequality (<ref>) implies

\begin{equation} \label{eq:sumPHC}
\sum_{p\in\L^*_+:|p|\geq \frac{1}{2L} N^{1-\delta_\mathrm{H}}}|\eta_p|^2 \lesssim L^6 N^{-3+\delta_\mathrm{H}}
\end{equation}

as well as

\begin{equation} \label{eq:sumPHCC}
\sum_{p\in P_\mathrm{H}^{\mathrm{c}}}|\eta_p| \lesssim L^3 N^{-\delta_\mathrm{H}},
\end{equation}

where $P_\mathrm{H}^{\mathrm{c}}$ denotes the complement of $P_\mathrm{H}$ in (<ref>).

§ BOGOLIUBOV FREE ENERGY

The goal of this section is to prove the following lemma.

We consider the limit $N \to \infty$, $\beta = \kappa \beta_{\mathrm{c}}$ with $\kappa \in (0,\infty)$ and $\beta_{\mathrm{c}}$ in (<ref>). Recall  definition (<ref>) for $\tilde{\epsilon}(p)$. There exists a constant $C>0$ such that

\begin{align}
\frac{1}{\beta} \sum_{p \in P_{\mathrm{B}}} \ln \left( 1- \exp(-\beta \tilde{\epsilon}(p)) \right) \leq& \frac{1}{\beta} \sum_{p \in P_{\mathrm{B}}} \ln \left( 1- \exp(-\beta (p^2-\mu_0)) \right) + 8 \pi \mathfrak{a}_N \varrho_0 \sum_{p \in P_{\mathrm{B}}} \frac{1}{\exp(\beta(p^2-\mu_0))-1} \label{eq:AppB1} \\
&- \frac{1}{2 \beta} \sum_{p \in \Lambda_+^*} \left[ \frac{16 \pi \mathfrak{a}_N \varrho_0}{p^2} - \ln\left( 1 + \frac{16 \pi \mathfrak{a}_N \varrho_0}{p^2} \right) \right] + \frac{CN_0^2}{N^2} \left[ \frac{N^{\delta_{\mathrm{B}}}}{L^2} + \frac{1}{\beta N^{\delta_{\mathrm{B}}}} + \frac{L^2}{\beta^2 N_0} \right]. \nonumber
\end{align}

We first assume $\mu_0 = 0$ and then comment on how to adjust the proof to $\mu_0 < 0$. Let us define the function

\begin{equation}
F(\alpha) = \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \ln\left(1-\exp\left(-|p| \sqrt{p^2 + \alpha}\right)\right).
\label{eq:AppB2}
\end{equation}

For $\alpha = 16 \pi \mathfrak{a}_N \varrho_0 \beta$ it equals $\beta$ times the l.h.s. of (<ref>). In the following we derive an asymptotic expansion of $F$ for small values of $\alpha$. We also define the functions

\begin{equation}
g_1(x) = \frac{1}{\exp(x)-1} \quad \text{ and } \quad g_2(x) = - \frac{1}{4 \sinh^2(x/2)}
\label{eq:AppB3}
\end{equation}

and note that the bounds

\begin{equation}
g_1(x) \geq \frac{1}{x} - C \quad \text{ and } \quad g_2(x) \leq \frac{-1}{x^2} + \frac{C}{x}
\label{eq:AppB3b}
\end{equation}

hold for $0 < x \leq 1$. The first and the second derivative of $F$ can be written in terms of $g_1$ and $g_2$ as

\begin{align}
F'(\alpha) &= \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} g_1\left( |p| \sqrt{p^2 + \alpha} \right) \frac{|p|}{2 \sqrt{p^2 + \alpha}}, \nonumber \\
F''(\alpha) &= \frac{1}{4} \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \left[ g_2\left( |p| \sqrt{p^2 + \alpha} \right) \frac{p^2}{p^2 + \alpha} - g_1\left( |p| \sqrt{p^2 + \alpha} \right) \frac{|p|}{(p^2 + \alpha)^{3/2}} \right],
\label{eq:AppB4}
\end{align}

and hence

\begin{equation}
F(\alpha) - F(0) - F'(0) \alpha = \frac{1}{4} \int_0^{\alpha} \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \left[ g_2\left( |p| \sqrt{p^2 + t} \right) \frac{p^2}{p^2 + t} - g_1\left( |p| \sqrt{p^2 + t} \right) \frac{|p|}{(p^2 + t)^{3/2}} \right] (\alpha-t) \de t.
\label{eq:AppB5}
\end{equation}

It remains to investigate the r.h.s. of this above identity.

Using the bounds in (<ref>), we see that it is bounded from above by

\begin{equation}
-\frac{1}{2} \int_0^{\alpha} \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \frac{(\alpha-t)}{(p^2 + t)^2} \de t + C \alpha \int_0^{\alpha} \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \frac{|p|}{(p^2+t)^{3/2}} \de t.
\label{eq:AppB6}
\end{equation}

Here, the integral in the second term is bounded by

\begin{equation}
\sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \int_0^{\alpha} \frac{1}{p^2+t} \de t = \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \ln\left( 1 + \frac{\alpha}{p^2} \right) \leq \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \frac{\alpha}{p^2} \lesssim \frac{L^2 \alpha N^{\delta_{\mathrm{B}}}}{\beta}.
\label{eq:AppB7}
\end{equation}

A straightforward computation also shows

\begin{equation}
\int_0^{\alpha} \frac{(\alpha-t)}{(p^2 + t)^2} \de t = \frac{\alpha}{p^2} - \ln\left( 1 + \frac{\alpha}{p^2} \right).
\label{eq:AppB8}
\end{equation}

In combination, (<ref>)–(<ref>) imply

\begin{equation}
F(\alpha) - F(0) - F'(0) \alpha \leq - \frac{1}{2} \sum_{p \in \beta^{1/2} P_{\mathrm{B}}} \left[  \frac{\alpha}{p^2} - \ln\left( 1 + \frac{\alpha}{p^2} \right) \right] + \frac{C L^2 \alpha^2 N^{\delta_{\mathrm{B}}}}{\beta}.
\label{eq:AppB9}
\end{equation}

Finally, using $\ln(1+x) \geq x - x^2/2$ for $x \geq 0$ we see that

\begin{equation}
\sum_{p \in \beta^{1/2} P^{\mathrm{c}}_{\mathrm{B}}} \left[  \frac{\alpha}{p^2} - \ln\left( 1 + \frac{\alpha}{p^2} \right) \right] \leq \frac{\alpha^2}{2 \beta^2} \sum_{p \in P^{\mathrm{c}}_{\mathrm{B}}} \frac{1}{p^4} \lesssim \frac{\alpha^2 L^4}{\beta^2 N^{\delta_{\mathrm{B}}}}.
\label{eq:AppB10}
\end{equation}

When we put our findings together, we obtain a proof of (<ref>) if $\mu_0=0$ (the last error term excluded).

If $\mu_0 < 0$ our proof applies without changes and we obtain the first term in the second line of (<ref>) with $p^2$ replaced by $p^2-\mu_0$. It is not difficult to check that the difference between these two terms is bounded by a constant times $N_0 L^2/(\beta^2 N^2)$, which proves the claim of the lemma.

§ PROPERTIES OF THE FREE ENERGY OF THE CONDENSATE

In this appendix we prove several statements concerning the effective condensate free energy in (<ref>), one of which is Proposition <ref>. The other statements are needed for the proof of Theorem <ref>. We start our discussion with a lemma that provides us with the asymptotic behavior of the chemical potential.

We consider the limit $N \to \infty$, $\beta = \kappa \beta_{\mathrm{c}}$ with $\kappa \in (0,\infty)$ and $\beta_{\mathrm{c}}$ in (<ref>). Let $g$ be the Gibbs distribution in (<ref>) and assume that $\int_{\mathbb{C}} |z|^2 g(z) \de z = M$. The chemical potential $\mu$ related to $g$ satisfies the following statements for a given $\epsilon > 0$:

* If $M \gtrsim N^{5/6 + \epsilon}$ then there exists a constant $c>0$ such that

\begin{equation}
\left| \mu - 8 \pi \mathfrak{a}_N M/L^3 \right| \lesssim L^{-2} \exp\left( - c N^{\epsilon} \right).
\label{eq:LemChemPotBEC1}
\end{equation}

* If $M \lesssim N^{5/6 - \epsilon}$ then we have

\begin{equation}
\left| \mu + \frac{1}{\beta M} \right| \lesssim \frac{N^{-2 \epsilon}}{\beta M}.
\label{eq:LemChemPotBEC2}
\end{equation}

We write the two-dimensional integration over $\mathbb{C}$ w.r.t. the measure $\de z = \de x \de y/ \pi$ in polar coordinates $(r,\varphi)$ and afterwards introduce the variable $x = r^2$. This allows us to write

\begin{equation}
M = \int_{\mathbb{C}} |z|^2 g(z) \de z = \frac{\int_0^{\infty} x \exp\left(- \beta \left( h x^2 - \mu x \right) \right) \de x}{\int_0^{\infty} \exp\left(- \beta \left( h x^2 - \mu x \right) \right) \de x},
\label{eq:AppBEC1}
\end{equation}

where $h = 4 \pi \mathfrak{a}_N /L^3 \sim L^{-2} N^{-1}$. A short computation shows the integral in the numerator equals

\begin{equation}
\frac{1}{2 \beta h} + \frac{\mu}{4 h} \sqrt{ \frac{\pi}{\beta h} } \exp\left( \frac{\beta \mu^2}{4 h} \right) \mathrm{erfc}\left( - \sqrt{ \frac{\beta}{h}} \frac{\mu}{2}  \right),
\label{eq:AppBEC2}
\end{equation}

where $\mathrm{erfc}(x) = (2/\sqrt{\pi}) \int_x^{\infty} \exp(-t^2) \de t$ denotes the complementary error function. For the integral in the denominator we find

\begin{equation}
\frac{1}{2} \sqrt{ \frac{\pi}{\beta h} } \exp\left( \frac{\beta \mu^2}{4 h} \right) \mathrm{erfc}\left( - \sqrt{ \frac{\beta}{h}} \frac{\mu}{2} \right).
\label{eq:AppBEC3}
\end{equation}

Let us introduce the notation $\eta =\mu \sqrt{\beta/(4h)}$. Using (<ref>) and (<ref>), we bring (<ref>) to the form

\begin{equation}
\sqrt{\pi \beta h} M = \frac{ 1 + \sqrt{\pi} \eta \exp(\eta^2) \mathrm{erfc}(- \eta) }{\exp(\eta^2) \mathrm{erfc}(- \eta)} \eqqcolon \Upsilon(\eta).
\label{eq:AppBEC4}
\end{equation}

The function $\Upsilon$ is strictly positive, strictly monotone increasing, and satisfies $\lim_{x \to - \infty} \Upsilon(x) = 0 $ as well as $\lim_{x \to \infty} \Upsilon(x) = + \infty$. In the following we study the asymptotic behavior of the (unique) solution to this equation. We start with the parameter regime $M \gtrsim N^{5/6 + \epsilon}$, which implies $\sqrt{\pi \beta h} M \gtrsim N^{\epsilon}$.

In this case the l.h.s. of (<ref>) diverges in the limit $N \to \infty$, and hence $\eta \to \infty$. From <cit.> we know that

\begin{equation}
\frac{1}{x + \sqrt{x^2 + 2}} < \exp\left( x^2 \right) \int_x^{\infty} \exp\left( -t^2 \right) \de t \leq \frac{1}{x + \sqrt{x^2 + 4/\pi}}
\label{eq:AppBEC5}
\end{equation}

holds for $x \geq 0$. In combination with $\mathrm{erfc}(-\eta) = 2 - \mathrm{erfc}(\eta)$, this implies

\begin{equation}
2 \exp\left( \eta^2 \right) - \frac{2}{ \sqrt{\pi} \left( \eta + \sqrt{\eta^2 + 4/\pi } \right)}  \leq \exp\left( \eta^2 \right) \mathrm{erfc}(-\eta) < 2 \exp\left( \eta^2 \right) - \frac{2}{ \sqrt{\pi} \left( \eta + \sqrt{\eta^2 + 2 } \right)}
\label{eq:AppBEC6}
\end{equation}

as well as

\begin{equation}
\sqrt{\pi \beta h} M = \frac{ 1 + \sqrt{\pi} \eta \left[2 \exp\left( \eta^2 \right) + O(1/\eta) \right] }{2 \exp\left( \eta^2 \right) + O(1/\eta)}.
\label{eq:AppBEC7}
\end{equation}

We already know that $\eta \gg 1$, and hence $\eta \simeq \sqrt{\beta h} M$. Using this and our assumption $M \gtrsim N^{5/6 + \epsilon}$, which implies $\eta \gtrsim N^{\epsilon}$, we easily check that (<ref>) holds. It remains to prove (<ref>).

If $M \lesssim N^{5/6 - \epsilon}$ the l.h.s. of (<ref>) satisfies $\sqrt{\pi \beta h} M \lesssim N^{-\epsilon}$ and we therefore have $\eta \to -\infty$. To obtain the leading order behavior of $\eta$, the approximation provided by (<ref>) is not sufficiently accurate. A more precise approximation is provided by <cit.>, which implies

\begin{equation}
\sqrt{\pi} \exp( x^2 ) \mathrm{erfc}(x) = \frac{1}{x} - \frac{1}{2 x^3} + Q(x), \quad \text{ where $Q$ satisfies } \quad |Q(x)| \leq \frac{3}{4 x^5}
\label{eq:AppBEC8}
\end{equation}

for $x \geq 0$. We use this approximation in (<ref>) and find

\begin{equation}
\sqrt{\beta h} M = \frac{1}{2 | \eta |} \left( 1 + O\left(\eta^{-2} \right) \right).
\label{eq:AppBEC9}
\end{equation}

Eq. (<ref>) is a direct consequence of (<ref>). This proves our claim.

We are now prepared to give the proof of Proposition <ref>. Because of technical reasons, we prove it in a slightly more general situation.

We consider the limit $N \to \infty$, $\beta = \kappa \beta_{\mathrm{c}}$ with $\kappa \in (0,\infty)$ and $\beta_{\mathrm{c}}$ in (<ref>). The following statements hold for given $\epsilon > 0$:

* Assume that $M \gtrsim N^{5/6 + \epsilon}$. There exists a constant $c>0$ such that

\begin{equation}
F^{\mathrm{BEC}}(\beta,M,L,\mathfrak{a}_N) = 4 \pi \mathfrak{a}_N L^{-3} M^2 + \frac{\ln \left( 4 \beta \mathfrak{a}_N/L^3 \right)}{2 \beta} + O\left( L^{-2} \exp\left(- c N^{\epsilon} \right) \right).
\label{eq:FreeEnergyBECInteractingLimit}
\end{equation}

* Assume that $M \lesssim N^{5/6 - \epsilon}$. Then

\begin{equation}
F^{\mathrm{BEC}}(\beta,M,L,\mathfrak{a}_N) = - \frac{1}{\beta} \ln(M) - \frac{1}{\beta} + O\left( L^{-2} N^{2/3 - 2 \epsilon} \right)
\label{eq:FreeEnergyBECNonInteractingLimit}
\end{equation}

holds. In particular, $F^{\mathrm{BEC}}(\beta,M,L,\mathfrak{a}_N)$ is independent of $\mathfrak{a}_N$ at the given level of accuracy.

The free energy $F^{\mathrm{BEC}}(\beta,M,L,\mathfrak{a}_N)$ in (<ref>) consists of two terms. In the following, we denote the first by $\Phi(\beta,M,L,\mathfrak{a}_N)$. When we apply the same coordinate transformations that led to (<ref>), we can write it as

\begin{equation}
\Phi(\beta,M,L,\mathfrak{a}_N) = - \frac{1}{\beta} \ln\left( \int_0^{\infty} \exp\left( -\beta \left( h x^2 - \mu x \right) \right) \de x \right) = - \frac{1}{\beta} \ln \left( \frac{1}{2} \sqrt{ \frac{\pi}{\beta h} } \exp\left( \frac{\beta \mu^2}{4h} \right) \mathrm{erfc}\left( - \sqrt{ \frac{\beta}{h} } \frac{\mu}{2} \right) \right),
\label{eq:AppBEC10}
\end{equation}

where the second identity follows from the fact that the denominator in (<ref>) is given by (<ref>).

We first consider the parameter regime $M \gtrsim N^{5/6 + \epsilon}$, where $\eta \simeq \sqrt{\beta h} M \gtrsim N^{\epsilon}$. An application of (<ref>) shows that the r.h.s. of (<ref>) equals

\begin{equation}
-\frac{1}{\beta} \ln\left( \sqrt{ \frac{\pi}{\beta h} } \exp( \eta^2 ) \left( 1 + O\left( \exp(-\eta^2)/\eta \right) \right) \right) = \frac{1}{2 \beta} \ln\left( \frac{4 \beta \mathfrak{a}_N}{L^3} \right) - 4 \pi \mathfrak{a}_N M^2 L^{-3} + O\left( L^{-2} \exp \left(-c N^{2 \epsilon} \right) \right).
\label{eq:AppBEC11}
\end{equation}

From Lemma <ref> we know that

\begin{equation}
\mu M = 8 \pi \mathfrak{a}_N M^2 L^{-3} + O\left( L^{-2} \exp\left( -c N^{\epsilon} \right) \right).
\label{eq:AppBEC12}
\end{equation}

In combination, these consideration show

\begin{equation}
F^{\mathrm{BEC}}(\beta,M,L,\mathfrak{a}_N) = \Phi(\beta,M,L,\mathfrak{a}_N) + \mu M =  \frac{1}{2 \beta} \ln\left( \frac{4 \beta \mathfrak{a}_N}{L^3} \right) + 4 \pi \mathfrak{a}_N M^2 L^{-3} + O\left( L^{-2} \exp \left(-c N^{\epsilon} \right) \right),
\label{eq:AppBEC13}
\end{equation}

which proves (<ref>).

Next, we consider the case $M \lesssim N^{5/6 - \epsilon}$, where $\eta \simeq -1/(2 \sqrt{\beta h} M) \lesssim -N^{\epsilon}$. We use (<ref>) to write $\Phi$ as

\begin{equation}
\Phi(\beta,M,L,\mathfrak{a}_N) = -\frac{1}{\beta} \ln\left( \sqrt{ \frac{1}{\beta h} } \frac{1}{2 | \eta |} \left( 1 + O \left( \eta^{-2} \right) \right) \right) = -\frac{ \ln( M ) }{\beta} + O\left( N^{-2 \epsilon} / \beta \right).
\label{eq:AppBEC14}
\end{equation}

To obtain the second equality we applied Lemma <ref>. Another application of the same lemma yields

\begin{equation}
\mu M = -\frac{1}{\beta} \left( 1 + O\left( N^{-2 \epsilon} \right) \right).
\label{eq:AppBEC15}
\end{equation}

In combination, (<ref>) and (<ref>) prove (<ref>).

The last lemma provides us with a large deviations bound as well as with bounds for the moments of the distribution $\zeta$. The large deviations bound is needed in the proof of Lemma <ref> in Appendix <ref>, whereas the moment bound finds application in Section <ref> in our proof of Lemma <ref>. We recall that $\zeta$ equals $g$ in (<ref>) except that the chemical potential $\widetilde{\mu}$ is chosen s.t. $\int_{\mathbb{C}} |z|^2 \zeta(z) \de z = \widetilde{N}_0$ holds with $\widetilde{N}_0$ in (<ref>).

We consider the limit $N \to \infty$, $\beta = \kappa \beta_{\mathrm{c}}$ with $\kappa \in (0,\infty)$ and $\beta_{\mathrm{c}}$ in (<ref>), and assume that $0 \leq \widetilde{N}_0 \lesssim N$ holds. Then there exist constants $c,\tilde{c}>0$ such that

\begin{equation}
\int_{\mathbb{C}} (1+ |z|^{2} ) \mathds{1}(|z|^2 \geq cN) \zeta(z) \de z \lesssim \exp(-\tilde{c} N^{1/3} ).
\label{eq:largeDeviationsZeta}
\end{equation}

\begin{equation}
\int_{\mathbb{C}} |z|^{2k} \zeta(z) \de z \lesssim_k N^k
\label{eq:momentBoundsZeta}
\end{equation}

holds for all $k \in \mathbb{N}$.

Let us again use the notation $h = 4 \pi \mathfrak{a}_N L^{-3}$. We first consider (<ref>) with $|z|^2$ replaced by $|z|^{2k}$, $k \in \mathbb{N}_0$ and $\lesssim$ replaced by $\lesssim_k$, that is, we need to derive a bound for

\begin{equation}
\int_{\mathbb{C}} |z|^{2k} \mathds{1}(|z|^2 \geq cN) \zeta(z) \de z = \frac{\int_{cN}^{\infty} x^k \exp( -\beta(h x^2 - \widetilde{\mu} x) ) \de x}{\int_{0}^{\infty} \exp( -\beta(h x^2 - \widetilde{\mu} x) ) \de x} = \frac{\int_{cN}^{\infty} x^k \exp( -\beta h( x - \widetilde{\mu}/(2h))^2 ) \de x}{\int_{0}^{\infty} \exp( -\beta h( x - \widetilde{\mu}/(2h))^2 ) \de x}.
\label{eq:AppBEC16}
\end{equation}

To obtain the first equality, we used the same coordinate transformations as above (<ref>). Inspection of (<ref>) shows that the chemical potential assumes its largest (positive) values when $\widetilde{N}_0 \sim N$. This follows from the fact that the l.h.s. of (<ref>) is strictly increasing in $M$ and that the two maps $\widetilde{\eta} \mapsto \Upsilon(\widetilde{\eta})$ with $\Upsilon$ in (<ref>) and $\widetilde{\mu} \mapsto \widetilde{\eta} =\widetilde{\mu} \sqrt{\beta/(4h)}$ are strictly increasing. Application of part (a) of Lemma <ref> and the bound $\widetilde{N}_0 \lesssim N$ therefore show that $\widetilde{\mu}$ can be bounded from above by a constant times $\mathfrak{a}_N N L^{-3} \lesssim L^{-2}$.

Using this, we see that for $c>0$ large enough and $x \geq cN$, we have $x - \mu/2h \geq x/2$. We insert this bound on the r.h.s. of (<ref>) and find

\begin{equation}
\int_{\mathbb{C}} |z|^{2k} \mathds{1}(|z|^2 \geq cN)  \zeta(z) \de z \lesssim_k \exp\left( -\frac{\beta h c^2 N^2}{16} \right) \frac{ \int_{0}^{\infty} \exp( -(\beta h/2)( x - \widetilde{\mu}/(2h))^2 ) \de x}{\int_{0}^{\infty} \exp( -\beta h( x - \widetilde{\mu}/(2h))^2 ) \de x},
\label{eq:AppBEC17}
\end{equation}

where the fraction on the r.h.s. equals $\sqrt{2}$ times

\begin{equation}
\frac{ \int_{-(\widetilde{\mu}/2) \sqrt{\beta/(2h)}}^{\infty} \exp( -x^2 ) \de x}{\int_{-(\widetilde{\mu}/2) \sqrt{\beta/h}}^{\infty} \exp( -x^2 ) \de x}.
\label{eq:AppBEC18}
\end{equation}

If $\widetilde{\mu} \geq 0 $ we obtain an upper bound when we replace the lower integration boundary in the numerator by $-\infty$ and that in the denominator by $0$. This yields an upper bound of order $1$. If $\widetilde{\mu}<0$ we apply (<ref>) on the r.h.s. of (<ref>) and see that it is bounded from above by a constant times $\exp(\beta \widetilde{\mu}^2/(8h))$. In combination, these considerations imply the bound

\begin{equation}
\int_{\mathbb{C}} |z|^{2k} \mathds{1}(|z|^2 \geq cN) \zeta(z) \de z \lesssim_k \exp\left( -\frac{\beta h c^2 N^2}{16} \right) \max\{ 1, \exp(\beta \widetilde{\mu}^2/(8h))  \}.
\label{eq:AppBEC19}
\end{equation}

When we assume that $0 > \widetilde{\mu} \geq -C/L^2$ for some $C>0$ and choose $c$ large enough, then (<ref>) proves our claim. It remains to consider the case $\widetilde{\mu}<-C/L^2$.

In this case we start with the term after the first equality sign in (<ref>). We pick $c_1>0$ and realize that it is bounded from above by

\begin{align}
\frac{\int_{cN}^{\infty} x^k \exp( \beta \widetilde{\mu} x ) \de x}{\int_{0}^{c_1 N} \exp( -\beta(h x^2 - \widetilde{\mu} x )) \de x} &\leq  \frac{ \exp(\beta h c_1^2 N^2) \int_{cN}^{\infty} x^k \exp( \beta \widetilde{\mu} x ) \de x}{  \int_{0}^{c_1N} \exp( \beta \widetilde{\mu} x )) \de x} \nonumber \\
&\lesssim_k \frac{\exp(\beta h c_1^2 N^2) (|\beta \widetilde{\mu}|^{-k-1} + N^k) \exp(\beta \widetilde{\mu} cN)}{1-\exp(\beta \widetilde{\mu} c_1 N)} \lesssim_k \exp(-\tilde{c} N^{1/3}). \label{eq:AppBEC20}
\end{align}

for some $\tilde{c}>0$. In the last step we used $\widetilde{\mu}<-C/L^2$ and that $c_1$ can be chosen as small as we wish. In combination with our previous considerations, this proves

\begin{equation}
\int_{\mathbb{C}} |z|^{2k} \mathds{1}(|z|^2 \geq cN) \zeta(z) \de z  \lesssim_k \exp(-\tilde{c} N^{1/3}) \label{eq:AppBEC21}
\end{equation}

for all $k \in \mathbb{N}_0$. In particular, (<ref>) holds.

Eq. (<ref>) follows from (<ref>) when we use the decomposition

\begin{equation}
\int_{\mathbb{C}} |z|^{2k} \zeta(z) \de z = \int_{\{ |z|^2 \leq cN \}} |z|^{2k} \zeta(z) \de z + \int_{\{ |z|^2 > cN \}} |z|^{2k} \zeta(z) \de z \lesssim_k N^k + \exp(-\tilde{c} N^{1/3} ).
\end{equation}

§ THE EXPECTED PARTICLE NUMBER IN THE TRIAL STATE

In this appendix we prove Lemma <ref>. An essential ingredient of the proof are large deviations bounds for $G_{\mathrm{B}}(z)$ in (<ref>) and $G_{\mathrm{free}}$ in (<ref>). Before we state them, we define

\begin{equation}\label{eq:defNBNI}
\cN_{\mathrm{B}}=\sum_{p\in P_\mathrm{B}}a^*_pa_p\qquad\text{and}\qquad \cN_{\mathrm{I}}=\sum_{p\in P_\mathrm{I}}a^*_pa_p.
\end{equation}

We consider the limit $N \to \infty$, $\beta = \kappa \beta_{\mathrm{c}}$ with $\kappa \in (0,\infty)$ and $\beta_{\mathrm{c}}$ in (<ref>). For any $c>0, r \in \mathbb{N}$ we have

\begin{equation}
\tr[(1+\cN_\mathrm{B})\mathds{1}(\cN_\mathrm{B}\geq cN)G_{\mathrm{B}}(z)] \lesssim_r N^{-r}\big(N^{2/3+\delta_B}\big)^{r+1}.
\label{eq:plybound}
\end{equation}

Moreover, there exist positive constants $c,\tilde{c}>0$ such that

\begin{equation}
\tr[(1+\cN_\mathrm{I})\mathds{1}(\cN_\mathrm{I}\geq cN)G_{\mathrm{free}}] \lesssim \exp(-\tilde{c} N^{1/3} ).
\label{eq:expgf}
\end{equation}

For the sake of simplicity, we give the proof of the first bound with $1+\mathcal{N}_{\mathrm{B}}$ replaced by $\mathcal{N}_{\mathrm{B}}$, and similarly for the second bound. For any $r\geq 1$, we have $\mathds{1}(\cN_\mathrm{B}\geq cN)\leq \cN_\mathrm{B}^r(cN)^{-r}$. Hence,

\begin{align}
\tr[\cN_\mathrm{B}\mathds{1}(\cN_\mathrm{B}\geq cN)G_{\mathrm{B}}(z)]&\lesssim (cN)^{-r}\tr[\cN_\mathrm{B}^{r+1}G_{\mathrm{B}}(z)] \nonumber \\
&= (cN)^{-r}\sum_{p_1, \dots, p_{r+1}\in P_\text{B}}\tr[a^*_{p_1}a_{p_1}\dots a^*_{p_{r+1}}a_{p_{r+1}} G_{\mathrm{B}}(z)].
\end{align}

After normal ordering and an application of Wick's theorem, we can use (<ref>) and (<ref>) to see that

\begin{align}
\tr[\cN_\mathrm{B}\mathds{1}(\cN_\mathrm{B}\geq cN)G_{\mathrm{B}}(z)]&\lesssim_r N^{-r}\big(N^{2/3+\delta_B}\big)^{r+1}
\end{align}

holds. This proves the first bound in (<ref>)

Next, we prove the second bound. Let $0<k\leq L^{-2}\beta$ and observe that

\begin{align}
k\tr[\cN_\mathrm{I}\mathds{1}(\cN_\mathrm{I}\geq cN)G_{\mathrm{free}}]&\leq  k\tr[(\cN_\mathrm{I}-cN)\mathds{1}(\cN_\mathrm{I}\geq cN)G_{\mathrm{free}}]+ckN \tr[\mathds{1}(\cN_\mathrm{I}\geq cN)G_{\mathrm{free}}]\nonumber\\
&\leq (1+ckN)\tr[\exp(k(\cN_\mathrm{I}-cN))G_{\mathrm{free}}].\label{eq:d1}
\end{align}

The trace on the r.h.s. can be written as

\begin{align}
\tr[\exp(k(\cN-cN)G_{\mathrm{free}})] &= \exp(-k cN)\frac{\tr_{\mathscr{F}_{\mathrm{I}}}[ \exp(-\beta \de \Gamma( \mathds{1}( -\mathrm{i} \nabla \in P_{\mathrm{I}})(- \Delta - \mu_0-k\beta^{-1} )))]}{\tr_{\mathscr{F}_{\mathrm{I}}}[ \exp(-\beta \de \Gamma( \mathds{1}( -\mathrm{i} \nabla \in P_{\mathrm{I}})(- \Delta - \mu_0 )))]} \nonumber \\
&= \exp(-k cN) \exp\left(\beta\,\big(\Phi(\mu_0)-\Phi(\mu_0+k\beta^{-1})\big)\right),\label{eq:d2}
\end{align}

\begin{equation}
\Phi(\mu)=-\beta^{-1}\ln\tr_{\mathscr{F}_{\mathrm{I}}}\Big[\exp\Big(-\beta \sum_{p\in P_{\mathrm{I}}}(p^2-\mu)a^*_pa_p)\Big)\Big] = \beta^{-1} \sum_{p\in P_{\mathrm{I}}}\Big[\ln\big(1-\exp(-\beta(p^2-\mu)) \big) \Big].
\end{equation}

Using that $\Phi(\mu)$ is a concave and monotone decreasing function of $\mu$, we obtain the lower bound

\begin{align}
\beta \Phi(\mu_0+k\beta^{-1}) & \geq \beta \Phi(\mu_0) - \sum_{p\in P_{\mathrm{I}}}\frac{k}{\exp(\beta(p^2-\mu_0-k\beta^{-1}))-1} \nonumber \\
&\geq \beta \Phi(\mu_0) - CN k
\label{eq:d3}
\end{align}

for some $C>0$. To come to the last line, we used used $0<k\leq L^{-2}\beta$ and applied Lemma <ref>. In combination, these considerations show

\begin{equation}\label{eq:expgb1}
\tr[e^{k(\cN-cN)}G_{\mathrm{free}}]\lesssim  e^{-k N (c - C) }.
\end{equation}

When we choose $k= L^{-2} \beta$ and $c>C$ in the above equation this proves the second bound in (<ref>).

We are now prepared to give the proof of Lemma <ref>.

We use $[\cN,B]=0$ to write

\begin{align}\label{eq:trng}
\tr [\cN \G] = \sum_{\alpha} \lambda_{\alpha}\frac{\braket{(1+B)\psi_{\alpha},\cN(1+B)\psi_{\alpha}}} {\braket{(1+B)\psi_{\alpha}, (1+B)\psi_{\alpha}}}= \sum_{\alpha} \lambda_{\alpha}\frac{\langle\psi_{\alpha},(\cN+B^*\cN B)\psi_{\alpha}\rangle} {\braket{(1+B)\psi_{\alpha}, (1+B)\psi_{\alpha}}},
\end{align}

and apply the lower bound $\|(1+B)\psi_{\alpha}\|^2=\langle\psi_\alpha, (1+B^*B)\,\psi_\alpha \rangle\geq 1$ to see that the r.h.s. is bounded from above by $\mathrm{Tr}[ (\cN+B^* B\cN) \Gamma_0 ]$. Using the definition of $B$ in (<ref>),  Lemma <ref> and the bound for $\eta_p$ in (<ref>), we find

\begin{align}\label{eq:NBB}
\text{Tr} [B^* B\cN \Gamma_0] &= \frac{1}{4L^6} \sum_{\substack{q\in\L^*,p_1,p_2 \in P_{\mathrm{H}},\\ \ u_1,v_1,u_2,v_2 \in P_{\mathrm{L}}}}\eta_{p_1}\eta_{p_2} \tr [ a_{u_1}^* a_{v_1}^* a_{u_1+p_1}a_{v_1-p_1} a_{u_2+p_2}^*a_{v_2-p_2}^* a_{u_2 }a_{v_2}a^*_q a_q \Gamma_0 ] \nonumber\\
&=\frac{1}{2 L^6} \sum_{\alpha} \lambda_{\alpha}\sum_{\substack{p_1 \in P_{\mathrm{H}},\\ \ q,u_1,v_1,u_2 \in P_{\mathrm{L}}}}\eta_{p_1}\eta_{p_1+u_1-u_2} \tr[ a^*_{u_1} a_{v_1}^* a_{u_2 }a_{v_1+u_1-u_2}a^*_qa_q \Gamma_0 ] \lesssim N^{\delta_H},
\end{align}

and therefore $\tr [\cN \G] - \tr [\cN \G_0]  \lesssim N^{\delta_H}$. It remains to prove a lower bound.

To that end, we note that $B^*\cN B \geq 0$ and $\mathds{1}(\cN_0\leq cN)\mathds{1}(\cN_{\text{B}}\leq cN)\mathds{1}(\cN_{\text{I}}\leq cN) \leq 1$ implies

\begin{equation}
\tr [\cN \G] \geq \sum_{\alpha} \lambda_{\alpha}\frac{\langle\psi_{\alpha},\cN\mathds{1}(\cN_0\leq cN)\mathds{1}(\cN_{\text{B}}\leq cN)\mathds{1}(\cN_{\text{I}}\leq cN)\psi_{\alpha}\rangle} {\braket{(1+B)\psi_{\alpha}, (1+B)\psi_{\alpha}}}
\label{eq:AppD2}
\end{equation}

for sufficiently large $c>0$. Let $\cN_0=a^*_0a_0$ and recall the definition of $\cN_{\mathrm{B}}$ and $\cN_{\mathrm{I}}$ in (<ref>). In combination, (<ref>), $1/(1+x) \geq 1 - x$ for $x\geq 0$ and the bound $\langle\psi_{\alpha},\cN\mathds{1}(\cN_0\leq cN)\mathds{1}(\cN_{\text{B}}\leq cN)\mathds{1}(\cN_{\text{I}}\leq cN)\psi_{\alpha}\rangle\leq 3cN$ allow us to show that

\begin{align}
\tr [\cN \G] &\geq \sum_{\alpha} \lambda_{\alpha}\langle\psi_{\alpha},\cN\mathds{1}(\cN_0\leq cN)\mathds{1}(\cN_{\text{B}}\leq cN)\mathds{1}(\cN_{\text{I}}\leq cN)\psi_{\alpha}\rangle(1-\langle\psi_\alpha, B^*B\,\psi_\alpha \rangle) \nonumber \\
&\geq \tr [\cN \G_0]- \tr\left[ \cN\big\{\mathds{1}(\cN_0\geq cN)+\mathds{1}(\cN_\textsc{B}\geq cN)+\mathds{1}(\cN_\textsc{I}\geq cN)\big\} \Gamma_0 \right] -3cN \tr [ B^*B\, \Gamma_0 ]. \label{eq:lbng}
\end{align}

The last contribution can be bounded from below by $-CN^{\delta_H}$ (this can be seen similarly as for (<ref>)).

To obtain a bound for the second term on the r.h.s. we apply Lemmas <ref> and <ref>: we have

\begin{equation}
\tr\left[ \cN\big[\mathds{1}(\cN_\textsc{B}\geq cN)+\mathds{1}(\cN_\textsc{I}\geq cN)\big] \Gamma_0 \right] \lesssim_r N^{1-r}\big(N^{2/3+\delta_{\mathrm{B}}}\big)^{r+1} + \exp(-\tilde{c} N^{1/3} ) \lesssim 1,
\label{eq:AppD5b}
\end{equation}

provided $c>1$ and $r \in \mathbb{N}$ are chosen large enough. To obtain the second bound we also used the assumption $\delta_{\mathrm{B}} < 1/3$. Next, we consider

\begin{equation}
\int_{\mathbb{C}} \tr[ (1 + \mathcal{N}_0) \mathds{1}(\mathcal{N}_0 \geq cN) |z \rangle \langle z | ] \zeta(z) \de z = \int_{\mathbb{C}} (1 + |z|^2) \tr[ \mathds{1}(\mathcal{N}_0 + 1 \geq cN ) |z \rangle \langle z | ] \zeta(z) \de z,
\label{eq:AppD3}
\end{equation}

where we used $a_0 \mathds{1}(\mathcal{N}_0 \geq cN) = \mathds{1}(\mathcal{N}_0 + 1 \geq cN) a_0$.

Pick $c'>0$. An application of Lemma <ref> shows that the term on the r.h.s. of (<ref>) is bounded from above by

\begin{align}
&\int_{\{ |z|^2 \leq c' N \} } (1 + |z|^2) \tr[ \mathds{1}(\mathcal{N}_0 + 1 \geq cN) |z \rangle \langle z | ] \zeta(z) \de z + \int_{\{ |z|^2 > c' N \}} (1+|z|^2) \zeta(z) \de z \nonumber \\
&\hspace{4cm}\leq (1 + c' N) \int_{\{ |z|^2 \leq c' N \} } \tr[ \mathds{1}(\mathcal{N}_0 + 1 \geq cN) |z \rangle \langle z | ] \zeta(z) \de z + C \exp(- \tilde{c} N^{1/3})
\label{eq:AppD4}
\end{align}

for some $\tilde{c} > 0$ as long as $c>0$ is chosen large enough. It remains to consider the first term on the r.h.s. of (<ref>).

We evaluate the trace in the eigenbasis $\{ |n \rangle \}_{n \in \mathcal{N}_0}$ of $\mathcal{N}_0$:

\begin{align}
\int_{\{ |z|^2 \leq c' N \} } \tr[ \mathds{1}(\mathcal{N}_0 + 1 \geq cN) |z \rangle \langle z | ] \zeta(z) \de z &= \int_{\{ |z|^2 \leq c' N \} } \sum_{n=0}^{\infty} \mathds{1}(n \geq cN - 1) | \langle z,n \rangle |^2 \zeta(z) \de z \nonumber \\
&=\int_{\{ |z|^2 \leq c' N \} } \exp(-|z|^2) \left[ \sum_{n \geq cN-1} \frac{|z|^{2n}}{n!} \right] \zeta(z) \de z.
\label{eq:AppD5}
\end{align}

To come to the second line, we also used the identity $| \langle z,n \rangle |^2 = \exp(-|z|^2) |z|^{2n} / n!$. An application of Taylor's theorem with an explicit form of the remainder allows us to see that the series on the r.h.s. is bounded from above by $\exp(|z|^2) |z|^{2 cN}/(cN !)$ (here we assume for the sake of simplicity that $c \in \mathbb{N}$). An application of Stirling's approximation formula therefore shows

\begin{equation}
\int_{\{ |z|^2 \leq c' N \} } \tr[ \mathds{1}(\mathcal{N}_0 + 1 \geq cN) |z \rangle \langle z | ] \zeta(z) \de z \leq \frac{(c'N)^{cN}}{cN!} \lesssim N^{-1/2} \left( \frac{c' e}{c} \right)^{cN},
\label{eq:AppD6}
\end{equation}

which is exponentially small in $N$ as long as $c' e/c < 1$ holds. When we put (<ref>)–(<ref>) and (<ref>) together, we find the bound

\begin{equation}
\tr\left[ \cN\big[\mathds{1}(\cN_0\geq cN)+\mathds{1}(\cN_\textsc{B}\geq cN)+\mathds{1}(\cN_\textsc{I}\geq cN)\big] \Gamma_0 \right] \lesssim 1.
\label{eq:AppD7}
\end{equation}

In combination with (<ref>) and the assumption $\delta_{\mathrm{H}} > 0$ this proves $\tr [\cN \G] - \tr [\cN \G_0]  \gtrsim -N^{\delta_{\mathrm{H}} }$. We put this result and (<ref>) together and obtain a proof of (<ref>).

[1] M. Abramowitz, I. A. Stegun, eds., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Applied Mathematics Series. Vol. 55, Washington D.C.; New York: United States Department of Commerce, National Bureau of Standards; Dover Publications (1972)

[2] A. Adhikari, C. Brennecke, B. Schlein, Bose-Einstein Condensation Beyond the Gross-Pitaevskii Regime, Ann. Henri Poincaré 22, 1163 (2021)

[3] J. Agerskov, R. Reuvers, J. P. Solovej, Ground state energy of dilute Bose gases in 1d, 	arXiv:2203.17183 [math-ph] (2022)

[4] M.H. Anderson, J.R. Ensher, M.R. Matthews, C.E. Wieman, E.A. Cornell, Observation of bose-einstein condensation in a dilute atomic vapor, Science 269, 198 (1995)

[5] G. Basti, S. Cenatiempo, A. Olgiati, G. Pasqualetti, B. Schlein, A second order upper bound for the ground state energy of a hard-sphere gas in the Gross-Pitaevskii regime, Commun. Math. Phys. 399, 1 (2023)

[6] G. Basti, S. Cenatiempo, B. Schlein, A new second order upper bound for the ground state energy of dilute Bose gases, Forum Math. Sigma 9, e74 (2021)

[7] N. Benedikter, G. de Oliveira, B. Schlein, Quantitative Derivation of the Gross-Pitaevskii Equation, Commun. Pure Appl. Math. 68, 1399 (2015)

[8] N. Benedikter, M. Porta, B. Schlein, Effective Evolution Equations from Quantum Dynamics, Springer, Berlin (2016)

[9] F.A. Berezin, Covariant and contravariant symbols of operators, Izv. Akad. Nauk, Ser. Mat. 36, 1134–1167 (1972); English translation: USSR Izv. 6, 1117 (1973). F.A. Berezin, General concept of quantization, Commun. Math. Phys. 40, 153 (1975)

[10] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein, Bogoliubov Theory in the Gross-Pitaevskii Limit, Acta Math. 222, 219 (2019)

[11] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein, The Excitation Spectrum of Bose Gases Interacting Through Singular Potentials,
J. Eur. Math. Soc. 22, 2331 (2020)

[12] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein, Optimal Rate for Bose–Einstein Condensation in the Gross–Pitaevskii Regime, Commun. Math. Phys. 376, 1311 (2020)

[13] C. Boccato, R. Seiringer, The Bose Gas in a Box with Neumann Boundary Conditions, Ann. Henri Poincaré 24, 1505 (2023)

[14] N. N. Bogoliubov, On the theory of superfluidity, lzv. Akad. Nauk USSR 11, 77 (1947). Eng. Trans. J. Phys. (USSR) 11, 23 (1947)

[15] C. Brennecke, M. Caporaletti, B. Schlein, Excitation Spectrum for Bose Gases beyond the Gross-Pitaevskii Regime, Rev. Math. Phys. 34, No. 09, 2250027 (2022)

[16] C. Brennecke, B. Schlein, Gross–Pitaevskii dynamics for Bose–Einstein condensates, Anal. PDE 12, No. 6, 1513 (2019)

[17] C. Brennecke, B. Schlein, S. Schraven, Bose-Einstein Condensation with Optimal Rate for Trapped Bosons in the Gross-Pitaevskii Regime, Math. Phys. Anal. Geom. 25, 12 (2022)

[18] C. Brennecke, B. Schlein, S. Schraven, Bogoliubov Theory for Trapped Bosons in the Gross-Pitaevskii Regime, Ann. Henri Poincaré 23, 1583 (2022)

[19] M. Caporaletti, A. Deuchert and B. Schlein, Dynamics of mean-field bosons at positive temperature, Ann. Inst. H. Poincaré Anal. Non Linéaire (2023)

[20] C. Caraci, S. Cenatiempo, B. Schlein, Bose–Einstein Condensation for Two Dimensional Bosons in the Gross–Pitaevskii Regime, J. Stat. Phys. 183, 39 (2021)

[21] C. Caraci, S. Cenatiempo, B. Schlein, The Excitation Spectrum of Two-Dimensional Bose Gases in the Gross–Pitaevskii Regime, Ann. Henri Poincaré, Online first (2023)

[22] S. Chatterjee, P. Diaconis, Fluctuations of the Bose–Einstein condensate, J. Phys. A: Math. Theor. 47, 085201 (2014)

[23] K. B. Davis, M.-O. Mewes, M. R. Andrews, N.J. van Druten, D.S. Durfee, D.M. Kurn, W. Ketterle, Bose-Einstein Condensation in a Gas of Sodium Atoms, Phys. Rev. Lett. 75, 3969 (1995)

[24] A. Deuchert, S. Mayer, R. Seiringer, The free energy of the two-dimensional dilute Bose gas. I. Lower bound, Forum Math. Sigma 8, e20 (2020)

[25] A. Deuchert, R. Seiringer, Gross-Pitaevskii Limit of a Homogeneous Bose Gas at Positive Temperature, Arch. Ration. Mech. 236(3), 1217 (2020)

[26] A. Deuchert, R. Seiringer, J. Yngvason, Bose-Einstein Condensation for a Dilute Trapped Gas at Positive Temperature, Commun. Math. Phys. 368, 723 (2019)

[27] F. J. Dyson, Ground-State Energy of a Hard-Sphere Gas, Phys. Rev. 106, 20 (1957)

[29] L. Erdös, B. Schlein, H.-T. Yau, Rigorous derivation of the Gross-Pitaevskii equation with a large interaction potential, J. Amer. Math. Soc. 22, 1099 (2009)

[30] L. Erdös, B. Schlein, H.-T. Yau, Derivation of the Gross-Pitaevskii equation for the dynamics of Bose-Einstein condensate, Ann. of Math. 172, 291 (2010)

[31] S. Fournais, T. Girardot, L. Junge, L. Morin, M. Olivieri, The Ground State Energy of a Two-Dimensional Bose Gas, arXiv:2206.11100 [math-ph] (2022)

[32] S. Fournais, J.P. Solovej, The energy of dilute Bose gases, Ann. Math. 192, No. 3, 893 (2020)

[33] S. Fournais, J.P. Solovej, The energy of dilute Bose gases II: The general case, Invent. math. 232, 863 (2023)

[34] P. Grech, R. Seiringer, The excitation spectrum for weakly interacting bosons in a trap, Commun. Math. Phys. 322 (2), 559 (2013)

[35] F. Haberberger, C. Hainzl, P. T. Nam, R. Seiringer, A. Triay, The free energy of dilute Bose gases at low temperatures, arXiv:2304.02405 [math-ph] (2023)

[36] C. Hainzl, Another proof of BEC in the GP-limit, J. Math. Phys. 62, 051901 (2021)

[37] C. Hainzl, B. Schlein, A. Triay, Bogoliubov Theory in the Gross-Pitaevskii Limit: a Simplified Approach, Forum Math. Sigma 10, e90 (2022)

[38] M. Jeblick, N. Leopold, P. Pickl, Derivation of the Time Dependent Gross-Pitaevskii Equation in Two Dimensions, Commun. Math. Phys. 372, 1 (2019)

[39] T. D. Lee, K. Huang, C. N. Yang, Eigenvalues and eigenfunctions of a bose system of hard spheres and its low-temperature properties, Phys. Rev. 106, 1135 (1957)

[40] T. D. Lee, C. N. Yang, Low-Temperature Behavior of a Dilute Bose System of Hard Spheres. I. Equilibrium Properties, Phys. Rev. 112, No. 5, 1419 (1958)

[41] E.H. Lieb, The classical limit of quantum spin systems, Commun. Math. Phys. 31, 327 (1973)

[42] E.H. Lieb, R. Seiringer, Proof of Bose-Einstein Condensation for Dilute Trapped Gases, Phys. Rev. Lett. 88, 170409 (2002)

[43] E.H. Lieb, R. Seiringer, Derivation of the Gross-Pitaevskii equation for rotating Bose gases, Commun. Math. Phys. 264, 505 (2006)

[44] E.H. Lieb, R. Seiringer, J. Yngvason, Superfluidity in Dilute Trapped Bose Gases, Phys. Rev. B 66, 134529 (2002)

[45] E.H. Lieb, R. Seiringer, J. P. Solovej, J. Yngvason, The Mathematics of the Bose Gas and its Condensation, Birkhäuser, Basel (2005)

[46] E.H. Lieb, R. Seiringer, J. Yngvason, Bosons in a trap: A rigorous derivation of the Gross-Pitaevskii energy functional, Phys. Rev. A 61, 043602 (2000)

[47] E.H. Lieb, R. Seiringer, J. Yngvason, Justification of c-Number Substitutions in Bosonic Hamiltonians, Phys. Rev. Lett. 94, 080401 (2005).

[48] E.H. Lieb, J. Yngvason, Ground State Energy of the Low Density Bose Gas, Phys. Rev. Lett. 80, 2504 (1998)

[49] E. H. Lieb, J. Yngvason, The ground state energy of a dilute two-dimensional Bose gas, J. Stat. Phys. 103, 509 (2001)

[50] S. Mayer, R. Seiringer, The free energy of the two-dimensional dilute Bose gas. II. Upper bound, J. Math. Phys. 61, 061901 (2020)

[51]
P. T. Nam, M. Napiórkowski, J. Ricaud, A. Triay, Optimal rate of condensation for trapped bosons in the Gross-Pitaevskii regime, Anal. PDE 15, 1585 (2022)

[52] P. T. Nam, N. Rougerie, R. Seiringer, Ground states of large bosonic systems: The Gross-Pitaevskii limit revisited, Anal. PDE 9, 459 (2016)

[53] P. T. Nam, A. Triay, Bogoliubov excitation spectrum of trapped Bose gases in the Gross-Pitaevskii regime, arXiv:2106.11949 [math-ph] (2021)

[54] M. Napiórkowski, Dynamics of interacting bosons: a compact review, Lecture Notes Series, Institute for Mathematical Sciences, National University of Singapore, Density Functionals for Many-Particle Systems, pp. 117 (2023)

[55] M. Napiórkowski, R. Reuvers, J. P. Solovej, The Bogoliubov free energy functional II: the dilute limit, Commun. Math. Phys. 360, 347 (2017)

[56] P. Pickl, Derivation of the time dependent Gross Pitaevskii equation with external fields, Rev. Math. Phys. 27, 1550003 (2015)

[58] N. Rougerie, De Finetti theorems, mean-field limits and Bose-Einstein condensation, Lecture notes, preprint arXiv:1506.05263 (2015)

[59] J. Schmitt, T. Damm, D. Dung, F. Vewinger, J. Klaers, M. Weitz, Observation of Grand-Canonical Number Statistics in a Photon Bose-Einstein Condensate, Phys. Rev. Lett. 112, 030401 (2014)

[60] R. Seiringer, Ground state asymptotics of a dilute, rotating gas, J. Phys. A: Math. Gen. 36, 9755 (2003).

[61] R. Seiringer: The thermodynamic pressure of a dilute Fermi gas, Commun. Math. Phys. 261, 729 (2006)

[62] R. Seiringer, Free Energy of a Dilute Bose Gas: Lower Bound, Commun. Math. Phys. 279, 595 (2008)

[63] R. Seiringer, The excitation spectrum for weakly interacting bosons, Commun. Math. Phys. 306,  565 (2011)

[64] E. C. I. van der Wurff, A.-W. de Leeuw, R. A. Duine, H. T. C. Stoof, Interaction Effects on Number Fluctuations in a Bose-Einstein Condensate of Light, Phys. Rev. Lett. 113, 135301 (2014)

[65] H.-T. Yau, J. Yin, The Second Order Upper Bound for the Ground Energy of a Bose Gas, J. Stat. Phys. 136, 453 (2009)

[66] J. Yin, Free Energies of Dilute Bose Gases: Upper Bound, J. Stat. Phys. 141, 683 (2010)

(Chiara Boccato) Università degli Studi di Milano

Via Saldini 50, 20133 Milano, Italy

E-mail address:

(Andreas Deuchert) Institute of Mathematics, University of Zurich

Winterthurerstrasse 190, 8057 Zurich, Switzerland

E-mail address:

(David Stocker)

E-mail address:


# Super-localisation of a point-like emitter in a resonant environment :
correction of the mirage effect

Lorenzo Baldassari Department of Computational Applied Mathematics and
Operations Research, Rice University, Houston, TX 77005 USA. Pierre Millien
ESPCI Paris, PSL Research University, CNRS, Institut Langevin, France  Alice
L. Vanel Department of Mathematics, ETH Zürich, Rämistrasse 101, CH-8092
Zürich, Switzerland.

(August 27, 2024)

###### Abstract

In this paper, we show that it is possible to overcome one of the fundamental
limitations of super-resolution microscopy techniques: the necessity to be in
an _optically homogeneous_ environment. Using recent modal approximation
results from [10, 7] we show as a proof of concept that it is possible to
recover the position of a single point-like emitter in a _known resonant
environment_ from far-field measurements with a precision two orders of
magnitude below the classical Rayleigh limit. The procedure does not involve
solving any partial differential equation, is computationally light
(optimisation in $\mathbb{R}^{d}$ with $d$ of the order of $10$) and therefore
suited for the recovery of a very large number of single emitters.

Mathematics Subject Classification (MSC2000). 35R30, 35C20.

Keywords. plasmonic resonance, biomedical imaging, mirage effect, super
resolution, super localisation

##  1 Introduction

###  1.1 Context and position of the problem

Super-resolution imaging [20, 19] has had a huge impact on biological imaging
(its contribution was recognised by the Nobel Prize in $2014$ [15]). With
these techniques, biological structures can be imaged with a recovery of
details that are two order of magnitude below the classical diffraction limit
of microscopy. The technique is now being used in other areas of material
science such as polymer dynamics study [17, 29] or plasmonics [11]. We refer
the reader to the review paper [30] and references therein for an overview of
the recent research on plasmonics using super-resolved microscopy. Super
resolution microscopy’s principle is to excite a single point-like emitter in
a diffraction-limited spot, and to find the position of the emitter by fitting
a point spread function over the spot. The procedure is iterated several
thousand times to get a precise image, similar to the technique of
_pointillism_ in painting. There are several types of microscopy depending on
whether the excitation of the emitters is tailored on their emissive
properties. The latter is what can be called localisation-based super-
resolution imaging, and is known under several acronyms such as STORM, PALM,
PAINT, …

Single emitters such as quantum dots or fluorescent molecules interact with
their environment in several ways.

It has been known that the lifetime emission of single fluorescent molecules
depends on the local density of states (LDOS) and therefore they are excellent
probes for the near-field around a nano-structure. In [11] the authors
experimentally measured the lifetime of single emitters around a plasmonic rod
while simultaneously super-localising the emitters. They were able to
construct a map of the fluorescence lifetime with a resolution approaching
$20$ nm.

However, the single emitter can also excite the plasmonic resonances of the
nano-structure nearby if its emission frequency is close to a plasmonic
resonant frequency. In this case, the far-field image does not only consist of
the point spread function of the single emitter but also contains artefacts
due to the electromagnetic field scattered by the resonant nanoparticle. It is
demonstrated in [27] that coupling between the plasmonic particle and the
fluorescent molecule leads to an incorrect molecular localisation of up to
$30$ nm, namely a single-molecule mirage. We refer to the review [30] and
references therein for more details.

The only current solution to this problem is to perform a least-square
optimisation on the position of the emitter, simulating the measurements at
each step by solving the full Maxwell system with a different source position
at every step. This solution could give satisfying results for the
localisation of one emitter, but the computational cost makes this method
impossible to use to get a real super-resolution image, as the single emitter
localisation procedure needs to be done tens of thousands times to get an
image.

The challenge is therefore to find a model that is accurate enough to
precisely describe the interaction between a single emitter and a resonant
structure, yet computationally light enough to be able to retrieve the
positions of tens of thousands of emitters.

###  1.2 Mathematical challenges and recent advances

A natural approach to tackle this problem is to model the response of the
plasmonic system with excitation independant quantities. A natural way to do
this is to use modal analysis. The idea is to express the measured field as a
sum of the field generated by the single point-wise emitter and the
contribution from the nanoparticle described by a few modes. The obvious
strength of this description is that the modes can be pre-computed, and the
procedure of solving the forward problem can be reduced to computing the modal
excitation coefficients.

The main challenge is that the mirage phenomenon happens in a regime where the
extinction cross section of the nanoparticle is dominated by the scattering
cross section. This happens for larger nanoparticles, of typical size larger
than $50$ nm (see [6, section 4] for more details).

While there has been plenty of work on the mathematical description of
plasmonic resonances as an eigenvalue problem [18, 8, 1], most of these works
rely on an electro-static approximation that does not describe well the
behaviour of metallic nanoparticles when their size becomes _small but
comparable to the wavelength_ of the electromagnetic field. For instance, when
the size of the nanoparticle is larger than one tenth of the wavelength
(around $50$ nm in practice), the quasi-static theory breaks down [21, 23]. It
is only recently, via the introduction of a perturbative spectral method, that
the field scattered by larger nanoparticles has been precisely described [5,
10, 7]. In this paper, using these recent perturbed modal decompositions, we
tackle the problem of the super-localisation of a single emitter near a
plasmonic structure in the regime that is relevant for applications.

###  1.3 Contribution of the paper

Using a simple scalar _toy model_ that still encompasses all the mathematical
complexity of plasmonic resonances, we show that it is possible to recover the
position of a single point-like emitter near a resonant structure up to a
precision that is two orders of magnitudes below the diffraction limit with a
low computational complexity.

The paper is structured as followed: in section 2 we present the mathematical
modelling of the physical system and justify its relevance. Then, in section
3, we give the modal approximation of the field. In section 4, we move on to
introduce the imaging functionals that allow the reconstruction of the
dipole’s position and orientation. Finally, we show numerical experiments that
corroborate our results in section 5 and give concluding remarks in section 6.

##  2 Model

###  2.1 Notations

We consider a system composed of a single point-like dipole source (modelling
a fluorescent molecule) and a plasmonic particle embedded in a homogeneous
medium. The plasmonic particle occupies a bounded simply connected domain
$D\subset\mathbb{R}^{2}$ of class $C^{1,\alpha}$ for some $0<\alpha<1$. It has
characteristic size $\delta$. The point-like source is modelled by an ideal
electric dipole of centre $z^{*}\in\mathbb{R}^{2}\setminus{\overline{D}}$ and
direction $\mathbf{p^{*}}\in\mathbb{S}^{1}$. We denote the permittivity of the
plasmonic particle by $\varepsilon_{D}$. The permittivity of the background
medium is denoted by $\varepsilon_{m}$. In other words, the permittivity
distribution $\varepsilon$ is given by

$\varepsilon=\varepsilon_{D}\chi\left(D\right)+\varepsilon_{m}\chi\left(\mathbb{R}^{2}\setminus\overline{D}\right).$

Note that $\varepsilon_{D}$ depends on the frequency $\omega$. The model used
for $\varepsilon_{D}$ depends on the type of metal. It is nonrestrictive to
assume that $\varepsilon_{D}$ is described by a Drude model [28]:

$\displaystyle\varepsilon_{D}(\omega)=\varepsilon_{0}\left(1-\frac{\omega_{p}^{2}}{\omega\left(\omega+i\tau^{-1}\right)}\right),$

where $\omega_{p}$ is the plasma frequency of the bulk material and $\tau$ is
the bulk electron relaxation time. Any similar model can be used [26]. The
permeability $\mu$ is constant everywhere, $\mu=\mu_{m}$. We define the
wavenumbers $k_{D}=\omega\sqrt{\varepsilon_{D}\mu_{m}}$ and
$k_{m}=\omega\sqrt{\varepsilon_{m}\mu_{m}}$. We denote by $z_{D}$ the centre
of the particle $D$. Up to a translation it is nonrestrictive to assume that
$z_{D}=0$.

Consider $\Omega\supset D$ a domain of characteristic size $R\gg{k_{m}}^{-1}$
and $d(D,\partial\Omega)\gg{k_{m}}^{-1}$. Our goal is to reconstruct
$\mathbf{p^{*}}$ and $z^{*}$ from the knowledge of far-field data, _i.e._ ,
the knowledge on $\partial\Omega$ of $u$ the solution of

$\displaystyle\nabla\cdot\left(\frac{1}{\varepsilon(x)}\nabla
u\right)+\omega^{2}\mu\ u=\mathbf{p^{*}}\cdot\nabla\delta_{z^{*}}\qquad\text{
in }\mathbb{R}^{2},$ (1)

satisfying the outgoing Sommerfeld radiation condition

$\left|\frac{\partial
u}{\partial|x|}-ik_{m}u\right|=\mathcal{O}\left(|x|^{-3/2}\right),\qquad\mbox{as
}|x|\rightarrow\infty,$

uniformly in $x/|x|$, for $\Re{k_{m}}>0$. For simplicity we assume that
$\Omega=B_{R}$ the disk centred at the origin of radius $R$. The setting is
sketched in Figure 1.

Figure 1: Setting.

######  Proposition 2.1.

If $\Im\varepsilon_{D}\neq 0$, then equation (1) is well posed.

###### Proof.

This is a classical result. A recent treatment of this exact case can be found
in [6].

∎

######  Remark 2.1.

The condition $\Im\varepsilon_{D}\neq 0$ is not sharp, but it is sufficient.
The full treatment of the case $\Im\varepsilon_{D}\rightarrow 0$ is done in
[25].

###  2.2 Relevance of the scalar model

The scalar Helmholtz equation is the correct way to describe the
electromagnetic field propagation in a two-dimensional environment [22, Remark
2.1]. It is relevant to model the propagation of an electromagnetic wave in a
medium that is invariant along a coordinate.

Two cases can be modelled. One case is when the electric field is
perpendicular to the plane of propagation. It has the form $(0,0,E_{z})$,
while the magnetic field has the form $(H_{x},H_{y},0)$. This case is often
called the _transverse electric_ case. The other case is when the magnetic
field is transverse to the plane of propagation. It is called the _transverse
magnetic_ case and the magnetic field and electric field are $(0,0,H_{z})$ and
$(E_{x},E_{y},0)$, respectively.

For a non-magnetic medium described by a space-dependent permittivity
distribution $\varepsilon(x)$ with $x\in\mathbb{R}^{2}$ and a constant
magnetic permeability, the transverse magnetic case is the most mathematically
interesting case, as the spectral properties of the transmission problem for a
small obstacle (of characteristic size small compared to the wavelength) are
very similar to the ones of the Maxwell transmission problem. Indeed, in the
_static regime_ , the field $\mathbf{E}$ can be written as
$\mathbf{E}:=\nabla\tilde{u}:=(E_{x},E_{y},0):=(\nabla u,0)$. Then the
Lippman-Schwinger equation associated with problem (1), _i.e._ , the equation
solved by the electric field inside the particle $\nabla u=(E_{x},E_{y})$ is
[5]:

$\displaystyle\eta\nabla u(x)-\nabla\int_{D}\nabla\Gamma^{0}(x-y)\cdot\nabla
u(y)\mathrm{d}y=\frac{1}{1-\varepsilon_{r}}D^{2}\Gamma^{0}(x-z^{*})\mathbf{p}^{*},\qquad
x\in D,$ (2)

where $\Gamma^{k}$, $k\geq 0$, is a fundamental outgoing solution to the
Helmholtz operator $\Delta+k^{2}$ given by

$\Gamma^{k}(x)=\begin{dcases}\frac{1}{2\pi}\ln|x|&\mbox{if }k=0,\\\
-\frac{i}{4}H_{0}^{(1)}{(k|x|)}&\mbox{if }k>0,\end{dcases}$ (3)

and $\eta=\varepsilon_{m}/(\varepsilon_{m}-\varepsilon_{D})$ is a contrast
parameter depending on the relative permittivity of the obstacle with respect
to the surrounding medium.

In the three-dimensional case, the integral equation solved by the electric
field $\mathbf{E}$ inside the particle is very similar (see for instance [13]
or [5, equation (3.1)])

$\displaystyle\eta\mathbf{E}(x)-\nabla\int_{D}\nabla\tilde{\Gamma}^{0}(x-y)\cdot\mathbf{E}(y)\mathrm{d}y=\frac{1}{1-\varepsilon_{r}}D^{2}\tilde{\Gamma}^{0}(x-z^{*})\mathbf{p}^{*},\qquad
x\in D,$ (4)

where $\tilde{\Gamma}^{k}$ is the fundamental solution of the Helmholtz
equation in three dimensions.

We can notice a strong similarity between equations (2) and (4), and therefore
the scalar model is a good starting point. The full static and low frequency
spectral analysis of equation (4) is performed in [10], while the Maxwell case
can be found in [7]. In this paper we only focus on the two-dimensional case.

###  2.3 Limits of the model

The mathematical results presented here hold when the size of the nanoparticle
becomes very small compared to the wavelength or when the dipole source gets
very close to the particle. Nevertheless, the model used to describe the
interaction between the fluorescent molecule and the nanoparticle is incorrect
when the different scales in the problem become too small, and fails to
properly describe the physics of the problem, especially any quantum effect
that may appear. We can mention the _quenching_ phenomenon that occurs when
the molecule gets very close to the nanoparticle (a few nanometers) [12, 14].
When the nanoparticle becomes very small (of characteristic size less than
$10$nm), the size of the fluorescent molecule can become comparable to the
size of the nanoparticle itself. Therefore it is not justified to model the
fluorescent molecule as a point dipole, and a more sophisticated (quantum)
model has to be considered [24]. However, as mentioned in subsection 1.2, the
novelty of the method proposed here is that it can be applied to _larger_
nanoparticles, of size smaller but comparable to the wavelength of the
electromagnetic field emitted by the fluorescent molecule.

##  3 Modal Expansion of the solution

In the following section, we recall the results from [10] that are needed for
the localisation procedure.

We will use standard _layer potentials_ notations. The definitions and
symmetrisation process of the Neumann-Poincaré operator on which the theorem
relies are recalled in appendix A. We also refer to the book [4] and
references therein.

######  Definition 3.1.

We define the _modes_ of the system by

$e^{\omega}_{n}(x):=\mathcal{S}_{D}^{k_{m}}\left[\varphi_{n}\right](x),\qquad
x\in\mathbb{R}^{2}\setminus\overline{D},\leavevmode\nobreak\
n\in\mathbb{N}^{*}.$

$\mathcal{S}_{D}^{k_{m}}$ is the single-layer potential and $\varphi_{n}$ are
the eigenvectors of the Neumann-Poincaré operator $\mathcal{K}_{D}^{*}$.

######  Theorem 3.1.

The unique solution $u$ of (1) can be written as:

$\displaystyle
u(x)=\nabla\Gamma^{k_{m}}(x,z^{*})\cdot\mathbf{p^{*}}+\sum_{n=1}^{N}\alpha_{n}e^{\omega}_{n}(x)+\mathcal{E}_{N}(x),\qquad
x\in\mathbb{R}^{2}\setminus\overline{D},N\in\mathbb{N}^{*},$

where $\alpha_{n}$ are coupling coefficients between the dipole source and the
modes that depend on $z^{*}$, $\mathbf{p}^{*}$ and $\omega$.
$\mathcal{E}_{N}(x)$ is an error term that depends on the position of the
source, the shape of the particle $D$, the number of modes considered, and the
observation position $x$. The coupling coefficients $\alpha_{n}$ have a
(rather lengthy) explicit formula that is given in the proof.

###### Proof.

The theorem was proved in [10, proposition 5]. We recall here for the sake of
completeness the expression of the coupling coefficients. Define $F$ as

$F=-\frac{1}{\varepsilon_{m}}\bm{\nu}(x)^{\top}\mathrm{D}^{2}\Gamma^{k_{m}}(x,z^{*})\mathbf{p^{*}}-\frac{1}{\varepsilon_{D}}\left(\frac{1}{2}I-\mathcal{K}^{k_{D},*}_{D}\right)\left(\mathcal{S}^{k_{D}}_{D}\right)^{-1}[\nabla\Gamma^{k_{m}}(x,z^{*})\cdot\mathbf{p^{*}}],\quad
x\in\partial D,$

where $\bm{\nu}(x)$ is the normal vector on $\partial D$. We denote by
$\lambda_{n}$ the eigenvalues of the Neumann-Poincaré operator on $D$ (see
appendix A) and introduce the contrast parameter:

$\lambda(\omega):=\dfrac{\varepsilon_{m}+\varepsilon_{D}(\omega)}{2\left(\varepsilon_{m}-\varepsilon_{D}(\omega)\right)},$

as well as

$\tau_{n}(\omega)=\left(\frac{1}{\varepsilon_{D}}-\frac{1}{\varepsilon_{m}}\right)\left(\lambda(\omega)-\lambda_{n}\right)+\left(\omega\delta
c^{-1}\right)^{2}\log{\left(\omega\delta c^{-1}\right)}\tau_{n,1},$

where $\tau_{n,1}$ is a perturbative spectral parameter defined in definition
A.1. Then [10, proposition 5] gives us the following expansion for $u$:

$\displaystyle
u(x)=\nabla\Gamma^{k_{m}}(x,z^{*})\cdot\mathbf{p^{*}}+\sum_{n=1}^{N}\frac{\left\langle
F,\varphi_{n}\right\rangle_{\mathcal{H}^{*}}}{\tau_{n}(\omega)}\mathcal{S}_{D}^{k_{m}}\left[\varphi_{n}\right](x)+\mathcal{E}_{N}(x),\qquad
x\in\mathbb{R}^{2}\setminus\overline{D},N\in\mathbb{N}^{*}.$

Setting $\alpha_{n}:=\left\langle
F,\varphi_{n}\right\rangle_{\mathcal{H}^{*}}/(\tau_{n}(\omega))$ we get the
result. ∎

######  Remark 3.1.

The coupling coefficients

$\displaystyle\alpha_{n}=\frac{\left\langle
F,\varphi_{n}\right\rangle_{\mathcal{H}^{*}}}{\tau_{n}(\omega)}$

can become very large if $\tau_{n}(\omega)\ll 1$. This will happen if $\omega$
is close to a plasmonic resonance frequency. It is then easy to see that in
the case that $\left\langle F,\varphi_{n}\right\rangle_{H^{*}}\not=0$ and
$\tau_{n}(\omega)\ll 1$, the field scattered by the nanoparticle can become
very large and contribute to the far field as much as the radiating dipole,
making it very difficult to localise the dipole, as the imaging system will be
_blinded_ by the nanoparticle.

######  Remark 3.2.

The number of modes that have to be considered to obtain a good approximation
of the field is a complicated question because the modal approximation is not
uniform, neither with respect to the observation point $x$, nor with the
frequency $\omega$ or the parameters of the source ($\mathbf{p}^{*},z^{*}$).
However, in practice, only a few modes are needed to correctly describe the
scattered field. This is due to the fact that:

  1. 1.

the _coupling coefficient_ of the mode $n$ to the source

$\displaystyle\alpha_{n}=\frac{\left\langle
F,\varphi_{n}\right\rangle_{\mathcal{H}^{*}}}{\tau_{n}(\omega)}$

have superpolynomial decay with respect to $n$ for a fixed source and
frequency [10, proposition 3],

  2. 2.

high order modes do not contribute to the far field. It has been shown in [9]
that, in the static case ($\omega=0$), the contribution of the modes outside
the particle decay exponentially with the order of the mode. It has not been
demonstrated in the dynamic case $\omega\neq 0$ but there are numerical
indications that the result probably holds for $\omega\neq 0$.

So even though it is very difficult to get a quantitative estimate on the
error term in theorem 3.1, in practical situations it is a good approximation.

##  4 Localisation

We now introduce a method that recovers a dipole source position and
orientation in the presence of a resonant nanoparticle at a low computational
cost. It is based on the modification of a back-propagation/holography
procedure, but the principle of the method works with any linear imaging
functional.

In the absence of a resonant structure, super-localisation is possible through
the modelling of a _point spread function_ (PSF) of the imaging system, that
maps a single point-like source to an image. Then the location of the point
source is obtained by finding the centre of the PSF, which can be done with a
precision that depends only on the signal-to-noise ratio. In general, it is
achieved with a least-square procedure, i.e. by fitting the centre of the PSF
on the image.

In the presence of the resonant structure, the image will not be the PSF, but
the superposition of the PSF at the location of the dipole with the image of
the nanoparticle. The current state-of-the-art method to recover the correct
position of an emitter near a nano-particle requires solving the forward
problem (the full-wave Maxwell equations) for each position of the source
[16]. This is very costly in terms of computation, making it impossible to
apply in practical situations where thousands of isolated single emitters have
to be localised quickly.

The idea behind our method is that even though the field scattered by the
nanoparticle cannot be known in advance (as it depends on both the position
and the orientation of the dipole), it can be expressed as a linear
combination of pre-computable quantities (the modes). This suggests that the
distortion of the PSF due to the presence of the nanoparticle is simply a
linear combination of pre-computable artefacts. Therefore, instead of
performing the least-square procedure on the position of the centre of the
PSF, we optimise on both the centre of the PSF and the coefficients of the
linear combination of artefacts.

###  4.1 Back-propagation in the absence of a resonant structure

Super-localisation of the fluorescent molecule is obtained by back-propagating
$u(x)$, for $|x|$ large enough, and maximising the back-propagation function.
The precision of the localisation is limited by the signal-to-noise ratio
(SNR).

We introduce the following vector-valued imaging functional.

######  Definition 4.1.

$\displaystyle\mathcal{I}(z)=\int_{\partial
B_{R}}\overline{\nabla_{z}\Gamma^{k_{m}}(x,z)}u(x)\mathrm{d}\sigma(x),\quad
z\in\Omega,$

where the imaging domain $\Omega$ is such that

$\displaystyle D\Subset\Omega\Subset B(0,R).$

We then include a lemma from [3] with its proof for sake of completeness.

######  Lemma 4.1.

Let

$\displaystyle\mathcal{R}(y,z)=\int_{\partial
B_{R}}\overline{\nabla_{y}\Gamma^{k_{m}}(x,y)}\nabla_{z}\Gamma^{k_{m}}(x,z)^{\top}\mathrm{d}\sigma(x),\qquad(y,z)\in\mathbb{R}^{2}.$

Then

$\displaystyle\lim_{R\rightarrow\infty}\mathcal{R}(z,z)=\frac{k_{m}}{8}I_{2},$

where $I_{2}$ is the identity matrix.

###### Proof.

We start from the Helmholtz-Kirchhoff theorem [2]

$\displaystyle\lim_{R\rightarrow\infty}\mathcal{R}(y,z)=-\frac{1}{k_{m}}\nabla_{y}\nabla_{z}\Im\,\Gamma^{k_{m}}(y,z),$

where $\Im$ is the imaginary part. We can then compute an approximation of
$\mathcal{R}$:

$\frac{1}{k_{m}}\nabla_{y}\nabla_{z}\Im\,\Gamma^{k_{m}}(y,z)=\frac{1}{4}\Bigg{[}k_{m}J_{0}\left(k_{m}|y-z|\right)\frac{\left(y-z\right)\left(y-z\right)^{\top}}{|y-z|^{2}}\\\
-2\frac{J_{1}\left(k_{m}|y-z|\right)}{|y-z|}\frac{\left(y-z\right)\left(y-z\right)^{\top}}{|y-z|^{2}}+\frac{J_{1}\left(k_{m}|y-z|\right)}{|y-z|}I_{2}\Bigg{]}.$

One can see that $\mathcal{R}(y,z)$ decays as $|y-z|^{-1/2}$ and the imaging
functional has a peak at $y=z$. Evaluating at $y=z$ gives the result. ∎

It is easy to see from lemma 4.1 that in the absence of the resonant
structure, $\mathcal{I}(z)$ would have an amplitude peak at $z=z^{*}$ the
position of the dipole source and its orientation would be parallel to
$\mathbf{p^{*}}$. To get $z^{*}$ and $\mathbf{p^{*}}$ one needs to solve the
following optimisation problem in $\mathbb{R}^{4}$:

$\displaystyle(z^{*},\mathbf{p^{*}})=\mathrm{argmin}_{z,\mathbf{p}}\left\|\mathcal{I}(\cdot)-\mathcal{R}(\cdot,z)\mathbf{p}\right\|_{L^{2}(\Omega)}.$
(5)

###  4.2 Back-propagation in the presence of the resonant structure

######  Proposition 4.1.

In the presence of the nanoparticle, the imaging functional has the form:

$\displaystyle\mathcal{I}(z)=\mathcal{R}(z,z^{*})\mathbf{p^{*}}+\sum_{n=1}^{N}\alpha_{n}I_{e_{n}}(z)+I_{\mathcal{E}_{N}}(z),$

where $\alpha_{n}$ are the (unknown) coupling coefficients defined in theorem
3.1, $I_{e_{n}}$ the (known) image by the imaging functional of the mode $n$
and $I_{\mathcal{E}_{n}}(z)$ the (unknown) image of the error term.

###### Proof.

The result is a direct consequence of theorem 3.1 and the linearity of the
imaging functional. Let

$\displaystyle I_{e_{n}}(z)=\int_{\partial
B_{R}}\overline{\nabla_{z}\Gamma^{k_{m}}(x,z)}e_{n}(x)\mathrm{d}\sigma(x),$

and

$\displaystyle I_{\mathcal{E}_{N}}(z)=\int_{\partial
B_{R}}\overline{\nabla_{z}\Gamma^{k_{m}}(x,z)}\mathcal{E}_{N}(x)\mathrm{d}\sigma(x).$

We obtain the result by applying lemma 4.1. ∎

To recover $z^{*}$ and $\mathbf{p^{*}}$ in the presence of the resonant
structure, we now solve an optimisation problem in $\mathbb{R}^{4+N}$:

$\displaystyle\left(z^{*},\mathbf{p^{*}},(\alpha_{n}^{*})_{n=1}^{N}\right)=\mathrm{argmin}_{z,\mathbf{p},(\alpha_{n})_{n=1}^{N}}\left\|\mathcal{I}(\cdot)-\mathcal{R}(\cdot,z)\mathbf{p}-\sum_{n=1}^{N}\alpha_{n}I_{e_{n}}(z)\right\|^{2}_{L^{2}(\widetilde{\Omega})}.$
(6)

##  5 Numerical simulations

In this section we illustrate the method’s efficiency with numerical
simulations. We perform a virtual experiment where a dipole source is placed
near a metallic nanoparticle whose permittivity is described by a Drude model.
We consider nanoparticles with various shapes and subject to noise to
illustrate the robustness of the method.

###  5.1 Geometry and physical parameters

Throughout this section, we consider the three domains sketched on Figure 2 to
illustrate our results. The five-petal flower (a) is defined by
$\varrho=\delta(2+0.6\cos(5\theta))$ in polar coordinates. The rounded diamond
(b) is defined by the parametric curve
$\zeta(\theta)=2\delta\left(e^{i\theta}+0.066e^{-3i\theta}\right)$, for
$\theta\in[0,2\pi]$. The narrow ellipse (c) semi-axes are on the $x_{1}$\- and
$x_{2}$\- axes and are of length $a=1\delta$ and $b=5\delta$, respectively.

Figure 2: Sketch of the three reference domains: the five-petal flower (a),
the rounded diamond (b) and the narrow ellipse (c).

We recall the permittivity inside the particle is described by a Drude model

$\displaystyle\varepsilon_{D}(\omega)=\varepsilon_{0}\left(1-\frac{\omega_{p}^{2}}{\omega\left(\omega+i\tau^{-1}\right)}\right).$

The physical parameters are chosen as $\omega_{p}=2\cdot 10^{15}$ Hz,
$T=10^{-14}$s, $\varepsilon_{0}=8.854187128\cdot 10^{-12}$ Fm-1,
$\mu_{0}=4\pi\cdot 10^{-7}$ Hm-1 and $\delta=10^{-8}$ m.

###  5.2 Modes and resonances

For each shape, we compute a list of $N$ plasmonic resonant frequencies
$(\omega_{n})_{n\in[1,N]}$ and the associated modes. We give an illustration
for the rounded diamond.

Figure 3: First $6$ modes for the rounded diamond, corresponding to
$\omega=1.5050\cdot 10^{15}$ Hz.

###  5.3 Imaging experiment

For each shape we pick a resonant frequency $\omega_{n}$ and place an
oscillating dipole at frequency $\omega_{n}$ near the nanoparticle. We choose
the frequency that corresponds to a mode that has a dipole radiation pattern,
to maximise the effect of the nanoparticle’s presence on the far-field image
and make the localisation of the dipole source as hard as possible.

The scattered field (solution of (1)) is then computed with a system of
boundary integral equations (that uses layer potential techniques) detailed in
appendix A. The boundary of the nanoparticle is discretised with $2^{8}$
points.

The field (amplitude and phase) is then collected in the far field on a disk
of radius 3000$\delta$. The measured field is then back-propagated using the
imaging functional introduced in definition 4.1.

###  5.4 Localisation by optimisation and mirage

The localisation (without correction) of the particle is done by fitting a
focal spot on the imaging functional, as described in equation (5). The
optimisation is done using the _fmincon_ MATLAB® function. We constrain the
search to a $300\times 300\text{nm}^{2}$ square. The position of the initial
guess does not influence the result.

For the corrected functional, we use the _fmincon_ function to search for a
solution of equation (6). The number of modes to use for the correction will
be discussed later in this section.

We show, on Figure 4, examples of the results of the optimisation procedure
for the three different domains.

Figure 4: Mirage and optimisation on the three domains: for the flower (a),
the localisation error is $29$nm before correction and $8$nm. For the diamond
(b), the error is $37$nm before correction and $5.6$nm. The dipole direction
is $(-1,1)$. For the ellipse (c), the localisation error is $71$nm before
correction and $8$nm after and the dipole direction is $(0,1)$.

###  5.5 Sensitivity analysis

On Figure 5 we plot the error on the localisation (a) and orientation (b) of
the dipole with respect to the distance between the emitter and the
nanoparticle. We initially place the dipole at $z^{*}=[18.65,16.65]$nm and fix
its orientation ($\mathbf{p}^{*}=(-1,1)$). We gradually increase the distance
to the nanoparticle by translating the dipole along a radial component of a
disk centred at the origin. As expected, the uncorrected functional does not
behave well when the source is close to the nanoparticle. Note that the error
is still one order of magnitude below the usual diffraction limit. The
localisation error on the corrected functional is one order of magnitude lower
than the uncorrected one. The correction is a massive improvement if the
emitter is close enough to excite the nanoparticle. When the emitter is far
away and the coupling is weaker, both procedures have similar performances.

Figure 5: (a) Absolute error on the position $z^{*}$ in nanometers with
(orange crosses) and without (purple circles) correction against the distance
between the emitter and the rounded diamond boundary. (b) Absolute error on
the orientation $\mathbf{p}^{*}$ in degrees with (orange crosses) and without
(blue triangles) correction against the distance between the emitter and the
rounded diamond boundary. The initial position is $z^{*}=[18.65,16.65]$nm. The
dipole direction is $(-1,1)$.

###  5.6 Stability with respect to measurement noise

We compute the imaging functional with a set of data perturbed by Gaussian
white noise. We consider the case in which the measured field $u(x)$ is
corrupted by an additive noise $\nu(x)$, $x\in\partial B_{R}$. We assume that
$\partial B_{R}$ is covered uniformly with sensors and that the additive
noises $\nu(x)$ have independent and identically distributed
$\mathcal{N}(0,\sigma_{\text{noise}})$ entries. Hence, the entries of $\nu(x)$
are independent Gaussian random variables with mean zero and variance

$\sigma_{\text{noise}}=\sigma_{0}\|u\|_{F}/\sqrt{N},$

where $||\leavevmode\nobreak\ ||_{F}$ is the Frobenius norm, $\sigma_{0}$ the
noise percentage and $N$ the number of sensors.

For each noise level, we average the results over $100$ realisations. Figure 6
presents the results of computational experiments in the case in which the
nanoparticle is the diamond. It shows that the corrected imaging functional
performs well at high levels of noise, which means that the corrected imaging
functional is robust with respect to additive measurement noise.

Figure 6: Boxplots showing the error on the position (a) and orientation (p)
for noise levels $\sigma_{0}$ varying from $40\%$ to $150\%$. The boxplots
show the median, 25- and 75-quartiles, the notches correspond to the 95
confidence interval for 100 realisations. The outliers are shown as black
symbols.

###  5.7 Number of modes

For the numerical experiments of this section we have assumed to know in
advance the pre-computable artefacts that describe the distortion of the PSF
due to the presence of the nanoparticle, _i.e._ , the modes excited by the
fluorescent molecule. In practice, this may not the case. However, we know
that only a few modes are excited and contribute significantly to approximate
the scattered field. One can ask whether optimising over excited _and_
unexcited modes will give results as good as the ones obtained by optimising
over the excited modes only.

For these simulations we consider the rounded diamond. The dipole position is
$[18.65,16.65]$nm. The initial direction is $(-1,1)$. In this case, we can
check numerically that only modes 4 and 6 suffice to approximate the field.

On Table 1, we show the error on the position and orientation for the first
$N$ modes. As expected, the error becomes small as soon as $N\geq 6$, _i.e._ ,
modes 4 and 6 are considered, regardless of the number of unexcited modes used
for the optimisation procedure.

$N\leavevmode\nobreak\ \text{modes}$ | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---
$||\bf{p}_{N}-\bf{p}^{*}||\leavevmode\nobreak\ \text{(in deg)}$ | 0.57 | 89.2 | 13.9 | 45.0 | 0.57 | 0.30
$||z_{N}-z^{*}||\leavevmode\nobreak\ \text{(in nm)}$ | 116 | 145 | 116 | 166 | 4.74 | 4.22
Table 1: Error on the position and orientation for the first $N$ modes,
$N=2..7$.

##  6 Concluding remarks

In this paper we have shown a new method to recover the localization of a
single isolated source in a resonant environment with a low computational
cost. Using the results of [10, 7] and pre-computing the modes of the resonant
structure, we transform the _inverse source problem for the wave equation_
into a simple minimization in $\mathbb{R}^{d}$ with $d$ being the number of
modes to consider, usually less than $10$. We can recover the position of the
source with a precision that is two orders of magnitude below the usual
diffraction limit. Since most of the computational cost is the pre-computation
of the modes, the position of many single sources can be recovered very
quickly, as is needed in real _PALM_ or _STORM_ experiments, where the
position of thousands of emitters are needed to produce an image. This paper
acts like a proof of principle for this method. In a forthcoming paper, we
will adapt the numerics for the $3D$ Maxwell case that is more suited for the
applications, and consider a more complex structured environment as well as
amplitude only measurements.

##  Appendix A Layer potentials

The results presented in this section are all classical and can be found in
the books [4, 2].

###  A.1 Operator

For $k\geq 0$, a fundamental solution to the Helmholtz operator $\Delta+k^{2}$
is given by

$\Gamma^{k}(x)=\begin{dcases}\frac{1}{2\pi}\ln|x|&\mbox{if }k=0,\\\
-\frac{i}{4}H_{0}^{(1)}{(k|x|)}&\mbox{if }k>0,\end{dcases}$

for $x\neq 0$, where $H_{0}^{(1)}$ is the well-known Hankel function of the
first kind and order zero. Let $\mathcal{S}^{k}_{D}$ be the single layer
potential, defined by

$\mathcal{S}^{k}_{D}[\phi](x)=\int_{\partial
D}\Gamma^{k}(x,y)\phi(y)\mathrm{d}\sigma(y),\qquad x\in\mathbb{R}^{2},$ (7)

for $\phi\in L^{2}(\partial D)$. We also define the Neumann-Poincaré operator
by

$\mathcal{K}^{k,*}_{D}[\phi](x)=\int_{\partial
D}\frac{\partial\Gamma^{k}(x,y)}{\partial\nu(x)}\phi(y)\mathrm{d}\sigma(y),\qquad
x\in\partial D,$ (8)

for $\phi\in L^{2}(\partial D)$, where $\partial/\partial\nu(x)$ denotes the
outward normal derivative at $x\in\partial D$. When $k=0$, we omit the
superscript and write $\mathcal{S}_{D}$ and $\mathcal{K}^{*}_{D}$ for
simplicity. Let $\left\langle\cdot,\cdot\right\rangle_{-1/2,1/2}$ be the
duality pairing between $H^{-\frac{1}{2}}(\partial D)$ and
$H^{\frac{1}{2}}(\partial D)$, where $H^{\frac{1}{2}}(\partial D)$ is the
Sobolev space of order 1/2. The single layer potential is, in general, not
invertible in $L^{2}(\partial D)$. Let us introduce

$\widetilde{\mathcal{S}}_{D}[v]=\begin{cases}\mathcal{S}_{D}[v]&\text{if
}\langle v,\chi(\partial D)\rangle_{-\frac{1}{2},\frac{1}{2}}=0,\\\
-\chi(\partial D)&\text{if }v=\varphi_{0},\end{cases}$

with $\varphi_{0}$ being the unique (in the case of a single particle)
eigenfunction of $\mathcal{K}^{*}_{D}$ associated with eigenvalue $1/2$ such
that $\langle\varphi_{0},\chi(\partial
D)\rangle_{-\frac{1}{2},\frac{1}{2}}=1$. The operator
$\mathcal{K}^{*}_{D}:H^{-1/2}(\partial D)\rightarrow H^{-1/2}(\partial D)$ is
compact and the following Plemelj’s symmetrisation principle identity (also
known as Calderón) holds on $H^{-1/2}(\partial D)$:

$\widetilde{\mathcal{S}}_{D}\mathcal{K}^{*}_{D}=\mathcal{K}_{D}\widetilde{\mathcal{S}}_{D}.$

Let $\mathcal{H}^{*}(\partial D)$ be the space $H^{-1/2}(\partial D)$ equipped
with the following inner product:

$<u,v>_{\mathcal{H}^{*}(\partial
D)}=-<\widetilde{\mathcal{S}}_{D}[v],u>_{1/2,-1/2}.$

The Neumann-Poincaré operator $\mathcal{K}^{*}_{D}$ is self-adjoint in the
Hilbert space $\mathcal{H}^{*}(\partial D)$. Let
$(\lambda_{n},\varphi_{n})_{n\in\mathbb{N}}$ be the eigenvalue and normalised
eigenfunction pair of $\mathcal{K}^{*}_{D}$ in $\mathcal{H}^{*}(\partial D)$.
Then $\lambda_{0}=1/2$, $-1/2<\lambda_{n}<1/2$ for $n\geq 1$ and
$\lambda_{n}\rightarrow 0$ as $n\rightarrow+\infty$.

###  A.2 Solution for the wave equation

Let $H^{1/2}(\partial D)$ be the usual Sobolev space and let
$H^{-1/2}(\partial D)$ be its dual space with respect to the duality pairing
$\left\langle\cdot,\cdot\right\rangle_{\frac{1}{2},-\frac{1}{2}}$. The field
$u$ can be represented using the single-layer potentials
$\mathcal{S}^{k_{c}}_{D}$ and $\mathcal{S}^{k_{m}}_{D}$, introduced in
equation (7), as follows:

$u(x)=\begin{dcases}\mathcal{S}^{k_{c}}_{D}[\Phi](x),&x\in D,\\\
u^{\text{in}}(x)+\mathcal{S}^{k_{m}}_{D}[\Psi](x),&x\in\mathbb{R}^{d}\setminus\overline{D},\end{dcases}$
(9)

where the pair $(\Phi,\Psi)\in H^{-\frac{1}{2}}(\partial D)\times
H^{-\frac{1}{2}}(\partial D)$ is the unique solution to

$\begin{dcases}\mathcal{S}^{k_{m}}_{D}[\Psi](x)-\mathcal{S}^{k_{c}}_{D}[\Phi](x)=F_{1},&x\in\partial
D,\\\
\frac{1}{\varepsilon_{m}}\left(\frac{1}{2}I+\mathcal{K}^{k_{m},*}_{D}\right)[\Psi](x)+\frac{1}{\varepsilon_{c}}\left(\frac{1}{2}I-\mathcal{K}^{k_{c},*}_{D}\right)[\Phi](x)=F_{2},&x\in\partial
D,\end{dcases}$ (10)

and

$F_{1}=-u^{\text{in}}(x),\qquad F_{2}=-\frac{1}{\varepsilon_{m}}\frac{\partial
u^{\text{in}}(x)}{\partial\nu},\qquad x\in\partial D,$

where $\mathcal{K}^{k_{m},*}_{D}$ is the Neumann-Poincaré operator introduced
in equation (8) and
$u^{\text{in}}=\mathbf{p}^{*}\cdot\nabla\Gamma^{k}(\cdot,z^{*})$.

###  A.3 Definition of $\tau_{n,1}$

######  Lemma A.1.

For $k$ small enough, the two-dimensional boundary operator
$\widehat{\mathcal{S}}_{D}^{k}:\mathcal{H}^{*}(\partial
D)\rightarrow\mathcal{H}^{*}(\partial D)$ defined as

$\widehat{\mathcal{S}}_{D}^{k}[\phi](x)=\mathcal{S}_{D}^{0}[\phi](x)+\eta_{k}\int_{\partial
D}\phi(y)\mathrm{d}\sigma(y),$ (11)

is invertible and

$\left(\widehat{\mathcal{S}}_{D}^{k}\right)^{-1}=\widetilde{\mathcal{S}}_{D}^{-1}-\left\langle\widetilde{\mathcal{S}}_{D}^{-1}[\cdot],\varphi_{0}\right\rangle_{\mathcal{H}^{*}(\partial
D)}\varphi_{0}-\mathcal{U}_{k},$ (12)

where

$\mathcal{U}_{k}=\dfrac{\left\langle\widetilde{\mathcal{S}}_{D}^{-1}[\cdot],\varphi_{0}\right\rangle_{\mathcal{H}^{*}(\partial
D)}}{\mathcal{S}_{D}[\varphi_{0}]+\eta_{k}}\varphi_{0}$

and $\eta_{k}=(1/2\pi)(\log{k}+\gamma-\log{2})-i/4,$ with the constant
$\gamma$ being the Euler constant. Note that
$\mathcal{U}_{k}=\mathcal{O}(1/\log k)$.

######  Definition A.1.

$\tau_{n,1}=\left\langle\mathcal{A}_{D,1}\varphi_{n},\varphi_{n}\right\rangle_{\mathcal{H}^{*}(\partial
D)}$

with

$\mathcal{A}_{D,1}=\frac{1}{\varepsilon_{m}}\mathcal{K}^{(1)}_{D,1}(I-\mathcal{P}_{\mathcal{H}^{*}_{0}})+\left(\frac{1}{2}I-\mathcal{K}^{*}_{D}\right)\widetilde{\mathcal{S}}_{D}^{-1}\mathcal{S}_{D,1}^{(1)}\left(\frac{1}{\varepsilon_{D}}I-\frac{1}{\varepsilon_{m}}\mathcal{P}_{\mathcal{H}^{*}_{0}}\right),$

## References

  * [1] H. Ammari, Y. Deng, and P. Millien, Surface plasmon resonance of nanoparticles and applications in imaging, Archive for Rational Mechanics and Analysis, 220 (2016), pp. 109–153.
  * [2] H. Ammari, J. Garnier, W. Jing, H. Kang, M. Lim, K. Sølna, and H. Wang, Mathematical and statistical methods for multistatic imaging, vol. 2098, Springer, 2013.
  * [3] H. Ammari, J. Garnier, and P. Millien, Backpropagation imaging in nonlinear harmonic holography in the presence of measurement and medium noises, SIAM Journal on Imaging Sciences, 7 (2014), pp. 239–276.
  * [4] H. Ammari, H. Kang, and H. Lee, Layer Potential Techniques in Spectral Analysis, vol. 153 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, 2009.
  * [5] H. Ammari and P. Millien, Shape and size dependence of dipolar plasmonic resonance of nanoparticles, Journal de Mathématiques Pures et Appliquées, 129 (2019), pp. 242 – 265.
  * [6] H. Ammari, P. Millien, M. Ruiz, and H. Zhang, Mathematical analysis of plasmonic nanoparticles: The scalar case, Arch. Ration. Mech. Anal., 224 (2017), pp. 597–658.
  * [7] H. Ammari, P. Millien, and A. L. Vanel, Modal approximation for strictly convex plasmonic resonators in the time domain: The maxwell’s equations, Journal of Differential Equations, 309 (2022), pp. 676–703.
  * [8] K. Ando and H. Kang, Analysis of plasmon resonance on smooth domains using spectral properties of the Neumann-Poincaré operator, Journal of Mathematical Analysis and Applications, 435 (2016), pp. 162–178.
  * [9] K. Ando, H. Kang, Y. Miyanishi, and T. Nakazawa, Surface localization of plasmons in three dimensions and convexity, SIAM Journal on Applied Mathematics, 81 (2021), pp. 1020–1033.
  * [10] L. Baldassari, P. Millien, and A. L. Vanel, Modal approximation for plasmonic resonators in the time domain: the scalar case, Partial Differential Equations and Applications, 2 (2021), pp. 1–40.
  * [11] D. Bouchet, J. Scholler, G. Blanquer, Y. De Wilde, I. Izeddin, and V. Krachmalnicoff, Probing near-field light–matter interactions with single-molecule lifetime imaging, Optica, 6 (2019), pp. 135–138.
  * [12] E. Castanié, M. Boffety, and R. Carminati, Fluorescence quenching by a metal nanoparticle in the extreme near-field regime, Optics letters, 35 (2010), pp. 291–293.
  * [13] M. Costabel, E. Darrigrand, and H. Sakly, The essential spectrum of the volume integral operator in electromagnetic scattering by a homogeneous body, Comptes Rendus Mathematique, 350 (2012), pp. 193–197.
  * [14] A. Delga, J. Feist, J. Bravo-Abad, and F. Garcia-Vidal, Quantum emitters near a metal nanoparticle: strong coupling and quenching, Physical review letters, 112 (2014), p. 253601.
  * [15] M. Ehrenberg, The nobel prize in chemistry 2014 (press release), 2014\.
  * [16] B. Fu, B. P. Isaacoff, and J. S. Biteen, Super-resolving the actual position of single fluorescent molecules coupled to a plasmonic nanoantenna, ACS nano, 11 (2017), pp. 8978–8987.
  * [17] M. Gramlich, J. Bae, R. Hayward, and J. Ross, Fluorescence imaging of nanoscale domains in polymer blends using stochastic optical reconstruction microscopy (storm), Optics express, 22 (2014), pp. 8438–8450.
  * [18] D. Grieser, The plasmonic eigenvalue problem, Reviews in Mathematical Physics, 26 (2014), p. 1450005.
  * [19] S. W. Hell, S. J. Sahl, M. Bates, X. Zhuang, R. Heintzmann, M. J. Booth, J. Bewersdorf, G. Shtengel, H. Hess, P. Tinnefeld, et al., The 2015 super-resolution microscopy roadmap, Journal of Physics D: Applied Physics, 48 (2015), p. 443001.
  * [20] B. Huang, M. Bates, and X. Zhuang, Super-resolution fluorescence microscopy, Annual review of biochemistry, 78 (2009), pp. 993–1016.
  * [21] K. L. Kelly, E. Coronado, L. L. Zhao, and G. C. Schatz, The optical properties of metal nanoparticles: the influence of size, shape, and dielectric environment, 2003.
  * [22] A. Moiola and E. A. Spence, Acoustic transmission problems: wavenumber-explicit bounds and resonance-free regions, Mathematical Models and Methods in Applied Sciences, 29 (2019), pp. 317–354.
  * [23] A. Moroz, Depolarization field of spheroidal particles, JOSA B, 26 (2009), pp. 517–527.
  * [24] T. Neuman, R. Esteban, D. Casanova, F. J. García-Vidal, and J. Aizpurua, Coupling of molecular emitters and plasmonic cavities beyond the point-dipole approximation, Nano letters, 18 (2018), pp. 2358–2364.
  * [25] H.-M. Nguyen, Limiting absorption principle and well-posedness for the helmholtz equation with sign changing coefficients, Journal de Mathématiques Pures et Appliquées, 106 (2016), pp. 342–374.
  * [26] M. A. Ordal, L. Long, R. Bell, S. Bell, R. Bell, R. Alexander, and C. Ward, Optical properties of the metals al, co, cu, au, fe, pb, ni, pd, pt, ag, ti, and w in the infrared and far infrared, Applied optics, 22 (1983), pp. 1099–1119.
  * [27] M. Raab, C. Vietz, F. D. Stefani, G. P. Acuna, and P. Tinnefeld, Shifting molecular localization by plasmonic coupling in a single-molecule mirage, Nature Communications, 8 (2017), p. 13966.
  * [28] D. Sarid and W. A. Challener, Modern introduction to surface plasmons: theory, Mathematica modeling, and applications, Cambridge University Press, 2010.
  * [29] B. E. Urban, B. Dong, T.-Q. Nguyen, V. Backman, C. Sun, and H. F. Zhang, Subsurface super-resolution imaging of unstained polymer nanostructures, Scientific reports, 6 (2016), pp. 1–9.
  * [30] K. A. Willets, A. J. Wilson, V. Sundaresan, and P. B. Joshi, Super-resolution imaging and plasmonics, Chemical Reviews, 117 (2017), pp. 7538–7582.