text:
# Renormalization Scheme and Higher Loop Stability in Hadronic 𝜏 Decay within Analytic Perturbation Theory

## I  Introduction

The ratio of hadronic to leptonic widths for the inclusive decay of the $`\tau `$-lepton, $`R_\tau =\mathrm{\Gamma }(\tau ^{}hadrons\nu _\tau )/\mathrm{\Gamma }(\tau ^{}\mathrm{}\overline{\nu }_{\mathrm{}}\nu _\tau )`$, gives important information about the QCD running coupling at relatively small energy scales. The theoretical analysis of the hadronic decay of a heavy lepton was performed in  before the experimental discovery of the $`\tau `$-lepton in 1975. Since then, the properties of the $`\tau `$ have been studied very intensively. Numerous publications are devoted to the QCD description of the inclusive decay of the $`\tau `$-lepton and determination of the QCD running coupling $`\alpha _s`$ at the $`\tau `$ mass scale. A detailed consideration of this subject has been given in . Recently, an updated QCD analysis has been performed by the ALEPH  and OPAL  collaborations, where applications of different theoretical approaches to the $`\tau `$-decay have been analyzed.

At present, the $`R_\tau `$-ratio is known experimentally to high accuracy, $`0.5\%`$. Nevertheless, the value of $`\alpha _s`$ extracted from the data has a rather large error, in which theoretical uncertainties are dominant. For example, the ALEPH Collaboration result is $`\alpha _s(M_\tau =1.777\text{GeV})=0.334\pm 0.007_{\mathrm{expt}}\pm 0.021_{\mathrm{theor}}`$ . It should be emphasized that nonperturbative terms, the values of which are not well known, do not dominate these uncertainties, because their contribution is rather small . The main difficulty is associated with the perturbative description.

The original theoretical expression for the width $`\mathrm{\Gamma }(\tau ^{}hadrons\nu _\tau )`$ involves integration over small values of timelike momentum . The perturbative description with the standard running coupling, which has unphysical singularities, becomes ill-defined in this region and some additional ansatz has to be applied to get a finite result for the hadronic width. To this end, one usually transforms to a contour representation for $`R_\tau `$ , which allows one to give meaning to the initial expression and, in principle, perform calculations in the framework of perturbative QCD. Assuming the validity of this transformation it is possible to present results in the form of a truncated power series with $`\alpha _s(M_\tau )`$ as the expansion parameter . There are also other approaches to evaluating the contour integral. The Le Diberder and Pich prescription  allows one to improve the convergence properties of the approximate series and reduce the renormalization scheme (RS) dependence of theoretical predictions. The possibility of using different approaches in the perturbative description of $`\tau `$-decay leads to an uncertainty in the value of $`\alpha _s(M_\tau )`$ extracted from the same experimental data. Moreover, any perturbative description is based on this contour representation, i.e., on the possibility of converting the initial expression involving integration over timelike momenta into a contour integral in the complex momentum plane. To carry out this transition by using Cauchy’s theorem requires certain analytic properties of the hadronic correlator or of the corresponding Adler function. However, the required analytic properties are not automatically maintained in perturbative QCD resummed by the renormalization group. It is well known that at the one-loop level the so-called ghost pole occurs in the invariant charge. Higher-loop corrections do not solve this problem, but merely add some unphysical branch points. The occurrence of incorrect analytic properties in the conventional perturbative approximation makes it impossible to exploit Cauchy’s theorem in this manner and therefore prevents rewriting the initial expression for $`R_\tau `$ in the form of a contour integral in the complex momentum-plane.

In this paper we will use the analytic approach proposed in  (see also  for details). Being inspired by Käll$`\stackrel{´}{\mathrm{e}}`$n–Lehmann analyticity, which is based on general principles of quantum field theory, this method ensures that the running coupling possesses the correct analytic properties, leads to a self-consistent definition of the effective charge in the timelike region  (which cannot be a symmetrical reflection of the spacelike one ), and provides equality between the initial $`R_\tau `$-expression and the corresponding contour representation . A distinguishing feature of the analytic approach is the existence of a universal infrared limiting value of the analytic running coupling at $`q^2=0`$ which is independent of both the QCD scale parameter $`\mathrm{\Lambda }`$ and the choice of renormalization scheme. This limiting value is defined by the general structure of the Lagrangian and turns out to be stable with respect to higher-loop corrections in contrast to the corresponding quantity in conventional perturbation theory (PT). The higher-loop stability of the analytic perturbation theory (APT) holds also for physical observables .

However, it is not sufficient to study the stability with respect to higher-loop corrections; one must also investigate the stability with respect to choice of renormalization scheme. This is also essential in order to estimate the uncertainty of the results obtained. The theoretical ambiguity which is connected with higher-loop corrections and with RS dependence becomes considerable at low energy scales (see, e.g., ). The APT method, as an invariant analytical version of perturbative QCD , improves the situation and gives very stable results over a wide range of renormalization schemes. This has been demonstrated for the $`e^+e^{}`$ annihilation ratio  and for the Bjorken  and Gross–Llewellyn Smith  deep inelastic scattering sum rules.

The main aim of the paper is a study of the RS dependence which appears in the description of the inclusive $`\tau `$ decay within the APT approach. We will consider the $`R_\tau `$-ratio at the next-to-next-to-leading order (NNLO) and the next-to-leading order (NLO) and compare results obtained with those of standard perturbation theory.

## II  QCD parametrization of $`R_\tau `$

The ratio of hadronic to leptonic $`\tau `$-decay widths can be written as

$$R_\tau =3S_{\mathrm{EW}}(|V_{ud}|^2+|V_{us}|^2)(1+\delta _{\mathrm{QCD}}),$$
(1)
where $`S_{\mathrm{EW}}=1.0194\pm 0.0040`$  is the electroweak factor, $`|V_{ud}|=0.9752\pm 0.0007`$ and $`|V_{us}|=0.2218\pm 0.0016`$  are the CKM matrix elements, and $`\delta _{\mathrm{QCD}}`$ is the QCD correction (see  for details).

We first introduce some definitions: $`\mathrm{Im}\mathrm{\Pi }1+r`$ for the hadronic correlator $`\mathrm{\Pi }(q^2)`$ and $`D1+d`$ for the Adler function $`D(q^2)`$. Then for massless quarks one can write $`\delta _{\mathrm{QCD}}`$ as an integral over timelike momentum $`s`$:

$$\delta _{\mathrm{QCD}}=2_0^{M_\tau ^2}\frac{ds}{M_\tau ^2}\left(1\frac{s}{M_\tau ^2}\right)^2\left(1+2\frac{s}{M_\tau ^2}\right)r(s).$$
(2)
Within the conventional perturbative approximation of $`r(s)`$ this integral is ill-defined due to unphysical singularities of the running coupling lying in the range of integration. The most useful trick to rescue the situation is to appeal to analytic properties of the hadronic correlator $`\mathrm{\Pi }(q^2)`$. This opens up the possibility of exploiting Cauchy’s theorem by rewriting the integral in the form of a contour integral in the complex $`q^2`$-plane with the contour being a circle of radius $`M_\tau ^2`$:

$$\delta _{\mathrm{QCD}}=\frac{1}{2\pi i}_{|z|=M_\tau ^2}\frac{dz}{z}\left(1\frac{z}{M_\tau ^2}\right)^3\left(1+\frac{z}{M_\tau ^2}\right)d(z).$$
(3)
Starting from the contour representation (3) the PT description can be developed in the following two ways (see, e.g., ). One is Braaten’s approach  in which the quantity (3) is represented in the form of truncated power series with the expansion parameter $`\alpha _s(M_\tau ^2)`$. The NNLO representation for $`\delta _{\mathrm{QCD}}`$ is written as follows

$$\delta _{\mathrm{QCD}}^{\mathrm{Br}}=a_\tau +r_1a_\tau ^2+r_2a_\tau ^3,$$
(4)
where $`a_\tau \alpha _s(M_\tau ^2)/\pi `$. The coefficients $`r_1`$ and $`r_2`$ in the $`\overline{\mathrm{MS}}`$ scheme with three active flavors are $`r_1=5.2023`$ and $`r_2=26.366`$ . The running coupling satisfies the renormalization group equation:

$$\mu ^2\frac{da}{d\mu ^2}=\frac{b}{2}a^2(1+c_1a+c_2a^2),$$
(5)
where $`b`$, $`c_1`$ and $`c_2`$ are the $`\beta `$-function coefficients. For three active flavors $`b=9/2`$, $`c_1=16/9`$ and $`c_2^{\overline{\mathrm{MS}}}=3863/864`$.

In the approach of Le Diberder and Pich (LP) , the PT expansion is applied to the $`d`$-function***We use the definition $`q^2<0`$ in the Euclidean region. We have made a few changes in notation from that given in : now $`a=\alpha _s/\pi `$, and consequently $`d_1`$ and $`d_2`$ are what we called $`d_2`$ and $`d_3`$ previously.

$$d(q^2)=a(q^2)+d_1a^2(q^2)+d_2a^3(q^2),$$
(6)
where in the $`\overline{\mathrm{MS}}`$-scheme $`d_1^{\overline{\mathrm{MS}}}=1.6398`$ and $`d_2^{\overline{\mathrm{MS}}}=6.3710`$  for three active quarks. Substituting Eq. (6) into Eq. (3) leads to the following expansion, which is not a power series in $`a`$,

$$\delta _{\mathrm{QCD}}^{\mathrm{LP}}=A^{(1)}(a)+d_1A^{(2)}(a)+d_2A^{(3)}(a)$$
(7)
with

$$A^{(n)}(a)=\frac{1}{2\pi i}_{|z|=M_\tau ^2}\frac{dz}{z}\left(1\frac{z}{M_\tau ^2}\right)^3\left(1+\frac{z}{M_\tau ^2}\right)a^n(z).$$
(8)
As noted above, transition to the contour representation (3) requires certain analytic properties of the correlator, namely, that it must be an analytic function in the complex $`q^2`$-plane with a cut along the positive real axis. The correlator parametrized by the PT running coupling does not have this virtue . Moreover, the conventional renormalization group method determines the running coupling in the spacelike region, whereas the initial expression (2) for $`R_\tau `$ contains an integration over timelike momentum. Thus, we are in need of some method of continuing the running coupling from the spacelike to the timelike region that takes into account the proper analytic properties of the running coupling . Because of this failure of analyticity, Eqs. (2) and (3) are not equivalent in the framework of PT  and if one remains within PT, nothing can be said about the errors introduced by this transition.

The analytic approach may eliminate these problems. To make our analysis more transparent and to demonstrate clearly the differences between the consequences of the PT and APT methods, we restrict our consideration here to massless NNLO. The NNLO analysis can be performed in a more rigorous way without model assumptions that allows us to avoid minor details and exhibit the principal features of the APT approach. Thus, other effects, such as nonperturbative terms, higher-loop corrections, and renormalon contributions lie outside of the purpose of this paper. Note, however, that the NNLO approximation is adequate to describe the actual physical situation because numerically the corresponding terms give the principal contribution to the $`R_\tau `$-ratio.

The function $`d(q^2)`$, which is analytic in the cut $`q^2`$-plane, can be expressed in terms of the effective spectral function $`\rho (\sigma )`$, the basic quantity in the APT method,

$$d(q^2)=\frac{1}{\pi }_0^{\mathrm{}}\frac{d\sigma }{\sigma q^2}\rho (\sigma ).$$
(9)
The connection between the QCD corrections to the $`D`$\- and $`R`$-functions can be written down in the form of the dispersion integral

$$d(q^2)=q^2_0^{\mathrm{}}\frac{ds}{(sq^2)^2}r(s),$$
(10)
which is inverted by the following formula

$$r(s)=\frac{1}{2\pi i}_{siϵ}^{s+iϵ}\frac{dz}{z}d(z).$$
(11)
Here, the contour lies in the region of analyticity of the $`D`$-function. In terms of $`\rho (\sigma )`$ the function $`r(s)`$ defined for timelike momenta can be expressed as follows :

$$r(s)=\frac{1}{\pi }_s^{\mathrm{}}\frac{d\sigma }{\sigma }\rho (\sigma ).$$
(12)
Eqs. (9) and (12) determine the QCD corrections $`d(q^2)`$, which is defined in the Euclidean (spacelike) region of momenta, and $`r(s)`$ defined for the Minkowskian (timelike) argument, in terms of the spectral function $`\rho (\sigma )`$. For $`\delta _{\mathrm{QCD}}`$, using Eq. (2) or equivalently Eq. (3), in terms of $`\rho (\sigma )`$, we find To distinguish APT and PT cases, we will use subscripts “an” and “pt”.

$$\delta _{\mathrm{an}}=\frac{1}{\pi }_0^{\mathrm{}}\frac{d\sigma }{\sigma }\rho (\sigma )\frac{1}{\pi }_0^{M_\tau ^2}\frac{d\sigma }{\sigma }\left(1\frac{\sigma }{M_\tau ^2}\right)^3\left(1+\frac{\sigma }{M_\tau ^2}\right)\rho (\sigma ).$$
(13)
In the APT approach, the spectral function is defined as the imaginary part of the perturbative approximation to $`d_{\mathrm{pt}}(q^2)`$ on the physical cut:

$$\rho (\sigma )=\varrho _0(\sigma )+d_1\varrho _1(\sigma )+d_2\varrho _2(\sigma ),$$
(14)
where

$$\varrho _n(\sigma )=\mathrm{Im}[a^{n+1}(\sigma +iϵ)].$$
(15)
Substituting Eq. (14) into Eq. (13), we can rewrite $`\delta _{\mathrm{an}}`$ in the form of the APT expansion

$$\delta _{\mathrm{an}}=\delta ^{(0)}+d_1\delta ^{(1)}+d_2\delta ^{(2)}.$$
(16)
Note that the APT representations of the $`d`$-function and the QCD correction $`\delta _{\mathrm{QCD}}`$ are not in the form of power series.

