text:
# New algorithms and technologies for the un-supervised reduction of Optical/IR images

## 1 Introduction

The ESO Imaging Survey (EIS)  project at ESO is an ongoing effort to conduct public imaging surveys and to develop the required infrastructure and software tools to support the large volume of optical/infrared data expected from such surveys, using dedicated imaging telescopes such as VST and VISTA.

Since 1997 EIS  has conducted a variety of medium-size optical/infrared surveys. The challenge has been to cope with the variety of imagers and observing strategies and, more recently, with the large volume of data coming from the wide-field imager (WFI) mounted on the 2.2 m telescope at La Silla. In the early phases of this project the main effort was to adapt existing software to the project’s particular needs leading to the use of different packages (e.g., IRAF, ECLIPSE). This made it difficult to develop a common environment for operations. Therefore, for the past two years an effort has been made to develop a common, end-to-end, integrated system capable of supporting un-supervised reduction. From the EIS survey system it is possible to prepare observations, retrieve raw data and convert these to science grade products in the form of fully calibrated pixel maps, catalogs and other derived products . This effort has progressed in several fronts and includes the development of an administrative wrapper using Python, graphic user interfaces (Tcl/Tk), the integration of a supporting database (Sybase), a suite of survey tools and an image processing package.

Coping with the large volume of data expected in the future, and the multi-wavelength nature of the public imaging surveys currently being conducted requires the image reduction software to have a high-throughput and to be instrument-independent. This has been accomplished using a new, C-based library which integrates XML technology and a host of new algorithms to deal with a variety of situations encountered in observations covering the optical and infrared domains, while using single- and multi-chip instruments. Here a brief review of the main concepts and algorithms is presented.

## 2 Preparation Steps for Reduction

As stated earlier, the software developed is capable of supporting single-CCD or multi-CCD instruments, and a variety of electronic layouts (several over-scan areas and read-out ports). To support an instrument the system requires the following information. First, a list of basic keywords (seven in total) in the FITS header of the original must be translated into those recognized by the system. All others are re-computed as required to conform to the internal convention and ensure consistency throughout the processing phase. Second, information about the camera used such as: the number of CCDs, their geometric layout in the array, and the electronic layout of each CCD. The electronic layout is defined by the location of the *overscan area*s, and the number, location and properties (read-out noise and gain) of different read-out ports. These define different quadrants on the image, which must be known for a proper reduction. Third, the definition of Reduction Modes (RM; see section 4) specific to that instrument.

Once this information is available and the raw data are on disk the process can start. The first step is to check the number of fits extensions in the raw images. If this number corresponds to the number of CCD chips expected for the instrument, the original raw image is split into individual ones. The second step is to translate the keywords in the fits header of the raw images, using the dictionary mentioned above, and to check for their existence and integrity. The third step is to organize the raw images by recognizing their types (science, standards, calibration) and to create reduction blocks (RBs). These consist of groups of images taken sequentially and sharing the same filter and/or position (in the case of science images). Calibration frames are also sorted by their type namely bias, dark, dome flats and twilight flats. Dark exposures with the same integration time are grouped together. These calibration RBs can be further combined according to a period of validity defined by the user. Finally, each RB is assigned a RM which tells the system how these images should be reduced and the parameters to use for each sub-process.

The reduction mode carries information about the instrument (number of chips, chip size, overscan region, read-out ports, gain, readout-noise) and attributes such as the type of image (dark, bias, flat, standard, science), type of survey (extragalactic or galactic fields), and passband (optical/infrared). Each of these attributes may require a special sequence of processes and/or suitable parameters. The information contained in the RB is used to define the sequence of processes and the optimal parameters for each instrument and passband. The present system is capable of sorting out these various pieces of information and use the specific modules required with the appropriate set of parameters. To a large extent this is the major advantage of an integrated system as the one described here.

For science frames, besides the RB, one has also to provide additional information such as the astrometric reference catalog to be used in the registration of the images and the definition of a reference grid onto which the final reduced images will be mapped. This reference grid is defined by a reference sky position, a pixel scale and orientation (e.g., north-up; east-left) and a projection (e.g., TAN, COE). The use of a reference grid is not mandatory but is recommended for deep observations and mosaics since it is easier to co-add images sharing the same grid convention.

In the survey system version, the Python wrapper can use the original groups and apply additional constraints to filter out images. For instance, it is possible to restrict images based on the variation of the seeing (currently based on the values of the seeing monitor), the time or airmass interval covered, and enforce that the reduction block corresponds to the original observation blocks (OB), as defined by ESO’s data flow system.

Before the actual reduction of the data, some additional checks are required. The first check is on the accuracy of the reference pixel stored in the header. This is important in order to detect possible large offsets from the expected location which may occur due to changes in the pointing model at the telescope. This is also critical given the convention adopted by ESO for the CRVALs and CRPIXs values in the FITS header of the raw images (especially SUSI2and WFI). This check is carried out using a cross-correlation technique between two images. The first is a low-resolution (under-sampling by a factor of 8 the original pixel scale of the CCD considered) version of a mock image created from the reference catalog drawn from a much larger area than one CCD (ranging from 4 to 36 times the original field-of-view). The other is constructed from the science image to the same resolution as the one adopted earlier using wavelet transform, padded with zeros to match the mock image size. In preparing the low-resolution science image, care must be taken to remove cosmic rays and the effect that a few bright stars may have on the results. A more detailed description of the method will be presented in a forthcoming paper. The peak of the cross-correlation yields the offset that must be applied to the values of CRPIXs in the header to correspond to the values of the CRVALs. This correction is then applied to all frames (individual chips) belonging to the same night or an arbitrary period of validity specified by the user. This procedure is possible thanks to the speed of the algorithm which works with low-resolution images and benefits from the Fast Fourier Transform (FFT). One advantage of carrying out this reference pixel check at the preparation phase is to assure that during reduction only a refinement of the astrometric solution is required. Note, however, that the method is only valid for correcting translations. When dealing with multi-chip instruments some algorithms (see section 3) require information about the exact location of each chip and this must be checked by the system. Since it is not advisable to rely in the order in which they appear in the FITS file, which may (have) change in time, this must be computed. Once new CRPIX values are available these are used to define the geometry of the array and to cross-check if the computed CRPIXs are consistent.

## 3 Handling Calibration Frames

Normally, the first RBs to be reduced are those associated with the calibration frames. However, it is also possible to assign external master calibrations, in which case this step can be skipped. Each calibration RB of the “closed” type (dark and bias) containing at least three frames a median value is computed and a master bias and dark are created following the validity period specified. Several master darks may exist with different exposure times. In general, darks are only used for IR reduction, since their contribution for modern optical CCDs is negligible. The RB for flats (dome or twilight) contains only frames within a specified range of counts, the others being rejected in the preparation phase. Each frame is then subtracted by the master bias or, in the case of IR, by the corresponding (same exposure time) master darks. Each frame is then divided by their mean to provide normalized images and a median value is computed to produce the master flat. The same procedure is adopted for dome and twilight flats, but currently for test reductions the twilight flat is normally used leading to good results. One exception is the IR images taken with SOFI at La Silla. SOFI uses the so called special-flats, requiring a special reduction procedure which is currently not yet available. Instead the flats prepared on the mountain are used.

Once a master flat is available a bad-pixel (hot/cold/dead pixels) map is created by flagging pixels with intensities greater or smaller than a threshold level on an image storing the difference between the master flat and a smooth version of it produced by computing a median in a 10$`\times `$10 box. This is done automatically and on-the-fly for each reduction with no need for saving this map. At this point the information about the read-out ports is important in order to avoid smoothing over the border of different quadrants. Otherwise these borders will be flagged as bad-pixels.

The master flat also characterizes the noise level of the pixels which can be used to define the notion of weight map (pixels with a low efficiency should have a smaller signal-to-noise). This notion is critical in the co-addition of images and for a homogeneous source extraction at a fixed signal-to-noise level. The weight map is equal to the master flat. It can also be generalized by multiplying it by the bad-pixel map which gives zero weight to the rejected pixels. This weight map is then associated to all raw science frames within the validity period of the master flat. It is important to emphasize that the weight map is equivalent to an effective exposure map.

The procedure described above also applies to multi-CCD instruments, with the system splitting the raw images into individual CCDs. The administration of these several files is carried out internally by the system and it is transparent to the user. However, for CCD arrays an additional step is necessary, namely the correction for variations of the gain from chip-to-chip along the array so as to enforce that all the master flats of the different CCDs have the same gain.

This correction is computed separately and for each filter using long exposures and assuming that all CCDs in the array have the same gain. Hence, after they are corrected for instrumental effects (bias, dark, flat), they should share the same background value. Under this assumption the background is computed at the four borders (left, top, right, bottom) of each CCD in the array. If differences in the counts are found, a flux-scaling factor for each CCD is computed such as to minimize these differences. The flats are then multiplied by the derived flux scale. For the instruments considered (e.g., SUSI2, WFI) these differences are typically less than 10%. It is important to note that this calculation does not take into account nor corrects for illumination effects. The derivation of the illumination correction needs to be done separately and it is not yet fully implemented.

In principle, this procedure prevents the need for the calculation of independent photometric calibrations for each chip. However, the robustness of the procedure must be independently monitored by computing photometric solutions chip by chip. Such a facility is provided by the survey system.

## 4 Handling Science Frames

### 4.1 Reduction

The pipeline offers different ways to process a reduction block of raw images. As mentioned earlier these are set by the *Reduction Mode* (RM) assigned to it. In the present framework four main modules are available as schematically shown in Figure 1. Also indicated at the bottom are the external dependencies (e.g., master calibration images and chip layout, weight map, and reference catalogs for the absolute astrometric solution) of each of the modules. The current system consists of the following steps:

* Removal of instrument signature \- For each optical image the master bias is overscan corrected quadrant by quadrant and subtracted from the science images. Next the science images are divided by the master flat. A similar procedure is adopted for the IR images except that instead of the master bias, the master dark with the appropriate exposure time is used. An option is available to interpolate across bad pixels. In general, however, since multiple exposures are available this is taken into account by the weight map during the co-addition process.
* Removal of local background \- this module is normally used to deal with images observed in passbands where both low and high-frequency features such as fringes are present in the background. This is typically the case for the red passbands in the optical (R and more conspicuously in I) and in the infrared domains. A first crude attempt to estimate the local background is done by assuming that the observed features are stable in pixel coordinates. A first estimate of the local background level is computed as the median of the the intensity of all images in the RB at that pixel. Next all images in the RB are subtracted by this estimated background. However, as shown in sub-section 5.1 this estimate is bias and greatly affects the flux of faint objects. To improve the background estimate it is necessary to mask out objects. In order to do that the images are registered using one of them as reference, and the values are added and divided by the number of images contributing at the position of the pixel in the new grid. Objects are extracted from this image and a mask image is created. This mask is then applied to the original images in the RB at their original pixel coordinates. Pixels inside the masks are then discarded in the calculation of the mean for a second estimate of the background. This second estimate is then used to subtract the background of each image. The resulting RB is then used by the next module. In cases where there are not enough exposures or for very crowded fields the system uses an external fringing map as discussed in sub-section 5.1.
* Removal of low-frequency background \- this module computes and removes the background of each image. This procedure assumes that the images have a smooth background (using a median of a top-hat filter) and that the smoothing scale is larger than the largest object in the image. For each image a mini-background is computed (according to the smoothing scale adopted) and an estimate of the background, for each quadrant, is computed merging a user specified number of mini-backgrounds. Once this estimate is available the mini-backgrounds are re-sampled back to the original image scale and subtracted from all images in the RB. In the case of IR observations carried out in nights of variable conditions, it is recommended to deal with each image separately. At this point any residual imprint of the quadrant structure in the CCD is removed.
* RB stacking & Astrometric calibration \- the final module of the reduction part is responsible for producing a final image which combines all the images in the RB, normalized to one second integration time (and computing an effective gain), by calculating the mean value at each location, and astrometrically calibrating the final image. The noise of the normalized images in the RB are computed and for RBs with at least three images the noise distribution is examined and outliers are removed. In case the individual images have sufficient high S/N objects, the astrometric reference catalog specified in the previous section is used to create a mock image (as discussed above) using the reference grid specified at the preparation level, to determine the absolute astrometric registration for each image which are warped and then combined. This is normally the case used for optical images with exception perhaps of usually shallower $`U`$-band exposures. For the latter and for all IR images a relative astrometric solution is found, images are then warped and combined. This final image is then astrometrically calibrated using the reference catalog mapped (warped) to the final reference grid, exactly in the same way as above. Once the images in the RB are registered (absolutely or relatively) cosmic rays are removed using a sigma-clipping technique. At this point, the Hough transform (see sub-section 5.3) is applied to identify possible satellite tracks. If found, pixels along a strip surrounding the track are rejected. Finally, the weight map is used to mask out bad pixels. This is done before warping the images. In general, this module deals with sky-subtracted images coming out from the previous modules. However, this is not required since in some cases it is of interest to keep the smooth part of the original image background. On the other hand, this module can only be used after fringes have been removed, since they greatly impact the wavelet decomposition, specially in the medium-resolution scales which is used for the first estimate of the astrometric calibration. In the case of images for which the background has not been removed, offsets in the background level, due to variations in sky brightness during the time interval of the RB, are corrected for using one of the images as reference. The final stacked image is trimmed according to a user specified minimum number of contributing images. The weight map of the final image is the sum of the original weights properly normalized and warped, taking into account the exposure times of the original images.

