text:
# Normal-Branch Quasi-Periodic Oscillations during the High Intensity State of Cygnus X-2

## 1 Introduction

Although the low-mass X-ray binary (LMXB) and Z source (Hasinger & van der Klis 1989) Cygnus X-2 is one of the best studied LMXBs, its complex behaviour is far from understood. On time-scales of hours to days, the source traces out its characteristic Z shape pattern in the X-ray colour-colour diagram (CD). The branches of this Z are called, from top to bottom, the horizontal branch (HB), the normal branch (NB), and the flaring branch (FB). On time-scales of weeks to months, the overall intensity varies smoothly by a factor of about 4 (see, e.g., Fig. 1). The exact morphology of the Z track and the position of this track in the CD (and hence the overall hardness of the X-ray spectrum) varies significantly when the overall intensity changes (see, e.g., Kuulkers, van der Klis, & Vaughan 1996 and Wijnands et al. 1997 for overviews of these changes). The first few months of the Rossi X-ray Timing Explorer (RXTE; Bradt, Rothschild, & Swank 1993) All Sky Monitor (ASM) data of Cygnus X-2 suggested an $``$78 day period in the long-term variations (Wijnands, Kuulkers, & Smale 1996). Although this period is still one of the strongest periods, analyses of all RXTE/ASM data now available (up to 2000 October 9) show that the long-term variations are more complex than only one single periodicity (see also Kong, Charles, & Kuulkers 1998 and Paul, Kitamoto, & Makino 2000).

The physical mechanism behind the long-term variations in Cyg X-2 is unknown. Several mechanisms have been proposed, such as variations in the mass accretion rate onto the neutron star, a precessing accretion disc, or a precessing neutron star (see Wijnands et al. 1997 and references therein for a discussion of these possible mechanisms). Although variations in the mass accretion rate are unlikely to cause the long-term variability (these variations are thought to produce the motion of Cyg X-2 along the Z track; see, e.g., Hasinger & van der Klis 1989, Kuulkers et al. 1996, and Wijnands et al. 1997), a precessing accretion disc or a precessing neutron star are still possible. Besides affecting the X-ray count rate and the X-ray spectrum, both mechanisms would most likely also affect the rapid X-ray variability, although it is unclear how exactly the variability would be affected. A detailed study of the rapid X-ray variability during different overall intensities might give new insights into the physics behind the long-term variability of Cyg X-2. However, such a detailed study has not yet been performed. Below we briefly list the results so far available.

On the HB, quasi-periodic oscillations (QPOs) between 15 and 55 Hz are observed (the horizontal branch QPOs or HBOs). These QPOs are also often observed on the upper part of the NB. A detailed comparison of the HBO properties during different overall intensities (hereafter referred to as different intensity states) has not yet been made. Although Wijnands et al. (1997) reported that the HBO is always observed on the HB regardless of the overall intensity of Cygnus X-2, the HBO properties could differ in detail between different intensity states.

On the NB, QPOs with frequencies near 5 Hz can also be observed (the normal-branch oscillations or NBOs; sometimes observed simultaneous with the HBO). Also the NBO properties have not yet been compared between different intensity states, although Wijnands et al. (1997) did not detect the NBO (on the lower NB; with typical upper limits of $`<`$1% rms on the NBO amplitude) during the high intensity states observed using Ginga data of Cyg X-2. The NBOs were observed during lower intensity states.

Another type of QPO (with frequencies between 300 and 900 Hz) was recently discovered (Wijnands et al. 1998) in Cygnus X-2 when this source was on the HB during a medium intensity state: the kHz QPOs. These QPOs were not observable during a higher intensity state (Smale 1998) with upper limits which were smaller than the measured values.

Here, we report the first detection of a NBO in Cyg X-2 during a high intensity state demonstrating that this type of QPO is so far always observable regardless of the overall intensity of the source. The energy dependence of the rms amplitude and the time lags of this NBO are very similar to those observed during lower intensity states.

## 2 Observations and Analysis

Cygnus X-2 was observed by the RXTE proportional counter array (PCA; Jahoda et al. 1996) on several occasions between 1998 June 2–7 for a total of $``$42 ksec (see Table 1 for a log of the observations). Data were obtained in 129 photon energy channels (covering 2–60 keV) with a time resolution of 16 s. Simultaneous data were obtained in two “Single Bit” modes (one energy channel each; time resolution of 128 $`\mu `$s; energy range 2–5.0 and 5.0–13.0 keV), one “Event” mode (64 channels; 16 $`\mu `$s; 13.0–60 keV), and one “Binned” mode (16 channels; 2 ms; 2–13.0 keV). We used the 16 s data to create the CD and the hardness-intensity diagram (HID; see the caption of Fig. 3 for the energy bands used to calculate the colours). We used the Single Bit modes and the Event modes to create 1/16–2048 Hz power spectra to study the rapid X-ray variability above 100 Hz (i.e., to search for kHz QPOs); we used the Binned and the Event modes to create 1/16–256 Hz power spectra and cross spectra to study the rapid X-ray variability below 100 Hz (i.e., to study the NBO and HBO). The strengths of the power spectral components are all given for the energy range 2–60 keV, unless otherwise noted.

## 3 Results

Figure 1 shows the RXTE/ASM light curve of Cyg X-2 clearly showing the long-term X-ray variations. During the observations presented here, Cygnus X-2 was at high count rates. The RXTE/ASM count rates were approximately 60 counts s<sup>-1</sup> (Fig. 2). The upper envelope of the ASM light curves in Figures 1 and 2 corresponds to epochs during which Cyg X-2 was on the NB; the short-lasting “drop-outs” from this upper envelope are excursions into the HB and the FB. These excursions do not define a lower envelope demonstrating that the overall intensity of Cyg X-2 is mainly defined by the count rate on the NB.

The CD and HID of the data are shown in Figure 3. In both diagrams, clear horizontal and normal branches are visible. The FB is not clearly visible in the CD (only as a broadening of the lower part of the normal branch) but in the HID an extended FB is present. The count rate decreased when Cygnus X-2 entered the FB. We selected the power spectra based on the position of the source on the track in the HID. The HID was used instead of the CD because by using the HID it was possible to disentangle the power spectra corresponding to the lower part of the NB and those corresponding to the FB.

In Figure 4, typical 2–60 keV power spectra are shown for different positions of the source on the Z track. On the HB, clear HBOs, their second harmonics, and strong (8%–10% rms) band-limited noise (the Z source low-frequency noise or LFN; Fig. 4a) were present. The frequency of the HBO was $``$30 Hz at the left end of the HB and increased smoothly to 49 Hz at the right end (see also Fig. 3b). The rms amplitude of the HBO at the left end of the HB was 5% and decreased to 4% when the source moved slightly to the right on the HB. The amplitude of the HBO on the rest of the HB stayed approximately constant at 4% rms. The FWHM of the HBO was 5 Hz at the left end of the HB and increased to 20 Hz when the source moved to the right end. Together with the HBO, its second harmonic could be detected from the left end up to halfway the HB, with an rms amplitude decreasing from 4.4% to 2.6% and a FWHM increasing from 24 to 36 Hz.

On the upper part of the observed normal branch (Fig. 4b), the HBO was still visible around 55 Hz (with an rms amplitude of 2%–3% and a FWHM of 20–30 Hz) together with a noise component below 1 Hz following a power law (the Z source very-low frequency noise or VLFN) and a peaked noise component near 5 Hz. Further down the NB and at the beginning of the FB (Fig. 4c), no QPOs were present, but the VLFN still was, together with a weak (several percent) noise component following a cutoff power law (cutoff frequency near 10 Hz; probably the Z source high-frequency noise or HFN). On the left end of the FB only the VLFN remained (Fig. 4d).

So, when using power spectra created for the total RXTE/PCA photon energy range (2–60 keV), no NBO were detected on the NB. However, when we examined the peaked noise component on the upper part of the observed NB (Fig. 4b) in more detail and in different energy ranges, a clear QPO around 5 Hz emerged. Figure 5 shows the power spectra corresponding to the upper normal branch but at different energies. The power spectrum in Figure 5a is the same as the one presented in Figure 4b, although now plotted with the power axis in linear scale in order to show the QPO more clearly. Figures 5b to d show the corresponding power spectra at different energies. At energies above 7.9 keV a very significant (29 $`\sigma `$) QPO is clearly visible with a frequency of 5.37$`{}_{0.03}{}^{}{}_{}{}^{+0.08}`$ Hz, a FWHM of 4.4$`\pm `$0.2 Hz, and an rms amplitude (7.9–60 keV) of 6.7%$`\pm `$0.1%. Note that the FWHM of the QPO is such that according to the usual criterion for a noise component to be called a QPO (the ratio of the FWHM to the frequency is $`<`$ 0.5), this NBO was not truly a QPO. However, below we will use the term QPO for this peaked noise feature. By comparing Figures 5a and d, it is clear that this QPO is only visible as a broad peaked noise component near 5 Hz when combining the data of the whole RXTE/PCA energy range. The QPO was also visible at low ($`<`$5 keV) photon energies (Fig. 5b), but not at energies between 5.0 and 6.4 keV.

After discovering the NBO at high photon energies, we reexamined all the power spectra corresponding to the different locations of the source on the Z track. We used the 7.9–60 keV energy range instead of the 2–60 keV one. The 7.9–60 keV power spectra on the horizontal branch were very similar to those in the 2–60 keV band. As said above, on the upper part of the observed normal branch a broad NBO was now clearly visible. The frequency (5.4 Hz), the FWHM (4 Hz), and the rms amplitude (7%) of the NBO did not change significantly when the source moved down the NB. Halfway down the NB, the QPO became much broader (with a FWHM of 7–13 Hz) and the NBO decreased slightly in amplitude (to 5%–6% rms). On the lower part of the NB (near the NB-FB vertex), the NBO could not be detected anymore with upper limits on the rms amplitude of 3%–4% (see also Fig. 3b for the region where the NBO could be detected). At the NB-FB vertex and on the FB the 7.9–60 keV power spectra were very similar to the 2–60 keV power spectra.

Figure 6a shows the count rate in the energy bands used to calculate the rms amplitude of the NBO. The rms amplitude of the NBO as a function of photon energy is plotted in Figure 6b. In order to calculate these rms amplitudes, we only used the data of the upper part of the observed NB for which the NBO had a FWHM which was lower then 5 Hz (see also Fig. 3b). The NBO decreased from about 3% rms at the lowest energies to $`<`$2% rms between 5 and 6.4 keV. Above 6.4 keV, the strength of the NBO rapidly increased to 10% at 15 keV. At higher energies the strength seemed to level of, although a firm conclusion cannot be made.

For the same data, we determined the time lags of the NBO between the different energy bands. In order to calculate the time lags of the NBO, we calculated the lags in the 4.4 Hz interval (corresponding to the FWHM of the NBO) centred on the peak frequency of the QPO (5.37 Hz). We subtracted the average cross-vector between 80 and 250 Hz in order to minimize the dead-time effects on the lags induced by Poisson fluctuations (see van der Klis et al. 1987). The resulted time lags are presented in Figure 6c. We used the 2–3.1 keV energy band as reference band. The time lags between this reference band and the bands below 5 keV were consistent with zero. Above 5 keV, the time lags increase rapidly to $``$70 ms for energies between 8 and 9 keV. At higher energies, the time lags remained constant near 70 ms. This means that the hard photons lag the soft photons by about 70 ms (2.4 radians or 140°).

