text:
# The DIRECT project: Catalogs of stellar objects in nearby galaxies. II. Eastern arm and NGC 206 in M31.Tables 1-5 are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strsbg.fr (130.79.125.5) or via http://cdsweb.u-strsbg.fr/Abstract.html They are also available through anonymous ftp on cfa-ftp.harvard.edu, in pub/kstanek/DIRECT directory.

## 1 Introduction

Starting in 1996 we undertook a long term project, DIRECT, to obtain the distances to two important galaxies in the cosmological distance ladder, M31 and M33. These “direct” distances will be obtained by determining the distance of Cepheids using the Baade-Wesselink method and by measuring the distance to detached eclipsing binaries (DEBs).

As the first step of the DIRECT project we have searched for DEBs and new Cepheids in M31 and M33. In the M31 galaxy we have analyzed five $`11\mathrm{}\times 11\mathrm{}`$ fields, A$``$D and F (Kaluzny et al. pap1 (1998), pap4 (1999); Mochejska et al. pap5 (1999); Stanek et al. pap2 (1998), pap3 (1999); hereafter Papers I, IV, V, II, III). A total of 410 variables, mostly new, were found: 48 eclipsing binaries, 206 Cepheids and 156 other periodic, possible long-period or non-periodic variables. In the first paper of the series of stellar catalogs we have presented the catalog of stars in the central part of M33 (Macri et al. lmm (2001)). In this research note we present the equatorial coordinates and $`BVI`$ photometry for 26 712 stars detected in the M31 galaxy within the monitored area of 557.8 arcmin2.

## 2 Observations and data reduction

We have observed five fields, A$``$D, concentrated on the rich spiral arm, and F, containing the giant star formation region known as NGC 206 (Fig. 1). The center ($`\alpha _{2000}`$, $`\delta _{2000}`$) coordinates of the fields were: A (11$`\stackrel{}{.}`$34, 41$`\stackrel{}{.}`$73); B (11$`\stackrel{}{.}`$20, 41$`\stackrel{}{.}`$59); C (11$`\stackrel{}{.}`$10, 41$`\stackrel{}{.}`$42); D (11$`\stackrel{}{.}`$03, 41$`\stackrel{}{.}`$27); F (10$`\stackrel{}{.}`$10, 40$`\stackrel{}{.}`$72);

M31 was primarily observed in 1996 with the 1.3 m McGraw-Hill Telescope at the Michigan-Dartmouth-MIT (MDM) Observatory. We used the front-illuminated, Loral $`2048^2`$ CCD “Wilbur” (Metzger, Tonry & Luppino metz (1993)). Data for M31 were also obtained, mostly in 1997, with the 1.2 m telescope at the F. L. Whipple Observatory (FLWO), where we used “AndyCam” (Szentgyorgyi et al. sze (1999)), with a thinned, back-illuminated, AR coated Loral $`2048^2`$ pixel CCD. The pixel scale was essentially the same at both telescopes, $`0.32`$ arcsec pixel-1, giving a field of view of roughly $`11\times 11`$ arcmin2.

For the full description of the applied data reduction, calibration and astrometry procedures the reader is referred to Paper I. Here we present only a very brief summary. Stellar profile photometry was extracted using the Daophot/Allstar package (Stetson ste1 (1987)). The transformation of instrumental magnitudes to the standard system was based on 18 standard stars (Landolt lan (1992)) observed on the night of 1996 September 14/15 at MDM. The residuals in $`V`$, $`VI`$ and $`BV`$ showed no overall offsets and no dependence on color (Fig. 2 in Paper I). A comparison with the Magnier et al. (mag (1992)) photometry showed very good agreement in $`V`$ (average $`VV_{Ma92}`$=0.013 for stars with $`V<20`$) and a strong trend in $`VI`$ residuals with the $`VI`$ color (Fig. 4 in Paper I). We have obtained independent calibrations at MDM on the night of 1996 October 2/3 with the Charlotte 10242 CCD (field B) and at FLWO on the night of 1997 October 9/10 (fields C, D, F). The offsets in $`V`$ and $`VI`$, respectively, were 0.040 and 0.016 in field B, 0.012 and 0.024 in field C, -0.014 and 0.047 in field D, -0.020 and 0.057 in field F. Apart from the offsets, we did not see anything resembling the strong trend in the $`VI`$ residuals, seen in the comparison with the photometry of Magnier et al. (mag (1992)). This discrepancy certainly deserves further attention. To check the internal consistency of our calibration, we have compared the photometry in the overlapping regions between the fields. The offsets in $`V`$ and $`I`$, respectively, were 0.022 and 0.018 between fields A and B, 0.034 and 0.024 between B and C, and -0.063 and -0.040 between C and D. The offset in $`B`$ between fields C and D was 0.007.

The transformation from rectangular to equatorial coordinates was derived using stars from the list published by Magnier et al. (mag (1992)) for fields A$``$D, and the USNO-A2 catalog (Monet et al. mon (1996)) for field F.

## 3 The catalog

In Figure 2 we plot the $`V/VI`$ and $`V/BV`$ color-magnitude diagrams (CMDs) from the combined A$``$D and F field catalogs. In the $`V/VI`$ CMD (left panel), stars near $`V22`$ mag and $`VI1.8`$ represent the top of the evolved red giant population. The vertical strip of stars with $`0.6<VI<1.2`$ and $`V<20`$ are Galactic foreground stars. Stars bluer than $`VI<0.6`$ are the upper main sequence stars in M31. In the $`V/BVCMD`$ (right panel), the most prominent feature is the upper main sequence at $`BV0`$. The Galactic foreground stars are also present, between 0.4 and 1.0 in $`BV`$.

In Tables 1$``$5, we present the catalogs for fields A$``$D and F, with the equatorial coordinates and photometry in $`VI`$ (fields A and B) and $`BVI`$ (C, D, F). For each star we list its ID, $`\alpha _{2000}`$, $`\delta _{2000}`$ equatorial coordinates, standard V, I and B magnitudes with their respective errors, and the Stetson variability index JS (Stetson ste2 (1996)). The IDs, based on the equatorial coordinates, are in the format D31Jhhmmss.s+ddmmss.s. The first three correspond to $`\alpha `$, expressed in hours (hhmmss.s), the last three to $`\delta `$ in degrees (ddmmss.s).

###### Acknowledgements.

BJM and JK were supported by the Polish KBN grant 2P03D003.17. BJM was also supported by the Foundation for Polish Science stipend for young scientists and the Polish KBN grant 2P03D025.19 and JK by the NSF grant AST-9819787. DDS acknowledges support from the Alfred P. Sloan Foundation and from NSF grant No. AST-9970812.


format: text

text:
# Scaling of thermal conductivity of helium confined in pores

## I Introduction

The superfluid transition of liquid $`{}_{}{}^{4}He`$ offers a unique opportunity for testing the finite-size scaling theory of static and dynamic critical phenomena. Recently, a sophisticated experimental study was carried out in microgravity environment, the so-called confined helium experiment (CHeX). Lipa et al.chex  measured the specific heat of helium confined in a parallel plate geometry with a spectacular nanokelvin resolution, thus, providing experimental results within a few nanokelvin of $`T_\lambda `$. When the critical temperature is approached, the bulk correlation length $`\xi `$ of the fluid can become of the order of the confining length. CHeX approached so close to the lambda point that the correlation length became macroscopic in size. In this case the whole fluid acts in a correlated way and this changes the value of global properties, such as the specific heat, relative to their bulk values. In a parallel approach, Mehta and Gasparinimetha1 ; metha2  have also reported earth-bound measurements on samples with smaller plate spacing $`L`$. The size of $`L`$ used is these measurements is smaller so that the results are not significantly influenced by the change in $`T_\lambda `$ between the top and bottom of the film because of hydrostatic pressure difference which exists due to the earth’s gravitational field.

The finite-size scaling (FSS) theoryfinite  and the renormalization-group theory (RGT)rgt  were expected to describe the behavior of the system at temperature close to $`T_\lambda `$. A testable implication of this theory is that very close to the lambda point, in a confined system with a confining length of size $`H`$, a dimensionless quantity or the ratio of two quantities having the same dimensions, is only a function of the ratio of $`\xi /H`$. Therefore the values of a given observable $`O(t,H)`$, for various values of $`H`$ and of the reduced temperature $`t=|1T/T_\lambda |`$, divided by its bulk value of $`O(t,H=\mathrm{})`$ should be a dimensionless scaling function $`f(x)`$, where $`x=\xi (t)/H`$. The results of CHeX were in remarkable agreement with predictions which were available prior to the experiment based on scaling functions obtained from renormalization-group theorythe1  and those obtained by combining FSS and the results obtained from large-scale simulationssm95 .

A second equally important step toward understanding the FSS theory is to study dynamical and transport properties near a critical point. A well-suited candidate problem for this study is the thermal conductivity $`\lambda `$ of $`{}_{}{}^{4}He`$ near $`T_\lambda `$. When $`T_\lambda `$ is approached from above, the thermal conductivity of the fluid divergesahlers1 ; kerrisk . The precise behavior of bulk $`\lambda `$ as a function of $`t`$ was studied in great detail both experimentallytam ; dingus ; tam2  and theoreticallyhal .

There are several recent theoretical studies of dynamical critical phenomena and dynamical scaling. Koch, Dohm, and Staufferkoch1  presented field-theoretical and numerical studies of the validity of dynamic finite-size scaling for relaxational dynamics in cubic geometry with periodic boundary conditions above and below $`T_c`$. Quantitative agreement between theory and Monte Carlo data was obtained by them. Koch and Dohmkoch2  have provided a prediction for the dynamic finite-size scaling function for the effective diffusion constant of model $`C`$ of Hohenberg and Halperinhohenberg . Bhattacharjeebha  derived an approximate form of the scaling function for the thermal conductivity using a decoupled-mode approximation and model $`E`$. Krech and Landaukrech  calculated the transport coefficient of the out-of-plane magnetization component at the critical point, which is related to the thermal conductivity of liquid $`{}_{}{}^{4}He`$ using Monte Carlo spin dynamics simulations of the $`XY`$ model in three dimensions on a simple cubic lattice with periodic boundary conditions. They determined the critical exponent characterizing the thermal conductivity.

Accurate experimental studies have been carried out not only for dynamic bulk phenomena with improved resolution but also dynamic properties in confined geometries deeply in the critical regionlipa ; ahlers2 . Rather recently, Kahn and Ahlerskahn  measured the thermal conductivity of liquid $`{}_{}{}^{4}He`$ confined in a glass capillary array of thickness 3 $`mm`$ with holes $`2\mu m`$ in diameter. Their results show that long cylindrical samples have a transition from three-dimensional to one-dimensional behavior and there is no phase transition in the one-dimensional system. However, as measurements over a wide range of hole-diameter are required in order to test the finite-size scaling theory for transport properties, further experimental studies are plannedbest  in order to reveal dynamical exponents near the critical point and to study the finite-size scaling behavior of the thermal conductivity in such confining geometries. To avoid the limitations imposed by the earth’s gravity, this experimental effortbest  will be carried out under microgravity conditions on the Low Temperature Microgravity Facility on International Space Station.