Typical examples of reduction modes are: 1) deep IR images and optical images in I-band all modules are active; 2) in the case of the bluer passbands, the local background module can be bypassed; 3) in the case of crowded fields, such as those in the Pre-Flames survey, only the first and the last modules are used.

### 4.2 Final Products

Often the observations are carried out as a sequence of OBs. Therefore, at the end of the reductions several reduced RBs are available all sharing the same reference grid. At this point these intermediate images can be combined using the weight maps and the noise, as weight. At this stage a further rejection of cosmic rays is applied leading to the removal of lower S/N features not detected in the original frames. In the case of multi-chip instruments at the end of the process the RBs of the individuals CCDs are merged together into the final image as specified by the reference grid. This assumes that flux scaling correction of the flats described earlier is adequate. Otherwise, an external correction factor for each CCD has to be provided and an external module has to be used to combine and merge the different chips.

## 5 Algorithms

### 5.1 Removal of Fringes

As mentioned in the previous section fringes are a common phenomenom present on near-infrared images. Fringes are an additive signal stemming from sky emission lines, making a non-negligible contribution to the total signal, typically of the order of 10% of the background level. Hence, subtracting its contribution is essential for a proper reduction of the image. While in the optical, the fringing pattern is relatively stable (but exceptions exist) over the time scale of a night, in the infrared the fringes vary very quickly. Fortunately, the fringing pattern is largely independent of the pointing and consecutive images share the same pattern. This allows them to be separated from the astronomical objects by shifting the pointing observations. While the objects will fall in different pixel locations the fringes will not.

Therefore, a first crude estimate of the fringing pattern is to compute the median at every pixel using the estimator $`E`$ given by

$$B_{i,j}=E(P_{n_{i,j}}\mu _n,n)$$
(1)
where $`B_{i,j}`$ is the computed background at the location $`(i,j)`$, $`P_{n_{i,j}}`$ is the counts of the pixel $`(i,j)`$ of the $`n^{th}`$ image, and, $`\mu _n`$ is the mean of the $`n^{th}`$ image. Subtraction of $`\mu _n`$ is important because it may vary from one image to another in the case of non-photometric nights. The estimator apply on the RB is shown in Figure 2 on the left panel.

Unfortunately, this estimator is biased due to the presence of faint astronomical objects which cannot be detected in the individual images due to their low S/N. However, if all images are combined these objects can have a significant contribution. For instance, for an RB with 60 images objects detected in the final image with a S/N $``$ 2 have a S/N $``$ 0.25 on the individual images. Neglecting their contribution in estimating the background leads to an overestimate of the background, and consequently to an underestimate of the flux of the objects. This effect can be as large as 0.2 mag near the limiting magnitude.

In order to improve the estimator $`E`$ in equation (1) it is necessary to distinguish pixels which have the “*background value*” from those contaminated by the light of objects as illustrated in the right side of Figure 2. In order to do that one has to use the first estimate of the background, combine the sky-subtracted images, detect objects, and create masks which encompass surrounding pixels containing the object. These masks are then mapped back to the individual images and used to discard pixels in a second estimate of the background. Typically, the mask size is taken to be about twice the estimated size of the objects. The estimator apply on the RB with the mask is shown in Figure 2 on the right pannel.

If the background is evolving during the total exposure time, it can be estimated only from some images observed consecutively. The algorithm allows to specify a range of images to be used to estimated the background. In the IR the background is evolving faster, and a range of typically 15 images is used to compute the background, corresponding to 15 minutes of observation.

One limitation of the method is that, for large objects and very dense regions it may be impossible to “see” the sky, in which case holes will appear in the final image. In order to avoid this case a proper choice of the amplitude of the jitter pattern is critical.

### 5.2 Astrometric registration using wavelet transform

In the pipeline astrometric solutions are computed either in a relative or in an absolute way. For instance, relative astrometry is done for the stacking of an RB using one of the images as reference, while an absolute astrometric calibration requires an astrometric reference catalogue and a grid. In both cases the solution is found using the concept of multi-resolution (reference), by comparing low-resolution images obtained from a wavelet decomposition. In the case of a reference catalog this is done by first creating a mock image to map the objects listed in the catalog onto a reference image according to the grid convention. This method is referred to as pixel-based image registration .

The calculation of the geometrical distortion between the two images requires: 1) the identification of a set of well defined astronomical objects (control points, hereafter CP) in both images; 2) the definition of a functional form for the deformation model which maps a point $`(x,y)`$ in the input image to a point $`(X,Y)`$ on the reference image; 3) the creation of the distorted or corrected image, which requires the adoption of an interpolation or warping method.

The distortion model adopted by the pipeline is a pair of bivariate polynomials given by:

