The random fractional matching problem

Random combinatorial optimization problems (Rcops) are optimization problems in which some parameters of the cost function are random variables. Rcops live therefore at the intersection of mathematical optimization theory and probability theory, and typically the average properties of the solution of the problem are of primary interest. Surprisingly enough, methods developed for the study of disordered systems, such as the replica approach and the cavity method, have proven to be powerful tools for the study of Rcops, they have been able to provide precise analytical predictions [2–6], and they inspired new algorithms for their solution [5–9]. Indeed, the search for an optimal solution in a Rcop can be thought as the search for the ground state of a constrained disordered system, whose energy coincides with the cost function appearing in the original Rcop. Despite the fact that many of the results obtained by this analogy were originally non-rigorously derived, they were later confirmed by rigorous mathematical proofs [1, 10–14].

In this paper we will consider an example of Rcop that belongs to a well studied class of combinatorial optimization problems, namely the class of matching problems. Minimum matching problems are combinatorial optimization problems that can be defined on a graph $\newcommand{\mE}{\mathcal{E}} \newcommand{\mV}{\mathcal{V}} {\mathcal G}=(\mV, \mE)$ having vertex set ${\mathcal V}$ and edge set ${\mathcal E}\subseteq{\mathcal V}\times{\mathcal V}$, both of finite cardinality. The graph is said to be weighted if a weight $\newcommand{\R}{\mathds{R}} w_e\in\R$ is associated to each edge $\newcommand{\mE}{\mathcal{E}} e\in\mE$. A (minimum) matching is therefore a subset $\newcommand{\mE}{\mathcal{E}} {\mathcal M}\subseteq\mE$ of edges that do not share a common vertex, such that the cardinality $|{\mathcal M}|$ is maximal and the cost of the matching, defined as $\sum_{e\in{\mathcal M}}m_ew_e$, is minimum. If ${\mathcal M}$ matches all vertices in $\newcommand{\mV}{\mathcal{V}} \mV$, the matching is said to be perfect [15, 16]. Once the weighted graph is assigned, the solution of a specific instance of the matching problem on a graph can be found quite easily using efficient algorithms available in the literature [5, 17–20].

Random-link matching problems have been the first class of optimization problems to be studied by statistical physics techniques, along with the random-link travelling salesman problem (Rtsp) [3, 21]. The random-link minimum matching problem, simply called random matching problem (Rmp), is defined on the complete graph ${\mathcal K}_{2N}$, i.e. a graph having 2N vertices such that each of them is connected to all the others. In this case, the optimal matching is always perfect. Each edge e of the graph is associated to a random weight w_e, to be extracted, independently from all the others, according to a given probability density $\varrho$, that is assumed to be the same for all edges (see figure 1, left). In the Rmp we search for the set of 'edge occupations' $\boldsymbol{\mathsf m}=\{m_e\}_e$ such that the cost (or energy) function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]:= \sum_{e\in{\mathcal E}}m_e w_e, \nonumber \end{align} \tag{ 1a }$$

is minimized, with the constraints

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_e\in\{0,1\}\quad \forall e\in{\mathcal E},\quad \sum_{e\to v}m_e=1\quad \forall v\in{\mathcal V}. \nonumber \end{align} \tag{ 1b }$$

Zoom In Zoom Out Reset image size

In the previous expression, the sum $\sum_{e\to v}$ runs over all edges having $v$ as an endpoint. Denoting by $\overline{\bullet}$ the average over all instances of the problem, it has been shown that, if $\lim_{w\to 0}\varrho(w)=1$ [3, 14]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{N\to+\infty}\overline{\min_{\boldsymbol{\mathsf m}} E[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]}=\frac{\pi^2}{12}. \label{asintE} \nonumber \end{align} \tag{ 1c }$$

The coefficient of the ${1}/{N}$ correction to this asymptotic result has been derived in [22–24], but its final expression is given in terms of an infinite series, whose terms have to be evaluated numerically, except for the first ones, for which an analytic expression is available.

In the bipartite version of the Rmp, called random assignment problem (Rap), the complete graph ${\mathcal K}_{2N}$ is replaced by a bipartite complete graph, ${\mathcal K}_{N, N}$, in which two sets of vertices, each one of cardinality N, are present, and each vertex of one set is connected to all vertices of the other set only. Both the asymptotic properties [3, 11] and the finite-size corrections to the average optimal cost (aoc) [23, 25] of the Rap have been extensively studied in the literature, and a closed formula is available for the aoc at any value of N in the case of exponentially distributed weights [12, 13, 26], namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\min_{\boldsymbol{\mathsf m}} E[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]}=\sum_{n=1}^{N}\frac{1}{n^2}=\frac{\pi^2}{6}-\frac{1}{N}+o\left(\frac{1}{N}\right).\nonumber \end{align} \tag{ 2 }$$

Both the Rmp and the Rap play the role of mean-field versions of matching problems embedded in the Euclidean space, called random Euclidean matching problems, in which the vertices of the considered graphs correspond to random points in an Euclidean domain and the weight associated to an edge is a function of the Euclidean distance of the endpoints. This fact introduces correlation amongst the edge weights that has been treated as a correction to the purely random case [27, 28], and some analytical results have been obtained for the aoc or its finite-size corrections, using a direct field-theoretical approach [29–31] or, in one dimension, purely probabilistic arguments [32]. Analytical results regarding sparse instances of the matching problem [33–35] have been obtained through the statistical physics cavity method [5].

In the present paper, we will apply the replica formalism to a different type of matching problem, namely the random fractional matching problem (Rfmp), the linear relaxation of the matching problem. We start considering a weighted complete graph ${\mathcal K}_{2N}$. The weights $\boldsymbol{\mathsf w}:= \{w_e\}_e$ associated to the edges $\{e\}_e=:{\mathcal E}$ are non-negative independent random variables, identically distributed according to a probability density $\varrho(w)$, exactly as in the Rmp case. In the Rfmp we search for the set of quantities $\boldsymbol{\mathsf m}:= \{m_e\}_e$ that minimize the cost

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]:= \sum_{e\in{\mathcal E}}m_e w_e, \nonumber \end{align} \tag{ 3a }$$

with the additional constraints

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{costraintsm} m_e\in[0,1]\quad \forall e\in{\mathcal E},\quad \sum_{e\to v}m_e=1\quad \forall v\in{\mathcal V}. \nonumber \end{align} \tag{ 3b }$$

It is easy to show for general graphs that the problem has semi-integer solutions, i.e. optimal configurations ${\mathcal M}$ with $m_e\in\{0, {1}/{2}, 1\}$. ${\mathcal M}$ contains only (odd) cycles and edges that do not share their endpoints (see figure 1, center). It is expected that the aoc of the Rfmp is less than or equal to the aoc of the Rmp obtained with the same weight probability density $\varrho$, due to the fact that all matching configurations feasible for the Rmp, are also feasible for the Rfmp.

If defined on bipartite graphs, the relaxed problem above has integer solutions only, due to the absence of odd cycles, and therefore it is equivalent to the standard Rap. Noticeably, the Belief Propagation algorithm, deeply related to the cavity method, is able to recover the optimal solution of an assignment problem when the solution is unique, and its computational complexity is on par with that one of the best alternative solvers [36]. The statistical physics of linear or convex relaxations of discrete optimization problems has been actively investigated in recent years [37, 38].

In an interesting variation of the Rfmp, that we will call 'loopy Rfmp', an additional non-negative weight $w_{v}$ is associated to each vertex $v\in{\mathcal V}$ of the graph. Each weight $w_v$ is a random variable extracted independently from all other weights in the problem with the same distribution of the edge weights $\varrho(w)$. The loopy Rfmp corresponds therefore to a Rfmp defined on a graph $\hat{{\mathcal K}}_{2N}$ obtained allowing self-loops on ${\mathcal K}_{2N}$ (see figure 1, right). The cost is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]:= \sum_{e\in{\mathcal E}}m_e w_e+2\sum_{v\in{\mathcal V}}m_v w_v, \nonumber \end{align} \tag{ 4a }$$

with the constraints

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{costraintsml} m_e\in[0,1]\quad \forall e\in{\mathcal E},\quad m_v\in[0,1]\quad \forall v\in{\mathcal V},\quad \sum_{e\to v}m_e+2m_v=1\quad \forall v\in{\mathcal V}. \nonumber \end{align} \tag{ 4b }$$

Wästlund proved that, in the loopy Rfmp, in the optimal configuration $m_{e}\in\{0, {1}/{2}, 1\}$ $\forall e$ and $m_v\in\{0, {1}/{2}\}$ $\forall v$ [1]. Remarkably, he also obtained the expression for the aoc for any N assuming $\varrho(w)={\rm e}^{-w}\theta(w)$, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{wastformula} \overline{\min_{\mathsf m}E[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]}=\sum_{n=1}^{2N}\frac{(-1)^{n-1}}{n^2}=\frac{\pi^2}{12}-\frac{1}{8N^2}+o\left(\frac{1}{N^2}\right). \nonumber \end{align} \tag{ 4c }$$

Again, since the loopy Rfmp is a relaxed version of the Rfmp, that in turn is a relaxed version of the Rmp, it follows from equation (4c) that all three Rcops have the same asymptotic aoc, equal to ${\pi^2}/{12}$. The presence of cycles in the Rfmp and cycles and loops in the loopy Rfmp does not affect, therefore, the value of the aoc, but the finite-size corrections only.

In the present paper, we will study, using the replica approach, both the Rfmp and the loopy Rfmp, and their finite-size corrections. We shall consider weight probability densities $\varrho(w)$ with non-negative support and such that $\varrho(w)=1-\mu w+o(w)$ for $w\to 0^+$. In particular, if the weights are exponentially distributed we have $\mu=1$, whereas $\mu=0$ if they are uniformly distributed. We will show, in particular, that the Rmp saddle-point solution naturally appears as asymptotic solution of the problem, alongside with another saddle-point solution of higher cost corresponding to the solution to the Rtsp. Moreover, we will evaluate the finite-size corrections to the aoc on the matching saddle-point, obtaining as the main result of the paper the closed formula

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{costofinaleintro} E_{\varepsilon}(\mu)=\frac{\zeta(2)}{2}+\frac{1}{2N}\left[(\mu-1)\zeta(3)+\frac{1-\varepsilon}{4}\zeta(2)\right]+o\left(\frac{1}{N}\right). \nonumber \end{align} \tag{ 5 }$$

Here $\zeta(z)$ is the Riemann zeta function and, since $\zeta(2)={\pi^2}/{6}$, we obtain the result of equation (1c) at the leading order. The parameter ε takes value  +1 if self-loops are allowed in the model, −1 otherwise, i.e. when m_ii is not present. The analytic predictions from the replica calculation are supported by the numerical simulations in the last section of the paper.

