text:
# The entrance dynamics of polymers into a gel

## Abstract

We use computer simulations to study the (driven) dynamics of a charged polymer in a disordered medium, thus mimicking the setting used in gel electrophoresis. In agreement with experiments, we find that inside the gel the mobility of the polymer is only a weak function of its length $`N`$. In contrast to this, the mean entrance time into the gel, is a very strong function of $`N`$, $`\tau _eN^{1.3}`$, and does not show any sign of saturation with increasing $`N`$. We rationalize this effect by means of a simple model and propose an experimental setup that should allow to increase the separation ability of gel electrophoresis significantly.

PACS numbers: 36.20.E, 87.15.T, 61.20.J, 36.20.E

Understanding the dynamics of macromolecules in disordered porous media is important in many situations of practical interest such as filtration, oil recovery, transport through membranes, and gel electrophoresis , the problem we treat in this letter. Gel electrophoresis is a simple but powerful technique used to separate polyelectrolyte strands with different length. One of its main fields of application is molecular biology, where it is used to sequence DNA. In its standard setting, fragments of various lengths are injected in a thin layer of gel, composed by, e.g., agarose powder and a buffer solution. The latter imparts a charge to each base-pair. If an external electric field is applied, the fragments migrate in its direction with shorter strands moving faster than longer ones. This, in principle, allows one to separate polymers with different length $`N`$. However, for large $`N`$ the velocity depends only weakly on $`N`$ and therefore it is not feasible to use gel electrophoresis as an efficient separation technique for long molecules.

In the usual experimental setup, the molecules are driven through the interface between the buffer solution and the gel by the external field. Very recently this process has called the attention of experimentalists who tried to determine if it can be used as a method to improve the efficiency of separation.

Using fluorescence microscopy, Oana et al.  concluded that DNA molecules stay at the interface for a time that decreases with increasing pore size. Unfortunately, due to the difficulty with this technique to investigate many such entrance events, it is hard to extract quantitative information, such as the dependence of the delay time on the size of the molecules or the applied field. In addition, the limited spatial resolution of $`0.1\mu m`$ prevents a detailed knowledge of what the polyelectrolytes actually do at the interface.

Klepárník et al.  investigated the injection process by measuring the time that a polymer needs to propagate a macroscopic distance. From its extrapolation to vanishing distance, they inferred that the molecules are stuck at the interface during an interval that increases with length and speculated that it will grow with decreasing field. They attributed the delay to an electro-osmotic flow in the opposite direction of the field which hinders some of the molecules from entering the gel.

In the past, numerical simulations have proven to be a very useful tool to understand the motion of polymers in disordered media . In the context of electrophoresis, some of the most significant results are the ones by Deutsch et al.  who solved numerically the Langevin dynamics of a driven chain moving in a two dimensional random array of obstacles, allowing, however, the self-intersection of the polymer. For rather strong fields, the numerical results showed the existence of U-like conformations inside the bulk of the gel, a result which was later confirmed experimentally by video microscopy .

Here we study the details of the entrance process of charged polymers into the gel by means of simulations. In particular we analyze how the time it takes the polymer to enter depends on its size and the strength of the field. Since simulations are not restricted by any resolution limit, they allow us to identify the conformations of the molecule during the entrance in the gel, which in turn will help us to understand the details of this process. By using the gathered information, we propose a setup that should significantly improve the separation-efficiency in constant field electrophoresis.

The polymer is modeled by a chain of $`N`$ beads connected to each other by an unbreakable spring . The interaction between two monomers $`i`$ and $`j`$, that are a distance $`r_{ij}`$ apart, is

$`V(r_{ij})`$ $`=`$ $`4ϵ[(\sigma /r_{ij})^{12}(\sigma /r_{ij})^6]K_0R_0^2/2\mathrm{ln}[1(r_{ij}/R_0)^2](\delta _{j,i+1}+\delta _{j,i1}),`$ (1)
where $`\delta _{i,j}`$ is the Kronecker $`\delta `$. The first term is strongly repulsive at short distances and represents the excluded volume. For computational efficiency it is truncated and shifted at the location of its minimum, $`r_m=2^{1/6}\sigma `$. The second term represents the bond between neighboring monomers. In the following we use reduced units and measure length, energy and time in units of $`\sigma `$, $`ϵ`$, and $`\sqrt{\sigma m/ϵ}`$, where $`m`$ is the mass of one monomer and the Boltzmann constant $`k_B`$ is set to 1. The constants $`K_0`$ and $`R_0`$ are 30 and 1.5, respectively. We mimic the charging effect by assigning a charge $`q=1`$ to each monomer. (Note that we neglect the Coulomb interaction between monomers since it can be assumed that the solvent leads to strong screening.) If we assume that the polymer is a DNA strand with charge $`150e150\times \mathrm{1.6\hspace{0.17em}10}^{19}`$C per base pair, these numbers translate into physical units of around 200Å for the unit of length and 2500 V/cm for the unit of field. The gel is modeled by a disordered 2d array of point particles obtained by equilibrating a system of soft disks interacting via a $`r^{12}`$ potential at $`T=1`$ . A snapshot of this system was expanded by a factor of 5 and the resulting, permanently fixed, configuration with a liquid-like structure at short range, is our 2d gel. The typical distance between neighboring particles is $`7`$ which gives pore sizes of the order of $`515`$. The size of the gel perpendicular to the field is 160, and is sufficiently large to avoid finite size effects if periodic boundary conditions are used. The interaction of the polymer with the gel is given by a truncated and shifted Lennard-Jones potential with the form of the first term in Eq. (1). The solvent imparts a stochastic force on each monomer, that is mimicked by substituting after each time step the velocities of the monomers by new ones drawn from a Maxwell-distribution at $`T=1`$. The positions of the monomers are then updated by means of the velocity version of the Verlet algorithm, using a time-step $`=0.01`$.

At the beginning of each simulation we equilibrate the polymer by taking into account only the monomer-monomer interactions. Typically, the polymer acquires an open random-walk like conformation with a Kuhn length $`b25`$ (estimated from visual inspection), which is of the order of the pore size and corresponds to $`4001000`$Å. Once the polymer is equilibrated we slide it in the direction of the field until it is at a distance $`r_m`$ from the closest particle in the gel. This instant defines the initial time $`t=0`$ when the polymer starts interacting with the gel. We then turn on the field, driving the polymer into the gel.

We have studied polymer sizes $`N=25`$, 50, 100, 200, and 400 and fields $`E=0.0312`$, 0.0625, 0.125, 0.25, 0.5, and 1.0. In all runs the temperature is such that $`k_BT=1`$, which we identify with room temperature $`T=300K`$. This leads to an adimensional parameter $`qEb/k_BT`$ between $`0.1`$ and $`10`$. To improve the statistics we average all results over at least 30 independent runs.

One possibility to characterize the dynamics inside the gel is to monitor the time dependence of the averaged center-of-mass position in the direction of the applied field, $`y_{cm}(t)`$, shown in Fig. 1 for various $`N`$. After a transient, the curves become straight lines with slope $`v(t)`$. In Fig. 2 we plot the mobility at large times, $`\mu _{\mathrm{}}=lim_t\mathrm{}v(t)/E`$, as a function of $`N`$. We see that to a good approximation $`\mu _{\mathrm{}}`$ is independent of $`E`$ and that it is only a weak function of $`N`$. (A change of $`N`$ by a factor of 100 changes $`\mu _{\mathrm{}}`$ only by a factor of less than two.) The reason for this is that for long polymers not only the force due to the field is proportional to $`N`$ but also the effective friction due to the gel is linear in $`N`$. For large $`N`$ the two forces cancel and the resulting velocity has only a weak $`N`$-dependence . Therefore it is very hard to separate strands with similar lengths.

From Fig. 1 we also recognize that at short times the mean velocities have a stronger $`N`$-dependence than at long times, when the polymer is in the bulk. In Fig. 3 we show a series of snapshots that help us understanding this result. At early times, Fig. 3a, the field drives the (relatively open) equilibrium configuration in an unhindered fall towards the gel. The polymer hits the gel and its mobility decreases slightly. After a short interval some monomers, either at the end or in the middle of the polymer, find their way into the pores while others are still pinned at the boundary, Fig. 3b. The configuration acquires little hernias that penetrate the gel, Fig. 3c. Subsequently, the little hernias are eaten up by the longer ones  and a “staple-like” configuration is reached with two elongated arms that penetrate the gel, Fig. 3d . For much longer times the longer arm will finally drag the polymer into the gel. The transient dynamics due to the entering process has ended and the asymptotic mean velocity is reached. For large $`N`$, the time it takes to resolve this staple configuration is determined by the difference in length between the arms. Since this time is significantly longer that the time it takes the polymer to reach this configuration, it determines the entrance delay almost completely.

We mention that such staple-like configurations also occur inside the gel, though with two important differences. i) Unless the field is extremely weak, their width is generally smaller, since the pinning is due to only one or two obstacles. In the entering case instead, the pinning occurs at more sites (therefore we call the configuration “staple-like” and not hernias) since outside the gel the polymer is open with a typical cross-section that scales like the radius of gyration. ii) Inside the gel the generation of such pinned configurations is a non-synchronized stochastic process. At the interface most strands are trapped almost immediately after touching the gel and hence the trapping is synchronized. These two points explain why the motion of the polymers close to the surface is so much slower than inside the gel. We also mention that none of the polymers we generated bounced back into the solvent when they encountered the interface, i.e. moved against the direction of the field. This shows that there can be a significant delay also without an electro-osmotic flow, in contrast to the explanation put forward by Klepárník et al.  of their experimental results.

In order to analyze how the entrance time depends on length and field, we define a time $`\tau _e`$ such that $`y_{cm}(\tau _e)=\alpha N`$, with $`\alpha =0.75`$. (The precise value of $`\alpha `$ is irrelevant as long as the probability for the polymer to have entered more than $`\alpha N`$ without having been entangled is very small.) In Fig. 4 we show $`\tau _e`$ multiplied by $`E`$ versus $`N`$ for different values of $`E`$. From this plot we conclude that the entering time is inversely proportional to the field. Much more important, however, is the strong dependence of $`\tau _e`$ on $`N`$ (note the logarithmic scale on the $`\tau _e`$ axis!). This dependence does not saturate with $`N`$ and it is approximately given by the law $`N^{1.3}`$ (bold solid line). Therefore it should be possible to use this entering process as a tool to increase the selectivity of gel electrophoresis.