In this paper we wish to study the thermal conductivity $`\lambda `$ of confined helium and to calculate the scaling function associated with $`\lambda `$ for a fixed geometry. Since there are already experimental resultskahn  for the scaling function of $`\lambda `$ for the pore-like geometry, in this paper we will focus our attention to this geometry because we hope to compare with the experiment. We will examine the FSS theory for the thermal conductivity of helium confined in a bar-like geometry i.e., on an $`H^2\times L`$ lattice with $`L>>H`$. This confining geometry is similar to that of Kahn and Ahlerskahn  because two of the dimensions of a pore used in their experimental studies are confining as is the case of the bar-like geometry. We will consider the limit in which our results are independent of the bar length $`L`$. We will apply periodic boundary conditions (BC) in the $`L`$ direction because these BC approach the bulk limit faster. In the other two directions which are kept finite we will apply open boundary conditions. We will use the dynamics of planar-magnet model and Monte Carlo simulation to study $`\lambda (t,H)`$. We find that $`\lambda (t,H)H^{\pi /\nu }`$ plotted as a function of $`x=tH^{\frac{1}{\nu }}`$ fall on the same curve for a wide range of values of $`H`$ and $`t`$, using the known values of $`\nu `$ and $`\pi `$. This demonstrates that finite-size is also valid for dynamical critical properties. In addition we obtain the scaling function which fits very well the experimental data of Kahn and Ahlerskahn  using the scale of temperature and the thermal conductivity scale as free parameters.

## II The Method

We will first briefly describe the model and the numerical method used and show how the thermal conductivity is computed in our model. To describe the dynamics of a superfluid, we will use the planar magnet model which is classified as model $`F`$ (or $`E`$ in the absence of an external field) by Hohenberg and Halperinhohenberg . Matsubara and Matsudamat  has proposed model $`F`$ to explain the properties of liquid $`{}_{}{}^{4}He`$. The problem of hard core bosons can be described by a lattice gas model which can be mapped to the quantum antiferromagnet in which the superfluid order parameter corresponds to $`S_xiS_y`$ while the density of the boson system corresponds to $`1/2S_z`$. In order to study equilibrium critical properties of a superfluid one uses the $`XY`$ modelsm95 ; sm  because the planar magnet model and the $`XY`$ model belong to the same universality classnho . For critical dynamics of a superfluid, however, one needs to use the full planar magnet model in which the role of the third component of the pseudospin is crucial.hohenberg

In the pseudospin notation, the planar magnet model takes the following form:

$$H=J\underset{<ij>}{}(S_i^xS_j^x+S_i^yS_j^y),$$
(1)
where the summation is over all nearest neighbors, $`\stackrel{}{S}_i=(S_i^x,S_i^y,S_i^z)`$, and $`J`$ sets the energy scale. In the usual $`XY`$ model the two-component pseudo-magnetization corresponds to the superfluid order parameter. In the planar magnet model, the third component corresponds to the particle density and it is necessary in order to study the dynamics.

In our calculations, we use a bar-like geometry, i.e. a $`H\times H\times L`$ lattice with $`LH`$. This geometry is chosen in order to mimic the pore geometry used in experimental studies. In our calculations, open boundary conditions are imposed in the $`H`$ direction, and in the $`L`$ direction we applied periodic boundary conditions.

We use a hybrid Monte Carlo procedurekrech  which consists of a combination of steps using the Metropolis update, the Cluster updatewolff , and the over-relaxation algorithmbrown . Using this hybrid algorithm, we generate approximately 3,000-10,000 uncorrelated configurations from the equilibrium canonical ensemble at a given temperature. Each configuration is evolved using the equations of motion for the planar magnet model which are given as followskrech ; hohenberg

$$\frac{d}{dt}\stackrel{}{S}_i=\frac{H}{\stackrel{}{S}_i}\times \stackrel{}{S}_i.$$
(2)
Starting from a particular initial spin configuration, we perform numerical integration of these equations of motion. Following Ref.krech2  we use a recently developed decomposition methodde  where the integration is carried out to a maximum time $`t_{max}`$ (typically of the order of $`t_{max}`$=400) with a time step $`\delta t`$=0.05. We made sure that this way we determined the real-time history of every configuration within a sufficiently long interval of time (0$`tt_{max}`$). Finally, we compute the thermal average of a time-dependent observable (such as the thermal current-current correlation function) by averaging over all the values of the observable obtained by evolving all the independent initial equilibrium configurations generated via the hybrid Monte Carlo procedure.

Compared to calculating static critical properties, the computation of dynamical properties is far more CPU time intensive and demands large computational resources. The computations described here were carried out on a dedicated massively parallel cluster of 64 CPUs which was designed by our group to achieve high performance to cost ratio.

We computed the thermal conductivity on $`H\times H\times L`$ lattices, where $`H`$=6,8,10,12,14,20 and $`L=5H`$. The thermal conductivity $`\lambda `$ of liquid $`{}_{}{}^{4}He`$ at a given temperature $`T`$ can be calculated using the dynamic current-current correlation functionkrech :

$$\lambda =\frac{1}{k_BT\chi _{zz}}\frac{2}{\pi }_0^{\mathrm{}}𝑑t\underset{i}{}<j_0^z(0)j_i^z(t)>,$$
(3)
where the out-of-plane static susceptibility

$$\chi _{zz}=<M^{z2}>/(k_BTL^3)$$
(4)
is needed for normalization. The $`z`$-component $`j_i^z`$ of the current density $`\stackrel{}{j}_i`$ associated with the lattice point $`i`$ is defined by

$$j_i^z=J(S_i^yS_{i+e_z}^xS_i^xS_{i+e_z}^y),$$
(5)
where the notation $`i+e_z`$ denotes the nearest neighbor of the lattice site $`i`$ in the $`z`$ lattice direction.

Now, we would like to examine the finite-size scaling hypothesis for the thermal resistivity $`R(t,H)=1/\lambda (t,H)`$, and to compare our results with the existing experimental resultskahn . The dependence upon $`t`$ of the bulk thermal resistivity can be described by the power law

$$R(t)=R_0t^\pi ,$$
(6)
where $`\pi `$ is a dynamic critical exponent. Using Eq. 6, the finite-size scaling expression for the thermal resistivity $`R(t,H)`$ is given by

$$R(t,H)H^{\pi /\nu }=f(tH^{1/\nu }),$$
(7)
where the function $`f(x)`$ is universal and $`\nu `$=0.6705 is the critical exponent of the correlation lengthgold .

## III Results

In this section, we calculate the thermal resistivity, we examine its scaling behavior with respect to $`H`$ and then we compare the scaling function with the experimental results. To calculate these observables with small statistical errors even with our utilization of the most recent numerical advances and with using the 64-node dedicated cluster, it requires signinificant amount of time of high-throughput computation.

Fig. 1 shows some of our results for the thermal resistivity $`R(T,H)`$ as a function of temperature $`T`$ for various lattice sizes with open boundary conditions in the $`H`$ direction. Our results for $`R(T,H)`$ for several values of $`H`$ and $`T`$ are given in Table 1 and 2. We wish to make $`L`$ large enough so that finite-size effects due to $`L`$ are smaller than our statistical errors. We have found that taking $`L5H`$ and applying periodic boundary conditions along the direction of $`L`$ introduces insignificant finite size effects due to the finite size of $`L`$ for the temperature range studied here. Since we wish to remain in the 3D critical region, it appears that keeping $`L=5H`$ introduces small finite-size effects due to $`L`$ in this region. Lowering the temperature further, when the correlation length becomes comparable to $`L`$, the value of $`R`$ is influenced by the finite size of $`L`$. In this region the system exits the 3D critical region and it behaves as a one-dimensional system.

Notice in Fig. 1 that the thermal resistivity feels strong finite-size effects due to the bar thickness $`H`$. The arrow shows the bulk transition temperature $`T_\lambda =1.5518`$ obtained from Monte Carlo simulation using the planar magnet modelnho . In bulk helium $`R(t)`$ approaches zero as the bulk transition temperature $`T_\lambda `$ is approached from above.

We wish to avoid using any adjustable parameters to obtain scaling of our results. Thus, we need to examine if our results obey scaling using the known values of the critical exponents $`\nu `$ and $`\pi `$. The value of $`\nu `$ is accurately known from theoretical and experimental studies of static critical properties and we shall use the value $`\nu =0.6705`$ as determined by Goldner and Ahlersgold . There is less agreement between theory and experiment on the actual value of the dynamical critical exponent $`\pi `$. Ahlersahlers3  used a power law fit to the data of Tam and Ahlerstam  for their “Cell $`F`$” and he found the value $`\pi =0.4397`$. However, the dynamic scaling theoryds  had predicted a divergence in $`\lambda `$ with a critical exponent given by $`\pi =\nu /20.335`$.

Fig. 2 shows a scaling plot of the thermal resistivity scaling function $`f(x)=R(t,H)t^{\pi /\nu }`$ versus the scaled reduced temperature parameter $`x=tH^{\frac{1}{\nu }}`$, where the reduced temperature is taken relative to the bulk transition temperature $`T_\lambda `$. Our Monte Carlo data collapse onto a universal curve using the value of $`\pi 0.44`$ determined by Ahlersahlers3 . In Fig. 2 we compare our universal function $`f(x)`$ with the experimental data obtained by Kahn and Ahlerskahn  represented by a solid line. In order to do this, we used two multiplicative constants as free fitting parameters, one multiplying the scale of $`x`$ axis and another the scale of $`y`$. The agreement between Monte Carlo simulation and experiment is quite satisfactory. In the past it has been demonstratedsm ; mc5  that the boundary conditions play a significant role in defining the universal function $`f(x)`$. We believe that if we use more realistic boundary conditions, such as Dirichlet boundary conditions, along the $`H`$-direction we can reduce the number of fitting parameters to only one.

Using the theoretical value of $`\pi `$, the results of our simulation also collapse on a different scaling function given in Fig. 3 . However, if we attempt to fit the scaling curve with the experimental resistivity of Kahn and Ahlers we obtained a lower quality fit than that of Fig. 2.

In summary we have calculated the thermal resistivity $`R(t,H)`$ of liquid $`{}_{}{}^{4}He`$ in a pore-like geometry (on a $`H\times H\times L`$ lattice) applying open boundary conditions in the $`H`$ direction. We have been able to demonstrate the validity of finite-size scaling theory and we obtained the thermal resistivity scaling function $`f(x)`$ using known values for the critical exponents and no adjustable parameters. In addition, the scaling function $`f(x)`$ for $`R(t,H)`$ agrees rather well with experimental data using the temperature scale and thermal resistivity scale as free parameters.

## IV acknowledgments

This work was supported by the National Aeronautics and Space Administration under grant no. NAG3-1841.


text:
# The X-ray Spectrum of the Rapid Burster using the Chandra HETGS