In the spirit of the seminal works of Orland [21] and Mézard and Parisi [3], let us first write down the partition function for the Rfmp, both in its usual and in its loopy version. As anticipated, we will consider a random-link formulation of the problem on the complete graph with 2N vertices $\newcommand{\mV}{\mathcal{V}} \mV=\{i\}_{i=1, \dots, 2N}$, where the weights $\{w_{ij}\}_{ij}$ are non-negative independent random variables identically distributed with distribution $\varrho(w)$. As anticipated, in the following, we will consider a particular class of probability densities, i.e. we will assume that $\varrho(w)=1-\mu w+o(w)$ for $w\to 0^+$. The uniform distribution on the unit interval, $\varrho(w)=\theta(w)\theta(1-w)$, and the exponential distribution on the positive real axis, $\varrho(w)={\rm e}^{-w}\theta(w)$, belong to this class, with $\mu=0$ and $\mu=1$ respectively. Due to the fact that all distributions in this class have the same limit for $w\to 0^+$, i.e. $\lim_{w\to 0^+}\varrho(w)=1$, we expect that the asymptotic aoc will be the same for all of them, as it happens in the Rmp. Indeed, in the general framework of the analysis of the Rmp performed in [23, 25, 39] it is easily seen that, whereas the asymptotic cost only depends on the behavior of the first term in the Maclaurin expansion of $\varrho(w)\simeq\varrho_0 w^r$, the $O({1}/{N})$ finite-size corrections depend on the expansion up to, at least, the second term [23], and the power r also affects the scaling of the corrections themselves [25].

Let us start observing that, in the Rfmp, the occupation number $m_{ij}=m_{ji}$ of the edge $(i, j)$ between the node i and the node j can assume the values

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_{ij}\in\{0,1,2\},\quad {\rm{with~the~constraint~}} \sum_{j=1}^{2N}m_{ij}+\varepsilon m_{ii}=2,\quad 1\leqslant i\leqslant 2N, \nonumber \end{align} \tag{ 6 }$$

where $\varepsilon=+1$ if loops are allowed, $\varepsilon=-1$ otherwise. The parameter ε, therefore, allows us to switch between the two variations of the model described in the Introduction. The cost of a given matching configuration is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{costoE} E_{\varepsilon}[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]=\frac{1}{2}\sum_{i\leqslant j}\left(1+\varepsilon\delta_{ij}\right)m_{ij}w_{ij}. \nonumber \end{align} \tag{ 7 }$$

For calculation convenience we consider, for each edge e, $m_{e}\in\{0, 1, 2\}$ and not $m_{e}\in\{0, {1}/{2}, 1\}$ as in the Introduction. This fact does not affect the results apart from the necessary rescaling of the cost that indeed we have introduced in equation (7). The partition function for a given instance of the problem can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_\varepsilon(\beta):= \sum_{m_{ij}\in\{0,1,2\}}\left[\prod_{i=1}^{2N}{\mathbb I}\left(\sum_{j=1}^{2N}m_{ij}+\varepsilon m_{ii}=2\right)\right]{\rm e}^{-2\beta N E_{\varepsilon}[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]}, \nonumber \end{align} \tag{ 8 }$$

where ${\mathbb I}(\bullet)$ is an indicator function that is equal to one when the condition in the brackets is satisfied, and zero otherwise. The aoc is recovered as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\varepsilon}(\mu):= \overline{\min_{\boldsymbol{\mathsf m}}E_{\varepsilon}[\boldsymbol{\mathsf m},\boldsymbol{\mathsf w}]}=-\lim_{\beta\to+\infty}\frac{1}{2N}\frac{\partial\overline{\ln Z(\beta)}}{\partial\beta}. \nonumber \end{align} \tag{ 9 }$$

Note that we have made explicit the dependence of the aoc on the value of $\mu=-\left.\partial_w\varrho\right|_{w=0}$. From the results in [1], we know that $E_{\varepsilon}(\mu)=O(1)$, i.e. the aoc is not extensive, and the cost density scales as $O\left({1}/{N}\right)$: this is indeed expected, due to the fact that the shortest link amongst N scales as ${1}/{N}$. We have therefore rescaled β in the exponent of the partition function accordingly, in such a way that a finite thermodynamical limit at fixed β can be obtained, and an extensive functional $2NE_{\varepsilon}(\mu)=O(N)$ appears in the exponent in the low-temperature regime.

To average over the disorder, we use the following integral representation of the Kronecker delta,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \delta_{a,0}=\int_{0}^{2\pi}{\rm e}^{{\rm i}\lambda a}\frac{\dd\lambda}{2\pi}, \nonumber \end{align} \tag{ 10 }$$

and we apply, as usual, the replica trick. The average replicated partition function for the fractional matching problem can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \label{Zpm} \overline{Z^n_\varepsilon} =\left[\prod_{a=1}^n\prod_{i=1}^{2N}\int_0^{2\pi}\frac{{\rm e}^{-2{\rm i}\lambda_i^a}\dd\lambda_i^a}{2\pi}\right]\prod_{i<j}\left(1+\frac{T_{ij}}{N}\right)\prod_{i=1}^{2N}\left(1+\frac{\varepsilon+1}{2}\frac{R_{i}}{N}\right). \nonumber \end{align} \tag{ 11 }$$

In equation (11) we have introduced the quantities

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1+\frac{T_{ij}}{N}&:= \overline{\prod_{a=1}^n\left[1+\exp\left({\rm i}\lambda_i^a+{\rm i}\lambda_j^a-\beta N w_{ij}\right)+\exp\left(2{\rm i}\lambda_i^a+2{\rm i}\lambda_j^a-2\beta N w_{ij}\right)\right]}\nonumber \\ &=1+\frac{1}{N}\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}\hat g_{|\alpha|+2|\beta|}\exp\left[{\rm i}\sum_{a\in\alpha}\left(\lambda_i^a + \lambda_j^a\right) + 2{\rm i}\sum_{b\in\beta}\left(\lambda_i^b + \lambda_j^b\right)\right],\nonumber \end{align} \tag{ 12a }$$

and, in the presence of loops, the on-site contribution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1+\frac{R_{i}}{N}&:= \overline{\prod_{a=1}^n\left[1+\exp\left(2{\rm i}\lambda_i^a-2\beta N w_{ii}\right)\right]}\nonumber \\ &=1+\frac{1}{N}\sum_{{\alpha\neq\emptyset}}\hat g_{2|\alpha|}\exp\left[2{\rm i}\sum_{a\in\alpha}\lambda_i^a\!\right].\nonumber \end{align} \tag{ 12b }$$

In equation (12) the sums run over the elements of ${\mathcal P}([n])$, set of subsets of $[n]:= \{1, \dots, n\}$, and we have denoted the cardinality of $\alpha\in{\mathcal P}([n])$ by $|\alpha|$. We have also introduced the quantity

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \label{gp} \hat g_p&:= N\int_0^{+\infty}{\rm e}^{-\beta p N w}\varrho(w)\dd w\nonumber \\ &=g_p\left[1-\mu\frac{g_p}{N}+o\left(\frac{1}{N}\right)\right],\quad {\rm{where}}\ g_p:= \frac{1}{\beta p}. \nonumber \end{align} \tag{ 13 }$$

Using the previous equation, and in order to evaluate the aoc and its first finite-size correction, we use the fact that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 2R_i=&2\sum_{{\alpha\neq\emptyset}}\hat g_{2|\alpha|}\exp\left[2{\rm i}\sum_{a\in\alpha}\lambda_i^a\!\right]\nonumber \\ & =\sum_{{\alpha\neq\emptyset}}\left[g_{|\alpha|}-\frac{2g_{2|\alpha|}^2}{N}+o\left(\frac{1}{N}\right)\right]\exp\left[2{\rm i}\sum_{a\in\alpha}\lambda_i^a\!\right]\sim T_{ii}+O\left(\frac{1}{N}\right),\nonumber \end{align} \tag{ 14 }$$

since $2g_{2p}=g_{p}$. In the last step we have used the fact that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{ii}=\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}\hat g_{|\alpha|+2|\beta|}\exp\left[2{\rm i}\sum_{a\in\alpha}\lambda_i^a + 4{\rm i}\sum_{b\in\beta}\lambda_i^b\right] \nonumber \end{align} \tag{ 15 }$$

gives zero contribution, unless $\newcommand{\e}{{\rm e}} \beta=\emptyset$, due to the overall constraint imposed by the integration on $\{\lambda_i^a\}$. Neglecting $O({1}/{N})$ terms in the exponent (i.e. $O({1}/{N^2})$ to the cost), we can write the partition function as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \overline{Z^n_\varepsilon} =\left[\prod_{a=1}^n\prod_{i=1}^{2N}\int_0^{2\pi}\frac{{\rm e}^{-2{\rm i}\lambda_i^a}\dd\lambda_i^a}{2\pi}\right] \exp\left[\frac{1}{2N}\sum_{i,j}\left(T_{ij}-\frac{T_{ij}^2}{2N}\right)+\frac{\varepsilon-1}{4N}\sum_{i=1}^{2N}T_{ii}+O\left(\frac{1}{N}\right)\right]. \nonumber \end{align} \tag{ 16 }$$

We introduce now the placeholders

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_{\alpha,\beta}:= {\mathbb I}(\alpha\cap\beta=\emptyset) \sum_{i=1}^{2N}\exp\left({\rm i}\sum_{a\in\alpha}\lambda_i^a+2{\rm i}\sum_{b\in\beta}\lambda_i^b\right), \nonumber \end{align} \tag{ 17 }$$

