\documentclass[fleqn,a4paper]{article}

\usepackage[utf8]{inputenc} % why not type "Bézout" with unicode?
\usepackage[T1]{fontenc} % vector fonts plz
\usepackage{fullpage,amsmath,amssymb,latexsym,graphicx}
\usepackage{newtxtext,newtxmath} % Times for PR
\usepackage{appendix}
\usepackage[dvipsnames]{xcolor}
\usepackage[
  colorlinks=true,
  urlcolor=MidnightBlue,
  citecolor=MidnightBlue,
  filecolor=MidnightBlue,
  linkcolor=MidnightBlue
]{hyperref} % ref and cite links with pretty colors
\usepackage[
  style=phys,
  eprint=true,
  maxnames = 100
]{biblatex}
\usepackage{anyfontsize,authblk}

\usepackage{tikz}

\addbibresource{marginal.bib}

\begin{document}

\title{
  None yet
}

\author{Jaron Kent-Dobias}
\affil{Istituto Nazionale di Fisica Nucleare, Sezione di Roma I}

%\maketitle
%\begin{abstract}
%\end{abstract}

\section{Introduction}

Systems with rugged landscapes are important across many disciplines, from the
physics to glasses and spin-glasses to the statistical inference problems. The
behavior of these systems is best understood when equilibrium or optimal
solutions are studied and averages can be taken statically over all possible
configurations. However, such systems are also infamous for their tendency to
defy equilibrium and optimal expectations in practice, due to the presence of
dynamic transitions or crossovers that leave physical or algorithmic dynamics
stuck exploring only a subset of configurations.

In some simple models of such landscapes, it was recently found that marginal
minima are significant as the attractors of gradient descent dynamics
\cite{Folena_2020_Rethinking, Folena_2023_On}. This extends to more novel
algorithms, like message passing \cite{} \textbf{Find out if this is true}.
\textbf{Think of other examples.}
While it is still not known how to predict which marginal minima will be
attractors, this ubiquity of behavior suggests that cartography of marginal
minima is a useful step in bounding out-of-equilibrium dynamical behavior.

In the traditional methods for analyzing the geometric structure of rugged
landscapes, it is not necessarily straightforward to condition an analysis on
the marginality of minima. Using the method of a Legendre transformation of the
Parisi parameter corresponding to a set of real replicas, one can force the
result to be marginal by restricting the value of that parameter, but this
results in only the marginal minima at the energy level at which they are the
majority of stationary points \cite{Monasson_1995_Structural}. It is now
understood that out-of-equilibrium dynamics usually goes to marginal minima at
other energy levels \cite{Folena_2023_On}.

The alternative, used to great success in the spherical models, is to start by
making a detailing understanding of the Hessian matrix at stationary points.
Then, one can condition the analysis on whatever properties of the Hessian are
necessary to lead to marginal minima. This strategy is so successful in the
spherical models because it is very straightforward to implement: a natural
parameter in the analysis of these models linearly shifts the spectrum of the
Hessian, and so fixing this parameter by whatever means naturally allows one to
require that the Hessian spectrum have a pseudogap.
Unfortunately this strategy is less straightforward to generalize. Many models
of interest, especially in inference problems, have Hessian statistics that are
poorly understood.

Here, we introduce a generic method for conditioning the statistics of
stationary points on their marginality. The technique makes use of a novel way
to condition an integral over parameters to select only those that result in a
certain value of the smallest eigenvalue of a matrix that is a function of
those parameters. By requiring that the smallest eigenvalue of the Hessian at
stationary points be zero, we restrict to marginal minima, either those with a
pseudogap in their bulk spectrum or those with outlying eigenvectors. We
provide a heuristic to distinguish these two cases. We demonstrate the method
on the spherical models, where it is unnecessary but instructive, and on
extensions of the spherical models with non-GOE Hessians where the technique is
more useful.

\section{How to condition on the smallest eigenvalue}

