format: text

text:
# A Submillimeter Study of the Star-Forming Region NGC 7129

## 1 Introduction

NGC 7129 is a reflection nebula seen against a molecular cloud and estimated by Shevchenko et al. (Shev89 (1989)) to be at a distance of 1.25 kpc. There are many signs in this region which support the idea of triggered and possibly ongoing star formation: the association of a very young cluster formed by the B3-type stars BD +65 1637 and BD +65 1638, and the B5/7-type star LkH$`\alpha `$ 234, several embedded infrared sources (Harvey et al. Harvey84 (1984); Weintraub et al. Wein94 (1994); Cabrit et al. Cabrit97 (1997)) as well as reflection nebulae and Herbig-Haro (HH) objects (Hartigan & Lada Hart85 (1985); Miranda et al. Miranda94 (1994)). There are also several molecular outflows in the region: one associated with an optical jet near LkH$`\alpha `$ 234 (Ray et al., Ray90 (1990); Edwards & Snell Edwards83 (1983)), another bipolar outflow near the far-infrared source NGC7129 FIRS 2 and one may be driven by the T Tauri star V350 Cep (Hartigan & Lada Hart85 (1985); Goodrich Goodrich86 (1986); Miranda et al. Miranda94 (1994)).

A molecular cavity is revealed in both molecular and infrared maps (Bechis et al. Bechis78 (1978); Bertout Bertout87 (1987)), but its origin is still unknown. The location of the two bright stars, BD +65 1637 and BD +65 1638, inside the cavity, suggests three mechanisms by which the molecular cavity could have been created: (i) through the advancement of a photon-dominated region (PDR) around those stars, in which the molecular gas was dissociated by their intense UV radiation, (ii) by an expanding shell associated with stellar mass loss, or (iii) by radiation pressure acting on grains, together with gas-grain drag (Bechis et al. Bechis78 (1978)). The two mid-B type stars in the cavity, BD +65 1637 and BD +65 1638, are also the oldest in the cluster. Therefore stellar mass loss may have been very effective at earlier times and the winds from those two stars may have triggered the formation of other bright stars in NGC 7129, including LkH$`\alpha `$ 234 and SVS 13. As for what mechanism is dominating at present time, the question remains open. PDR models can marginally explain the line fluxes of \[OI\] and \[CII\] emission in the region, although the C-type shock models cannot be ruled out (Lorenzetti et al. Loren99 (1999); Giannini et al. Giannini99 (1999)).

This paper employs new, high resolution submillimeter observations of the molecular ridge in dust continuum emission, 13CO line emission, and C18O line emission. One goal of this study is to improve the present picture of the structure of the ridge and the ridge-cavity boundary. Another goal is to obtain the physical properties of the compact far-infrared sources in the ridge. A third aim is to obtain the gas/dust ratio, thus determining whether CO isotopomers are a good surrogate for molecular hydrogen, or whether CO has substantially frozen out onto dust grains. The final goal is to search for variations in the dust opacity with position through the field. The mapped region includes the dense molecular cloud, a photodissociation region, and a powerful outflow. We might expect processes such as grain coagulation and mantle evaporation to have a strong spatial dependency here and for this to be reflected in spatial variations of the dust emissivity index, $`\beta `$.

## 2 Observations and data reduction

### 2.1 $`{}_{}{}^{13}CO`$ and C18O observations

Observations of 13CO and C18O J=3$``$2 were made with the James Clerk Maxwell Telescope, on Mauna Kea, Hawaii, in August 1994 and November 1997 respectively, with the receiver B3i and the Digital Autocorrelation Spectrometer (DAS). The observations were made by position switching, with the reference position 10 East of LkH$`\alpha `$ 234, which was determined to be free of significant CO emission. The main beam efficiency ($`\eta _{bm}`$) was 0.58 and the half-power beam (HPBW) was $``$ 14′′ for both 13CO and C18O observations.

A grid map with the spacing of 7$`^{\prime \prime }.`$5 x 7$`^{\prime \prime }.`$5 was taken in the 13CO J=3$``$2 (330.58 GHz) rotational transition. The map was centered near the star LkH$`\alpha `$ 234111The astrometric position for LkH$`\alpha `$ 234 given by Clements & Argyle (Clements84 (1984)) is $`\alpha `$(2000) = 21h 43m 06s.802 and $`\delta (2000)`$ = 66 06 54$`^{\prime \prime }.`$5., at coordinates $`\alpha (2000)`$ = 21h 43m 06s.08 and $`\delta (2000)`$ = 66 06 56$`^{\prime \prime }.`$29, and extending from $``$100′′ to +100′′ in $`\alpha `$ and from $``$150′′ to 75′′ in $`\delta `$ offsets. The instrumental velocity resolution was 0.28 km s-1. C18O J=3$``$2 (329.33 GHz) spectra were taken at nine positions around the star LkH$`\alpha `$ 234 and near the molecular ridge. Figure 1 shows the 13CO integrated intensity ($`T_B𝑑v`$) contour map superimposed on a K image of the NGC 7129 region from Hodapp (Hodapp94 (1994)). The nine C18O offsets with respect to the star LkH$`\alpha `$ 234 are shown with cross symbols.

### 2.2 Submillimeter continuum observations

Two deep 850 $`\mu `$m and 450 $`\mu `$m–maps of LkH$`\alpha `$ 234 and the bright-rimmed cloud ridge southwest of LkH$`\alpha `$ 234 have been obtained using the Submillimeter Common User Bolometer Array, SCUBA (Holland et al. Holland99 (1999)), on the 15 m James Clerk Maxwell Telescope, on Mauna Kea, Hawaii, on October 15, 1997. The weather conditions were dry and stable with an 850 $`\mu `$m zenith optical depth of $``$ 0.19. The two overlapping fields were obtained in jiggle-map mode with a chop-throw of 150′′ in Right Ascension with a total integration time of 85 minutes/field. For a description of SCUBA and its observing modes, see Holland et al. (Holland99 (1999)). Unfortunately the maps were taken without usable calibration and pointing observations and we have therefore reduced the data by determining the HPBW of the telescope from observations of Uranus from other nights in September and October with similar sky conditions. The same nights were also used to determine gain conversion factors for the 850 $`\mu `$m and 450 $`\mu `$m filters. The HPBW was found to be $``$ 15$`^{\prime \prime }.`$6 $`\times `$ 13$`^{\prime \prime }.`$5 at 850 $`\mu `$m and 9$`^{\prime \prime }.`$1 $`\times `$ 7$`^{\prime \prime }.`$8 at 450 $`\mu `$m with the beam broadened in the chop direction. The submillimeter position of LkH$`\alpha `$ 234 (Table 3), which has an accuracy of $``$ 1′′, was taken from three short-integration maps interlaced with pointing observations in SCUBA commissioning time during Spring 1997. These maps have not been added in the final data set, because they were taken with a 120′′ Azimuth chop and there is clear evidence that the array was chopping onto emission. However, if we compare the total flux density of LkH$`\alpha `$ 234 (i.e. background subtracted flux density) it agrees very well with data taken in October.

We also include a map of the H2O maser and bipolar outflow source NGC 7129 FIRS2, which was used as a test source during SCUBA commissioning time (Sandell Sandell97 (1997)). A composite map at 850 $`\mu `$m, including all the three overlapping fields is shown in Figure 2.

The basic data reduction was done with SURF (Jenness & Lightfoot Jenness99 (1999)). Even with a 150′′-chop it is more than likely that we may have chopped onto some emission in both fields, but it does not appear to affect the morphology of the strong dust ridge in which LkH$`\alpha `$ 234 is embedded. Some negative emission is detected west of LkH$`\alpha `$ 234, which confirms that there is some emission in the off-position. It is therefore likely that emission east of the ridge is also affected, which is why the 850 and 450$`\mu `$m maps appear somewhat different. To map a cloud ridge in beam switch mode is not very advisable, but at the time SCUBA was only available in jiggle map mode. Although chopping into the cloud cannot be avoided, as long as the off-source emission is faint and relatively uniform, the morphology of the cloud ridge should not be affected. However, since the beam size and error lobe contribution differ between 850$`\mu `$m and 450$`\mu `$m, this leads to systematic errors in the integrated intensities (Section 3.2).

The peak fluxes are uncertain to about 10$`\%`$ for 850$`\mu `$m and 20$`\%`$ for 450$`\mu `$m. The final pointing corrected and calibrated maps have an rms noise of $``$ 15 mJy/beam and $``$ 120 mJy/beam for 850 and 450$`\mu `$m, respectively. These maps were written out as FITS files and read into the MIRIAD package (Sault, Teuben, & Wright Sault95 (1995)) for further analysis.

## 3 Analysis

### 3.1 13CO and C18O Emission

A sharply bounded molecular cavity is seen in 13CO J=3$``$2 emission westward of LkH$`\alpha `$ 234, most probably excavated by the B3-type star BD +65 1637. The lack of molecular emission and the ridge morphology were also noted in the CO J=1$``$0 and NH3 observations (Bechis et al. Bechis78 (1978); Fuente et al. Fuente98 (1998); Güsten & Marcaide Gusten86 (1986)). However, our high resolution (half-beam sampling with a 14′′ beamsize) observations resolve the ridge structure in much more detail.

There are four 13CO peaks along the ridge, visible in Figure 1 as contour peaks. These 13CO structures are assumed to be physical clumps and named NGC7129 CO1, NGC 7129 CO2, NGC 7129 CO3, and NGC 7129 CO4. The positions of these four 13CO peaks, together with their peak antenna temperatures, line widths, and line center velocities for the 13CO emission are given in Table 1. NGC7129 CO1 coincides within errors with LkH$`\alpha `$ 234 (position 5 in Table 2), NGC 7129 CO4 with position 7 and the other two clumps, NGC7129 CO2 and NGC7129 CO3 are in located close to the positions 6 and 3, respectively.

Perhaps the most striking feature of the 13CO map is the sharpness of the transition between the molecular cavity to the west and the molecular cloud to the east. 13CO line strengths drop by factors of 20 or more across the boundary. The narrowness of the transition region between atomic and molecular gas has implications for the three dimensional topology of the region. There seem to be three possibilities: Firstly, the boundary region may be a sheet which happens to be perpendicular to the line-of-sight. Secondly, the molecular cavity may be a bowl, which is surrounded by molecular gas. Thirdly, the molecular cloud may take the form of a rather narrow ridge, which is partially surrounded by atomic gas. The first possibility is improbable, while the second should lead to a broader transition region (molecular to atomic). We favor the third picture, in which the dense portion of the molecular cloud consists of a clumpy ridge with a line-of sight extent comparable to its extent in the plane of the sky (roughly 0.5 pc).

Another remarkable feature which can be seen in Figure 1 is the narrow filament of emission seen in the K image (Hodapp 1994). It extends roughly south from LkH$`\alpha `$ 234 and then curves to the west, keeping a nearly constant surface brightness. In an unpublished observations using an IR camera with a circularly variable filter on the CFHT, Mitchell and Nadeau found this filament to be also a source of strong emission in the S(1) v= 1$``$0 line of H2. The optical filament follows closely the outermost 13CO contours, which separate the cavity from the molecular cloud. Can this close vicinity be explained by standard PDR models? The presence of the photodissociation region is a natural explanation for both the sharp molecular boundary and the line-emitting filament. The H2 S(1) v= 1$``$0 line emission is then due to fluorescence. The PDR is very likely produced by UV radiation from the early B star BD+65 1638. However, the detection of all three types of emission: optical, CO and H2 v= 1$``$0, at the same location, is in contrast with predictions of steady state stationary PDR models (Hollenbach & Tielens Hollenbach97 (1997)). These models predict a more gradual succession of emission layers as one moves away from the UV ionization front, and where the H2 v= 1$``$0 transition layer would be separated by about 10′′ from the molecular CO emission. On the other hand, the observed morphology may be explained in the framework of nonstationary PDR models (Bertoldi & Draine Bertoldi96 (1996)), which show that the inhomogeneities in the molecular cloud surface can lead to the merging of the ionization and the dissociation front.