that allow us to write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{i,j}T_{ij}+\frac{\varepsilon-1}{2}\sum_i T_{ii} =\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}\hat g_{|\alpha|+2|\beta|} q_{\alpha,\beta}^2+\frac{\varepsilon-1}{2}\sum_{\alpha\neq\emptyset}\hat g_{|\alpha|}q_{\emptyset,\alpha}, \nonumber \end{align} \tag{ 18a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \sum_{i,j}T_{ij}^2 =\sideset{}{'}\sum_{{\alpha,\beta|\hat\alpha,\hat\beta}}\hat g_{|\alpha|+2|\beta|}\hat g_{|\hat \alpha|+2|\hat \beta|} q_{\alpha\triangle \hat \alpha,\beta\cup\hat\beta\cup(\alpha\cap\hat\alpha)}^2.\label{dst} \nonumber \end{align} \tag{ 18b }$$

In the previous equations $\newcommand{\tr}{{\rm tr}} \alpha\triangle\beta:= (\alpha\setminus\beta)\cup(\beta\setminus\alpha)$ is the symmetric difference of the sets α and β, we have denoted by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sideset{}{'}\sum_{\alpha,\beta|\hat\alpha,\hat\beta}=\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}\quad \sum_{{{\hat\alpha\cap\hat \beta=\emptyset}\atop{\hat\alpha\cup\hat \beta\neq\emptyset}}}{\mathbb I}(\beta\cap\hat\beta=\emptyset){\mathbb I}((\alpha\cup\hat\alpha)\cap(\beta\cup\hat\beta)=\emptyset). \nonumber \end{align} \tag{ 19 }$$

Denoting by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \varphi_{\alpha,\beta}:= \sideset{}{'}\sum_{\sigma,\rho|\hat\sigma,\hat\rho}\hat g_{|\rho|+2|\sigma|}\hat g_{|\hat \rho|+2|\hat \sigma|}{\mathbb I}\left(\rho\triangle\hat\rho=\alpha\right){\mathbb I}\left(\sigma\cup\hat\sigma\cup(\rho\cap\hat\rho)=\beta\right), \nonumber \end{align} \tag{ 20 }$$

and using a Hubbard–Stratonovich transformation in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle &\exp\left(\frac{\hat g_{|\alpha|+2|\beta|} -{1}/{2N}\,\varphi_{\alpha,\beta}}{2N}q_{\alpha,\beta}^2\right)\nonumber \\ &=\sqrt{\frac{N}{2\pi \left(\hat g_{|\alpha|+2|\beta|} - {\varphi_{\alpha,\beta}}/{2N}\right)}}\int_{-\infty}^{+\infty}\exp\left(-\frac{NQ_{\alpha,\beta}^2}{2\hat g_{|\alpha|+2|\beta|} - {\varphi_{\alpha,\beta}}/{N}} +Q_{\alpha,\beta}q_{\alpha,\beta}\right)\!\dd Q_{\alpha,\beta}\nonumber \\ &\simeq \sqrt{\frac{N}{2\pi \hat g_{|\alpha|+2|\beta|}}}\int_{-\infty}^{+\infty}\exp\left(-\frac{NQ_{\alpha,\beta}^2}{2\hat g_{|\alpha|+2|\beta|}} -\frac{\varphi_{\alpha,\beta}Q^2_{\alpha,\beta}}{4\hat g^2_{|\alpha|+2|\beta|}}+Q_{\alpha,\beta}q_{\alpha,\beta}+O\left(\frac{1}{N}\right)\right)\dd Q_{\alpha,\beta} \nonumber \end{align} \tag{ 21 }$$

we can finally introduce the order parameters $Q_{\alpha, \beta}$ as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{quisopra} &\left[\prod_{a=1}^n\prod_{i=1}^{2N}\int_0^{2\pi}\frac{{\rm e}^{-2{\rm i}\lambda_i^a}\dd\lambda_i^a}{2\pi}\right]\exp\left[\frac{1}{2N}\sum_{i,j}\left(T_{ij}-\frac{T_{ij}^2}{2N}\right)+\frac{\varepsilon-1}{4N}\sum_{i=1}^{2N}T_{ii}\right]\nonumber \\ \simeq&\left[\prod_{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}\int\dd Q_{\alpha,\beta}\sqrt{\frac{N}{2\pi \hat g_{|\alpha|+2|\beta|}}}\right]\nonumber \\ &\times\exp\left[-N\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}\ \frac{Q_{\alpha,\beta}^2}{2\hat g_{|\alpha| + 2|\beta|}} + 2N\ln z[\bQ]\! -\sideset{}{'}\sum_{{\alpha,\beta|\hat\alpha,\hat\beta}}\frac{g_{|\alpha| + 2|\beta|}g_{|\hat\alpha| + 2|\hat\beta|}}{4 g^2_{|\alpha| + |\hat\alpha| + 2|\beta| + 2|\hat\beta|}}Q_{\alpha\triangle\hat\alpha,\beta\cup\hat\beta\cup(\alpha\cap\hat\alpha)}^2 \right]. \nonumber \end{align} \tag{ 22 }$$

The expression of $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} \ln z[\bQ]$, in which we have exponentiated the integration on $\{\lambda_i^a\}$, is given in equation (A.1) in appendix A. Equation (22) generalizes the equivalent expression for the partition function obtained for the Rmp in [3, 22, 23], that is recovered imposing that $\newcommand{\e}{{\rm e}} Q_{\alpha, \beta}\equiv Q_{\beta}\delta_{|\alpha|, 0}$.

Observe that in equation (22) $\hat g_p$ appears, a quantity defined in equation (13). We have that $\lim_{N\to +\infty}\hat g_p=g_p$. However, expanding for large N the quantities $\hat g_p$ in equation (22), new ${1}/{N}$ finite-size corrections to the cost will appear⁶. In appendix A we show that the replicated action can be finally written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{ints} \overline{Z^n} \simeq\left[\prod_{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}\int\dd Q_{\alpha,\beta}\sqrt{\frac{N}{2\pi \hat g_{|\alpha|+2|\beta|}}}\right]{\rm e}^{-N{\mathcal S}[\bQ]}\nonumber \\ {\mathcal S}[\bQ]:= S[\bQ]+\Delta S^{\rm T}[\bQ]+\Delta S^{\varrho}[\bQ]+o\left({1}/{N}\right).\nonumber \end{align} \tag{ 23a }$$

The three contributions appearing in the previous expression are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle S[\bQ] := \sum_{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}\frac{Q_{\alpha,\beta}^2}{2g_{|\alpha|+2|\beta|}}-2\ln z_0[\bQ],\label{s0} \nonumber \end{align} \tag{ 23b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \Delta S^{\rm T}[\bQ] =\frac{1}{N}\sideset{}{'}\sum_{{\alpha,\beta|\hat\alpha,\hat\beta}}\frac{g_{|\alpha|+2|\beta|}g_{|\hat\alpha|+2|\hat \beta|}}{4 g^2_{|\alpha\triangle \hat\alpha|+2|\beta|+2|\hat\beta|+2|\alpha\cap\hat\alpha|}}Q_{\alpha\triangle\hat \alpha,\beta\cup\hat\beta\cup(\alpha\cap\hat\alpha)}^2-\frac{\varepsilon-1}{2N}\sum_{\alpha\neq\emptyset}g_{|\alpha|}\frac{\partial \ln z_0[\bQ]}{\partial Q_{0,\alpha}}, \nonumber \end{align} \tag{ 23c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \Delta S^{\varrho}[\bQ] := \frac{\mu}{2\beta N}\sum_{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}Q_{\alpha,\beta}^2, \nonumber \end{align} \tag{ 23d }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{z0} z_0[\bQ]:= \lim_{{N\to\infty}} z[\bQ]= \left[\prod_{a=1}^n\int_{0}^{2\pi}\frac{{\rm e}^{-2{\rm i}\lambda^a}\dd\lambda^a}{2\pi}\right] \exp\left[\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}Q_{\alpha,\beta}\exp\left({\rm i}\sum_{a\in\alpha}\lambda^a+2{\rm i}\sum_{b\in\beta}\lambda^b\right)\right] \nonumber \end{align} \tag{ 23e }$$

is the (leading) one-site partition function. It follows that S is the leading term in ${\mathcal S}$. The $\Delta S^{\rm T}$ term contains the finite-size correction due to the re-exponentiation and to the one-site partition function z. Finally, the $\Delta S^{\varrho}$ term contains an additional contribution due to the finite-size corrections to g_p appearing in equation (13). Observe, once again, that this contribution is absent in the case of flat distribution.

The integral in equation (23a) can now be evaluated using the saddle-point method. The saddle-point equation for $Q_{\alpha, \beta}$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{saddlepointQ} \frac{Q_{\alpha,\beta}}{g_{|\alpha|+2|\beta|}}=2\frac{\partial \ln z_0[\bQ]}{\partial Q_{\alpha,\beta}}=2\left\langle\exp\left({\rm i}\sum_{a\in\alpha}\lambda^a+2{\rm i}\sum_{b\in\beta}\lambda^b\right)\right\rangle_{z_0}, \nonumber \end{align} \tag{ 24 }$$

where $\langle\bullet\rangle_{z_0}$ is the average performed respect to the one-site partition function z₀. Denoting by $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} \bQ^{\rm sp}$ the solution of equation (24), the replicated partition function becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{sfinale} \overline{Z^n}\simeq \exp\left(-N{\mathcal S}[\bQ^{\rm sp}]-\frac{\ln\det\boldsymbol{\mathsf \Omega}[\bQ^{\rm sp}]}{2}\right), \nonumber \end{align} \tag{ 25 }$$

where $\boldsymbol{\mathsf \Omega}$ is the Hessian matrix of $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} S[\bQ]$ evaluated on the saddle-point $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} \bQ^{\rm sp}$, solution of equation (24), i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{Omega} \Omega_{\alpha\beta,\hat\alpha\hat\beta}[\bQ^{\rm sp}]:= \sqrt{g_{|\alpha|+2|\beta|}g_{|\hat\alpha|+2|\hat\beta|}}\left.\frac{\partial^2 S[\bQ]}{\partial Q_{\alpha,\beta}\partial Q_{\hat\alpha,\hat \beta}}\right|_{\bQ=\bQ^{\rm sp}}. \nonumber \end{align} \tag{ 26 }$$

The additional term $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} -{1}/{2}\ln\det\boldsymbol{\mathsf \Omega}[\bQ^{\rm sp}]$ provides, in general, a nontrivial finite-size correction to the leading term [22, 23].

2.1. Replica symmetric ansatz and matching saddle-point

To proceed further, let us assume that a replica-symmetric ansatz holds, i.e. we search for a solution of the saddle-point equation in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q_{\alpha,\beta}\equiv Q_{|\alpha|,|\beta|}. \nonumber \end{align} \tag{ 27 }$$

This is a common and successful assumption in the study of random combinatorial optimization problems [3, 27, 39] that greatly simplifies the calculation. The stability of a replica symmetric solution is however not obvious and, in general, a replica symmetry breaking might occur [5]. In the present paper, this assumption will be justified a posteriori, on the basis of the agreement between the analytical predictions and the numerical computation. Using the replica symmetric hypothesis, the leading term in equation (23b) becomes

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{Snsmall} &S[\bQ]=\sum_{p+q\geqslant 1}\binom{n}{p\ q}\frac{Q_{p,q}^2}{2g_{p+2q}}-2\ln z_0[\bQ]\nonumber \\ &\qquad \qquad \qquad \qquad \xrightarrow{n\to 0} n\sum_{{p+q\geqslant 1}}\frac{(-1)^{\,p+q-1}}{p+q}\binom{p+q}{q}\frac{Q_{p,q}^2}{2g_{p+2q}}-2n\lim_{n\to 0}\frac{\ln z_0[\bQ]}{n}.\nonumber \end{align} \tag{ 28 }$$

In the previous expression, and in the following, we will adopt the notation

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \binom{a}{b_1\cdots b_s}:= \frac{\Gamma(a+1)}{\Gamma\left(a+1-\sum_{i=1}^sb_i\right)\prod_{i=1}^s\Gamma(b_i+1)}. \nonumber \end{align} \tag{ 29 }$$

Even under the replica symmetric hypothesis, the evaluation of $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} \lim_{n\to 0}\frac{1}{n}\ln z_0[\bQ]$ remains nontrivial. However, in appendix A we show that a special replica symmetric saddle-point solution exists, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{matchingsp} Q_{p,q}^{\rm sp}=\delta_{p,0}Q_{0,q}\equiv \delta_{p,0}Q_q, \nonumber \end{align} \tag{ 30 }$$

corresponding to the replica-symmetric saddle-point solution of the Rmp. This fact is not surprising: as anticipated in the Introduction, the aoc of the Rfmp coincides with the aoc of the Rmp in the $N\to+\infty$ limit, and, indeed, the evaluation of $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} z_0[\bQ]$ on the matching saddle-point can be performed exactly, and coincides with the one of the Rmp [3, 22, 23]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{zmatchingsp1}\ln z_0[\bQ^{\rm sp}] =n\int_{-\infty}^{+\infty}\left({\rm e}^{-{\rm e}^x}-{\rm e}^{-G(x)}\right)\dd x, \nonumber \end{align} \tag{ 31a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{zmatchingsp2}G(x) := \sum_{k=1}^\infty \frac{(-1)^{k-1}}{k!}Q_k{\rm e}^{xk}. \nonumber \end{align} \tag{ 31b }$$

The saddle-point equations become

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle Q_{p,q}=\delta_{p,0}Q_q=\frac{\delta_{p,0}}{\beta}\int_{-\infty}^{+\infty}\frac{{\rm e}^{qy-G(y)}}{q!}\dd y, \nonumber \end{align} \tag{ 32 }$$

implying the self-consistent equation for G given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle G(x)= \frac{1}{\beta}\int_{-\infty}^{+\infty}B(x+y){\rm e}^{-G(y)}\dd y, \nonumber \end{align} \tag{ 33 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle B(x):= \sum_{k=1}^\infty (-1)^{k-1} \frac{{\rm e}^{kx}}{\Gamma^2(k+1)}. \nonumber \end{align} \tag{ 34 }$$

As expected, the leading contribution corresponds therefore to the Rmp free-energy,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \frac{1}{n}S[\bQ^{\rm sp}]&=\beta\sum_{k=1}^\infty (-1)^{k-1} Q_{k}^2 -2 \int_{-\infty}^{+\infty}\left({\rm e}^{-y}-{\rm e}^{-G(y)}\right)\dd y \nonumber \\ & =\int_{-\infty}^{+\infty}G(y){\rm e}^{-G(y)}\dd y-2\int_{-\infty}^{+\infty}\left({\rm e}^{-y}-{\rm e}^{-G(y)}\right)\dd y. \nonumber \end{align} \tag{ 35 }$$

In the $\beta\to+\infty$ limit we can introduce

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat G(x):= G\left(\beta x\right) \nonumber \end{align} \tag{ 36 }$$

and we can use the fact that [3]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{\beta\to+\infty} B(\beta x)=\theta(x). \nonumber \end{align} \tag{ 37 }$$

We can solve for $\hat G$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \label{saddleG} \hat G(x)=\int_{-x}^{+\infty}{\rm e}^{-\hat G(y)}\dd y\Longrightarrow \hat G(x)=\ln(1+{\rm e}^x). \nonumber \end{align} \tag{ 38 }$$

We finally obtain that, in both the considered formulations, the aoc of the Rfmp is equal to the aoc of the Rmp, as anticipated in the Introduction,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \lim_{N\to+\infty}E_\varepsilon(\mu)=\frac{1}{2}\int_{-\infty}^{+\infty}\hat G(l){\rm e}^{-\hat G(l)}\dd l=\frac{\pi^2}{12}. \nonumber \end{align} \tag{ 39 }$$

2.2. Finite-size corrections

In equation (25) three contributions to the finite-size corrections appear, namely $\Delta S^\varrho$, $\Delta S^{\rm T}$, and ${1}/{2}\ln\det\boldsymbol{\mathsf\Omega}$. The first contribution can be evaluated straightforwardly as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \Delta S^{\varrho}[\bQ^{\rm sp}]&=\frac{\mu}{2N}\sum_{\beta\neq \emptyset} Q_{|\beta|}^2 =\frac{\mu}{2N}\sum_{q=1}^\infty\binom{n}{q} Q_{q}^2=\frac{n\mu}{2N}\sum_{q=1}^\infty\frac{(-1)^{q-1}}{q} Q_{q}^2+o(n) \nonumber \\ &=\frac{n\mu}{2N\beta}\int_{-\infty}^{+\infty}\dd y\int_0^{+\infty}\dd x\,G(y-x){\rm e}^{-G(y)}+o(n). \nonumber \end{align} \tag{ 40 }$$

In the $\beta\to+\infty$ limit we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \lim_{\beta\to+\infty}\lim_{n\to 0}\frac{\Delta S^{\varrho}[\bQ^{\rm sp}]}{n\beta}=\frac{\mu}{2N}\int_{-\infty}^{+\infty}\dd y\int_0^{+\infty}\dd x\,\hat G(y-x){\rm e}^{-\hat G(y)}=\frac{\mu}{N}\zeta(3). \nonumber \end{align} \tag{ 41 }$$

The last integration was performed in [22]. The $\Delta S^{\rm T}$ contribution depends on ε, i.e. on the presence or not of loops. In appendix B we show that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \lim_{\beta\to+\infty}\lim_{n\to 0}\frac{\Delta S^{\rm T}[\bQ^{\rm sp}]}{n\beta}=-\frac{\varepsilon-1}{4N}\zeta(2)-\frac{1}{N}\zeta(3). \nonumber \end{align} \tag{ 42 }$$

Finally, the fluctuation contribution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\Delta E^{\boldsymbol{\mathsf \Omega}}}{N}:= \lim_{\beta\to+\infty}\lim_{n\to 0}\frac{\ln\det\boldsymbol{\mathsf \Omega}[\bf Q^{\rm sp}]}{4n\beta}. \nonumber \end{align} \tag{ 43 }$$

The expression of this fluctuation term is more involved than the corresponding one for the Rmp, that has been studied in [22, 23]. As it happens in that case, an exact evaluation through the replica formalism is quite complicated (we give some details in appendix C). However, observe that $\Delta E^{\boldsymbol{\mathsf \Omega}}$ depends on neither ε nor μ. We can therefore avoid the complex, direct evaluation and extract its value from Wästlund's formula in equation (4c) for the aoc in the loopy Rfmp with exponentially distributed weights, comparing our result in this specific case with Wästlund's one. In particular, equation (4c) predicts no ${1}/{N}$ corrections in the loopy Rfmp, implying the simple result

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta E^{\boldsymbol{\mathsf \Omega}}=0. \nonumber \end{align} \tag{ 44 }$$

Moreover, in the spirit of the analysis in [40], the fact that $\Delta E^{\boldsymbol{\mathsf \Omega}}$ is a finite and well-defined quantity also suggests that the Hessian $\boldsymbol{\mathsf \Omega}$ remains positive definite within the replica symmetric ansatz for $\beta\to+\infty$, and, therefore, that the replica symmetric solution remains stable.

Collecting all contributions, we can finally write down a general expression for the aoc for the Rfmp, and its finite-size corrections, as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{asintoticafinale} E_{\varepsilon}(\mu)=\frac{\zeta(2)}{2}+\frac{1}{2N}\left[(\mu-1)\zeta(3)+\frac{1-\varepsilon}{4}\zeta(2)\right]+o\left(\frac{1}{N}\right). \nonumber \end{align} \tag{ 45 }$$

The analytic results in equation (45) have been verified numerically using the LEMON graph library [41]. For each one of the considered models, the results have been obtained averaging over at least $3\cdot 10^6$ instances for each value of N. In figure 2 we show our results, both for the case of uniform weight distribution ($\mu=0$) and the case of exponential weight distribution ($\mu=1$). In both cases it is evident that the aoc of the Rmp is greater than the corresponding aoc of the Rfmp, and similarly $E_{-1}(\mu)\geqslant E_{+1}(\mu)$, as expected. The asymptotic formula, equation (45), has been verified for all cases. We performed moreover a parametric fit of our data using the fitting function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{fitrfmp} E=\frac{\pi^2}{12}+\frac{a}{2N}+\frac{b}{4N^2}+\frac{c}{8N^3} \nonumber \end{align} \tag{ 46a }$$

to verify our predictions for the Rfmp in two ways, i.e. either assumiming all parameters free, or fixing a equal to our analytical prediction (when available) to improve our evaluation of b and c. We also performed, for comparison, a similar fit for the Rmp, using the fitting function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{fitrmp} E=\frac{\pi^2}{12}+\frac{a}{2N}+\frac{\hat c}{(2N)^{{3}/{2}}}+\frac{b}{4N^2}, \nonumber \end{align} \tag{ 46b }$$

to take into account the anomalous scaling of the corrections in this case [28]. The results are summarized in table 1. It is remarkable that the $N^{-{3}/{2}}$ correction introduced in equation (46a) for the Rmp cost, that is numerically present in agreement with the results in [28], is absent in all considered variants of the Rfmp, as we numerically verified and as analytically predicted in equation (4c) for the loopy Rfmp with exponential weight distribution.

Zoom In Zoom Out Reset image size

Table 1. Results of a fitting procedure of the aoc obtained numerically compared with the theoretical predictions. Wästlund's formula predicts a  =  0 and $b=-c=-{1}/{2}$ for the aoc of the Rfmp with $\mu=1$ and $\newcommand{\e}{{\rm e}} \epsilon=+1$.

	Problem	Theoretical a	a, b, c/$\hat c$ free			b, c free
			a	b	c or $\hat c$	b	c
$\mu=0$	Rmp	—	$-0.049(4)$	$0.03(9)$	$-1.54(4)$	—	—
	Rfmp $\varepsilon=+1$	$-1.202\, 0569\dots$	$-1.204(1)$	$1.26(5)$	$-1.0(4)$	$1.19(2)$	$-0.5(2)$
	Rfmp $\varepsilon=-1$	$-0.379\, 5898\dots$	$-0.381(1)$	$-2.43(5)$	$4.7(4)$	$-2.50(2)$	$5.2(2)$
$\mu=1$	Rmp	–	$1.131(3)$	$-0.03(8)$	$-1.13(3)$	—	—
	Rfmp $\varepsilon=+1$	0	$-0.001(1)$	$-0.48(5)$	$0.5(5)$	$-0.53(2)$	$0.9(2)$
	Rfmp $\varepsilon=-1$	$0.822\, 4670\dots$	$0.821(1)$	$-1.15(6)$	$0.5(5)$	$-1.19(3)$	$0.9(3)$

Equation (45) allows us to make predictions about differences of aoc for different types of models, due to the fact that we have isolated the different contributions depending on the presence of loops, or on the chosen distribution $\varrho$. For example, we expect that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{dif\!floops} \delta E_\ell:= E_{-1}(\mu)-E_{+1}(\mu)=\frac{\zeta(2)}{4N}+o\left(\frac{1}{N}\right). \nonumber \end{align} \tag{ 47a }$$

Similarly, we have that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{dif\!frho} \delta E_\varrho:= E_{\varepsilon}(1)-E_{\varepsilon}(0)=\frac{\zeta(3)}{2N}+o\left(\frac{1}{N}\right). \nonumber \end{align} \tag{ 47b }$$

Both the relations above have been verified numerically. Our results are shown in figure 3 and they are in agreeement with equation (47).

Zoom In Zoom Out Reset image size

The Rfmp on the complete graph ${\mathcal K}_{2N}$ generalizes the more famous Rmp, allowing cycles in the optimal solution. Here we have studied, using the replica formalism, its aoc and we have compared it with the aoc of the Rmp. We have considered the case of random weights on the graph edges distributed according to $\varrho(w)=1-\mu w+o(w)$, and we have evaluated the asymptotic expression of the aoc on the matching saddle-point, obtaining for the Rfmp the same asymptotic aoc of the Rmp. Remarkably enough, another saddle-point solution naturally appears in the calculation, corresponding to the Rtsp solution, whose average optimal cost is however higher.

We have also explicitly obtained the finite-size corrections to the asymptotic aoc in two variations of the model, the standard Rfmp and the loopy Rfmp, in which each node can be matched to itself (loop). For the latter model, in particular, an explicit formula for the aoc had been found by Wästlund for any value of N in the case of exponentially distributed weights [1]. The two models have different, non-trivial finite-size corrections.

By means of our approach, we have been able to separate the different contributions in the finite size corrections coming from different details of the problem, namely a first one due to the possibility of having loops, a second one due to the chosen weight distribution and a third one that is independent from the aforementioned characteristics. The first and the second contributions have been explicitly derived, whereas the third one has been obtained comparing Wästlund's result with ours in the corresponding specific case. We have been able to write down an explicit expression for the aoc for all the considered cases up to $o({1}/{N})$ corrections, given in equation (5). We finally verified our results in all the considered cases comparing them with the values for the aoc obtained numerically. Going beyond the ${1}/{N}$ corrections, numerical results also suggest that the anomalous scaling ${1}/{N^{{3}/{2}}}$ in higher order corrections in the Rmp, discussed in [28] and verified in the present work, does not appear in the Rfmp.

Despite the fact that an explicit expression for the aoc has been found for the Rfmp in all the considered formulations of it, a number of open problems still remains. For example, a complete replica calculation of the contribution $\Delta E^{\boldsymbol{\mathsf \Omega}}$ deriving from $\ln\det\boldsymbol{\mathsf \Omega}$ is still missing. We have been able to estimate this contribution relying on Wästlund's results. However, an explicit replica calculation is still interesting for a series of reasons. One of them is that, as shown in appendix C, $\Delta E^{\boldsymbol{\mathsf \Omega}}$ appears to be the sum of two terms, the first one identical to a corresponding one appearing in the Rmp that is known to be nonzero, the second one that we expect, from Wästlund's formula, to be exactly opposite to the former one, in such a way that $\Delta E^{\boldsymbol{\mathsf \Omega}}=0$. The presence of this cancellation suggests that the contribution of cycles plays an important role in the anomalous scaling of higher corrections in the Rmp, and in its absence in the Rfmp. Remembering also that an explicit formula for the ${1}/{N}$ correction in the Rmp is still missing, further investigations in this direction are in order, to shed some light on the problem of finite-size corrections both in the Rfmp and in the Rmp even beyond the $O({1}/{N})$ terms.

The authors are grateful to Johan Wästlund for useful discussions and correspondence. GP and GS acknowledge the financial support of the Simons foundation (Grant No. 454949, Giorgio Parisi). The work presented in this article was supported by the project 'Meccanica statistica e complessità', a research grant funded by PRIN 2015 (Agreement no. 2015K7KK8L).

The evaluation of the one-site partition function in the Rfmp follows the same type of arguments adopted in the literature for the Rmp [3, 22, 23] and for the Rtsp [39]. The one-site partition function z in the Rfmp is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{zone} z[\bQ]:=& \prod_{a=1}^n \left[\int_{0}^{2\pi}\frac{{\rm e}^{-2{\rm i}\lambda^a}\dd\lambda^a}{2\pi}\right] \nonumber \\ &\times \exp\left[\sum_{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}Q_{\alpha,\beta}\exp\left({\rm i}\sum_{a\in\alpha}\lambda^a+2{\rm i}\sum_{b\in\beta}\lambda^b\right)+\frac{\varepsilon-1}{4N}\sum_{\alpha\neq\emptyset}g_{|\alpha|}\exp\left(2{\rm i}\sum_{a\in\alpha}\lambda^a\right) \right]\nonumber \\ \equiv &~ z_0[\bQ]+\frac{\varepsilon-1}{2N}\sum_{\alpha\neq\emptyset}g_{|\alpha|}\frac{\partial z_0[\bQ]}{\partial Q_{\emptyset,\alpha}}. \nonumber \end{align} \tag{ A.1 }$$

In the previous equation, we have introduced $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} z_0[\bQ]$, that coincides with the expression given in equation (23e). It follows that, for $N\gg 1$, $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} \ln z[\bQ]$ provides a contribution both to the leading term $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} S[\bQ]$ and to the finite-size corrections $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} \Delta S^{\bf T}$$\newcommand{\bQ}{\boldsymbol{\mathsf Q}}[\bQ]$, namely, up to $o\left({1}/{N}\right)$ terms,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle 2\ln z[\bQ]=2\ln z_0[Q]-\frac{\varepsilon-1}{2N}\sum_{\alpha\neq\emptyset}\frac{g_{|\alpha|}}{z_0[\bQ]}\frac{\partial z_0[\bQ]}{\partial Q_{\emptyset,\alpha}}, \nonumber \end{align} \tag{ A.2 }$$

