text:
# Evidence for prolonged main sequence stellar evolution of F stars in close binariesBased on the data from the Hipparcos astrometry satellite (European Space Agency)

## 1 Introduction

Many of the F stars cataloged in the HIPPARCOS as single were recently argued to be in fact unidentified binaries (Suchkov & McMaster 1999). The criterion used to isolate binary candidates involves the difference between the absolute magnitude $`M_V`$ based on the HIPPARCOS parallax and the absolute magnitude $`M_{c_0}`$ derived from the uvby luminosity index $`\mathrm{\Delta }c_0`$, i.e., $`\mathrm{\Delta }M_{c_0}=M_{c_0}M_V`$. For the best known nearby single F stars ($`d<25`$ pc), the two magnitudes coincide within $`\pm 0.15`$ mag. However, most of the unresolved binary stars are anomalously bright for their uvby colors, i.e., their absolute magnitude, $`M_V`$, is substantially brighter than the magnitude derived from $`\mathrm{\Delta }_{c_0}`$: $`\mathrm{\Delta }M_{c_0}>0.15`$ (for simplicity, we will often be referring to $`\mathrm{\Delta }M_{c_0}`$ as brightness anomaly and consider a star as anomalously bright if $`\mathrm{\Delta }M_{c_0}>0.15`$). Unlike the nearby single stars, many of the distant single F stars, ($`d>25`$ pc) turned out to be anomalously bright as well, similar to the known binary F stars. The analysis showed that in this case, the brightness anomaly is likely associated with the presence of undetected companions (Suchkov & McMaster 1999). These binary candidates were called C binaries (the follow-up radial-velocity survey of northern C binaries with $`\mathrm{\Delta }M_{c_0}0.5`$ and $`d<80`$ pc conducted by Griffin & Suchkov 2001 has already confirmed that more than $`35\%`$ of the sample stars actually are binaries; orbits for many of these binaries have been determined).

For unresolved binaries, $`M_{c_0}`$ is derived from the color indices originating from a double star combined spectrum. Because of that $`\mathrm{\Delta }M_{c_0}`$ underestimates the luminosity of that star if its components are comparably bright, by up to $`2.5\mathrm{log}20.75`$ mag in the case of identical components (this obviously imposes an upper limit of 0.75 on $`\mathrm{\Delta }M_{c_0}`$). However, this seems to be not the only, and perhaps not even the main effect leading to the discrepancy between $`M_{c_0}`$ and $`M_V`$. For one thing, a substantial fraction of binary stars have brightness anomaly well above the indicated upper limit of 0.75 mag. The stars with $`\mathrm{\Delta }M_{c_0}>0.75`$, which cannot be explained in terms of the combined flux of the unresolved binary components, indicate that along with the combined luminosity something else may impact the absolute magnitude of unresolved binary stars.

More evidence favoring the existence of an additional source for brightness anomaly has come from the study of the relationship between $`\mathrm{\Delta }M_{c_0}`$ and stellar age (Suchkov 2000). It turns out that in terms of kinematics and metallicity, anomalously bright binary candidates, C binaries, are, on average, substantially older than the “truly single” F stars, i.e., the stars with normal brightness. If the anomaly were caused entirely by the luminosity contribution from the undetected secondary, one would have to conclude that binaries with the components having about the same mass live longer or tend to form earlier in the galactic evolution than binaries with disparate components. But such a scenario does not seem very likely. A more plausible interpretation of the relationship between $`\mathrm{\Delta }M_{c_0}`$ and age might be that at least some part of the star’s brightness anomaly is associated with a peculiar stellar evolution in the binary system, so that the primary F star in that binary stay within the main sequence longer than a normal single F star. Given the obviously important implications of such a possibility, we have examined a number of relationships between $`\mathrm{\Delta }M_{c_0}`$ and other parameters of F stars in order to get a deeper insight into the nature of brightness anomaly.

## 2 Data

The star sample used in this study originates from the list of $`10,000`$ HIPPARCOS stars that have measured uvby colors, with effective temperatures corresponding to the spectral range occupied by F and early G stars (Suchkov & McMaster 1999). The uvby data are from Hauck & Mermilliod (1998). We have computed tangential velocities of these stars from the HIPPARCOS parallaxes and proper motions. The HIPPARCOS parallax has also been utilized to derive absolute magnitude $`M_V`$ from the Johnson $`V`$ magnitude and absolute magnitude $`M_{Hp}`$ from the HIPPARCOS $`Hp`$ magnitude. Absolute magnitude $`M_{c_0}`$ has been calculated from the uvby data with the algorithm published by Moon (1985; see also Moon & Dworetsky 1985); these calculations involve the dereddened uvby luminosity index, $`\mathrm{\Delta }c_0`$, along with the dereddened metallicity and temperature indices. The algorithm in Moon (1985) has also been used to compute standard $`\beta `$ from the blanketing-corrected color index $`(by)_{cor}`$, from which effective temperature has been obtained.

Metallicity has been calculated from the dereddened uvby color indices, utilizing the calibration given by Schuster & Nissen (1989) for F stars. This calibration is somewhat different from that in Carlberg et al. (1985), which was used in Suchkov (2000) and Suchkov & McMaster (1999). This difference, however, has been found not to impact the results of both the previous and present studies. The Schuster & Nissen (1989) calibration is commonly utilized in current studies, so we will use it in this paper.