$$\{\begin{array}{ccc}X_i& =\underset{p=0}{\overset{N}{}}\underset{q=0}{\overset{Np}{}}a_{pq}x_i^py_j^q=& P(x_i,y_i)\\ Y_i& =\underset{p=0}{\overset{N}{}}\underset{q=0}{\overset{Np}{}}b_{pq}x_i^py_i^q=& Q(x_i,y_i)\end{array}$$
(2)
Where $`\text{(}X_\text{i},Y_i)`$ are the coordinates of the $`i^{th}`$ CP in the reference image and $`(x_i,y_i)`$ the corresponding CP from the input image. $`N`$ is the degree of the polynomial. Usually, the degree of the polynomial is one or two for instruments with a small field-of-view and up to four otherwise. Then follows the computation of the unknown parameters ($`\frac{(N+1)(N+2)}{2}`$ for each polynomial) using the least mean square estimator (LMSE).

The most difficult step of the process is to associate the CP from the input image with the CP from the reference image. CP are randomly spread and can be very dense. The associations cannot be done using any *window search* because multiple associations will occure and will bias the distortion estimation. The idea of the wavelet transform, is to work with different resolutions and therefore give a constraint for the minimum distance between CP.

Let $`I_r`$ be the reference image and $`I_i`$ be the input image. On each image the wavelet transform with the so called *algorithm à trous*  is performed up to the scale $`L`$ producing two wavelet images $`W_{il}`$ and $`W_{rl}`$ with: $`l(1,L)`$ from the two images $`I_i`$ and $`I_r`$. $`L`$ being the initial dyadic step, it corresponds to the large details of the images. A detection procedure on $`W_{iL}`$ is achieved and only the local maxima with amplitude greater than three times the sigma of the noise are kept. Theses local maxima play the role of CP. They correspond to significant image patterns, and have to be found in the corresponding wavelet plan of the reference image $`W_{rL}`$. The association is done with a window search. The search radius is set to $`2^{L1}`$ pixels, and corresponds to the half of the maximum wavelength (Shannon constraint) appearing at the scale $`L`$. This radius assures no multiple associations. At this scale identical CP can be matched with confidence, and therefore the relationship between the coordinates of the different images can be determined. This provides a rough estimation of the distortion model. Considering now the wavelet images of order $`L1`$, a new set of maxima in each wavelet plan is detected. The coordinates of each maximum detected are then transformed in the reference image using the previous estimation of the distortion model. This follows an iterative process until the last scale corresponds to the wavelet plan of the best resolution, establishing thus the best goemetrical correspondance. A first degree polynomial is used for the scale $`L`$, but the degree can increased for the last scales.

For each input pixel at location $`(i,j)`$ is computed $`(k,l)`$, where $`k=P(i,j)`$ and $`l=Q(i,j)`$. The coordinates of the distorted pixel $`(k,l)`$ are generally not integers, so an interpolation must be done to calculate the intensity value for the output pixel. Nearest-neighbor, bi-linear, bi-cubic, splines, and Lancoz interpolations are the most widely used . Currently, only the nearest-neighbor and Lancoz are used with the latter providing the best results.

The advantage of the multi-resolution approach is that it is insensitive to the presence of spurious objects in the reference catalog (stars with large proper motion, objects that do not appear on the observed image due to their color). This multi-resolution approach has also proven to be very efficient in very crowded areas such as the Small Magellanic Cloud. Exhaustive tests have also been performed on images of the open cluster Messier 67 observed with wide-field imager .

### 5.3 Satellite track detection

A frequent problem encountered in imaging surveys using wide-field imagers is the presence of a significant number of satellite tracks. This is a non-trivial problem for survey work since these features, if not identified in the original exposures, will appear in the final co-added image. The presence of these tracks are not only unpleasant from the cosmetic point of view but also have a considerable impact on the derived products (source extraction). Once propagated to the final image, none of the available options are adequate. One can either mask the final image, which may continue to affect the source extraction, or to go back to the individual exposures, mask the visible tracks, and stack the images again. While this manual process can be done for a few images it is certainly not feasible for large data flows given the high incidence of these features. Therefore, a robust system must be able to identify and remove these unwanted features automatically during the processing phase which precedes the stacking .

Fortunately, this problem can be tackled using the Fast Hough Transform (FHT)  since a track consists of a set of points aligned in an otherwise quasi-random distribution of points. The method is based on the fact that the FHT transforms points to the original space into lines (using Bresenham’s procedure) in the Hough space. In this space the intersection of the lines correspond to the parameters of a straight line in the original space.

The Hough space is a two-dimensional space in which each point $`x,y`$ in the image space is transformed into a line given by $`y=xX+y`$. The amplitude in the Hough space is the sum of the number of points contributing at each position in this space. The alignement of $`N`$ points in real space following the relatiion $`y_i=ax_i+b`$, where $`i=1,N`$, will imply a peak in the intensity of the amplitude in Hough space at the point ($`X=a,Y=b`$) given by $`b=x_ia+y_i`$. This transformation is illustrated in Figure 4 where two tracks in real space are shown in the left panel and the Hough transform of all points along these two lines are represented in the Hough space in the right panel. As discussed above the intersection of the point-transformed leads to an intersection in the Hough space yielding a local maximum. The identification of these high-amplitude peaks yields the coefficients $`a`$ and $`b`$ that correspond to the slope and the offset of the track in real space.

Note, however, that the case shown is valid only for slopes in the interval $`[1,+1]`$. In order to deal with slopes outside this range and limit the size of the FHT space, the original catalog is flipped to satisfy the same conditions. Therefore, the actual implementation of the algorithm consists of two steps.

Using the method described above a catalog is extracted from single exposures considering significant local maxima. This catalog is FHTed, and a high-pass filter is used to remove the constant response of FHT. Finally, possible tracks are identified on the Hough filtered image applying a fixed threshold given by $`\frac{\sqrt{N}}{k}`$, where $`N`$ is the number of points in catalog, and, $`k`$ is the sensitivity factor. Tests show that $`k=2`$ is adequate for very faint detections without false positives. To illustrate the method Figure 5 shows a portion of an image with a faint satellite track (left panel), where the dashed line represents the line detected by the procedure on the Hough filtered image (right panel). This methodology has now been implemented in the EIS pipeline and it is routinely used with very good results.

## 6 General Comments about the code and Performance

The software package presented here is an evolution of the so called multi-resolution visual model or MVM, which consists of basic functions for operations on images, wavelet transform and image registration, complemented by open source libraries such as CFITSIO, FFTW and LIBXML2, among others. The MVM library is written in C++ and runs on different platforms: SUN-OS, TRUE-64 on an ALPHA computer, and PC-Linux. Using these basic blocks the package described in this paper has been developed over the past two years to deal with the specific needs of astronomical images and consists of 90.000 lines of code. The various processes are bound together by employing the Extensible Mark-up Language (XML) technology which is used to input parameters, as a way of communicating between different modules, as an internal database and to store the history of the process. The latter can be much more comprehensive than it would be possible in the FITS header of the image. Furthermore, the associated XML files can be translated to any format and make their content available in a human-readable form using suitable style-sheets (XSLT protocol), which can be processed to produce and display HTML pages. An example is shown in figure 6 where the parameters characterizing a particular reduction mode are shown. The figure illustrates the use of a particular style-sheet which produces a concise description for internal use. A more comprehensive style-sheet is used for other users, which gives more detailed information about each parameter.

By using a common library approach, instead of a concatenation of stand-alone programs, the package avoids unnecessary access to disk and the entire reduction is carried out in memory, thereby enhancing significantly the data rate of the process. This is paramount for dealing with large data volumes, and is true both for the stand-alone implementation and the plug-in version for the survey system. In the latter case, the Python wrapper is only responsible for the preparation and the administration of the RBs, thereafter launching the reduction process. The communication between Python and the library is only done at the beginning and at the end of the process via the XML configuration and log files produced.

The stand-alone version of the package has been extensively used to reduce images from a variety of past and ongoing public surveys allowing thus a fairly robust estimate of its performance. Some examples are given in Table 1, where the survey name, type, instrument, detector size, number of images and total time required for the reduction are listed.

The performance refers to an Alpha DS20E running at 666 MHz with two CPUs. The values listed above yield a typical reduction data-rate of $`0.5`$ Mega-pixels per second for optical images, and $`0.4`$ Mega-pixels per second for IR images (which requires the removal of fringes). The overhead time for the preparation and the reduction of the calibration images depends on the instrument, and it is typically $`50`$ percent of the reduction time. The pipeline has also been tested successfully on other multi-CCD instruments such as: the CFHT 12k camera, and the SUSI2 instrument mounted on the ESO-NTT telescope.

## 7 Conclusion

In this paper a new image reduction system, incorporating new technologies and algorithms is presented. It is part of the effort of the EIS team at ESO to develop an integrated, end-to-end survey system capable of converting raw images to science-grade products from multi-wavelength, multi-instrument public surveys. The system incorporates a number of procedures design to cope with a variety of particular situations encountered in real observations covering the optical/infrared domain. The system has a high-throughput ($``$ 0.5 Mega-pixels/sec) and handles the entire suite of imagers available at ESO. It also makes extensive use of XML technology to interface with the user and to provide the full history of the process. The code is written in C++ and it is platform independent. The development of new algorithms is still ongoing (e.g., illumination correction) as well as work aiming at the full integration of the library into the survey system. Efforts are also made to prepare the required documentation and user interface to make the stand-alone version of the image reduction pipeline publicly available.

###### Acknowledgements.

I would like to express my thanks to A. Bijoui for useful discussions over the years and to the EIS project and team members for their support and contribution during the development of this work.


format: text

text:
# Replica bounds for optimization problems and diluted spin systems

## I Introduction

The replica method , originally devised as a trick to compute thermodynamical quantities of physical systems in presence of quenched disorder, has found applications in the analysis of systems of very different nature, as Neural Networks, Combinatorial optimization problems , Error Correction Codes  etc.

Although many physicists believe that the method, within the Replica Symmetry Breaking scheme of Parisi , is able to potentially give the exact solution of any problem treatable as a mean field theory, the necessary mathematical foundation of the theory is still lacking, after more then 20 years from its introduction in theoretical physics. The last times have seen a growing interest of the mathematical community in the method, leading to important but still partial results, confirming in certain cases the replica analysis, with more conventional and well established techniques . Apart the remarkable exception of the analysis of the fully connected $`p`$-spin model in ref.  and the rigorous analysis of Random Energy Models , the analysis of the mathematicians has been, as far as we know, restricted to the high temperature regions and/or to problem of replica symmetric nature.

Very welcomed have been the techniques recently introduced by Guerra and Toninelli  which allow rigorous analysis not relying on the assumption of high temperature, and valid even in problems with replica symmetry breaking. Along these lines, an important step towards the rigorous comprehension of the replica method, has been undertaken in , where it has been shown how in the case of the Sherrington-Kirkpatrick model, and its $`p`$-spin generalizations, the replica free-energies with arbitrary number of replica symmetry breaking steps constitute variational lower bounds to the true free-energy of the model. As stated in that paper, the analysis is restricted to fully-connected models, whose replica mean field theory can be formulated in terms of a single $`n\times n`$ matrix. However, in recent times, many of the more interesting problems analyzed with replica theory pertain to the so called “diluted models” where each degree of freedom interacts with a finite number of neighbors. The introduction of a “population dynamics algorithm”  has allowed to treat in full generality -within statistical precision- complicated sets of probabilistic functional equations appearing in the one step symmetry broken framework of diluted models. The same algorithm has been used as a starting point of a generalized “belief propagation” algorithm for optimization problems . Furthermore, at the analytic level, simplifications due to graph homogeneities in some cases , and to the vanishing temperature limit in some other cases  have led to supposedly exact solutions of the ground state properties of diluted models, culminated in the resolution of the random XOR-SAT on uniform graphs in  and the random K-SAT problem in  within the framework of “one-step replica symmetry breaking” (1RSB).

The aim of this paper, is to show that the replica analysis of diluted models provides lower bounds for the exact free-energy density, and ground state energy density. We analyze in detail the cases of the diluted $`p`$-spin model on the Poissonian degree hyper-graphs also known as random XOR-SAT problem and the random K-SAT problems. We expect that along similar lines free-energy lower bounds can be found for many other diluted cases.

The Guerra method we use sheds some light on the meaning of the replica mean field theory. The physical idea behind the method is that within mean field theory one can modify the original Hamiltonian weakening the strength of the interaction couplings or removing them partially or totally, and compensate this removal by some auxiliary external fields. In disordered systems these fields should be random fields, taken from appropriate probability distributions and possibly correlated with the original values of the quenched variables eliminated from the systems. One is then led to consider Hamiltonians interpolating between the original model and a pure paramagnet in a random field, and by means of these models achieving free-energy lower bounds. We will see that the RS case corresponds to assuming independence between the random fields and the quenched disorder. The Parisi RSB scheme, assumes at each breaking level a peculiar kind of correlations, and gives free-energy bounds improving the RS one.

Our paper is organized in this way: in section 2 we introduce some notations that will be extensively used in the following sections. In section 3 we introduce the general strategy to get the replica bounds; we then specialize to the replica symmetric and the one step replica symmetry broken bounds, giving the results in the $`p`$-spin and the $`K`$-SAT cases. Conclusions are drawn in section 4. In the appendices some details of the calculations in both the $`p`$-spin and the $`K`$-SAT cases are shown.

Our results will be issue of explicit calculations. Although at the end we will get bounds, formalizable as mathematical theorems, the style and most of the notations of the paper will be the ones of theoretical physics.

## II Notations

The spin models we will consider in this work are defined by a collection of $`N`$ Ising $`\pm 1`$ spins $`𝐒=\{S_1,\mathrm{},S_N\}`$, interacting through Hamiltonians of the kind

$$^{(\alpha )}(𝐒,𝐉)=\underset{\mu =1}{\overset{M}{}}H_{J^{(\mu )}}(S_{i_1^\mu },\mathrm{},S_{i_p^\mu })$$
(1)
where the indices $`i_l^\mu `$ are i.i.d. quenched random variables chosen uniformly in $`\{1,\mathrm{},N\}`$. We will call each term $`H_{J^{(\mu )}}`$ a clause. The subscript $`J^{(\mu )}`$ in the clauses indicates the dependence on a single or a set of quenched random variables, as it will be soon clear. The number of clauses $`M`$ will be taken to be proportional to $`N`$. For convenience we will choose it to be for each sample a Poissonian number with distribution $`\pi (M,\alpha N)=\mathrm{e}^{\alpha N}\frac{(\alpha N)^M}{M!}`$. The fluctuations of $`M`$ will not affect the free-energy in the thermodynamic limit, and this choice, which slightly simplify the analysis, will be equivalent to choosing a fixed value of $`M`$ equal to $`\alpha N`$. The clauses themselves will be random. The $`p`$-spin model has clauses of the form

$$H_{J^{(\mu )}}(S_{i_1^\mu },\mathrm{},S_{i_p^\mu })=J^\mu S_{i_1^\mu }\mathrm{}S_{i_p^\mu }.$$
(2)
This form reduces to $`H_{J^{(\mu )}}(S_{i_1^\mu },S_{i_2^\mu })=J^\mu S_{i_1^\mu }S_{i_2^\mu }`$ in the case of the Viana-Bray spin glass $`p=2`$. In both cases the $`J^\mu `$ will be taken as i.i.d. random variable with regular symmetric distribution $`\mu (J)=\mu (J)`$. Notice that for $`\mu (J)=1/2[\delta (J+1)+\delta (J1)]`$ the model reduces to the random XOR-SAT problem  of computer science. The random K-SAT clauses have the form

$$H_{J^{(\mu )}}(S_{i_1^\mu },\mathrm{},S_{i_p^\mu })=\underset{l=1}{\overset{p}{}}\frac{1+J_{i_l^\mu }^\mu S_{i_l^\mu }}{2},$$
(3)
where the $`J_{i_l^\mu }^\mu =\pm 1`$ are i.i.d. with symmetric probability. (The number $`p`$ of spin appearing in a clause is usually called $`K`$ in the K-SAT problem, for uniformity of notation we will deviate from this convention). Notice that in all cases, on average each spin participate to $`\alpha =\frac{M}{N}`$ clauses, and that the set of spins and interactions defines a random diluted hyper-graph of uniform rank $`p`$ and random local degree with Poissonian statistics in the thermodynamic limit. At high enough temperature, the existence of the free-energy in the thermodynamic limit for models of this kind has been proved in by Talagrand in , together with the validity of the RS solution. A proof valid at all temperature based on the ideas presented in this paper, can be obtained for even $`p`$ in analogy of the analysis in  for long range models. We sketch it in appendix C in the case of the $`p`$-spin model.

In establishing the free-energy bounds we will need several kind of averages:

* The Boltzmann-Gibbs average for fixed quenched disorder: given an observable $`A(𝐒)`$

$$\omega (A)=\frac{_𝐒A(𝐒)\mathrm{exp}(\beta (𝐒,𝐉))}{Z}$$
(4)
where $`Z=_𝐒\mathrm{exp}(\beta (𝐒,𝐉))`$ and $`\beta `$ is the inverse temperature.

Obviously, $`\omega (A)`$, as well as $`Z`$ will be functions of the quenched variables, the size of the system and the temperature. This dependence will be made explicit only when needed.
* The disorder average: given an observable quantity $`B`$ dependent on the quenched variables appearing in the Hamiltonian, we will denote as $`E(B)`$ its average. This will include the average with respect to the $`J`$ variables and the choice of the random indices in the clauses as well as with respect to other quenched variables to be introduced later.
* We will need in several occasion the “replica measure”

$$\mathrm{\Omega }(A_1,\mathrm{},A_n)=E(\omega (A_1)\mathrm{}\omega (A_n))$$
(5)
and some generalizations that we will specify later.
* We will occasionally use other kinds of averages, as well as other notations, for which we will use an angular bracket notation, with a subscript indicating the variable(s) over which the average is performed. e.g. an average over a random variable $`u`$ with probability distribution $`Q(u)`$ will be denoted equivalently as $`𝑑uQ(u)()𝑑Q(u)()_u`$. Analogously, averages over distribution families of $`Q(u)`$ will be denoted as $`𝑑Q𝒬(Q)()𝒟𝒬(Q)()_Q`$. Subscripts will be omitted whenever confusion is not possible.
* Another notation we will have the occasion to use in the one for the overlaps among $`l`$ spin configurations $`\{S_i^{a_1},\mathrm{},S_i^{a_l}\}`$, out of a population of $`n`$ $`\{S_i^1,\mathrm{},S_i^n\}`$:

$$q^{(a_1,\mathrm{},a_l)}=\frac{1}{N}\underset{i=1}{\overset{N}{}}S_i^{a_1}\mathrm{}S_i^{a_l}(1a_rnr),$$
(6)
and in particular

$$q^{(n)}=q^{(1,\mathrm{},n)}=\frac{1}{N}\underset{i=1}{\overset{N}{}}S_i^1\mathrm{}S_i^n,$$
(7)
This notation will be extended to multi-overlaps in the 1RSB case, as we will specify in section III B.

In the following we will need to consider averages where some of the variables are excluded, e.g. the averages when a variable $`u_i^{k_i}`$ is erased. These average will be denoted with a subscript $`u_i^{k_i}`$ e.g. if an $`\omega `$ average is concerned the notation will be $`\omega ()_{u_i^{k_i}}`$. Other notations will be defined later in the text whenever needed.

Our interest will be confined to bounds to the free-energy density $`F_N=\frac{1}{\beta N}E\mathrm{log}Z`$ and the ground state energy density $`U_{GS}=lim_N\mathrm{}1/NE\left[min(U_N)\right]`$ valid in the thermodynamic limit, so that $`O(1/N)`$ will be often implicitly neglected in our calculations.

## III The general strategy

The strategy to get the replica bound is a generalization of the one introduced by Guerra in the case of fully connected models . We will consider models which will interpolate between the original ones we want to analyze and pure paramagnet in random fields with suitably chosen distribution. The underlying idea is that, given the mean field nature of the models involved, if one was able to reconstruct the real local fields acting on a given spin variable via a given hyper-edge, and to introduce auxiliary fields acting on that variable in such a way to energetically balance the deletion of the hyper-edge, then it would be possible to have an exact expression for the free-energy in terms of such auxiliary fields even when the whole edge set was emptied. However, if the replacement is done with some approximate form of the auxiliary fields distribution function, the real free-energy will be the one calculated using the approximate fields plus an excess term at every step of the graph deletion process. The proof of the definite sign of this excess term gives a way to determine bounds for the thermodynamic quantities.

We will prove the existence of replica lower bounds to the free-energy density of the p-spin model and the random K-SAT problem. In this last case our result proves that the recent replica solution of  gives a lower bound to the ground state energy and therefore an upper bound for the satisfiability threshold. The proofs will strictly hold in the $`N\mathrm{}`$ limit, due to the presence of corrections of order $`1/N`$ in the calculated expressions for any finite size graph. Moreover, our proofs will be restricted to the $`p`$-spin model the the K-SAT with even $`p`$. In the cases of odd $`p`$ the same bound would hold if one could rely on some physically reasonable assumptions on the overlap distribution (see below).

Our analysis will start from the TAP equations for the models , and their probabilistic solutions implied by the cavity, or equivalently the replica method at various degrees of approximation. We will consider in particular the replica symmetric (RS) and one step replica symmetry broken solutions, but it should be clear from our analysis how to generalize to more steps of replica symmetry breaking. In the TAP/cavity equations one singles out the contribution of the clauses and the sites to the free-energy and defines cavity fields $`h_i^{(\mu )}`$ and $`u_\mu ^{(i)}`$ respectively as the local field acting on the spin $`i`$ in absence of the clause $`\mu `$ and the local field acting on $`i`$ due to the presence of the clause $`\mu `$ only. If we define $`Z_N[S_i]`$ as the partition function of a given sample with $`N`$ spins where all but the spin $`i`$ are integrated, $`F_{N,i}`$ the free-energy of the corresponding systems where the spin $`S_i`$ and all the clauses to which it belongs are removed, we can write,

$`Z_N[S_i]`$ $`=`$ $`\mathrm{e}^{\beta F_{N,i}}{\displaystyle \underset{\mu T_i}{}}{\displaystyle \underset{S_{i_2^\mu },\mathrm{},S_{i_p^\mu }}{}}\mathrm{e}^{\beta H_{J^{(\mu )}}(S_{i^\mu },S_{i_2^\mu },\mathrm{},S_{i_p^\mu })+_{l=2}^ph_{i_l^\mu }^{(\mu )}S_{i_l^\mu }}`$ (8)
$`=`$ $`\mathrm{e}^{\beta F_{N,i}}{\displaystyle \underset{\mu T_i}{}}B_\mu ^{(i)}\mathrm{e}^{\beta u_\mu ^{(i)}S_i}`$ (9)
where $`T_i`$ is the set of clauses containing the spin $`i`$, and the constant $`B_\mu ^{(i)}=\mathrm{e}^{\beta \mathrm{\Delta }F_\mu ^{(i)}}`$ can be interpreted as suitable shifts in the free-energy due to the contribution of the clause $`\mu `$ for fixed value of the spin $`i`$. We notice that denoting $`J^\mu `$ as $`J`$, and renaming the fields in (9) into $`h_1,\mathrm{},h_{p1}`$, Eq. (9) defines functions

$$u_J(h_1,\mathrm{},h_{p1})\mathrm{and}B_J(h_1,\mathrm{},h_{p1}).$$
(10)
The equation are closed by the self-consistent condition:

$$h_i^{(\mu )}=\underset{\nu \{T_i\mu \}}{}u_\nu ^{(i)}$$
(11)
These equations are at the basis of iterative algorithms such as the “belief propagation” or “sum-product” know for a long time in statistical inference  and coding theory  and the more recently proposed algorithm of “survey propagation”. Conditions (9) and (11) can be diagrammatically represented as in fig.(1).

The cavity fields solutions of (9,11) are random variables which fluctuate for two reasons . First, they differ from sample to sample. Second, within the same sample the equations can have several solutions which can level-cross. The cavity/replica method provides under certain assumption probabilistic solutions. In the RS approximation, one just supposes a single solution to give the relevant contribution in a given sample. The sample to sample fluctuation induce probability distributions $`P(h)`$ and $`Q(u)`$ whose relations implied by (9,11) are:

$`P(h)`$ $`=`$ $`{\displaystyle \underset{k}{}}\mathrm{e}^{\alpha p}{\displaystyle \frac{(\alpha p)^k}{k!}}{\displaystyle 𝑑u_1Q(u_1)\mathrm{}𝑑u_kQ(u_k)\delta (h\underset{i=1}{\overset{k}{}}u_k)}`$ (12)
$`Q(u)`$ $`=`$ $`{\displaystyle 𝑑h_1P(h_1)\mathrm{}𝑑h_{p1}P(h_{p1})\delta (uu_J(h_1,\mathrm{},h_{p1}))_J}`$ (13)
where $`_J`$ denotes the average over the random variables appearing in a clause. In addition to sample to sample fluctuations, the 1RSB solution assumes fluctuations of the fields from solution to solution of the equations, so that the functions $`P(h)`$ and $`Q(u)`$ will be themselves randomly distributed according to some functional probability distributions $`𝒫(P)`$ and $`𝒬(Q)`$ related by the self-consistency equations

$$𝒬(Q)=𝒟P_1𝒫(P_1)\mathrm{}𝒟P_{p1}𝒫(P_{p1})\delta (Q()Q(|P_1,\mathrm{},P_{p1},J))_J$$
(14)
$$𝒫(P)=\underset{k=0}{\overset{\mathrm{}}{}}\mathrm{e}^{\alpha p}\frac{(\alpha p)^k}{k!}\underset{l=1}{\overset{k}{}}𝒟Q_l𝒬(Q_l)\delta (P()P(|Q_1,\mathrm{},Q_k))$$
(15)
where:

$`Q(u|P_1,\mathrm{},P_{p1},H)=𝒩_P[P_1,\mathrm{},P_{p1}]{\displaystyle 𝑑h_1P_1(h_1)\mathrm{}𝑑h_{p1}P_1(h_{p1})B_𝐉(h_1,\mathrm{},h_{p1})^m\delta (uu_𝐉(h_1,\mathrm{},h_{p1}))}`$ (16)
$`P(h|Q_1,\mathrm{},Q_k)=𝒩_{Q,k}[Q_1,\mathrm{},Q_k](2\mathrm{cosh}(\beta h))^m{\displaystyle \underset{l=1}{\overset{k}{}}du_l\frac{Q_l(u_l)}{(2\mathrm{cosh}(\beta u_l))^m}\delta (h\underset{l=1}{\overset{k}{}}u_l)}`$ (17)
where $`𝒩_{Q,k}[Q_1,\mathrm{},Q_k]`$ and $`𝒩_G[G_1,\mathrm{},G_{p1}]`$ insure normalization and $`B_𝐉(g_1,\mathrm{},g_{p1})`$ is a rescaling term of the form (10) that can be reabsorbed in the normalization in the case of the $`p`$-spin model. Its form for the $`K`$-SAT case is given in the appendix. $`m`$ is a number in the interval $`(0,1]`$, which within the formalism selects families of solutions at different free-energy levels. The physical free-energy is estimated maximizing over $`m`$.

The interpretation of these equations has been discussed many times in the literature . We will show here, that such choices in the field distributions result in lower bounds for the free-energy analogous to the ones first proved by Guerra in fully connected models. In order to prove these bounds, we will have to consider auxiliary models where the number of clauses $`\alpha N`$ will be reduced to $`\alpha tN`$ ($`0t1`$), while this reduction will be compensated in average by some external field terms of the kind:

$$_{ext}^{(t)}=\underset{i}{}\underset{l_i=1}{\overset{k_i}{}}u_i^{l_i}S_i$$
(18)
where the numbers $`k_i`$ will be i.i.d. Poissonian variables with average $`\alpha p(1t)`$. As the notation suggests, the fields $`u_i^l`$ will play the role of the cavity fields $`u_\mu ^{(i)}`$ of the TAP approach, and they will be i.i.d. random variables with suitable distribution. Indeed, for each field $`u_i^{l_i}`$ we will chose in an independent way $`p1`$ primary fields $`g_i^{l_i,n}`$ ($`n=1,\mathrm{},p1`$) and clause variables $`J_i^{l_i,n}`$ such that the relation

$$u_i^{l_i}=u_{J_i^{l_i,n}}(g_i^{l_i,1}\mathrm{}g_i^{l_i,p1})$$
(19)
is verified. Notice that the compound Hamiltonian

$$_{tot}^{(t)}[𝐒]=^{(\alpha t)}[𝐒]+_{ext}^{(t)}[𝐒]$$
(20)
will constitute a sample with the original distribution for $`t=1`$, while it will consist in a system of non interacting spins for $`t=0`$. The key step of the procedure, consists in the choice of the distribution of the primary fields $`g_i^{l_i}`$. We will also find useful to define fields $`h_i`$ verifying

$$h_i=\underset{l=1}{\overset{k_i}{}}u_i^l.$$
(21)
The field $`u`$ are related to the $`g`$’s by a relation similar to (9), while the $`h`$’s are related to the $`u`$’s by a relation similar to (11). Of course, the statistics of the fields $`h`$ and the $`g`$’s do coincide in the TAP approach. It is interesting to note that the bounds we will get, are optimized precisely when their statistical ensemble coincide. As we mentioned, various Replica bounds are obtained assuming for the fields $`g_i^{l_i}`$ the type of statistics implied by the different replica solution. So, the Replica Symmetric bound is got just supposing the field as quenched variables completely independent of the quenched disorder and with distribution $`G(g)`$. For the one-step RSB bound on the other hand the distribution $`G`$ will itself be considered as random, subject to a functional probability distribution $`𝒢[G]`$. More complicated RSB estimates, not considered in this paper, can be obtained along the same lines. The case of the fully connected models considered by Guerra can be formalized in this way where the various field distributions involved are Gaussian.

### A The RS bound

We consider in this case i.i.d. fields $`u`$ and $`h`$ distributed according probabilities $`Q(u)`$ and $`P(h)`$ verifying the following relation with the distribution $`Q(g)`$ of the primary fields.

$`Q(u)`$ $`=`$ $`{\displaystyle 𝑑g_1G(g_1)\mathrm{}𝑑g_{p1}G(g_{p1})\delta (uu_J(g_1,\mathrm{},g_{p1}))_J}`$ (22)
$`P(h)`$ $`=`$ $`P(h|k)\pi (k,\alpha p(1t))`$ (23)
$`P(h|k)`$ $`=`$ $`{\displaystyle 𝑑u_1Q(u_1)\mathrm{}𝑑u_kQ(u_k)\delta (h\underset{i=1}{\overset{k}{}}u_k)}`$ (24)
The distribution $`G(g)`$ will be chosen to be symmetric under change of sign of $`g`$, and regular enough for all the expression below to make sense. The RS bound can now be obtained following a procedure to the one of Guerra for the SK model, and considering the $`t`$ dependent free-energy; with obvious notation:

$$F(t)=\underset{N\mathrm{}}{lim}F_N(t)=\underset{N\mathrm{}}{lim}\frac{1}{\beta N}E\mathrm{log}Z_N(t)$$
(25)
where $`E`$ represents the average over all the quenched variables, the one defining the clauses and the external fields. We then consider the $`t`$ derivative of $`F_N`$

$$\frac{d}{dt}F_N(t)=\frac{1}{N\beta }\frac{d}{dt}E(\mathrm{log}Z_N).$$
(26)
As in  we will then write

$$F(1)=F(0)+_0^1𝑑t\frac{d}{dt}F(t)$$
(27)
and show, by an explicit computation, that, up to $`O(1/N)`$ terms that will be systematically neglected, the expression coincides with the variational RS free-energy plus a remainder. In fortunate cases this term will have negative sign and neglecting it will immediately result in a lower bound for the free-energy. This happens in the Viana-Bray model, the $`p`$-spin and the K-SAT for even $`p`$. In the cases of odd $`p`$ we were not able to prove the sign definiteness of the remainder, although as we will discuss we believe this to be the case on a physical basis.

The time derivative of $`F`$ take contributions from the derivative of the distribution of the number of clauses $`M`$

$$\frac{d\pi (M,\alpha tN)}{dt}=N\alpha (\pi (M,\alpha tN)\pi (M1,\alpha tN))$$
(28)
and the distribution of the number of $`u`$ fields on each site

$$\frac{d\pi (k_i,\alpha p(1t))}{dt}=\alpha p(\pi (k_i,\alpha p(1t))\pi (k_i1,\alpha p(1t)))$$
(29)
so that:

$`{\displaystyle \frac{d}{dt}}E\mathrm{log}Z(t)=`$ $`N\alpha {\displaystyle \underset{M}{}}(\pi (M,\alpha tN)\pi (M1,\alpha tN))E^{}\mathrm{log}Z(t)`$ (31)
$`+\alpha p{\displaystyle \underset{i}{}}{\displaystyle \underset{k_i}{}}(\pi (k_i,\alpha p(1t))\pi (k_i1,\alpha p(1t)))E_i^{\prime \prime }\mathrm{log}Z(t)`$
where we have denoted as $`E^{}`$ the average with respect to all the quenched variables except $`M`$ and with $`E_i^{\prime \prime }`$ the average with respect to all the quenched variables except $`k_i`$, and simply $`Z(t)`$ the partition function of the $`N`$ spin system $`Z_N(t)`$.

In the first term of (31) we can single out the $`M`$-th clause, and write $`Z(t)=Z_M(t)\omega (\mathrm{e}^{\beta H_M(S_{i_1^M},\mathrm{},S_{i_p^M}})_M`$, where by $`Z_M(t)`$ we denote the partition function of the system in absence of the $`M`$-th clause, and $`\omega ()_M`$ is the canonical average in absence of the $`M`$-th clause. In the following terms we single out the $`k_i`$-th field $`u`$ term, $`Z(t)=Z_{u_i^{k_i}}(t)\omega (\mathrm{e}^{\beta u_i^{k_i}S_i})_{u_i^{k_i}}`$, where $`Z_{u_i^{k_i}}(t)`$ is the partition function in absence of the field $`u_i^{k_i}`$ and analogously for the average $`\omega ()_{u_i^{k_i}}`$. Finally, rearranging all terms we find

$`{\displaystyle \frac{d}{dt}}E\mathrm{log}Z(t)=`$ $`N\alpha {\displaystyle \underset{M}{}}(\pi (M1,\alpha tN))E^{}\mathrm{log}[\omega (\mathrm{e}^{\beta H_{J^{(M)}}(S_{i_1^M},\mathrm{},S_{i_p^M})})_M)]`$ (33)
$`p\alpha {\displaystyle \underset{i}{}}{\displaystyle \underset{k_i}{}}\pi (k_i1,\alpha p(1t))E_i^{\prime \prime }\mathrm{log}[\omega (\mathrm{e}^{\beta u_i^{k_i}S_i})_{u_i^{k_i}}].`$
where we have used $`_M\pi (M1,\alpha tN)E^{}\mathrm{log}Z_M=_{k_i}\pi (k_i1,\alpha p(1t))E_i^{\prime \prime }\mathrm{log}Z_{u_i^{k_i}}=E\mathrm{log}Z`$. We notice at this point that the statistical ensemble defined by $`\pi (M1,\alpha tN))E^{}`$ can be substituted with the original one $`E`$ and the average of the variables appearing in the clause we have singled out. To be more precise, we remark that the average $`\omega ()`$ depends on the quenched variables $`D=\{𝐉,𝐮\}`$ appearing in the Hamiltonian. Writing explicitly this dependence as $`\omega (|D)`$, and denoting as $`D_M`$ all the quenched variables except the ones appearing in the $`M`$-th clause, our statement is that thanks to the Poissonian distribution of $`M`$ and the uniform choice of the indices of each clause,

$$\underset{M}{}(\pi (M1,\alpha tN))E^{}\mathrm{log}[\omega (\mathrm{e}^{\beta H_{J^{(M)}}(S_{i_1^M},\mathrm{},S_{i_p^M})}|D_M)]=E\left(\frac{1}{N^p}\underset{i_1,\mathrm{},i_p}{}\mathrm{log}[\omega (\mathrm{e}^{\beta H_J(S_{i_1},\mathrm{},S_{i_p})}|D)]_J\right).$$
(34)
where by $`_J`$ we denote the average with respect to the random variables appearing in the clause. This is a crucial step in our analysis, in fact, similar considerations apply to the term in the second line of (33), which can be written as

$`{\displaystyle \underset{k_i}{}}\pi (k_i1,\alpha p(1t))E_i^{\prime \prime }\mathrm{log}[\omega (\mathrm{e}^{\beta u_i^{k_i}S_i})_{u_i^{k_i}}]=`$ $`E\mathrm{log}\omega \left(e^{\beta uS_i}\right)_u.`$ (35)
The same kind of averages $`E`$ and $`\omega `$ appear in the two terms which can be therefore directly compared as we will do in the next section. This property, linked to the Poissonian character of the graph defined by the model would not hold for other ensembles of random graphs and the analysis would be technically more involved. Substituting in (33) we find:

$$\frac{1}{N}\frac{d}{dt}E\mathrm{log}Z(t)=\alpha E\left[\frac{1}{N^p}\underset{i_1,\mathrm{},i_p}{}\mathrm{log}[\omega (\mathrm{e}^{\beta H(S_{i_1},\mathrm{},S_{i_p})})]_H\frac{p}{N}\underset{i}{}\mathrm{log}\omega (e^{\beta uS_i})_u\right]$$
(36)
Rearranging terms and using (27) we finally find that the free-energy $`F_N`$ can be written as

$$F_N=F_{var}[G]+_0^1𝑑tR_{RS}[G,t]+O(1/N)$$
(37)
where $`F_{var}[G]`$ coincides the expression of the variational free-energy in the replica treatment under condition $`G[h]=P[h]`$ $`h`$ at $`t=0`$ and $`_0^1𝑑tR_{RS}[G,t]`$ is a remainder term. Instead of writing the formulae for general clauses, in order to keep the notations within reasonable simplicity, we specialize now to the specific cases of the $`p`$-spin model and the K-SAT. Notice that in all models

$$F[0]=\frac{1}{\beta }\mathrm{log}(2\mathrm{cosh}(\beta h))_h|_{t=0}$$
(38)
#### 1 $`p`$-spin

In the case of the $`p`$-spin $`H_J(S_{i_1},\mathrm{},S_{i_p})=JS_{i_1}\mathrm{}S_{i_p}`$. Substituting in eq.(36) and rearranging terms one immediately finds:

$`F_{var}^{pspin}[G]`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}[\alpha (p\mathrm{log}(\mathrm{cosh}\beta u)_u\mathrm{log}(\mathrm{cosh}\beta J)_J)\mathrm{log}(2\mathrm{cosh}\beta h)_h+`$ (40)
$`\alpha (p1)\mathrm{log}(1+\mathrm{tanh}(\beta J){\displaystyle \underset{t=1}{\overset{p}{}}}\mathrm{tanh}(\beta g_t))_{\{g_t\},J}]`$
while the remainder is the $`t`$ integral of

$`R_{RS}^{pspin}[G,t]=`$ $`{\displaystyle \frac{\alpha }{\beta }}[{\displaystyle \frac{1}{N^p}}{\displaystyle \underset{i_1,\mathrm{},i_p}{}}E\mathrm{log}(1+\mathrm{tanh}(\beta J)\omega (S_{i_1}\mathrm{}S_{i_p}))_JpE\mathrm{log}(1+\mathrm{tanh}(\beta u)\omega (S_i))_u+`$ (42)
$`(p1)E\mathrm{log}(1+\mathrm{tanh}(\beta J){\displaystyle \underset{t=1}{\overset{p}{}}}\mathrm{tanh}(\beta g_p))_{\{g_t\},J}].`$
The expression for $`F_{var}^{pspin}[G]`$ coincides with the RS free energy once extremized over the variational space of probability distributions $`G`$ . Terms have been properly added and subtracted in order to get a remainder which equal to zero if maximization over $`G`$ is taken, and the temperature is high enough for replica symmetry to be exact . As we will see, the remainder turns out to be positive. $`F_{var}^{pspin}[G]`$ is therefore, for all $`G`$ for which its expression makes sense, a lower bound to the free-energy. At saturation the condition

$$G[h]=P[h]|_{t=0}h$$
(43)
should hold, which is simply the self-consistency RS equation.

By using equation

$$E\mathrm{log}(1+\mathrm{tanh}(\beta u)\omega (S_i))_u=E\mathrm{log}(1+\mathrm{tanh}(\beta J)\underset{t=1}{\overset{p1}{}}\mathrm{tanh}(\beta g_r)\omega (S_i))_{\{g_t\},J}$$
(44)
we can establish that the remainder is positive for even $`p`$. We expand the logarithm of the three terms in (absolutely converging) series of $`\mathrm{tanh}(\beta J)`$, and notice that thanks to the parity of the $`J`$ and the $`g`$ distributions, they will just involve negative terms. We can then take the expected value of each terms and write

$$R_{RS}^{pspin}[G,t]=\frac{1}{\beta }\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{tanh}^{2n}\beta J_J\frac{1}{n}\mathrm{\Omega }\left[(q^{(2n)})^ppq^{(2n)}\mathrm{tanh}^{2n}\beta g_g^{p1}+(p1)\mathrm{tanh}^{2n}\beta g_g^p\right]$$
(45)
where we have introduced the overlap $`q^{(l)}`$ and the replica measure $`\mathrm{\Omega }`$ defined in section 2. The series in (45) is an average of positive terms in the case of the Viana-Bray model $`p=2`$, where we get perfect squares, and more in general for all even $`p`$, as we can easily, starting from the observation that in this case $`x^ppxy^{p1}+(p1)y^p`$ is positive or zero for all $`x=q^{(2n)}`$ ,$`y=\mathrm{tanh}^{2n}\beta J_J`$ real.

In the case of $`p`$ odd, the same term is positive only if $`x`$ is itself positive or zero. The bound of the free-energy would therefore be established if we were able to prove that the probability distributions of the $`q^{(2n)}`$ has support on the positives.<sup>*</sup><sup>*</sup>*A different sufficient condition for the series to have positive terms is that $`|\mathrm{\Omega }(q^{(2n)})|\mathrm{tanh}(\beta g)^{2n}_g`$, but it is not clear its physical meaning. This property, which tells that anti-correlated states are not possible, is physically very sound whenever the Hamiltonian is not symmetric under change of sign of all spins. In fact, one expects the probability of negative values of the overlaps to be exponentially small in the size of the system for large $`N`$. Unfortunately however we have not been able to prove this property in full generality. Notice that upon maximization on $`G`$, the results of  imply that the remainder is exactly equal to zero if the temperature is high enough for replica symmetry to hold.

#### 2 K-SAT

In the case of the $`K`$-SAT, using def.(3) for the clause $`H`$, we find relation:

$$u_𝐉(h_1,\mathrm{},h_{p1})u_J(\{J_t\},\{h_t\})=\frac{J}{\beta }\mathrm{tanh}^1\left[\frac{\frac{\xi }{2}_{t=1}^{p1}\left(\frac{1+J_t\mathrm{tanh}(\beta h_t)}{2}\right)}{1+\frac{\xi }{2}_{t=1}^{p1}\left(\frac{1+J_t\mathrm{tanh}(\beta h_t)}{2}\right)}\right],$$
(46)
where $`\xi e^\beta 1<0`$. Via direct inspection, the variational free-energy coincides with the RS expression

$`F_{var}^{KSAT}[G]`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}[\alpha (p1)\mathrm{log}(1+(e^\beta 1){\displaystyle \underset{t=1}{\overset{p}{}}}\left({\displaystyle \frac{1+\mathrm{tanh}(\beta g_t)}{2}}\right))_{\{g_t\},\{J_t\}}\mathrm{log}(2\mathrm{cosh}(\beta h))_h+`$ (48)
$`\alpha p\mathrm{log}(2\mathrm{cosh}(\beta u))_u\alpha p\mathrm{log}(1+{\displaystyle \frac{(e^\beta 1)}{2}}{\displaystyle \underset{t=1}{\overset{p1}{}}}\left({\displaystyle \frac{1+\mathrm{tanh}(\beta g_t)}{2}}\right))_{\{g_t\},\{J_t\}}]`$
while the remainder is the $`t`$ integral of