to be compared with the terms appearing in equation (23). The evaluation of $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} z_0[\bQ]$ is nontrivial in general. In the replica symmetric hypothesis, $\newcommand{\bQ}{\boldsymbol{\mathsf Q}} z_0[\bQ]$ can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle z_0[\bQ]=\prod_{a=1}^n \left[\int_{0}^{2\pi}\frac{{\rm e}^{-2{\rm i}\lambda^a}\dd\lambda^a}{2\pi}\right] \exp\left[\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}Q_{|\alpha|,|\beta|}\exp\left({\rm i}\sum_{a\in\alpha}\lambda^a+2{\rm i}\sum_{b\in\beta}\lambda^b\right)\right]. \nonumber \end{align} \tag{ A.3 }$$

Denoting now by $\xi_a:= {\rm e}^{{\rm i}\lambda_a}$, observe now that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\sum_{{{\alpha\cap\beta=\emptyset}\atop{\alpha\cup\beta\neq\emptyset}}}Q_{|\alpha|,|\beta|}\prod_{a\in\alpha}\xi_a\prod_{b\in\beta}\xi_b^2=\sum_{{p+q\geqslant 1}} Q_{p,q}\sum_{|\beta|=q}\prod_{b\in\beta}\xi_b^2\sum_{{{|\alpha|=p}\atop{\alpha\cap\beta=\emptyset}}}\ \prod_{a\in\alpha}\xi_a \nonumber \\ &\to \sum_{{p+q\geqslant 1}} Q_{p,q}\sum_{|\beta|=q}\prod_{b\in\beta}\xi_b^2\sum_{{|\alpha|=p}}\ \prod_{a\in\alpha}\xi_a \to \sum_{{p+q\geqslant 1}}\frac{Q_{p,q}}{q!} \left(\sum_{|\alpha|=p}\prod_{a\in\alpha}\xi_a\right)\left(\sum_b\xi_b^2\right)^q,\nonumber \end{align} \tag{ A.4 }$$