We searched for kHz QPOs in the total energy band and in the energy range 5.0–60 keV (as used by Wijnands et al. 1998 when kHz QPOs were present in Cyg X-2). None were found with typical (assuming a FWHM of 150 Hz) 95% confidence upper limits of 2%–3% (1%–2%) on the HB, 2%–3% (1%–1.5%) on the NB, and 2.5%–3.5% (1.5%–2.0%) on the FB in the energy range 5.0–60 keV (the limits enclosed within parentheses are for the total RXTE/PCA energy range). Two X-ray bursts were observed. We searched for (nearly) coherent oscillations, as reported in several other LMXBs (e.g., Strohmayer et al. 1996; Smith, Morgan, & Bradt 1997), but none were found with upper limits of 4%–7% (2–60 keV).

## 4 Discussion

We have detected the NBO on the upper part of our observed normal branch of Cygnus X-2 when this source was at relative high overall intensities (Fig. 1). The shape of the Z track in the CD and the HID was very similar to those observed with Ginga for Cygnus X-2 in the high intensity state (Wijnands et al. 1997; i.e., compare their Fig. 1 with our Fig. 3). Therefore, we conclude that Cygnus X-2 was in a high overall intensity state during the observations reported here. This conclusion is strengthened by the fact that our power spectra on the lower normal branch resemble those obtained by Wijnands et al. (1997) in the high intensity state, especially the non-detection of the NBO on the lower part of the NB (compare their Figs. 8g–i with our Figs. 4c and d). Wijnands et al. (1997) were unable to search for the NBO further upwards on the NB (the middle and upper part of the NB) due to the fact that they did not have enough high time resolution data at those positions on the Z track. Our RXTE/PCA data are the first high time resolution data available on the middle and upper part of the NB of Cygnus X-2 during the high intensity state.

Our first clear detection of the NBO in the high intensity state of Cygnus X-2 shows that all power spectral components which are observed below 100 Hz in Cygnus X-2 are present at high overall intensities. So far, the only clear difference in the timing properties of Cygnus X-2 during different intensity states are the detection of the kHz QPOs on the HB during a medium intensity state (Wijnands et al. 1997) but not during high intensity states (Smale 1998 and this paper). More detections of the kHz QPOs in Cygnus X-2 are needed in order to determine the exact reason(s) why the kHz QPOs are only present during preferred intensity states.

Wijnands et al. (1997) pointed out that a precessing accretion disc, one of the proposed explanations for the long-term X-ray variations of Cyg X-2 (e.g., Kuulkers et al. 1996; Wijnands et al. 1996), could not easily account for the non-detection of the NBO during the high intensity state. Our detection of the NBO during this state now shows that the arguments used by Wijnands et al. (1997) are not valid any more when only considering the NBO. However, when also considering the kHz QPOs, it is clear that the precessing accretion disc model has still difficulties explaining the difference in rapid X-ray timing variability (i.e., the kHz QPOs) when Cyg X-2 is at different intensity states. In the precessing accretion disc scenario, it is possible that more matter is blocking the emission region during low intensity states (hence the lower overall intensity) than during high intensity states. If this the case then during these low states, the amplitudes of the kHz QPOs should be significantly lower than in the high states due to more scattering of the radiation. However, we see exactly the opposite: no kHz QPOs are detected during the high intensity state. It is also possible that the precession of the accretion disc causes changes in the projected area of the inner disc which could result in the different intensity states. However, then the kHz QPOs would be expected in all intensity states, i.e., also during the high intensity state. The non-detections of the kHz QPOs during the high intensity state constitute a serious problem for the precessing accretion disc model.

We now compare our results obtained for the NBO to those which have been published previously. The best cases to compare our results with are the Ginga June 1987 observations of Cygnus X-2 reported by Mitsuda & Dotani (1989; see also Wijnands et al. 1997) and the EXOSAT November 14/15 1985 observations which are reported by Dieters et al. (2000; see also Hasinger et al. 1985,1987). Wijnands et al. (1997) classified the intensity state of Cyg X-2 during the Ginga observations as intermediate between the high and the medium intensity state. Kuulkers et al (1996) classified the state during the EXOSAT observations as the medium state (see Wijnands et al. 1997 for the classification of the different states of Cygnus X-2).

During the Ginga observations (Mitsuda & Dotani 1989; Wijnands et al. 1997) and the EXOSAT observations (Dieters et al. 2000), the NBO was best visible on the middle part of the NB. The NBO was much less clear on the upper and lower parts of the NB (although excess power on top of a power-law function was present near 5 Hz). The exact position on the Z track where the NBO is present during those observations is very similar to the positions on the track where we detect the NBO. Although we detected the NBO up to the upper most part of our observed NB and Mitsuda & Dotani (1989) and Dieters et al. (2000) did not detect the NBO on the upper part of their observed NB, this difference is most likely due to the fact that we did not observe the NB all the way to the HB/NB vertex (see Fig. 3b) and most likely they did. If we would have observed the NB all the way to the HB/NB vertex most likely we would not have detected the NBO either on this part of the Z track. Wijnands et al. (1997) also reported that the NBO was best visible at the middle of the NB during the medium state and that it tends to get broader when the source moves further down the NB. The energy dependence and time lags of our NBO are also nearly identical to those obtained for the NBO by Mitsuda & Dotani (1989) and Dieters et al. (2000). This indicates that the physical mechanism behind the long-term variation of the X-ray flux of Cygnus X-2 does not significantly affects the properties of the NBO.

From the presently available results, we conclude that so far no significant difference is present in the NBO properties during the high intensity state of Cygnus X-2 and during its intermediate and medium intensity state. Recently, Kuulkers, Wijnands, & van der Klis (1999) reported on RXTE observations of Cyg X-2 when this source was in a very low intensity state. No NBOs were detected. However, no complete Z track was traced out during these observations but most likely only the lower part of the NB and part of the FB. It is possible that when a complete Z track would have been traced out, the NBO would have been present at the middle part or upper part of the NB.

To explain the NBO in Z sources, it has been proposed that, when those sources are accreting close to the Eddington accretion limit, part of the accretion will occur in an approximately spherical radial inflow. In this radial flow a radiation pressure feed-back loop can be set up that causes 5–7 Hz oscillations in the optical depth of the flow (Fortner et al. 1989; Lamb 1989; Miller & Lamb 1992). This radiation-hydrodynamic model is able to explain the combination of the minimum in the rms amplitude spectrum of the NBO and the $``$140–150° phase shift in the phase lags at the same photon energy (Miller & Lamb 1992). Miller & Lamb (1992) used the results reported by Mitsuda & Dotani (1989) to constrain their model. Because our results are so similar to those published by Mitsuda & Dotani (1989), we are unable to constrain this radiation-hydrodynamic model further. In order to constrain this model, data are needed with a higher energy resolution than our RXTE/PCA data and/or data have to be obtained for the energy range below 3 keV, a region which we are unable to probe with our RXTE/PCA data. However, it has already been showed that this model cannot easily explain the differences in the NBO properties observed for Cyg X-2 and those observed for GX 5–1 and Sco X-1 (see, e.g., Dieters et al. 2000).

The recent discovery of a QPO near 7 Hz in the atoll source 4U 1820–30 (Wijnands, van der Klis, & Rijkhorst 1999) might even impose a more serious problem to all NBO models. The properties of this QPO in 4U 1820–30 (i.e., its frequency and its presence only at the highest observed mass accretion rates in this source) are similar to the NBO in the Z sources. However, if indeed these QPOs are the same phenomenon, then the models explaining the NBO in Z sources which require near-Eddington mass accretion rate, will not hold because the highest mass accretion rate observed in 4U 1820–30 (i.e., the accretion rate at the time of the 7 Hz QPO) is significantly lower than the critical Eddington mass accretion rate. This is not only true for the radiation-hydrodynamic model of Fortner et al. (1989), but also for alternative models, such as the sound wave model of Alpar et al. (1992). Clearly, the NBO models have to be significantly adjusted in order to explain all NBO properties and to explain the similarities between the NBO and the 7 Hz QPO in 4U 1820–30.

## Acknowledgments

This work was supported in part by the Netherlands Foundation for Research in Astronomy (ASTRON) grant 781-76-017, by the Netherlands Research School for Astronomy (NOVA), and the NWO Spinoza grant 08-0 to E. P. J. van den Heuvel. Support for this work was provided by the NASA through the Chandra Postdoctoral Fellowship grant number PF9-10010 awarded by the Chandra X-ray Center, which is operated by the Smitsonian Astrophysical Observatory for NASA under contract NAS8-39073. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.


format: text

text:
# The Physics of Blazar Optical Emission Regions I: Alignment of Optical Polarization and the VLBI Jet

## 1. Introduction

Polarization observations have long been a very important probe of the internal structure of blazar jets. Bright spatially resolved knots often show radio polarization (E vector) aligned with the projected jet direction, indicating a perpendicular magnetic field. This suggests that shocks are responsible for compressing the jet magnetic field and accelerating the synchrotron-emitting electrons (e.g. Aller, Aller & Hughes 1985).

The relation between optical polarization and VLBI structure provides a unique tool for investigating the regions of jet formation on $`<`$ parsec scales. While previous statistical investigations have shown a tendency for optical polarization to be aligned with the jet (Impey et al. 1991, Rusk & Seaquist 1985), the interpretation of optical polarization is less clear because the emitting regions in blazars are not resolved, blazars often show violent short-term optical polarization variability, and the old radio observations did not have sufficient angular resolution to probe the region near the optically emitting core. A few quasi-simultaneous optical-VLBI observations indicate that the optical polarization is aligned with the direction of newly-ejected blobs at the highest VLBI resolutions (Gabuzda & Sitko 1996, Lister & Smith 2000). The optical polarization may originate in shocks at the base of the jet.

We re-address the question of optical alignment, taking advantage of a more extensive optical polarization database, and more and improved VLBI maps.

## 2. Data and Derived Parameters

Our sample consists of 31 BL Lac type objects (BLLs) and 20 high polarization QSOs (simply called QSOs) with extensive optical linear polarization data and high quality VLBI maps taken from the literature. Optical polarization data from observations spanning 20–30 years were collected from the literature and McDonald Observatory archives.

We determined the following parameters for each blazar:

1. The probability that the measured optical polarization vectors are drawn from an isotropic distribution.
2. The preferred optical polarization position angle. This is the angle of the vector average of the unit vectors corresponding to each polarization measurement. We calculate this for objects with an isotropic distribution probability $`<`$ 15%. In these cases, our data are consistent with a single preferred angle (for an exception with two preferred angles, see the paper by Cross & Wills, these Proceedings).
3. The position angle for the VLBI inner structure. While some blazar jets are straight, many are curved even on very small scales (Gomez et al. 1999, Kellermann et al. 1998). Therefore, we measure position angles on both 1–3mas and 5–20mas scales, and determine a jet bending angle (the difference between them).

## 3. Results

1. Most BLLs show long-term, preferred optical-polarization angles despite their violent short-term variability (Figure 1–left). The probability of this distribution arising by chance is $`<<10^4`$ for BLLs alone, and for BLLs and QSOs combined. The QSOs’ distribution is significantly different from the BLLs’ (0.5% chance for the two to arise from the same underlying distribution) and consistent with an isotropic angle distribution.
2. When we look only at the objects with preferred optical polarization position angles most BLLs have preferred optical polarization aligned with the VLBI 1–3mas jet. For BLLs, or BLLs and QSOs combined, the probability that Figure 1–right represents an isotropic distribution of angles is $`<`$0.1%.
3. The objects with preferred optical polarization angles show a strong tendency to have straight VLBI jets (bending angle $`<15^{}`$) compared with objects with no preferred optical polarization angles (Figure 2). The probability for the objects with preferred optical angles to have the same VLBI bending angle distribution as the ones with no preferred angles, is less than 1%. Objects with preferred optical polarization angle and small VLBI bending are mostly BLLs.