The function $`\varrho _0(\sigma )`$ in Eq. (14) defines the analytic spacelike running coupling

$$a_{\mathrm{an}}(q^2)=\frac{1}{\pi }_0^{\mathrm{}}\frac{d\sigma }{\sigma q^2}\varrho _0(\sigma ).$$
(17)
In the one-loop approximation it leads to

$$a_{\mathrm{an}}(q^2)=a_{\mathrm{pt}}(q^2)+\frac{2}{b}\frac{\mathrm{\Lambda }^2}{\mathrm{\Lambda }^2+q^2}.$$
(18)
Unlike the one-loop PT running coupling, $`a_{\mathrm{pt}}(q^2)=2/b\mathrm{ln}(q^2/\mathrm{\Lambda }^2)`$, the analytic running coupling (18) has no unphysical ghost pole and, therefore, possesses the correct analytic properties, arising from Källén-Lehmann analyticity that reflects the general principles of the theory. The nonperturbative (non-logarithmic) term, which appears in the analytic running coupling, does not change the ultraviolet limit of the theory and thus the APT and the PT approaches coincide with each other in the asymptotic region of high energies.

Thus, the APT approach provides a self-consistent description of the hadronic $`\tau `$ decay. This description can be equivalently phrased either on the basis of the original expression (2), which involves the Minkowskian quantity $`r(s)`$, or on the contour representation (3), which involves the Euclidean quantity $`d(q^2)`$.

An important feature of the APT approach is the fact that $`d_{\mathrm{an}}(q^2)`$ and $`a_{\mathrm{an}}(q^2)`$ have a universal limit at the point $`q^2=0`$. This limiting value, generally, is independent of both the scale parameter $`\mathrm{\Lambda }`$ and the order of the loop expansion being considered. Because $`d_{\mathrm{an}}(0)`$ and $`a_{\mathrm{an}}(0)`$ are equal to the reciprocal of the first coefficient of the QCD $`\beta `$-function, they are also RS invariant (we consider only gauge- and mass-independent RSs). The existence of this fixed point plays a decisive role in the improved convergence properties relative to PT and in the very weak RS dependence of our results.

To find the analytic function $`d(q^2)`$ involved in Eq. (3), we solve the transcendental equation for the running coupling

$$\frac{b}{2}\mathrm{ln}\left(\frac{q^2}{\mathrm{\Lambda }_{\overline{\mathrm{MS}}}^2}\right)i\pi \frac{b}{2}=d_1^{\overline{\mathrm{MS}}}d_1+\frac{1}{a}+c_1\mathrm{ln}\left(\frac{b}{2c_1}\right)+F^{(l)}(a),$$
(19)
where at the NLO

$$F^{(2)}(a)=c_1\mathrm{ln}\left(\frac{c_1a}{1+c_1a}\right),$$
(20)
and at the NNLO

$$F^{(3)}(a)=F^{(2)}(a)+c_2_0^a\frac{dx}{(1+c_1x)(1+c_1x+c_2x^2)},$$
(21)
on the physical cut lying along the positive real axis in the complex $`q^2`$-plane and then use Eqs. (14), (15) and (9). Eq. (19) holds in any MS-like renormalization scheme and allows us to normalize the results obtained by using the scale parameter $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$. Having found $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$, we can study how $`\delta _{\mathrm{an}}`$ varies with a change of renormalization scheme. To do that one has to select parameters which determine the RS. The function $`d(q^2)`$ in Eq. (3) is parametrized by a set of RS-dependent parameters. There are RS invariant combinations which constrain these parameters . At the NNLO there are two RS-invariant quantities; the first of them expresses an energy dependence, the second is just a number

$$\omega _2=c_2+d_2c_1d_1d_{1}^{}{}_{}{}^{2},$$
(22)
which in our case equals 5.2378. Here, $`c_1`$ is RS invariant and we can choose $`d_1`$ and $`c_2`$ as independent variables, which define some RS.

There are no fundamental principles upon which one can choose one or another preferable RS. Nevertheless, a natural way of studying the RS dependence is to supplement results in a certain scheme with an estimate of the variability of the predictions over a range of a priori acceptable schemes specified by some criterion. In  it was proposed to consider the class of ‘natural’ RSs, which obey the condition

$$|c_2|+|d_2|+c_1|d_1|+d_{1}^{}{}_{}{}^{2}C|\omega _2|.$$
(23)
This inequality is called the “cancellation index criterion” which means that the degree of cancellation in the second RS invariant (22) should not be too large. To define a boundary of ‘acceptable’ schemes which is defined by the value of the cancellation index $`C`$, we will require no more cancellation than that which occurs in the scheme obeying the principle of minimal sensitivity (PMS) , which leads to $`C2`$.

## III APT: convergence properties and RS dependence

For various physical quantities, the APT approach allows one to construct a series that has improved convergence properties as compared to a perturbative expansion. To demonstrate this fact for the hadronic $`\tau `$-decay, we compare the convergence properties of the PT expansions (4) and (7) on the one hand, and the APT approach given by Eq. (16) on the other hand. For our calculation we take as input the TAU’98 conference value: $`R_\tau =3.642\pm 0.019`$ , which is consistent with the PDG’98 fit $`R_\tau =3.642\pm 0.024`$ . In Table I we present NNLO results obtained by the methods mentioned above for the central experimental value in the $`\overline{\mathrm{MS}}`$ scheme. The relative contributions of higher-order terms depends on the method which is applied. The convergence properties of the APT expansion seem to be much improved compared to those of the PT expansions.

The values of the scale parameter $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ and the coupling $`\alpha _s(M_\tau ^2)`$ obtained from above PT and APT expansions are noticeably different from each other. The corresponding numerical estimations are given in Table II, in which, in order to clarify the situation concerning higher-loop stability of different expansions, we also present the NLO result. This table demonstrates that the theoretical ambiguity, which associated with different versions of the perturbative description, leads to a rather large uncertainty, $`\alpha _{\mathrm{PT}(\mathrm{Br})}^{\mathrm{NNLO}}\alpha _{\mathrm{PT}(\mathrm{LP})}^{\mathrm{NNLO}}=0.012`$. At the same time the experimental error is $`\mathrm{\Delta }\alpha _{\mathrm{expt}}=0.007`$$`0.009`$ . The distinction between NLO and NLLO running coupling values is $`12\%`$ for PT (Br) and $`5\%`$ for PT (LP) approaches, while for the APT approach it is less than $`0.5\%`$.

The non-logarithmic terms, which ensure the correct analytic properties and allow a self-consistent description of $`\tau `$ decay, turn out to be very important for the numerical analysis and influence in an essential way the value of $`\mathrm{\Lambda }`$ parameter extracted from the data. Indeed, at the one-loop level one can write a simple relation: $`\delta _{\mathrm{an}}(\mathrm{\Lambda })\delta _{\mathrm{pt}}^{\mathrm{LP}}(\mathrm{\Lambda })(2/b)\mathrm{\Lambda }^2/M_\tau ^2`$. The second term, which is ‘invisible’ in the perturbative expansion, turns out to be numerically important  (see the detailed discussion in ). Note also that there is a difference between the shapes of the analytic and perturbative running couplings, for example, $`\alpha _{\mathrm{an}}(\mathrm{\Lambda }=907\text{MeV})=0.403`$, while at the same scale, the value of the perturbative coupling much larger, $`\alpha _{\mathrm{pt}}(\mathrm{\Lambda }=907\text{MeV})=0.796`$. Here the question may arise, how is the large APT value of $`\mathrm{\Lambda }`$ consistent with high energy experimental data? We have estimated the ratio of hadronic to leptonic Z-decay widths, $`R_\mathrm{Z}`$, using the above value of $`\mathrm{\Lambda }_{\mathrm{an}}`$ and the matching procedure proposed in . We obtained the value $`R_\mathrm{Z}=20.82`$, which lies within the range of experimental errors; for example, the PDG’98 average is $`\mathrm{R}_\mathrm{Z}=20.77\pm 0.07`$ . This fact can be understood if one takes into account that there are differences between the shapes of the analytic and perturbative running couplings and also in the terms of the corresponding series.

We found that the value of $`\delta _{\mathrm{an}}`$ depends so slightly on $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ that a $`0.9\%`$ error in $`R_\tau `$ gives $`18\%`$ error in the value of $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$. (This is the reason why the errors in the values of $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ and $`\alpha (M_\tau ^2)`$ given by APT are larger than those in PT.) We illustrate this feature in Table III. According to the table, when we change $`M_\tau ^2/\mathrm{\Lambda }^2`$ from 2.0 to 6.5 (corresponding to a variation of $`\mathrm{\Lambda }`$ from 1.256 GeV to 0.697 GeV), $`\delta _{\mathrm{an}}`$ is only altered by about $`20\%`$. The sensitivity to $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ increases as $`M_\tau ^2/\mathrm{\Lambda }^2`$ gets smaller.

Consider now the RS dependence of the APT result and compare it with the perturbative LP approach,The LP approach is often called the contour-improved fixed-order PT (CIPT) . which of the two PT schemes is more preferable from the viewpoint of sensitivity to RS dependence. In the framework of the PT, the RS dependence of $`\delta _{\mathrm{QCD}}`$ has been discussed in detail in .

In the $`\overline{\mathrm{MS}}`$ scheme we adopt $`\delta _{\mathrm{an},\mathrm{pt}}=0.1906`$ and consider some RS belonging to the domain described above \[see Eq. (23)\]. Take two schemes, $`A`$ and $`B`$, located at the two lower corners of the boundary of the domain (see Fig. 1), i.e., they have the same cancellation index as does the PMS scheme, with $`A=(1.6183,0)`$ and $`B=(0.9575,0)`$, where the first coordinate is $`d_1`$ and the second is $`c_2`$. Then for the PT case in NNLO we get $`\delta _{\mathrm{pt}}(A)=0.2025`$ and $`\delta _{\mathrm{pt}}(B)=0.1911`$. Therefore, even for this sufficiently narrow class of RS the perturbative approach gave a $`6\%`$ deviation in $`\delta _{\mathrm{QCD}}`$ that corresponds to a RS uncertainty for the running coupling value in the $`\overline{\mathrm{MS}}`$-scheme of $`\mathrm{\Delta }\alpha _{\mathrm{pt}}^{\mathrm{RS}}=0.0153`$. The difference between APT results is much smaller: $`\delta _{\mathrm{an}}(A)=0.1890`$ and $`\delta _{\mathrm{an}}(B)=0.1905`$, and we have only $`0.8\%`$ deviation, which corresponds in the $`\overline{\mathrm{MS}}`$-scheme to a RS uncertainty of $`\mathrm{\Delta }\alpha _{\mathrm{an}}^{\mathrm{RS}}=0.0035`$. The similar RS stability holds also at the two-loop level: one has a $`5\%`$ deviation in the PT case and only a $`0.4\%`$ for the APT. We display our NNLO results in the form of a contour plot, in Fig. 1.

It is worthwhile to analyze some schemes lying outside the domain considered above with the relatively small value of the cancellation index $`C2`$. Among them there is, for instance, the commonly used $`\overline{\mathrm{MS}}`$ scheme which does not belong this domain. In  it was shown that the so-called $`V`$ scheme  lies very far from the domain described above and gives so a large value of $`\delta _{\mathrm{pt}}`$ that it cannot be used at this low energy. For the $`V`$ scheme we have $`d_1=0.109`$ and $`c_2=26.200`$. The three-loop perturbative result is $`\delta _{\mathrm{pt}}(V)=0.3060`$ that corresponds to about a $`61\%`$ deviation from the $`\overline{\mathrm{MS}}`$ scheme. On the other hand, if we turn to APT we have $`\delta _{\mathrm{an}}=0.1902`$, i.e., only about a $`0.2\%`$ deviation from the $`\overline{\mathrm{MS}}`$ scheme. So the $`V`$ scheme is still useful at this energy in APT.

In PT at high energies the weak RS dependence is a consequence of the small value of the coupling constant. At lower energies the uncertainty increases. In APT, at high energies, the situation is the same. However, at low energies the theory has a universal RS-invariant infrared limiting value $`d_{\mathrm{an}}(0)`$, which restricts the RS ambiguity over a very wide range of momentum. Another way to illustrate the remarkable stability of APT is to calculate the spectral functions $`\varrho _n(\sigma )`$ given by Eq. (15); one sees that $`\varrho _1(\sigma )`$ is much smaller than $`\varrho _0(\sigma )`$ over the whole spectral region. The same statement is true for the relationship between $`\varrho _1(\sigma )`$ and $`\varrho _2(\sigma )`$. This monotonically decreasing behavior reduces the RS dependence strongly, since the perturbative coefficients $`d_1`$ and $`d_2`$ in expression (14) for $`\rho (\sigma )`$ are multiplied by these functions. For the $`\overline{\mathrm{MS}}`$ scheme, this situation is demonstrated in Fig. 2.

## IV Conclusion

We have considered inclusive $`\tau `$-decay in three-loop order within analytic perturbation theory concentrating on the analysis of theoretical uncertainties coming from the perturbative short distance part of the QCD correction to the $`R_\tau `$-ratio, which defines the principal contribution to this physical quantity. For the low energy $`\tau `$-mass scale, the main source of theoretical uncertainties results from the inevitable truncation of the perturbative series, which leads to the essential RS dependence and higher loop sensitivity of the theoretical predictions. In order to resolve this problem within the conventional perturbative approach it is possible to try, in principle, to compute higher loop contributions. However, even if this were to be done, one has to keep in mind that from the rigorous point of view it will hardly be sufficient because the series is asymptotic, and, in any finite order, the analytic properties of the hadronic correlator, which arise from general principles of the theory, are violated. Thus, to resolve this problem one has to use a modification of the perturbative expansion at low energy scales.