Age estimates are based on isochrone fitting in the $`\mathrm{log}T_eM_V`$ diagram using Yale 1996 isochrones (Demarque et al. 1996111Available at http://shemesh.gsfc.nasa.gov/astronomy.html). Age has been used only in terms of mean values obtained by averaging over large groups of stars, and only in differential analysis involving relative ages. Therefore, age accuracy, both in terms of random and systematic errors, is not of much concern for the present study.

The sample has been constrained as follows: $`0.22(by)0.39`$; $`5800T_e7500`$, K; $`0.6[\mathrm{Fe}/\mathrm{H}]\text{ }0.5`$; $`0.02e_{c1}`$ ($`e_{c1}`$ is the mean error on $`c_1`$ for an individual star as given in Hauck & Mermilliod 1998; the sample mean for $`e_{c1}`$ is $`10`$ times smaller); $`0.015e_{m1}`$ ($`e_{m1}`$ is the mean error on $`m_1`$ for an individual star); $`3e_\pi `$, mas (error on parallax); $`3e_{\mu _\alpha }`$, mas $`\mathrm{yr}^1`$(error on proper motion in RA); $`3e_{\mu _\delta }`$, mas $`\mathrm{yr}^1`$(error on proper motion in DEC).

## 3 Results and discussion

### 3.1 Brightness anomaly distribution

Brightness anomaly is very conspicuous in unresolved binary stars (Fig. 1). As shown in Suchkov & McMaster (1999), the distribution of $`\mathrm{\Delta }M_{c_0}`$ for the best studied single stars within 25 pc is centered at $`\mathrm{\Delta }M_{c_0}=0`$ and is rather narrow, with standard deviation of only 0.15 mag. The corresponding gaussian ($`\sigma _{\mathrm{\Delta }M_{c_0}}=0.15`$) is shown in Fig. 1. For the unresolved binaries, the maximum of the $`\mathrm{\Delta }M_{c_0}`$ distribution is not only offset from zero by $`0.5`$ mag but is also much broader. Within 200 pc, the stars cataloged in the HIPPARCOS as single have a broad distribution as well, but with the maximum only slightly shifted ($`0.1`$ mag) from that of the single stars within 25 pc.

Note that the left tail of the distribution, albeit resembling a segment of a gaussian, is far above the level suggested by $`\sigma _{\mathrm{\Delta }M_{c_0}}=0.15`$, meaning that it does probably not represent the error distribution for $`\mathrm{\Delta }M_{c_0}`$. Therefore, the population of anomalously faint stars, $`\mathrm{\Delta }M_{c_0}<0.15`$, should be as real as that of anomalously bright stars. As seen below, its distinctive identity is reflected in lowest velocities and highest metallicities among the F stars, which makes it the youngest group of F stars. Therefore, we conclude that in general the distributions of both binary and single stars in Fig. 1 are most likely dominated by some stellar population physics rather than random or systematic errors.

The most striking feature in Fig. 1 is the positive brightness anomaly of binary stars. Therefore, in the rest of the paper, we will focus primarily on the nature of anomalously bright binary stars.

### 3.2 Correlation between $`\mathrm{\Delta }M_{c_0}`$ and binary components separation

If the combined flux of unresolved binary components is not the only reason for the enhanced brightness of a binary star, and there is a noticeable contribution from a physical anomaly in the binary’s primary star caused by its interaction with the companion, one may expect the brightness anomaly to depend on binary components separation, because the effect of the interaction, whatever it is, is probably stronger in tighter pairs. We have tested this hypothesis by correlating brightness anomaly and projected separation for two very different, non-overlapping samples of binary F stars within 125 pc, for which the HIPPARCOS provides angular separation, $`\rho `$. Only the stars with $`\mathrm{\Delta }M_{c_0}`$ in the range $`0.5`$ to 1.0 have been included, so that the potentially less reliable extreme values of $`\mathrm{\Delta }M_{c_0}`$ do not impact the results. The chosen distance limit rejects the most distant stars, thus ensuring better data quality for the remainder of the sample; at the same time it retains enough stars for statistical tests we are interested in.

The stars in the first sample are binaries discovered by HIPPARCOS (see Mignard et al. 1992, Söderhjelm et al. 1992 for the procedures to detect and measure binary stars on the basis of HIPPARCOS observations). Only the binaries with angular separation $`\rho 1`$ arcsec have been included; this constraint rejects only a very small fraction of the HIPPARCOS binaries and is not of great significance for the test below. The sample has also been constrained with respect to angular separation error, $`ϵ_\rho `$, to include only the stars with $`ϵ_\rho /\rho <0.1`$.

The second sample comprises unresolved binaries known prior to HIPPARCOS. With respect to angular separations, it has been constrained similar to the sample of the HIPPARCOS binaries, $`\rho 1`$ arcsec. The pre-HIPPARCOS binaries are much more numerous, which allows us to impose a more stringent constraint on the angular separation accuracy: $`ϵ_\rho /\rho <0.01`$. Additionally, this sample has been constrained to include only hot stars, $`65007500`$ K. This should improve chances to detect any separation-dependent evolutionary effects in $`\mathrm{\Delta }M_{c_0}`$ just because the fraction of evolved stars increases toward higher temperatures (at cooler temperatures, brightness anomaly in the sample stars may be dominated by the separation-independent effect of the combined flux of binary’s unresolved components).

Projected separation for both samples has been derived from the parallax, $`\pi `$, and the angular separation, $`\rho `$, given in the HIPPARCOS catalogue.

The stars of the two samples are presented in Fig. 2 in Fig. 3. As one can see, they do indeed show a correlation between separation and brightness anomaly, with the slope of the linear regression being non-zero well above the error in both cases (the slope error is given as a percent of the slope value).

It is to be noted that the correlation remains the same when absolute magnitude $`M_V`$ in the definition of $`\mathrm{\Delta }M_{c_0}`$ is replaced with absolute magnitude $`M_{Hp}`$ obtained from the HIPPARCOS $`Hp`$ magnitude (for the details of HIPPARCOS photometry see Mignard, Froschle, & Falin 1992; Evans et al. 1992). The corresponding regression slopes are: $`0.005\pm 25.2`$% for the HIPPARCOS binaries, and $`0.011\pm 21.6`$% for the pre-HIPPARCOS binaries.

The correlation in Fig. 2 and Fig. 3 reveals that $`\mathrm{\Delta }M_{c_0}`$ is, on average, larger at smaller component separations. This suggests that brightness anomaly is somehow associated with the components interaction in a tight binary pair. The question is then: what is the physics behind this anomaly?

### 3.3 Relationship between $`\mathrm{\Delta }M_{c_0}`$ and stellar kinematics

To answer the above question, we have explored whether there is any dependence in $`\mathrm{\Delta }M_{c_0}`$ on age, because the effects caused by interaction in a close binary may accumulate with time. First of all, we have correlated $`\mathrm{\Delta }M_{c_0}`$ and tangential velocity, employing the fact that the statistics from stellar kinematics, such as mean tangential velocity and velocity dispersion used here, are age-sensitive.

As seen in Fig. 4, there is indeed a significant correlation between $`<v_{\mathrm{tangent}}>`$ and $`\mathrm{\Delta }M_{c_0}`$. It implies that $`\mathrm{\Delta }M_{c_0}`$ is larger, on average, at older ages. Given that the HIPPARCOS single stars with large $`\mathrm{\Delta }M_{Vc_0}`$ are probable close binaries, we conclude that this result is consistent with increase in brightness anomaly for a tight binary pair as the primary star evolves.

Fig. 4 also reveals that the kinematics of known binary stars is “younger” than that of the single stars at the same $`\mathrm{\Delta }M_{c_0}`$. This reflects the generally younger age of the population of the known binary F stars, which was discussed in Suchkov (2000).

### 3.4 Relationship between $`\mathrm{\Delta }M_{c_0}`$ and metallicity

Mean metallicity is another age-dependent statistics that can be used to check if there is any relationship between $`\mathrm{\Delta }M_{c_0}`$ and age. The correlation between $`[\mathrm{Fe}/\mathrm{H}]`$  and $`\mathrm{\Delta }M_{c_0}`$ is shown in Fig. 5 for the same groups of stars as in Fig. 4. Here, the linear regression slope is non-zero within more than one sigma for both the binary and single star samples, indicating that $`\mathrm{\Delta }M_{c_0}`$ is larger at lower metallicities. Since lower metallicity is indicative of older age, the correlation in Fig. 5 is consistent with the inference from the preceding subsection that the stars with larger $`\mathrm{\Delta }M_{c_0}`$ are, on average, older.

The correlation in Fig. 5 could have been expected on the basis of the correlation in Fig. 4 because metallicity is known to correlate with kinematics. At the same time this correlation is obviously weaker than that in Fig. 4. Since metallicity is also known to have a large spread at any given age, a possible interpretation of weaker correlation may be that metallicity is coupled with age more loosely than kinematics.

The weaker correlation in Fig. 5, whatever its reason, argues against the possibility that the correlation in Fig. 4 results from a metallicity-dependent bias in estimates of $`M_{c_0}`$, hence $`\mathrm{\Delta }M_{c_0}`$. If such a bias did exist (say, due to flaws in the algorithm that computes $`M_{c_0}`$) and cause all the correlation in Fig. 5, one would have expected it to translate into a weaker rather than stronger correlation between $`\mathrm{\Delta }M_{c_0}`$ and kinematics in Fig. 4. It is to be added that the aforementioned hypothetical bias seems unlikely anyway, because the available evidence, including the results from the binary components separation data and the age–velocity relation (see below), argues against any role of this kind of systematic errors in brightness anomaly.

Thus, the interpretation of both Fig. 4 and Fig. 5 involving the dependence of brightness anomaly on age appears to be reasonably justified.

### 3.5 Age–velocity relation: discrepancy between normal and anomalously bright stars

So far, we have been referring to brightness anomaly in the sense that the absolute magnitude of a star with this anomaly is brighter than that of a normal single star with the same uvby colors. But it would be equally correct to say that the star has the uvby colors that are anomalous for a normal star at a given effective temperature and luminosity. The anomaly in this case refers to the color index $`c_0`$.

For the main sequence F stars, $`c_0`$ provides a measure of Balmer discontinuity (see Crawford 1975 and references therein). In these stars, smaller Balmer discontinuity at a given effective temperature (hence smaller $`c_0`$) means a denser stellar atmosphere, i.e., larger surface gravity. In this way, $`c_0`$ actually measures stellar surface gravity (at a given $`T_e`$), through which it is related to the star’s luminosity. Anomalous $`c_0`$ in an F star then means anomalous surface gravity, $`\mathrm{log}g`$. Specifically, the anomaly resulting in a positive $`\mathrm{\Delta }M_{c_0}`$ corresponds to the surface gravity that is too high for a given effective temperature, $`T_e`$, and luminosity, $`L`$. Consequently, the stellar mass of such a star is larger than the mass of a normal single star at the same $`T_e`$ and $`L`$. (Alternatively, one can say that the luminosity of an anomalously bright star is too high for a given $`T_e`$ and $`\mathrm{log}g`$). This obviously requires that stellar evolution of a normal star is different from that of an anomalously bright star. So, if $`c_0`$ of the anomalous stars in our sample is dominated by the effect of enhanced surface gravity, these stars can be predicted to evolve in a non-standard way.

We can verify the latter prediction as follows. Let us take the sample of the HIPPARCOS “single” stars and split it into two subsamples: the first one including only anomalously bright stars, $`\mathrm{\Delta }M_{c_0}>0.15`$ (C binaries), and the second subsample that includes the stars with $`\mathrm{\Delta }M_{c_0}<0.15`$; for simplicity, the latter stars will be referred to as “truly single”, although we expect a fraction of them to be unevolved and/or wide binaries with disparate components. Now, we can derive and compare the age–velocity relations (AVR) for these two group of stars. If the same stellar evolution model (meaning the same set of isochrones used to derive age) is applicable to all stars, we should obviously get the same AVR for both groups, simply because stellar kinematics and stellar evolution are entirely unrelated things. If this proves not to be the case, and the age–velocity relations are discrepant, the two groups of stars evolve probably differently.

We have derived ages for both groups of stars based on isochrones from Demarque et al. (1996). The isochrones as well as the stars used to calculate the age–velocity relations are shown in Fig. 6. The stars have been selected in a narrow metallicity range, $`0.3<[\mathrm{Fe}/\mathrm{H}]\text{ }<0.1`$, that matches the chemical composition of the isochrones, $`Z=1.0\times 10^2,Y=0.27`$. Additionally they have been constrained to the distance range $`d125`$ pc. They have been grouped into 1 to 2 Gyr age bins, for which mean tangential velocity has been calculated. Because of large age errors for the young stars near the ZAMS, ages younger 1 Gyr have not been considered. The resulting age–velocity relations for the age range from 1 to 8 Gyr are given in Fig. 7.

Fig. 7 reveals a conspicuous discrepancy between the two age-velocity relations, which will be referred to as the AVR discrepancy. The possibility that kinematic differences between the normal and anomalously bright stars arise from different distance sampling (the latter stars are, on average, more distant) was checked and ruled out in Suchkov (2000). We have also looked into the possibility that the bias associated with the brighter, on average, absolute magnitudes of anomalously bright stars (see Fig. 6) plays a role in the discrepancy (for example, this may happen if the isochrones underestimate the age of far evolved stars). To this end, we have compared the kinematics of the normal and anomalously bright stars with $`0.6<[\mathrm{Fe}/\mathrm{H}]\text{ }<0.3`$ in the two regions of the $`\mathrm{log}T_eM_V`$ diagram, above and below the line corresponding to 1.0 mag above the ZAMS shown in Fig. 6. As seen in Table 1, anomalously bright stars close to the ZAMS have velocities higher, on average, than those of the normal stars far above the ZAMS. This rules out the indicated bias as the cause of the AVR discrepancy. Therefore, the discrepancy must be real, reflecting some inconsistency in the derived ages. Assuming that the isochrones in Fig. 6 adequately represent the age of the normal stars (truly single), we must then conclude that the age of the stars with large brightness anomaly (C binaries) is underestimated. The amount of the AVR discrepancy suggests that the underestimation is, on average, as large as 2 Gyr or more. Thus, anomalously bright stars are consistent with being actually older than normal stars at the same positions in the $`\mathrm{log}T_eM_V`$ diagram.

There are at least two reasons why isochrone fitting could be inadequate for anomalously bright stars and result in incorrect ages: (i) combined flux of unidentified binaries with comparably bright components, and (ii) non-standard stellar evolution.

If anomalously bright stars are unidentified binaries with comparably bright components, isochrone fitting results in either overestimated or underestimated age, depending on the position of the star in the $`\mathrm{log}T_eM_V`$ diagram and the contribution of the secondary into the star’s $`V`$ magnitude. For the stars cataloged in the HIPPARCOS as single, age would be, on average, underestimated (Suchkov & McMaster 1999).

The following, however, argues against the hypothesis that such binaries dominate our sample of anomalously bright stars.

First, the sample stars are much older (on average) than the known unresolved binaries (Suchkov, 2000; see also Fig. 4 and Fig. 5). This means that most of these stars are different from normal binaries. Even if they are double systems, in most cases their secondary component is hardly a regular star whose emission is comparable to that of the primary.

Second, as mentioned before, even among known binary stars a substantial fraction have brightness anomaly exceeding the maximum value of 0.75 mag, hence it cannot be attributed to the combined flux of the binary ’s two components only (see Fig. 1). Statistically, the excess is quite significant, so there is little doubt that in many cases the cause of brightness anomaly is different from the light contribution from the star’s unresolved companion.

So, it seems unlikely that the AVR discrepancy in Fig. 7 is due only to binaries with comparably bright components.

The above isochrone fitting can also be inadequate if for some reasons the evolution of a star is different from the normal evolution of a single star.

In such a case the position of a star in the $`\mathrm{log}T_eM_V`$ diagram would correspond to an age different from the one predicted by standard isochrones. The age–velocity relation based on ages from the standard isochrones would be incorrect for these stars, and this may be the reason for the AVR discrepancy revealed by Fig. 7.

As discussed below, anomalous stellar evolution, unlike the effect of combined flux, offers a way to explain the relationship between age and brightness anomaly; in particular, it explains why C binaries are, on average, older than the truly single stars. Also anomalous evolution appears to be in line with the recent developments in studies of eclipsing binaries. This makes it a stronger candidate for being the dominant source of the AVR discrepancy in Fig. 7. The combined evidence discussed in this paper suggests then that anomalous stellar evolution occurs apparently in tight binary pairs where the evolution of the primary star is affected by the presence of the companion.

## 4 Discussion and summary

The older ages of anomalously bright stars (Suchkov 2000) seem to indicate that these stars stay within the main sequence longer than normal single stars. Such a prolonged main sequence evolution is known to occur in stars with enhanced central mixing in the convective core, and is apparently quite ubiquitous among close binary systems at stellar masses slightly above solar (see, e.g., Chiosi 1999, Andersen & Nordström 1999, Gimenéz et al. 1999). This suggests that anomalous brightness discussed in this paper reflects non-standard stellar evolution that involves enhanced core mixing. Extreme manifestation of brightness anomaly among known unresolved binary stars seems to indicate that the anomaly is somehow related to binarity. Given that $`\mathrm{\Delta }M_{c_0}`$ correlates with binary components separation, one may infer that interaction in a tight binary pair somehow induces extra mixing in the convective core of the binary’s primary component, additional to whatever mixing occurs in a similar single star. Of course, there are other mechanisms that may impact stellar evolution in a tight binary pair. One may recall, for example, modified convection in the outer convective envelope of a star in a tight binary pair or mass loss/transfer in very close binary systems. These mechanisms were discussed in the literature as related to the existence of blue stragglers and other unusual stars in stellar clusters (Pols & Marinus 1994), apparent age discrepancy between the components in low mass binaries (Clausen et al. 1999), young age of the known unresolved binary stars (Suchkov 2000), etc. However, enhanced convective core mixing seems to be the only mechanism capable of extending the main sequence evolution by a significant amount.

It is to be noted that despite the potentially longer than normal lifetime, the actual average main sequence lifetime of the known unresolved binary stars is shorter than that of the single stars (Suchkov 2000). As hypothesized in the latter paper, this may be related to violent interaction of the binary’s stellar components, which effectively removes the binary’s primary in the color–magnitude diagram from the region occupied by the main sequence F stars; mechanisms of such an interaction were considered, for instance, in Pols & Marinus (1994). On the other hand, C binaries were found to be indeed significantly older, on average, than the truly single stars, possibly indicating that many of C binaries have substellar secondaries incapable of causing too much damage to the primary F star.

Unidentified close binaries may represent a large fraction of the local population of F stars. Assuming that “single” F stars with $`\mathrm{\Delta }M_{c_0}>0.15`$ are in fact tight binary pairs, the fraction of such binaries is at least $`40`$ %. Thus, close binaries may heavily impact stellar population statistics used to probe the formation and evolution of the Galaxy. Therefore, these stars are important not only for better understanding of stellar evolution, but also for developing the adequate picture of the history of the Galaxy.

In conclusion, we can summarize our results as follows. (i). The unresolved binary F stars from the HIPPARCOS are typically much brighter than predicted on the basis of their uvby colors. (ii). This discrepancy, called here brightness anomaly, inversely correlates with the binary components separations, suggesting that at least part of it is associated with binary components interaction in tight pairs. (iii). Brightness anomaly correlates with stellar kinematics in a way that implies higher, on average, anomaly values at older ages. (iv). With ages from the same set of isochrones, the age–velocity relations of the normal single stars and binary candidates with brightness anomaly (C binaries) are discrepant in a way implying that the actual ages of C binaries are older than predicted.

These results argue that the main sequence evolution of the primary F star in a tight binary pair is different from that of a similar single star, evidently because of the interaction with the secondary. The difference appears to be consistent with extra mixing in the convective core of the primary, presumably induced by that interaction; however, other mechanisms need to be carefully studied before a definitive conclusion can be reached.

###### Acknowledgements.

I would like to thank the anonymous referee for helpful comments and suggestions. It is my pleasure to thank Claus Leitherer, Stefano Casertano, and Nikolai Piskunov for stimulating discussions. I am grateful to S. Lyapustina and G. Galas for careful reading the paper and numerous revision suggestions.


format: text

text:
# Surface effects in multiband superconductors. Application to MgB2.

## I Introduction

Recent experiments report the existence of superconductivity at nearly 40K in MgB2. Its origin is not completely elucidated. The material shows a pronounced isotope effect, and the density of states is well approximated by the BCS theory. Tunneling experiments suggest that the superconducting properties at the surface of the material differ from the expected bulk behavior. In addition, photoemission results suggest the existence of an s-like gap, $`\mathrm{\Delta }`$, such that $`\mathrm{\Delta }3k_BT_c`$. A possible explanation of this result is that the measured $`\mathrm{\Delta }`$ is the average of different gaps. The present work is motivated by the persistent discrepancy between the gap values measured in different experiments, and, particularly, the excellent fit to a BCS gap too low to explain the value of the critical temperature observed in tunneling experiments reported in.

Band structure calculations suggest that there are, at least, two types of bands at the Fermi surface: a hole band, built up of boron $`\sigma `$ orbitals, with a weak dispersion in the direction perpendicular to the boron planes, and a broader band, built up mainly of $`\pi `$ boron orbitals, which shows a significant dispersion in the direction perpendicular to the boron planes. Theoretical arguments favor, as the origin of the superconductivity, the hole like $`\sigma `$ band, or the $`\pi `$ band. The existence of two bands with different physical properties is assumed in other models for the superconducting properties of MgB2. It has been argued that the upper critical field can be best modelled if the superconducting properties depend on the specific band at the Fermi level. On general grounds, it is reasonable to assume that the $`\sigma `$ and $`\pi `$ bands in MgB2 will have different contributions to the superconducting phases, and that the superconducting gap needs not be the same in the two bands.

The existence of many bands at the Fermi level, with very different physical properties, is probably a generic feature of intermetallic superconductors. In these materials, it can be expected that the pairing interaction which gives rise to the superconductivity will depend on the details of each band. If this is the case, there is not a uniform gap at the Fermi level. The superconducting state resembles, in this respect, that of an anisotropic superconductor. The effects of interband scattering on the bulk properties of a superconductor with two different bands at the Fermi level was studied in.

Interband scattering can be induced by any defect which breaks the translational symmetry of the lattice, including lattice imperfections, phonons, and surfaces or grain boundaries. Thus, in materials with many bands crossing the Fermi level, one can expect the superconducting properties near a surface to differ from those at the bulk. This effect should be more pronounced in systems with many irregular interfaces, such as ceramic and granular materials.

In the present work, we study pair breaking effects at the surfaces of many band superconductors. The main difference with the theory developed in is that, when interband scattering is restricted to the surface region, the superconducting gaps are inhomogeneous, because the value of the gaps within the bulk of the system is not affected by the presence of a surface. In order to take this effect into account, we develop a method which differs from the standard treatment of pair breaking effects in conventional or unconventional superconductors. We use the simplest possible model where non trivial effects are expected, presented in the next section. Section III presents the main results, and the main conclusions which can be drawn from them are discussed in section IV.

## II The model

### A The hamiltonian.

Superconductors with different bands at the Fermi level, where each of these bands have different pairing interactions have already been discussed in the literature. The simplest model contains two bands, with two different densities of states and pairing interactions. We will consider the hamiltonian:

$``$ $`=`$ $`_0+_{int}+_{pb}`$ (1)
$`_0`$ $`=`$ $`{\displaystyle \underset{i=1,2;\stackrel{}{𝐤},s}{}}ϵ_{i,\stackrel{}{𝐤}}c_{i,\stackrel{}{𝐤},s}^{}c_{i,\stackrel{}{𝐤},s}`$ (2)
$`_{int}`$ $`=`$ $`{\displaystyle \underset{i=1,2;\stackrel{}{𝐤}\stackrel{}{𝐤}^{}}{}}g_ic_{i,\stackrel{}{𝐤}}^{}c_{i,\stackrel{}{𝐤}}^{}c_{i,\stackrel{}{𝐤}^{}}c_{i,\stackrel{}{𝐤}^{}}`$ (3)
$``$ $`{\displaystyle \underset{\stackrel{}{𝐤}\stackrel{}{𝐤}^{}}{}}g^{}c_{1,\stackrel{}{𝐤}}^{}c_{1,\stackrel{}{𝐤}}^{}c_{2,\stackrel{}{𝐤}^{}}c_{2,\stackrel{}{𝐤}^{}}`$ (4)
$`_{pb}`$ $`=`$ $`{\displaystyle \underset{s}{}}{\displaystyle d^3\stackrel{}{𝐫}Vf(z)\psi _{1,s}^{}(\stackrel{}{𝐫})\psi _{2,s}(\stackrel{}{𝐫})}`$ (5)
Electrons within each band experience a different pairing interaction, $`g_i`$, leading to two superconducting gaps. The two bands are coupled by the interaction $`g^{}`$. Otherwise, in the absence of interband scattering the two gaps open at different temperatures, $`T_{c,i}\omega _0\mathrm{exp}(W_i/g_i)`$, where $`W_i`$ is the bandwidth, and $`\omega _0`$ is a cutoff related to the pairing mechanism. Specific heat measurements seem to exclude this possibility in MgB2. We neglect intraband scattering, which does not give rise to pair breaking effects, at least to lowest order. Finally, we assume a constant interband scattering term, localized near the surface. The function $`f(z)`$ is peaked near the surface, located at $`z=0`$. The width of $`f(z)`$ is of the order of a lattice spacing. With these restrictions, the model described in eq.(5) includes six parameters, with dimensions of energy: the density of states at the Fermi level of each band, $`N_i(ϵ_F)`$, the pairing interactions, $`g_i`$ and $`g^{}`$, and the scattering potential, $`V`$.

### B Pair breaking effects.

The lack of translational symmetry induced by the surface makes it convenient to solve directly the the Bogoliubov-de Gennes equations derived from (5). We use a discrete tight binding model for this purpose. We assume that each band can be described by a single orbital per site, and that there are local attractive interactions with induce the pairing. The model reduces, in the absence of interband scattering, to two coupled negative $`U`$ Hubbard models. In this basis, interband scattering can be included by allowing for hopping from one orbital to the other at any given lattice site. Our discretized model in real space becomes:

$``$ $`=`$ $`_0+_{int}+_{pb}`$ (6)
$`_0`$ $`=`$ $`{\displaystyle \underset{l;ijs}{}}t_lc_{l,i,s}^{}c_{l,j,s}+\mathrm{h}.\mathrm{c}.`$ (7)
$`_{int}`$ $`=`$ $`{\displaystyle \underset{l,i}{}}U_lc_{l,i,}^{}c_{l,i,}c_{l,i,}^{}c_{l,i,}`$ (8)
$``$ $`{\displaystyle \underset{i}{}}U^{}c_{1,i,}^{}c_{1,i,}c_{2,i,}^{}c_{2,i,}`$ (9)
$`_{pb}`$ $`=`$ $`{\displaystyle \underset{s,ii_s}{}}Vc_{1,i,s}^{}c_{2,i,s}`$ (10)
We assume that the lattice is a semiinfinite chain, and that interband scattering is restricted to the outermost site, as schematically depicted in Fig..

The BCS equations are solved using the standard equivalence of the attractive Hubbard model in a bipartite lattice to the repulsive Hubbard model. The semiinfinite model is solved using transfer matrix techniques, which are described in the Appendix. The gaps are calculated selfconsistently in a layer of $`m`$ sites, which are coupled to an homogeneous chain where the gaps take the bulk values, where the gaps are defined as:

$$\mathrm{\Delta }_i=\frac{U_i}{2}c_i^{}c_i^{}+\frac{U^{}}{2}c_j^{}c_j^{}$$
(11)
We present results obtained with $`m=32`$. This idealized one dimensional model can be viewed as an approximation to the inhomogeneous layered structure expected near the surface. Note that the BCS equations exclude the possibility of one dimensional fluctuations, so that the solutions to be discussed below do not show unphysical one dimensional features.

## III Results

We use as unit of energy one of the hopping parameters in eq.(10), $`t_1=1`$. We assume that band 1 represents the total effect of the two $`\sigma `$ bands in MgB2. The density of states in the $`\pi `$ band is roughly one third that of the two $`\sigma `$ bands, and we set $`t_2=3t_1`$. We assume that the superconductivity is mostly due to the $`\sigma `$ band, and we choose the parameters $`U_1,U_2`$ and $`U^{}`$ so that the smallest gap is less than one half the largest one, at zero temperature. A reasonable combination is $`U_1=0.5t_1,U_2=0`$ and $`U^{}=0.2t_1`$ so that $`\mathrm{\Delta }_1=0.0212t_1`$ and $`\mathrm{\Delta }_2=0.0079t_1`$. The critical temperature in these units is $`T_c=0.011t_1`$. Using the experimental value $`T_c40`$K, the value of the unit of energy, $`t_1`$ is equal to 0.3 eV, and the density of states of the $`\sigma `$ band is 0.5 states/eV-cell, which has the right order of magnitude. Note that, with these parameters, the model is well into the BCS weak coupling regime. Finally, we must determine the strength of the interband scattering at the surface, $`V`$. The mixing of the boron $`\pi `$ and $`\sigma `$ bands is strongly suppressed in the bulk. This effect is weakened if the relative angles between nearest neighbor B ions is modified at the surface, so that the $`\sigma `$ bands are not built up of sp2 orbitals only. The interband hopping associated with these deformations will be a fraction of the bulk hopping terms. We take in the following $`V=0.4t_10.12`$eV, and restrict this interband hopping to the outermost layer. The gaps near the surface, at zero temperature, are plotted in the inset in Fig.. The gaps are almost uniform, suggesting that the influence of the interband scattering on the value of the gaps is very small.

Interband scattering has a stronger influence on the density of states near the surface. The results for the outermost site are plotted in Fig.. The coupling between the two bands shows in the existence of peaks in the density of states at the two gap positions, while, in the bulk, each band displays a single peak at the value of the corresponding gap.

The overall features in the density of states remain the same at relatively high temperatures, as shown in fig., where the results at T=36K are shown. The scale at which the smallest gap closes is determined by the largest gap.

In some experiments, like point contact spectroscopy, the probe can be an additional source of interband scattering, at the position where the measurement is being made. In the presence of strong interband scattering at the surface, the perturbation in the densities of states is more pronounced, as shown in Fig., calculated using $`V=t_10.3`$eV. A single smeared gap will be observed in an experiment of this type.

It is interesting to note that, if the value of the gaps in the bulk were of opposite sign, a midgap state, induced by Andreev reflections, should arise. The problem considered here maps onto that of a dimerized chain. If the value of the dimerization changes sign, a gap state always arises, as extensively discussed in connexion to solitons in polyacetilene.

## IV Conclusions

We have analyzed the effects of interband scattering at the surface, in a superconductor with a many bands at the Fermi level, and different pairing strengths for each band.

Surface scattering, for reasonable values of the parameters, is inneffective in changing the gaps near the surface. This is consistent with the fact that, in a conventional weak coupling BCS superconductor, the coherence length is much greater than the surface layer where strong scattering is expected.

The influence of the surface scattering induces significant changes in the density of states near the surface. The density of states of each band shows peaks at the positions of the gaps associated with the other bands. There is, however, always a minimum gap, below which the density of states vanishes. This gap, $`\mathrm{\Delta }_{min}`$, has a weak temperature dependence until close to the bulk critical temperature.

We have used a set of parameters appropiate for MgB2. We assume that there is a wide, delocalized band, and a narrower and more localized band, which determines mostly the superconducting properties. The wide band (derived from the $`\pi `$ orbitals) has the smallest gap, and is weakly influenced by the narrow band. The narrow band shows stronger features at both the small and large gaps. Tunneling experiments are highly sensitive to the delocalization of the wavefunction. Hence, it is possible that they measure the density of states associated to the wider band ($`\pi `$ orbitals). If this is the case, the observations will show a single gap, well approximated by the BCS expression. This gap should have a weak temperature dependence, until temperatures comparable to the bulk critical temperature.

The strength of the interband scattering at the surface studied here can depend on the experimental setup, and it may be enhanced in some experiments. If that is the case, a single smeared gap will be observed, of magnitude comparable to 1.7 Tc.

Finally, it is interesting to note that, if the gaps in the two bands were of opposite signs, as expected from electron-electron pairing mechanisms, a surface state near the center of the gap should appear.

## V Acknowledgements

We are thankful to S. Vieira, H. Suderow and G. Rubio-Bollinger for many helpful discussions. Financial support from CICyT (Spain) through grant PB96/0875 is gratefully acknowledged.

## VI Appendix

The problem defined by the hamiltonian in eq.(10) with the geometry shown in Fig. reduces to the calculation of the density of states in $`2m`$ sites, with gap values $`\mathrm{\Delta }_i`$ and connected by hopping terms which can be either $`t_1,t_2`$ or $`V`$. The attractive Hubbard model, $`U_i<0`$, can be mapped onto the repulsive Hubbard model, in a bipartite lattice, by the transformation:

$`c_i^{}`$ $``$ $`d_i^{}`$ (12)
$`c_i`$ $``$ $`d_i`$ (13)
$`c_i^{}`$ $``$ $`(1)^id_i`$ (14)
$`c_i`$ $``$ $`(1)^id_i^{}`$ (15)
The superconducting order parameter, $`\mathrm{\Psi }_i=c_i^{}c_i^{}`$, is mapped onto a staggered magnetization in the transverse direction. At half filling, there is an additional symmetry which allows us to rotate this magnetization to the z-axis. Then, $`\mathrm{\Psi }_i=(1)^i(c_i^{}c_ic_i^{}c_i)`$. In this representation, the hamiltonian does not mix the spins, and it can be decomposed into two boxes, one for each spin direction. The problem is reduced to the calculation of the density of states in a tight binding chain with variable hoppings, $`t_{i,i\pm 1}`$, and energy levels, $`ϵ_i`$, which are related to the local value of the gaps. The gaps must be determined selfconsistency from the values of the $`\mathrm{\Psi }_i`$’s.

The fractions:

$$g_{n,n\pm 1}(\omega )=t_{n,n\pm 1}\frac{G_{n,n^{}}(\omega )}{G_{n\pm 1,n^{}}(\omega )}$$
(16)
are independent of $`n^{}`$, and satisfy:

$$g_{n,n\pm 1}(\omega )=\frac{t_{n,n\pm 1}^2}{\omega ϵ_{n,n\pm 1}t_{n\pm 1,n\pm 2}g_{n\pm 1,n\pm 2}(\omega )}$$
(17)
and:

$$G_{n,n}(\omega )=\frac{1}{\omega ϵ_nt_{n,n1}g_{n,n1}(\omega )t_{n,n+1}g_{n,n+1}(\omega )}$$
(18)
Thus, the problem can be solved by iteration from the boundaries, provided that one knows the values of $`g_{\pm m,\pm m\pm 1}`$. These values can be easily be calculated, if one assumes that the values of the $`ϵ_i`$’s are constant beyond position $`m`$. Finally, the selfconsistency requirement for the values of $`ϵ_i,i=1,2m`$ must be satisfied.


text:
# On the thin-string limit of the 6d stringlike defect model

## Proof of Eqn. (15)

It remains to be shown that taking constant $`R_0`$ implies (15). The proof is straightforward. For the sake of clarity of notation, let us denote $`\alpha =\sigma ^{3/4}\gamma ^{1/2}`$. In terms of $`\alpha `$, the boundary conditions at $`\rho =0`$ and the junction condition read:

$$\alpha |_{\rho =0}=0,\alpha ^{}|_{\rho =0}=A^{3/4},\mathrm{and}\alpha |_{\rho =ϵ}=R_0e^{5/4cϵ}$$
(17)
Assuming that $`\alpha `$ is a continuously differentiable function, it is clear that there must exist a $`\rho ^{}<ϵ`$ such that:

$$\alpha ^{}(\rho ^{})=\frac{\alpha (ϵ)\alpha (0)}{ϵ}=\frac{R_0e^{5/4cϵ}}{ϵ}.$$
(18)
This implies, assuming further that $`\alpha ^{}`$ is continuously differentiable in $`(0,\rho ^{})`$, that there exists a $`\rho ^{}<\rho ^{}`$ such that:

$`\alpha ^{\prime \prime }(\rho ^{})`$ $`=`$ $`{\displaystyle \frac{\alpha ^{}(\rho ^{})\alpha ^{}(0)}{\rho ^{}}}={\displaystyle \frac{R_0e^{5/4cϵ}ϵA^{3/4}}{\rho ^{}ϵ}}`$ (19)
$``$ $`{\displaystyle \frac{R_0e^{5/4cϵ}ϵA^{3/4}}{ϵ^2}}.`$ (20)
Suppose first that

$$\alpha (\rho ^{})<R_0e^{5/4cϵ}.$$
(21)
The equation (19) implies:

$$\frac{\alpha ^{\prime \prime }(\rho ^{})}{\alpha (\rho ^{})}>\frac{1}{ϵ^2}\left[1ϵ\frac{A^{3/4}}{R_0}e^{5/4cϵ}\right].$$
(22)
For $`ϵ0`$ and $`A,R_0𝒪(1)`$, the first term in this inequality becomes dominant and therefore:

$$\frac{\alpha ^{\prime \prime }(\rho ^{})}{\alpha (\rho ^{})}𝒪\left(\frac{1}{ϵ^2}\right).$$

Suppose now that

$$\alpha (\rho ^{})R_0e^{5/4cϵ}.$$

Then we can repeat the same reasoning to find an another point, now inside the interval $`(0,\rho ^{})`$, satisfying an inequality analogous to (19). If this new point still does not satisfy (21), we can repeat this procedure over and over until we find a point satisfying both (19) and (21). Given that $`\alpha `$ is continuously differentiable and that $`\alpha (0)=0`$, this iteration must eventually stop and such a point must exist.

In conclusion, there exists a point $`\rho ^\mathrm{\#}(0,ϵ)`$ where, for $`ϵ`$ sufficiently small:

$$\frac{\alpha ^{\prime \prime }(\rho ^\mathrm{\#})}{\alpha (\rho ^\mathrm{\#})}𝒪\left(\frac{1}{ϵ^2}\right).$$

(In fact, given the regularity properties of $`\alpha `$, it will be also valid for the points in a small neighborhood of $`\rho ^\mathrm{\#}`$.)


text:
# Multispin Ising spin glasses with ferromagnetic interactions

## 1 Introduction

Competition between quenched randomness, frustration, and a tendency to order uniformly causes non-trivial behaviour in many-body systems. Typical examples include spin glasses with ferromagnetic bias, neural networks, proteins, and error-correcting codes. When the quenched randomness dominates, these systems exhibit randomly frozen (spin glass) states. When the ordering tendency dominates, they exhibit uniformly ordered (non-glassy ferromagnetic) states. In the intermediate regime they can exhibit mixed (glassy ferromagnetic) states. The phase diagrams around this crossover often have a rich structure.

One of the principal purposes of this paper is to present a new method to construct the phase diagrams of such systems from an understanding of the related system with only the random term and an external field term. We apply this method to an Ising spin glass with a $`p`$-spin Gaussian random interaction and an $`r`$-spin ordering interaction, both of infinite range. It has been shown using replica theory that for $`p>2`$, the system with neither ordering nor field terms has two phase transitions: a discontinuous one step replica symmetry breaking, followed by a full replica symmetry breaking at lower temperature. Subsequent work investigated the one step transition with the addition of either a field or an ordering interaction with $`r=p`$. We extend the results to the case of general $`r`$ and also examine the transition from one step to full replica symmetry breaking. This work follows on from a previous paper in which we dealt not with an Ising spin system, but with a technically simpler spherical spin system.

The order of the paper is as follows. In §2 we present the general method. In §3 we review the replica theory of the Ising $`p`$-spin glass. In §4 we discuss the phase diagrams of the Ising $`p`$-spin glass with $`r`$-spin ferromagnetism. In §5 we present some exact results using a gauge transformation for the case where $`r=p`$. Finally, we make concluding remarks and present technical appendices.

## 2 The constrained magnetization approach to mean field ferromagnetism

We consider a general system of scalar spins with a Hamiltonian

$$_r(J_0)=_{\mathrm{base}}\frac{J_0(r1)!}{N^{r1}}\underset{i_1<i_2\mathrm{}<i_r}{}\sigma _{i_1}\mathrm{}\sigma _{i_r}$$
(1)
where the first term is the Hamiltonian of any system (referred to as the base system), and the second term is an infinite-range $`r`$-spin ferromagnetic interaction. To leading order in $`1/N`$ the partition function is

$`Z_r(J_0)`$ $`={\displaystyle \underset{\{\sigma _i\}}{}}\mathrm{exp}\left[\beta _{\mathrm{base}}+N\beta J_0{\displaystyle \frac{1}{r}}\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\sigma _i\right)^r\right]`$ (2)
$`={\displaystyle 𝑑M\mathrm{exp}}N\beta \left(f_{\mathrm{base}}(M){\displaystyle \frac{1}{r}}J_0M^r\right)`$ (3)
where
$`f_{\mathrm{base}}(M)`$ $`={\displaystyle \frac{1}{N\beta }}\mathrm{ln}Z_{\mathrm{base}}(M),`$ (4)
$`Z_{\mathrm{base}}(M)`$ $`={\displaystyle \underset{\{\sigma _i\}}{}}\delta \left(M{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\sigma _i\right)\mathrm{exp}\beta _{\mathrm{base}}.`$ (5)
We recognize $`f_{\mathrm{base}}(M)`$ as the free energy per site of the base system when the magnetization per spin is constrained (via the $`\delta `$-function in the partition function) to be $`M`$. Given a knowledge of this function, we can carry out the integral over $`M`$ in (3) to leading order in $`1/N`$ by the saddle point method, giving a free energy per spin of

$$f_r(J_0)=\frac{1}{N\beta }\mathrm{ln}Z_r(J_0)=f_{\mathrm{base}}(M_r)\frac{1}{r}J_0M_r^r$$
(6)
where the equilibrium magnetization $`M_r(J_0)`$ minimizes the right hand side, and so is a solution to the equation

$$f_{\mathrm{base}}^{}(M_r)=J_0M_r^{r1}.$$
(7)
We note further that if we take $`r=1`$, the second term of (1) represents an external magnetic field $`h=J_0`$, and so the free energy of the constrained base system is related simply to that of the unconstrained base system in an external magnetic field by

$$f_{\mathrm{field}}(h)=f_{\mathrm{base}}(M_{\mathrm{field}}(h))hM_{\mathrm{field}}(h),$$
(8)
$$f_{\mathrm{base}}^{}(M_{\mathrm{field}}(h))=h.$$
(9)
Therefore a knowledge of the free energy of either system is enough to determine that of the ferromagnetic system via (6) and (7): the ferromagnetic term acts as an effective field

$$h=J_0M^{r1}$$
(10)
which must be determined self-consistently.

## 3 The replica theory of the Ising $`p`$-spin glass

In the remainder of this paper we shall concentrate on the example of an infinite-ranged Ising $`p`$-spin glass model. The Hamiltonian we use is

$$_r=\underset{i_1<i_2\mathrm{}<i_p}{}J_{i_1\mathrm{}i_p}\sigma _{i_1}\mathrm{}\sigma _{i_p}\frac{J_0(r1)!}{N^{r1}}\underset{i_1<i_2\mathrm{}<i_r}{}\sigma _{i_1}\mathrm{}\sigma _{i_r}$$
(11)
where $`J_{i_1\mathrm{}i_p}`$ are independent quenched random couplings given by a Gaussian distribution with zero mean and variance $`p!J^2/2N^{p1}`$, and $`\sigma _i=\pm 1`$ are Ising spins. For $`p=2`$ we recover a generalized Sherrington–Kirkpatrick model which has only full replica symmetry breaking. In the present paper we consider $`p>2`$ so that there also is a region of one step replica symmetry breaking. In the limit $`p\mathrm{}`$, we recover the random energy model (REM) which has no full replica symmetry breaking: this will be discussed further in §4.

In the spirit of the previous section, we consider first the simpler

$$_{\mathrm{field}}=\underset{i_1<i_2\mathrm{}<i_p}{}J_{i_1\mathrm{}i_p}\sigma _{i_1}\mathrm{}\sigma _{i_p}h\underset{i}{}\sigma _i.$$
(12)
Previous studies used the replica method to obtain the free energy per site $`f_{\mathrm{field}}`$ given by

$$\begin{array}{c}\beta f_{\mathrm{field}}=\underset{n0}{lim}\frac{1}{n}\left[\frac{1}{4}(p1)\beta ^2J^2\underset{ab}{}q_{ab}^p\mathrm{ln}\underset{\{\sigma ^a\}}{}\mathrm{exp}\left(\frac{1}{4}p\beta ^2J^2\underset{ab}{}q_{ab}^{p1}\sigma ^a\sigma ^b+\beta h\underset{a}{}\sigma ^a\right)\right]\frac{\beta ^2J^2}{4}\hfill \end{array}$$
(13)
which must be extremized for the replica overlaps $`q_{ab}=\frac{1}{N}_i\sigma _i^a\sigma _i^b`$ for $`1abn`$. In the replica symmetric (RS) ansatz ($`q_{ab}=q`$ for all $`ab`$) this is

$$\begin{array}{c}\beta f_{\mathrm{field}}^{(\mathrm{RS})}=\frac{1}{4}(p1)\beta ^2J^2q^p\frac{1}{4}\beta ^2J^2+\frac{1}{4}p\beta ^2J^2q^{p1}\hfill \\ \hfill \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{ln}2\mathrm{cosh}\left(\sqrt{\frac{1}{2}p\beta ^2J^2q^{p1}}z+\beta h\right)\end{array}$$
(14)
where the overlap $`q`$ is given by the self-consistency equation

$$q=\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2\left(\sqrt{\frac{1}{2}p\beta ^2J^2q^{p1}}z+\beta h\right).$$
(15)
This RS solution is stable against small replica symmetry breaking Almeida–Thouless fluctuations if

$$\frac{1}{2}p(p1)\beta ^2J^2q^{p2}\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sech}^4\left(\sqrt{\frac{1}{2}p\beta ^2J^2q^{p1}}z+\beta h\right)<1.$$
(16)
The magnetization of the RS solution is

$$M=\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}\left(\sqrt{\frac{1}{2}p\beta ^2J^2q^{p1}}z+\beta h\right).$$
(17)
In the one step replica symmetry breaking (1RSB) ansatz the free energy per site is $`f_{\mathrm{field}}^{(1\mathrm{R}\mathrm{S}\mathrm{B})}`$ given by

$$\begin{array}{c}\beta f_{\mathrm{field}}^{(1\mathrm{R}\mathrm{S}\mathrm{B})}=\frac{1}{4}(p1)\beta ^2J^2\left[(1x)q_1^p+xq_0^p\right]\frac{1}{4}\beta ^2J^2+\frac{1}{4}p\beta ^2J^2q_1^{p1}\hfill \\ \hfill \frac{1}{x}\frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\mathrm{ln}2^x\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG\end{array}$$
(18)
where

$$G\sqrt{\frac{1}{2}p\beta ^2J^2q_0^{p1}}z_0+\sqrt{\frac{1}{2}p\beta ^2J^2(q_1^{p1}q_0^{p1})}z_1+\beta h,$$
(19)
and the mutual overlap $`q_0`$, the self overlap $`q_1`$, and the break point $`x`$ (i.e. the weight of $`q_0`$ in the averaged overlap distribution function) are given by the self-consistency equations

$`q_0`$ $`={\displaystyle \frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\left(\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG\mathrm{tanh}G}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG}\right)^2},`$ (20a)
$`q_1`$ $`={\displaystyle \frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG\mathrm{tanh}^2G}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG}},`$ (20b)
$`{\displaystyle \frac{1}{4}}(p1)\beta ^2J^2(q_1^pq_0^p)`$ $`={\displaystyle \frac{1}{x^2}}{\displaystyle \frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\mathrm{ln}\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG}`$
$`+{\displaystyle \frac{1}{x}}{\displaystyle \frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG\mathrm{ln}\mathrm{cosh}G}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG}}.`$ (20c)
This 1RSB solution is stable against further replica symmetry breaking fluctuations if

$$\frac{1}{2}p(p1)\beta ^2J^2q_1^{p2}\frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^{x4}G}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG}<1.$$
(21)
Where the 1RSB solution becomes unstable, we expect a continuous transition to a solution given by the full Parisi ansatz, which we do not consider in detail here. The magnetization of the 1RSB solution is