$`R_{RS}^{KSAT}[G,t]`$ $`=`$ $`{\displaystyle \frac{\alpha }{\beta }}E[{\displaystyle \frac{1}{N^p}}{\displaystyle \underset{i_1,\mathrm{},i_p}{}}\mathrm{log}(1+(\mathrm{e}^\beta 1)\omega ({\displaystyle \underset{t=1}{\overset{p}{}}}{\displaystyle \frac{1+J_tS_{i_t}}{2}}))_{\{J_t\}}`$ (51)
$`{\displaystyle \frac{p}{N}}{\displaystyle \underset{i}{}}\mathrm{log}\left(1+\xi \omega \left({\displaystyle \frac{1+JS_i}{2}}{\displaystyle \underset{t=1}{\overset{p1}{}}}{\displaystyle \frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}}\right)\right)_{\{g_t\},J,\{J_t\}}+`$
$`(p1)\mathrm{log}(1+\xi {\displaystyle \underset{t=1}{\overset{p}{}}}{\displaystyle \frac{1+J_r\mathrm{tanh}(\beta g_t)}{2}})_{\{g_t\},\{J_t\}}].`$
Considerations analogous to the case of the $`p`$-spin, have led us to add and subtract terms from eq.(36) to single out the proper remainder term. Expanding in series the logarithms, exploiting the symmetry of the probabilities distribution functions and taking the expectation of each term of the absolutely convergent series we finally obtain:

$$R_{RS}^{KSAT}[G,t]=\frac{\alpha }{\beta }\underset{n1}{}\frac{(1)^n}{n}(\xi ^{})^n\mathrm{\Omega }\left[(1+Q_n)^pp(1+Q_n)(1+J\mathrm{tanh}(\beta g))^n_{J,g}^{p1}+(p1)(1+J\mathrm{tanh}(\beta g))^n_{J,g}^p\right]$$
(52)
where we have defined $`\xi ^{}\xi /(2^p)<0`$ and $`Q_n_{l=1}^nJ^l_J_{a_1<\mathrm{}<a_l}^{1,n}q^{a_1\mathrm{}a_l}`$. Detailed calculations are given in the appendix. As in the $`p`$-spin case, the previous sum is obviously positive for $`p`$ even. For $`p`$ odd we should again rely on the physical wisdom that all $`q^{(a_1,\mathrm{},a_l)}`$ have positive support and so have the functions $`1+Q_n0`$. Again, the variational free-energy coincides with the RS expression once extremized over $`G`$ at the condition $`P=G`$ at $`t=0`$.

### B The 1RSB Bound

We establish here a more complex estimate, in a larger variational space of functional probability distributions. The general strategy will be here to consider the same form for the auxiliary Hamiltonian, but now with a more involved choice for the fields distribution. The fields on different sites or different index $`l_i`$ will be still independent, but each site field distribution $`G_i^{l_i}(g_i^{l_i})`$ will be itself random i.i.d., chosen with a probability density functional $`𝒢[G]`$, with support on symmetric distributions $`G(g)=G(g)`$. It will be assumed that $`𝒢`$ is such that all the expressions below make sense. In this case, the variational approximation for the free-energy will be obtained from an estimate of