Because 13CO is likely to have an appreciable optical depth, C18O 3$``$2 spectra were obtained at several positions. The positions for which we have C18O spectra are marked by the cross symbols in Figure 1. A comparison between the 13CO and C18O spectra at the nine offset positions is shown in Figure 3. Positions are given as offsets in arcseconds from LkH$`\alpha `$ 234. The line intensities are given in antenna temperatures, $`T_{A}^{}{}_{}{}^{}`$, after corrections for losses due to the Earth atmosphere and the telescope. For all spectra only the linear baselines were removed, using the SPECX software package (Padman Padman92 (1992)).

Near LkH$`\alpha `$ 234, the lines are broadened due to contributions from the outflow. Position 1 (Table 2) is near the peak of CO emission in the redshifted outflow lobe. Both the 13CO and C18O spectra are double peaked at this position. The value of the line ratio shows that C18O is optically thin. The two peaks therefore represent two real kinematic components rather than a single line suffering self absorption. The peak at $``$10 km s-1 is the ambient gas at this position, while the peak at $``$7 km s-1 is the postshock gas.

Several line parameters have been calculated at each offset, assuming local thermodynamical equilibrium (LTE) conditions. A value of 10 was used for the abundance ratio $`N_{{}_{}{}^{13}CO}`$/$`N_{C{}_{}{}^{18}O}`$, estimated from the partial ratios $`N_{C{}_{}{}^{18}O}/N_{H_2}`$ 2 $`\times `$ 10-7 (Frerking, Langer $`\&`$ Wilson Frerking82 (1982)) and the standard value $`N_{{}_{}{}^{12}CO}/N_{{}_{}{}^{13}CO}`$ 56 (Wilson and Rood, Wilson94 (1994)). The intensity ratio, 13CO/C18O, was used to find the optical depth, assuming a common Tex for the two transitions. The value of Tex was computed from the LTE expression relating line strength (given by the radiation temperature), line optical depth $`\tau `$, the excitation and radiation temperatures:

$$T_{ex}=\frac{h\nu }{k}\frac{1}{\mathrm{ln}[1+(1e^\tau )h\nu /kT_R]}.$$
(1)
In the above relation, the radiation temperature, $`T_R`$, is related to the antenna temperature, $`T_{}^{}{}_{A}{}^{}`$, by:

$$T_R=\frac{T_{}^{}{}_{A}{}^{}}{\eta _{bm}f},$$
(2)
where the beam efficiency $`\eta _{bm}`$ was 0.58 and the filling factor f was assumed to be 1. We will further assume that the gas kinetic temperature is equal to the excitation temperature, but caution that a density greater than 104 cm-3 is required for this to be true. The results of the LTE analysis are shown in Table 2. For LkH$`\alpha `$ 234, the excitation temperature is found to be about 30 K, in agreement with previous values of 20 $``$ 30 K obtained from lower resolution (70′′) CO observations (Bechis et al. Bechis78 (1978)) and with the 20 $``$ 25 K obtained from NH3 (40′′) observations (Güsten & Marcaide Gusten86 (1986)).

At position 6 (Table 2), the gas has an excitation temperature of $``$ 34 K, which is higher than the typical value for a dark core, but lower than the previous value of 63 K obtained by Mitchell and Matthews (Mitch94 (1994)) from 12CO and 13CO data. These two temperatures can be reconciled, taking into account that the isotopomers from which they have been derived (12CO for 63K and 13CO for 34 K) have different optical depths and hence probe different regions of the gas. With its lower optical depth, 13CO will probe deeper into the PDR, where the temperature is lower.

The three last entries in Table 2 represent positions increasing in distance from the cloud boundary. The excitation temperature is found to decrease with distance from the boundary, from 38 K through 31 K to 21 K. This decrease is expected as the stellar radiation is increasingly attenuated with distance into the molecular cloud. At position 1 (Table 2), which is farthest from the cloud boundary, the excitation temperature of the ambient gas is lowest, $``$ 13 K.

Since the clumps have approximately the same size as our beam and since the beam does not always include a clump of gas, it is reasonable to calculate the masses of the gas within the beam (14′′). The results, including the fractional He abundance and assuming a distance of 1.25 kpc to the source, are shown in Table 2. The NGC 7129 CO1 clump which coincides with LkH$`\alpha `$ 234, contains 10.7 M. The masses within the beam, at the two offsets adjacent to the clump matching the pre-stellar source LkH$`\alpha `$ 234 SMM 2 (positions 3 and 6 in Table 2) are 12.8 and 7.2 M, respectively. However, those values should be regarded as upper limits since the temperatures are probably lower, as suggested by the more optically thin NH3 data.

These results also depend on the accuracy in the distance determination. The distance estimate of 1kpc to NGC7129, reported by Racine (Racine68 (1968)), has been most commonly used. In our analysis we have adopted the value of 1.25 kpc, determined by Shevchenko et al. (Shev89 (1989)) on an improved photometry and spectral classification of 35 stars. However, if the PDR is associated with to the Cepheus Bubble, as suggested by Ábrahám , Balázs & Kun (Abraham99 (1999)), then the distance to could be as low as 500$``$700 pc, and of course, the derived masses and luminosities would be significantly reduced.

The masses obtained in this analysis will be discussed and compared to masses obtained from the dust emission, in Section 4 below.

### 3.2 Dust Emission

In order to compare the 850 and 450$`\mu `$m dust emission with the molecular line observations, the 450$`\mu `$m images need to be convolved to the same HPBW as the rest of the observations. However, since the JCMT telescope is far from perfect at 450$`\mu `$m, one first needs to remove the error beam. A model beam was therefore constructed, by fitting three Gaussians to representative beam maps of Uranus. The main beam and the near error lobe were fitted with elliptical Gaussians, while the far error lobe was fitted with a circularly symmetric beam. The relative amplitude of the main beam and the near and far error lobes were found to be 0.925:0.07:0.005. The HPBW of the near error lobe was 37′′$`\times `$ 26′′, with the major axis aligned in the chop direction similarly to the main beam (c.f. Section 2.2), and the far error lobe gave a HPBW $``$ 120′′. The deconvolution of images was tested with both MAXEN (MIRIAD’s maximum entropy task) and CLEAN. Both worked reasonably well, but CLEAN is better in preserving flux, and was therefore adopted. We also constructed a model beam for 850$`\mu `$m and and ran CLEAN on the 850$`\mu `$m–images as well. The deconvolved images were restored with a symmetric Gaussian to a 14′′ beam. For 450$`\mu `$m we also made maps restored to a HPBW of 8′′. Figure 4 shows the high resolution maps at 450$`\mu `$m and 850$`\mu `$m.

Even CLEAN, however, is not ideal for recovering smooth extended emission, which could cause to underestimate the 450 $`\mu `$m emission more than at 850$`\mu `$m, since the 450 $`\mu `$m–map is restored from an 8$`^{\prime \prime }.`$5 resolution to 14′′, while the resolution in the 850$`\mu `$m–maps remains roughly the same. Emission in the off position will have an even more severe systematic effect. Since the error beam has a much higher amplitude at 450$`\mu `$m than at 850$`\mu `$m, any emission in the off position will affect the baselevel of the 450 $`\mu `$m–map much more than that at 850$`\mu `$m. This leads to a systematic underestimate of the flux densities at 450$`\mu `$m. This was tested by adding a constant flux level of 5% of the peak value at 450$`\mu `$m to the “raw” 450$`\mu `$m–map of LkH$`\alpha `$ 234. The baselevel corrected map was run through CLEAN and restored to 14′′. When the flux densities in the 14′′ beam are compared for positions more than 40′′ away from LkH$`\alpha `$ 234, it is found that adding a zero-level of 5% of the peak flux density increases the flux level by more than a factor of two compared to the uncorrected map. Since the base level can easily have an error of 5%, this means that the flux densities for faint extended emission are uncertain by more than a factor of two. In this case the 450$`\mu `$m flux densities are more than likely to be systematically underestimated.

Two-dimensional Gaussian fits were then used to determine the size and total flux of the two “compact” sources in the field, i.e. LkH$`\alpha `$ 234 SMM 1 and the protostellar source, LkH$`\alpha `$ 234 SMM 2, $``$ 28′′ North-West of LkH$`\alpha `$ 234 SMM 1. Here the flux estimates are much more accurate, since the sources are compact and the emission from the surrounding cloud ridge can be substracted. LkH$`\alpha `$ 234 SMM 1 emission is found to be extended, with a size of $``$ 5′′ (Table 3). This agrees well with earlier 800$`\mu `$m and 450$`\mu `$m maps (Sandell & Weintraub San\_Wei94 (1994)) as well as with 1.3 mm maps by Fuente et al. (Fuente98 (1998)) and Henning et al. (Henning98 (1998)). The source LkH$`\alpha `$ 234 SMM 2 is not obvious in the 1.3 mm data of Fuente et al. (Fuente98 (1998)) nor in the map published by Henning et al. (Henning98 (1998)), but both 1.3 mm maps show an extension towards SMM 2 in good agreement with our SCUBA maps. By using the size derived from our SCUBA maps, we derive an integrated 1.3 mm flux of 0.2 $`\pm `$ 0.04 Jy for the map published by Fuente et al. (Fuente98 (1998)). We are interested to derive the $`\beta `$-index (where $`\beta `$ reflects the change in dust emissivity or, equivalently, the change in optical depth, with frequency, $`\kappa \tau \nu ^\beta `$) and the total (gas + dust) mass. For this, we have fitted the 1.3 mm data point and the two SCUBA data points with a simple isothermal model (see e.g. Sandell Sandell00 (2000)). Since both IRAS and Kuiper Airborne Observatory (KAO) data (Bechis et al. Bechis78 (1978); Harvey, Wilking & Joy Harvey84 (1984); Di Francesco et al. diFrancesco98 (1998)) employ large beam sizes that are likely to include hot dust from the surrounding reflection nebulosity as well, we have not included FIR data in our least squares fit, although we constrained our LkH$`\alpha `$ 234 SMM 1 model to the peak dust temperature derived by Harvey, Wilking & Joy (Harvey84 (1984)). The results of the fit are given in Table 4 and also shown in Figure 5. LkH$`\alpha `$ 234 SMM 2 is even less constrained. However, by looking at gas temperatures, either from our own CO observations or at the temperature estimates from the NH3 map by Güsten & Marcaide (Gusten86 (1986)), we find the gas temperatures to be $``$ 30 K. Therefore the fit presented in Table 4 was constrained to 30 K. The mass we derive for the dust disk or envelope surrounding LkH$`\alpha `$ 234 SMM 1 (Dent et al, Dent89 (1989)), is about half of that derived by Henning et al. (Henning98 (1998)) after it is adjusted to the same distance. The difference is largely due to the difference in the adopted mass opacity for the dust.

In order to compare the gas mass estimates derived from our 13CO and C18O data with dust mass estimates, submillimeter flux densities have been derived from our submillimeter maps convolved with a 14′′ beam. The results are given in Table 5. If one uses the baseline un-adjusted 450$`\mu `$m data and ignores the differences in dust properties and dust temperatures between the two embedded stars by assuming that the dust is at the same temperature as the CO, the $`\beta `$-values appear unplausibly low. If instead the baseline adjusted 450$`\mu `$m data were to be used one would obtain a $`\beta `$ of $``$ 1.7. However, since we cannot accurately correct for the emission in the off position, we prefer to assume a $`\beta `$-index and only use the flux density measurements at 850$`\mu `$m, which are much less affected by emission in the off position. The CO temperatures are likely to be dominated by the gas heated by the PDR, and they therefore give only an upper limit to the mass averaged gas temperatures. The dust emission, however, is more optically thin and the submillimeter dust emission contains also emission from the colder dust inside the dense cloud ridge. This regime is more accurately measured by NH3, which shows kinetic temperatures of $``$ 25 K (Güsten & Marcaide Gusten86 (1986)) along the cloud ridge. In the second entry for each offset lower gas temperatures are assumed, which are typically smaller than the CO temperatures. The two compact sources in the ridge have been substracted from the map in order to separate the dust cloud from the disk/envelope emission. However, the total mass values in Table 5 include also the contribution from LkH$`\alpha `$ 234 SMM 1 and LkH$`\alpha `$ 234 SMM 2.