where each substitution is justified because of the overall constraint that allows us to neglect powers $\xi_b^{k}$ with $k\geqslant 3$. It can be seen that [39, 42]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\He}{{\rm He}} \displaystyle \sum_{|\alpha|=p}\prod_{a\in\alpha}\xi_a=\frac{\left(\sum_{a}\xi_a^2\right)^\frac{p}{2}}{p!}\He_p\left(\frac{\sum_a\xi_a}{\sqrt{\sum_a\xi_a^2}}\right). \nonumber \end{align} \tag{ A.5 }$$

Here $\newcommand{\He}{{\rm He}} \He_p(x)$ is the probabilists' Hermite polynomial [42]. Substituting the previous identity in the expression for z₀, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\He}{{\rm He}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle z_0[\bQ] &=\int_{-\infty}^{+\infty}\frac{\dd x\dd k_x \dd y\dd k_y}{(2\pi)^2} \exp\left[\sum_{p+q\geqslant 1}\frac{Q_{p,q}}{p!q!}y^{\frac{p}{2}+q}\He_p\left(\frac{x}{\sqrt{y}}\right)+{\rm i}k_x x+{\rm i}k_y y\right]\Phi(k_x,k_y), \nonumber \end{align} \tag{ A.6 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \Phi(k_x,k_y)=\left[\int_{0}^{2\pi}\frac{\dd\lambda}{2\pi}\exp\left(-2{\rm i}\lambda-{\rm i}k_x {\rm e}^{{\rm i}\lambda}-{\rm i}k_y {\rm e}^{2{\rm i}\lambda}\right)\right]^n=\left(-\frac{k_x^2}{2}-{\rm i}k_y\right)^n. \nonumber \end{align} \tag{ A.7 }$$

