%~Mouliné par MaN_auto v.0.40.4 (9e9f6e09) 2026-08-27 10:08:59 \documentclass[AHL,Unicode,longabstracts]{cedram} \usepackage{dsfont} \usepackage{amssymb} \usepackage{mathtools} \AtBeginDocument{ \newtheorem{claim}[cdrthm]{Claim} \theoremstyle{plain} \newtheorem{claiml}{Claim} \renewcommand*\theclaiml{\Alph{claiml}} } \newcommand{\norm}[1]{\lVert#1\rVert} \newcommand{\ntwo}[1]{\left\lVert#1\right\rVert} \newcommand{\ovlineM}[1]{\bar{#1}^M} \newcommand{\rest}[1]{\bar{#1}^{\mathrm{rest}}} \newcommand{\1}{\mathbf{1}} \newcommand{\ind}{\mathds{1}} \newcommand{\E}{\mathbb{E}} \renewcommand{\P}{\mathbb{P}} \newcommand{\Ss}{\mathbb{S}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}} \newcommand{\T}{\mathbb{T}} \newcommand{\N}{\mathbb{N}} \newcommand{\F}{\mathcal{F}} \newcommand{\Mm}{\mathcal{M}} \newcommand{\cO}{\mathcal{O}} \newcommand{\cN}{\mathcal{N}} \newcommand{\cB}{\mathcal{B}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfun}{\mathbf{1}} \renewcommand{\d}{\mathrm{d}} \DeclareMathOperator{\Var}{Var} \DeclareMathOperator{\Unif}{Unif} \DeclareMathOperator{\Pois}{Pois} \DeclareMathOperator{\diag}{diag} \DeclareMathOperator{\Trace}{Trace} \DeclareMathOperator{\Id}{Id} \DeclareMathOperator{\Lip}{Lip} \DeclareMathOperator{\Law}{Law} \DeclareMathOperator{\avg}{\text{-}avg.} \newcommand{\odd}{\mathrm{odd}} \newcommand{\law}{\mathrm{law}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Better widebar \makeatletter \let\save@mathaccent\mathaccent \newcommand*\if@single[3]{% \setbox0\hbox{${\mathaccent"0362{#1}}^H$}% \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}% \ifdim\ht0=\ht2 #3\else #2\fi } %The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath: \newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna} %If there's a superscript following the bar, then no negative kern may follow the bar; %an additional {} makes sure that the superscript is high enough in this case: \newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}} %Use a separate algorithm for single symbols: \newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}} \newcommand*\wide@bar@[3]{% \begingroup \def\mathaccent##1##2{% %Enable nesting of accents: \let\mathaccent\save@mathaccent %If there's more than a single symbol, use the first character instead (see below): \if#32 \let\macc@nucleus\first@char \fi %Determine the italic correction: \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}% \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}% \dimen@\wd\tw@ \advance\dimen@-\wd\z@ %Now \dimen@ is the italic correction of the symbol. \divide\dimen@ 3 \@tempdima\wd\tw@ \advance\@tempdima-\scriptspace %Now \@tempdima is the width of the symbol. \divide\@tempdima 10 \advance\dimen@-\@tempdima %Now \dimen@ = (italic correction / 3) - (Breite / 10) \ifdim\dimen@>\z@ \dimen@0pt\fi %The bar will be shortened in the case \dimen@<0 ! \rel@kern{0.6}\kern-\dimen@ \if#31 \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}% \advance\dimen@0.4\dimexpr\macc@kerna %Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows: \let\final@kern#2% \ifdim\dimen@<\z@ \let\final@kern1\fi \if\final@kern1 \kern-\dimen@\fi \else \overline{\rel@kern{-0.6}\kern\dimen@#1}% \fi }% \macc@depth\@ne \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}% \macc@set@skewchar\relax \let\mathaccentV\macc@nested@a %The following initialises \macc@kerna and calls \mathaccent: \if#31 \macc@nested@a\relax111{#1}% \else %If the argument consists of more than one symbol, and if the first token is %a letter, use that letter for the computations: \def\gobble@till@marker##1\endmarker{}% \futurelet\first@char\gobble@till@marker#1\endmarker \ifcat\noexpand\first@char A\else \def\first@char{}% \fi \macc@nested@a\relax111{\first@char}% \fi \endgroup } \makeatother \let\oldbar\bar \renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \graphicspath{{./figures/}} \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} \hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red} \newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}} \newcommand*{\claimenumi}{\renewcommand*{\theenumi}{{\textbf{\arabic{enumi}}}}\relabel} \newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel} \newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel} \newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel} \newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel} \let\oldtilde\tilde \renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}} \let\oldhat\hat \renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}} \let\oldexists\exists \renewcommand*{\exists}{\mathrel{\oldexists}} \let\oldforall\forall \renewcommand*{\forall}{\mathrel{\oldforall}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \renewcommand*{\to}{\mathchoice{\longrightarrow}{\rightarrow}{\rightarrow}{\rightarrow}} \let\oldmapsto\mapsto \renewcommand*{\mapsto}{\mathchoice{\longmapsto}{\oldmapsto}{\oldmapsto}{\oldmapsto}} %~ This fixes spacing issues around \left( and \right) and other similar symbols \let\originalleft\left \let\originalright\right \renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft} \renewcommand{\right}{\aftergroup\egroup\originalright} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \title[The Spherical~model and the large-$N$ limit of the Spin~\texorpdfstring{$O(N)$}{O(N)}~model via GFF]{The~Spherical~model and the large-$N$ limit of the Spin~\texorpdfstring{$O(N)$}{O(N)}~model via the Gaussian free field} \alttitle{Le modèle sphérique et la limite pour $N$ grand du modèle de spins \texorpdfstring{$O(N)$}{O(N)} via le champ libre gaussien} \subjclass{60F05, 82B20} \keywords{spherical model, Spin O(N) model, Gaussian free field, zero-average Green's function} \author[\initial{J.} \lastname{Aru}]{\firstname{Juhan} \lastname{Aru}} \address{EPFL SB MATH RGM\\ MA B2 467 (Bâtiment MA)\\ Station 8\\ 1015 Lausanne (Switzerland)} \email{juhan.aru@epfl.ch} \author[\initial{A.} \lastname{Korzhenkova}]{\firstname{Aleksandra} \lastname{Korzhenkova}} \address{EPFL SB MATH RGM\\ MA A2 407 (Bâtiment MA)\\ Station 8\\ 1015 Lausanne (Switzerland)} \email{aleksandra.korzhenkova@gmail.com} \thanks{We would like to thank A.~Prévost for interesting discussions {and the referee for their careful reading and helpful comments and suggestions.} Both authors are supported by Eccellenza grant 194648 of the Swiss National Science Foundation and are members of NCCR Swissmap.} \begin{abstract} Using probabilistic limit theorems and concentration of measure, we prove that the infinite-volume limit of the spherical model on a $d$-dimensional torus is a massive GFF in the high-temperature regime, a standard GFF at the critical temperature, and a standard GFF plus a Rademacher random constant in the low-temperature regime. The proof at the critical temperature appears to be new and relies on a fine analysis of the zero-average Green's function on the torus. We study the spin $O(N)$ model in the double limit of spin dimensionality and torus size. Sending $N \to \infty$ first, and then the torus size to infinity, we show that the different spin coordinates become i.i.d. fields, distributed as a massive GFF in the high-temperature regime, a standard GFF at the critical temperature, and a standard GFF plus a Gaussian random constant in the low-temperature regime. \end{abstract} \begin{altabstract} En utilisant des théorèmes limites probabilistes et des résultats de concentration de la mesure, nous démontrons que la limite en volume infini du modèle sphérique sur un tore de dimension $d$ est un champ libre gaussien massif dans le régime de haute température, un champ libre gaussien standard à la température critique, et un champ libre gaussien standard auquel s'ajoute une constante aléatoire de Rademacher dans le régime de basse température. La preuve à la température critique semble être nouvelle et repose sur une analyse fine de la~fonction de Green de moyenne nulle sur le tore. Nous étudions le modèle de spins $O(N)$ dans la double limite de la dimension des spins et de la taille du tore. En faisant d’abord tendre $N$ vers l’infini, puis la taille du tore vers l’infini, nous montrons que les différentes coordonnées des spins deviennent des champs i.i.d., distribués comme un champ libre gaussien massif dans le régime de haute température, un~champ libre gaussien standard à la température critique, et un champ libre gaussien standard auquel s’ajoute une constante aléatoire gaussienne dans le régime de basse température. \end{altabstract} \datereceived{2024-07-21} \daterevised{2026-03-10} \dateaccepted{2026-05-19} \editors{S. Gou\"ezel and H. Lacoin} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} \maketitle \section{Introduction}\label{sec:intro} \subsection{Spherical model} Let $\T_n^d$ denote the $d$-dimensional discrete torus of side length $n$. To each vertex $x \in \T_n^d$, we associate a continuous ``spin'' $\theta_x \in \R$ in such a way that the full configuration $\theta = (\theta_x)_{x \in \T^d_n}$ lies on the $(n^d-1)$-dimensional sphere of radius $\sqrt{n^d}$, that is, $\norm{\theta}^2 \coloneqq \sum_x \theta_x^2 = n^d$. We say that $\theta$ is a configuration of the \emph{spherical model} on $\T^d_n$ at inverse temperature $\beta \geq 0$ if \begin{align*} \theta \sim \nu_{\T^d_n, \beta}[\d \theta] = \frac{1}{Z_{\T^d_n,\beta}} \exp \biggl(\frac{\beta}{2} \sum_{x \sim y} \theta_x \theta_y\biggr) \Unif_{\sqrt{n^d} \Ss^{n^d-1}}[\d \theta], \end{align*} where $x\sim y$ denotes the ordered pairs of neighboring vertices. The spherical model was first introduced by Berlin and Kac~\cite{SpherModel} in 1952 as a simplification of the Ising model that still exhibits a phase transition in dimensions $d \geq 3$, but allows for an~analytic treatment of the free energy in the infinite-volume limit. In their original work, Berlin and Kac used variations of the steepest-descent method to compute the limiting free energy, identify the phase transition in the three-dimensional lattice case, and argue that the limiting distribution of the model below the critical temperature is non-Gaussian. In a short note, Molchanov and Sudarev~\cite{MS75} stated that this method yields further results on the limiting distribution of the spherical model, but the proofs never appeared. In this article, we rigorously prove a subset of these results using probabilistic tools (cf. Theorem~\ref{thm:spherical}), which we believe also helps clarify the conceptual picture. An alternative approach is the so-called \emph{mean-spherical} one, in which one adds a~mass term $(s \sum_x \theta_x^2)$ to the Hamiltonian rather than conditioning the spin configuration. The parameter $s$ is then chosen such that the average value of $\norm{\theta}^2$ equals $n^d$. This method was proposed by Lewis and Wannier~\cite{LW52} shortly after the original work~\cite{SpherModel}. They observed that various quantities, such as correlations and the free energy, are significantly simpler to compute in this formulation. While the method works well in the high-temperature regime, it was soon realized~\cite{Lax55,LW53_corr} that it cannot be directly applied in the zero-magnetic-field setting below the critical temperature, as discrepancies with the original spherical model arise. However, in the presence of a nonzero external magnetic field $h \neq 0$, it was shown~\cite{YW65} that the thermodynamic limits of the spherical and mean-spherical models agree, and that the zero-magnetic-field case of the spherical model can be recovered by first taking the thermodynamic limit and then sending $h \to 0$. Since then, this method has appeared in several works~\cite{BD87,KT77, Shcherbina1988}, although none seem to provide a fully rigorous proof of the critical and low-temperature regimes for the zero-magnetic-field case. {In particular, \cite{BD87} proves the results of Theorem~\ref{thm:spherical} at the non-critical temperature under unverified technical assumptions.} Finally, the low-temperature regime of the spherical model can be also viewed as a toy model for condensation phenomena, {where the random constant term that appears in the limit (see Theorem~\ref{thm:spherical}) plays the role of a condensate}. In this broader context and motivated by the phenomena of Bose--Einstein condensation, \cite{Lukkarinen2020} investigates a complexified version of the spherical model below the critical temperature and {studies under which conditions the condensate could also consist of more than one component. Relevantly to the current article,} the author establishes convergence of local correlation functions, which, combined with tightness, implies local convergence in law in the aforementioned regime. While this proof also applies to the usual spherical model, it does not suffice to prove convergence at the critical temperature. {The key method in~\cite{Lukkarinen2020} differs from ours and is based on an interesting quantitative coupling, but, like our approach, it uses Fourier expansion}. Let us now state our main result about the limiting distribution of the spherical model and explain the reasoning leading to the proof. \begin{theo}[Infinite-volume limit of the spherical model]\label{thm:spherical} Let $d\geq 2$. {There exists $\beta_c = \beta_c(d) \in (0,\infty]$ such that} the spherical model on $\T^d_n$ at inverse temperature $\beta > 0$, $\theta = (\theta_x)_{x \in \T^d_n}$, converges in law, uniformly over compact subsets of $\Z^d$ (viewed as subsets of $\T^d_n$ for each $n$), as $n \to \infty$ to: \begin{enumerate} \item\label{theo1.1.1} $\beta < \beta_c$: a massive Gaussian free field (GFF) on $\Z^d$ scaled by $1/\sqrt{\beta}$ with the mass $m^2> 0$ depending on $\beta$ and $d$ in a specific way; \item\label{theo1.1.2} $\beta = \beta_c$: a GFF on $\Z^d$ scaled by $1/\sqrt{\beta}$; \item\label{theo1.1.3} $\beta > \beta_c$: a GFF on $\Z^d$ scaled by $1/\sqrt{\beta}$ plus an independent constant random drift $\sqrt{\frac{\beta - \beta_c}{\beta}} X$ with $X$ being a Rademacher random variable. \end{enumerate} Furthermore, all local correlations of spins converge and the critical temperature is given by $\beta_c = G_{\Z^d}(0,0)$. \end{theo} Here, as mentioned, the constant random drift appearing in the statement can be compared to the condensate in Bose--Einstein condensation phenomena~\cite{Lukkarinen2020}. To explain where the result stems from, we start by observing that for any $n$ and $m^2 > 0$, the constraint $\sum_x \theta_x^2 = n^d$ implies that \begin{align*} \nu_{\T^d_n, \beta}[\d \theta] \propto \exp \left(-\frac{\beta}{2} \left\langle \theta, \left(-\Delta_{\T^d_n} + m^2\right)\theta\right\rangle\right) \Unif_{\sqrt{n^d} \Ss^{n^d-1}}[\d \theta]. \end{align*} \pagebreak \noindent Hence, the law of $\sqrt{\beta}\theta$ is related to a massive GFF $\phi$ on $\T^d_n$ (for an arbitrary mass $m^2 > 0$) through conditioning on its norm to be equal to $\sqrt{\beta n^d}$, i.e., $\norm{\phi}^2 \coloneqq \sum_{x \in \T^d_n} \phi_x^2 = \beta n^d$ (see Proposition~\ref{prop:rel_GFF_spher}). Now the idea is somewhat similar to the one of Lewis and Wannier~\cite{LW52}, but from a probabilistic point of view: namely, we choose the mass $m_n^2$ in such a way that $\E[\norm{\phi}^2] = \beta n^d$ (or equivalently, $G_{\T^d_n, m_n^2}(0,0) = \beta$), which is natural due to the concentration properties of Gaussian measures. For large $n$ and fixed $m > 0$, $G_{\T^d_n, m^2}$ is roughly comparable to $G_{\Z^d, m^2}$; thus, the question of the existence of a critical point of the model (that is $\beta_c \in (0, \infty)$) boils down to whether there exists $m^2 \geq 0$ such that $G_{\Z^d, m^2}(0,0) = \beta$. Notice that in $d = 2$, since $G_{\Z^d}(0,0)$ is infinite, by choosing the mass small enough, one could get an arbitrarily large value of the variance $G_{\Z^d, m^2}(0,0)$, suggesting $\beta_c = \infty$; whereas in $d \geq 3$, there is a maximal possible value $\beta_c \coloneqq G_{\Z^d}(0,0) < \infty$. {In particular, the model exhibits a phase transition in dimensions three and higher, while it remains subcritical at all temperatures in dimension two.} More precisely, we observe that in the high-temperature regime $\beta < \beta_c$, the aforementioned sequence of masses $(m_n^2)_n$ converges to a positive number $m^2>0$ in such a way that \[ G_{\T^d_n, m_n^2}(0,0) \to G_{\Z^d, m^2}(0,0) = \beta. \] In this case, the concentration of measure for $\norm{\phi}^2$ is sufficient to make the conditioning disappear in the limit, and we obtain a purely Gaussian limit. Somewhat refined versions of classical limiting theorems help us handle this phase rather directly. The low-temperature phase $\beta > \beta_c$ is already trickier. One can check that \[ G_{\T^d_n, m_n^2}(x,y) \sim G_{\T^d_n}^{0\avg}(x,y) + \beta - \beta_c \overset{n \to \infty}{\sim} G_{\Z^d}(x,y) + \beta - \beta_c, \] where $G_{\T^d_n}^{0\avg}$ is the correlation function of the zero-average GFF on torus. This heuristically allows us to restate the problem in terms of the zero-average GFF $\gamma$ plus an independent zero mode (a constant-in-space Gaussian drift of the form $Z \bfun_{\T^d_n}$) conditioned on the norm of the sum being $\sqrt{\beta n^d}$. Since the drift is constant across all sites, almost surely $\frac{1}{n^d}\norm{\gamma + Z\bfun_{\T^d_n}}^2 = \frac{1}{n^d}\norm{\gamma}^2 + \lvert Z\rvert^2$. Due to high concentration of $\frac{1}{n^d}\norm{\gamma}^2$ around $\beta_c$, conditioning on the norm of the sum forces $\lvert Z\rvert$ to become constant in the limit; and by symmetry, we obtain a Rademacher random variable for the zero mode. Making this precise is slightly more subtle; nevertheless, the proof follows from a combination of basic concentration of measure results with relatively direct density bounds. Finally, in the critical case, we observe that $G_{\T^d_n}^{0\avg}(x,y) \sim G_{\Z^d}(x,y)$. However, further refinement of this relation is necessary and is granted by Theorem~\ref{thm:improved}: it shows that in $d \geq 3$, \[ G_{\T^d_n}^{0\avg}(x,y) = G_{\Z^d}(x,y) + O(n^{2-d}); \] and in $d=3$, where we further need to determine the sign of the error on the diagonal, that \[ G_{\T^d_n}^{0\avg}(x,x) < G_{\Z^d}(x,x) \] uniformly for large $n$. Obtaining the sign was surprisingly tricky and is related to the so-called Madelung constant, arising in the electrostatic potential of certain ionic crystals—a quantity of interest in chemical physics since its introduction in the early $20^{\rm th}$ century~\cite{Madelung}. See Section~\ref{subsec:further_questions} for a discussion of this connection. To deduce convergence of local spin correlations from convergence in law (Corollary~\ref{cor:spher_correlations}), we prove the existence of moments of single spherical spins (see Lemma~\ref{lemma:spher_moments}), inspired by an analogous result for complexified spherical spins in~\cite{Lukkarinen2020}. \subsection{Spin \texorpdfstring{$O(N)$}{O(N)} model} The spin $O(N)$ model, introduced (for general $N$) by Stanley in 1968~\cite{StanleyON}, is a fundamental model in statistical mechanics: for $N = 1$, one recovers the Ising model; for $N = 2$, the XY model; and for $N = 3$, the classical Heisenberg model. Here, we consider the infinite spin-dimensionality limit, $N \to \infty$, and its relation to the spherical model. A connection between the spin $O(N)$ model in the $N \to \infty$ limit and the spherical model was first suggested by Stanley~\cite{Stanley68} and later proved for any fixed non-critical temperature by Kac and Thompson~\cite{SpherO(N)Model}; see also~\cite{Shcherbina1988} and~\cite{gough1993}, which together cover the case of the critical temperature. More precisely, combining the steepest-descent method with an approach similar to the aforementioned mean-spherical one, they showed the equality of the free energies per site of the two models: if $S \in \prod_{x \in \T^d_n} (\sqrt{N} \Ss^{N-1})$ is a configuration of the \emph{spin $O(N)$ model} on $\T^d_n$ at inverse temperature $\beta \geq 0$, i.e., \begin{align*} S \sim \frac{1}{Z_{\T^d_n,N,\beta}} \exp\bigg(\frac{\beta}{2} \sum_{x \sim y} {S}_x \cdot {S}_y\bigg) \prod_{x \in \T^d_n} \Unif_{\sqrt{N} \Ss^{N-1}} [\d {S}_x], \end{align*} then (for $\beta \neq \beta_c$ in~\cite{SpherO(N)Model}), \begin{align*} \lim_{n, N \to \infty} \frac{1}{n^d N} \log Z_{\T^d_n,N,\beta} = \lim_{n \to \infty} \frac{1}{n^d} \log Z_{\T^d_n,\beta}. \end{align*} for any ordering of the limits $N, n \to \infty$. In this article, we take a closer look at this connection and observe that it actually passes by the Gaussian free field. In the following theorem, we describe the local limit of the spin $O(N)$ model when we first take the infinite spin-dimensionality limit $N \to \infty$ and then the thermodynamic limit $\T_n^d \nearrow \Z^d$. \begin{theo}[Infinite spin-dimensionality thermodynamic limit of the spin \texorpdfstring{$O(N)$}{O(N)} model]\label{thm:spinO(N)} Let $d\geq 2$. Any projection on the finitely many ($M \leq N$) coordinates of the spin $O(N)$ model on $\T^d_n$ at inverse temperature $\beta > 0$, $\ovlineM{S} \coloneqq (S^i_x)_{x \in \T^d_n}^{i = 1, \dots, M}$, converges in law as $N \to \infty$ to that of an $M$-vectorial massive GFF scaled by $1/\sqrt{\beta}$ with the mass $m_n^2$ depending on $\beta, d$ and $n$ in a specific way. {Furthermore, there exists $\beta_c = \beta_c(d) \in (0,\infty]$ such that} the consequent local (uniform on compacts of $\Z^d$) infinite-volume distributional limit ($n \to \infty$) has the law of: \begin{enumerate} \item\label{theo1.2.1} $\beta < \beta_c$: an $M$-vectorial massive GFF on $\Z^d$ scaled by $1/\sqrt{\beta}$ with the mass $m^2$ depending on $\beta$ and $d$ in a specific way; \item\label{theo1.2.2} $\beta = \beta_c$: an $M$-vectorial GFF on $\Z^d$ scaled by $1/\sqrt{\beta}$; \item\label{theo1.2.3} $\beta > \beta_c$: an $M$-vectorial GFF on $\Z^d$ scaled by $1/\sqrt{\beta}$ plus an independent constant random drift $\sqrt{\frac{\beta - \beta_c}{\beta}} \ovlineM{Z}$ with $\ovlineM{Z}$ being an $M$-dimensional standard normal vector. \end{enumerate} More precisely, $\beta_c = G_{\Z^d}(0,0)$. \end{theo} As a corollary (maybe a little surprisingly) we see that at least in the order of limits we take, although the free energies per site of the two models agree at all temperatures, their limiting laws are, in fact, different in the low-temperature regime. We expect this to hold regardless of the order in which the limits are taken. The a~posteriori explanation of this appearing discrepancy in the low-temperature regime is not particularly hard: differences in the law of a single degree of freedom (such as the drift part of the limiting field) vanish when computing the limiting free energy per site. \begin{coro}\label{intro:cor:sphr_spin} In the high-temperature regime and at the critical temperature $\beta_c$, the local distributional limit of the spherical model and the corresponding limiting law (obtained by first taking $N \to \infty$ and then the volume of the torus to infinity) of a single spin-coordinate process of the spin $O(N)$ model agree; while in the low-temperature regime, they differ in the law of the global zero mode. \end{coro} The proof of Theorem~\ref{thm:spinO(N)} consists of two steps: first we show that as the spin-dimensionality $N$ tends to infinity, the law converges locally (in spin-dimensionality) to a vector-valued massive GFF with mass $m_n^2$ chosen in such a way that $G_{\T^d_n, m_n^2}(0,0)\linebreak = \beta$; and then take the infinite-volume limit of the Gaussian measure of $m_n^2$-massive GFF on torus. The first step follows from a multidimensional version of the classical local limit theorem, similarly to our proof of the high-temperature regime of the spherical model, the proof of the latter step is quite simple and requires only a few well-known results about the Green's functions on $\T^d_n$ and $\Z^d$. We believe that the theorem should also hold in the opposite order of limits. In fact in the high-temperature regime this would follow directly from the work of A.~Kupiainen~\cite{o(n)expansion}; and we are hopeful that by combining some of his ideas with our approach, we would be able to extend this interchange of limits to all temperatures. The key issues in exchanging limits in our set-up are as follows. First, we cannot use similar density arguments in the whole space as on finite boxes. A~natural alternative would be to rather work with convergence of mixed moments. However, exchanging limits in that case requires uniformity in one of the directions. Understanding uniformity in $N$ for infinite-volume limits currently seems out of reach (at least at all temperatures), so one may instead hope to obtain uniformity in the box size for the $N \to \infty$ limits. This amounts to being able to say that, independently of the box size, local correlations converge uniformly. This is very reasonable to believe and is indeed shown in the high-temperature regime in~\cite{o(n)expansion} using reflection positivity. Trying to make use of reflection positivity in our set-up to obtain such uniformity seems to be an interesting direction. \subsection{Further questions and wider context} \label{subsec:further_questions} To complete the probabilistic picture of the spherical model, which at the very least provides a statistical physics model that can be fully analyzed in all dimensions and phases, one would need to investigate its behavior across a variety of domains, under different boundary conditions, and in the presence of a magnetic field, ideally extending the analysis to more general graphs and interaction structures. The qualitative effects of adding boundary conditions or a constant magnetic field are not difficult to anticipate, and our proofs can be adapted to cover these cases {with some additional technical work}. Similarly, one should be able to extend the proofs to long-range interactions, where the limiting Gaussian fields would given by fractional Gaussian fields; this problem has been considered in the physics literature~\cite{Joyce}, and the non-critical case has been treated in a recent master's thesis~\cite{master_thesis}. Extending our methods to more general (also weighted) graphs may, however, present new challenges {due to weaker control on the spectral/Fourier side}. As mentioned before, an interesting aspect of the critical case of the spherical model was the study of the zero-average Green's function on the torus and the question of whether it is smaller than the lattice Green's function. To this end, we analyzed the zero-average Green's function on the ``boundary'' (of the preimage of the torus in $\Z^d$ under the canonical projection), which is closely related to the so-called Madelung constant in chemical physics and mathematical chemistry~\cite{Mbook,Madelung}. Let us explain this connection: one way to study the sign of $G^{0\avg}_{\T^d_n}(0, y)$ for $y \in (\partial[-\lfloor n/2\rfloor, \lfloor n/2\rfloor]^d) \cap \Z^d$ is to take an appropriate scaling limit and examine the zero-average continuum Green's function $G_{\T^d}(0, x)$ for $x \in \partial[-1/2, 1/2]^d$, viewed as points on the torus $\T^d$ of side-length one. These quantities provide versions of the Madelung constant---for example, $d(2)$, $c(2)$, and $b(2)$ in~\cite[(1.3.29) \& Table~1.4]{Mbook} correspond precisely to $(2\pi)^2 G_{\T^3}(0, x)$ with $x = (1/2, 1/2, 1/2)$, $x = (0, 1/2, 1/2)$, and $x = (0, 0, 1/2)$, respectively. Determining or approximating the Madelung constants is a well-known problem. While explicit expressions exist (for instance, via Fourier series) the corresponding expansions are not absolutely convergent and may even diverge under natural truncations. In our case, we found a way to prove an inequality sufficient for our purposes directly on the lattice, although our bound is far from optimal. This naturally leads to several interesting questions of potential analysis. For example, it would be interesting to better understand the negativity region of the zero-average Green's function in all dimensions $d \geq 3$. Another perspective on the spherical model comes from its role as a simple framework for studying condensation phenomena~\cite{Lukkarinen2020}. The connection between the spherical model and the ideal Bose gas (via their shared critical behavior~\cite{PhysRev.166.152} and the emergence of Bose--Einstein condensation~\cite{gough1993}) was first noted in the late 1960s and early 1990s, and has attracted renewed interest in recent years~\cite{PhysRevResearch.1.023022, Lukkarinen2020}. It could be interesting to revisit this link, and to study other models that exhibit condensation, using the kind of probabilistic limit theorems and concentration-of-measure techniques we apply here, to see what extra insight such an approach might offer. We now turn to questions concerning the spin $O(N)$ model. To complete the picture of its limiting distribution, the first objective would be to find a way of interchanging the order of the infinite spin-dimensionality and infinite-volume limits at all temperatures. As mentioned above, we hope this to be within reach by combining our methods with those of~\cite{o(n)expansion}. In a broader context, an intriguing question is whether our approach could provide new insights into Polyakov's conjecture~\cite{POLYAKOV} on the exponential decay of correlation functions in 2D spin $O(N)$ models with $N \geq 3$. While it is difficult to see how our techniques might help to resolve the conjecture, one could hope to extract useful quantitative bounds on the critical temperatures in the 2D spin $O(N)$ model as $N \to \infty$ (cf.~\cite{o(n)expansion}). Finally, our results for the large-$N$ limit of the spin $O(N)$ model (particularly in light of Proposition~\ref{prop:dens_conv}) can also be placed in the context of propagation of chaos~\cite{SznitmanChaos}. The classical example is the Poincaré lemma, which states that the marginals of the spherical measure in high dimensions converge to independent Gaussians. In our setting, we observe that across all temperature regimes, the spin $O(N)$ measures are chaotic: as $N \to \infty$, the projection of the measure to any finite set of coordinates converges to a product measure. In principle, one could extend this analysis to spin models with additional non-quadratic interaction terms (such as quartic terms in the Hamiltonian) though it is not yet clear whether such generalizations would reveal anything beyond the behavior already observed. \subsection{Outline} Section~\ref{sec:setup} defines the precise setting of the spherical and spin $O(N)$ models, provides the necessary background on various versions of the GFF, and explains their connection to the two models. It also contains new results on the zero-average Green's function (Theorem~\ref{thm:improved}), which are proved in Section~\ref{sec:zero-avgG}. Section~\ref{sec:sphr} focuses on the spherical model: in particular, we prove Theorem~\ref{thm:spherical} and establish convergence of correlations for any finite number of its spins. Section~\ref{sec:O(N)} turns to the spin $O(N)$ model and proves Theorem~\ref{thm:spinO(N)}. To keep the article concise, some proofs less central to the main narrative have been moved to Appendix~\ref{sec:Appendix}. \section{Setup and preliminaries}\label{sec:setup} In this section, we introduce the notation used throughout the article, along with several variants of the Gaussian free field that are closely related to our models of interest, namely the spherical and spin $O(N)$ models. We also recall some properties of these GFFs and their corresponding Green’s functions (their correlation structures), which will be relevant for our proofs. The only new contribution stated in this section is {Theorem~\ref{thm:improved}}, which provides refined bounds on the zero-average Green's function on the discrete torus and compares it with the Green's function on $\Z^d$. Its proof is postponed to Section~\ref{sec:zero-avgG}. \subsection{Notation, definitions, and setup of our two models} Let $d \geq 2$, $n \in \N$, and $G$ be either the discrete $d$-dimensional torus of side length $n$, i.e., $G = \T^d_n \coloneqq (\Z/n\Z)^d$, or $\Z^d$, both viewed as a graph with the canonical nearest-neighbor structure. For $x,y \in G$, we write $\d(x,y)$ for the graph distance and $\lvert x-y \rvert$ for the Euclidean distance between the vertices. $x\sim y$ denotes two neighboring sites of $G$, and sums over $x\sim y$ run over pairs $(x, y)$ such that $\{x,y\}$ is an edge in $G$. Given $A\subset G$, we write $A^c = G\setminus A$ for the complement of $A$ in $G$, $\bfun_A$ for the vector $(1)_{x \in A}$, and, analogously, for any given $w\in \R^{G}$ we set $w_A \coloneqq (w_x)_{x \in A}$. If $G = \T^d_n$, to highlight the dependence on $n$, we often write $\Lambda_n$ instead. When considering $A \subset \Z^d$, we also view it as a subset of $\Lambda_n$ for each $n$ (by canonically projecting it onto $\Lambda_n$). Conversely, for $U \subset \Lambda_n$, we denote by $\hat{U}$ the unique preimage of $U$ under the canonical projection contained in $[-n/2, n/2)^d$—in particular, $\hat{\Lambda}_n = [-n/2, n/2)^d \cap \Z^d$. Moreover, let $\langle\,\cdot\,,\,\cdot\,\rangle$ be the canonical inner product on $\R^{A}$ and $\norm{\,\cdot\,}$ be the associated norm, in order to distinguish them from the Euclidean inner product and norm on $\R^k$ for $k \in \N$, which we denote by $\cdot$ and $|\cdot|$, respectively. For the reader's convenience, we recall the definitions of the spherical and spin $O(N)$ models in the setting of interest. \begin{defi}[Spherical model]\label{def:spher} Let $\beta > 0, n \in \N$. We call $\theta \coloneqq (\theta_x)_{x\in \Lambda_n} \in \sqrt{n^d} \;\Ss^{n^d-1}$ a configuration of the {spherical model} on $\Lambda_n$ at inverse temperature $\beta$ if $\theta \sim \nu_{\Lambda_n,\beta}$ with \begin{align*} \nu_{\Lambda_n,\beta} [\d \theta] = \frac{1}{Z_{\Lambda_n, \beta}} \exp\biggl(\frac{\beta}{2} \sum_{x\sim y} \theta_x\theta_y \biggr) \Unif_{\sqrt{n^d}\;\Ss^{n^d-1}}[\d \theta], \end{align*} where $Z_{\Lambda_n, \beta}$ is the normalizing constant. \end{defi} \begin{defi}[Spin $O(N)$ model]\label{def:spinO(N)} Let $\beta > 0, n, N \in \N$. We call $S \coloneqq (S_x)_{x\in \Lambda_n} \in (\sqrt{N} \: \Ss^{N-1})^{n^d}$ a configuration of the {spin $O(N)$ model} on $\Lambda_n$ at inverse temperature $\beta$ if $S \sim \mu_{\Lambda_n, N, \beta}$ with \begin{align*} \mu_{\Lambda_n, N,\beta} [\d S] = \frac{1}{Z_{\Lambda_n, N, \beta}} \exp\biggl(\frac{\beta}{2} \sum_{x\sim y} S_x \cdot S_y \biggr) \prod_{x \in \Lambda} \Unif_{\sqrt{N}\;\Ss^{N-1}}[\d S_x], \end{align*} where $Z_{\Lambda_n, N, \beta}$ is the normalizing constant. \end{defi} In the sequel we will always refer to a configuration of the spherical model as $\theta$ and of the spin $O(N)$ model as $S$. \subsection{Massive GFF, its vectorial descendant and zero-average GFF}\label{subsec:GFF} In the present subsection, we review several variants of the GFF, which, as we will see later, are closely related to our models of interest both in the finite setting and in the distributional limit. We begin by defining the vector-valued (massive) Gaussian free field on a discrete torus with Dirichlet boundary conditions. \begin{defi}[$m^2$-massive $N$-vectorial GFF]\label{def:GFF} Let $m \geq 0, d\geq 2$, $N, n \in \N$, and $G = \T^d_n$. Let $U \subset G$ be an arbitrary subset (possibly empty if $m>0$). We call a~random function $\phi: G \to \R^N$ an $m^2$-massive $N$-vectorial Gaussian free field (GFF) on $G$ with Dirichlet boundary condition on $U$ if $\phi_x = 0$ for all $x \in U$ and \begin{align}\label{eq:gff_def} \P[\d \phi] \propto \exp \biggl(- \frac{1}{4}\sum_{x\sim y} |\phi_x - \phi_y|^2 - \frac{m^2}{2} \sum_{x \in G} |\phi_x|^2\biggr) \prod_{x \in G\setminus U} \d \phi_x. \end{align} If $m^2 = 0$, we drop the term ``massive'' from the name; if $N=1$, the ``$N$-vectorial'' part. If $U = \emptyset$, we drop the ``Dirichlet'' part altogether. \end{defi} \begin{rema} From the expression in~\eqref{eq:gff_def}, it is evident that the coordinate processes of $\phi$ are i.i.d., each being an $m^2$-massive GFF on $G$ with Dirichlet boundary condition on $U$. We therefore restrict to the scalar case when studying the covariance structure of the field and when discussing component-wise properties of $\phi$. \end{rema} When $N = 1$, we can rewrite the expression in the exponent of~\eqref{eq:gff_def} as follows \begin{align}\label{eq:discrGaussform} \frac{1}{2}\sum_{x\sim y} |\phi_x - \phi_y|^2 + m^2\sum_{x} |\phi_x|^2 = \bigl\langle\phi, (-\Delta_{G} + m^2) \phi\bigr\rangle, \end{align} where $\langle\,\cdot\,,\,\cdot\,\rangle$ stands for the inner product on $\R^{G}$ and \begin{align*} -\Delta_{G} f (x) = 2d f(x) - \sum_{y \in G: y\sim x} f(y) \; \text{ for any } \,x \in G, \quad f: G \to \R \end{align*} is the discrete Laplacian on $G$. When restricted to the set $F_U$ of functions $f: G \to \R$ with $f(x) = 0$ for all $x \in U$, we denote it by $-\Delta_{G\setminus U}$ (or $-\Delta_{U^c}$) and call it the discrete Laplacian with Dirichlet boundary condition on $U$. As a straightforward consequence of~\eqref{eq:discrGaussform}, we see that $-\Delta_{U^c}$ is invertible on $F_U$ (note that we assume that $U \neq \emptyset$ or $m^2 >0$) and its inverse $(G_{G\setminus U, m^2}(x,y))_{x,y \in G}$, which is then also the covariance matrix of $\phi$, is the massive (if $m^2>0$) Green's function on $G$ with Dirichlet boundary condition on $U$. If $m^2 = 0$, we often drop the corresponding subscript and write simply $G_{G\setminus U}$. It is a well-known fact (see for instance~\cite[Section~1]{Rod17}) that $G_{U^c, m^2}$ (where $U^c = G\setminus U$) is related to random walks in the following way: for all $x,y \in G$, $\alpha \coloneqq \frac{m^2}{2d+m^2}$ \begin{equation}\label{eq:Gfct_RW} \begin{aligned}[t] G_{U^c, m^2}(x,y) \,\MK &\!\Mk:= \frac{1-\alpha}{2d} \sum_{k \geq 0}\P^x_{\alpha, U}[X_k = y]\\ &= \frac{1-\alpha}{2d} \sum_{k \geq 0} (1-\alpha)^k \P^x_{0}[X_k = y, k< H_U], \end{aligned} \end{equation} where $\P_{\alpha,U}^x$ denotes the law of a random walk $X$ on the graph $G\cup \{\dag\}$ (obtained from $G$ by adding an edge connecting each vertex of $G$ to $\dag$) starting at $x \in G$ with transition probabilities \begin{align*} p_{u,v} = \frac{1-\alpha}{2d} \ind_{\{u\sim v\}}, \quad p_{u,\dag} = \alpha, \quad p_{w,w} = 1, \quad \text{for all } u \in U^c,\; v\in G, \;w\in U \cup \{\dag\}, \end{align*} $H_U\coloneqq \inf\{k \geq 0: X_k \in U\}$ is the first hitting time of $U$, and $\P^x_{0}$ is the law of a simple random walk (SRW) on $G$ started at $x \in G$. Observe that the right-hand side of~\eqref{eq:Gfct_RW} remains well-defined (finite) if $G = \Z^d$ for $d\geq 2$ (we denote it by $G_{\Z^d\setminus U, m^2}$) under the assumption that $U\neq \emptyset$ or $m^2 > 0$ if $d=2$, since a SRW on $\Z^d$ is transient for $d\geq 3$, while it is recurrent for $d=2$. This also explains our assumption on $m^2, U$ if $G = \T^d_n$, since a SRW on $\T^d_n$ is recurrent. Under this assumption on $U$ and $m^2$, we can now define an \emph{$m^2$-massive $N$-vectorial GFF on $\Z^d$ with Dirichlet boundary condition on $U \subset \Z^d$} as a vector of $N$ i.i.d. centered Gaussian processes with covariance matrix $G_{\Z^d\setminus U, m^2}$ each. Note furthermore that from~\eqref{eq:Gfct_RW}, if $G = \Lambda_n$, we get that for any fixed finite $U \subset \Z^d$ and $m^2 \geq 0$ ($m^2 > 0$ if $U=\emptyset$), as $\Lambda_n\nearrow \Z^d$, an $m^2$-massive $N$-vectorial GFF $\phi$ on $\Lambda_n$ with Dirichlet boundary condition on $U$ converges in law, uniformly on compact subsets of $\Z^d$, to an $m^2$-massive $N$-vectorial GFF $\psi$ on $\Z^d$ with Dirichlet boundary condition on $U$. Moreover, we can extend the argument to a sequence of $m_n^2$-massive GFFs, each defined on $\Lambda_n$, under the assumption that $m_n^2\to m^2 > 0$ as $n\to \infty$. Indeed, \eqref{eq:Gfct_RW} guarantees the existence of $C>0$, uniform in $n$, such that $\lvert G_{\Lambda_n\setminus U, m_n^2} - G_{\Lambda_n \setminus U, m^2}\rvert \leq C \lvert \alpha - \alpha_n\rvert G_{\Lambda_n \setminus U, m^2}$. The latter, in turn, converges to zero uniformly over any compact subset of $\Z^d$ as $n$ tends to infinity. Let us end the discussion of a massive GFF in this subsection by stating one more of its canonical properties that will be frequently used in the sequel. Its proof may be found in~\cite[Lemma~1.1]{Rod17} (see also~\cite[Lemma~3.1]{Biskup_notes} for the massless case). \begin{prop}[Domain Markov property of the (massive) GFF]\label{prop:domainMP} Let $G$ be either $\T^d_n$ or $\Z^d$, $U \subset G$, and $m^2 \geq 0$ ($m^2 > 0$ if $U=\emptyset$ and $G = \T^d_n$ or $G = \Z^2$). Let $\phi$ be an $m^2$-massive GFF on $G$ with Dirichlet boundary condition on $U$. For $K \subset G$, set \begin{align*} h^{\phi_K}_x = \E[\phi_x |\F_K] \quad \text{for all } x \in G, \end{align*} where $\F_K \coloneqq \sigma(\phi_x: x \in K)$. Then $h^{\phi_K} = (h^{\phi_K}_x)_{x\in G}$ is discrete massive-harmonic in $G\setminus (K \cup U)$ (i.e., $(-\Delta_{G} + m^2) h^{\phi_K} (y) = 0$ for any $y \in G\setminus (U\cup K)$) with boundary values determined by $h^{\phi_K}_x = \phi_x$ for any $x \in K\cup U$. In terms of the random walk $X \sim \P_{\alpha,U}^x$ as in~\eqref{eq:Gfct_RW}, it is given by \begin{align*} h^{\phi_K}_x = \sum_{y \in K} \P^x_{\alpha, U}[H_K < \infty, X_{H_K} = y] \:\phi_y. \end{align*} Moreover, the field $\phi^{G\setminus K} \coloneqq (\phi_x - h^{\phi_K}_x)_{x\in G}$ is independent of $\F_K$ and has the law of an $m^2$-massive GFF on $G$ with Dirichlet boundary condition on $K \cup U$. \end{prop} We conclude this subsection by introducing another version of a (massless) GFF on the discrete torus $\Lambda_n = \T^d_n$ in dimension $d \geq 3$: the so-called zero-average GFF $\gamma$. The name reflects the constraint imposed to reduce the dimensionality of the space on which the GFF measure is defined—namely, that $\sum_x \gamma_x = 0$. \begin{defi}[Zero-average GFF on torus]\label{def:zeroavGFF} Let $n \in \N, d \geq 3$. We call a random function $\gamma: \Lambda_n \to \R$ a {zero-average Gaussian free field on $\Lambda_n$} if it is distributed according to \begin{align}\label{eq:0-avf_GFF_def} \P[\d \gamma] \propto \left(- \frac{1}{4}\sum_{x\sim y} (\gamma_x - \gamma_y)^2 \right) \delta_0 \left(\sum_x \gamma_x \right) \prod_{x\in \Lambda_n} \d \gamma_x. \end{align} \end{defi} To find the covariance matrix of $\gamma$ while avoiding direct computation with the Dirac measure, let $(\eta_w)_{w \in [0,n)^d \cap \Z^d}$ be the eigenvalues of $-\Delta_{\Lambda_n}$ and $(q^w)_w$ be the corresponding orthogonal eigenvectors (cf. Lemma~\ref{lemma:massive_EVs}); in particular, $q^0 = \frac{1}{\sqrt{n^d}} \1_{\Lambda_n}$. Set $Q \coloneqq (q^0, (q^w)_{w \neq 0})$ and observe that \begin{align*} \frac{1}{2} \sum_{x\sim y} (\gamma_x - \gamma_y)^2 = \left\langle Q^T \gamma, \,\diag\left(0, (\eta_w)_{w \neq 0}\right) Q^T\gamma \right\rangle; && \frac{1}{\sqrt{n^d}} \sum_{x \in \T^d_n} \gamma_x = \langle q^0, \gamma\rangle = (Q^T \gamma)_0. \end{align*} This allows us to conclude that the zero-average GFF has the same law as the orthogonal transformation of a vector of $n^d-1$ independent centered normal random variables $Z_w$ with variance $1/\eta_w$ for $w \in [0,n)^d \cap\Z^d \setminus \{0\}$: \begin{align}\label{eq:gamma_law_iid} \gamma \overset{\law}{=} Q \bigl(0, (Z_w)_{w\neq 0}\bigr)^T, \end{align} which in turn implies that the covariance matrix of $\gamma$, known as the \emph{zero-average Green's function}, is given by \begin{align*} \left(G_{\Lambda_n}^{0\avg}(x,y)\right)_{x,y \in \Lambda_n} \coloneqq Q \,\diag\left(0, \left(\frac{1}{\eta_k} \right)_{w \in [0,n)^d \cap\Z^d \setminus \{0\}} \right) \,Q^T. \end{align*} We recall and prove several of its properties in Section~\ref{subsec:Greens_fcts}. \subsection{Massive GFF, its vectorial descendant, and the zero-average GFF: relation to the spherical and spin \texorpdfstring{$O(N)$}{O(N)} models}\label{subsec:rel_GFF_spher_spinO(N)} Let $d \geq 2$, $n$ and $N \in \N$. Recall that $\Lambda_n = \T^d_n$ is the $d$-dimensional discrete torus. Let $\phi$ be an $m^2$-massive GFF on $\Lambda_n$ for some $m^2 >0$. Then, using polar coordinates, it is easy to obtain an expression for the conditional density of $\frac{\phi}{\text{\small{$\sqrt{\beta}$}}}$ given $\norm{\phi} = \sqrt{\text{\small{$\beta n^d$}}}$ at $v \in \R^{\Lambda_n}$ with $\norm{v} = \sqrt{n^d}$: \begin{align*} f_{\frac{\phi}{\sqrt{\beta}}\big\vert \ntwo{\phi} = \sqrt{\beta n^d}} (v) \d v &\propto \exp\left(-\frac{\beta}{2} \left\langle v, \left(-\Delta_{\Lambda_n} + m^2\right) v \right\rangle\right) \Unif_{\sqrt{n^d} \: \Ss^{n^d-1}} (\d v) \\ &\propto \exp\left(\frac{\beta}{2} \sum_{x\sim y} v_x v_y \right) \Unif_{\sqrt{n^d} \: \Ss^{n^d-1}} (\d v) \propto \nu_{\Lambda_n, \beta}. \end{align*} In the last line, we used that $\langle v, v \rangle = n^d$. If $d\geq 3$, consider a zero-average GFF $\gamma$ on $\Lambda_n$, and let and $c$ is an ``independent'' Lebesgue constant, that is, $(\gamma, c)$ is distributed as $\cN (0, G_{\Lambda_n}^{0\avg}) \times \d c$. We are interested in the conditional law of the sum of two ($\gamma + c\bfun_{\Lambda_n}$) given that it equals $\sqrt{\beta n^d}$. Note that upon conditioning, we indeed have a probability measure. To better understand this, let us instead consider the orthogonal transformation of the sum (this does not affect the conditioning). More precisely, let $Q$ be the orthogonal matrix diagonalizing $G_{\Lambda_n}^{0\avg}$, that is, $Q^T G_{\Lambda_n}^{0\avg} Q$ is diagonal. It is clear that the ``law'' of $\gamma_{c, Q} \coloneqq Q^T(\gamma + c \bfun_{\Lambda_n})$ is given by $\frac{1}{\sqrt{n^d}} \d c \times \bigtimes_{k=2}^{n^d} \cN (0, \eta_k^{-1})$, where $(\eta_k)_{k \geq 2}$ are the non-zero eigenvalues of $-\Delta_{\Lambda_n}$. Restricted to a finite ball and normalized accordingly, it is a probability measure that we denote by $\rho$. For our purposes we can consider, e.g., $\cB^{n^d}(0, 2\sqrt{\beta n^d})$. Then, analogously to the computation in the massive case, for any Borel-measurable $B \subset \R^{\Lambda_n}$, \begin{multline}\label{eq:cond_prob_meas} \rho\Bigl[B \Bigm\vert~\ntwo{\gamma_{c, Q}} = \sqrt{\beta n^d}\Bigr] \propto \int \ind_{B} (v) \prod_{k=2}^{n^d} \exp\left(-\frac{\beta}{2} \eta_k v_k^2 \right) \Unif_{\sqrt{n^d} \: \Ss^{n^d-1}} (\d v) \\ = \int \ind_{Q B} (\sigma) \exp\left(\frac{\beta}{2} \sum_{x\sim y} \sigma_x \sigma_y \right) \Unif_{\sqrt{n^d} \: \Ss^{n^d-1}} (\d{\sigma}) \propto \nu_{\Lambda_n, \beta} [Q B] \end{multline} In the second line we used the substitution $v = Q^T \sigma$ so that $\sum_{k=2}^{n^d} \eta_k v_k^2 = \langle \sigma, -\Delta_{\Lambda_n} \sigma\rangle$, and the fact that the uniform measure on the sphere is invariant under orthogonal transformations. Let us summarize our findings in the following proposition. \begin{prop}[Massive and zero-average GFFs: relation to the spherical model]\label{prop:rel_GFF_spher} Let $\nu_{\Lambda_n, \beta}$ be the spherical model on $\Lambda_n$ at inverse temperature $\beta > 0$. Then, for any $m^2>0$, the conditional law of an $m^2$-massive GFF $\phi$ on $\Lambda_n$ given $\norm{\phi} = \sqrt{\beta n^d}$ is exactly $\nu_{\Lambda_n, \beta}$. That is, \begin{align*} \nu_{\Lambda_n, \beta} = \Law \left(\frac{\phi}{\sqrt{\beta}}\,\middle|\,\ntwo{\phi} = \sqrt{\beta n^d} \right), \end{align*} Moreover, if $d\geq 3$, $\gamma$ is a zero-average GFF on $\Lambda_n$, and $c$ is an ``independent'' Lebesgue constant, i.e., the distribution of $(\gamma, c)$ is $\cN (0, G_{\Lambda_n}^{0\avg}) \times \d c$, then \begin{align*} \nu_{\Lambda_n, \beta} = \Law \left(\frac{\gamma + c \bfun_{\Lambda_n}}{\sqrt{\beta}}\,\middle|\ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d} \right) \end{align*} in the sense of~\eqref{eq:cond_prob_meas}. \end{prop} We also have a similar statement for the spin $O(N)$ model, which was implicitly used already in~\cite{o(n)expansion} and proven in a slightly different setting in~\cite[Proposition~2.3]{AGS22}. We omit its proof as it is fully analogous to the proof of the first statement of the previous proposition and can also be easily adjusted from the latter reference. \begin{prop}[Relation between an $N$-vectorial massive GFF and the spin \texorpdfstring{$O(N)$}{O(N)} model]\label{prop:mGFF-O(N)} Let $m^2> 0$ and $\phi$ be an $N$-vectorial $m^2$-massive GFF on $\Lambda_n$, $\beta > 0$. The conditional law of $\phi/\sqrt{\beta}$ given $|\phi_x| = \sqrt{\beta N}$ for all $x\in \Lambda_n$ is that of the spin $O(N)$ model on $\Lambda_n$ at inverse temperature $\beta$. \end{prop} \subsection{Preliminary estimates on the Green's functions}\label{subsec:Greens_fcts} This subsection mainly recalls but also proves some further results about Green's functions and their eigenvalues. We start off with the massive Green's function and then turn to the zero-average Green's function. \subsubsection{Massive Green's function} We begin with some canonical results, whose proofs can be found, for example, in~\cite[Section~1]{Rod17} and~\cite[Proposition~8.30]{friedli_velenik_2017}. \begin{lemm}[{Properties of massive Green's function}]\label{lemma:prop:lapl,grfct} Let $n \in \N$, $G = \T^d_n$ or $\Z^d$, and $U$ be a subset of $G$, possibly $U = \emptyset$. Let $m^2\geq 0$ if $U \neq \emptyset$ or $G = \Z^d$ for $d\geq 3$ and $m^2 > 0$ otherwise. \begin{enumerate} \item\label{lemm2.9.1} For any finite subset $K \subset G\setminus U \eqqcolon U^c$, \begin{align}\label{eq:Gfct_decomp} G_{U^c, m^2}(x,y) = G_{(U\cup K)^c, m^2} + \E_{\alpha, U}^x\left[\ind_{\{H_K < \infty\}} G_{U^c, m^2}(X_{H_K},y)\right] \end{align} with the notation of~\eqref{eq:Gfct_RW}. \item\label{lemm2.9.2} If $m^2>0$, there exist two constants $c, C> 0$ depending only on $m^2$ and $d$ such that \begin{align}\label{eq:Gfct_expdec} G_{U^c, m^2}(x,y) \leq C e^{-c \lvert x-y\rvert}, \quad\text{for all } U\subset G, x, y \in G. \end{align} \end{enumerate} \end{lemm} We will frequently find it helpful to calculate with the basis of the eigenfunctions of the massive Green's function on $\T^d_n$. \begin{lemm}[{Eigenvalues of massive Green's function}]\label{lemma:massive_EVs} Let $n \in \N, m^2 > 0$, and $U$ be a finite subset of $\Z^d$ viewed as a subset of $\Lambda_n$ for each $n$. Set $U^c \coloneqq \Lambda_n \setminus U$ and $u_n \coloneqq \lvert U^c\rvert$. \begin{enumerate} \item\label{lemm2.10.1} $-\Delta_{\Lambda_n}$ has one zero eigenvalue $\eta_0 = 0$ and $(n^d-1)$ positive eigenvalues given by \begin{align*} \eta_{w} = 2 \sum_{i=1}^d \left(1-\cos\left(2\pi \frac{w_i}{n} \right)\right) \quad\text{for } w=(w_i)_{i=1}^d \in [0, n)^d \cap \Z^d \setminus \{0\}. \end{align*} The corresponding eigenvectors $(q^w)_{w\in [0,n)^d \cap \Z^d}$ generating an orthonormal system are given by \begin{align*} q^w_x = \frac{1}{\sqrt{n^d}} \prod_{i=1}^d \left(\cos\left(2\pi \frac{x_i w_i}{n}\right) + \sin\left(2\pi \frac{x_i w_i}{n}\right)\right)\quad \text{for all } x \in [0, n)^d \cap \Z^d. \end{align*} \item\label{lemm2.10.2} Let $(\eta_k)_{k \geq 2}$\footnote{Here and in the sequel, we slightly abuse notation and write both $(\eta_k)_{k=1}^{n^d}$ (when ordered non-decreasingly) and $(\eta_x)_{x \in [0, n)^d \cap \Z^d}$ (when indexed by the domain) for the eigenvalues of $-\Delta_{\Lambda_n}$\label{footnote:abuse_EVs}. } be the positive eigenvalues of $-\Delta_{\Lambda_n}$ ordered non-decreasingly. Then, \begin{equation}\label{eq:EV_asympt} \eta_{k} = \Theta\left(\frac{k^{2/d}}{n^2}\right). \end{equation} \item\label{lemm2.10.3} The eigenvalues of $G_{\Lambda_n, m^2}$ are given by $(1/(m^2 + \eta_x))_{x \in [0,n)^d \cap \Z^d}$, while the eigenvalues $(\tau_i)_{i=1}^{u_n}$ of $G_{U^c, m^2}\vert_{U^c \times U^c}$ ordered non-increasingly satisfy \begin{align}\label{eq:unif_bdsEVs} \tau_i \in \left[\frac{1}{m^2 + 4d}, \frac{1}{m^2}\right] \quad \text{for all } 1\leq i \leq u_n. \end{align} \end{enumerate} \end{lemm} \begin{proof} The eigenvalues of $-\Delta_{\Lambda_n}$ can be found, for instance, in~\cite[Example~5.17]{aldous-fill-2014} or easily derived from~\cite[Appendix~A]{SpherModel}; the eigenvectors are also classical. Using explicit formulae for $(\eta_x)_{x \in [0, n/2]^d \cap \Z^d}$ (any remaining $x$ gives an eigenvalue coinciding with one indexed by some point in $[0, n/2]^d \cap \Z^d$), we obtain $\eta_x = \Theta(\lvert x\rvert^2/n^2)$. This together with symmetries of the torus and properties of $\cos$-function implies that $\eta_{k} = \Theta (k^{2/d}/n^2 )$. The third point follows immediately since $G_{\Lambda_n, m^2}$ is the inverse of $-\Delta_{\Lambda_n} + m^2$, while $G_{U^c, m^2}\vert_{U^c \times U^c}$ is the inverse of $(-\Delta_{U^c} + m^2)\vert_{U^c \times U^c}$---see discussion following~\eqref{eq:discrGaussform}. From~\eqref{eq:discrGaussform}, one furthermore sees that the smallest eigenvalue of $(-\Delta_{U^c} + m^2)\vert_{U^c \times U^c}$ is at least $m^2$, while the largest is upper bounded by $4d + m^2$. \end{proof} \begin{rema}\label{rem:EV_indexing} Note that equivalently the eigenvalues and the corresponding eigenfunctions as well as their domains can be indexed by $[-n/2, n/2)^d\cap \Z^d$ instead of $[0,n)^d\cap \Z^d$ due to the invariance of cosine functions under translation of the argument by $2\pi \Z$. More precisely, for $w,x \in [0,n)^d\cap \Z^d$, $\eta_w = \eta_{\tilde w}$ and $q^w_x = q^{\tilde w}_{\tilde x}$, where $\tilde y_i = y_i$ if $y_i < n/2$ and $y_i - n$ otherwise for $y = x, w$. \end{rema} \subsubsection{The zero-average Green's function} We now proceed to the results on the zero-average Green's functions corresponding to the covariance structure of our zero-average GFF on $\Lambda_n$. We start by recalling some known results stemming from~\cite{Aba17}: \begin{prop}[Properties of the zero-average Green's function]\label{prop:props:zeroAvGfct} Let $d \geq 3$, $n \in \N$. The zero-average Green's function $G_{\Lambda_n}^{0\avg}$ satisfies the following: \begin{enumerate} \item\label{prop2.12.1} For any $y \in \Lambda_n$, $\sum_{x \in \Lambda_n} G_{\Lambda_n}^{0\avg}(x,y) = 0$. \item\label{prop2.12.2} Let $\bar{X} = (\bar{X}_t)_{t \geq 0}$ be a continuous-time simple random walk on $\Lambda_n$. We denote its law when started at $x$ by $\P_{\Lambda_n}^x$. Then, for all $x,y \in \Lambda_n$, \begin{equation*} G_{\Lambda_n}^{0\avg}(x,y) = \frac{1}{2d} \int_0^\infty \left(\P^x_{\Lambda_n}[\bar{X}_t = y] - \frac{1}{n^d} \right) \d t. \end{equation*} \item\label{prop2.12.3} For any $U \subsetneq \Lambda_n$, $x,y \in \Lambda_n$, \begin{align}\label{eq:0-avGfct_decomp} G_{\Lambda_n}^{0\avg}(x,y) = G_{U^c, 0}(x,y) + \E_{0}^x \left[G_{\Lambda_n}^{0\avg}(X_{H_U},y)\right] - \frac{1}{2d n^d} \E_{0}^x[H_U], \end{align} where $H_U= \inf\{n \geq 0: X_n \in U\}$ is the first hitting time of $U$ by a simple random walk $(X_k)_{k \in \N_0} \sim \P^x_{0}$ on $\Lambda_n$ started at $X_0 = x$. \item\label{prop2.12.4} $G_{\Lambda_n}^{0\avg}$ converges to $G_{\Z^d}$ uniformly over compact subsets of $\Z^d$ (viewed also as subsets of $\Lambda_n$). \item\label{prop2.12.5} There exists $C>0$ such that for all $x, y \in \Lambda_n$ and $n \in \N$ sufficiently large, \begin{align}\label{eq:zero-avGfct_polydec} \bigl\lvert G_{\Lambda_n}^{0\avg}(x,y) \bigr\rvert \leq C (\log n)^{3d/2} \d(x,y)^{2-d}. \end{align} Recall that $\d(\,\cdot\,,\,\cdot\,)$ stands for the graph-distance on $\Lambda_n$. \end{enumerate} \end{prop} For most of the results of this paper, the bound~\eqref{eq:zero-avGfct_polydec} would have been sufficient, however, not for the critical regime of the spherical model in dimensions~$3$ and $4$. Therefore, we establish the improved estimates, whose proof is given in Section~\ref{sec:zero-avgG}. \begin{theo}[Fine estimates for the zero-average Green's function]\label{thm:improved} Let $d \geq 3, n \in \N$. The zero-average Green's function $G_{\Lambda_n}^{0\avg}$ satisfies the following: \begin{enumerate} \item\label{theo2.13.1} There exists $C = C(d)>0$ such that for all $n$ sufficiently large and $y \in [-n/2, n/2)^d \cap \Z^d$, \begin{align}\label{eq:0avgGfct_polydec} \left\lvert G_{\Lambda_n}^{0\avg}(0,y) - G_{\Z^d}(0,y) \right\rvert \leq C n^{2-d}. \end{align} In particular, for $\lvert y\rvert = \Theta(n)$ as above, $G_{\Lambda_n}^{0\avg}(0,y) = \cO(n^{2-d})$. \item\label{theo2.13.2} For $d=3$ and any $x \in \Z^d$, \begin{align}\label{eq:negativity_Greens_diff} G^{0\avg}_{\Lambda_n}(x,x) - G_{\Z^d}(x,x) = \Theta(n^{2-d}) < 0 \quad \text{uniformly for large $n$}. \end{align} \end{enumerate} \end{theo} Finally, we discuss the eigenvalues of the submatrices of $G^{0\avg}_{\Lambda_n}$. But first recall that its own eigenvalues are given by $0$ and $(1/\eta_k)_{k= 2, \dots, n^d}$, where $0 = \eta_1 < \eta_2 \leq \eta_3\leq \dots \leq \eta_{n^d}$ are the eigenvalues of $-\Delta_{\Lambda_n}$. \begin{lemm}[Eigenvalues of zero-average Green's function] Let $d\geq 3$, $n \in \N$, and $U$ be a finite non-empty subset of $\Z^d$ viewed as a subset of $\Lambda_n$ for each $n$. Set $u_n \coloneqq \lvert U^c\rvert$ and let $(\mu_k)_{k\leq u_n}$ ordered non-increasingly denote the eigenvalues of \[ G_{\Lambda_n}^{0\avg}\vert_{U^c\times U^c} - G_{\Lambda_n}^{0\avg}\vert_{U^c\times U} G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1} G_{\Lambda_n}^{0\avg}\vert_{U\times U^c}. \] Then, for all $k \leq u_n - \lvert U\rvert$, \begin{equation}\label{eq:auxiliary_estim} \mu_k \in \left[\frac{1}{\eta_{n^d-u_n + k + \lvert U\rvert +1}}, \frac{1}{\eta_{k+1}}\right] \quad \left(=\left[0, \frac{1}{\eta_{k+1}}\right] \;\text{for } k=u_n - \lvert U\rvert\right). \end{equation} \end{lemm} \begin{proof} Let $(\zeta_k)_{k\leq u_n}$ and $(\upsilon_k)_{k\leq u_n}$ ordered non-increasingly denote the eigenvalues of the matrices $G_{\Lambda_n}^{0\avg}\vert_{U^c\times U^c}$ and $G_{\Lambda_n}^{0\avg}\vert_{U^c\times U} G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1} G_{\Lambda_n}^{0\avg}\vert_{U\times U^c}$, respectively. Since $G_{\Lambda_n}^{0\avg}\vert_{U^c\times U^c}$ is a principal submatrix of $G_{\Lambda_n}^{0\avg}$, by~\cite[Theorem~4.3.28]{horn2012matrix}, for all $1 \leq k \leq u_n - 1$ \begin{align*} \frac{1}{\eta_{n^d-u_n + k + 1}} \leq \zeta_k \leq \frac{1}{\eta_{k+1}}, \qquad 0 \leq \zeta_{u_n} \leq \frac{1}{\eta_{u_n+1}}. \end{align*} Furthermore, since $G_{\Lambda_n}^{0\avg}\vert_{U^c\times U} G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1} G_{\Lambda_n}^{0\avg}\vert_{U\times U^c}$ is positive semi-definite and has rank at most $\lvert U\rvert$, $\upsilon_k = 0$ for all $k > \lvert U\rvert$. Thus, by Weyl's inequality~\cite[Theorem~4.3.1]{horn2012matrix}, $\mu_k \in [\zeta_{k+ \lvert U\rvert}, \zeta_k]$ for all $k \leq u_n - \lvert U\rvert$. \eqref{eq:auxiliary_estim} follows. \end{proof} \section{The infinite-volume limit of the Spherical model}\label{sec:sphr} The main objective of this section is to characterize the infinite-volume limit of the spherical model, in particular, prove Theorem~\ref{thm:spherical}. We first establish convergence in law: \begin{theo}[Infinite-volume limit of the spherical model]\label{thm:spher_detailed} Let $d\geq 2$, $\Lambda_n = \T^d_n$, and $\theta = (\theta_x)_{x \in \Lambda_n}$ be a configuration of the spherical model on $\Lambda_n$ at inverse temperature $\beta > 0$. Set $\beta_c \coloneqq G_{\Z^d}(0,0)$ (in particular, $\beta_c = \infty$ if $d= 2$ and\pagebreak{} $\beta_c \in (0,\infty)$ otherwise). For any finite set $U \subset \Z^d$ (considered as a subset of $\Lambda_n$ for each $n$), $\theta_{U} \coloneqq (\theta_x)_{x \in U}$ converges in law as $n\to \infty$ to: \begin{enumerate} \item\label{theo3.1.1} $\beta < \beta_c$: an $m^2$-massive GFF on $\Z^d$ restricted to $U$ scaled by $1/\sqrt{\beta}$ with the mass $m^2$ depending on $\beta$ and $d$ in such a way that \begin{align*} G_{\Z^d, m^2}(x,x) = \beta \quad \text{for all } x\in \Z^d; \end{align*} \item\label{theo3.1.2} $\beta = \beta_c$: a GFF on $\Z^d$ restricted to $U$ scaled by $1/\sqrt{\beta}$; \item\label{theo3.1.3} $\beta > \beta_c$: a GFF on $\Z^d$ restricted to $U$ scaled by $1/\sqrt{\beta}$ plus an independent constant random drift of the form $\sqrt{\frac{\beta - \beta_c}{\beta}} X \bfun_U$ with $X$ being a Rademacher random variable, i.e., $\P[X=1] = \P[X=-1] = \frac{1}{2}$. \end{enumerate} \end{theo} The proof is divided into three parts, corresponding to the different temperature regimes, which are discussed separately in the following three subsections. Besides convergence in law, we further obtain convergence of local covariance functions of the spherical model. More precisely, in Section~\ref{subsec:spher_correlations} we prove the following, which together with Theorem~\ref{thm:spher_detailed} yields Theorem~\ref{thm:spherical}. \begin{coro}\label{cor:spher_correlations} Let $\Lambda_n$, $\theta = (\theta_x)_{x \in \Lambda_n}$, and $U \subset \Z^d$ be as above, $(i_x)_{x \in U} \in \N_0^{U}$ be an arbitrary but fixed vector of non-negative integers. Then, the following holds: \begin{align*} \nu_{\Lambda_n, \beta} \left[\prod_{x \in U} \theta_x^{i_x}\right] \xrightarrow{n \to \infty} \E\left[\prod_{x \in U} \alpha_x^{i_x}\right], \end{align*} where $\sqrt{\beta} \alpha \coloneqq \sqrt{\beta}(\alpha_x)_{x \in \Z^d}$ in correspondence with Theorem~\ref{thm:spher_detailed} is either an $m^2$-massive GFF on $\Z^d$ if $\beta < \beta_c$, a GFF on $\Z^d$ if $\beta = \beta_c$, or a GFF on $\Z^d$ plus an independent drift $\sqrt{\beta - \beta_c} X \bfun_{\Z^d}$, where $X$ is a Rademacher random variable, if $\beta > \beta_c$. \end{coro} \subsection{The high-temperature regime}\label{subsec:sph_Thigh} This section presents a proof of $\beta < \beta_c$ case of Theorem~\ref{thm:spher_detailed}. We start by recalling that by Proposition~\ref{prop:rel_GFF_spher}, \begin{align*} \Law\left(\theta \right) = \Law\left(\frac{\phi}{\sqrt{\beta}}\,\middle|\,\ntwo{\phi} = \sqrt{\beta n^d}\right) \end{align*} for any $m_n^2$-massive GFF $\phi$ on $\Lambda_n$ with $m_n^2>0$. Thus, for any finite set $U \subset \Z^d$, we have the following equality of Lebesgue densities on $\R^{U}$: \begin{align}\label{eq:dens_spher} f_{\theta_{U}}(v) &= f_{\frac{\phi_U}{\sqrt{\beta}} \big\vert \ntwo{\phi} =\sqrt{\beta n^d}} (v) = f_{\frac{\phi_U}{\sqrt{\beta}}} \left(v \right) \frac{f_{\|{\phi}^{U^c} + h^{\sqrt{\beta} v}\|^2} (\beta n^d)}{f_{\ntwo{\phi}^2} (\beta n^d)} \eqqcolon f_{\frac{\phi_U}{\sqrt{\beta}}} \left(v \right) {J_n}(v), \end{align} where $v \coloneqq (v_x)_{x \in U} \in \R^{U}$. {In the second equality, we applied the domain Markov property of $\phi$ (Proposition~\ref{prop:domainMP}).} Provided that $m_n^2 \to m^2 > 0$, the vector $\phi_U$ converges in distribution to an $m^2$-massive GFF on $\Z^d$ restricted to $U$. This in particular implies pointwise convergence of the corresponding Gaussian densities. To further ensure that the ratio $J_n(v)$ in~\eqref{eq:dens_spher} converges to $1$, we need to get rid of conditioning on the norm of $\phi$ in the limit. Since $\E[\norm{\phi}^2] = G_{\Lambda_n, m_n^2}(0,0) n^d$, a natural choice would be to select $m_n^2, m^2 > 0$ such that $\lim_n G_{\Lambda_n, m_n^2}(0,0) = G_{\Z^d, m^2}(0,0) = \beta$, which is possible since $\beta < \beta_c = G_{\Z^d}(0,0)$. We may then pick a sequence $(m_n^2)_n$ such that $G_{\Lambda_n, m_n^2}(0,0) = \beta$ for all $n$ sufficiently large. To conclude the proof in this regime, we show that with the above choice of $(m^2_n)_n$, the ratio $J_n(v)$ converges to one. To this end, let $P \in \R^{U^c\times U^c}, Q \in \R^{\Lambda_n\times \Lambda_n}$ be the orthonormal matrices diagonalizing $G_{U^c, m_n^2}\vert_{U^c \times U^c}$ and $G_{\Lambda_n, m_n^2}$, respectively. In particular, \begin{align*} P^T \left(G_{U^c, m_n^2}\vert_{U^c \times U^c}\right) P &= \diag\left((\tau_i)_{i=1}^{\lvert U^c\rvert} \right), \\ Q^T G_{\Lambda_n, m_n^2} Q &= \diag \left((\lambda_i \coloneqq (\eta_i + m_n^2)^{-1})_{i=1}^{n^d}\right) \end{align*} where $\tau_1 \geq \dots \geq \tau_{\lvert U^c\rvert} > 0$ are the eigenvalues of $G_{U^c, m_n^2}\vert_{U^c \times U^c}$, \[ 0 = \eta_1 < \eta_2 \leq \dots \leq \eta_{n^d} \] are those of $-\Delta_{\Lambda_n}$ (cf. Lemma~\ref{lemma:massive_EVs}). Set $u_n \coloneqq \lvert U^c\rvert$, $\tilde h \coloneqq P^T h^{\sqrt{\beta} v}_{U^c}$, and for\linebreak a given triangle array $(Y_{i, u_n})_{i \leq u_n, u_n \in \N}$ of i.i.d. standard normal random variables, define $Z_{i, u_n} \coloneqq (\tilde h_i + \sqrt{\tau_i} Y_{i, u_n})^2 - \tau_i - \tilde h_i^2$. In terms of these variables, we can then rewrite $J_n(v)$ as follows: \begin{align*} J_n(v) &= \frac{f_{\left\|\tilde h + \cN\left(0, \diag(\tau_j)_{j=1}^{u_n}\right)\right\|^2} \left(\beta n^d - \beta \ntwo{v}^2\right)}{f_{\left\|\cN\left(0, \diag(\lambda_j)_{j=1}^{n^d}\right)\right\|^2} (\beta n^d)}\\ &= \frac{f_{\sum_{i=1}^{u_n} Z_{i, u_n}} \left(\beta n^d - \sum_{i=1}^{u_n} \tau_i - \ntwo{\tilde h}^2 - \beta \ntwo{v}^2\right)}{f_{\sum_{i=1}^{n^d} \lambda_i \left(Y_{i,n^d}^2-1\right)} (0)}. \end{align*} Here, we used that $\beta n^d = \Trace(G_{\Lambda_n, m_n^2}) = \sum_{i=1}^{n^d} \lambda_i$. Note that by the canonical properties of the massive Green's function (see Section~\ref{subsec:Greens_fcts}, particularly~\eqref{eq:Gfct_decomp}), we have \begin{multline}\label{eq:trace_comparison} \sum_{i=1}^{u_n} \tau_i + \ntwo{\tilde h}^2\\ \begin{aligned} &= \Trace \bigl(G_{\Lambda_n, m_n^2}\vert_{U^c \times U^c} \bigr) - \sum_{x \in U^c} \E_{\alpha_n}^x\left[\ind_{\{H_{U} < \infty\}} G_{\Lambda_n, m_n^2}(X_{H_{U}},x)\right] + \left\lVert h^{\sqrt{\beta} v}_{U^c}\right\rVert^2\\ &= u_n \beta - \sum_{y \in U} \sum_{x \in U^c} \P_{\alpha_n}^x[H_{U} < \infty, X_{H_{U}} = y] G_{\Lambda_n, m_n^2}(y,x) + \left\|h^{\sqrt{\beta} v}_{U^c}\right\|^2\\ &\eqqcolon n^d \beta - R_{U, n}, \end{aligned} \end{multline} where $\P_{\alpha_n}^x$ denotes the law of a simple random walk $X$ on $\Lambda_n \cup \{\dagger\}$ started at $x$ and killed uniformly with probability $\alpha_n = \frac{m^2_n}{2d + m^2_n}$ at every step; $H_U$ and $\tilde H_U$ are the first hitting time and the first return time to $U$ by this walk, respectively. Furthermore,\pagebreak{} since by the weak Markov property of the (massive) random walk $X$ (here $\omega_k$ stands for the shift of the time by $k$) and~\eqref{eq:Gfct_expdec}, \begin{align*} \P^{x}_{\alpha_n}[H_{U} < \infty] &= \sum_{k \geq 0} \sum_{z \in U} \P^{x}_{\alpha_n}\left[X_k = z,\; \tilde H_{U} \circ \omega_k = \infty\right] \\ &= \sum_{z \in U} \P^{z}_{\alpha_n} \left[\tilde H_{U} = \infty\right] \sum_{k \geq 0} \P^{x}_{\alpha_n}[X_k = z]\\ &\leq \sum_{z \in U} \frac{2d}{1-\alpha_n} G_{\Lambda_n, m^2}(x,z) \leq (2d + m_n^2) |U| \sup_{z \in U} G_{\Lambda_n, m^2}(x,z)\\ &\leq c(m^2, d, U) e^{-c'(m^2, d) \min_{z \in U}|x-z|} \end{align*} for appropriate constants $c, c' > 0$, we obtain \begin{equation}\label{eq:h_norm_bound} \bigl\|h^{\sqrt{\beta} v}_{U^c}\bigr\|^2 \leq \beta \ntwo{v}^2 \sum_{x \in U^c} \P^{x}_{\alpha_n}[H_{U} < \infty]^2 \leq c^2 \beta \ntwo{v}^2 \sum_{x \in U^c} e^{-2c'\min_{z \in U}|x-z|}. \end{equation} And therefore, $0 \leq |R_{U,n}| \leq \beta |U| + \sum_{y \in U} \sum_{x \in U^c} c e^{-c' |x-y|} + \norm{h^{\sqrt{\beta} v}_{U^c}}^2 \leq M$ for an absolute constant $M>0$. Next, we study the variances of the two sums $\sum_{i=1}^{u_n} Z_{i, u_n}$ and $\sum_{i=1}^{n^d} \lambda_i (Y_{i,n^d}^2-1)$, denoted $s_n^2$ and $\sigma_n^2$, respectively. Using Lemma~\ref{lemma:massive_EVs}, similarly to~\eqref{eq:trace_comparison}, we deduce that \begin{align*} \sigma_n^2 = 2 \sum_{k=1}^{n^d} \lambda_k^2 = \Theta(n^d);\qquad s_n^2 = 2\sum_{k=1}^{u_n} \tau_k^2 + \underbrace{4\sum_{j=1}^{u_n} \tau_j \tilde{h}_j^2}_{\leq 4\tau_1 M = \Theta(1)} = \sigma_n^2(1 - o(1)). \end{align*} This allows us to rewrite the ratio in question as \begin{align*} J_n(u) = (1 + o(1)) \frac{s_n f_{\sum_{i=1}^{u_n} Z_{i, u_n}} \left(s_n \frac{R_{U, n} - \beta \ntwo{v}^2}{s_n}\right)}{\sigma_n f_{\sum_{i=1}^{n^d} \lambda_i \left(Y_{i,n^d}^2-1\right)} (0)}. \end{align*} To conclude the proof, we use the following result (stated in a more general form to also cover the critical case), proven in Appendix~\ref{A:proof_densconv_1D}. \begin{prop}\label{prop:dens_conv_1D} Let $(X_{i,n})_{i\leq n, n \in \N}$ be a triangular array of independent centered random variables with the probability density functions $f_{i,n}$ (with respect to Lebesgue measure on $\R$). Assume that \begin{enumerate} \item\label{prop3.3.1} $f_{i,n} \in L^r(\R)$ for some $r \in (1,2]$ (independent of $i,n$) such that $\sup_{i, n} \norm{f_{i,n}}_{L^r} \leq M$ for some $M> 0$; \item\label{prop3.3.2} For all $i \leq n, n \in \N$, $\sigma_{i,n}^2 \coloneqq \Var[X_{i,n}] < \infty$ ordered such that $\sigma_{1,n}^2 \geq \dots \geq \sigma_{n,n}^2$; \item\label{prop3.3.3} Lindeberg's condition is satisfied, that is, for any $\varepsilon > 0$, \begin{align*} \frac{\sum_{i \geq 1} \E\left[X^2_{i,n} \ind_{\{ \lvert X_{i,n}\rvert > \varepsilon s_n\}}\right]}{s_n^2} \xrightarrow{n \to \infty} 0, \end{align*} where $s_n^2 \coloneqq \sum_{i=1}^n \sigma_{i,n}^2$. \item\label{prop3.3.4} There exist $\delta > 0$ (uniform), $K(n) \geq 1$, $l_*(n) \geq 1, n \geq l^*(n) \geq 2 \lceil\frac{r}{r-1}\rceil$ such that for all $n$ sufficiently large \begin{align} \tag{$a$}\label{prop3.3.4.a} \frac{\sum_{i \geq l_*(n)} \sigma_{i,n}^2}{\sum_{i=1}^n \sigma_{i,n}^2} &\geq \delta;\\ \tag{$b$}\label{prop3.3.4.b} \frac{\sum_{i \geq l_*(n)} \E\left[X^2_{i,n} \ind_{\{ \lvert X_{i,n}\rvert > K(n)\}}\right]}{\sum_{i \geq l_*(n)} \sigma_{i,n}^2} &\leq \frac{1}{8}; \\ \tag{$c$}\label{prop3.3.4.c} \frac{n - l^*(n)}{\sigma_{l^*, n}^2 \vee K(n)^2} &\gg \log \left(\sum_{i=1}^n \sigma_{i,n}^2\right). \end{align} \end{enumerate} Then the relation \begin{align*} s_n f^{(n)} (s_n x) \xrightarrow{n\to \infty} \frac{1}{\sqrt{2\pi}} \: e^{-\frac{x^2}{2}} \end{align*} holds uniformly over $\R$. Here, $f^{(n)}$ is the convolution of $f_{1,n}, \dots, f_{n,n}$, as well as the density function of $\sum_{i=1}^n X_{i,n}$. Note that $s_n f^{(n)} (s_n x)$ is the density of $\frac{\sum_{i=1}^n X_{i,n}}{s_n}$. \end{prop} Indeed, if the arrays $(Z_{i, u_n})_{i,u_n}$ and $(\lambda_i (Y_{i,n^d}^2-1))_{i,n^d}$ satisfy the assumptions of the proposition, then the densities \[ s_n f_{\sum_{i=1}^{v_n} Z_{i, v_n}} (s_n x)\quad\text{and}\quad\sigma_n f_{\sum_{i=1}^{n^d} \lambda_i \left(Y_{i,n^d}^2-1\right)} (\sigma_n x) \] both converge to the standard normal density $f_{\cN(0,1)}(x)$, uniformly in $x\in \R$. In particular, this implies $J_n(v) \to 1$, as desired. We now proceed to the remaining step, namely to verify that, up to reordering, the array $(Z_{i, u_n} = (\tilde h_i + \sqrt{\tau_i} Y_{i, u_n})^2 - \tau_i - \tilde h_i^2)_{i,u_n}$ (and fully analogously $(\lambda_i (Y_{i,n^d}^2-1))_{i,n^d}$)\linebreak satisfies the assumptions of Proposition~\ref{prop:dens_conv_1D} with $l_* = 1$, $l^* = 2 \lceil \frac{r}{r-1}\rceil$, and $K$ sufficiently large but independent of $n$. \begin{proof}[Proof that $(Z_{i, u_n})_{i,u_n}$ satisfies the assumptions of Proposition~\ref{prop:dens_conv_1D}] It is clear that the variables $Z_{i, u_n}$ are centered and independent, with variances \[ \sigma^2_{i,u_n} \coloneqq \Var[Z_{i,u_n}] = 2 \tau_i^2 + 4 \tilde h_i^2 \tau_i \in (0, \infty). \] Let $f_{i,u_n}$ denote the density of $Z_{i,u_n}$. It is given by \begin{align*} f_{i,u_n}(t) &= \frac{\left(f_Y \left(\frac{\sqrt{t + \tau_i + \tilde h_i^2} - \tilde h_i}{\sqrt{\tau_i}}\right) + f_Y \left(\frac{\sqrt{t + \tau_i + \tilde h_i^2} + \tilde h_i}{\sqrt{\tau_i}} \right) \right)}{2\sqrt{\tau_i} \sqrt{t + \tau_i + \tilde h_i^2}} \ind_{\left\{t > - \tau_i - \tilde h_i^2\right\}}, \end{align*} where $Y$ is a standard normal random variable. For any $r \in (1,2)$, we estimate (and thus verify condition~\eqref{prop3.3.1}). \begin{align*} \ntwo{f_{i,u_n}}_{L^r}^r &= \int_0^\infty \frac{\left(f_Y \left(\frac{\sqrt{s} - \tilde h_i}{\sqrt{\tau_i}}\right) + f_Y \left(\frac{\sqrt{s} + \tilde h_i}{\sqrt{\tau_i}} \right) \right)^r}{2^r \tau_i^{r/2} s^{r/2}} \d s \leq \int_0^\infty \frac{f^r_Y \left(\frac{\sqrt{s} - \left\lvert \tilde h_i\right\rvert}{\sqrt{\tau_i}} \right)}{\tau_i^{r/2} s^{r/2}} \d s\\ &= \int_0^\infty \frac{f_Y^r \left(\sqrt{s} - \left\lvert\frac{\tilde h_i}{\sqrt{\tau_i}}\right\rvert\right)}{\tau_i^{r-1} s^{r/2}} \d s = 2\int_0^\infty \frac{f_Y^r \left(t - \left\lvert\frac{\tilde h_i}{\sqrt{\tau_i}}\right\rvert\right)}{\tau_i^{r-1} t^{r-1}} \d t\\ &\leq \frac{2}{\tau_i^{r-1}} \left(\int_{0}^1 t^{1-r} \d t + \int_1^{\infty} f_Y \left(t - \left\lvert\frac{\tilde h_i}{\sqrt{\tau_i}} \right\rvert \right) \d t\right)\\ &\leq \frac{2}{\tau_i^{r-1}} \left(\frac{1}{2-r} + 1 \right) \eqqcolon c(r) \frac{1}{\tau_i^{r-1}}. \end{align*} Hence, $\sup_{i, n} \norm{f_{i,u_n}}_{L^r}^r \leq c(r) (4d + 2 m^2)^{r-1}$, by uniform bounds on $(\tau_i)_i$—$\tau_i \in [L_*, L^*] \coloneqq [(2m^2+4d)^{-1}, 2/m^2]$—implied by~\eqref{eq:unif_bdsEVs} and the convergence $m_n^2 \xrightarrow{n\to \infty} m^2>0$. Since by~\eqref{eq:h_norm_bound}, \begin{align*} \tilde h_i^2 \leq \bigl\|h^{\sqrt{\beta} v}_{U^c}\bigr\|^2 \leq c(m^2, d, U) \beta \ntwo{v}^2, \end{align*} we further get that \[ 0 < 2 L_*^2 \leq \sigma^2_{i,u_n} \leq 2(L^*)^2 + 4c \beta\ntwo{v}^2 L^* < \infty \] uniformly over $(i, n)$ (condition~\eqref{prop3.3.2} up to reordering). Consequently, $s_n^2 = \Theta(u_n)$. We now check condition~\eqref{prop3.3.4} (and~\eqref{prop3.3.3}): \eqref{prop3.3.4.a} and~\eqref{prop3.3.4.c} are immediate with our choice of $l_*, l^*$, and large $K$ independent of $n$. For all $K$ sufficiently large (e.g., $K > 4(L^* + c \beta\norm{v}^2)$), we have \begin{align*} \bigl\{ \lvert Z_{i,u_n}\rvert > K \bigr\} &\subset \Bigl\{\left(\tilde h_i + \sqrt{\tau_i} Y_{i,u_n}\right)^2 > K\Bigr\} \\ &\subset \Bigl\{\left|\sqrt{\tau_i} Y_{i,u_n}\right| > \sqrt{K} - |\tilde h_i| \Bigr\} \subset \left\{ \lvert Y_{i,u_n} \rvert > \frac{\sqrt{K}}{2 \sqrt{L^*}} \right\}. \end{align*} Furthermore, $Z^2_{i,u_n} \leq \tilde c(d, m^2, U, L^*, v) (Y_{i, u_n}^4 + \lvert Y_{i, u_n}\rvert^3 + Y_{i, u_n}^2 + \rvert Y_{i, u_n}\rvert + 1)$ for an appropriate $\tilde c > 0$, so \begin{align*} \frac{\E\left[Z^2_{i,u_n} \ind_{\left\{ \left\lvert Z_{i,u_n}\right\rvert > K\right\}}\right]}{ \sigma_{i,u_n}^2} &\leq \frac{5\tilde c}{2L^2_*} \E\left[Y_{i, u_n}^4 \ind_{\left\{ \lvert Y \rvert > \frac{\sqrt{K}}{2 \sqrt{L^*}} \right\}}\right] = \hat{c} \E\left[Y^4 \ind_{\left\{ \lvert Y \rvert > \frac{\sqrt{K}}{2 \sqrt{L^*}} \right\}}\right]. \end{align*} This in turn yields, \begin{align*} \frac{\sum_{i=1}^{u_n} \E[Z^2_{i,u_n} \ind_{\{ \lvert Z_{i,u_n}\rvert > K\}}]}{s_n^2} \leq \hat{c}\E\biggl[Y^4 \ind_{\bigl\{ \lvert Y \rvert > \frac{\sqrt{K}}{2 \sqrt{L^*}} \bigr\}}\biggr], \end{align*} which can be made arbitrarily small by choosing $K$ large. The Lindeberg condition follows with $\varepsilon s_n$ for $\varepsilon > 0$ in place of $K$. \end{proof} \begin{rema} In the high-temperature regime, one could alternatively deduce the uniform convergence of the numerator and denominator in $J_n$ to the standard normal density using results from Malliavin calculus~\cite[Corollary~6.6, Theorem~7.3, Lemma~7.1]{HU2014814}. However, these results are not applicable in the supercritical regime and are of limited use in the critical case (helpful only for $d \geq 4$, see Remark~\ref{rem:Malliavin}). We therefore opt for a more direct approach. \end{rema} \subsection{The low-temperature regime}\label{subsec:sph_Tlow} This section proves the $\beta > \beta_c$ case of Theorem~\ref{thm:spher_detailed} and lays the groundwork for the critical regime. Throughout, we assume $d \geq 3$ and $\beta \geq \beta_c$, explicitly indicating when the analysis is restricted to $\beta > \beta_c$. Recall that by Proposition~\ref{prop:rel_GFF_spher}, \begin{align*} \Law\left(\theta \right) = \Law\left(\frac{\gamma + c \bfun_{\Lambda_n}}{\sqrt{\beta}} \,\middle|\, \ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d}\right), \end{align*} where $\gamma$ is a zero-average GFF on $\Lambda_n$ and $c$ is an independent Lebesgue constant. Our strategy is to analyze the conditional joint law of $(\gamma_U, c)$ given $\norm{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d}$, for a finite set $U \subset \Z^d$ (viewed as a subset of $\Lambda_n$). More precisely, we show that it converges to the product measure of the GFF $\psi$ on $\Z^d$ restricted to $U$ and a Rademacher random variable $X$ scaled by $\sqrt{\beta - \beta_c}$. This immediately yields the desired local convergence in distribution of $\theta$ towards $\frac{1}{\sqrt{\beta}} \psi + \sqrt{\frac{\beta - \beta_c}{\beta}} X \bfun_{\Z^d}$. Let $B \subset \R^{U}$ be a Borel set and $I \subset [0,\sqrt{\beta}]$ be a proper interval with the closure $\bar{I} = [a,b]$ and such that $a,b \notin \{\pm \sqrt{\beta - \beta_c}\}$. Observe that since $\gamma$ is a zero-average GFF, $\langle\gamma, \bfun_{\Lambda_n}\rangle = 0$ almost surely. Therefore, \begin{multline*} \P \Bigl[\gamma_U \in B, c \in I\Bigm|\ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d}\Bigr] = \P \Bigl[\gamma_U \in B, c \in I \Bigm| \ntwo{\gamma} ^2 + c^2 n^d = \beta n^d\Bigr]\\ \begin{aligned} &= \int_B \frac{\int_I f_{\gamma_U, \ntwo{\gamma}^2}\left(y, n^d (\beta - c_0^2)\right) \d c_0 }{\int_{-\sqrt{\beta}}^{\sqrt{\beta}} f_{\ntwo{\gamma}^2} \left(n^d(\beta - \tilde c^2)\right) \d \tilde c} \d y= \int_B f_{\gamma_U}(y) \frac{\int_I f_{\ntwo{\hat{\gamma}^y}^2}\left(n^d (\beta - c_0^2)\right) d c_0 }{2\int_{0}^{\sqrt{\beta}} f_{\ntwo{\gamma}^2} \left(n^d(\beta - \tilde c^2)\right) \d \tilde c} \d y \\ &=\mkern-5mu \int_B f_{\gamma_U}(y) \frac{\int_{a^2}^{b^2} \frac{1}{2 \sqrt{u}} n^d f_{\ntwo{\hat{\gamma}^y}^2 }\left(n^d (\beta - u)\right) \d u}{2 \int_{0}^{\beta} \frac{1}{2 \sqrt{\tilde u}} n^d f_{\ntwo{\gamma}^2} \left(n^d(\beta - \tilde u)\right) \d \tilde u} \d y =\mkern-5mu \int_B f_{\gamma_U}(y) \frac{\int_{a^2}^{b^2} \frac{1}{\sqrt{u}} f_{\frac{\ntwo{\hat{\gamma}^y}^2}{n^d}}(\beta - u) \d u}{2 \int_{0}^{\beta} \frac{1}{\sqrt{\tilde u}} f_{\frac{ \ntwo{\gamma}^2}{n^d}} (\beta - \tilde u) \d \tilde u} \d y \\ &= \int_B f_{\gamma_U}(y) \frac{\E \left[\frac{1}{\sqrt{\beta - \ntwo{\hat{\gamma}^y}^2/ n^d}} \ind_{(\beta - b^2, \beta - a^2)} \left(\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} \right)\right]}{2\E \left[\frac{1}{\sqrt{\beta - \ntwo{\gamma}^2/n^d}} \ind_{[0, \beta)} \left(\frac{\ntwo{\gamma}^2}{n^d} \right)\right]}\d y =: \int_B f_{\gamma_U}(y) R_n(y,I) \d y, \end{aligned} \end{multline*} where $\hat{\gamma}^y$ is defined by $\hat{\gamma}^y_U \equiv y$ and \begin{multline}\label{eq:cond_GFFlaw} \hat{\gamma}^y_{U^c}\sim \cN \left( \begin{aligned} &G_{\Lambda_n}^{0\avg}\vert_{U^c\times U} \;G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1} y, \; G_{\Lambda_n}^{0\avg}\vert_{U^c\times U^c}\\ &\mkern50mu - G_{\Lambda_n}^{0\avg}\vert_{U^c\times U}\; G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1}\; G_{\Lambda_n}^{0\avg}\vert_{U\times U^c} \end{aligned} \right)\\ \eqqcolon \cN (\nu(y), \Xi) = \Law (\gamma_{U^c}\mid\gamma_U = y). \end{multline} Note that we may assume $B$ is bounded. Indeed, by the symmetries of the torus, \begin{align*} \E\left[\ntwo{\gamma_U}^2\,\middle|\,\ntwo{\gamma}^2 + c^2 n^d = \beta n^d \right] = \lvert U\rvert \:\E\left[\gamma_0^2\,\middle|\,\ntwo{\gamma}^2 + c^2 n^d = \beta n^d \right] \leq \lvert U\rvert \beta. \end{align*} Therefore, by choosing $t > 0$ large enough, the probability \[ \P \left[\gamma_U \notin B^{\lvert U\rvert}(0,t)\,\middle|\,\ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d} \right] \] can be made arbitrarily small. To conclude the proof of the low-temperature regime, it suffices to verify the following claim. \begin{enonce}{Claim} \label{claim:ratio_conv} {For any bounded Borel set $B\subset \R^U$, \begin{align}\label{eq:ratio_bound} \limsup_{n\to \infty} \;\sup_{y\in B} R_n(y,I) < \infty, \end{align} and, moreover, for any $y\in B$, } \begin{equation}\label{eq:ratio_convergence} R_n(y,I) = \frac{\E \left[\frac{1}{\sqrt{\beta - \ntwo{\hat{\gamma}^y}^2/ n^d}} \ind_{(\beta - b^2, \beta - a^2)} \left(\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} \right)\right]}{2\E \left[\frac{1}{\sqrt{\beta - \ntwo{\gamma}^2/n^d}} \ind_{[0, \beta)} \left(\frac{\ntwo{\gamma}^2}{n^d} \right)\right]} \xrightarrow{n \to \infty} \frac{1}{2} \ind_{I}\left(\sqrt{\beta - \beta_c}\right). \end{equation} \end{enonce} Indeed, let $s \in \R \setminus\{\pm \sqrt{\beta - \beta_c}\}$\footnote{These are the continuity points of the cumulative distribution function of $\sqrt{\beta - \beta_c} X$.} and set $I_+ = (-\infty, s] \cap (0,\sqrt{\beta}]$, $I_- = (-\infty, s] \cap [-\sqrt{\beta}, 0]$. Since $\gamma_U$ converges in law towards $\psi_U$, which in particular implies uniform convergence of their densities, Claim~\ref{claim:ratio_conv}, together with dominated convergence theorem, yields \begin{multline}\label{eq:concl_converg_prod_meas} \P \left[\gamma_U \in B, c \leq s \,\middle|\,\ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d}\right] \\ \begin{aligned} &\overset{\phantom{{n \to \infty}}}{=} \P \left[\gamma_U \in B, c \in I_+\,\middle|\,\ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d}\right] \ind_{\{I_+ \neq \emptyset\}} \\ &\qquad+ \P \left[\gamma_U \in B, c \in -I_-\,\middle|\,\ntwo{\gamma + c \bfun_{\Lambda_n}} = \sqrt{\beta n^d}\right]\ind_{\{I_- \neq \emptyset\}}\\ &= \int_B f_{\gamma_U}(y) (R_n(y,I_+) + R_n(y,-I_-))\d y \\ &\xrightarrow{n \to \infty} \frac{1}{2}\left(\ind_{I_+}\left(\sqrt{\beta - \beta_c}\right) + \ind_{-I_-}\left(\sqrt{\beta - \beta_c}\right)\right) \P[\psi_U \in B]\\ &= \P\left[\sqrt{\beta - \beta_c} X \leq s\right] \P[\psi_U \in B], \end{aligned} \end{multline} where $\ind_\emptyset \equiv 0$. It now only remains to prove Claim~\ref{claim:ratio_conv} for $\beta >\beta_c$. To this end, we establish the following auxiliary results (two of which also hold for $\beta = \beta_c$): \begin{enonce*}[defi]{Claims} \leavevmode %\advance\leftmargini by 2.4\normalparindent \begin{enumerate}\claimenumi \item\label{claim1} For any $y \in B$, $\frac{\norm{\gamma}^2}{n^d}$ and $\frac{\norm{\hat{\gamma}^y}^2}{n^d}$ converge almost surely to $\beta_c$. \item\label{claim2} If $\beta > \beta_c$, there exist $c = c(\beta, d), \hat{c} = \hat{c}(\beta, d, B) > 0$ such that for any $0 < \varepsilon < (\beta - \beta_c)/4$ and $y \in B$, \begin{align*} \P\left[\frac{\ntwo{\gamma}^2}{n^d} \geq \beta - \varepsilon\right] \leq e^{-c n^{d-2}}, \qquad \P\left[\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} \geq \beta - \varepsilon \right] \leq e^{-\hat{c} n^{d-2}}. \end{align*} \item\label{claim3} There exists an absolute constant $K> 0$ such that the density functions of $\norm{\gamma}^2$ and $\norm{\hat{\gamma}^y}^2$ for any $y\in \R^U$ are uniformly bounded by $K$, i.e., for all $n$ sufficiently large, \begin{align*} \sup_{y\in \R^U}\left\|f_{\ntwo{\hat{\gamma}^y}^2}\right\|_{L^\infty} \vee \left\|f_{\ntwo{\gamma}^2}\right\|_{L^\infty} \leq K. \end{align*} \end{enumerate} \end{enonce*} To motivate these claims, observe that in view of~\eqref{eq:ratio_convergence}, the convergence of $R_n(y,I)$ crucially depends on the behavior of the squared norms $\norm{\gamma}^2/n^d$ and $\norm{\hat{\gamma}^y}^2/n^d$. Claim~\ref{claim1} captures the expected concentration of both quantities around $\beta_c$, thereby providing a heuristic justification for the limiting behavior of the ratio $R_n(y,I)$. To make this precise, however, we must control the contributions from the upper end of the integration domain, specifically the interval $(\beta - \varepsilon, \beta)$, where the integrand may become singular due to the factor $1/\sqrt{\beta - x}$. A direct route would have been to prove uniform boundedness of the density of $\norm{\hat{\gamma}^y}^2/n^d$ on this interval, but this is a highly non-trivial and technically involved fact. Instead, we split $(\beta - \varepsilon, \beta)$ into two subintervals: one where the contribution is negligible due to exponential concentration away from $\beta$ (provided by Claim~\ref{claim2}), and another, sufficiently small neighborhood of $\beta$, where we can appeal to the integrability of $x \mapsto 1/\sqrt{\beta - x}$ and use boundedness of the density function of $\norm{\hat{\gamma}^y}^2$ (ensured by Claim~\ref{claim3}). \begin{proof}[Proof of Claim~\ref{claim:ratio_conv} assuming Claims~\ref{claim1},\ref{claim2} and \ref{claim3}] Let $\varepsilon > 0$ be small enough (as in Claim~2). Since $\R_+ \ni x \mapsto 1/\sqrt{x}$ away from zero is bounded, by dominated convergence theorem and Claim~1 we get \begin{equation}\label{eq:aux_conv_ratioR_n} \begin{aligned}[t] &\E \left[\frac{1}{\sqrt{\beta - \ntwo{\hat{\gamma}^y}^2/ n^d}} \ind_{(\beta - b^2, (\beta - a^2) \wedge (\beta - \varepsilon))} \left(\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} \right)\right]\\ &\mkern260mu\xrightarrow{n \to \infty} \frac{1}{\sqrt{\beta - \beta_c}} \ind_{(a,b)}\left(\sqrt{\beta - \beta_c}\right); \\ &\E \left[\frac{1}{\sqrt{\beta - \ntwo{\gamma}^2/n^d}} \ind_{[0, \beta - \varepsilon)} \left(\frac{\ntwo{\gamma}^2}{n^d} \right)\right]\xrightarrow{n \to \infty} \frac{1}{\sqrt{\beta - \beta_c}}. \end{aligned} \end{equation} %\goodbreak We now consider the remaining bits and show that they are vanishing in the limit. Let $\delta > 0$ be small, \begingroup \allowdisplaybreaks \begin{align*} 0 &\leq \E \left[\frac{1}{\sqrt{\beta - \ntwo{\gamma}^2/n^d}} \ind_{[\beta - \varepsilon, \beta)} \left(\frac{\ntwo{\gamma}^2}{n^d} \right)\right] \eqqcolon r_\varepsilon(\gamma)\\ &\leq \frac{n^d}{\sqrt{\delta}} \P \left[\frac{\ntwo{\gamma}^2}{n^d} \in \left[\beta - \varepsilon, \beta - \delta/n^{2d}\right) \right] +\mkern-5mu \sqrt{n^d} \E \left[\frac{1}{\sqrt{\beta n^d - \ntwo{\gamma}^2}} \ind_{\beta n^d + [- \delta/n^d, 0)} (\ntwo{\gamma}^2)\right]\\ &\leq \frac{n^d}{\sqrt{\delta}} e^{-c n^{d-2}} + \sqrt{n^d}\ntwo{f_{\ntwo{\gamma}^2}}_{L^\infty} \int_0^{\delta/n^d} \frac{\d t}{\sqrt{t}}\leq \frac{n^d}{\sqrt{\delta}} e^{-c n^{d-2}} + 2K \sqrt{\delta}. \end{align*} \endgroup The two inequalities in the last line follow by Claims~\ref{claim2} and~\ref{claim3}, respectively. By taking limit $n \to \infty$, we eliminate the first term in the upper bound; the second one can be made arbitrarily small by adjusting $\delta$. By replacing $\gamma$ with $\hat{\gamma}^y$, we recover the same bound for $r_\varepsilon(\hat{\gamma}^y)$. The above observations(\eqref{eq:aux_conv_ratioR_n} and that $r_\varepsilon(\hat{\gamma}^y)$ and $r_\varepsilon(\gamma)$ can be made arbitrarily small) readily yield~\eqref{eq:ratio_convergence}. Moreover, note that \begin{align*} R_n(y,I) \leq \frac{\frac{1}{\sqrt{\varepsilon}} + r_\varepsilon(\hat{\gamma}^y)}{\frac{2}{\sqrt{\beta}} \P \left[\frac{\ntwo{\gamma}^2}{n^d} < \beta\right]} \leq \sqrt{\beta} \left(\frac{1}{\sqrt{\varepsilon}} + \frac{n^d}{\sqrt{\delta}} e^{-\hat c n^{d-2}} + 2K \sqrt{\delta}\right), \end{align*} where we used Claim~\ref{claim1} in the last inequality to estimate the denominator. The bound~\eqref{eq:ratio_bound} follows. \end{proof} We conclude this subsection and the proof of the low-temperature regime by verifying the remaining claims. \begin{proof}[{Proof of Claims~\ref{claim1} and~\ref{claim2}}] We will show that there exists $c>0$ such that for any $t > 0$ and for all $n$ sufficiently large, \begin{equation} \label{eq:concentration_of_measure} \P[S_n > t] \vee \P[S_n < -t] \leq e^{-c n^{d-2} \left(\sqrt{\beta_c + t/2} - \sqrt{\beta_c}\right)^2}, \end{equation} where $S_n \in \left\{\frac{\norm{\gamma}^2 - \E[\norm{\gamma}^2]}{n^d}, \frac{\norm{\hat{\gamma}^y}^2 - \E[\norm{\hat{\gamma}^y}^2]}{n^d}\right\}$. Together with the observations that \begin{multline}\label{eq:observation} G^{0\avg}_{\Lambda_n}(0,0) \xrightarrow{n \to \infty} \beta_c \\ \text{and}\quad \frac{1}{n^d} \E[\ntwo{\hat{\gamma}^y_{U^c}}^2] = \frac{1}{n^d}\Trace(\Xi) + \frac{1}{n^d} \ntwo{\nu(y)}^2 \xrightarrow{n \to \infty} \beta_c + 0, \end{multline} the inequality~\eqref{eq:concentration_of_measure} implies Claim~\ref{claim1} by Borel--Cantelli lemma (as the estimates are clearly summable in $n$) and Claim~\ref{claim2} by setting $t = \beta - \beta_c - \varepsilon \pm o(1)$. As for~\eqref{eq:observation}, note that by Proposition~\ref{prop:props:zeroAvGfct}, $G_{\Lambda_n}^{0\avg}(x,y)$ converges uniformly over $U^2$ to $G_{\Z^d}(x,y)$, and so do the corresponding eigenvalues. The latter matrix $G_{\Z^d}\vert_{U \times U}$ is positive-definite, thus, there exists $c>0$ (depending only on $U$) such that \[ 00$, \[ \P[f(\bar{Z}) \geq \E[f(\bar{Z})] + x] \leq 2 \frac{1}{\sqrt{2\pi}} \int_{x/\Lip(f)}^\infty e^{-u^2/2} \d u \leq e^{-\frac{1}{2} ({x}/{\Lip(f)})^2}. \] With $\hat{f}^y$ in place of $f$, we obtain an analogous bound for the choice $S_n = \frac{\norm{\hat{\gamma}^y}^2 - \E[\norm{\hat{\gamma}^y}^2]}{n^d}$. To estimate $\P[S_n \leq -t]$, consider $-f$ and $-\hat{f}^y$. \end{proof} \begin{proof}[Proof of Claim~\ref{claim3}] Since $\Xi$ (cf.~\eqref{eq:cond_GFFlaw}) is a covariance matrix, it is positive semi-definite, and we denote its eigenvalues by $\mu_1 \geq \mu \geq\dots \geq \mu_{\lvert U^c\rvert} \geq 0$. We set $u_n \coloneqq \lvert U^c\rvert$ and observe that $\mu_{u_n} = 0$ with the unit eigenvector $p^{u_n} = \frac{1}{\sqrt{u_n}} \bfun_{U^c}$. Let $P = (p^1, \dots, p^{u_n})$ and $Q$ be the orthonormal matrices diagonalizing $\Xi$ and $G_{\Lambda_n}^{0\avg}$, respectively, i.e., \begin{align*} P^T \Xi P = \diag \bigl((\mu_k)_{k=1}^{u_n}\bigr) \quad \text{and}\quad Q^T G_{\Lambda_n}^{0\avg} Q = \diag \bigl(0, (1/\eta_k)_{k=2}^{n^d}\bigr), \end{align*} where $0=\eta_1 < \eta \leq \dots\leq \eta_{n^d}$ are the eigenvalues of $-\Delta_{\Lambda_n}$. Let $(Y_i)_i$ be i.i.d. standard normal random variables. Then, if we set $h(y) \coloneqq P^T \nu(y)$, we obtain the following equalities in law: \begin{equation}\label{eq:EV_representation} \begin{aligned} \ntwo{\gamma} ^2 \overset{\law}{=} \sum_{k=2}^{n^d} \frac{1}{\eta_k} Y_k^2;\qquad \ntwo{\hat{\gamma}^y_{U^c}} ^2 \overset{\law}{=} \sum_{k=1}^{u_n-1} (h_k(y) + \sqrt{\mu_k} Y_k)^2 + h_{u_n}(y)^2. \end{aligned} \end{equation} By Young's convolution inequality and the fact that $\norm{f}_{L^1} = 1$ for any probability density function $f$, we get that \begin{align*} \ntwo{f_{\ntwo{\hat{\gamma}^y}^2}}_{L^\infty} \leq \ntwo{f_{\left(h_l(y) + \sqrt{\mu_l} Y_l\right)^2 + \left(h_i(y) + \sqrt{\mu_i} Y_i\right)^2}}_{L^\infty}, && \ntwo{f_{\ntwo{\gamma} ^2}}_{L^\infty} \leq \ntwo{f_{Y_k^2/\eta_k + Y_j^2/\eta_j}}_{L^\infty} \end{align*} for any $k\neq j, l \neq i$ (the choice will be specified later). To complete the proof, it now suffices to show that there exists an absolute (uniform in $n$ and $y \in \R^U$) constant $K > 0$ such that \[ \ntwo{f_{\left(h_l(y) + \sqrt{\mu_l} Y_l\right)^2 + \left(h_i(y) + \sqrt{\mu_i} Y_i\right)^2}}_{L^\infty} \leq K, \] as the remaining case can be easily deduced from this one by setting $U = \emptyset$. To this end, fix $x \in \R_+$ and, for simplicity, write $(h_l(y) + \sqrt{\mu_l} Y_l)^2 + (h_i(y) + \sqrt{\mu_i} Y_i)^2 = (a + c Y)^2 + (b + d Z)^2$ (with $Y, Z$ independent standard normal random variables). The density of $(a/c + Y)^2$ at $x$ is given by \[ \frac{1}{2 \sqrt{x}} \bigl(f_Y\left(\sqrt{x} - \lvert a/c\rvert\right) + f_Y\left(\sqrt{x} + \lvert a/c\rvert\right)\bigr) \leq \frac{1}{\sqrt{x}}. \] Therefore, \begin{align*} f_{(a + c Y)^2 + (b + d Z)^2}(x) &= \int_0^x f_{c^2(a/c + Y)^2} (z) f_{d^2(b/d + Z)^2} (x-z) \d z \\ &\leq \frac{1}{c d} \int_0^1 \frac{1}{\sqrt{u (1-u)}} \d u = \frac{\pi}{c d}. \end{align*} If $c, d > 0$ are uniformly bounded away from zero, the claim follows. In our setting we can, for example, pick $(l,i) = (\lceil u_{n}/2\rceil, \lceil u_{n}/2\rceil + 1) = (k, j)$, then by~\eqref{eq:EV_asympt} and~\eqref{eq:auxiliary_estim}: $\mu_{l}, \mu_i, 1/ \eta_k, 1/\eta_j = \Theta(1)$ uniformly in $n$. \end{proof} \subsection{At the {critical temperature}}\label{subsec:spher_crit} In this subsection, we analyze the model at the critical temperature $\beta = \beta_c$. Unlike in the other two regimes, the case $d=3$ requires a slightly different treatment compared to higher dimensions. Our proof for $d \geq 4$ combines arguments from both non-critical regimes and makes use of Proposition~\ref{prop:dens_conv_1D} along with~\eqref{eq:0avgGfct_polydec}. For $d=3$, the argument relies heavily on the fact that $\beta_c - G_{\Lambda_n}^{0\avg}(0,0) = \Theta(1/n) > 0$ (see~\eqref{eq:negativity_Greens_diff}), and on the observation that the variance of $\norm{\gamma}^2$ is of the same order as the square of the largest eigenvalue of $G_{\Lambda_n}^{0\avg}(0,0)$. We start by recalling from Section~\ref{subsec:sph_Tlow} that, to prove Theorem~\ref{thm:spher_detailed} at the critical temperature, it suffices to show that the conditional law of $(\gamma_U, c)$, given $\norm{\gamma + c \bfun_{\Lambda_n}}^2 = \beta_c n^d$, converges to $(\psi_U, 0)$, where $\psi$ is the GFF on $\Z^d$. By the argument in~\eqref{eq:concl_converg_prod_meas}, this is in turn implied by the following refined version of Claim~\ref{claim:ratio_conv}, which additionally takes into account intervals $I$ whose boundary contains the discontinuity point $0$ of the constant zero-variable. \begin{enonce}{Claim}\label{claim:refined} Let $I\subset [0,\sqrt{\beta}]$ be a proper interval with $\bar I \eqqcolon [a,b]$, $B \subset \R^U$ be a bounded Borel set, then for any $y\in B$, \begin{align}\label{eq:ratio_critical_convergence} R_n(y,I) = \frac{\E \Big[\frac{1}{\sqrt{\beta_c - \ntwo{\hat{\gamma}^y}^2/ n^d}} \ind_{(\beta_c - b^2, \beta_c - a^2)} \Big(\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} \Big)\Big]}{2\E \Big[\frac{1}{\sqrt{\beta_c - \ntwo{\gamma}^2/n^d}} \ind_{[0, \beta_c)} \Big(\frac{\ntwo{\gamma}^2}{n^d} \Big)\Big]} \xrightarrow{n\to \infty} \frac{1}{2} \ind_{\{a=0\}}, \end{align} and~\eqref{eq:ratio_bound} holds for $R_n(y,I)$. \end{enonce} We first discuss the case $a>0$, in which we can work with the original expression for the ratio $R_n(y,I)$, while relying on two auxiliary results proven first below. Namely, that $\P[\norm{\gamma}^2 < \beta_c n^d]$\footnote{Note that $\P[\norm{\gamma}^2 < \beta_c n^d] = \P[X_n < T_n] \geq \P[X < T]/2$ for all $n$ sufficiently large by Claim~\ref{claimA} and since $T = \lim_n T_n$.} is uniformly bounded away from zero (see the proof of the case $a=0$), i.e., there exists $\delta >0$ such that for all $n$ sufficiently large, $\P[\norm{\gamma}^2 < \beta_c n^d] \geq \delta$, and that \begin{equation}\label{eq:var_error} \begin{aligned}[t] \Var[\ntwo{\hat{\gamma}^y}^2] &= \Var[\ntwo{{\gamma}}^2] + \cO \Bigl(n^2 \bigl(n \ind_{\{d=3\}} + \log n \ind_{\{d= 4\}} + \ind_{\{d \geq 5\}}\bigr) \Bigr)\\ &= (1+o(1)) \Var[\ntwo{{\gamma}}^2] \\ &= \Theta \bigl(n^4 \ind_{\{d=3\}} + n^d \log n \ind_{\{d=4\}} + n^d \ind_{\{d \geq 5\}}\bigr), \end{aligned} \end{equation} \begin{proof}[Proof of Claim~\ref{claim:refined} for $a>0$ assuming~\eqref{eq:var_error} and {$\liminf_n \P[\norm{\gamma}^2 < \beta_c n^d] \geq \delta$}] Observe that \begin{align} \limsup_{n\to\infty} \sup_{y\in B} R_n(y,I) &\leq \frac{\sqrt{\beta_c} \limsup_{n} \sup_{y\in B}\P\bigl[\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} < \beta_c - a^2\bigr]}{2a \liminf_{n} \P\bigl[\frac{\ntwo{{\gamma}}^2}{n^d} < \beta_c - a^2\bigr]} \\ &\leq \frac{\sqrt{\beta_c}}{2a \delta} \limsup_{n\to\infty} \sup_{y\in B} \P\Bigl[\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} < \beta_c - a^2\Bigr]. \end{align} Furthermore, similarly to~\eqref{eq:bound_shift}, we obtain that \begin{align}\label{eq:expect_gamma^y_sqrd} \E[\ntwo{\hat{\gamma}^y}^2] = G^{0\avg}_{\Lambda_n}(0,0) u_n + o(n^{d/2}) \end{align} uniformly in $y\in B$; recall that $u_n = |U^c|$. Therefore, by Chebyshev's inequality and~\eqref{eq:var_error}, we further have \begin{align*} \limsup_n \sup_{y\in B} \P\left[\frac{\ntwo{\hat{\gamma}^y}^2}{n^d} < \beta_c - a^2\right] \leq \limsup_n \sup_{y\in B} \frac{\Var[{\ntwo{\hat{\gamma}^y}^2}]}{n^{2d} (a^2/2)^2} = 0. \end{align*} This proves both~\eqref{eq:ratio_bound} and~\eqref{eq:ratio_critical_convergence}, and thus, completes the proof in the case $a>0$. \end{proof} Observe that the $a>0$ case implies that value of $b \in (a, \sqrt{\beta_c}]$ is irrelevant, thus, without loss of generality, we can assume that $b^2 = \beta_c$. So, hereinafter let $a=0, b^2 = \beta_c$, and write simply $R_n(y)$ for this choice of the interval $I$. Next, we define the following variables: \begin{alignat*}{3} X_n &\coloneqq \frac{\ntwo{\gamma}^2 - \E[\ntwo{\gamma}^2]}{\sqrt{\Var[\ntwo{\gamma}^2]}}, \quad &\hat{X}_n &\coloneqq \frac{\ntwo{\hat{\gamma}^y}^2 - \E[\ntwo{\hat{\gamma}^y}^2]}{\sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]}};\\ T_n &\coloneqq \frac{\beta_c n^d - \E[\ntwo{\gamma}^2]}{\sqrt{\Var[\ntwo{\gamma}^2]}}, \qquad\: &\hat{T}_n &\coloneqq \frac{\beta_c n^d - \E[\ntwo{\hat{\gamma}^y}^2]}{\sqrt{\Var[\ntwo{{\gamma}^y}^2]}}. \end{alignat*} Using them, we can rewrite $R_n(y)$ in view of~\eqref{eq:var_error} as \begin{align*} R_n(y) = \frac{1 + o(1)}{2} \frac{\E \left[\frac{1}{\sqrt{\hat{T}_n - \hat{X}_n}} \ind_{\{\hat{X}_n < \hat{T}_n\}}\right]}{\E \left[\frac{1}{\sqrt{T_n - X_n}} \ind_{\{X_n < T_n\}}\right]}. \end{align*} Moreover, note that~\eqref{eq:expect_gamma^y_sqrd} combined with~\eqref{eq:0avgGfct_polydec} and~\eqref{eq:negativity_Greens_diff} implies that both $T_n$ and $\hat{T}_n$ converge to the same finite limit $T$ (depending on $d$): for $d \geq 4$, $T = 0$; for $d = 3$, $T > 0$. As for $\hat{X}_n$ and $X_n$, we have the following. %\begin{enonce*}[defi]{Claims} %\leavevmode %\begin{enumerate}\Alphenumi \begin{claiml}\label{claimA} $X_n$ and $\hat{X}_n$ converge in law to the same non-trivial real random variable $X$. For $d \geq 4$, $X$ is a standard normal random variable. \end{claiml} \begin{claiml}\label{claimB} There exists a constant $M = M(B,d)>0$ such that, for all sufficiently large $n$ and all $y\in B$, the densities of $\hat{X}_n$ and $X_n$ are bounded by $M$. \end{claiml} \begin{proof}[Proof of Claim~\ref{claim:refined} for $a=0$ assuming~\eqref{eq:var_error} and Claims~\ref{claimA} and~\ref{claimB}]% \!\! Using Claims~\ref{claimA}--\ref{claimB} and applying Portmanteau theorem to the upper/lower semicontinuous functions $x \mapsto \frac{1}{\sqrt{x}} \ind_{\{x \geq \varepsilon\}}$ and $x \mapsto \frac{1}{\sqrt{x}} \ind_{\{x > 0\}}$, respectively, for $\varepsilon > 0$ arbitrary small such that $\varepsilon < T(d=3)/2$, we obtain that \begin{align*} \limsup_{n \to \infty} \E \left[\frac{1}{\sqrt{\hat{T}_n - \hat{X}_n}} \ind_{\{\hat{X}_n < \hat{T}_n\}}\right] &\leq 2M \sqrt{\varepsilon} + \E\left[\frac{1}{\sqrt{T - X}} \ind_{\{X < T - \varepsilon\}}\right];\\ {\limsup_{n \to \infty} \sup_{y\in B} \E \left[\frac{1}{\sqrt{\hat{T}_n - \hat{X}_n}} \ind_{\{\hat{X}_n < \hat{T}_n\}}\right]} &\leq 2M \sqrt{\varepsilon} + \frac{1}{\sqrt{\varepsilon}};\\ \liminf_{n \to \infty} \E \left[\frac{1}{\sqrt{T_n - X_n}} \ind_{\{X_n < T_n\}}\right] &\geq \E\left[\frac{1}{\sqrt{T - X}} \ind_{\{X < T\}} \right]\\ &\geq \E\left[\frac{1}{\sqrt{T - X}} \ind_{\{X < T - \varepsilon\}}\right]. \end{align*} Since $X$ is non-trivial and of mean zero, the probability of it being non-positive is strictly positive. Hence, if $d=3$, since $T - \varepsilon > 0$ in this case, there exists an absolute (independent of $\varepsilon < T/2$) constant $p>0$ such that \[ \E\left[\frac{1}{\sqrt{T - X}} \ind_{\{X < T - \varepsilon\}}\right] \wedge \P[X < T] \geq p. \] Existence of such a constant for $d\geq 4$ even though $T(d\geq 4) = 0$ is clear since $X$ is standard Gaussian. Thus, for all $n$ sufficiently large, \begin{align*} R_n(y) &\leq \frac{1+o(1)}{2} \left(1 + \frac{2M\sqrt{\varepsilon}}{p}\right); \\ {\sup_{y\in B} R_n(y)} &\leq \frac{2M \sqrt{\varepsilon} + \frac{1}{\sqrt{\varepsilon}}}{p}. \end{align*} This proves~\eqref{eq:ratio_bound}; and by taking $\varepsilon>0$ much smaller than $p$, we conclude that $R_n(y)$ is less or equal to $1/2$ in the limit. By reversing the roles of the numerator and denominator, we derive that $R_n(y)$ is bounded from below by $1/2$. \end{proof} The proof of the critical-temperature regime is thus complete, provided~\eqref{eq:var_error} and Claims~\ref{claimA}--\ref{claimB} are verified. We so proceed to discuss the remaining proofs of~\eqref{eq:var_error} and Claims~\ref{claimA}--\ref{claimB}. Our arguments for the latter will significantly differ for $d \geq 4$ and $d=3$. \begin{proof}[Proof of~\eqref{eq:var_error}] Using the eigenvalue representation~\eqref{eq:EV_representation}, we easily observe that $\Var[\norm{\hat{\gamma}^y}^2] = 2 \sum_{k=1}^{u_n-1} (\mu_k^2 + 2 h_k(y)^2 \mu_k)$. By~\eqref{eq:EV_asympt}, \eqref{eq:auxiliary_estim}, and~\eqref{eq:bound_shift}, \[ \mu_1 \ntwo{h(y)}^2 \ll \Theta \bigl(n^4 \ind_{\{d=3\}} + n^d \log n \ind_{\{d=4\}} + n^d \ind_{\{d \geq 5\}}\bigr) = 2 \sum_{k=1}^{u_n-1} \mu_k^2 \] uniformly over $y$ in a fixed bounded set. Combined with~\eqref{eq:0avgGfct_polydec}, this further implies that \begin{align*} \frac{1}{2} \Var[\ntwo{\hat{\gamma}^y}^2] &= \sum_{k=1}^{u_n-1} \mu_k^2 + \cO \bigl(n^2 \bigl(n \ind_{\{d=3\}} + \log n \ind_{\{d= 4\}} + \ind_{\{d \geq 5\}}\bigr) \bigr);\\ \sum_{k=1}^{u_n-1} \mu_k^2 &= \Trace(\Xi^2) = \frac{1}{2} \Var[\ntwo{{\gamma}}^2] + A_1 + A_2 + A_3, \end{align*} where $\hat{\Lambda}_n = [-n/2, n/2)^d \cap \Z^d$ and \begingroup \allowdisplaybreaks \begin{align*} A_1 \,\MK &\!\Mk:= \sum_{x, y \in U^c} G_{\Lambda_n}^{0\avg}(x,y)^2 - \frac{1}{2} \Var[\ntwo{{\gamma}}^2] = - 2 \sum_{\substack{x \in U^c, \\ z \in U}} G_{\Lambda_n}^{0\avg}(x,z)^2 \\ &= \cO \bigl(\bigl(n \ind_{\{d=3\}} + \log n \ind_{\{d= 4\}} + \ind_{\{d \geq 5\}}\bigr) \bigr); \\%HERE (\log n)^{3d} A_2 \,\MK &\!\Mk:= \sum_{\substack{z, z',\\ w, w' \in U\\ x,y \in U^c}} \left( \begin{aligned} &G_{\Lambda_n}^{0\avg}(x,z)\; G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1}(z, z')\; G_{\Lambda_n}^{0\avg}(z',y)\\ &\qquad \times G_{\Lambda_n}^{0\avg}(y,w)\; G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1}(w, w')\; G_{\Lambda_n}^{0\avg}(w',x) \end{aligned} \right)\\ &=\cO\left(\left(\int_{1}^n r^{4-2d} r^{d-1} \d r\right)^2\right) = \cO \bigl(n^2 \ind_{\{d=3\}} + \log^2 n \ind_{\{d= 4\}} + \ind_{\{d \geq 5\}}\bigr);\\ A_3 \,\MK &\!\Mk:= - 2\sum_{\substack{z, z' \in U\\ x,y \in U^c}}\; G_{\Lambda_n}^{0\avg} \vert_{U^2}^{-1}(z, z')\; G_{\Lambda_n}^{0\avg}(x,y) \;G_{\Lambda_n}^{0\avg}(y,z) \;G_{\Lambda_n}^{0\avg}(z',x) \\ &= \cO \left(\sum_{\substack{x\neq y \in \hat{\Lambda}_n \setminus\{0\}}} \lvert x\rvert^{2-d} \lvert y \rvert^{2-d} \d(x,y)^{2-d} \right)\\ &= \cO \left(2 \sum_{x \in \hat{\Lambda}_n \setminus \{0\}} \lvert x\rvert^{4-2d} \sum_{\substack{y \in \hat{\Lambda}_n \setminus \{0,x\}\\ \lvert y\rvert \geq \lvert x\rvert}} \d(x,y)^{2-d} \right) \\ &= \cO \left(n^2 \sum_{x \in \hat{\Lambda}_n \setminus \{0\}} \lvert x\rvert^{4-2d} \right) = \cO \bigl(n^2 \bigl(n \ind_{\{d=3\}} + \log n \ind_{\{d= 4\}} + \ind_{\{d \geq 5\}}\bigr) \bigr). \end{align*} \endgroup In conclusion, as desired \begin{align*} \Var[\ntwo{\hat{\gamma}^y}^2] &= \Var[\ntwo{{\gamma}}^2] + \cO \Bigl(n^2 \bigl(n \ind_{\{d=3\}} + \log n \ind_{\{d= 4\}} + \ind_{\{d \geq 5\}}\bigr) \Bigr)\\ &= (1+o(1)) \Var[\ntwo{{\gamma}}^2].\qedhere \end{align*} \end{proof} \subsubsection*{Proof of Claims~\ref{claimA} and~\ref{claimB} for \texorpdfstring{$d \geq 4$}{d>=4}} The key observation in this case is the following. \begin{claiml}\label{claimC} The densities of $X_n$ and $\hat{X}_n$\footnote{up to removing the last $\lvert U\rvert$ terms in the eigenvalue decomposition~\eqref{eq:EV_representation} of $\hat{X}_n$ containing $\mu_{u_n - q}$ for $1 \leq q \leq \lvert U\rvert$} converge uniformly towards the density function of a standard normal random variable. \end{claiml} The proof of this result is based on Proposition~\ref{prop:dens_conv_1D} (note that for $d=3$ even Lindeberg's condition fails) and the eigenvalue representation~\eqref{eq:EV_representation}. Before proving Claim~\ref{claimC}, we explain how it yields Claims~\ref{claimA}--\ref{claimB}. Note that by the eigenvalue decomposition we can write $\hat{X}_n$ as a sum of two independent random variables $P_n + L_n$, where $P_n$ is the ``good'' part treated in Claim~\ref{claimC} and $L_n$ is the sum of the remaining $\lvert U\rvert$ terms. Then, by Claim~\ref{claimC}, $P_n$ converges in law to a standard normal random variable and its density is uniformly bounded. Note that \[ L_n = \frac{\sum_{u_n - \lvert U\rvert}^{u_n} \mu_k (Y_k^2 - 1) + 2h_k(y) \sqrt{\mu_k} Y_k}{\sqrt{\Var[\norm{\hat{\gamma}^y}^2]}}. \] Since these $\mu_k$'s are uniformly bounded and $\norm{h(y)}$ by~\eqref{eq:bound_shift} is growing slower than any polynomial, it is straightforward to verify that $L_n$ admits a density function and converges almost surely to zero. This implies that $L_n + P_n$ converges in law to a standard Gaussian variable (Claim~\ref{claimA}) and that the density of $\hat{X}_n$ is uniformly bounded as a convolution of a bounded function and an integrable non-negative (with unit $L^1$-norm) function (Claim~\ref{claimB}). \begin{proof}[Proof of Claim~\ref{claimC}] We check that the required assumptions of Proposition~\ref{prop:dens_conv_1D} are satisfied by $\bfY_n \coloneqq ((h_k(y) + \sqrt{\mu_k} Y_k)^2 - \mu_k - h_k(y)^2)_{k=1}^{u_n - \lvert U\rvert - 1}$ contributing to~$\hat X_n$, and thus, also by $(\frac{1}{\eta_k}(Y_k^2 - 1))_{k=2}^{n^d}$ contributing to $X_n$, since the latter can be recovered from the former by setting $h(y) \equiv 0$ and $U = \emptyset$. More precisely, we consider independent families $(\bfY_n)_n$ that form a triangular array. For better readability, we do not use double-indexing for the elements of each $\bfY_n$. The first assumption of the proposition can be directly deduced from the\linebreak corresponding verification in the high-temperature regime since the eigenvalues $(\mu_k)_{k=1}^{u_n-\lvert U\rvert - 1}$ of $\Xi$ and $(1/\eta_k)_{k=2}^{n^d}$ of $G_{\Lambda_n}^{0\avg}$ are uniformly bounded away from zero. The second condition (up to reordering which will be done for the verification of the fourth criterion) is clear since $\sigma^2_{k} = 2 \mu_k^2 + 4 h_k(y)^2 \mu_k < \infty$. From this point on, the index $k$ refers to the original ordering, while $\tilde k$ denotes the ordering by decreasing variance. Recall that Lindeberg's condition (the third condition of the proposition) follows from Lyapunov's condition, which in our setting for $\delta = 2$ appears as \begin{align*} \frac{1}{s_n^4} \sum_{k=1}^{u_n-\lvert U\rvert - 1} \E\left[\left(\left(h_k(y) + \sqrt{\mu_k} Y_k\right)^2 - \mu_k - h_k(y)^2\right)^4 \right] \xrightarrow{n \to \infty} 0, \end{align*} where $s_n^2 = \sum_{k=1}^{u_n-\lvert U\rvert - 1} (2 \mu_k^2 + 4 h_k(y)^2 \mu_k) = 2(1+o(1)) \sum_{k=1}^{u_n-\lvert U\rvert - 1} \mu_k^2$ by~\eqref{eq:bound_shift}. We observe that this convergence indeed takes place since \begin{multline*} \sum_{k=1}^{u_n - \lvert U\rvert - 1}\E\left[\left(\left(h_k(y) + \sqrt{\mu_k} Y_k\right)^2 - \mu_k - h_k(y)^2\right)^4 \right] \\ \begin{aligned} &\leq c \sum_{k=1}^{u_n - \lvert U\rvert - 1} \left(\mu_k^4 + \mu_k^2 h_k(y)^4 + \mu_k^3 h_k(y)^2\right) \\ &\leq c \left(\sum_{k=1}^{u_n - \lvert U\rvert - 1} \mu_k^4 + \mu_1^2 \ntwo{h(y)}^4 + \mu_1^3 \ntwo{h}^2\right), \end{aligned} \end{multline*} and by~\eqref{eq:EV_asympt}, \eqref{eq:auxiliary_estim}, and~\eqref{eq:bound_shift}, for $d \geq 4$, this is of order $\cO (n^8\ind_{\{d \leq 7\}} + n^d \log n \ind_{\{d = 8\}} + n^d \ind_{\{d \geq 9\}})$, while $s_n^4 = \Theta(n^8 (\log n)^2 \ind_{\{d = 4\}} + n^{2d} \ind_{\{d \geq 5\}})$. For the fourth criterion, consider the set $J = \{\sigma_k^2: k \leq u_n - \lvert U\rvert -1$ such that $\lvert h_k(y)\rvert > 1/\log n\}$ (there are at most $\log n \norm{h(y)}^2$ such indices $k$). Notice that since $(\mu_k)_{k \leq u_n - \lvert U\rvert - 1}$ are uniformly bounded away from zero and $\sigma_{k}^2 = 2 \mu_k^2 + 4\mu_k h_k(y)^2$, the order of the variances on the complement of $J$ will be determined by the order of $\mu_k$'s. Furthermore, elements of $J$ after reordering can only be sent forward or remain at their place. Hence, for any $L \geq 1$, \begin{align*} \left\{\sigma_k^2: k \geq L\right\} \setminus J \subset \left\{\sigma_{\tilde k}^2: \tilde k \geq L\right\} \setminus J \subset \left\{\sigma_k^2: k \geq L - 2\lvert J\rvert\right\} \setminus J. \end{align*} Let $\varepsilon \in (0,1)$ be fixed but arbitrarily small. We choose the constants as follows: \begin{itemize} \item for $d \geq 5$, let $l^* = l_* \coloneqq \inf\{l \geq \varepsilon n^d: l - 2\lvert J\rvert \notin J, l \in \Z\}$, $K = K(\varepsilon)$ a large constant depending only on $\varepsilon$; \item for $d=4$, let $l^*\coloneqq \inf\{l \geq \varepsilon n^d: l - 2\lvert J\rvert \notin J, l \in \Z\}$, $l_*\coloneqq \inf\{l \geq n^{d-\varepsilon}: l - 2\lvert J\rvert \notin J, l \in \Z\}$, $K(n)= n^\varepsilon$. \end{itemize} Note that $l^*, l_*$ can at most be equal to the reference value (either $\varepsilon n^d$ or $n^{d-\varepsilon}$) plus $\lvert J\rvert = o(n^\delta)$ for any $\delta > 0$. Moreover, for $L \in \{l_*, l^*\}$, \begin{align*} \left\{\sigma_k^2: k \geq L\right\} \setminus J \subset \left\{\sigma_{\tilde k}^2: \tilde k \geq L\right\} \subset \left\{\sigma_k^2: k \geq L - 2\lvert J\rvert\right\}. \end{align*} This implies that the $l^{*\rm th}$ element of the reordered sequence $(\sigma^2_{\tilde k})_{\tilde k}$ is smaller or equal to $\sigma^2_{l^* - 2\lvert J\rvert}$ (of the original order). We first verify the assumptions of the fourth criterion for $d\geq 5$. By~\eqref{eq:EV_asympt} and~\eqref{eq:auxiliary_estim}, for $u_n-\lvert U\rvert - 1 \geq k \geq l_* - 2\lvert J\rvert, k \notin J$, $0 < c_* \leq \mu_k = \Theta(n^2/k^{2/d}) = \Theta(1) \leq \mu_{l_* - 2\lvert J\rvert} \leq c^* < \infty$, $\sigma_{k}^2 = 2\mu_k^2 + 4 \mu_k h_k(y)^2 = \Theta(1)$ uniformly. Thus, for $Z_{k} = (h_k(y) + \sqrt{\mu_k} Y_k)^2 - \mu_k - h_k(y)^2$, $Z \sim \cN(0,1)$, \begingroup \allowdisplaybreaks \begin{align*} \frac{\sum_{\tilde k \geq l_*}^{u_n-\lvert U\rvert - 1} \E\left[Z_{\tilde k}^2 \ind_{\{\lvert Z_{\tilde k}\rvert > K\}}\right]}{\sum_{\tilde k \geq l_*}^{u_n-\lvert U\rvert - 1} \sigma_{\tilde k}^2} &\leq \frac{\sum_{k \geq l_* - 2\lvert J\rvert, k \notin J}^{u_n-\lvert U\rvert - 1} \E\left[Z_k^2 \ind_{\{\lvert Z_k\rvert > K\}}\right] + \lvert J\rvert \sigma^2_{l_* - 2\lvert J\rvert}}{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \sigma^2_{k}}\\ &\leq \frac{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \E\left[Z_k^2 \ind_{\{\lvert Z_k\rvert > K\}}\right] + 3\lvert J\rvert \sigma^2_{l_* - 2\lvert J\rvert}}{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \sigma^2_{k}}\\ &\leq \frac{\E\left[(Z^2 - 1)^2 \ind_{\left\{\mu_{l_* - 2\lvert J\rvert} Z^2 + \frac{2 \sqrt{\mu_{l_* - 2\lvert J\rvert}}}{\log n} \lvert Z\rvert > K \right\}}\right] \sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \mu_k^2}{2 \sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \mu_k^2} \\ &\qquad+\frac{4\mu_1 \ntwo{h(y)}^2 + \frac{1}{\log n} 4 \sqrt{3} \sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \mu_k^{3/2} + o(n)}{2 \sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \mu_k^2}\\ &\leq \frac{1}{2} \E\left[(Z^2 - 1)^2 \ind_{\{ Z^2 + \lvert Z\rvert > K/(2 c^*)\}}\right] + o(n^{3-d}) + \Theta \left(\frac{1}{\log n}\right) \end{align*} \endgroup can be made arbitrarily small by choosing $K > c^*$ sufficiently large. Furthermore, \begin{align*} \frac{\sum_{\tilde k \geq l_*}^{u_n-\lvert U\rvert - 1} \sigma_{\tilde k}^2}{s_n^2} \geq \frac{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \sigma_{k}^2}{s_n^2} &\geq \frac{\sum_{k = \varepsilon n^d + n}^{u_n-\lvert U\rvert - 1} 2 \mu_k^2}{\Theta (n^d)} \\ &\geq \frac{\Theta \left(n^4 \sum_{k = \varepsilon n^d + n^{d/2}}^{u_n-\lvert U\rvert - 1} k^{-\frac{4}{d}} \right)}{\Theta (n^d)} = \left(1 - \varepsilon^{\frac{d-4}{d}}\right) \Theta(1); \\ \frac{u_n-\lvert U\rvert - 1 - l^*}{\sigma^2_{l^*, \text{reorder.}} \vee K^2} &\geq \frac{u_n-\lvert U\rvert - 1 - l^*}{\sigma^2_{l^* - 2\lvert J\rvert} \vee K^2} = \Theta(n^d) \gg \log(n^d). \end{align*} This completes the verification of all required conditions for $d \geq 5$. For $d= 4$, the argument proceeds analogously. Using~\eqref{eq:EV_asympt} and~\eqref{eq:auxiliary_estim}, we obtain: \begingroup \allowdisplaybreaks \begin{align*} \frac{\sum_{\tilde k \geq l_*}^{u_n-\lvert U\rvert - 1} \sigma_{\tilde k}^2}{s_n^2} \geq \frac{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \sigma_{k}^2}{s_n^2} &\geq \frac{\Theta\left(\sum_{k = 2n^{d-\varepsilon}}^{u_n-\lvert U\rvert - 1} n^4/k \right)}{\Theta (n^4 \log n)} = \frac{\Theta \left(\log \left(\frac{u_n}{n^{d-\varepsilon}}\right)\right)}{\Theta (\log n)} = \varepsilon \Theta(1);\\ \frac{u_n-\lvert U\rvert - 1 - l^*}{\sigma^2_{l^*, \text{reorder.}} \vee K(n)^2} &\geq \frac{u_n-\lvert U\rvert - 1 - l^*}{\sigma^2_{l^* - 2\lvert J\rvert} \vee K(n)^2} = \Theta \left(\frac{n^d}{(n^4/(\varepsilon n^d)^{4/d}) \vee n^{2\varepsilon}}\right)\\ &= \Theta(n^{d-2 \varepsilon}) \gg \log(n^4 \log n); \\ \frac{\sum_{\tilde k \geq l_*}^{u_n-\lvert U\rvert - 1} \E\left[Z_{\tilde k}^2 \ind_{\{\lvert Z_{\tilde k}\rvert > K(n)\}}\right]}{\sum_{\tilde k \geq l_*}^{u_n-\lvert U\rvert - 1} \sigma_{\tilde k}^2} &\leq \frac{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \E\left[Z_k^2 \ind_{\{\lvert Z_k\rvert > K(n)\}}\right] + 3\lvert J\rvert \sigma^2_{l_*-2\lvert J\rvert}}{\sum_{k \geq l_*, k \notin J}^{u_n-\lvert U\rvert - 1} \sigma_{k}^2} \\ &\leq \frac{1}{2} \E\left[(Z^2 - 1)^2 \ind_{\left\{ Z^2 + \lvert Z\rvert > \frac{K(n)}{2 \mu_{l_* - 2\lvert J\rvert}}\right\}}\right] \\ &\qquad+\frac{\Theta \left(n^2 \ntwo{h(y)}^2 + \frac{\sum_{k \geq l_*}^{u_n} \mu_k^{3/2}}{\log n} + \lvert J\rvert \sigma^2_{l_* - 2\lvert J\rvert}\right)}{\varepsilon \Theta (n^4 \log n)}\\ &= \frac{1}{2} \E\left[(Z^2 - 1)^2 \ind_{\left\{ Z^2 + \lvert Z\rvert > \frac{n^\varepsilon}{\Theta(n^{\varepsilon/2})}\right\}}\right] + \Theta \left(\frac{1}{\log^2 n}\right).\qedhere \end{align*} \endgroup The latter is arbitrarily small. \end{proof} \subsubsection*{Proof of Claims~\ref{claimA} and~\ref{claimB} for \texorpdfstring{$d=3$}{d=3}} \begin{proof}[Proof of Claim~\ref{claimB}] The statement can be verified fully analogously to Claim~\ref{claim3} in the low-temperature regime (see Section~\ref{subsec:sph_Tlow}) using eigenvalue representation~\eqref{eq:EV_representation} once we notice that for $\mu \in \{\mu_1, \mu_2, \eta_2^{-1}, \eta_3^{-1}\}$ and $v \in \{\Var[\norm{\hat{\gamma}^y}^2], \Var[\norm{{\gamma}}^2]\}$: $\mu/\sqrt{v} = \Theta(1)$ uniformly in $n$. The latter is due to the above observations that $\mu, \eta^{-1} = \Theta(n^2)$ and the variances are of order $\Theta(n^4)$. \end{proof} \begin{proof}[Proof of Claim~\ref{claimA}] Recall that if a sequence of random variables $(X_n)_n$ possesses moment generating functions (MGF) $(M_n)$ that are finite on a common open interval around zero, $(-a,a)$ for $a>0$, i.e., $M_n(u) < \infty$ for all $n \in \N, u \in (-a,a)$, and there exists $0 0$ as $n\to \infty$. Thus, a suitable $b > 0$ can be chosen as required. Next, we find explicit expressions for the MGFs of $\hat{X}_n$ and $X_n$. To this end, let $t \in [-b,b]$ and recall that the MGF of $(Z+\mu)^2$ for a standard normal variable $Z$ is given by $\exp (\frac{\mu^2 x}{1-2x} - \frac{1}{2} \log(1-2x))$ for all $x< 1/2$. Combined with the eigenvalue representation~\eqref{eq:EV_representation} of $\hat{X}_n$ and the fact that $\Var[\norm{\hat{\gamma}^y}^2] = 2(1+o(n^{\varepsilon - 1})) \sum_{k=1}^{u_n-1} \mu_k^2$, this yields \begingroup \allowdisplaybreaks \begin{align*} \hat{M}_n(t) &= \exp \left(-t \sum_{k=1}^{u_n-1}\frac{\mu_k}{\sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]}} - \frac{1}{2} \sum_{k=1}^{u_n-1} \log\left(1 - \frac{2t\mu_k}{\sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]}}\right) \right) \\ &\qquad\times \exp \Biggl(\underbrace{t \sum_{k=1}^{u_n-1}\frac{h_k(y)^2}{\sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]}} \biggl(\underbrace{ \biggl(1 - \frac{2t \mu_k}{\sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]}} \biggr)^{-1} - 1}_{\lvert \cdot\rvert \in [0,3)}\biggr)}_{\lvert \cdot\rvert \in \bigl[0, 3b \ntwo{h(y)}^2 / \sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]} \bigr) \overset{\eqref{eq:bound_shift}}{\subset} (0, o(n^{\varepsilon - 1}))} \Biggr) \\ &= (1 + o(n^{\varepsilon - 1})) \exp \Bigg(\sum_{l=2}^\infty \underbrace{\frac{1}{2l} \biggl(\frac{\sqrt{2} t}{1+o(n^{\varepsilon - 1})}\biggr)^l \sum_{k=1}^{u_n-1} \biggl(\frac{\mu_k^2}{ \sum_{p=1}^{u_n-1} \mu_p^2}\biggr)^{l/2}}_{\eqqcolon f_{n}(l)} \Biggr) \end{align*} \endgroup for any $\varepsilon \in (0,1)$. Here, we additionally used that \[ 0 < \frac{2\lvert t\rvert\mu_k}{\sqrt{\Var[\ntwo{\hat{\gamma}^y}^2]}} \leq \frac{(1+o(1)) 2 b\mu_1}{\sqrt{\Var[\ntwo{{\gamma}}^2]}} < \frac{1+o(1)}{2} \eta_2 \mu_1 \leq \frac{1+o(1)}{2} < \frac{3}{4} \] for all $n$ sufficiently large, which follows by~\eqref{eq:auxiliary_estim} and~\eqref{eq:EV_asympt}. In order to show that the limits of the MGFs $\hat{M}_n$ and $M_n$ exist and coincide, we first observe that the sum in $f_n(l)$ equals $1$ when $l=2$. Since moreover $\mu_k^2/ \sum_{p=1}^{u_n-1} \mu_p^2 < 1$ for all $k$, it remains bounded by $1$ for all $l > 2$. Therefore, uniformly in $t \in [-b,b]$, \[ \lvert f_n(l)\rvert \leq \frac{1}{2l} \left(\frac{3}{4}\right)^l. \] Since the latter bound is summable, by dominated convergence theorem, we obtain \begin{align*} \lim_{n\to\infty} \hat{M}_n(t) = \exp \left(\sum_{l=2}^\infty \frac{1}{2l} (\sqrt{2} t)^l \lim_{n \to \infty} \sum_{k=1}^{u_n-1} \left(\frac{\mu_k^2}{ \sum_{p=1}^{u_n-1} \mu_p^2}\right)^{l/2} \right). \end{align*} And analogously, \begin{align*} \lim_{n\to\infty} M_n(t) = \exp \left(\sum_{l=2}^\infty \frac{1}{2l} (\sqrt{2} t)^l \lim_{n \to \infty} \sum_{k=2}^{n^d} \left(\frac{1/\eta_k^2}{ \sum_{p=2}^{n^d} 1/\eta_p^2}\right)^{l/2} \right). \end{align*} It thus remains to show that for each fixed $l \geq 3$, \[ \lim_{n} \sum_k \left(\frac{\mu_k^2}{ \sum_p \mu_p^2}\right)^{l/2} = \lim_{n} \sum_k \left(\frac{1/\eta_k^2}{ \sum_p 1/\eta_p^2}\right)^{l/2}, \] and that these limits are indeed well-defined. But before turning to proving this, observe that for $l=2$, both sides of the latter expression are equal to $1$, which then suffices to conclude that the limiting distribution of $(\hat{X}_n)$ and $(X_n)_n$ is non-trivial. We start by verifying the well-definedness of \[ \lim_{n} \sum_k \left(\frac{1/\eta_k^2}{ \sum_p 1/\eta_p^2}\right)^{l/2}. \] First note that \begin{align*} \sum_{k=2}^{n^d} \left(\frac{1/\eta_k^2}{ \sum_{p=2}^{n^d} 1/\eta_p^2}\right)^{l/2} &=\sum_{x \in [0,n) \cap \Z^d \setminus \{0\}}\hspace{-0.2cm} \left(\frac{1/\eta_x^2}{ \sum_{y \in [0,n) \cap \Z^d \setminus \{0\}} 1/\eta_y^2}\right)^{l/2} \\ &= \sum_{x \in \hat{\Lambda}_n \setminus \{0\}} \left(\frac{1/\eta_x^2}{\sum_{y \in \hat{\Lambda}_n \setminus \{0\}} 1/\eta_y^2}\right)^{l/2}, \end{align*} where $\eta_x = 2 \sum_{i=1}^3 (1-\cos(2\pi x_i/n))$ and $\hat{\Lambda}_n = [-n/2, n/2)^3 \cap \Z^3$. Let $f(n)$ be an arbitrary function slowly growing to infinity such that $f(n) = o(n)$, e.g., $f(n) = \log(n)$. Then, using explicit form of $\eta_x$ and Taylor expanding it, we easily get that \begin{align*} \sum_{y \in \hat{\Lambda}_n \setminus \{0\}} \frac{1}{n^4 \eta_y^2} = \bigl(1 + \cO\bigl(f(n)^{-1}\bigr)\bigr) \sum_{\substack{y \in \hat{\Lambda}_n \setminus \{0\}\\ \lvert y\rvert \leq f(n)}} \frac{1}{(2\pi)^4 \lvert y\rvert^4} + o(1) \overset{n \to \infty}{\sim} \sum_{y \in \Z^3\setminus \{0\}} \frac{1}{(2\pi)^4 \lvert y\rvert^4}. \end{align*} Fully analogously, \begin{align*} \sum_{x \in \hat{\Lambda}_n \setminus \{0\}} \left(\frac{1}{n^2 \eta_x}\right)^l \overset{n \to \infty}{\sim} \sum_{x \in \Z^3\setminus \{0\}} \frac{1}{(2\pi)^{2l} \lvert x\rvert^{2l}}, \end{align*} which then yields \begin{align*} \lim_{n\to \infty} \sum_k \left(\frac{1/\eta_k^2}{ \sum_p 1/\eta_p^2}\right)^{l/2} = \sum_{x \in \Z^3\setminus \{0\}} \frac{1}{\lvert x\rvert^{2l}}\,\Bigg/\,\left(\sum_{y \in \Z^3\setminus \{0\}} \frac{1}{\lvert y\rvert^4} \right)^{l/2} < \infty. \end{align*} To conclude the proof of the claim, it only remains to show that the ratio \[ \sum_k \left(\frac{\mu_k^2}{ \sum_p \mu_p^2}\right)^{l/2}\,\Bigg/\,\sum_k \left(\frac{1/\eta_k^2}{ \sum_p 1/\eta_p^2}\right)^{l/2} \xrightarrow{n\to \infty} 1. \] The latter can further be reduced to showing $T_n \coloneqq \sum_k \mu_k^l/\sum_k 1/\eta_k^l \nearrow 1$ as $n$ tends to infinity, which is due to the fact that $\sum_p \mu_p^2 = \frac{1}{2} \Var[\norm{\gamma}^2] + o(1) = \sum_p 1/\eta_p^2 +o(1)$ and $\sum_k \mu_k^l/\sum_k 1/\eta_k^l \leq 1$. The former observation follows from the proof of~\eqref{eq:var_error}, while the latter is a consequence of~\eqref{eq:EV_asympt} and~\eqref{eq:auxiliary_estim} (since the contribution of the smallest $\lvert U\rvert$ eigenvalues $\mu_k$ to the sum is irrelevant). Our objective now is to find an alternative expression for $T_n$ in terms of traces of the zero-average Green's function and its appropriate principal submatrix. To this end, note that for $k \ll n^d$, by the proof of~\eqref{eq:auxiliary_estim}, we get that $\mu_k \geq \max(\zeta_k -\upsilon_1, \zeta_{k+\lvert U\rvert})$, where $\zeta_k$ is the $k^{\rm th}$ largest eigenvalue of $G' \coloneqq G^{0\avg}_{\Lambda_n}\vert_{U^c\times U^c}$ and $\upsilon_1$ is the largest eigenvalue of $A \coloneqq G^{0\avg}_{\Lambda_n}\vert_{U^c\times U}\,G^{0\avg}_{\Lambda_n}\vert_{U\times U}^{-1} \,G^{0\avg}_{\Lambda_n}\vert_{U\times U^c}$. We now show that for $k \leq \log^3 n$: $\mu_k = \zeta_k(1-o(1))$ with $o(1)$ uniform in this range of $k$. First observe that by~\eqref{eq:EV_asympt}, $\zeta_k \geq 1/\eta_{k+\lvert U\rvert + 1} = \Theta(n^2/k^{2/3})$. A sufficiently good upper bound for $\upsilon_1 > 0$ follows from~\eqref{eq:0avgGfct_polydec}: \[ \upsilon_1 \leq \Trace(A) \leq c(U) \sum_{x \in U^c, y \in U} G^{0\avg}_{\Lambda_n}(x,y)^2 = \cO(n). \] Analogously to the above, this implies that \begin{align*} \frac{\sum_{k=1}^{u_n-1} (\mu_k/n^2)^l}{\sum_{k=2}^{n^d} 1/(n^2\eta_k)^l} {\sim} \frac{\sum_{k=1}^{\log^3 n} (\mu_k/n^2)^l}{\sum_{k = 2}^{\log^3 n} 1/(n^2\eta_k)^l} &\sim \frac{\sum_{k=1}^{\log^3 n} (\zeta_k/n^2)^l}{\sum_{k = 2}^{\log^3 n} 1/(n^2\eta_k)^l} \\ &\sim \frac{\sum_{k=1}^{u_n} (\zeta_k/n^2)^l}{\sum_{k=2}^{n^d} 1/(n^2\eta_k)^l} = \frac{\Trace((G')^l)}{\Trace((G^{0\avg}_{\Lambda_n})^l)} \end{align*} as $n \to \infty$, and that $\Trace((G')^l) = \Theta(n^{2l})$. To show that the latter ratio converges to $1$, observe first that~\eqref{eq:0avgGfct_polydec} yields that there is an absolute constant $K>0$ such that for all $x,y \in \hat{\Lambda}_n$, $|G^{0\avg}_{\Lambda_n}(x,y)| \leq K G_{\Z^d}(x,y)$. With the convention that $x_{l+1} = x_1$, we have \pagebreak \begin{multline*} \left|\Trace\left(\left(G^{0\avg}_{\Lambda_n}\right)^l\right) - \Trace((G')^l)\right|\\ \begin{aligned} &\leq \sum_{k=1}^l \sum_{\substack{I \subset \{1, \dots, l\}\\ \lvert I\rvert = k}} \sum_{\substack{x_i \in U: i \in I;\\ x_j \in U^c: j \notin I}} \prod_{i=1}^l \left|G^{0\avg}_{\Lambda_n} (x_i, x_{i+1})\right|\\ &\leq \sum_{k=1}^l K^l \sum_{\substack{I \subset \{1, \dots, l\}\\ \lvert I\rvert = k}} \sum_{\substack{x_i \in U: i \in I;\\ x_j \in U^c: j \notin I}} \prod_{i=1}^l G_{\Z^d} (x_i, x_{i+1})\\ &\leq \sum_{k=1}^l K^l \binom{l}{k} \sum_{u \in U} \sum_{\substack{x_i \in \hat{\Lambda}_n: i \neq 1}} G_{\Z^d}(u, x_2) \left(\prod_{i=2}^{l-1} G_{\Z^d} (x_i, x_{i+1})\right) G_{\Z^d}(x_l, u)\\ &\leq (2^l - 1) K^l \sum_{u \in U} \lambda_{\max} \left(G_{\Z^d}\vert_{\hat{\Lambda}_n \times \hat{\Lambda}_n}^{l-2}\right) \sum_{\substack{x\in \hat{\Lambda}_n}} G_{\Z^d}(u, x)^2\\ &\leq c n \lambda_{\max} \left(G_{\Z^d}\vert_{\hat{\Lambda}_n \times \hat{\Lambda}_n}\right)^{l-2} \end{aligned} \end{multline*} for an appropriate $c = c(l, U)>0$. Furthermore, by Cauchy--Schwarz inequality, \begin{align*} \lambda_{\max} \left(G_{\Z^d}\vert_{\hat{\Lambda}_n \times \hat{\Lambda}_n}\right) \leq \sup_{h \in \R^{\hat{\Lambda}_n}: \ntwo{h} = 1} \sum_{x, y \in \hat{\Lambda}_n} h_x G_{\Z^d}(x,y) h_y \leq \sup_{x \in \hat{\Lambda}_n} \sum_{y \in \hat{\Lambda}_n} G_{\Z^d}(x,y) \leq c' n^2. \end{align*} Thus, the ratio of interest $\Trace((G')^l)/ \Trace((G^{0\avg}_{\Lambda_n})^l) = 1 + \cO(n^{-3})$ converges to $1$ as $n$ tends to infinity. \end{proof} \begin{rema}\label{rem:Malliavin} As mentioned above, instead of relying on Proposition~\ref{prop:dens_conv_1D} and the estimate~\eqref{eq:0avgGfct_polydec}, which together yield uniform convergence of the probability density functions of the variables $X_n$ and $\hat{X}_n$ to the standard normal density for $d \geq 4$, one could consider an alternative route via Malliavin calculus~\cite{HU2014814}. More precisely, by slightly modifying the proof of Lemma~7.1 in that reference, one can verify the assumptions of their Corollary~6.6 (cf. Theorem~7.3) for $\hat{X}_n$ and $X_n$ in dimensions $d \geq 5$. For $d = 4$, however, checking the fifth assumption would require an additional argument to replace Lemma~7.1, which we expect would resemble the proof of our Proposition~\ref{prop:dens_conv_1D}. Finally, this alternative approach does not extend to $d = 3$, as the variables $\hat{X}_n$ and $X_n$ are not Gaussian in the limit. \end{rema} \subsection{Local correlation functions of the Spherical model}\label{subsec:spher_correlations} The key objective of this section is to prove Corollary~\ref{cor:spher_correlations}, namely to show that \begin{equation}\label{eq:cov_converg} \nu_{\Lambda_n, \beta} \left[\prod_{x \in U} \theta_x^{i_x}\right] \xrightarrow{n \to \infty} \E\left[\prod_{x \in U} \alpha_x^{i_x}\right], \end{equation} where $\sqrt{\beta} \alpha$, as described in Theorem~\ref{thm:spher_detailed}, is either an $m^2$-massive GFF on $\Z^d$ if $\beta < \beta_c$, a GFF on $\Z^d$ if $\beta = \beta_c$ or a GFF on $\Z^d$ plus an independent drift $\sqrt{\beta - \beta_c} X \bfun_{\Z^d}$ with a Rademacher random variable $X$ if $\beta > \beta_c$. Since $\theta_U$ already converges in law to $\alpha_U$, it suffices to verify the uniform integrability of the variables $\prod_{x \in U} \theta_x^{i_x}$ in order to conclude~\eqref{eq:cov_converg} (see, e.g., \cite[Theorem~25.12]{billingsley1995probability}). The latter follows instantaneously from the following: \begin{lemm}[Moments of spins of the spherical model]\label{lemma:spher_moments} Let $\theta$ be a configuration of the spherical model on $\Lambda_n$ at inverse temperature $\beta > 0$. For any $p \in \N_0$, there exists a constant $C_p > 0$, independent of $n$, such that for all $n$, \begin{align*} \nu_{\Lambda_n, \beta}\left[\lvert\theta_0\rvert^{2p}\right] \leq C_p. \end{align*} \end{lemm} Indeed, since for a sequence of random variables $(X_n)$ to be uniformly integrable, it suffices to satisfy the condition $\sup_n \E[\lvert X_n\rvert^2] < \infty$; and by generalized H\"older's inequality, we have \begin{align*} \nu_{\Lambda_n, \beta} \left[\prod_{x \in U} \lvert\theta_x\rvert^{2i_x}\right] \leq \prod_{x \in U} \nu_{\Lambda_n, \beta}\left[\lvert\theta_x\rvert^{2i_x \lvert U\rvert}\right]^{1/\lvert U\rvert} \leq \nu_{\Lambda_n, \beta}\left[\lvert\theta_0\rvert^{2 \lvert U\rvert \max_x(i_x)}\right] \leq C_{|U| \max_x (i_x)}. \end{align*} The second to last inequality follows from the fact that spins of the spherical model on torus are identically distributed. We now proceed to the proof of Lemma~\ref{lemma:spher_moments}, which is motivated by~\cite[Lemma~3]{Lukkarinen2020}. \begin{proof} Recall that by Proposition~\ref{prop:rel_GFF_spher}, $\sqrt{\beta} \theta$ has the same law as any massive ($m^2>0$) GFF $\phi$ on $\Lambda_n$ conditioned on $\norm{\phi}^2 = \beta n^d$. Therefore, \begin{align*} \beta^p \nu_{\Lambda_n, \beta}\left[\theta_0^{2p}\right] = \E\left[\phi_0^{2p}\,\middle|\,\ntwo{\phi}^2 = \beta n^d\right]. \end{align*} Let $Q = (q_{x}^w)_{x, w \in [0, n)^d \cap \Z^d}$ be the orthonormal matrix diagonalizing $G_{\Lambda_n, m^2}$, whose columns are the eigenvectors $(q^w)_w$; in particular, we have $Q^T G_{\Lambda_n, m^2} Q = \diag(((m^2+\eta_w)^{-1})_{w})$, where $(\eta_w)_{w \in [0, n)^d \cap \Z^d}$ are the eigenvalues of $-\Delta_{\Lambda_n}$. Let $(Z_w)_w$ be i.i.d. standard normal random variables. Then, \begin{align*} \phi = (\phi_x)_{x \in \Lambda_n} &\overset{\law}{=} \left(\sum_{w \in [0, n)^d \cap \Z^d} q_{x}^w \frac{1}{\sqrt{m^2+\eta_w}} Z_w\right)_{x \in \Lambda_n}; \\ \ntwo{\phi}^2 &\overset{\law}{=} \sum_{w \in [0, n)^d \cap \Z^d} \frac{1}{\eta_w + m^2} Z_w^2. \end{align*} Using this, we obtain on the one hand that \begin{align*} \beta^{p} n^{dp} = \E\left[\ntwo{\phi}^{2p}\,\middle|\,\ntwo{\phi}^2 = \beta n^d\right] = \sum_{\substack{w_1, \dots, w_{p}\\ \in[0, n)^d \cap \Z^d }} \E\left[\prod_{j=1}^{p} \frac{Z_{w_j}^2}{m^2+ \eta_{w_j}}\,\middle|\,\sum_{w} \frac{Z_w^2}{\eta_w + m^2} = \beta n^d\right]. \end{align*} And on the other hand, by Lemma~\ref{lemma:massive_EVs}, \begin{align*} \E\left[\phi_0^{2p}\,\middle|\,\ntwo{\phi}^2 = \beta n^d\right] &= \sum_{\substack{w_1, \dots, w_{2p}\\ \in[0, n)^d \cap \Z^d }} \left(\prod_{j=1}^{2p} q_{0}^{w_j}\right) \E\left[\prod_{j=1}^{2p} \frac{Z_{w_j}}{\sqrt{m^2+ \eta_{w_j}}}\,\middle|\,\sum_{w} \frac{Z_w^2}{\eta_w + m^2} = \beta n^d\right]\\ &= \frac{1}{n^{dp}} \sum_{\substack{w_1, \dots, w_{2p}\\ \in[0, n)^d \cap \Z^d }}\E\left[\prod_{j=1}^{2p} \frac{Z_{w_j}}{\sqrt{m^2+ \eta_{w_j}}}\,\middle|\,\sum_{w} \frac{Z_w^2}{\eta_w + m^2} = \beta n^d\right]. \end{align*} Since the conditional law of $(Z_{w}/\sqrt{m^2+ \eta_{w}})_w$, given $\sum_{w \in [0, n)^d \cap \Z^d} \frac{1}{\eta_w + m^2} Z_w^2 = \beta n^d$, is symmetric in each coordinate, it is invariant under sign flips. Therefore, any product involving an odd power of some $Z_w$ vanishes under the conditional expectation, and hence, \begin{align*} \E\left[\phi_0^{2p}\,\middle|\,\ntwo{\phi}^2 = \beta n^d\right] &\leq \frac{(2p)!}{n^{dp}} \sum_{\substack{w_1, \dots, w_{p}\\ \in[0, n)^d \cap \Z^d }}\E\left[\prod_{j=1}^{p} \frac{Z^2_{w_j}}{m^2+ \eta_{w_j}}\,\middle|\,\sum_{w} \frac{Z_w^2}{\eta_w + m^2} = \beta n^d\right] \\ &= (2p)! \beta^p. \end{align*} This concludes the proof, establishing the uniform moment bound $\nu_{\Lambda_n, \beta}[\theta_0^{2p}] \leq (2p)! = C_p$. \end{proof} \section{The infinite spin-dimensionality distributional limit of the spin \texorpdfstring{$O(N)$}{O(N)} model}\label{sec:O(N)} This section studies the ``local'' infinite spin-dimensionality limit ($N \to \infty$) of the spin $O(N)$ model both on the finite domain (the discrete $d$-dimensional torus $\Lambda_n$ for $d \geq 2$) and in the subsequent infinite-volume limit ($n \to \infty$). More precisely, we prove the following two theorems corresponding to these respective cases that combined with Scheffé's lemma yield Theorem~\ref{thm:spinO(N)} in the introduction. \begin{theo}[$N\to \infty$ limit of the spin $O(N)$ model on a finite torus]\label{thm:spinO(N)_torus} Let $\Lambda_n = \T^d_n$ and $S = (S_x)_{x \in \Lambda_n}$ be a configuration of the spin $O(N)$ model at inverse temperature $\beta > 0$. For any $M \in \N$, the density function of $\ovlineM{S} \coloneqq (S^i_x)_{x \in \Lambda_n}^{i = 1, \dots, M}$ converges as $N \to \infty$ to that of an $M$-vectorial massive GFF on $\Lambda_n$ scaled by $1/\sqrt{\beta}$ with the mass $m_n^2$ depending on $\beta, d$, and $n$ in such a way that \begin{align}\label{eq:choice_of_mass} G_{\Lambda_n, m_n^2}(x,x) = \beta \quad \text{for all } x\in \Lambda_n. \end{align} or more explicitly\footnote{This follows from the fact that the trace of a matrix equals the sum of its eigenvalues and Lemma~\ref{lemma:massive_EVs}.}, $m^2$ is the unique solution to \begin{align}\label{eq:mass_spinO(N)} \frac{1}{m_n^2 n^d} + \frac{1}{n^d} \sum_{x \in [0,n)^d \cap \Z^d \setminus\{0\}} \frac{1}{m_n^2 + 2\sum_{k=1}^d (1 - \cos(2\pi x_k/n))} = \beta. \end{align} \end{theo} \begin{theo}[Infinite-volume $N\to \infty$ limit of the spin $O(N)$ model]\label{thm:spinO(N)_infVol} For each $n, N \in \N$, let $S = (S_x)_{x \in \Lambda_n}$ be a configuration of the spin $O(N)$ model on $\Lambda_n$ at inverse temperature $\beta > 0$. Set $\beta_c \coloneqq G_{\Z^d}(0,0)$. Then, for any $M \in \N$ and a finite set $U \subset \Z^d$ (considered as a subset of $\Lambda_n$ for each $n$), $\ovlineM{S}_{U} \coloneqq (S^i_x)_{x \in U}^{i = 1, \dots, M}$ converges in law as $N \to \infty$ followed by $n \to \infty$ towards: \begin{enumerate} \item\label{theo4.2.1} $\beta < \beta_c$: an $M$-vectorial massive GFF on $\Z^d$ restricted to $U$ scaled by $1/\sqrt{\beta}$ with the mass $m^2$ depending on $\beta$ and $d$ in such a way that \begin{align*} G_{\Z^d, m^2}(x,x) = \beta \quad \text{for all } x\in \Z^d; \end{align*} \item\label{theo4.2.2} $\beta = \beta_c$: an $M$-vectorial GFF on $\Z^d$ restricted to $U$ scaled by $1/\sqrt{\beta}$; \item\label{theo4.2.3} $\beta > \beta_c$: an $M$-vectorial GFF on $\Z^d$ restricted to $U$ scaled by $1/\sqrt{\beta}$ plus an independent constant random drift of the form $\sqrt{\frac{\beta - \beta_c}{\beta}} \ovlineM{Z} \bfun_U$ with $\ovlineM{Z}$ being an $M$-dimensional standard normal vector. \end{enumerate} \end{theo} \subsection{The infinite spin-dimensionality limit on the torus}\label{subsec:spinO(N)_torus} This subsection is concerned with the infinite spin-dimensionality limit of the spin $O(N)$ model on the discrete torus, and in particular proves Theorem~\ref{thm:spinO(N)_torus}. Since $n \in \N$ is fixed throughout this section, we write $\Lambda \coloneqq \T^d_n$ for brevity. Before proceeding to the proof, let us quickly discuss the strategy. First, by Proposition~\ref{prop:mGFF-O(N)}, we know that for any $m^2$-massive $N$-vectorial GFF $\phi$ with $m^2>0$, \begin{align*} \Law \bigl(\ovlineM{S}\bigr) = \Law \left(\frac{\ovlineM{\phi}}{\sqrt{\beta}}\,\middle|\,|\phi_x| = \sqrt{\beta N}, \,\forall x \in \Lambda \right), \end{align*} where $\ovlineM{\phi}$ denotes the first $M$ components and $\rest{\phi}$ the remaining $N-M$ components of $\phi$. Thus, for any $\xi \coloneqq (\xi_x)_{x\in \Lambda} \in (\R^M)^{|\Lambda|}$, we have the following equality of densities with respect to Lebesgue measure on $\R^{M |\Lambda|}$: \begin{equation}\label{eq:dens_O(N)} f_{\ovlineM{S}}(\xi) = f_{\frac{\ovlineM{\phi}}{\sqrt{\beta}}\big\vert (|\phi_x|)_x = \left(\sqrt{\beta N}\right)_x}(\xi) = f_{\frac{\ovlineM{\phi}}{\sqrt{\beta}}} \left(\xi \right) \frac{f_{(|\rest{\phi_x}|^2)_x} (({\beta N} - \beta |\xi_x|^2)_x)}{f_{(|\phi_x|^2)_x} (({\beta N})_x)}. \end{equation} Since by the central limit theorem, \begin{equation*} \left(\frac{|{\phi_x}|^2- \E[|{\phi_x}|^2]}{\sqrt{N}}\right)_{x \in \Lambda} = \left(\frac{|{\phi_x}|^2- N G_{\Lambda, m^2}(0,0)}{\sqrt{N}}\right)_{x \in \Lambda} \end{equation*} (and similarly for $\rest{\phi}$) converges to a centered $n^d$-dimensional Gaussian vector, we expect that, by choosing $m^2 = m^2_n$ as in~\eqref{eq:choice_of_mass}, the above ratio converges to one, which would complete the proof. Motivated by this discussion, we state and prove (see Appendix~\ref{A:proof_densconv}) the following result, which can be seen as a multidimensional generalization of~\cite[Chapter~8, Section~46, Theorem~1]{gnedenko1968limit}. \begin{prop}\label{prop:dens_conv} Let $(X^i)_i$ be a sequence of i.i.d. {square-integrable}, centered random vectors with probability density function (with respect to Lebesgue measure on $\R^k$) $f: \R^k \to \R$. Assume that \begin{enumerate} \item\label{prop4.3.1} $f \in L^r(\R^k)$ for some $r \in (1,2]$; \item\label{prop4.3.2} The covariance matrix $C$ of $X^1 = (X^1_j)_{j=1}^k$ is positive definite. \end{enumerate} Then the relation \begin{align*}%\label{prop4.3.3}\tag{$3$} n^{k/2} f^{(n)} \left(\sqrt{n} v\right) \xrightarrow[]{n\to \infty} \frac{1}{(2\pi)^{k/2}\sqrt{\det C}} \: e^{-\frac{1}{2} \left(v, C^{-1}v\right)} \end{align*} holds uniformly over $v \in \R^k$. Here, $f^{(n)}$ is the $n$-fold convolution of $f$, as well as the density function of $\sum_{i=1}^n {X}^i$. \end{prop} \begin{proof}[Proof of Theorem~\ref{thm:spinO(N)_torus}] Choose $m^2> 0$ as in~\eqref{eq:mass_spinO(N)}. Then, \[ \E[(\phi^i_x)^2] = G_{\Lambda, m_n^2}(x,x) = \beta. \] We start by checking that the i.i.d. vectors $((\phi^i)^2 - \beta \coloneqq ((\phi^i_x)^2 - \beta)_{x\in \Lambda})_{i}$ satisfy assumptions of Proposition~\ref{prop:dens_conv}. As for property~\eqref{prop4.3.2}, it is clear that the entries of $C$ are well-defined since we work with squares of Gaussians and that $C$ is positive semi-definite as a covariance matrix. Furthermore, $v^T Cv = 0$ if and only if $\langle v, (\phi^i)^2-\beta\rangle = \sum_{x\in \Lambda} v_x ((\phi_x^i)^2-\beta) = 0$ almost surely, which in turn is true only for $v = 0 \in \R^{\Lambda}$. To verify~\eqref{prop4.3.1}, we derive a formula for the density function $f$ of $((\phi^1_x)^2 - \beta)_{x\in \Lambda}$ and show that it belongs to $L^r(\R^{n^d})$ for any $r \in [1,2)$. Let $f_{\phi}$ be the density function of an $m_n^2$-massive GFF on $\Lambda$, then for any $t = (t_x)_{x \in \Lambda} \in \R^\Lambda$ \begin{align*} f(t) = \left(\prod_{x \in \Lambda} \frac{1}{2\sqrt{t_x + \beta}} \right) \sum_{(k_x)_x \in \{0,1\}^{n^d}} f_{\phi} \left(\left((-1)^{k_x} \sqrt{t_x + \beta}\right)_x\right) \ind_{\{\forall x\colon t_x + \beta > 0\}}. \end{align*} Recall that the eigenvalues of $(-\Delta_\Lambda + m_n^2)$ belong to $[m_n^2, m_n^2 + 4d]$ (see Lemma~\ref{lemma:massive_EVs}). Therefore, \begin{align*} f(t) \leq \frac{2^{n^d}}{(2\pi)^{n^d/2} \sqrt{\det G_{\Lambda, m_n^2}}} \prod_{x\in \Lambda} \left(\frac{1}{2\sqrt{t_x + \beta}} e^{-\frac{1}{2} m_n^2 (t_x + \beta)} \ind_{\{t_x + \beta > 0\}}\right). \end{align*} From this explicit formula we directly see that $f \in L^r(\R^{n^d})$ for any $r \in [1, 2)$. Having verified the assumptions of Proposition~\ref{prop:dens_conv}, we now apply it to the ratio in~\eqref{eq:dens_O(N)}. \begin{multline*} \frac{f_{(|\rest{\phi_x}|^2)_x} \left(({\beta N} - \beta |\xi_x|^2)_x \right)}{f_{(|\phi_x|^2)_x} (({\beta N})_x)}\\ = \left({\frac{N}{N-M}}\right)^{{n^d}/{2}} \times \frac{(N-M)^{{n^d}/{2}} f_{(|\rest{\phi_x}|^2 - (N-M)\beta)_x} \left(\sqrt{N-M}\left(\frac{\beta M - \beta |\xi_x|^2}{\sqrt{N-M}}\right)_x \right)}{N^{n^d/2} f_{(|\phi_x|^2 - N\beta)_x} (\sqrt{N}({0})_x)}\\ \xrightarrow{N \to \infty} 1. \end{multline*} \end{proof} \subsection{The infinite-volume limit: \texorpdfstring{$N\to \infty$}{N goes to infty} followed by \texorpdfstring{$n\to \infty$}{n goes to infty}}\label{subsubsec:infdim,infvol_limit} In this subsection we discuss the subsequent infinite-volume limit of the model and prove Theorem~\ref{thm:spinO(N)_infVol}. Taking the limit $N \to \infty$ first significantly simplifies the problem. Indeed, by Theorem~\ref{thm:spinO(N)_torus}, we know that for any fixed sufficiently large $n \in \N$ and finite $U \subset \Z^d$ (viewed as a subset of $\Lambda_n$), the marginal $\ovlineM{S}_U \coloneqq (S^{i}_x)_{x \in U}^{i = 1, \dots, M}$ converges in law as $N\to \infty$ to \[ \frac{1}{\sqrt{\beta}} \ovlineM{\phi}_U \coloneqq \frac{1}{\sqrt{\beta}}(\phi^i_x)_{x\in U}^{i= 1,\dots, M}, \] where $\ovlineM{\phi}_U$ is an $m^2_n$-massive $M$-vectorial GFF on $\Lambda_n$ restricted to $U$, with $m_n^2$ determined by~\eqref{eq:mass_spinO(N)}. Since the fields $(\phi^i)_i$ are i.i.d. centered Gaussian vectors, it suffices to verify that the covariance matrix of a single coordinate process converges to that of the limiting field. \begin{proof}[Proof of Theorem~\ref{thm:spinO(N)_infVol} (Convergence of the covariance structure)] Recall that the critical inverse temperature $\beta_c \in (0, \infty]$ corresponds to $G_{\Z^d}(0,0)$. In particular, $\beta_c = \infty$ if $d=2$, $\beta_c \in (0,\infty)$ for $d\geq 3$. Let us have a closer look at~\eqref{eq:mass_spinO(N)}: using Riemann sum approximation, for $m^2> 0$ if $d=2$, or $m^2 \geq 0$ if $d\geq 3$, we get \begin{multline}\label{eq:asympt_Gdiag} \frac{1}{n^d} \sum_{x \in [0,n)^d \cap \Z^d \setminus\{0\}} \frac{1}{m^2 + 2\sum_{k=1}^d (1 - \cos(2\pi x_k/n))} \\ \overset{n\to \infty}{\sim} \int_{(0,1)^d} \frac{1}{m^2 + 2\sum_{k=1}^d (1 - \cos(2\pi x_k))} \d x = G_{\Z^d, m^2}(0,0). \end{multline} If $\beta < \beta_c$, this implies that $(m^2_n)_n$—solutions to~\eqref{eq:mass_spinO(N)}—converge to $m^2> 0$, the unique solution to $G_{\Z^d, m^2}(0,0) = \beta$. Combined with the fact that $G_{\Lambda_n, m^2}(x,y) \to G_{\Z^d, m^2}(x,y)$ for all $x,y \in U$, this completes the proof in this case. Suppose that $\beta > \beta_c$ ($d\geq 3$). By~\eqref{eq:asympt_Gdiag}, we see that in the limit the sum in~\eqref{eq:mass_spinO(N)} can be at most $\beta_c$. Hence, $m_n^2$ must be of order $\Theta(1/n^d)$. Since all the non-zero eigenvalues of $-\Delta_{\Lambda_n}$ are at least of order $\Omega(1/n^2)$, \begin{align*} \frac{1}{n^d} \sum_{x \in [0,n)^d \cap \Z^d \setminus\{0\}} \frac{1}{m_n^2 + \eta_x} = \frac{1-\cO(n^{2-d})}{n^d} \sum_{x \in [0,n)^d \cap \Z^d \setminus\{0\}} \frac{1}{\eta_x} \overset{n \to \infty}{\sim} \beta_c. \end{align*} Therefore, $m_n^2 = \frac{1 + o(1)}{(\beta - \beta_c) n^d}$, where the $o(1)$ term is uniform in large $n$. Let\linebreak $(q^w)_{w \in [0,n)^d \cap \Z^d}$ be the orthonormal eigenvectors of $-\Delta_{\Lambda_n}$ corresponding to the eigenvalues $(\eta_w)_{w \in [0,n)^d \cap \Z^d}$. In particular, $q^0_x = \frac{1}{\sqrt{n^d}}$ for any $x \in [0,n)^d \cap \Z^d$. For any $x, y \in U$, we~have \begin{equation}\label{eq:Gfct_zeroavGfct} \begin{aligned}[t] G_{\Lambda_n, m_n^2}(x,y) &= \sum_{w \in [0,n)^d \cap \Z^d \setminus\{0\}} \frac{1}{m_n^2+\eta_w} q^w_x q^w_y + q^0_x q^0_y \frac{1}{m^2_n} \\ &= (1-o(1)) \left(\sum_{w \in [0,n)^d \cap \Z^d \setminus\{0\}} \frac{1}{\eta_w} q^w_x q^w_y + (\beta - \beta_c)\right)\\ &= (1-o(1)) \left(G_{\Lambda_n}^{0\avg}(x,y) + (\beta - \beta_c) \right) \\ &\xrightarrow[\text{unif. in $x,y\in U$}]{n \to \infty} G_{\Z^d}(x,y) + (\beta - \beta_c), \end{aligned} \end{equation} where $G_{\Lambda_n}^{0\avg}$ is the zero-average Green's function on $\Lambda_n$. We conclude by observing that $(G_{\Z^d}(x,y) + (\beta - \beta_c))_{x,y \in U}$ is the covariance matrix of the sum of the GFF on $\Z^d$ restricted to $U$ and a Gaussian field of the form $(\sqrt{\beta - \beta_c} Z)_{x \in U}$, where $Z \sim \cN(0,1)$ is independent of the GFF. Let $\beta = \beta_c$. By~\eqref{eq:asympt_Gdiag}, we see that for any $\varepsilon > 0$, there exists $n_0$ such that for all $n \geq n_0$, $0 < m_n^2 < \varepsilon$. Combined with~\eqref{eq:Gfct_RW}, this implies that for any $\delta \in (0,1)$, \begin{align*} G_{\Z^d, \varepsilon}(x,y) \xleftarrow[\text{unif. over $U$}]{n \to \infty} G_{\Lambda_n, \varepsilon}(x,y) \leq G_{\Lambda_n, m_n^2}(x,y) \leq G_{\Lambda_n,\, \min\left(m_n^2, n^{-2 - \delta}\right)}(x,y). \end{align*} As $\varepsilon \to 0$, $G_{\Z^d, \varepsilon}(x,y)$ converges to $G_{\Z^d}(x,y)$ uniformly in $x,y \in U$. Moreover, analogously to~\eqref{eq:Gfct_zeroavGfct}, by~\eqref{eq:mass_spinO(N)} \begin{align*} G_{\Lambda_n, \min(m_n^2, n^{-2 - \delta})}(x,y) &= (1 - o(1)) G_{\Lambda_n}^{0\avg}(x,y) + \frac{1}{\min(m_n^2, n^{-2 - \delta}) n^d} \\ &= (1-o(1)) G_{\Lambda_n}^{0\avg}(x,y) + \cO(n^{2+ \delta - d})\\ &\qquad+ \left(\beta_c - \frac{1-\cO(m_n^2 n^2)}{n^d} G_{\Lambda_n}^{0\avg}(0,0)\right) \ind_{\{m_n^2 \leq n^{-2 - \delta}\}} \\ &\xrightarrow[\text{unif. in $x,y\in U$}]{n \to \infty} G_{\Z^d}(x,y) \end{align*} Hence, $G_{\Lambda_n, m_n^2}$ converges uniformly over $U\times U$ to $G_{\Z^d}$ as desired. \end{proof} \section{Fine estimates for the zero-average Green's function}\label{sec:zero-avgG} This section is dedicated to the proof of Theorem~\ref{thm:improved}. For reader's convenience, we restate it here. \begin{theo}\label{thm:improved_copy} Let $d \geq 3, n \in \N$, and recall that $\Lambda_n = \T^d_n$. The zero-average Green's function $G_{\Lambda_n}^{0\avg}$ satisfies the following: \begin{enumerate} \item\label{theo5.1.1} There exists $C = C(d)>0$ such that for all $n$ sufficiently large and $y \in [-n/2, n/2)^d \cap \Z^d$, \begin{equation}\label{eq:0avgGfct_polydec_copy} \left\lvert G_{\Lambda_n}^{0\avg}(0,y) - G_{\Z^d}(0,y) \right\rvert \leq C n^{2-d}. \end{equation} In particular, for $\lvert y\rvert = \Theta(n)$ as above, $G_{\Lambda_n}^{0\avg}(0,y) = \cO(n^{2-d})$. \item\label{theo5.1.2} For $d=3$ and any $x \in \Z^d$, \begin{equation}\label{eq:negativity_Greens_diff_copy} G^{0\avg}_{\Lambda_n}(x,x) - G_{\Z^d}(x,x) = \Theta(n^{2-d}) < 0 \quad \text{uniformly for large $n$}. \end{equation} \end{enumerate} \end{theo} The proof of Theorem~\ref{thm:improved_copy} is divided into two parts, each verifying one of the above statements. The proof of the estimate~\eqref{eq:0avgGfct_polydec_copy} is given in Subsection~\ref{subsec:zero_avg_polydecay} and mainly follows the argument verifying~\eqref{eq:zero-avGfct_polydec} in~\cite[Proposition~1.5]{Aba17}, but with much higher precision. As for~\eqref{eq:negativity_Greens_diff_copy}, we first obtain the following key ingredient (proven in Section~\ref{subsec:3D_bound}), which is then used to derive the desired result in the proof below. \goodbreak \begin{prop}\label{prop:0avg_bdry_bound} Let $d=3$. Consider $\hat U \coloneqq (\partial [-\lfloor n/2\rfloor, \lfloor n/2\rfloor]^d) \cap \Z^d$, and let $U$ be its canonical projection onto $\Lambda_n$. There exists an absolute constant $c_* < \frac{3^{2/3}}{2 (4\pi)^{2/3}}$\footnote{We believe that the optimal constant is non-positive.} such that \begin{align} \sup_{y \in U} G^{0\avg}_{\Lambda_n}(0,y) \leq \frac{c_*}{n} \label{eq:bdry_bound} \end{align} for all $n$ sufficiently large. \end{prop} \begin{proof}[Proof of~\eqref{eq:negativity_Greens_diff_copy} assuming Proposition~\ref{prop:0avg_bdry_bound}] Let $x = 0$, $\hat U$ and $U$ be as in Proposition~\ref{prop:0avg_bdry_bound}. Then, since the random walk on torus killed upon entering $U$ and the one on $\Z^d$ (started in the interior of $\hat U$) killed on $\hat{U}$ have the same law, by~\eqref{eq:0-avGfct_decomp} and~\eqref{eq:Gfct_decomp}: \begin{align} G^{0\avg}_{\Lambda_n} (0,0) &= G_{U^c,0}(0,0) + \E_{0}^0\left[G_{\Lambda_n}^{0\avg}(X_{H_U},0)\right] - \frac{1}{2d n^d} \E_{0}^0 [H_U] \\ &= G_{\Z^d}(0,0) - \E_{0}^{0,\Z^d} \left[G_{\Z^d}(X_{H_U},0) \ind_{\{H_U < \infty\}}\right]\nonumber \\ &\qquad+\E_{0}^0\left[G_{\Lambda_n}^{0\avg}(X_{H_U},0)\right] - \frac{1}{2d n^d} \E_{0}^0 [H_U] \nonumber\\ &= G_{\Z^d}(0,0) - \E_{0}^{0,\Z^d} [G_{\Z^d}(X_{H_U},0)]\label{eq:aux_eq_bnd}\\ &\qquad+\E_{0}^0\left[G_{\Lambda_n}^{0\avg}(X_{H_U},0)\right] - \frac{1}{2d n^d} \E_{0}^{0,\Z^d} [H_U] \nonumber \end{align} (cf.~\cite[(15)]{das2018extremal}). By the optional sampling theorem, using that under $\P^{0, \Z^d}_0$, $\lvert X_k\rvert^2 - k$ is a martingale and $H_U$ is integrable, \begin{equation}\label{eq:opt_sampl} \E_{0}^{0,\Z^d} [H_U] = \Theta(n^2) = \E_{0}^{0,\Z^d} \bigl[|X_{H_U}|^2 \bigr] \end{equation} (c.f. \cite[(1.21)]{lawler2012intersections}). By~\cite[Theorem~1.5.4]{lawler2012intersections}, we further know that for $d\geq 3$, \begin{equation}\label{eq:asympt_Greensfct} G_{\Z^d}(0,x) \sim C(d) \lvert x\rvert^{2-d} \quad \text{with} \quad C(d) = \frac{1}{2d} \frac{d}{2} \Gamma\left(\frac{d}{2} - 1\right) \pi^{-d/2}. \end{equation} By plugging in~\eqref{eq:opt_sampl}, \eqref{eq:asympt_Greensfct}, and Proposition~\ref{prop:0avg_bdry_bound}, into~\eqref{eq:aux_eq_bnd}, and noticing that $|X_{H_U}| \in [\lfloor n/2\rfloor, {\sqrt{3} n}/{2}]$, we obtain \begin{align*} G^{0\avg}_{\Lambda_n}(0,0) - G_{\Z^3}(0,0) &\leq - \inf_{r \in \left[\left\lfloor \frac{n}{2}\right\rfloor, \frac{\sqrt{3} n}{2}\right]} \left(\frac{C(3)}{r} + \frac{1}{6n^3} r^2\right) + \frac{c_*}{n} \\ &= \left(c_* - \frac{(3 C(3))^{2/3}}{2}\right) \frac{1}{n} = -\frac{1}{n} \left(\frac{3^{2/3}}{2 (4\pi)^{2/3}} - c_* \right)\eqqcolon - \frac{c_{**}}{n} \end{align*} with an absolute constant $c_{**} > 0$, as desired. The infimum is attained at \[ r = (3 C(3))^{1/3} n = \left(\frac{3}{4\pi}\right)^{1/3} n \in \left[\frac{n}{2}, \frac{\sqrt{3} n}{2}\right].\qedhere \] \end{proof} \subsection{Proof of~\eqref{eq:0avgGfct_polydec_copy}: Polynomial decay of zero-average Green's function}\label{subsec:zero_avg_polydecay} This subsection proves~\eqref{eq:0avgGfct_polydec_copy}, namely that for $d\geq 3$, there exists $C = C(d)>0$ such that for all $n$ sufficiently large and $y \in \hat \Lambda \coloneqq [-n/2, n/2)^d \cap \Z^d$, \[ \left\lvert G_{\Lambda}^{0\avg}(0,y) - G_{\Z^d}(0,y) \right\rvert \leq C n^{2-d}. \] Here and in what follows, we write $\Lambda \coloneqq \T^d_n$ for brevity. \begin{proof}[Proof of~\eqref{eq:0avgGfct_polydec_copy}] Recall that (see Proposition~\ref{prop:props:zeroAvGfct}), \begin{align*} 2d G_{\Lambda}^{0\avg}(0,y) = \int_0^\infty \Big(\P^0_{\Lambda}[\bar{X}_t = y] - \frac{1}{n^d} \Big) \d t, \end{align*} where $\bar{X} = (\bar{X}_t)_{t \geq 0}$ is the continuous-time simple random walk on $\Lambda$ viewed as a~graph. To obtain the desired estimate we follow the strategy of the proof of~\cite[Proposition~1.5]{Aba17} corresponding to our~\eqref{eq:zero-avGfct_polydec} and split the above integral approximately at the time $\bar{{X}}$ reaches equilibrium (the uniform distribution). More precisely, as in the reference (with $N = n$ in our setting), for \[ \lambda_* \coloneqq \frac{1}{d}(1-\cos(2\pi/n)) \overset{n \to \infty}{\sim} \frac{2 \pi^2}{d n^2} \] (the spectral gap of $\bar{{X}}$) and $t_* \coloneqq \frac{\log(n^d)}{\lambda_*}$ by~\cite[Theorem~20.6]{levin2017markov}, \begin{align*} \int_{t_*}^\infty \left\lvert \P^0_{\Lambda}[\bar{X}_t = y] - \frac{1}{n^d} \right\rvert \d t \leq \int_{t_*}^\infty e^{-\lambda_* t} \d t = \frac{1}{\lambda_* n^d} = \Theta(n^{2-d}). \end{align*} The remaining part of the integral will be treated similarly to the referenced proof but with higher precision tailored to our setting. Note that \begin{align*} &\int_0^{t_*} \left(\P^0_{\Lambda}[\bar{X}_t = y] - \frac{1}{n^d} \right) \d t\\ &= \int_0^{t_*} \sum_{k=0}^\infty \P^0_{\Lambda} [N_t = k] \P^0_0 [X_k = y] \d t - \frac{t_*}{n^d} \in\mkern-15mu \sum_{k=0}^{t_* + \sqrt{4 t_* \log(n^d)}} \P^0_0 [X_k = y] \int_0^{t_*} \frac{t^k e^{-t}}{k!} \d t - \frac{t_*}{n^d}\\ &\mkern320mu + \left[0, \int_0^{t_*} \P^0_{\Lambda} [N_t > t_* + \sqrt{4 t_* \log(n^d)}] \d t\right]. \end{align*} Here $N_t$ is the number of jumps of $(\bar{X}_s)_s$ up to time $t$. We now show that $\int_0^{t_*} \P^0_{\Lambda} [N_t > t_* + \sqrt{4 t_* \log(n^d)}] \d t = o(n^{2-d})$ as $n \to \infty$. For $t \in (0, t_*)$ by exponential Markov's inequality, \begin{equation}\label{eq:Pois_tails} \begin{aligned}[t] \P^0_{\Lambda} \left[N_t > t_* + \sqrt{4 t_* \log(n^d)}\right] &\leq \inf_{\lambda > 0} \exp \left(t(e^\lambda - 1) - \lambda \left(t_* + \sqrt{4 t_* \log(n^d)}\right) \right)\\ &\leq \exp \left(t_* \inf_{\lambda > 0} \left[(e^\lambda - 1) - \lambda \left(1 + \sqrt{4 \lambda_*}\right) \right] \right) \\ &= \exp \biggl(t_* \biggl[\sqrt{4 \lambda_*} - \underbrace{\log \bigl(1 + \sqrt{4 \lambda_*}\bigr)}_{\geq \sqrt{4 \lambda_*} - 4\lambda_*/2} \bigl(1 + \sqrt{4 \lambda_*}\bigr) \biggr] \biggr)\\ &\leq \exp \bigl(- 2 t_* \lambda_*(1 - o(1))\bigr) \leq (n^d)^{-3/2}. \end{aligned} \end{equation} And thus, $\int_0^{t_*} \P^0_{\Lambda} [N_t > t_* + \sqrt{4 t_* \log(n^d)}] \d t \leq t_* n^{-3d/2} = \Theta(n^{2-3d/2} \log(n)) = o(n^{2-d})$ as desired. Furthermore, for $0 \leq k \leq t_* + \sqrt{4 t_* \log(n^d)} \eqqcolon t_*'$, \begin{align*} \P^0_0 [X_k = y] = \sum_{\substack{v \in \Z^d:\\ \d_{\Z^d}(0, y+nv) \leq t_*'}} \P^{0, \Z^d}_0 [X_k = y + nv]. \end{align*} We split the latter sum into two parts: \[ \sqrt{t_*'} C_d \sqrt{\log(n^d)} \leq \d_{\Z^d}(0, y+nv) \leq t_*'\quad \text{and}\quad\d_{\Z^d}(0, y+nv) \leq \sqrt{t_*'} C_d \sqrt{\log(n^d)} \] for some $C_d > 0$ to be determined later. We show that the contribution of the former fragment of the sum to $\sum_{k=0}^{t_*'} \P^0_0 [X_k = y] \int_0^{t_*} \frac{t^k e^{-t}}{k!} \d t$ can be made as small as $o(n^{2-d})$ if $C_d$ is chosen appropriately (the lower bound is clearly zero). For some appropriate absolute (depending only on the dimension) constants $c, c', r_0 > 0, R_0 < \infty$, \begin{multline*} \underbrace{\sum_{k=0}^{t_*'} \int_0^{t_*} \frac{t^k e^{-t}}{k!} \d t}_{\leq t_*} \underbrace{\sum_{\substack{v \in \Z^d:\sqrt{t_*'} C_d \sqrt{\log(n^d)} \leq \\ \d_{\Z^d}(0, y+nv) \leq t_*'}}}_{\leq c (t_*'/n)^d} \P^{0, \Z^d}_0 [X_k = y + nv] \\ \begin{aligned} &\leq \P^{0, \Z^d}_0 \left[\max_{k=0}^{\lfloor t_*'\rfloor} \d_{\Z^d}(0, X_k) \geq \sqrt{t_*'} C_d \sqrt{\log(n^d)}\right] \underbrace{t_* c(t_*'/n)^d}_{= \Theta\left(n^{2+d} \log^{d+1}(n)\right)}\\ &\leq c'n^{2+d} \log^{d+1}(n) \P^{0, \Z^d}_0 \leq \P^{0, \Z^d}_0 \left[\max_{k=0}^{\lfloor t_*'\rfloor} \lvert X_k\rvert \geq (C_d/d) \sqrt{t_*' \log(n^d)}\right]\\ &\leq c'n^{2+d} \log^{d+1}(n) R_0 e^{-r_0(C_d/d)^2 \log(n^d)} = \Theta\left(n^{2+d - r_0C_d^2/d} \log^{d+1}(n)\right). \end{aligned} \end{multline*} The last inequality (and in particular, existence of $r_0, R_0$ as above) follows from~\cite[Proposition~2.1.2$\MK$(b)]{lawler_limic_2010}. By choosing $C_d$ such that $C_d > d\sqrt{2/r_0}$, say, $C_d = (2d/\sqrt{r_0}) \vee 1$, we get that the latter expression is $o(n^{2-d})$. It remains to bound \begin{align*} I_n \coloneqq \sum_{k=0}^{t_*'} \sum_{\substack{v \in \Z^d: \d_{\Z^d}(0, y+nv) \\\qquad \leq C_d \sqrt{t_*'\log(n^d)}}} \P^{0, \Z^d}_0 [X_k = y + nv] \int_0^{t_*} \frac{t^k e^{-t}}{k!} \d t - \frac{t_*}{n^d}. \end{align*} First observe that $\int_0^{t_*} \frac{t^k e^{-t}}{k!} \d t = \P[\Pois(t_*) \geq k+1]$ (by definition of the upper incomplete gamma function and explicit expression of Poisson cumulative distribution function in terms of it). Therefore, \begin{align*} I_n \leq \sum_{k=0}^{t_*'} \sum_{\substack{v \in \Z^d: \d_{\Z^d}(0, y+nv) \\\qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \P^{0, \Z^d}_0 [X_k = y + nv] - \frac{t_*}{n^d}, \end{align*} but also \begingroup \allowdisplaybreaks \begin{align*} I_n &\geq \P\left[\Pois(t_*) \geq t_*-\sqrt{4 t_* \log(n^d)}\right] \\ &\quad\times \sum_{k=0}^{t_*-\sqrt{4 t_* \log(n^d)}} \sum_{\substack{v \in \Z^d: \d_{\Z^d}(0, y+nv) \\ \qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \P^{0, \Z^d}_0 [X_k = y + nv] - \frac{t_*}{n^d}\\ &\geq (1- n^{-3d/2}) \sum_{k=0}^{t_*-\sqrt{4 t_* \log(n^d)}} \sum_{\substack{v \in \Z^d: \d_{\Z^d}(0, y+nv) \\\qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \P^{0, \Z^d}_0 [X_k = y + nv] - \frac{t_*}{n^d}. \end{align*} \endgroup In the latter inequality we have used that fully analogously to~\eqref{eq:Pois_tails}, $\P[\Pois(t_*) < t_*-\sqrt{4 t_* \log(n^d)}] \leq n^{-3d/2}$. Define $p_k: \R^d \to \R$ by $p_k(x) \coloneqq 2 (\frac{d}{2\pi k})^{d/2} e^{-\frac{d\lvert x\rvert^2}{2k}}$ and write $x \leftrightarrow k$ if $k$ and $x$ have the same parity, i.e., if $k + x_1 + \dots + x_d$ is even. Set further $E_k(x)\coloneqq (\P^{0, \Z^d}_0[X_k = x] - p_k(x))\ind_{\{k \leftrightarrow x\}}$. By~\cite[Lemma~1.5.2]{lawler2012intersections}, $\lvert x\rvert^\alpha \sum_{k \geq 0} \lvert E_k(x)\rvert \xrightarrow{\lvert x\rvert \to \infty} 0$ for every $\alpha < d$. Note that for $y \in \hat{\Lambda}$ and $v \in \Z^d \setminus\{0\}$, $\lvert y+nv\rvert = \Omega(n) \geq n/2$. Altogether, this implies that for $t \in \{t_*', t_* - \sqrt{4 t_* \log(n^d)}\}$ and $\alpha \in (d-2, d)$, \begin{multline*} \sum_{k=0}^{t} \sum_{\substack{v \in \Z^d: \d_{\Z^d}(0, y+nv) \\\qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \P^{0, \Z^d}_0 [X_k = y + nv] \\ = \underbrace{\sum_{k=1}^{t} \sum_{\substack{v \in \Z^d\setminus \{0\}: \d_{\Z^d}(0, y+nv) \\\qquad\leq C_d \sqrt{t_*'\log(n^d)}}} p_k(y + nv) \ind_{\{k \leftrightarrow y+nv\}}}_{\eqqcolon J_n(t)} \\ + \sum_{k=0}^t \P^{0, \Z^d}_0 [X_k = y] \pm \underbrace{\Theta((\log n)^d) o(n^{-\alpha})}_{= o(n^{2-d})}. \end{multline*} \begin{enonce*}[plain]{Claims} \leavevmode \begin{enumerate}\claimenumi \item\label{claimbis1} $0 \leq 2d G_{\Z^d}(0,y) - \sum_{k=0}^t \P^{0, \Z^d}_0 [X_k = y] \leq o(n^{2-d})$. \item\label{claimbis2} $\lvert J_n(t) - \frac{t_*}{n^d}\rvert = \cO(n^{2-d})$. \end{enumerate} \end{enonce*} Note that combined with the above arguments these claims would imply the desired result. Indeed, by the above and since $G_{\Z^d}(0,y)$ is finite for all $y \in \Z^d$, \begingroup \allowdisplaybreaks \begin{align*} -\cO(n^{2-d}) + \underbrace{(1 - n^{-3d/2})\Big(J_n + \sum_{k=0}^t \P^{0, \Z^d}_0 [X_k = y]\Big) - \frac{t_*}{n^d}}_{\geq 2d G_{\Z^d}(0,y) - \cO(n^{2-d})} \leq 2d G^{0\avg}(0,y) \\ \leq \cO(n^{2-d}) + \underbrace{J_n - \frac{t_*}{n^d} + \sum_{k=0}^t \P^{0, \Z^d}_0 [X_k = y]}_{\leq 2d G_{\Z^d}(0,y) + \cO(n^{2-d})}. \end{align*} \endgroup Hence, $\lvert G^{0\avg}(0,y) - G_{\Z^d}(0,y)\rvert = \cO(n^{2-d})$. \pagebreak \begin{proof}[{Proof of Claim~\ref{claimbis1}}] Recall that $2d G_{\Z^d}(0,y) = \sum_{k\geq 0} \P^{0, \Z^d}_0 [X_k = y]$. Therefore, it only remains to check that the tail ($k\geq t$) of this sum decays at rate $o(n^{2-d})$. This follows almost immediately from the following three observations: $p_k(y) \leq C(d) r^{-\frac{d}{2}}$, $\lvert E_k(y)\rvert \leq \cO(k^{-\frac{d}{2}-1})$ by~\cite[Theorem~1.2.1$\MK$(1.10)]{lawler2012intersections} and $t = \Omega(n^2 \log n)$, which together imply that the tail is bounded by $\int_{C n^2 \log n}^\infty (r^{-\frac{d}{2}-1} + r^{-\frac{d}{2}}) \leq \cO(n^{2-d}/(\log n)^{\frac{d}{2} - 1}) = o(n^{2-d})$ as $d\geq 3$. \end{proof} \begin{proof}[{Proof of Claim~\ref{claimbis2}}] We first show that the contribution of $k < \varepsilon n^2$ terms (for any $0 < \varepsilon < \frac{1}{4}$) is $\cO(n^{2-d})$. Indeed, since $\R_+ \ni r \mapsto p_r(z)$ is increasing on $(0, \lvert z\rvert^2)$ for any $z\in \R^d$ and $\lvert y + n v\rvert^2 \geq \frac{n^2}{4}$ for $v \in \Z^d \setminus\{0\}$, \begin{align*} \sum_{k=1}^{\varepsilon n^2} \sum_{\substack{v \in \Z^d\setminus \{0\}\,:\, \d_{\Z^d}(0, y+nv) \\ \qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \hspace{-0.6cm} p_k(y + nv) &\leq \hspace{-0.6cm} \sum_{\substack{v \in \Z^d\setminus \{0\}\,:\, \d_{\Z^d}(0, y+nv) \\ \qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \int_{1}^{2\varepsilon n^2} 2 \left(\frac{d}{2\pi r}\right)^{d/2} e^{-\frac{d\lvert y + n v\rvert^2}{2r}} \d r\\ &\leq \sum_{\substack{v \in \Z^d\setminus \{0\}\,:\, \d_{\Z^d}(0, y+nv) \\\qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \hspace{-0.6cm} c_d \lvert y + n v\rvert^{2-d} \int_{\frac{d \lvert y + n v\rvert^2}{4 \varepsilon n^2}}^\infty s^{\frac{d}{2} - 2} e^{-s} \d s\\ &\leq \tilde c_d n^{2-d} \varepsilon^{2 - \frac{d}{2}} \hspace{-0.6cm} \sum_{\substack{v \in \Z^d\setminus\{0\}\,:\, \d_{\Z^d}(0, y+nv) \\\qquad\leq C_d \sqrt{t_*'\log(n^d)}}} \left\lvert \frac{y}{n} + v\right\rvert^{-2} e^{- \frac{d}{4\varepsilon} \lvert \frac{y}{n} + v\rvert^2}\\ &\leq c_d' n^{2-d} \varepsilon^{2 - \frac{d}{2}} \prod_{i=1}^d \left(\sum_{v_i = - 2\tilde C_d \log n}^{2 \tilde C_d \log n} e^{- \frac{d}{4\varepsilon} \left(\frac{y_i}{n} + v_i\right)^2} \right)\\ &\leq c_d'' n^{2-d} \varepsilon^{2 - \frac{d}{2}} \left(\int_{\R} e^{- \frac{d}{4\varepsilon} x^2} \d x + 1\right)^d = \hat{c}_d n^{2-d} \varepsilon^{2 - \frac{d}{2}}. \end{align*} In the third line we used that for all $x > 1$, \begin{align*} \int_{x}^\infty s^{d/2} e^{-s}\d s \leq x^{d/2-2} e^{-x} \left(1 + \sum_{k=2}^{\lfloor d/2\rfloor} \prod_{j=2}^k \left(\frac{d}{2} -j\right)\right). \end{align*} Note also that there exists an absolute constant $c > 0$ such that for all $\lvert z\rvert \leq 2 C_d \sqrt{t_*' \log(n^d)}, t_*' \geq s \geq \varepsilon n^2$: \[ p_{\lfloor s\rfloor}(z) / p_{s}(z) \in \left(1 - c \frac{\log^2 n}{\varepsilon^2 n^2}, 1 + c \frac{\log^2 n}{\varepsilon^2 n^2}\right). \] Since $\lvert y+nv \rvert \leq \d_{\Z^d}(0, y+nv) \leq d \lvert y+nv \rvert$, and $\frac{1}{2} p_k$ is the density function of $\cN(0, \frac{k}{d} \Id_d)$, this yields that \begingroup \allowdisplaybreaks \begin{align*} J_n(t) &\lessgtr_a \left(1 \pm c \frac{\log^2 n}{\varepsilon^2 n^2}\right) \int_{\varepsilon n^2}^{t} \d s \int_{\cB^d\left(\frac{-y}{n}, a C_d \sqrt{\frac{t_*' \log(n^d)}{n^2}} \right) \setminus \left[-\frac{1}{2}, \frac{1}{2}\right)^d} \\ &\qquad\frac{p_s\left(n \left(\frac{y}{n} + \lfloor v\rfloor \right)\right)}{2} \d v + \cO (n^{2-d}),\\ &= \frac{n^2}{n^d} \left(1 \pm c \frac{\log^2 n}{\varepsilon^2 n^2}\right) \int_{\varepsilon}^{\frac{t}{n^2}} \d r \int_{\cB^d\left(\frac{-y}{n}, aC_d \sqrt{\frac{t_*' \log(n^d)}{n^2}} \right) \setminus \left[-\frac{1}{2}, \frac{1}{2}\right)^d} \\ &\qquad \frac{p_{r}\left(\frac{y}{n} + \lfloor v\rfloor \right)}{2} \d v + \cO(n^{2-d}) \end{align*} \endgroup where $\lfloor v\rfloor \in \Z^d$ stands for the point on the grid $\Z^d$ such that $v \in \lfloor v\rfloor + [-1/2,1/2)^d$ and $\lessgtr_a$ here means that $J_n(t) \leq$ r.h.s. with $a= 2$, non-negative $\cO(n^{2-d})$ and $+$ in $1 \pm c \frac{\log^2 n}{\varepsilon^2 n^2}$ and $J_n(t) \geq $ r.h.s. for $a= 1/(2d)$ with $\cO(n^{2-d}) = 0$ and $-$ in $1 \pm c \frac{\log^2 n}{\varepsilon^2 n^2}$. Observe that \begin{align*} 0 \leq \frac{n^2}{n^d} \int_{\varepsilon}^{\frac{t}{n^2}} \d r \frac{p_{r}\bigl(\frac{y}{n}\bigr)}{2} \leq C_d \frac{n^2}{n^d} \int_{\varepsilon}^{\frac{t}{n^2}} r^{-\frac{d}{2}} \d r \leq \cO(n^{2-d}). \end{align*} Therefore, we can add $\bigl[-\frac{1}{2}, \frac{1}{2}\bigr)^d$ back to the integral by adding $-\cO(n^{2-d})$ to the lower bound of $J_n(t)$, for the upper bound other than that nothing changes. Using first less precise bounds, we further show that \begin{equation}\label{eq:J_n_term} J_n(t) = \frac{n^2}{n^d} \int_{\varepsilon}^{\frac{t_*}{n^2}} \d r \int_{\R^d} \frac{p_{r}\bigl(\frac{y}{n} + \lfloor v\rfloor \bigr)}{2} \d v + \cO(n^{2-d}). \end{equation} By the mean-value theorem, for any $v \in \cB^d(\frac{-y}{n}, 2C_d \sqrt{\frac{t_*' \log(n^d)}{n^2}} )$, there exists $c_v \in [0,1]$ such that \begin{multline*} \left\lvert p_{r} \left(\frac{y}{n} + \lfloor v\rfloor \right) - p_{r} \left(\frac{y}{n} + v \right)\right\rvert \\ \begin{aligned} &\leq p_{r} \left(\frac{y}{n} + (1-c_v) v + c_v \lfloor v\rfloor \right) \frac{d \left\lvert\frac{y}{n} + (1-c_v) v + c_v \lfloor v\rfloor \right\rvert}{r} \lvert v - \lfloor v\rfloor\rvert \\ &\leq e^{\frac{d}{2r}c_v^2 \lvert v - \lfloor v\rfloor\rvert^2} p_{r} \left(\frac{1}{\sqrt{2}} \left(\frac{y}{n} + v\right) \right) d^{3/2} \frac{\left\lvert\frac{y}{n} + v \right\rvert + \sqrt{d}}{r} \\ &\leq c_d(\varepsilon) p_{r}\left(\frac{1}{\sqrt{2}} \left(\frac{y}{n} + v\right) \right) \frac{\left\lvert\frac{y}{n} + v \right\rvert + 1}{r}. \end{aligned} \end{multline*} Furthermore, for any $C > 0$, $0 < s \leq t_*'$, \begin{align*} \P\left[\ntwo{\cN (0, (s/d) \Id_d)} \geq C \sqrt{t_*' \log(n^d)}\right] &\leq \P\left[\ntwo{\cN(0, \Id_d)} \geq C \sqrt{d\log(n^d)}\right]\\ &\leq 1 - \P\left[\lvert \cN(0,1)\rvert \leq C \sqrt{\log(n^d)}\right]^d \\ &\leq 1 - \left(1 - n^{-d C^2}\right)^d = \cO(n^{-d C^2}). \end{align*} Hence, for $R \coloneqq a C_d \sqrt{\frac{t_*' \log(n^d)}{n^2}}$, \begin{multline*} \int_{\cB^d\left(\frac{-y}{n}, R\right)} \frac{p_{r}\left(\frac{y}{n} + \lfloor v\rfloor \right)}{2} \d v\\ \begin{aligned} &= \int_{\cB^d(0, R )} \frac{p_{r} (u)}{2} \d u \pm \frac{c_d(\varepsilon)}{r} \cO \left(\int_{\cB^d\left(0, R \right)} \frac{ p_{r} (u/\sqrt{2}) (\lvert u\rvert + 1)}{2} \d u \right)\\ &= \P[\ntwo{\cN(0, (r/d) \Id_d)} \leq R] \pm \frac{c_d(\varepsilon)}{r} \cO \left(1 + \E[\ntwo{\cN(0, (r/d) \Id_d)}]\right) \\ &= 1 - \cO\left(n^{-d C_d^2}\right) \pm \frac{1}{\sqrt{r}} \cO(1) \end{aligned} \end{multline*} and analogously, \begin{multline*} \int_{\cB^d\left(\frac{-y}{n}, R\right)^c} \frac{p_{r}\left(\frac{y}{n} + \lfloor v\rfloor \right)}{2} \d v\\ \begin{aligned} &= \P[\ntwo{\cN(0, (r/d) \Id_d)} \geq R] (1+ \cO(1/r)) \\ &\quad\pm \frac{1}{r} \cO \left(\E\left[\ntwo{\cN(0, (r/d) \Id_d)}^2\right]^{1/2} \P\left[\ntwo{\cN(0, (r/d) \Id_d)} \geq R\right]^{1/2}\right) \\ &= \cO\left(n^{-d C_d^2/2}\right) \left(1 \pm \frac{1}{\sqrt{r}}\right) \end{aligned} \end{multline*} where we additionally used Cauchy inequality. This in turn implies that \begin{alignat*}{2} \frac{\log^2 n}{n^d} \int_{\varepsilon}^{\frac{t}{n^2}} \d r \int_{\cB^d\left(\frac{-y}{n}, R\right)} \frac{p_{r}\left(\frac{y}{n} + \lfloor v\rfloor \right)}{2} \d v &= \cO\left(\frac{t}{n^2} \frac{\log^2 n}{n^d} \right) &= o(n^{2-d}); \\ \frac{n^2}{n^d} \int_{\varepsilon}^{\frac{t}{n^2}} \d r \int_{\cB^d\left(\frac{-y}{n}, R\right)^c} \frac{p_{r}\left(\frac{y}{n} + \lfloor v\rfloor \right)}{2} \d v &= \cO\left(n^{-d C_d^2/2} \frac{t}{n^d}\right) &= o(n^{2-d}); \\ \frac{n^2}{n^d} \int_{\frac{t_*}{n^2}}^{\frac{t}{n^2}} \d r \int_{\cB^d\left(\frac{-y}{n}, R\right)} \frac{p_{r}\left(\frac{y}{n} + \lfloor v\rfloor \right)}{2} \d v &= \cO\left(\frac{t - t_*}{n^d}\right) &= o(n^{2-d}). \end{alignat*} Hence, \eqref{eq:J_n_term}. Note that the above arguments as well imply that $\lvert J_n(t) - \frac{t_*}{n^d}\rvert = \cO (n^{2-d} \sqrt{\frac{t_*}{n^2}}) = \cO(n^{2-d} \sqrt{\log n})$. We now prove that $\lvert J_n(t) - \frac{t_*}{n^d}\rvert = \cO(n^{2-d})$. By Taylor's theorem for multivariate functions, for $k \in \Z^d, x \coloneqq \frac{y}{n} + k$ and $u, w \in B \coloneqq [-1/2, 1/2)^d$: \begin{align*} p_{r} (x) = p_{r} (x + u) - u^T\nabla p_r (x) - u^T \int_0^1 (1-t) D^2 p_r (x + tu) \d t \: u. \end{align*} Note that the second term on the right-hand side is symmetric w.r.t. the origin as a function in $u$. Hence, if we integrate both sides w.r.t. $u$ and $w$ over $B$, \begin{multline*} p_r(x) - \int_{B} p_{r} (x + u) \d u + \int_B \int_{B}\int_0^1 (1-t) u^T D^2 p_r (x + tu + w) \:u \:\d t \d u \d w \\ = \int_B \int_{B} \int_0^1 (1-t) u^T \bigl(D^2 p_r (x + tu + w) - D^2 p_r (x + tu)\bigr) u\: \d t \d u \d w. \end{multline*} since $\lvert B\rvert = 1$. Observe that $(D^2 p_r (z))_{i,j} = (\frac{d}{r})^2 p_r(z) (z_i z_j - \frac{r}{d} \delta_{i,j})$ with the maximal eigenvalue $(\frac{d}{r})^2 p_r(z) (\lvert z\rvert^2 - \frac{r}{d})$ and the remaining ones equal to $- \frac{d}{r} p_r(z)$. Thus, since $r > \varepsilon > 0$, \begingroup \allowdisplaybreaks \begin{multline*} \left\lvert \int_{B}\int_0^1 (1-t) u^T D^2 p_r (z + tu) \:u \:\d t \d u \right\rvert \\ \begin{aligned} &\leq \int_{B}\int_0^1 \lvert u\rvert^2 \left(\frac{d}{r}\right)^2 p_r(z+ tu) \left(2\lvert z\rvert^2 + 2\lvert u\rvert^2 + \frac{r}{d}\right) \d t \d u\\* &\leq \left(\frac{d}{r}\right)^2 \frac{d}{4} \left(2\lvert z\rvert^2 + \frac{d}{2} + \frac{r}{d}\right) \int_{B}\int_0^1 e^{\frac{d}{2r} t^2\lvert u\rvert^2} p_r\left(\frac{z}{\sqrt{2}}\right) \d t \d u\\ &\leq \left(\frac{d}{r}\right)^2 \left(2\lvert z\rvert^2 + \frac{d}{2} + \frac{r}{d}\right) c_d(\varepsilon) p_r\left(\frac{z}{\sqrt{2}}\right). \end{aligned} \end{multline*} \endgroup The latter as a function in $z$ is clearly integrable, hence, by Fubini's theorem \begin{multline*} \sum_{x \in \frac{y}{n} + \Z^d} \left(\int_{B} p_{r} (x + u) \d u - \int_B \int_{B}\int_0^1 (1-t) u^T D^2 p_r (x + tu + w) \:u \:\d t \d u \d w \right) \\ = \int_{\R^d} p_r(z) \d z - \int_{B}\int_0^1 (1-t) u^T \int_{\R^d} D^2 p_r (z+ tu) \d z\:u \:\d t \d u =2. \end{multline*} Here we also used that $\frac{1}{2} p_r(z)$ is a density function of $\cN(0, \frac{r}{d} \Id_d)$ and so, \[ \int_{\R^d} p_r(z) \left(z_i z_j - \frac{r}{d} \delta_{i,j}\right) = 2\left(\frac{r}{d} \delta_{i,j} - \frac{r}{d} \delta_{i,j}\right) = 0. \] Moreover, by the mean-value theorem, for any $i, j \leq d$, $z \in \R^d$, \begingroup \allowdisplaybreaks \begin{align*} &\left\lvert \left(D^2 p_r (z + w) - D^2 p_r (z)\right)_{i,j}\right\rvert \leq \sup_{c \in [0,1]} \left\lvert \nabla (D^2 p_r)_{i,j} (z+ cw)\right\rvert \lvert w\rvert \\ &= \left(\frac{d}{r}\right)^3 \sup_{c \in [0,1]} p_r(z+ cw) \left\lvert \left(-\prod_{l \in \{i,j,k\}}(z+ cw)_l + \frac{r}{d}\sum_{l \in \{i,j,k\}} (z + cw)_l \delta_{\{i,j,k\} \setminus \{l\}} \right)_k \right\rvert \lvert w\rvert \\ &\leq \left(\frac{d}{r}\right)^3 p_r \left(\frac{z+ w}{\sqrt{2}}\right) \frac{\sqrt{d}}{2} \sup_{c \in [0,1]} e^{\frac{d}{2r} (1-c)^2 \lvert w\rvert^2} \left(\lvert z+ cw\rvert^3 + 3 \sqrt{d} \frac{r}{d}\lvert z+ cw\rvert \right) \\ &\leq C(d, \varepsilon) \left(\frac{d}{r}\right)^3 p_r \left(\frac{z+ w}{\sqrt{2}}\right) \left(\lvert z + w\rvert^3 + \lvert z + w\rvert^2 + \frac{r}{d} \lvert z + w\rvert + \lvert z + w\rvert + \frac{r}{d} + 1 \right)\\ &\leq \tilde C(d, \varepsilon) \left(\frac{d}{r}\right)^3 p_r \left(\frac{z+ w}{\sqrt{2}}\right) \left(\lvert z + w\rvert^3 + \frac{r}{d} \lvert z + w\rvert + \frac{r}{d}\right). \end{align*} \endgroup For us, of interest is $z = x + tu$. Note that by analogous estimation procedures we can further conclude that \begin{multline*} p_r \left(\frac{z+ w}{\sqrt{2}}\right) \left(\lvert z + w\rvert^3 + \frac{r}{d} \lvert z + w\rvert + \frac{r}{d}\right) \\ \leq C(d, \varepsilon) p_r \left(\frac{x+ w}{2}\right) \left(\lvert x + w\rvert^3 + \frac{r}{d} \lvert x + w\rvert + \frac{r}{d}\right) \end{multline*} This then yields that \begin{multline*} \sum_{x \in \frac{y}{n} + \Z^d} \int_B \left\lvert \int_{B} \int_0^1 (1-t) u^T \left(D^2 p_r (x + tu + w) - D^2 p_r (x + tu)\right) u\: \d t \d u \right\rvert \d w\\ \begin{aligned} &\leq C'(d, \varepsilon) \left(\frac{d}{r}\right)^3 \sum_{x \in \frac{y}{n} + \Z^d} \int_B p_r \left(\frac{x+ w}{2}\right) \left(\lvert x + w\rvert^3 + \frac{r}{d} \lvert x + w\rvert + \frac{r}{d}\right) \d w \\ &\leq \hat C(d, \varepsilon) \left(\frac{d}{r}\right)^3 \left(\E\left[\left\lVert\cN\left(0, \frac{r}{d} \Id_d \right)\right\rVert^3 \right] + \frac{r}{d} \E\left[\left\lVert\cN\left(0, \frac{r}{d} \Id_d \right)\right\rVert \right] + \frac{r}{d} \right) \\ &\leq \Bar C(d, \varepsilon) r^{-3/2}. \end{aligned} \end{multline*} Altogether, this part combined with~\eqref{eq:J_n_term} show that \begin{align*} J_n(t) = \frac{n^2}{n^d} \int_{\varepsilon}^{\frac{t_*}{n^2}} \bigl(1 \pm \cO(1) r^{-3/2} \bigr) \d r + \cO(n^{2-d}) = \frac{t_*}{n^d} \pm \cO (n^{2-d}) \end{align*} with $\cO(1)$ uniform in $r > \varepsilon$. This concludes the proof of Claim~\ref{claimbis2} as desired. \end{proof}\let\qed\relax \end{proof} \subsection{Proof of Proposition~\ref{prop:0avg_bdry_bound}: {``}Boundary{''} estimate for the zero-average Green's function in 3D}\label{subsec:3D_bound} In this subsection, we prove Proposition~\ref{prop:0avg_bdry_bound}, namely that for $d=3$, \begin{align*} \exists \, c_* < \frac{3^{2/3}}{2 (4\pi)^{2/3}} \,\text{ such that }\, G^{0\avg}_{\Lambda}(0,y) \leq \frac{c_*}{n} \, \text{ uniformly over } \, y \in U, \end{align*} where $U$ is the canonical projection of $\hat U = (\partial [-\lfloor n/2\rfloor, \lfloor n/2\rfloor]^d) \cap \Z^d$ onto $\Lambda \coloneqq \T^d_n$. \begin{proof}[Proof of Proposition~\ref{prop:0avg_bdry_bound}] Without loss of generality, due to the symmetries of the torus, we may assume that $y_1 \geq \dots \geq y_d \geq 0$. In particular, since $y \in \hat{U}$, $y_1 = \left\lfloor \frac{n}{2}\right\rfloor$. In Section~\ref{subsec:zero_avg_polydecay} we showed that \begin{align*} 2d \: n^{d-2} G^{0\avg}_\Lambda(0,y) &\leq \frac{d}{2\pi^2} + o(1) + c_d \varepsilon^{2-\frac{d}{2}} + \int_{\varepsilon}^{\frac{t_*}{n^2}} \d r \left(\sum_{\substack{k \in \Z^3 \\ \forall i: \lvert k_i\rvert \leq C_d \log n}} q_{r}\left(\frac{y}{n} + k \right) - 1\right), \end{align*} where $t_* = \frac{d \log(n^d)}{2\pi^2} n^2$ and $q_r(x) = (\frac{d}{2\pi r})^{d/2} e^{-\frac{d\lvert x\rvert^2}{2r}}$. We begin by verifying that even if $n$ is odd, we can set $y_1 = n/2$, which would lead to an error of order at most $o(1)$. Note that since $\lvert k_1\rvert \leq C_d \log n$, $r >\varepsilon > 0$, \begin{align*} q_{r}\left(\frac{y}{n} + \left(\frac{1}{2n}, 0, \dots, 0\right) + k \right) &= q_{r}\left(\frac{y}{n} + k \right) e^{-\frac{d}{2r}\left(\frac{1}{(2n)^2} + \left(k_1+ \frac{y_1}{n}\right) \frac{1}{n}\right)} \\ &= q_{r}\left(\frac{y}{n} + k \right) \left(1 + \cO\left(\frac{\log n}{\varepsilon n}\right)\right). \end{align*} Since we have already proven in Section~\ref{subsec:zero_avg_polydecay} that $\sum_{\substack{k \in \Z^3}} q_{r}(\frac{y}{n} + k) \leq 1 + \cO(1) r^{-3/2} \leq c(\varepsilon)$, the above implies that \begin{equation*} \cO\left(\frac{\log n}{\varepsilon n}\right) \int_{\varepsilon}^{\frac{t_*}{n^2}} \d r \sum_{\substack{k \in \Z^3 \\ \forall i\,:\, \lvert k_i\rvert \leq C_d \log n}} q_{r}\left(\frac{y}{n} + k \right) = \cO \left(\frac{\log^2 n}{n} \right) = o(1). \end{equation*} Therefore, \begin{equation}\label{eq:bound_on_0avgG_bdry} 2d \: n^{d-2} G^{0\avg}_\Lambda(0,y) \leq \frac{d}{2\pi^2} + o(1) + c_d \varepsilon^{2-\frac{d}{2}} + \int_{\varepsilon}^{\frac{t_*}{n^2}} \d r \left(\sum_{k \in \Z^3} q_{r}\left(\frac{y^*}{n} + k \right) - 1\right), \end{equation} where $y^*_1 = n/2$ and the remaining coordinates are the same as of $y$. We next show that $\sum_{k \in \Z^3} q_{r}(\frac{y^*}{n} + k )$ under our previous assumptions on $y^*$ is maximized at $(\frac{n}{2}, 0, \dots, 0)$. For this notice that $q_r(x) = \prod_{i=1}^d f_r(x_i)$ with $f_r$ being\linebreak a~density function of a centered normally distributed random variable of variance~$r/d$. Hence, \begin{align*} \sum_{k \in \Z^3} q_{r}\left(\frac{y^*}{n} + k \right) = \prod_{i=1}^d \left(\sum_{l \in \Z} f_{r}\left(\frac{y^*_i}{n} + l \right)\right) \eqqcolon \prod_{i=1}^d F_r\left(\frac{y^*_i}{n}\right). \end{align*} Note that $F_r$ is a function on a one-dimensional continuous torus $\T = [-1/2,1/2]/\sim$ (or alternatively, a periodic function on $\R$). Since furthermore it is in $L^2(\T)$, by Carleson's theorem it coincides with its Fourier series almost everywhere, and thus by continuity, everywhere. We have \begin{align*} F_r\left(\frac{y^*_i}{n}\right) = \sum_{p \in \Z} e^{- \frac{r}{2d} (2\pi)^2 p^2} \cos\left(2\pi p \frac{y^*_i}{n}\right) \leq \sum_{p \in \Z} e^{- \frac{r}{2d} (2\pi)^2 p^2} \end{align*} with equality reached at $y^*_i = 0$. So, from now on let $y^* = (\frac{n}{2}, 0, \dots, 0)$ (independently of the value of $y \in \hat{U}$). Let us further investigate the last summand of~\eqref{eq:bound_on_0avgG_bdry}. Using the aforementioned Fourier series and noticing that all of them are absolutely convergent and integrable as functions of $r$ on the interval $(\varepsilon, \infty)$, we get \begin{align*} I \,\MK &\!\Mk:= \int_{\varepsilon}^{\frac{t_*}{n^2}} \d r \biggl(\sum_{k \in \Z^3} q_{r}\Bigl(\frac{y^*}{n} + k \Bigr) - 1\biggr) \\ &= \int_{\varepsilon}^{\frac{t_*}{n^2}} \hspace{-0.1cm} \d r \Biggl[\sum_{w \in \Z^2 \setminus\{0\}} \hspace{-0.25cm} e^{-\frac{r}{2d}(2\pi)^2 \lvert w\rvert^2} \biggl(\underbrace{1 - 2 \sum_{l=1}^\infty (-1)^{l-1} e^{-\frac{r}{2d} (2\pi)^2 l^2}}_{\in [0,1] \text{ unif. in } r > \varepsilon} \biggr)\\ &\mkern300mu - 2 \sum_{l=1}^\infty (-1)^{l-1} e^{-\frac{r}{2d} (2\pi)^2 l^2} \Biggr]\\ &=: \int_{\varepsilon}^{\frac{t_*}{n^2}} h(r) \d r. \end{align*} Observe that $\int_{t_*/n^2}^\infty \sum_{l=1}^\infty e^{-\frac{r}{2d} (2\pi)^2 l^2} \d r = \cO(n^{-1})$. This implies that \begin{align*} I \leq \int_{\varepsilon}^{\infty} h(r) \d r + \cO(n^{-1}), \end{align*} where the integral on the right-hand side by Fubini's theorem, in turn, equals to \begin{multline*} \sum_{w \in \Z^2 \setminus\{0\}} \int_{\varepsilon}^{\infty} \d r e^{-\frac{r}{2d}(2\pi)^2 \lvert w\rvert^2} \biggl(1 - 2 \sum_{l=1}^\infty (-1)^{l-1} e^{-\frac{r}{2d} (2\pi)^2 l^2}\biggr) \\- 2 \sum_{l=1}^\infty (-1)^{l-1} \frac{2d}{(2\pi)^2 l^2} e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}. \end{multline*} Note furthermore that for a fixed constant $M > 0$ (to be determined, but independent of $\varepsilon$), \begingroup \allowdisplaybreaks \begin{multline*} \sum_{\substack{w \in \Z^2 \setminus\{0\}\\ \lvert w\rvert > M/\sqrt{\varepsilon}}} \int_{\varepsilon}^{\infty} \d r e^{-\frac{r}{2d}(2\pi)^2 \lvert w\rvert^2} \biggl(\underbrace{1 - 2 \sum_{l=1}^\infty (-1)^{l-1} e^{-\frac{r}{2d} (2\pi)^2 l^2}}_{\leq 1}\biggr)\\ \begin{aligned} &\leq \frac{2d}{(2\pi)^2} \sum_{\substack{w \in \Z^2 \setminus\{0\}\\ \lvert w\rvert > M/\sqrt{\varepsilon}}} \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 \lvert w\rvert^2}}{\lvert w\rvert^2} \\ &\leq \frac{2d}{(2\pi)^2} \left(2\pi \int_{\frac{M}{\sqrt{\varepsilon}}}^\infty \d s \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 s^2}}{s} + 4 \int_{\frac{M}{\sqrt{\varepsilon}}}^\infty \d s \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 s^2}}{s^2} \right)\\ &\leq \frac{(2d)^2 e^{-\frac{M^2}{2d} (2\pi)^2}}{(2\pi)^4 M^2} \left(\pi + 2\frac{\sqrt{\varepsilon}}{M}\right). \end{aligned} \end{multline*} \endgroup Altogether, \begin{multline*} I \leq \cO\left(\frac{1}{n}\right) + \frac{(2d)^2 e^{-\frac{M^2}{2d} (2\pi)^2}}{(2\pi)^4 M^2} \left(\pi + 2\frac{\sqrt{\varepsilon}}{M}\right) \\ \begin{aligned} &+ \frac{2d}{(2\pi)^2} \sum_{\substack{w \in \Z^2 \setminus\{0\}\\ \forall i: \lvert w_i\rvert \leq M/\sqrt{\varepsilon}}} e^{-\frac{\varepsilon}{2d} (2\pi)^2 \lvert w\rvert^2} \left[\frac{1}{\lvert w\rvert^2} - 2 \sum_{\substack{l = 1\\ \odd}}^{M/\sqrt{\varepsilon}} \left(\frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{\lvert w\rvert^2 + l^2} - \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}}{\lvert w\rvert^2 + (l+1)^2}\right) \right]\\ &- 2 \frac{2d}{(2\pi)^2}\sum_{\substack{l = 1\\ \odd}}^{\infty} \left(\frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{l^2} - \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}}{(l+1)^2}\right). \end{aligned} \end{multline*} We can write \begin{align*} \frac{1}{\lvert w\rvert^2} &= 2 \frac{2 + \pi \lvert w\rvert \tanh \bigl(\frac{\pi \lvert w\rvert}{2}\bigr) - \pi \lvert w\rvert \coth \bigl(\frac{\pi \lvert w\rvert}{2}\bigr)}{4 \lvert w\rvert^2} + \pi \frac{\coth \bigl(\frac{\pi \lvert w\rvert}{2}\bigr)- \tanh \bigl(\frac{\pi \lvert w\rvert}{2}\bigr)}{2 \lvert w\rvert}\\ &= 2\sum_{k=1}^{\infty} \frac{4k - 1}{(\lvert w\rvert^2 + (2k-1)^2)(\lvert w\rvert^2 + 4k^2)} + \pi \frac{1}{\lvert w\rvert \sinh (\pi \lvert w\rvert)}\\ &= 2\sum_{\substack{l=1\\ \odd}}^{\infty} \biggl(\frac{1}{\lvert w\rvert^2 + l^2} - \frac{1}{\lvert w\rvert^2 + (l+1)^2}\biggr) + \pi \frac{1}{\lvert w\rvert \sinh (\pi \lvert w\rvert)} \end{align*} using the following series' representations {(cf.~\cite[1.217 1, 1.421 2, 1.421 4]{Series_book} for the formulas and~\cite{Dun09} for their derivation)} for the hyperbolic trigonometric functions \begin{align*} \tanh(x) &= 8x \sum_{k=1}^\infty \frac{1}{\pi^2 (2k-1)^2+ 4x^2} \quad \forall \; x \in \R;\\ \coth(x) &= 2x \sum_{k=1}^\infty \frac{1}{\pi^2 k^2+x^2} + \frac{1}{x} \quad \forall \; x \in \R\setminus \{0\}. \end{align*} Thus far, we have shown that \begingroup \allowdisplaybreaks \begin{multline*} I \leq C_M \sqrt{\varepsilon} + \frac{(2d)^2 e^{-\frac{M^2}{2d} (2\pi)^2}}{(2\pi)^4 M^2} \pi\\ \begin{aligned} &+ \frac{2d}{(2\pi)^2}\sum_{\substack{w \in \Z^2 \setminus\{0\}\\ \lvert w\rvert \leq M/\sqrt{\varepsilon}}} e^{-\frac{\varepsilon}{2d} (2\pi)^2 \lvert w\rvert^2} 2 \left[\sum_{\substack{l = 1\\ \odd}}^{M/ \sqrt{\varepsilon}} \biggl(\frac{1 -e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{\lvert w\rvert^2 + l^2} - \frac{1 - e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}}{\lvert w\rvert^2 + (l+1)^2}\biggr)\right. \\ &\mkern210mu\left. + \sum_{\substack{l=1 + M/\sqrt{\varepsilon}\\ \odd}}^{\infty} \biggl(\frac{1}{\lvert w\rvert^2 + l^2} - \frac{1}{\lvert w\rvert^2 + (l+1)^2}\biggr) \right]\eqqcolon A_1 \\ &+ \frac{2d}{(2\pi)^2} \sum_{\substack{w \in \Z^2 \setminus\{0\}\\ \lvert w\rvert \leq M/\sqrt{\varepsilon}}} e^{-\frac{\varepsilon}{2d} (2\pi)^2 \lvert w\rvert^2} \pi \frac{1}{\lvert w\rvert \sinh (\pi \lvert w\rvert)}\eqqcolon A_2\\ &- 2 \frac{2d}{(2\pi)^2}\sum_{\substack{l = 1\\ \odd}}^{\infty} \biggl(\frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{l^2} - \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}}{(l+1)^2}\biggr) \eqqcolon A_3. \end{aligned} \end{multline*} \endgroup for an appropriate $C_M > 0$. For the term $A_3$, notice that by Fatou's lemma \begin{align*} \liminf_{\varepsilon \to 0} \sum_{\substack{l = 1\\ \odd}}^{\infty} \left(\frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{l^2} - \frac{e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}}{(l+1)^2}\right) \geq \sum_{\substack{l = 1\\ \odd}}^{\infty} \left(\frac{1}{l^2} - \frac{1}{(l+1)^2}\right) = \frac{\pi^2}{12}. \end{align*} Therefore, for all $\varepsilon > 0$ very small, there exists $c(\varepsilon) > 0$ which converges to $0$ as $\varepsilon$ does, such that $A_3 \leq -\frac{2d}{(2\pi)^2} \frac{\pi^2}{6} + c(\varepsilon)$. As for $A_2$, \begin{align*} \frac{(2\pi)^2}{2d} A_2 &\leq \sum_{w \in \Z^2 \setminus\{0\}} \frac{\pi}{\lvert w\rvert \sinh (\pi \lvert w\rvert)} \leq 2\pi\int_{2}^\infty \d s \frac{\pi s}{s \sinh (\pi s)} + \sum_{k=1}^\infty \frac{4\pi}{k \sinh (\pi k)} \\ &\quad+ \frac{4\pi}{\lvert (1,1)\rvert \sinh (\pi \lvert (1,1)\rvert)} + \frac{4\pi}{\lvert (2,2)\rvert \sinh (\pi \lvert (2,2)\rvert)} + \frac{8\pi}{\lvert (1,2)\rvert \sinh (\pi \lvert (1,2)\rvert)}\\ &\leq -2\pi \log(\tanh(\pi)) + \frac{4\pi}{\sinh(\pi)} + \int_1^\infty \d s\frac{2\pi}{\sinh (\pi s)} \\ &\quad+ 4\pi \left(\frac{1}{\sqrt{2} \sinh (\pi \sqrt{2})} + \frac{1}{2\sqrt{2} \sinh (\pi 2\sqrt{2})} + \frac{2}{\sqrt{5} \sinh (\pi \sqrt{5})}\right)\\ &= -2\pi \log(\tanh(\pi)) - 2\log(\tanh(\pi/2))\\ &\quad+ 4\pi \left(\frac{1}{\sinh(\pi)} + \frac{1}{\sqrt{2} \sinh (\pi \sqrt{2})} + \frac{1}{2\sqrt{2} \sinh (\pi 2\sqrt{2})} + \frac{2}{\sqrt{5} \sinh (\pi \sqrt{5})}\right). \end{align*} We have $\frac{(2\pi)^2}{2d} (A_2 + A_3) \leq c'(\varepsilon) + 2\pi(0.2411 - \pi/12) \leq c'(\varepsilon) - 0.04\pi$. It remains to estimate $A_1$. We begin by analyzing $l \mapsto \frac{1 -e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{\lvert w\rvert^2 + l^2} \eqqcolon g(l)$ (of interest to us are its regions of monotonicity). For $l \geq 1$, \begin{align*} \frac{(\lvert w\rvert^2 + l^2)}{2l} g'(l) &= e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2} \left(1 + \frac{\varepsilon}{2d} (2\pi)^2(l^2 + \lvert w\rvert^2) \right) - 1\\ &\geq \left(1- \frac{\varepsilon}{2d} (2\pi)^2 l^2\right) \left(1 + \frac{\varepsilon}{2d} (2\pi)^2(l^2 + \lvert w\rvert^2) \right) - 1\\ &= -\left(\frac{\varepsilon (2\pi)^2}{2d} \right)^2 \left(l^4 + l^2 \lvert w\rvert^2 - \frac{2d}{(2\pi)^2\varepsilon} \lvert w\rvert^2 \right). \end{align*} In particular, the subinterval of $[1, \infty)$ on which $g(l)$ is decreasing ($\{l \geq 1: g'(l) < 0\}$) is therefore included in the range of solutions of the following inequality \begin{align*} l^4 + l^2 \lvert w\rvert^2 &\geq \frac{2d}{(2\pi)^2\varepsilon} \lvert w\rvert^2\\ &\Updownarrow\\ l^2 + \frac{\lvert w\rvert^2}{2} &\geq \lvert w\rvert^2 \sqrt{ \frac{2d}{(2\pi)^2\varepsilon} \frac{1}{\lvert w\rvert^2} + \frac{1}{4}}\\ &\Updownarrow l\geq 1\\ l &\geq \lvert w\rvert \sqrt{\sqrt{ \frac{2d}{(2\pi)^2\varepsilon} \frac{1}{\lvert w\rvert^2} + \frac{1}{4}} - \frac{1}{2}} \eqqcolon w_*. \end{align*} Furthermore, observe that \begin{align*} g(l) - g(l+1) &= \frac{(2l+1)\left(1- e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}\right) - (\lvert w\rvert^2 + l^2) \left(e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2} - e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}\right)}{(\lvert w\rvert^2 + l^2)(\lvert w\rvert^2 + (l+1)^2)}\\ &\leq \frac{2l+1}{(\lvert w\rvert^2 + l^2)(\lvert w\rvert^2 + (l+1)^2)}. \end{align*} Hence, since on $\{l \geq 1: g'(l) \geq 0\}, \quad g(l) - g(l+1) \leq 0$ (up to a borderline value of~$l$), \begin{multline*} 2\mkern-8mu\sum_{\substack{l = 1\\ \odd}}^{M/ \sqrt{\varepsilon}} \biggl(\frac{1 -e^{-\frac{\varepsilon}{2d} (2\pi)^2 l^2}}{\lvert w\rvert^2 + l^2} - \frac{1 - e^{-\frac{\varepsilon}{2d} (2\pi)^2 (l+1)^2}}{\lvert w\rvert^2 + (l+1)^2}\biggr) + 2\mkern-15mu\sum_{\substack{l=1 + M/\sqrt{\varepsilon}\\ \odd}}^{\infty} \biggl(\frac{1}{\lvert w\rvert^2 + l^2} - \frac{1}{\lvert w\rvert^2 + (l+1)^2}\biggr)\\ \begin{aligned} &\leq 2\sum_{\substack{l=w_* - 1\\ \odd}}^{\infty} \frac{2l+1}{(\lvert w\rvert^2 + l^2)(\lvert w\rvert^2 + (l+1)^2)} \\ &\leq \int_{w_* - 1}^\infty \frac{(2x+1) \d x}{(\lvert w\rvert^2 + x^2)(\lvert w\rvert^2 + (x+1)^2)} + \frac{C}{\lvert w\rvert^3} \ind_{\bigl\{\lvert w\rvert \geq \frac{c}{\sqrt{\varepsilon}}\bigr\}}\\ &= \frac{1}{\lvert w\rvert} \biggl[\arctan \Bigl(\frac{w_*}{\lvert w\rvert}\Bigr) - \arctan \Bigl(\frac{w_*-1}{\lvert w\rvert}\Bigr)\biggr] + \frac{C}{\lvert w\rvert^3} \ind_{\bigl\{\lvert w\rvert \geq \frac{c}{\sqrt{\varepsilon}}\bigr\}}\\ &\leq \frac{1}{\lvert w\rvert^2 + (w_*-1)^2} + \frac{C}{\lvert w\rvert^3} \ind_{\bigl\{\lvert w\rvert \geq \frac{c}{\sqrt{\varepsilon}}\bigr\}} \end{aligned} \end{multline*} for appropriate constants $c, C>0$. The second summand in the third line comes from the evaluation of the integrand at its local maximum point (which clearly has the leading term proportional to $\lvert w\rvert$\footnote{An explicit computation shows that the pre-coefficient of the leading term of the local maximum point is $1/\sqrt{3}$.}) if it is greater or equal to $w_*-1$. Note further that $w_* - 1 \geq w_* (1 - \sqrt{c} \varepsilon^{1/4})$ for some $c>0$ and \begin{align*} \lvert w\rvert^2 + w_*^2 \left(1 - \sqrt{c} \varepsilon^{1/4}\right)^2 &\geq \lvert w\rvert^2 \left(1 + (1 - c \sqrt{\varepsilon}) \left(\sqrt{\frac{1}{4}+ \frac{2d}{(2\pi)^2\varepsilon} \frac{1}{\lvert w\rvert^2}} - \frac{1}{2}\right) \right) \\ &\geq \frac{1}{2} \left(\lvert w\rvert^2 + 2 \frac{1-c\sqrt{\varepsilon}}{1+c\sqrt{\varepsilon}} \sqrt{\frac{2d}{(2\pi)^2\varepsilon}} \lvert w\rvert \right). \end{align*} For simplicity set $a \coloneqq 2 \frac{1-c\sqrt{\varepsilon}}{1+c\sqrt{\varepsilon}} \sqrt{\frac{2d}{(2\pi)^2\varepsilon}} \approx 2 \sqrt{\frac{2d}{(2\pi)^2\varepsilon}}$. It is now possible to estimate $A_1$, \begingroup \allowdisplaybreaks \begin{align*} \frac{(2\pi)^2}{2d} A_1 &\leq \sum_{\substack{w \in \Z^2 \setminus\{0\}\\ \lvert w\rvert \leq M/\sqrt{\varepsilon}}} \left[\frac{2}{\lvert w\rvert^2 + a \lvert w\rvert} + \frac{C}{\lvert w\rvert^3} \ind_{\left\{\lvert w\rvert \geq \frac{c}{\sqrt{\varepsilon}}\right\}}\right] \\ &\leq 2\pi \int_{1}^{\frac{M}{\sqrt{\varepsilon}}} \frac{2r \d r}{r^2 + a r} + 4 \int_{1}^{\infty} \frac{2\d r}{r^2 + a r} + C' \int_{\frac{c}{\sqrt{\varepsilon}}}^\infty \frac{r\d r}{r^3} + \frac{8}{1 + a}\\ &\leq C\sqrt{\varepsilon} + 4\pi \log \left(\frac{a + M/\sqrt{\varepsilon}}{1 + a}\right) + 8 \frac{\log(a+1)}{a}\\* &\leq 4\pi \log \left(1 + \frac{2\pi M}{2 \sqrt{2d}}\right) + c(\varepsilon) \end{align*} \endgroup with $c(\varepsilon) \to 0$ as $\varepsilon \to 0$. To sum up, for all sufficiently large $n$, \begin{align*} n^{d-2} &G^{0\avg}_{\Lambda}(0,y) \\ &\leq \frac{1 - 0.04\pi + 4\pi \log \left(1 + \frac{2\pi M}{2 \sqrt{2d}}\right) + \pi \frac{2d}{(2\pi)^2 M^2} e^{-\frac{M^2}{2d} (2\pi)^2}}{(2\pi)^2} + c(\varepsilon) + c_d \varepsilon^{2-\frac{d}{2}} \end{align*} with $c(\varepsilon)>0$ that can be made arbitrarily small by choosing $\varepsilon > 0$ sufficiently small (but non-zero). Note that since $d= 3$, the last summand is as well arbitrary small. Choose $M> 0$ such that $\frac{2\pi M}{\sqrt{2d}} = 1$. Then, \begin{equation}\label{eq:final_bound_G0avg} n G^{0\avg}_{\Lambda}(0,y) \leq \frac{1 - 0.04\pi + 4\pi \log(3/2) + \pi e^{-1}}{(2\pi)^2} + c(\varepsilon) \eqqcolon c_* < \frac{3^{2/3}}{2 (4\pi)^{2/3}}.\qedhere \end{equation} \end{proof} \appendix \section{}\label{sec:Appendix} \subsection{Multivariate local central limit theorem: Proof of Proposition~\ref{prop:dens_conv}} \label{A:proof_densconv} In this subsection, we prove a local central limit theorem for the densities of i.i.d. random vectors used in the study of the infinite spin-dimensionality limit of the spin $O(N)$ model. \begin{prop} Let $(X^i)_i$ be a sequence of i.i.d. square-integrable, centered random vectors with probability density function (with respect to Lebesgue measure on $\R^k$) $f: \R^k \to \R$. Assume that \begin{enumerate} \item\label{propA.1.1} $f \in L^r(\R^k)$ for some $r \in (1,2]$; \item\label{propA.1.2} The covariance matrix $C$ of $X^1 = (X^1_j)_{j=1}^k$ is positive definite. \end{enumerate} Then the relation \begin{align*} n^{k/2} f^{(n)} (\sqrt{n} v) \xrightarrow{n\to \infty} \frac{1}{(2\pi)^{k/2}\sqrt{\det C}} \: e^{-\frac{1}{2} (v, C^{-1}v)} \end{align*} holds uniformly over $v \in \R^k$. Here, $f^{(n)}$ is the $n$-fold convolution of $f$, as well as the density function of $\sum_{i=1}^n {X}^i$. \end{prop} The following proof follows the lines of the one for the case $k=1$ in~\cite[Chapter~8, Section~46, Theorem~1]{gnedenko1968limit}. Note that, deviating from the notation used in the main body of the paper, we write $(x,t)$ to denote the canonical inner product on $\R^k$. \begin{proof} Let $\psi$ denote the characteristic function of $X^1$, and $\psi_n$ the characteristic function of $\sum_{i=1}^n X^i$. In particular, \begin{align*} \psi({t}) = \int_{\R^k} e^{i (x,t)} f({x}) \d{x} \quad \text{and} \quad \psi_n({t}) = \int_{\R^k} e^{i (x,t)} f^{(n)}({x}) \d{x} = (\psi({t}))^n. \end{align*} By assumption, $f$ belongs to $L^r(\R^k)$ and to $L^1(\R^k)$, as it is a density function. Therefore, also $f \in L^u(\R^k)$ for any $u \in (1,r)$. This follows from H\"older's inequality, as we now explain. Let $\theta = \frac{r}{r-1} \frac{u-1}{u} \in (0,1)$, $p = \frac{r-1}{u-1} \in (1,\infty)$ and $p'= p/(p-1)$ such that $\frac{1}{p} + \frac{1}{p'} = 1$. Then, by H\"older's inequality, \begin{align*} \int_{\R^k} \lvert f\rvert^{u} \d {x} \leq \ntwo{ \lvert f\rvert^{u\theta}}_{L^p} \ntwo{\lvert f\rvert^{u(1-\theta)}}_{L^{p'}} = \ntwo{f}_{L^r}^{r/p} \ntwo{f}_{L^1}^{1/{p'}} <\infty. \end{align*} Let $\hat{f}({t}) \coloneqq \int_{\R^k} e^{-2\pi i(x,t)} f({x}) \d{x}$ be the Fourier transform of $f$. Then, by Hausdorff--Young inequality (see, e.g., \cite[Proposition~2.2.16]{Grafakos_2014}) we know that if $f \in L^{p}(\R^k)$ for some $p \in [1,2]$, then $\hat{f} \in L^{p'}(\R^k)$ for $p'=\frac{p}{p-1}$. Thus, the above observations about $f$ and the fact that $\hat{f}({t}) = \psi(-2\pi {t})$ yield that $\psi$ is in $L^{q}(\R^k)$ for any $q \geq \frac{r}{r-1}$. In particular, we get that for any $n \geq \frac{r}{r-1}$, $\psi_n$ is integrable. Hence, by Fourier inversion, \begin{align*} f^{(n)}({x}) &= \frac{1}{(2\pi)^k} \int_{\R^k} e^{- i (x,t)} \psi_n({t}) \d{t}, \intertext{and thus,} n^{k/2} f^{(n)}(\sqrt{n}{x}) &= \frac{1}{(2\pi)^k} \int_{\R^k} e^{- i (x,t)} \psi_n({t}/\sqrt{n}) \d{t}. \end{align*} On the other hand, we know that \begin{align*} \frac{1}{(2\pi)^{k/2}\sqrt{\det(C)}} \: e^{-\frac{1}{2} ({x}, C^{-1}{x})} = \frac{1}{(2\pi)^k} \int_{\R^k} e^{- i (x,t) - \frac{1}{2} ({t}, C{t})} \d{t}. \end{align*} Therefore, to conclude the proof it suffices to show that \begin{align*} R_n({x}) \coloneqq \int_{\R^k} e^{- i (x,t)} \left(\psi_n({t}/\sqrt{n}) - e^{- \frac{1}{2} ({t}, C{t})}\right)\d{t} \end{align*} converges to zero as $n$ tends to infinity uniformly over ${x} \in \R^k$. By central limit theorem (multivariate version) and dominated convergence theorem~($\norm{\psi_n}_{L^\infty} \leq 1$), for any $r> 0$ fixed, \begin{align*} \left\lvert \int_{B^k(0,r)} e^{- i (x,t)} \left(\psi_n({t}/\sqrt{n}) - e^{- \frac{1}{2} ({t}, C{t})}\right)\d{t} \right\rvert \xrightarrow{n\to \infty} 0 \quad\text{uniformly over } \:{x}. \end{align*} Moreover, by taking $r$ sufficiently large we can make $\left\lvert \int_{\lvert {t}\rvert > r} e^{- i (x,t)- \frac{1}{2} ({t}, C{t})}\d{t} \right\rvert$ arbitrarily small. It remains to treat the term \begin{align*} I \coloneqq \int_{\lvert {t}\rvert > r} e^{- i (x,t)} \psi_n({t}/\sqrt{n}) \d{t}, \end{align*} which we further split into two integrals $I_1$ and $I_2$ over $r < \lvert {t}\rvert \leq \varepsilon \sqrt{n}$ and $\lvert {t}\rvert > \varepsilon\sqrt{n}$, respectively, for some $\varepsilon>0$ small. We start with $I_1$. By our assumptions, $\psi$ is twice differentiable with \begin{align*} \nabla\psi ({0}) = {0}, \quad D^2 \psi({0}) = -C. \end{align*} Hence, by Taylor's theorem, in the neighbourhood of ${0}$, \begin{align*} \psi({t}) - 1 = (\nabla\psi ({0}))^T {t} + \frac{1}{2} {t}^T (D^2\psi({0})) {t} + o(\lvert {t}\rvert^2) = -\frac{1}{2} {t}^T C {t} + o(\lvert {t}\rvert^2). \end{align*} By assumption, $C$ is positive definite, so it is possible to choose $\varepsilon>0$ sufficiently small such that for all $\lvert{t}\rvert \leq \varepsilon$, the remainder of the expansion $o(\lvert {t}\rvert^2)$ is smaller than $\frac{1}{4} {t}^T C {t}$. Then, $\lvert \psi({t})\rvert \leq 1 - \frac{1}{4} {t}^T C {t} \leq e^{-\frac{1}{4} {t}^T C {t}}$ and \begin{align*} \lvert I_1\rvert \leq \int_{r < \lvert {t}\rvert \leq \varepsilon \sqrt{n}} \lvert\psi_n({t}/\sqrt{n}) \rvert\d{t} \leq \int_{r < \lvert {t}\rvert \leq \varepsilon \sqrt{n}} e^{-n\frac{1}{4n} {t}^T C {t}} \d{t} \leq \int_{r < \lvert {t}\rvert} e^{-\frac{1}{4} {t}^T C {t}} \d{t}. \end{align*} Thus, by taking $r$ sufficiently large we can make $\lvert I_1\rvert$ arbitrarily small. Finally, we treat $I_2$. Recall that $\hat{f}({t}) = \psi(-2\pi {t})$, and thus by Riemann-Lebesgue lemma (see~\cite[Proposition~2.2.17]{Grafakos_2014}), $\lvert \psi({t}) \rvert \to 0$ as $\lvert{t}\rvert \to \infty$. Moreover, $\lvert \psi({t})\rvert < 1$ for ${t} \neq{0}$. Indeed, suppose first that there exists a non-zero ${t}$ such that $\psi({t}) = 1$. Then the following must hold: \begin{align*} \int \underbrace{[1- \cos({(x,t)})]}_{\geq 0} f({x}) \d{x} = \int \underbrace{[1-\sin({(x,t)})]}_{\geq 0} f({x}) \d{x} = 0. \end{align*} And hence, $1- \cos({(x,t)}) = 0$ for $f({x}) \d{x}$-a.e. ${x}$, which is clearly impossible unless ${t}$ is zero. Suppose now that $\lvert\psi({t})\rvert = 1$, and consider $Z \coloneqq X^1 - X^2$. Then $\psi_Z({t}) = 1$, and we get that for $f_Z({x}) \d{x}$-a.e. ${x}$, $1- \cos({(x,t)}) = 0$. This again implies that ${t} = 0$. Altogether, we proved that there exists $c>0$ such that $\lvert\psi({t})\rvert < e^{-c}$ for all $\lvert {t}\rvert > \varepsilon$. Thus, for a fixed $q > \frac{r}{r-1}$, \[ \lvert I_2\rvert \leq e^{-c(n-q)} \int_{\lvert {t}\rvert > \varepsilon \sqrt{n}} \lvert\psi({t}/\sqrt{n}) \rvert^q \d{t} = \underbrace{e^{-c(n-q)} n^{k/2}}_{\xrightarrow{n\to\infty} 0} \underbrace{\int_{\lvert {t}\rvert > \varepsilon} \lvert\psi({t}) \rvert^q \d{t}}_{<\infty}.\qedhere \] \end{proof} \subsection{Local CLT for triangular arrays of independent random variables: Proof of Proposition~\ref{prop:dens_conv_1D}} \label{A:proof_densconv_1D} This subsection establishes a local CLT for the normalized sum of independent, non-identically distributed random variables. \begin{prop}\label{propA.2} Let $(X_{i,n})_{i\leq n, n \in \N}$ be a triangular array of independent centered random variables with the probability density functions $f_{i,n}$ (w.r.t. Lebesgue measure on $\R$). Assume that \begin{enumerate} \item\label{propA.2.1} $f_{i,n} \in L^r(\R)$ for some $r \in (1,2]$ (independent of $i,n$) such that $\sup_{i, n} \norm{f_{i,n}}_{L^r} \leq M$ for some $M> 0$; \item\label{propA.2.2} For all $i \leq n, n \in \N$, $\sigma_{i,n}^2 \coloneqq \Var[X_{i,n}] < \infty$ ordered such that $\sigma_{1,n}^2 \geq \dots \geq \sigma_{n,n}^2$; \item\label{propA.2.3} Lindeberg's condition is satisfied, that is, for any $\varepsilon > 0$, \begin{align*} \frac{\sum_{i \geq 1} \E\left[X^2_{i,n} \ind_{\{ \lvert X_{i,n}\rvert > \varepsilon s_n\}}\right]}{s_n^2} \xrightarrow{n \to \infty} 0, \end{align*} where $s_n^2 \coloneqq \sum_{i=1}^n \sigma_{i,n}^2$. \item\label{propA.2.4} There exist $\delta > 0$ (uniform), $K(n) \geq 1$, $l_*(n) \geq 1, n \geq l^*(n) \geq 2 \lceil\frac{r}{r-1}\rceil$ such that for all $n$ sufficiently large \begin{gather} \tag{$a$} \frac{\sum_{i \geq l_*(n)} \sigma_{i,n}^2}{\sum_{i=1}^n \sigma_{i,n}^2} \geq \delta;\label{propA.2.4.a}\\ \tag{$b$} \frac{\sum_{i \geq l_*(n)} \E\left[X^2_{i,n} \ind_{\{ \lvert X_{i,n}\rvert > K(n)\}}\right]}{\sum_{i \geq l_*(n)} \sigma_{i,n}^2} \leq \frac{1}{8}; \label{propA.2.4.b}\\ \tag{$c$} \frac{n - l^*(n)}{\sigma_{l^*, n}^2 \vee K(n)^2} \gg \log (s_n^2).\label{propA.2.4.c} \end{gather} \end{enumerate} Then the relation \begin{align*} s_n f^{(n)} (s_n x) \xrightarrow{n\to \infty} \frac{1}{\sqrt{2\pi}} \: e^{-\frac{x^2}{2}} \end{align*} holds uniformly over $\R$. Here, $f^{(n)}$ is the convolution of $f_{1,n}, \dots, f_{n,n}$, as well as the density function of $\sum_{i=1}^n X_{i,n}$. Note that $s_n f^{(n)}(s_n x)$ is the density of the normalized sum $\frac{1}{s_n} \sum_{i=1}^n X_{i,n}$. \end{prop} \begin{proof} Let $\psi_{i,n}$ denote the characteristic function of $X_{i,n}$, and $\psi^{(n)}$ the characteristic function of $\sum_{i=1}^n X_{i,n}$, $\hat{g}(t) \coloneqq \int_{\R} e^{-2\pi i x t} g(x) \d x$ be the Fourier transform of $g: \R \to \R$ (for an appropriate $g$). By the same argument as in Section~\ref{A:proof_densconv}, $\hat f_{i,n}, \psi_{i,n} \in L^{q}(\R)$ for any $q \geq \frac{r}{r-1}$. In particular, by H\"older's inequality, for any fixed $q \geq \frac{r}{r-1}$ and $n \geq q$, \begin{align*} \ntwo{\psi^{(n)}}_{L^1} = \left\lVert\prod_{i=1}^n \psi_{i,n}\right\rVert_{L^1} \leq \prod_{i=1}^n \ntwo{\psi_{i,n}}_{L^n} < \infty. \end{align*} Hence, by inversion formula theorem, $f^{(n)}(x) = \frac{1}{2\pi} \int_{\R} e^{- i x t} \psi^{(n)}(t) \d t$. To conclude the proof, it suffices to show that \begin{align*} D_n({x}) \coloneqq \int_{\R} e^{- i x t} \left(\psi^{(n)} (t/s_n) - e^{- \frac{t^2}{2}}\right)\d t \end{align*} converges to zero as $n$ tends to infinity uniformly over $x \in \R$. We first split $D_n$ into two integrals over $\lvert t \rvert \leq r$ and $\rvert t \rvert > r$ for some $r>0$ (to be specified later). By Lindeberg--Feller central limit theorem and dominated convergence theorem~($\norm{\psi_n}_{L^\infty} \leq 1$), the former converges to zero as $n$ tends to infinity uniformly in $x \in \R$. By taking $r$ sufficiently large we can also make $\lvert \int_{\lvert t\rvert > r} e^{- i x t- \frac{t^2}{2}}\d {t} \rvert$ arbitrarily small. The remaining part $\int_{\lvert t\rvert > r} e^{- i xt} \psi^{(n)} (t/s_n) \d t$ we split into two integrals $I_1$ and $I_2$ over $r < \lvert t\rvert \leq \delta_n s_n$ and $\lvert t\rvert > \delta_n s_n$, respectively, for some $\delta_n >0$ small (to be chosen later in dependence of $K(n)$). By our assumptions, $\psi_{i,n}' (0) = 0$ and $\psi_{i,n}''(0) = - \sigma_{i,n}^2$ for all $i,n$. Hence, by Taylor's theorem, in the neighborhood of $0$ ($\lvert t \rvert \leq \delta_n$), \begin{align*} \psi_{i,n}(t) = 1 -\frac{\sigma_{i,n}^2 t^2}{2} + R_{i,n}(t). \end{align*} Since $e^{ix} = 1 + ix -\frac{x^2}{2} + R(x)$ with $\lvert R(x)\rvert \leq \min(\lvert x\rvert^3/6, x^2)$, the remainder term $R_{i,n}(t)$ in the above expansion of $\psi_{i,n}$ satisfies \begin{align*} \lvert R_{i,n}(t)\rvert &\leq \frac{1}{6} \E\left[\lvert t X_{i,n}\rvert^3 \ind_{\{\lvert X_{i,n}\rvert \leq K\}}\right] + \E\left[\lvert t X_{i,n}\rvert^2 \ind_{\{\lvert X_{i,n}\rvert \geq K\}}\right] \\ &\leq \frac{K \delta_n}{6} t^2 \sigma_{i,n}^2 + t^2 \E\left[\lvert X_{i,n}\rvert^2 \ind_{\{\lvert X_{i,n}\rvert \geq K\}}\right] \end{align*} for any $K > 0$. In particular, also for $K(n)$ as in the fourth assumption. By choosing $\delta_n = \frac{3}{4 K(n)} > 0$, we obtain that $\lvert \psi_{i,n}(t)\rvert \leq 1 - \frac{3}{8} \sigma_{i,n}^2 t^2 + t^2 \E[\lvert X_{i,n}\rvert^2 \ind_{\{\lvert X_{i,n}\rvert \geq K(n)\}}]\leq \exp(-\frac{3}{8} \sigma_{i,n}^2 t^2 + t^2 \E[\lvert X_{i,n}\rvert^2 \ind_{\{\lvert X_{i,n}\rvert \geq K(n)\}}])$ for all $\lvert t\rvert \leq \delta_n$. Hence, \begin{align*} \lvert I_1\rvert &\leq \int_{r < \lvert t\rvert \leq \delta_n s_n} \lvert\psi^{(n)} (t/s_n) \rvert \d t \leq \int_{r < \lvert t\rvert \leq \delta_n s_n} \prod_{i=l_*(n)}^n e^{-\frac{3}{8} \frac{\sigma_{i,n}^2}{s_n^2} t^2 + t^2 \frac{\E\left[\lvert X_{i,n}\rvert^2 \ind_{\{\lvert X_{i,n}\rvert \geq K(n)\}}\right]}{s_n^2}} \d t\\ &\leq \int_{r < \lvert t\rvert \leq \delta_n s_n} \exp \Biggl(-t^2 \underbrace{\frac{\sum_{i \geq l_*(n)} \sigma_{i,n}^2}{s_n^2}}_{\geq \delta} \underbrace{\Biggl[\frac{3}{8} - \frac{\sum_{i \geq l_*(n)}\E\left[\lvert X_{i,n}\rvert^2 \ind_{\{\lvert X_{i,n}\rvert \geq K(n)\}}\right]}{\sum_{i \geq l_*(n)} \sigma_{i,n}^2} \Biggr]}_{\geq 1/4} \Biggr) \d t\\ &\leq \int_{r < \lvert t\rvert} \exp (- t^2 \delta/ 4) \d t \end{align*} By further increasing $r$, we can make $\lvert I_1\rvert$ arbitrarily small. We now turn to the remaining $I_2$. Let $\Mm_{i,n}$ be the median of $f_{i,n}(X_{i,n})$, i.e., $\P[f_{i,n}(X_{i,n}) \geq \Mm_{i,n}] \geq 1/2$ and $\P[f_{i,n}(X_{i,n}) \leq \Mm_{i,n}] \geq 1/2$. By~\cite[Theorem~2]{BCG12}, there exist two absolute constants $c_1, c_2 > 0$ such that \begin{align*} \lvert \psi_{i,n}(t)\rvert < \begin{cases} 1 - c_1/(\Mm_{i,n}^2 \sigma_{i,n}^2), & \lvert t\rvert \geq {\pi}/(4 \sigma_{i,n});\\ 1 - c_2 t^2/\Mm_{i,n}^2, & 0 < \lvert t\rvert < {\pi}/({4 \sigma_{i,n}}). \end{cases} \end{align*} Since $f_{i,n}(X_{i,n})$ is non-negative, Markov's inequality implies that for any $\lambda > 0$, \begin{align*} \P[f_{i,n}(X_{i,n}) \geq \lambda] \leq \frac{\E[f_{i,n}(X_{i,n})^{r-1}]}{\lambda^{r-1}} = \frac{\int f_{i,n}^{r}(x) \d x}{\lambda^{r-1}} \leq \frac{M^r}{\lambda^{r-1}}, \end{align*} which directly implies that $0< \Mm_{i,n} \leq (2 M^r)^{1/(r-1)} \eqqcolon C^{1/2}$. Thus, for all $i \leq n,\linebreak n \in \N$: \begin{align*} \lvert \psi_{i,n}(t)\rvert < \begin{cases} 1 - c_1/(C \sigma_{i,n}^2), & \lvert t\rvert \geq {\pi}/(4 \sigma_{i,n});\\ 1 - c_2 \delta_n^2/C, & \delta_n < \lvert t\rvert < {\pi}/({4 \sigma_{i,n}}); \end{cases} \end{align*} and $\sup_{i \geq l^*(n)} \lvert \psi_{i,n}(t)\rvert < 1 - \min(c_1/(C \sigma_{l^*,n}^2), c_2 \delta_n^2/C) \leq e^{-c' \min(\sigma_{l^*,n}^{-2}, \delta_n^2)}$ for some $c'>0$ and all $\lvert t \rvert \geq \delta_n$. Thus, for any fixed integer $2\lceil \frac{r}{r-1}\rceil > q > \frac{r}{r-1}$, \begingroup\allowdisplaybreaks \begin{align*} \lvert I_2\rvert &\leq s_n \int_{\lvert t\rvert > \delta_n} \prod_{i=1}^n \lvert \psi_{i,n}(t) \rvert \d t \leq e^{-(n-l^*) c' \min\left(\sigma_{l^*,n}^{-2}, \delta_n^2\right)} s_n \int_{\lvert t\rvert > \delta_n} \prod_{i=1}^q \lvert \psi_{i,n}(t) \rvert \d t \\ &\leq e^{-(n-l^*) c' \min\left(\sigma_{l^*,n}^{-2}, \delta_n^2\right)} s_n \prod_{i=1}^q \ntwo{\psi_{i,n}}_{L^q} \\&\leq e^{-(n-l^*) c' \min\left(\sigma_{l^*,n}^{-2}, \delta_n^2\right)} s_n \prod_{i=1}^q (2\pi)^{1/q} \ntwo{f_{i,n}}_{L^{q/(q-1)}} \\ &\leq e^{-(n-l^*) c' \min\left(\sigma_{l^*,n}^{-2}, \delta_n^2\right)} s_n 2 \pi \prod_{i=1}^q \left[\ntwo{f_{i,n}}_{L^r}^{r/(q(r-1))} \ntwo{f_{i,n}}_{L^1}^{1-r/(q(r-1))}\right]\\ &\leq 2 \pi e^{-(n-l^*) c' \min\left(\sigma_{l^*,n}^{-2}, \delta_n^2\right)} s_n M^{r/(r-1)} \xrightarrow[\text{unif. in }x]{n \to \infty} 0, \end{align*} \endgroup where we used H\"older's and Hausdorff--Young inequalities and the third assumption ($\delta_n = c/K(n)$). \end{proof} \bibliography{20260519_aru} \end{document}