$$M=\frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG\mathrm{tanh}G}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}^xG}.$$
(22)
## 4 Full results and phase diagrams

The phase diagrams for the Hamiltonian (11) were found in two stages. They are shown for $`p=5`$ and $`r=1,\mathrm{},6`$ in figure 1.

Firstly, the free energy and magnetization of the system in a field but without ferromagnetism, as given by (12), were found. The RS solution was found by numerical solution of (15), the 1RSB by numerical solution of (20). No attempt was made to determine the solution with full replica symmetry breaking (FRSB) in the Parisi scheme, but the onset of unstable modes in the 1RSB solution, as given by (21), is taken to signify a continuous transition to FRSB.

The resulting phase diagram for $`p=5`$ is plotted in figure 1(a). This Hamiltonian is equivalent to (11) with $`p=5`$, $`r=1`$, and $`J_0=h`$. There are three phases: an RS paramagnet, a 1RSB spin glass, and an FRSB spin glass. The replica symmetry breaking transition curve has been found previously. Where the field is less than a critical value, there is a discontinuous one step RSB transition (D1RSB) where $`x=1`$, and the transition temperature rises with $`h`$. Above the critical value, the transition is continuous, first continuous one step (C1RSB) where $`q_0=q_1`$, and then FRSB; and the transition temperature falls with $`h`$, vanishing asymptotically as $`h\mathrm{}`$. (Note that the D1RSB transition is continuous in the thermodynamic sense.) The numerical calculation of the 1RSB/FRSB curve is troublesome near the C1RSB transition line: this is because the difference between the left and right hand sides of (21) is very small, of the order of $`(q_1q_0)^4`$. Consequently we do not have reliable results in this region, and a small gap is left in the plot to indicate this. However, the point where this curve meets the C1RSB line (which we shall call the RSB triple point, since RS, 1RSB, and FRSB phases meet at it) has been determined accurately using the RS solution and a perturbative calculation, details of which are given in Appendix A. The magnetization $`M`$ was calculated with (17) and (22), and curves of $`h`$ as a function of $`M`$ at fixed temperature for $`p=5`$ are shown in figure 2. We indicate where the solution is of each type (RS, 1RSB, FRSB), but note that $`h`$ is not exact in the FRSB region, since it is approximated by the 1RSB value.