### 3.3 NGC 7129 FIRS 2

NGC 7129 FIRS 2 was first discovered by Bechis et al. (Bechis78 (1978)) as a cold FIR source associated with a 13CO column density peak south of LkH$`\alpha `$ 234. It was found to coincide with an H2O maser (Cesarsky et al. Cesarsky78 (1978); Rodríguez et al. Rodriguez80 (1980); Sandell & Olofsson Sandell81 (1981)), suggesting that it is a young protostellar source, yet the source is invisible in the optical and near-IR. Additional FIR mapping by Harvey, Wilking & Joy (Harvey84 (1984)) confirmed the discovery by Bechis et al. (Bechis78 (1978)) and detected the source at both 50 and 100 $`\mu `$m, but not at 20 $`\mu `$m or shorter wavelengths. Edwards & Snell (Edwards83 (1983)) found that FIRS 2 was associated with a CO outflow, which further strengthened its association with a young, intermediate PMS star. Jenness, Padman & Scott (Jenness95 (1995)) included FIRS 2 in their survey of FIR cores in the vicinity of H2O masers and reported a flux density of 5.2 Jy at 800 $`\mu `$m and 30 Jy at 450 $`\mu `$m. Eiróa, Palacios & Casali (Eiroa98 (1998)) proposed that FIRS 2 is an intermediate mass Class 0 object, based on submillimeter photometry at 2 mm, 1.3 mm, 1.1 mm and 800 $`\mu `$m as well as from high resolution reconstructed IRAS data. They detected FIRS 2 at 25, 60 and 100 $`\mu `$m, but not at 12 $`\mu `$m. They fitted the SED of FIRS 2 with a standard isothermal graybody fit and found a $`\beta `$-index of 0.9 and a dust temperature of 35 K corresponding to a total mass of 6 M.

Our SCUBA data give a more accurate position for FIRS 2 (Table 3), which agrees well with the H2O position (uncertainty $`\pm `$ 10′′). The source is resolved with a FWHM of $``$ 3.7′′. The surrounding cloud core is also seen, but because most observations were made using relatively short chop throws ($``$ 100′′), the emission of the cloud core is not well mapped. The data is supplemented with maps and photometry made by us using UKT14 on JCMT. Our UKT14 flux densities are generally lower than those found by Jenness, Padman & Scott (Jenness95 (1995)) and Eiróa, Palacios and Casali (Eiroa98 (1998)), but agree much better with the well-calibrated SCUBA data. We find flux densities of 0.48 $`\pm `$ 0.10 Jy, 1.18 $`\pm `$ 0.03 Jy, 1.69 $`\pm `$ 0.05 Jy and 3.66 $`\pm `$ 0.22 Jy for 2 mm, 1.3 mm, 1.1 mm and 800 $`\mu `$m, respectively. The flux densities at 2 mm and 1.3 mm are based on photometry and therefore also include a small contribution from the surrounding cloud. The 1.1 mm and 800 $`\mu `$m values are derived from Gaussian fits to calibrated maps, i.e. essentially the same way SCUBA data were analyzed.

The isothermal fit (Fig. 5) predicts a $`\beta `$ $``$ 1 for dust temperature of 42 K. The fit underestimates slightly the FIR data points, which are observed with a large beam size and therefore also include emission from the surrounding cloud core. Our $`\beta `$-index and mass (Table 4) agree rather well with the value derived by Eiróa, Palacios & Casali (Eiroa98 (1998)) and demonstrate that FIRS 2 is a cold protostellar source with a luminosity of $``$ 400 L.

## 4 Discussion

### 4.1 A Gas and Dust Comparison

The 13CO emission (Fig. 1) and dust emission (Fig. 4) show broadly similar distributions. The molecular cavity is also a dust cavity, although the boundary appears broader in dust continuum emission than in 13CO emission. The presence of a dust cavity suggests that the interior is not filled with HI, but has been substantially cleared of dust and gas. Bechis et al. (1978) showed that radiation pressure from older stars in the cavity is capable of doing this. Moreover, recent 21cm HI observations by H. E. Matthews et al. (2001, in preparation) in the NGC 7129 region support this hypothesis, by showing that HI in the cavity has a low number density, of 85 cm-3, and it is in pressure equilibrium with the molecular gas in the ridge (both having P/k $``$ 2x105 cm-3K).

In both 13CO and dust, the mapped region is seen to contain two main clumps, one surrounding the star LkH$`\alpha `$ 234 and extending to the northwest, the second being about 2 to the south. Although the dust and gas maps are consistent on large scales, there are significant differences on smaller scale. The southern clump, centered on the peak NGC 7129 CO4, has an extension to the east with no counterpart in the continuum maps. LkH$`\alpha `$ 234 SMM 1 (identified as NGC 7129 CO1) dominates the emission in submillimeter, while NGC 7129 CO2 (SMM 2 in Table 3) is almost equally bright in 13CO and much fainter in dust continuum. The peak NGC 7129 CO3 is not present in the dust emission maps. The two 13CO peaks, NGC 7129 CO2 and CO3, have essentially equal intensity, so the difference in continuum emission between the two is surprising. The answer may lie in the proximity of the star SVS 13 to SMM 2. Radiation from SVS 13 may be responsible for maintaining the dust temperature of $``$ 30 K, which we find for this source. The projected distance of SMM 2 from SVS 13 is about three times that of NGC 7129 CO2 and therefore the stellar radiation field is less intense. However, the absence of dust emission at NGC 7129 CO3 cannot be explained by temperature alone. It must have a lower density than NGC 7129 CO2.

A quantitative comparison of gas and dust is provided by the masses listed in Tables 2 and 5. The mass from the CO isotopomers was obtained from LTE expressions (Section 3.1), while the mass from the 850 $`\mu `$m flux required values of $`\beta `$, Td and a gas/dust mass ratio of 100 (Hildebrand Hilde83 (1983)):

$$M(M_{})=1.88\times 10^2\times F\times (\frac{\lambda }{0.25})^{\beta +3}\times [e^{\frac{14.4}{\lambda T_d}}1]\times d^2,$$
(3)
where F is the flux (in Jy) at the wavelength $`\lambda `$ (in mm), d is the distance to the source (in kpc) and the mass opacity coefficient is 10 cm2/g at 250$`\mu `$m. Two values for dust masses are given in Table 5, for different values of T and $`\beta `$. We will refer below to those found with $`\beta `$ = 2, which is the typical value expected from theoretical grain models (Hildebrand Hilde83 (1983)). In view of all the assumptions that have been made, we consider the mass values derived from the CO and from the dust to be in good agreement. In the worst cases, the difference is a factor of three. In all these cases, the agreement could be improved by making defensible adjustments to the dust temperature. The results are also consistent with little or no CO depletion in the NGC 7129 ridge, perhaps due to the enhanced radiation field, which maintains grain temperatures above the 16 K CO freeze-out temperature.

### 4.2 The nature of the submillimeter sources LkHa 234 SMM 1 and SMM 2

Early studies of NGC 7129 (Bechis et al. Bechis78 (1978); Harvey, Wilking & Joy Harvey84 (1984)) indicated that the FIR emission could be explained by dust heated by the B stars that illuminate the reflection nebulae. The FIR emission peaked on LkH$`\alpha `$ 234, which was also known to be a weak radio continuum source (Bertout & Thum Bertout82 (1982)), associated with an H2O maser (Cesarsky et al. Cesarsky78 (1978)) and driving a high velocity CO outflow (Edwards & Snell Edwards83 (1983)). The second FIR peak, FIRS 2 was also associated with an H2O maser and a high velocity CO outflow, but lacked an optical counterpart (Section 3.3). Since LkH$`\alpha `$ 234 was associated with HH objects (Hartigan & Lada Hart85 (1985)), coincided with a free-free emission source at 15 and 22 GHz, with compact dust emission at 3 mm (Wilking, Mundy & Schwartz Wilking86 (1986)) and drove a highly collimated optical jet (Ray et al. Ray90 (1990)), this identification appeared rather secure.

However, Skinner, Brown & Steward (Skinner93 (1993)) noticed that the free-free emission was offset by $``$ 1$`^{\prime \prime }.`$7 from the optical position of LkH$`\alpha `$ 234 if they adopted the more accurate optical position given by Herbig & Bell (Herbig88 (1988)) or by 2$`^{\prime \prime }.`$2 when compared with the astrometric position of Clements & Argyle (Clements84 (1984)). Near-IR polarimetric imaging by Weintraub, Kastner & Mahesh (Wein94 (1994)) showed that there is a deeply embedded companion $``$ 3′′ North-West of LkH$`\alpha `$ 234, which was invisible in total intensity maps at K band. High resolution mid-IR imaging by Cabrit et al. (Cabrit97 (1997)) identified the embedded companion star, IRS 6, at 2$`^{\prime \prime }.`$7 NW of LkH$`\alpha `$ 234 and showed that it has a steeply rising spectrum, making IRS 6 as bright as LkH$`\alpha `$ 234 at 17$`\mu `$m. IRS 6 therefore coincides with the radio free-free emission source, with two H2O maser spots, and with the 3 mm dust emission (Wilking, Mundy & Schwartz Wilking86 (1986); Sandell, unpublished data). Cabrit et al. (Cabrit97 (1997)) show that IRS 6 illuminates an arc-shaped reflection nebula with very red colors and that it drives an H2 jet at a position angle (p.a.) of 226. The p.a. of the H2-jet differs from that of the optical \[SII\]-jet, which has a p.a. of 252 (Ray et al. Ray90 (1990)). Whether the much larger \[SII\] jet also originates from IRS 6 is not clear (Cabrit et al. Cabrit97 (1997)). Although Mitchell & Matthews (Mitch94 (1994)) report that the p.a. of the CO outflow agrees with the optical jet, this is not true if the p.a. is determined from the symmetry axis of the large scale red CO outflow, which is $``$ 230. Neither does the p.a. of the optical jet agree with the location of the two HH objects GGD 32 and HH 103 SW of LkH$`\alpha `$ 234 and IRS 6, which lie in the blue-shifted lobe of the CO outflow.

The submillimeter position for LkH$`\alpha `$ 234 SMM 1 determined from our 850 $`\mu `$m map (Table 3) falls between LkH$`\alpha `$ 234 and IRS 6, but since our positional accuracy is $``$ 1′′, the emission could peak on either star. However, since the dust emission at 3 mm peaks on IRS 6, it is plausible to assume that the submillimeter emission also peaks on IRS 6. It therefore appears that it is the cold, heavily obscured IRS 6 that drives the CO outflow and which is associated with the massive dust disk or envelope we see in the submillimeter. This disk-like structure is orthogonal to the near-IR H2 jet to within a few degrees. From our simple isothermal modeling (Section 3.2) we derive a luminosity of $``$ 9 102 L, which would make IRS 6 an early B star, i.e. a young Herbig Be star still heavily enshrouded by dust. It is clear that LkH$`\alpha `$ 234 as well as BD$`+`$651638 and even BD$`+`$651637 and SVS 13 contribute to the 50 and 100 $`\mu `$m emission, but the submillimeter emission is totally dominated by IRS 6 and the second submillimeter source LkH$`\alpha `$ 234 SMM 2. There is some dust emission in the molecular cavity as well, especially in the vicinity of BD$`+`$651638. This dust is hotter than the dust in the cloud ridge, which is also seen from the spectral index plot (Fig. 6).