A very simple model allows us to understand on a qualitative basis the origin of this strong $`N`$-dependence: Let us assume, as suggested by the simulations, that we can neglect the time needed to reach the staple-like conformation and hence can approximate $`\tau _e`$ by the time necessary to resolve it. The problem can now be modelled by a massive rope with uniform linear density $`q`$ that hangs under gravity, with strength $`E`$, from a pulley and whose motion is slowed down due to friction that is proportional to its velocity. The arm carrying more mass makes the rope slide in its direction. If we call $`\delta `$ the difference in length of the two arms, normalized by the total length of the rope, it is easy to calculate the time needed for the rope to fall from the pulley and one finds

$$\tau _e=\frac{\xi N}{2qE}\mathrm{ln}\delta ,$$
(2)
where $`\xi `$ is the friction constant per unit length. In our problem $`\delta `$ fluctuates from sample to sample. Therefore an average over samples translates into an average over $`\delta `$. If we assume that the distribution of $`\delta `$ is independent of $`\delta `$ we obtain $`\tau _e\xi N/E`$. The stronger $`N`$-dependence shown in Fig. 4 can be accounted for by assuming that the distribution of $`\delta `$ has a somewhat larger weight at small $`\delta `$.

Equipped with the results of the simulation we propose an experimental setup that should allow to increase the ability of gel electrophoresis to separate polymers of different length. In this setup we have an alternating sequence of strips made out of gel and strips with just solvent. If the field is pointing upwards, these strips are arranged vertically and have “heights” $`g_1`$, $`s_1`$, $`g_2`$, $`s_2`$, $`g_3`$, etc. The basic idea is as follows. If the polymer mixture is placed below the first strip all polymers enter this strip, leading to the very strong separation studied here. The polymers later leave the first strip and enter a strip of solvent. The polymer changes its conformation from the elongated structure acquired in the gel to a random walk-like structure since in the pure solvent the elongated structure is not stable. After some time it hits the next solvent-gel interface and the process starts again. In order to allow a strong separability at the interface, the width of the gel has to be at least on the order of the longest polymers. To ensure that long polymers open up again during their fall in the solvent, this strip also needs to be sufficiently long. Therefore it is advisable to increase progressively the width of the strips i.e. $`g_1<g_2<g_3<\mathrm{}`$ and $`s_1<s_2<s_3\mathrm{}`$.

We conclude by stating that these results, together with the ones in Refs. , make us expect that the use of interfaces can lead to a strong improvement of the separability in gel electrophoresis.

We thank S. Bottani, M. Picco, and J.-L. Viovy for very useful discussions. Part of this work has been supported by the DFG through grant SFB 262, the ITP in Santa Barbara, and the SPHINX collaboration of the ESF.


format: text

text:
# DIFFRACTIVE STRUCTURE FUNCTION AT VERY SMALL 𝛽 AND UNITARITY CORRECTIONS IN THE COLOR DIPOLE REGGE APPROACH

## References


text:
# From the quantum Zeno to the inverse quantum Zeno effect

## Acknowledgements

We thank L.S. Schulman for interesting comments.


text:
# The IRAS view of the Local Universe

## 1. Introduction

Redshift surveys of galaxies can be used to infer the dynamics of the Local Universe. Comparison with observed peculiar velocities allows determination of the mass associated with galaxies (usually parameterised by the the quantity $`\beta \mathrm{\Omega }^0.6/b`$), and also reconstruction of the initial conditions. Because gravity is extensive, accurate predictions formally require uniform surveys with infinite depth and complete sky coverage; failing that, we need a depth much larger than the region of interest, and with coverage gaps small enough to interpolate over. IRAS surveyed 96% of the sky with negligible extinction, tracing the large-scale structure to a depth of at least $`20,000\mathrm{km}\mathrm{s}^1`$. For these reasons, IRAS-based surveys have dominated this field for the last decade. The IRAS dipole has remained astonishingly stable in direction over this period, despite the ever-increasing depth and sky coverage of the follow-up redshift surveys. However, the amplitude has increased continuously, with extra contributions uncovered at larger distances and lower latitudes.

The Point Source Calatog Redshift Survey (PSCz) consists of essentially all galaxies in the PSC, down to its nominal flux limit of $`0.6\mathrm{Jy}`$, across the 84% of the sky where extinction was estimated to be low enough ($`A_B<2^m`$) to allow reasonably complete optical identifications and spectroscopy with 2-metre class telescopes. The PSCz data is now available from http://www-astro.physics.ox.ac.uk/$``$wjs/pscz.html; also there available is the full, labelled 3D density distribution, using the Fourier interpolation method presented in Saunders and Ballinger (2000). The catalogue and main science papers to date are summarised in Saunders et al. (2000a).

The PSC itself is useable much deeper into the Plane. The Behind The Plane Survey extends the PSCz to cover the entire 93% of the sky with useable PSC data, and is now almost completed. The remaining size of the Zone of Avoidance is just 4% of the sky, and this is hopefully small enough to interpolate over with confidence. The sky distribution for both surveys is shown in Figure 1.

## 2. The PSCz velocity field and the Path Interchange Zeldovich Approximation

A linear theory reconstruction of the velocity field and comparison with the Mark III dataset was been performed by Branchini et al. (1999). However, the heavy smoothing needed to make linear theory applicable seriously compromises the detail with which such comparisons can be done.

The PIZA algorithm, originally presented by Croft and Gaztañaga (1997), is a quasi-linear method for reconstructing the velocity field from galaxy surveys, allowing the full resolution of the data to be kept. The basic idea is wonderful and very simple: to Zeldovich order, velocities and trajectories depend in a fixed way on the scale factor; so minimising the action is equivalent to minimising the rms trajectories of the galaxies. Given a set of uniform initial positions and the current galaxy positions, the least-action solution is found by setting up arbitrary trajectories and swapping pairs to reduce the action. The formalism can be readily adapted to minimise (for a given $`\beta `$) the real-space trajectories from redshift-space information. With more complexity, the changes in selection function for individual galaxies caused by peculiar velocities can also be included, though we have found doing this leads to instabilities without any improvement in accuracy, and this is not done here. This and edge effects lead to a systematic velocity underestimation of $`20\%`$ in comparison with simulations. Further details are presented in Valentine, Saunders and Taylor (2000). The resulting velocity field in the Supergalactic Plane is shown in Figure 2. The flow field is dominated by infall into the Supergalactic Plane, and coherent flow along it involving the Local Supercluster, Centaurus and the Shapley Concentration.

## 3. Comparison with observed peculiar velocities

The velocity field predicted from PIZA can be compared with observations. The derived trajectory of the Local Group is robust, with weak dependence on how the mask is filled or the input $`\beta `$. Comparison of this trajectory with the CMB dipole gives $`\beta =0.53\pm 0.1`$. A comparison of the bulk flow as a function of radius with the MarkIII dataset implies a value for $`\beta `$ of $`0.5\pm 0.15`$.

We have also undertaken a new comparison with the SFI dataset of peculiar velocities of Spiral Galaxies (Giovanelli et al. 1998). We have calculated the average peculiar velocity with respect to the Local Group, of galaxies in redshift slices, so as to be as closely as possible comparable to Giovanelli et al. . We find reasonable agreement in amplitude, though there are significant directional differences, for $`\beta =0.55\pm 0.1`$ (Table 1). We agree with Giovanelli et al. that most of the dipole is generated within $`5000\mathrm{km}\mathrm{s}^1`$. However, we do find substantial contribution to the dipole beyond this distance (see below), even though our derived velocity dipole at this distance is in excellent agreement with the CMB. The sum derived by Giovanelli et al. is a galaxy-weighted vector sum, rather than volume-weighted. Therefore, a good agreement between velocity vector and CMB dipole does not necessarily imply that a cosmic rest frame has been reached.

## 4. PSCz dipole

In the course of the BTP survey, we also took several hundred new reshifts for low-latitude PSCz galaxies, giving an overall redshift incompleteness for bright galaxies of 70 galaxies (0.5%). A further 230 galaxies (1.5%) without redshifts are assumed from their optical and/or infrared appearance to be beyond $`30,000\mathrm{km}\mathrm{s}^1`$. In the original, publicly available survey, there was known incompleteness at lower redshifts towards the Galactic anti-centre as discussed by Saunders et al. (2000a) and whose effect on the dipole is discussed in Rowan-Robinson et al. (2000). It is appropriate then to recalculate the dipole analysed by Rowan-Robinson et al. . We have corrected for peculiar velocities using PIZA as above, with $`\beta =0.5`$ for self-consistency between gravity and velocity dipoles. The mask is filled in via the Fourier interpolation method of Saunders and Ballinger (2000). The revised dipole is presented in Figure 3. The dipole is remarkably stable between $`125180h^1\mathrm{Mpc}`$, However, it shows a significant shift at $`180220h^1\mathrm{Mpc}`$ in the X and Y components. This shift is seen in the analysis of Rowan-Robinson et al. , and is too large to be attributed to the small remaining redshift incompleteness. It remains when the interpolated galaxies in the Plane are replaced by randomly distributed galaxies at the mean density, so it is not due to the interpolation either. It does not appear to be associated with any single particular structure.

In linear theory, the contribution to the velocity dipole per galaxy is

$$𝐕(𝐫)=\frac{\beta H_0}{4\pi }\frac{𝐫}{r^3\psi (𝐫)}$$
(1)
where $`\psi `$ is the selection function. In Figure 4, we plot the amplitude of this velocity as a function of distance for the PSCz. The dipole contribution per galaxy is remarkably constant in the range $`20100h^1\mathrm{Mpc}`$. Nearby, the contribution is large because of the $`1/r^2`$ term, but at least the numbers of galaxies involved are small and the structures they reside in are real and well-sampled. At large distances, the contribution blows up because $`\psi `$ falls off much faster than $`1/r^2`$, and the dipole becomes increasingly dominated by shot noise (e.g. Taylor and Valentine 1999). Strauss et al. (1992) derived an optimal weighting scheme, in the sense of producing the minimum variance ensemble difference between true and estimated gravity vectors. At large distances, this amounts to giving each galaxy a ‘$`J_3`$ weighting’ i.e.