Secondly, the methods of §2 were used to calculate the transitions in the full Hamiltonian, essentially by finding the stable solutions of (10) with the $`h(M)`$ found previously. (In the FRSB region, the use of the 1RSB values introduces a small error in the ferromagnetic transition lines.) The resulting phase diagrams for $`p=5`$ and $`r=2,\mathrm{},6`$ are plotted in figures 1(b)1(f).

For $`r=2`$ and $`p>2`$, the phase diagram contains six regions. For small $`J_0`$, there are three $`M=0`$ phases: on cooling, the RS paramagnet undergoes first a D1RSB transition, and then an FRSB (not shown explicitly in figure 1). Increasing $`J_0`$, the $`M=0`$ phases become unstable, and a continuous transition to $`M0`$ occurs. There are three ferromagnetic phases: RS, 1RSB, and FRSB. Within the ferromagnetic region, the RS/1RSB curve has D1RSB and C1RSB parts, as for $`r=1`$. (The RS/1RSB transition in the $`M=0`$ region is always D1RSB, since the effective field $`J_0M^{r1}`$ vanishes.)

For $`r>2`$, the ferromagnetic transition is first order: on increasing $`J_0`$, there is first a spinodal transition, at which a metastable $`M0`$ solution appears, followed by a thermodynamic transition, at which the comparative free energies make it the favoured state.