The submillimeter source LkH$`\alpha `$ 234 SMM 2 does not have a known optical or IR counterpart. It is $``$ 12′′ to the east of SVS 13, which does not appear to be associated with any dust emission (Fig. 4). SMM 2 does not coincide with the CO hotspot seen by Mitchell & Matthews (Mitch94 (1994)) either, which is $``$ 10′′ \- 15′′ SE of the submillimeter source. SMM 2 coincides within errors with the NGC 7129 CO2 peak (Fig. 1) and is also apparent in the CS J=3$``$2 map presented by Fuente et al. (Fuente98 (1998)) as a dense gas condensation. Two of our C18O J=3$``$2 spectra, at offsets 3 and 6 in Table 2, are close to SMM 2, i.e. $``$ 5′′ north and south of the submillimeter peak. They show slightly elevated excitation temperatures (see Table 2), but the higher kinetic temperature is more likely due to the PDR rim than to any embedded source. The linewidths, especially in C18O, are narrow and there is no evidence for outflow. SMM 2 therefore has the appearance of a pre-stellar core. The virial mass is $``$ 8.2 M, assuming a linewidth of 1 km s-1, and a diameter of $``$ 14′′ i.e. several times larger than the mass we derive from our submillimeter observations. However, it is clear that the the hot gas in the cavity has compressed the whole ridge in which SMM 2 is embedded. We also note that our derived value for $`\beta `$ as well as our assumed dust temperature is rather uncertain and our mass estimate is at most accurate to a factor of two. SMM 2 could therefore either be a pre-stellar core, or a protostar.

### 4.3 Dust properties

To investigate dust properties across the region, a map of the spectral index $`\alpha `$ ($`F_\nu \nu ^\alpha `$) was made, using the 850 and 450$`\mu `$m deconvolved images smoothed to 14′′ resolution. Explicitly,

$$\alpha =\mathrm{log}(F_{\lambda _1}/F_{\lambda _2})/\mathrm{log}(\lambda _2/\lambda _1),$$
(4)
where $`\lambda _1`$ and $`\lambda _2`$ are the two wavelengths used in our observations.

We know that the maps are distorted by chopping into the cloud emission (Section 3.2). This will affect the 450 $`\mu `$m more than the 850 $`\mu `$m flux, because the error lobe is larger at 450 $`\mu `$m. The 450 $`\mu `$m map was therefore corrected by adding a constant level of 5% of the peak value to the map, before convolving it to a resolution of 14′′. The resulting map is shown in Figure 6.

The spectral index map has two trends. Firstly, a comparison with the 450 $`\mu `$m or 850 $`\mu `$m map shows a general anticorrelation between submillimeter continuum flux and spectral index. The highest flux positions are minima of $`\alpha `$, and the value of $`\alpha `$ increases towards regions of weaker continuum emission. The $`\alpha `$ minima are seen at locations correlating with dense gas clumps, LkH$`\alpha `$ 234 SMM 1, SMM 2 and the flux peak near the offset ($``$45$`^{\prime \prime }.`$0,$``$85$`^{\prime \prime }.`$0). The curving ridge of 450 $`\mu `$m emission (Fig. 4) can be followed throughout its length in the $`\alpha `$ map. Secondly, the map also shows a steep rise in $`\alpha `$ towards the molecular cavity. This large gradient in $`\alpha `$ occurs along the entire transition region.

It is important to determine whether the observed increase in $`\alpha `$ is due to the change in temperature through the PDR, or in the dust opacity index, $`\beta `$. Because Fν $``$ $`\kappa _\nu `$βBν(Td), the spectral index, $`\alpha `$, and the dust opacity, $`\beta `$, can be related as follows:

$$\alpha =\beta +2+\gamma ,$$
(5)
where $`\gamma `$ can be thought of as a correction to the Rayleigh-Jeans form of the Planck function and is given by (Visser et al. 1998)

$$\gamma =1\mathrm{log}[(e^{h\nu _{850}/kT}1)/(e^{h\nu _{450}/kT}1)]/\mathrm{log}(\nu _{850}/\nu _{450}).$$
(6)
An increase in $`\alpha `$ may be due to an increase in the dust temperature, T, and/or in $`\beta `$. The temperature is expected to be higher in the cavity because of the more intense radiation field. Our analysis in Section 3.1 showed that the gas temperature is higher towards the ridge. Is it possible to understand the larger values of $`\alpha `$ only as due to higher dust temperatures? The answer is no because $`\gamma `$ approaches zero at large T so that the maximum value of $`\alpha `$ is $`\beta `$ \+ 2. Values of $`\alpha `$ larger than 4 require $`\beta `$$``$2, even if T is large. While the expected increase in Td toward the cavity will result in higher values of $`\alpha `$, the increase in Td must be accompanied by an increase in $`\beta `$.

Ideally, the above relations can yield the exact value of $`\beta `$ if one knows the spectral index $`\alpha `$ and the temperature of the dust. In practice, the determination of the true value of $`\beta `$ is restricted by the current uncertainties in the submillimeter fluxes, especially at short wavelengths, which will reflect into uncertainties in $`\alpha `$. The dust will be more optically thick at 450 $`\mu `$m than at 850 $`\mu `$m (since the optical depth changes with frequency as $`\tau \nu ^\beta `$). Thus at 450 $`\mu `$m we do not see as deeply into the densest hot cores as one does at 850 $`\mu `$m. Another more severe restriction is our ability to accurately measure the submillimeter fluxes, especially at short wavelengths. As mentioned before, we tried to remedy this situation by adding a constant zero-level to the 450 $`\mu `$m flux, but this is only a first order correction. In reality, the dust structure will not be flat, but will have its own substructure, similar to that of the molecular cloud. Therefore, the numerical values of $`\alpha `$ in regions of low flux are highly uncertain.

In summary, we have confidence in the values of $`\alpha `$ $``$ 3 that we find at high flux levels. We also believe in the trend of increasing $`\alpha `$ towards lower flux values. In the dense cores we find $`\beta `$ $`<`$ 2, as suggested by the minima in the spectral index plot as well as our dust fits in Figure 5. Since small grains produce a high $`\beta `$-index and large grains lead to smaller values of $`\beta `$, our results are consistent with the growth of grains in regions of high density.

These present observations of NGC 7129 support the conclusions of other similar studies that high density clumps or cores have smaller values of $`\beta `$ than does the surrounding gas (Visser et al. Viss98 (1998); Johnstone $`\&`$ Bally Johnstone99 (1999); Sandell et al. Sandell99 (1999)). Values of $`\beta `$ $``$ 2 outside the flux peaks have been also reported in reflection nebulae surrounding Herbig Ae/Be stars (Whitcomb et al. Whitcomb81 (1981)). Our $`\alpha `$ map goes further, however, in showing a general correspondence between flux and $`\alpha `$ over an extended area, and in showing a steep increase in $`\alpha `$ (and, therefore, in $`\beta `$) along an extended PDR.

## 5 SUMMARY

New high resolution maps of the NGC 7129 ridge have been obtained, both in submillimeter continuum emission and in 13CO J=3$``$2 line emission. The main results of this study are:

1. A sharp boundary between the molecular ridge and cavity is clearly seen in the 13CO map. A previously detected narrow filament of optical and near infrared emission follows closely the cloud-cavity boundary. Strong emission in the v=1$``$0 S(1) line of H2 from the filament, if due to fluorescence, supports the interpretation of the boundary as a PDR.

2. The 13CO emission consists of two main cloud cores connected by a bridge of gas. The northern cloud core has three 13CO peaks, while the southern one has a single extended peak.

3. The 850 $`\mu `$m and 450 $`\mu `$m maps show a very good morphological agreement with the 13CO map. Both data sets show a sharp boundary (cloud rim) toward the cavity where the two young stars BD +65 1637 and BD +65 1638 are located.

4. Three compact sources are seen in the continuum maps (SMM 1, SMM 2, and FIRS 2 in Table 3). Isothermal graybody fits to the measured spectral energy distributions have been made, by combining our SCUBA data with earlier observations. The resulting masses, dust temperatures, luminosities and dust emissivities are given in Table 4. The SMM 1 continuum peak is identified with the deeply embedded infrared source IRS 6, which according to these present fits is a young Herbig Be star.

5. For nine offsets around the LkH$`\alpha `$ 234 star, the 13CO/C18O intensity ratio was used to obtain the mass of gas in an 14′′ beam. The masses derived from CO are found to be consistent with masses derived from the dust continuum emission. Within our analysis uncertainties the results are consistent with little or no CO depletion in this region.

6. A map of the 850$`\mu `$m/450$`\mu `$m spectral index, $`\alpha `$, shows a general correlation between the $`\alpha `$ minima and the submillimeter continuum flux peaks. The low values of $`\alpha `$ in the peaks also correlate with low values of dust opacity indices, $`\beta `$, in the range 1 to 1.6, derived for the compact source, suggesting grain growth. The spectral index $`\alpha `$ rises steeply towards the molecular cavity, along the entire length of the boundary region. Consequently, the dust opacity index in this region is also high, $`\beta `$ $``$ 2. This suggests very hot, small size grains in the PDR ridge, whose mantles have been evaporated by the intense UV radiation.

The James Clerk Maxwell Telescope is operated on a joint basis between the United Kingdom Particle Physics and Astronomy Research Council (PPARC), the Netherlands Organization for the Advancement of Pure Research (ZWO), the Canadian National Research Council (NRC), and the University of Hawaii (UH). This research was partly supported by an operating grant from the Natural Sciences and Engineering Research Council of Canada.


text:
# Superconductivity in the Pseudogap State due to Fluctuations of Short – Range Order.

## I Introduction.

The pseudogap state, observed mainly in the underdoped region of the phase diagram of high – temperature superconducting (HTSC) cuprates, leads to a wide range of anomalies of physical properties of these system both in normal and superconducting states . There are two main theoretical scenarios to explain these anomalies. The first is based upon the model of Cooper pairs formation already above the temperature of superconducting transition  with phase coherence appearing only for $`T<T_c`$. The second assumes that the pseudogap state is induced by fluctuations of antiferromagnetic (AFM) short – range order which exist in the underdoped region of the phase diagram . Recently a number of experiments appeared which rather convincingly favor the second scenario .

Most theoretical papers in the field are dedicated to the study of the models of the pseudogap state in normal phase for $`T>T_c`$. In Refs.  we proposed a very simple (toy) exactly solvable model of the pseudogap state, which assumes the existence of certain “hot” (flat) patches on the Fermi surface. Within this model we derived the Ginzburg – Landau expansion for different types of Cooper pairing  and studied the anomalies of superconducting state for $`T<T_c`$  induced by fluctuations of AFM short – range order. In these papers we dealt with (over) simplified model of static Gaussian fluctuations of short – range order with infinite correlation length, which allowed us to obtain an exact analytical solution for the pseudogap state. In real systems the correlation length of AFM fluctuations is finite and relatively short . The main aim of the current paper is to generalize the main results of Refs.  to the case of finite correlation lengths of fluctuations of AFM short – range order and to analyze the dependence of the main physical characteristics of the superconducting state on this correlation length, as well as on the effective width of the pseudogap.

## II Model of the Pseudogap State.

We consider the simplified model of the pseudogap  based upon the picture of well – developed fluctuations of AFM short – range order and is in some respects close to the “hot – spot” model of Ref. . We assume the Fermi surface of two – dimensional electronic system as shown in Fig.1. In fact such a Fermi surface was observed in a number of ARPES – experiments on HTSC – cuprates . Note that the assumption of the flatness of these patches is non crucial for our model, but significantly simplify calculations, which, in principle, can be done also in more realistic “hot-spot” approach with apparently similar results. The model Fermi surface as in Fig.1 has already been applied to HTSC – cuprates in Refs. , where the details of microscopic criteria for the existence of antiferromagnetic and superconducting phases were analyzed. Here we just assume a kind of phenomenological model with static Gaussian fluctuations of short – range order with correlation function (structure factor) of the following form :

$$S(𝐪)=\frac{1}{\pi ^2}\frac{\xi ^1}{(q_xQ_x)^2+\xi ^2}\frac{\xi ^1}{(q_yQ_y)^2+\xi ^2}$$
(1)
where $`\xi `$ – is correlation length of fluctuations, while the scattering vector is taken to be $`Q_x=\pm 2k_F`$$`Q_y=0`$ or $`Q_y=\pm 2k_F`$$`Q_x=0`$, anticipating that these fluctuations are incommensurate. The factorized form of (1) was first introduced in Ref.  and leads to great simplification of all calculations. At the same time we stress that in practice it practically coincides with the usual isotropic Lorentzian in the most important region of $`|𝐪𝐐|<\xi ^1`$ .