## 1 Introduction

The Rapid Burster (RB) is a recurrent transient low-mass X-ray binary with an X–ray burst phenomenology that makes it unique in the Galaxy. The many idiosyncrasies of the RB are overwhelming, and it is fair to say that we do not understand why the RB, and only the RB, behaves so differently from all other low mass X-ray binaries (LMXBs). For a detailed review see Lewin, Van Paradijs & Taam (1993), hereafter LVT. The RB (MXB 1730$``$335) is the only low-mass X-ray binary (LMXB) known to produce two types of X-ray bursts (Hoffman, Marshall, & Lewin, 1978). Type I bursts, which are observed from $``$40 LMXBs, are due to thermonuclear flashes on the surface of an accreting neutron star. They typically recur on time scales of hours to tens of hours (LVT). Type II bursts are due to accretion instabilities (gravitational potential energy). In the RB, their durations range from $``$2 to $``$680 sec, and thousands of short bursts can occur in a single day. The behavior of the type II bursts from the RB is like that of a relaxation oscillator: the type II burst fluence $`E`$ is roughly proportional to the time interval, $`\mathrm{\Delta }t`$, to the following burst (the “$`E`$-$`\mathrm{\Delta }t`$” relation). The time averaged luminosity of the RB is $``$1037 erg s-1 during its active periods which typically last for several weeks. The type II burst luminosities at burst maximum range from $``$5$`\times `$1037 to $``$4$`\times `$1038 erg s-1 (LVT). Henceforward, whenever we use the word “burst(s)”, we mean type II burst(s) unless we specify otherwise.

Among the LMXBs, the RB stands alone in its complex QPO behavior, none of which is understood (LVT). With the Rossi X-ray Timing Explorer (RXTE) we may have detected ($`>98`$% significance) the signature of near coherent 306.5 Hz pulsations in the average power-density spectrum of the first second of 31 Type I bursts (Fox, et al., 2001).

Though we are largely at a loss in our understanding of this system, it is important to remember that there have been occasions that, when the RB became active, it behaved like an ordinary LMXB and then transitioned to rapid bursting. Starting August 5, 1983, no type II bursts were detected but instead a strong persistent X-ray flux and type I bursts were observed for several days Kunieda, et al. (1984); Barr, et al. (1987). In the period November 6–17, 1996, and July 7–10, 1997 (during a total of 18 ksec of RXTE observations) the RB again produced only persistent emission and type I bursts at the beginning of its activity cycle.

We are reporting results from the first high resolution X-ray spectral observation of the RB, using the Chandra High Energy Transmission Grating Spectrometer (HETGS). A description of the HETGS and its performance are given by Canizares et al. (in preparation). Our objective was to search for specral emission and absorption features that might distinguish the bursts from the persistent emission. The type II bursts, due to their E-$`\mathrm{\Delta }`$t character, are believed to be emission from the surface of the NS. Discovery of any emission or absorption features might then provide a measure of the gravitational redshift at the NS photosphere.

## 2 Observations and Data Reduction

The RB was observed with the Chandra HETGS on 24-25 October 1999 (MJD = JD - 2400000.5 = 51476.2-6.6). The Advanced CCD Imaging Spectrometer (ACIS-S) detector was used for readout and was operated in the so-called “continuous clocking” (CC) mode. This mode is mentioned in the Chandra Proposer’s Observatory Guide (p. 101). In this mode, the charge is shifted row by row to the frame storage buffer at 2.84996 ms intervals so the X-ray events are timed to an accuracy of about 3 ms by locating events within the storage buffer, assuming that the event is associated with the target. The total exposure time was 29562.7 s, during which a total of 702 bursts were detected. A few samples of the light curve from zeroth order are shown in Fig.1, which shows a type I burst about 11580 s into the observation. The type I burst was fairly typical of those discovered by Hoffman, Marshall, & Lewin (1978), having a similar peak flux to that of the type II bursts but with a longer decay time constant.

The median peak burst count rate was about 45 count s-1 but the average burst count rate is about a factor of two smaller due to the nearly triangular shapes of the burst profiles. Average burst spectra were determined by selecting events where the count rate over a 1 s interval exceeded 5 count s-1, eliminating only 10% of the total counts in type II bursts while simultaneously reducing the background (see 2.2) by more than a factor of three. A “peak” burst spectrum was also derived by restricting the selection to intervals where the count rate exceeded 30 count s-1, which gave half as many counts as the average spectrum. After removing 82 s containing the type I burst, the exposure time was 4326 s in the average spectrum and 1349 s in the peak spectrum. The average count rate was 23.4 count s-1 for the average spectrum and 39.7 count s-1 in the peak spectrum. In addition to the type I burst there were 18 “mini-bursts” where the peak flux is 10-40% of the peak flux from the remaining bursts which are mostly included in the average spectrum but entirely excluded from the peak spectrum. The first seven of these mini-bursts are shown in Fig.1.

The CC mode introduces complications into the spectral processing. Imaging information along the readout direction is lost in favor of improved event timing. In the timed exposure (TE) mode, the imaging information is used to select the HEG and MEG data for spectral extraction and for the selection of background. In CC mode, only the ACIS pulse height data are used. The MEG and HEG first orders are well separated in energy when at any specific dispersion distance although there is the ambiguity that an event dispersed into 2nd order by the MEG gratings will have the same energy as an HEG photon and dispersed by the same distance. In practice, however, the efficiency of the MEG in 2nd order is a few percent of the HEG efficiency in the 1-10 keV range, so this ambiguity is not an issue.

Background selection for spectral analysis is not simple, so this instrument mode is best suited for bright targets such as the RB during bursts. We cannot define the background spatially so we instead select regions of pulse height as a function of dispersion distance with the same width, $`\mathrm{\Delta }E`$, as used for selecting the source but displaced in energy (see step 10 of the spectral processing). Due to the difficulty in processing and analyzing the CC mode data, we took great care to check our methods using data from other sources observed with this mode, such as the Crab pulsar (observation ID 170) and Cyg X-1 (observation ID 1511).

### 2.1 Imaging

As stated, there is limited spatial information; we obtain $``$1″ resolution along one dimension of the detector array. In this dimension, when the spacecraft dither and the thermal drifts are removed, the result is a very narrow image of the HETGS zeroth order. It is well fit with a Gaussian of 0.77″ FWHM, which is consistent with the superb imaging of the Chandra mirror assembly if the image is integrated along one dimension.111Although the one dimensional response function is not well documented, one may examine results posted on the web by Dan Dewey: http://space.mit.edu/HETG/technotes/zo\_1d/focus\_trends.010227.html. In order to properly remove the telescope dither, however, we required advancing the times assigned to the events by about 5 s. Other instruments on Chandra do not require this timing offset, so we conclude that the offset is an unaccounted fixed timing offset. It was also found to be required for dither correction in all other CC mode observations, including an observation of the Crab pulsar and the pulsar B0540-69 (Kaaret, et al., 2001).222Note that Kaaret, et al. (2001) misstated the sign of the correction. A constant offset of about 5 s must be subtracted from the times of events in order to match the time of the aspect system. Of five CC mode observations we have processed, the sign and magnitude of the offset is the same. We estimate that the offset time is accurate to about 0.5s.

### 2.2 Spectra

The spectral data were reduced starting from level 1 data provided by the Chandra X-ray Center (CXC) using IDL custom processing scripts; the method is quite similar to processing data from timed exposure (TE) mode but has two wrinkles that make it more difficult. This procedure was successfully verified using CC mode data from Cyg X-1 and is more thoroughly described by Marshall, et al. (in preparation) Briefly, the procedure was to: 1) restrict the event list to the nominal grade set (0,2,3,4, and 6), 2) remove the event streaks in the S4 chip (which are not related to the readout streak and appear in CC and TE mode data alike), 3) compute a projected dispersion distance ($`x_p`$, in mm) by correcting for telescope dither along the dispersion direction after applying the timing offset determined using the dither pattern of the zeroth order events along the readout direction (as described in the last section), 4) determine the location of zeroth order, $`x_{p0}`$, by fitting a Gaussian to the one dimensional profile, 5) estimate the event wavelength $`m\lambda =P(x_px_{p0})/(R\mathrm{cos}\theta )`$ (where $`P`$ is the grating period, $`m`$ is the dispersion order, $`R`$ is the Rowland distance in mm, and $`\theta `$ is the clocking angle of the MEG or HEG), 6) correct event energies for detector node-to-node gain variations, 7) select $`\pm 1`$ order events in order to minimize background using $`|E_{ACIS}m\lambda /(hc)|1<\mathrm{\Delta }`$ where $`E_{ACIS}`$ is the event energy inferred from the ACIS pulse height, and $`\mathrm{\Delta }`$ is a function of $`\lambda `$, ranging from 0.025 at 1Å to 0.15-0.20 for $`\lambda >6`$ Å, 8) eliminate events in bad columns where the counts in an histogram deviate by 5$`\sigma `$ from a 50 pixel running median,333There were many bad columns, especially on the S2 and S4 chips. The cause may be a poor on-board bias map calculation (Peter Ford, private communication). These mostly affect the dispersed spectrum below 1 keV so they have a negligible effect on our spectral results. 10) select background events by applying a $`E_{ACIS}`$ criterion similar to that applied to select source events but accepting events in the range $`\mathrm{\Delta }`$ to $`2\mathrm{\Delta }`$, 11) bin MEG (HEG) events at 0.01Å (0.005Å), 12) eliminate data affected by detector gaps, and 13) generate and apply an instrument effective area based on the pre-flight calibration data and a set of correction polynomials derived from the CC mode observation of Cyg X-1 (Marshall et al., in preparation).

The time averaged spectra of the type II bursts are shown in Figs. 2 and 3. The data have been adaptively rebinned to give a nearly constant signal/noise ratio better than 10 except in the bin at the lowest energy. A blackbody model has been fit to the data and overplotted in the upper panel. A 5% systematic error was added to the uncertainties in quadrature so that the $`\chi ^2`$ came out to 1.00. This model of the uncertainties is not precisely correct because there are correlations between the systematic error in adjacent bins. One systematic error is due to the difference between the background spectra derived from events with $`E_{ACIS}`$ larger than $`(hc)/(m\lambda )`$ compared to those events where $`E_{ACIS}`$ is smaller than $`(hc)/(m\lambda )`$. The background is less than 10% of the signal for $`1.5<E<8.0`$ keV. The systematic errors are small in this energy band so that their true dependence on energy will not have a significant effect on our results.

