Instance Independence of Single Layer Quantum Approximate Optimization Algorithm on Mixed-Spin Models at Infinite Size

The mixed-spin SK model (often called the mixed $p$-spin model, although we will not use this terminology to avoid confusion with the $p$ of qaoa) is given by [15, 16]

$\displaystyle\begin{split}H&=\sum_{j=1}^{n}J_{j}Z_{j}+\frac{1}{\sqrt{n}}\sum_{1\leq j<k\leq n}J_{jk}Z_{j}Z_{k}+\frac{1}{n}\sum_{1\leq j<k<\ell\leq n}J_{jk\ell}Z_{j}Z_{k}Z_{\ell}+\cdots+\frac{1}{n^{\frac{d-1}{2}}}\sum_{1\leq j_{1}<\cdots j_{d}\leq n}J_{j_{1}\cdots j_{d}}Z_{j_{1}}\cdots Z_{j_{d}}\\ &=\sum_{q=1}^{d}n^{\frac{1-q}{2}}\sum_{\begin{subarray}{c}S\subseteq\{1,\dots,n\}\\ |S|=q\end{subarray}}J_{S}Z_{S}\end{split}$

(5)

where in the last line we denote the product $\prod_{i\in S}Z_{i}$ by $Z_{S}$. We assume each $J_{S}$ is sampled independently from a Gaussian distribution $N(0,\sigma^{2}_{|S|})$ that only depends on $|S|$, and we will let $\sigma_{q}$ be the standard deviation of the coupling constants $J_{S}$ with $|S|=q$, for $q=1,\dots,d$.

The ground state of the mixed-spin SK Hamiltonian is known to have a fixed energy per spin in the limit $n\rightarrow\infty$ [18]. In fact we can also allow arbitrarily high orders in the mixed-spin SK model ($d=\infty$), provided the variances decrease quickly enough to make $\sum_{q}2^{q}\sigma_{q}^{2}$ finite, and the ground state model will still have a fixed energy per spin as $n\rightarrow\infty$ [18]. However, for simplicity we consider some finite bound $d$ on the degree.

Our main result is as follows. Define the $n$-spin model by Eq. 5 with $J_{S}\sim N(0,\sigma^{2}_{|S|})$. Then, using depth $1$ qaoa with angles $\beta,\gamma$, the expectation of the energy per spin equals in the large $n$ limit is given by

$\displaystyle\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]$

$\displaystyle=2\gamma\sum_{q=1}^{d}\sigma_{q}^{2}\sum_{\begin{subarray}{c}a\text{ odd}\\ a\leq q\end{subarray}}\frac{(-1)^{\frac{a-1}{2}}}{a!(q-a)!}\exp\left(-2\gamma^{2}a\sum_{q^{\prime}=1}^{d}\frac{\sigma_{q^{\prime}}^{2}}{(q^{\prime}-1)!}\right)\sin^{a}(2\beta)\cos^{q-a}(2\beta),.$

(6)

If we define the variables $c_{q}:=\sigma_{q}/\sqrt{q!}$ and the polynomial $\xi(x):=\sum_{q=1}^{d}c_{q}^{2}x^{q}$, as is common in discussions of the mixed-spin model (see Section 3.2), we can also write this as

$\displaystyle\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]=2\gamma\Im\left[\xi\Big{(}\cos(2\beta)+i\sin(2\beta)\exp\left({-2\gamma^{2}\xi^{\prime}(1)}\right)\Big{)}\right].$

(7)

Furthermore, the second moment of the energy per spin equals

$\displaystyle\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{2}\rangle\right]$

$\displaystyle=\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]^{2}.$

(8)

Our second result, Eq. 8, allows us to prove that $p=1$ qaoa applied to the mixed-spin SK model has both concentration of measure and instance independence. To see this, we note that we may write

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{2}\rangle\right]-\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]^{2}$