Increasing $`r`$ at fixed $`p`$, we identify two qualitative changes, illustrated in figure 3. For $`r`$ slightly greater than $`2`$, the diagram has the form of figure 3(a). On increasing $`r`$, the peak of the RS/1RSB curve approaches the spinodal ferromagnetic transition line, and the D1RSB curve between them shortens. The peak reaches the spinodal line for $`r=r_1`$, and there is no D1RSB curve in the ferromagnetic region for $`rr_1`$, as in figure 3(b). The RSB triple point approaches the spinodal ferromagnetic transition line on increasing $`r`$, and the 1RSB region shrinks. The triple point reaches the spinodal line at $`r=r_2`$, and there is no ferromagnetic 1RSB region for $`rr_2`$, as in figure 3(c). These special values of $`r`$ may be evaluated for given $`p`$ with some accuracy, as discussed in Appendix B. Some results are shown in table 1. They are in accordance with an exact result, to be proved in §5, that $`r_1`$ cannot exceed $`p`$.

Of particular interest is the behaviour of the model in the case $`r=p`$; this is relevant to error-correcting codes, and is amenable to an exact calculation on the Nishimori line. We observe that there is a first order transition to a ferromagnetism which is glassy for small enough values of $`J_0`$ and $`T`$, and that there is a triple point where the spinodal ferromagnetic transition line and the glassy transition line within the ferromagnet meet with infinite slope. (Exact results in the $`r=p`$ case, including the fact that the solution is RS on the line $`p\beta J^2=2J_0`$, are given in §5.) The RSB region shrinks as $`r=p`$ increases. The existence of an RS ferromagnetic phase persists in the limit $`r=p\mathrm{}`$. Our results therefore contradict conclusions of Dorlas and Wedagedera, who claim that the ferromagnetic region in the model with $`r=p\mathrm{}`$ is entirely glassy. This discrepancy is due to the order in which the limits $`rp`$ and $`p\mathrm{}`$ are taken. The approach of Dorlas and Wedagedera is to consider the REM with general $`r`$, effectively taking the limit $`p\mathrm{}`$ in the Ising model, and to subsequently take the limit $`r\mathrm{}`$. In doing this, they obtain results which would apply to $`pr\mathrm{}`$. In our approach we set $`r=p`$ and then take $`p\mathrm{}`$, obtaining results which apply to $`r=p\mathrm{}`$, the appropriate case when one considers error-correcting codes. This difference may be illustrated using figures 1(c) and 1(e): the former is for $`r<p`$, and the ferromagnetic phase is dominated by the 1RSB region in the limit $`p\mathrm{}`$ with $`r<p`$; the latter is for $`r=p`$, and the ferromagnetic phase is dominated by the RS region in the limit $`r=p\mathrm{}`$.

## 5 Special features on the Nishimori line

Some of our conclusions about the phase diagram for $`r=p`$ are reinforced from a different point of view. We show in the present section that on the Nishimori line, the distribution function of the spin glass order parameter consists only of two delta functions. This immediately implies that the ferromagnetic phase is not glassy on the Nishimori line for any $`r=p`$.

We start our argument from the definition of the distribution function of magnetization under the Hamiltonian (11) with $`r=p`$:

$$P_m(y)=\left[\frac{_{\{\sigma _i\}}\delta (yN^1_i\sigma _i)\mathrm{exp}\left(\beta _p\right)}{_{\{\sigma _i\}}\mathrm{exp}\left(\beta _p\right)}\right]_J.$$
(23)
Here the outer brackets denote the configurational average over the Gaussian distribution of exchange interactions. To define the distribution function of the spin glass order parameter, it is convenient to introduce two replicas of the same system, with spins $`\{\sigma _i^{(1)}\}`$ and $`\{\sigma _i^{(2)}\}`$:

$$P_q(y)=\left[\frac{_{\{\sigma _i^{(1)}\}}_{\{\sigma _i^{(2)}\}}\delta (yN^1_i\sigma _i^{(1)}\sigma _i^{(2)})\mathrm{exp}\left(\beta _p(\{\sigma _i^{(1)}\})\beta _p(\{\sigma _i^{(2)}\})\right)}{_{\{\sigma _i^{(1)}\}}_{\{\sigma _i^{(1)}\}}\mathrm{exp}\left(\beta _p(\{\sigma _i^{(1)}\})\beta _p(\{\sigma _i^{(2)}\})\right)}\right]_J.$$
(24)
Our main result is the identity