The burst peak spectra are well fit by a blackbody with $`kT=1.590\pm 0.016`$ keV with a radius of $`8.88\pm 0.14`$ km for a distance of 8.6 kpc, and the interstellar column density, $`N_H`$, was measured to be $`1.81\pm 0.03\times 10^{22}`$ cm-2. The fit to the burst average spectra gave similar values for the temperature and $`N_H`$, $`1.577\pm 0.012`$ keV and $`1.74\pm 0.03\times 10^{22}`$ cm-2, but the fitted radius was significantly smaller, $`6.54\pm 0.08`$ km. The uncertainties determined in the fit are indicative of the high counting statistics ($`>10^5`$ counts) but would increase as one takes systematic errors into account; we estimate that the uncertainties should be regarded as good to a factor of 2. In principle, we would be able to verify the fitted value of $`N_H`$ by measuring the ionization edges due to Mg and Si at 1.30 and 1.84 keV, respectively. With large $`N_H`$ and relatively low source flux below the spectral peak of the blackbody emission, the flux is insufficient to measure the optical depths at the edges to better than a factor of 2, so we cannot obtain a useful independent constraint on $`N_H`$.

We searched the time-averaged spectra for narrow absorption and emission features and found several candidates. All candidates were found to be spurious, however, by comparing the HEG and MEG spectra and the $`+1`$ orders to the $`1`$ orders. The burst peak spectra also showed no features. In Fig. 4, we show estimates of the 3$`\sigma `$ limits that can be placed on any emission lines whose widths (FWHMs) are comparable to the instrument resolution. The computation uses the model fitted to the time-averaged HEG and MEG data and the effective areas. We assume that candidate features are only 2 bins wide. The curves are rather smooth, except for locations of chip gaps, so one may derive limits on broader features using these curves until the scale of the feature becomes comparable to that of the instrument calibration uncertainties. Limits on absorption lines are identical when there are a lot of counts but are systematically larger. The 3 $`\sigma `$ upper limits to the equivalent widths of any features are $`<`$ 10 eV in the 1.1-7.0 keV band and as small as 1.5 eV near 1.7 keV.

## 3 Discussion

We find that the color temperature at the peak of the type II bursts is the same as the average color temperature to within the very small statistical uncertainties even though the count rate drops by a factor of two from the average of the burst peaks to the average of all type II bursts. This finding confirms a result that is well established for significantly longer bursts – that the color temperature is not observed to vary significantly during the bursts (Marshall, et al., 1979). Thus, the profile of a type II burst seems to result from variation of the size of the emission region with a maximum of about 9 km. If all emission takes place on the surface of a neutron star, then any smaller apparent radius must result from nonuniformities over the stellar surface; in the simplest geometry, the emission is assumed to emanate from a spot on the surface. During the fainter parts of normal type II bursts and near the peaks of minibursts the observed intensity drops to 10% of the peak so the apparent fraction of the surface, $`f_a`$, is about 0.1. Marshall (1982) computed $`f_a`$ of an emission region spot covering a true fraction $`f`$ of the surface for various values of the ratio of the stellar radius to its Schwarzschild radius. For a 1.4 $`M_{\mathrm{}}`$ neutron star with a 7 km radius which is uniformly radiating across its surface, the apparent radius (measured at infinity) is 11 km, so we would infer that the type II emission region covers between 5 and 20% of the surface when faint and an average of 50-85% at burst peak. The endpoints of the allowed range depend on whether the spot is centered on the nearest or farthest points of the surface. Gravitational bending of light trajectories ensures that all bursts can be detected as long as at least 5-10% of the surface radiates, for relatively compact neutron stars.

The main objective of this observation was to determine if there would be emission lines or absorption features that could be used to determine the gravitational redshift of the emission region. At the estimated temperature of the emission region, Fe will be highly ionized, so a significant fraction of the Fe will be He-like Fe XXV, for which the strongest resonance line is at 6.70 keV. If the radiating region is spherically symmetric, then the emission region is a shell at a radius of 7 km and the Fe XXV line will be gravitationally redshifted to about 4.3 keV, where our 3$`\sigma `$ upper limit to absorption features is 3 eV. The instrumental resolution is 0.01Å, which is 15 eV at this energy. Pressure broadening of the resonance line due to the Stark effect in a highly ionized, hydrogen-dominated plasma is estimated to be about 4 eV in the rest frame of the gas while fine-structure splitting is of order 21 eV (Paerels,, 1997). There is a 30% narrowing to the observed frame due to the gravitational redshift so the line FWHM is expected to be comparable to the instrumental resolution. The thermal line FWHM is smaller, $`1`$ eV. Thus, the upper limit on the column density of Fe XXV is about $`8\times 10^{15}`$ cm-2. For a simple model atmosphere dominated by electron scattering, the column density of Fe to one optical depth is of order $`X_{Fe}/\sigma _T=5\times 10^{19}`$ cm-2, where $`X_{Fe}`$ is the abundance of Fe and $`\sigma _T`$ is the Thomson cross section. The discrepancy of $`10^4`$ leads us to suggest that Compton scattering dominates the spectrum, destroying possible line features.

At the measured color temperature, Comptonization is expected to shift low energy photons to higher energies. The color temperature is systematically larger than the effective temperature of the plasma. The magnitude of this bias is not well determined but can be a factor of 1.5 (cf. LVT). A decrease of the temperature requires a compensating increase of the inferred radius of the emission region, yielding a radius of approximately 20 km. This radius is somewhat larger than observed in other burst sources but still feasible from theoretical standpoint.

This work has been supported in part under NASA grant NAG5-3239, and NASA contracts NAS8-38249, NAS8-39073, and SAO SV1-61010.


text:
# Localized Random Lasing Modes and a New Path for Observing Localization

Localization theory and laser theory were both developed in the sixties. The propagation of quantum and classical waves in disordered media has been well understood , while the laser physics has been well established  at the same time. It was always assumed that disorder was detrimental to lasing action. However, Letokhov  theoretically predicted the possibility of lasing in a random system, called “random laser”. But only after the experimental observations by Lawandy et al. , the random laser systems have been studied intensively . Since then, many experiments were carried out that showed a drastic spectral narrowing  and a narrowing of the coherent backscattering peak . Recently, new experiments  showed random laser action with sharp lasing peaks. To fully explain such an unusual behavior of stimulated emission in random systems, many theoretical models were constructed. John and Pang  studied the random lasing system by combining the electron number equations of energy levels with the diffusion equation. Berger et al.  obtained the spectral and spatial evolution of emission from random lasers by using a Monte Carlo simulation. Very recently, Jiang and Soukoulis , by combining a FDTD method with laser theory , were able to study the interplay of localization and amplification. They obtained  the field pattern and the spectral peaks of localized lasing modes inside the system. They were able to explain  the multi-peaks and the non-isotropic properties in the emission spectra, seen experimentally . Finally, the mode repulsion property which gives a saturation of the number of lasing modes in a given random laser system was predicted. This prediction was checked experimentally by Cao et al. .

One very interesting question that has not been answered by previous studies is what is the form of the wavefunction in a random laser system. How does the wavefunction in a random system change as one introduces gain? Does the wavefunction retain its shape in the presence of gain? Another very interesting point is if we can predict a priori, which mode will lase first? What will be its emission wavelength? If we understand these issues, we will be able to design random lasers with the desired emission wavelengths. In addition, we will be able to use the random laser as a new tool to study the localization properties of random systems.

In this paper, we explore the evolution of the wavefunction without and with gain, by the time-independent transfer matrix theory , as well as by time-dependent theory . The emission spectra can be also obtained. In addition, we can also use the semi-classical theory of lasing  to obtain analytical results for the threshold of lasing, as well as which mode will lase first. This depends on the gain profile, as well as the quality factor $`Q`$ of the modes before gain is introduced.

Our system is essentially a one-dimensional simplification of the real experiments . It consists of many dielectric layers of real dielectric constant $`\epsilon _2=2.56\epsilon _0`$ ($`\epsilon _0`$ is the dielectric permeability of free space) of fixed thickness ($`b_0=100`$ nm), sandwiched between two surfaces, with the spacing between the dielectric layer assumed to be a random variable $`a_n=a_0(1+W)`$, where $`a_0=200`$ nm and $`W`$ has a random value in the range of \[-0.7, 0.7\]. We choose a 30-cell random system, as the first system of our numerical study. In Fig. 1a, we present the results of the logarithm of the transmission coefficient as a function of frequency $`f`$. These results were obtained by using the transfer matrix techniques introduced in Ref. 15. Notice that we have three typical resonance peaks (denoted $`P_1`$, $`P_2`$ and $`P_3`$) in the frequency range of 600 to 660 THz. As one can see from Fig. 1a, the linewidths of the three modes are different. $`P_3`$ has the smallest linewidth and therefore the largest $`Q`$, while $`P_1`$ has the largest linewidth and therefore the smallest $`Q`$. We have also numerically calculated the wavefunctions corresponding to these three peaks and indeed find out that the more localized wavefunction is the one with the larger $`Q`$. All the results above were obtained for the case without gain.

According to the semiclassical theory of laser physics , we generally use a polarization due to gain $`P_{gain}=\epsilon _0\chi (\omega )E=\epsilon _0(\chi ^{}(\omega )+i\chi ^{\prime \prime }(\omega ))E`$ to introduce amplifying medium effects. Both $`\chi ^{}(\omega )`$ and $`\chi ^{\prime \prime }(\omega )`$ are proportional to the outside pumping rate $`P_r`$ and can be expressed by the parameters of the gain material .

To determine which peak will lase first, we can again use the time-independent transfer matrix method (see Eq. (4) of Ref. 15) with a frequency-$`independent`$ gain, which means the width of the gain profile is very large. It is well understood that time-independent theory  for random lasers can be used to obtain the threshold for lasing. At threshold, the transmission coefficient $`T`$ goes to infinity. In Fig. 1b, we plot the $`log_{10}(T)`$ versus $`\chi ^{\prime \prime }`$ for the three peaks $`P_1`$, $`P_2`$ and $`P_3`$. In Fig. 1b, the $`P_3`$ which has the largest $`Q`$ has the smallest threshold for lasing. So, in the frequency-independent gain case, the transfer matrix method gives that the mode with largest $`Q`$ will lase first. We have also used the transfer-matrix method with a frequency-dependent gain profile, given as a dotted line in Fig. 1a. We choose the central frequency as $`f_a=\omega _a/2\pi =618.56`$ THz, which is exactly between the $`P_1`$ and $`P_2`$, and its width to be $`\mathrm{\Delta }f_a=\mathrm{\Delta }\omega _a/2\pi =15`$ THz. In this case, when we increase the pumping rate $`P_r`$, we find the $`P_2`$ will lase first, then $`P_1`$. So two important conditions determine which mode will lase first, the first is the quality factor of the mode and the second is the gain profile. Experimentally, it is quite often that only part of the random medium has been pumped. Then a third factor comes in, i.e. the spatial overlap of the mode function and the gain region.

The next issue we address is the shape of the wave function, as one introduces gain. In Fig. 2a, we present the amplitude of the electric field versus distance for $`f=626.6`$ THz, which corresponds to the peak $`P_2`$ of Fig. 1a. In Fig. 2b, we show the wavefunction with near-threshold gain at the exact high-peak frequency of $`P_2`$. Both Fig. 2a and 2b are obtained by the time-independent transfer matrix method. Notice that in the presence of gain, the shape of the wavefunction in Fig. 2b is almost same as that without gain in Fig. 2a. The only change is its amplitude which increases uniformly. Actually, by keeping the incident amplitude the same, when we increase the gain from zero to the threshold value, we find that the amplitude of the wavefunction increases from a small to a very large value, but its shape remains almost same.