$\displaystyle=\Big{(}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{2}\rangle-\langle H/n\rangle^{2}\right]\Big{)}+\Big{(}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle^{2}\right]-\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]^{2}\Big{)}$

(9)

When applying qaoa to mixed-spin SK models, the measured value of $(H/n)$ varies for two reason. First, it varies because the bonds $J_{S}$ vary from instance to instance (the $\operatorname*{\mathbb{E}}_{J}\left[\cdot\right]$ expectation). Second, it varies because the qaoa state $|\gamma,\beta\rangle$ is not an eigenstate of the $Z$-operators, so that the measurement outcomes have randomness even for fixed $J_{S}$ (the $\langle\cdot\rangle$ expectation). The left hand side of Eq. 9 represents the total variance in $(H/n)$ due to both sources of randomness. The right hand side demonstrates that the total variance can be decomposed into two terms. The first term is the average over $J_{S}$ of the variance due only to the measurement randomness. The second term is the variance in the expected value $\langle H/n\rangle$ due to the randomness in the bonds $J_{S}$. Since both of these terms are non-negative, they must both tend to zero as $n\rightarrow\infty$ as well.

The fact that the first variance approaches zero gives us concentration of measure: it says that for typical couplings $J_{S}$ the variance in the measurement outcomes vanishes, and thus we always measure a string with energy per spin equal to the expectation value (note that the term inside the $\operatorname*{\mathbb{E}}\left[\cdot\right]$ is always positive, so that for the average over $J_{S}$ to go to zero the magnitude of the typical value must also go to zero). The second term approaching zero clearly gives instance independence of the expectation value, since it shows that the variance in the expectation value due to different couplings vanishes.

Finally, we note that the methods we use also suffice to derive a formula for all higher moments of the energy per spin

$\displaystyle\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{m}\rangle\right]$

$\displaystyle=\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]^{m},$

(10)

although unlike the $m=2$ result, the $m>2$ result does not have any obvious implication for the performance of the qaoa.

As an example of using our Eqs. 6 and 7, we can compute the optimal angles for pure $d$-spin models given by $\sigma_{q}=\delta_{d,q}\sqrt{d!/2}$. In this case, our central equation reduces to

$\displaystyle\begin{split}\lim_{n\rightarrow\infty}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]&=\gamma d!\sum_{\begin{subarray}{c}a\text{ odd}\\ a\leq d\end{subarray}}\frac{(-1)^{\frac{a-1}{2}}}{a!(d-a)!}\exp\left(-2\gamma^{2}ad\right)\sin^{a}(2\beta)\cos^{d-a}(2\beta)\\ &=\gamma\Im\left[\Big{(}\cos(2\beta)+i\sin(2\beta)\exp\left({-2\gamma^{2}d}\right)\Big{)}^{d}\right].\end{split}$

(11)

Examples of this energy landscape for small $d$ are plotted in Fig. 1, and numerically optimized angles and corresponding energy per spin are plotted in Fig. 2.

As another example application of our formula, we explore the instance independence at finite $n$ to demonstrate convergence as $n\rightarrow\infty$. In Fig. 3a-e, we plot the expected energy per spin $\langle H/n\rangle$ for five randomly generated mixed-spin optimization problems generated from a distribution with $(\sigma_{1},\sigma_{2},\sigma_{3})=(1/3,1/2,1)$ (and all higher-order terms zero). We see that as we increase $n$, the energy landscape of any given problem instance quickly approaches the $n\rightarrow\infty$ landscape (Fig. 3f), as predicted.

To put Eq. 6 in perspective we explicitly calculate what it implies for two basic cases $q=2$ and $q=3$ and compare it with the literature on mixed-spin models.

The standard Sherrington-Kirkpatrick model has $q=2$, $\sigma_{2}=1$, and $\sigma_{q}=0$ for all other $q\neq 2$. In this case our Eq. 6 tells us that

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H_{SK}/n\rangle\right]$