Using now the identity $a^n=\left.\partial_t^n{\rm e}^{at}\right|_{t=0}$, we can write [3, 28, 39]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \newcommand{\He}{{\rm He}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \label{zfractional} &\frac{\ln z_0[\bQ]}{n}=\frac{1}{n}\ln\left[\left.\frac{\partial^n}{\partial t^n}\int_{-\infty}^{+\infty}\frac{{\rm e}^{-\frac{z^2}{2}}\dd z}{\sqrt{2\pi}} \exp\left(\sum_{p+q\geqslant 1}\frac{Q_{p,q}}{p!q!}t^{\frac{p}{2}+q}\He_p(z)\right)\right|_{t=0}\right]\nonumber \\ &=\int_0^{+\infty}\frac{\dd t}{t}\left[{\rm e}^{-t}-\int_{-\infty}^{+\infty}\frac{{\rm e}^{-\frac{z^2}{2}}\dd z}{\sqrt{2\pi}} \exp\left(\sum_{p+q\geqslant 1}\frac{Q_{p,q}}{p!q!}(-t)^{\frac{p}{2}+q}\He_p(z)\right)\right]+O(n). \nonumber \end{align} \tag{ A.8 }$$

At this point, some considerations are in order. Let us first observe that, on the matching saddle-point in equation (30), $Q_{p, q}^{\rm sp}=\delta_{p, 0}Q_q$, the previous expression becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \lim_{n\to 0}\frac{\ln z_0[Q]}{n}=\int_0^{+\infty}\frac{\dd t}{t}\left[{\rm e}^{-t}- \exp\left(\sum_{q=1}^\infty\frac{Q_{q}}{q!}(-t)^{q}\right)\right]. \nonumber \end{align} \tag{ A.9 }$$

The expression above coincides with the one-site partition function in the Rmp and it has been discussed in details in [3, 22, 23, 25]. In particular, equation (31a) can be obtained introducing the function G in equation (31b). On the other hand, if we consider the Rtsp saddle-point solution, $\tilde Q^{\rm sp}_{q, p}=\delta_{q, 0}\tilde Q_p$, the Rtsp one-site partition function is recovered, as it can be easily seen comparing equation (A.8) with the results in [39].

To evaluate the $\Delta S^{\rm T}$ contribution on the matching saddle-point, we follow the approach in [22, 23]. In particular,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle &\frac{1}{n}\Delta S^{\rm T}[\bQ^{\rm sp}]=\frac{1}{Nn}\sideset{}{'}\sum_{{\alpha,\beta|\alpha,\hat\beta}}\frac{g_{|\alpha|+2|\beta|}g_{|\alpha|+2|\hat \beta|}}{4 g^2_{2|\beta|+2|\hat\beta|+2|\alpha|}}Q_{|\beta|+|\hat\beta|+|\alpha|}^2-\frac{\varepsilon-1}{4nN}\sum_{\alpha\neq\emptyset}\frac{g_{|\alpha|}}{g_{2|\alpha|}}Q_{0,|\alpha|}\nonumber \\ &\qquad =\frac{1}{N}\sum_{{{s+p\geqslant 1}\atop{s+q\geqslant 1}}}\frac{(-1)^{\,p+q+s-1}\Gamma(s+p+q)}{p!q!s!}\frac{g_{s+2p}g_{s+2q}}{4 g^2_{2p+2q+2s}}Q_{p+q+s}^2 -\frac{\varepsilon-1}{4N}\sum_{p=1}^\infty\frac{(-1)^{\,p-1}}{p}\frac{g_{p}}{g_{2p}}Q_{p}. \nonumber \end{align} \tag{ B.1 }$$

The asymptotic value of the last sum is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \lim_{\beta\to +\infty}&\frac{\varepsilon-1}{4N\beta}\sum_{p=1}^\infty\frac{(-1)^{\,p-1}}{p}\frac{g_{p}}{g_{2p}}Q_{p} =\lim_{\beta\to +\infty}\frac{\varepsilon-1}{2N\beta}\sum_{p=1}^\infty\frac{(-1)^{\,p-1}}{p p!}\int_{-\infty}^{+\infty}{\rm e}^{\beta p x-\hat G(x)}\dd x \nonumber \\ &=\frac{\varepsilon-1}{2N}\lim_{\beta\to +\infty}\int_{-\infty}^{+\infty}\frac{\gamma_{\rm E}+\beta x+\Gamma(0,{\rm e}^{\beta x})}{\beta}{\rm e}^{-\hat G(x)}\dd x\nonumber \\ & =\frac{\varepsilon-1}{2N}\int_{0}^{+\infty}x{\rm e}^{-\hat G(x)}\dd x=\frac{\varepsilon-1}{4N}\zeta(2). \nonumber \end{align} \tag{ B.2 }$$

In the previous expression, $\gamma_E$ is the Euler–Mascheroni constant, whereas $\newcommand{\dd}{{\rm d}} \Gamma(a, z):= \int_z^{\infty}{\rm e}^{-t}t^{a-1}\dd t$ is the incomplete gamma function [42]. In the contribution from $\Delta S^{\rm T}$ we also have

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle &\frac{1}{n}\sum_{{s+p\geqslant 1}\atop{s+q\geqslant 1}}^\infty\binom{n}{s\ p\ q}\frac{g_{s+2p}g_{s+2q}}{4 g_{2p+2q+2s}}Q_{p+q+s}^2\nonumber \\ &=\sum_{{{s+p\geqslant 1}\atop{s+q\geqslant 1}}}\frac{(-1)^{\,p+q+s-1}\Gamma(s+p+q)}{p!q!s!}\frac{g_{s+2p}g_{s+2q}}{4 g^2_{2p+2q+2s}}Q_{p+q+s}^2\nonumber \\ & =\sum_{{{s+p\geqslant 1}\atop{s+q\geqslant 1}}}\frac{(-1)^{\,p+q+s-1}g_{s+2p}g_{s+2q}}{2p!q!s!g_{2p+2q+2s}}Q_{p+q+s}\int_{-\infty}^{+\infty}{\rm e}^{(\,p+q+s) x-G(x)}\dd y =\frac{1}{\beta}\int_{-\infty}^{+\infty}\dd x\, {\rm e}^{-G(x)}\nonumber \\ &\times \sum_{k=1}^\infty (-1)^{k-1}{\rm e}^{k x}Q_k\left[\sum_{p=1}^{k-1}\frac{1}{2pp!(k-p)!}+\sum_{s=1}^k\sum_{p=0}^{k-s}\frac{1}{s!(s+2p)\,p!(k-s-p)!}\right]. \nonumber \end{align} \tag{ B.3 }$$

In [22] it has been proved that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \lim_{\beta\to+\infty}\frac{1}{\beta^2}\int_{-\infty}^{+\infty}\dd x\, {\rm e}^{-G(x)} \sum_{k=1}^\infty (-1)^{k-1}{\rm e}^{k x}Q_k \sum_{p=1}^{k-1}\frac{1}{2pp!(k-p)!}=-\frac{\zeta(3)}{N}. \nonumber \end{align} \tag{ B.4 }$$

The remaining contribution is zero. To prove this fact, let us start from

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle & \sum_{s=1}^k\sum_{p=0}^{k-s}\frac{1}{s!(s+2p)\,p!(k-s-p)!}=\sum_{s=1}^k\frac{\rm i}{2\pi s!}\int_0^{+\infty}\dd t\oint_{\gamma_{\mathrm H}}\dd z \sum_{p=0}^\infty(-z)^{\,p+s-k-1}\frac{{\rm e}^{-(s+2p)t-z}}{p!} \nonumber \\ & =\sum_{s=1}^k\frac{\rm i}{2\pi s!}\int_0^{+\infty}{\rm e}^{-st}\dd t\oint_{\gamma_{\mathrm H}}{\rm e}^{-(1+{\rm e}^{-2t})z}(-z)^{s-k-1}\dd z =\sum_{s=1}^k\int_0^{+\infty}\frac{\left(1+{\rm e}^{-2t}\right)^{k-s-1}{\rm e}^{-st}}{\Gamma(k-s+1) s!}\dd t \nonumber \\ &\qquad =\int_0^{1}\sum_{s=1}^k\frac{\left(1+\tau\right)^{k-s}\tau^{{s}/{2}-1}}{2s!(k-s)!}\dd\tau=\int_0^1\dd \tau\frac{\left(\tau+\sqrt \tau+1\right)^k-(1+\tau)^k}{2k!\tau}. \nonumber \end{align} \tag{ B.5 }$$

In the expression above, $\gamma_{\mathrm H}$ is the Hankel contour in the complex plane. Summing over k we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle &\frac{1}{2\beta}\int_{-\infty}^{+\infty}{\rm e}^{-G(x)}\dd x \int_0^1\dd \tau\frac{G\left(x+\ln\left(1+\sqrt \tau+\tau\right)\right)-G\left(x+\ln(1+\tau)\right)}{\tau} \nonumber \\ &=\frac{\beta}{2}\int_{-\infty}^{+\infty}{\rm e}^{-\hat G(x)}\dd x \int_0^{+\infty}\dd u\left[\hat G\left(x + \frac{\ln\left(1 + {\rm e}^{-\frac{\beta u}{2}} + {\rm e}^{-\beta u}\right)}{\beta}\right) - \hat G\left(x + \frac{\ln(1 + {\rm e}^{-\beta u})}{\beta}\right)\right], \nonumber \end{align} \tag{ B.6 }$$

whose corresponding contribution goes to zero as $\beta\to\infty$. We can finally write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bQ}{\boldsymbol{\mathsf Q}} \displaystyle \lim_{\beta\to+\infty}\lim_{n\to 0}\frac{\Delta S^{\rm T}[\bQ^{\rm sp}]}{n\beta}=-\frac{1}{N}\left(\frac{\varepsilon-1}{4}\zeta(2)+\zeta(3)\right). \nonumber \end{align} \tag{ B.7 }$$

In this appendix we will give some details about the evaluation of the logarithm of the determinant of the Hessian matrix $\boldsymbol\Omega$ on the matching saddle-point using the replica approach, showing that it is equal to the Rmp contribution, plus an additional contribution that we expect to be opposite from Wästlund's formula. We start from its general expression in equation (26),

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle & \Omega_{\alpha\beta,\hat\alpha\hat\beta}[\boldsymbol{\mathsf Q}^{\rm sp}]:= \sqrt{g_{|\alpha|+2|\beta|}g_{|\hat\alpha|+2|\hat\beta|}}\left.\frac{\partial^2 S[\boldsymbol{\mathsf Q}]}{\partial Q_{\alpha,\beta}\partial Q_{\hat \alpha,\hat \beta}}\right|_{\boldsymbol{\mathsf Q}=\boldsymbol{\mathsf Q}^{\rm sp}}\nonumber \\ &\qquad =\delta_{\alpha,\hat\alpha}\delta_{\beta,\hat\beta}{\mathbb I}(\alpha\cap\beta =\emptyset) -2\sqrt{g_{|\alpha|+2|\beta|}g_{|\hat\alpha|+2|\hat\beta|}}\left\langle\exp\left({\rm i}\sum_{{a\in\alpha\cup\hat\alpha}}\lambda^a+2{\rm i}\sum_{{b\in\beta\cup\hat\beta}}\lambda^b\right)\right\rangle_{z_0}\nonumber \\ & +2\sqrt{g_{|\alpha|+2|\beta|}g_{|\hat\alpha|+2|\hat\beta|}}\left\langle\exp\left({\rm i}\sum_ {{a\in\alpha}}\lambda^a+2{\rm i}\sum_{b\in\beta}\lambda^b\right)\right\rangle_{z_0}\left\langle\exp\left({\rm i}\sum_{a\in\hat\alpha}\lambda^a+2{\rm i}\sum_{b\in\hat\beta}\lambda^b\right)\right\rangle_{z_0}\nonumber \\ &\qquad \qquad \qquad \qquad =\delta_{\alpha,\hat\alpha}\delta_{\beta,\hat\beta}+\frac{1}{2}\frac{Q_{\alpha,\beta}}{\sqrt{g_{|\alpha|+2|\beta|}}}\frac{Q_{\hat\alpha,\hat\beta}}{\sqrt{g_{|\hat\alpha|+2|\hat\beta|}}}\nonumber \\ & -\frac{\sqrt{g_{|\alpha|+2|\beta|}g_{|\hat\alpha|+2|\hat\beta|}}Q_{\alpha\triangle\hat\alpha,\beta\cup\hat\beta\cup(\alpha\cap\hat\alpha)}}{g_{|\alpha\triangle\hat\alpha|+2|\beta\cup\hat\beta\cup(\alpha\cap\hat\alpha)|}}{\mathbb I}(\beta\cap\hat\beta=\emptyset){\mathbb I}((\alpha\cup\hat\alpha)\cap(\beta\cup\hat\beta)=\emptyset). \nonumber \end{align} \tag{ C.1 }$$

In the replica symmetric ansatz on the matching saddle-point solution the expression greatly simplifies, becoming

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega_{\alpha\beta,\hat\alpha\hat\beta}=&~\delta_{\alpha,\hat\alpha}\delta_{\beta,\hat\beta}{\mathbb I}(\beta\cap \alpha=\emptyset)+\frac{\delta_{|\alpha|,0}\delta_{|\hat \alpha|,0}}{2}\frac{Q_{|\beta|}}{\sqrt{g_{2|\beta|}}}\frac{Q_{|\hat\beta|}}{\sqrt{g_{2|\hat\beta|}}} \nonumber \\ &-\delta_{\alpha,\hat\alpha}\frac{Q_{|\beta\cup\hat\beta\cup\alpha|}\sqrt{g_{|\alpha|+2|\beta|}g_{|\alpha|+2|\hat\beta|}}{\mathbb I}(\beta\cap\hat\beta=\emptyset){\mathbb I}((\beta\cup\hat\beta)\cap \alpha=\emptyset)}{g_{2|\beta\cup\hat\beta\cup\alpha|}}.\nonumber \end{align} \tag{ C.2 }$$

We follow the approach of [22, 23]. Observe first that the quantity above is diagonal respect to the index α. In particular

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \ln\det\boldsymbol{\mathsf \Omega}[\boldsymbol{\mathsf Q}^{\rm sp}]&=\ln\det\boldsymbol{\mathsf \Omega}^{(0)}[\boldsymbol{\mathsf Q}^{\rm sp}]+\sum_{|\alpha|=1}^{\infty}\binom{n}{|\alpha|}\ln\det\boldsymbol{\mathsf \Omega}^{(|\alpha|)}[\boldsymbol{\mathsf Q}^{\rm sp}]\nonumber \\ & =\ln\det\boldsymbol{\mathsf \Omega}^{(0)}[\boldsymbol{\mathsf Q}^{\rm sp}]+n\sum_{s=1}^{\infty}\frac{(-1)^{s-1}}{s}\ln\det\boldsymbol{\mathsf \Omega}^{(s)}[\boldsymbol{\mathsf Q}^{\rm sp}]+o(n), \nonumber \end{align} \tag{ C.3 }$$

where we have separated the contributions for different values of $s=|\alpha|$ and introduced

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega^{(s)}_{\beta\hat\beta}:= \delta_{\beta,\hat\beta}-Q_{|\beta|+|\hat\beta|+s}{\mathbb I}(\beta\cap\hat\beta=\emptyset)\frac{\sqrt{g_{2|\beta|+s}g_{2|\hat\beta|+s}}}{g_{2|\beta|+2|\hat\beta|+2s}}+\frac{\delta_{s,0}}{2}\frac{Q_{|\beta|}Q_{|\hat\beta|}}{\sqrt{g_{2|\beta|}g_{2|\hat\beta|}}}, \nonumber \end{align} \tag{ C.4 }$$

where β is a set of n  −  s replica indices. A vector $\boldsymbol{\mathsf q}$ is eigenvector of $\boldsymbol{\mathsf \Omega}^{(s)}$ with eigenvalue λ if

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eigeneq} q_\beta-\sum_{{\hat\beta\colon \beta\cap\hat\beta=\emptyset}}\frac{\sqrt{g_{s+2|\beta|}g_{s+2|\hat\beta|}}}{g_{2|\beta|+2|\hat\beta|+2s}}Q_{|\beta|+|\hat\beta|+s}q_{\hat \beta} +\frac{\delta_{s,0}}{2}\sum_{|\hat\beta|\neq\emptyset}\frac{Q_{|\beta|}Q_{|\hat\beta|}q_{\hat\beta}}{\sqrt{g_{2|\beta|}g_{2|\hat\beta|}}}=\lambda q_\beta. \nonumber \end{align} \tag{ C.5 }$$