To check if indeed this amazing property is also obeyed by the full time-dependent theory  of the semiclassical laser theories with Maxwell-equations, we repeat our theory for this case too. While in the time-independent theory, every mode is independent and amplification is not saturated, this is not true for a real lasing system. In a real lasing system, modes will compete with each other. As was discussed in Ref. 12, in a $`short`$ system, the first lasing mode will suppress all other modes, so we observe only one lasing mode even for a large pumping rate. This is exactly the case when we use the FDTD method to simulate our 30-cell system with an over-threshold gain whose gain profile is the same as given in Fig.1a. At first, the electric field is a random one due to spontaneous emission, then a strong lasing mode evolves from the noisy background and a sharp peak appears in the emission spectrum after Fourier transform. The stable field profile of the FDTD calculation is given in Fig. 2c, where the wavefunction is quite the same as in Fig. 2a and 2b, but the amplitude is very large. The emission spectrum of this case is obtained, and gives a sharp peak very close to the resonant peak $`P_2`$ of Fig. 1a. We also checked the wavefunctions as well as the lasing threshold as we shift the gain profile. If the central frequency of the gain profile is near $`P_3`$, the $`P_3`$ mode will lase first, and the shape of wavefunction of $`P_3`$ remains unchanged. We have also checked the above ideas for a larger ($`L\xi `$, where $`\xi `$ is the localization length) system. In such strong localized case, the form of the wavefunction remains unchanged as one introduces gain.

The numerical results which are presented above can be explained by the semi-classical Lamb theory  of laser physics. According to Lamb theory, the Maxwell equation of the random laser system can be written as

$`^2E(x,t)+\mu _0\sigma {\displaystyle \frac{E(x,t)}{t}}+\mu _0\epsilon (x){\displaystyle \frac{^2E(x,t)}{t^2}}=\mu _0{\displaystyle \frac{^2P_{gain}(x,t)}{t^2}}`$ (1)
Where $`\mu _0`$ is the permeability of free space, $`E(x,t)`$ is the electric field , and the dielectric constant $`\epsilon (x)`$ is determined by the random configuration of the system. $`\sigma `$ is not just the common conductivity loss, but can be interpreted by the total mode loss including the surface loss of the system by radiation. The polarization due to gain $`P_{gain}(x,t)`$  is the same as the one defined above.

We assume that the system satisfies the slow varying approximation (it is always satisfied if we care about the stable lasing state). So $`E(x,t)=E_m(t)\varphi _m(x)exp(i\omega t)`$ (our later discussion shows that the separation of spatial and time parts of wavefunction is reasonable), where $`E_m(t)`$ is the field amplitude, $`\varphi _m(x)`$ is the normalized wavefunction, and $`\omega `$ is the frequency of the field. The surface loss of the mode is $`\sigma _m=\epsilon _0\omega /Q_m`$, where the $`Q_m`$ is the quality factor of the mode. Then, we get two equations  for the real and the imaginary terms of Eq. (1)

$`^2\varphi _m(x)+\mu _0(\epsilon (x)+\epsilon _0\chi ^{}(x,\omega ))\omega ^2\varphi _m(x)=0`$ (2)
$`{\displaystyle \frac{E_m(t)}{t}}=\left(\overline{\chi ^{\prime \prime }(\omega )}{\displaystyle \frac{1}{Q_m}}\right){\displaystyle \frac{\epsilon _0\omega E_m(t)}{2\overline{\epsilon }}}`$ (3)
where $`\overline{\epsilon }=_0^L\epsilon (x)𝑑x/L`$ is the spatial averaged dielectric constant inside the system, and $`\overline{\chi ^{\prime \prime }(\omega )}=_0^L\varphi _m(x)^{}\chi ^{\prime \prime }(x,\omega )\varphi _m(x)𝑑x/L`$ is the spatial averaged gain. The last integral is done to take into account of the overlap between the wavefunction and the spatial region of gain.

Equation (2) determines the field distribution, quality factor $`Q_m`$, and the vibration frequency $`\omega `$ of the lasing mode. The term $`\epsilon _0\chi ^{}(x,\omega )`$ will cause the vibration frequency to shift away from the original eigenfrequency of the mode, called pulling effect. For a well-defined mode, generally $`Q_m1`$, we need a very small gain to lase. Then $`\chi ^{}(x,\omega )1`$, the pulling effect is very weak, and $`\omega \mathrm{\Omega }_m`$, where $`\mathrm{\Omega }_m`$ is the eigenfrequency of the mode. So the wavefunction of the lasing mode should be similar to the eigenfunction of the mode. Theoretically, we can use the perturbation method to obtain $`\omega `$ and the wavefunction.

Equation (3) is the time-dependent amplitude equation. First we can use it to determines the threshold condition when $`\overline{\chi }^{\prime \prime }(\omega )=1/Q_m`$. For our system, with homogeneous pumping, we have

$`P_r^c=C_0{\displaystyle \frac{1+\frac{4(\omega _a\mathrm{\Omega }_m)^2}{\mathrm{\Delta }\omega _a^2}}{Q_m}}`$ (4)
where $`C_0=\frac{a_0+b_0}{a_0}\frac{\gamma _cm_e\epsilon _0\omega _a\mathrm{\Delta }\omega _a}{\gamma _rN_0\tau _{21}e^2}`$ is a constant. Eq. (4) indeed gives that the threshold value of lasing is inversely proportional to the quality factor $`Q`$.

Second, Eq. (3) give us the stable amplitude of field when the gain is over threshold. Actually, the gain is saturable, $`\chi ^{\prime \prime }(\omega )\mathrm{\Delta }N1/(1+C_2|E_m|^2)`$ , in real systems and in our FDTD model . With an over-threshold gain, $`E_m(t)`$ will increase and the gain parameter $`\chi ^{\prime \prime }(\omega ,E_m)`$ will decrease until $`\overline{\chi }^{\prime \prime }(\omega ,E_m)=1/Q_m`$, then the field is stable. So Eq. (3) also determines the amplitude of the stable field for over-threshold pumping cases.

Our numerical and analytical results clearly suggest that states of the random system with gain can easily lase, provided their $`Q`$-factor is large. These findings provide a new path for observing localization of light. Since localized states have large $`Q`$ values, will lase with a small pumping rate. On the other hand, strongly fluctuating extended states have smaller $`Q`$ values because of the radiation loss on the surface of the system and can lase only after a stronger pumping. In a real experiment, even in the presence of absorption, if the gain profile is close to the mobility edge, there is going to be a discontinuity in the critical pumping rate needed for lasing. Localized states will lase first at a low pumping rate, while extended states need a high pumping rate. In Fig. 3, we present the results of the DOS vs frequency $`f`$ for a 1d quasi-disordered system. On the insert two eigenfunctions are given, one localized at $`f=468`$ THz and the one “extended” at $`f=485`$ THz. In addition the emission intensity $`I_e`$ vs $`f`$ is also shown with gain, $`P_r=10^5`$ $`1/s`$, which demonstrate that the localized modes will lase first. In a more realistic 3D case, we expect $`I_e`$ vs frequency to have many peaks for the localized region and no peaks in the extended region for a give pumping rate . It would be very interesting if this discontinuity can be observed experimentally.

In summary, we have used the time-independent transfer matrix method and the FDTD method to show that all lasing modes are coming from the eigenstates of the random system. A knowledge of the eigenstates and the density-of-states of the random system, in conjunction with the frequency profile of the gain, can accurately predict the mode that will lase first, as well as its critical pumping rate. Our detailed numerical results clearly demonstrate that the shape of the wavefunction remains unchanged as gain is introduced into the system. The role of the gain is to just increase uniformly the amplitude of the wavefunction without changing its shape. These results can be understood by generalizing the semi-classical Lamb theory of lasing in random systems. These findings can help us unravel the conditions of observing the localization of light, as well as manufacturing random laser with specific emission wavelengths.

Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was supported by the director for Energy Research, Office of Basic Energy Sciences.


text:
# The interaction of radio sources and cooling flows

## 1. Introduction

The X-ray surface brightness continues to rise towards the centre in two thirds of all clusters of galaxies. The radiative cooling time of the gas there continues to decrease, dropping below 10 Gyr at 70–150 kpc and 1 Gyr at 5–30 kpc. The smaller radius values apply to low luminosity clusters such as that in Virgo, the larger ones to the high luminosity ones such as that in Perseus. Within this inner region the gas temperature also drops, from 3–10 keV down to 1–3 keV, as determined by simple single-temperature spectral fits. The gross appearance of these regions fits that of a cooling flow (see Fabian 1994 for a review), in which radiative cooling of the central dense parts of a hydrostatic atmosphere causes it to subsonically slump inward under the influence of gravity and the weight of overlying gas. About 10 per cent of clusters, such as the Coma cluster, do not show the central peak in X-ray emission and have no cooling flow. Another 10–20 per cent of clusters, such as A754, have complex structure, probably due to a cluster-cluster merger.

The inferred rate at which gas is cooling out in cooling flows ranges from $`<10\mathrm{M}_{}\mathrm{yr}^1`$ in poor clusters to $`>1000\mathrm{M}_{}\mathrm{yr}^1`$ in massive rich clusters. These values are uncertain and depend on an assumed age for the flow if determined only from x-ray imaging. Only when X-ray spectra show the cooling components can the rate be clearly defined. In that case an estimate for the age of the flow is obtained as well. Work using ROSAT and ASCA data (Allen & Fabian 1997; Allen et al 1999) indicates that the ages may be about 5 Gyr. Many clusters appear to have rates of a few $`100\mathrm{M}_{}\mathrm{yr}^1`$, which suggests total cooled gas masses of $`10^{11}10^{12}\mathrm{M}_{}`$. Blue light from massive stars and emission line nebulosities are common in such regions (see Crawford et al 1999 and references therein) but normal star formation accounts for only about 10 per cent of the cooling rate.

What happens to the rest is not clear. It has been suggested that heating balances cooling so that little or no gas actually cools from the X-ray emitting hot gas component (Rosner & Tucker 1989; Binney & Tabor 1995; Soker et al 2000). This seems plausible when it is noted that most cooling flows have a radio source in the host central cluster galaxy. If somehow the accretion power from a small amount of gas accreted onto a central black hole can be distributed into the cooler and cooling gas then perhaps some equilibrium may be obtained. Thermal conduction too must be worried about since we are discussing cooler gas components embedded in hotter ones.

The problem here is that the X-ray data show little sign that the radio sources are actually heating the bulk of the gas. Of course there is something going on in the radio lobes, but these occupy only a small fraction of the flow. What is needed to stem a cooling flow is a lot of energy into the cooler gas. It has to hit the gas with the shortest cooling times. The heating needs to be complete in the sense that it shouldn’t just cause the gas to hang up at some temperature below the outer cluster temperature or it would easily be detected. It must also allow gas to cool (or give the appearance of cooling) down to 1–3 keV. This is not trivial.