## 4. Discussion

A natural explanation for the result that optical polarization tends to align with the jet, is that the optical synchrotron emission arises from a shock front in which the jet magnetic field has been compressed, on average, perpendicular to the jet. The large scatter in the optical polarization angles for a given blazar suggests that, in the inner jet region, the compressed magnetic field changes direction with time. Possible explanations are that the inner jet is internally unstable, or shocks may form via interaction with gas surrounding the central engine. The jets of QSOs may be affected by gas in the NLR and BLR, present in QSOs but absent in BLLs. This may explain why QSOs show preferred optical polarization angles less frequently. The variations may be enhanced by the effects of foreshortening and relativistic beaming.

The tendency that objects with well-defined preferred optical polarization directions also have very small VLBI scale jet bending indicates that a well-behaved straight jet on parsec to Kpc scales corresponds to a well-behaved jet on sub-parsec (optical) scales. Large curvature is likely to be the effect of projection of small jet curvature at very small viewing angles (Gower et al. 1982). Possible causes of jet curvature are  an interaction with the environment, or  an apparent curvature. In the first case, how does the base of the jet know about the environment on much larger scales? The angular resolution of optical observations is at best a factor of 100 worse than VLBI, often $`>100`$mas. So we do not have direct evidence to test our assumption that the more energetic optical photons arise near the base of the jet. The optical emission could arise in the same shocks giving rise to cm-wavelength emission. The observation of rapid polarization variation at cm wavelengths, outside the core, gives credence to this idea (Gabuzda et al. 2000). In the second case, the direction of particle ejection may vary with time, for example, via a precession jet (e.g. Hummel et al. 1997). Present data are inadequate to address changes in optical polarization position angles on precession time scales.

## References

Aller, H. D., Aller, M. F. & Hughes, P. A. 1985, ApJ, 298, 296

Gabuzda, D. C., Kochenov, P. et al. 2000, MNRAS, 313, 627

Gabuzda, D. C., Sitko, M. L. & Smith, P. S. 1996, AJ, 112, 1877

Gomez, J., Marscher, A. P. et al. 1999, ApJ, 519, 642

Gower, A. C., Gregory, P. C., Unruh, W. G. & Hutchings, J. B. 1982, ApJ, 262, 478

Hummel, C. A., Krichbaum, T. P. et al. 1997, A&A, 324 857

Impey, C. D., Lawrence, C. R. & Tapia, S. 1991, ApJ, 375, 46

Lister, M. L. & Smith, P. S. 2000, ApJ, 541, 66

Kellermann, K. I., Vermeulen, R. C. et al., 1998, AJ, 115, 1295

Rusk, R. & Seaquist, E. R. 1985, AJ, 90, 30


text:
# Non-Equilibrium Duality

\[

## Abstract

The dual ’worldline’, or ’string’, description of adiabatic Ginzburg-Landau field theory near a phase transition in terms of quasi-Brownian strings and loops is well understood. In reality, the implementation of a transition is intrinsically non-equilibrium. We sketch how time-dependent Ginzburg-Landau theory leads to a modified dual string picture in which a causal bound on the growth of unstable string prevents the uncontrolled proliferation of string (the Shockley-Hagedorn transition) suggested by the adiabatic approximation.

preprint: cond-mat/??? Imperial/TP/??
\]

From the early days of lattice models it has been appreciated that time-independent Ginzburg-Landau (TIGL) theory has a dual ’worldline’ representation in terms of near-Brownian strings and loops. More details are given in Refs. .

The main ingredient in all these calculations is that they are performed for systems in equilibrium, in which the temperature $`T`$ is fixed at values close to the critical temperature $`T_c`$. For the case that interests us, that of continuous transitions, the instabilities that characterise the transitions are, in this adiabatic dual picture, a consequence of the uncontrolled proliferation of string.

In practice, transitions occur in a finite, often short, time in which the temperature $`T`$ crosses $`T_c`$ from its initial value. The adiabatic approximation is only valid away from the transition, when the field can keep in step with the changing environment. Close to the transition the adiabatic regime is replaced by an impulse regime, in which the field falls out of step with what its equilibrium behaviour would be, and freezes in. Just as the correlation length, in reality, is unable to diverge because there is not enough time, so we would expect a realistic dual worldline description in terms of strings to show a limit to their proliferation as we pass through a transition.

The question that we shall begin to address here is how the intrinsically non-adiabatic transition has a dual representation in terms of strings and loops.

However, we begin with a brief recapitulation of equilibrium dual theory. The simplest condensed matter TIGL theory is that of a single complex field $`\varphi `$, with Ginzburg-Landau free energy

$`F(T)={\displaystyle d^3x\left(\frac{\mathrm{}^2}{2m}|\varphi |^2+\alpha (T)|\varphi |^2+\beta |\varphi |^4\right)}`$in which the chemical potential $`\alpha (T)=\alpha _0ϵ(T)`$, where $`ϵ(T)=(T/T_c1)`$, vanishes at $`T_c`$. Such an energy provides a reasonable description of superfluid $`{}_{}{}^{4}He`$, a simplified model for $`{}_{}{}^{3}He`$ and a good description of the scalar sector of low-$`T_c`$ superconductors.

It is convenient to work in spatial units of $`\xi _0=\sqrt{\mathrm{}^2/2m\alpha _0}`$ when, on rescaling the field, self-coupling and temperature,

$`\overline{F}(T)={\displaystyle d^3x\left(|\varphi |^2+ϵ(T)|\varphi |^2+\overline{\beta }|\varphi |^4\right)}.`$At temperature $`T_0>T_c`$ the free-field correlation function that follows from

$`\overline{F}_0(T_0)={\displaystyle d^3x\left(|\varphi |^2+ϵ(T_0)|\varphi |^2\right)}`$is

$`\varphi (𝐱)\varphi ^{}(\mathrm{𝟎})=G_0(r)={\displaystyle d/^3ke^{i𝐤.𝐱}P(k)},`$($`r=|𝐱|`$) in which the power spectrum

$`P(k)={\displaystyle \frac{1}{𝐤^2+ϵ(T_0)}}={\displaystyle _0^{\mathrm{}}}𝑑\tau e^{\tau k^2}e^{\tau ϵ(T_0)}`$has the usual representation in terms of the Schwinger proper-time $`\tau `$. In turn, this gives

$$G_0(r)=_0^{\mathrm{}}𝑑\tau \left(\frac{1}{4\pi \tau }\right)^{3/2}e^{r^2/4\tau }e^{\tau ϵ(T_0)}.$$
(1)
$`G_0(0)`$, necessary for loops, diverges from the UV singularities at $`\tau =0`$, but these can be regulated with a fixed cutoff. We choose a cutoff of order unity in units of $`\xi _0`$, which will be implicit throughout.

The dual picture is obtained by observing that $`(1/4\pi \tau )^{3/2}e^{r^2/4\tau }`$ is the probability distribution for a Brownian ’worldline’ or, more usefully, a ’polymer’, path $`𝐱(\tau )`$, parametrised by $`\tau `$, having its endpoints separated by a distance $`r`$. Specifically, it is expressible as the sum over such paths

$`\left({\displaystyle \frac{1}{4\pi \tau }}\right)^{3/2}e^{r^2/4\tau }={\displaystyle _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}}𝒟𝐱\mathrm{exp}\left[{\displaystyle _0^\tau }𝑑\tau ^{}{\displaystyle \frac{1}{4}}({\displaystyle \frac{d𝐱}{d\tau ^{}}})^2\right].`$If we think of the path as having step length unity (in units of $`\xi _0`$) then $`\tau `$ is proportional to the length of the path and we shall use $`\tau `$ and path length synonymously.

Then $`G_0(r)`$ is the sum over string paths of all lengths,

$`G_0(r)={\displaystyle _0^{\mathrm{}}}𝑑\tau {\displaystyle _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}}𝒟𝐱e^{S_{eq}[𝐱;\tau ,ϵ(T_0)]},`$where $`S_{eq}[𝐱;\tau ,ϵ(T_0)]`$ is the equilibrium Euclidean action

$`S_{eq}[𝐱;\tau ,ϵ(T_0)]={\displaystyle _0^\tau }𝑑\tau ^{}\left({\displaystyle \frac{1}{4}}({\displaystyle \frac{d𝐱}{d\tau ^{}}})^2+ϵ(T_0)\right).`$

The free-field partition function

$`𝒵={\displaystyle 𝒟\varphi 𝒟\varphi ^{}e^{\overline{F}_0}}=\mathrm{exp}[tr\mathrm{ln}(^2+ϵ(T_0))]`$is just that of a gas of free orientable polymer or string loops, whose lengths are parametrised by $`\tau `$, with one-loop partition function

$$\mathrm{ln}𝒵=_0^{\mathrm{}}\frac{d\tau }{\tau }\left(\frac{1}{4\pi \tau }\right)^{3/2}e^{\tau ϵ(T_0)}.$$
(2)
$`ϵ(T_0)`$ is understood as the energy/unit length or tension of the string. The factor $`\tau ^{5/2}`$ in (2) is the length distribution for Brownian strings of length $`\tau `$. As we drive $`T_0T_c+`$ it is easier to create string. If, formally, we take $`T_0`$ below $`T_c`$ then the onset of ’negative’ tension makes long loops overwhelmingly favoured and the condensation of loops that follows is understood by condensed matter physicists as a Shockley-Feynman transition, and by quantum field theorists as a Hagedorn transition.

The average loop length is

$$\tau =\frac{d\mathrm{ln}𝒵}{dϵ(T_0)}=\frac{_0^{\mathrm{}}𝑑\tau \left(\frac{1}{4\pi \tau }\right)^{3/2}e^{\tau ϵ(T_0)}}{_0^{\mathrm{}}\frac{d\tau }{\tau }\left(\frac{1}{4\pi \tau }\right)^{3/2}e^{\tau ϵ(T_0)}}.$$
(3)
Yet again, forcing $`T`$ below $`T_c`$ makes $`\tau `$ diverge from the IR divergence of the integrals at large $`\tau `$. A formal IR cutoff at $`\tau =\tau _{max}`$ leads to $`\tau =O(\tau _{max})`$.

The incorporation of the $`\overline{\beta }|\varphi |^4`$ field self-interaction into equilibrium dual string theory is implemented by the introduction of a repulsive steric string interaction at the points where paths or strings cross. Although this reduces the weight of string configurations, qualitatively the transition still occurs in this adiabatic picture because of the proliferation of strings.

In practice, a change of phase is enforced by changing the sign of $`ϵ(t)`$, either by reducing $`T`$, or by varying the critical temperature $`T_c`$. Experiments with superconductors and (the neutron bombardment of) $`{}_{}{}^{3}He`$ do the former, while pressure quenches of $`{}_{}{}^{4}He`$ do the latter. If we write

$`ϵ(t)=ϵ(T(t))={\displaystyle \frac{T(t)}{T_c(t)}}1,`$the adiabatic approximation for the ’free’ correlation function $`G_0(r,t)`$ at time $`t`$ is