An arbitrary function $g$ of the minimum eigenvalue of a matrix $A$ can be expressed as
\begin{equation}
  g(\lambda_\textrm{min}(A))
  =\lim_{\beta\to\infty}\int\frac{d\mathbf x\,\delta(N-\mathbf x^T\mathbf x)e^{-\beta\mathbf x^TA\mathbf x}}{\int d\mathbf x'\,\delta(N-\mathbf x'^T\mathbf x')e^{-\beta\mathbf x'^TA\mathbf x'}}g\left(\frac{\mathbf x^TA\mathbf x}N\right)
\end{equation}
\begin{equation}
  \lim_{\beta\to\infty}\frac{d\mathbf x\,\delta(N-\mathbf x^T\mathbf x)e^{-\beta\mathbf x^TA\mathbf x}}{\int d\mathbf x'\,\delta(N-\mathbf x'^T\mathbf x')e^{-\beta\mathbf x'^TA\mathbf x'}}
  =d\mathbf x\,\frac12\left[\delta(\mathbf x_\mathrm{min}(A)-\mathbf x)+\delta(\mathbf x_\mathrm{min}(A)+\mathbf x)\right]
\end{equation}
The first equality makes use of the normalized eigenvector $x_\mathrm{min}(A)$
associated with the minimum eigenvalue. By definition,
$x_\mathrm{min}(A)^TAx_\mathrm{min}(A)=x_\mathrm{min}(A)^Tx_\mathrm{min}(A)\lambda_\mathrm{min}(A)=N\lambda_\mathrm{min}(A)$
assuming the normalization is $\|x_\mathrm{min}(A)\|^2=N$. The second equality
extends a technique first introduced in \cite{Ikeda_2023_Bose-Einstein-like}
and used in \cite{Kent-Dobias_2024_Arrangement}. A Boltzmann distribution is introduced over a spherical
model whose Hamiltonian is quadratic with interaction matrix given by $A$. In
the limit of zero temperature, the measure will concentrate on the ground
states of the model, which correspond with the eigenvectors $\pm x_\mathrm{min}$
associated with the minimal eigenvalue $\lambda_\mathrm{min}$.

Consider a matrix $A=B+\omega I$ for $B$ a GOE matrix with entries whose variance is $\sigma^2/N$. As an example, we compute
\begin{equation}
  e^{NG_\sigma(\omega)}=\overline{\lim_{\beta\to\infty}\int\frac{d\mathbf x\,\delta(N-\mathbf x^T\mathbf x)e^{-\beta\mathbf x^T(B+\omega I)\mathbf x}}{\int d\mathbf x'\,\delta(N-\mathbf x'^T\mathbf x')e^{-\beta\mathbf x'^T(B+\omega I)\mathbf x'}}\,\delta\big(\mathbf x^T(B+\omega I)\mathbf x\big)}
\end{equation}
where the overline is the average over $B$. Using replicas to treat the
denominator and transforming the $\delta$-function to its Fourier
representation, we have
\begin{equation}
  e^{NG_\sigma(\omega)}=\overline{\lim_{\beta\to\infty}\lim_{n\to0}\int d\lambda\prod_{a=1}^n\left[d\mathbf x_a\,\delta(N-\mathbf x_a^T\mathbf x_a)\right]
  \exp\left\{-\beta\sum_{a=1}^n\mathbf x_a^T(B+\omega I)\mathbf x_a+\lambda\mathbf x_1^T(B+\omega I)\mathbf x_1\right\}}
\end{equation}
Taking the average over $B$, we have
\begin{equation}
  e^{NG_\sigma(\omega)}
  =\lim_{\beta\to\infty}\lim_{n\to0}\int d\lambda\prod_{a=1}^n\left[d\mathbf x_a\,\delta(N-\mathbf x_a^T\mathbf x_a)\right]
  \exp\left\{-Nn\beta\omega+N\lambda\omega+\frac{\sigma^2}{N}\left[\beta^2\sum_{ab}^n(\mathbf x_a^T\mathbf x_b)^2
    -2\beta\lambda\sum_a^n(\mathbf x_a^T\mathbf x_1)^2
    +\lambda^2N^2
  \right]\right\}
\end{equation}
We make the Hubbard--Stratonovich transformation to the matrix field $Q_{ab}=\frac1N\mathbf x_a^T\mathbf x_b$. This gives
\begin{equation}
  e^{NG_\sigma(\omega)}
  =\lim_{\beta\to\infty}\lim_{n\to0}\int d\lambda\,dQ\,
  \exp N\left\{
    -n\beta\omega+\lambda\omega+\sigma^2\left[\beta^2\sum_{ab}^nQ_{ab}^2
      +-\beta\lambda\sum_a^nQ_{1a}^2
    +\lambda^2
  \right]+\frac12\log\det Q\right\}
\end{equation}
where $Q_{aa}=1$ because of the spherical constraint. We can evaluate this
integral using the saddle point method. We make a replica symmetric ansatz for
$Q$, because this is a 2-spin model, but with the first row singled out because
of its unique coupling with $\lambda$. This gives
\begin{equation}
  Q=\begin{bmatrix}
    1&\tilde q_0&\tilde q_0&\cdots&\tilde q_0\\
    \tilde q_0&1&q_0&\cdots&q_0\\
    \tilde q_0&q_0&1&\ddots&q_0\\
    \vdots&\vdots&\ddots&\ddots&\vdots\\
    \tilde q_0&q_0&q_0&\cdots&q_0
  \end{bmatrix}
\end{equation}
with $\sum_{ab}Q_{ab}^2=n+2(n-1)\tilde q_0^2+(n-1)(n-2)q_0^2$, $\sum_aQ_{1a}^2=1+(n-1)\tilde q_0^2$,
and
\begin{equation}
  \log\det Q=(n-2)\log(1-q_0)+\log(1+(n-2)q_0-(n-1)\tilde q_0^2)
\end{equation}
Inserting these expressions and taking the limit of $n$ to zero, we find
\begin{equation}
  \mathcal S(q_0,\tilde q_0,\lambda)=\lambda\omega+\sigma^2\left[
    2\beta^2(q_0^2-\tilde q_0^2)-2\beta\lambda(1-\tilde q_0^2)+\lambda^2
  \right]-\log(1-q_0)+\frac12\log(1-2q_0+\tilde q_0^2)
\end{equation}
The integral is then given by its value at the stationary point of this
expression with respect to its three arguments.
The extremal conditions are
\begin{align}
  0&=\frac{\partial\mathcal S}{\partial q_0}
  =\frac1{1-q_0}-\frac1{1-2q_0+\tilde q_0^2}+4\beta^2\sigma^2q_0 \\
  0&=\frac{\partial\mathcal S}{\partial \tilde q_0}
  =\frac{\tilde q_0}{1-2q_0+\tilde q_0^2}-4\sigma^2(\beta^2-\beta\lambda)\tilde q_0 \\
  0&=\frac{\partial\mathcal S}{\partial\lambda}
  =\omega+2\sigma^2\big(\lambda-\beta(1-\tilde q_0^2)\big)
\end{align}
We expect the overlaps to concentrate on one as $\beta$ goes to infinity. We therefore take
\begin{align}
  q_0=1-y\beta^{-1}-z\beta^{-2}+O(\beta^{-3})
  &&
  \tilde q_0=1-\tilde y\beta^{-1}-\tilde z\beta^{-2}+O(\beta^{-3})
\end{align}
The first equations expanded in $\beta$ give
\begin{align}
  &0=4\sigma^2\beta^2+\bigg(\frac1{y}-\frac12\frac1{y-\tilde y}-4y\sigma^2\bigg)\beta+O(\beta^0) \\
  &0=-4\sigma^2\beta^2+\bigg(\frac12\frac1{y-\tilde y}+4\sigma^2(\lambda+\tilde y)\bigg)\beta+O(\beta^0)
\end{align}
One cannot satisfy this equation order-by-order in $\beta$. However, a solution
suggests itself: the expansion is singular when $\tilde y=y$. Making this
identification, we find instead
\begin{align}
  &0=\left(4\sigma^2-\frac1{y^2+2(z-\tilde z)}\right)\beta^2+\left(\frac1y+\frac{2y\tilde z}{(y^2+2(z-\tilde z))^2}-4\sigma^2y\right)\beta+O(\beta^0) \\
  &0=\left(-4\sigma^2+\frac1{y^2+2(z-\tilde z)}\right)\beta^2+\left(-\frac y{y^2+2(z-\tilde z)}-\frac{2y\tilde z}{(y^2+2(z-\tilde z))^2}+4\sigma^2(y+\lambda)\right)\beta+O(\beta^0)
\end{align}
Finally, expanding the equation for $\lambda$ to lowest order, we have
\begin{equation}
  0=\omega+2\sigma^2(\lambda-2y)+O(\beta^{-1})
\end{equation}
Simultaneously solving these five equations stemming from the coefficients of $\beta$ for $y$, $z$, $\tilde z$, and $\lambda$, we have
\begin{align}
  \lambda=-\frac1\sigma\sqrt{\frac{\omega^2}{(2\sigma)^2}-1}
  &&
  y=\frac1{2\sigma}\left(\frac{\omega}{2\sigma}-\sqrt{\frac{\omega^2}{(2\sigma)^2}-1}\right)
  \\
  z=\frac1{2\sigma^2}\left(1-\frac{\omega^2}{(2\sigma)^2}\right)
  &&
  \tilde z=\frac1{4\sigma^2}\left(1-\frac{\omega}{2\sigma}\left(\frac\omega{2\sigma}+\sqrt{\frac{\omega^2}{(2\sigma)^2}-1}\right)\right)
\end{align}
Inserting this solution into $\mathcal S$ and taking the limit of $\beta$ to zero, we find
\begin{equation}
  G_\sigma(\omega)=-\frac{\omega}{2\sigma}\sqrt{\frac{\omega^2}{(2\sigma)^2}-1}
  +\log\left[
    \frac{\omega}{2\sigma}+\sqrt{\frac{\omega^2}{(2\sigma)^2}-1}
  \right]
\end{equation}
This function is plotted in Fig. For $\omega<2\sigma$ $G_\sigma(\omega)$ has an
imaginary part, which makes any additional integral over $\omega$ highly
oscillatory. This indicates that the existence of a marginal minimum for this
parameter value corresponds with a large deviation that grows faster than $N$,
rather like $N^2$, since in this regime the bulk of the average spectrum is
over zero and therefore extensively many eigenvalues have to have large
deviations in order for the smallest eigenvalue to be zero. For
$\omega\geq2\sigma$ this function gives the large deviation function for the
probability of seeing a zero eigenvalue given the shift $\omega$.
$\omega=2\sigma$ is the maximum of the function with a real value, and
corresponds to the intersection of the average spectrum with zero.


\begin{equation}
  H(\mathbf s)+\sum_{i=1}^r\omega_ig_i(\mathbf s)
\end{equation}
\begin{align}
  \nabla H(\mathbf s,\pmb\omega)=\partial H(\mathbf s)+\sum_{i=1}^r\omega_i\partial g_i(\mathbf s)
  &&
  \operatorname{Hess}H(\mathbf s,\pmb\omega)=\partial\partial H(\mathbf s)+\sum_{i=1}^r\omega_i\partial\partial g_i(\mathbf s)
\end{align}

\begin{equation}
  d\mu_H(\mathbf s,\pmb\omega)=d\mathbf s\,d\pmb\omega\,\delta\big(\nabla H(\mathbf s,\pmb\omega)\big)\,\delta\big(\mathbf g(\mathbf s)\big)\,\big|\det\operatorname{Hess}H(\mathbf s,\pmb\omega)\big|
\end{equation}
\begin{equation}
  d\mu_H(\mathbf s,\pmb\omega\mid E)=d\mu_H(\mathbf s,\pmb\omega)\,\delta\big(NE-H(\mathbf s)\big)
\end{equation}

\begin{equation}
  \begin{aligned}
    &\mathcal N_\text{marginal}(E)
    =\int d\mu_H(\mathbf s,\pmb\omega\mid E)\,\delta\big(N\lambda_\mathrm{min}(\operatorname{Hess}H(\mathbf s,\pmb\omega))\big) \\
    &=\lim_{\beta\to\infty}\int d\mu_H(\mathbf s,\pmb\omega\mid E)
    \frac{d\mathbf x\,\delta(N-\mathbf x^T\mathbf x)\delta(\mathbf x^T\partial\mathbf g(\mathbf s))e^{\beta\mathbf x^T\operatorname{Hess}H(\mathbf s,\pmb\omega)\mathbf x}}
    {\int d\mathbf x'\,\delta(N-\mathbf x'^T\mathbf x')\delta(\mathbf x'^T\partial\mathbf g(\mathbf s))e^{\beta\mathbf x'^T\operatorname{Hess}H(\mathbf s,\pmb\omega)\mathbf x'}}
    \delta\big(\mathbf x^T\operatorname{Hess}H(\mathbf s,\pmb\omega)\mathbf x\big)
  \end{aligned}
\end{equation}
where the $\delta$-functions
\begin{equation}
  \delta(\mathbf x^T\partial\mathbf g(\mathbf s))
  =\prod_{s=1}^r\delta(\mathbf x^T\partial g_i(\mathbf s))
\end{equation}
ensure that the integrals are constrained to the tangent space of the configuration manifold at the point $\mathbf s$.

\begin{equation}
  \begin{aligned}
    &\Sigma_\text{marginal}(E)
    =\frac1N\overline{\log\mathcal N_\text{marginal}(E)} \\
    &=\lim_{\beta\to\infty}\lim_{n\to0}\frac\partial{\partial n}\int\prod_{a=1}^n\left[d\mu_H(\mathbf s_a,\pmb\omega_a\mid E)\lim_{m_a\to0}
    \left(\prod_{b=1}^{m_a} d\mathbf x_a^b\,\delta(N-(\mathbf x_a^b)^T\mathbf x_a^b)\delta((\mathbf x_a^b)^T\partial\mathbf g(\mathbf s_a))e^{\beta(\mathbf x_a^b)^T\operatorname{Hess}H(\mathbf s_a,\pmb\omega_a)\mathbf x_a^b}\right)\,\delta\big((\mathbf x_a^1)^T\operatorname{Hess}H(\mathbf s_a,\pmb\omega_a)\mathbf x_a^1\big)\right]
  \end{aligned}
\end{equation}

\section{Application to the spherical models}

\begin{align}
  C_{ab}=\frac1N\mathbf s_a\cdot\mathbf s_b
  &&
  R_{ab}=-i\frac1N\mathbf s_a\cdot\hat{\mathbf s}_b
  &&
  D_{ab}=\frac1N\hat{\mathbf s}_a\cdot\hat{\mathbf s}_b
  \\
  A_{ab}^{cd}=\frac1N\mathbf x_a^c\cdot\mathbf x_b^d
  &&
  X^c_{ab}=\frac1N\mathbf s_a\cdot\mathbf x_b^c
  &&
  \hat X^c_{ab}=\frac1N\hat{\mathbf s}_a\cdot\mathbf x_b^c
\end{align}

\begin{equation}
  \begin{aligned}
    &\sum_{ab}^n\left[\beta\omega A_{aa}^{bb}+\hat x\omega A_{aa}^{11}+\beta^2f''(1)\sum_{cd}^m(A_{ab}^{cd})^2+\hat x^2f''(1)(A_{ab}^{11})^2+\beta\hat xf''(1)\sum_c^m A_{ab}^{1c}\right]\\
    &+\hat\beta^2f(C_{ab})+(2\hat\beta(R_{ab}-F_{ab})-D_{ab})f'(C_{ab})+(R_{ab}^2-F_{ab}^2)f''(C_{ab})
    +\log\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix}-\log\det F
  \end{aligned}
\end{equation}

$X^a$ is $n\times m_a$, and $A^{ab}$ is $m_a\times m_b$.
\begin{equation}
  \begin{bmatrix}
    C&iR&X^1&\cdots&X^n \\
    iR&D&i\hat X^1&\cdots&i\hat X^m\\
    (X^1)^T&i(\hat X^1)^T&A^{11}&\cdots&A^{1n}\\
    \vdots&\vdots&\vdots&\ddots&\vdots\\
    (X^n)^T&i(\hat X^n)^T&A^{n1}&\cdots&A^{nn}
  \end{bmatrix}
\end{equation}
$X_{ab}^c$ will be nonzero if the lowest eigenvector of the hessian at the
point $\mathbf s_c$ are correlated with the direction of the point $\mathbf
s_a$. Since the eigenvector problem is always expected to be replica symmetric,
we expect no $b$-dependence of this matrix. $A^{aa}$ is the usual
replica-symmetric overlap matrix of the spherical two-spin problem. $A^{ab}$
describes overlaps between eigenvectors at different stationary points and should be a constant $m_a\times m_b$ matrix.

We will discuss at the end of this paper when these order parameters can be expected to be nonzero, but in this and most isotropic problems all of the $X$s, $\hat X$s, and $A^{ab}$ for $a\neq b$ are zero.

\begin{equation}
  \Sigma_\textrm{marginal}(E)
  =\operatorname{max}_\omega\big[\Sigma(E,\omega)+G_{\sqrt{f''(1)}}(\omega)\big]
\end{equation}
where the maximum over $\omega$ needs to lie at a real value.

\section{Twin spherical model}

$\Omega=S^{N-1}\times S^{N-1}$
\begin{equation}
  H(\pmb\sigma)=H_1(\pmb\sigma^{(1)})+H_2(\pmb\sigma^{(2)})+\epsilon\pmb\sigma^{(1)}\cdot\pmb\sigma^{(2)}
\end{equation}
\begin{equation}
  \overline{H_s(\pmb\sigma_1)H_s(\pmb\sigma_2)}
  =Nf_s\left(\frac{\pmb\sigma_1\cdot\pmb\sigma_2}N\right)
\end{equation}

\begin{equation}
  \mathcal S(C,R,D,W,\hat\beta,\omega)
  =\frac12\frac1n
    \sum_{ab}\left(
      \hat\beta^2f(C_{ab})+(2\hat\beta R_{ab}-D_{ab})f'(C_{ab})+(R_{ab}^2-W_{ab}^2)f''(C_{ab})
    \right)
\end{equation}

\begin{equation}
  \mathcal S(C^{11},R^{11},D^{11},W^{11},\hat\beta)+\mathcal S(C^{22},R^{22},D^{22},W^{22},\hat\beta)
  -\epsilon(r_{12}+r_{21})-\omega_1(r^{11}_d-w^{11}_d)-\omega_2(r^{22}_d-w^{22}_d)+\hat\beta E
  +\frac12\log\det\begin{bmatrix}C^{11}&iR^{11}\\iR^{11}&D^{11}\end{bmatrix}
  +\frac12\log\det\left(
  \begin{bmatrix}C^{22}-q_{12}^2C^{11}&iR^{22}\\iR^{22}&D^{22}\end{bmatrix}
  \right)
  -\log\det(W^{11}W^{22}+W^{12}W^{21})
\end{equation}

\begin{equation}
  \begin{aligned}
    &\sum_a^n\left[\hat q_a(Q^{11}_{aa}+Q^{22}_{aa}-1)-\beta(\omega_1Q^{11}_{aa}+\omega_2Q^{22}_{aa}+2\epsilon Q^{12}_{aa})\right]
    +\lambda(\omega_1Q^{11}_{11}+\omega_2Q^{22}_{11}+2\epsilon Q^{12}_{11}) \\
    &+\sum_{i=1,2}f_i''(1)\left[\beta^2\sum_{ab}^n(Q^{ii}_{ab})^2-2\beta\lambda\sum_a^n(Q^{ii}_{1a})^2+\lambda^2(Q^{ii}_{11})^2\right]
    +\frac12\log\det\begin{bmatrix}
      Q^{11}&Q^{12}\\
      Q^{12}&Q^{22}
    \end{bmatrix}
  \end{aligned}
\end{equation}
\begin{equation}
  \log\det\begin{bmatrix}
    Q^{11}&Q^{12}\\
    Q^{12}&Q^{22}
  \end{bmatrix}
  +\log\det(Q^{11}Q^{22}-Q^{12}Q^{12})
\end{equation}

\section{Multi-species spherical model}

We consider models whose configuration space consists of the product of $r$
spheres, each with its own dimension $N_s$, or
$\Omega=S^{N_1-1}\times\cdots\times S^{N_r-1}$. Coordinates on this space we
will typically denote
$\pmb\sigma=(\pmb\sigma^{(1)},\ldots,\pmb\sigma^{(r)})\in\Omega$, with
$\pmb\sigma^{(s)}\in S^{N_s-1}$ denoting the subcomponent restricted to a
specific subsphere. The model can be thought of as consisting of centered
random functions $H:\Omega\to\mathbb R$ with covariance
\begin{equation}
  \overline{
    H(\pmb\sigma_1)H(\pmb\sigma_2)
  }
  =f\left(
      \frac{\pmb\sigma^{(1)}_1\cdot\pmb\sigma^{(1)}_2}{N_1},
      \ldots,
      \frac{\pmb\sigma^{(r)}_1\cdot\pmb\sigma^{(r)}_2}{N_r}
    \right)
\end{equation}
where $f:[-1,1]^r\to\mathbb R$ is an $r$-component function of the overlaps that defines the model.

\printbibliography

\end{document}