Here I review the recent results from Chandra which resolve the radio – X-ray interaction regions in some luminous clusters. The simple result is that the coolest gas is seen to be closest to the radio lobes. Not what a heating model would be expected to show if the lobes are responsible for the heating.

## 2. Chandra results

### 2.1. Hydra A and 3C295

Early data from Chandra showed clear evidence for ‘holes’ in the X-ray in the vicinity of the radio lobes in the Hydra A cluster (McNamara et al 2000). The central radio source here is fairly powerful and the lobes occupy one of the larger fractions of the cooling region. Later work suggests that a small cooling flow may operate only near the centre of the region (David et al 2000) at a rate compatible with the star formation seen.

The radio lobes in the powerful FRII source 3C295 are seen in X-ray emission, presumable as a consequence of inverse-Compton scattering (Harris et al 2000). The surrounding emission does however appear to be reasonably undisturbed and undergoing a cooling flow (Allen et al 2000).

### 2.2. The Perseus cluster and 3C84

The core of the Perseus cluster around the central galaxy NGC1275 and its radio source 3C84 has long been known to have X-ray ‘holes’ at the positions of the radio lobes (Bohringer et al 1993). ROSAT data showed good agreement between the holes and inner lobes. Chandra imaging now shows this very clearly (Fabian et al 2000a; Fig. 1).

The exact 3D geometry is not obvious (Fig. 2). Presumably the S jet is pointing toward us at some angle. There are also several outer holes seen (Fabian et al 2000a) which may represent older, buoyant lobes (Churazov et al 2000a). The energetic electrons may now have cooled in these regions but the nonthermal pressure from the protons, magnetic field and cooled electrons may keep them inflated and less dense than the surroundings. A spur of radio emission toward the W outer hole is seen in the 74 MHz map of Blundell et al (2000) which supports this idea.

The Chandra data (Fabian et al 2000a) covers the 0.3–7 keV range and can be divided into colours from which absorption and temperature maps can be constructed. Significant excess absorption is seen in the region of the ‘high velocity system’, an irregular galaxy falling into the cluster core close to our line of sight. The temperature map (Fig. 3) shows that the gas is cooler toward the centre with the coolest gas along the rims of the inner radio lobes and in the bright E blob. The surface brightness can be used to obtain the gas density and thus pressure and radiative cooling time. The cooling times are also shortest close to the inner radio lobes.

We see no sign that the lobes have heated the gas beyond the lobes. Something has happened within the lobes. They have either been cleared of gas, or perhaps just cleared, of cooler gas. The intracluster medium there may be multiphase with cooler, cooling, blobs embedded within a hotter phase. The lobes may have only pushed aside the cooler blobs. A simple rearrangement of the cooler gas is in agreement with the observed surface brightness excess around the rims of the lobes. There is certainly no sign of the strong shocks suggested by Heinz et al (1998).

### 2.3. A1795 and A2199

Other clusters with central radio sources which we have observed are A1795 (Fabian et al 2000) and A2199 (Johnstone et al 2001, in preparation). A1795 shows a filament of soft X-ray emission (Fig. 4) which coincides with an optical H$`\alpha `$ filament discovered by Cowie et al in 1983. This could be a cooling wake due to motion of the central cluster galaxy. The velocity map of the H$`\alpha `$ emission (Hu et al 1985) shows that the wake has the same velocity as the bulk of the cluster (Oegerele & Hill 1994), there is a velocity shift of $`+150\mathrm{km}\mathrm{s}^1`$ at the cD galaxy. The small radio source associated with the cD (Ge & Owen 1993) shows some interaction with the brighter local X-ray emission. A possible X-ray ‘shadow’, a fainter X-ray region, parallel to but further out from the radio lobes, is also seen.

A2199 has large ‘holes’ conicident with its more diffuse outer radio structures (Fig. 5; radio data from Giovannini et al 1998, see also Owen & Eilek 1998) but there is little agreement between the inner radio and X-ray structures. Analysis of the Chandra data on A2199 is in progress (Johnstone et al 2001).

### 2.4. The Virgo cluster and M87

Chandra images of the inner parts of the Virgo cluster around M87 have been made in the ACIS-I detector by Fabian et al (2001 in preparation) and in the ACIS-S detector by Wilson et al (2001). The excess emission along the line of the radio structures, seen in ROSAT data (Bohringer et al 1995), is clearly detected (Fig. 6). Only the X-ray structure to the E shows any good correspondence with the radio images (Owen et al 2000). The excess emission there is thermal and cooler than the surrounding gas (Bohringer et al 1995) and possibly due to uplift of cooler gas from nearer the centre by a buoyant radio plume (Bohringer et al 1995; Churazov et al 2000b).

## 3. Summary

The Chandra data so far show little evidence for any widespread steady heating by the radio source. Perhaps there are sporadic upheavals due to large outbursts by the central engine, or a pervasive leakage of energetic particles. Or perhaps cooling dominates.

Conduction may be highly suppressed in cooling flows. Chandran et al (1999) have discussed theoretical possiblilities including the mirroring of particles. The Chandra observations of ‘cold fronts’ in some clusters (Markevich et al 2000; Vikhlinin et al 2000) where 5 keV gas lies a few kpc away from 10 keV gas suggests observationally that conduction may in some intracluster regions be highly suppressed (Ettori & Fabian 2000).

## 4. Post-conference appendix

In between the conference and the writing of this paper, Reflection-Grating Spectrometer (RGS) data from XMM have appeared. These have a strong bearing on the appearance of cooling flows so a brief outline of the issues is now given to steer the reader to the relevant literature.

The RGS data of the strong, hot, cooling flow cluster A1835 (Peterson et al 2000) show none of the strong emission lines expected from gas cooling below 3 keV. It seems like the gas cools from about 9 keV down to about 3 keV then vanishes. Similarly there are no lines from gas below about 2 keV from A1795 (Tamura et al 2000) or A1101 (Kaastra et al 2000). The temperature profiles of the gas drop inward in the manner expected from a cooling flow but the spectral signature of the final temperature plunge is missing.

There are many possible explanations, discussed by Peterson et al (2000) and Fabian et al (2000c). Perhaps there is heating. But why then does the gas apparently cool most of the way and is then all heated back up to its original temperature? (It cannot be heated to some other intermediate temperature or it would be detected as accumulating there.) Perhaps the ages of the flows are all only about one Gyr. Perhaps resonance scattering and some absorption operate on the inner, cooler gas making it undetectable. Perhaps the metallicity of the gas is bimodal. A recent paper on M87 by Bohringer et al (2000) shows that the situation is complex, and metallicity, resonance scattering or continuum dilution may be factors. They find that the metallicity gradient in the Virgo cluster in which the iron (and other element) abundance rises inwards reverses in the innermost arcmin resulting in the metallicity going to zero! Not only are no lines seen from cooling gas, but there are no lines from anything else. Clearly further observations, and thought, are required.

## References

Allen S.W., Fabian A.C., 1997, MNRAS, 286, 583

Allen S.W., Fabian A.C., Johnstone R.M., Nulsen P.E.J., Arnaud K.A., 1999, MNRAS, in press, astro-ph/9910188

Allen S.W. et al 2000, MNRAS in press

Binney J., Tabor G., 1995, MNRAS, 276, 663

Blundell K., Kassim N., Perley R., 2000, astro-ph004005

Bohringer H., Voges W., Fabian A.C., Edge A.C., Neumann D.M., 1993, MNRAS, 264, L25

Böhringer H., Nulsen P.E.J., Braun R., Fabian A.C., 1995, MNRAS, 274, L67

Böhringer H. et al 2000, A&A, in press, astro-ph 0011459

Bregman J.N., David L.P., 1988, ApJ, 326, 639

Chandran B., Cowley S., 1998, Physical Review Lett, 80, 3077 (1998)

Churazov E., Forman W., Jones C., Bohringer H., 2000, A&A, 356, 788

Churazov E., Bruggen M., Kaiser C., Bohringer H., Forman W., 2000, ApJ submitted (astro-ph 0008215)

Cowie L.L., Hi E.M., Jenkins E.B., York D.G., 1983 ApJ, 272, 29

Crawford C.S., Allen S.W., Ebeling H., Edge A.C., Fabian A.C., 1999, MNRAS, 306, 875

David L, Nulsen P.E.J., McNamara B.R., Forman, W., Jones C., Robertson B., Wise M., 2000 , ApJ, submitted, astro-ph/0010224

Ebeling H., Edge A.C., Fabian A.C., Allen S.W., Crawford C.S., Böhringer H., 1997, ApJ, 479, L101

Ettori S., Fabian A.C., 2000, MNRAS, 317, L57

Fabian A.C., 1994, ARAA, 32, 277

Fabian A.C., et al 2000a, MNRAS, 318, L65

Fabian A.C., et al 2000b, MNRAS in press (astro-ph/0011547)

Fabian A.C., Mushotzky R.F., Nulsen P.E.J., Peterson J., 2000 MNRAS in press (astro-ph/0010509)

Giovannini G., Cotton W.D., Feretti L., Lara L., Venturi T., 1998, ApJ, 493, 632

Ge J.P., Owen F.N., 1993, AJ, 105, 778

Harris D., et al, 2000, ApJ, 530, L81

Heinz S., Reynolds C.S., Begelman M.C., 1998, ApJ, 501, 126

Hu E.M., Cowie L.L., Wang Z., 1985, ApJS, 59, 447

Kaastra J., et al, 2000, A&A in press (asttro-ph/0010423)

Markevitch M. et al, 2000, ApJ, 541, 542

McNamara B.R. et al 2000, ApJ submitted, astro-ph/0001402

McNamara B. et al, 2000a, ApJ, 534, L135

Oegerle, W.R., Hill, J.M., 1994, AJ, 107, 857

Owen F.N., Eilek J.A., Kassim N.E., 2000, ApJ, 543, 611

Owen F.N., Eilek J.A., 1998, ApJ, 493, 730

Peterson J.A. et al 2000, A&A in press (astro-ph 0010658)

Rosner R., Tucker W., 1989, ApJ, 338, 761

Soker N., White R.E., David L.P., McNamara B.R., 2000, ApJ, submitted, astro-ph/009173

Tamura T. et al 2000, A&A in press, astro-ph/0010362

Tucker W.H., Rosner R., 1983, 267, 547

Tucker W., David L., 1997, ApJ, 484, 602

Vikhlinin A., Markevitch M., Murray S.S., 2000, ApJ, submitted, astro-ph/0008496

Wu K.K.S., Fabian A.C., Nulsen P.E.J., 2000, MNRAS, 318, 889


text:
# 1 Introduction

## 1 Introduction