$$P_m(y)=P_q(y)\text{wherever}p\beta J^2=2J_0.$$
(25)
The condition $`p\beta J^2=2J_0`$ defines a curve in the phase diagram called the Nishimori line. This line passes through the paramagnetic and ferromagnetic phases as shown by the dot-dashed curve in figure 4. The result is proved in Appendix C.1.

It is well established that the distribution of magnetization $`P_m(y)`$ has a simple structure with two delta functions at symmetric positions. Then (25) proves that the distribution function of the spin glass order parameter $`P_q(y)`$ is also simple on the Nishimori line, which is not the case if the ferromagnetic phase is glassy.

In Appendix C.2 we demonstrate three further facts about the Nishimori line: it passes through the point where the spin glass transition line meets the thermodynamic ferromagnetic transition line (the lower-right point marked in figure 4); at this point, the $`q_1`$ of the spin glass and the $`q`$ of the ferromagnet are equal; and it passes through the point where the continuous glassy transition line (or, equivalently, the RS Almeida–Thouless line) within the ferromagnetic phase meets the spinodal ferromagnetic transition line (the upper-left point marked in figure 4). This last result demonstrates that there is no D1RSB transition line within the ferromagnetic region at $`r=p`$, as the continuous glassy transition reaches the spinodal line which marks the edge of the ferromagnetic phases. It follows that $`r_1p`$, as claimed in §4.

## 6 Conclusions

In this paper, we have shown that the behaviour of a spin glass with an $`r`$-spin ferromagnetic interaction of infinite range is simple to obtain from a knowledge of the free energy either of the base spin glass with a constrained magnetization, or equivalently of the base spin glass in an external field (although these latter problems are often numerically subtle); and we have applied this method to an Ising glass with a $`p`$-spin Gaussian random interaction. We have found that, for $`p>2`$ and $`r2`$, there are both glassy and non-glassy ferromagnetic phases; and that for $`r`$ less than a $`p`$-dependent critical value the glassy ferromagnet has regions of both one step and full replica symmetry breaking, while for higher $`r`$ it has only full replica symmetry breaking. We have also proved several exact results for the case where $`r=p`$, which is relevant to error-correcting codes: in particular, we have explicitly shown the existence of a curve of simple form in the phase diagram, passing through both paramagnetic and ferromagnetic phases, on which the solution is certainly non-glassy. We hope that the present analysis will stimulate investigation of the phase structures of related systems.

## Appendix A Results from perturbation theory close to the Almeida–Thouless line

In this appendix we derive expression for $`x_\mathrm{c}`$, the value of $`x`$ along the C1RSB curve, and the condition to determine the RSB triple point where the RS, 1RSB, and FRSB ferromagnetic phases meet. We use the definitions

$`𝒯_n`$ $`{\displaystyle \frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\mathrm{tanh}^n\left(\sqrt{\frac{1}{2}p\beta ^2J^2q_0^{p1}}z_0+\beta h\right)},`$ (26)
$`𝒮_n`$ $`{\displaystyle \frac{dz_0}{\sqrt{2\pi }}e^{z_0^2/2}\mathrm{sech}^n\left(\sqrt{\frac{1}{2}p\beta ^2J^2q_0^{p1}}z_0+\beta h\right)}.`$ (27)
A continuous one step replica symmetry breaking (C1RSB) transition occurs on the line where $`q_0=q_1`$ in the 1RSB ansatz. Details of the 1RSB solution close to this line, using only the numerical solution to the RS equation, are obtained by an expansion of these equations in small $`ϵ(q_1q_0)`$. Since the hyperbolic functions of $`G`$ become polynomials in $`z_1`$, the integrals over $`z_1`$ can be carried out immediately, simplifying matters considerably. In particular, one can find the RSB triple point using only the numerical RS data.

To leading order, both (20a) and (20b) are

$$q_0=𝒯_2+O(ϵ)$$
(28)
and so $`q_0=q+O(ϵ)`$ where $`q`$ is the RS order parameter satisfying (15). To obtain a second piece of information from these two equations, $`[(\text{20b})(\text{20a})]`$ must be expanded to first order, giving

$$ϵ=\frac{1}{2}p(p1)\beta ^2J^2q_0^{p2}𝒮_4ϵ+O(ϵ)^2.$$
(29)
This shows that (16) holds as an equality to leading order in $`ϵ`$, and hence that the C1RSB transition and the onset of RS AT instability coincide. This is to be expected: where the replica symmetry breaking mode of the RS solution softens, a continuous transition to an RSB solution occurs.

(20c) is trivially satisfied for $`q_0=q_1`$. Expanding to first order gives (28) multiplied by $`k_1\frac{1}{4}p(p1)\beta ^2J^2q_0^{p2}ϵ`$, which gives no new information. Expanding $`[(\text{20c})k_1(\text{20a})]`$ to second order gives (28) multiplied by $`k_2\frac{1}{8}p(p1)(p2)\beta ^2J^2q_0^{p3}ϵ^2`$ plus (29) multiplied by $`\frac{1}{2}k_1`$, which again gives no new information. Expanding $`[(\text{20c})(k_1+k_2)(\text{20a})\frac{1}{2}k_1((\text{20b})(\text{20a}))]`$ to third order gives

$$\frac{1}{12}\left[\frac{1}{2}p(p1)\beta ^2J^2q_0^{p2}\right]^3𝒮_6(xx_\mathrm{c})ϵ^3+O(ϵ)^4=0$$
(30)
$$x=x_\mathrm{c}+O(ϵ)$$
(31)
where
$$x_\mathrm{c}\frac{2(p2)q_0^3+4p(p1)\beta ^2J^2q_0^{p+2}2[p(p1)\beta ^2J^2q_0^p]^2𝒮_6}{[p(p1)\beta ^2J^2q_0^p]^2𝒮_6}.$$
(32)
This gives a perturbative expression for $`x`$ on the C1RSB line using only the RS value of $`q`$.<sup>1</sup><sup>1</sup>1Our (32) for $`x_\mathrm{c}`$ differs from that given by (36) of de Oliveira and Fontanari, due to their omission of one of these higher order corrections. Further details may be found in our comment on that paper.

The calculation of the 1RSB stability of the C1RSB transition line follows a similar method. Expanding (21) in powers of $`ϵ`$, everything up to third order vanishes by the self-consistent equations. Taking a linear combination of the inequality (21) and the equations (20) such that all terms to third order cancel, the fourth order term gives

$$\begin{array}{c}4(p2)(8p15)q+48p(p1)(p2)\beta ^2J^2q^p\hfill \\ \hfill 4p^2(p1)^2\beta ^4J^4q^{2(p1)}\left[4q+3(p2)(2+x_\mathrm{c})𝒮_6\right]\\ \hfill +3p^3(p1)^3\beta ^6J^6q^{3(p1)}\left[8(1+x_\mathrm{c})𝒮_6(6+8x_\mathrm{c}+x_\mathrm{c}^2)𝒮_8\right]>0.\end{array}$$
(33)
At high temperatures, this is satisfied and the continuous replica symmetry breaking transition is to 1RSB. Below a point (called the RSB triple point) found by solving (15), and (16) and (33) as equalities, the transition is to FRSB. Within the RSB phase there is a line marking a transition between 1RSB and FRSB which meets the RS AT line at this triple point.

## Appendix B An equation for the spinodal transition line

In this appendix we use the definition

$$t\mathrm{tanh}\left(\sqrt{\frac{1}{2}p\beta ^2J^2q^{p1}}z+\beta h\right).$$
(34)
The spinodal transition to ferromagnetism occurs at the minimum $`J_0`$ for which an $`M0`$ solution is possible. In the RS region, the ferromagnetic solution is given by (15), (17), and (10). Differentiating these three with respect to $`h`$, and using integration by parts,

$`{\displaystyle \frac{q}{h}}`$ $`={\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\left[(13t^2)(1t^2)\frac{1}{2}p(p1)\beta ^2J^2q^{p2}\frac{q}{h}+2t(1t^2)\beta \right]},`$ (35)
$`{\displaystyle \frac{M}{h}}`$ $`={\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\left[t(1t^2)\frac{1}{2}p(p1)\beta ^2J^2q^{p2}\frac{q}{h}+(1t^2)\beta \right]},`$ (36)
$`1`$ $`={\displaystyle \frac{J_0}{h}}M^{r1}+(r1)J_0M^{r2}{\displaystyle \frac{M}{h}}.`$ (37)
Asserting the spinodal condition $`J_0/h=0`$ and eliminating the unwanted derivatives gives

$$M=(r1)\beta h\left[1\frac{dz}{\sqrt{2\pi }}e^{z^2/2}t^2\frac{2\left[\frac{dz}{\sqrt{2\pi }}e^{z^2/2}t(1t^2)\right]^2}{\frac{2}{p(p1)\beta ^2J^2q^{p2}}\frac{dz}{\sqrt{2\pi }}e^{z^2/2}(1t^2)(13t^2)}\right].$$
(38)
We are often interested in the intersection between the spinodal line and the C1RSB line. At this point, (16) holds as an equality, and (38) simplifies to

$$M=(r1)\beta h\left[1\frac{dz}{\sqrt{2\pi }}e^{z^2/2}t^2\frac{\left[\frac{dz}{\sqrt{2\pi }}e^{z^2/2}t(1t^2)\right]^2}{\frac{dz}{\sqrt{2\pi }}e^{z^2/2}t^2(1t^2)}\right].$$
(39)
This result can be used as discussed in §4. The peak of the 1RSB transition in the ferromagnetic region is given by the solution of the RS self-consistency equation (15), the RS AT condition (16) as an equality, and the D1RSB equation $`x_\mathrm{c}=1`$ with $`x_\mathrm{c}`$ given by (32). When this peak hits the spinodal line, the D1RSB curve in the ferromagnetic region vanishes. This occurs at $`r=r_1`$, obtained for given $`p`$ by solving these three equations and (39). Similarly, the RSB triple point is given by (15), (16) as an equality, and the perturbative 1RSB AT condition (33) as an equality. When this point hits the spinodal line, the 1RSB ferromagnetic region vanishes. This occurs at $`r=r_2`$, obtained by solving these three equations and (39). Some values are shown in table 1.

## Appendix C Proofs for §5

### Appendix C.1 Proof from gauge transformation

In this subsection, we prove (25) using a gauge transformation technique. We explicitly write the configurational average appearing in (23):

