%~Mouliné par MaN_auto v.0.32.4 2024-08-19 12:29:01
\documentclass[CRPHYS, Unicode, XML, screen, biblatex]{cedram}

\addbibresource{crphys20240530.bib}


\newcommand{\bfS}{\mathbf{S}}
\newcommand{\bfJ}{\mathbf{J}}
\newcommand{\rPro}{\mathrm{Pro}}
%\newcommand{\bbR}{\mathbb{R}}
%\newcommand{\bbC}{\mathbb{C}}
%\newcommand{\bbZ}{\mathbb{Z}}
%\newcommand{\cal}[1]{\mathcal{#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}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\CDRsetmeta{articletype}{review}
%
%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Special issue dedicated to Gérard Toulouse}
%\def\cdr@specialtitle@french{Numéro spécial dédié à Gérard Toulouse}
%\makeatother{}



\title{Generalizations of Parisi's replica symmetry breaking and overlaps in random energy models}

\alttitle{Généralisations de la brisure de symétrie des répliques de Parisi et des overlaps dans les modèles d’énergies aléatoires}


\author{\firstname{Bernard} \lastname{Derrida}\IsCorresp \CDRorcid{0000-0001-6994-0226}}
\address{Collège de France, 11 place Marcelin Berthelot, 75005 Paris, France}
\address{Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France}
\email[B. Derrida]{bernard.derrida@college-de-france.fr}

\author{\firstname{Peter} \lastname{Mottishaw}\CDRorcid{0000-0002-0091-4094}}
\address{SUPA, School of Physics and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9 3FD, United Kingdom}
\email[P. Mottishaw]{peter.mottishaw@ed.ac.uk}

\keywords{Disordered systems, Spin glasses, Replica symmetry breaking, Random Energy Model}

\altkeywords{Systèmes désordonnés, Verres de spin, Brisure de symétrie des répliques, Modèle d'énergies aléatoires}

\begin{abstract}
The random energy model (REM) is the simplest spin glass model which exhibits replica symmetry breaking. It is well known since the 80's that its overlaps are non-selfaveraging and that their statistics satisfy the predictions of the replica theory. All these statistical properties can be understood by considering that the low energy levels are the points generated by a Poisson process with an exponential density. Here we first show how, by replacing the exponential density by a sum of two exponentials, the overlaps statistics are modified. One way to reconcile these results with the replica theory is to allow the blocks in the Parisi matrix to fluctuate. Other examples where the sizes of these blocks should fluctuate include the finite size corrections of the REM, the case of discrete energies and the overlaps between two temperatures. In all these cases, the block sizes not only fluctuate but need to take complex values if one wishes to reproduce the results of our replica-free calculations.
\end{abstract}

\begin{altabstract}
Le modèle d'énergies aléatoires (REM) est le modèle de verre de spin le plus simple qui présente une brisure de symétrie des répliques. Il est bien connu depuis les années 80 que ses overlaps ne sont pas automoyennants et que leurs statistiques sont celles prédites par la méthode des répliques. Ces propriétés statistiques peuvent être comprises en considérant que les niveaux d'énergie les plus bas sont les points générés par un processus de Poisson de densité exponentielle. Nous montrons ici dans un premier temps comment ces statistiques d'overlaps sont modifiées lorsqu'on remplace la densité exponentielle par une somme de deux exponentielles. Une fa\c{c}on de concilier ces résultats avec la théorie des répliques est de permettre aux blocs de la matrice de Parisi de fluctuer. D'autres exemples où la taille de ces blocs doit fluctuer incluent les corrections de taille finie du REM, le cas des énergies discrètes et les overlaps entre deux températures. Dans tous ces cas, non seulement la taille des blocs fluctue mais elle doit prendre des valeurs complexes si l'on souhaite reproduire nos résultats obtenus directement, c'est à dire sans utiliser la méthode des répliques.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\maketitle

\section{Introduction}


The introduction by Edwards and Anderson~\cite{Edwards_1975_Theory} in 1975 of the replica method and of the overlaps as an order parameter was a big step in the theory of spin glasses~\cite{Parisi_2004_OVERLAP}. For a spin glass model with $N$ Ising spins, the overlap $q ({\bf S,S'}) $ between two spin configurations $\bfS \equiv \{S_i \pm 1\}$ and ${\bf S'} \equiv \{S_i' \pm 1\}$ is defined as
\begin{equation}\label{A1}
q\left({\bf S,S'}\right) =\frac{1}{N} \sum_{i=1,N} S_i \, S_i'
\end{equation}
Qualitatively, the idea was that, because spin glasses have a rugged energy landscape, two typical spin configurations at equilibrium have, at low temperature, the tendency of being trapped in the same valley giving rise to a non-zero overlap. Quantitatively, for a spin glass model with quenched pair interactions $\bfJ\equiv \{J_{i,j}\}$
\begin{equation}\label{A2}
E_{\bfJ }(\bfS) = - \sum_{i,j} J_{i,j} S_i S_j
\end{equation}
the overlaps can be characterized by their probability distribution
\begin{equation}\label{A3}
P_{\bfJ }(Q) = \frac{ \sum_\bfS \sum_{\bf S'} \delta\Big(q(\bfS,{\bf S'})-Q \Big)
\ \exp\left[-\beta \Big(E_\bfJ(\bfS)+ E_\bfJ(\bfS') \Big) \right]}{
\sum_\bfS \sum_{\bf S'}
\ \exp\left[-\beta \Big(E_\bfJ(\bfS)+ E_\bfJ(\bfS') \Big) \right] }
\end{equation}
where $\beta$ is the inverse temperature. As long as the system is finite ($N < \infty$) this distribution $P_{\bfJ }(Q)$ is a broad function of $Q$ and depends on the sample $\bfJ$. However in the thermodynamic limit, it was initially expected that, in the spin-glass phase and for almost all samples $\bfJ$, it becomes a sum of two delta functions
\begin{equation}\label{A4}
\lim_{N\,\to\,\infty} P_{\bfJ}(Q) = \frac{1}{2} \left[\delta\left(Q-q_{EA}\right) + \delta\left(Q+q_{EA}\right) \right]
\end{equation}
and in~\cite{Edwards_1975_Theory} a mean field theory was developed to determine the value of the Edwards--Anderson order parameter $q_{EA}$. (if some odd interactions were added to~\eqref{A2} the limit~\eqref{A4} would reduce to a single delta function $\lim_{N\,\to\,\infty} P_{\bfJ}(Q) = \delta(Q-q_{EA}) $).

Soon after the Edwards Anderson 1975 paper~\cite{Edwards_1975_Theory}, Sherrington and Kirkpatrick~\cite{Sherrington_1975_Solvable} considered the infinite range version of the model~\eqref{A2}, for which the mean field approximation was expected to become exact. Using the replica approach, they could obtain explicit expressions of $q_{EA}$ and of the average free energy. However they realized that their analytic solution could not be correct because, at low temperature, it leads to a negative entropy~\cite{Sherrington_1975_Solvable} and to a negative variance of the free-energy~\cite{Toulouse_1981_Free}

It was only in 1979 that Parisi was able to overcome these difficulties with a Replica Symmetry Breaking (RSB) scheme which turned out to give the correct free energy for the Sherrington Kirkpatrick model. Initially, the solution was based on unconventional mathematics such as using matrices of non-integer size or replacing maxima by minima. It took then a few more years before the physical interpretation of Parisi's solution was understood~\cite{Parisi_1983_Order,Mezard_1984_Nature,Mezard_1984_Replica,Mezard_1987_Spin} and even longer before it was confirmed by a series of rigorous mathematical proofs~\cite{Guerra_2003_Broken,Talagrand_2006_Parisi,Arguin_2009_structure,Panchenko_2013_Parisi}. See~\cite{Charbonneau_2023_Spin} for a recent review of the theory of spin glasses and replica symmetry breaking. Besides an analytic way of calculating the exact free energy of the Sherrington--Kirkpatrick model, Parisi's solution led to a number of surprising predictions. One of them was that, instead of~\eqref{A4}, the distribution $P_{\bfJ}(Q)$ remains broad and sample dependent even in the thermodynamic limit ($N \to \infty$). Moreover it gave an explicit way of calculating these fluctuations~\cite{Mezard_1984_Nature,Mezard_1984_Replica,Mezard_1987_Spin}.

One way of describing these fluctuations is to consider the cumulative function
\begin{equation}\label{A5}
Y(Q)= \lim_{N\,\to\,\infty}\int_Q^1 P_{\bfJ}\left(Q'\right) dQ'
\end{equation}
which represents the probability that, at equilibrium, two configurations in the same energy landscape have an overlap larger than $Q$. According to the RSB theory (see Appendix A), $Y(Q)$ remains sample dependent even in the thermodynamic limit and its probability distribution $\Pi_\mu(Y)$ is universal for all the models which can be solved by the RSB theory. This distribution $\Pi_\mu(Y)$ is indexed by a single parameter $0\leq \mu \leq 1 $ which depends on $Q$, on the temperature, on the magnetic field and on all the other parameters which may characterize a specific model. As a consequence if $ \langle Y \rangle $ is known for a given system (where $\langle. \rangle$ denotes the average over the samples, i.e. over the $\bfJ$'s) all the moments of $Y$ are known. For example if $ \langle Y \rangle=1-\mu $, one has
\begin{equation}\label{A6}
\langle Y \rangle = 1-\mu\; ; \quad
\left\langle Y^2 \right\rangle = \frac{(1-\mu) (3-2 \mu)}{3} \; ; \quad
\left\langle Y^3 \right\rangle = \frac{(1-\mu)\left(15-17 \mu + 5 \mu^2\right)}{15}
\end{equation}
A typical shape of the distribution $\Pi_\mu(Y)$ is shown in Figure~\ref{figure1}~\cite{Mezard_1984_Replica,Derrida_1987_Statistical}


\begin{figure}[h]
\centerline{\includegraphics[width=10.5cm]{Figure-1-bitmap}}
\vspace*{-15pt}
\caption{The distribution $\Pi_\mu (Y)$ as a function of $Y$ in the case $\mu=1/3$.}\label{figure1}
\end{figure}
More generally one can consider, for a given sample $\bfJ$, the probability $Y_k(Q)$ that $k$ configurations have all their overlaps larger than $Q$, or the probability $Y_{k,k'}(Q)$ that, for two groups of $k$ and $k'$ configurations, all the overlaps between pairs of configurations inside each group are larger than $Q$ but all pairs between two groups have an overlap less than $Q$. According to the RSB theory (see Appendix A), the averages of these $Y_k$ or $Y_{k,k'} $ over the samples have also explicit expressions in terms of the parameter $\mu$
\begin{equation}\label{A7}
\langle Y_k \rangle = \frac{\Gamma(k-\mu)}{\Gamma(k) \ \Gamma(1-\mu) }
\; ; \quad \left\langle Y_{k,k'} \right\rangle = \frac{\mu \ \Gamma(k-\mu) \ \Gamma(k'-\mu)}{\Gamma(k+k') \ \Gamma(1-\mu)^2 }
\end{equation}
These satisfy particular cases of the Ghirlanda-Guerra relations~\cite{Ghirlanda_1998_General}: for example it is easy to check that
\[
\left\langle Y_{k+1}\right\rangle = \frac{\langle Y_{2} \rangle \,
\langle Y_{k} \rangle + (k-1)\langle Y_{k} \rangle}{k}
\]

All these statistical properties of the $Y_k$'s were first derived in 1984 using the replica method (see Appendix~\ref{appA}). In a joint work~\cite{Derrida_1985_Sample} with Gérard Toulouse in 1985, they were confirmed by a replica-free calculation for the Random Energy Model (REM).


The goal of the present paper is to exhibit several simple models where these statistical properties are no longer valid and have to be modified. As explained in Section~\ref{sect3}, these models are cousins of the Random Energy Model, in the sense that the energies of the configurations are independent random variables and that their properties can be obtained using replica-free methods. Based on some exact expressions \eqref{D8}-\eqref{D9} and \eqref{D11}-\eqref{D12} which replace~\eqref{A7} for these models, we will discuss possible generalizations of the Parisi matrix (see the Appendix~\ref{appA}) where the sizes of the blocks have to fluctuate (sometimes with complex sizes !).


The paper is organized as follows. In Section~\ref{sect2}, we recall some well known properties of the REM and why, in the low temperature phase, the energies can be generated by a Poisson process with an exponential density. In Section~\ref{sect3}, we establish some general expressions allowing us to calculate non-integer moments of the partition function as well as average overlaps $\langle Y_k\rangle, \langle Y_{k,k'} \rangle $ for a Poisson process with an arbitrary density. In the case of an exponential density, this allows one to recover~\eqref{A7}. In Section~\ref{sect4} we consider the case where the density of the Poisson process is the sum of two exponentials and show how the expressions~\eqref{A7} are modified. We will then suggest a way of adapting the RSB approach to reproduce these new expressions. This is probably the simplest case for which the size of the blocks in the Parisi matrix have to fluctuate. A straightforward generalization to the case of a sum of an arbitrary number of exponentials will allow us to recover, in a unified way, earlier results on finite size corrections and in the case of discrete energies. In Section~\ref{sect5} we generalize the approach to consider the overlaps between two configurations at different temperatures, and show that, again, the size of the blocks of the Parisi matrix have to fluctuate. Finally, in Section~\ref{sect6}, we discuss the replica method for a Poisson REM with an arbitrary density of energy, and see how the replica scheme has to be modified in this case.



\section{Short review of known results on the REM}\label{sect2}

In this section, we first recall a few known features of the REM, in particular the fact that, in the thermodynamic limit, the low temperature properties can be understood by considering that the energies are the realization of a Poisson process with an exponential density~\eqref{B15}.


\subsection{The REM with Gaussian energies}

The Random Energy model was introduced~\cite{Derrida_1980_Random,Derrida_1981_Random} in 1980 as a simple spin glass model which exhibits a spin glass transition and which can be solved without resorting to the replica trick. As for other Ising spin glass models like the Edwards Anderson or the Sherrington--Kirkpatrick models, there are $2^N $ spin configurations $\bfS$ whose energies are Gaussian random variables
\begin{equation}\label{B1}
P(E) = \frac{1 }{ \sqrt{N \pi} \, J} \exp\left[-\frac{ E^2 }{ N \, J^2} \right]
\end{equation}
We will set $J=1$, for simplicity, in the rest of the paper. The only difference with other spin glass models is that, in the REM, the energies of different configurations are independent random variables. Therefore, instead of the interactions $\bfJ$ in~\eqref{A2}, a given sample $\bfJ$ in the REM is specified by $2^N$ random quenched variables $E(\bfS)$ distributed according to~\eqref{B1}. Then the partition function $Z(\beta)$ becomes a sum of $2^N$ independent random variables making many calculations much easier than for other spin-glass models
\begin{equation}\label{B2}
Z(\beta)=\sum_{\bfS=1}^{2^N} e^{- \beta E(\bfS)}
\end{equation}
Still the REM exhibits a phase transition~\cite{Derrida_1980_Random,Derrida_1981_Random} at some $\beta_c$.
\begin{equation}\label{B3}
\beta_c = 2 \sqrt{\log2}
\end{equation}
and, in the thermodynamic limit, the average free is given by
\begin{equation}\label{B3a}
\lim_{N\,\to\,\infty} \frac{\left\langle \log Z(\beta)\right\rangle }{ N} =
\left\{
\begin{array}{lll} \log 2 + \frac{\beta^2 J^2 }{ 4} & \text{ for} & \beta< \beta_c
\\ & & \\
\beta \, \sqrt{\log2} & \text{for} & \beta > \beta_c
\end{array} \right.
\end{equation}
where $\langle.\rangle$ denotes an average over the energies $E(\bfS)$. These properties have been confirmed in a number of rigorous works~\cite{Eisele_1983_third,Galves_1989_Fluctuations,Olivieri_1984_existence,Bovier_2002_Fluctuations}; for a review see~\cite{Bolthausen_2007_Random} or~\cite{Kistler_2015_Derridas}. Note that the REM is an example of sums of exponentials of random i.i.d. variables which have been considered in several other contexts~\cite{Pastur_1989_limit,BenArous_2005_Limit}.

Many other properties of the REM have been calculated (see~\cite{Derrida_2023_Random} for a review) like, finite size corrections~\cite{Derrida_1981_Random,Campellone_1995_Some,Campellone_2009_Replica,Derrida_2015_Finite}, the effect of a magnetic field~\cite{Derrida_1981_Random}, the location of the zeroes in the complex plane of $\beta$~\cite{Moukarzel_1991_Numerical,Derrida_1991_zeroes,
Saakian_2000_Random} and the fluctuations of the spectral form factor of chaotic systems~\cite{Bunin_2023_Fisher}, the effect of discrete energies~\cite{Ogure_2005_exact,Jana_2007_Contributions,Derrida_2022_discrete} or the integer and the non-integer moments of the partition function~\cite{Gardner_1989_probability}.

For example, in the low temperature phase (i.e. $\beta > \beta_c$), the negative moments are given~\cite{Gardner_1989_probability} for large $N$ by
\begin{equation}\label{B4}
\left\langle Z(\beta)^n \right\rangle = \left(\frac{A}{\beta_c} \ \Gamma(1-\mu) \right)^\frac{n}{\mu} \frac{\Gamma\left(1-\frac{n}{\mu}\right)}{\Gamma(1-n)} \quad \text{for} \quad n <0
\end{equation}
where
\begin{equation}\label{B5}
A = \frac{e^{N \beta_c \sqrt{\log2}}}{\sqrt{ \pi N} }
\end{equation}
and
\begin{equation}\label{B6}
\mu=\frac{\beta_c}{\beta}
\end{equation}

In the low temperature phase (see~\eqref{B3a}) the extensive part of the energy is constant. This is because the only configurations $\bfS$ which really contribute to the free energy are those whose energies are very close to the ground state. These configurations are likely to be very scattered in phase space and so to have zero overlap between themselves. Therefore, in the large $N$ limit for almost all samples, $P_\bfJ(Q)$ has the form
\begin{equation}\label{B7}
P_\bfJ (Q)= (1-Y_2)\, \delta(Q) + Y_2 \, \delta(Q-1)
\end{equation}
where $Y_2$ is the probability of finding at equilibrium two copies of the system in the same configuration
\begin{equation}\label{B8}
Y_2=
\frac{\sum_\bfS e^{-2 \beta E(\bfS)}}{\left(\sum_\bfS e^{- \beta E(\bfS)} \right)^2}
\end{equation}
Similarly
\begin{equation}\label{B9}
Y_k=
\frac{\sum_\bfS e^{-k \beta E(\bfS)}}{\left(\sum_\bfS e^{- \beta E(\bfS)} \right)^k}
\; ;
\quad Y_{k,k'}=
\frac{\sum_{\bfS\neq\bfS'} e^{-k \beta E(\bfS) -k' \beta E(\bfS')} }{
\left(\sum_\bfS e^{- \beta E(\bfS)} \right)^{k+k'} }
\end{equation}


These quantities depend on the sample (i.e. on the realization of the $2^N$ energies $E(\bfS)$). For the REM, their sample averages have been calculated in the low temperature phase~\cite{Derrida_1985_Sample} (see Section~\ref{sect3} for a simple derivation)
\begin{equation}\label{B10}
\frac{\left\langle Y_k \, Z(\beta)^n \right\rangle}{\left\langle Z(\beta)^n \right\rangle} = \frac{\Gamma(k-\mu)}{\Gamma(1-\mu)} \ \frac{\Gamma(1-n)}{\Gamma (k-n) }= \prod_{j=1}^{k-1} \frac{j-\mu}{j-n}
\end{equation}
\begin{equation}\label{B11}
\frac{\left\langle Y_{k,k'} \ Z(\beta)^n \right\rangle}{\left\langle Z(\beta)^n \right\rangle} = (\mu-n) \frac{\Gamma(k-\mu)}{\Gamma(1-\mu)} \ \frac{\Gamma(k'-\mu)}{\Gamma(1-\mu)} \ \frac{\Gamma(1-n)}{\Gamma (k+k'-n) }
\end{equation}
and they agree with the predictions of the replica theory~\cite{Mezard_1984_Nature,
Mezard_1984_Replica,Derrida_1997_random,
Derrida_2015_Finite}.

Note that setting $n=0$, the expressions~\eqref{B10} and~\eqref{B11} reduce to~\eqref{A7}. Also by expanding~\eqref{B10} in powers of $n$ one can see that the free-energy and the overlaps are correlated
\begin{equation}\label{B12}
\left\langle Y_k \log Z(\beta) \right\rangle - \langle Y_k \rangle \, \left\langle \log Z(\beta) \right\rangle = \frac{\Gamma(k-\mu)}{\Gamma(k) \, \Gamma(1 - \mu)} \left(\sum_{q=1}^{k-1} \frac{1}{q} \right)
\end{equation}

\subsection{Other distributions of energies}

All the above expressions \eqref{B4}, \eqref{B6}-\eqref{B12} remain valid for more general distributions of energies. For example if the distribution~\eqref{B1} of energies $E(\bfS)$ is, for large $N$, of the form
\begin{equation}\label{B13}
P(E) \simeq \sqrt\frac{G''\left(\frac{E}{N}\right)}{2 \pi N} \exp \left[- N G\left(\frac{E}{N} \right)\right]
\end{equation}
where $G$ is a convex function (this would be the case if each energy was the sum of $N$ i.i.d. random variables distributed according to a continuous distribution), the only change being that the value of $\beta_c$ and the constant $A$ in~\eqref{B3} and~\eqref{B5} would become
\begin{equation}\label{B14}
\beta_c=-G'(\epsilon_c)
\; ; \quad A= \sqrt\frac{ G''(\epsilon_c)}{2 \pi N} \, e^{-N \beta_c \epsilon_c}
\end{equation}
where $\epsilon_c$ is the minimal solution of $\log(2) - G(\epsilon_c)=0$.

Note, as we will see in Section~\ref{sect4}, that if the distribution of energies is not of the form~\eqref{B13}, in particular when the energies $E(\bfS)$ take only discrete integer values, the expressions of the overlaps will be quite different.

\subsection{The Poisson process}

Overlaps are non-zero only in the low temperature phase of the REM. This is why, in this whole paper, we limit our discussion to this low temperature phase where only the configurations whose energies are close to the ground state energy matter. The energies of these configurations for a given sample can be generated by a Poisson process of density~\cite{Derrida_2015_Finite}
\begin{equation}\label{B15}
\rho(E)= A \, e^{\beta_c E}
\end{equation}
which approximates $2^N P(E)$ in the neighborhood of the ground state energy. Depending on the choice of $P(E)$ in~\eqref{B1} or~\eqref{B13}, the amplitude $A$ and $\beta_c$ are given by \eqref{B3}-\eqref{B5} or~\eqref{B14}. In the next section we will see that~\eqref{B15} allows one to recover the predictions~\eqref{A7} of the replica approach.

One can notice that, for the density~\eqref{B15}, the average partition function $\langle Z(\beta)\rangle= \int \rho(E) d E \ \exp[-\beta E] $ is infinite for all values of $\beta$. However, in the low temperature phase $\beta > \beta_c$, (which is the only temperature range where the overlaps are non-zero), $Z(\beta)$ is finite with probability 1. (In fact it is easy to prove that for all $C >0$
\begin{equation}\label{B15a}
\rPro\Big(Z(\beta)< C\Big) \ > \
\rPro(Z_1 < C)) \ \rPro(Z_2=0) \ > \
\left(1-\frac{\langle Z_1 \rangle}{C} \right) \ \rPro(Z_2=0)
\end{equation}
where $Z(\beta) =Z_1+Z_2$ and $Z_1$ represents the contributions to $Z$ of the energies $E > \Lambda$ and $Z_2$ the contributions of the energies $E<\Lambda$. One has
\[
\langle Z_1 \rangle = A\, \frac{e^{(\beta_c-\beta)\Lambda}}{\beta-\beta_c} \quad \text{and}
\quad \rPro (Z_2=0)= \exp\left[-A\, \frac{e^{\beta_c \Lambda}}{\beta_c} \right]
\]
By choosing $\Lambda$ sufficiently negative and $C$ much larger than $\langle Z_1 \rangle$, one can make the r.h.s. of~\eqref{B15a} as close as needed to 1.) We will see in Section~\ref{sect3} that the Poisson process with the density~\eqref{B15} allows one to recover the known expressions of the negative moments~\eqref{B4} of the REM.


\subsection{Link with sums of random variables with a heavy tail}

In the low temperature phase, the partition function $Z(\beta)$ of the REM can be viewed as a sum of i.i.d. random variables distributed according to a heavy tailed distribution~\cite{Derrida_1997_random}. Let $E_1 < E_2 < \dots E_p$ be the $p$ lowest energies of a realization of the Poisson process with density~\eqref{B15}. The probability distribution of these energies is
\[
\rPro\left(E_1,\,\cdots E_p\right) =A^p \exp\left[\beta_c\left(E_1+ E_2 + \cdots E_p\right) - A \frac{e^{\beta_c E_p }}{\beta_c} \right]
\]
On the other hand if one considers a large number $M$ of i.i.d. random positive variables $x_i$ distributed according to a distribution $\rho(x)$ with an heavy tail
\[
\rho(x) \sim B \, x^{-1-\mu} \quad \text{for large}\ \ x
\]
and if one orders these $x_i$'s
the distribution of $x_1>x_2 \cdots x_p$ is given for large $M$ by
\[
P\left(x_1,\,\cdots x_p\right) =M^p \rho(x_1) \rho(x_2) \cdots \rho\left(x_p\right) \exp\left[-M\int_{x_p}^\infty \rho(x) dx \right]
\sim \frac{(M B)^p}{\left(x_1 \, x_2 \cdots x_p\right)^{1+\mu}} \exp\left[-\frac{M B}{\mu \, x_p^{\mu} }\right]
\]
We see that the two previous distributions are identical through the change of variables
\[
x_i= \left(\frac{M B \beta}{A}\right)^\frac{1}{\mu} e^{-\beta E_i} \quad \text{with} \quad \mu=\frac{\beta_c}{\beta}
\]
This shows that the partition function $Z(\beta)= \sum_i e^{-\beta E_i}$ is up to a rescaling a sum of i.i.d. random variables with an heavy tail.


\subsection{The \texorpdfstring{$p$}{p}-spin Ising spin glass}

The REM is the large $p$ limit of the $p$-spin model~\cite{Derrida_1980_Random,
Derrida_1981_Random} which is a generalisation of the Sherrington--Kirkpatrick model with an energy given by
\begin{equation}\label{p-spin}
E_{\bf A} (\bfS) = -\sum_{i_1\,\le\,i_2 \cdots i_p} A_{i_1\cdots i_p} \ S_{i_1} \ S_{i_2} \cdots S_{i_p}
\end{equation}
where the $p$-spin interactions $A_{i_1\cdots i_p}$ are quenched random variables distributed according to
\begin{equation}\label{p-spin1}
\rho\left(A_{i_1\cdots i_p}\right)= \sqrt{\frac{N^{p-1}}{\pi \, p!}} \exp\left[- \frac{\left(A_{i_1\cdots i_p}\right)^2 \, N^{p-1} }{ p! } \right]
\end{equation}
It is easy to see from \eqref{p-spin},\eqref{p-spin1} that the energies $E(\bfS)$ are still distributed according to~\eqref{B1} and that their covariances are given (for large $N$) by
\begin{equation}\label{p-spin2}
\left\langle E(\bfS) \, E({\bf S'})\right\rangle \simeq \frac{N}{2 } \big[q(\bfS,{\bf S'}) \big]^p
\end{equation}
Clearly this covariance vanishes in the large $p$ limit when $ |q(\bfS,{\bf S'})| < 1$ and it was argued in~\cite{Derrida_1980_Random,Derrida_1981_Random} that in the large $p$ limit, the $p$-spin model should be equivalent to the REM. Indeed, using the RSB approach Gross and Mézard~\cite{Gross_1984_simplest} were able to calculate the average free-energy as well as $\langle Y_2 \rangle $ in the large $p$ limit and their predictions agree with~\eqref{B3} and~\eqref{A7} with $\mu$ given by~\eqref{B6}. A RSB solution for finite $p$ was then developed by Gardner in 1985~\cite{Gardner_1985_Spin}.


\section{The general framework of a Poisson REM}\label{sect3}

In this section we give a few general formulas \eqref{C2}-\eqref{C6} allowing to calculate the moments of the partition and the average overlaps for a REM whose energies are the points of a Poisson process of an arbitrary density $\rho(E)$.

\subsection{An arbitrary density}

For a Poisson REM of density $\rho(E)$, the energies $E(\bfS)$ are the points of a Poisson process of density $\rho(E)$ meaning that in each infinitesimal interval $(E,E+dE)$ with $dE \ll 1$ there is one configuration with probability $\rho(E) dE$ and with probability $1 - \rho(E) dE$ that there is no configuration at that energy. (If the integral of $\rho(E)$ is infinite, then the number of configurations is infinite with probability 1. Then the partition function at some low value of $\beta$ i.e. at high enough temperature could become infinite. However, in all cases we consider here, only the lowest energy levels matter in the low temperature phase implying that the negative moments of the partition function are finite and the overlaps have non trivial values.)

For each realization of the process, one has a sequence of energies $E_1, E_2, \cdots E_p, \cdots$ and the partition function is by definition
\begin{equation}\label{AA1}
Z(\beta)= \sum_p e^{-\beta E_p}
\end{equation}
Given $\rho(E)$ one can first show that
\begin{equation}\label{C1}
\left\langle e^{-t Z(\beta)}\right\rangle = e^{\phi(t)}
\end{equation}
where
\begin{equation}\label{C2}
\phi(t)= \int \rho(E) dE \left(\exp\left[-t \, e^{-\beta E}\right] -1\right)
\end{equation}
To do so one simply writes that
\begin{align*}
\left\langle e^{- t Z(\beta)}\right\rangle & =
\prod_E \Big(1- \rho(E) dE +e^{-t e^{-\beta E}} \, \rho(E) dE \Big)
\end{align*}
and because $dE \ll 1$ one can exponentiate to obtain \eqref{C1},\eqref{C2}. Knowing $\phi(t)$, one can then write the expressions of integer and non integer moments of the partition function. For example for negative $n$ one has
\begin{equation}\label{C3}
\left\langle Z(\beta)^n\right\rangle = \frac{1}{\Gamma(-n)} \int_0^\infty t^{-n-1} dt \, e^{\phi(t)} \quad \text{for} \; n<0
\end{equation}
In the same manner (see the proof below) one can show (see~\eqref{B9}) that for $n < k$
\begin{equation}\label{C4}
\left\langle Y_k Z(\beta)^n\right\rangle = \left\langle Z(k \beta) \, Z(\beta)^{n-k} \right\rangle = \frac{1}{\Gamma(k-n)} \int_0^\infty t^{k-n-1} dt \, \Phi_k(t)\ e^{\phi(t)}
\end{equation}
or
\begin{equation}\label{C5}
\left\langle Y_{k,k'} Z(\beta)^n\right\rangle = \frac{1}{\Gamma(k+k'-n)} \int_0^\infty t^{k+k'-n-1} dt \, \Phi_k(t)\ \Phi_{k'}(t) \ e^{\phi(t)}
\end{equation}
where
\begin{equation}\label{C6}
\Phi_k(t) = \int \rho(E) dE \, \exp\left[-k \beta E -t \, e^{-\beta E}\right]
\end{equation}


\begin{proof}
For example to obtain~\eqref{C4} one can use (see\eqref{B9}) the fact that
\[
Y_k = \frac{Z(k \beta)}{Z(\beta)^k}
\]
and that
\[
\left\langle Z(k\beta) e^{-t Z(\beta)}\right\rangle = \sum_E \rho(E) dE e^{-k \beta E-t e^{-\beta E}} \
\prod_{E'\ne E} \left(1- \rho(E') dE' +e^{-t e^{-\beta E'}} \, \rho(E') dE' \right)
\]
Then one can replace the sum over $E$ by an integral and the product as the exponential of an integral (in fact the condition $E\neq E'$ can be forgotten). A similar reasoning leads to~\eqref{C5}.
\end{proof}

\subsection{The case of the exponential density~\eqref{B15}}

In this section we show how to recover \eqref{B4}, \eqref{B10}, \eqref{B11} using \eqref{C1}-\eqref{C6} when the density $\rho(E)$ is exponential. For the density given by~\eqref{B15}, the functions $\phi(t)$ and $\Phi_k(t)$ defined in~\eqref{C2} and~\eqref{C6} have explicit exact expressions in the low temperature phase ($\beta> \beta_c$).
\begin{equation}\label{C7}
\phi(t)=\frac{A}{\beta}\ \Gamma\left(-\frac{\beta_c}{\beta} \right) \, t^\frac{\beta_c}{\beta}
\; ; \quad
\Phi_k(t)=\frac{A}{\beta}\ \Gamma\left(k-\frac{\beta_c}{\beta} \right) \, t^{\frac{\beta_c}{\beta}-k}
\end{equation}
(In the high temperature phase, i.e. for $\beta< \beta_c$, it turns out that $\phi(t)$ is infinite for the density~\eqref{B15}. There it is no longer appropriate to replace the original REM by a Poisson process with the exponential density~\eqref{B15}.) Using~\eqref{C7} the integrals in \eqref{C3}-\eqref{C5} can be computed exactly. One finds
\begin{equation}\label{C8}
\left\langle Z(\beta)^n \right\rangle =\left(- \frac{A}{\beta}\, \Gamma\left(-\frac{\beta_c}{\beta} \right)\right)^{n \frac{\beta}{\beta_c}} \ \frac{\Gamma\left(1- n \frac{\beta}{\beta_c}\right)}{\Gamma(1-n)}
\end{equation}
which is identical to~\eqref{B4} with the definition~\eqref{B6} of $\mu$. One also recovers that way from \eqref{C4}, \eqref{C5} the expressions \eqref{B10}, \eqref{B11} of the overlaps.



\section{The double exponential and its consequences}\label{sect4}

One simple case for which the expressions~\eqref{B10} and~\eqref{B11} of the overlaps are no longer valid is when the density $\rho(E)$ is a sum of two exponentials. Here we show that one needs to let fluctuate the parameter $\mu$ in these expressions. The generalization to the sum of an arbitrary number of exponentials will be straightforward. This will allow us to calculate to recover, in a much easier way, earlier results on finite size corrections or on the effect of discrete energies

\subsection{The double exponential case}

In this case, the density $\rho(E)$ is
\begin{equation}\label{D1}
\rho(E) = A_1 \, e^{\beta_1 \, E} + A_2 \, e^{\beta_2 \, E}
\end{equation}
Then~\eqref{C2} and~\eqref{C6} become
\begin{equation}\label{D2}
\phi(t) = B_1 \, t^{\mu_1} \, \Gamma(-\mu_1) + B_2 \, t^{\mu_2} \, \Gamma(-\mu_2)
\end{equation}
\begin{equation}\label{D3}
\Phi_k(t) = B_1 \, t^{\mu_1-k} \, \Gamma(k-\mu_1) + B_2 \, t^{\mu_2-k} \, \Gamma(k-\mu_2)
\end{equation}
with
\begin{equation}\label{D4}
B_i= \frac{A_i}{\beta} \quad \text{and} \quad \mu_i=\frac{\beta_i}{\beta}
\end{equation}

The relation~\eqref{C3} (after an integration by parts) can be written as
\begin{equation}\label{D5}
\left\langle Z(\beta)^n\right\rangle = \frac{-1}{\Gamma(1-n)} \, \int_0^\infty t^{-n} dt \, \phi'(t) \, e^{\phi(t)}
\end{equation}
Then by replacing $\phi'(t)$ using~\eqref{D2} and using again~\eqref{C3}, one gets for $n <0$
\[
\left\langle Z(\beta)^n \right\rangle = B_1 \, \frac{\Gamma(1-\mu_1)\, \Gamma(\mu_1-n)
}{\Gamma(1-n)}
\, \left\langle Z(\beta)^{n-\mu_1}\right\rangle + B_2 \, \frac{\Gamma(1-\mu_2)\, \Gamma(\mu_2-n)
}{\Gamma(1-n)}
\, \left\langle Z(\beta)^{n-\mu_2} \right\rangle
\]
which can be written as
\begin{equation}\label{D6}
\left\langle Z(\beta)^n \right\rangle = \sum_{i=1,2} B_i \ \frac{\Gamma(1-\mu_i)\, \Gamma(\mu_i-n)
}{\Gamma(1-n)}
\, \left\langle Z(\beta)^{n-\mu_i}\right\rangle
\end{equation}

For the overlaps, starting from~\eqref{C4} with $\Phi_k(t)$ given by~\eqref{D3} one can repeat the same procedure to get
\begin{equation}\label{D6a}
\left\langle Y_k\ Z(\beta)^n \right\rangle =
\int_0^\infty t^{k-n-1} dt
\left[\sum_{i=1,2} B_i \ \frac{\Gamma(k-\mu_i)}{\Gamma(k-n)}
\, t^{\mu_i-k}
\, e^{\phi(t)} \right]
\end{equation}
which gives
\begin{equation}\label{D7}
\left\langle Y_k\ Z(\beta)^n \right\rangle = \sum_{i=1,2} B_i \, \frac{ \Gamma(k-\mu_i)
\ \Gamma(\mu_i-n)
}{\Gamma(k-n)}
\, \left\langle Z(\beta)^{n-\mu_i}\right\rangle
\end{equation}
Using~\eqref{D6} one can rewrite~\eqref{D7} as
\begin{equation}\label{D8}
\frac{\left\langle Y_k\ Z(\beta)^n \right\rangle}{\left\langle Z(\beta)^n \right\rangle} =
\sum_{i=1,2} \frac{
\Gamma(k-\mu_i)}{
\Gamma(1-\mu_i)} \ \frac{\Gamma(1-n)}{\Gamma(k-n) }\ W_i
\end{equation}
where the weights $W_i$ are given by
\begin{equation}\label{D9}
W_i = B_i \ \frac{\Gamma(1-\mu_i)\, \Gamma(\mu_i-n)}{\Gamma(1-n)}
\,
\frac{\left\langle Z(\beta)^{n-\mu_i}\right\rangle}
{\left\langle Z(\beta)^n\right\rangle }
\end{equation}
We see that the average overlaps in~\eqref{D8} have the same expression as in~\eqref{B10} except that now, for the double exponential~\eqref{D1}, the parameter $\mu$ fluctuates between two values $\mu_1$ or $\mu_2$.

The averages of other overlap functions can be derived from~\eqref{C5} in the same way: for example using the fact that by iterating~\eqref{D6} one has
\begin{equation}\label{D10}
\left\langle Z(\beta)^n\right\rangle =
\sum_{i=1,2} \sum_{j=1,2} B_i \, B_j \, \frac{ \Gamma(1-\mu_i) \Gamma(1-\mu_j) \, \Gamma(\mu_i+ \mu_j-n)}{(\mu_i-n) \, \Gamma(1-n)} \left\langle Z(\beta)^{n-\mu_i-\mu_j}\right\rangle
\end{equation}
and one gets, using~\eqref{D3} twice in~\eqref{C5}
\begin{equation}\label{D11}
\frac{ \left\langle Y_{k,k'} Z(\beta)^n \right\rangle}{\left\langle Z(\beta)^n \right\rangle} =
\sum_{i=1,2}
\sum_{j=1,2} (\mu_i-n) \frac{ \Gamma(k-\mu_i)}{\Gamma(1-\mu_i) } \frac{ \Gamma\left(k'-\mu_j\right)}{\Gamma\left(1-\mu_j\right) } \frac{\Gamma(1-n)}{\Gamma(k+k'-n)} W_{i,j}
\end{equation}
with the weights $W_{i,j}$ given by
\begin{equation}\label{D12}
W_{i,j}= B_i \, B_j \, \frac{ \Gamma(1-\mu_i) \Gamma(1-\mu_j) \, \Gamma(\mu_i+ \mu_j-n)}{(\mu_i-n) \, \Gamma(1-n)} \frac{ \left\langle Z(\beta)^{n-\mu_i-\mu_j}\right\rangle}{\left\langle Z(\beta)^n \right\rangle}
\end{equation}
Clearly~\eqref{D11} is a generalization of~\eqref{B11} where the $\mu$'s fluctuate. We see in~\eqref{D12} that the $\mu_i$'s are correlated (as $W_{i,j} \ne W_i \, W_j$). Note that although the $W_{i,j}$ is not symmetric under the exchange $\mu_i \leftrightarrow \mu_j$, the symmetry is restored in the sum~\eqref{D11}.



\begin{rema}
One possible realization of the double exponential density~\eqref{D1} would be to consider a REM with a total of $\alpha_1^N + \alpha_2^N $ configurations, the first $\alpha_1^N$ configurations having energies distributed according to $P_1(E) = \exp[-E^2 /(N a_1)] / \sqrt{\pi a_1}$ and the last $\alpha_2^N$ according to $P_2(E) = \exp[-E^2 /(N a_2)] / \sqrt{\pi a_2}$. When the parameters $\alpha_1,\alpha_2,a_1,a_2$ satisfy the following condition
\[
\sqrt{a_1 \log \alpha_1}= \sqrt{a_2 \log \alpha_2} \equiv -\epsilon
\]
the density
$\rho(E)$ near the ground state energy becomes
\[
\rho(E) =
\frac{1}{\sqrt{\pi N a_1} } \exp\left[2 \sqrt\frac{\log \alpha_1}{a_1}(E- N \epsilon) \right] + \frac{1}{\sqrt{\pi N a_2} } \exp\left[2 \sqrt\frac{\log \alpha_2}{a_2}(E- N \epsilon) \right]
\]
which is indeed of the form~\eqref{D1}. For more general spin glass models, one can imagine that the double exponential density~\eqref{D1} could also be relevant near a first order transition between two low temperature phases that exhibit a one step RSB.
\end{rema}

\subsection{The replica approach in the case of the double exponential}

In the appendix A, we recall how to obtain the expressions \eqref{A7}, \eqref{B10}, \eqref{B11} using the Parisi matrix shown in Figure~\ref{peter-fig1}. We are now going to see that in order to recover the above expressions \eqref{D8}, \eqref{D11} one needs to consider matrices of the form shown in Figure~\ref{peter-fig2} and to average over matrices of this shape by letting the number $n_i$ of blocks of size $\mu_i$ to fluctuate.


\begin{figure}[!h]
\centerline{\includegraphics[width=6.cm]{peter-fig2}}
\caption{In the case of the double exponential, the overlap matrix of Figure~\ref{peter-fig1} shoud be replaced by a matrix whose block sizes fluctuate.}
\label{peter-fig2}
\end{figure}

Let us imagine that a matrix of overlaps is a $n \times n$ matrix, as in Figure~\ref{peter-fig2}, where the diagonal blocks have varying sizes $\nu_\alpha$ so that
\begin{equation}\label{D13}
n= \sum_\alpha \nu_\alpha
\end{equation}
The probability that $k$ different replicas chosen among $n$ replicas belong to the same block is
\begin{equation}\label{D14}
Y_k= \sum_{\alpha} \frac{\nu_\alpha !}{(\nu_\alpha-k) !} \, \frac{(n-k)!}{n!}
\end{equation}
These expressions depend on the number $n$ of replicas and on the sizes $\nu_\alpha$ of the blocks. Now if we allow $n$ and the $\nu_\alpha$'s to take non-integer values, as in the original RSB scheme~\cite{Parisi_1979_Infinite} this can be written as
\begin{equation}\label{D15}
Y_k= \sum_{\alpha} \frac{\Gamma(k-\nu_\alpha)}{\Gamma(-\nu_\alpha)} \, \frac{\Gamma(-n)}{\Gamma(k-n)}
\end{equation}
Similarly
\begin{equation}\label{D16}
Y_{k,k'}= \sum_{\alpha} \sum_{\alpha' \neq \alpha} \frac{\Gamma(k-\nu_\alpha)}{\Gamma(-\nu_\alpha) } \, \frac{\Gamma(k'-\nu_{\alpha'})}{\Gamma(-\nu_{\alpha'})} \, \frac{\Gamma(-n)}{\Gamma(k+k'-n)}
\end{equation}



Let us imagine that the $\nu_\alpha$ can take possible values $\mu_i$ (in the case of two exponentials there are two possible values $\mu_1$ and $\mu_2$ defined in~\eqref{D4}) and let $n_i$ be the number of $\nu_\alpha$'s taking the value $\mu_i$. Then~\eqref{D13} becomes
\begin{equation}\label{D16a}
\sum_i n_i \, \mu_i=n
\end{equation}
and going from one matrix to the other the values of the $n_i$ fluctuate (while the $\mu_i$'s don't) keeping this sum constant.

Averaging over the $n_i$'s one gets
\[
\frac{\left\langle Y_k \,
Z(\beta)^n \right\rangle}{\left\langle Z(\beta)^n\right\rangle } = \sum_i \langle n_i \rangle \, \frac{\Gamma(k-\mu_i)}{\Gamma(-\mu_i) } \, \frac{\Gamma(-n)}{\Gamma(k-n)} = \sum_i \frac{\langle n_i \rangle \, \mu_i}{n} \ \frac{\Gamma(k-\mu_i)}{\Gamma(1-\mu_i) } \, \frac{\Gamma(1-n)}{\Gamma(k-n)}
\]
which has the form~\eqref{D8} if one identifies
\begin{equation}\label{D17}
W_i =\frac{\langle n_i \,\rangle \mu_i}{n}
\end{equation}
One can interpret the weight $W_i$ as the probability that one among $n$ replicas belongs to a block of size $\mu_i$.

The matrix in Figure~\ref{peter-fig2} is characterized by the numbers $n_1$ and $n_2$ of blocks of size $\mu_1$ and $\mu_2$. The sizes $\mu_1$ and $\mu_2$ of the blocks are fixed. Only $n_1$ and $n_2$ fluctuate with the constraint~\eqref{D16a} and to recover \eqref{D8}, \eqref{D11} one needs to average over $n_1$ and $n_2$. Note that in this picture, although $\mu_1$ and $\mu_2$ are simply given by~\eqref{D4}, the expressions of the weights $W_i$ and $W_{i,j}$ in \eqref{D9}, \eqref{D12} are less explicit as they require the knowledge of negative moments of the partition function.


Similarly
\begin{multline}\label{D17a}
\frac{\left\langle Y_{k,k'} \, Z(\beta)^n \right\rangle}{\left\langle Z(\beta)^n \right\rangle }\\
=\left[\sum_i \left\langle n_i^2- n_i\right\rangle \ \frac{\Gamma(k-\mu_i)}{\Gamma(-\mu_i) } \frac{\Gamma(k'-\mu_i)}{\Gamma(-\mu_i) } + \sum_i \sum_{j \neq i} \left\langle n_i n_j\right\rangle \ \frac{\Gamma(k-\mu_i)}{\Gamma(-\mu_i) } \frac{\Gamma(k'-\mu_j)}{\Gamma(-\mu_j) } \,
\right] \frac{\Gamma(-n)}{\Gamma(k+k'-n)}
\end{multline}
which has the form~\eqref{D11} if one identifies
\begin{equation}\label{D18}
W_{i,i} =\frac{\langle n_i (n_i-1)\rangle\, \mu_i^2}{n(n-\mu_i)} \quad \text{ and for } i\neq j \quad W_{i,j} =\frac{\langle n_i n_j\rangle \, \mu_i \mu_j}{n(n-\mu_i)}
\end{equation}

In all cases i.e. whether $i=j$ or $i\neq j$, one can interpret $W_{i,j}$ as the probability that choosing 2 replicas belonging to different blocks, the first one is in a block of size $\mu_i$ and the second one in a block of size $\mu_j$.


\subsection{The finite size corrections to the REM}

It is straightforward to generalize the above results for the double exponential~\eqref{D1} to the sum of an arbitrary number of exponentials
\begin{equation}\label{D19}
\rho(E) =\sum_j A_j \ e^{\beta_j E}
\end{equation}
The density of energies~\eqref{B1} of the REM near the ground state energy can be written as
\begin{equation}\label{D20}
\frac{2^N}{\sqrt{\pi N}} \exp \left[-\frac{E^2}{N} \right] = \int \frac{dy}{\pi \, N } e^{- N y^2} \exp\left[\left(\beta_c + 2 i y\right)\left(E+ N \sqrt{\log 2}\right)\right] 
\end{equation}
which is of the form~\eqref{D19} if one allows the $\beta_j$'s to take small imaginary values $\beta_c + 2 iy $. Therefore as in the double exponential case, this implies that the block sizes $\mu_j $ (see~\eqref{D4}) can become complex.

Because the possible values of the block sizes $\mu_j= \frac{\beta_c + 2 i y}{\beta}$ are now complex (see~\eqref{D4}), the weights $W_i$ and the numbers $n_i$ may be complex too. Using the fact that $y$ is distributed as in~\eqref{D20}, one gets that for large $N$
\[
\left\langle \mu_j^2\right\rangle - \langle \mu_j \rangle^2 \sim -\frac{2}{N \beta^2}
\]
recovering the prediction in equation (55) of~\cite{Derrida_2015_Finite} that the block size have a negative variance. (Note that the variance is negative simply because the box sizes $\mu_i$ are complex).

\subsection{The case of discrete energies}

The case where the energies take only discrete values can also be treated the same way, for example if the energies $E(\bfS)$ were sums of $N$ random numbers taking only integer values~\cite{Derrida_2022_discrete,Ogure_2005_exact,
Bouchaud_2003_Energy}. In that case the density $\rho(E)$ near the ground state energy could be written as
\begin{equation}\label{D21}
\rho(E) =A \ e^{\beta_c E} \sum_{n=-\infty}^\infty \delta(E-n)= \sum_{p=-\infty}^\infty A \ e^{\left(\beta_c + 2 i \pi p\right)E}
\end{equation}

In this case too, the density $\rho(E)$ is a particular case of~\eqref{D19}. Therefore one can use the above formulas allowing the $\beta_j$ to take complex values $\beta_c + 2 i \pi p$. The block sizes $\mu_i$'s then take complex values. Therefore according to \eqref{D9}, \eqref{D17} the replica interpretation is that the ``probabilities'' or rather the $W_p$ (which do not need here to be positive or even real) that a replica belongs to a block of size $\mu_p=\frac{\beta_c +2i \pi p}{\beta}$ are given by
\begin{equation}\label{D22}
W_p = \frac{A}{\mu_p} \ \frac{\Gamma(1-\mu_p)\, \Gamma(\mu_p-n)
}{\Gamma(1-n)}
\, \frac{\left\langle Z(\beta)^{n-\mu_p}\right\rangle
}{\langle Z(\beta)^n \rangle }
\quad \text{with} \quad \mu_p=\frac{\beta_c +2i \pi p}{\beta}
\end{equation}

One way of illustrating that the statistics of the overlaps are strongly modified is to look at the distribution $\Pi_\mu(Y)$ drawn in Figure~\ref{peter-fig4} and to compare it to the shape of Figure~\ref{figure1} for continuous energies.

\begin{figure}[h]
\centering
\includegraphics[width=10.5cm]{Figure-3-bitmap}
\vspace*{-10pt}
\caption{The distribution $\Pi_\mu(Y)$ when the energies take only integer values. Here like in Figure~\ref{figure1}, $\beta=3 \beta_c$. }\label{peter-fig4}
\end{figure}

It is easy to see that changing the amplitude $A$ in~\eqref{D21} to $A e^{\beta_c}$ has the effect of shifting all the energies by $1$ implying that $\Pi_\mu$ and the overlaps remain unchanged. Other changes of $A$ (for example~$A \to A e^{\beta_c x }$ when $x$ is not an integer change the shape of the distribution $\Pi_\mu (Y)$ and the average overlaps $Y_k$ with a periodic dependence on the parameter $x$ (as observed in~\cite{Derrida_2022_discrete}).

\section{The overlaps of the REM at two temperatures}\label{sect5}

Considering the overlaps between configurations at different temperatures in the same energy landscape is a way to study temperature chaos~\cite{Bray_1987_Chaotic} (for a review see~\cite{Rizzo_2009_Chaos}). It is well known that the REM as well as other models, such as directed polymers in a random medium (in their tree version), do not exhibit temperature chaos~\cite{Sales_2001_Rejuvenation,Krzakala_2002_Chaotic,Kurkova_2003_Temperature,Pain_2021_Two}. Still it is interesting to understand how the replica approach works to predict these overlaps. A detailed analysis of this question is given in~\cite{Derrida_2021_One}. In this section we provide an alternative analysis extending the replica-free calculation of the previous section to the two temperature case and obtaining an expression for the overlaps that gives an explicit formula in terms of the probabilities of fluctuating block sizes. As in the case of finite size corrections or of discrete energies (see Section~\ref{sect4}) we will see that, in the two temperature case, the block sizes take complex values.

Let us consider the overlaps between $k$ configurations at inverse temperature $ \beta$ and $k'$ configurations at inverse temperature $\beta'$. For a given sample, the probability $U_{k,k'}$ that all these $k+k'$ configurations occupy the same energy level is
\begin{equation}\label{F0}
U_{k,k'}= \frac{\sum_\bfS e^{-\left(k\beta +k'\beta'\right) E(\bfS)}}{Z(\beta)^k \ Z(\beta')^{k'}} =\frac{Z\left(k\beta +k'\beta'\right)}{Z(\beta)^k \ Z(\beta')^{k'}}
\end{equation}
and our goal is to determine the following sample averages $\langle Z(\beta)^n Z(\beta')^{n'} \rangle$ and $\langle U_{k,k'} Z(\beta)^n Z(\beta')^{n'}\rangle $ in the low temperature phase. Although there is a perfect symmetry under the exchange $(\beta,k,n) \leftrightarrow (\beta',k',n')$, this symmetry is not apparent in some expressions below, but in fact it is indeed respected.


For negative $n$ and $n'$, one can repeat what we did in \eqref{C1}-\eqref{C3} and write
\begin{equation}\label{F1}
\left\langle Z(\beta)^n \, Z(\beta')^{n'} \right\rangle = \frac{1}{\Gamma(-n) \, \Gamma(-n')} \int_0^\infty t^{-n-1} dt \int_0^\infty u^{-n'-1} du \ e^{\varphi(t,u)}
\end{equation}
where
\begin{equation}\label{F2}
\varphi(t,u) = \int \rho(E) dE \left(\exp\left[-t e^{-\beta E} - u e^{-\beta' E} \right] -1 \right)
\end{equation}
One can rewrite~\eqref{F1} as
\begin{equation}\label{F3}
\left\langle Z(\beta)^n \, Z(\beta')^{n'} \right\rangle = \frac{-1}{\Gamma(1-n) \, \Gamma(-n')} \int_0^\infty t^{-n} dt \int_0^\infty u^{-n'-1} du \ e^{\varphi(t,u)} \ \frac{d \varphi(t,u)}{dt}
\end{equation}
We are now going to use the following complex integral representation (called the Cahen--Mellin representation, see for example~\cite{Paris_2001_Asymptotics}) of the exponential
\begin{equation}\label{F4}
e^{-t} = \int_{-\infty}^\infty \frac{dy}{2 \pi} \ \Gamma(\nu+i y) \ t^{-\nu-i y} 
\end{equation}
where $\nu$ is real and $ \nu > 0$. (In~\eqref{F4} the integral does not depend on $\nu$ as long as $\nu >0$ simply because the integral path can be deformed without crossing any singulatity.)

For the exponential density~\eqref{B15} and transforming $\exp[- t e^{-\beta E}]$ according to~\eqref{F4} one gets, (note that in the low temperature phase $\beta>\beta_c$ and $\beta'>\beta_c$)
\[
\frac{d \varphi(t,u)}{dt}
= -A \int dE e^{(\beta_c -\beta)E}
\exp\left[-u e^{-\beta' E} \right] \int_{-\infty}^\infty \frac{dy}{2 \pi} \ \Gamma(\nu+i y) \ t^{-\nu-i y} e^{(\nu + i y)\beta E}
\]
which, after integrating over $E$, becomes
\begin{equation}\label{F5}
\frac{d \varphi(t,u)}{dt} = - \frac{A}{\beta'}
\int_{-\infty}^\infty \frac{dy}{2 \pi} \ \Gamma(\nu+i y) \
\Gamma\left(\frac{\beta(1-\nu-i y) -\beta_c}{\beta'}\right) t^{-\nu-i y} \ u^\frac{ \beta_c -\beta (1-\nu - i y)}{\beta'}
\end{equation}
Then one gets from~\eqref{F3} using~\eqref{F1}
\begin{multline*}
\left\langle Z(\beta)^n \, Z(\beta')^{n'} \right\rangle = \frac{A}{\beta'}
\int_{-\infty}^\infty \frac{dy}{2 \pi} \ \frac{\Gamma(\nu+i y) \
\Gamma\left(\frac{\beta(1-\nu-i y) -\beta_c}{\beta'}\right)
}{\Gamma(1-n) \, \Gamma(-n')}
\Gamma(1-n-\nu - i y) \\
\times
\Gamma\left({-n'+\frac{\beta_c -\beta(1-\nu-i y}{\beta'}}\right)
\left\langle Z(\beta)^{n-1+\nu+i y} Z(\beta')^{n'-\frac{\beta_c}{\beta'} +\frac{\beta(1-\nu-iy)}{\beta'}}
\right\rangle
\end{multline*}
in other words
\begin{multline}\label{F6}
\left\langle Z(\beta)^n \, Z(\beta')^{n'} \right\rangle\\
 = \frac{A}{\beta'}
\int_{-\infty}^\infty \frac{dy}{2 \pi} \ \frac{\Gamma(1-\mu_y) \
\Gamma\left(-\mu'_y\right)
}{\Gamma(1-n) \, \Gamma(-n')}
\Gamma(\mu_y-n)
\Gamma\left(\mu'_y-n'\right)
\left\langle Z(\beta)^{n-\mu_y} Z(\beta')^{n'-\mu'_y}
\right\rangle
\end{multline}
where
\begin{equation}\label{F7}
\mu_y= 1-\nu - i y
\quad \text{and} \quad
\mu'_y = \frac{\beta_c- \beta \mu_y}{\beta'}
\end{equation}
(We have checked~\eqref{F6} against the explicit expression in~\cite{Derrida_2021_One} and we have also checked that when $\beta=\beta'$ we recover~\eqref{C8}.) Comparing with~\eqref{D6} we see that here the $n$ replicas at inverse temperature $\beta$ and they are grouped in blocks of sizes $\mu_y$ and the $n'$ at inverse temperature $\beta'$ in blocks of sizes $\mu'_y$. Each block of size $\mu_y$ is associated to a block of size $\mu'_y$ which satisfy
\begin{equation}\label{F8}
\beta \mu_y + \beta' \mu'_y= \beta_c
\end{equation}
So the $\mu_y$ and $\mu'_y$ fluctuate but the relation~\eqref{F8} remains fixed. The imaginary part of the block sizes fluctuates between $-i \,\infty$ and $+i \,\infty$, but the real part can be arbitrarily fixed within the constraints
\begin{equation}\label{F9}
n<\Re(\mu_y)<1 \;, \; n'<\Re(\mu'_y)<0
\end{equation}
which is a result of the requirement for $\nu>0$ in~\eqref{F4} and that the gamma functions in~\eqref{F6} must have positive arguments. Combining~\eqref{F8} and~\eqref{F9} we obtain
\begin{equation}\label{F10}
\frac{\beta_c}{\beta} < \Re(\mu_y) < \min \left(1, \frac{\beta_c-\beta' n'}{\beta}\right).
\end{equation}
(Remember that $n' <0$).

We now turn to the two temperature overlaps defined in~\eqref{F0}. With negative $n,n'$ and positive integer $k,k'$ we can use the same approach as~\eqref{C4} to obtain
\begin{equation}\label{F11}
\left\langle U_{k,k'} \, Z(\beta)^n \, Z(\beta')^{n'}\right\rangle = \frac{1}{\Gamma(k-n) \, \Gamma(k'-n')} \int_0^\infty t^{k-n-1} dt \int_0^\infty u^{k'-n'-1} du \, \Phi_{k,k'}(t,u) \ e^{\varphi(t,u)}
\end{equation}
where
\begin{equation}\label{F12}
\Phi_{k,k'}(t,u) = \int \rho(E) dE \exp\left[-\beta k E -\beta' k' E -t e^{-\beta E} - u e^{-\beta' E} \right].
\end{equation}
Using the Cahen--Mellin representation~\eqref{F4} we can write this as a contour integral and for exponential disorder~\eqref{B15} this simplifies to
\begin{equation}\label{F13}
\Phi_{k,k'}(t,u) = \frac{A}{\beta'}
\int_{-\infty}^\infty \frac{dy}{2 \pi} \ \Gamma(\nu+i y) \
\Gamma\left(k'+ \frac{\beta(k-\nu-i y) -\beta_c}{\beta'}\right) t^{-\nu-i y} \ u^{ -k'- \frac{\beta(k-\nu-i y) -\beta_c}{\beta'}}.
\end{equation}
Substituting in~\eqref{F11} and using~\eqref{F1} gives
\begin{multline}\label{F15}
\left\langle U_{k,k'} \, Z(\beta)^n \, Z(\beta')^{n'} \right\rangle\\
= \frac{A}{\beta'}
\int_{-\infty}^\infty \frac{dy}{2 \pi} \
\frac{\Gamma(k-\mu_y) \ \Gamma\left(k'-\mu'_y\right)}{\Gamma(k-n) \ \Gamma(k'-n')} 
\times \Gamma(\mu_y-n) \ \Gamma\left(\mu'_y -n'\right)
\left\langle Z(\beta)^{n-\mu_y} Z(\beta')^{n'-\mu'_y} \right\rangle
\end{multline}
where $\mu_y$ and $\mu'_y$ are defined in~\eqref{F7} and must satisfy~\eqref{F8} and~\eqref{F10}.

We can express the overlaps in terms of a normalised weight function $W(y)$ similar in concept to formula~\eqref{D11}
\begin{equation}\label{F17}
\frac{\left\langle U_{k,k'} \, Z(\beta)^n \, Z(\beta')^{n'}\right\rangle}{\left\langle Z(\beta)^n \, Z(\beta')^{n'} \right\rangle}=
\int_{-\infty}^\infty dy \
\frac{\Gamma(k-\mu_y) \ \Gamma\left(k'-\mu'_y\right)}{\Gamma(1-\mu_y) \ \Gamma\left(-\mu'_y\right)}
\frac{\Gamma(1-n) \ \Gamma(-n')}{\Gamma(k-n) \ \Gamma(k'-n')} \, W(y)
\end{equation}
with
\begin{equation}\label{F18}
W(y)= \frac{A}{2 \pi \,\beta'}\frac{\Gamma(1-\mu_y) \
\Gamma\left(-\mu'_y\right)
}{\Gamma(1-n) \, \Gamma(-n')}
\Gamma(\mu_y-n)
\Gamma\left(\mu'_y-n'\right)
\frac{\left\langle Z(\beta)^{n-\mu_y} Z(\beta')^{n'-\mu'_y} \right\rangle}{ \langle Z(\beta)^n \, Z(\beta')^{n'} \rangle}.
\end{equation}

The weights $W(y)$ can be thought of as the probability that, for a given block, the block size is $\mu_y$ at inverse temperature $\beta$ and $\mu'_y$ at inverse temperature $\beta'$ subject to the condition~\eqref{F8}. However, they are complex and there is some freedom, due the choice of $\nu$ in~\eqref{F4}, in the choice of where the contour crosses the real axis. Despite this lack of uniqueness we expect that the different choices of weights will lead to the same moments for $\mu_y$ and $\mu'_y$.

The above expressions \eqref{F6}-\eqref{F10} are clearly non-symmetric under the exchange $(n,\beta,\mu_y) \leftrightarrow (n',\beta',\mu_y')$. We believe that very much like for the choice of the real part of $\mu_y$, there are many different choices of the distribution of the $\mu_y$ and of the $\mu_y'$ which lead to exactly the same predictions for the moments $\langle Z^n(\beta) Z^{n'}(\beta') \rangle $ or of the average overlaps even though the expressions look different.

\begin{figure}[h]
\centerline{\includegraphics[width=9cm]{replic-bis.pdf}}
\caption{The Parisi overlap matrix at two temperatures. The single temperature form of Figure~\ref{peter-fig1} must be replaced by an $(n+n')\times(n+n')$ matrix which has fluctuating block sizes. In addition it has block diagonal parts that represent overlaps between replicas at the same temperature and off diagonal parts that represent overlaps between the two different temperatures. }
\label{replic-bis}
\end{figure}

It is difficult to visualise block sizes with real and imaginary parts: in Figure~\ref{replic-bis} one should imagine that the fluctuations of the block sizes represent the fluctuations of their imaginary part while the real parts of these block sizes can remain fixed.

\section{The case of a general density \texorpdfstring{$\rho(E)$}{} and the replica approach}\label{sect6}

We consider now the case of a Poisson REM with a general density $\rho(E)$. (For technical reasons, we will assume that $\int \rho(E) dE= \infty$ to ensure that $Z(\beta)\neq 0$ with probability $1$. We will also assume that the average partition function $\langle Z(\beta) \rangle$ is finite at least for large enough $\beta$ and of the form $\langle Z(\beta) \rangle = e^{N \Psi(\beta)}$ in contrast to the exponential densities we considered before for which $\langle Z(\beta) \rangle =\infty$.

The replica approach consists usually in computing various quantities such as moments of the partition function or overlaps for an integer number $n$ of replicas and then to try to extend these expressions to non-integer values of $n$. Here we will see that both for an integer and a non-integer number $n$ of replicas, one obtains very similar expressions to those obtained in the previous section for the double exponential density.

\subsection{The replica approach for integer numbers of replicas}

For integer values of the number of replicas one can establish the following relations (see Appendix~\ref{appB}) which are consequences of \eqref{C1}-\eqref{C6}
\begin{itemize}
\item For integer $n >0$
\begin{equation}\label{E1}
\left\langle Z(\beta)^n \right\rangle = \sum_{\mu =1}^n \frac{(n-1)!}{(\mu-1)! \, (n-\mu)!} \ \langle Z(\mu \, \beta) \rangle \left\langle Z(\beta)^{n-\mu} \right\rangle
\end{equation}
\item and integer $n \ge k \ge 1 $
\begin{equation}\label{E2}
\left\langle Y_k\, Z(\beta)^n\right\rangle = \sum_{\mu =1}^n \frac{(n-k)!}{(\mu-k)! \, (n -\mu)!} \ \langle Z(\mu \, \beta) \rangle \left\langle Z(\beta)^{n-\mu}\right\rangle
\end{equation}
In~\eqref{E2}, as well as everywhere else in this paper, we use the convention that $(-n)!=\infty$ when $n$ is a positive integer.
\item This approach can be extended to obtain similar, but more complicated expressions for $\langle Y_{k,k'} Z(\beta)^n \rangle$.
\end{itemize}
One can notice the similarities between~\eqref{E1}, \eqref{E2} and~\eqref{D6}, \eqref{D7}. Then as in~\eqref{D8}, \eqref{D9}, \eqref{D11}, \eqref{D12} we see, by comparing~\eqref{E1} and~\eqref{E2}, that one can write
\begin{equation}\label{E5}
\frac{\left\langle Y_k\, Z(\beta)^n\right\rangle
}{\left\langle Z(\beta)^n\right\rangle }= \sum_{\mu=1}^n \frac{(n-k)! \ (\mu-1)!}{(n-1)! \, (\mu-k)! }
\ W(\mu)
\end{equation}
where the weights $W(\mu)$ are given by
\begin{equation}\label{E6}
W(\mu)= \frac{(n-1)!}{(\mu-1)! \, (n-\mu)!} \ \frac{\langle Z(\mu \, \beta) \rangle \ \left\langle Z(\beta)^{n-\mu} \right\rangle
}{\left\langle Z(\beta)^n\right\rangle}
\end{equation}
A replica interpretation of $W(\mu)$ is that if we choose one replica at random out of the $n$ replicas the probability that it is in a block of size $\mu$ is $W(\mu)$. The combinatorial factor in~\eqref{E5} is then given by the probability that the remaining $k-1$ replicas are also chosen from the same block (see Appendix~\ref{appA}). Already, at the level of an integer number $n$ of replicas, we see that, as in \eqref{D8}, \eqref{D11}, the $\mu's$ can fluctuate.

%and from the expression for $\langle Y_{k,k'} \, Z(\beta)^n \rangle$ we have found that they are correlated as they were in the replica calculation~\eqref{D12} above for the double exponential.


\subsection{When the number \texorpdfstring{$n$}{n} of replicas is not an integer}

We are now going to write \emph{exact} expressions which generalize \eqref{E1}-\eqref{E2} to non-integer values of $n$. As some expressions may change depending on the range of $n$, we will consider only the case $n<0$. As shown in the Appendix B, our derivation is based on the Cahen-Mellin representation~\eqref{F4} of the exponential.
\begin{itemize}
\item one gets for $n<0$ and $n<\mu < 1$
\begin{equation}\label{E10}
\left\langle Z(\beta)^n\right\rangle =
\int \frac{dy}{2\pi} \frac{\Gamma(1- \mu - i y) \ \Gamma(\mu-n + i y)}{\Gamma(1-n)} \left\langle Z\Big(\left(\mu + i y\right)\beta \Big) \right\rangle \left\langle Z(\beta)^{n-\mu-i y}\right\rangle
\end{equation}
\item For $ n < \mu<k$
\begin{equation}\label{E11}
\left\langle Y_k \, Z(\beta)^n\right\rangle =
\int \frac{dy}{2\pi} \frac{\Gamma(k- \mu - i y) \ \Gamma(\mu-n + i y)}{\Gamma(k-n)} \left\langle Z\Big(\left(\mu + i y\right)\beta \Big) \right\rangle \left\langle Z(\beta)^{n-\mu-i y} \right\rangle
\end{equation}
\end{itemize}

One can notice the similarity between~\eqref{E1}, \eqref{E2} and~\eqref{E10}, \eqref{E11}: the ratios of factorials become ratios of Gamma functions and the sums over $ \mu$ become complex integrals over~$y$.


For large $N$ (here $N$ is the system size in, for example, \eqref{B13}), the integrals~\eqref{E10}-\eqref{E11} are dominated by a single saddle point at $y=0$ and $\mu=\beta_c/\beta$ when the distribution of energies is of the form~\eqref{B13} and one recovers~\eqref{B10}, \eqref{B11}. One can in principe calculate all the finite size corrections by analysing the neighbourhood of this saddle point. On the other hand, when the energies take only integer values, there are several saddle points at the same height which contribute at $\mu=(\beta_c+ 2\pi i p)/\beta$ and $\mu'=(\beta_c+ 2\pi i p')/\beta$. By analogy with the case of double exponential of Section~\ref{sect4}, one can say that the sizes of the blocks in the Parisi matrix fluctuate and take complex values. In fact it is easy to see that the prefactors corresponding to these different saddle points in~\eqref{E10}, \eqref{E11}  coincide with the weights $W_p$ in~\eqref{D22}.

\section{Conclusion}

In this paper we have discussed a number of examples of systems exhibiting one step replica symmetry breaking but for which the original Parisi scheme has to be modified. The simplest case is a Poisson REM where the density of energies is a sum of two exponentials (see Section~\ref{sect4}). Our approach is based on developing a recursion relation between the negative moments of the partition function (see for example~\eqref{D6}). The overlaps can also be expressed in terms of the same negative moments of the partition function~\eqref{D7}. Then taking the ratio we obtain expressions \eqref{D8},\eqref{D9} for the overlaps as a weighted sum of two overlap expressions of the form~\eqref{B10} for the original REM. We have shown that the replica interpretation of this is that the block size of the matrix of overlaps is no longer fixed but each block can choose one of two possible block sizes, $\mu_1$ and $\mu_2$. The Parisi matrix is still block diagonal, but the block sizes fluctuate, taking either size $\mu_1$ or size $\mu_{2}$ as illustrated in Figure~\ref{peter-fig2}. The weights $W_i$ in formula~\eqref{D9} can be interpreted as the probability of a replica chosen at random being in a block of size $\mu_i$.

This approach is easily extended to other densities of energies that are sums of an arbitrary number of exponentials. This allows us to address REMs that have a disorder distribution that can be represented as a sum of exponentials. An example is the case of finite-size corrections to the REM where we show that the block sizes fluctuate and take complex values. As a consequence, the variance of the block sizes is negative and the corresponding weights are complex. Another example is the REM with discrete energies. Again in the replica approach the block sizes fluctuate and take complex values. We show (see Section~\ref{sect4}) that this gives a very different form for the distribution $\Pi_\mu(Y)$ (see Figure~\ref{peter-fig4}) from the case of Gaussian disorder (see Figure~\ref{figure1}).

In Section~\ref{sect5} we addressed the case of overlaps between two different temperatures. A recursion relation between negative moments of the partition functions allows us to obtain an expression~\eqref{F17} for the overlaps. In the replica approach this can be interpreted as fluctuating block sizes at the two different temperatures but with a relationship~\eqref{F8} between the block sizes. The fluctuations in the block sizes are continuous and imaginary. They are also not unique, with some freedom to chose the real part. However, we expect the overlap expressions and the physical predictions to be the same for each of the choices.

The REM can be viewed as the paradigmatic model for the one step replica symmetry breaking transition which is known to occur in other spin glass models like the $p$-spin Ising model~\cite{Gross_1984_simplest,Gardner_1985_Spin} or the p-spin spherical spin glass~\cite{Crisanti_1992_spherical}. The obvious advantage of the REM is that there are exact approaches that do not depend on the replica method. In the original random energy model, with Gaussian disorder, both the replica method and the exact methods produce the same expression~\eqref{B10} for the overlap in the thermodynamic limit. In this case the replica method uses Parisi's replica symmetry breaking scheme with a single step of replica symmetry breaking where the $n$ replicas are partitioned into blocks of fixed size $\mu$.

Here we saw that rather simple changes in the energy distribution lead to fluctuations of the block sizes of the overlap matrix. It would be interesting to see whether similar effects can be seen in more complex models. We think that the best candidates are systems with discrete energies for which the ground state can be degenerate. These include the directed polymers problem on a tree~\cite{Derrida_1988_Polymers} when the the energies on each bond take integer values, diluted mean field spin glass models with interactions $J_{ij}= \pm 1$~\cite{Obuchi_2009_Complex}, the binary perceptron~\cite{Gardner_1988_Optimal,Gardner_1989_Three,Krauth_1989_Storage,Huang_2014_Origin,Ding_2019_Capacity} and the K-sat~\cite{Monasson_1996_Entropy,Monasson_1997_Statistical,Mezard_2009_Information}. In all these cases, because the ground state has a non-zero probability of being degenerate, the overlaps $\langle Y_k \rangle \neq 1$ at zero temperature in contrast with what one would expect from~\eqref{A7} in the limit $\mu \to 0$, leading to a $\Pi_\mu(Y)$ consisting of a sum of delta-peaks at values $Y=\frac{1}{p}$ (where $p$ is the degeneracy of the ground state). One can then wonder whether this kind of shape would survive at least at low temperature as it does in our Figure~\ref{peter-fig4} or whether it takes the universal shape of Figure~\ref{figure1} as soon as the temperature is non-zero.

In all the versions of the REM that we considered here, the matrix of overlaps described a one-step RSB. There is no doubt that one can generalize what we have seen here, i.e. overlap matrices with fluctuating block sizes, to systems with several steps of RSB, for example by looking at a generalized random energy model GREM~\cite{Derrida_1985_generalization} with discrete energies. Clearly the statistics predicted by the Ruelle cascade~\cite{Ruelle_1987_mathematical} would then be modified. A challenging question would be to see whether the same phenomenon is present in other models with full RSB.

\subsection*{Gérard Toulouse}


\emph{A la fin des années 70 Gérard Toulouse, en comprenant l'importance du problème des verres de spins, a entrainé dans son sillage plusieurs jeunes théoriciens qui sont devenus des leaders mondialement connus. A nouveau au milieu des années 80, il a été l'un des premiers en France à encourager ses plus jeunes collègues à se lancer dans le domaine des sciences cognitives. Il a ainsi eu une influence majeure sur toute une génération de jeunes physiciens.}

(At the end of the 1970s, Gérard Toulouse, recognised the importance of the spin glass problem and inspired the interest of several young theorists, who have since become world-renowned physicists. Again in the mid-1980s, he was one of the first in France to encourage his younger colleagues to enter the field of cognitive science. He has thus had a major influence on an entire generation of young physicists.)

\appendix
\section{The RSB way of computing overlaps}\label{appA}

This appendix gives the derivation of~\eqref{A7}, \eqref{B10}, \eqref{B11} based on replicas. In the RSB approach~\cite{Mezard_1984_Nature,Mezard_1984_Replica,Mezard_1987_Spin}, one considers $n$ configurations called replicas, and given an overlap $Q$, one assumes that these replicas are organized into $\frac{n}{\mu}$ blocks of $\mu$ replicas. All pairs of replicas have an overlap larger than $Q$ if they belong to the same block and an overlap less than $Q$ if they belong to different blocks. This structure can be represented by the famous $n \times n$ Parisi's matrices shown in Figure~\ref{peter-fig1}.


\begin{figure}[h]
\centerline{\includegraphics[width=6cm]{peter-fig1}}
\caption{Parisi's matrix of overlaps between the $n$ replicas: there are $\frac{n}{\mu} $ blocks of $\mu $ replicas. Pairs of replicas inside a block have an overlap larger than $Q$ and pairs of replicas in different blocks have an overlaps less than $Q$.}\label{peter-fig1}
\end{figure}

According the the RSB scheme~\cite{Mezard_1984_Replica}, the average $\langle Y_k \rangle $ is, in the limit $n \to 0$, the probability that, $k$ different replicas chosen among the $n$ replicas belong all to the same block
\begin{equation}\label{a1}
\left\langle Y_k\right\rangle = \lim_{n\,\to\,0} \left[\frac{n}{\mu} \, \frac{\mu(\mu-1) (\mu-2) \cdots (\mu-k+1)}{n (n-1)(n-2) \cdots (n-k+1)} \right]
\end{equation}
while $\langle Y_{k,k'} \rangle$ is the probability that, choosing $k+k'$ different replicas, all the first $k$ ones belong to one block and the last $k'$ ones belong to a different block
\begin{equation}\label{a2}
\left\langle Y_{k,k'}\right\rangle = \lim_{n\,\to\,0} \left[\frac{n}{\mu} \, \left(\frac{n}{\mu}-1 \right) \, \frac{\Big(\mu(\mu-1) \cdots (\mu-k+1) \Big)
\, \Big(\mu(\mu-1) \cdots (\mu-k'+1) \Big)
}{n (n-1) \cdots (n-k-k'+1)} \right]
\end{equation}
and more generally for $p$ groups
\begin{equation}\label{a3}
\left\langle Y_{k_1, \cdots k_p} \right\rangle = \lim_{n\,\to\,0} \left[(-)^p \frac{\Gamma\left(p-\frac{n}{\mu} \right)}{\Gamma\left(-\frac{n}{\mu}\right)}\, \frac{\Gamma(-n)}{\Gamma\left(k_1 + \cdots k_p-n\right)} \; \prod_{i=1}^p \frac{\Gamma(k_i-\mu)}{\Gamma(- \mu)} \right]
\end{equation}
One can then obtain the expressions~\eqref{A6} using the fact that
\begin{equation}\label{a4}
\left\langle Y^2 \right\rangle= \langle Y_4 \rangle + \langle Y_{2,2} \rangle \quad ; \quad
\left\langle Y^3 \right\rangle= \langle Y_6 \rangle + 3 \langle Y_{4,2} \rangle + \langle Y_{2,2,2} \rangle
\end{equation}
(which can be understood for $\langle Y^2 \rangle$ by noticing that for replicas 1,2 to be in the same group and replicas 3,4 to be also in the same group, either 1,2,3,4 are all in the same group or the pair 1,2 and the pair 3,4 belong to different groups. This reasoning can easily be generalized to establish the expression of $\langle Y^3 \rangle$ in~\eqref{a4}.)


\section{Derivation of various identities of Sections~\ref{sect5} and~\ref{sect6}}\label{appB}

\subsection{ An integer number of replicas at a single temperature: derivation of \eqref{E1},\eqref{E2}}

One way to establish~\eqref{E1} is to expand (as a formal series) in powers of $t$ the following relation (see~\eqref{C1})
\[
\phi'(t) \left\langle e^{-t Z(\beta)} \right\rangle = \frac{d}{dt} \left\langle e^{-t Z(\beta)} \right\rangle
\]
given that (see~\eqref{C2})
\[
\phi'(t) = \sum_{\mu=1}^n \frac{ (-)^\mu \ t^{\mu-1}}{(\mu-1)!} \left\langle Z(\beta \mu) \right\rangle + O\left(t^{n}\right)
\]
Similarly, using the fact that
\[
\left\langle Y_k \, Z(\beta)^k e^{-t Z(\beta)}\right\rangle = \left\langle Z(k \beta) e^{-t Z(\beta)} \right\rangle = \Phi_k(t) \left\langle e^{-t Z(\beta) } \right\rangle
\]
and expanding~\eqref{C6} and $ \langle e^{-t Z(\beta) } \rangle $ in powers of $t$ one gets for integer $n \ge 0$ :
\[
\left\langle Y_k \, Z(\beta)^{k+n} \right\rangle = \sum_{\mu=0}^n \frac{n!}{\mu! \, (n-\mu)!} \left\langle Z\Big((k+\mu)\beta\Big) \right\rangle \, \left\langle Z(\beta)^{n-\mu} \right\rangle
\]
which is the same as~\eqref{E2} if one makes the changes of variables $n\to n-k, \mu \to \mu-k$ and one uses the convention that $(\mu-k)!= \infty$ for $\mu < k$.


\subsection{ A non-integer negative number of replicas at a single temperature: derivation of \eqref{E10}-\eqref{E11}}

One can rewrite~\eqref{C3} as
\[
\left\langle Z(\beta)^n \right\rangle = \frac{-1}{\Gamma(1-n)} \int_0^\infty t^{-n} dt \, \phi'(t) \, \left\langle e^{-t Z(\beta)} \right\rangle \quad \text{for} \; n<0
\]
Using first the identity~\eqref{F4} valid for real $\nu>0$
\[
\left\langle Z(\beta)^n\right\rangle = \frac{-1}{\Gamma(1-n)}\int_{-\infty}^\infty \frac{dy}{2\pi} \int_0^\infty t^{-n} dt \, \phi'(t) \, \Gamma(\nu+iy)t^{-\nu-iy} \left\langle Z(\beta)^{-\nu-iy}\right\rangle
\]
then an expression of $\phi'(t)$ coming from~\eqref{C2} and then
choosing
\[
0< \nu <1-n
\]
to ensure the convergence of the integral over $t$, one gets after integrating over $t$ and then over $E$
\[
\left\langle Z(\beta)^n \right\rangle =
\int_{-\infty}^\infty \frac{dy}{2\pi} \frac{\Gamma(-n-\nu -iy +1)\, \Gamma(\nu+iy)
}{\Gamma(1-n)}
\left\langle Z(n+\nu +iy)\beta)\right\rangle \, \left\langle Z(\beta)^{-\nu-iy}\right\rangle
\]
which reduces to~\eqref{E10} in terms of $\mu=n+\nu$.

Similarly from \eqref{C1},\eqref{C4} one has
\[
\left\langle Y_k Z(\beta)^n\right\rangle = \frac{1}{\Gamma(k-n)} \int_0^\infty t^{k-n-1} dt \, \Phi_k(t)\left\langle e^{-t Z(\beta)} \right\rangle
\]
which becomes using~\eqref{C6} and~\eqref{F4} for $0<\nu <k-n$
\[
\left\langle Y_k Z(\beta)^n\right\rangle =
\int_{-\infty}^\infty \frac{dy}{2\pi} \frac{\Gamma(k-n-\nu -iy +1)\, \Gamma(\nu+iy)
}{\Gamma(k-n)}
\left\langle Z\Big((n+\nu +iy)\beta\Big)
\right\rangle \, \left\langle Z(\beta)^{-\nu-iy}\right\rangle
\]
and this leads to~\eqref{E11} in terms of $\mu=n+\nu$.



\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography
\end{document}