$\displaystyle=2\gamma e^{-2\gamma^{2}}\sin(2\beta)\cos(2\beta)=\gamma e^{-2\gamma^{2}}\sin(4\beta),$

(12)

which agrees with previous work in [2] and [8]. This expectation is minimized by the parameters $\beta_{*}=\pm\pi/8\approx\pm 0.392699$ and $\gamma_{*}=\mp 0.5$ for which $\operatorname*{\mathbb{E}}\left[\langle H_{SK}/n\rangle\right]=-1/\sqrt{4e}\approx-0.303265$. This is the same result we found numerically for $d=2$ in Fig. 2.

For the more general case $q\geq 2$ we will briefly review the notation used by the articles [5, 4, 18, 17]. In the current paper the mixed-spin model is described by the Hamiltonian

$\displaystyle H$

$\displaystyle=\sum_{q=1}^{d}\frac{1}{\sqrt{n^{q-1}}}\sum_{1\leq j_{1}<\cdots<j_{q}\leq n}{J_{j_{1}\cdots j_{q}}Z_{j_{1}}\cdots Z_{j_{q}}},$

(13)

with $J_{j_{1}\cdots j_{q}}\sim N(0,\sigma_{q}^{2})$ and for a fixed $q$ the summation has $\binom{n}{q}$ terms. In [5] and elsewhere a different notation is used where the Hamiltonian is defined by

$\displaystyle H^{\prime}$

$\displaystyle=\sum_{q=1}^{d}\frac{c_{q}}{\sqrt{n^{q-1}}}\sum_{1\leq j_{1},\dots,j_{q}\leq n}{J_{j_{1}\cdots j_{q}}Z_{j_{1}}\cdots Z_{j_{q}}},$

(14)

with $J_{j_{1}\cdots j_{q}}\sim N(0,1)$. Crucially, now the indices $j_{1},\dots,j_{q}$ can have repeated values and permutations are treated distinctly, hence for each $q$ the summation involves $n^{q}$ terms.

In the limit of large $n$ one can show that the $q$-plets $(j_{1},\dots j_{q})$ with repeated values can be ignored as their relative contribution decreases as a function $n$. For distinct indices there are $q!$ permutations to consider in $H^{\prime}$, hence we have a sum of $q!$ standard normal distributions: $J^{\prime}_{j_{1}\dots j_{q}}+\cdots+J^{\prime}_{j_{q}\cdots j_{1}}$, which is identical to a single normal distribution with variance $q!$. As a result we re-express $H^{\prime}$ as

$\displaystyle H^{\prime}$

$\displaystyle=\sum_{q=1}^{d}\frac{c_{q}\sqrt{q!}}{\sqrt{n^{q-1}}}\sum_{1\leq j_{1}<\dots<j_{q}\leq n}{J_{j_{1}\cdots j_{q}}Z_{j_{1}}\cdots Z_{j_{q}}},$

(15)

with $J\sim N(0,1)$. Thus when $\sigma_{q}=c_{q}\sqrt{q!}$ we have that $H$ and $H^{\prime}$ describe the same model. It is common to capture the different $c_{q}$ coefficients by the mixture function $\xi(x)=\sum_{q}c_{q}^{2}x^{q}$, such that the standard SK model has $\xi(x)=x^{2}/2$ (i.e. $c_{2}=1/\sqrt{2}$ and hence $\sigma_{2}=1$).

For the $3$-body model $H_{3}$, [4] analyzes the mixture $\xi(x)=c_{3}^{2}x^{3}=x^{3}/2$, such that $c_{3}=1/\sqrt{2}$. It was shown that the expected ground state energy-per-spin of this model equals $-0.8132\pm 0.0001$ In the notation of the current article this model is equivalent to $\sigma_{3}=\sqrt{3}$ and with this our Eq. 6 gives the expectation

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H_{3}/n\rangle\right]$