$`P_m(y)`$ $`={\displaystyle \underset{i_1<\mathrm{}i_p}{}P_J(J_{i_1\mathrm{}i_p})dJ_{i_1\mathrm{}i_p}\frac{1}{Z_J(\beta )}\underset{\{\sigma _i\}}{}\delta \left(yN^1\underset{i}{}\sigma _i\right)\mathrm{exp}\left(\beta _p(\{\sigma _i\})\right)}`$ (40)
$`={\displaystyle 𝒟J\mathrm{exp}\left(\frac{N^{p1}}{p!J^2}\underset{i_1<\mathrm{}i_p}{}J_{i_1\mathrm{}i_p}^2\right)}`$
$`\times {\displaystyle \frac{1}{Z_J(\beta )}}{\displaystyle \underset{\{\sigma _i\}}{}}\delta \left(y\frac{1}{N}{\displaystyle \underset{i}{}}\sigma _i\right)\mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}\left(J_{i_1\mathrm{}i_p}{\displaystyle \frac{J_0(p1)!}{N^{p1}}}\right)\sigma _{i_1}\mathrm{}\sigma _{i_p}\right]`$
$`\times {\displaystyle \frac{1}{Z_J(\beta )}}{\displaystyle \underset{\{\tau _i\}}{}}\mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}\left(J_{i_1\mathrm{}i_p}{\displaystyle \frac{J_0(p1)!}{N^{p1}}}\right)\tau _{i_1}\mathrm{}\tau _{i_p}\right]`$ (41)
where we have introduced a factor of unity at the end, defined $`𝒟J_{i_1<\mathrm{}i_p}\sqrt{\frac{N^{p1}}{\pi pJ^2}}dJ_{i_1\mathrm{}i_p}`$, and $`Z_J(\beta )`$ is the partition function for a given choice of the couplings $`J_{i_1\mathrm{}i_p}`$:
$`Z_J(\beta )`$ $`={\displaystyle \underset{\{\sigma _i\}}{}}\mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}\left(J_{i_1\mathrm{}i_p}{\displaystyle \frac{J_0(p1)!}{N^{p1}}}\right)\sigma _{i_1}\mathrm{}\sigma _{i_p}\right].`$ (42)
Shifting the integration variable by $`J_0(p1)!/N^{p1}`$ and reordering trivially, we obtain
$`P_m(y)`$ $`={\displaystyle \underset{\{\tau _i\}}{}}{\displaystyle 𝒟J\mathrm{exp}\left(\frac{N^{p1}}{p!J^2}\underset{i_1<\mathrm{}i_p}{}J_{i_1\mathrm{}i_p}^2\underset{i_1<\mathrm{}i_p}{}\frac{2J_0}{pJ^2}J_{i_1\mathrm{}i_p}\underset{i_1<\mathrm{}i_p}{}\frac{J_0^2(p1)!}{pN^{p1}J^2}\right)}`$
$`\times {\displaystyle \frac{1}{Z_J(\beta )^2}}{\displaystyle \underset{\{\sigma _i\}}{}}\delta \left(y\frac{1}{N}{\displaystyle \underset{i}{}}\sigma _i\right)\mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}J_{i_1\mathrm{}i_p}\sigma _{i_1}\mathrm{}\sigma _{i_p}\right]`$
$`\times \mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}J_{i_1\mathrm{}i_p}\tau _{i_1}\mathrm{}\tau _{i_p}\right].`$ (43)
We now perform a gauge transformation with $`\sigma _i\sigma _i\tau _i`$ and $`J_{i_1\mathrm{}i_p}J_{i_1\mathrm{}i_p}\tau _{i_1}\mathrm{}\tau _{i_p}`$, and obtain
$`P_m(y)`$ $`={\displaystyle \underset{\{\tau _i\}}{}}{\displaystyle 𝒟J\mathrm{exp}\left(\frac{N^{p1}}{p!J^2}\underset{i_1<\mathrm{}i_p}{}J_{i_1\mathrm{}i_p}^2\underset{i_1<\mathrm{}i_p}{}\frac{2J_0}{pJ^2}J_{i_1\mathrm{}i_p}\tau _{i_1}\mathrm{}\tau _{i_p}\underset{i_1<\mathrm{}i_p}{}\frac{J_0^2(p1)!}{pN^{p1}J^2}\right)}`$
$`\times {\displaystyle \frac{1}{Z_J(\beta )^2}}{\displaystyle \underset{\{\sigma _i\}}{}}\delta \left(y\frac{1}{N}{\displaystyle \underset{i}{}}\sigma _i\tau _i\right)\mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}J_{i_1\mathrm{}i_p}\sigma _{i_1}\mathrm{}\sigma _{i_p}\right]`$
$`\times \mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}J_{i_1\mathrm{}i_p}\right].`$ (44)
If $`2J_0/pJ^2=\beta `$ then exchanging two terms gives
$`P_m(y)`$ $`={\displaystyle \underset{\{\tau _i\}}{}}{\displaystyle 𝒟J\mathrm{exp}\left(\frac{N^{p1}}{p!J^2}\underset{i_1<\mathrm{}i_p}{}J_{i_1\mathrm{}i_p}^2\underset{i_1<\mathrm{}i_p}{}\frac{2J_0}{pJ^2}J_{i_1\mathrm{}i_p}\underset{i_1<\mathrm{}i_p}{}\frac{J_0^2(p1)!}{pN^{p1}J^2}\right)}`$
$`\times {\displaystyle \frac{1}{Z_J(\beta )^2}}{\displaystyle \underset{\{\sigma _i\}}{}}\delta \left(y\frac{1}{N}{\displaystyle \underset{i}{}}\sigma _i\tau _i\right)\mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}J_{i_1\mathrm{}i_p}\sigma _{i_1}\mathrm{}\sigma _{i_p}\right]`$
$`\times \mathrm{exp}\left[\beta {\displaystyle \underset{i_1<\mathrm{}i_p}{}}J_{i_1\mathrm{}i_p}\tau _{i_1}\mathrm{}\tau _{i_p}\right].`$ (45)
Finally, shifting the integration variable back, we obtain

$$\begin{array}{c}P_m(y)=\hfill \\ \hfill \underset{i_1<\mathrm{}i_p}{}P_J(J_{i_1\mathrm{}i_p})dJ_{i_1\mathrm{}i_p}\frac{1}{Z_J(\beta )^2}\underset{\{\sigma _i\}}{}\underset{\{\tau _i\}}{}\delta \left(y\frac{1}{N}\underset{i}{}\sigma _i\tau _i\right)\mathrm{exp}\left[\beta _p(\{\sigma _i\})\beta _p(\{\tau _i\})\right]\\ \hfill =P_q(y)\end{array}$$
(46)
which completes the proof.

### Appendix C.2 Other proofs concerning the Nishimori line

In this subsection, we prove the results concerning the Nishimori line and two triple points mentioned in §5. (The demonstration of several mathematical identities used here is left to Appendix C.3.) We first notice that (25) means not only that the ferromagnet on the Nishimori line ($`p\beta J^2=2J_0`$) is replica symmetric, but that it has $`q=M`$. It follows that $`\frac{1}{2}p\beta ^2J^2q^{p1}=\beta J_0M^{r1}`$, and with (10) this is $`\beta h`$. We may verify $`q=M`$ independently of the gauge result, as using (15) and (17) this becomes

$$\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2\left(\sqrt{\beta h}z+\beta h\right)=\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}\left(\sqrt{\beta h}z+\beta h\right),$$
(47)
an identity proved in Appendix C.3.

The spin glass transition line satisfies (20) with $`x=1`$ and $`q_0=h=0`$ (as $`M=0`$). Matters are simplified by the independence of $`G`$ from $`z_0`$, which means that the integrals over $`z_0`$ become trivial. The assertion that the $`q_1`$ here is the same as the $`q`$ of the RS ferromagnet on the Nishimori line becomes, with (20b),

$$\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2\left(\sqrt{\beta h}z+\beta h\right)=\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)\mathrm{tanh}^2\left(\sqrt{\beta h}z_1\right)}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)},$$
(48)
which is also proved in Appendix C.3. The other non-trivial self-consistent equation, (20c), becomes

$$\begin{array}{c}\frac{1}{4}(p1)\beta ^2J^2q^p=\mathrm{ln}\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)\hfill \\ \hfill +\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)\mathrm{ln}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)}.\end{array}$$
(49)
The thermodynamic transition line is located by equating the free energy of the ferromagnetic solution and that of the paramagnetic or spin glass solution. The free energy for $`r=p`$ is related to that in a field via (6) and (8). Where the thermodynamic and spin glass transition lines cross, the free energy of the spin glass is $`f_{\mathrm{SG}}=f_{\mathrm{field}}^{(1\mathrm{R}\mathrm{S}\mathrm{B})}`$ with $`x=1`$, $`q_0=h=0`$, and $`f_{\mathrm{field}}^{(1\mathrm{R}\mathrm{S}\mathrm{B})}`$ given by (18). If the Nishimori line passes through this point then the previous results apply. The free energy of the ferromagnet is therefore

$$f_{\mathrm{FM}}=f_{\mathrm{field}}^{(\mathrm{RS})}+hM\frac{1}{r}J_0M^r=f_{\mathrm{field}}^{(\mathrm{RS})}+\frac{1}{2}(p1)\beta J^2q^p$$
(50)
with $`f_{\mathrm{field}}^{(\mathrm{RS})}`$ given by (14). The thermodynamic transition condition $`f_{\mathrm{FM}}=f_{\mathrm{SG}}`$ becomes

$$\frac{1}{4}(p1)\beta ^2J^2q^p\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{ln}\mathrm{cosh}\left(\sqrt{\beta h}z+\beta h\right)=\mathrm{ln}\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)$$
(51)
or, using (49),
$$\frac{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)\mathrm{ln}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)}{\frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}\mathrm{cosh}\left(\sqrt{\beta h}z_1\right)}\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{ln}\mathrm{cosh}\left(\sqrt{\beta h}z+\beta h\right)=0,$$
(52)
which is also proved in Appendix C.3. Thus, the Nishimori line does pass through the thermodynamic triple point where the paramagnet, spin glass, and RS ferromagnets meet.

The equations for the intersection of spinodal line and the C1RSB transition line are discussed in Appendix B. If this point lies on the Nishimori line then the above results apply. Setting $`p=r`$, $`M=q`$, and $`\frac{1}{2}p\beta ^2J^2q^{p1}=\beta h`$ in (39) and in (16) as an equality, dividing the former by the latter, and dividing the result by $`\frac{1}{2}p(p1)\beta ^2J^2q^{p1}`$, gives a condition for the consistency of this hypothesis which may be written

$$\begin{array}{c}\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sech}^4\left(\sqrt{\beta h}z+\beta h\right)1+\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2\left(\sqrt{\beta h}z+\beta h\right)\hfill \\ \hfill +\frac{\left[\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}\left(\sqrt{\beta h}z+\beta h\right)\mathrm{sech}^2\left(\sqrt{\beta h}z+\beta h\right)\right]^2}{\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2\left(\sqrt{\beta h}z+\beta h\right)\mathrm{sech}^2\left(\sqrt{\beta h}z+\beta h\right)}=0,\end{array}$$
(53)
which is proved in Appendix C.3. Thus, the Nishimori line also passes through the spinodal triple point where the paramagnet, RS ferromagnet, and RSB ferromagnet meet.

### Appendix C.3 Proofs of some mathematical identities

In this section we prove several identities used in Appendix C.2. These proofs will all rely on the lemmas

$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)f(z)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}f_\mathrm{E}(z+a)}`$ (54)
$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sinh}(az)f(z)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}f_\mathrm{O}(z+a)}`$ (55)
where $`f_\mathrm{E}`$ and $`f_\mathrm{O}`$ are the even and odd parts of the function $`f`$. The proof of these runs as follows:

$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\left\{\genfrac{}{}{0pt}{}{\mathrm{cosh}(az)}{\mathrm{sinh}(az)}\right\}f(z)}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\left(e^{az}\pm e^{az}\right)f(z)}`$ (56)
$`={\displaystyle \frac{1}{2}}e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}\left(e^{(za)^2/2}\pm e^{(z+a)^2/2}\right)f(z)}`$ (57)
$`={\displaystyle \frac{1}{2}}e^{a^2/2}\left({\displaystyle \frac{dz_1}{\sqrt{2\pi }}e^{z_1^2/2}f(z_1+a)}\pm {\displaystyle \frac{dz_2}{\sqrt{2\pi }}e^{z_2^2/2}f(z_2a)}\right)`$ (58)
$`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\left\{\genfrac{}{}{0pt}{}{f_\mathrm{E}(z+a)}{f_\mathrm{O}(z+a)}\right\}}.`$ (59)
The proof of (47) now runs as follows:

$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2(az+a^2)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{tanh}^2(az)}`$ using (54) (60)
$`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sinh}(az)\mathrm{tanh}(az)}`$ (61)
$`={\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}(az+a^2)}`$ using (55). (62)
The proof of (48) runs as follows:

$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}}`$ using (54) (63)
$`=e^{a^2/2};`$ (64)
$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{tanh}^2(az)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2(az+a^2)};`$ using (54) (65)
$`{\displaystyle \frac{\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{tanh}^2(az)}{\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)}}`$ $`={\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2(az+a^2)},`$ (66)
and the proof of (52) is identical with $`\mathrm{tanh}^2`$ replaced by $`\mathrm{ln}\mathrm{cosh}`$.

The proof of (53) runs as follows: using (54),

$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sech}^4(az+a^2)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{sech}^4(az)},`$ (67)
$`1`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)},`$ (68)
$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2(az+a^2)}`$ $`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{tanh}^2(az)};`$ (69)
using (55),
$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}}`$ $`\mathrm{tanh}(az+a^2)\mathrm{sech}^2(az+a^2)`$
$`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sinh}(az)\mathrm{tanh}(az)\mathrm{sech}^2(az)}`$ (70)
$`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{tanh}^2(az)\mathrm{sech}^2(az)};`$ (71)
and using (54) again
$`{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}}`$ $`\mathrm{tanh}^2(az+a^2)\mathrm{sech}^2(az+a^2)`$
$`=e^{a^2/2}{\displaystyle \frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\mathrm{tanh}^2(az)\mathrm{sech}^2(az)};`$ (72)
thus