$$𝐕=\frac{\beta H_0}{4\pi }\underset{gals}{}\frac{J_3(r)\psi (𝐫)}{1+J_3(r)\psi (𝐫)}\frac{𝐫}{r^3\psi (r)}$$
(2)
Where $`J_3(r)=4\pi _0^rr^2\xi (r^{})𝑑r^{}`$ is the integral over the correlation function. $`J_3(r)`$ is not well-determined on large scales, but the precise value used is not critical as any error only means the weighting is quite optimal, while leaving the dipole estimate unbiased in an ensemble sense. Based on the PSCz power spectrum (Sutherland et al. 1999), we have simply assigned $`J_3(r)`$ a constant value $`J_3=20,000h^1\mathrm{Mpc}`$ for all $`r`$. The resulting velocity dipole per galaxy is shown in Figure 4; there is now a sharp attenuation above $`200h^1\mathrm{Mpc}`$.

Based on Lineweaver et al. (1996) and Yahil, Sandage and Tamman (1977), we have assumed a value for the LG motion with respect to the CMB of $`(276,30,627\mathrm{km}\mathrm{s}^1)`$. The misalignment of CMB and IRAS dipoles remains $`1525^{}`$, depending on chosen upper cutoff radius for the gravity sum, precise treatment of the mask etc. The inferred value of $`\beta `$ is $`\beta =0.54\pm 0.10`$ from the overall amplitude; however, the Z-component, which is very stable and also less sensitive to the inferred structure behind the Plane, gives $`\beta =0.40\pm 0.07`$.

The dipole and peculiar velocity results are largely independent, since the latter is sensitive only to tidal effects from large scales. The completely independent analysis of Tadros et al. (1999), using the statistical distortion of the clustering in redshift space, also gives a consistent, low-$`\beta `$ result $`\beta =0.47\pm 0.16`$.

## 5. The Behind The Plane survey

The 84% sky coverage of the PSCz survey is effectively limited by the need to get, for every galaxy, an optical identification from sky survey plates. The fractional incompleteness in sky coverage translates directly into a much larger uncertainty in predicting the gravity dipole on the Local Group: suppose the fractional masked area is $`f`$, and also that this is large and compact enough that we have essentially no information about the distribution behind the mask. Then we can think of the sky as consisting of $`1/f`$ independent contributions to the dipole, one of which is unknown, giving a fractional uncertainty in the measured dipole vector of $`\sqrt{f/(1f)}`$. Thus even 10% missing sky coverage can lead to a error of 1/3 in the dipole. The bulk of the dipole is generated around $`40h^1\mathrm{Mpc}`$, at which distance the correlation length $`r_0`$ subtends just $`5^{}`$, so gaps much larger than this are too large to interpolate usefully and lead to large errors in the dipole.

The IRAS PSC data itself is reliable to much lower latitudes, although genuine galaxies are outnumbered by Galactic sources with similar IRAS properties. Previous attempts to go further into the Plane have either been restricted to the Arecibo declination range, or have relied on optical identifications from Sky Survey Plates. Because the extinction may be several magnitudes or more, they have inevitably suffered from progressive and unquantifiable incompleteness as a function of latitude. Hence in 1994 we embarked on a program, parallel with the PSCz survey, to systematically identify and get redshifts for low latitude IRAS galaxies wherever the PSC data allowed - essentially limited only by IRAS coverage and confusion, to 93% of the sky. A detailed description is given in the associated poster paper (Saunders et al. 2000b). The identification program is complete; and we have found over 200 galaxies which were not included in the 1.2Jy survey but should have been. The BTP sky distribution is shown along with the PSCz in Figure 1. At low latitudes towards the Galactic centre, the PSC is not complete to our nominal flux limit of $`0.6\mathrm{Jy}`$, and the density of identified BTP galaxies reflect this. This incompleteness is interpreted, quantified and corrected for in the poster paper. Notwithstanding this, the huge overdensity of galaxies in the Great Attractor region $`(325,\pm 5)`$ is obvious; as is the Puppis supercluster $`(250,0)`$ and the extension from Perseus (150,-10) and Pegasus (95,-10) into the Plane.

The overall surface density of IRAS galaxies in the BTP area, when corrected for known incompleteness, is 78% of the PSCz average. Unquantified incompletenesses could conceivably push this up to the average value, but clearly we can rule out large overall overdensity behind the Plane. We have used the Fourier interpolation technique to find the projected surface density of the PSCz+BTP galaxy distribution; this is shown in Figure 5. In addition to the above features, we see the Ophiucus Spercluster clearly at (0,10) and the local void behind the Galactic Centre.

## 6. 3D distribution

The spectroscopy is now completed in the south, and useably completed in the north. In addition to BTP galaxies, we also systematicall pursued faint, low-latitude PSCz galaxies where there was any likelihood that these were within $`30,000\mathrm{km}\mathrm{s}^1`$. There are c.370 very faint galaxies in the combined PSCz+BTP with unknown redshift, but these are likely beyond $`30,000\mathrm{km}\mathrm{s}^1`$. There are a further 240 for which we would like a redshift but do not currently have one, though for half of these we have unreduced data. As a temporary fix, for these latter 240 galaxies we have simply assigned the redshift of a nearby galaxy with redshift. The N(z) distribution, both with and without this fix, is shown in Figure 6, together with the prediction from the PSCz survey. There is some evidence that we are missing galaxies at most or all redshifts, but also that there is no cutoff in the redshift completeness with distance till at least $`25,000\mathrm{km}\mathrm{s}^1`$ \- much better than our original target of $`15,000\mathrm{km}\mathrm{s}^1`$.

We can then make a 3-dimensional map of the galaxy distribution. The full, animated 3D map is available from http://www.roe.ac.uk/willwww. Figure 7 shows the peak we find in the GA region at $`(325,5,3500\mathrm{km}\mathrm{s}^1)`$. The overdensity is quite modest compared with other superclusters in the region, and the peak overdensity is closer than the velocity of the Norma cluster at $`4844\mathrm{km}\mathrm{s}^1`$ (Woudt et al. 2000). However, we find the GA to be very extended in redshift space, from $`30005500\mathrm{km}\mathrm{s}^1`$. More surprising is Figure 8, showing an apparent ‘Greater Attractor’, directly behind the GA at $`(326,3,12500\mathrm{km}\mathrm{s}^1)`$. Because of the incompleteness weighting, and because the overdensity is right on the edge on the mask, the finding represents only a dozen measured redshifts and is at this stage preliminary.

## 7. The PSCz+BTP Dipole

The combined PSCz+BTP sample is now complete enough to attempt a revised dipole. One of the main motivations for this is that current low $`\mathrm{\Omega }`$ models predict small misalignments between gravity and velocity vectors, especially as the Hubble flow in the immediate vicinity of the Local Group (where non-linear effects are most likely to arise) is so cold: a misalignment as large as the $`2025^{}`$ between IRAS and CMB dipoles would present problems for current models. As above, we have assigned a nearest neighbour redshift to 240 galaxies. We have used our interpolated 3D maps to generate random ‘galaxies’ to add to the catalogue when determining the dipole. However, because the missing region is so much smaller, the results are now much less sensitive to this interpolation.

To get real-space positions, we used the PIZA reconstruction for the PSCz galaxies as above, with an assumed $`\beta =0.5`$. For the BTP sources we just used the redshift-space position. The resulting dipole is shown in Figure 9. Overall, the dipole looks much the same as for the PSCz alone. However, the misalignment with the CMB is now very much smaller, being $`10^{}`$ at $`250h^1\mathrm{Mpc}`$ and $`5^{}`$ at $`300h^1\mathrm{Mpc}`$, with the shot noise error being $`10^{}`$. The misalignment is $`13^{}`$ for the $`J_3`$-weighted sum. The directional agreement also appears to be robust: if the galaxies with unknown redshift are dropped altogether, the dipole only changes by $`5^{}`$ in direction and 7% in amplitude, while if random sources are put uniformly in the mask at the average density instead of interpolated, the dipole only alters by a few degrees in direction and 10% in amplitude.

The shift at $`180220h^1\mathrm{Mpc}`$ seen in the PSCz remains. At this stage, it remains possible that variable redshift completeness is responsible. Also, at this distance, very slight systematic effects in the PSC itself could have large effects on the dipole, but these would not be expected to be restricted to a particular velocity range as seen in Figure 8. Beyond that, there is a drift in the X component, but this too is sensitive to variable redshift incompleteness around the Galactic Plane and is not secure at this stage.

Peacock (1992) noted that apparent convergeance in the dipole may be illusory. The apparent convergeance of the IRAS dipole between $`80180h^1\mathrm{Mpc}`$, in a direction only $`20^{}`$ away from the CMB, was seen in the QDOT, 1.2Jy and PSCz surveys, and was in each case interpreted as the ‘true’ value for the determination of $`\beta `$, with no significant contribution from larger distances. We now see that there may indeed be a significant contribution from beyond $`180h^1\mathrm{Mpc}`$. Of course the same argument can in principal be applied this work, in that we have not ruled out contributions from beyond $`300h^1\mathrm{Mpc}`$. However, we can hope that the excellent directional agreement we find, combined with the implausible size of structures required, make this unlikely.

## References

Branchini E. et al. 1999. MNRAS 308 1.

Croft,R.A.C. & Gaztañaga,E., 1997. MNRAS , 285 793.

Giovanelli,R., Haynes,M.P., Freudling,W., da Costa,L.N., Salzer,J.J., Wegner,G. 1998. ApJ 505 L91.

Lineweaver,C.H., Tenorio,L., Smoot,G.F., Keegstra,P., Banday,A.J., Lubin,P. 1996. ApJ 470 38.

Peacock,J.A. 1992. MNRAS 258 581.

Rowan-Robinson,M. et al. 2000. MNRAS in press (astro-ph/9912223).