$\displaystyle=3\gamma e^{-3\gamma^{2}}\sin(2\beta)\cos^{2}(2\beta)-\gamma e^{-9\gamma^{2}}\sin^{3}(2\beta).$

(16)

This expectation is minimized to $\operatorname*{\mathbb{E}}_{J}[\langle H_{3}/n\rangle]\approx-0.270638$ by the angles $(\beta_{*},\gamma_{*})\approx(\pm 0.290003,\mp 0.430091)$, which are the solutions to

	$\displaystyle\beta_{*}$	$\displaystyle=-\frac{1}{4}\arccos\left(1-\frac{1}{9\gamma_{*}^{2}}\right)$		(17)
	$\displaystyle\exp(-6\gamma_{*}^{2})$	$\displaystyle=18\gamma_{*}^{2}-3$		(18)
	$\displaystyle\operatorname*{\mathbb{E}}_{J}[\langle H_{3}/n\rangle]$	$\displaystyle=\sqrt{4\gamma_{*}^{2}-\frac{2}{3}}.$		(19)

This is the same result we found numerically for $d=3$ in Fig. 2. We thus see how depth $p=1$ qaoa can approximate the ground state energy by a factor of $0.332806$. This $0.33$ approximation factor had been reported earlier by Zhou et al.[19]

Following [2], to evaluate $\operatorname*{\mathbb{E}}_{J}[\langle H/n\rangle]$ and $\operatorname*{\mathbb{E}}_{J}[\langle(H/n)^{2}\rangle]$ we use the moment-generating function $\operatorname*{\mathbb{E}}_{J}[\langle e^{i\lambda H/n}\rangle]$. From the moment-generating function, we can find the moments via

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{m}\rangle\right]$

$\displaystyle=\left.(-i\partial_{\lambda})^{m}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle e^{i\lambda H/n}\rangle\right]\right|_{\lambda=0}$

(21)

We can simplify the expectation inside the moment-generating function by inserting three complete sets of states:

$\displaystyle\begin{split}\langle e^{i\lambda H/n}\rangle&=\langle\Sigma|e^{i\gamma H}e^{i\beta B}e^{i\lambda H/n}e^{-i\beta B}e^{-i\gamma H}|\Sigma\rangle\\ &=\sum_{z,z^{\prime},z^{\prime\prime}\in\{\pm 1\}^{n}}\langle\Sigma|z\rangle\langle z|e^{i\gamma H}e^{i\beta B}|z^{\prime\prime}\rangle\langle z^{\prime\prime}|e^{i\lambda H/n}e^{-i\beta B}e^{-i\gamma H}|z^{\prime}\rangle\langle z^{\prime}|\Sigma\rangle\\ &=\frac{1}{2^{n}}\sum_{z,z^{\prime},z^{\prime\prime}\in\{\pm 1\}^{n}}\exp{\left(i\gamma\left[H(z)-H(z^{\prime})\right]+i\frac{\lambda H(z^{\prime\prime})}{n}\right)}\langle z|e^{i\beta B}|z^{\prime\prime}\rangle\langle z^{\prime\prime}|e^{-i\beta B}|z^{\prime}\rangle\\ &=\frac{1}{2^{n}}\sum_{z,z^{\prime},z^{\prime\prime}\in\{\pm 1\}^{n}}\exp{\left(i\sum_{q=1}^{d}n^{\frac{1-q}{2}}\sum_{|S|=q}J_{S}\left[\gamma(z_{S}-z^{\prime}_{S})+\frac{\lambda z^{\prime\prime}_{S}}{n}\right]\right)}\langle z|e^{i\beta B}|z^{\prime\prime}\rangle\langle z^{\prime\prime}|e^{-i\beta B}|z^{\prime}\rangle\\ &=\frac{1}{2^{n}}\sum_{z,z^{\prime},z^{\prime\prime}\in\{\pm 1\}^{n}}\exp{\left(i\sum_{q=1}^{d}n^{\frac{1-q}{2}}\sum_{|S|=q}z^{\prime\prime}_{S}J_{S}\left[\gamma(z_{S}-z^{\prime}_{S})+\frac{\lambda}{n}\right]\right)}\langle z|e^{i\beta B}|(+1)^{n}\rangle\langle(+1)^{n}|e^{-i\beta B}|z^{\prime}\rangle,\end{split}$