$$G_0^{ad}(r,t)=_0^{\mathrm{}}𝑑\tau _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}𝒟𝐱e^{S_{eq}[𝐱,\tau ,ϵ(t)]},$$
(4)
in which we just make a straightforward substitution of the equilibrium $`ϵ(T_0)`$ with $`ϵ(t)`$ in $`S_{eq}[𝐱,\tau ,ϵ(T_0)]`$. This is treating the equilibrium pictures as a series of snapshots that can be run together as a continuous film. In this film string lengths increase uncontrollably as we cross the transition. In particular, $`G_0^{ad}(r,t)`$ of (1) is simply calculable as the Yukawa correlator

$`G_0^{ad}(r)={\displaystyle \frac{1}{4\pi r}}e^{r/\xi _{ad}(t)},`$where, on rescaling, $`\xi _{ad}(t)=\xi _0/\sqrt{ϵ(t)}`$ diverges as $`T(t)T_c`$. This cannot be the case, since causality alone prevents the correlation length diverging in a finite time. In terms of the dual picture this implies that the production of an infinity of string in a finite time is equally prohibited. The question is, how?

To answer this, we adopt the time-dependent Landau-Ginzburg (TDLG) equation for $`\overline{F}`$,

$$\frac{1}{\mathrm{\Gamma }}\frac{\varphi }{t}=\frac{\delta \overline{F}}{\delta \varphi }+\eta ,$$
(5)
where $`\eta `$ is Gaussian thermal noise, satisfying

$$\eta (𝐱,t)\eta ^{}(𝐲,t^{})=2T(t)\mathrm{\Gamma }\delta (𝐱𝐲)\delta (tt^{}).$$
(6)
Let us continue to consider the free-field case, $`\overline{F}=\overline{F}_0`$. In (5) the natural unit of time is $`\tau _0=1/\alpha _0\mathrm{\Gamma }`$ and, in units of $`\tau _0`$ and $`\xi _0`$, Eq.(6) becomes

$$\dot{\varphi }(𝐱,t)=[^2+ϵ(t)]\varphi (𝐱,t)+\overline{\eta }(𝐱,t).$$
(7)
where $`\overline{\eta }`$ is the renormalised noise. It is straightforward to show that the equal-time correlation function is now

$`\varphi (𝐱,t)\varphi ^{}(\mathrm{𝟎},t)=G_0(r,t)={\displaystyle d/^3ke^{i𝐤.𝐱}P(k,t)}.`$in which the power spectrum $`P(k,t)`$ has a representation in terms of the Schwinger proper-time $`\tau `$ as

$`P(k,t)={\displaystyle _0^{\mathrm{}}}𝑑\tau \overline{T}(t\tau /2)e^{\tau k^2}e^{_0^\tau 𝑑\tau ^{}ϵ(t\tau ^{}/2)},`$where $`\overline{T}`$ is the renormalised temperature. In turn, this gives

$`G_0(r,t)=`$

$$_0^{\mathrm{}}𝑑\tau \overline{T}(t\tau /2)\left(\frac{1}{4\pi \tau }\right)^{3/2}e^{r^2/4\tau }e^{_0^\tau 𝑑\tau ^{}ϵ(t\tau ^{}/2)}.$$
(8)
This differs significantly from (4), which would replace $`ϵ(t\tau ^{}/2)`$ in (8)with $`ϵ(t)`$. Further, because of $`\overline{T}(t\tau /2)`$, the paths are no longer Brownian in their length distribution.

However, $`G_0(r,t)`$ can still be expressed in terms of paths, as

$`G_0(r,t)={\displaystyle _0^{\mathrm{}}}𝑑\tau \overline{T}(t\tau /2){\displaystyle _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}}𝒟𝐱e^{S[𝐱,\tau ,t]}.`$$`S[𝐱,\tau ,t]`$ is not the equilibrium action $`S_{eq}[𝐱,\tau ,ϵ(t)]`$, but

$`S[𝐱,\tau ,t]={\displaystyle _0^\tau }𝑑\tau ^{}\left({\displaystyle \frac{1}{4}}({\displaystyle \frac{d𝐱}{d\tau ^{}}})^2+ϵ(t\tau ^{}/2)\right).`$Unlike the case for the adiabatic approximation (4), the local tension $`ϵ(t\tau /2)`$ is dependent on both the time $`t`$ at which we examine the strings and on the position $`\tau `$ along them. As a result the integrated tension $`_0^\tau 𝑑\tau ^{}ϵ(t\tau ^{}/2)`$ does not vanish at the transition time. Since it is the uniform vanishing of tension which triggers the avalanche of string production we see already that it will not happen in this case.

As a simple, if unrealistic, demonstration that transitions implemented in a finite time do not impose singular dual behaviour we consider the case of an instantaneous quench at time $`t=0`$ from a temperature $`T_0`$ above the critical temperature $`T_c`$ to absolute zero. That is, $`ϵ(t)=ϵ(T_0)\theta (t)`$.

Simple calculation shows that $`G_0(r,t)`$ takes the form of (1) for $`t<0`$, whereas for $`t>0`$

$`G_0(r,t)`$ $`=`$ $`e^{2tT_0/T_c}{\displaystyle _{2t}^{\mathrm{}}}𝑑\tau {\displaystyle _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}}𝒟𝐱e^{S_{eq}[𝐱,\tau ,ϵ(T_0)]}`$ (9)
$`=`$ $`e^{2tT_0/T_c}{\displaystyle _{2t}^{\mathrm{}}}𝑑\tau \left({\displaystyle \frac{1}{4\pi \tau }}\right)^{3/2}e^{r^2/4\tau }e^{\tau ϵ(T_0)}.`$ (10)
After the transition we have a representation in terms of positive tension paths at the initial temperature $`T_0`$. The unstable (negative tension) paths, for which $`\tau <2t`$, are totally excluded. Instead, the instabilities are encoded in the non-singular exponential prefactor.

This lack of IR singular behaviour is made even more explicit in a saddle-point approximation for $`G_0(r,t)`$. For $`r^2>4ϵ(T_0)t`$ we find

$`G_0(r,t)e^{2tT_0/T_c}{\displaystyle \frac{1}{4\pi r}}e^{r/\xi },`$where $`\xi =\xi (T_0)`$ remains frozen in at its initial equilibrium value. Unlike $`\xi _{ad}(t)`$ there is no divergence of $`\xi `$ as we cross $`T=T_c`$ in this abrupt way.

For the equilibrium theory, $`\tau `$ of (3) can also be expressed as

$`\tau ={\displaystyle \frac{G_0(0)}{2G_0^{\prime \prime }(0)}},`$where the prime denotes differentiation with respect to $`r`$. If we adopt the same definition out of equilibrium the benign effect of the quench to $`T_f=0`$ is again apparent in that $`\tau _t`$ agrees with (3) for $`t<0`$, but is

$$\tau _t=\frac{_{2t}^{\mathrm{}}𝑑\tau \left(\frac{1}{4\pi \tau }\right)^{3/2}e^{\tau ϵ(T_0)}}{_{2t}^{\mathrm{}}\frac{d\tau }{\tau }\left(\frac{1}{4\pi \tau }\right)^{3/2}e^{\tau ϵ(T_0)}}$$
(11)
for $`t>0`$. The exponential prefactors have cancelled to reproduce the equilibrium result (3), again at the initial temperature $`T=T_0`$, but for the absence of unstable loops with length less than $`2t`$, to which there is now no reference. Because of the exponential damping of long string the dominant string length is $`\tau =2t`$.

There is no question of an IR divergence of loop length as naively suggested by the equilibrium theory. We understand the stability of loops with length $`\tau >2t`$ in (11) as a causal bound. In our units, negative tension can only propagate at speed $`c=1`$, which happens to be the cold speed of sound in the $`\varphi `$ field.

This is reinforced by an extension to non-zero final temperature $`T_f`$

$`ϵ(t)`$ $`=`$ $`ϵ(T_0)>0\text{if }t<0\text{,}`$ (12)
$`=`$ $`ϵ(T_f)<0\text{if }t>0.`$ (13)
For $`t>0`$ we now find

$`G_0(r,t)`$ $`=`$ $`e^{2t(T_0T_f)/T_c}{\displaystyle _{2t}^{\mathrm{}}}𝑑\tau {\displaystyle _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}}𝒟𝐱e^{S_{eq}[𝐱,\tau ,ϵ(T_0)]}`$ (14)
$`+`$ $`{\displaystyle \frac{T_f}{T_0}}{\displaystyle _0^{2t}}𝑑\tau {\displaystyle _{𝐱(0)=\mathrm{𝟎}}^{𝐱(\tau )=𝐱}}𝒟𝐱e^{S_{eq}[𝐱,\tau ,ϵ(T_f)]}.`$ (15)
We see explicitly how the limited length $`\tau <2t`$ of string with negative tension $`ϵ(T_f)<0`$ prevents such string from giving a divergent contribution. Small unstable loops are no longer precluded, but their contribution is finite. The dominant length remains $`\tau =2t`$.

Although an instantaneous quench is impossible, more general quenches show similar qualitative behaviour. Specifically, suppose that $`ϵ(t)`$ decreases monotonically, with a single zero $`ϵ(0)=0`$ at time $`t=0`$. As before, only strings of limited length $`\tau <2t`$, for which $`ϵ(t\tau /2)<0`$, have negative tension. However, for $`\tau >2t`$, for which $`ϵ(t\tau /2)>0`$, strings have segments with negative and positive tension. There is a causal bound $`c=1`$ on the speed at which instability can propagate along a string. The string with negative tension, with a contribution that is independent of $`\tau `$, gives a prefactor growing at least exponentially. The effect is to leave only stable string of length $`\tau >2t`$ in the $`\tau `$ integral for $`\tau >2t`$. This boundedness on unstable string prevents the unlimited production of string suggested by the adiabatic approximation.

The exponential growth of $`G_0(r,t)`$ with time in (10), and more generally, can only be accommodated for a very short period, since $`|\varphi |^2=G_0(0,t)`$ must be constrained by the value of the order parameter $`|\varphi |`$ after the transition, equal to $`\sqrt{\overline{\beta }^1/2}`$ in the absence of corrections. A rough guide to the maximum time $`\overline{t}`$ for which the free-field approximation is valid is that $`G_0(0,\overline{t})=\overline{\beta }^1/2`$. As $`t`$ approaches $`\overline{t}`$, and thereafter, the reaction of the field with itself will cut off the exponential growth.

It is difficult to see how the dual string picture will survive at later times without further approximation. One indication is through a mean-field (or large-N) approach in which the linear nature of the theory is maintained self-consistently. In this approach (7) is replaced by

$$\dot{\varphi }(𝐱,t)=[^2+ϵ_{eff}(t)]\varphi (𝐱,t)+\overline{\eta }(𝐱,t).$$
(16)
with

$`ϵ_{eff}(t)`$ $`=`$ $`ϵ(t)+p\overline{\beta }|\varphi |^2`$ (17)
$`=`$ $`ϵ(t)+p\overline{\beta }G(0,t),`$ (18)
where $`G(0,t)`$ is determined self-consistently from (16). The coefficient $`p`$ depends on whether we adopt a mean-field (Hartree) approximation or a large-N limit for $`N=2`$.

The assumed single zero of $`ϵ(t)`$ at $`t=0`$ will lead to a zero of $`ϵ_{eff}(t)`$ at $`t0`$, and the previous analysis applies; with only segments no longer than $`2t`$ with negative tension, there is no singular behaviour as we cross the transition.