Kinetic simulations of protein folding are notoriously difficult. Thermodynamic simulations may use unphysical moves and are therefore potentially easier, but existing methods need improvement. Three properties that a successful thermodynamic algorithm must possess are as follows. First and foremost, it must be able to alleviate the multiple-minima problem. Methods like the multicanonical algorithm  and simulated tempering  try to do so by the use of generalized ensembles. Second, it must provide an efficient evolution of large-scale properties of unfolded chains. The simple pivot method  does remarkably well  in that respect. Third, it must be able to alter local properties of folded chains without causing too drastic changes in their global structure. This paper is concerned with the third problem, which is important if the backbone potentials are stiff and especially if the mobility is restricted to the biologically most relevant torsional degrees of freedom.

An update that rearranges a restricted section of the chain without affecting the remainder is local. For chains with flexible or semiflexible backbones, there exists a variety of local updates, ranging from simple single-site moves to more elaborate methods  where inner sections are removed and then regrown site by site in a configurational-bias manner . However, these methods break down if bond lengths and bond angles are completely rigid.

The problem of generating local deformations of chains with only torsional degrees of freedom was analyzed in a classic paper by Gō and Scheraga . Based on this analysis, Dodd et al.  devised the first proper Monte Carlo algorithm of this type, the concerted-rotation method. This method works with seven adjacent torsion angles along the chain. One of these angles is turned by a random amount. Possible values of the remaining six angles are then determined by numerically solving a set of equations that guarantee that the move is local. The new conformation is finally drawn from the set of all possible solutions to this so-called rebridging problem. Variations and generalizations of this method have been discussed by several groups . There are also methods  that combine elements of the configurational-bias and concerted-rotation approaches. One of these methods  uses an analytical rebridging scheme, inspired by the solution for a similar problem in robotic control .

The concerted-rotation approach is a powerful method that can generate large local deformations by finding the discrete solutions to the rebridging problem. However, the method is not easy to implement and large local deformations may be difficult to accomplish if, for example, the chain is folded and has bulky side groups. Hence, there are situations where this method is not the obvious choice.

In this paper, we discuss a different and less sophisticated type of Monte Carlo move in torsion angle space. This algorithm is by nature a “small-step” algorithm so large local deformations cannot take place. Drastic global changes would still occur if the steps were random. To avoid that, a biasing probability is introduced. The method becomes approximately local if the bias is made strong. Compared to a strictly local update, this method has the disadvantage that a much smaller part of the energy function is left unchanged, so the CPU time per update is larger. However, this problem is not too severe for moderate chain lengths. Moreover, both our method and strictly local ones are typically combined with some truly nonlocal update like pivot, and such an update is not faster than ours.

The algorithm proceeds as follows. We consider $`n`$ torsion angles $`\varphi _i`$, where $`n=8`$ in our calculations. To update these angles, we introduce a conformation-dependent $`n\times n`$ matrix $`𝐆`$ such that $`\delta \overline{\varphi }^T𝐆\delta \overline{\varphi }0`$ for changes $`\delta \overline{\varphi }=(\delta \varphi _1,\mathrm{},\delta \varphi _n)`$ that correspond to local deformations. The steps $`\delta \overline{\varphi }`$ are then drawn from the Gaussian distribution

$$P(\delta \overline{\varphi })\mathrm{exp}\left[\frac{a}{2}\delta \overline{\varphi }^T(\mathrm{𝟏}+b𝐆)\delta \overline{\varphi }\right],$$
(1)
where $`\mathrm{𝟏}`$ denotes the $`n\times n`$ unit matrix and $`a`$ and $`b`$ are tunable parameters. The parameter $`a`$ controls the acceptance rate, whereas $`b`$ sets the degree of bias. The new conformation is finally subject to an accept/reject step. Important to the implementation of the algorithm is that the matrix $`𝐆`$ is non-negative and symmetric. Hence, it is possible to take the “square root” of $`\mathrm{𝟏}+b𝐆`$, which facilitates the calculations.

This method, which is quite general, is tested on a reduced model protein  with 54 amino acids and the Ramachandran torsion angles as its only degrees of freedom. This chain forms a three-helix bundle in its native state and exhibits an abrupt collapse transition that coincides with its folding transition. The performance of the method is studied both above and below the folding temperature, for different values of the parameters $`a`$ and $`b`$. For a suitable choice of $`b`$, we find that the effective step size can be increased by a factor of three in the folded phase, compared to the unbiased $`b=0`$ case. The optimal value of $`b`$ corresponds to a relatively strong bias, that is an approximately local update.

## 2 The Model

In our calculations, we consider a reduced protein model  where each amino acid is represented by five or six atoms. The three backbone atoms N, Cα and C are all included, whereas the side chain is represented by a single atom, Cβ. The Cβ atom can be hydrophobic, polar or absent, which means that there are three different types of amino acids in the model. For a schematic illustration of the chain representation, see Fig. 1.

All bond lengths, bond angles and peptide torsion angles ($`180^{}`$) are held fixed, which leaves us with two degrees of freedom per amino acid, the Ramachandran torsion angles (see Fig. 1).

The energy function

$$E=E_{\text{loc}}+E_{\text{sa}}+E_{\text{hb}}+E_{\text{AA}}$$
(2)
is composed of four terms. The local potential $`E_{\text{loc}}`$ has a standard form with threefold symmetry,

$$E_{\text{loc}}=\frac{ϵ_{\text{loc}}}{2}\underset{i}{}(1+\mathrm{cos}3\varphi _i).$$
(3)
The self-avoidance term $`E_{\text{sa}}`$ is given by a hard-sphere potential of the form

$$E_{\text{sa}}=ϵ_{\text{sa}}\stackrel{}{}_{i<j}\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12},$$
(4)
where the sum runs over all possible atom pairs except those consisting of two hydrophobic Cβ. The hydrogen-bond term $`E_{\text{hb}}`$ is given by

$$E_{\text{hb}}=ϵ_{\text{hb}}\underset{ij}{}u(r_{ij})v(\alpha _{ij},\beta _{ij}),$$
(5)
where $`i`$ and $`j`$ represent H and O atoms (see Fig. 1), respectively, and