(22)

where in the last step we used that $\langle z|e^{i\beta B}|z^{\prime}\rangle=\langle zz^{\prime}|e^{i\beta B}|(+1)^{n}\rangle$ and made the replacements $z\mapsto zz^{\prime\prime}$ and $z^{\prime}\mapsto z^{\prime}z^{\prime\prime}$.

We will now treat the $J_{S}$ couplings as a random variable and consider the expectation $\operatorname*{\mathbb{E}}_{J}$ of the energy. We assume that the distribution is symmetric with respect to $J\leftrightarrow-J$, such that we can replace $z^{\prime\prime}_{S}J_{S}$ by $J_{S}$ to get

$\displaystyle\begin{split}\operatorname*{\mathbb{E}}_{J}\left[\langle e^{i\lambda H/n}\rangle\right]&=\frac{1}{2^{n}}\sum_{z,z^{\prime},z^{\prime\prime}\in\{\pm 1\}^{n}}\operatorname*{\mathbb{E}}_{J\sim N}\left[\exp{\left(i\sum_{q=1}^{d}n^{\frac{1-q}{2}}\sum_{|S|=q}z^{\prime\prime}_{S}J_{S}\left[\gamma(z_{S}-z^{\prime}_{S})+\frac{\lambda}{n}\right]\right)}\right]\langle z|e^{i\beta B}|(+1)^{n}\rangle\langle(+1)^{n}|e^{-i\beta B}|z^{\prime}\rangle\\ &=\sum_{z,z^{\prime}\in\{\pm 1\}^{n}}\operatorname*{\mathbb{E}}_{J\sim N}\left[\exp{\left(i\sum_{q=1}^{d}n^{\frac{1-q}{2}}\sum_{|S|=q}J_{S}\left[\gamma(z_{S}-z^{\prime}_{S})+\frac{\lambda}{n}\right]\right)}\right]\langle z|e^{i\beta B}|(+1)^{n}\rangle\langle(+1)^{n}|e^{-i\beta B}|z^{\prime}\rangle.\end{split}$

(23)

When we further assume that the $J_{S}$ variables are independent between different $S$ we can continue by

$\displaystyle\operatorname*{\mathbb{E}}_{J}\left[\langle e^{i\lambda H/n}\rangle\right]$

$\displaystyle=\sum_{z,z^{\prime}\in\{\pm 1\}^{n}}\prod_{q=1}^{d}\prod_{|S|=q}\operatorname*{\mathbb{E}}_{J_{S}}\left[\exp{\left(i\cdot n^{\frac{1-q}{2}}J_{S}\left[\gamma(z_{S}-z^{\prime}_{S})+\frac{\lambda}{n}\right]\right)}\right]\langle z|e^{i\beta B}|(+1)^{n}\rangle\langle(+1)^{n}|e^{-i\beta B}|z^{\prime}\rangle.$

(24)

Next we assume that the $J_{S}$ are normally distributed with a standard deviation that is the same for all sets $S$ of the same size, that is $J_{S}\sim N(0,\sigma^{2}_{|S|})$. We note that taking the expectation value of a Gaussian random variable $J$ with standard deviation $\sigma$ gives

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N(0,\sigma^{2})}\left[e^{icJ}\right]=\frac{1}{\sigma\sqrt{2\pi}}\int e^{-J^{2}/2\sigma^{2}+icJ}=e^{-c^{2}\sigma^{2}/2}$

(25)

