text:
# Limits on the infrared photometric monitoring of brown dwarfs

## 1 Introduction

Thermodynamical considerations show that dust particles form in the atmospheres of ultra cool M, L and T dwarfs (e.g. Allard et al. allard01 (2001); Burrows & Sharp burrows99 (1999); Helling helling01 (2001); Lodders lodders99 (1999); Lunine et al. lunine89 (1989); Sharp & Huebner sharp90 (1990); Tsuji et al. tsuji96 (1996)). This has been confirmed in many cases from the modelling of their optical and near infrared spectral energy distributions (e.g. Leggett et al. leggett01 (2001); Marley et al. marley02 (2002); Schweitzer et al. schweitzer01 (2001)schweitzer02 (2002)). Physical models suggest that as the effective temperature drops below about 2500 K, the photosphere becomes increasingly dusty due to increased gas condensation, whereas at even lower temperatures the dust opacity is reduced due to gravitational settling of the dust below the photosphere (Ackerman & Marley ackerman01 (2001); Tsuji tsuji01 (2001)).

However, such models assume uniform horizontal distributions of this dust, whereas recent observational work has shown a significant fraction of ultra cool dwarfs to be photometrically variable at the level of 0.01–0.08 mag (Bailer-Jones & Mundt bjm99 (1999)bjm01 (2001); Clarke et al. 2002a 2002b ; Gelino et al. gelino02 (2002); Martín et al. martin01 (2001); Nakajima et al. nakajima00 (2000); Tinney & Tolley tinney99 (1999)). Bailer-Jones & Mundt bjm01 (2001) found evidence for a rapid evolution of features on the surfaces of a few field L dwarfs. Being field objects, they are presumably part of an older population which have long lost their disks (age $`>10`$ Myr), thus ruling out accretion hot spots or other star–disk activity as a cause of the variability. This leaves magnetic star spots or dust clouds as the most plausible candidates, although as yet there is no direct evidence for the physical nature of these surface features. Theoretical arguments from Gelino et al. gelino02 (2002) and Subhanjoy et al. subhanjoy02 (2002) argue against the existence of spots in ultra cool dwarfs due to the low ionization fractions. On the other hand, a non-monotonic variation with spectral type in the strength of the FeH band at 9896 Å (which would be effected by iron condensation) noted by Burgasser et al. burgasser02 (2002) may be indicative of an inhomogeneous dust cloud coverage, at least around the L/T transition.

If the variability is due to dust, then understanding its physical, chemical and dynamical nature is central to determining the fundamental physical parameters of ultra cool dwarfs ($`\mathrm{T}_{\mathrm{eff}}`$, abundances, ages etc.) and hence their formation mechanisms. Using the dusty and condensed dust models of Allard et al. allard01 (2001), one of us has recently predicted the observational signatures of different dust cloud and spot scenarios (Bailer-Jones bj02 (2002)). This shows that there would be a distinct signature in the J and K bands for these phenomena, depending on effective temperature and feature size.

In this paper, we present the results of a pilot project to look for variability in ultra cool dwarfs through near-simultaneous J and K band monitoring. A comparison of the amplitudes and correlation of any variability can be compared with surface feature scenarios from various models.

## 2 Data acquisition and reduction

Details of the three targets monitored are given in Table 1. They were chosen from a list of sufficiently bright objects with appropriate RA/Dec for the observatory/season to represent a range of L spectral types. The targets were observed over 16 nights in March/April 2001 from the 1.23m telescope at Calar Alto, Spain, with the MAGIC camera, equipped with a 256$`\times `$256 NICMOS3 array. The pixel scale was 1.15′′/pix providing a field-of-view of 5$`{}_{}{}^{}\times `$5. The objects were observed in both the J (50% transmission points at 1.12 and 1.40$`\mu `$m) and K (50% transmission points at 1.93 and 2.27$`\mu `$m) filters. Each object was monitored continuously for several hours each night (alternating between J and K), weather conditions permitting. Table 2 shows the log of observations. The seeing varied between about 1.2′′ and 2.2′′ in J (poorer seeing images are not retained in the final light curves). Weather conditions varied across the run.

A problem with the NICMOS arrays is the variation of the quantum efficiency (QE) across a single pixel (intra-pixel QE variation). On account of the coarse pixel scale and undersampled photometry, this could result in a variation in the detected flux of a non-variable star as its image moves relative to the center of a pixel. For this reason, each epoch (point in the final light curves) is derived from 45 frames taken with a dithering macro with non-integer pixel offsets: A field is observed at each of the nine points in a 3$`\times `$3 square grid, with x/y distance 12′′ (10.4 pixels) between grid points. Five consecutive 5s integrations (repeats) are taken at each grid point, read out in double correlated readout (DCR) mode. The macro required 4.7 minutes to execute. Due to a variable bias level, the first repeat at each grid point could not be used, leaving a total of 36 usable frames. A clipped average of these points provides a single point in the final light curve (see section 3). This procedure is necessary to achieve high precision photometry with infrared arrays.

It was further found that the bias levels for repeats 2–5 were all different, although stable in time. Thus the flat fielding and sky subtraction had to be performed separately for each repeat.

Flat fields were created on each night from the differences between illuminated and non-illuminated dome flats. For each epoch, a sky image was then created from a median combination of the nine frame positions for that repeat stage (each offset to have a common flux zero point to accommodate changes in the overall background). This process removed all sources, leaving an image of the 2D background (predominantly sky brightness) variations which was subtracted from each frame.

No correction was made for nonlinearity in the flux response of the pixels. In principle, a nonlinearity can introduce apparent variability if the stars do not remain exactly in the same place on the detector or if the stars or background undergo a uniform change in brightness (e.g. due to changes in sky brightness or atmospheric extinction). However, we have calculated that such effects introduce errors in the relative photometry of much less than 1% under good conditions, and no more than 1% in the very worst case. Note that bright reference stars were avoided to reduce the influence of nonlinearity.

## 3 Photometry

A reference image was established and the positions of the target and reference stars in this propagated to all other images of the field. This avoids shifting and interpolating images, which can only degrade the data quality. Photometry was carried out using the CCDCAP aperture photometry program written by K. Mighell. CCDCAP was developed to do accurate aperture photometry with small, undersampled apertures, using a bilinear pixel interpolation algorithm (see appendix B of Mighell & Rich mighell95 (1995)). After detailed investigation of the magnitude error (see below) as a function of aperture size, aperture radii of 2.4 and 2.0 pixels were used for J and K respectively. A ‘hardness’ of 1.0 was used in CCDCAP, which corresponds to the maximum subdivision of a pixel. A ‘drift’ of 1.5 pixels in CCDCAP allowed for optimal recentering of the aperture.

Variability monitoring requires a reliable assessment of the photometric errors. Error determinations based on readout noise and Poisson statistics in the source and sky generally underestimate the true error because they ignore various systematic errors. The grid observing method outlined above is ideally suited for determining an empirical error measure based on the 36 frames taken in rapid succession. However, due to numerous bad pixels – up to four of the nine grid positions could be contaminated – simply taking the mean and standard error of these measures frequently gave poor results, because the bad pixels could bias the photometry by several magnitudes. As the PSF is undersampled, interpolation over bad pixels (and cosmic rays) was not deemed reliable (and this can only be done to make reduction convenient: it does not add information). Instead, an iterative rejection scheme was developed to remove bad photometry. First the median and its standard deviation, $`\sigma `$, were determined from all 36 measures. If $`\sigma `$ exceeded a threshold, the highest 8 and lowest 8 measures were rejected, typically corresponding to all measures at four grid positions. The median is then re-established. Iterative clipping of points more than n-$`\sigma `$ about the median was then done to remove outliers (each iteration starting with all 36 measures so that re-inclusion was possible). The initial threshold stage is required because the median of 36 points containing just a few bad measures can result in an unrepresentative median and hence large amounts of data being clipped. A threshold of 0.2 mag and n$`=`$5 were used, although the results are not sensitive to these values. The iterative clipping converged after 0–3 iterations. This procedure was found to be very robust, both excluding obvious bad pixels but resulting in more than 32 (of 36) measures being retained in the majority of cases.

The mean of the remaining $`N`$ measurements was taken to be the final magnitude for that epoch. As these are all independent measurements of the same thing,111This assumes that the star does not vary over 4.7 minutes. If it does, then this error is an overestimate of the true error. the appropriate error, $`\delta m`$, in this measure is the standard error in this mean, i.e. $`\sigma /\sqrt{N}`$. We found that the theoretical errors produced by IRAF’s apphot task (based on source Poisson noise and detector read-out noise) – using the same aperture sizes – were as much as ten times smaller, i.e. grossly underestimated. This demonstrates the importance of using empirical error measures in infrared photometry, at least with NICMOS3 arrays.

Photometry was obtained in this manner for the target and reference stars in all frames. Differential photometry was performed as described in Bailer-Jones & Mundt bjm01 (2001) to produce a light curve for the target relative to the average flux of the reference stars. This removes first order (wavelength independent) variations in the sky extinction between epochs.222This is true provided the extinction is constant across the 5 image field size averaged across the $`70\mathrm{cos}\delta `$ arcmins on the sky through which the field moves during the 4.7 min macro execution time. Light curves were similarly constructed for each reference star relative to all the other reference stars. 6, 8 and 5 reference stars were retained for 2M1112, SDSS1203 and 2M1439 respectively.

In what follows, light curve shall always mean the relative magnitude light curve, i.e. the magnitude of the star relative to its particular set of reference stars. By reference level we mean the time series of the reference magitude for the L dwarf target, i.e. the magnitude formed from the average of the fluxes of the target’s reference stars (see Bailer-Jones & Mundt bjm01 (2001) section 3.3).

## 4 Light curve analysis

The light curves for the target star and reference stars were plotted and examined visually for features. The $`\chi ^2`$ test as used by Bailer-Jones & Mundt bjm01 (2001) was also applied to look for variability across the whole data set and within individual nights. When interpreted as probabilities of variability on the assumption of Gaussian errors, these $`\chi ^2`$ values often indicated significant variability ($`p<<0.01`$) in both the target and some of the reference stars. Assuming that not all of these really are intrinsically variable at this level, this implies one or more of the following: 1. the photometric errors have been underestimated; 2. the errors are non-Gaussian; 3. the relative photometry is not representative of the intrinsic brightness of the stars. Point (1) seems unlikely given the thorough testing of the method described in section 3 to evaluate the errors. Point (2) is generally true due to outliers: this will inflate all $`\chi ^2`$ values for all stars, meaning that one should adopt a conservative threshold for flagging variability (i.e. a small value of the probability, $`p`$, or equivalently a higher $`\chi ^2`$ value). However, this would not be magnitude dependent so would not lead to a differential effect between the target and reference stars. Thus the relative values of $`\chi ^2`$ can still be used as an indication of variability. Point (3) could be problematic on account of colour effects in the extinction variation, as will be discussed in section 6.

As we have near-simultaneous observations in two bands, J and K, we introduce a parameter to look for correlated changes in the relative magnitudes. For a given star at epoch number $`t`$, this parameter is defined as