Here, we have applied the analytic approach which is not inconsistent with the general principles of quantum field theory and which opens up the possibility of reducing the theoretical uncertainties associated with short distance contributions mentioned above. Let us summarize the important features of this method: (i) the method maintains the correct analytic properties and leads to a self-consistent definition of the procedure of analytic continuation from the spacelike to the timelike region; (ii) the APT approach has much improved convergence properties and turns out to be stable with respect to higher-loop corrections; (iii) the RS dependence of the results obtained is reduced drastically. For example, the $`V`$ scheme, which gives a very large discrepancy in standard perturbation theory, can be used in analytic perturbation theory without any difficulty and the APT predictions are practically RS independent over a wide region of RS parameters.

The nonperturbative power corrections coming from the operator product expansion (in this connection see a discussion in ), renormalon and other effects are beyond the scope of our present consideration. Note, however, that the process of enforcing analyticity modifies the perturbative contributions by incorporating some nonperturbative terms. The form of the APT running coupling and the non-power structure of the APT expansion are essentially different from the PT ones. Numerically, this difference becomes very important in the region less than a few GeV and in order to get the same physical quantity the contribution of power corrections should also be changed.

The value of $`\mathrm{\Lambda }_{\mathrm{APT}}`$ is very sensitive to the experimental value of $`R_\tau `$. For example, as has been demonstrated in  the use of the value of $`R_\tau `$ obtained by the CLEO collaboration  gives a value of the scale parameter some 30% smaller than that found here. Note also that the renormalon contribution influences the value of $`\mathrm{\Lambda }`$ extracted from the $`\tau `$ data (see  for a review). Within the usual approach, renormalons reduce the value of $`\alpha _s(M_\tau ^2)`$ by about $`15\%`$. A similar situation holds also in APT and for the nonperturbative $`a`$-expansion approach , which allows one, as in APT, to maintain the required analyticity . These two analytic approaches often lead to rather similar consequences. For example, they allow one to get a good description of experimental data corresponding to the Euclidean and Minkowskian characteristics of the process of $`e^+e^{}`$ annihilation into hadrons down to the lowest energy scale .

Pure massless APT analysis, which has been performed here, leads to an unusually large value of the QCD scale parameter $`\mathrm{\Lambda }`$ as compared to the conventional PT value. This is connected with the presence of nonperturbative contributions that appear in the APT method which have a negative relative sign. The effects mentioned above can change the value of the scale parameter extracted from the $`\tau `$ data. However, this fact is not relevant for the essential conclusion which we have claimed in this paper, that the APT method provides predictions which are stable with respect to the choice of renormalization scheme and to the inclusion of higher loop corrections. Thus, the analytic approach discussed here is not in conflict with the general principles of the theory and allows one to reduce the uncertainties of theoretical predictions drastically.

## Acknowledgements

The authors would like to thank D.V. Shirkov for interest in this work. Partial support of the work by the US National Science Foundation, grant PHY-9600421, and by the US Department of Energy, grant DE-FG-03-98ER41066, and by the RFBR, grants 99-01-00091, 99-02-17727 is gratefully acknowledged. The work of ILS and OPS is also supported in part by the University of Oklahoma, through its College of Arts and Science, the Vice President for Research, and the Department of Physics and Astronomy.


format: text

text:
# The Statistical Physics of Regular Low-Density Parity-Check Error-Correcting Codes

## I Introduction

Error-correcting codes are commonly used for reliable data transmission through noisy media, especially in the case of memoryless communication where corrupted messages cannot be repeatedly sent. These techniques play an important role in a wide range of applications from memory devices to deep space explorations, and are expected to become even more important due to the rapid development in mobile phones and satellite-based communication.

In a general scenario, the sender encodes an $`N`$ dimensional Boolean message vector $`𝝃`$, where $`\xi _i(0,1),i`$, to an $`M(>N)`$ dimensional Boolean codeword $`𝒛_0`$, which is then being transmitted through a noisy communication channel. Noise corruption during transmission can be modelled by the noise vector $`𝜻`$, where corrupted bits are marked by the value 1 and all other bits are zero, such that the received corrupted codeword takes the form $`𝒛=𝒛_0+𝜻\text{(mod 2)}`$. The received corrupted message is then decoded by the receiver for retrieving the original message $`𝝃`$.

The error-correcting ability comes at the expense of information redundancy. Shannon showed in his seminal work that error-free communication is theoretically possible if the code rate, representing the fraction of informative bits in the transmitted codeword, is below the channel capacity; in the case of unbiased messages transmitted through a Binary Symmetric Channel (BSC), which we will focus on here, $`R=N/M`$ satisfies

$$R<1+p\mathrm{log}_2p+(1p)\mathrm{log}_2(1p).$$
(1)
The expression on the right is termed Shannon’s bound. However, Shannon’s derivation is non-constructive and the quest for codes which saturate Eq.(1) has been one of the central topics of information theory ever since.

In this paper we examine the efficiency and limitations of Gallager-type error-correcting code, which attracted much interest recently among researchers in this field. This code was discovered almost forty years ago by Gallager but was abandoned shortly after its invention due to the computational limitations of the time. Since their recent rediscovery by MacKay and Neal, different variations of Gallager-type codes have been developed  attempting to get as close as possible to saturating Shannon’s bound.

In this paper we will examine the typical properties of a family of codes based on one variation, the MN code , using the established methods of statistical physics, to provide a theoretical study based on the typical performance of codes rather on the worst case analysis.

This paper is organised as follows: In the next two sections, we introduce Gallager-type error-correcting codes in detail and link them to the statistical mechanics framework. We then examine the equilibrium properties of various members of this family of codes using the replica method (section IV) and compare the bit-error rate below criticality. In section V, we examine the relation between Belief-Propagation (BP) decoding and the Thouless-Anderson-Palmer (TAP) approach to diluted spin systems; we then use it for comparing the analytical results obtained via the replica method to those obtained from simulations using BP decoding. In section VI we show a computationally efficient construction for one of the more practical constructions. Finally, we present conclusions for the current work and suggest future research directions.

## II Gallager-type error-correcting codes

There are several variations in Gallager-type error-correcting codes. The one discussed in this paper is termed the MN code, recently introduced by MacKay and Neal. In these codes, a Boolean message $`𝝃`$ is encoded into a codeword $`𝒛_0`$ using two randomly constructed Boolean sparse matrices $`C_s`$ and $`C_n`$, which are characterised in the following manner.

The rectangular sparse matrix $`C_s`$ is of dimensionality $`M\times N`$, having randomly chosen $`K`$ non-zero unit elements per row and $`C`$ per column. The matrix $`C_n`$ is an $`M\times M`$ (mod 2)-invertible matrix having randomly chosen $`L`$ non-zero elements per row and column. These matrices are shared by the sender and the receiver.

Using these matrices, one can encode a message $`𝝃`$ into a codeword $`𝒛_0`$ in the following manner

$$𝒛_0=C_n^1C_s𝝃\text{(mod 2)},$$
(2)
which is then transmitted via a noisy channel. Note that all matrix and vector components are Boolean $`(0,1)`$, and all summations are carried out in this field. For simplicity, the noise process is modelled hereafter by a binary symmetric channel (BSC), where each bit is independently flipped with probability $`p`$. Extending the code presented here to other types of noise is straightforward.

During transmission, a noise vector $`𝜻`$ is added to $`𝒛_0`$ and a corrupted codeword $`𝒛=𝒛_0+𝜻`$ (mod 2) is received at the other end of the channel. Decoding is then carried out by taking the product of the matrix $`C_n`$ and the received codeword $`𝒛`$, which results in $`C_s𝝃+C_n𝜻=C_n𝒛𝑱`$. The equation

$$C_s𝑺+C_n𝝉=𝑱\text{(mod 2)},$$
(3)
is solved via the iterative methods of belief propagation (BP) to obtain the most probable Boolean vectors $`𝑺`$ and $`𝝉`$. BP methods in this context have recently been shown to be identical to a TAP based solution of a similar physical system.

## III A Statistical Mechanics Perspective

Sourlas was the first to point out that error-correcting codes of this type have a similarity to Ising spin systems in statistical physics; he demonstrated this using a simple version of the same nature. His work, that focused on extensively connected systems, was recently extended to finitely connected systems. We follow a similar approach in the current investigation; preliminary results have already been presented in.

To facilitate the statistical physics analysis, we first employ the binary representation $`(\pm 1)`$ of the dynamical variables $`𝑺`$ and $`𝝉`$ and of the check vector $`𝑱`$ rather than the Boolean one $`(0,1)`$. The $`\mu `$-th component of Eq.(3) is then rewritten as

$$\underset{i_s(\mu )}{}S_i\underset{j_n(\mu )}{}\tau _j=J_\mu ,$$
(4)
where $`_s(\mu )`$ and $`_n(\mu )`$ are the sets of all indices of non-zero elements in row $`\mu `$ of the sparse matrices $`C_s`$ and $`C_n`$, respectively. The check $`\mu `$ is given by message $`𝝃`$ and noise $`𝜻`$ as $`J_\mu =_{i_s(\mu )}\xi _i_{j_n(\mu )}\zeta _j`$; it should be emphasised that the message vector $`𝝃`$ and the noise vector $`𝜻`$ themselves are not known to the receiver.

An interesting link can now be formulated between the Bayesian framework of MN codes and Ising spin systems. Rewriting Kronecker’s delta for binary variables $`x`$ and $`y`$ as $`\delta [x;y]=(1+xy)/2=lim_\beta \mathrm{}\mathrm{exp}(\beta \delta [1;xy])`$, one may argue that, using it as a likelihood, equation (4) gives rise to the conditional probability of the check $`𝑱`$ for given $`𝑺`$, $`𝝉`$, $`C_s`$ and $`C_n`$

$`𝒫(𝑱|𝑺,𝝉,C_s,C_n)=\underset{\beta \mathrm{}}{lim}\mathrm{exp}\left(\beta {\displaystyle \underset{\mu =1}{\overset{M}{}}}\delta [1;J_\mu {\displaystyle \underset{i_s(\mu )}{}}S_i{\displaystyle \underset{j_n(\mu )}{}}\tau _j]\right).`$ (5)
Prior knowledge about possibly biased message and noise is represented by the prior distributions

$`𝒫_s(𝑺)={\displaystyle \frac{\mathrm{exp}\left(F_s\underset{i=1}{\overset{N}{}}S_i\right)}{\left(2\mathrm{cosh}F_s\right)^N}},𝒫_n(𝝉)={\displaystyle \frac{\mathrm{exp}\left(F_n\underset{j=1}{\overset{M}{}}\tau _j\right)}{\left(2\mathrm{cosh}F_n\right)^M}},`$ (6)
respectively. Non-zero field $`F_s`$ is introduced for biased message and $`F_n`$ is determined by flip rate $`p`$ of channel noise as $`F_n=(1/2)\mathrm{ln}\left((1p)/p\right)`$. Using equations (5) and (6), the posterior distribution of $`𝑺`$ and $`𝝉`$ for given check $`𝑱`$ and matrices $`C_s`$ and $`C_n`$ is of the form

$`𝒫(𝑺,𝝉|𝑱,C_s,C_n)`$ $`=`$ $`{\displaystyle \frac{𝒫(𝑱|𝑺,𝝉,C_s,C_n)𝒫_s(𝑺)𝒫_n(𝝉)}{𝒫\left(𝑱|C_s,C_n\right)}}`$ (7)
$`=`$ $`\underset{\beta \mathrm{}}{lim}{\displaystyle \frac{\mathrm{exp}\left(\beta (𝑺,𝝉|𝓙,𝒟)\right)}{𝒵(𝓙,𝒟)}},`$ (8)
where $`𝒫(𝑱|C_s,C_n)=_{\{𝑺,𝝉\}}𝒫(𝑱|𝑺,𝝉,C_s,C_n)𝒫_s(𝑺)𝒫_n(𝝉)`$,

$`(𝑺,𝝉|𝓙,𝒟)`$ $`=`$ $`{\displaystyle \underset{\mu =1}{\overset{M}{}}}\delta [1;J_\mu {\displaystyle \underset{i_s(\mu )}{}}S_i{\displaystyle \underset{j_n(\mu )}{}}\tau _j]{\displaystyle \frac{F_s}{\beta }}{\displaystyle \underset{i=1}{\overset{N}{}}}S_i{\displaystyle \frac{F_n}{\beta }}{\displaystyle \underset{j=1}{\overset{M}{}}}\tau _j`$ (9)
$`=`$ $`{\displaystyle \underset{<i_1,..,i_K;j_1,..,j_L>}{}}𝒟_{<i_1,..,i_K;j_1,..,j_L>}\delta [1;𝒥_{<i_1,..,i_K;j_1,..,j_L>}S_{i_1}\mathrm{}S_{i_K}\tau _{j_1}\mathrm{}\tau _{j_L}]`$ (10)
$`{\displaystyle \frac{F_s}{\beta }}{\displaystyle \underset{i=1}{\overset{N}{}}}S_i{\displaystyle \frac{F_n}{\beta }}{\displaystyle \underset{j=1}{\overset{M}{}}}\tau _j,`$ (11)
and

$`𝒵(𝓙,𝒟)`$ $`=`$ $`\underset{\beta \mathrm{}}{lim}{\displaystyle \underset{\{𝑺,𝝉\}}{}}\mathrm{exp}\left(\beta (𝑺,𝝉|𝓙,𝒟)\right)`$ (12)
$`=`$ $`{\displaystyle \underset{\{𝑺,𝝉\}}{}}{\displaystyle \underset{<i_1,..,i_K;j_1,..,j_L>}{}}\left[1+{\displaystyle \frac{1}{2}}𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}\left(𝒥_{<i_1,..,i_K;j_1,..,j_L>}S_{i_1}\mathrm{}S_{i_K}\tau _{j_1}\mathrm{}\tau _{j_L}1\right)\right]`$ (13)
$`\times `$ $`\mathrm{exp}\left(F_s{\displaystyle \underset{i=1}{\overset{N}{}}}S_i+F_n{\displaystyle \underset{j=1}{\overset{M}{}}}\tau _j\right).`$ (14)
The final form of posterior distribution (8) implies that the MN code is identical to an Ising spin system defined by the Hamiltonian (11) in the zero temperature limit $`T=\beta ^10`$. In equations (11) and (14), we introduced the sparse connectivity tensor $`𝒟_{<i_1,..,j_L>}`$ which takes the value 1 if the corresponding indices of both message and noise are chosen (i.e., if all corresponding indices of the matrices $`C_s`$ and $`C_n`$ are 1) and 0 otherwise, and coupling $`𝒥_{<i_1,..,i_K;j_1,..,j_L>}=\xi _{i_1}\xi _{i_2}\mathrm{}\xi _{i_K}\zeta _{j_1}\zeta _{j_2}\mathrm{}\zeta _{j_L}`$. These come to isolate the disorder in choosing the matrix connections, embedded in $`𝒟_{<i_1,..,j_L>}`$, and to simplify the notation.