so that our overall expression becomes

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle e^{i\lambda H/n}\rangle\right]$

$\displaystyle=\smashoperator[r]{\sum_{z,z^{\prime}\in\{\pm 1\}^{n}}^{}}\exp{\left(-\sum_{q=1}^{d}\frac{\sigma_{q}^{2}}{2n^{q-1}}\sum_{|S|=q}\left[\gamma^{2}(z_{S}-z^{\prime}_{S})^{2}+\frac{2\gamma\lambda(z_{S}-z^{\prime}_{S})}{n}+\frac{\lambda^{2}}{n^{2}}\right]\right)}\langle z|e^{i\beta B}|(+1)^{n}\rangle\langle(+1)^{n}|e^{-i\beta B}|z^{\prime}\rangle.$

(26)

To do the sum over $z,z^{\prime}$, we claim that the summand does not depend on all $2n$ spin values of $z$ and $z^{\prime}$. Instead, it is only a function of the four integer values $(n_{++},n_{+-},n_{-+},n_{--})$, where $n_{ss^{\prime}}$ is defined to be the number of positions $k\in\{1,\dots,n\}$ with $(z_{k},z^{\prime}_{k})=(s,s^{\prime})$. Note that only three of these variables are actually independent, as we always have $n_{++}+n_{+-}+n_{-+}+n_{--}=n$. As these numbers summarize the crucial information of the strings, we will refer to $(n_{++},n_{+-},n_{-+},n_{--})$ as the sketch of $(z,z^{\prime})$. Writing the summand in terms of the sketch rather than $(z,z^{\prime})$ was introduced in [2]; here we establish that we can still write the summand in terms of the sketch for the mixed-spin SK model. To start, it is straightforward to verify that

$\displaystyle\langle z|e^{i\beta B}|(+1)^{n}\rangle\langle(+1)^{n}|e^{-i\beta B}|z^{\prime}\rangle$

$\displaystyle=\smashoperator[r]{\prod_{ss^{\prime}\in\{\pm\}^{2}}^{}}Q_{ss^{\prime}}^{n_{ss^{\prime}}},\text{\leavevmode\nobreak\ with\leavevmode\nobreak\ }\begin{cases}Q_{++}:=\cos^{2}(\beta)\\ Q_{--}:=\sin^{2}(\beta)\\ Q_{-+}:=i\sin(\beta)\cos(\beta)\\ Q_{+-}:=-i\sin(\beta)\cos(\beta).\end{cases}$

(27)

We can also write explicit combinatorial formulas for the sums in the exponential:

	$\displaystyle\begin{split}\sum_{\begin{subarray}{c}S\subseteq\{1,\dots,n\}\\ |S|=q\end{subarray}}(z_{S}-z_{S}^{\prime})&=\sum_{k=0}^{q}(-1)^{q-k}\left[\binom{n_{++}+n_{+-}}{k}\binom{n_{-+}+n_{--}}{q-k}-\binom{n_{++}+n_{-+}}{k}\binom{n_{+-}+n_{--}}{q-k}\right]\\ &=\sum_{k=0}^{q}(-1)^{q-k}\left[\binom{\frac{n+(n_{+-}-n_{-+})+(n_{++}-n_{--})}{2}}{k}\binom{\frac{n-(n_{+-}-n_{-+})-(n_{++}-n_{--})}{2}}{q-k}\right.\\ &\qquad\qquad\qquad\qquad-\left.\binom{\frac{n-(n_{+-}-n_{-+})+(n_{++}-n_{--})}{2}}{k}\binom{\frac{n+(n_{+-}-n_{-+})-(n_{++}-n_{--})}{2}}{q-k}\right]\\ &=:f_{q}(n_{++},n_{+-},n_{-+},n_{--})\end{split}$			(28)
	$\displaystyle\begin{split}\sum_{\begin{subarray}{c}S\subseteq\{1,\dots,n\}\\ |S|=q\end{subarray}}(z_{S}-z_{S}^{\prime})^{2}&=4\sum_{k\text{ odd}}\binom{n_{+-}+n_{-+}}{k}\binom{n_{++}+n_{--}}{q-k}\\ &=4\sum_{k\text{ odd}}\binom{n_{+-}+n_{-+}}{k}\binom{n-(n_{+-}+n_{-+})}{q-k}\\ &=:g_{q}(n_{++},n_{+-},n_{-+},n_{--}).\end{split}$			(29)