To have a quantitative estimate of the effect of back-reaction we make a Gaussian approximation for $`G(0,t)`$ by expanding about the zero of $`ϵ_{eff}(t)`$. Assuming that the quench is not too rapid, we find that, for $`t\overline{t}`$, the average loop length is

$`\tau _t{\displaystyle \frac{_0^{\mathrm{}}𝑑\tau \frac{\overline{T}(t\tau /2)}{(4\pi \tau )^{3/2}}e^{(\tau 2t)^2|ϵ^{}(0)|/4}}{_0^{\mathrm{}}\frac{d\tau }{\tau }\frac{\overline{T}(t\tau /2)}{(4\pi \tau )^{3/2}}e^{(\tau 2t)^2|ϵ^{}(0)|/4}}},`$As before, the prefactors encoding the unstable strings cancel approximately, and we need no information about the self-consistent mass, but for the fact that $`ϵ_{eff}(t)0`$ at large times to stop $`G(0,t)`$ from growing. We see the peak at $`\tau =2t`$ in the length distribution moving clear of the UV endpoint behaviour. This continues for later times and $`\tau _t=O(t)`$ once the temperature has become low enough that $`\overline{T}(t)`$ has suppressed the UV singular behaviour. Details are given in our earlier work.

We stress that this discussion has been restricted to the ’first quantised’ dual representation of the Ginzburg-Landau theory. However, it is a familiar result from a different viewpoint. For a linear system like (16) it can be shown that the mean loop length $`\tau _t`$ satisfies

$$\tau _t=\frac{1}{4\pi n(t)},$$
(19)
where $`n(t)`$ is the density of line zeroes of the complex field $`\varphi `$. The relevance of this is that, at later times, the global vortices of this $`U(1)`$ scalar theory can be identified by the line zeroes of their cores. The linear growth of dual loop lengths with time corresponds to a $`t^{1/2}`$ behaviour for $`\varphi `$-field line-zero separation and, when vortices are well-defined, vortex separation. This scaling behaviour can be justified for the decay of vortices of $`{}_{}{}^{4}He`$, and is used to determine initial vortex densities at the $`{}_{}{}^{4}He`$ transition.

This work is the result of a network supported by the European Science Foundation. I thank Adriaan Schakel for helpful comments.


text:
# Progress using generalized lattice Dirac operators to parametrize the Fixed-Point QCD action

## 1 Motivation for improved fermionic actions

There are many reasons to seek an improved lattice Dirac operator. The most commonly used one, the Wilson Dirac operator, has a number of problems. Chiral symmetry, an essential element of QCD, is broken explicitly by a term which removes the fermion-doublers and so the pion is not massless for zero quark mass. The pion can be tuned massless by varying the hopping parameter. However, this does not remove other problems, such as the mixing of operators in different chiral representations, which makes the extraction of chiral matrix elements difficult. There are also problems measuring chiral quantities such as the quark condensate and zero modes of the Dirac operator and connecting these to the topology of the gauge fields, due to, among other things, topological defects and exceptional configurations. In addition, the Wilson Dirac operator has large lattice artifacts.

A solution to problems of chirality has been proposed by Ginsparg and Wilson : a lattice Dirac operator which satisfies

$$\{\gamma _5,D^1\}=2aR\gamma _5,$$
(1)
where $`R`$ is any local function which commutes with $`\gamma _5`$ and $`a`$ is the lattice spacing, was expected to define a lattice theory with the same chiral properties as the continuum theory. Unfortunately, no solution to this relation was found and the idea was abandoned. The first lattice regularization of fermions with chiral symmetry, the Domain Wall fermions  and the related Overlap construction , followed a different path and seemed to be unrelated to the Ginsparg-Wilson (GW) relation.

Following the observation that the Fixed-Point Dirac operator satisfies the GW relation , the interest turned to this general formulation again. The GW relation turned out to be a powerful theoretical tool: it implies an index theorem on the lattice  and all the chiral properties of the formal continuum theory . This last feature became trivial after it was recognized that an exact chiral symmetry transformation follows from eq.(1) . The observation that the Overlap Dirac operator satisfies the GW relation connected the Domain Wall approach to this general formulation . The progress in vector chiral symmetry was followed by a breakthrough in the formulation of regularized chiral gauge theories also .

All these developments have led to much recent activity on the theoretical understanding of lattice chiral symmetry  and on constructing improved lattice Dirac operators which approximately satisfy the GW relation . At the same time, detailed tests and the first large-scale simulations with the simplest version of Domain Wall and Overlap Dirac operators demonstrated the power and difficulties of working with chirally symmetric QCD actions .

## 2 General lattice Dirac operators

Before we explain our approach to the parametrization of the Fixed-Point Dirac operator we briefly discuss the structure of general Dirac operators which satisfy the basic symmetries on the lattice. This will make clear what kind of lattice operators can be used in the parametrization of the Fixed-Point Dirac operator or in any other parametrization of a Dirac operator. A more detailed discussion of these issues is given in . A general gauge covariant lattice operator with color, space and Dirac indices can be written as:

$$D=\underset{A}{}\mathrm{\Gamma }_A\underset{l}{}c(\mathrm{\Gamma }_A,l)\widehat{U}(l),$$
(2)
where the $`\mathrm{\Gamma }_A`$ are elements of the Clifford algebra basis which we define as $`\mathrm{\Gamma }=1,\gamma _\mu ,i\sigma _{\mu \nu },\gamma _5,\gamma _\mu \gamma _5`$, $`l`$ denotes a path, $`\widehat{U}(l)`$ the product of the link matrices along this path and $`c(\mathrm{\Gamma }_A,l)`$ is the coupling for the given path and Clifford algebra element.

The basic symmetries of the Dirac operator, which are gauge symmetry, $`\gamma _5`$-hermiticity, charge conjugation, hypercubic rotations and reflections, impose the following restrictions on eq.(2).

* Translation invariance requires that the Dirac operator $`D(n,n+r)`$ connecting lattice sites $`n`$ and $`n+r`$ depends on $`n`$ only through the $`n`$-dependence of the gauge field, i.e. the couplings $`c(\mathrm{\Gamma }_A,l)`$ are constants or gauge invariant functions of gauge fields, respecting locality and invariance under the symmetry transformations.
* $`\gamma _5`$-hermiticity and charge conjugation together imply that the couplings $`c(\mathrm{\Gamma }_A,l)`$ for our choice of the Clifford algebra basis are real. Further, from hermiticity it follows that the path $`l`$ and the opposite path $`\overline{l}`$ (or equivalently, $`\widehat{U}(l)`$ and $`\widehat{U}(l)^{}`$) should enter in the combination

$$\mathrm{\Gamma }\left(\widehat{U}(l)+ϵ_\mathrm{\Gamma }\widehat{U}(l)^{}\right),$$
(3)
where the sign $`ϵ_\mathrm{\Gamma }`$ is defined by $`\gamma ^5\mathrm{\Gamma }^{}\gamma ^5=ϵ_\mathrm{\Gamma }\mathrm{\Gamma }`$.
* Permutations and reflections of the coordinate axes (hypercubic rotations) imply that for a given path $`l_0`$ a whole class of paths belongs to the Dirac operator. These paths are related to $`l_0`$ by all the $`16\times 24=384`$ reflections and permutations of the coordinate axes. Under such a symmetry transformation $`\alpha =1,\mathrm{},384`$ the Clifford algebra element $`\mathrm{\Gamma }_0`$ associated with $`l_0`$ generally transforms into $`\mathrm{\Gamma }^{(\alpha )}`$. Furthermore the sign of the couplings may change, whereas their absolute value remains unchanged.

A Dirac operator satisfying all the basic symmetries can be written as