The posterior distribution (8) can be used for decoding. One can show that expectation of the overlap between original message $`𝝃`$ and retrieved one $`\widehat{𝝃}`$

$$m=\frac{1}{N}\underset{i=1}{\overset{N}{}}\xi _i\widehat{\xi }_i,$$
(15)
is maximised by setting $`\widehat{𝝃}`$ to its Bayes-optimal estimator

$$\widehat{\xi }_i^B=\text{sign}(m_i^S),m_i^S=\underset{\{𝑺,𝝉\}}{}S_i𝒫(𝑺,𝝉|𝑱,C_s,C_n).$$
(16)
It is worth while noting that this optimal decoding is realized at zero temperature rather than at a finite temperature as in . The reason is that the true likelihood term (5) corresponds to the ground state of the first term of the Hamiltonian (11) due to the existence of more degrees of freedom, in the form of the dynamical variables $`𝝉`$, which do not appear in other systems. Introducing the additional variables $`𝝉`$, the degrees of freedom in the spin system increase from $`N`$ to $`N+M`$, while the number of constraints from the checks $`𝑱`$ remains $`M`$. This implies that in spite of the existence of quenched disorder caused by $`𝓙`$ and $`𝒟`$, the system is free from frustration even in the low temperature limit, which is useful for practical decoding using local search algorithms. The last two terms in Eq.(11) scale with $`\beta `$ remain finite even in the zero temperature limit $`\beta \mathrm{}`$ representing the true prior distributions, which dominates the statistical properties of the system, while the first term vanishes to satisfy the parity check condition (4).

## IV Equilibrium Properties: The Replica Method

As we use the methods of statistical mechanics, we concentrate on the case of long messages, in the of limit of $`N,M\mathrm{}`$ while keeping code rate $`R=N/M=K/C`$ finite. This limit is quite reasonable for this particular problem since Gallager-type codes are usually used in the transmission of long ($`10^410^5`$) messages, where finite size corrections are likely to be negligible.

Since the first part of the Hamiltonian (11) is invariant under the gauge transformation $`S_i\xi _iS_i`$, $`\tau _j\zeta _j\tau _j`$ and $`𝒥_{i_1,\mathrm{},j_L}1`$, it is useful to decouple the correlation between the vectors $`𝑺`$, $`𝝉`$ and $`𝝃`$, $`𝜻`$. Rewriting the Hamiltonian using this gauge, one obtains a similar expression to Eq.(11) apart from the second terms which become $`F_s/\beta _{i=1}\xi _iS_i`$ and $`F_n/\beta _{j=1}\zeta _j\tau _j`$.

Due to the existence of several types of quenched disorder in the system, it is natural to resort to replica method for investigating the typical properties in equilibrium. More specifically, we calculate expectation values of $`n`$-th power of partition function (14) with respect to the quenched variables $`𝝃`$, $`𝜻`$ and $`𝒟`$ and take the limit $`n0`$.

Carrying out the calculation in the zero temperature limit $`\beta \mathrm{}`$ gives rise to a set of order parameters

$`q_{\alpha ,\beta ,..,\gamma }={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_iS_i^\alpha S_i^\beta ,..,S_i^\gamma _\beta \mathrm{},r_{\alpha ,\beta ,..,\gamma }={\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}Y_j\tau _j^\alpha \tau _j^\beta ,..,\tau _j^\gamma _\beta \mathrm{}`$ (17)
where $`\alpha `$, $`\beta ,..`$ represent replica indices, and the variables $`Z_i`$ and $`Y_j`$ come from enforcing the restriction of $`C`$ and $`L`$ connections per index, respectively:

$$\delta \left(\underset{i_2,..,i_K}{}𝒟_{<i,i_2,..,j_L>}C\right)=_0^{2\pi }\frac{dZ}{2\pi }Z^{_{i_2,..,i_K}𝒟_{<i,i_2,..,j_L>}(C+1)},$$
(18)
and similarly for the restriction on the $`j`$ indices.

To proceed further, it is necessary to make an assumption about the order parameters symmetry. The assumption made here is that of replica symmetry in both the order parameters and the related conjugate variables

$`q_{\alpha ,\beta ..\gamma }`$ $`=`$ $`a_q{\displaystyle 𝑑x\pi (x)x^l},\widehat{q}_{\alpha ,\beta ..\gamma }=a_{\widehat{q}}{\displaystyle 𝑑\widehat{x}\widehat{\pi }(\widehat{x})\widehat{x}^l}`$ (19)
$`r_{\alpha ,\beta ..\gamma }`$ $`=`$ $`a_r{\displaystyle 𝑑y\rho (y)y^l},\widehat{r}_{\alpha ,\beta ..\gamma }=a_{\widehat{r}}{\displaystyle 𝑑\widehat{y}\widehat{\rho }(\widehat{y})\widehat{y}^l},`$ (20)
where $`l`$ is the number of replica indices, $`a_{}`$ are normalisation coefficients, and $`\pi (x),\widehat{\pi }(\widehat{x}),\rho (y)`$ and $`\widehat{\rho }(\widehat{y})`$ represent probability distributions. Unspecified integrals are over the range $`[1,+1]`$. This ansatz is supported by the facts that (i) the current system is free of frustration and (ii) there has never been observed replica symmetry breaking at Nishimori’s condition which corresponds to using correct priors $`F_s`$ and $`F_n`$ in our case. The results obtained hereafter also support this ansatz. Extremizing the partition function with respect to distributions $`\pi ()`$, $`\widehat{\pi }()`$, $`\rho ()`$ and $`\widehat{\rho }()`$, one then obtains the free energy per spin

$`f={\displaystyle \frac{1}{N}}\mathrm{ln}𝒵_{\xi ,\zeta ,𝒟}`$ (21)
$`=`$ $`\text{extr}_{\{\pi ,\widehat{\pi },\rho ,\rho \}}\{{\displaystyle \frac{C}{K}}\mathrm{ln}2+C{\displaystyle }dxd\widehat{x}\pi (x)\widehat{\pi }(\widehat{x})\mathrm{ln}(1+x\widehat{x})+{\displaystyle \frac{CL}{K}}{\displaystyle }dyd\widehat{y}\rho (y)\widehat{\rho }(\widehat{y})\mathrm{ln}(1+y\widehat{y})`$ (22)
$`{\displaystyle \frac{C}{K}}{\displaystyle \left[\underset{k=1}{\overset{K}{}}dx_k\pi (x_k)\right]\left[\underset{l=1}{\overset{L}{}}dy_l\rho (y_l)\right]\mathrm{ln}\left[1+\underset{k=1}{\overset{K}{}}x_k\underset{l=1}{\overset{L}{}}y_l\right]}`$ (23)
$`{\displaystyle \left[\underset{k=1}{\overset{C}{}}d\widehat{x}_k\widehat{\pi }(\widehat{x}_k)\right]\mathrm{ln}\left[e^{F_s\xi }\underset{k=1}{\overset{C}{}}(1+\widehat{x}_k)+e^{F_s\xi }\underset{k=1}{\overset{C}{}}(1\widehat{x}_k)\right]_\xi }`$ (24)
$`{\displaystyle \frac{C}{K}}{\displaystyle }\left[{\displaystyle \underset{l=1}{\overset{C}{}}}d\widehat{y}_l\widehat{\rho }(\widehat{y}_l)\right]\mathrm{ln}[e^{F_n\zeta }{\displaystyle \underset{l=1}{\overset{L}{}}}(1+\widehat{y}_l)+e^{F_n\zeta }{\displaystyle \underset{l=1}{\overset{L}{}}}(1\widehat{y}_l)]_\zeta \},`$ (25)
where angled brackets with subscript $`𝝃`$, $`𝜻`$ and $`𝒟`$ denote averages over the message and noise distributions respectively, and sparse connectivity tensor $`𝒟`$. Message averages take the form

$$\mathrm{}_\xi =\underset{\xi =\pm 1}{}\frac{1+\xi \mathrm{tanh}F_s}{2}\left(\mathrm{}\right)$$
(26)
and similarly for $`\mathrm{}_\zeta `$. Details of the derivation are given in Appendix A.

Taking the functional variation of $`f`$ with respect to the distributions $`\pi `$, $`\widehat{\pi }`$, $`\rho `$ and $`\widehat{\rho }`$, one obtains the following saddle point equations

$`\pi (x)`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{C1}{}}d\widehat{x}_l\widehat{\pi }(\widehat{x}_l)\delta \left(x\mathrm{tanh}\left(\xi F_s+\underset{l=1}{\overset{C1}{}}\mathrm{tanh}^1\widehat{x}_l\right)\right)_\xi },`$ (27)
$`\widehat{\pi }(\widehat{x})`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{K1}{}}dx_l\pi (x_l)\underset{l=1}{\overset{L}{}}dy_l\rho (y_l)\delta \left(\widehat{x}\underset{l=1}{\overset{K1}{}}x_l\underset{l=1}{\overset{L}{}}y_l\right)},`$ (28)
$`\rho (y)`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{L1}{}}d\widehat{y}_l\widehat{\rho }(\widehat{y}_l)\delta \left(y\mathrm{tanh}\left(\zeta F_n+\underset{l=1}{\overset{L1}{}}\mathrm{tanh}^1\widehat{y}_l\right)\right)_\zeta },`$ (29)
$`\widehat{\rho }(\widehat{y})`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{K}{}}dx_l\pi (x_l)\underset{l=1}{\overset{L1}{}}dy_l\rho (y_l)\delta \left(\widehat{y}\underset{l=1}{\overset{K}{}}x_l\underset{l=1}{\overset{L1}{}}y_l\right)}.`$ (30)
After solving these equations, the expectation of the overlap between the message $`𝝃`$ and the Bayesian optimal estimator (16), which serves as a performance measure, can be evaluated as

$`m={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\xi _i\text{sign}S_i_\beta \mathrm{}_{𝝃,𝜻,𝒟}={\displaystyle 𝑑z\varphi (z)\text{sign}(z)},`$ (31)
where

$`\varphi (z)={\displaystyle \left[\underset{l=1}{\overset{C}{}}d\widehat{x}_l\widehat{\pi }(\widehat{x}_l)\right]\delta \left(z\mathrm{tanh}\left(F_s\xi +\underset{i=1}{\overset{C}{}}\mathrm{tanh}^1\widehat{x}_i\right)\right)_\xi }.`$ (32)
The derivation of Eqs.(31) and (32) is given in Appendix B.

Examining the physical properties of the solutions for various connectivity values exposes significant differences between the various cases. In particular, these solutions fall into three different categories: the cases of $`K=1`$ and general $`L`$ value, the case of $`K=L=2`$ and all other parameter values where either $`K3`$ or $`L3`$ (and $`K>1`$). We describe the results obtained for each one of these cases separately.

### A Analytical solution - the case of $`K3`$ or $`L3`$, $`K>1`$

Results for the cases of $`K3`$ or $`L3`$, $`K>1`$ can be obtained analytically and have a simple and transparent interpretation; we will therefore focus first on this simple case. For unbiased messages (with $`F_s=0`$), one can easily verify that the ferromagnetic phase, characterised by $`m=1`$, and the probability distributions

$`\pi (x)=\delta (x1),\widehat{\pi }(\widehat{x})=\delta (\widehat{x}1),\rho (y)=\delta (y1),\widehat{\rho }(\widehat{y})=\delta (\widehat{y}1);`$ (33)
and the paramagnetic state of $`m=0`$ with the probability distributions

$`\pi (x)`$ $`=`$ $`\delta (x),\widehat{\pi }(\widehat{x})=\delta (\widehat{x}),\widehat{\rho }(\widehat{y})=\delta (\widehat{y}),`$ (34)
$`\rho (y)`$ $`=`$ $`{\displaystyle \frac{1+\mathrm{tanh}F_n}{2}}\delta (y\mathrm{tanh}F_n)+{\displaystyle \frac{1\mathrm{tanh}F_n}{2}}\delta (y+\mathrm{tanh}F_n),`$ (35)
satisfy saddle point equations (30). Other solutions may be obtained numerically; here we have represented the distributions by $`10^310^4`$ bins and iterated Eqs.(30) $`100500`$ times with $`10^5`$ Monte Carlo sampling steps for each iteration. No solutions other than the above two have been discovered.

The thermodynamically-dominant state is found by evaluating the free energy of the two solutions using Eq.(25), which yields

$`f_{\text{ferro}}={\displaystyle \frac{C}{K}}F_n\mathrm{tanh}F_n={\displaystyle \frac{1}{R}}F_n\mathrm{tanh}F_n,`$ (36)
for the ferromagnetic solution and

$`f_{\text{para}}={\displaystyle \frac{C}{K}}\mathrm{ln}2\mathrm{ln}2{\displaystyle \frac{C}{K}}\mathrm{ln}2\mathrm{cosh}F_n={\displaystyle \frac{1}{R}}\mathrm{ln}2\mathrm{ln}2{\displaystyle \frac{1}{R}}\mathrm{ln}2\mathrm{cosh}F_n,`$ (37)
for the paramagnetic solution.