Therefore, the summand indeed depends only on $(n_{++},n_{+-},n_{-+},n_{--})$ and the number of ways to assign the $n$ positions into four groups of these sizes is the multinomial $n!/(n_{++}!n_{+-}!n_{-+}!n_{--}!)$. To condense our notation we will use $\{n_{*}\}$ to denote the set of sketches $n_{*}=(n_{++},n_{+-},n_{-+},n_{--})$, allowing us to use the shorthand

$\displaystyle\sum_{\{n_{*}\}}\binom{n}{n_{*}}F(n_{*})$

$\displaystyle:=\sum_{\begin{subarray}{c}(n_{++},n_{+-},n_{-+},n_{--})\in\mathbb{N}^{4}\\ n_{++}+n_{+-}+n_{-+}+n_{--}=n\end{subarray}}\binom{n}{n_{++},n_{+-},n_{-+},n_{--}}F(n_{++},n_{+-},n_{-+},n_{--}).$

(30)

Note that this summation has $\binom{n+3}{3}$ terms. We thus have

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle e^{i\lambda H/n}\rangle\right]$

$\displaystyle=\sum_{\{n_{*}\}}\binom{n}{n_{*}}\exp\left({-\sum_{q=1}^{d}\left[\gamma^{2}g_{q}(n_{*})+2\gamma\lambda\frac{f_{q}(n_{*})}{n}+\lambda^{2}\frac{\left(\begin{smallmatrix}n\\ q\end{smallmatrix}\right)}{n^{2}}\right]\frac{\sigma_{q}^{2}}{2n^{q-1}}}\right)\smashoperator[r]{\prod_{ss^{\prime}\in\{\pm\}^{2}}^{}}Q_{ss^{\prime}}^{n_{ss^{\prime}}}.$

(31)

Eq. 31 is the form of the moment-generating function we will use to evaluate $\operatorname*{\mathbb{E}}_{J}[\langle H/n\rangle]$ and $\operatorname*{\mathbb{E}}_{J}[\langle(H/n)^{2}\rangle]$.

Using Eq. 21 combined with the form of the moment-generating function given in Eq. 31, we can write the first moment as:

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]$

$\displaystyle=i\gamma\sum_{q=1}^{d}\frac{\sigma_{q}^{2}}{n^{q}}\sum_{\{n_{*}\}}\binom{n}{n_{*}}\exp\left({-\sum_{q^{\prime}=1}^{d}\gamma^{2}g_{q^{\prime}}(n_{*})\frac{\sigma_{q^{\prime}}^{2}}{2n^{q^{\prime}-1}}}\right)f_{q}(n_{*})\smashoperator[r]{\prod_{ss^{\prime}\in\{\pm\}^{2}}^{}}Q_{ss^{\prime}}^{n_{ss^{\prime}}},$

(32)

and the second moment as:

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{2}\rangle\right]=$

$\displaystyle\sum_{q=1}^{d}\binom{n}{q}\frac{\sigma_{q}^{2}}{n^{{q}+1}}-\gamma^{2}\sum_{q,q^{\prime}=1}^{d}\frac{\sigma_{q}^{2}\sigma_{q^{\prime}}^{2}}{n^{q+q^{\prime}}}\sum_{\{n_{*}\}}\binom{n}{n_{*}}\exp\left({-\sum_{q^{\prime\prime}=1}^{d}\gamma^{2}g_{q^{\prime\prime}}(n_{*})\frac{\sigma_{q^{\prime\prime}}^{2}}{2n^{q^{\prime\prime}-1}}}\right)f_{q}(n_{*})f_{q^{\prime}}(n_{*})\smashoperator[r]{\prod_{ss^{\prime}\in\{\pm\}^{2}}^{}}Q_{ss^{\prime}}^{n_{ss^{\prime}}}.$