$`D`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Gamma }_0,l_0}{}}c(\mathrm{\Gamma }_0,l_0){\displaystyle \underset{\alpha }{}}O^\alpha `$
$`O^\alpha `$ $`=`$ $`\mathrm{\Gamma }^{(\alpha )}\left(\widehat{U}(l^{(\alpha )})+ϵ_\mathrm{\Gamma }\widehat{U}(l^{(\alpha )})^{}\right).`$ (4)
where the sum runs over a set of reference paths defined by $`\mathrm{\Gamma }_0`$ and $`l_0`$ as well as over all the symmetry transformations $`\alpha `$ defined by the group of permutations and reflections of the coordinate axes.

It is quite natural to include the full Clifford algebra in the Dirac operator. The Wilson operator already contains $`1`$ and $`\gamma _\mu `$, the $`O(a)`$ Symanzik condition involves $`\sigma _{\mu \nu }`$ and the topological charge is proportional to $`\text{Tr}(\gamma _5D)`$, which is zero unless $`D`$ contains $`\gamma _5`$. Furthermore, for $`\gamma _5`$-hermitian Dirac operators, the Ginsparg-Wilson relation (1) can be written as

$$D+D^{}=2aD^{}RD,$$
(5)
which can only be satisfied if $`D`$ contains the complete Clifford algebra.

## 3 Efficient implementation of general Dirac operators

On first appearances one might think that it is not feasible to calculate such a general Dirac operator which has many different couplings, where every coupling might contain as many as $`768`$ paths. But one has to keep in mind that the calculation of propagators for small quark masses needs several hundreds or even thousands of conjugate gradient steps and therefore one can afford to spend some time to precalculate and store the whole Dirac operator before one starts the calculation of the propagator. On top of this there are two reasons why the building of general operators can be done in a very efficient way:

* There are a lot of paths which are invariant under certain subgroups of the reflections and permutations which reduces the number of terms significantly and in some case even leads to a cancellation of certain terms because they have relative minus signs.
* A less trivial fact is that the sum of paths for many couplings can be factorized in an efficient way, i.e. that such large sums of many paths can be written as a product of smaller sums of fewer paths.

As an example we consider a nearest neighbour coupling with $`\mathrm{\Gamma }_0=\gamma _5`$ and generating path $`l_0=[2,1,2,3,4,3,4]`$, where a path $`l`$ is a sequence of steps in the $`\pm \mu `$ direction. All the paths of this coupling can be written in the following compact way:

$$\gamma _5\underset{\mu \nu \rho \sigma }{}ϵ_{\mu \nu \rho \sigma }\left(S_{\mu \nu }P_{\rho \sigma }+P_{\rho \sigma }S_{\mu \nu }+\text{h.c.}\right),$$
(6)
where the color matrices $`S_{\mu \nu }`$ and $`P_{\rho \sigma }`$ are certain combinations of staples or plaquettes, respectively. When all the plaquettes and staples and the most frequent combinations, e.g. $`P_{\rho \sigma }`$, are precalculated then operators as in eq.(6) can be calculated very quickly. As an illustration of this we consider one of our parametrizations of the Fixed-Point Dirac operator. It has 39 couplings in total, at least one per offset on the 1-hypercube and per type of Clifford algebra element. Building this operator takes roughly 15 times as long as to multiply the operator with a vector and therefore it is a very small fraction of the time used to perform a calculation of one propagator.

The basic operation required to calculate propagators or eigenvalues is the multiplication of the Dirac operator times a vector. For a Dirac operator with hypercubic fermion offsets which contains the complete Clifford algebra, the matrix times vector multiplication is roughly 40 times more expensive than for the Wilson Dirac operator.

In , we have included a MAPLE code which, given an input path $`l`$ and any Clifford algebra element $`\mathrm{\Gamma }`$, gives all of the paths $`l^{(\alpha )}`$ and algebra elements $`\mathrm{\Gamma }^{(\alpha )}`$ generated by all permutations and reflections of the coordinate axes. This can be used by anyone who wishes to use a very general lattice Dirac operator, independently of how this operator is optimized.

## 4 Fixed-Point improved actions

The Fixed-Point method of improving lattice actions is inspired by the Renormalization Group flow of asymptotically-free theories . Consider a lattice action which contains all possible interactions. The Renormalization Group is defined by some blocking function which averages over the fields to produce a new action on a coarser lattice with fewer fields. The new action typically has different couplings from the original action and so we can imagine the blocking step as a flow in the coupling parameter space. Infinitely many blocking steps give the complete Renormalization Group (RG) trajectory. In Figure 1, we show the RG trajectory for QCD with massless quarks. The Fixed-Point (FP) has the property that the couplings are reproduced after a blocking step. For asymptotically-free theories, the Fixed-Point is on the surface of vanishing coupling $`g=0`$.

If one starts on this surface, the RG trajectory flows quickly to the Fixed-Point. If one starts close to this surface at some small coupling $`g`$ (i.e. small lattice spacing $`a`$), the RG trajectory flows quickly towards the Fixed-Point and then flows away from it. Starting on the RG trajectory with arbitrarily fine lattices with arbitrarily small lattice artifacts, one can reach any point on the trajectory by making sufficiently many blocking steps. All actions on the trajectory describe the same physics. The physical observables of the continuum quantum theory are identical to those of any lattice quantum theory on the RG trajectory, independently of the lattice coarseness. Such lattice actions are said to be quantum perfect. The Fixed-Point action is an approximation to the RG trajectory for small couplings $`g`$ and is classically perfect, i.e. it completely describes the continuum classical theory without discretization errors.

Fixed-Point actions have many desirable features. By closely approximating the RG trajectory, they are expected to have much reduced quantum lattice artifacts. They can be optimized for locality. The Fixed-Point Dirac operator satisfies the Ginsparg-Wilson relation and so has good chiral behavior. The Fixed-Point QCD action has well-defined topology and satisfies the index theorem on the lattice. Although more expensive than standard actions, studies of other models have shown that practical Fixed-Point actions can be constructed and used .

## 5 Free Fixed-Point fermions

We first consider free massless fermions. Because the fermionic action is quadratic in the fermionic fields, the Renormalization Group step for the fermionic fields amounts to Gaussian integration, which can be done exactly. The blocking step connects the Dirac operator on the fine (F) and coarse (C) lattices by

$$D_\mathrm{C}^1=\frac{1}{\kappa }+\mathrm{\Omega }D_\mathrm{F}^1\mathrm{\Omega }^{},$$
(7)
provided $`D_\mathrm{F}`$ has no zero modes, where $`\kappa `$ is an optimizable free parameter of the blocking and $`\mathrm{\Omega }`$ is the blocking function used to integrate out the fine fields to produce the coarse fields. The FP Dirac operator is reproduced under the blocking step i.e. $`D_\mathrm{C}=D_\mathrm{F}=D_{\mathrm{FP}}`$. For free fermions, this equation can be solved exactly analytically. The FP Dirac operator is local, i.e. the couplings fall off exponentially with distance, and the rate of fall off can be maximized by varying $`\kappa `$. However, $`D_{\mathrm{FP}}`$ contains infinitely many couplings. For practicality, $`D_{\mathrm{FP}}`$ is approximated with an ultra-local operator, for which each point is only coupled to its neighbours on the hypercube. In Figure 2, we see the exact $`D_{\mathrm{FP}}`$ reproduces exactly the continuum dispersion relation, while the approximate FP Dirac operator has much smaller discretization errors than the Wilson operator. The unphysical branches in the approximate $`D_{\mathrm{FP}}`$ occur at large momenta and only have a small effect.

The Ginsparg-Wilson relation eq.(1) is a constraint on $`\{\gamma _5,D^1\}`$. The FP equation (7) connects the propagators $`D^1`$ on the coarse and fine lattices. Combining these equations gives

$$R_\mathrm{C}=\frac{1}{\kappa }+\mathrm{\Omega }R_\mathrm{F}\mathrm{\Omega }^{},$$
(8)
and at the Fixed-Point, $`R_\mathrm{C}=R_\mathrm{F}=R_{\mathrm{FP}}`$. For free fermions, this equation can also be solved exactly analytically. Choosing a symmetric overlapping block transformation $`\mathrm{\Omega }`$ with a scale factor 2 averaging over hypercubes , the exact $`R_{\mathrm{FP}}`$ is ultra-local and has only hypercubic couplings — no approximation for practicality is required. The block transformation $`\mathrm{\Omega }`$ determines $`R`$. With other methods to build a Dirac operator satisfying the Ginsparg-Wilson relation, for example the Overlap construction, $`R`$ is unconstrained and typically $`R=\frac{1}{2}`$.

Defining a rescaled Dirac operator $`d=\sqrt{2R}D\sqrt{2R}`$, the Ginsparg-Wilson relation can be written (with lattice spacing $`a=1`$)

$$d+d^{}=d^{}d,$$
(9)
i.e. the eigenvalues of $`d`$ lie on a circle of radius 1 centred at $`(1,0)`$. Using $`R_{\mathrm{FP}}`$ and the hypercubic approximation of $`D_{\mathrm{FP}}`$, we show the eigenvalue spectrum of $`d`$ in Figure 3. It lies almostly exactly on the circle, indicating that the hypercubic truncation has only slightly affected the Ginsparg-Wilson property.

## 6 Parametrization of the QCD FP Dirac operator

To parametrize the QCD FP Dirac operator $`D_{\mathrm{FP}}`$ for massless quarks, we use a general Dirac operator as defined in eq.(2) with all the couplings of the $`1`$-hypercube and all elements of the Clifford algebra. The gauge invariant functions we use as couplings of $`D_{\mathrm{FP}}`$ are polynomials in local fluctuations of the gauge fields. Furthermore we include the possibility to smear the gauge fields and to project them back to SU(3), i.e. we are using so called fat links. We use the new parametrization of the FP gauge action  which also makes use of fat gauge links.

The QCD FP action is also quadratic in the fermion fields and the Renormalization Group step can again be done analytically. The QCD FP equation is the generalization of the free equation (7) given by including gauge fields. In case $`D_\mathrm{F}`$ has zero modes, the QCD FP equation is most conveniently written as

$`D_\mathrm{C}(V)`$ $`=`$ $`\kappa \mathbf{\hspace{0.17em}1}`$
$`\kappa ^2\mathrm{\Omega }(U)[D_\mathrm{F}(U)+\kappa \mathrm{\Omega }^{}(U)\mathrm{\Omega }(U)]^1\mathrm{\Omega }^{}(U),`$
where $`\kappa `$ is an optimizable free parameter of the block transformation and $`U`$ and $`V`$ are the gauge fields on the fine and coarse lattice, respectively. They are related through the FP equation of the pure SU(3) gauge theory

$$S^{\mathrm{FP}}(V)=\underset{\{U\}}{\mathrm{min}}\left(S^{\mathrm{FP}}(U)+T(U,V)\right),$$
(11)
where $`S^{\mathrm{FP}}`$ is the FP action of the pure SU(3) gauge theory and $`T(U,V)`$ is the blocking kernel of the block transformation. As the fluctuations of the gauge fields on the fine lattice are much smaller than those of the coarse gauge fields, the Dirac operator used on the fine lattice has to have good chiral properties on gauge fields with small fluctuations. This however is by far easier than to have a Dirac operator with good chiral properties on gauge fields with large fluctuations. An important fact for the parametrization of the FP Dirac operator is that eq.(6) can also be given in terms of the propagators

$$D_\mathrm{C}^1(V)=\frac{\mathrm{𝟏}}{\kappa }+\mathrm{\Omega }(U)D_\mathrm{F}^1(U)\mathrm{\Omega }^{}(U),$$
(12)
as long as $`D_\mathrm{F}`$ has no zero modes. In contrast to eq.(6) the equation for the propagator gives much more weight to the small physical modes of the Dirac operator and can therefore be used to improve the properties of the small modes of the parametrized FP Dirac operator.

The parametrization is an iterative procedure. We first start at a large value of $`\beta `$, generate thermal coarse gauge configurations $`V`$ with the FP gauge action and determine the corresponding fine configurations $`U`$ via minimization as in eq.(11). Using the free FP Dirac operator on the fine configurations (which have very small fluctuations), $`D_\mathrm{C}(V)`$ and $`D_\mathrm{C}^1(V)`$ are calculated from eqs.(6) and (12). The couplings of the parametrized Dirac operator $`D_{\mathrm{par}}`$ are determined by minimizing the following $`\chi ^2`$-function:

$`\chi ^2`$ $`=`$ $`{\displaystyle D_{\mathrm{par}}vD_\mathrm{C}v^2}+`$ (13)
$`\lambda {\displaystyle D_{\mathrm{par}}^1vD_\mathrm{C}^1v^2},`$
where the sum runs over a number of random vectors $`v`$ on different configurations and $`\lambda `$ is a weighting factor. The use of vectors for the calculation of a $`\chi ^2`$ function for $`D_{\mathrm{par}}`$ is mandatory because the definition of $`D_\mathrm{C}`$ requires a matrix inversion which we can afford only for a limited number of vectors. The minimization of the $`\chi ^2`$ yields a parametrized FP Dirac operator $`D_{\mathrm{par}}(V)`$ which has good chiral properties over a larger range of gauge couplings than the initial truncated free FP Dirac operator.

The fitted parametrized operator $`D_{\mathrm{par}}`$ is now used on fine configurations $`U^{}`$, determined via minimization from coarse configurations $`V^{}`$ generated thermally at a smaller value of $`\beta `$. Minimizing the $`\chi ^2`$ in eq.(13) again gives $`D_{\mathrm{par}}(V^{})`$ which performs well on an even larger range of gauge couplings. The whole procedure is repeated until we reach $`\beta =3.0`$ which corresponds to $`\beta _{\mathrm{Wilson}}5.75`$. During this procedure we keep leading terms in the naive continuum limit fixed such that the tree level mass is zero, the $`O(a)`$ Symanzik condition is fulfilled, the dispersion relation starts with slope $`1`$ and the normalization of the topological charge is correct . Furthermore we fix the free field limit such that we recover the truncated free FP operator on the trivial gauge configuration.

In order to parametrize the operator $`R`$ in eq.(8) we proceed in a similar way as for the Dirac operator. We also use a general operator with fat links and fluctuation polynomials. The parametrization of $`R`$ is however simpler as it is trivial in Dirac space and therefore contains a smaller number of operators. In contrast to the equation for the block transformation of $`D_{\mathrm{FP}}`$ the corresponding equation for $`R`$ eq.(8) contains no inversion and therefore the $`\chi ^2`$-function which we minimize in order to find the optimal parametrization of $`R_{\mathrm{par}}`$ can be defined as

$$\chi ^2=R_{\mathrm{par}}R_\mathrm{C}^2,$$
(14)
where the norm here is the matrix norm $`A^2=_{i,j}a_{ij}^2`$.

## 7 Preliminary results

We have some preliminary results for the parametrized QCD FP Dirac operator. In Figure 4, we show the spectrum of exact eigenvalues of $`d=\sqrt{2R}D\sqrt{2R}`$ at $`\beta =3.0`$, i.e. $`\beta _{\mathrm{Wilson}}5.75`$, corresponding to a lattice spacing of $`a0.16\mathrm{fm}`$. We show all the eigenvalues on a $`4^4`$ volume (crosses) and those closest to the origin on an $`8^4`$ volume (circles). We see that the spectrum lies quite close to the Ginsparg-Wilson circle, with a small additive mass renormalization. This is evidence that the parametrization is sufficiently rich to reproduce well the chiral Fixed-Point properties.

In order to measure hadron masses, one needs to calculate the quark propagators. The quark propagators are determined by solving the equation $`X=D^1\eta `$ for some source vector $`\eta `$. We have measured the accuracy $`ϵ_k`$ in calculating the inverse, $`ϵ_k=D(D^1\eta )_k\eta `$, after $`k`$ conjugate gradient steps, i.e. how many times $`\eta `$ is multiplied by $`D`$ in calculating $`D^1\eta `$. To achieve a given accuracy, we find that the parametrized FP Dirac operator requires roughly 3 times fewer steps than the unpreconditioned Wilson operator corresponding to the same effective pion mass for the same gauge configuration. This reduces the overhead of $`D_{\mathrm{par}}`$ in production runs.

The chiral properties of a Dirac operator are improved by using the Overlap construction. Given some input Dirac operator $`D_0`$, a Dirac operator which satisfies the Ginsparg-Wilson relation is given by

$$D_{\mathrm{OV}}=1\frac{A}{\sqrt{A^{}A}},A=1D_0.$$
(15)
If the input operator already satisfies the Ginsparg-Wilson relation, then $`A^{}A=1`$ and $`D_{\mathrm{OV}}=D_0`$. The inverse square root is calculated by some expansion in $`A^{}A`$. How quickly the expansion converges is determined by the smallest and largest eigenvalues of $`A^{}A`$, which tells us how far the input operator is from satisfying the Ginsparg-Wilson relation.

In Figure 5, we plot the 10 smallest eigenvalues of $`A^{}A`$ for 30 configurations of size $`8^4`$ at $`\beta =3.4`$, i.e. $`\beta _{\mathrm{Wilson}}5.95`$, using $`d`$ (squares) and $`D_{\mathrm{Wilson}}`$ (crosses) as input. The largest eigenvalues of $`A^{}A`$ are $`1.5`$ and $`41`$ for $`d`$ and $`D_{\mathrm{Wilson}}`$ respectively. The parametrized FP operator $`d`$ approximately satisfies the Ginsparg-Wilson relation and the inverse square root can be determined accurately by, for example, a Legendre expansion with only a few terms, which would very much improve the already good chiral behavior of $`d`$. The expansion converges faster by projecting out the smallest eigenvalues, which are treated exactly. Using $`d`$, we see $`A^{}A`$ has few small eigenvalues and so the projection can be done quite cheaply. In contrast, using $`D_{\mathrm{Wilson}}`$ to construct $`A`$, most of the eigenvalues of $`A^{}A`$ are far from $`1`$ and very expensive expansions with up to hundreds of terms are required to determine the inverse square root . There are also more very small eigenvalues to be projected out to improve the convergence, which would also be quite expensive.

In Figure 6, we plot a histogram of the number of $`A^{}A`$ multiplications required to find the 10 smallest eigenvalues of $`A^{}A`$ using $`d`$ and $`D_{\mathrm{Wilson}}`$ as input which were plotted in Figure 5. Using $`d`$ in the Overlap construction, fewer small eigenvalues need be treated exactly and these can be determined with fewer multiplications than are necessary using $`D_{\mathrm{Wilson}}`$.

## 8 Conclusions

The project of parametrizing the QCD FP Dirac operator is ongoing. Preliminary results show that the parametrization achieved so far gives a Dirac operator which satisfies the Ginsparg-Wilson relation quite well with a small additive mass renormalization. Although a matrix vector multiplication with $`D_{\mathrm{par}}`$ is roughly 40 times more expensive than with $`D_{\mathrm{Wilson}}`$, we see that calculating eigenvalues and propagators converges much faster. The chiral properties of a Dirac operator are greatly improved by using it as input for the Overlap construction. As $`D_{\mathrm{par}}`$ is already quite close to satisfying the Ginsparg-Wilson relation, the Overlap construction using $`D_{\mathrm{par}}`$ is comparable in cost to using $`D_{\mathrm{Wilson}}`$.

## 9 Acknowledgements

This work has been supported in part by the Schweizerischer Nationalfonds.


text:
# Energetic Components of Cooperative Protein Folding

## Abstract

A new lattice protein model with a four-helix bundle ground state is analyzed by a parameter-space Monte Carlo histogram technique to evaluate the effects of an extensive variety of model potentials on folding thermodynamics. Cooperative helical formation and contact energies based on a 5-letter alphabet are found to be insufficient to satisfy calorimetric and other experimental criteria for two-state folding. Such proteinlike behaviors are predicted, however, by models with polypeptide-like local conformational restrictions and environment-dependent hydrogen bonding-like interactions.

PACS Numbers: 87.15.Aa, 87.15.Cc, 87.15.He, 87.15.By

Proteins are complex systems. Insight into their behaviors has been gained by simplified models of generic proteins \[1-6\]. To serve as stepping stones towards an elucidation of the physics of real proteins, however, these models must be subjected to rigorous evaluations against experiments. Recently, we found that a number of popular lattice protein models with 2, 3, and 20 residue types and pairwise additive contact energies do not satisfy the experimental criteria for two-state thermodynamic cooperativity which, among other conditions, requires a protein’s van’t Hoff to calorimetric enthalpy ratio $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}1`$ . While certain Gō models are cooperative, they do not address the physical nature of protein interactions since their potentials are teleological  in the sense that the energetic favorability of the native conformation as a whole is presupposed.