A first step is to diagonalize $\boldsymbol{\mathsf \Omega}^{(s)}$ according to the irreducible representations of the permutation group [43] in the space of n  −  s replica indices. In the spirit of the strategy of De Almeida and Thouless [40], and strictly following [22, 23], we observe that an eigenvector $\boldsymbol{\mathsf q}^{(c)}$ with c distinguished replicas, is such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_\beta^{(c)}:= \left\{\begin{array}{@{}clc@{}} 0&{\rm{if~}}|\beta|<c, \nonumber \\ \omega_{|\beta|}^i&{\rm{if~}} \beta {\rm{~contains~}} c-i {\rm{~of~the~}} c {\rm{~distinguished~indices,~with~}} i=0,1,\dots,c. \end{array}\right. \nonumber \end{align} \tag{ C.6 }$$

For c  =  0, than we only have n  −  s possible eigenvectors in the form $\newcommand{\e}{{\rm e}} q_\beta^{(0)}\equiv q^{(0)}_{|\beta|}$. The eigenvalue equation can be written down for the $(n-s)\times (n-s)$ matrix $\boldsymbol{\mathsf N}^{(s, 0)}$ given by

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle N^{(s,0)}_{pq}=\delta_{pq}-\binom{n-s-p}{q}\frac{\sqrt{g_{2p+s}g_{2q+s}}Q_{p+q+s}}{g_{2p+2q+2s}} +\frac{\delta_{s,0}}{2}\binom{n}{q}\frac{Q_pQ_q}{\sqrt{g_{2p}g_{2q}}}, \nonumber \end{align} \tag{ C.7 }$$

whose eigenvalues have multiplicity 1 in the set of eigenvalues of $\boldsymbol{\mathsf \Omega}^{(s)}$. For $c\geqslant 1$, imposing the ortogonality relation between $\boldsymbol{\mathsf q}^{(c)}$ and $\boldsymbol{\mathsf q}^{(c-1)}$, we obtain

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle &\sum_{\beta}q^{(c)}_\beta q^{(c-1)}_\beta=\sum_{|\beta|\geqslant c}q^{(c)}_\beta q^{(c-1)}_\beta \nonumber \\ & =\sum_{p\geqslant c}\sum_{j=0}^{c-1}\sum_{r=0}^{c-j}\binom{n-s-c}{p-j-r}\binom{c}{j}\binom{c-j}{r}\omega_{p}^{c-(r+j)}\omega_{p}^{c-1-j}=0 \nonumber \\ &\Longrightarrow \sum_{r=0}^{c-j}\binom{n-s-c}{p-j-r}\binom{c-j}{r}\omega_{p}^{c-(r+j)}=0{\rm{~for~}} p\geqslant c {\rm{~and~}} j=0,1,\dots,c-1,\label{orto} \nonumber \end{align} \tag{ C.8a }$$

that for $n\to 0$ becomes, for $p\geqslant c$ and $j=0, 1, \dots, c-1$,

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \sum_{r=0}^{c-j}(-1)^{r}\binom{c-j}{r}\frac{\Gamma(s+c+p-r-j)}{\Gamma(\,p-r-j+1)}\omega_{p}^{c-(r+j)}=0,\label{orto2} \nonumber \end{align} \tag{ C.8b }$$

where we have used the property

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \lim_{n\to 0}\binom{n-a}{b}=\frac{(-1)^b\Gamma(a+b)}{\Gamma(a)\Gamma(b+1)}. \nonumber \end{align} \tag{ C.9 }$$

Equation (C.8a) allows us to keep $\omega_{p}^0$ as independent only. In particular, for c  =  1 we have $p\omega_{p}^0+(n-s-p)\omega_{p}^1=0$ and therefore the diagonalization of $\boldsymbol{\mathsf \Omega}$ in the subspace c  =  1 can be reduced to the diagonalization of the $(n-s-1)\times(n-s-1)$ matrix

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle N^{(s,1)}_{pq}=&~\delta_{pq} -\binom{n-s-p}{q}\frac{q}{q+s-n}\frac{\sqrt{g_{2p+s}g_{2q+s}}Q_{p+q+s}}{g_{2p+2q+2s}} \nonumber \\ & +\frac{\delta_{s,0}}{2}\left[\binom{n-1}{q}\frac{q}{q-n}+\binom{n-1}{q-1}\right]\frac{Q_pQ_q}{\sqrt{g_{2p}g_{2q}}}, \nonumber \end{align} \tag{ C.10 }$$

with eigenvalue multiplicity n  −  s  −  1 respect to the original matrix $\boldsymbol{\mathsf \Omega}^{(s)}$. Before proceeding further, some considerations are in order. As observed in [22, 23], the matrices $\boldsymbol{\mathsf N}^{(0, 0)}$ and $\boldsymbol{\mathsf N}^{(0, 1)}$ have the same limit as $n\to 0$, in particular $\lim_{n\to 0}\boldsymbol{\mathsf N}^{(0, 1)}=\boldsymbol{\mathsf N}^{(0, 0)}$. The calculation of the contribution of these two matrices requires some care, but it can be proved that it is eventually zero for $\beta\to+\infty$ [22, 23]. For $c\geqslant 1$ equation (C.8b) implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\frac{\omega_{p}^0}{\prod_{u=0}^{c-1}(s+p+u)}=\frac{\omega_{p}^1}{(\,p-c+1)\prod_{u=1}^{c-1}(s+p+u)} =\cdots \nonumber \\ & =\cdots=\frac{\omega_{p}^l}{\prod_{v=1}^{l}(\,p-c+v)\prod_{u=l}^{c-1}(s+p+u)} =\cdots =\frac{\omega_{p}^{c}}{\prod_{v=1}^{c}(\,p-c+v)} \nonumber \\ & \Longrightarrow \omega^l_p=\frac{\Gamma(\,p+l-c+1)\Gamma(s+p)}{\Gamma(\,p-c+1)\Gamma(s+p+l)}\omega_p^0. \nonumber \end{align} \tag{ C.11 }$$

Using the previous result, equation (C.5) for the $|\beta|=0$ component of an eigenvector with $c\geqslant1$ becomes

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \lambda\omega_p^0=&~\frac{\delta_{s,0}}{2}\sum_{q=1}^\infty\sum_{i=0}^c\binom{n-c}{q+{\rm i}-c}\binom{c}{c-i}\frac{Q_pQ_q\omega_q^i}{\sqrt{g_{2p}g_{2q}}} \nonumber \\ &+\omega_p^0-\sum_{q=1}^\infty\binom{n-s-p}{q}\frac{\sqrt{g_{s+2p}g_{s+2q}}}{g_{2s+2p+2q}}Q_{p+q+s}\omega_{q}^c \nonumber \\ & =\omega_p^0-\sum_{q=1}^\infty (-1)^{q}\frac{\Gamma(s+q)\Gamma(s+p+q)\sqrt{g_{s+2p}g_{s+2q}}}{\Gamma(s+p)\Gamma(s+q+c)\Gamma(q-c+1)}\frac{Q_{p+q+s}}{g_{2s+2p+2q}}\omega_{q}^0+o(n) \nonumber \end{align} \tag{ C.12 }$$

(observe that the quadratic term is zero because of equation (C.8a)). We can finally write for each value of s

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \ln\det\boldsymbol{\mathsf \Omega}^{(s)}=\sum_{c=0}^\infty\left[\binom{n-s}{c}-\binom{n-s}{c-1}\right]\ln\det\boldsymbol{\mathsf N}^{(s,c)} \nonumber \end{align} \tag{ C.13 }$$

being the $(n-s-c)\times(n-s-c)$ matrix $\boldsymbol{\mathsf N}^{(s, c)}$ in the $n\to 0$ limit

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N^{(s,c)}_{pq}=\delta_{pq} -(-1)^q\frac{\Gamma(s+q)\Gamma(\,p+q+s)\sqrt{g_{2p+s}g_{2q+s}}}{\Gamma(s+p)\Gamma(q-c+1)\Gamma(s+q+c)}\frac{Q_{p+q+s}}{g_{2p+2q+2s}}. \nonumber \end{align} \tag{ C.14 }$$

Shifting by c the indices in the expression above, transposing and then multiplying by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (-1)^{\,p+q}\sqrt{\frac{g_{2q+2c+s}}{g_{2p+2c+s}}}\frac{\Gamma(q+c+s)\Gamma(\,p+1)}{\Gamma(\,p+c+s)\Gamma(q+1)} \nonumber \end{align*}$$

we get a new matrix $\boldsymbol{\mathsf M}^{(s, c)}$ with the same spectrum of $\boldsymbol{\mathsf N}^{(s, c)}$, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M^{(s,c)}_{pq}=\delta_{pq} -(-1)^{q+c}\frac{\Gamma(\,p+q+s+2c)}{\Gamma(\,p+s+2c)\Gamma(q+1)}\frac{g_{2q+2c+s}}{g_{2(\,p+q+s+2c)}}Q_{p+q+s+2c}. \nonumber \end{align} \tag{ C.15 }$$

To evaluate the correction to the free-energy we need to evaluate therefore

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{omegamsc} &\ln\det\boldsymbol{\mathsf \Omega}\nonumber \\ &=\sum_{c=2}^\infty\left[\binom{n}{c}-\binom{n}{c-1}\right]\ln\det\boldsymbol{\mathsf M}^{(0,c)}+\sum_{s=1}^\infty\binom{n}{s}\sum_{c=0}^\infty\left[\binom{n-s}{c} - \binom{n-s}{c-1}\right]\ln\det\boldsymbol{\mathsf M}^{(s,c)}\nonumber \\ &= n\sum_{c=2}^\infty(-1)^{c-1}\frac{2c-1}{c(c-1)}\ln\det\boldsymbol{\mathsf M}^{(0,c)}\nonumber \\ &\quad -n\sum_{s=1}^\infty \sum_{c=0}^\infty (-1)^{c+s}(2c+s-1)\frac{\Gamma(s+c-1)}{s!c!}\ln\det\boldsymbol{\mathsf M}^{(s,c)}+o(n). \nonumber \end{align} \tag{ C.16 }$$

At this point is important to observe that, in the expression above, the s  =  0 contribution exactly coincides with the fluctuation contribution appearing in the finite-size corrections of Rmp. Its evaluation has been performed in [22, 23], but no closed formula is known for it. The $s\geqslant1$ contribution is instead absent in the Rmp. Generalizing therefore the analysis of [22, 23], we note that eigenvalues of $\boldsymbol{\mathsf M}^{(s, c)}$ are the same as the eigenvalues of the operator

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal M}^{(s,c)}(x,y)=\delta(x-y)-(-1)^c{\mathcal A}^{(c+{s}/{2})}(x,y) \nonumber \end{align} \tag{ C.17 }$$