$$\beta F_N[m,t]=\frac{1}{mN}E_1\mathrm{log}E_2(Z^m(t))$$
(53)
where we have denoted with:

* $`E_2`$ the average w.r.t. $`g_i^{l_i,n}`$ for fixed distributions $`G_i^{l_i,n}`$ according to the measure

$$C\underset{i=1}{\overset{N}{}}\underset{l_i=1}{\overset{k_i}{}}\underset{n=1}{\overset{p1}{}}dg_i^{l_i,n}G_i^{l_i,n}(g_i^{l_i,n})\left(\frac{B_{J_i^{l_i,n}}(g_i^{l_i,1}\mathrm{}g_i^{l_i,p1})}{2\mathrm{cosh}(\beta u_{J_i^{l_i,n}}(g_i^{l_i,1}\mathrm{}g_i^{l_i,p1}))}\right)^m$$
(54)
where $`C`$ ensures the normalization.
* $`E_1`$ the average with respect to the quenched clause variable, distributions the $`G_i^{l_i}`$’s and the Poissonian variables $`k_i`$’s, which will be i.i.d. with probabilities $`\mu (J)`$, $`𝒢(G_i^{l_i})`$ and $`\pi (k_i,(1t)\alpha )`$ respectively.

The number $`m`$ is real in the interval (0,1\]. The statistical ensemble of the auxiliary fields $`u`$ and $`h`$ will be now related to the one of the $`g`$ by:

$$𝒬(Q)=𝒟G_1𝒢(G_1)\mathrm{}𝒟G_{p1}𝒢(G_{p1})\delta (Q()Q(|G_1,\mathrm{},G_{p1},J))_J$$
(55)
$$𝒫(P)=\underset{k=0}{\overset{\mathrm{}}{}}\mathrm{e}^{\alpha p(1t)}\frac{(\alpha p(1t))^k}{k!}\underset{l=1}{\overset{k}{}}𝒟Q_l𝒬(Q_l)\delta (P()P(|Q_1,\mathrm{},Q_k))$$
(56)
where:

$`Q(u|G_1,\mathrm{},G_{p1},J)=𝒩_G[G_1,\mathrm{},G_{p1}]{\displaystyle 𝑑g_1G_1(g_1)\mathrm{}𝑑g_{p1}G_1(g_{p1})B_J(g_1,\mathrm{},g_{p1})^m\delta (uu_J(g_1,\mathrm{},g_{p1}))}`$ (57)
$`G(g|Q_1,\mathrm{},Q_k)=𝒩_{Q,k}[Q_1,\mathrm{},Q_k](2\mathrm{cosh}(\beta g))^m{\displaystyle \underset{l=1}{\overset{k}{}}du_l\frac{Q_l(u_l)}{(2\mathrm{cosh}(\beta u_l))^m}\delta (g\underset{l=1}{\overset{k}{}}u_l)}`$ (58)
where $`𝒩_{Q,k}[Q_1,\mathrm{},Q_k]`$, $`𝒩_G[G_1,\mathrm{},G_{p1}]`$ and $`B_J(g_1,\mathrm{},g_{p1})`$ have been previously defined. With notations similar to the ones of the RS case, we can write

$`{\displaystyle \frac{d}{dt}}(\beta F_N[m,t])=`$ $`\alpha {\displaystyle \underset{M}{}}(\pi (M,\alpha tN)\pi (M1,\alpha tN))E_1^{}{\displaystyle \frac{1}{Nm}}\mathrm{log}E_2Z(t)^m`$ (60)
$`+{\displaystyle \frac{\alpha p}{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{k_i}{}}(\pi (k_i,\alpha p(1t))\pi (k_i1,\alpha p(1t)))E_{1,i}^{\prime \prime }{\displaystyle \frac{1}{Nm}}\mathrm{log}E_2Z(t)^m`$
which, extracting explicitly the contribution from the $`M`$-th close in the first term and the $`k_i`$-th field $`u`$ in the second, following considerations similar to the RS case we find:

$`{\displaystyle \frac{d}{dt}}(\beta F_N[m,t])=`$ $`\alpha {\displaystyle \underset{M}{}}(\pi (M1,\alpha tN)){\displaystyle \frac{1}{m}}E_1^{}\mathrm{log}\left[{\displaystyle \frac{E_2Z_M^m\omega (\mathrm{e}^{\beta H_{J^{(\mu )}}(S_{i_1^M},\mathrm{},S_{i_p^M})})_M^m}{E_2Z_M^m}}\right]`$ (62)
$`{\displaystyle \frac{p\alpha }{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{k_i}{}}\pi (k_i1,\alpha p(1t)){\displaystyle \frac{1}{m}}E_{1,i}^{\prime \prime }\mathrm{log}\left[{\displaystyle \frac{E_2Z_{u_i^{k_i}}^m\omega (\mathrm{e}^{\beta u_i^{k_i}S_i})_{u_i^{k_i}}^m}{E_2Z_{u_i^{k_i}}^m}}\right].`$
Again it can be recognized that the primed averages coincide with the averages over the original ensembles plus the averages on the variables appearing in the terms we extracted. Finally we get:

$$\frac{d}{dt}(\beta F_N[m,t])=\frac{\alpha }{m}E_1\left[\frac{1}{N^p}\underset{i_1,\mathrm{},i_p}{}\mathrm{log}\left(\frac{E_2Z^m\omega (e^{\beta H_J(S_{i_1},\mathrm{},S_{i_p})})^m}{E_2Z^m}\right)_J\frac{p}{N}\underset{i}{}\mathrm{log}\left(\frac{E_2Z^m\omega (e^{\beta uS_i})^m_u}{E_2Z^m}\right)_Q\right].$$
(63)
Rearranging all terms one finds the estimate:

$$F_N=F_{var}[𝒢]+_0^1𝑑tR_{1RSB}[𝒢,t]+O(1/N)$$
(64)
where this time $`F_{var}[𝒢]`$ coincides with $`F_{1RSB}[𝒢]`$, the expression of the variational free-energy in the 1RSB treatment at the saddle point $`𝒢=𝒫`$ at $`t=0`$, and $`_0^1𝑑tR_{1RSB}[𝒢,t]`$ is the remainder. Notice that the derivation immediately suggests how to generalize the analysis to more steps of replica symmetry breaking. Let us now specialize the formulae for the $`p`$-spin model and the K-SAT. Again, in this case we will need the expression for $`F[0]`$:

$$F[0]=\frac{1}{\beta m}\left[\mathrm{log}\left(\frac{1}{2\mathrm{cosh}(\beta h)}\right)^m_h_P\right]_{|_{t=0}}.$$
(65)
#### 1 $`p`$-spin

In this case, plugging def.(2) in eq.(63) rearranging, adding and subtracting terms one finds:

$`F_{var}^{pspin}[𝒢]`$ $`=`$ $`{\displaystyle \frac{1}{\beta m}}[\mathrm{log}\left({\displaystyle \frac{1}{2\mathrm{cosh}(\beta h)}}\right)^m_h_P\alpha m\mathrm{log}(2\mathrm{cosh}(\beta J))_J\alpha p\mathrm{log}\left({\displaystyle \frac{1}{2\mathrm{cosh}(\beta u)}}\right)^m_u_Q`$ (67)
$`+\alpha (p1)\mathrm{log}(1+\mathrm{tanh}(\beta J)\mathrm{tanh}(\beta g_1)\mathrm{}\mathrm{tanh}(\beta g_p))^m_{g_1,\mathrm{},g_p}_{G_1,\mathrm{},G_p;J}]`$
while the remainder is the $`t`$ integral of

$`R_{1RSB}^{pspin}[𝒢,t]=`$ $`{\displaystyle \frac{\alpha }{\beta m}}E_1[{\displaystyle \frac{1}{N^p}}{\displaystyle \underset{i_1,\mathrm{},i_p}{}}\mathrm{log}\left({\displaystyle \frac{E_2Z^m(1+\omega (S_{i_1}\mathrm{}S_{i_p})\mathrm{tanh}(\beta J))^m}{E_2Z^m}}\right)_J`$ (70)
$`{\displaystyle \frac{p}{N}}{\displaystyle \underset{i}{}}\mathrm{log}\left({\displaystyle \frac{E_2Z^m(1+\omega (S_i)\mathrm{tanh}(\beta u))^m_u}{E_2Z^m}}\right)_Q+`$
$`(p1)\mathrm{log}(1+\mathrm{tanh}(\beta J)\mathrm{tanh}(\beta g_1)\mathrm{}\mathrm{tanh}(\beta g_p))^m_{g_1,\mathrm{},g_p}_{G_1,\mathrm{},G_p;J}]`$
The expression for $`F_{var}^{pspin}[𝒢]`$ coincides with the 1RSB free-energy  once maximized over the variational space of probability distribution functionals $`𝒢`$. The maximization condition reads:

$$𝒢[P]=𝒫[P]|_{t=0}P,$$
(71)
which is simply the self consistency 1RSB condition. For even $`p`$ (and in particular for $`p=2`$ that corresponds to the Viana-Bray case), one can check that the remainder is positive just expanding the logarithm in series and exploiting the parity of the $`J`$ and the $`g`$ distributions. As this is considerably more involved then in the RS case, we relegate this check to appendix A.

#### 2 $`K`$-SAT

In the $`K`$-SAT case the expression for function $`B_𝐉(h_1,\mathrm{},h_{p1})`$ reads:

$$B_𝐉(h_1,\mathrm{},h_{p1})B(\{J_t\},\{h_t\})=1+\frac{\xi }{2}\underset{t=1}{\overset{p1}{}}\left(\frac{1+J_t\mathrm{tanh}(\beta h_t)}{2}\right),$$
(72)
while the corresponding one for $`u_𝐉(h_1,\mathrm{},h_{p1})`$ is the same as in the RS case. The corresponding replica free-energy and remainder read

$`F_{var}^{KSAT}[𝒢]`$ $`=`$ $`{\displaystyle \frac{1}{m\beta }}[\alpha (p1)\mathrm{log}(1+\xi {\displaystyle \underset{t=1}{\overset{p}{}}}\left({\displaystyle \frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}}\right))^m_{\{g_t\}}_{\{G_t\},\{J_t\}}`$ (74)
$`\alpha p\mathrm{log}\left({\displaystyle \frac{B(\{J_t\},\{g_t\})}{2\mathrm{cosh}(\beta u_J(\{J_t\},\{g_t\}))}}\right)^m_{\{g_t\}}_{\{G_t\},\{J_t\},J}+\mathrm{log}\left({\displaystyle \frac{1}{2\mathrm{cosh}(\beta h)}}\right)^m_h_P]`$
The remainder is the $`t`$ integral of

$`R_{1RSB}^{KSAT}[𝒢,t]`$ $`=`$ $`{\displaystyle \frac{\alpha }{\beta m}}E_1[{\displaystyle \frac{1}{N^p}}{\displaystyle \underset{i_1,\mathrm{},i_p}{}}\mathrm{log}\left({\displaystyle \frac{E_2Z^m\left(1+\xi \omega \left(_{t=1}^p\frac{1+J_tS_{i_t}}{2}\right)\right)^m}{E_2Z^m}}\right)_{\{J_t\}}`$ (77)
$`{\displaystyle \frac{p}{N}}{\displaystyle \underset{i}{}}\mathrm{log}\left({\displaystyle \frac{E_2Z^m\left(1+\xi \frac{1+J\omega (S_i)}{2}_{t=1}^{p1}\frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}\right)^m_{\{g_t\}}}{E_2Z^m}}\right)_{\{G_t\},\{J_t\},J}+`$
$`(p1)\mathrm{log}(1+\xi {\displaystyle \underset{t=1}{\overset{p}{}}}\left({\displaystyle \frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}}\right))^m_{\{g_t\}}_{\{G_t\},\{J_t\}}]`$
The expression for $`F_{var}^{KSAT}[𝒢]`$ coincides with the 1RSB free energy once extremized under condition (71), with the corresponding $`K`$-SAT probability distribution functionals. Notice that The proof of the positivity of (77) for even $`p`$ is again dove via series expansion, all the detail are explained in Appendix B.

At this point we can take the zero temperature limit, finding that the resulting expression gives us a lower bound for the ground-state energy of the system, i.e. the minimal number of unsatisfied clauses. Notice that the $`T0`$ limit of the replica free-energy is not trivial. The necessary assumptions on the field distributions to get it correct are well known in the physical literature, and have been recently reviewed in . Recently Mezard, Parisi and Zecchina  have worked out the K-SAT 1RSB solution for $`p=3`$ predicting a non zero ground-state energy for values of $`\alpha `$ above a satisfiability threshold of $`\alpha _c=4.256`$, very well in agreement with the numerical simulations. Our results, together with the additional hypothesis of positivity of the support of the overlap functions imply that this value is an upper bound to the true threshold.