Less justified physically is the assumption of the static nature of fluctuations which can be reasonable only at high enough temperatures . At low temperatures, including the superconducting phase, spin dynamics can be very important, e.g. for the microscopics of Cooper pairing according to the picture of “nearly – antiferromagnetic” Fermi liquid . However, we assume that our static approximation can be sufficient for qualitative analysis of the influence of superconductivity upon superconductivity, which will be described below within the standard BCS – like approximation.

Let us write down the interaction of electrons with AFM fluctuations in the following form:

$$V_{eff}=(2\pi )^2W^2S(𝐪)$$
(2)
where $`W`$ is determining the energy scale (width) of the pseudogap. We assume that only electrons from the flat (“hot”) of the Fermi surface interact with AFM fluctuations, so that $`W`$ is effectively non zero only for these electrons . Note that we completely neglect the spin structure of interaction which can be rather easily taken into account , but makes calculations more cumbersome. In this sense, strictly speaking our discussion is more appropriate for the case of electrons interacting with fluctuations short – range order of CDW rather than SDW (AFM) type. We also assume that this simplification is relatively unimportant for the analysis of qualitative effects of pseudogap on superconductivity.

The factorized form of the correlator (1) and of interaction (2) leads to the one – dimensional nature of scattering on fluctuations. In the limit of infinite correlation length $`\xi \mathrm{}`$ this model acquires an exact solution . For finite $`\xi `$ we can construct “nearly exact” solution , directly generalizing the one – dimensional approach proposed in Ref. . In this case we can (approximately) sum the whole diagram series for one – particle Green’s function of electrons from the flat parts of the Fermi surface (where the “nesting” condition for the electronic spectrum $`\xi _{𝐩\pm 𝐐}=\xi _𝐩`$ is satisfied).

For the contribution of an arbitrary diagram for the self – energy of an electron of the $`N`$-th order in (2) we use the following Ansatz :

$`\mathrm{\Sigma }^{(N)}(\epsilon _n𝐩)=W^{2N}{\displaystyle \underset{j=1}{\overset{2N1}{}}}G_{0k_j}(\epsilon 𝐩),`$ (3)
$`G_{0k_j}(\epsilon _n𝐩)={\displaystyle \frac{1}{i\epsilon _n(1)^j\xi _𝐩+ik_j\kappa }}`$ (4)
where $`\kappa =v_F\xi ^1`$ ($`v_F`$ – Fermi velocity), $`k_j`$ – is the number of interaction lines surrounding the $`j`$-th (from the beginning) electronic line in diagram, $`\epsilon _n=2\pi T(n+1/2)`$ (in the following we write expressions for the case of $`\epsilon _n>0`$). Thus the contribution of an arbitrary diagram is determined, in fact, only by some set of integer numbers $`k_j`$. An arbitrary diagram with intersections of interaction lines is seen to be equal to some diagram without intersections and the contribution of all diagrams with intersection can be taken into account with the help of some combinatorial factors $`v(k_j)`$, attributed to interaction lines on diagrams with no intersections . In the model of incommensurate fluctuations which we only consider:

$$v(k)=\{\begin{array}{cc}\frac{k+1}{2}& \text{for odd }k\\ \frac{k}{2}& \text{for even }k\end{array}$$
(5)
As a result we obtain the following recursion procedure (continuous fraction representation) to calculate the one – particle Green’s function $`G(\epsilon _n𝐩)`$ of electrons from “hot patches” :

$$G_k(\epsilon _n𝐩)=\frac{1}{i\epsilon _n(1)^k\xi _𝐩+ik\kappa W^2v(k+1)G_{k+1}(\epsilon _n𝐩)};G(\epsilon _n𝐩)G_0(\epsilon _n𝐩)$$
(6)
Diagrammatically this procedure is shown in Fig.2.

Ansatz (4) for the contribution of an arbitrary diagram of the $`N`$-th order, in general case, is not exact . However, in two – dimensional case we can explicitly demonstrate the topologies of the Fermi surface when (4) becomes exact , for the general case we can show  that this representation in some sense overestimates the role of finiteness of correlation length $`\xi `$ in the given order of perturbation theory. For one – dimensional case, when this problem is particularly serious , it appears that for incommensurate fluctuations calculations of the density of states using (4) quantitatively almost ideally reproduce  the results of an exact numerical solution given in Refs.  ***In the case of one – dimensional problem with commensurate fluctuations Ansatz (4) does not describe a certain weak Dyson – type singularity of the density of states near the center of the pseudogap , while outside the region of this singularity it also gives rather good quantitative description of exact results. Note that in two -dimensional case Dyson – singularity, most probably, is just absent.. In the limit of $`\xi \mathrm{}`$ Ansatz (4) reduces to an exact solution of Ref. , while in the limit of $`\xi 0`$ for the fixed value of $`W`$ it describes the physically correct limit of free electrons.

Outside “hot patches” electrons in our model do not interact with fluctuations and the Green’s function remains a free one:

$$G(\epsilon _n𝐩)=G_{00}(\epsilon _n𝐩)=\frac{1}{i\epsilon _n\xi _𝐩}$$
(7)
This model leads to non Fermi – liquid like behavior of the spectral density on “hot patches” and to the smeared pseudogap in the density of states (cf. similar results in “hot spots” model ). On “cold” parts of the Fermi surface we just get the usual Fermi – liquid (free electron) behavior.

## III Gorkov Equations for a Superconductor with Pseudogap.

In Refs.  we have analyzed the anomalies of superconducting state in an exactly solvable model of the pseudogap state induced by AFM short – range order fluctuations with infinite correlation length ($`\xi \mathrm{}`$). In particular, in Ref.  it was shown that AFM fluctuations may lead to strong fluctuations of superconducting order parameter (energy gap $`\mathrm{\Delta }`$) breaking the standard assumption of self – averaging gap  which allows independent averaging (over configurations of the random field of static short – range order fluctuations) of the order parameter $`\mathrm{\Delta }`$ and different combinations of electronic Green’s functions, entering the basic equations of the theory. Usually, the possibility of such independent averaging is supported by the following argument : the value of $`\mathrm{\Delta }`$ changes on characteristic scale of the order of superconducting coherence length $`\xi _0v_F/\mathrm{\Delta }_0`$ of BCS – theory, while Green’s functions typically vary on much shorter lengths of the order of interatomic distances. Naturally, the last assumption becomes wrong if a new length scale $`\xi \mathrm{}`$ appears in electronic system. However, in case of AFM correlation length $`\xi \xi _0`$ (i.e. when AFM fluctuations correlate on distances much smaller than the size of Cooper pairs) the self – averaging property of $`\mathrm{\Delta }`$ remains, being broken only for $`\xi >\xi _0`$. For this reason all our analysis below will be done assuming the self – averaging nature of the energy gap, allowing us to use the standard approach of the theory of disordered superconductors (mean – field approach in terms of Ref. ). Thus, we temporarily put aside the very interesting problem of superconductivity in the absence of self – averageness of the order parameter. Note that in real HTSC – systems we, in fact, usually have $`\xi \xi _0`$, so that these materials are in the region of parameters most difficult for the theory.

Similarly to Refs.  we assume that superconducting pairing is due to the following simplest attraction potential:

$$V(𝐩,𝐩^{})=V(\varphi ,\varphi ^{})=Ve(\varphi )e(\varphi ^{}),$$
(8)
where $`\varphi `$ is the polar angle determining the direction of electronic momentum $`𝐩`$ in highly conducting plane, while for $`e(\varphi )`$ we take the simple model dependence:

$$e(\varphi )=\{\begin{array}{cc}1\hfill & (\text{ }s\text{-wave pairing})\hfill \\ \sqrt{2}\mathrm{cos}(2\varphi )\hfill & (\text{ }d\text{-wave pairing})\hfill \end{array}.$$
(9)
The attraction constant $`V`$ is, as usual, non zero in some shell of the width of $`2\omega _c`$ around the Fermi level ($`\omega _c`$ – is characteristic frequency of the quanta responsible for electron attraction). In this case superconducting gap takes the form:

$$\mathrm{\Delta }(𝐩)\mathrm{\Delta }(\varphi )=\mathrm{\Delta }e(\varphi ).$$
(10)
In the following for brevity we shall denote gap simply as $`\mathrm{\Delta }`$ instead of $`\mathrm{\Delta }(\varphi )`$ writing down explicit angular dependence only where it is necessary.

In the superconducting state perturbation theory over the interaction with AFM fluctuations (1) has to be built on “free” normal and anomalous Green’s functions of a superconductor:

$$G_{00}(\epsilon _n𝐩)=\frac{i\epsilon _n+\xi _𝐩}{\epsilon _n^2+\xi _𝐩^2+|\mathrm{\Delta }|^2};F_{00}^+(\epsilon _n𝐩)=\frac{\mathrm{\Delta }^{}}{\epsilon _n^2+\xi _𝐩^2+|\mathrm{\Delta }|^2}$$
(11)
In the model with flat patches on the Fermi surface the electronic spectrum on the patches orthogonal to $`p_x`$ is: $`\xi _𝐩=v_F(|p_x|p_F)`$ as electron velocity $`𝐯`$ is perpendicular to $`p_y`$ (everything is symmetric for patches orthogonal to $`p_y`$). Thus, in the case of $`s`$-wave pairing, when $`\mathrm{\Delta }`$ is independent of the direction of electronic momentum, for the model interaction (1), (2) the problem becomes one – dimensional. In case of $`d`$-wave pairing situation is more difficult as even on the flat patches orthogonal to $`p_x`$ the value of $`\mathrm{\Delta }(\varphi )`$ depends on $`p_y`$ (and symmetrically on patches orthogonal to $`p_y`$). For this reason, in the analysis of $`d`$-wave pairing it is convenient to introduce instead of (1) the fluctuation correlator of the form

$$S(𝐪)=\frac{1}{\pi }\left\{\frac{\xi ^1}{(q_x2p_F)^2+\xi ^2}\delta (q_y)+\frac{\xi ^1}{(q_y2p_F)^2+\xi ^2}\delta (q_x)\right\}$$
(12)
In this case interaction does not change $`p_y`$ or $`p_x`$ on flat patches orthogonal to $`p_x`$ or $`p_y`$ and the problem again reduces to one – dimensional.

Now we can formulate an analogue of approximation (4) also for superconducting state. Some details of derivation are given in the Appendix A. The contribution of an arbitrary diagram of $`N`$-th order over interaction (2) to normal or anomalous Green’s function is given by the product of $`N+1`$ “free” normal $`G_{0k_j}`$ or anomalous $`F_{0k_j}^+`$ Green’s functions with renormalized frequencies and gaps (cf. below). Here $`k_j`$ – is the number of interaction lines surrounding the $`j`$-th (from the beginning of the diagram) electronic line. As in normal phase, the contribution of an arbitrary diagram is defined by the set of integers $`k_j`$ and each diagram with intersection of interacting lines is equal to a certain diagram of the same order without intersections. Thus again we can use only diagrams without intersections of interaction lines accounting diagrams with intersections by the same combinatorial factors $`v(k)`$ (attributed to interaction lines) as in the normal phase. As a result we obtain diagrammatic analogue of Gorkov equations  shown in Fig.3. Accordingly we get a system of two recursion relations for normal and anomalous Green’s functions:

$`G_k=G_{0k}+G_{0k}\stackrel{~}{G}G_kG_{0k}\stackrel{~}{F}F_k^+F_{0k}\stackrel{~}{G}^{}F_k^+F_{0k}\stackrel{~}{F}^+G_k`$ (13)
$`F_k^+=F_{0k}^++F_{0k}^+\stackrel{~}{G}G_kF_{0k}^+\stackrel{~}{F}F_k^++G_{0k}^{}\stackrel{~}{G}^{}F_k^++G_{0k}^{}\stackrel{~}{F}^+G_k`$ (14)
where

$$\stackrel{~}{G}=W^2v(k+1)G_{k+1};\stackrel{~}{F}^+=W^2v(k+1)F_{k+1}^+$$
(15)
$$G_{0k}(\epsilon _n𝐩)=\frac{i\epsilon _n+(1)^k\xi _𝐩}{\stackrel{~}{\epsilon }_n^2+\xi _𝐩^2+|\stackrel{~}{\mathrm{\Delta }}|^2};F_{0k}^+(\epsilon _n𝐩)=\frac{\stackrel{~}{\mathrm{\Delta }}^{}}{\stackrel{~}{\epsilon }_n^2+\xi _𝐩^2+|\stackrel{~}{\mathrm{\Delta }}|^2}$$
(16)
and we introduced renormalized frequency $`\stackrel{~}{\epsilon }`$ and gap $`\stackrel{~}{\mathrm{\Delta }}`$ as:

$$\stackrel{~}{\epsilon }_n=\eta _k\epsilon _n;\stackrel{~}{\mathrm{\Delta }}=\eta _k\mathrm{\Delta };\eta _k=1+\frac{k\kappa }{\sqrt{\epsilon _n^2+|\mathrm{\Delta }|^2}}$$
(17)
similar to those appearing for superconductors with impurities .

From (14)-(17) it is easy to obtain the system of recursion relations for real and imaginary parts of normal Green’s function, as well as for anomalous Green’s function:

$`ImG_k={\displaystyle \frac{\stackrel{~}{\epsilon }Im\stackrel{~}{G}}{(\stackrel{~}{\epsilon }Im\stackrel{~}{G})^2+((1)^k\xi _𝐩+Re\stackrel{~}{G})^2+|\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{F}|^2}}`$ (18)
$`ReG_k={\displaystyle \frac{(1)^k\xi _𝐩+Re\stackrel{~}{G}}{(\stackrel{~}{\epsilon }Im\stackrel{~}{G})^2+((1)^k\xi _𝐩+Re\stackrel{~}{G})^2+|\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{F}|^2}}`$ (19)
$`F_k^+={\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^{}+\stackrel{~}{F}^+}{(\stackrel{~}{\epsilon }Im\stackrel{~}{G})^2+((1)^k\xi _𝐩+Re\stackrel{~}{G})^2+|\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{F}|^2}}`$ (20)
Let us introduce the following notations:

$$ImG_k=\epsilon _nJ_k;ReG_k=(1)^k\xi _𝐩R_k;F_k^+=\mathrm{\Delta }^{}f_k$$
(21)
Then we can see that recursion relations for $`J_k`$ and $`f_k`$ just coincide, so that $`J_k=f_k`$. Finally we obtain the following system of recursion relations for $`J_k`$ and $`R_k`$:

$`J_k={\displaystyle \frac{\eta _k+W^2v(k+1)J_{k+1}}{(\epsilon _n^2+\mathrm{\Delta }^2)(\eta _k+W^2v(k+1)J_{k+1})^2+\xi _𝐩^2(1+W^2v(k+1)R_{k+1})^2}}`$ (22)
$`R_k={\displaystyle \frac{1+W^2v(k+1)R_{k+1}}{(\epsilon _n^2+\mathrm{\Delta }^2)(\eta _k+W^2v(k+1)J_{k+1})^2+\xi _𝐩^2(1+W^2v(k+1)R_{k+1})^2}}`$ (23)
Then the normal and anomalous Green’s functions of a superconductor are determined through $`R_0`$ and $`J_0`$:

$$ImG=\epsilon _nJ_0;ReG=\xi _𝐩R_0;F^+=\mathrm{\Delta }^{}J_0$$
(24)
and represent the complete sum of perturbation series for electron in a superconductor interacting with AFM short – range order fluctuations.

## IV Critical Temperature and Temperature Dependence of the Gap.

Energy gap of a superconductor is determined by the equation:

$$\mathrm{\Delta }(𝐩)=T\underset{𝐩^{}}{}\underset{\epsilon _n}{}V_{sc}(𝐩,𝐩^{})F(\epsilon _n𝐩^{})$$
(25)
On the flat parts of the Fermi surface the anomalous Green’s function is defined by (24) and recursion relations (23). On the rest (“cold” part) of the Fermi surface the scattering on AFM fluctuations is absent (in our model), so that the anomalous Green’s there function is given by (11). As a result, for the case of $`s`$-wave pairing, with the account of (9), the gap equation (25) takes the form:

$$1=\lambda \left\{\stackrel{~}{\alpha }T\underset{\epsilon _n}{}_{\omega _c}^{\omega _c}𝑑\xi J_0(\epsilon _n\xi )+(1\stackrel{~}{\alpha })_0^{\omega _c}𝑑\xi \frac{th\frac{\sqrt{\xi ^2+\mathrm{\Delta }^2}}{2T}}{\sqrt{\xi ^2+\mathrm{\Delta }^2}}\right\}$$
(26)
where $`\lambda =VN_0(0)`$ is dimensionless coupling constant of pairing interaction ($`N_0(0)`$ – free – electron density of states at the Fermi level ), $`\stackrel{~}{\alpha }=4\alpha /\pi `$, where $`\alpha `$ is the angular size of a flat patch on the Fermi surface (cf. Fig.1). In numerical calculations below we, rather arbitrarily, choose $`\stackrel{~}{\alpha }=2/3`$, i.e. $`\alpha =\pi /6`$, which is close e.g. to the data of Ref. .

In case of $`d`$-wave pairing we have to take into account the angular dependence of the gap (10), so that Eq. (25) becomes:

$$1=\lambda \frac{4}{\pi }\left\{T_0^\alpha 𝑑\varphi e^2(\varphi )\underset{\epsilon _n}{}_{\omega _c}^{\omega _c}𝑑\xi J_0(\epsilon _n\xi )+_\alpha ^{\pi /4}𝑑\varphi e^2(\varphi )_0^{\omega _c}𝑑\xi \frac{th\frac{\sqrt{\xi ^2+\mathrm{\Delta }^2e^2(\varphi )}}{2T}}{\sqrt{\xi ^2+\mathrm{\Delta }^2e^2(\varphi )}}\right\}$$
(27)
In Fig.4 we show the calculated (from Eq. (26)) temperature dependences of the gap for the case of $`s`$-wave pairing and for different values of correlation length (parameter $`\kappa =v_F\xi ^1`$) of fluctuations. For the case of $`d`$-wave pairing results are qualitatively similar.

Equation for superconducting critical temperature $`T_c`$ immediately follows from (26), (27) as $`\mathrm{\Delta }0`$. In this case $`J_0(\mathrm{\Delta }0)`$ is independent of $`\varphi `$ and is the same both for $`s`$ and $`d`$-wave pairing. Then $`T_c`$ – equation takes the following form:

$$1=\lambda \left\{\alpha _{eff}T_c\underset{\epsilon _n}{}_{\omega _c}^{\omega _c}𝑑\xi J_0(\epsilon _n\xi ;\mathrm{\Delta }0)+(1\alpha _{eff})_0^{\omega _c}𝑑\xi \frac{th\frac{\xi }{2T_c}}{\xi }\right\}$$
(28)
where an “effective” size of flat patches on the Fermi surface is defined as:

$$\alpha _{eff}=\{\begin{array}{cc}\stackrel{~}{\alpha }\hfill & (\text{ }s\text{-wave pairing})\hfill \\ \stackrel{~}{\alpha }+\frac{1}{\pi }\mathrm{sin}(\pi \stackrel{~}{\alpha })\hfill & (\text{ }d\text{-wave pairing})\hfill \end{array}.$$
(29)
Calculated dependences of $`T_c`$ on the width of the pseudogap $`W`$ and correlation length (parameter $`\kappa =v_F\xi ^1`$) are shown in Fig.5 ($`T_{c0}`$ – transition temperature in the absence of pseudogap).

The general qualitative conclusion is the same as in Refs. : pseudogap suppresses superconductivity due to a partial “dielectrization” of electronic spectrum on “hot” parts of the Fermi surface. This suppression effect is maximal for $`\kappa =0`$ (infinite correlation length of AFM fluctuations)  and weakens as correlation length becomes shorter. These results are in full accordance with experimental phase diagram of HTSC – cuprates.

Let us stress once again that all our results are valid in case of self – averaging superconducting order parameter (mean – field approach of Ref. ), which is valid for not very large correlation lengths $`\xi <\xi _0`$, where $`\xi _0`$ – is superconducting coherence length (the size of Cooper pairs at $`T=0`$). For $`\xi \xi _0`$ important effects due to non self – averaging gap fluctuations appear, leading e.g to characteristic “tails” in temperature dependence of the average gap for $`T_c<T<T_{co}`$ .

## V Cooper Instability. Recurrence Procedure for the Vertex Part.

It is well known that the critical temperature can also be determined from Cooper instability of the normal phase:

$$1V\chi (0,0)=0$$
(30)
where the generalized Cooper susceptibility is defined by the graph shown in Fig.6. Here we have to calculate the “triangular” vertex part accounting for interaction with AFM fluctuations. For the similar one – dimensional problem (and for real frequencies $`T=0`$) the appropriate recurrence procedure was formulated in Ref. . For our two – dimensional model this procedure was used in Ref.  to calculate optical conductivity. Generalization to Matsubara frequencies is rather straightforward. Below, for definiteness, we assume $`\epsilon _n>0`$. Then we obtain:

$`\mathrm{\Gamma }_{k1}(\epsilon _n,\epsilon _n,𝐪)=1+`$ (31)
$`+W^2v(k)G_k\overline{G}_k\left\{1+{\displaystyle \frac{2ik\kappa }{2i\epsilon _n(1)^kv_FqW^2v(k+1)(G_{k+1}\overline{G}_{k+1})}}\right\}\mathrm{\Gamma }_k(\epsilon _n,\epsilon _n,𝐪)`$ (32)
$`\mathrm{\Gamma }(\epsilon _n,\epsilon _n,𝐪)\mathrm{\Gamma }_0(\epsilon _n,\epsilon _n,𝐪)`$ (33)
where $`G_k=G_k(\epsilon _n𝐩+𝐪)`$ and $`\overline{G}_k=G_k(\epsilon _n,𝐩)`$ are calculated from (6).

To find $`T_c`$ we have to know the vertex part at $`𝐪=0`$. Then $`\overline{G}_k=G_k^{}`$ and the vertex $`\mathrm{\Gamma }_k`$ becomes real, which considerably simplifies (33). Using notations similar to (21), from (6) and (33) we get:

$$\mathrm{\Gamma }_{k1}=1+W^2v(k)\frac{J_k}{1+W^2v(k+1)J_{k+1}}\mathrm{\Gamma }_k$$
(34)
while for $`R_k`$ and $`J_k`$ we have recursion relations given by (23) with $`\mathrm{\Delta }=0`$.

There exists the following exact relation similar to the Ward identity (the proof of this relation will be given below):

$$G(\epsilon _n𝐩)G(\epsilon _n𝐩)\mathrm{\Gamma }(\epsilon _n,\epsilon _n,0)=(\xi _𝐩^2R_0^2(\epsilon _n\xi _𝐩)+\epsilon _n^2J_0^2(\epsilon _n\xi _𝐩))\mathrm{\Gamma }_0(\epsilon _n,\epsilon _n,0)J_0(\epsilon _n\xi _𝐩)=\frac{ImG(\epsilon _n𝐩)}{\epsilon _n}$$
(35)
Numerical analysis fully confirms this relation, demonstrating the self – consistency of our recursion relations for one – particle Green’s function and vertex part Note that an analytic proof of this relation from direct comparison of recursion procedures for Green’s function and vertex part is non obvious, to say the least.. As $`J_0(\mathrm{\Delta }0)`$ coincides with $`J_0`$ in normal phase, the relation (35) leads to $`T_c`$ – equation obtained from Cooper instability (30):

$$1=\lambda \left\{\alpha _{eff}T_c\underset{\epsilon _n}{}_{\omega _c}^{\omega _c}𝑑\xi (\xi _𝐩^2R_0^2(\epsilon _n\xi _𝐩)+\epsilon _n^2J_0^2(\epsilon _n\xi _𝐩))\mathrm{\Gamma }_0(\epsilon _n,\epsilon _n,0)+(1\alpha _{eff})_0^{\omega _c}𝑑\xi \frac{th\frac{\xi }{2T_c}}{\xi }\right\}$$
(36)
being the same as Eq. (28), obtained by linearization of the gap equation, despite seemingly different recursion procedures used to obtain these equations.

## VI Ginzburg – Landau Expansion.

In Ref.  we derived the Ginzburg – Landau expansion in exactly solvable model of the pseudogap with infinite correlation length of AFM fluctuations. Here we generalize these results to the case of finite correlation lengths.

Let us write the Ginzburg – Landau expansion for the difference of free energies of superconducting and normal state in the following form:

$$F_sF_n=A|\mathrm{\Delta }_𝐪|^2+q^2C|\mathrm{\Delta }_𝐪|^2+\frac{B}{2}|\mathrm{\Delta }_𝐪|^4,$$
(37)
where $`\mathrm{\Delta }_𝐪`$ – is the amplitude of the Fourier – component of the order parameter

$$\mathrm{\Delta }(\varphi ,𝐪)=\mathrm{\Delta }_qe(\varphi ).$$
(38)
Now (37) is defined by diagrams of loop – expansion for the free energy of an electron in the field of fluctuations of the order parameter with small wave – vector $`𝐪`$ .

We express the coefficients of Ginzburg – Landau expansion as:

$$A=A_0K_A;C=C_0K_C;B=B_0K_B,$$
(39)
where $`A_0`$, $`C_0`$ and $`B_0`$ denote the standard expressions for these coefficients in the case of isotropic $`s`$-wave pairing:

$$A_0=N_0(0)\frac{TT_c}{T_c};C_0=N_0(0)\frac{7\zeta (3)}{32\pi ^2}\frac{v_F^2}{T_c^2};B_0=N_0(0)\frac{7\zeta (3)}{8\pi ^2T_c^2},$$
(40)
Then all the anomalies of the model under consideration, connected with the appearance of the pseudogap, are contained in dimensionless coefficients $`K_A`$, $`K_C`$ and $`K_B`$. In the absence of pseudogap all these coefficients are equal to 1, only in case of $`d`$-wave pairing we have $`K_B=3/2`$. Thus for $`d`$\- wave pairing we shall appropriately normalize $`K_B`$, giving numerical results for $`\stackrel{~}{K}_B=2/3K_B`$.

Consider the generalized Cooper susceptibility shown in Fig.6.

$$\chi (𝐪0;T)=T\underset{\epsilon _n}{}\underset{𝐩}{}G(\epsilon _n𝐩+𝐪)G(\epsilon _n𝐩)e^2(\varphi )\mathrm{\Gamma }(\epsilon _n,\epsilon _n,𝐪)$$
(41)
Using (35) we can express the coefficients $`K_A`$ and $`K_C`$ as:

$`K_A={\displaystyle \frac{\chi (𝐪0;T)\chi (00;T_c)}{A_0}}=`$ (42)
$`=\alpha _{eff}{\displaystyle \frac{T_c}{TT_c}}\left\{T{\displaystyle \underset{\epsilon _n=\pi T(2n+1)}{}}{\displaystyle _{\omega _c}^{\omega _c}}𝑑\xi J_0(\epsilon _n\xi )T_c{\displaystyle \underset{\epsilon =\pi T_c(2n+1)}{}}{\displaystyle _{\omega _c}^{\omega _c}}𝑑\xi J_0(\epsilon _n\xi )\right\}+1\alpha _{eff}`$ (43)
$`K_C=\underset{q0}{lim}{\displaystyle \frac{\chi (𝐪0;T_c)\chi (00;T_c)}{q^2C_0}}=`$ (44)
$`={\displaystyle \frac{32\pi ^2T_c^3}{7\zeta (3)v_Fq^2}}\alpha _{eff}\{{\displaystyle \underset{\epsilon _n=\pi T_c(2n+1)}{}}{\displaystyle _{\omega _c}^{\omega _c}}d\xi J_0(\epsilon _n\xi )`$ (45)
$`{\displaystyle \underset{\epsilon =\pi T_c(2n+1)}{}}{\displaystyle _{\omega _c}^{\omega _c}}d\xi G(\epsilon _n,\xi +{\displaystyle \frac{1}{2}}v_Fq)\mathrm{\Gamma }(\epsilon _n,\epsilon _n,q)G(\epsilon _n,\xi {\displaystyle \frac{1}{2}}v_Fq)\}+1\alpha _{eff}`$ (46)
Situation with coefficient $`B`$, in general case, is more complicated. Important simplifications appear if we limit ourselves in the order of $`|\mathrm{\Delta }_q|^4`$, as is usually done, by considering only the case of $`q=0`$. Then the coefficient $`B`$ can be determined directly from the anomalous Green’s function $`F`$, for which we already have the recurrence procedure (23), (24). Let us consider the diagrammatic expansion of the anomalous Green’s function, shown in Fig.7(a). From this it becomes clear that:

$$\underset{\mathrm{\Delta }0}{lim}\frac{F(\epsilon _n𝐩)}{\mathrm{\Delta }}=G(\epsilon _n𝐩)G(\epsilon _n𝐩)+\mathrm{}=G(\epsilon _n𝐩)G(\epsilon _n𝐩)\mathrm{\Gamma }(\epsilon _n,\epsilon _n,0)$$
(47)
which, by the way, immediately proves (35) (use also (24)). Thus, for the two – particle loop $`\chi (0,0)`$ we get:

$$\chi (0,0)=T\underset{𝐩}{}\underset{\epsilon _n}{}\underset{\mathrm{\Delta }0}{lim}\frac{F(\epsilon _n𝐩)}{\mathrm{\Delta }}=T\underset{𝐩}{}\underset{\epsilon _n}{}J_0(\mathrm{\Delta }=0)$$
(48)
For the “four – tail” diagram of Fig.7(b), determining the coefficient $`B`$, in the same way we obtain:

$$T\underset{𝐩}{}\underset{\epsilon _n}{}\underset{\mathrm{\Delta }0}{lim}\frac{\frac{F(\epsilon _n𝐩)}{\mathrm{\Delta }}lim_{\mathrm{\Delta }0}\frac{F(\epsilon _n𝐩)}{\mathrm{\Delta }}}{|\mathrm{\Delta }|^2}=T\underset{𝐩}{}\underset{\epsilon _n}{}\underset{\mathrm{\Delta }0}{lim}\frac{J_0(\mathrm{\Delta })J_0(\mathrm{\Delta }=0)}{|\mathrm{\Delta }|^2}$$
(49)
where $`J_0(\mathrm{\Delta })`$ is defined by the recursion procedure (23). Finally, for the dimensionless coefficient $`K_B`$ we have:

$$K_B=\alpha _B\frac{8\pi ^2T_c^3}{7\zeta (3)}\underset{\epsilon _n}{}_{\omega _c}^{\omega _c}𝑑\xi \underset{\mathrm{\Delta }0}{lim}\frac{J_0(\mathrm{\Delta }=0)J_0(\mathrm{\Delta })}{|\mathrm{\Delta }|^2}+1\alpha _B$$
(50)
where

$$\alpha _B=\{\begin{array}{cc}\stackrel{~}{\alpha }\hfill & (\text{ }s\text{-wave pairing})\hfill \\ \stackrel{~}{\alpha }+\frac{4}{3\pi }\mathrm{sin}(\pi \stackrel{~}{\alpha })+\frac{1}{6\pi }\mathrm{sin}(2\pi \stackrel{~}{\alpha })\hfill & (\text{ }d\text{-wave pairing})\hfill \end{array}.$$
(51)
These expressions allow direct numerical calculations of the coefficients $`K_A,K_C,K_B`$. In Fig.8, for example, we present the calculated dependence of $`K_C`$ on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F/xi^1`$). The appropriate dependences of $`K_A`$ and $`K_B`$ are qualitatively quite similar. In particular, for $`\kappa =0`$ we have just $`K_B=K_C`$ .

## VII Physical Properties of Superconductors with Pseudogap.

Ginzburg – Landau expansion defines two characteristic lengths of superconductors: the coherence length and penetration depth of magnetic field.

The coherence length for given temperature $`\xi (T)`$ determines the characteristic scale of inhomogeneities of the order parameter $`\mathrm{\Delta }`$:

$$\xi ^2(T)=\frac{C}{A}.$$
(52)
In the absence of the pseudogap:

$`\xi _{BCS}^2(T)={\displaystyle \frac{C_0}{A_0}},`$ (53)
$`\xi _{BCS}(T)0.74{\displaystyle \frac{\xi _0}{\sqrt{1T/T_c}}},`$ (54)
where $`\xi _0=0.18v_F/T_c`$. In our model:

$$\frac{\xi ^2(T)}{\xi _{BCS}^2(T)}=\frac{K_C}{K_A}.$$
(55)
The dependences of $`\xi ^2(T)/\xi _{BCS}^2(T)`$ on the width of the pseudogap $`W`$ and correlation length of fluctuations (parameter $`\kappa `$) for the case $`d`$-wave pairing are shown in Fig.9. Note that the changes of coherence length are relatively small.

For the penetration depth of a superconductor without the pseudogap we have:

$$\lambda _{BCS}(T)=\frac{1}{\sqrt{2}}\frac{\lambda _0}{\sqrt{1T/T_c}},$$
(56)
where $`\lambda _0^2=\frac{mc^2}{4\pi ne^2}`$ is penetration depth at $`T=0`$. In general case:

$$\lambda ^2(T)=\frac{c^2}{32\pi e^2}\frac{B}{AC}.$$
(57)
Then for our model:

$$\frac{\lambda (T)}{\lambda _{BCS}(T)}=\left(\frac{K_B}{K_AK_C}\right)^{1/2}.$$
(58)
Graphical dependences of penetration depth for the case of $`d`$-wave pairing are shown in Fig.10.

Near $`T_c`$ the upper critical magnetic field $`H_{c2}`$ is defined via Ginzburg – Landau coefficients as:

$$H_{c2}=\frac{\varphi _0}{2\pi \xi ^2(T)}=\frac{\varphi _0}{2\pi }\frac{A}{C},$$
(59)
where $`\varphi _0=c\pi /e`$ is magnetic flux quantum. Then the slope of the upper critical field close to $`T_c`$ is defined as:

$$\left|\frac{dH_{c2}}{dT}\right|_{T_c}=\frac{24\pi \varphi _0}{7\zeta (3)v_F^2}T_c\frac{K_A}{K_C}.$$
(60)
Graphic dependences of the slope of the upper critical field $`\left|\frac{dH_{c2}}{dT}\right|_{T_c}`$, normalized to the slope at $`T_{c0}`$, on the effective width of the pseudogap $`W`$ and parameter of correlation length $`\kappa `$ for the case of $`d`$-wave pairing are shown in Fig.11. It is seen that the slope for large enough correlation lengths rapidly drops with the width of the pseudogap. However, for short enough correlation lengths we can observe even some weak growth of this parameter for small values of the pseudogap width. For fixed pseudogap width the slope of $`H_{c2}`$ significantly grows as correlation length becomes smaller.

Finally, let us consider the discontinuity of specific heat at the transition point:

$$\frac{C_sC_n}{\mathrm{\Omega }}=\frac{T_c}{B}\left(\frac{A}{TT_c}\right)^2,$$
(61)
where $`C_s,C_n`$ are specific heats of superconducting and normal states, $`\mathrm{\Omega }`$ – sample volume. For $`T=T_{c0}`$ (in the absence of pseudogap, $`W=0`$):

$$\left(\frac{C_sC_n}{\mathrm{\Omega }}\right)_{T_{c0}}=N(0)\frac{8\pi ^2T_{c0}}{7\zeta (3)}.$$
(62)
Then the normalized discontinuity of specific heat in our model can be expressed as:

$$\frac{(C_sC_n)_{T_c}}{(C_sC_n)_{T_{c0}}}=\frac{T_c}{T_{c0}}\frac{K_A^2}{K_B}.$$
(63)
Appropriate dependences on effective width of the pseudogap $`W`$ and parameter of correlation length $`\kappa `$ for the case of $`d`$-wave pairing are shown in Fig.12. It is seen that specific heat discontinuity rapidly drops with the growth of the pseudogap width and grows as correlation length of AFM fluctuations becomes smaller.

For $`s`$-wave superconductor the dependences of physical properties are quite similar, the only change is in larger scale of $`W`$ for which the appropriate changes appear, corresponding to larger stability of isotropic superconductors to partial “dielectrization” of electronic spectrum due to pseudogap formation on “hot patches” of the Fermi surface .

Among the physical characteristics, analyzed above, relatively detailed experimental data are available for specific heat discontinuity . In complete agreement with our conclusions, specific heat discontinuity for $`Bi2212`$ rapidly drops as the system moves to the underdoped region, where the width of the pseudogap grows. According to Ref.  the width of the pseudogap (our parameter $`2W`$) changes from the values of the order of $`700K`$ for hole concentration $`p=0.05`$ to the values of the order of $`T_c100K`$ near optimal concentration $`p=0.16`$ and drops to zero for $`p=0.19`$. A relation between the drop of specific heat discontinuity and the growth of the pseudogap width is clearly observed. Unfortunately, we do not know detailed enough data on concentration dependence of correlation length of fluctuations and appropriate dependences of physical characteristics of superconductors in pseudogap state. Qualitatively it is obvious that correlation length grows as system moves to the underdoped region, so that the drop of specific heat discontinuity is natural also from this point of view.

## VIII Conclusion.

In this paper we continued the study of anomalies of superconducting state in the framework of rather crude model of the pseudogap state of two – dimensional electronic system , which is, however, in qualitative agreement with a number of observed anomalies of electronic structure of underdoped HTSC – cuprates. In Refs.  we considered rather unrealistic limit of infinite correlation length of AFM short – range order fluctuations, which allowed us to obtain an exact analytic solution. Here we generalized our model to realistic case of finite correlation lengths, taking into account, as in Refs. , all diagrams of perturbation theory on electron interaction with fluctuations of short – range order. Our analysis was performed in a standard (mean – field in terms of Ref. ) approach, assuming self – averaging property of superconducting order – parameter over fluctuations of the random field of AFM fluctuations. In Ref.  we have shown that this assumption is not justified in the limit of $`\xi \mathrm{}`$. At the same time this assumption is apparently well justified for the case of $`\xi \xi _0`$ (where $`\xi _0`$ is the coherence length of a superconductor at $`T=0`$, i.e. the size of Cooper pairs). Thus, we are left with rather complicated task of the accounting of non self – averaging effects for $`\xi >\xi _0`$. We have already mentioned that in real HTSC systems in most cases $`\xi \xi _0`$, so that effects of non self – averaging superconducting gap, similar to those considered in Ref. , may be very important, leading e.g. to characteristic “tails” in the temperature dependence of the average gap for $`T>T_c`$ and the physical picture of superconducting “drops” of Ref. .

Another serious simplification of our model is the assumption of static (and Gaussian) nature of short – range order of fluctuations. This assumption may be justified only for high enough temperatures $`T\omega _{sf}`$ (where $`\omega _{sf}`$ – is characteristic scale of spin fluctuations) . Thus, the use of static approximation in superconducting state for $`T<T_c`$ is rather doubtful. However, we think that our simplified treatment allows us to describe most important effects of the changes of the electronic spectrum (due to pseudogap formation on “hot patches” of the Fermi surface) upon superconductivity. The account of spin dynamics inevitably requires to drop the simple phenomenology of BCS model and consider the microscopic nature of pairing interaction. It is doubtful that such a program can be realized in near future. In particular, the problem of summation of all perturbation theory diagrams for the interaction with dynamical spin fluctuations seems absolutely hopeless.

This work was partially supported by grants 99-02-16285 of the Russian Foundation of Basic Research and REC-005 of CRDF, as well as by the projects 108-11(00-$`\mathrm{\Pi }`$) of the State Program “Statistical Physics” and 96-051 of the State Program on HTSC.

## A Coordinate representation. Normal and anomalous Green’s functions.

Let us consider some technical details of derivation of recursion relations for Gorkov equations (14) – (17). It is sufficient to limit the analysis to consideration of two flat parts of the Fermi surface, orthogonal to $`p_x`$ – axis, which are connected by the scattering vector $`𝐐=(\pm 2p_F,0)`$. Then the problem becomes purely one – dimensional as the velocity projection $`v_y=0`$ and the linearized electronic spectrum $`\xi _{p_xp_F}=\pm v_Fp_x`$ does not depend on $`y`$-component of electronic momentum. For brevity in the following we just put $`v_F=1`$.

Calculations simplify in coordinate representation , considering the electron propagation if the field of Gaussian AFM fluctuations $`W(x)W^{}(x)`$ (incommensurate case) with correlator:

$$<W^{}(x)W(x^{})>=W^2e^{\kappa |xx^{}|}$$
(A1)
Then electron propagators corresponding to normal and anomalous Green’s functions of a superconductor (11), take the form:

$`G_{00}(x)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dp_x}{2\pi }}e^{ip_xx}G_{00}(p_x)={\displaystyle \frac{i}{2}}\left({\displaystyle \frac{\epsilon _n}{\sqrt{\epsilon _n^2+|\mathrm{\Delta }|^2}}}+\sigma _3sign(x)\right)e^{\sqrt{\epsilon _n^2+|\mathrm{\Delta }|^2}|x|}`$ (A2)
$`F_{00}(x)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dp_x}{2\pi }}e^{ip_xx}F_{00}^+(p_x)={\displaystyle \frac{\mathrm{\Delta }^{}}{\sqrt{\epsilon _n^2+|\mathrm{\Delta }|^2}}}e^{\sqrt{\epsilon _n^2+|\mathrm{\Delta }|^2}|x|}`$ (A3)
where $`\sigma _3=1`$ for right moving particles, and $`\sigma _3=1`$ for left moving particles. Scattering by fluctuations transforms “right” particles into “left” and vice versa. From (A3) it is seen that the particle moving along the path of the length $`l`$ produces the factor $`e^{\sqrt{\epsilon _n^2+|\mathrm{\Delta }|^2}l}`$.

During calculation of specific diagrams it is convenient  to change integration variables from coordinates of interaction vertices $`x_k`$ to lengths of the paths $`l_k`$ traversed by the particle between separate scatterings, fixing the total displacement $`xx^{}`$. With interaction line connecting the vertices $`m`$ and $`n`$ on electronic line we have to associate the factor:

$$W^2exp(\kappa |x_mx_n|)=W^2exp(\kappa |\underset{k=m}{\overset{n1}{}}(1)^kl_k|)$$
(A4)
Integration over all $`l_k`$ is performed from $`0`$ to $`\mathrm{}`$.

It is seen that the finite correlation length of fluctuations in a given diagram leads to some “damping” of the appropriate transition amplitude with distance traversed by an electron. The exact treatment of this effect is difficult, but in Ref.  we used an obvious inequality:

$$exp\left(\kappa |\underset{k=m}{\overset{n1}{}}(1)^kl_k|\right)>exp\left(\kappa \underset{k=m}{\overset{n1}{}}l_k\right)$$
(A5)
and replaced the exponential of (A4) by the exponential from the r.h.s. of (A5). This is equivalent to the replacement of correlator of random fields (A1) by similar expression, where in the exponent we just replace the distance $`|xx^{}|`$ by the total length of the path traversed by the particle between scattering acts at $`x`$ and $`x^{}`$. According to (A5) this procedure somehow overestimates the effect of damping $`\kappa `$ in each diagram of perturbation theory. After such a replacement the diagrams of all orders are easily calculated and precisely reproduce the Ansatz of (4) for the normal phase . We have already mentioned that results obtained in this way, e.g. for the density of states, are in good quantitative agreement with exact numerical simulation of the one – dimensional problem , which gives additional support for our approximation strengthening qualitative estimates of Ref. .

Let us use the same approximation during the analysis of diagrams of perturbation theory in superconducting phase, which are built upon propagators (A3). In this case the role of interaction reduces just to the appearance of additional factor of $`e^{\kappa l_k}`$ in each normal or anomalous Green’s function (A3), surrounded by the given interaction line, or (which is the same) to the addition of $`\kappa `$ to $`\sqrt{\epsilon _n+|\mathrm{\Delta }|^2}`$ in the exponential of each Green’s function. Making transformation back to the momentum representation it is easily seen that the contribution of an arbitrary diagram of the higher order of perturbation theory is determined by the product of the appropriate number of normal and anomalous Green’s functions of the form:

$$G_{0k}(p)=\frac{i\epsilon _n\frac{\epsilon _k}{\sqrt{\epsilon _n+|\mathrm{\Delta }|^2}}+(1)^k\xi _p}{\epsilon _k^2+\xi _p^2};F_{0k}^+(p)=\frac{\mathrm{\Delta }^{}\frac{\epsilon _k}{\sqrt{\epsilon _n+|\mathrm{\Delta }|^2}}}{\epsilon _k^2+\xi _p^2};$$
(A6)
where $`\epsilon _k=\sqrt{\epsilon _n+|\mathrm{\Delta }|^2}+k\kappa `$, while $`k`$ – is the number of interaction lines, surrounding the given Green’s function. The factor of $`(1)^k`$ is due to scattering transforming “right” particles into “left” and vice versa. Introducing the renormalized frequency and gap as in (17), we can see that (A6) reduces to the standard form (16), which completes the justification of our recursion procedure (14), (17).

Figure Captions:

Fig.1. Fermi surface of two – dimensional system. “Hot patches” are shown by thick lines of the width $`\xi ^1`$.

Fig.2. Diagrammatic representation of recursion relation for one – particle Green’s function.

Fig.3. Diagrammatic representation of recursion relations for Gorkov equations.

Fig.4. Temperature dependence of superconducting energy gap for the case of $`s`$-wave pairing and different values of correlation length (parameter $`\kappa =v_F\xi ^1`$) of AFM fluctuations, calculated for $`\lambda =0.4,\frac{\omega _c}{W}=3`$:

$`\frac{\kappa }{W}`$=0 (1); 1.0 (2); 10.0 (3).

Dashed line — $`\mathrm{\Delta }(T)`$ in the absence of the pseudogap.

Fig.5. Dependence of superconducting transition temperature on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F\xi ^1`$):