It has since been proposed that a cooperative interplay between local conformational preferences and nonlocal interactions could give rise to proteinlike thermodynamics . To evaluate the viability of this scenario, we introduce in this Letter a 55-mer chain model with a four-helix bundle ground state (Fig. 1) that shares common features with the corresponding protein motif . The cubic-lattice helices are right-handed (called type (i) in Ref. ), as are most $`\alpha `$-helices. The model has the following energetic components. (In the analysis below all energies are dimensionless and temperature independent.)

Contact energies. A reduced 5-letter alphabet (Fig. 1) identical to the one recently optimized  is used for nearest-neighbor interactions, with energy parameters from Table III of Kolinski et al. . While these energies are proteinlike to some extent because hydrophobic groups are placed in the native core (Fig. 1), they do not represent the full interactions between amino acid residues . Here they are adopted to capture heterogeneous aspects of intraprotein interactions .

Unfavorable local conformations. Two types of local non-proteinlike bond geometries are discouraged: the initiation of a left-handed (lh) helix (Fig. 2a), and one end of a helix taking a sharp turn (st) to fold back onto itself (Fig. 2b) are assigned unfavorable ($`>0`$) energies to take into account that in real proteins left-handed $`\alpha `$-helices are sterically disfavored, and that polypeptides are stiffer than a fully flexible chain .

Environment-dependent hydrogen bonding. The favorable many-body interactions ($`_{\mathrm{Hb}}<0`$) in Fig. 2c,d are introduced to explore an idea, suggested by experiments , that the collective strength of hydrogen bonds is stronger when they are buried in the core of a protein than when they are exposed to water. Hence we focus mainly on $`b>1`$ cases below. Analogous interactions have been used before .

Cooperative helical propagation. An extra favorable energy is assigned to every two consecutive helical turns (Fig. 2e) to encourage helix elongation. Such an effect may arise from dipole-dipole interactions between amide groups in real $`\alpha `$-helices .

The total energy of a conformation from these contributions is

$$E=E_{\mathrm{contact}}+\gamma _{\mathrm{lh}}N_{\mathrm{lh}}+\gamma _{\mathrm{st}}N_{\mathrm{st}}+_{\mathrm{Hb}}N_{\mathrm{Hb}}^{(6)}+b_{\mathrm{Hb}}N_{\mathrm{Hb}}^{(8)}+_{\mathrm{Helix}}N_{\mathrm{Helix}},$$
$`(1)`$
where $`E_{\mathrm{contact}}`$ is the sum of 5-letter contact energies, $`N_{\mathrm{lh}}`$ and $`N_{\mathrm{st}}`$ are respectively the numbers of all incidences of Fig. 2a,b. In the present analysis, $`N_{\mathrm{Helix}}`$ only counts those helices (Fig. 2e) that are parts of the four helices in the ground-state conformations (monomers 1–12, 15–26, 30–41, and 44–55; see Fig. 1); and the hydrogen bonding pairs counted by $`N_{\mathrm{Hb}}^{(6)}`$ (Fig. 2c) and $`N_{\mathrm{Hb}}^{(8)}`$ (Fig. 2d) are the 36 $`(i,i+3)`$ or $`(i,i+5)`$ contacting monomer pairs in the native helices. A first-principle treatment would have assigned hydrogen bonds and helical segments in a manner that do not require knowledge of the native structure. However, progress can nonetheless be made by the approach taken here, which seeks, as a first step in the inquiry, to ascertain the consequence of presupposing local native preferences, a presupposition that is notably less dependent on a priori knowledge than the assumption of global native preference in Gō potentials.

To efficiently determine how thermodynamic properties vary with the model energetic parameters, we use a generalization  of the standard Metropolis Monte Carlo histogram technique  to eliminate the need to perform separate direct simulations for every parameter set of interest. Typical simulations are carried out with $`_{\mathrm{Hb}}`$ $`=`$ $`_{\mathrm{Helix}}=0`$ at a certain temperature $`T^{}`$, during which numbers $`P`$ of sampled conformations are binned into a multiple-dimensional array (histogram) according to $`(E^{},N_{\mathrm{Hb}}^{(6)},N_{\mathrm{Hb}}^{(8)},N_{\mathrm{Helix}})`$, where $`E^{}`$ is the energy of the conformation. Thus the density of states of the simulated system $`g(E^{},N_{\mathrm{Hb}}^{(6)},N_{\mathrm{Hb}}^{(8)},N_{\mathrm{Helix}})`$ $`=`$ $`P(E^{},N_{\mathrm{Hb}}^{(6)},N_{\mathrm{Hb}}^{(8)},N_{\mathrm{Helix}})e^{E^{}/k_BT^{}}`$, where $`k_B`$ is Boltzmann’s constant. It follows that the partition function for any $`_{\mathrm{Hb}}`$, $`_{\mathrm{Helix}}`$ at any temperature $`T`$ is given by $`g(E^{},N_{\mathrm{Hb}}^{(6)},N_{\mathrm{Hb}}^{(8)},N_{\mathrm{Helix}})e^{E/k_BT}`$, where the summation is over $`E^{}`$, $`N_{\mathrm{Hb}}^{(6)}`$, $`N_{\mathrm{Hb}}^{(8)}`$, and $`N_{\mathrm{Helix}}`$, and $`E=E^{}+_{\mathrm{Hb}}(N_{\mathrm{Hb}}^{(6)}+bN_{\mathrm{Hb}}^{(8)})+_{\mathrm{Helix}}N_{\mathrm{Helix}}`$. A similar procedure is used to study the effects of $`\gamma _{\mathrm{lh}}`$ and $`\gamma _{\mathrm{st}}`$.