## IV Summary and conclusions

In this paper we have established that the free-energy of some families of diluted random spin models can be written as the sum of a term identical to the ones got in the cavity/replica plus an error term. Both the replica term and the remainder are different in different replica scheme, corresponding to the choice of statistical ensemble of the cavity fields. We believe that the sign of the remainder is in general negative in the model we have considered, although we have been able to prove that only in the case of even $`p`$. For odd $`p`$ our belief is supported by the physical wisdom that the overlap distributions are supported on the positives in the large $`N`$ limit.

We have considered the cases of replica symmetry and one step of replica symmetry breaking. It is clear that the analysis could be extended to further levels of replica symmetry breaking, although the complexity of the analysis would greatly increase. The 1RSB level is thought to give the exact scheme to treat the $`p`$ spin model and the K-SAT problem for $`p3`$. For the Viana-Bray model on the other hand it is believed that no finite RSB scheme furnish the exact solution, and one needs to consider the limit of infinite number of replica symmetry breaking. It is not clear to us how to generalize the analysis to this case.

Our analysis of the diluted models underlines a strong link between the Guerra method and the cavity method which remained rather hidden in the fully connected case. In the cavity approach one considers incomplete graphs in which either sites or clauses are removed from the complete graph. Then, with the aid of precise physical hypothesis, consistency equations are written that allow to compute the free-energy from the comparison between the site and clause contributions. In the approach presented in this paper the removal of clauses is compensated in average by the addition of some external fields which have precisely the statistics which is assumed with cavity. The novelty of the approach is that it gives some control on the approximation involved, and proves the variational nature of the replica free-energies. Of course a complete control on the remainder in various situations would result in rigorous solutions.

Although in this paper we have mainly worked at finite temperature, the zero temperature limit can be considered without harm. This is particularly relevant in random satisfiability problem, where it is typically found a SAT-UNSAT transition where the ground state energy passes from zero to non zero values. Our analysis shows that, the replica estimates for many of the models considered in the literature, and possibly some of the ones to be obtained in the future with the same method provide upper bounds for the satisfiability thresholds.

In this paper we have confined ourselves to spin models on graphs with Poissonian connectivity. The extension to more general diluted graph models will be presented in a forthcoming paper.

Finally we would like to remark that despite the heavy formalism, our proofs to the bounds are conceptually simple. They are issue of explicit computations and elementary positivity consideration. We hope that this contributes to illustrate the elegance of the construction first introduced in .

## V Acknowledgments

We thank M. Isopi for pointing us out the importance of ref. . We thank G. Gaeta for discussions.

## A p-spin

### 1 Check of the positive sign of $`R_{1RSB}^{pspin}`$

In this appendix we will explicitly show that expression (70) has positive sign. We proceed expanding in absolutely convergent series each of the three addend in (70) and showing, taking the expectation of each term that the resulting series is positive semidefinite.

The first term writes:

$`E_1\mathrm{log}{\displaystyle \frac{E_2Z^m(1+\mathrm{tanh}(\beta J)\omega (S_{i_1}\mathrm{}S_{i_p}))^m}{E_2Z^m}}_J=`$ (A1)
$`{\displaystyle \underset{l1}{}}{\displaystyle \frac{(1)^{l+1}}{l}}{\displaystyle \underset{k_1,\mathrm{},k_l}{\overset{1,\mathrm{}}{}}}\left({\displaystyle \frac{m(1)^{k_11}}{k_1!}}{\displaystyle \underset{r_1=1}{\overset{k_11}{}}}(r_1m)\right)\mathrm{}\left({\displaystyle \frac{m(1)^{k_l1}}{k_l!}}{\displaystyle \underset{r_l=1}{\overset{k_l1}{}}}(r_lm)\right)(\mathrm{tanh}(\beta J)^{_{s=1}^lk_s}_J`$ (A2)
$`E_1\left({\displaystyle \underset{s=1}{\overset{l}{}}}{\displaystyle \frac{E_2(Z^m\omega (S_{i_1}\mathrm{}S_{i_p})^{k_s})}{E_2(Z^m)}}\right)`$ (A3)
where the term $`E_1()`$ in the last line of eq.(A3) can be written as

$$E_1\left(\frac{E_2^{(1)}\mathrm{}E_2^{(l)}Z_{(1)}^m\mathrm{}Z_{(l)}^m\omega _{(1)}(S_{i_1}^{1,1}\mathrm{}S_{i_p}^{1,1}\mathrm{}S_{i_1}^{k_1,1}\mathrm{}S_{i_p}^{k_1,1})\mathrm{}\omega _{(l)}(S_{i_1}^{1,l}\mathrm{}S_{i_p}^{1,l}\mathrm{}S_{i_1}^{k_l,l}\mathrm{}S_{i_p}^{k_l,l})}{(E_2Z^m)^l}\right)=\mathrm{\Omega }^{(l)}\left[(q^{(k_1,\mathrm{},k_l)})^p\right].$$
(A4)
where each $`\omega _{(s)}`$ ($`s=1,\mathrm{},l`$) is a product of $`k_s`$ Gibbs measure with independent fields (variables appearing in the $`E_2^{(s)}`$ averages), and same fields distributions and quenched disorder (variables appearing in $`E_1`$). The quantities $`q^{(k_1,\mathrm{},k_l)}`$ have been defined as:

$$q^{(k_1,\mathrm{},k_l)}=\frac{1}{N}\underset{i}{}S_i^{1,1}\mathrm{}S_i^{k_1,1}\mathrm{}S_i^{1,l}\mathrm{}S_i^{k_l,l}$$
(A5)
and in this case the averages are performed using a a generalized replica measure, defined as:

$$\mathrm{\Omega }^{(l)}[(q^{(k_1,\mathrm{},k_l)})^n]=E_1\left[\frac{_{s=1}^lE_2^{(s)}Z_{(s)}^m\omega _{(s)}(S_{i_1}\mathrm{}S_{i_n})^{k_s}}{(E_2Z^m)^l}\right]$$
(A6)
for any integer $`n`$. The average over $`J`$ selects the terms with even $`_{s=1}^lk_l`$ in (A3) so that we finally find

$$\underset{l1}{}\frac{m^l}{l}\underset{\genfrac{}{}{0pt}{}{k_1,\mathrm{},k_l}{_{s=1}^lk_s\mathrm{even}}}{\overset{1,\mathrm{}}{}}\underset{s=1}{\overset{l}{}}\left(\frac{_{r_s=1}^{k_s1}(r_sm)}{k_s!}\right)(\mathrm{tanh}(\beta J))^{_{s=1}^lk_s}_J\mathrm{\Omega }^{(l)}\left[(q^{(k_1,\mathrm{},k_l)})^p\right]$$
(A7)
notice that $`(r_sm)0`$ $``$ integer $`r_s>0`$ only in the current hypothesis that $`m[0,1]`$. Analogously, the term

$$E_1\mathrm{log}\frac{E_2Z^m(1+\mathrm{tanh}(\beta u)\omega (S_i))^m_u}{E_2Z^m}_Q$$
(A8)
writes

$$\underset{l1}{}\frac{m^l}{l}\underset{\genfrac{}{}{0pt}{}{k_1,\mathrm{},k_l}{_{s=1}^lk_s\mathrm{even}}}{\overset{1,\mathrm{}}{}}\underset{s=1}{\overset{l}{}}\left(\frac{_{r=1}^{k_s1}(rm)}{k_s!}\right)\underset{s=1}{\overset{l}{}}\mathrm{tanh}(\beta u)^{k_s}_u_Q\mathrm{\Omega }^{(l)}\left[(q^{(k_1,\mathrm{},k_l)})\right]$$
(A9)
or, making use of the definition of $`G(g)`$,

$$\underset{l1}{}\frac{1}{l}\underset{\genfrac{}{}{0pt}{}{k_1,\mathrm{},k_l}{_{s=1}^lk_s\mathrm{even}}}{\overset{1,\mathrm{}}{}}\underset{s=1}{\overset{l}{}}\left(\frac{_{r=1}^{k_s1}(rm)}{k_s!}\right)\underset{s=1}{\overset{l}{}}(\mathrm{tanh}(\beta g))^{k_s}_g_G^{p1}(\mathrm{tanh}(\beta J))^{_{s=1}^lk_s}_J\mathrm{\Omega }^{(l)}\left[(q^{(k_1,\mathrm{},k_l)})\right]$$
(A10)
Eventually, following analogous manipulations, the last term

$$\mathrm{log}(1+\mathrm{tanh}(\beta J)\underset{t=1}{\overset{p}{}}\mathrm{tanh}(\beta g_t))^m_{\{g_t\}}_{J,\{G_t\}}$$
(A11)
can be written as

$$\underset{l1}{}\frac{m^l}{l}\underset{\genfrac{}{}{0pt}{}{k_1,\mathrm{},k_l}{_{s=1}^lk_s\mathrm{even}}}{\overset{1,\mathrm{}}{}}\underset{s=1}{\overset{l}{}}\left(\frac{_{r=1}^{k_s1}(rm)}{k_s!}\right)\underset{s=1}{\overset{l}{}}(\mathrm{tanh}(\beta g))^{k_s}_g_G^p(\mathrm{tanh}(\beta J))^{_{s=1}^lk_s}_J.$$
(A12)
Invoking 55 and collecting all

$`R_{1RSB}^{pspin}[𝒢,t]`$ $`=`$ $`{\displaystyle \frac{\alpha }{\beta m}}{\displaystyle \underset{l1}{}}{\displaystyle \frac{m^l}{l}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{k_1,\mathrm{},k_l}{_{s=1}^lk_s\mathrm{even}}}{\overset{1,\mathrm{}}{}}}{\displaystyle \underset{s=1}{\overset{l}{}}}\left({\displaystyle \frac{_{r=1}^{k_s1}(rm)}{k_s!}}\right)\left(\mathrm{tanh}(\beta J)\right)^{_{s=1}^lk_s}_J`$ (A14)
$`\mathrm{\Omega }^{(l)}\left[(q^{(k_1,\mathrm{},k_l)})^ppA(k_1,\mathrm{},k_l)^{p1}(q^{(k_1,\mathrm{},k_l)})+(p1)A(k_1,\mathrm{},k_l)l^p\right]`$
where we have defined:

$$A(k_1,\mathrm{},k_l)\underset{s=1}{\overset{l}{}}(\mathrm{tanh}(\beta g))^{k_s}_g_G$$
(A15)
Each inner term of the series (A14)

$$\mathrm{\Omega }^{(l)}\left[(q^{(k_1,\mathrm{},k_l)})^ppA(k_1,\mathrm{},k_l)^{p1}(q^{(k_1,\mathrm{},k_l)})+(p1)A(k_1,\mathrm{},k_l)^p\right]$$
(A16)
is always positive semidefinite for $`p`$ even while we need the condition conditions $`q^{(k_1,\mathrm{},k_l)}0`$ for $`p`$ odd. For $`p=2`$ one retrieves the Viana-Bray result where (A16) is a perfect square. As in the RS case, one can now integrate eq.(A14) and recognize that once more the total true free-energy can be written as variational term plus a positive extra one. The variational term coincides with the 1RSB free-energy at stationarity and under condition

$`𝒢(P)`$ $`=`$ $`𝒫(P)|_{t=0}P.`$ (A17)
## B K-SAT

### 1 Check of the positive sign of $`R_{RS}^{KSAT}`$

The aim of this appendix is to show that the expression for the remainder $`R_{RS}[G,t]`$ in (37) for the $`K`$-SAT model case as positive sign. For the $`K`$-SAT $`R_{RS}[G,t]`$ specializes to:The sum of the site indices has been eliminated by symmetry.

$`R_{RS}^{KSAT}[G,t]`$ $`=`$ $`{\displaystyle \frac{\alpha }{\beta }}E[\mathrm{log}\left(\omega \left(\mathrm{exp}^{\beta _{r=1}^p\frac{1+J_rS_r}{2}}\right)\right)_{\{J_t\}}`$ (B3)
$`p\mathrm{log}(1+\omega (S)\mathrm{tanh}(\beta u))_up\mathrm{log}(1+{\displaystyle \frac{\xi }{2}}{\displaystyle \underset{t=1}{\overset{p1}{}}}\left({\displaystyle \frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}}\right))_{\{g_t\},\{J_t\}}+`$
$`(p1)\mathrm{log}(1+\xi {\displaystyle \underset{t=1}{\overset{p}{}}}{\displaystyle \frac{1+J_r\mathrm{tanh}(\beta g_t)}{2}})_{\{g_t\},\{J_t\}}]`$
which thanks to the relation between $`Q(u)`$ and $`G(g)`$, rewrites as