$$Q(t)=\underset{t^{}=1}{\overset{t^{}=t}{}}\frac{m_\mathrm{J}(t^{})m_\mathrm{J}(t^{}1)}{\delta m_\mathrm{J}}.\frac{m_\mathrm{K}^{}(t^{})m_\mathrm{K}^{}(t^{}1)}{\delta m_\mathrm{K}^{}}$$

where $`m_\mathrm{J}(t^{})`$ is the relative J magnitude of the star at epoch $`t^{}`$ and similarly for $`m_\mathrm{K}^{}(t^{})`$. $`\delta m_\mathrm{J}`$ is the error in $`m_\mathrm{J}(t^{})m_\mathrm{J}(t^{}1)`$, obtained from the quadrature sum of the empirical error measures, $`\delta m`$, for each epoch (see section 3). Each of the two terms on the right hand side of this equation gives the change in relative magnitude from one epoch to the next in the units of the random error; significant changes have a modulus larger than unity. As the two terms, one for J and one for K, are multiplied together, correlated changes in the two bands give a positive contribution to the sum; anticorrelated changes a negative contribution. Thus if the star shows a series of correlated changes, $`Q`$ will become more positive; if they are anticorrelated it becomes more negative. If changes are random, i.e. sometimes correlated and sometimes anticorrelated, $`Q`$ will do a random walk; in particular it has an expectation value of zero and variance of $`t`$ (the number of epochs). By comparing the variation in $`Q`$ for the target with that for its reference stars, we can see whether the target tends to show a greater level of colour correlated (or colour anticorrelated) changes than the reference stars. This is summarized by $`Q_\mathrm{s}`$, the difference between the maximum and minimum values of $`Q`$ across all epochs.

As cool spots and dust features effect the radiative flux over a wide wavelength range in the optical and infrared, the evolution of such features will produce some kind of correlated or anticorrelated change in the J and K bands and hence some pattern in the Q parameter. This is the motivation for the Q parameter. The specific relative amplitudes in these bands depends on the physical mechanism (spots or clouds or something else) and of course on the specific atmospheric models used. For example, making predictions based on the models by Allard et al. allard01 (2001), we see that dust variations in early L dwarfs produce anticorrelated J/K variations, whereas star spots produce correlated ones (Bailer-Jones bj02 (2002)). Different relative amplitudes may be predicted when making different assumptions in the models, e.g. with a different dust grain size distribution or treatment of convection.

Our J and K observations are not strictly simultaneous. However, for SDSS1203 and 2M1439, a strict pairing of K observations immediately after J observations was maintained in the final light curves, and $`Q`$ has been evaluated using this pairing. For various reasons this was not possible with 2M1112, so $`Q`$ values have not been calculated for this field. A summary of the relative photometry, errors and $`Q_\mathrm{s}`$ is given in Table 3.

## 5 Results

SDSS1203 Calculated across the entire data set, SDSS1203 gives $`\chi ^2`$ values of 200 and 306 for J and K respectively. These formally correspond to very small $`p`$ values (a test with 103 degrees of freedom) but, as noted in section 4, we may expect $`\chi ^2`$ to be large due to outliers across this large number of observations. Indeed, the reference stars also give large $`\chi ^2`$ values: ranging from 114 to 158 (for J) and 98 to 155 (for K). Ignoring the formal $`p`$ values and just looking at the relative $`\chi ^2`$ values, we see that SDSS1203 is indeed more variable than any of the reference stars in both filters.

Fig. 1 shows the J and K light curves across the four nights of observation. There are several points which deviate from the mean by more than 3$`\sigma `$, although generally the variability is quite weak ($`<0.05`$ mag deviation). On night 13, there is a dip in K at around AJD 2001.4. A similar but opposite effect is seen in reference star 2 (although not ref. 1) and the egress from this dip is accompanied by a rise in the reference level, so this is probably a telluric effect. The 4$`\sigma `$ jump in J at AJD 2001.51 is not accompanied by any similar feature in the reference stars or reference level so could be intrinsic to SDSS1203, although it is neither large nor that statistically significant. The J dip at AJD 2002.63, on the other hand, is accompanied by a dip in the reference level (but not the reference stars). This could be due to a colour dependent variation in the extinction (see section 6). We also note that SDSS1203 appeared to be brighter in K at the beginning of night 17, an effect which is marginally significant and cannot obviously be reduced to the behaviour of the reference stars or reference level.

The $`Q`$ parameter in Fig. 1 shows no long lasting colour correlated variations. Moreover, the range of $`Q`$, given by $`Q_\mathrm{s}`$, is no larger than for the reference stars (Table 3). The J/K correlation coefficient, $`\rho `$, is slightly larger for SDSS1203 than the other stars ($`\rho `$=0.34 against $`\rho `$=0.07–0.24 for the reference stars), although it is still small.

We conclude that there is no good evidence for variability in SDSS1203 from these data. SDSS1203 was found to be variable in the I-band by Bailer-Jones & Mundt bjm01 (2001), primarily due to an apparent brightening of the object over a duration of one to two hours. As far as we could tell, those data were taken under good observing conditions near culmination of the object (i.e. very small airmass changes) so we do not believe that result was an artifact due to the Earth’s atmosphere.

2M1439 As with SDSS1203, this object shows more significant variability than any of the reference stars for both J and K according to the $`\chi ^2`$ test applied to all data. However, its $`Q`$ values, as well as the correlations between J and K and J and J$``$K are no larger than for some of the reference stars.

It can be seen in Fig. 2 that the K light curve of 2M1439 dims by 0.14 mag over a period of five hours from AJD 1992.4 on night 4. Although no such trend (or opposite trend) is seen in the light curve of the reference stars individually, the opposite trend is seen in the reference level with a scale of 0.4 magnitudes, indicating that all of the reference stars increased in brightness (Fig. 3). This coincides with the field rising – from an airmass of 2.6 to 1.05 – so is presumably just due to the usual decrease in extinction with decreasing airmass. However, while the reference stars brighten by 0.4 mag, 2M1439 brightens by only 0.26 mag. This is possible if the extinction coefficient has a strong gradient across the K band and if 2M1439 has a different flux gradient from the reference stars. In this case, second order (colour) effects in the terrestrial extinction may become important (see section 6). The same effect is seen in the J band (reference level brightens by 0.35 mag, 2M1439 dims relative to this by 0.07 mag).

However, on night 7 we see the same brightening of the reference level due to decreasing airmass (0.65 mag in both bands from airmass 2.8 to 1.05) but without any trend or significant variability in the light curves of either the reference stars or 2M1439 (see right hand panel in Fig. 3). Clearly, different atmospheric conditions must have prevailed on this night than on night 4.

Correlated changes of up to 1 magnitude in J and K are found in the reference level at the beginning of night 6, caused by thin clouds and humidity. Correlated changes with amplitude up to 0.1 mag are also seen in the light curves of the reference stars and 2M1439, indicating that these large extinction variations have not cancelled out in the relative photometry. This is presumably caused by water vapour changes associated with the clouds (see section 6).

In conclusion, most of the features which account for the higher variability in 2M1439 (compared to the reference stars) appear to be related to terrestrial atmospheric phenomena.

2M1439 has previously been monitored in the I-band by Bailer-Jones & Mundt bjm01 (2001) and by Gelino et al. gelino02 (2002). Neither found good evidence for variability above about 0.01 mag, the former from temporally-dense observations (48 observations spread over four nights) and the latter over longer timescales (a few observations per night over 40 nights).

2M1112 This light curve shows some variations which at first sight would be of significance, but on closer inspection appear to be due to variations in the reference level. These are not qualitatively different from the “problematic” variations seen in the other two targets, discussed above.

We performed a frequency-domain search for periodic variations using two methods, the CLEAN and Lomb-Scargle periodograms. These methods yielded consistent results and showed no evidence for periodic variations in SDSS1203 or 2M1439. As implemented, these methods are known to be effective at detecting even low signal-to-noise periodic signals (Lamm et al. lamm03 (2003)).

## 6 Second order extinction effects

The method of relative photometry employed in this program (and all those referred to in section 1) makes the implicit assumption that the Earth’s atmospheric extinction is the same for the target star as for the reference stars. However, the extinction coefficient varies with wavelength, so if the target and reference stars have different spectral energy distributions (SEDs), their wavelength integrated extinctions are also different. This is relevant when monitoring ultra cool dwarfs, as the reference stars will usually be significantly hotter and hence bluer in the optical and near infrared (see Table 3). This would not matter for relative photometry if the extinction did not vary in time or space. But extinction does vary polychromatically with airmass and atmospheric conditions, thus introducing a change in the relative magnitude between stars with different SEDs. Such effects are referred to as second order extinction.

In the visible, the extinction is dominated by molecular (Rayleigh) and aerosol scattering and generally shows a smooth decrease with increasing wavelength. In this case, colour effects are either small or a broad band colour term can be used to correct relative photometry (also as a function of airmass). In the near infrared, however, the extinction is dominated by molecular absorption, particularly of H2O, CO2 and O3, which show much sharper variations with wavelength. Precipitable water vapour and ozone concentrations in particular can vary rapidly in time and space. Thus even under apparently good observing conditions, second order extinction effects can have a significant impact on high precision, differential, broad band infrared photometry. Note that terrestial water clouds and fog are dominated by large ($`>>`$ 1 $`\mu `$m) water droplets, so these cause a scattering of visible and infrared light with only a weak wavelength dependence ($`\lambda ^\alpha `$, where $`0<\alpha <1`$). Nonetheless, some absorption also occurs in the liquid drops themselves which does show rapid changes in some wavelength ranges, although in the optical at least the colour dependence of terrestrial cloud extinction appears to be small (Honeycutt honeycutt71 (1971); Serkowski serkowski70 (1970)). However, where there are clouds there is presumably also a high water vapour column density, with the result that clouds would contribute to second order extinction.

Unfortuntately, we cannot quantify the effects of second order extinction on our data, because we have no independent measure of the atmospheric constituents – in particular the column density of the precipitable water vapour – and we would further require an appropriate atmospheric model matched to the Calar Alto conditions at the time. Moreover we do not know the SEDs of the reference stars: the broad band J$``$K colour is too undersampled to determine the integrated extinction from a water spectrum (Young young89 (1989)).

An idea of the scale of the problem can be found from reference to the models of J and K extinction at Kitt Peak for a cool giant ($`\mathrm{T}_{\mathrm{eff}}`$=4000 K) and Vega ($`\mathrm{T}_{\mathrm{eff}}`$=9650 K) from Manduca & Bell (1979). At low water column density, the differential extinction between the two stars (E$`{}_{4000}{}^{}`$E9650) is 0.049 mag at J and 0.002 mag at K. Increasing the water column density by a factor of 35 changes these to 0.067 mag at J and 0.001 mag at K. Thus increasing the precipitable water density increases the differential extinction (so decreases the apparent brightness of a cool star relative to a hotter star) by approximately 0.02 mag at J and $``$0.001 mag at K.333Derived from Tables IIa and IIb of Manduca & Bell for Johnson filters at an airmass of 1.5. The actual precipitable water column density experienced during our observations could be larger, so it is not implausible that water vapour changes cause the apparent variability we see in the L dwarfs. Contrary to the figures given above, however, we generally see similar amplitude variations in J as in K. But our K filter extends further into the blueward telluric absorption band than does the Johnson K band of Manduca & Bell and these figures are quite sensitive to the band profile. Also, our target (and probably reference) stars have much lower temperatures than 4000 K and 9650 K.