$`\frac{\kappa }{W}`$=0.1 (1); 1.0 (2); 10.0 (3).

Dashed line — $`\kappa =0`$ .

At the insert: dependence of $`T_c`$ on $`\kappa `$ for $`\frac{W}{T_{c0}}=5`$.

Fig.6. Diagram for the generalized Cooper susceptibility.

Fig.7. (a) – Diagram series for anomalous Green’s function, dashed lines – AFM fluctuations. (b) – Diagram determining $`K_B`$.

Fig.8. Dependence of $`K_C`$ on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F\xi ^1`$):

$`\frac{\kappa }{W}`$=0.1 (1); 1.0 (2); 10.0 (3).

Dashed line — $`\kappa =0`$ .

At the insert: dependence of $`K_C`$ on $`\kappa `$ for $`\frac{W}{T_{c0}}=5`$.

Fig.9. Dependence of supercoducting coherence length on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F\xi ^1`$):

$`\frac{\kappa }{W}`$=0.1 (1); 1.0 (2); 10.0 (3).

Dashed line — $`\kappa =0`$ .

At the insert: dependence of coherence length on $`\kappa `$ for $`\frac{W}{T_{c0}}=5`$.

Fig.10. Dependence of penetration depth on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F\xi ^1`$):

$`\frac{\kappa }{W}`$=0.1 (1); 1.0 (2); 10.0 (3).

At the insert: dependence of penetration length on $`\kappa `$ for $`\frac{W}{T_{c0}}=5`$.

Fig.11. Dependence of the slope of the upper critical field on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F\xi ^1`$):

$`\frac{\kappa }{W}`$=0.1 (1); 1.0 (2); 10.0 (3).

Dashed line — $`\kappa =0`$ .

At the insert: dependence of the slope of $`H_{c2}`$ on $`\kappa `$ for $`\frac{W}{T_{c0}}=5`$.

Fig.12. Dependence of specific heat discontinuity on the width of the pseudogap $`W`$ and correlation length of AFM fluctuations (parameter $`\kappa =v_F\xi ^1`$):

$`\frac{\kappa }{W}`$=0.1 (1); 1.0 (2); 10.0 (3).

Dashed line — $`\kappa =0`$ .

At the insert: dependence of specific heat discontinuity on $`\kappa `$ for $`\frac{W}{T_{c0}}=5`$.