$`R_{RS}^{KSAT}[G,t]`$ $`=`$ $`{\displaystyle \frac{\alpha }{\beta }}E[\mathrm{log}(1+(\mathrm{e}^\beta 1)\omega ({\displaystyle \underset{t=1}{\overset{p}{}}}{\displaystyle \frac{1+J_tS_t}{2}}))_{\{J_t\}}`$ (B6)
$`p\mathrm{log}\left(1+\xi \omega \left({\displaystyle \frac{1+JS}{2}}{\displaystyle \underset{t=1}{\overset{p1}{}}}{\displaystyle \frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}}\right)\right)_{\{g_t\},J,\{J_t\}}+`$
$`(p1)\mathrm{log}(1+\xi {\displaystyle \underset{t=1}{\overset{p}{}}}{\displaystyle \frac{1+J_r\mathrm{tanh}(\beta g_t)}{2}})_{\{g_t\},\{J_t\}}]`$
The last term has been added and subtracted from eq.(37) in order to extract a remainder that would vanish if replica symmetry holds, and maximization is performed on $`G(g)`$. As in the $`p`$-spin case, we will proceed in a Taylor expansion of expression (B6) in powers of $`\xi `$, and rely on absolute convergence to average each term of the series.

Expanding the first term in (B6) we can write

$`E\left[\mathrm{log}\left(1+\xi \omega ({\displaystyle \underset{t=1}{\overset{p}{}}}{\displaystyle \frac{1+J_tS_t}{2}})\right)_{\{J_t\}}\right]=`$ (B7)
$`{\displaystyle \underset{n1}{}}{\displaystyle \frac{(1)^{n+1}}{n}}(\xi ^{})^nE\left[\omega \left({\displaystyle \underset{t=1}{\overset{p}{}}}(1+J_tS_t)\right)^n_{\{J_t\}}\right]=`$ (B8)
$`{\displaystyle \underset{n1}{}}{\displaystyle \frac{(1)^{n+1}}{n}}(\xi ^{})^n\mathrm{\Omega }\left[{\displaystyle \underset{t=1}{\overset{p}{}}}\left(1+{\displaystyle \underset{l=1}{\overset{n}{}}}J_t^l_{J_t}{\displaystyle \underset{a_1<\mathrm{}<a_l}{\overset{1,n}{}}}S_t^{a_1}\mathrm{}S_t^{a_l}\right)\right]=`$ (B9)
$`{\displaystyle \underset{n1}{}}{\displaystyle \frac{(1)^{n+1}}{n}}(\xi ^{})^n\mathrm{\Omega }\left[{\displaystyle \underset{t=1}{\overset{p}{}}}\left(1+{\displaystyle \underset{l=1}{\overset{n}{}}}J_t^l_{J_t}{\displaystyle \underset{a_1<\mathrm{}<a_l}{\overset{1,n}{}}}q^{a_1\mathrm{}a_l}\right)\right]=`$ (B10)
$`{\displaystyle \underset{n1}{}}{\displaystyle \frac{(1)^{n+1}}{n}}(\xi ^{})^n\mathrm{\Omega }[(1+Q_n)^p]`$ (B11)
where we have defined $`\xi ^{}(e^\beta 1)/(2^p)`$ and $`_{l=1}^nJ^l_J_{a_1<\mathrm{}<a_l}^{1,n}q^{a_1\mathrm{}a_l}Q_n`$. Notice that due to the negative sign of $`\xi ^{}`$, the coefficients $`(1)^{n+1}(\xi ^{})^n`$ are all negative.

The analogous expansion of the second term is:

$`E\left[\mathrm{log}\left(1+\xi \omega \left({\displaystyle \frac{1+JS}{2}}{\displaystyle \underset{t=1}{\overset{p1}{}}}{\displaystyle \frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}}\right)\right)_{\{J_t\},J,\{g_t\}}\right]=`$ (B12)
$`{\displaystyle \underset{n1}{}}{\displaystyle \frac{(1)^{n+1}}{n}}(\xi ^{})^n\mathrm{\Omega }\left[\left(1+{\displaystyle \underset{l=1}{\overset{n}{}}}J^l_J{\displaystyle \underset{a_1<\mathrm{}<a_l}{\overset{1,n}{}}}q^{a_1\mathrm{}a_l}\right){\displaystyle \underset{t=1}{\overset{p1}{}}}{\displaystyle \underset{l=1}{\overset{n}{}}}\left(1+J_t\mathrm{tanh}(\beta g_t)\right)_{\{J_t\},\{g_t\}}\right]=`$ (B13)
$`{\displaystyle \underset{n1}{}}{\displaystyle \frac{(1)^{n+1}}{n}}(\xi ^{})^n\mathrm{\Omega }\left[(1+Q_n)(1+J\mathrm{tanh}(\beta g))^n_{J,g}^{p1}\right]`$ (B14)
Finally, the third terms in eq.(B6) immediately reads

$$\mathrm{log}\left(1+\xi \underset{t=1}{\overset{p}{}}\frac{1+J_t\mathrm{tanh}(\beta g_t)}{2}\right)_{\{J_t\},\{g_t\}}=\underset{n1}{}\frac{(1)^{n+1}}{n}(\xi ^{})^n(1+J\mathrm{tanh}(\beta g))^n_{J,g}^p$$
(B15)
The sum of the three pieces in eq.(B6) gives:

$$R_{RS}^{KSAT}[G,t]=\frac{\alpha }{\beta }\underset{n1}{}\frac{(1)^n}{n}(\xi ^{})^n\mathrm{\Omega }\left[(1+Q_n)^pp(1+Q_n)(1+J\mathrm{tanh}(\beta g))^n_{J,g}^{p1}+(p1)(1+J\mathrm{tanh}(\beta g))^n_{J,g}^p\right]$$
(B16)
The previous sum is always positive semidefinite for $`p`$ even while we need $`1+Q_n0`$ for $`p`$ odd.

### 2 …and of $`R_{1RSB}^{KSAT}`$

We proceed in the same way as in the $`p`$-spin case. The algebra is elementary but more tedious and involved, therefore we will only list the final results of the calculation. Starting from eq.(77), we again expand in series the first term, getting, with a treatment similar to the RS case:

$`R_{1RSB}^{KSAT}[𝒢,t]={\displaystyle \underset{l1}{}}{\displaystyle \frac{m^l}{l}}{\displaystyle \underset{k_1,\mathrm{},k_l}{\overset{1,\mathrm{}}{}}}(\xi ^{})^{_{s=1}^lk_s}{\displaystyle \underset{s=1}{\overset{l}{}}}\left({\displaystyle \frac{_{r=1}^{k_s1}(rm)}{k_s!}}\right)\mathrm{\Omega }^{(l)}\left[(1+𝐐(k_1,\mathrm{},k_l))^p\right]`$ (B17)
where we have defined:

$`𝐐(k_1,\mathrm{},k_l)`$ $``$ $`{\displaystyle \underset{s=1}{\overset{l}{}}}{\displaystyle \underset{r_1,\mathrm{},r_s}{\overset{k_1,\mathrm{},k_s}{}}}J^{(r_1+\mathrm{}+r_s)}_J{\displaystyle \underset{t=1}{\overset{s}{}}}{\displaystyle \underset{a_1<\mathrm{}<a_{r_t}=1}{\overset{k_1,\mathrm{},k_s}{}}}q^{(a_{r_1},\mathrm{},a_{r_s})}`$ (B18)
Analogous steps give for the second term in eq.(77)

$`{\displaystyle \underset{l1}{}}{\displaystyle \frac{m^l}{l}}{\displaystyle \underset{k_1,\mathrm{},k_l}{\overset{1,\mathrm{}}{}}}(\xi ^{})^{_{s=1}^lk_s}{\displaystyle \underset{s=1}{\overset{l}{}}}\left({\displaystyle \frac{_{r=1}^{k_s1}(rm)}{k_s!}}\right)\mathrm{\Omega }^{(l)}\left[1+𝐐(k_1,\mathrm{},k_l)\right]{\displaystyle \underset{s=1}{\overset{l}{}}}\left(1+J\mathrm{tanh}(\beta g)\right)^{k_l}_g_{G,J}^{p1}`$ (B19)
and for the third term

$`{\displaystyle \underset{l1}{}}{\displaystyle \frac{m^l}{l}}{\displaystyle \underset{k_1,\mathrm{},k_l}{\overset{1,\mathrm{}}{}}}(\xi ^{})^{_{s=1}^lk_s}{\displaystyle \underset{s=1}{\overset{l}{}}}\left({\displaystyle \frac{_{r=1}^{k_s1}(rm)}{k_s!}}\right){\displaystyle \underset{s=1}{\overset{l}{}}}\left(1+J\mathrm{tanh}(\beta g)\right)^{k_l}_g_{G,J}^p,`$ (B20)
where in the last two terms we can further expand

$$\underset{s=1}{\overset{l}{}}\left(1+J\mathrm{tanh}(\beta g)\right)^{k_l}_g_{G,J}^n=\left(\underset{r_1,\mathrm{},r_l=1}{\overset{k_1,\mathrm{},k_l}{}}\underset{s=1}{\overset{l}{}}\left(\genfrac{}{}{0pt}{}{k_s}{r_s}\right)J^{(r_1+\mathrm{}+r_l)}_J\underset{s=1}{\overset{l}{}}\left(\mathrm{tanh}(\beta g)\right)^{r_s}_g_G^n\right)$$
(B21)
with $`n`$ equal to $`p1`$ and $`p`$ respectively. Since $`\xi ^{}<0`$ it is easy to see how only positive terms of the series survive.

Collecting all, we eventually find the complete power expansion for $`R_{1RSB}^{KSAT}`$:

$`{\displaystyle \frac{\alpha }{\beta m}}{\displaystyle \underset{l1}{}}{\displaystyle \frac{m^l}{l}}{\displaystyle \underset{k_1,\mathrm{},k_l}{\overset{1,\mathrm{}}{}}}(\xi ^{})^{_{s=1}^lk_s}{\displaystyle \underset{s=1}{\overset{l}{}}}\left({\displaystyle \frac{_{r=1}^{k_s1}(rm)}{k_s!}}\right)\times `$ (B22)
$`\mathrm{\Omega }^{(l)}\left[(1+𝐐(k_1,\mathrm{},k_l))^pp(1+𝐐(k_1,\mathrm{},k_l))𝐀(k_1,\mathrm{},k_l)^{p1}+(p1)𝐀(k_1,\mathrm{},k_l)^p\right]`$ (B23)
where we have defined

$$𝐀(k_1,\mathrm{},k_l)\underset{s=1}{\overset{l}{}}\left(1+J\mathrm{tanh}(\beta g)\right)^{k_l}_g_G$$
(B24)
Again, every term of the expansion is positive for even $`p`$ and for $`p`$ odd under condition $`1+𝐐(k_1,\mathrm{},k_l)0`$.

## C Appendix C

Let us briefly sketch the proof of the existence of the thermodynamic limit of free-energy of the $`p`$ spin model for $`p`$. Let us define a model which interpolates between two non interacting systems with $`N_1`$ and $`N_2`$ spins respectively, and a system of $`N=N_1+N_2`$ spins. Each clause $`\mu =1,\mathrm{},M`$ will belong to the total system with probability $`t`$, to the first subsystem with probability $`N_1/N(1t)`$ and to the second subsystem with probability $`N_2/N(1t)`$. We chose the indices $`i_1^\mu ,\mathrm{},i_p^\mu `$ in the following way: for each clause the indices will be i.i.d. with probability $`t`$, the indices will be chosen uniformly in the set $`\{1,\mathrm{},N\}`$, with probability $`(1t)N_1/N`$ the indices will be chosen in $`\{1,\mathrm{},N_1\}`$ and with probability $`(1t)N_2/N`$ in the set $`\{N_1+1,\mathrm{},N\}`$. Let us consider the free-energy $`F_N(t)=\frac{1}{N\beta }\mathrm{log}Z(t)`$. A direct calculation of its $`t`$-derivative

$$\frac{dF_N(t)}{dt}=\frac{1}{\beta }\left[\frac{1}{N^p}\underset{i_1,\mathrm{},1_p}{\overset{1,N}{}}+\frac{N_1}{N}\frac{1}{N_1^p}\underset{i_1,\mathrm{},1_p}{\overset{1,N_1}{}}+\frac{N_2}{N}\frac{1}{N_2^p}\underset{i_1,\mathrm{},1_p}{\overset{N_1+1,N}{}}\right]E\mathrm{log}(1+\mathrm{tanh}(\beta J)\omega (S_{i_1}\mathrm{}S_{i_p}))_J.$$
(C1)
Expanding the logarithm in series, observing that thanks to the symmetry of the $`J`$ distribution the odd term vanish, introducing the replica measure and using the convexity of the function $`x^p`$ for even $`p`$ one proves that $`\frac{dF_N(t)}{dt}0`$ which implies sub-additivity $`F_N\frac{N_1}{N}F_{N_1}+\frac{N_2}{N}F_{N_2}`$; this is in turn is a sufficient condition to the existence of the free-energy density. The same prove applies to the even $`p`$ random K-SAT model. For odd $`p`$ we face a difficulty similar to the one in the replica bounds. We can not prove sub-additivity due to the need to consider negative values of the overlaps, and non convexity of $`x^p`$ for negative $`x`$.