Figure 1(a) describes schematically the nature of the solutions for this case, in terms of the free energy and the magnetisation obtained, for various flip rate probabilities. The difference between the free energies of Eqs.(36) and (37)

$`f_{\text{ferro}}f_{\text{para}}={\displaystyle \frac{\mathrm{ln}2}{R}}\left[R1+H_2(p)\right],`$ (38)
vanishes in the boundary between the two phase

$`R_c=1H_2(p)=1+p\mathrm{log}_2(p)+(1p)\mathrm{log}_2(1p),`$ (39)
which coincides with Shannon’s channel capacity.

Equation (39) indicates that all constructions with either $`K3`$ or $`L3`$ (and $`K>1`$) can potentially realize error-free data transmission for $`R<R_c`$ in the limit where both message and codeword lengths $`N`$ and $`M`$ become infinite, thus saturating Shannon’s bound.

### B The case of $`K=L=2`$

All codes with either $`K=3`$ or $`L=3`$, $`K>1`$ potentially saturate Shannon’s bound and are characterised by a first order phase transition between the ferromagnetic and paramagnetic solutions. On the other hand, numerical investigation based on Monte Carlo methods indicates of significantly different physical characteristics for $`K=L=2`$ codes shown in Fig.1(b).

At the highest noise level, the paramagnetic solution (35) gives the unique extremum of the free energy until noise level reaches the first critical point $`p_1`$, at which the ferromagnetic solution (33) of higher free energy appears to be locally stable. As the noise level decreases, a second critical point $`p_2`$ appears, where the paramagnetic solution becomes unstable and a sub-optimal ferromagnetic solution and its mirror image emerge. Those solutions have lower free energy than the ferromagnetic solution until the noise level reaches the third critical point $`p_3`$. Below $`p_3`$, the ferromagnetic solution becomes the global minimum of the free energy, while the sub-optimal ferromagnetic solutions still remain locally stable. However, the sub-optimal solutions disappear at the spinodal point $`p_s`$ and the ferromagnetic solution (and its mirror image) becomes the unique stable solution of the saddle point Eqs.(30) as shown by the numerical investigation for all $`p<p_s`$.

The analysis implies that $`p_3`$, the critical noise level below which the ferromagnetic solution becomes thermodynamically dominant, is lower than $`p_c=H_2^1(1R)`$ which corresponds to Shannon’s bound. Namely, $`K=L=2`$ does not saturate Shannon’s bound in contrast to $`K3`$ codes even if optimally decoded. Nevertheless, it turns out that the free energy landscape, for noise levels $`0<p<p_s`$, offers significant advantages in the decoding dynamic comparing to that of other codes ($`K3`$ or $`L3`$, $`K>1`$).

### C General $`L`$ codes with $`K=1`$

The particular choice of $`K=1`$, independently of the value chosen for $`L`$, exhibits a different behaviour presented schematically in Fig.1(c); also in this case there are no simple analytical solutions and all solutions in this scenario, except for the ferromagnetic solution, have been obtained numerically. The first important difference to be noted is that the paramagnetic state (35) is no longer a solution of the saddle point equations (30) and is being replaced by a sub-optimal ferromagnetic state. Convergence to the perfect solution of $`m=1`$ can only be guaranteed for corruption rates smaller than that of the spinodal point, marking the maximal noise level for which only the ferromagnetic solution exists, $`p<p_s`$.

The $`K=1`$ codes do not saturate Shannon’s bound in general; however, we have found that at rates $`R<1/3`$ they outperform the $`K=L=2`$ code while offering slightly improved dynamical (decoding) properties. Studying the free energy in this case shows that as the corruption rate increases, sub-optimal ferromagnetic solutions (stable and unstable) emerge at the spinodal point $`p_s`$. When the noise increases further this sub-optimal state becomes the global minimum at $`p_1`$, dominating the system’s thermodynamics. The transition at $`p_1`$ must occur at noise levels lower or equal to the value predicted by Shannon’s bound. In Fig.2 we show free energy values computed for a given code rate and several values of $`L`$, marking Shannon’s bound by a dashed line; it is clear that the thermodynamical transition observed numerically (i.e. the point where the ferromagnetic free energy equals the sub-optimal ferromagnetic free energy) is bellow, but very close, to the channel capacity. It implies that these codes also do not quite saturate Shannon’s bound if optimally decoded but get quite close to it.

## V Decoding: Belief propagation/TAP approach

The Bayesian message estimate (16) potentially provides the optimal retrieval of the original messages. However, it is computationally difficult to follow the prescription exactly as it requires a sum over $`𝒪(2^N)`$ terms. Belief propagation (BP) can be used for obtaining an approximate estimate. It was recently shown that the BP algorithm can be derived, at least in the current context, from the TAP approach to diluted systems in statistical mechanics.

Both algorithms (BP/TAP) are iterative methods which effectively calculate the marginal posterior probabilities $`𝒫(S_i|𝑱,C_s,C_n)=_{\{\{S_{ki}\},𝝉\}}𝒫(𝑺,𝝉|𝑱,C_s,C_n)`$ and $`𝒫(\tau _j|𝑱,C_s,C_n)=_{\{𝑺,\{\tau _{kj}\}\}}𝒫(𝑺,𝝉|𝑱,C_s,C_n)`$ based on the following three assumptions:

1. The posterior distribution is factorizable with respect to dynamical variables $`S_{i=1,\mathrm{},N}`$ and $`\tau _{j=1,\mathrm{},M}`$.
2. The influence of check $`J_{\mu =1,\mathrm{},M}`$ on a specific site $`S_i`$ (or $`\tau _j`$) is also factorizable.
3. The contribution of a single variables $`S_{i=1,\mathrm{},N}`$, $`\tau _{j=1,\mathrm{},M}`$ and $`J_{\mu =1,\mathrm{},M}`$ to the macroscopic variables is small and can be isolated.

Parameterising pseudo-marginal posteriors and marginalized conditional probabilities as

$`𝒫(S_i|\{J_{\nu \mu }\},C_s,C_n)`$ $`=`$ $`{\displaystyle \frac{1+m_{\mu i}^SS_i}{2}},𝒫(\tau _j|\{J_{\nu \mu }\},C_s,C_n)={\displaystyle \frac{1+m_{\mu j}^n\tau _j}{2}},`$ (40)
$`𝒫(J_\mu |S_i,\{J_{\nu \mu }\},C_s,C_n)`$ $``$ $`{\displaystyle \frac{1+\widehat{m}_{\mu i}^SS_i}{2}},𝒫(J_\mu |\tau _j,\{J_{\nu \mu }\},C_s,C_n){\displaystyle \frac{1+\widehat{m}_{\mu j}^n\tau _j}{2}},`$ (41)
the above assumptions provide a set of self-consistent equations

$`m_{\mu l}^S=\mathrm{tanh}\left(F_s+{\displaystyle \underset{\nu _S(l)/\mu }{}}\mathrm{tanh}^1(\widehat{m}_{\nu l}^S)\right),m_{\mu l}^n=\mathrm{tanh}\left(F_n+{\displaystyle \underset{\nu _n(l)/\mu }{}}\mathrm{tanh}^1(\widehat{m}_{\nu l}^n)\right).`$ (42)
and

$`\widehat{m}_{\mu l}^S=J_\mu {\displaystyle \underset{k_S(\mu )/l}{}}m_{\mu k}^S{\displaystyle \underset{j_n(\mu )}{}}m_{\mu j}^n,\widehat{m}_{\mu l}^n=J_\mu {\displaystyle \underset{k_S(\mu )}{}}m_{\mu k}^S{\displaystyle \underset{j_n(\mu )/l}{}}m_{\mu j}^n.`$ (43)
Here, $`_s(l)`$ and $`_n(l)`$ indicate the set of all indices of non-zero components in the $`l`$-th column of the sparse matrices $`C_s`$ and $`C_n`$, respectively. Similarly, $`_s(\mu )`$ and $`_n(\mu )`$ denote the set of all indices of non-zero components in $`\mu `$-th row of the sparse matrices $`C_s`$ and $`C_n`$, respectively. The notation $`_s(\mu )/l`$ represents the set of all indices belonging to $`_s(\mu )`$ except the index $`l`$.

Equations (42) and (43) are solved iteratively using the appropriate initial conditions. After obtaining a solution to all $`m_{\mu l}`$ and $`\widehat{m}_{\mu l}`$, an approximated posterior mean can be calculated as

$`m_i^S=\mathrm{tanh}\left(F_s+{\displaystyle \underset{\mu _S(l)}{}}\mathrm{tanh}^1(\widehat{m}_{\mu i}^S)\right),`$ (44)
which provides an approximation to the Bayes-optimal estimator (16) in the form of $`\widehat{\xi }^B=\text{sign}(m_i^S)`$.

Notice that the rather vague meaning of the fields distributions introduced in the previous section becomes clear by introducing the new variables $`x=\xi _im_{\mu i}^S`$, $`\widehat{x}=\xi _i\widehat{m}_{\mu i}^S`$, $`y=\zeta _jm_{\mu j}^n`$ and $`\widehat{y}=\zeta _j\widehat{m}_{\mu j}^n`$. If one considers that these variables are independently drawn from the distributions $`\pi (x)`$, $`\widehat{\pi }(\widehat{x})`$, $`\rho (y)`$ and $`\widehat{\rho }(\widehat{y})`$, the replica symmetric saddle point equations (30) are recovered from the BP/TAP equations (42) and (43). This connection can be extended to the free energy as equations (42) and (43) extremize the TAP free energy