Strauss,M., Yahil,A., Davis,M., Huchra,J.P., Fisher,K. 1992. ApJ 397 395.

Saunders,W. and Ballinger 2000. These proceedings.

Saunders,W., Sutherland,W.J., Maddox,S.J., Keeble,O., Oliver,S.J., Rowan-Robinson,M., McMahon,R.G., Efstathiou,G.P., Tadros,H., White,S.D.M., Frenk,C.S., Carramiñana,A., Hawkins,M.R.S. 2000a. MNRAS , in press (astro-ph/0001117).

Saunders,W. et al. 2000b. These proceedings.

Sutherland W.J.et al. 1999. MNRAS 308 289.

Tadros H. et al. 1999. MNRAS 305 527.

Valentine,H.E.M. and Taylor,A.N. 1999. MNRAS 306 491.

Valentine,H.E.M., Saunders,W., and Taylor,A.N. 2000. MNRAS , submitted.

Woudt,P.A., Kraan-Korteweg,R.C., Fairall,A.P. These proceedings.

Yahil,A., Tammann,G., Sandage,A. 1977. ApJ 217 903.


text:
# X-ray spectral features from GRBs: Predictions of progenitor models

## 1 Introduction

With the advent of the new generation of X-ray telescopes, such as Chandra and XMM-Newton, the detection of X-ray spectral signatures from the environments of cosmological gamma-ray bursts (GRBs) has become a realistic prospect. The marginal detection of a redshifted Fe K$`\alpha `$ emission line in the afterglow of GRB 970508 (Piro et al. (1999)) with the BeppoSAX NFI has stimulated a vital discussion about the possible origin of this line feature (Ghisellini et al. (1999); Lazzati at al. (1999); Böttcher et al. 1999a ; Böttcher, Dermer, & Liang 1999b; Vietri et al. (1999); Weth et al. (2000); Paerels et al. (2000); Böttcher (2000)). So far, this discussion has concentrated on determining the required / inferred physical conditions in the vicinity of GRB 970508 needed to explain the iron-line feature in the afterglow (assuming this feature is real).

An essential assumption in these papers was that the material responsible for potential emission line features is illuminated (and partially or completely ionized) by burst emission which is qualitatively similar to the observed GRB and afterglow radiation. Also, these papers assume that the blast wave which interacts with the line-emitting material is relativistic, and qualitatively similar to the one which is associated with the GRB and its afterglow. The conditions inferred under these assumptions are rather extreme, requiring large amounts ($`10^4M_{}`$) of iron to be concentrated anisotropically in a small ($`R10^3`$ pc) region around the central engine of the GRB (e.g., Lazzati at al. (1999), Weth et al. (2000), Böttcher (2000)). Paerels et al. (2000) have recently suggested that high-resolution X-ray spectroscopy may help to distinguish between different line production mechanisms as a way to unveal the underlying physical scenario. They argue that, in particular, photoionization can be distinguished from collisional excitation of the line by virtue of the lower plasma temperature resulting in photoionization scenarios. In their detailed analysis of the BeppoSAX MECS spectrum of GRB 970508 during the time segment in which the $`3.6`$ keV emission line was marginally detected, they found that the line width and apparent strength of the radiative recombination continuum blueward of the emission line, in combination with the measured redshift of the GRB at $`z=0.835`$, are inconsistent with a photoionization scenario in which the source of the line emission is in photoionization equilibrium with the afterglow radiation.

In addition to this evidence, most of the currently popular gamma-ray burst models may be hard-pressed to produce the conditions inferred from photoionization-based scenarios to explain the observed iron-line feature in GRB 970508. These models are not ruled out, both because the iron line feature in GRB 970508 is marginal, and even if it is real, it may appear only in one small subclass of GRBs. However, Vietri et al. (1999) have used this line to argue for their alternative supranova model (Stella & Vietri (1998)) which invokes the collapse of a neutron star as it spins down. Unfortunately for this model, realistic calculations of collapsing neutron stars (Ruffert, Janka, & Schäfer 1996; Fryer & Woosley 1998) find that these collapses eject too much baryonic material and have too little energy to produce GRBs. Because these models do not include magnetic fields, this does not rule out the supranova model. Even so, the supranova model would most likely produce a short burst, and hence can not explain the line feature seen in GRB 970508, since that burst was a long-duration burst with $`t_\gamma 25`$ s.

Instead of trying to explain GRB 970508, in this paper we look at the more viable “black-hole accretion disk” class of GRB engines (see Fryer, Woosley, & Hartmann 1999 for a review) and estimate potential X-ray spectral features in GRB afterglows that these models produce. We make the initial assumption that these features are produced by the collision of the GRB ejecta with the environment produced during the formation process of the GRB engine. In the Collapsar and he-merger engines, the GRB is surrounded by a disk of material which is produced during their progenitor evolution. MacFadyen & Woosley (1999) found that collapsar (and probably he-merger) GRBs are beamed, forming a jet along the angular momentum axis of the accreting disk. They also found that the explosion along the equator is likely to be very baryon loaded ($`>`$1 M) moving at velocities much less than the speed of light (less than $`10^9`$ cm s-1). In §2, we discuss the formation of these nonrelativistic disks around black-hole accretion disk GRBs, estimating the structure of each disk. We then use these structures in §§3 and 4 to predict the X-ray afterglow signature of these GRBs. We conclude with a discussion of the implications these results have on current and future X-ray missions in §5.

We note that the predictions concerning potentially time-variable absorption features as derived by Ghisellini et al. (1999) and Böttcher et al. (1999a) are valid, independent of the detailed GRB mechanism. Obviously, material responsible for such X-ray absorption features is located along the line of sight to the burst, so that the observed GRB and afterglow radiation and the assumption of a relativistic blast wave may be used.

## 2 Disk formation

Nearly all of the formation scenarios of black hole accretion disk GRBs invoke a “common-envelope” phase. Binaries are said to be in a common-envelope phase when the hydrogen envelope of one star engulfs its companion (this usually occurs when the star expands into a giant or supergiant). Friction and/or tidal forces cause the companion to spiral in towards the giant’s helium core, ejecting the hydrogen envelope. The evolution of two stars that enter a common envelope phase is one of the major uncertainties in binary population synthesis and, despite years of effort, still remains an open question (see Sandquist et al. 1998 and references therein). This limits our ability to make any strong quantitative predictions of this phase, but we can use the latest simulations to guide our estimates. Current simulations (e.g. Sandquist et al. 1998) suggest that the companion inspiral occurs on a timescale of 1 – 10 orbital periods ($`t_{\mathrm{orbit}}`$). The outcome of this inspiral is either: a) a close binary system if the companion is able to eject the hydrogen envelope before it merges with the star’s helium core or b) a merged object if the stars merge before the ejection of the hydrogen envelope.

The hydrogen envelope carries away much of the orbital angular momentum and is preferentially ejected in the orbital plane (Sandquist et al. (1998)). It is this material which forms the disks around GRB engines. Double NS, BH-NS, and BH-WD mergers go through a common envelope phase long before they actually merge, and the progenitors of these GRBs travel far from their formation sites (and their common-envelope disks) before producing GRBs (e.g. Fryer et al. (1999), Bulik & Belczyński (1999) and references therein). Collapsar and he-merger GRBs, on the other hand, occur shortly after their common-envelope phase. In this paper, we study the interaction of the ejecta from the GRB explosion with the common-envelope disks produced by collapsar and he-merger progenitors and determine the spectral features of this interaction. But first we must estimate the structure and size of these disks.

### 2.1 He-merger disks

To understand the characteristics of the disks produced in he-mergers, we must first understand the formation process of he-mergers (for more details, see Fryer et al. (1999)). The progenitor of a he-merger GRB is a binary system with two massive stars (both stars have masses in excess of $`8`$ M). The more-massive star (primary) evolves through its life, collapsing to form a compact remnant (either a neutron star or black hole). During the initial expansion, the primary may transfer mass to its companion (secondary star), and this mass transfer (or even a common-envelope phase) may tighten the orbital separation of these binaries. For some binaries, asymmetries in the primary’s supernova explosion also lead to a tighter binary system as the compact remnant is “kicked” into a closer orbit with its companion. In addition, the ejecta from the supernova explosion may enrich the envelope of the secondary (see Israelian et al. (1999)) with r-process elements.

When the secondary evolves off the main sequence, it envelops the primary’s compact remnant, and the binary goes into a common-envelope phase. The compact remnant ejects the hydrogen envelope in a disk-like structure, but not before merging with the secondary’s helium core. A he-merger occurs after the compact remnant has settled into the helium core. The inspiral process spins up the helium core, producing a disk around the compact remnant. In addition, the compact remnant accretes $`1`$ – 3 Mduring the inspiral and disk formation process (Zhang & Fryer (2000)), causing it to collapse to a black hole if it is not one already. This black-hole accretion disk system is surrounded by a disk formed from the enriched hydrogen envelope of the secondary.

There are no simulations of the common-envelope evolution of a compact remnant into a massive star which reliably predict the ejecta from the inspiral process. However, we may extrapolate from simulations such as Sandquist et al. (1998) to obtain a rough estimate of the characteristics of the hydrogen disk. The ejection velocity ($`v_{\mathrm{ejection}}`$) is roughly the escape velocity at any given radius in the star:

$$v_{\mathrm{ejection}}=\sqrt{\frac{2GM(r)_{\mathrm{secondary}}}{r}}.$$
(1)
The time between GRB outburst and the ejection of matter is equal to the inspiral time (since the outburst occurs as soon as the neutron star spirals into the center of the secondary):

$$t_{\mathrm{GRB}}10t_{\mathrm{orbit}}=\frac{20\pi r^{1.5}}{G^{0.5}M(r)_{\mathrm{secondary}}^{0.5}}.$$
(2)
At the time of the GRB outburst, the location of any layer of star with radius $`r`$ is simply:

$$D_{\mathrm{ejecta}}v_{\mathrm{ejection}}\times t_{\mathrm{GRB}}=20\sqrt{2}\pi r.$$
(3)
Clearly, the density structure of any he-merger disk depends upon the density structure of the secondary at the start of the common envelope phase. This structure is a function of both the size of the star (and hence, the orbital separation of the binary) at the beginning of the common-envelope phase and the mass distribution of the companion star. Using the binary population synthesis code developed in Fryer et al. (1999), we can calculate the distribution of orbital separations and companion masses for he-merger progenitors (Figs. 1 and 2). Fig. 3 shows the range of hydrogen disk structures formed by he-mergers with a 15 Mcompanion star for a series of orbital separations. Note that the outer disk radius is most sensitive to the orbital separation prior to the common envelope evolution. Fig. 4 shows the range of disk structures formed by 15, 25, and 40 Mstars assuming the orbital separation is set to the maximum radius of the companion. Although our rough estimates of the disk formation and the uncertainties in binary population synthesis and stellar evolution (see Fryer et al. (1999) for a discussion) make it difficult to predict the density profile of these disks accurately, these figures give a flavor of the range of possible disk structures.

### 2.2 Collapsar disks

The ring around Supernova 1987A proves that at least some massive stars have disks. The progenitor of Supernova 1987A was probably a binary system in which the more massive primary engulfed its companion, causing the companion to spiral into the primary during a common-envelope phase (Podsiadlowski (1992)). The companion was unable to eject the entire hydrogen envelope and it merged with the primary’s helium core. This process could lead to the formation of an outward moving disk which, in the case of supernova 1987A, was lit up 10,000 – 100,000 years later by the supernova to reveal a “ring” (see for example, Collins et al. (1999)).

The majority of Collapsar progenitors also follow an evolutionary path where a binary system goes through a common envelope phase although the companion often does not merge with the primary (Fryer et al. (1999)). The “classical” Collapsar model requires a massive helium star (without a hydrogen envelope) to avoid baryon contamination (MacFadyen & Woosley (1999)), and a common-envelope phase is required to eject most of the hydrogen envelope (Fryer et al. (1999)). Unlike the he-merger disks, the common envelope phase can occur more than 100,000 years before the eventual GRB outburst. Although the Collapsar progenitor is still likely to be surrounded by this ejecta, disks formed in Collapsar progenitors will be much further away at outburst than those disks formed in he-mergers. Assuming an ejection velocity equal to the escape velocity (eq. 1) and setting the time from disk ejection to GRB outburst ($`t_{\mathrm{GRB}}`$) to 100,000 y, we find that the inner edge of most Collapsar disks exceeds $`10^{17}`$ cm.

However, if the binary does not go into a common-envelope until just before the collapse of the primary, a much more compact disk may be formed. Recall that a common-envelope phase occurs when the radius of a star expands enough to engulf its companion. Most stellar models reveal that a star actually contracts during the last 10,000 – 100,000 y of its life (Schaerer et al. (1993), Woosley & Weaver (1995)). If the binary does not enter a common-envelope phase before this contraction, it is unlikely that it ever will. So with the current stellar models, $`t_{\mathrm{GRB}}`$ must be greater than 10,000 – 100,000 y. However, stellar models do not accurately predict the radii of massive stars, and there is a growing concern that these radii may be drastically incorrect (Fryer et al. (1999), Wellstein & Langer (1999), Fryer & Kalogera (2000)). It is possible that massive stars reach their maximum size just before collapse. If the common-envelope phase occurs 10 – 100 years before collapse, the inner edge of the disk could be less than $`10^{15}`$ cm. Until accurate stellar models are produced, we can not refine our estimates further. However, even with these rough estimates of the disk structure, we can now estimate the expected X-ray afterglow spectral features from both of these objects.

X-rays may also be produced when the GRB ejecta strikes the Collapsar’s companion star. Recall that the companion generally does not merge with the collapsar during the common envelope phase. The ejecta will hit this companion star, causing it to heat and expand, producing an X-ray emitting nebula (P. Pinto - private communication). The magnitude and spectra of X-rays under this mechanism is difficult to predict quantitatively, and we will delay further discussion of this emission for a later paper.

## 3 Predicted X-ray spectral features

### 3.1 Analytic estimate of the maximum iron line luminosity

In this section, we present a simple analytic estimate for an upper limit to the total, isotropic luminosity in the two major constituents of the iron K$`\alpha `$ line blend from a hot, highly ionized plasma, namely the Fe XXV He$`\alpha `$ (2p1s $``$ 1s2 resonance) and the Fe XXVI H$`\alpha `$ (2p $``$ 1s) transitions. Apart from energy conservation constraints, the luminosity in the resonant emission lines considered here is restricted by several line destruction mechanisms, such as Compton scattering, photoelectric absorption by lighter elements, and collisional de-excitation. Considering Compton scattering and photoelectric absorption — both of which processes would remove a line photon from the line — one can define an effectively emitting volume given by the fraction of the total disk volume through which the optical depth $`\tau _\mathrm{L}\tau _\mathrm{T}+\tau _{\mathrm{pa}}\tau _\mathrm{T}(1+f_{\mathrm{pa}})`$ due to electron scattering and photoelectric absorption is $`1`$. In a neutral plasma of cosmic element abundance, we have $`f_{\mathrm{pa}}4`$ at the energies of the $`n=21`$ transitions of Fe XXV and Fe XXVI. In an ionized plasma, this number obviously becomes smaller, and for our simple estimate we assume $`f_{\mathrm{pa}}2`$. The effect of collisional de-excitation becomes relevant if the density of the line-emitting material becomes comparable to or greater than the critical density $`n_{\mathrm{crit}}`$ at which the collisional de-excitation rate equals the spontaneous decay rate. Using van Regemorter’s (van Regemorter (1962)) $`\overline{g}`$ approximation, the critical density can be roughly approximated as

$$n_{\mathrm{crit}}^{\mathrm{c}.\mathrm{d}.}2.4\times 10^{11}\frac{\sqrt{T_e/\mathrm{K}}}{\lambda _{\mathrm{nm}}^3}\mathrm{cm}^33.8\times 10^{17}T_8^{1/2}\mathrm{cm}^3,$$
(4)
where $`T_8=T/(10^8\mathrm{K})`$.

For the ease of computation, we assume that the emitting region can be geometrically represented by a torus, located at a distance $`r=10^{13}r_{13}`$ cm from the center of the burst source, with a cross-sectional radius of $`a=10^{12}a_{12}`$ cm, containing a total mass of $`M_\mathrm{T}=m_\mathrm{T}M_{}`$, the volume of the torus will be given by $`V_\mathrm{T}=2\times 10^{38}a_{12}^2r_{13}`$ cm3. The average density of the torus material is then $`n_\mathrm{T}=6\times 10^{18}[m_\mathrm{T}/(a_{12}^2r_{13})]\mathrm{cm}^3`$. If the observer is looking along a line of sight close to the symmetry axis of the torus, then the Thomson depth through the torus is reasonably well approximated by

$$\tau _\mathrm{T}=an_\mathrm{T}\sigma _\mathrm{T}=4\times 10^6\left(\frac{m_\mathrm{T}}{a_{12}r_{13}}\right).$$
(5)
We define a critical density for Thomson scattering, $`n_{\mathrm{crit}}^\mathrm{T}`$ as the density at which the Thomson depth equals 1, so that for densities $`n>n_{\mathrm{crit}}^\mathrm{T}`$ line photons are likely to be scattered out of the line before leaving the emitting volume. This critical density is given by $`n_{\mathrm{crit}}^\mathrm{T}1.510^{12}a_{12}^1`$ cm-3. In the situations of interest here, we always find $`n_{\mathrm{crit}}^\mathrm{T}n_{\mathrm{crit}}^{\mathrm{c}.\mathrm{d}.}`$, indicating that line destruction by Compton scattering and photoelectric absorption is dominant over collisional de-excitation.

The luminosity in the emission lines may thus be estimated as $`L_{\mathrm{line}}=j_\mathrm{L}V_\mathrm{T}/\mathrm{max}\{\tau _\mathrm{L},\mathrm{\hspace{0.17em}1}\}`$, where $`j_L`$ is the emissivity in the line. In order to parametrize the line emissivities by an emissivity parameter $`x`$, we use the notation of Raymond & Smith (1977), where $`j_\mathrm{L}=n_\mathrm{H}n_e10^{23x}`$ erg cm-3 s-1. The emissivity parameters are inferred from runs of the XSTAR code (Kallman & McCray (1982)) with a negligibly small ionization parameter of $`\xi =10^8`$ ergs cm s-1 and constant, pre-specified plasma temperature. This yields

$$L_{\mathrm{line}}610^{45x}\frac{m_\mathrm{T}}{a_{12}}\mathrm{ergs}\mathrm{s}^1.$$
(6)
Some representative values of the emissivity parameters $`x`$ and the upper limits on the line luminosities are listed in Table 1, and the estimated luminosity upper limits are plotted as a function of plasma temperature in Fig. 5. Although these are very crude estimates, they indicate that Fe K$`\alpha `$ luminosities in excess of $`10^{44}`$ ergs s-1 are very well possible in a shock-heated-disk scenario.

### 3.2 Simulations of the shockwave / disk interaction

In order to get a realistic prediction of the expected X-ray line and continuum emission from the shock-heated disk, we simulate the shock-wave evolution and thermal history of the shocked material as the supernova ejecta, associated with the GRB explosion, interact with the pre-ejected disk. We assume that in the course of the supernova/GRB a total mass of $`M_01M_{}`$ is ejected quasi-isotropically at a speed of $`v_s=10^9v_9`$ cm s-1, initiating a shock-wave when interacting with the disk of pre-ejected material.

For the case of an intermediate-mass (secondary) progenitor (15, 25 $`M_{}`$), we parametrize the density profile of the disk by $`\rho _\mathrm{d}(r)=\rho _0(r/r_{\mathrm{in}})^{2.5}`$ for $`r_{\mathrm{in}}rr_{\mathrm{out}}`$ and its geometry by a constant $`h/r0.1`$. For more massive stars ($`40M_{}`$) the density profile may be approximated by a broken power-law with $`\rho _\mathrm{d}(r)=\rho _0(r/r_{\mathrm{in}})^4`$ for $`r_{\mathrm{in}}rr_{\mathrm{br}}`$, and $`\rho _\mathrm{d}(r)=\rho _0(r_{\mathrm{br}}/r_{\mathrm{in}})^4(r/r_{\mathrm{br}})^1`$ for $`r_{\mathrm{br}}rr_{\mathrm{out}}`$ (see Fig. 4). Writing $`r_{\mathrm{in}}=10^xr_{i,x}`$ cm and $`n_0=\rho _0/(\overline{A}m_p)=10^{17}n_{17}`$ cm-3, where $`\overline{A}`$ is the average atomic weight of the disk material, the total mass in the pre-ejected disk is