Changes in airmass under stable conditions can also change the relative photometry. For example, Manduca & Bell show that the J band magnitude of a cool giant relative to Vega would be at least 0.01 mag higher when measured at airmass 2.0 than at airmass 1.0. The effect is much smaller in the K band, however, so this may not be cause of the dimming of 2M1439 on night 4 discussed in section 5 (see Fig. 3), as there the dimming is larger in K than J. However, this may again be due to the differences between the K and K bands and the SEDs of the stars. Recall that no such dimming was seen on night 7 for the same airmass change, perhaps indicating different causes of extinction on these nights.

The fact that both SDSS1203 and 2M1439 showed larger $`\chi ^2`$ values for J and K than the reference stars can be attributed to the likelihood that the reference stars have more similar SEDs to each other than to L dwarfs, so that differential extinction effects between any one reference star and all the other reference stars (which form its reference level) are much smaller. Hence the reference stars show less variability.

Second order extinction is less of a problem in the optical ($`\lambda <0.65`$$`\mu `$m), but here most brown dwarfs are far too faint to be monitored. The I band ($``$ 0.78–0.92 $`\mu `$m) is a good compromise. This too is intersected by a water absorption band centered at 0.81 $`\mu `$m, but this band is much weaker than those effecting the J and K bands (e.g. McCord & Clark mccord79 (1979)).

## 7 Conclusions

The prime requirements of any search for variability are (1) a measurement of the scale of the variability intrinsic to the source, and (2) an accurate determination of the observational errors in these measurements. However the analysis is done, a variability detection relies on the former being significantly larger than the latter.

We have developed an observation and data reduction technique which gives an accurate and reliable determination of the photometric errors, arising from the source, background, detector and effect of the data processing (see section 3).

One of the main problems with ground-based monitoring is variations in the extinction (scattering plus absorption) of the Earth’s atmosphere. The standard procedure of differential photometry was used to give a measure of intrinsic variability. Using this, some evidence for variability was found in the observed L dwarfs. However, closer analysis has shown that this measure is probably “contaminated” by second order extinction effects in the Earth’s atmosphere (see section 6). Although we have no direct evidence for this, it seems the most plausible cause of the observed variability given the features we have described and the elimination of other potential causes. We therefore conclude that we have no good evidence for intrinsic variability in any of the three L dwarfs monitored. None of these objects showed significant correlations in their J and K light curves, as evidenced by correlation plots and our $`Q`$ parameter. Some such (anti)correlation would be expected if the variability is caused by cool spots or dust features. Upper limits on the peak-to-peak amplitude of persistent variability are set by the scatter in the light curves under the most stable Earth atmospheric conditions. For both J and K, these limits are set at 0.04 mag for SDSS1203 and 2M1439 and 0.08 for 2M1112.

As late M, L and T dwarfs are at their brightest in the near infrared (0.9–2.5 $`\mu `$m) and show spectral signatures which could discriminate between various surface feature scenarios (Bailer-Jones bj02 (2002)), this is a desirable wavelength range for variability monitoring. We have shown that high precision differential photometry (random errors of less than two percent) can be achieved on L dwarfs with infrared arrays. However, as a measure of the intrinsic variability in L dwarfs, broad band differential infrared photometry appears to be limited in accuracy to a few percent by variations in terrestrial molecular extinction, in particular precipitable water vapour. These second order extinction effects could be reduced in one or more ways:

1. use of passbands which avoid molecular extinction, in particular use of better designed (narrower) J and K filters to avoid the strong H2O bands which generally intersect these “standard” filters444There is significant variance between the profiles of filters called I, Z, J, H and K at different observatories, so there is no such thing as a “standard” infrared filter set, although most Z, J and K filters are intersected by significant water absorption bands.;
2. use of narrower band filters, which reduces the effect of second order extinction555A variation on this method has been used by Bailer-Jones bj02 (2002), in which an L dwarf was monitored spectrophotometrically relative to another star in the spectrograph slit. The relative photometry was then performed independently in each narrow wavelength bin (0.002 $`\mu `$m), reducing effects of second order extinction to well below other error sources. (this would obviously require a larger telescope or longer integration times; the latter will be unacceptable for detecting short-term variability);
3. use of specially designed filters to monitor the dominant molecular absorption (in particular H2O) from which the broad band relative photometry could, in principle, be corrected (Angione angione87 (1987)). This could be done using a multi-channel camera. However, additional spectroscopic observations are required to calibrate the measures and so apply the correction. Nonetheless, some kind of time-dependent monitoring of extinction would at least allow one to recognise when observing conditions were variable;
4. observing from very dry sites under atmospherically stable conditions. Space-based observations remove all problems related to second order extinction. A strong discriminant between candidate surface features in brown dwarfs is the variability signatures in the water absorption bands at 1.35–1.45$`\mu `$m and 1.80–2.00$`\mu `$m (e.g. Bailer-Jones bj02 (2002). These bands cannot be reliably monitored from the ground, providing another argument for space-based observations.

## Acknowledgements

We are grateful to Ken Mighell for use of his CCDCAP algorithm and for useful information on its application and interpretation.


format: text

text:
# On the formal equivalence of the TAP and thermodynamic methods in the SK model

## I Introduction

The static properties of mean-field spin-glasses have been investigated in the past mainly by means of two different approaches: standard thermodynamics, which through the replica method ea ; sk ; sk2 ; 1rsb ; 2rsb  aims to compute the partition function and the equilibrium free energy density of the system, and the Thouless-Anderson-Palmer (TAP) method tap ; bm ; bm1 ; bm2 ; dd1 ; innocent2 ; ddy , which introduces a mean-field free energy $`F_{TAP}`$, function of the set $`m`$ of local magnetizations, $`m\{m_i\}_{i=1}^N`$. The local minima of $`F_{TAP}(m)`$ are identified with the metastable states of the system, which at the mean-field level are well defined. At low temperatures the total number $`𝒩`$ of minima of the TAP free energy becomes exponentially large bm ; bm1 ; dd1 , and the density of states is given by,

$$\rho (f)e^{N\mathrm{\Sigma }(f)}.$$
(1)
In this expression $`\mathrm{\Sigma }(f)`$ is the complexity of the TAP states with free energy density between $`f`$ and $`f+df`$.

The existence of two apparently different methods for the study of the static properties of mean-field spin-glasses clearly poses a problem of consistency. Both these methods are in principle correct, and the equivalence between the two approaches should hold in general. A direct argument in its favour is the cavity approach MPV ; cavity : on the one hand the cavity approach is known to be mathematically equivalent to the replica approach; on the other hand the first step of the cavity approach consists in clustering the configurations into states which are nothing but the TAP states. From the cavity analysis, one can expect a priori the existence of a formal relationship between the replica results computed at a given order $`k`$ of replica symmetry breaking, and the results of TAP solutions, correctly weighted and with an average computed at the order $`k1`$. Unfortunately, when considering the detailed results published so far, there seemed to exist a discrepancy between the two approaches. For instance the number of TAP solutions computed directly at the order $`k=0`$ (annealed average) in bm , disagrees with the one obtained indirectly through the replica method with one step of replica symmetry breaking (i.e. a $`k=1`$ computation) cavity . This apparent discrepancy has led us to give a careful look at this consistency problem. We show here that the TAP approach must be reconsidered: using a supersymmetry, one finds that the two approaches are in fact fully consistent.

From a general point of view, the problem of consistency between the TAP approach and standard thermodynamics can be posed in a twofold way. First, the two methods must agree in the calculation of the thermodynamic quantities that can be computed in both frameworks. This requirement is a strong form of consistency. At a weaker level, regarding objects which are inherently defined only within the TAP method, we may still expect to find some consistency, and even some formal connections, with the standard thermodynamic approach.

The strong consistency basically asks that the equilibrium free energy computed within the two approaches must be the same. Given that $`F_{TAP}`$ is a function of the local magnetization, it is not obvious how to compute the partition function in this framework. More specifically, due to the large degeneracy of TAP minima, there is the problem of how to weight correctly these states. This problem was solved by De Dominicis and Young ddy  (DDY), who proposed to weight each TAP state $`\alpha `$ with,

$$w_\alpha =\mathrm{exp}[\beta F_{TAP}(m^\alpha )]/Z_{TAP}.$$
(2)
The strong equivalence of the TAP approach to the standard thermodynamic approach is therefore encoded in the equality (which should hold in the thermodynamic limit) $`\frac{1}{N}\mathrm{log}Z=\frac{1}{N}\mathrm{log}Z_{TAP}`$, that is,

$$\frac{1}{N}\mathrm{log}\left(\underset{\sigma }{}\mathrm{exp}[\beta H(\sigma )]\right)=\frac{1}{N}\mathrm{log}\left(\underset{\alpha =1}{\overset{𝒩}{}}\mathrm{exp}[\beta F_{TAP}(m^\alpha )]\right).$$
(3)
Proving relation (3) in general is highly nontrivial, especially because in spin-glasses we can only compare averages over the disorder, and thus, from a practical point of view, we need to compare the average of the replicated partition functions, which will typically require some tricky integration over complicated order parameters. This program has been carried out by DDY for the Sherrington-Kirkpatrick (SK) model sk . Their conclusion is that relation (3) is indeed verified, provided that the complexity of the dominant TAP states at any temperature is zero. This restrictive hypothesis is satisfied in the SK model, but not in other models where the TAP approach is nevertheless used. Thus, it may seem that either condition (3) is not broadly valid, or that a general proof of it must give up the hypothesis used by DDY.

In the present paper we show that the proof of (3) given by DDY for the SK model is valid with no extra hypothesis. The key point is that the TAP partition function can be written in an integral form, with an action which is invariant under the Becchi-Rouet-Stora-Tyutin (BRST) supersymmetry brs ; tito . This symmetry provides the mathematical relations needed to prove equation (3), without the need to invoke any further hypothesis on the complexity of the states.

Concerning the weak consistency of the two methods, the key object is the complexity $`\mathrm{\Sigma }(f)`$. Even though there is no obvious way to compute $`\mathrm{\Sigma }`$ in the standard thermodynamic framework, we may still impose an important consistency condition: at low temperatures (below the static transition) we expect the equilibrium thermodynamic states to have the same free energy density as the lowest TAP states. In other words, we require that at low enough temperatures there must be no static contribution of the metastable states. This condition is encoded in the equation,

$$f_{eq}=f_0,$$
(4)
where $`f_{eq}`$ is the equilibrium free energy density computed in the standard thermodynamic framework, and $`f_0`$ is defined by the TAP relation $`\mathrm{\Sigma }(f_0)=0`$. Apart from this minimal requirement, one may ask whether in the TAP calculation of $`\mathrm{\Sigma }`$ emerge some deeper formal connections with the standard calculation of the free energy. This fact, of course, seems more than likely if the strong consistency condition (3) holds, even though it is in general not obvious. In the $`p`$-spin spherical model crisa1 , relation (4) is satisfied kpz ; crisatap ; monasson ; franzparisi . Moreover, for this same model it has been shown in juanpe  that by means of the BRST supersymmetry the saddle point equations involved in the TAP calculation of $`\mathrm{\Sigma }`$ at $`f=f_0`$ become identical to the static equations in the replica approach. In the SK model, however, relation (4) is not easy to check bm ; bm1 ; bm2 ; dd1 ; innocent2 ; potters , since it requires a full replica symmetry breaking (RSB) calculation of the complexity. Moreover, up to now there was no formal connection between the TAP complexity and equilibrium free energy in the SK model.

The fact that the BRST supersymmetry plays a crucial role in the $`p`$-spin, suggests that even in the SK model this symmetry may help to prove (4), and discover other formal connections between the two methods. This is what we show in Section IV, where we revisit the annealed calculation of the SK complexity of Bray and Moore bm  (BM), and thanks to the BRST supersymmetry show that at $`f=f_0`$ this calculation becomes equal to the standard static calculation of the free energy, at one step level of replica symmetry breaking. This result comes hardly as a surprise, once the strong consistency (3) of the two methods is proved. However, there is a subtle point: it is not clear a priori whether a formal connection between static and TAP approaches is valid only for the correct full RSB solution of the SK model, or if such a connection is preserved at each finite (though approximated) level of replica symmetry breaking. Our result shows that this second scenario is the correct one.

A further, highly nontrivial consistency check of the static and TAP approaches, comes from a comparison of the TAP complexity with the complexity computed by means of constrained thermodynamics, which results in the Legendre transform construction of monasson . Within this method the complexity is given by the Legendre transform of the thermodynamic free energy of the system, thus establishing a connection between TAP and static methods at a generic value of $`f`$. We will prove that, as in the $`p`$-spin spherical model, in the SK model the TAP complexity, once the BRST relations are considered, coincides with the complexity of monasson .

In Section II we introduce the BRST supersymmetry in the context of the TAP approach, and derive the BRST relations that will be used in the rest of the paper. In Section II we show how the BRST relations enforce the DDY proof of (3), without the need of any extra hypothesis. In Section IV we present the BM calculation of the SK complexity, exploiting the BRST supersymmetry, and in Section V we show the connection between our results and the Legendre transform method for the computation of the complexity. Finally, we draw our conclusion in Section VI.

## II The BRST supersymmetry in the TAP context

In this Section we show why the BRST supersymmetry is helpful in the context of the TAP approach. This is just a specific example of the application of the supersymmetric formalism to statistical mechanics, and in particular to the field of disordered systems ps1 ; ps2 ; zinnjustin ; jorge1 ; jorge2 . As we have seen, the TAP free energy $`F_{TAP}`$ is a function of the local magnetizations $`m_i`$, and the minima of this function are identified with the metastable states of the system (TAP states). In the rhs of equation (3) we have a typical example of a sum over different TAP states, as often needed in this approach. In addition to the partition function $`Z_{TAP}`$, we may also compute the density of TAP states (1),

$$\rho (f)=\underset{\alpha =1}{\overset{𝒩}{}}\delta [F_{TAP}(m^\alpha )Nf],$$
(5)
which is the quantity one needs for the computation of the complexity $`\mathrm{\Sigma }(f)`$. In (5), as in (3), TAP states are labeled by $`\alpha `$, and $`m^\alpha `$ indicates the corresponding set of local magnetizations. More generally, in the TAP context we always have to deal with expressions of the form,

$$R=\underset{\alpha =1}{\overset{𝒩}{}}r[F_{TAP}(m^\alpha )],$$
(6)
$`r`$ being a generic function of $`F_{TAP}`$. The quantity $`R`$ can be written as,

$$R=\underset{\alpha =1}{\overset{𝒩}{}}\underset{i}{}dm_i\delta (m_im_i^\alpha )r[F_{TAP}(m)]=\underset{i}{}dm_i\delta (_iF_{TAP}(m))|det(_i_jF_{TAP}(m))|r[F_{TAP}(m)],$$
(7)
where the TAP states have been identified with solutions of the TAP equations $`_iF_{TAP}(m)=0`$. In (7) the modulus of the determinant is quite hard to handle, and due to this it is disregarded in most supersymmetric calculations discuss . This approximation is a priori unjustified, since without the modulus we are weighting each TAP solution with the sign of the Hessian determinant, with the risk of uncontrolled cancellations. The situation may slightly improve if the function $`r(F_{TAP})`$ is peaked only on low values of the free energy $`F_{TAP}`$, as it happens in the calculation of $`Z_{TAP}`$ at low temperatures, and of $`\rho (f)`$ at low free energies. In these cases, we may hope that at low temperatures (free energies) only minima of $`F_{TAP}`$ will dominate the integral in (7). Minima have positive-defined Hessian, and thus the modulus becomes redundant. Moreover, the very identification of TAP states with TAP solution is only sensible if these solutions are minima, rather than generic saddles. Yet, the validity of this scenario should be checked carefully, as it can be done for the $`p`$-spin spherical model, where the effect of disregarding the modulus is completely under control, due to the fact that the free energy distribution of minima and saddles is completely known mod . In the SK model we do not have such a precise information about the TAP free energy saddle points, and therefore dropping the modulus is quite risky. However, our aim here is to show that all the classic TAP calculations (where the modulus was disregarded) are formally consistent with the standard static approach, thus justifying a posteriori this approximation.

We can use an exponential representation for the delta function and the determinant,

$`{\displaystyle \underset{i}{}}\delta (_iF_{TAP})`$ $`=`$ $`{\displaystyle _i\mathrm{}^{+i\mathrm{}}}{\displaystyle \underset{i}{}}{\displaystyle \frac{dx_i}{2\pi i}}\mathrm{exp}\left({\displaystyle \underset{i}{}}x_i_iF_{TAP}(m)\right)`$ (8)
$`det(_i_jF_{TAP})`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \underset{i}{}}d\overline{\psi }_id\psi _i\mathrm{exp}\left[{\displaystyle \underset{ij}{}}\overline{\psi }_i\psi _j_i_jF_{TAP}(m)\right],`$ (9)
where $`\{\overline{\psi },\psi \}`$ are anti-commuting Grassmann variables. In this way we can write zinnjustin ,

$$R=𝒟m𝒟x𝒟\overline{\psi }𝒟\psi e^{S(m,x,\overline{\psi },\psi )},$$
(10)
where the action $`S`$ is given by,

$$S(m,x,\overline{\psi },\psi )=\underset{i}{}x_i_iF_{TAP}(m)+\underset{ij}{}\overline{\psi }_i\psi _j_i_jF_{TAP}(m)+w[F_{TAP}(m)].$$
(11)
In the case $`r`$ is an ordinary function we have $`w=\mathrm{log}r`$, while if $`r(F_{TAP})=\delta (F_{TAP})`$ then $`w=uF_{TAP}`$, where $`u`$ is an imaginary integration variable to implement the $`\delta `$-function (see next Section). The measures in (10) include the sum over the indices and the constant prefactors. A key property of action (11) is its invariance under a generalization of the Becchi-Rouet-Stora-Tyutin (BRST) supersymmetry brs ; tito  (see also zinnjustin ): if $`ϵ`$ is an infinitesimal Grassmann parameter, it is straightforward to verify that (11) is invariant under the following transformation,

$$\delta m_i=ϵ\psi _i\delta x_i=ϵw^{}\psi _i\delta \overline{\psi }_i=ϵx_i\delta \psi _i=0\delta S=0.$$
(12)
This generalization of the standard BRST supersymmetry to the case where $`w0`$ has been first introduced in juanpe , in the context of a TAP calculation for the $`p`$-spin spherical model. The BRST invariance does not depend on the explicit form of the function $`F_{TAP}(m)`$, but simply on the formal structure of action (11). In particular, it is essential that $`w`$ is a function of the magnetizations $`m_i`$ only through $`F_{TAP}(m)`$. The fact that $`\delta S=0`$ under the BRST supersymmetry implies that the average of any observable of the same variables performed with this action must be invariant too. For an observable $`𝒪`$,

$$𝒪(\mathrm{\Lambda })=𝒟\mathrm{\Lambda }𝒪(\mathrm{\Lambda })e^{S(\mathrm{\Lambda })},$$
(13)
with $`\mathrm{\Lambda }=\{m,x,\overline{\psi },\psi \}`$, we have

$$𝒪(\mathrm{\Lambda })𝒪(\mathrm{\Lambda }\delta \mathrm{\Lambda })=\delta 𝒪(\mathrm{\Lambda })=0.$$
(14)
This property can be used to generate some useful Ward identities. A fruitful choice is $`𝒪=m_i^k\overline{\psi }_i`$ and $`𝒪=x_i^k\overline{\psi }_i`$, whose variation gives,

$`km_i^{k1}\overline{\psi }_i\psi _i`$ $`=`$ $`m_i^kx_i`$ (15)
$`kx_i^{k1}w^{}\overline{\psi }_i\psi _i`$ $`=`$ $`x_i^kx_i.`$ (16)
In particular, the case $`k=1`$ gives the BRST equations already used in juanpe ,

$`\overline{\psi }_i\psi _i`$ $`=`$ $`m_ix_i[\mathrm{BRST1}]`$ (17)
$`w^{}\overline{\psi }_i\psi _i`$ $`=`$ $`x_ix_i[\mathrm{BRST2}].`$ (18)
The BRST relations are crucial. Thanks to them we will be able to prove in general that any TAP average must give the same result as the equivalent thermodynamic average, and more specifically to reduce the TAP calculation of the complexity to the standard replica calculation of the free energy.

## III The De Dominicis-Young calculation revisited

In 1983, De Dominicis and Young ddy  (DDY) introduced the statistical weight (2) for the TAP states and explicitly proved that under a key hypothesis the average of the replicated TAP partition function is equal to the average of the replicated standard partition function, i.e.

$$𝒟JP(J)Z^n(J)=𝒟JP(J)Z_{TAP}^n(J),$$
(19)
which is the practical way to write equation (3). This result is very important, because it proves the identity of the TAP and static approaches before any Ansatz for the overlap matrix is done. The hypothesis that DDY invoked to justify their computation is that the complexity of the dominant TAP states at any temperature is zero. More precisely, we can write,

$$Z_{TAP}=𝑑f\rho (f)\mathrm{exp}[\beta Nf]=𝑑f\mathrm{exp}[\beta N(fT\mathrm{\Sigma }(f))]=\mathrm{exp}[\beta N(f^{}T\mathrm{\Sigma }(f^{}))],$$
(20)
where the dominant free energy $`f^{}`$ is solution of the equation,

$$\beta =\frac{\mathrm{\Sigma }(f^{})}{f}.$$
(21)
DDY assumed as a necessary condition that

$$\mathrm{\Sigma }(f^{})=0.$$
(22)
In the present Section we argue that, although relation (22) is in fact verified in the SK model, it is not a necessary condition for the self-consistency of the TAP approach and for its equivalence with the statics.

### III.1 The problem and its solution in the DDY formulation

Let us briefly summarize the arguments of DDY leading to condition (22) (in order to keep the notation as clear as possible we treat $`m`$ as a one-dimensional variable in this Section; things do not change with the full $`N`$-dimensional representation). First, we introduce an auxiliary magnetic field $`h`$ both in $`F_{TAP}(m)`$ through the term $`hm`$, and in the Hamiltonian $`H(\sigma )`$ with $`h\sigma `$. From the very definition of equilibrium magnetization and energy we trivially have,

$`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{Z}{h}}`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}\sigma \mathrm{exp}[\beta H(\sigma )]=\sigma Z,`$ (23)
$`{\displaystyle \frac{Z}{\beta }}`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}H(\sigma )\mathrm{exp}[\beta H(\sigma )]=HZ.`$ (24)
On the other hand, from (3) we get,

$`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{Z_{TAP}}{h}}`$ $`=`$ $`X{\displaystyle \underset{\alpha }{\overset{𝒩(h,\beta )}{}}}{\displaystyle \frac{F_{TAP}(m_\alpha )}{h}}\mathrm{exp}[\beta F_{TAP}(m_\alpha )]=X+\sigma Z`$ (25)
$`{\displaystyle \frac{Z_{TAP}}{\beta }}`$ $`=`$ $`Y+{\displaystyle \underset{\alpha }{\overset{𝒩(h,\beta )}{}}}{\displaystyle \frac{F_{TAP}(m_\alpha )}{\beta }}\mathrm{exp}[\beta F_{TAP}(m_\alpha )]=Y+HZ.`$ (26)
To obtain this result we used the following relation,

$$\frac{F_{TAP}(m_\alpha )}{h}=\underset{\sigma \alpha }{}\sigma \mathrm{exp}[\beta H(\sigma )],$$
(27)
where the sum over $`\sigma `$ is restricted to those configurations belonging to state $`\alpha `$ (an analogous relation is valid for $`F_{TAP}/\beta `$). On the other hand, in equations (25)-(26) $`X`$ and $`Y`$ are the contributions coming, respectively, from the dependence on $`h`$ and $`\beta `$ of the total number of TAP states, $`𝒩(h,\beta )`$. DDY correctly noted that $`X=0`$ and $`Y=0`$ are necessary conditions for the consistency of the TAP approach. Moreover, by using (20), they noted that if $`\mathrm{\Sigma }(f^{})=0`$, these extra contributions are indeed vanishing in the thermodynamic limit, and thus consistency is recovered. Thus, DDY assumed $`\mathrm{\Sigma }(f^{})=0`$ and accordingly imposed $`X=Y=0`$ in their computation. More precisely, after averaging over the disorder, the relations $`X=Y=0`$ become two equations for the auxiliary fields (eqs.(32) and (42) of ddy ), which DDY used to prove (3). Our point, however, is that the equations $`X=Y=0`$ are always satisfied, irrespective of the value of $`\mathrm{\Sigma }(f^{})`$.

### III.2 Role of the supersymmetry

We now show that the mathematical conditions on the auxiliary variables that DDY impose are nothing else than the BRST relations, which are valid in general, and contain no information on $`\mathrm{\Sigma }(f^{})`$. Let us define,

$$G(m,h)=\frac{(\beta F_{TAP})}{m},r(m,h)=e^{\beta F_{TAP}(m,h)}.$$
(28)
We consider relation (23) (a similar reasoning can be easily done also for relation (24)) and use the integral representation of Section II for $`Z_{TAP}`$, to obtain,

$$\frac{Z_{TAP}}{h}=𝑑m𝑑xxe^{xG(m,h)}\frac{G}{m}r(m,h)+𝑑m𝑑xe^{xG(m,h)}\frac{G}{m}\frac{r}{h}.$$
(29)
In this expression the first term is $`X`$, whereas the second term is just $`\sigma Z`$. Integrating by parts we have,

$$X=𝑑m\delta (G(m,h))\frac{r}{m}.$$
(30)
At this point note that for a generic function $`r(m,h)`$ there is no need for $`X`$ to be zero, in much the same way as action (11) is not BRST invariant if $`r`$ is not a function of $`F_{TAP}`$. If, however, $`r(h,m)=r(F_{TAP}(m,h))`$, as in the case under examination, we have,

$$X=𝑑m\delta (G(m,h))G(m,h)\frac{r}{F_{TAP}}=0.$$
(31)
Thus, the extra term coming from the differentiation of $`𝒩(h,\beta )`$ with respect to $`h`$ vanishes if the TAP states are weighted with a function of the TAP free energy. A similar conclusion can be easily drawn for differentiation with respect to $`\beta `$, which leads to $`Y=0`$. In particular, these relations hold in the case $`r(F_{TAP})=\mathrm{exp}(\beta F_{TAP})`$, which is the one analyzed by DDY.

The important fact is that from a mathematical point of view, $`X=0,Y=0`$ are consequences of the particular form of action (11), and of its symmetry properties. Indeed, if we express (23) and (24) in the supersymmetric representation (13), we get,

$`X`$ $`=`$ $`x`$
$`Y`$ $`=`$ $`\left[x+\overline{\psi }\psi {\displaystyle \frac{}{m}}\right]{\displaystyle \frac{G}{\beta }}={\displaystyle \underset{k}{}}c_km^kx+km^{k1}\overline{\psi }\psi ,`$ (32)
with,

$$c_k=\frac{1}{k!}\frac{^k}{m^k}\frac{G}{\beta }.$$
(33)
A comparison between (32) and (15) clearly shows that the relations $`X=0`$ and $`Y=0`$ are a direct consequence of the BRST symmetry. Thus, the relations imposed by DDY in the calculation of $`Z_{TAP}^n`$, and in their demonstration of the equivalence between TAP and static averages, are a consequence of the BRST supersymmetric form of the TAP action (11). These relations are therefore not the expression for the absence of an extensive complexity of the equilibrium states, i.e. they do not imply $`\mathrm{\Sigma }(f^{})=0`$. For this reason, the result of DDY is more general than what originally thought. This is consistent with the fact that a similar formal connection between TAP and static approaches is valid in the $`p`$-spin spherical model, where, in a certain range of temperature, $`\mathrm{\Sigma }(f^{})>0`$.

### III.3 Disregarding the modulus

Before closing this Section, a word of caution on the modulus of the determinant ps2 ; jorge1 ; mod . As already noted above, disregarding the modulus can be very risky, especially if the function $`r(F_{TAP})`$ weighting the TAP states in (6) is not peaked on low free energies, and thus minima. As an extreme illustration of this risk we consider the function $`r(F_{TAP})=1`$. Clearly, from (6) we have,

$$R(h)=𝒩(h)\frac{dR}{dh}=\frac{d𝒩}{dh}0.$$
(34)
However, if we use for $`R`$ the integral representation (7) and we disregard the modulus, we obtain,

$$\frac{dR}{dh}=0,$$
(35)
since the BRST supersymmetry is trivially satisfied with $`r(F_{TAP})=1`$. The problem here is that by disregarding the modulus, and with a flat weight $`r`$, we are summing over all stationary points of $`F_{TAP}`$, each multiplied by the sign of the determinant. The Morse theorem states that this quantity must be a topological constant, only dependent on the manifold over which $`F_{TAP}`$ is define and on the boundary conditions on $`F_{TAP}`$. Result (35) is therefore correct, but the quantity in this equation is not the same $`R(h)`$ as in (34). This trivial example shows how important is to weight TAP solutions with a function peaked as much as possible on low free energies when the modulus is disregarded.

## IV The Bray-Moore calculation revisited

In the present Section we will compute the annealed complexity of the TAP states for the SK model, following closely the classic calculation of Bray and Moore (BM) bm . Our new contribution will be to exploit the BRST relations (17) and (18) in order to simplify the resulting saddle point equations. In this way we will prove that the complexity is intimately connected to the 1RSB static free energy.

### IV.1 General definitions

The TAP free energy for the SK model is given by tap ,

$$F_{TAP}(m)=\frac{1}{2}\underset{ij}{}J_{ij}m_im_j+\frac{1}{\beta }\underset{i}{}\varphi _0(q,m_i),$$
(36)
with,

$`\varphi _0(q,m)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+m)\mathrm{log}\left[{\displaystyle \frac{1}{2}}(1+m)\right]+{\displaystyle \frac{1}{2}}(1m)\mathrm{log}\left[{\displaystyle \frac{1}{2}}(1m)\right]{\displaystyle \frac{\beta ^2}{4}}(1q)^2`$ (37)
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{log}(1m^2)+m\mathrm{tanh}^1(m)\mathrm{log}2{\displaystyle \frac{\beta ^2}{4}}(1q)^2.`$
The variables $`m_i`$ are the local magnetizations, and $`q`$ is the self-overlap of the TAP states,