(33)

The explicit expression for $f_{q}$ (Eq. 28) shows that $f_{q}$ is a degree-$q$ polynomial in the variables $(n_{+-}-n_{-+})$, $(n_{++}-n_{--})$, and $n$, so we can expand $f_{q}$ as

$$f_{q}(n_{*})=\sum_{a+b+c\leq q}f_{q}^{abc}(n_{+-}-n_{-+})^{a}(n_{++}-n_{--})^{b}n^{c}$$

(34)

where the $f^{abc}_{q}$ are constants independent of $n$. In terms of this expansion, we have

$\displaystyle\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle H/n\rangle\right]$

$\displaystyle=i\gamma\sum_{q=1}^{d}\sigma_{q}^{2}\sum_{a+b+c\leq q}\frac{f_{q}^{abc}}{n^{q-c}}\sum_{\{n_{*}\}}\binom{n}{n_{*}}\exp\left(-\sum_{q^{\prime}=1}^{d}g_{q^{\prime}}(n_{*})\frac{\gamma^{2}\sigma_{q^{\prime}}^{2}}{2n^{q^{\prime}-1}}\right)(n_{+-}-n_{-+})^{a}(n_{++}-n_{--})^{b}\smashoperator[r]{\prod_{ss^{\prime}\in\{\pm\}^{2}}^{}}Q_{ss^{\prime}}^{n_{ss^{\prime}}}$

(35)

and

$\displaystyle\begin{split}\operatorname*{\mathbb{E}}_{J\sim N}\left[\langle(H/n)^{2}\rangle\right]&=\sum_{q=1}^{d}\binom{n}{q}\frac{\sigma_{q}^{2}}{n^{q+1}}-\gamma^{2}\sum_{q,q^{\prime}=1}^{d}\sigma_{q}^{2}\sigma_{q^{\prime}}^{2}\sum_{\begin{subarray}{c}a+b+c\leq q\\ a^{\prime}+b^{\prime}+c^{\prime}\leq q^{\prime}\end{subarray}}\frac{f^{abc}_{q}}{n^{q-c}}\frac{f^{a^{\prime}b^{\prime}c^{\prime}}_{q^{\prime}}}{n^{q^{\prime}-c^{\prime}}}\sum_{\{n_{*}\}}\binom{n}{n_{*}}\exp\left(-\sum_{q^{\prime\prime}=1}^{d}g_{q^{\prime\prime}}(n_{*})\frac{\gamma^{2}\sigma_{q^{\prime\prime}}^{2}}{2n^{q^{\prime\prime}-1}}\right)\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times(n_{+-}-n_{-+})^{a+a^{\prime}}(n_{++}-n_{--})^{b+b^{\prime}}\smashoperator[r]{\prod_{ss^{\prime}\in\{\pm\}^{2}}^{}}Q_{ss^{\prime}}^{n_{ss^{\prime}}}.\end{split}$

(36)

Ref. [2] could explicitly evaluate these terms for the small values of $a$ and $b$ relevant for the two-body SK model, using concise expressions for $f_{2}$ and $g_{2}$. However, to get explicit formulas beyond $q=2$ requires carefully counting powers of $n$ to establish which terms survive in the $n\rightarrow\infty$ limit, and using the general expressions for $f_{q}$ and $g_{q}$ (Eqs. 28 and 29) to derive explicit forms of the leading-order terms. To tame this sum our derivation thus uses techniques that go beyond a simple generalization of [2].