$$M_\mathrm{d}=0.1\left(\frac{h/r}{0.1}\right)n_{17}r_{i,13}^3\left(\sqrt{\frac{r_{\mathrm{out}}}{r_{\mathrm{in}}}}1\right)M_{}$$
(7)
for intermediate-mass progenitors, and

$$M_\mathrm{d}=0.5\left(\frac{h/r}{0.1}\right)n_{21}r_{i,12}^3\left(1+\frac{r_{\mathrm{in}}r_{\mathrm{out}}^2}{2r_{\mathrm{br}}^3}\right)M_{}$$
(8)
for high-mass progenitors.

The deceleration of the non-relativistic shock wave and the heating of the ejecta and and the swept-up material, are calculated by numerically solving simultaneously for the energy and momentum equations,

$$\frac{d}{dt}\left(Mv_\mathrm{s}\right)=4\pi R_\mathrm{s}^2\overline{P}_\mathrm{s},$$
(9)
$$\frac{dM}{dt}=4\pi R_\mathrm{s}^2\rho _\mathrm{d}(R_s)v_\mathrm{s},$$
(10)
$$\frac{dE}{dt}=\overline{P}_\mathrm{s}\frac{dV}{dt}+\dot{E}_{\mathrm{rad}},$$
(11)
where

$$E=\frac{1}{2}\overline{\rho }_\mathrm{s}v_\mathrm{s}^2V+\frac{\overline{P}_\mathrm{s}V}{\gamma 1},$$
(12)
$`R_\mathrm{S}`$ is the shock radius, $`\overline{P}_\mathrm{s}`$ and $`\overline{\rho }_\mathrm{s}`$ are the volume-averaged pressure and density of the shocked material, $`V`$ is the volume occupied by the shocked material, $`\gamma `$ is its adiabatic index, and $`\dot{E}_{\mathrm{rad}}`$ is the radiative cooling term. Throughout this paper, we assume $`\gamma =5/3`$. In Eq. 9, we have neglected the pressure of the disk material.

Due to the high densities involved and to the fact that the shock wave is non-relativistic, we may assume that at any point in time the shocked material is in approximate thermal and collisional ionization equilibrium. Photoionization precursors are not expected to play an important role since the Compton scattering depth $`\lambda _C=(n\sigma _\mathrm{T})^11.5\times 10^9n_{15}^1`$ cm is much smaller than any characteristic size scale of the system, so that most of the fluorescence photons produced in the photoionization precursor will be absorbed within the disk and thus be unobservable. The emission from the shocked region can therefore be represented by pure thermal plasma emission from an optically thick plasma at the temperature of the shocked material. At any given time, we calculate the emission from the shocked material using XSTAR (Kallman & McCray 1982) in a constant-temperature, purely thermal ionization mode (i. e. with very small photoionization parameter).

## 4 Results

We have done a series of simulations for a variety of parameters representative of both he-merger and collapsar/hypernova disks. In Fig. 6 we show the temporal evolution of the temperature of material behind the shock, resulting sample X-ray spectra at different times after the onset of the shock wave / disk interaction, and the light curve of the Fe K$`\alpha `$ line luminosity, for a disk with inner radius at $`r_{\mathrm{in}}=10^{13}`$ cm, representative for the he-merger case. In those simulations, we have assumed a disk density profile appropriate for a $`25M_{}`$ progenitor, and a mass of $`1M_{}`$ for the ejected material. The ejecta are assumed to have an initial velocity of $`10^9`$ cm s-1, corresponding to a kinetic energy of $`10^{51}`$ ergs in the ejecta. The figure shows that, especially in the later, decaying phase, $`2(1+z)\times 10^4`$ s after the onset of the shock wave / disk interaction, the thermal X-ray spectrum from the shocked disk material is strongly line dominated and might yield excellent prospects of detection by X-ray telescopes sensitive at $`1`$ keV.

Tab. 2 illustrates how the maximum Fe K$`\alpha `$ line luminosity and the decay time constant of the line emission depend on the mass and velocity of the ejecta shock-heating the pre-ejected disk from the $`25M_{}`$ progenitor. The time constant $`t_\mathrm{d}`$ is determined by fitting an $`\mathrm{exp}[(t/t_\mathrm{d})^2]`$ law to the decaying portion of the iron line light curves. We find that the shock wave / disk interaction expected for the he-merger scenario can very plausibly produce an Fe K$`\alpha `$ line of apparent quasi-isotropic luminosity $`L_{\mathrm{Fe}\mathrm{K}\alpha }10^{44}`$ ergs s-1 maintained over $`t_\mathrm{d}(1+z)\times 10^4`$ s after the ejecta begin to interact with the disk.

Fig. 7 shows the temperature evolution, X-ray spectra, and Fe K$`\alpha `$ light curve for a case representative of a collapsar / hypernova disk, if the system entered the common-envelope phase $`100`$ y before the GRB. The progenitor is assumed to be a $`25M_{}`$ star. In this case, the disk inner edge is expected to be located at $`r_{\mathrm{in}}10^{15}`$ cm. Results of simulations with different ejecta mass and velocity are summarized in Tab. 3. In this case, maximum iron line luminosities in excess of $`10^{42}`$ ergs s-1 are still possible, while the typical time delay between the GRB and the onset of the GRB is now $`(1+z)\times 10^6`$ s.

Comparing simulations with identical disk mass and ejecta mass and velocity, but different inner disk radii, we find an approximate scaling law for the maximum Fe K$`\alpha `$ line luminosity, $`L_{\mathrm{Fe}\mathrm{K}\alpha }r_{\mathrm{in}}^1`$.

## 5 Observational prospects

In the previous section, we found that in both the helium-merger and in the collapsar/hypernova scenarios, a quasi-thermal X-ray flash from the shock-heated disk may result. A critical and yet very uncertain parameter (especially for the collapsar/hypernova scenario) in the model simulations is the inner disk radius, which is primarily determined by the duration of the common-envelope phase prior to the GRB event. In the he-merger scenario, we expect typically $`r_{\mathrm{in}}10^{13}`$ cm, while in the collapsar/hypernova scenario, this parameter could have values $`10^{14}`$ cm $`r_{\mathrm{in}}10^{17}`$ cm. The onset and decay time scale of the resulting secondary X-ray flash scale as $`\mathrm{\Delta }t_Xr_{\mathrm{in}}`$, while the peak X-ray luminosity is approximately $`L_{\mathrm{Fe}\mathrm{K}\alpha }r_{\mathrm{in}}^1`$. Since the continuum X-ray afterglows from the (probably beamed) relativistic ejecta typically decay with temporal indices $`\chi 1.2`$, (if $`F_X(t)t^\chi `$), detection of the secondary, thermal X-ray outbursts predicted by the he-merger and collapsar/hypernova scenarios might be favored by larger disk radii, as long as the resulting X-ray flux remains above the detection threshold of currently operating X-ray telescopes.

Fig. 8 shows the absorbed 0.1 – 10 keV peak fluxes resulting from two of our simulations, as a function of redshift of the GRB source. Shown are representative cases for $`r_{\mathrm{in}}=10^{13}`$ cm and $`r_{\mathrm{in}}=10^{15}`$ cm (see also Figs. 6 and 7), corresponding to onset delays of $`(1+z)\times 10^4`$ s and $`(1+z)\times 10^6`$ s, respectively. Two plausible values of the Galactic neutral hydrogen column density, $`N_H`$, are used. Given that the nominal point source sensitivity for a 10 ksec exposure on Chandra’s ACIS detectors is $`4\times 10^{15}`$ ergs cm-2 s-1, while for the EPIC detectors on board XMM-Newton this limit is $`10^{14}`$ ergs cm-2 s-1, the predicted X-ray flashes from the shock wave / disk interaction may be detectable out to redshifts of at least $`z1`$ for an inner disk radius of $`r_{\mathrm{in}}10^{15}`$ cm. For a source at $`z=1`$ with $`r_{\mathrm{in}}=10^{15}`$ cm, we would expect the onset of the secondary X-ray burst $`3`$ weeks after the GRB.

## 6 Summary

During the progenitor evolution of collapsar/hypernova and he-merger GRBs, a hydrogen disk is formed around the central engine. These GRB engines produce a jet along the disk axis, and the relativistic outburst which produces the gamma-rays does not interact with this disk. However, both of these engines are likely to expel $`1`$ Malong the equator at lower velocities ($`10^9`$ cm s-1). The interaction of the explosion ejecta with the expelled disk may produce X-ray luminosities in excess of $`10^{44}(r_{\mathrm{in}}/10^{13}\mathrm{cm})^1`$ ergs s-1 with a delay of $`(r_{\mathrm{in}}/10^{13}\mathrm{cm})(1+z)\times 10^4`$ s after the GRB. For gamma-ray bursts with a redshift $`z1`$, this emission is well within the capabilities of the latest X-ray satellites (e.g. Chandra and XMM-Newton). Thus, long-term monitoring of X-ray afterglows over several weeks after the GRB (most notably, even after the direct afterglow radiation has faded to undetectable levels) may lead to the detection of these secondary X-ray flashes which would yield valuable information about the nature and pre-burst evolution of the GRB progenitor.

The work of M.B. is supported by NASA through Chandra Postdoctoral Fellowship Award Number PF 9-10007, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS 8-39073. C.L.F was supported by a Feynman Fellowship at LANL, NSF (AST-97-31569), and the US DOE ASCI Program (W-7405-ENG-48). It is a pleasure to thank P. Pinto for helpful advice and discussions.


text:
# Self-consistent description of simplified pseudospin-electron model

## Abstract