$$q=\frac{1}{N}\underset{i}{}m_i^2,$$
(38)
while the quenched couplings $`J`$ are random variables with Gaussian distribution,

$$P(J_{ij})=\sqrt{N/2\pi }\mathrm{exp}(NJ_{ij}^2/2).$$

The TAP equations and the Hessian of the free energy are respectively,

$`\beta _iF_{TAP}(m)`$ $`=`$ $`\beta {\displaystyle \underset{ji}{}}J_{ij}m_j+\varphi _1(q,m_i)=0,`$
$`\beta _i_jF_{TAP}(m)`$ $`=`$ $`\beta J_{ij}+\varphi _2(q,m_i)\delta _{ij},`$
with,

$`\varphi _1(q,m)`$ $`=`$ $`\beta ^2(1q)m+\mathrm{tanh}^1(m)`$
$`\varphi _2(q,m)`$ $`=`$ $`\beta ^2(1q)+{\displaystyle \frac{1}{1m^2}}+O(1/N).`$ (39)
The term of order $`1/N`$ in $`\varphi _2(q,m)`$ will be dropped in what follows. Following BM we perform an annealed calculation of the number of TAP states, i.e. we directly average $`\rho (f,\beta |J)`$ in (5) over the distribution of the quenched couplings $`J_{ij}`$,

$$𝒩(\beta ,f)=𝒟JP(J)\rho (f,\beta |J).$$

We use for $`\rho `$ the integral representation of eqs. (10) and (11), with

$$\delta (F_{TAP}Nf)=_i\mathrm{}^{+i\mathrm{}}\frac{du}{2\pi i}\mathrm{exp}[u(F_{TAP}Nf)],$$
(40)
and thus $`w(F_{TAP})=u(F_{TAP}Nf)`$ and $`w^{}=u`$. Thus, the average number of TAP states becomes,

$$𝒩(\beta ,f)=𝒟JP(J)𝒟m𝒟x𝒟\overline{\psi }𝒟\psi 𝑑ue^{\beta S(m,x,\overline{\psi },\psi ,u)},$$
(41)
where we have multiplied the action by $`\beta `$ in order to keep our calculation as close as possible to BM. Before proceeding we mention an important point: in the calculation of BM the $`J`$-dependent part of $`F_{TAP}(m)`$ in $`S`$ is eliminated by using the equations $`_iF_{TAP}(m)=0`$, which are enforced by the $`\delta `$-function. More specifically, BM use the equation,