$$\begin{array}{c}\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{sech}^4(az+a^2)1+\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2(az+a^2)\hfill \\ \hfill +\frac{\left[\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}(az+a^2)\mathrm{sech}^2(az+a^2)\right]^2}{\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{tanh}^2(az+a^2)\mathrm{sech}^2(az+a^2)}\\ \hfill =e^{a^2/2}\frac{dz}{\sqrt{2\pi }}e^{z^2/2}\mathrm{cosh}(az)\left[\mathrm{sech}^4(az)1+\mathrm{tanh}^2(az)+\mathrm{tanh}^2(az)\mathrm{sech}^2(az)\right];\end{array}$$
(73)
and this is identically zero since the term in square brackets vanishes using $`\mathrm{sech}^2(az)=1\mathrm{tanh}^2(az)`$.

## Acknowledgments

PG and DS would like to thank EPSRC (UK) for financial support, PG for research studentship 97304251 and DS for research grant GR/M04426. HN acknowledges the support of the Sumitomo Foundation. This work was also supported by the Anglo–Japanese Collaboration Programme between The Royal Society and the Japan Society for the Promotion of Science.


text:
# Lower Critical Field at Odds with A S-Wave Superconductivity in The New Superconductor 𝑀𝑔𝐵₂

\[

## Abstract

The lower critical field $`H_{c1}`$ has been carefully measured on a well shaped cylindrical sample of the new superconductor $`MgB_2`$ fabricated by high pressure synthesis. The penetration depth $`\lambda `$ is calculated from the $`H_{c1}`$ data. It is found that a linear relation of $`H_{c1}(T)`$ appears in whole temperature region below $`T_c`$. Furthermore a finite slope of $`dH_{c1}/dT`$ and $`d\lambda (T)/dT`$ remains down to the lowest temperature ( 2 K ). These are inconsistent with the expectation for a widely thought s-wave superconductivity in $`MgB_2`$.

\] Recently discovered new superconductor $`MgB_2`$ generates a new round of chasing in the field of superconductivity. Several thermodynamic parameters have already been derived, such as the upper critical field $`H_{c2}`$ = 15 - 20.4 T, the Ginzburg-Landau parameter $`\kappa `$ 26, and the superconducting critical current density $`j_c1.05\times 10^5A/cm^2`$ at 15 K and 1 T. Recently, Hall effect measurement has been done showing that the charge carriers are holes with very high density $`n_s=1.5\times 10^{23}cm^3`$. Measurement on Seebeck coefficient signals also positive charges in $`MgB_2`$. Theoretically there are two major diverse pictures accounting for the superconductivity. Based on band calculations, Kortus et al. proposed that the superconductivity is resulted from the strong electron-phonon interaction and the high phonon frequency associated with the light B element. This picture is supported by the recent isotope measurement on $`{}_{}{}^{10}B`$ and $`{}_{}{}^{11}B`$ yielding an isotope exponent $`\alpha `$ 0.26. In the second picture, Hirsch suggested an alternative explanation which conjectures the same mechanism of superconductivity in $`MgB_2`$ as in high temperature cuprate superconductors ( HTS ). Recently high pressure measurement shows that the superconducting transition temperature decreases linearly at a rate of $`dT_c/dP`$ = -1.6 K/GPa, being in contrast to that predicted by Hirsch’s theory. In order to know the superconducting mechanism, to probe the nature of superconducting gap and symmetry seems to be very important. Recent data from specific heat, tunneling spectroscopyand NMR study commonly suggest a phonon-mediated BCS picture with the s-wave superconducting energy gap $`\mathrm{\Delta }`$ ranging from 2 to 8 meV corresponding to the coupling strength from weak to strong. For a s-wave superconductor, it needs a minimum energy to excite the charges from the superconducting condensate to the excited states. Thus one can expect an exponentially activated temperature dependence in low temperature region for many thermodynamic parameters, such as the lower critical field $`H_{c1}`$ and penetration depth $`\lambda `$, etc.. Therefore one should observe a vanishing slope $`dH_{c1}(T)/dT`$ = 0 and $`d\lambda (T)/dT`$ = 0 in low temperature region. In this Letter, we present a precise determination of the lower critical field $`H_{c1}`$ in whole temperature region. Our data are clearly at odds with the prediction for a s-wave superconductivity in $`MgB_2`$.

Samples of $`MgB_2`$ investigated here were fabricated by high pressure synthesis ( P = 6 GPa at 950 C for 0.5 hour) described in detail elsewhere. The X-ray diffraction ( XRD ) analysis shows that it is nearly in a single phase with the second phase less than 1 wt.%. Our samples are very dense and hard with a metallic shiny surface after polishing. They look also like the hard $`YBa_2Cu_3O_7`$ bulks grown by melt-textured-growth technique, showing probably the trivial influence of possible weak links between grains ( if any ). Fig.1 shows that the onset $`T_c^{onset}`$ is 40 K and the transition width is less than 1 K. The inset of Fig. 1 shows a very sharp diamagnetic transition and a perfect diamagnetic signal, in consistent with the resistive result. To minimize the demagnetization factor, the sample was carefully cut and ground to a cylinder with a diameter of 1.58 mm and

length 3.50 mm. The DC magnetic measurements were carried out by a vibrating sample magnetometer ( VSM, Oxford 3001 ) and a superconducting quantum interference device ( SQUID, Quantum Design MPMS 5.5 T ). During the magnetic measurements, the sample axis is always aligned along the field direction so that the demagnetization effect can be neglected. After zero-field-cooled from 50 K to a desired temperature, the magnetic field is applied slowly up to 900 Oe, which is much higher than Hc1. Fig. 2(a) and 2(b) show the DC magnetization curves obtained by VSM and SQUID respectively. All curves show clearly a common linear dependence of the magnetization on field caused by Meissner effect at low fields.

The value of Hc1 is determined by examining the point of departure from linearity on the initial slope of the magnetization curve. To do this, a linear fit with points between 10 Oe and 100 Oe at 2 K is made. This line passes also through the data points measured at other temperatures in low field region. The fitted linear line is called

the ” Meissner line ” ( ML ), as the dotted line shown in Fig. 2, which represents the magnetization curve of Meissner state. The results of subtracting this ML from magnetization curves of SQUID measurements are plotted in Fig. 3 in a logarithmic scale. It is important to note that one should avoid to use a linear scale plot of $`\mathrm{\Delta }`$M vs. H to determine $`H_{c1}`$ since that will induce great uncertainty. The data from VSM bear close resemblance, thus they are not shown here. The value of Hc1 is easily obtained by choosing a proper criterion of $`\mathrm{\Delta }`$M. In Fig. 4(a), we plot the Hc1 acquired by using a criterion of $`\mathrm{\Delta }M=1\times 10^3emu`$. The Hc1 measured by VSM and SQUID coincide in the same line, showing that our results are independent on measuring device and technique. Extrapolating the Hc1 curve to zero temperature, we get Hc1(0) $``$ 295 Oe, being close to 320 Oe found by Takano et al. A striking feature here is that the curve $`H_{c1}(T)`$ shows a good linear behavior in whole temperature region as shown by the straight line in Fig. 4(a).

To obtain the penetration depth $`\lambda `$, we use the expression

$$H_{c1}=\frac{\mathrm{\Phi }_0}{4\pi \lambda ^2}ln\kappa ,$$
(1)
where $`\mathrm{\Phi }_0=hc/2e2\times 10^7Gcm^2`$ is the flux quantum, and $`\kappa `$ 26 is the Ginzburg-Landau parameter. Using Hc1(0) = 295 Oe obtained above, we find $`\lambda `$(0) $``$ 1325 $`\AA `$, being close to $`\lambda `$(0) = 1400 $`\AA `$ found by Finnemore et al.. The deviation of the penetration depth from its zero temperature value, $`\mathrm{\Delta }\lambda =\lambda (T)\lambda (0)`$, is shown in Fig. 4(b). It is remarkable that the dependence of $`\mathrm{\Delta }\lambda `$ on temperature below 20K ($`T_c/2`$) is almost linear, i.e., $`\mathrm{\Delta }\lambda T`$, which is more explict in the inset. A rough estimate gives the slope $`d\lambda /dT28\AA /K`$.

In order to know whether this linear behavior is induced by the uncertainty in using the criterion for $`\mathrm{\Delta }`$M, we have used different criterions, such as $`1\times 10^4emu`$, for the determination of $`H_{c1}`$, but the linear dependence is not changed. This linear dependence conflicts the prediction for a s-wave superconductivity that supported by some expriments mentioned in the first paragraph of this paper.

For a S-wave superconductor, the finite energy gap manifests itself with an exponentially activated temperature dependence of many thermodynamic properties. In the s-wave BCS theory, $`\mathrm{\Delta }\lambda `$ is given by,

$$\frac{\mathrm{\Delta }\lambda (\mathrm{T})}{\lambda (0)}3.3(\frac{\mathrm{T}_\mathrm{c}}{\mathrm{T}})^{\frac{1}{2}}exp(\frac{\mathrm{\Delta }}{\mathrm{k}_\mathrm{B}\mathrm{T}_\mathrm{c}}\frac{\mathrm{T}_\mathrm{c}}{\mathrm{T}}),$$
(2)
where $`\mathrm{\Delta }`$ is the energy gap. The exponential term causes

the very weak temperature dependence of $`\mathrm{\Delta }\lambda (T)`$, i.e., $`d\lambda (T)/dT0`$, as shown by the dotted line in the inset of Fig. 4(b). It is easy to see that the linear dependence of $`\mathrm{\Delta }\lambda (T)`$ versus T shown in Fig. 4(b) and its inset is clearly different from the prediction by the s-wave BCS theory. Even for a strong coupling superconductor with probably an anisotropic energy gap, it is still hard to understand the finite slope of $`H_{c1}(T)`$ and $`\lambda (T)`$ in low temperature region. There was an argument contributing this linearity to the classical superconducting phase fluctuation. This possibility can be however ruled out by the observation of very weak fluctuation effect in $`MgB_2`$ characterized by the non-rounded onset superconducting transitions measured at different fields. Another argument would be that our samples contain still grains and thus weak links. This is basically correct, but as mentioned already, our samples are prepared by high pressure synthesis with very high density and hardness together with rather high critical current density, therefore the couplings between grains ( if any ) are very strong. Furthermore, the zero temperature penetration depth $`\lambda (0)`$ determined from our high pressure synthesized cylindrical sample is very close to that found by Finnemore et al. based on a consideration on the intrinsic parameters of $`MgB_2`$. This may manifest that the linear $`H_{c1}(T)`$ determined in our present work is an intrinsic property.

This linear behavior is similar to that in HTS. In cuprate superconductors, it has been ascribed to the d-wave symmetry of energy gap based on the consideration as follows. In a pure d-wave superconductor, the energy gap along the node directions ( $`k_x=\pm k_y`$ ) vanishes and the spectrum $`N_s(E)`$ of quasiparticle excitations in the superconducting phase is gapless, which results in the linear dependence of $`N_s(E)`$ on E at low excitation energies. A finite temperature will excite certain quasiparticles leading to a linear dropping of the superfluid density $`\rho _s`$ with temperature. As to other p-wave superconductors, $`\mathrm{\Delta }\lambda `$ also shows non-exponential dependence on T because of the nodes of order parameter. The linear $`H_{c1}(T)`$ measured in $`MgB_2`$ may not imply a d-wave or p-wave superconductivity, but calls certainly for a new understanding. For example, according to the theory by Hirsch et al., it appears to be the same ” hole superconductivity ” in HTS and $`MgB_2`$. Further investigation may lead to an answer to the question whether the finite slope of $`H_{c2}(T)`$ and $`\mathrm{\Delta }\lambda (T)`$ at low temperatures for HTS and $`MgB_2`$ are due to the same cause. In addition, the penetration depth measurement by other techniques, such as $`\mu SR`$, microwave and a.c. susceptibility is very desirable in order to have a comparison to our data.

In conclusion, we have measured the M-H curve of a cylindrical $`MgB_2`$ sample and obtained its lower critical field $`H_{c1}(T)`$. A striking linear curve $`H_{c1}(T)`$ has been observed in the whole temperature region. Further calculation of the penetration depth shows also a finite slope in the low temperature limit. All these are at odds with the expectations for a s-wave superconductivity. The reason is unknown but warrants certainly further investigation in order to get a final consensus of the order parameter symmetry of $`MgB_2`$.

###### Acknowledgements.

This work is supported by the National Science Foundation of China (NSFC 19825111) and the Ministry of Science and Technology of China ( project: NKBRSF-G1999064602 ).