A method of the self-consistent calculation of the thermodynamical and correlation functions is presented. This approach is based on the GRPA (generalized random phase approximation) scheme with the inclusion of the mean field corrections. The numerical research shows that interaction between the electron and pseudospin subsystems leads to the possibility of the dipole (pseudospin) and charge-density instabilities as well as phase separation into the uniform and/or chess-board phases.

The aim of this article is to calculate the correlation functions, pseudospin and particle number mean values as well as the grand canonical potential for the simplified ($`U=0`$) pseudospin-electron model (PEM)  within the one self-consistent approach.

PEM includes the local interaction of the conducting electrons with some two level subsystem described by pseudospins $`S^z=\pm 1/2`$ (e.g anharmonic vibrations of the apex oxygen ions in YBaCuO-type crystals) . On the other hand, PEM can be transform into the binary alloy type model as well as into the Falikov-Kimball model (FKM) and they can be also studied analytically within the framework of the developed scheme.

The model Hamiltonian is the following:

$`H={\displaystyle \underset{i}{}}H_{i\mathrm{\hspace{0.17em}0}}+{\displaystyle \underset{ij\sigma }{}}t_{ij}c_{i\sigma }^+c_{j\sigma },`$ (1)
$`H_{i\mathrm{\hspace{0.17em}0}}=\mu {\displaystyle \underset{\sigma }{}}n_{i\sigma }+g{\displaystyle \underset{\sigma }{}}n_{i\sigma }S_i^zhS_i^z,`$
where an interaction with pseudospins ($`g`$-term) placed in some longitudinal field $`h`$ (chemical potential for ions in FKM) are included in the single-site part; $`\mu `$ is the chemical potential. The second term in the Hamiltonian describes an electron hopping from site to site.

The calculations are performed in the strong coupling case ($`gt`$) using single-site states as the basic one. A formalism of the electron annihilation (creation) operators $`a_{i\sigma }=c_{i\sigma }P_i^+,`$ $`\stackrel{~}{a}_{i\sigma }=c_{i\sigma }P_i^{}`$ ($`P_i^\pm =\frac{1}{2}\pm S_i^z`$) acting at a site with the certain pseudospin orientation is introduced.

Expansion of the calculated quantities in terms of the electron transfer leads to an infinite series of terms containing the averages of the $`T`$-products of the $`a_{i\sigma }`$, $`\stackrel{~}{a}_{i\sigma }`$ operators. The evaluation of such averages is made using the corresponding Wick’s theorem. The results are expressed in terms of the products of the nonperturbed Green’s functions and averages of the projection operators $`P_i^\pm `$ which are calculated by means of the semi-invariant expansion .

The calculation of the correlation functions is performed within a self-consistent scheme in the framework of the generalized random phase approximation  with the inclusion of the mean field type contributions coming from the effective pseudospin interaction via conducting electrons, i.e. all zero-order correlators (second order semi-invariants) as well as pseudospin mean value are calculated by the mean-field Hamiltonian:

$`H_{\mathrm{MF}}={\displaystyle \underset{i}{}}H_{i\mathrm{\hspace{0.17em}0}}+\alpha _1P_i^++\alpha _2P_i^{},`$ (2)
$`\alpha _1P_i^++\alpha _2P_i^{}={\displaystyle \frac{2}{N\beta }}{\displaystyle \underset{n,𝒌}{}}{\displaystyle \frac{t_𝒌^2}{g_n^1t_𝒌}}\left[{\displaystyle \frac{P_i^+}{\mathrm{i}\omega _n\epsilon _1}}+{\displaystyle \frac{P_i^{}}{\mathrm{i}\omega _n\epsilon _2}}\right]`$
$`g_n={\displaystyle \frac{P^+}{\mathrm{i}\omega _n\epsilon _1}}+{\displaystyle \frac{P^{}}{\mathrm{i}\omega _n\epsilon _2}},\epsilon _{1,2}=\mu \pm {\displaystyle \frac{g}{2}}.`$ (3)
Diagram equation on pseudospin correlator $`S^zS^z_𝒒`$ is following:

This equation differs from the one for the Ising model in MFA by the replacement of the exchange interaction by the electron loop (which describes an interaction between pseudospins via conducting electrons). In the analytical form its solution is equal

$`\chi ^{S^zS^z}(\omega _m,𝒒)={\displaystyle \frac{\delta (\omega _m)P^+P^{}}{T\mathrm{\Theta }(T,𝒒)}},`$ (4)
$`\mathrm{\Theta }(T,𝒒)={\displaystyle \frac{2}{\beta }}{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝒌}{}}\mathrm{\Lambda }_n^2\stackrel{~}{t}_n(𝒌)\stackrel{~}{t}_n(𝒌+𝒒),`$
$`\mathrm{\Lambda }_n={\displaystyle \frac{g\sqrt{P^+P^{}}}{(\mathrm{i}\omega +\mu )^2g^2/4}},\stackrel{~}{t}_n(𝒌)={\displaystyle \frac{t_𝒌}{(1g_nt_𝒌)}}.`$
Equation for pseudospin mean value in the uniform case ($`S_i^z=S^z`$) is following:

$$S^z=\frac{1}{2}\mathrm{tanh}\left\{\frac{\beta }{2}(h+\alpha _2\alpha _1)+\mathrm{ln}\frac{1+\mathrm{e}^{\beta \epsilon _1}}{1+\mathrm{e}^{\beta \epsilon _2}}\right\}$$

We also calculate particle number mean value and grand canonical potential . All quantities can be derived from the grand canonical potential by

$$\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}(h)}=S^z,\frac{\mathrm{d}S^z}{\mathrm{d}(\beta h)}=S^zS^z_{𝒒=0},$$

which show the thermodynamical consistence of the proposed approximations.

The analysis of the pseudospin correlator temperature behaviour shows that high temperature phase become unstable with respect to fluctuations with wave-vector $`𝒒=(\pi ,\pi )`$ (for some model parameters values) that indicates the possibility of the phase transition into a modulated (chess-board) phase. On the other hand, the phase transition between the uniform phases with different pseudospin mean values (bistability) is possible .

From the comparison of the grand canonical potential $`\mathrm{\Omega }`$ values for uniform and chess-board phases ($`\mu =\mathrm{const}`$ regime), the ($`\mu ,h`$) phase diagram is obtained (Fig. 1a).

Chess-board phase exists as intermediate one between the uniform phases with different $`S^z`$ and $`n`$ values. The transition between different uniform phases (bistability) is of the first order (Fig. 1a, dashed line), while the transition between the uniform and modulated ones is of the first (dotted line) or second (solid line) order.

On the other hand, the minimum of the free energy $`F=\mathrm{\Omega }+\mu N`$ is the equilibrium condition in the $`n=\mathrm{const}`$ regime. In this regime the first order phase transition with a jump of the pseudospin mean value accompanied by the change of electron concentration transforms into a phase separation into the regions with different phases (the uniform and chess-board ones) and with different electron concentrations and pseudospin mean values (Fig. 1b).

On the basis of the presented scheme the thermodynamics of phase transitions, possibility of the phase separations as well as appearance of the chess-board phase have been investigated.


text:
# Measurement of Newton’s Constant Using a Torsion Balance with Angular Acceleration Feedback

## Abstract

We measured Newton’s gravitational constant $`G`$ using a new torsion balance method. Our technique greatly reduces several sources of uncertainty compared to previous measurements: (1) it is insensitive to anelastic torsion fiber properties; (2) a flat plate pendulum minimizes the sensitivity due to the pendulum density distribution; (3) continuous attractor rotation reduces background noise. We obtain $`G=(6.674215\pm 0.000092)\times 10^{11}`$ m3 kg-1 s-2; the Earth’s mass is, therefore, $`M_{}=(5.972245\pm 0.000082)\times 10^{24}`$ kg and the Sun’s mass is $`M_{}=(1.988435\pm 0.000027)\times 10^{30}`$ kg.

The gravitational constant $`G`$ along with Planck’s constant $`\mathrm{}`$ and the speed of light $`c`$, is one of the most fundamental and universal constants in nature. Unlike most other physical constants, the value of $`G`$ is not precisely known due to the weakness and non-shieldability of gravity. Since the first laboratory measurement by Cavendish over 200 years ago, the reduction in uncertainty in $`G`$ has been only a factor of about ten per century. To make matters worse, measurements over the last ten years  all have larger error-bars than the value of Luther and Towler  on which the 1986 accepted value was based. In addition, some of these recent measurements lie far outside the $`\pm `$128 ppm uncertainty of the 1986 accepted value, bringing this accuracy of $`G`$ into question. Recognizing this situation, the currently recommended value is $`G=(6.673\pm 0.010)\times 10^{11}`$ m3 kg-1 s-2, corresponding to an uncertainty of $`\pm `$1500 ppm.

We developed a new technique  that offers several substantial advantages over previous methods. At the heart of the apparatus is a torsion balance placed on a turntable located between a set of attractor spheres. The turntable is first rotated at a constant rate so that the pendulum experiences a sinusoidal torque due to the gravitational interaction with the attractor masses. A feedback is then turned on that changes the rotation rate so as to minimize the torsion fiber twist. The resulting angular acceleration of the turntable, which is now equal to the gravitational angular acceleration of the pendulum , is determined from the second time-derivative of the turntable angle readout. Since the torsion fiber does not experience any appreciable deflection, this technique is independent of many torsion fiber properties including anelasticity, which may have led to a bias in previous measurements.

The gravitational angular acceleration $`\alpha `$ of a torsion pendulum expanded in spherical multipole moments is:

$$\alpha (\varphi )=\frac{4\pi G}{I}\underset{l=2}{\overset{\mathrm{}}{}}\frac{1}{2l+1}\underset{m=l}{\overset{+l}{}}mq_{lm}Q_{lm}e^{im\varphi },$$
(1)
where $`q_{lm}`$ are the multipole moments of the pendulum and the $`Q_{lm}`$multipole fields of the external mass distribution, $`\varphi `$ is the turntable angle with respect to the mass distribution, and $`I`$ is the pendulum moment of inertia about the torsion fiber. Since the series converges rapidly, the biggest term in Eq. 1 is the one with $`l=m=2`$:

$$\alpha _{22}(\varphi )=\frac{16\pi }{5}G\frac{q_{22}}{I}Q_{22}\mathrm{sin}2\varphi .$$
(2)
For an ideal, infinitely-thin two-dimensional vertical plate the quotient $`q_{22}/I`$ is independent of the pendulum mass distribution, $`\rho (\stackrel{}{r}_p)`$, i.e. mass, size, shape etc.:

$$\frac{q_{22}}{I}=\frac{\rho (\stackrel{}{r}_p)Y_{22}(\theta _p,\varphi _p)r_p^2d^3r_p}{\rho (\stackrel{}{r}_p)\mathrm{sin}^2\theta _pr_p^2d^3r_p}\stackrel{2D}{}\sqrt{\frac{15}{32\pi }}$$
(3)
and

$$\alpha (\varphi )\alpha _{22}(\varphi )=\sqrt{\frac{24\pi }{5}}GQ_{22}\mathrm{sin}2\varphi .$$
(4)
The flat pendulum geometry presents a significant advance over previous torsion balance measurements, where the biggest contribution to the systematic uncertainty was due to the pendulum’s mass distribution. A small and easily calculable modification can be applied to the ideal ratio for a thin rectangular plate with finite thickness $`t`$ and width $`w`$:

$$\frac{q_{22}}{I}=\frac{w^2t^2}{w^2+t^2}\sqrt{\frac{15}{32\pi }}.$$
(5)
All $`l>2,m=2`$ terms in Eq. 1 also contribute to the signal $`\mathrm{sin}(2\varphi )`$. However, these contributions are small since the expansion in Eq. 1 converges as $`(R_{pend}/R_{attr})^{l2}`$, where $`R_{pend}`$ and $`R_{attr}`$ are the characteristic radii of the pendulum and the source. In addition, all $`l=`$ odd pendulum moments and attractor fields are removed by designing the pendulum “up-down” symmetric. The $`l=4,m=2`$ contribution is also eliminated by design: the $`q_{42}`$ of a rectangular plate pendulum with height $`h`$ vanishes by choosing $`h^2=\frac{3}{10}(w^2+t^2)`$, and the $`Q_{42}`$ vanishes by using two spheres on each side with their centers spaced vertically by $`z=\sqrt{2/3}\rho `$, where $`\rho `$ is the radial distance from the pendulum axis. The lowest $`l>2,m=2`$ contribution which does not vanish by design is the small and easily calculable $`\alpha _{62}`$:

$$\frac{\alpha _{62}}{\alpha _{22}}=\frac{99}{7683200}\frac{213(w^4+t^4)+626w^2t^2}{\rho ^4}.$$
(6)
Table I contains the numeric values of the $`\alpha _{62}`$ and $`\alpha _{82}`$ corrections. We checked the accuracy of our multipole analysis with a full angular acceleration calculation using direct numeric integration.

The attractor spheres are located on a separate coaxial turntable which is rotated with angular velocity $`\omega _a(t)=\omega _d+\omega _i(t)`$, where $`\omega _i(t)`$ is the angular velocity of the torsion balance turntable. The difference angular velocity, $`\omega _d=\dot{\varphi }`$, is held constant. Rotation of the attractor masses allows us to cleanly remove gravitational interactions due to the environment. Furthermore, we are able to set the signal, $`\mathrm{sin}(2\omega _dt)`$, at a relatively high frequency to suppress the 1/$`f`$-noise characteristic of the torsion balance and the gravitational background.

A schematic of the apparatus is shown in Fig. 1. The torsion balance turntable consisted of an air-bearing, a precision angle encoder, and an eddy-current motor. The torsion pendulum was located in an aluminum vacuum chamber and was surrounded by a $`\mu `$-metal shield. The pendulum was hung from a 41.5 cm long, 17 $`\mu `$m diameter tungsten-fiber, which was attached to a swing damper. The pendulum was a 1.506 mm thick, 76 mm wide, and 41.6 mm tall Pyrex-glass plate with a thin gold-coating. The small pendulum deflection angle was sensed with an autocollimator using four reflections off the pendulum plate.

The centers of the spheres were located at $`\rho `$=16.76 cm on three stainless steel seats that were embedded in two cast-aluminum shelves. The shelves were supported by a turntable made from a steel bearing. The attractor spheres were machined from the same selected stock of ultrasonically tested #316 stainless steel. Their average diameter was 124.89 mm and their mass was $``$8.140 kg . A pressure-dependent air-density correction was applied to the $`G`$-measurements. The apparatus was located in the former cyclotron cave in the Nuclear Physics Lab at the University of Washington on a massive platform 3.5 m above the floor. The partially underground room was temperature stabilized. Temperature drifts and fluctuations were typically $`<0.05`$ K/day. The instrument itself was in a passive thermal enclosure made of Styrofoam. The apparatus-temperature was monitored and a correction was applied to compensate for the thermal expansion of the attractor mass assembly during the measurement.

A digital signal processor (DSP) recorded the data and controlled the experiment. The DSP’s loop frequency (2.5 KHz) was used for the timing, which was derived from a temperature-controlled quartz oscillator which was calibrated with a GPS receiver. The data were averaged by the DSP over exactly one second and uploaded to the host PC.

We recorded six data runs, each approximately three days long. A typical data segment is shown in Fig 2. After each data run the spheres were moved to locations differing by $`90^{}`$ in azimuthal angle on the attractor turntable shelves. In addition, the spheres from the upper shelf were placed on the lower shelf and the orientation of the spheres were changed in order to average over density fluctuations and non-sphericity. The position of the spheres were measured before and after each data run. We then repeated the entire measurement cycle using four different spheres.

Our largest systematic uncertainty was due to the attractor mass distance measurement. We used a specially fabricated micrometer tool made primarily of Invar for the horizontal distance measurement between the spheres. Before and after each distance measurement the tool was compared to an Invar ball bar standard that was calibrated at NIST to within 0.2 $`\mu `$m. The vertical spacing between the spheres was measured by inserting a small gauge block that was $``$10-20 $`\mu `$m thinner than the gap between the sphere surfaces. We inferred the spacing by optically measuring the angle through which the gauge block could be tilted. The temperature of the attractor mass was recorded during the distance measurements.

For most of the data $`\omega _i`$ was set to $``$5.3 mrad/s and $`\omega _d`$ was set to 20.01015 mrad/s so that the signal frequency occurred at $``$6.37 mHz. The rotation frequencies were chosen to be incommensurate. We tested a wide variety of $`\omega _i`$ and $`\omega _d`$ and found the results independent of these angular velocities.

The lab-fixed horizontal magnetic field at the apparatus was measured to be $``$100 mG. We ran tests by exaggerating the field at the center of the apparatus to $``$5 G and at the location of the spheres to $``$100 G. The observed acceleration difference due to the exaggerated fields was $`(6\pm 8)\times 10^{12}`$ rad/s2. Therefore, a 0.6 ppm error was attributed to magnetic accelerations.

To investigate the rotating temperature gradient sensitivity, two 5W-heaters were installed on the turntable. These rotating heaters exaggerated the normally observed temperature variation at $`2\omega _d`$ by greater than a factor of 250. The acceleration difference amplitude with the heaters activated was $`|\mathrm{\Delta }\alpha |=(21\pm 22)\times 10^{12}`$ rad/s2. Therefore the rotating temperature gradient coupled acceleration was $`<2\times 10^{13}`$ rad/s2 (0.4 ppm).

The torsion balance turntable angle was numerically differentiated twice to yield angular acceleration. All the data were divided into segments comprised of twenty sinusoidal cycles and fitted using a least squares method. The fitting function included the signal and its harmonics, the room background and its harmonics, an offset and a linear drift. The statistical error was derived from the scatter of the individual fit values.

Table 1 lists the corrections that were applied to the data. Since the DSP averaged over $`\tau =1`$ s intervals the true amplitudes are attenuated by a small amount: $`\frac{\mathrm{sin}(\omega _d\tau )}{\omega _d\tau }`$. Also the numerical derivative of $`\mathrm{sin}(2\omega _dt)`$ is smaller than the analytic derivative by $`\frac{\mathrm{sin}(\omega _d\mathrm{\Delta }t)}{\omega _d\mathrm{\Delta }t}`$, where $`\mathrm{\Delta }t`$ $`=10`$ s was the chosen time increment. Since this correction was substantial we tested a wide variety of increments and found that our results were independent of the choice of $`\mathrm{\Delta }t`$. We verified the accuracy of the data analysis by simulating the data numerically using various levels of drifts, 1/$`f`$-noise, and gravitational room background.

The overall linearity and insensitivity to the signals at $`4\omega _d`$ was tested using two pairs of spheres separated by $`45^{}`$ on each side. These 8-sphere results agreed with the regular 4-sphere results within the statistical errors.

Due to the finite gain of the acceleration feedback loop a small twist remained in the fiber. The twist angle at $`2\omega _d`$ was essentially unresolved so that a 0.35 ppm error was used instead of a correction.

We combined the results from the different attractor mass configurations into pairs that averaged out accelerations due to the attractor turntable itself and the combination of three pairs yielded our $`G`$-value (Fig. 3). After combining the two $`G`$-values obtained with different spheres our value for the gravitational constant is

$$G=(6.674215\pm 0.000092)\times 10^{11}\mathrm{m}^3\mathrm{kg}^1\mathrm{s}^2.$$
(7)
The systematic and statistical uncertainties are shown in Table 2. The total uncertainty is the quadrature sum of the individual uncertainties. Fig. 4 compares our results to other recent measurements.

Combining our value of $`G`$ with the geocentric value $`GM_{}`$  we determine the mass of the Earth to be

$$M_{}=(5.972245\pm 0.000082)\times 10^{24}\mathrm{kg}.$$
(8)
Likewise, using $`GM_{}`$ , the mass of the Sun is

$$M_{}=(1.988435\pm 0.000027)\times 10^{30}\mathrm{kg}.$$
(9)
We wish to thank the members of the Eöt-Wash group and Maximilian Schlosshauer for writing the data simulations. Many thanks go to NIST for the support of this work with precision measurement grant # 60NANB7D0053 and to NSF for grant # PHY-9602494 to the Eöt-Wash group. Further assistance was provided by the Physics Department of the University of Washington, and the Nuclear Physics Laboratory which is supported by the DOE.