$$\frac{1}{2}\underset{ij}{}J_{ij}m_im_j=\frac{1}{2\beta }\underset{i}{}m_i\varphi _1(q,m_i),$$
(42)
which is valid in the TAP states. This substitution simplifies considerably the calculation, but unfortunately the action obtained in this way is no longer BRST invariant. Thus, in the present calculation we must use the full form of $`F_{TAP}(m)`$, equation (36), and for this reason some parts of the calculation differ from BM.

### IV.2 The calculation

The $`J`$-dependent part of the action in (41) is given by,

$$\beta S_J=\beta \underset{ij}{}J_{ij}(x_im_j+\overline{\psi }_i\psi _j+\frac{1}{2}um_im_j),$$
(43)
and after averaging $`𝒩`$ over the disorder we obtain the new effective action,

$`\beta S`$ $`=`$ $`{\displaystyle \underset{i}{}}x_i\varphi _1(q,m_i)+{\displaystyle \underset{i}{}}\overline{\psi }_i\psi _i\varphi _2(q,m_i)+u{\displaystyle \underset{i}{}}\varphi _0(q,m_i)N\beta uf`$ (44)
$`+`$ $`{\displaystyle \frac{\beta ^2q}{2}}{\displaystyle \underset{i}{}}x_i^2+{\displaystyle \frac{\beta ^2}{2N}}{\displaystyle \underset{i}{}}(m_ix_i)^2{\displaystyle \frac{\beta ^2}{2N}}{\displaystyle \underset{i}{}}(\overline{\psi }_i\psi _i)^2+N{\displaystyle \frac{\beta ^2}{4}}u^2q^2+\beta ^2uq{\displaystyle \underset{i}{}}m_ix_i.`$
In order to linearize the quadratic terms, we introduce in (41) the following $`\delta `$-functions,

$`\delta (qN{\displaystyle \underset{i}{}}m_i^2)`$ $`=`$ $`{\displaystyle _i\mathrm{}^{+i\mathrm{}}}𝑑\lambda e^{\lambda qN+\lambda _im_i^2}`$ (45)
$`\delta (RN{\displaystyle \underset{i}{}}m_ix_i)`$ $`=`$ $`{\displaystyle _i\mathrm{}^{+i\mathrm{}}}𝑑re^{rRN+r_im_ix_i}`$ (46)
$`\delta (TN{\displaystyle \underset{i}{}}\overline{\psi }_i\psi _i)`$ $`=`$ $`{\displaystyle _i\mathrm{}^{+i\mathrm{}}}𝑑te^{tTN+t_i\overline{\psi }_i\psi _i},`$ (47)
and we integrate over $`q,R`$ and $`T`$. In this way the integrals in $`x_i`$ and $`(\overline{\psi }_i,\psi _i)`$ become Gaussian and can be performed explicitly. The effective action becomes,

$$\beta S(\mathrm{\Omega },m)=N\mathrm{\Sigma }_0(\mathrm{\Omega })+\underset{i}{}(\mathrm{\Omega },m_i),$$
(48)
where $`\mathrm{\Omega }=\{q,\lambda ,r,R,t,T,u\}`$, and,

$`\mathrm{\Sigma }_0(\mathrm{\Omega })`$ $`=`$ $`{\displaystyle \frac{\beta ^2}{2}}R^2{\displaystyle \frac{\beta ^2}{2}}T^2rRtT{\displaystyle \frac{1}{2}}\mathrm{log}(2\pi \beta ^2q)uf\beta \lambda q+{\displaystyle \frac{\beta ^2}{4}}u^2q^2+\beta ^2uqR,`$
$`(\mathrm{\Omega },m)`$ $`=`$ $`u\varphi _0(q,m){\displaystyle \frac{1}{2\beta ^2q}}[\varphi _1(q,m)+rm]^2+\mathrm{log}[\varphi _2(q,m)+t]+\lambda m^2.`$ (49)
We can now write the average number of TAP states as,

$$\overline{𝒩(\beta ,f)}=𝒟\mathrm{\Omega }e^{N\mathrm{\Sigma }_0(\mathrm{\Omega })}\underset{i}{}𝑑m_ie^{(\mathrm{\Omega },m_i)}=𝒟\mathrm{\Omega }e^{N\mathrm{\Sigma }_0(\mathrm{\Omega })+N\mathrm{log}{\scriptscriptstyle 𝑑me^{(\mathrm{\Omega },m)}}}.$$
(50)
Thanks to the prefactor $`N`$ in the exponential, the integral in $`𝒟\mathrm{\Omega }`$ can be performed with the steepest descent method. In this way we can write the complexity $`\mathrm{\Sigma }(\beta ,f)`$ as,

$$\mathrm{\Sigma }(\beta ,f)=\frac{1}{N}\mathrm{log}\overline{𝒩(\beta ,f)}=\mathrm{\Sigma }_0(\widehat{\mathrm{\Omega }})+\mathrm{log}𝑑me^{(\widehat{\mathrm{\Omega }},m)},$$
(51)
where $`\widehat{\mathrm{\Omega }}`$ is solution of the saddle point equations,

$$\frac{\mathrm{\Sigma }_0(\mathrm{\Omega })}{\mathrm{\Omega }}+\frac{(\mathrm{\Omega },m)}{\mathrm{\Omega }}=0,$$
(52)
with,

$$𝒪(m)=\frac{1}{𝑑me^{(\mathrm{\Omega },m)}}𝑑m𝒪(m)e^{(\mathrm{\Omega },m)}.$$
(53)
We can reduce the number of variables by directly solving the saddle point equations for $`R`$ and $`T`$,

$`{\displaystyle \frac{\mathrm{\Sigma }_0}{R}}`$ $`=`$ $`0R=r/\beta ^2qu`$
$`{\displaystyle \frac{\mathrm{\Sigma }_0}{T}}`$ $`=`$ $`0T=t/\beta ^2.`$ (54)
In order to have expressions as similar as possible to the ones of BM, we define,

$`B`$ $`=`$ $`\beta ^2(1q)+t`$
$`\mathrm{\Delta }`$ $`=`$ $`\beta ^2(1q)s.`$ (55)
Using the explicit forms of $`\varphi _1,\varphi _2`$ in (39), and relations (54) and (55), we can rewrite (49) as,

$`\mathrm{\Sigma }_0(\mathrm{\Omega })`$ $`=`$ $`\lambda q\beta uf(B+\mathrm{\Delta })(1q)+{\displaystyle \frac{(B^2\mathrm{\Delta }^2)}{2\beta ^2}}`$ (56)
$`{\displaystyle \frac{1}{2}}\mathrm{log}(2\pi \beta ^2q){\displaystyle \frac{\beta ^2}{4}}u^2q^2uq\mathrm{\Delta }\beta ^2uq(1q),`$
$`(\mathrm{\Omega },m)`$ $`=`$ $`\mathrm{log}\left({\displaystyle \frac{1}{1m^2}}+B\right){\displaystyle \frac{[\mathrm{tanh}^1(m)\mathrm{\Delta }m]^2}{2\beta ^2q}}+\lambda m^2+u\varphi _0(q,m),`$ (57)
where $`\varphi _0`$ is given in (37). The remaining saddle point equations (52) for the variable $`\{\lambda ,u,B,\mathrm{\Delta },q\}`$ are,

$`{\displaystyle \frac{\mathrm{\Sigma }}{\lambda }}=0`$ $``$ $`q=m^2`$ (58)
$`{\displaystyle \frac{\mathrm{\Sigma }}{u}}=0`$ $``$ $`\beta f=\varphi _0(q,m){\displaystyle \frac{\beta ^2}{2}}uq^2q\mathrm{\Delta }\beta ^2q(1q)`$ (59)
$`{\displaystyle \frac{\mathrm{\Sigma }}{B}}=0`$ $``$ $`B\left[1\beta ^2{\displaystyle \frac{(1m^2)^2}{1+B(1m^2)}}\right]=0`$ (60)
$`{\displaystyle \frac{\mathrm{\Sigma }}{\mathrm{\Delta }}}=0`$ $``$ $`\mathrm{\Delta }={\displaystyle \frac{\beta ^2}{2}}(1q)+{\displaystyle \frac{1}{2q}}m\mathrm{tanh}^1(m){\displaystyle \frac{\beta ^2}{2}}uq`$ (61)
$`{\displaystyle \frac{\mathrm{\Sigma }}{q}}=0`$ $``$ $`\lambda =B+\mathrm{\Delta }{\displaystyle \frac{1}{2q}}\left\{1{\displaystyle \frac{1}{\beta ^2q}}\left[\mathrm{tanh}^1(m)\mathrm{\Delta }m\right]^2\right\}`$ (62)
$`+u{\displaystyle \frac{\varphi _0(q,m)}{q}}\left[{\displaystyle \frac{\beta ^2}{2}}u^2q+u\mathrm{\Delta }+\beta ^2u(12q)\right].`$
These equations, together with relations (56) and (57) can be compared with equations (15), (16) and (17) of BM bm . The differences are due to the different representation of $`F_{TAP}(m)`$ we have taken, in order to preserve the BRST supersymmetry. However, the final result, i.e. the values of $`\mathrm{\Sigma }_0`$ and $``$ in the solution of the saddle point equations, is exactly the same (in the last of equations (17) of BM, however, there is a term $`uf/q`$ missing). Of course, if we set $`u=0`$, i.e. if we do not impose the constraint on the free energy, our expressions become formally identical to BM.

### IV.3 Role of the supersymmetry and connection with the statics

Solving the saddle point equations above, even numerically, is not a simple task. However, the problem becomes much easier if we make use of the two BRST relations (17) and (18). From (47), (54) and (55), we have,

$`\overline{\psi }_i\psi _i`$ $`=`$ $`T=t/\beta ^2=B/\beta ^2+(1q)`$
$`x_im_i`$ $`=`$ $`R=r/\beta ^2qu=\mathrm{\Delta }/\beta ^2(1q)qu,`$
and thus the first BRST equation becomes,

$$B+\mathrm{\Delta }=\beta ^2qu[\mathrm{BRST1}].$$
(63)
In order to write the second BRST equation we need $`x_ix_i`$, which can be computed from equation (44),

$$x_ix_i=\frac{1}{\beta ^2q}\left\{1\frac{1}{\beta ^2q}[\mathrm{tanh}^1(m)\mathrm{\Delta }m]^2\right\}.$$

In this way the second BRST relation becomes (we recall that $`w^{}=u`$),

$$\frac{1}{q}\left\{1\frac{1}{\beta ^2q}[\mathrm{tanh}^1(m)\mathrm{\Delta }m]^2\right\}=u[B\beta ^2(1q)][\mathrm{BRST2}].$$
(64)
Equation (60) admits the solution $`B=0`$, and following BM this is the solution we adopt. Setting $`B=0`$ into the first BRST relation we obtain,

$$\mathrm{\Delta }=\beta ^2qu,$$
(65)
while substituting equation (37) and the second BRST relation into (62), we find,

$$\lambda =\frac{1}{2}\beta ^2u^2q.$$
(66)
We now use these two results, eqs. (65) and (66), to rewrite the complexity $`\mathrm{\Sigma }(\beta ,f)`$ in (51). By making the change of variable $`mh=\mathrm{tanh}^1(m)`$ in the integral in (51), and using (56) and (57), we obtain,