where the operator ${\mathcal A}^{(k)}(x, y)$ was the one introduced in [23] and it is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal A}^{(k)}(x,y)=2{\rm e}^{-\frac{G(x)+G(y)}{2}}\sum_{q=0}^\infty \frac{(-1)^q{\rm e}^{(q+k)(x+y)}}{\Gamma(2k+q)\Gamma(q+1)}g_{2(q+k)}. \nonumber \end{align} \tag{ C.18 }$$

Indeed, if $\psi_p$ is an eigenvector of $\boldsymbol{\mathsf M}^{(s, c)}$ with corresponding eigenvalue λ, then, by a straightforward computation, it can be verified that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi(x):= {\rm e}^{\left(c+\frac{s}{2}\right)x-\frac{G(x)}{2}}\sum_{p=0}^{+\infty}\frac{(-1)^{\,p}}{p!}{g_{2p+2c+s}}\psi_p{\rm e}^{\,px} \nonumber \end{align} \tag{ C.19 }$$

is an eigenvector of ${\mathcal M}^{(s, c)}(x, y)$ with the same eigenvalue. To evaluate the $\beta\to\infty$ limit of ${\mathcal A}^{(k)}(x, y)$ we observe that it has the same eigenvalues of the operator

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal H}^{(k)}(x,y):= \beta {\mathcal A}^{(k)}(\beta x,\beta y) ={\rm e}^{-\frac{\hat G(x)+\hat G(y)}{2}}\sum_{q=0}^\infty \frac{(-1)^q{\rm e}^{(q+k)\beta (x+y)}}{\Gamma(2k+q)\Gamma(q+1)}\frac{1}{q+k}. \nonumber \end{align} \tag{ C.20 }$$

In [22, 23] it has been shown that, for $\beta\to+\infty$, if we impose $\ln k=\beta t$ with t fixed, the $\beta\to+\infty$ limit exists. In particular

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal H}^{(k)}(x,y)\xrightarrow[t {\rm{~fixed}}]{\ln k=\beta t} {\mathcal H}_{t}(x,y)={\rm e}^{-\frac{\hat G(x)+\hat G(y)}{2}}\theta\left(x+y-2t\right). \nonumber \end{align} \tag{ C.21 }$$

This result suggests that the evaluation of the sums in equation (C.16) must be performed scaling c, s and β in a proper way. Given the known result for the Rmp, we distinguish now between the contributions with $s\geqslant 1$ and the contribution obtained for s  =  0. We know indeed that we can obtain a finite limit for the s  =  0, that is the corresponding fluctuation correction to the aoc in the Rmp. We can write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}{{\rm d}} \displaystyle \label{detrmp} &n\sum_{c=2}^\infty(-1)^{c-1}\frac{2c-1}{c(c-1)}\ln\det{\bf M}^{(0,c)}\nonumber \\ &~=n\sum_{E=1}^\infty \frac{1}{E}\sum_{c=2}^\infty(-1)^{(E+1)c}\frac{2c-1}{c(c-1)}{\rm tr}\left[\left({\mathcal A}^{(c)}\right)^E\right]\nonumber \\ &~=n\sum_{E=1}^\infty\frac{1}{E} \sum_{c=1}^\infty\frac{(4c-1){\rm tr}\left[\left({\mathcal A}^{(2c)}\right)^E\right]}{2c(2c-1)} -n\sum_{E=1}^\infty\frac{(-1)^E}{E}\sum_{c=1}^\infty\frac{(4c+1){\rm tr}\left[\left({\mathcal A}^{(2c+1)}\right)^E\right]}{2c(2c+1)}\nonumber \\ &~\xrightarrow{\beta\to +\infty}2n\beta\sum_{E{\rm{~odd}}}\frac{1}{E}\int_0^{+\infty}{\rm tr}\left[{\mathcal H}^E_t\right]\dd t. \nonumber \end{align} \tag{ C.22 }$$

A numerical estimation of the quantity above can be found in [23, 28]. The evaluation of the $s\geqslant 1$ terms is more complicated and we will present here a non-rigorous treatment. We have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{settoriE} &n\sum_{s=1}^\infty\sum_{c=0}^\infty (-1)^{c+s-1}(2c+s-1)\frac{\Gamma(s+c-1)}{s!c!}\ln\det\boldsymbol{\mathsf M}^{(s,c)}\nonumber \\ &~=n\sum_{E=1}^\infty\frac{1}{E}\sum_{s=1}^\infty\sum_{c=0}^\infty (-1)^{(E+1)c+s}(2c+s-1)\frac{\Gamma(s+c-1)}{s!c!}{\rm tr}\left[\left({\mathcal A}^{(c+\frac{s}{2})}\right)^E\right]+o(n)\nonumber \\ &~= n\sum_{E=1}^\infty\frac{1}{E}\sum_{s=1}^\infty\sum_{c=0}^\infty (-1)^{s}(4c+s-1)\frac{\Gamma(s+2c-1){\rm tr}\left[\left({\mathcal A}^{(2c+\frac{s}{2})}\right)^E\right]}{s!(2c)!}\nonumber \\ &\quad -n(-1)^{E}\sum_{E=1}^\infty\frac{1}{E}\sum_{s=1}^\infty\sum_{c=0}^\infty (-1)^{s}(4c+s+1)\frac{\Gamma(s+2c){\rm tr}\left[\left({\mathcal A}^{(2c+1+\frac{s}{2})}\right)^E\right]}{s!(2c+1)!}+o(n). \nonumber \end{align} \tag{ C.23 }$$

Introducing 4c  +  s  =  z, the first sum in the last line of equation (C.23) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{s=1}^\infty\sum_{c=0}^\infty (-1)^{s}(4c+s-1)\frac{\Gamma(s+2c-1){\rm tr}\left[\left({\mathcal A}^{(2c+\frac{s}{2})}\right)^E\right]}{s!(2c)!}\nonumber \\ =\sum_{z=1}^\infty (z-1){\rm tr}\left[\left({\mathcal A}^{({z}/{2})}\right)^E\right]\sum_{s=1}^\infty\sum_{c=0}^\infty \frac{(-1)^{s}\Gamma\left(s+2c-1\right){\mathbb I}\left(s+4c=z\right)}{\Gamma(s+1)\Gamma(2c+1)} ={\rm tr}\left[\left({\mathcal A}^{(1)}\right)^E\right]\nonumber \\ +\sum_{z=3}^\infty (z-1){\rm tr}\left[\left({\mathcal A}^{({z}/{2})}\right)^E\right]\sum_{s=1}^\infty\sum_{c=0}^\infty \frac{(-1)^{s}\Gamma\left(s+2c-1\right){\mathbb I}\left(s+4c=z\right)}{\Gamma(s+1)\Gamma(2c+1)}. \nonumber \end{align} \tag{ C.24a }$$

Similarly, after the change of variable z  =  4c  +  s  +  2, the second sum becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle & \sum_{s=1}^\infty\sum_{c=0}^\infty (-1)^{s}(4c+s+1)\frac{\Gamma(s+2c)}{s!\Gamma(2c+2)}\tr\left[\left({\mathcal A}^{(2c+\frac{s}{2}+1)}\right)^E\right]\nonumber \\ &=\sum_{z=3}^\infty (z-1)\tr\left[\left({\mathcal A}^{({z}/{2})}\right)^E\right]\sum_{s=1}^\infty\sum_{c=0}^\infty \frac{(-1)^{s}\Gamma(s+2c){\mathbb I}\left(2+s+4c=z\right)}{\Gamma(s+1)\Gamma(2c+2)}. \nonumber \end{align} \tag{ C.24b }$$

In other words, we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \sum_{s=1}^\infty&\sum_{c=0}^\infty (-1)^{c+s-1}(2c+s-1)\frac{\Gamma(s+c-1)}{s!c!}\ln\det\boldsymbol{\mathsf M}^{(s,c)}\nonumber \\&= \sum_{E=1}^\infty\frac{1}{E}\sum_{z=3}^\infty (z-1)\tr\left[\left({\mathcal A}^{({z}/{2})}\right)^E\right]\left[h_1(z)-(-1)^E h_2(z)\right]. \nonumber \end{align} \tag{ C.25 }$$

The main difficulty in the evaluation of the quantities above is that the coefficients

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_1(z) := \sum_{s=1}^\infty\sum_{c=0}^\infty \frac{(-1)^{s}\Gamma\left(s+2c-1\right){\mathbb I}\left(s+4c=z\right)}{\Gamma(s+1)\Gamma(2c+1)} \nonumber \end{align} \tag{ C.26 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_2(z) := \sum_{s=1}^\infty\sum_{c=0}^\infty \frac{(-1)^{s}\Gamma(s+2c){\mathbb I}\left(2+s+4c=z\right)}{\Gamma(s+1)\Gamma(2c+2)} \nonumber \end{align} \tag{ C.27 }$$

are oscillating with diverging amplitude in z for $z\to +\infty$, and therefore the large z estimation is not straightforward. It is possible, however, that a different rearrangement of the contributions appearing in the sums might lead to a simpler asymptotic evaluation.

In the present work, we avoided this estimation using Wästlund's formula, but it is interesting to observe that equation (4c) implies that the contribution from equation (C.23) is equal and opposite to the one in equation (C.22), that coincides with the fluctuation finite-size correction to the aoc in the Rmp. A more comprehensive study of these quantities can be matter of future investigations.