$`u(r_{ij})`$ $`=`$ $`5\left({\displaystyle \frac{\sigma _{\text{hb}}}{r_{ij}}}\right)^{12}6\left({\displaystyle \frac{\sigma _{\text{hb}}}{r_{ij}}}\right)^{10}`$ (6)
$`v(\alpha _{ij},\beta _{ij})`$ $`=`$ $`\{\begin{array}{cc}\mathrm{cos}^2\alpha _{ij}\mathrm{cos}^2\beta _{ij}\hfill & \alpha _{ij},\beta _{ij}>90^{}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (9)
In these equations, $`r_{ij}`$ denotes the HO distance, $`\alpha _{ij}`$ the NHO angle, and $`\beta _{ij}`$ the HOC angle. Finally, the hydrophobicity term $`E_{\text{AA}}`$ has the form

$$E_{\text{AA}}=ϵ_{\text{AA}}\underset{i<j}{}\left[\left(\frac{\sigma _{\text{AA}}}{r_{ij}}\right)^{12}2\left(\frac{\sigma _{\text{AA}}}{r_{ij}}\right)^6\right],$$
(10)
where both $`i`$ and $`j`$ represent hydrophobic Cβ. In the following, $`kT`$ is given in dimensionless units, in which $`ϵ_{\text{hb}}=2.8`$ and $`ϵ_{\text{AA}}=2.2`$. Further details of the model, including numerical values of all the parameters, can be found in Ref. .

In this model, we study a designed three-helix-bundle protein with 54 amino acids. In Ref. , it was demonstrated that this sequence indeed forms a stable three-helix bundle, except for a twofold topological degeneracy, and that it has a first-order-like folding transition that coincides with the collapse transition. It should be noted that these properties are found without resorting to the widely used but drastic Gō approximation , where interactions that do not favor the desired structure are ignored.

## 3 The Algorithm

We now turn to the algorithm, which we describe assuming the particular chain geometry defined in Sec. 2. That this scheme can be easily generalized to other types of chains will be evident.

Consider a segment of four adjacent amino acids $`k`$, $`k+1`$, $`k+2`$ and $`k+3`$ along the chain, and let the corresponding eight Ramachandran angles (see Fig. 1) form a vector $`\overline{\varphi }=(\varphi _1,\mathrm{},\varphi _n)`$, where $`n=8`$. A change $`\delta \overline{\varphi }`$ of $`\overline{\varphi }`$ will, by construction, leave all amino acids $`k^{}<k`$, as well as the N, H and Cα atoms of amino acid $`k`$, unaffected. For all amino acids $`k^{}>k+3`$ to remain unaffected too, it is sufficient to require that the three atoms Cα, C and O of amino acid $`k+3`$ (see Fig. 1) do not move. If this condition is fulfilled, the deformation of the chain is local.

Denote the position vectors of the Cα, C and O atoms of amino acid $`k+3`$ by $`\overline{r}_I`$, $`I=1,2,3`$. A bias toward local deformations can be obtained by favoring changes $`\delta \overline{\varphi }`$ that correspond to small values of the quantity

$$\mathrm{\Delta }^2=\underset{I=1}{\overset{3}{}}(\delta \overline{r}_I)^2,$$
(11)
which for small $`\delta \varphi _i`$ can be written as

$$\mathrm{\Delta }^2\delta \overline{\varphi }^T𝐆\delta \overline{\varphi }=\underset{i,j=1}{\overset{n}{}}\delta \varphi _iG_{ij}\delta \varphi _j,$$
(12)
where

$$G_{ij}=\underset{I=1}{\overset{3}{}}\frac{\overline{r}_I}{\varphi _i}\frac{\overline{r}_I}{\varphi _j}.$$
(13)
Note that the three vectors $`\overline{r}_I`$ can be described in terms of six independent parameters, since bond lengths and angles are fixed. This implies that the $`n\times n`$ matrix $`𝐆`$, which by construction is non-negative and symmetric, has eigenvectors with eigenvalue zero for $`n=8>6`$. A bias toward small $`\mathrm{\Delta }^2`$ means that these soft modes are favored.

We can now define the update, which consists of the following two steps.

1. Draw a tentative new $`\overline{\varphi }`$, $`\overline{\varphi }^{}`$, from the Gaussian distribution

$$W(\overline{\varphi }\overline{\varphi }^{})=\frac{(det𝐀)^{1/2}}{\pi ^3}\mathrm{exp}\left[(\overline{\varphi }^{}\overline{\varphi })^T𝐀(\overline{\varphi }^{}\overline{\varphi })\right],$$
(14)
where the matrix

$$𝐀=\frac{a}{2}(\mathrm{𝟏}+b𝐆)$$
(15)
is a linear combination of the $`n\times n`$ unit matrix $`\mathrm{𝟏}`$ and the matrix $`𝐆`$ defined by Eq. 13. The shape of this distribution depends on the parameters $`a>0`$ and $`b0`$. The parameter $`b`$ sets the degree of bias toward small $`\mathrm{\Delta }^2`$. The bias is strong for large $`b`$ and disappears in the limit $`b0`$. The parameter $`a`$ is a direction-independent scale factor that is needed to control the acceptance rate. Larger $`a`$ means higher acceptance rate, for fixed $`b`$. If $`b=0`$, then the components $`\delta \varphi _i`$ are independent Gaussian random numbers with zero mean and variance $`a^1`$. Note that $`W(\overline{\varphi }\overline{\varphi }^{})W(\overline{\varphi }^{}\overline{\varphi })`$ since the matrix $`𝐆`$ is conformation dependent.
2. Accept/reject $`\overline{\varphi }^{}`$ with probability

$$P_{\text{acc}}=\mathrm{min}(1,\frac{W(\overline{\varphi }^{}\overline{\varphi })}{W(\overline{\varphi }\overline{\varphi }^{})}\mathrm{exp}[(E^{}E)/kT])$$
(16)
for acceptance. The factor $`W(\overline{\varphi }^{}\overline{\varphi })/W(\overline{\varphi }\overline{\varphi }^{})`$ is needed for detailed balance to be fulfilled, since $`W`$ is asymmetric.

It should be stressed that this scheme is quite flexible. For example, it can be immediately applied to chains with nonplanar peptide torsion angles. The use of the concerted-rotation method for simulations of such chains has recently been discussed .

A convenient and efficient implementation of the algorithm can be obtained if one takes the “square root” of the matrix $`𝐀`$, which can be done because $`𝐀`$ is symmetric and positive definite. More precisely, it is possible to find a lower triangular matrix $`𝐋`$ (with nonzero elements only on the diagonal and below) such that

$$𝐀=\mathrm{𝐋𝐋}^T.$$
(17)
An efficient routine for this so-called Cholesky decomposition can be found in .

### 3.1 Implementing step 1

Given the Cholesky decomposition of the matrix $`𝐀`$, the first step of the algorithm can be implemented as follows.

* Draw a $`\overline{\psi }=(\psi _1,\mathrm{},\psi _n)`$ from the distribution $`P(\overline{\psi })\mathrm{exp}(\overline{\psi }^T\overline{\psi })`$. The components $`\psi _i`$ are independent Gaussian random numbers and can be generated, for example, by using the Box-Muller method

$$\psi _i=(\mathrm{ln}R_1)^{1/2}\mathrm{cos}2\pi R_2,$$
(18)
where $`R_1`$ and $`R_2`$ are uniformly distributed random numbers between 0 and 1.
* Given $`\overline{\psi }`$, solve the triangular system of equations

$$𝐋^T\delta \overline{\varphi }=\overline{\psi }$$
(19)
for $`\delta \overline{\varphi }`$. It can be readily verified that the $`\delta \overline{\varphi }=\overline{\varphi }^{}\overline{\varphi }`$ obtained this way has the desired distribution Eq. 14.

### 3.2 Implementing step 2

The Cholesky decomposition is also useful when calculating the acceptance probability in the second step of the algorithm. The factor $`W(\overline{\varphi }\overline{\varphi }^{})`$ can be easily computed by using that

$$(det𝐀)^{1/2}=\underset{i=1}{\overset{n}{}}L_{ii}$$
(20)
and that $`\mathrm{exp}[(\overline{\varphi }^{}\overline{\varphi })^T𝐀(\overline{\varphi }^{}\overline{\varphi })]=\mathrm{exp}(\overline{\psi }^T\overline{\psi })`$. The reverse probability $`W(\overline{\varphi }^{}\overline{\varphi })`$ depends on $`𝐀(\overline{\varphi }^{})`$ and can be obtained in a similar way, if one makes a Cholesky decomposition of that matrix, too.

### 3.3 Pivot update

Previous simulations  of the model protein defined in Sec. 2 were carried out by using simulated tempering with pivot moves as the elementary conformation update. With this algorithm, the system was successfully studied down to temperatures just below the folding transition. However, the performance of the pivot update, where a single angle $`\varphi _i`$ is turned, deteriorates in the folded phase. What we hope is that the exploration of this phase can be made more efficient by alternating the pivot moves with moves of the type described previously.

## 4 Results

The character of the proposed update depends strongly on the bias parameter $`b`$. The suggested steps have a random direction if $`b=0`$. The distribution $`W(\overline{\varphi }\overline{\varphi }^{})`$ in Eq. 14 is, by contrast, highly asymmetric in the limit $`b\mathrm{}`$, with nonzero width only in directions corresponding to eigenvalue zero of the matrix $`𝐆`$. In particular, this implies that the reverse probability $`W(\overline{\varphi }^{}\overline{\varphi })`$ in the acceptance criterion Eq. 16 tends to be small for large $`b`$.

For the acceptance rate to be reasonable, it is necessary to use a very small step size if $`b`$ is small or large. The question is whether the step size can be increased by a better choice of $`b`$. To find that out, we performed a set of simulations of the three-helix-bundle protein defined in Sec. 2 for different $`a`$ and $`b`$. Two different temperatures were studied, $`kT=0.6`$ and 0.7, one on either side of the folding temperature $`kT_f0.66`$ .

In these runs, we monitored the step size $`S`$, where

$$S=|\delta \overline{\varphi }|=\left[\underset{i=1}{\overset{n}{}}(\delta \varphi _i)^2\right]^{1/2}$$
(21)
for accepted moves and $`S=0`$ for rejected ones. Measurements were taken only when the $`n=8`$ angles all were in the segment that makes the middle helix of the three. We focus on this segment because it is the most demanding part to update.

The average step size, $`S`$, depends strongly on $`b`$. A rough optimization of $`b`$ was carried out by maximizing $`S`$ as a function of $`a`$ for different fixed $`b=10^k`$ ($`k`$ integer). The best values found were $`b_{\mathrm{max}}=10(\mathrm{rad}/\mathrm{\AA })^2`$ and $`b_{\mathrm{max}}=0.1(\mathrm{rad}/\mathrm{\AA })^2`$ at $`kT=0.6`$ and $`kT=0.7`$, respectively. Note that the preferred degree of bias is higher in the folded phase.

In Fig. 2, we show $`S`$ against the average acceptance rate, $`P_{\text{acc}}`$, for $`b=0`$ and $`b=b_{\mathrm{max}}`$ at the two temperatures; $`P_{\text{acc}}`$ is an increasing function of $`a`$ for fixed $`b`$ and $`T`$.

Also shown are the corresponding results for the pivot update, where only one angle $`\varphi _i`$ is turned ($`S=|\delta \varphi _i|`$ if the change is accepted). At the higher temperature, we find that the $`b=b_{\mathrm{max}}`$ and $`b=0`$ updates show similar behaviors. The pivot update is somewhat better and has its maximum $`S`$ at low $`P_{\text{acc}}`$, where the proposed change $`\delta \varphi _i`$ is drawn from the uniform distribution between 0 and $`2\pi `$. This is consistent with the finding  that the pivot update is a very efficient method for self-avoiding walks, in spite of a low acceptance rate. The situation is different at the lower temperature, which is much harder to simulate. Here, the $`b=b_{\mathrm{max}}`$ update is the best. The maximum $`S`$ is approximately three times higher for this method than for the other two. This shows that the biasing probability Eq. 14 is indeed useful in the folded phase.

The $`b=0`$ update can be compared with the moves used by Shimada et al.  in a recent all-atom study of kinetics and thermodynamics for the protein crambin with 46 amino acids. These authors updated sets of two, four or six backbone torsion angles, using independent Gaussian steps with a standard deviation of $`2^{}`$. Our $`b=0`$ update has maximum $`S`$ at $`a6400(\mathrm{rad})^2`$ for $`kT=0.6`$, which corresponds to a standard deviation of $`0.7^{}`$. This value is in line with that used by Shimada et al., since we turn eight angles.

How local is the method for $`b=b_{\mathrm{max}}`$? To get an idea of that, we calculated the distribution of $`\mathrm{\Delta }^2`$ (see Eq. 11) for accepted moves, for $`b=b_{\mathrm{max}}`$ and $`b=0`$ at $`kT=0.6`$. As was previously the case, we restricted ourselves to angles in the middle helix. The two distributions are shown in Fig. 3 and we see that the one corresponding to $`b=b_{\mathrm{max}}`$ is sharply peaked near $`\mathrm{\Delta }^2=0`$. This shows that the $`b=b_{\mathrm{max}}`$ update is much more local than the unbiased $`b=0`$ update, although the average step size, $`S`$, is considerably larger for $`b=b_{\mathrm{max}}`$.

So far, we have discussed static (one-step) properties of the updates. We also estimated the dynamic autocorrelation function

$$C_i(t)=\frac{\mathrm{cos}\varphi _i(t)\mathrm{cos}\varphi _i(0)\mathrm{cos}\varphi _i(0)^2}{\mathrm{cos}^2\varphi _i(0)\mathrm{cos}\varphi _i(0)^2}$$
(22)
for different $`\varphi _i`$. This measurement is statistically very difficult at low temperatures. However, the sixteen most central angles $`\varphi _i`$ in the sequence, all belonging to the middle helix, were found to be effectively frozen at $`kT=0.6`$, and the time scale for the small fluctuations of these angles about their mean values was possible to estimate. In Fig. 4, we show the average $`C_i(t)`$ for these sixteen angles, denoted by $`C(t)`$, against Monte Carlo time $`t`$, for the $`b=0`$, $`b=b_{\mathrm{max}}`$ and pivot updates. One time unit corresponds to one elementary move, accepted or rejected, at a random position along the chain. We see that $`C(t)`$ decays most rapidly for the $`b=b_{\mathrm{max}}`$ update. So, the larger step size of this update does make the exploration of these degrees of freedom more efficient.

Let us finally comment on our choice to work with $`n=8`$ angles. This number can be easily altered and some calculations were done with $`n=6`$ and $`n=7`$, too. For $`n=6`$, the performance was worse, which is not unexpected because there are no soft modes available; there are not more variables than constraints. The results obtained for $`n=7`$ were, by contrast, comparable to or slightly better than the $`n=8`$ results.

## 5 Discussion

Straightforward Monte Carlo updates of torsional degrees of freedom tend to cause large changes in the global structure of the chains unless the step size is made very small, which is a problem in simulations of dense polymer systems. The strictly local concerted-rotation approach provides a solution to this problem but is rather complicated to implement. In this paper, we have discussed a method that may be less powerful but is much easier to implement, which suppresses rather than eliminates nonlocal deformations.

The method is flexible and not much harder to implement than simple unbiased updates. However, compared to such updates, it has two distinct advantages: the step size can be increased and the update becomes more local, as shown by our simulations of the three-helix-bundle protein in its folded phase.

Making the update more local is important in order to be able to increase the step size and thereby improve the efficiency. At the same time, it makes the dynamics more realistic; the proposed method is, in contrast to the other methods mentioned, tailored to avoid drastic deformations both locally and globally. Therefore, although this paper was focused on thermodynamic simulations, it should be noted that this method may be useful for kinetic studies, too.

### Acknowledgments

This work was supported in part by the Swedish Foundation for Strategic Research. G.F. acknowledges support from Universitá degli studi di Cagliari and the EU European Social Fund.