$`\mathrm{\Sigma }(\beta ,f)=u\left\{\mathrm{log}2{\displaystyle \frac{1}{2u}}\mathrm{log}(2\pi \beta ^2q)+{\displaystyle \frac{\beta ^2}{4}}[(u1)q^2+2q1]+{\displaystyle \frac{1}{u}}\mathrm{log}{\displaystyle 𝑑he^{(h;q,u)}}\right\}\beta uf,`$ (67)
with,

$$(h;q,u)=\frac{h^2}{2\beta ^2q}u\mathrm{log}\mathrm{cosh}h.$$
(68)
This form of the complexity can be fully appreciated by recalling the expression for the free energy of the SK model at the 1-step level of replica symmetry breaking (1RSB). The 1RSB free energy density is given by 1rsb ,

$$\beta F_{1RSB}(\beta ;q_1,x)=\mathrm{log}2+\frac{1}{2x}\mathrm{log}(2\pi \beta ^2q_1)+\frac{\beta ^2}{4}[(x1)q_1^2+2q_11]\frac{1}{x}\mathrm{log}𝑑he^{(h;q_1,x)},$$
(69)
where the self-overlap $`q_1`$ and the breaking point $`x`$ satisfy the saddle point equations $`F_{1RSB}/q_1=0`$, $`F_{1RSB}/x=0`$, and where we have set to zero the mutual overlap $`q_0`$. Clearly, we see that there is a striking formal correspondence between $`\mathrm{\Sigma }`$ and $`F_{1RSB}`$. Indeed we have,

$$\mathrm{\Sigma }(\beta ,f;q,u)=\beta u[F_{1RSB}(\beta ;q,u)f],$$
(70)
and,

$`0={\displaystyle \frac{\mathrm{\Sigma }}{q}}`$ $`=`$ $`\beta u{\displaystyle \frac{F_{1RSB}(\beta ;q,u)}{q}}`$ (71)
$`0={\displaystyle \frac{\mathrm{\Sigma }}{u}}`$ $`=`$ $`\beta u{\displaystyle \frac{F_{1RSB}(\beta ;q,u)}{u}}+\beta [F_{1RSB}(\beta ;q,u)f].`$ (72)
Therefore the saddle point equations for $`\mathrm{\Sigma }`$ and $`F_{1RSB}`$ coincide, provided that $`f=F_{1RSB}(\beta ;q_1,x)`$. Moreover, for $`f=F_{1RSB}`$ we trivially have

$$\mathrm{\Sigma }[\beta ,F_{1RSB}(\beta ;q_1,x)]=0,$$
(73)
and therefore $`f_0=F_{1RSB}`$: as expected, and anticipated in the introduction, the lowest TAP states have the static free energy density.

The fact that $`q_1=q`$ and $`x=u`$, at $`f_0`$ can also be obtained by directly setting $`f=f_0=F_{1RSB}`$ in the saddle point equations (58-62). From (59) and (61), after some algebra, we obtain,

$$\frac{1}{4}\beta ^2q^2u+\frac{1}{2u}\mathrm{log}(2\pi \beta ^2q)\frac{1}{u}\mathrm{log}𝑑he^{(h;q,u)}\mathrm{log}\mathrm{cosh}h=0,$$
(74)
while equation (58) becomes,

$$1+\beta ^2qu[(u+1)q1]\frac{1}{\beta ^2q}h^2=0,$$
(75)
where the averages $``$ are now performed with the distribution $`\mathrm{exp}[(h;q,u)]`$ of (68). Once we set $`u=x`$ and $`q=q_1`$ equations (74) and (75) coincide with the 1RSB saddle point equations, i.e. $`F_{1RSB}/x=0`$ and $`F_{1RSB}/q_1=0`$ 1rsb . Summarizing, the complexity $`\mathrm{\Sigma }`$ vanishes at the 1RSB free energy density, and at this point the two calculations are formally equivalent, since the saddle point parameters coincide, and $`\mathrm{\Sigma }`$ and $`F_{1RSB}`$ are related by equation (70).

Some considerations are in order. First, this identification of the saddle point parameters in the two calculations has a clear physical meaning, already discussed in juanpe . The self-overlap $`q_1`$ of the 1RSB states is just the same as the self-overlap of the TAP states at $`f=F_{1RSB}`$. Less obvious is the relation $`u=x`$, since the 1RSB breaking point $`x`$ does not seem trivially related to the parameter $`u`$. However, we must note that,

$$\frac{d\mathrm{\Sigma }(f)}{df}=\beta u(f),$$
(76)
and thus $`\beta u(F_{1RSB})`$ is the slope of the complexity at the lowest free energy, which is indeed equal to the static breaking point $`\beta x`$.

Secondly, we recall that we have set $`q_0=0`$ in the 1RSB calculation. This is due to the fact that we performed an annealed computation of $`\mathrm{\Sigma }`$, and thus we only had one value of the overlap. We believe that a quenched, but replica symmetric, calculation of $`\mathrm{\Sigma }`$ will be equivalent to the 1RSB static calculation with $`q_00`$. This brings us to the final point: the TAP method has one less step of replica symmetry breaking as compared to the standard static one. This is simply due to the fact that the elementary objects in the TAP approach are states, while in the static one are configurations. For this reason, a model as the $`p`$-spin spherical spin glass, which is exactly solved by a 1RSB Ansatz, has a complexity which is exactly replica symmetric. In the SK model, however, which needs a full RSB static solution, also the correct complexity will need to be computed at a full RSB level juanpe .

### IV.4 Numerical solution of the saddle point equations

We now consider the saddle point equations for $`\mathrm{\Sigma }`$ at a generic value of $`fF_{1RSB}`$. We are left with two unknown variable, $`q`$ and $`u`$, and with three unused equations, that is (58), (59) and (61). However, by using the BRST relations it is possible to show that actually equations (58) and (61) coincide. To prove this fact one must use the formula,

$$h^2=\beta ^2q\beta ^4q^2u+\beta ^4q^2u(u+1)\mathrm{tanh}^2h,$$
(77)
in (64), and then substitute this second BRST relation into either (58), or (61). Thus, the two remaining saddle point equations for $`q`$ and $`u`$, (58) and (59), become,

$`q`$ $`=`$ $`\mathrm{tanh}^2h`$
$`\beta f`$ $`+`$ $`\mathrm{log}\mathrm{cosh}h+{\displaystyle \frac{\beta ^2}{4}}[q^2(1+2u)2q+1]+\mathrm{log}2=0.`$ (78)
where the formula,

$$\frac{1}{2}\mathrm{log}(1m^2)=\mathrm{log}\mathrm{cosh}h,$$

was used.

It is not difficult to solve numerically the two equations above, and this shows how drastic is the simplification of the calculation due to the use of the BRST supersymmetry. In Fig.1 we plot $`\mathrm{\Sigma }`$ as a function of $`fF_{1RSB}(\beta )`$, at $`T=0.7T_c`$ (we recall that $`T_c=1`$). The system of equations (78) have solution only up to a threshold value $`f=f_{th}(\beta )`$, where the complexity takes its maximum value $`\mathrm{\Sigma }_{th}`$. The threshold free energy turns out to be quite close to $`F_{1RSB}(\beta )`$. An expansion of the equations close to $`T_c`$ shows that $`f_{th}=F_{1RSB}+O(ϵ^5)`$, where $`ϵ=T_cT`$, whereas $`q=O(ϵ)`$ and $`u=O(ϵ)`$.

An interesting comparison with the result of BM can be done by computing the total complexity of TAP solutions, $`\mathrm{\Sigma }_{tot}`$. BM obtain this quantity by removing the constraint given by the delta function on the free energy, that is by setting $`u=0`$ in the saddle point equations and in the complexity. The complexity $`\mathrm{\Sigma }_{tot}`$ that BM find in this way is nontrivial (Fig.1 of bm ). If we set $`u=0`$ in our equation (67), however, we obtain a trivial result, that is $`\mathrm{\Sigma }=0`$. This is consistent with our initial remarks on the modulus of the determinant: without the constraint on the free energy we are counting all the stationary points of the TAP free energy, each one weighted with the sign of the determinant. This quantity must be a topological invariant due to the Morse theorem, and thus the complexity must vanish. The fact that this happens in our calculation and not in the calculation of BM shows that the BRST symmetry selects the saddle point solution which preserves Morse’s invariance.

In order to obtain $`\mathrm{\Sigma }_{tot}`$ in our case we simply note that at each temperature the dominant value of the complexity is given by its maximal value, and therefore $`\mathrm{\Sigma }_{tot}(T)=\mathrm{\Sigma }_{th}(T)`$. Thus, the (annealed) average global number of solutions of the TAP equations is given by $`𝒩_{tot}(T)\mathrm{exp}(N\mathrm{\Sigma }_{th}(T))`$. The threshold values $`f_{th}`$ and $`\mathrm{\Sigma }_{th}`$ can be obtained by requiring that the system of equations (78) ceases to have solution, i.e. by imposing a marginality condition on their Hessian. The quantity $`\mathrm{\Sigma }_{th}(T)`$ is plotted in Fig.2: it turns out to be considerably smaller than the non-BRST symmetric value of bm . This might be one of the reasons why it is so difficult to find solutions of the TAP equations numerically. A numerical analysis of $`\mathrm{\Sigma }_{tot}(T)`$ shows that it goes to zero for $`ϵ=T_cT0`$ as,

$$\mathrm{\Sigma }_{tot}(T)=0.01ϵ^6+\mathrm{O}(ϵ^7),$$
(79)
We conclude this Section with an interesting remark on the comparison between the SK and the $`p`$-spin spherical model. In the $`p`$-spin model when we compute the complexity without the constraint on the free energy (i.e. setting $`u=0`$) and without the modulus, we find two saddle points: a BRST symmetric one, which gives $`\mathrm{\Sigma }=0`$, as in the SK case, and a non-BRST symmetric saddle point, which gives a nontrivial complexity, $`\mathrm{\Sigma }^{}0`$. In the $`p`$-spin spherical model this second value $`\mathrm{\Sigma }^{}`$ of the complexity coincides with $`\mathrm{\Sigma }_{th}`$, that is the maximal value of $`\mathrm{\Sigma }(f)`$ obtained by keeping $`u0`$ and which is BRST symmetric. The fact that our result, Fig.2, does not coincide with the result of BM (Fig.1 of bm ), shows that this identity does not hold in the SK model. Thus, the identity $`\mathrm{\Sigma }^{}=\mathrm{\Sigma }_{th}`$ seems to be a peculiarity of the $`p`$-spin spherical model.

### IV.5 The zero temperature limit of the total complexity

The largest value of the total complexity $`\mathrm{\Sigma }_{tot}(T)`$ is achieved at $`T=0`$. We find,

$$\underset{T0}{lim}\mathrm{\Sigma }_{tot}(T)=0.0073.$$
(80)
This value is remarkably smaller than the value $`\mathrm{\Sigma }_{tot}^0=0.199`$ obtained in tana ; dd1 ; bm  by directly solving the zero-temperature limit of the TAP equations. In this Section we want to discuss this difference and its implications.