Each Monte Carlo run consists of $`1.53\times 10^9`$ attempted moves, the first $`3\times 10^7`$ of which are excluded from data acquisition. To generate a multiple-dimensional histogram, a total of 20 runs at 10 different simulation $`T^{}`$s (around the transition region) are performed, with two different random initial conformations for each $`T^{}`$. Values of $`\kappa _2`$ estimated from different $`T^{}`$s agree well, with standard deviation $`4\%`$. We have conducted extensive comparisons with direct simulations to validate the method . No energy lower than that of the structure in Fig. 1 has been encountered.

As in experimental calorimetry , we characterize the thermodynamic cooperativity of a model protein by its specific heat capacity $`C_P`$ and $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$. When baseline subtraction is not applied to the $`C_P`$ function, $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ may be equated to $`\kappa _2`$ $`2T_{\mathrm{max}}\sqrt{k_BC_P(T_{\mathrm{max}})}/\mathrm{\Delta }H_{\mathrm{cal}}`$, where $`\mathrm{\Delta }H_{\mathrm{cal}}=_0^{\mathrm{}}𝑑TC_P(T)`$ is the calorimetric enthalpy, and $`C_P(T)`$ is maximum at $`T=T_{\mathrm{max}}`$ . Baseline subtractions amount to defining a multi-conformation native state and ignoring a part of the enthalpic variation in the denatured ensemble. This effectively reduces both the calorimetric enthalpy and the maximum heat capacity value, resulting in a modified enthalpy ratio $`\kappa _2^{(\mathrm{s})}`$ ($`>\kappa _2`$) .

Our main findings are summarized in Figs. 3 and 4. The apparent $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ $`=\kappa _2^{(\mathrm{s})}`$ after empirical baseline subtractions can often be close or equal to unity even when $`\kappa _2`$ is low (Fig. 4). However, as is the case for a 3-letter model, a large discrepancy between $`\kappa _2`$ and $`\kappa _2^{(\mathrm{s})}`$ is often symptomatic of non-proteinlike significant post-denaturational chain expansion at $`TT_{\mathrm{max}}`$ . Therefore, for model evaluation, proteinlike thermodynamic cooperativity requires both a small $`\kappa _2^{(\mathrm{s})}`$ $``$ $`\kappa _2`$ and $`\kappa _2^{(\mathrm{s})}1`$.

Fig. 3 compares three models by this criterion. The least cooperative is a flexible chain model with only pairwise additive contact energies. A second model with polypeptide-like local sterics has slightly enhanced cooperativity because populations of nonnative conformations with non-proteinlike local geometries that are hitherto favorable in a fully flexible chain model are reduced. However, both of these models are not proteinlike because of their significant post-denaturational chain expansions, as is evident from their thick denatured $`C_P`$ “tails” at $`TT_{\mathrm{max}}`$ (Fig. 3), which account for these models’ relatively large differences between $`\kappa _2^{(\mathrm{s})}`$ and $`\kappa _2`$ . On the other hand, the model that also incorporates environment-dependent hydrogen bonding has more proteinlike thermodynamics: Its native $`C_P`$ tail at $`TT_{\mathrm{max}}`$ is thin, with $`C_P`$ values lower on average than that of its denatured tail (Fig. 3). This implies that its native conformational diversity is limited, thus conforming better to NMR data  than a previously considered 20-letter model . As for real proteins , its average radius of gyration undergoes a sharp change around $`T_{\mathrm{max}}`$, but has no appreciable post-denaturational increase (data not shown).

Fig. 4 surveys a range of energetic parameters. Remarkably, local helical cooperativity has only a small impact on overall folding cooperativity, and proteinlike thermodynamics is possible at $`_{\mathrm{Helix}}=0`$. While $`_{\mathrm{Helix}}<0`$ stabilizes the native state, it also stabilizes denatured conformations with partially intact native helices. Therefore, its effect on calorimetric cooperativity is not substantial because it cannot widen the average enthalpy difference between native and denatured states significantly . In contrast, calorimetric cooperativity increases sharply with more negative $`_{\mathrm{Hb}}`$. For $`0.5b1.5`$, this effect is not very sensitive to $`b`$. For example, $`_{\mathrm{Hb}}=0.5`$, $`_{\mathrm{Helix}}=0`$, $`b=0.5`$ (and $`\gamma _{\mathrm{lh}}=6.0`$, $`\gamma _{\mathrm{st}}=5.0`$) lead to $`\kappa _2`$, $`\kappa _2^{(\mathrm{s})}=`$ $`0.85`$, $`0.99`$, which are only slightly lower than the $`\kappa _2`$, $`\kappa _2^{(\mathrm{s})}`$ values of $`0.90`$, $`1.0`$ for the same $`_{\mathrm{Hb}}`$ and $`_{\mathrm{Helix}}`$ in Fig. 4 for $`b=1.5`$.

In addition, we found that a lattice version of helix capping  has a slight attenuating effect on cooperativity . Consistent with experiments , in our model, tertiary interactions are essential in stablizing helices in native structures. When folded independently, the two 12-mer sequences for the native helices at 1–12 and 15–26 are much less stable ($`T_{\mathrm{max}}0.3`$), and their thermal transitions are not calorimetrically cooperative.

Figs. 3 and 4 suggest that a cooperative interplay between local conformational preferences and nonlocal contact interactions, as exemplified by the environment-dependent hydrogen bonding in the model, is a viable mechanism for proteinlike thermodynamics; and that the required cooperative effect ($`b_{\mathrm{Hb}}=0.8`$) need not be exceedingly strong relative to the pairwise contact energies (average magnitude $`=0.66`$). This observation is consistent with a previous high-coordination lattice model study , though the latter did not address the calorimetric criterion. The present approach did not consider non-native hydrogen bonding. Cooperativity in the present model would be reduced if such non-native conformations are favored. Further investigations using continuum models with polypeptide chain geometry are necessary to ascertain whether real proteins can have substantial number of non-native hydrogen bonds. While the present model should be regarded as tentative because it relies on local native information, its proteinlike features do not follow trivially from part of its interactions’ native-centric nature per se. Important physical principles have emerged from our analysis because not all native-centric interaction schemes can bring about comparable enhancements in thermodynamic cooperativity: (i) The fact that $`_{\mathrm{Helix}}<0`$ is neither necessary nor sufficient for proteinlike thermodynamics suggests that multi-body interactions that favors local native conformation irrespective of tertiary packing cannot account for calorimetric cooperativity. (ii) A Gō model for the ground-state conformations in Fig. 1, which exploits both local and nonlocal native information but is based exclusively on pairwise additive interactions, has $`\kappa _2=0.73`$ and is thus less cooperative than the model with $`\kappa _2=0.91`$ in Fig. 3. (iii) Proteinlike steric effects contribute to cooperativity. If non-proteinlike local conformations were not disfavored in the latter model (i.e., if $`\gamma _{\mathrm{lh}}=\gamma _{\mathrm{st}}=0`$), $`\kappa _2`$ would be reduced to $`0.85`$.

We have thus mapped out a general investigative strategy and established the viability of a folding scenario. It should be emphasized, however, that satisfying the requirements for thermodynamic cooperativity is clearly necessary but not sufficient for the validity of a scenario’s underlying physical mechanisms. Whether hydrogen bonding is favorable to native stability remains controversial . The present choice of $`_{\mathrm{Hb}}<0`$ is motivated by experiments . But there have been theoretical suggestions that hydrogen bonding disfavors the folded state of a protein . In the present modeling framework, that would be detrimental to calorimetric cooperativity as it corresponds to $`_{\mathrm{Hb}}>0`$, leading to $`\kappa _2`$s even smaller than the $`_{\mathrm{Hb}}=0`$ case (Fig. 3). For instance, when $`_{\mathrm{Hb}}=+0.1`$ and $`_{\mathrm{Helix}}=0`$, (and $`\gamma _{\mathrm{lh}}=6.0`$, $`\gamma _{\mathrm{st}}=5.0`$), $`\kappa _2`$, $`\kappa _2^{(\mathrm{s})}=`$ $`0.51`$, $`0.80`$ for $`b=0.5`$ and $`\kappa _2`$, $`\kappa _2^{(\mathrm{s})}=`$ $`0.45`$, $`0.70`$ for $`b=1.5`$, respectively. If that turns out to be the case, there would be added impetus to extend the present method to ascertain the role of other mechanisms such as sidechain packing  in protein calorimetric two-state cooperativity.

This work was supported by Medical Research Council of Canada grant MT-15323.

Figure Captions

Fig. 1. Ground state conformations of the model. Black beads denote nominally hydrophobic monomers. The numbers label selected sequence positions. Each of the 3 short loops at positions (13, 14), (27, 28, 29), and (42, 43) has two iso-energetic local conformations. Thus the ground state has 8 conformations, one of which is shown.

Fig. 2. Energetic components of the model. (a,b) Each contact marked by a double arrow is assigned an energy (as shown). The dotted lines in (b) depict an alternate path of the chain from monomer $`i3`$ to $`i`$ that has one instead of two unfavorable contacts. Similarly, an energy is associated with each buried hydrogen bond (c,d) and each occurrence of two consecutive turns (layers) of a right-handed lattice helix (e). In (c,d), hydrogen bonds are represented as ladders linking pairs of encircled monomers. Their burial requires occupation of at least 6 of their neighbor sites. The energy of a completely buried bond with 8 occupied neighbor sites (d) can be stronger ($`b>1`$) or weaker ($`b<1`$) than a partially buried bond with 6 or 7 occupied neighbor sites. The latter two cases have the same energy $`_{\mathrm{Hb}}`$ and are accounted for collectively by $`N_{\mathrm{Hb}}^{(6)}`$ in Eq. (1).

Fig. 3. Specific heat capacity functions. In all three cases shown, $`_{\mathrm{Helix}}=0`$. $`T_{\mathrm{max}}`$s are marked by vertical lines. ¿From left to right, the first model has only the 5-letter pairwise contact energies ($`\gamma _{\mathrm{lh}}=\gamma _{\mathrm{st}}=_{\mathrm{Hb}}=0`$). In addition to these, the second model incorporates the repulsive interactions in Fig. 2a,b, with $`\gamma _{\mathrm{lh}}=6.0`$ and $`\gamma _{\mathrm{st}}=5.0`$. The third model further adds the favorable hydrogen-bonding energies in Fig. 2c,d, with $`_{\mathrm{Hb}}=0.53`$, $`b_{\mathrm{Hb}}=0.8`$, and $`b=1.5`$. $`\kappa _2=0.55`$, $`0.62`$, and $`0.91`$, and after the plotted inclined baselines are subtracted, $`\kappa _2^{(\mathrm{s})}=0.93`$, $`0.93`$, and $`1.01`$, respectively.

Fig. 4. Calorimetric cooperativity as a function of local helical cooperativity and environment-dependent hydrogen bonding strength, with $`b=1.5`$, $`\gamma _{\mathrm{lh}}=6.0`$ and $`\gamma _{\mathrm{st}}=5.0`$. $`\kappa _2^{(\mathrm{s})}`$ and $`\kappa _2`$ are given by the upper and lower surfaces, respectively. The black dots connected by a vertical bar mark the parameters for the most cooperative case in Fig. 3. $`\kappa _2^{(\mathrm{s})}`$ are calculated using empirical baselines  constructed as tangents of the $`C_P`$ function at $`C_P^{\prime \prime }(T)/C_P^{\prime \prime }(T_{\mathrm{max}})=0.001`$, where $`C_P^{\prime \prime }d^2C_P/dT^2`$ (see Fig. 3).