$`f_{\text{TAP}}(\{𝒎\},\{\widehat{𝒎}\})`$ $`=`$ $`{\displaystyle \frac{M}{N}}\mathrm{ln}2+{\displaystyle \frac{1}{N}}{\displaystyle \underset{\mu =1}{\overset{M}{}}}{\displaystyle \underset{i_S(\mu )}{}}\mathrm{ln}\left(1+m_{\mu i}^S\widehat{m}_{\mu i}^S\right)+{\displaystyle \frac{1}{N}}{\displaystyle \underset{\mu =1}{\overset{M}{}}}{\displaystyle \underset{j_n(\mu )}{}}\mathrm{ln}\left(1+m_{\mu j}^n\widehat{m}_{\mu j}^n\right)`$ (45)
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{\mu =1}{\overset{M}{}}}\mathrm{ln}\left(1+J_\mu {\displaystyle \underset{i_S(\mu )}{}}m_{\mu i}^S{\displaystyle \underset{j_n(\mu )}{}}m_{\mu j}^n\right)`$ (46)
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}\left[e^{F_s}{\displaystyle \underset{\mu _S(i)}{}}\left(1+\widehat{m}_{\mu i}^S\right)+e^{F_s}{\displaystyle \underset{\mu _S(i)}{}}\left(1\widehat{m}_{\mu i}^S\right)\right]`$ (47)
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j=1}{\overset{M}{}}}\mathrm{ln}\left[e^{F_n}{\displaystyle \underset{\mu _n(j)}{}}\left(1+\widehat{m}_{\mu j}^n\right)+e^{F_n}{\displaystyle \underset{\mu _n(j)}{}}\left(1\widehat{m}_{\mu j}^n\right)\right].`$ (48)
This expression may be used for selecting the thermodynamically dominant state when Eqs.(42) and (43) have several solutions.

We have investigated the performance of the various codes using BP/TAP equations as the decoding algorithm. Solutions have been obtained by iterating the equations (42) and (43) $`100500`$ times under various initial conditions. Since the system is not frustrated, the dynamics converges within $`1030`$ updates in most cases except close to criticality. The numerical results mirror the behaviour predicted by the analytical solutions.

For either $`K3`$ or $`L3`$ ,$`K>1`$ codes, the ferromagnetic solution

$`m_{\mu i}^S=\xi _i,\widehat{m}_{\mu i}^S=\xi _i,m_{\mu j}^n=\zeta _j,\widehat{m}_{\mu j}^n=\zeta _j,`$ (49)
which provides perfect decoding ($`m=1`$) and the paramagnetic solution ($`m=0`$)

$`m_{\mu i}^S=0,\widehat{m}_{\mu i}^S=0,m_{\mu j}^n=\mathrm{tanh}F_n=12p,\widehat{m}_{\mu j}^n=0,`$ (50)
are obtained in various runs depending on the initial conditions (the message is assumed unbiased resulting in $`F_s=0`$). However, it is difficult to set the initial conditions within the basin of attraction of the ferromagnetic solution without prior knowledge about the transmitted message $`𝝃`$.

Biased coding is sometimes used for alleviating this difficulty . Using a redundant source of information (equivalent to the introduction of a non-zero field $`F_s`$ in the statistical physics description), one effectively increases the probability of the initial conditions being closer to the ferromagnetic solution. The main drawback of this method is that the information per transmitted bit is significantly reduced due to this redundancy. In order to investigate how the maximum performance is affected by transmitting biased messages, we have evaluated the critical information rate (i.e., code rate $`\times H_2(f_s=(1+\mathrm{tanh}F_s)/2)`$, the source redundancy), below which the ferromagnetic solution becomes thermodynamically dominant \[Fig.3(a)\]. The data were obtained by the BP/TAP method (diamonds) and numerical solutions of from replica framework (square); the dominant solution in the BP/TAP results, was selected by using the free energy (48). Numerical solutions have been obtained using $`10^310^4`$ bin models for each distribution and had been run for $`10^5`$ steps per noise level. The various results are highly consistent and practically saturate Shannon’s bound for the same noise level. However, it is important to point out that close to Shannon’s limit, prior knowledge on the original message is required for setting up appropriate initial conditions that ensure convergence to the ferromagnetic solution; such prior knowledge is not available in practice.

Although $`K,L3`$ codes seem to offer optimal performance when highly biased messages are transmitted, this seems to be of little relevance in most cases, characterised by the transmission of compressed unbiased messages or only slightly biased messages. In this sense, $`K=L=2`$ and $`K=1`$ codes can be considered more practical as the BP/TAP dynamics of these codes exhibit unique convergence to the ferromagnetic solution (or mirror image in the $`K=L=2`$ case) from any initial condition up to a certain noise level. This property results from the fact that the corresponding free energies have no local minima other than the ferromagnetic solution below $`p_s`$.

In figures 3(b) and (c) we show the value of $`p_s`$ for the cases of $`K=L=2`$ and $`K=1`$, $`L=2`$ respectively, evaluated by numerical solutions from the replica framework (diamonds) and using the BP/TAP method.

The case of $`K=L=2`$ shows consistent successful decoding for the code rates examined and up to noise levels slightly below, but close to, Shannon’s bound. It should be emphasised here that initial conditions are chosen almost randomly in the BP/TAP method, with a very slight bias of $`𝒪(10^{12})`$ in the initial magnetisation. This result suggests using $`K=L=2`$ codes (or similar), rather than $`K,L3`$ codes, although the latter may potentially have better equilibrium properties.

In Fig.3(c) we show that for code rates $`R<1/3`$, codes parametrised by $`K=1`$ and $`L=2`$ outperform $`K=L=2`$ codes with one additional advantage: Due to the absence of mirror symmetries these codes converge to the ferromagnetic state much faster, and there is no risk of convergence to the mirror solution. The difference in performance becomes even larger as the code rate decreases. Higher code rates will result in performance deterioration due to the low connectivity, eventually bringing the system below the percolation threshold.

In Fig.4 we examine the dependence of the noise level of the spinodal point $`p_s`$ on the value of $`L`$, and show that the choice of $`L=2`$ is optimal within this family. Codes with $`L=1`$ have very poor error-correction capabilities as their Hamiltonian (11) corresponds to the Mattis model, which is equivalent to a simple ferromagnet in a random field attaining magnetisation $`m=1`$ only in the noiseless case.

## VI Reducing Encoding Costs

The BP/TAP algorithm already offers an efficient decoding method, which requires $`𝒪(N)`$ operations; however, the current encoding scheme includes three costly processes: (a) The computational cost of constructing the generating matrix $`C_n^1C_s`$ requires $`𝒪(N^3)`$ operations for inverting the matrix $`C_n`$ and $`𝒪(N^2)`$ operations for the matrix multiplication. (b) The memory allocation for generating the matrix $`C_n^1C_s`$ scales as $`𝒪(N^2)`$ since this matrix is typically dense. (c) The encoding itself $`𝒛_0=C_n^1C_s𝝃`$ (mod $`2`$) requires $`𝒪(N^2)`$ operations.

These computational costs become significant when long messages $`N=10^410^5`$ are transmitted, which is typically the case for which Gallager-type codes are being used. This may require long encoding times and may delay the transmission.

These problems may be solved by utilising systematically constructed matrices instead of random ones, of some similarity to the constructions of . Here, we present a simple method to reduce the computational and memory costs to $`𝒪(N)`$ for $`K=L=2`$ and $`K=1`$, $`L=2`$ codes. Our proposal is mainly based on using a specific matrix for $`C_n`$,

$`\overline{C}_n=\left(\begin{array}{cccccccc}1& 0& 0& 0& \mathrm{}& 0& 0& \\ 1& 1& 0& 0& \mathrm{}& 0& 0& \\ 0& 1& 1& 0& \mathrm{}& 0& 0& \\ 0& 0& 1& 1& \mathrm{}& 0& 0& \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ 0& 0& 0& 0& \mathrm{}& 1& 1& \end{array}\right),`$ (57)
instead of a randomly-constructed one. For $`C_s`$, we use a random matrix of $`K=2`$ (or $`K=1`$) non-zero elements per row as before.

The inverse (mod 2) of $`\overline{C}_n^1`$ becomes the lower triangular matrix

$`\overline{C}_n^1=\left(\begin{array}{cccccccc}1& 0& 0& 0& \mathrm{}& 0& 0& \\ 1& 1& 0& 0& \mathrm{}& 0& 0& \\ 1& 1& 1& 0& \mathrm{}& 0& 0& \\ 1& 1& 1& 1& \mathrm{}& 0& 0& \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ 1& 1& 1& 1& \mathrm{}& 1& 1& \end{array}\right).`$ (64)
This suggests that encoding the message $`𝝃`$ into a codeword $`𝒛^0`$ would require only $`𝒪(N)`$ operations by carrying it out in two steps

$`t_\mu `$ $`=`$ $`(C_s𝝃)_\mu \text{(mod }2\text{)},\text{for }\mu =1,2,\mathrm{},M,`$ (65)
$`z_\mu ^0`$ $`=`$ $`(\overline{C}_n^1𝒕)_\mu =z_{\mu 1}^0+t_\mu \text{(mod }2\text{)},\text{for }\mu =2,\mathrm{},M,`$ (66)
with $`z_1^0=t_1`$. Both steps require $`𝒪(N)`$ operations due to the sparse nature of $`C_s`$. In addition, the required memory resources are also reduced to $`𝒪(N)`$ since only the sparse matrix $`C_s`$ should be stored.

The possible drawback of using the systematic matrix (57) is a deterioration in the error correction ability. We have examined numerically the performance of new construction to discover, to our surprise, that it is very similar to that of random matrix based codes as shown in Table I. Although our examination is only limited to BSC and i.i.d. messages, it seems to suggest that some deterministically constructed matrices may be implemented successfully in practice.

## VII Summary

In this paper, we have investigated the typical performance of the MN codes, a variation of Gallager-type error-correcting codes, by mapping them onto Ising spin models and making use of the established methods of statistical physics. We have discovered that for a certain choice of parameters, either $`K3`$ or $`L3`$, $`K>1`$ these codes potentially saturate the channel capacity, although this cannot be used efficiently in practice due to the decrease in the basin of attraction which typically diverts the decoding dynamics towards the undesired paramagnetic solution.

Codes with $`K=2`$ and $`L=2`$ show close to optimal performance while keeping a large basin of attraction, resulting in more practical codes. Constructions of the form $`K=1`$, $`L=2`$ outperform the $`K=L=2`$ codes for code rates $`R<1/3`$, having improved dynamical properties.

These results are complementary to those obtained so far by the information theory community and seem to indicate that worst-case analysis can be, in some situations, too pessimistic when compared to the typical performance results.

Beyond the theoretical aspects, we proposed an efficient method for reducing the computational costs and the required memory allocation by using a specific construction of the matrix $`C_n`$. These codes are highly attractive and provide lower computational costs for both encoding and decoding.

Various aspects that remain to be studied include a proper analysis of the finite size effects for rates below and above the channel capacity, which are of great practical relevance; and the use of statistical physics methods for optimising the matrix constructions.

## Acknowledgement

Support by the JSPS RFTF program (YK), The Royal Society and EPSRC grant GR/N00562 (DS) is acknowledged.

## A Replica Free Energy

The purpose of this appendix is to derive the averaged free energy per spin (25). Applying the gauge transformation

$`J_\mu `$ $``$ $`J_\mu {\displaystyle \underset{i_s(\mu )}{}}\xi _i{\displaystyle \underset{j_n(\mu )}{}}\zeta _j=1`$ (A1)
$`S_i`$ $``$ $`S_i\xi _i`$ (A2)
$`\tau _j`$ $``$ $`\tau _j\zeta _j,`$ (A3)
to eq. (14), one may rewrite the partition function in the form

$`𝒵(𝝃,𝜻,𝒟)={\displaystyle \underset{𝑺,𝝉}{}}\mathrm{exp}\left(F_s{\displaystyle \underset{i=1}{\overset{N}{}}}\xi _iS_i+F_n{\displaystyle \underset{j=1}{\overset{M}{}}}\zeta _j\tau _j\right)`$ (A4)
$`\times `$ $`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left[1𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}+𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{\displaystyle \frac{1}{2}}\left(1+S_{i_1}\mathrm{}S_{i_K}\tau _{j_1}\mathrm{}\tau _{j_L}\right)\right].`$ (A5)
Using the replica method, one calculates the quenched average of the $`n`$-th power of the partition function given by

$`𝒵(𝝃,𝜻,𝒟)^n_{𝝃,𝜻,𝒟}={\displaystyle \underset{𝑺^1\mathrm{}𝑺^n}{}}{\displaystyle \underset{𝝉^1\mathrm{}𝝉^n}{}}\mathrm{exp}\left(F_s{\displaystyle \underset{i=1}{\overset{N}{}}}\xi _i{\displaystyle \underset{\alpha =1}{\overset{n}{}}}S_i^\alpha \right)_𝝃\mathrm{exp}\left(F_n{\displaystyle \underset{j=1}{\overset{M}{}}}\zeta _j{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\tau _i^\alpha \right)_𝜻`$ (A6)
$`\times `$ $`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\left\{1+{\displaystyle \frac{1}{2}}𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}\left(S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha 1\right)\right\}_𝒟,`$ (A7)
where averages with respect to $`𝝃`$ can be easily performed

$`\mathrm{exp}\left(F_s{\displaystyle \underset{i=1}{\overset{N}{}}}\xi _i{\displaystyle \underset{\alpha =1}{\overset{n}{}}}S_i^\alpha \right)_𝝃`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left[\left({\displaystyle \frac{1+\mathrm{tanh}F_s}{2}}\right)e^{F_s_{\alpha =1}^nS_i^\alpha }+\left({\displaystyle \frac{1\mathrm{tanh}F_s}{2}}\right)e^{F_s_{\alpha =1}^nS_i^\alpha }\right]`$ (A8)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}\left(\xi F_s{\displaystyle \underset{a=1}{\overset{n}{}}}S_i^a\right)_\xi ,`$ (A9)
and similarly for $`\mathrm{}_𝜻`$. The main problem is in averages over the sparse tensor realisations $`𝒟`$, which have complicated constraints. Following the procedure introduced by Wong and Sherrington, it is being rewritten as

$`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\left[1+{\displaystyle \frac{1}{2}}𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}\left(S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha 1\right)\right]_𝒟`$ (A10)
$`=`$ $`𝒩^1{\displaystyle \underset{𝒟}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\delta \left({\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}𝒟_{i,i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}C\right){\displaystyle \underset{j=1}{\overset{M}{}}}\delta \left({\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}𝒟_{i_1,\mathrm{},i_K;j,j_2,\mathrm{},j_L}L\right)`$ (A12)
$`\times {\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}{\displaystyle \underset{\alpha =1}{\overset{n}{}}}[1+{\displaystyle \frac{1}{2}}𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}(S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha 1)],`$
where $`\delta (\mathrm{})`$ represents Dirac’s $`\delta `$-function and

$`𝒩={\displaystyle \underset{𝒟}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\delta \left({\displaystyle \underset{i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}𝒟_{i,i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}C\right){\displaystyle \underset{j=1}{\overset{M}{}}}\delta \left({\displaystyle \underset{i_1,\mathrm{},i_K;j_2,\mathrm{},j_L}{}}𝒟_{i_1,\mathrm{},i_K;j,j_2,\mathrm{},j_L}L\right)`$ (A13)
represents the normalisation constant.

We first evaluate this normalisation constant using the integral representation of the $`\delta `$-function and Eq.(A13), to obtain

$`𝒩`$ $`=`$ $`{\displaystyle \underset{𝒟}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\delta \left({\displaystyle \underset{i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}𝒟_{i,i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}C\right){\displaystyle \underset{j=1}{\overset{M}{}}}\delta \left({\displaystyle \underset{i_1,\mathrm{},i_K;j_2,\mathrm{},j_L}{}}𝒟_{i_1,\mathrm{},i_K;j,j_2,\mathrm{},j_L}L\right)`$ (A14)
$`=`$ $`{\displaystyle \underset{𝒟}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\lambda _i}{2\pi }}\mathrm{exp}\left[i\lambda _i\left({\displaystyle \underset{i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}𝒟_{i,i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}C\right)\right]\right\}`$ (A16)
$`\times {\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\lambda _j}{2\pi }}\mathrm{exp}\left[i\lambda _j({\displaystyle \underset{i_1,\mathrm{},i_K;j_2,\mathrm{},j_L}{}}𝒟_{i_1,\mathrm{},i_K;j,j_2,\mathrm{},j_L}L)\right]\right\}`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\lambda _i}{2\pi }}e^{iC\lambda _i}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\nu _j}{2\pi }}e^{iL\nu _j}\right\}`$ (A17)
$`\times {\displaystyle \underset{𝒟}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \underset{i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}e^{i\lambda _i𝒟_{i,i_2,\mathrm{},i_K;j_1,\mathrm{},j_L}}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle \underset{i_1,\mathrm{},i_K;j_2,\mathrm{},j_L}{}}e^{i\nu _j𝒟_{i_1,\mathrm{},i_K;j,j_2,\mathrm{},j_L}}\right\}`$ (A18)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\lambda _i}{2\pi }}e^{iC\lambda _i}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\nu _j}{2\pi }}e^{iL\nu _j}\right\}`$ (A19)
$`\times {\displaystyle \underset{𝒟}{}}{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left\{\left(e^{i\lambda _{i_1}}\mathrm{}e^{i\lambda _{i_K}}e^{i\nu _{j_1}}\mathrm{}e^{i\nu _{j_L}}\right)^{𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}}\right\}`$ (A20)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \frac{dZ_i}{2\pi i}Z_i^{(C+1)}}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle \frac{dY_j}{2\pi i}Y_j^{(L+1)}}\right\}{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left(1+Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right),`$ (A21)
where we made use of the transformations $`Z_i=e^{i\lambda _i},Y_j=e^{i\nu _j}`$, and carried out summations with respect to the realisation of $`𝒟`$. Expanding the product on the right hand side one obtains