The TAP equations at finite $`\beta `$ are given by,

$$m_i=\mathrm{tanh}\{\beta [h_i\beta (1q)m_i]\},$$
(81)
where the local field $`h_i`$ is,

$$h_i=\underset{ji}{\overset{N}{}}J_{ij}m_j.$$
(82)
For very large, but finite, $`\beta `$, we can write,

$$m_i=\mathrm{sign}[h_i\beta (1q)m_i],$$
(83)
where the magnetizations are now spin variable, $`m_i=\pm 1`$. In the limit $`\beta \mathrm{}`$ we have $`q1`$, and thus we must be careful with the expression in the square bracket. Let us write in general for $`\beta 1`$,

$$(1q)=h_0\beta ^\alpha ,\alpha >0,$$
(84)
such that the TAP equations become,

$$m_i=\mathrm{sign}\left[h_i\frac{h_0m_i}{\beta ^{\alpha 1}}\right].$$
(85)
According to the value of $`\alpha `$ we have different scenarios. For $`\alpha <1`$ there is no finite limit of the TAP equations (85). This case is therefore uninteresting. For $`\alpha =1`$, on the other hand, we have a well defined zero temperature limit of the TAP equations, i.e

$$m_i=\mathrm{sign}\left[h_ih_0m_i\right],$$
(86)
This is our case: in figure 3 we show that in the limit $`T0`$ the overlap $`q`$ corresponding to the threshold (i.e. total) complexity behaves like,

$$q=1h_0T,h_0=0.49,$$
(87)
so that $`\alpha =1`$. Thus, the zero-temperature limit of the BRST complexity, $`\mathrm{\Sigma }_{tot}^0=0.0073`$, is the complexity of the solutions of equations (86). Note that the $`T`$-dependence of the overlap may be different from (87) at values of the free energy different from the threshold.

In the case $`\alpha >1`$ the zero temperature limit of the TAP equations (85) yields,

$$m_i=\mathrm{sign}\left[h_i\right].$$
(88)
These are the equations considered in tana ; dd1 ; bm  in order to compute the zero temperature complexity (in particular, equations (88) were obtained by assuming that $`q=1O(T^2)`$, i.e. $`\alpha =2`$). Therefore, the value $`\mathrm{\Sigma }_{tot}^0=0.0073`$ we obtain as zero-$`T`$ limit of the BRST complexity and the value $`\mathrm{\Sigma }_{tot}^0=0.199`$ obtained in tana ; dd1 ; bm , refer to two different sets of equations, that is (86) and (88). The fact that solutions of (86) have a much smaller complexity than those of (88), is due to the fact that the former are actually a subset of the second. Indeed, from equations (86) we have,

$$|h_i|>h_0,$$
(89)
for all $`i`$, which is a very restrictive condition not necessarily satisfied by the solutions of (88).

We note that equations (85) imply,

$$m_ih_i>\frac{h_0}{\beta ^{\alpha 1}},$$
(90)
and thus, given a TAP solution, the change in energy $`\mathrm{\Delta }E(K)`$ when we flip $`K`$ spins, satisfies the relation,

$$\mathrm{\Delta }E(K)=\underset{i}{\overset{K}{}}m_ih_i\underset{ij}{\overset{K}{}}J_{ij}m_im_j\underset{i}{\overset{K}{}}m_ih_i\frac{K}{\sqrt{N}}>K\left(\frac{h_0}{\beta ^{\alpha 1}}\frac{1}{\sqrt{N}}\right).$$
(91)
The condition of stability under $`K`$ spin flips spinflip  requires $`\mathrm{\Delta }E(K)>0`$. A sufficient (but not necessary) condition for this to hold, is given by,

$$\frac{h_0}{\beta ^{\alpha 1}}\frac{1}{\sqrt{N}}>0.$$
(92)
In the case $`\alpha =1`$ this condition is trivially satisfied, while for $`\alpha >1`$ the order of the two limits $`\beta \mathrm{}`$ and $`N\mathrm{}`$ becomes important, and stability under $`K`$ spin flips may be trickier to prove. Nevertheless, a word of caution on the concept of stability is in order at this point. At finite temperature a stable solution of the TAP equations is a local minimum of the TAP free energy, that is a solution with positive Hessian. Unstable saddle points of the TAP free energy will exist, typically at high free energies. Unfortunately, at $`T=0`$ this topological stability is impossible to assess, because of the discrete nature of the model. On the other hand, the $`K`$-spin flips stability is easily defined at $`T=0`$. However, it is not straightforward to understand what is the relation between these two definitions of stability. We simply note that we have computed the total complexity by finding the threshold states: for free energies above the threshold value the saddle points equations become unstable. Results in continuous mean-field spin-glasses (as the $`p`$-spin spherical model) suggest that the threshold may indeed be the border between topologically stable solutions, that is TAP minima, and unstable saddle points. It is therefore plausible that our total complexity refers only to TAP minima. The connections we have discovered between complexity and thermodynamics support this hypothesis.

## V Connection with the Legendre transform approach

As discussed in the previous Section there is a precise correspondence between the TAP computation at the lower band edge $`f_0`$ and the computation of the thermodynamic free energy of the system. However, the relation between the TAP approach and the usual static one (which involves Boltzmann averages) is still deeper: as it can be appreciated in (70) there is a formal connection between the TAP complexity at generic $`f`$ and the free energy function $`F_{1RSB}`$. These equations also represent an important bridge between the TAP approach and some different methods to compute the complexity which do not rely on TAP equations, but rather on constrained thermodynamics monasson ; franzparisi . In monasson  it is argued that the complexity $`\mathrm{\Sigma }(\beta ,f)`$ of the ergodic components - i.e. states - present at low temperature is given by the Legendre Transform of the free energy $`F(\beta ,n)`$ of $`n`$ coupled real replicas of the original system. More precisely,

$$\mathrm{\Sigma }(\beta ,f)=\underset{n}{\mathrm{max}}[\beta nf\beta F(\beta ,n)],$$
(93)
where $`n`$ and $`f`$ are Legendre conjugated variables, i.e.

$$f=\frac{F(\beta ,n)}{n}.$$
(94)
If we assume, as it is generally accepted, that TAP minima correspond to low temperature states, then the complexity computed via eq.(93) and the complexity obtained from the TAP approach should be equal. In other words, we expect that the TAP complexity to satisfy eqs. (93) and (94).

To establish this point one has first to compute the free energy $`F(\beta ,n)`$ of the coupled system with the replica method and choose an appropriate Ansatz. This problem is discussed for the SK model in potters , where different kinds of Ansatz are analyzed: within a generalized two-group Ansatz for $`F(\beta ,n)`$, equation (93) is formally satisfied by the BM solution, as already noted by these authors in bm ; bm1 . Moreover, a simpler sub-class of two-group saddle point solutions exists, which corresponds to the following procedure: $`F(\beta ,n)`$ at the $`k`$-level of RSB is given by $`n`$ times the free energy of one system evaluated with $`k+1`$ steps of RSB, and breaking parameters $`nx_1,nx_2,\mathrm{}nx_k,n`$. Thus, at the lowest possible level of RSB ($`k=0`$) one has,

$$F(\beta ,n)=nF_{1RSB}(\beta ;q_1,q_0,n),$$
(95)
as also suggested in monasson . With this expression of $`F(\beta ,n)`$, eq. (93) has been verified within the TAP approach for models which are exactly solved by 1RSB replica solution potters2 ; juanpe . It has also been derived for the Bethe Lattice at zero temperature with the cavity method cavity , and seems therefore quite robust. Surprisingly, for the SK model this point was still unclear, since the BM solution bm , despite its formal correspondence to a more general two-group static solution, does not appear to satisfy relation (93) with the simpler Ansatz (95). The reason for this can be understood in the light of the previous Section: in bm  BM do not use the BRST supersymmetry and consequently consider a larger set of solutions than those physically relevant, much as the generalized two-group Ansatz has a larger set of saddle points than the simpler Ansatz (95). We therefore may expect that once the BRST relations are taken into account, as we have done in the previous Section, consistency must be recovered. Indeed, this is precisely what happens: as it can be easily verified, once (95) is implemented, eq.(70) is equal to (93), with $`q_0=0`$, $`q_1=q`$ and $`n=u`$. The parameters $`q_1`$ and $`n`$ are variationally fixed: (72) is equal to (94), while variation with respect to $`q_1`$ gives back equation (71).

Thus, also for the SK model, we have demonstrated within the TAP approach the validity of relation (93). We note that a crucial element to establish this point is the exploitation of the BRST relations. This is not a surprise since, as discussed in Section III for equilibrium averages, the BRST supersymmetry is what mathematically guarantees relations of physical significance as (23) and (24).

## VI Conclusions

In this paper we have discussed the equivalence between the TAP approach and the standard thermodynamic method for mean-field spin-glass systems, and we have shown that the use of the BRST symmetry of the TAP action is crucial to establish this equivalence.

We demonstrated in Section III that the BRST relations guarantee the mathematical consistency of the TAP approach and ensure the validity of physically relevant relations. We have shown that the DDY computation, which establishes for the SK model a formal connection between the TAP and the thermodynamic partition functions, is generally valid if the BRST symmetry is taken into consideration, and does not rely on any hypothesis about the nature of the equilibrium states. Furthermore, we have revisited the BM calculation of the complexity for the SK model, with the support of the BRST relations. Thanks to them, we have been able to solve the saddle point equations at given free energy, we have explicitly shown the equivalence of the lowest TAP solutions with the thermodynamic equilibrium states and we have verified the consistency of the TAP computation of the complexity with the Legendre Transform construction.

It should be noted that the BRST relations have a very delicate role due to the presence of quenched disorder. Indeed, once the average over the disorder is taken and the auxiliary variables are introduced, the BRST symmetry of the effective action is much less transparent. In the thermodynamic limit, the auxiliary variables become variational parameters which must satisfy the saddle point equations. The solutions of these equations correspond to averaged quantities (for example, in the BM calculation $`R=x_im_i`$) and should therefore satisfy the appropriate BRST relations. However, these relations are not automatically satisfied by all the saddle point solutions. In other words, the set of solutions of the saddle point equations in general includes also solutions which are not BRST invariant. Thus, the BRST relations must be explicitly imposed, providing some additional equations that must be satisfied by the variational parameters. Besides, as we have seen, this procedure provides a great simplification in the computation, since it drastically decreases the number of unknown variables.

Our BRST calculation gives a total complexity at $`T=0`$ which is much smaller than the value previously obtained by other classic calculations. We have shown that technically this is due to the fact that the BRST saddle point provides a set of TAP equations at $`T=0`$ which are different from the ones previously considered. In particular, the zero-temperature BRST-TAP solutions have local fields larger than a finite value $`h_0`$ at all sites, and this condition is responsible for the drastic decrease in the complexity. We have outlined what are the possible physical explanations of this difference, but more study on the $`T0`$ limit is necessary for a complete clarification of this point.

###### Acknowledgements.

We thank T. S. Grigera, A. Pagnani and F. Ricci-Tersenghi for numerous interesting discussions, and A. Bray and M. A. Moore for many important remarks on the manuscript.