$`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left[1+\left(Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right)\right]`$ (A22)
$`=`$ $`\mathrm{exp}\left[{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\mathrm{ln}\left\{1+\left(Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right)\right\}\right]`$ (A23)
$``$ $`\mathrm{exp}\left[{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left(Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right)\right]`$ (A24)
$``$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{K!}}\left({\displaystyle \underset{i=1}{\overset{N}{}}}Z_i\right)^K{\displaystyle \frac{1}{L!}}\left({\displaystyle \underset{j=1}{\overset{M}{}}}Y_j\right)^L\right],`$ (A25)
in the thermodynamic limit. Using the identities

$`1={\displaystyle 𝑑q\delta \left(\underset{i=1}{\overset{N}{}}Z_iq\right)},1={\displaystyle 𝑑r\delta \left(\underset{j=1}{\overset{M}{}}Y_jr\right)}`$ (A26)
Eq. (A21) becomes

$`𝒩`$ $`=`$ $`{\displaystyle 𝑑q\delta \left(\underset{i=1}{\overset{N}{}}Z_iq\right)𝑑r\delta \left(\underset{j=1}{\overset{M}{}}Y_jr\right)}`$ (A27)
$`\times {\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle }{\displaystyle \frac{dZ_i}{2\pi i}}Z_i^{(C+1)}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle }{\displaystyle \frac{dY_j}{2\pi i}}Y_j^{(L+1)}\right\}\mathrm{exp}\left({\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}\right)`$ (A28)
$`=`$ $`{\displaystyle 𝑑q\frac{d\widehat{q}}{2\pi i}\mathrm{exp}\left[\widehat{q}\left(\underset{i=1}{\overset{N}{}}Z_iq\right)\right]𝑑r\frac{d\widehat{r}}{2\pi i}\mathrm{exp}\left[\widehat{r}\left(\underset{j=1}{\overset{M}{}}Y_jr\right)\right]}`$ (A29)
$`\times {\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle }{\displaystyle \frac{dZ_i}{2\pi i}}Z_i^{(C+1)}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle }{\displaystyle \frac{dY_j}{2\pi i}}Y_j^{(L+1)}\right\}\mathrm{exp}\left({\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}\right)`$ (A30)
$`=`$ $`{\displaystyle 𝑑q\frac{d\widehat{q}}{2\pi i}𝑑r\frac{d\widehat{r}}{2\pi i}\mathrm{exp}\left(\frac{q^K}{K!}\frac{r^L}{L!}q\widehat{q}r\widehat{r}\right)}`$ (A31)
$`\times {\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle }{\displaystyle \frac{dZ_i}{2\pi i}}Z_i^{(C+1)}\mathrm{exp}\left(\widehat{q}Z_i\right)\right]{\displaystyle \underset{j=1}{\overset{M}{}}}\left[{\displaystyle }{\displaystyle \frac{dY_j}{2\pi i}}Y_j^{(L+1)}\mathrm{exp}\left(\widehat{r}Y_j\right)\right].`$ (A32)
The contour integrals provide the following constants

$`{\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{dZ_i}{2\pi i}Z_i^{(C+1)}\mathrm{exp}\left(\widehat{q}Z_i\right)}\right]=\left({\displaystyle \frac{\widehat{q}^C}{C!}}\right)^N,{\displaystyle \underset{j=1}{\overset{M}{}}}\left[{\displaystyle \frac{dY_j}{2\pi i}Y_j^{(L+1)}\mathrm{exp}\left(\widehat{r}Y_j\right)}\right]=\left({\displaystyle \frac{\widehat{r}^L}{L!}}\right)^M,`$ (A33)
respectively. Applying the saddle point method to the remaining integrals, one obtains

$`𝒩=\text{extr}_{\{q,\widehat{q},r,\widehat{r}\}}\left\{\mathrm{exp}\left[{\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}q\widehat{q}r\widehat{r}+NC\mathrm{ln}\widehat{q}N\mathrm{ln}(C!)+ML\mathrm{ln}\widehat{r}M\mathrm{ln}(L!)\right]\right\},`$ (A34)
which yields the following saddle point equations with respect to $`q`$, $`r`$, $`\widehat{q}`$ and $`\widehat{r}`$

$`q`$ $`=`$ $`{\displaystyle \frac{NC}{\widehat{q}}},r={\displaystyle \frac{ML}{\widehat{r}}}`$ (A35)
$`\widehat{q}`$ $`=`$ $`{\displaystyle \frac{q^{K1}}{(K1)!}}{\displaystyle \frac{r^L}{L!}},\widehat{r}={\displaystyle \frac{r^{L1}}{(L1)!}}{\displaystyle \frac{q^K}{K!}},`$ (A36)
providing the normalisation constant

$`𝒩=\left({\displaystyle \frac{\widehat{q}^C}{C!}}\right)^N\left({\displaystyle \frac{\widehat{r}^L}{L!}}\right)^M\mathrm{exp}\left({\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}q\widehat{q}r\widehat{r}\right).`$ (A37)
Equation (A12) can be evaluated similarly. Following a similar calculation to that of Eq.(A21) provides

$`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\left\{1+{\displaystyle \frac{1}{2}}𝒟_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}\left(S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha 1\right)\right\}_𝒟`$ (A38)
$`=`$ $`𝒩^1{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \frac{dZ_i}{2\pi i}Z_i^{(C+1)}}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle \frac{dY_j}{2\pi i}Y_j^{(L+1)}}\right\}`$ (A39)
$`\times {\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}[1+\left(Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right){\displaystyle \underset{\alpha =1}{\overset{n}{}}}{\displaystyle \frac{1}{2}}(1+S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha )].`$ (A40)
Using the expansion

$`{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\left(1+S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha \right)`$ (A41)
$`=`$ $`1+{\displaystyle \underset{\alpha =1}{\overset{n}{}}}S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha +{\displaystyle \underset{\alpha _1,\alpha _2}{}}(S_{i_1}^{\alpha _1}S_{i_1}^{\alpha _2})\mathrm{}(S_{i_K}^{\alpha _1}S_{i_K}^{\alpha _2})(\tau _{j_1}^{\alpha _1}\tau _{j_1}^{\alpha _2})\mathrm{}(\tau _{j_L}^{\alpha _1}\tau _{j_L}^{\alpha _2})`$ (A42)
$`+\mathrm{}+{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _n}{}}(S_{i_1}^{\alpha _1}\mathrm{}S_{i_1}^{\alpha _n})\mathrm{}(S_{i_K}^{\alpha _1}\mathrm{}S_{i_K}^{\alpha _n})(\tau _{j_1}^{\alpha _1}\mathrm{}\tau _{j_1}^{\alpha _n})\mathrm{}(\tau _{j_L}^{\alpha _1}\mathrm{}\tau _{j_L}^{\alpha _n})`$ (A43)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}(S_{i_1}^{\alpha _1}\mathrm{}S_{i_1}^{\alpha _m})\mathrm{}(S_{i_K}^{\alpha _1}\mathrm{}S_{i_K}^{\alpha _m})(\tau _{j_1}^{\alpha _1}\mathrm{}\tau _{j_1}^{\alpha _m})\mathrm{}(\tau _{j_L}^{\alpha _1}\mathrm{}\tau _{j_L}^{\alpha _m}),`$ (A44)
resulting in

$`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left[1+\left(Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right){\displaystyle \underset{\alpha =1}{\overset{n}{}}}{\displaystyle \frac{1}{2}}\left(1+S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha \right)\right]`$ (A45)
$``$ $`e^{_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}_{\alpha =1}^n\frac{1}{2}\left(1+S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha \right)}`$ (A46)
$`=`$ $`e^{\frac{1}{2^n}_{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}_{m=0}^n_{\alpha _1,\mathrm{},\alpha _m}(S_{i_1}^{\alpha _1}\mathrm{}S_{i_1}^{\alpha _m})\mathrm{}(S_{i_K}^{\alpha _1}\mathrm{}S_{i_K}^{\alpha _m})(\tau _{j_1}^{\alpha _1}\mathrm{}\tau _{j_1}^{\alpha _m})\mathrm{}(\tau _{j_L}^{\alpha _1}\mathrm{}\tau _{j_L}^{\alpha _m})}`$ (A47)
$`=`$ $`e^{\frac{1}{2^n}\left\{_{m=0}^n_{\alpha _1,\mathrm{},\alpha _m}_{i_1,\mathrm{},i_K}\left(S_{i_1}^{\alpha _1}\mathrm{}S_{i_1}^{\alpha _m}Z_{i_1}\right)\mathrm{}\left(S_{i_K}^{\alpha _1}\mathrm{}S_{i_K}^{\alpha _m}Z_{i_K}\right)_{j_1,\mathrm{},j_L}\left(\tau _{j_1}^{\alpha _1}\mathrm{}\tau _{j_1}^{\alpha _m}Y_{j_1}\right)\mathrm{}\left(\tau _{j_L}^{\alpha _1}\mathrm{}\tau _{j_K}^{\alpha _m}Y_{j_L}\right)\right\}}`$ (A48)
$``$ $`e^{\frac{1}{2^n}\left\{_{m=0}^n_{\alpha _1,\mathrm{},\alpha _m}\frac{1}{K!}\left(_{i=1}^NS_i^{\alpha _1}\mathrm{}S_i^{\alpha _m}Z_i\right)^K\frac{1}{L!}\left(\tau _j^{\alpha _1}\mathrm{}\tau _j^{\alpha _m}Y_j\right)^L\right\}}.`$ (A49)
Using the identities

$`1`$ $`=`$ $`{\displaystyle 𝑑q_{\alpha _1,\mathrm{},\alpha _m}\delta \left(\underset{i=1}{\overset{N}{}}S_i^{\alpha _1}\mathrm{}S_i^{\alpha _m}Z_iq_{\alpha _1,\mathrm{},\alpha _m}\right)},`$ (A50)
$`1`$ $`=`$ $`{\displaystyle 𝑑r_{\alpha _1,\mathrm{},\alpha _m}\delta \left(\underset{j=1}{\overset{M}{}}\tau _j^{\alpha _1}\mathrm{}\tau _j^{\alpha _m}Y_jr_{\alpha _1,\mathrm{},\alpha _m}\right)}`$ (A51)
and going through the same steps as in Eqs. (A26 \- A34), we arrive at

$`{\displaystyle \underset{i_1,\mathrm{},i_K;j_1,\mathrm{},j_L}{}}\left[1+\left(Z_{i_1}\mathrm{}Z_{i_K}Y_{j_1}\mathrm{}Y_{j_L}\right){\displaystyle \underset{\alpha =1}{\overset{n}{}}}{\displaystyle \frac{1}{2}}\left(1+S_{i_1}^\alpha \mathrm{}S_{i_K}^\alpha \tau _{j_1}^\alpha \mathrm{}\tau _{j_L}^\alpha \right)\right]`$ (A52)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}{\displaystyle 𝑑q_{\alpha _1,\mathrm{},\alpha _m}\delta \left(\underset{i=1}{\overset{N}{}}S_i^{\alpha _1}\mathrm{}S_i^{\alpha _m}Z_iq_{\alpha 1,\mathrm{},\alpha _m}\right)}`$ (A53)
$`\times {\displaystyle }dr_{\alpha _1,\mathrm{},\alpha _m}\delta ({\displaystyle \underset{j=1}{\overset{M}{}}}\tau _j^{\alpha _1}\mathrm{}\tau _j^{\alpha _m}Y_jr_{\alpha 1,\mathrm{},\alpha _m})\mathrm{exp}\left({\displaystyle \frac{1}{2^n}}\left\{{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}{\displaystyle \frac{q_{\alpha _1,\mathrm{},\alpha _m}^K}{K!}}{\displaystyle \frac{r_{\alpha _1,\mathrm{},\alpha _m}^L}{L!}}\right\}\right)`$ (A54)
$``$ $`\text{extr}_{\{𝒒,\widehat{𝒒},𝒓,\widehat{𝒓}\}}\{\mathrm{exp}[{\displaystyle \frac{1}{2^n}}\left\{{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}{\displaystyle \frac{q_{\alpha _1,\mathrm{},\alpha _m}^K}{K!}}{\displaystyle \frac{r_{\alpha _1,\mathrm{},\alpha _m}^L}{L!}}\right\}`$ (A55)
$`{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}q_{\alpha _1,\mathrm{},\alpha _m}\widehat{q}_{\alpha _1,\mathrm{},\alpha _m}{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}r_{\alpha _1,\mathrm{},\alpha _m}\widehat{r}_{\alpha _1,\mathrm{},\alpha _m}`$ (A56)
$`+{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}\widehat{q}_{\alpha _1,\mathrm{},\alpha _m}{\displaystyle \underset{i=1}{\overset{N}{}}}S_i^{\alpha _1}\mathrm{}S_i^{\alpha _m}Z_i+{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}\widehat{r}_{\alpha _1,\mathrm{},\alpha _m}{\displaystyle \underset{j=1}{\overset{M}{}}}\tau _j^{\alpha _1}\mathrm{}\tau _j^{\alpha _m}Y_j]\}.`$ (A57)
In order to proceed further, one has to make an assumption about the order parameter symmetry. We adopt here the replica symmetric ansatz for the order parameters $`q`$, $`r`$, $`\widehat{q}`$ and $`\widehat{r}`$. This implies that the order parameters do not depend on the explicit indices but only on their number. It is therefore convenient to represent them as moments of random variables defined over the interval $`[1,1]`$

$`q_{\alpha _1,\mathrm{},\alpha _l}`$ $`=`$ $`q{\displaystyle 𝑑x\pi (x)x^l},r_{\alpha _1,\mathrm{},\alpha _l}=r{\displaystyle 𝑑y\rho (y)y^l},`$ (A58)
$`\widehat{q}_{\alpha _1,\mathrm{},\alpha _l}`$ $`=`$ $`\widehat{q}{\displaystyle 𝑑\widehat{x}\widehat{\pi }(\widehat{x})\widehat{x}^l},\widehat{r}_{\alpha _1,\mathrm{},\alpha _l}=\widehat{r}{\displaystyle 𝑑\widehat{y}\widehat{\rho }(\widehat{y})\widehat{y}^l},`$ (A59)
Then, each term in Eq.(A57) takes the form

$`{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}{\displaystyle \frac{q_{\alpha _1,\mathrm{},\alpha _m}^K}{K!}}{\displaystyle \frac{r_{\alpha _1,\mathrm{},\alpha _m}^L}{L!}}`$ $`=`$ $`{\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}{\displaystyle \underset{m=0}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{m}}\right){\displaystyle \underset{k=1}{\overset{K}{}}dx_k\pi (x_k)x_k^m\underset{l=1}{\overset{L}{}}dy_l\rho (y_l)y_l^m}`$ (A60)
$`=`$ $`{\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}{\displaystyle \underset{k=1}{\overset{K}{}}dx_k\pi (x_k)\underset{l=1}{\overset{L}{}}dy_l\rho (y_l)\left(1+\underset{k=1}{\overset{K}{}}x_k\underset{l=1}{\overset{L}{}}y_l\right)^n}`$ (A61)
$`{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}q_{\alpha _1,\mathrm{},\alpha _m}\widehat{q}_{\alpha _1,\mathrm{},\alpha _m}`$ $`=`$ $`q\widehat{q}{\displaystyle \underset{m=0}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{m}}\right){\displaystyle 𝑑x𝑑\widehat{x}\pi (x)\widehat{\pi }(\widehat{x})x^m\widehat{x}^m}`$ (A62)
$`=`$ $`q\widehat{q}{\displaystyle 𝑑x𝑑\widehat{x}\pi (x)\widehat{\pi }(\widehat{x})\left(1+x\widehat{x}\right)^n}`$ (A63)
$`{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \underset{\alpha _1,\mathrm{},\alpha _m}{}}\widehat{q}_{\alpha _1,\mathrm{},\alpha _m}{\displaystyle \underset{i=1}{\overset{N}{}}}S_i^{\alpha _1}\mathrm{}S_i^{\alpha _m}Z_i`$ $`=`$ $`\widehat{q}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_i{\displaystyle 𝑑\widehat{x}\widehat{\pi }(\widehat{x})\underset{m=0}{\overset{n}{}}\widehat{x}^m\underset{\alpha _1,\mathrm{},\alpha _m}{}S_i^{\alpha _1}\mathrm{}S_i^{\alpha _m}}`$ (A64)
$`=`$ $`\widehat{q}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_i{\displaystyle 𝑑\widehat{x}\widehat{\pi }(\widehat{x})\underset{\alpha =1}{\overset{n}{}}\left(1+S_i^\alpha \widehat{x}\right)}.`$ (A65)
Substituting these into (A57), one obtains

$`Z(𝝃,𝜻,𝒟)^n_{𝝃,𝜻,𝒟}`$ (A66)
$`=`$ $`{\displaystyle \underset{𝑺^1\mathrm{}𝑺^n}{}}{\displaystyle \underset{𝝉^1\mathrm{}𝝉^n}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}\left(\xi F_s{\displaystyle \underset{\alpha =1}{\overset{n}{}}}S_i^\alpha \right)_\xi \times {\displaystyle \underset{j=1}{\overset{M}{}}}\mathrm{exp}\left(\zeta F_n{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\tau _j^\alpha \right)_\zeta `$ (A67)
$`\times 𝒩^1{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \frac{dZ_i}{2\pi i}Z_i^{(C+1)}}\right\}{\displaystyle \underset{j=1}{\overset{M}{}}}\left\{{\displaystyle \frac{dY_j}{2\pi i}Y_j^{(L+1)}}\right\}`$ (A68)
$`\times \text{extr}_{\{\pi ,\widehat{\pi },\rho ,\widehat{\rho }\}}\{\mathrm{exp}[{\displaystyle \frac{1}{2^n}}\left\{{\displaystyle \frac{q^K}{K!}}{\displaystyle \frac{r^L}{L!}}{\displaystyle }{\displaystyle \underset{l=1}{\overset{K}{}}}dx_l\pi (x_l){\displaystyle }{\displaystyle \underset{l=1}{\overset{L}{}}}dy_l\rho (y_l)(1+{\displaystyle \underset{l=1}{\overset{K}{}}}x_l{\displaystyle \underset{l=1}{\overset{L}{}}}y_l)^n\right\}`$ (A69)
$`q\widehat{q}{\displaystyle 𝑑x𝑑\widehat{x}\pi (x)\widehat{\pi }(\widehat{x})(1+x\widehat{x})^n}r\widehat{r}{\displaystyle 𝑑y𝑑\widehat{y}\rho (y)\widehat{\rho }(\widehat{y})(1+y\widehat{y})^n}`$ (A70)
$`+\widehat{q}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_i{\displaystyle }d\widehat{x}\widehat{\pi }(\widehat{x}){\displaystyle \underset{\alpha =1}{\overset{n}{}}}(1+S_i^\alpha \widehat{x})+\widehat{r}{\displaystyle \underset{j=1}{\overset{M}{}}}Y_j{\displaystyle }d\widehat{y}\widehat{\rho }(\widehat{y}){\displaystyle \underset{\alpha =1}{\overset{n}{}}}(1+\tau _j^\alpha \widehat{y})]\}.`$ (A71)
The term involving the spin variables $`S`$ is easily evaluated using the residue theorem

$`{\displaystyle \underset{𝑺^1\mathrm{}𝑺^n}{}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}\left(\xi F_s{\displaystyle \underset{\alpha =1}{\overset{n}{}}}S_i^\alpha \right)_\xi {\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \frac{dZ_i}{2\pi i}Z_i^{(C+1)}}\right\}\times \mathrm{exp}\left[\widehat{q}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_i{\displaystyle 𝑑\widehat{x}\widehat{\pi }(\widehat{x})\underset{\alpha =1}{\overset{n}{}}(1+S_i^\alpha \widehat{x})}\right]`$ (A72)
$`=`$ $`\left({\displaystyle \frac{\widehat{q}^C}{C!}}{\displaystyle \underset{l=1}{\overset{C}{}}d\widehat{x}_l\widehat{\pi }(\widehat{x}_l)\underset{\alpha =1}{\overset{n}{}}\left\{e^{\xi F_s}\underset{l=1}{\overset{C}{}}(1+\widehat{x}_l)+e^{\xi F_s}\underset{l=1}{\overset{C}{}}(1\widehat{x}_l)\right\}_\xi }\right)^N,`$ (A73)
and similarly for the term involving the variables $`\tau `$. Substituting these into Eq. (A71), one obtains the $`n`$-th moment of partition function

$`𝒵(𝝃,𝜻,𝒟)^n_{\xi ,\zeta ,𝒟}`$ (A74)
$`=`$ $`\text{extr}_{\{\pi ,\widehat{\pi },\rho ,\widehat{\rho }\}}\{\mathrm{exp}[NC\{{\displaystyle }dxd\widehat{x}\pi (x)\widehat{\pi }(\widehat{x})\mathrm{ln}(1+x\widehat{x})^n1\}`$ (A75)
$`ML\{{\displaystyle }dyd\widehat{y}\rho (y)\widehat{\rho }(\widehat{y})\mathrm{ln}(1+y\widehat{y})^n1\}`$ (A76)
$`+{\displaystyle \frac{1}{2^n}}\{{\displaystyle \frac{NC}{K}}{\displaystyle }\left[{\displaystyle \underset{k=1}{\overset{K}{}}}dx_k\pi (x_k)\right]\left[{\displaystyle \underset{l=1}{\overset{L}{}}}dy_l\rho (y_l)\right]\mathrm{ln}[1+{\displaystyle \underset{k=1}{\overset{K}{}}}x_k{\displaystyle \underset{l=1}{\overset{L}{}}}y_l]^n1\}]`$ (A77)
$`\times \left({\displaystyle \left[\underset{k=1}{\overset{C}{}}d\widehat{x}_k\widehat{\pi }(\widehat{x}_k)\right]\left(\left[e^{F_s\xi }\underset{k=1}{\overset{C}{}}(1+\widehat{x}_k)+e^{F_s\xi }\underset{k=1}{\overset{C}{}}(1\widehat{x}_k)\right]\right)^n_\xi }\right)^N`$ (A78)
$`\times \left({\displaystyle }\left[{\displaystyle \underset{l=1}{\overset{L}{}}}d\widehat{y}_l\widehat{\rho }(\widehat{y}_l)\right]\left([e^{F_n\zeta }{\displaystyle \underset{l=1}{\overset{L}{}}}(1+\widehat{y}_l)+e^{F_n\zeta }{\displaystyle \underset{l=1}{\overset{L}{}}}(1\widehat{y}_l)]\right)^n_\zeta \right)^M\}.`$ (A79)
Finally, in the limit $`n0`$ one obtains

$`{\displaystyle \frac{1}{N}}\mathrm{ln}𝒵(𝝃,𝜻,𝒟)_{\xi ,\zeta ,𝒟}=\underset{n0}{lim}{\displaystyle \frac{𝒵(𝝃,𝜻,𝒟)^n_{\xi ,\zeta ,𝒟}1}{nN}}`$ (A80)
$`=`$ $`\text{extr}_{\{\pi ,\widehat{\pi },\rho ,\widehat{\rho }\}}\{{\displaystyle \frac{C}{K}}\mathrm{ln}2C{\displaystyle }dxd\widehat{x}\pi (x)\widehat{\pi }(\widehat{x})\mathrm{ln}(1+x\widehat{x}){\displaystyle \frac{CL}{K}}{\displaystyle }dyd\widehat{y}\rho (y)\widehat{\rho }(\widehat{y})\mathrm{ln}(1+y\widehat{y})`$ (A81)
$`+{\displaystyle \frac{C}{K}}{\displaystyle \left[\underset{k=1}{\overset{K}{}}dx_k\pi (x_k)\right]\left[\underset{l=1}{\overset{L}{}}dy_l\rho (y_l)\right]\mathrm{ln}\left[1+\underset{k=1}{\overset{K}{}}x_k\underset{l=1}{\overset{L}{}}y_l\right]}`$ (A82)
$`+{\displaystyle \left[\underset{k=1}{\overset{C}{}}d\widehat{x}_k\widehat{\pi }(\widehat{x}_k)\right]\mathrm{ln}\left[e^{F_s\xi }\underset{k=1}{\overset{C}{}}(1+\widehat{x}_k)+e^{F_s\xi }\underset{k=1}{\overset{C}{}}(1\widehat{x}_k)\right]_\xi }`$ (A83)
$`+{\displaystyle \frac{C}{K}}{\displaystyle }\left[{\displaystyle \underset{l=1}{\overset{L}{}}}d\widehat{y}_l\widehat{\rho }(\widehat{y}_l)\right]\mathrm{ln}[e^{F_n\zeta }{\displaystyle \underset{l=1}{\overset{L}{}}}(1+\widehat{y}_l)+e^{F_n\zeta }{\displaystyle \underset{l=1}{\overset{L}{}}}(1\widehat{y}_l)]_\zeta \}.`$ (A84)
## B Evaluation of the Magnetisation

Here, we derive explicitly Eqs.(31) and (32). After using the gauge transformation $`S_i\xi _iS_i`$, the magnetisation can be written as

$$m=\frac{1}{N}\underset{i=1}{\overset{N}{}}\text{sign}(m_i)_{𝝃,𝜻,𝒟},$$
(B1)
introducing the notation $`m_i=S_i_\beta \mathrm{}`$ (gauged average).

For an arbitrary natural number $`p`$, one can compute $`p`$-th moment of $`m_i`$

$$m_{i}^{}{}_{}{}^{p}_{𝝃,𝜻,𝒟}=\underset{n0}{lim}\underset{\beta \mathrm{}}{lim}\underset{\{𝑺^1,𝝉^1\},\mathrm{},\{𝑺^n,𝝉^n\}}{}S_i^1S_i^2\mathrm{}S_i^pe^{\beta _{a=1}^n_a}_{𝝃,\zeta ,𝒟},$$
(B2)
where $`_a`$ denotes the gauged Hamiltonian of the $`a`$-th replica. Decoupling the dynamical variables and introducing auxiliary functions $`\pi ()`$, $`\widehat{\pi }()`$, $`\rho ()`$ and $`\widehat{\rho }()`$, of a similar form to Eq. (A59), one obtains

$$m_{i}^{}{}_{}{}^{p}_{𝝃,𝜻,𝒟}=\underset{l=1}{\overset{C}{}}d\widehat{x}_l\widehat{\pi }(\widehat{x}_l)\mathrm{tanh}^p\left(F_s\xi +\underset{k=1}{\overset{C}{}}\mathrm{tanh}^1\widehat{x}_k\right)_\xi ,$$
(B3)
using the saddle point solution of $`\widehat{\pi }()`$.

Employing the identity

$$\text{sign}(x)=1+2\underset{n\mathrm{}}{lim}\underset{m=0}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{2n}{m}\right)\left(\frac{1+x}{2}\right)^{2nm}\left(\frac{1x}{2}\right)^m$$
(B4)
which holds for any arbitrary real number $`x[1,1]`$ and Eqs.(B3) and (B4) one obtains

$`\text{sign}(m_i)_{𝝃,𝜻,𝒟}`$ $`=`$ $`1+2{\displaystyle 𝑑z\varphi (z)\underset{n\mathrm{}}{lim}\underset{m=0}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{2n}{m}\right)\left(\frac{1+z}{2}\right)^{2nm}\left(\frac{1z}{2}\right)^m}`$ (B5)
$`=`$ $`{\displaystyle 𝑑z\varphi (z)\text{sign}(z)},`$ (B6)
where we introduced a new notation for the distribution

$$\varphi (z)=\underset{l=1}{\overset{C}{}}d\widehat{x}_l\widehat{\pi }(\widehat{x}_l)\delta (zF_s\xi \underset{k=1}{\overset{C}{}}\mathrm{tanh}^1\widehat{x}_k)_\xi ,$$
(B7)
thus reproducing Eqs.(31) and (32).