text:
# Gravitational Collapse in Turbulent Molecular Clouds. II. Magnetohydrodynamical Turbulence

## 1 Introduction

All star formation takes place in molecular clouds. However, the star formation rate in these clouds is surprisingly low. From the Jeans argument one would expect that they should collapse within their free fall time

$$t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G\overline{\rho }}}(1.06\times 10^6\text{ yr})\left(\frac{n}{10^3\text{ cm}^3}\right)^{1/2},$$
(1)
where $`\overline{\rho }`$ is the mean mass density of the cloud, $`G`$ the gravitational constant and $`n=\overline{\rho }/\mu `$ the number density, with $`\mu =2.36m_H`$. A catastrophic collapse of a giant molecular cloud on this timescale would yield a single starburst event. However, molecular clouds have classically been thought to survive without global collapse for much longer than their free-fall time $`t_{\mathrm{ff}}`$ (Blitz & Shu 1980). Moreover, stars are not usually observed in nearby star-forming regions to form in such a catastrophic collapse. Instead, they form in localized regions dispersed through an apparently stable cloud.

Observations of spectral line widths in molecular clouds show that the gas moves at speeds exceeding the thermal velocities by up to an order of magnitude (Williams et al. 1999). These supersonic motions seem not to be ordered, so that turbulent support models have often been suggested, with the turbulence giving rise to an effectively isotropic turbulent pressure counteracting the gravitational forces.

However, such models have two problems. First, simulations of supersonic, compressible turbulence show that it typically decays in a time less than the cloud’s free-fall time $`t_{\mathrm{ff}}`$ (Gammie & Ostriker 1996, Mac Low et al. 1998, Mac Low 1999). So, in order to be really able to support the cloud, the turbulence would have to be constantly driven. Second, although hydrodynamical turbulence can prevent global collapse, it can never completely prevent local collapse except with unrealistically short driving length scale $`\lambda _D`$ (Klessen, Heitsch & Mac Low 2000; hereafter Paper I). The efficiency of local collapse depends on the wavelength and on the strength of the driving source. Long wavelength driving or no driving at all results in efficient, coherent star formation, with most collapsed regions forming near each other (Klessen & Burkert 2000). Strong, short wavelength driving, on the other hand, results in inefficient, incoherent star formation, with isolated collapsed regions randomly distributed throughout the cloud.

The model of molecular clouds being supported by turbulence has been widely discussed and investigated, as reviewed in Paper I and Vázquez-Semadeni et al. (2000). Recently, Ballesteros-Paredes et al. (1999), Hennebelle & Pérault (1999), and Elmegreen (2000) have suggested that molecular clouds might not have to be supported for these long time scales at all, but might be transient features caused by colliding flows in the interstellar medium. This would solve very naturally not only the problem of cloud support, but also, according to Ballesteros-Paredes et al. (1999), the pronounced lack of 5–10 million year old post-T Tauri-stars directly associated with star-forming molecular clouds.

Magnetic fields might alter the dynamical state of a molecular cloud sufficiently to prevent gravitationally unstable regions from collapsing (McKee 1999). They have been hypothesized to support molecular clouds either magnetostatically or dynamically through MHD waves.

Mouschovias & Spitzer (1976) derived an expression for the critical mass-to-flux ratio in the center of a cloud for magnetostatic support. Assuming ideal MHD, a self-gravitating cloud of mass $`M`$ permeated by a uniform flux $`\mathrm{\Phi }`$ is stable if the mass-to-flux ratio

$$\frac{M}{\mathrm{\Phi }}<\left(\frac{M}{\mathrm{\Phi }}\right)_{cr}\frac{c_\mathrm{\Phi }}{\sqrt{G}}.$$
(2)
with $`c_\mathrm{\Phi }`$ depending on the geometry and the field and density distribution of the cloud. A cloud is termed subcritical if it is magnetostatically stable and supercritical if it is not. Mouschovias & Spitzer (1976) determined that $`c_\mathrm{\Phi }=0.13`$ for their spherical cloud. Assuming a constant mass-to-flux ratio in a region results in $`c_\mathrm{\Phi }=1/(2\pi )0.16`$ (Nakano & Nakamura 1978). Without any other mechanism of support such as turbulence acting along the field lines, a magnetostatically supported cloud will collapse to a sheet which then will be supported against further collapse. Fiege & Pudritz (1999) discussed a sophisticated version of this magnetostatic support mechanism, in which poloidal and toroidal fields aligned in the right configuration could prevent a cloud filament from fragmenting and collapsing.

Investigation of the second alternative, support by MHD waves, concentrates mostly on the effect of Alfvén waves, as they (1) are not as subject to damping as magnetosonic waves and (2) can exert a force along the mean field, as shown by Dewar (1970) and Shu et al. (1987). This is because Alfvén waves are transverse waves, so they cause perturbations $`\delta 𝐁`$ perpendicular to the mean magnetic field $`𝐁`$. McKee & Zweibel (1994) argue that Alfvén waves can even lead to an isotropic pressure, assuming that the waves are neither damped nor driven. However, in order to support a region against self-gravity, the waves would have to propagate outwardly, rather than inwardly, which would only further compress the cloud. Thus, as Shu et al. (1987) comment, this mechanism requires a negative radial gradient in wave sources in the cloud.

Most molecular clouds show evidence of magnetic fields. However, the discussion of their relative strength is lively. Crutcher (1999) summarizes all 27 available Zeeman measurements of magnetic field strengths in molecular clouds. He concludes that (1) static magnetic fields are not strong enough to support the observed clouds alone, with typical ratios of the mass $`M`$ to the critical mass for the observed magnetic field $`M/M_{cr}2`$; (2) the ratio of thermal to magnetic pressure $`\beta =P_{th}/P_{mag}0.04`$; (3) internal motions are supersonic, with a velocity dispersion $`\sigma _vc_s`$ but approximately equal to the Alfvén speed, $`\sigma _vv_A`$; and (4) that the kinetic and magnetic energies are roughly equal, which he interprets as suggesting that static magnetic fields and MHD waves are equally important in cloud energetics. However, McKee (1999) remarks that (a) Crutcher’s data do not address the strength of the field on large scales (threading an entire GMC) and (b) that the data deal with dense regions in the clouds, so that ambipolar diffusion already might have altered the mass-to-flux ratio observed. Moreover, Crutcher’s fourth conclusion about the role of static fields and MHD waves is based on the assumption that the kinetic energy stems mainly from MHD waves.

Nakano (1998) made two further arguments against magnetostatic support of cloud cores. First, magnetically subcritical condensations cannot have column densities much higher than their surroundings. However, observed cloud cores have column densities significantly higher than the mean column density of the cloud, indicating that they are not magnetostatically supported. Second, if the cloud cores were magnetically supported and subcritical it would be difficult to maintain the observed non-thermal velocity dispersions for a significant fraction of their lifetime. Mac Low (1999) confirmed this by numerically determining the dissipation rate of supersonic, magnetohydrodynamic turbulence. He concludes that the typical decay time constant is far less than the free-fall time of the cloud.

Polarization measurements might give us a clue whether the fields are well ordered or in a turbulent state. However, up to now, most measurements refer to the highest density regions, thus giving information about the fields in small scale structures, but not about scales of the whole cloud. Hildebrand et al. (1999) present polarization measurements of cloud cores and envelopes. They find polarization degrees of at most $`10`$%. More recent observations of the molecular cloud filament OMC-3 by Matthews & Wilson (1999) suggest that the magnetic field is well ordered perpendicularly to the filament, but with a mean polarization degree of only $`4.2`$%.

Vázquez-Semadeni et al. (1996) performed three-dimensional simulations including self-gravity and MHD with a resolution of $`64^3`$. They found that hydrodynamical and supercritically magnetized turbulence can lead to gravitationally bound structures. Gammie & Ostriker (1996) did simulations in $`1\text{ }2/3`$ dimensions, while more recently $`2.5`$ dimensional models were presented by Ostriker et al. (1999). Mac Low et al. (1998), Stone et al. 1998, Padoan & Nordlund (1999) and Mac Low (1999) studied decaying magnetized turbulence and found short decay rates with as well as without magnetic fields.

We present the first high-resolution ($`256^3`$ zones) simulations of magnetized, self-gravitating, driven, supersonic turbulence to test the hypothesis that magnetic fields can contribute to the support of molecular clouds. The following section describes the technique and parameters used for the simulations. In section 3 we discuss requirements regarding the resolution needed for simulations of self-gravitating magnetized turbulence. In section 4 we present the results, and we summarize our conclusions in section 5.

## 2 Technique and Models

### 2.1 Technique

For our computations, we use ZEUS-3D, a well-tested, Eulerian, finite-difference code (Stone & Norman 1992a, b; Clarke 1994). It uses second-order advection and resolves shocks employing a von Neumann artificial viscosity. Self-gravity is implemented via an FFT-solver for cartesian coordinates (Burkert & Bodenheimer 1993). The magnetic forces are calculated via the constrained transport method (Evans & Hawley 1988) to ensure $`𝐁=0`$ to machine accuracy. In order to stably propagate shear Alfvén-waves, ZEUS uses the method of characteristics (Stone & Norman 1992b, Hawley & Stone 1995). This method evolves the propagation of Alfvén-waves as an intermediate step to compute time-advanced quantities for the evolution of the field components themselves to ensure that signals do not propagate upwind unphysically.

We use a three-dimensional, periodic, uniform, Cartesian grid for the models described here. This gives us equal resolution in all regions, and allows us to resolve shocks and magnetic field structures well everywhere. On the other hand, collapsing regions cannot be followed to scales less than a few grid zones.

We do not include ambipolar diffusion in this work, although numerical diffusion acts on scales in our simulations that, when translated to astronomical scales with typical parameters, correspond to the scales on which ambipolar diffusion begins to dissipate power (see Paper I). Moreover, we do not use a physical resistivity, relying on numerical dissipation at grid scale, due to averaging quantities over zones. The limited resolution in high-density regions can lead to excessive numerical diffusion of mass through the magnetic fields that must be accounted for when analyzing these or similar computations. In § 3 we derive the appropriate resolution criterion.

### 2.2 Models & Parameters

We employ two sets of parameters. The first one is the same as in Paper I and allows us to compare its results to the ones of this work. The second parameter set enables us to determine the numerical reliability of these results. All parameters are given in normalized units, where physical constants are scaled to unity and where we consider gas cubes with mass $`M1`$ and side length $`[1.0,1.0]`$. The system can be scaled to physical units using the Jeans mass $`M_J`$ and length scale $`\lambda _J`$. In Paper I we adopt a normalized sound speed of $`c_s=0.1`$, which yields $`M_J=0.0156`$ and $`\lambda _J=0.5`$ such that the computed box contains $`M=64M_J`$ and is $`L=4\lambda _J`$ in size. Our second parameter set is based on a sound speed of $`c_s=0.213`$, which in turn gives a ten times larger Jeans mass $`M_J=0.156`$ and yields $`\lambda _J=1.068`$, so that follows $`M=6.4M_J`$ and $`L=1.873\lambda _J`$.

We use the same driving mechanism described in Mac Low (1999). Each time step, a fixed velocity pattern is added to the actual velocity field thus, that the input energy rate is constant. The driving field pattern is derived from a Gaussian random field with a given spectrum. This allows us to drive the cloud on selected spatial scales. The rms Mach number for the first parameter set is $`_{rms}=10`$, for the second one it is $`_{rms}=5`$.

The MHD simulations start with a uniform magnetic field in the $`z`$-direction. As soon as the cube has evolved into a fully turbulent state, gravity is switched on. The critical field value according to equation 2 is $`B_{cr}=1.56`$ in code units. The initial field strength varies between $`B=0.191.77`$, corresponding to $`M/M_{cr}=0.48.3`$, covering the sub- and supercritical range. The ratio $`\beta =P_{th}/P_{mag}=0.014.04`$. For an overview of the parameters used see Table 1. There are four series of models: $`𝒟`$ (pure hydrodynamics, same runs as in Paper I); $``$ (MHD models with the same parameter set as the hydro models); $``$ (magnetostatic models for numerical tests); and $`𝒢`$ (MHD models with a reduced number of Jeans masses). The second letter in the model names stands for the resolution (low, $`64^3`$; intermediate, $`128^3`$; or high, $`256^3`$), followed by a digit denoting the driving scale. For the MHD models, a third letter gives the relative field strength (low, intermediate, moderate, strong).

The dynamical behavior of isothermal, self-gravitating gas depends only on the ratio between potential and kinetic energy. Therefore we can use the same scaling prescriptions as in Paper I, defining the physical time scale by the free-fall time (eq. 1), the length scale by the initial Jeans length

$$\lambda _J=c_\mathrm{s}(\pi /(G\overline{\rho }))^{1/2}(0.71\text{ pc})\widehat{\lambda }_J\left(\frac{c_\mathrm{s}}{0.2\text{ km s}^1}\right)\left(\frac{n}{10^3\text{ cm}^3}\right)^{1/2},$$
(3)
and the mass scale by the Jeans mass

$$M_J=\overline{\rho }\lambda _J^3(20.8M_{})\widehat{M}_J\left(\frac{c_\mathrm{s}}{0.2\text{ km s}^1}\right)\left(\frac{n}{10^3\text{ cm}^3}\right)^{1/2},$$
(4)
where $`\widehat{\lambda }_J`$ and $`\widehat{M}_J`$ are expressed in code units. We include the magnetic field $`B`$ via the pressure, using $`P=\rho c_\mathrm{s}^2=B^2/(8\pi )`$. This yields

$$B=(4.44\times 10^6G)\widehat{B}\left(\frac{c_\mathrm{s}}{0.2\text{ km s}^1}\right)\left(\frac{n}{10^3\text{ cm}^3}\right)^{1/2},$$
(5)
with $`\widehat{B}`$ again in code units.

## 3 Numerical Resolution Criterion

We must consider the resolution required to accurately follow gravitational collapse. To follow fragmentation in a grid-based, hydrodynamical simulation of self-gravitating gas, the criterion given by Truelove et al. (1997) holds. They studied fragmentation in self-gravitating, collapsing regions and found that the mass contained in one grid zone must remain significantly smaller than the local Jeans mass throughout the computation to accurately follow the fragmentation. Bate & Burkert (1997) found a similar criterion for particle methods. Applying it strictly would limit our simulations to the very first stages of collapse. We therefore only apply this criterion to the resolution of initial collapse. Thereafter, we only study the gross properties of collapsed cores, such as mass and location, but not their underresolved internal details.

The Truelove analysis does not include magnetic fields, which must also be sufficiently resolved to determine whether initial collapse occurs. Numerical diffusion can reduce the support provided by a static or dynamic magnetic field against gravitational collapse. Increasing the numerical resolution decreases the scale at which numerical diffusion acts. In this section we attempt to determine the resolution necessary to adequately resolve magnetic support against collapse in the presence of supersonic turbulence. Two regimes of field strength concern us. For strong, subcritical fields, the resolution should ensure that numerical diffusion remains unimportant even for the dense, shocked regions (subsection 3.1). For weaker, supercritical fields, the resolution should enable us to evolve MHD waves within the shocked regions (subsection 3.2).

### 3.1 Numerical Diffusion in Magnetostatic Configurations

We begin by considering magnetostatic support. Figure 1 demonstrates that numerical diffusion can dominate the behavior in this case. The left panel displays the peak density $`\rho _{max}`$ and magnetic field amplitude $`B_{max}`$ for the low and intermediate resolution undriven models $`l0s`$ and $`i0s`$, while the right panel contains the same quantities for the driven, but otherwise identical, models $`l1s`$ and $`i1s`$ (see Table 2). All these models have initial magnetic field strength sufficient to support the region, with $`M/M_{cr}=0.8`$.

Starting with a sinusoidal density perturbation in the undriven case, both the driven and undriven models first collapse into sheet structures. The dotted lines in the density plots show the density $`\rho _{\mathrm{sheet}}`$ corresponding to having all the mass in the box in a layer one zone thick. In a volume filled with otherwise unperturbed gas of initially uniform density, reaching peak densities $`\rho _{max}`$ higher than this threshold means that numerical diffusion across field lines must have begun, as can be seen in the top left panel of Figure 1. However, in driven, isothermal, supersonic turbulence, the peak densities in shocked regions scale with the Mach number $``$ as $`\rho _{pk}^2`$. Thus, $`\rho _{pk}`$ can reach values orders of magnitude higher than the mean density, easily exceeding $`\rho _{sheet}`$ even in the absence of diffusion.

The magnetic field provides an alternate diagnostic. In the absence of numerical diffusion, mass should be tied to the field lines running vertically through the cube, so the mass to flux ratio along any given field line should not change. In the lower panels, we plot the magnetic field strength $`B_{sup}(\rho )`$ required to support a region of density $`\rho _{max}`$ (thin lines) according to equation 2. The magnetic field starts out significantly stronger than $`B_{sup}`$, as the mass is insufficient for collapse to occur. If $`B_{max}`$ grows more slowly than $`B_{sup}`$, then density must be diffusing across field lines. If the two values cross, then numerical diffusion has allowed collapse to occur unphysically. This happens in the undriven models at $`t=4.7t_{\mathrm{ff}}`$, while in the driven models with their greater density contrasts, collapse already occurs at $`t=0.5t_{\mathrm{ff}}`$. We conclude that in these low and intermediate resolution models, the resolution is not sufficiently high for the magnetic field to follow the turbulence, especially in shocked regions.

We can use this example to derive a criterion for the resolution of a magnetostatically supported sheet. We can also confirm that we are using the correct numerical constant in the mass-to-flux criterion given by equation (2). We did a suite of models $``$ (see Table 1) varying the sound speed and the mass in the cube while holding the magnetic field strength constant, thus varying $`M/M_{cr}`$ and the number of zones in a Jeans length $`\lambda _J`$. For $`\lambda _J=1.0`$, we certainly cannot expect to get reliable results, as this Jeans length does not even fulfill the hydrodynamic criterion of Truelove et al. (1997). In Figure 2 we present the results. We find unphysical collapse occurring for physically supported regions until $`\lambda _J4.0`$ zones. We also find that at this resolution, a model with $`M/M_{cr}=1.1`$ collapses (thin line in first panel) while a model with $`M/M_{cr}=0.88`$ does not collapse, confirming equation 2. We conclude that for a self-gravitating magnetostatic sheet to be well resolved, its Jeans length must exceed four zones.

### 3.2 MHD Waves in High Density Regions

As we want to investigate whether MHD turbulence can prevent gravitationally unstable regions from collapsing, we have to be sure to resolve the MHD waves, which are the main agent in this mechanism. The same argument holds as for the criterion to prevent local collapse in the pure hydro case (Paper I): If the wavelength of the turbulent perturbation is smaller than the local Jeans length $`\lambda _J`$, stabilization should be possible, at least in principle. To sample a sinusoidal wave on a grid, we need at least four zones. Polygonal interpolation of a sine wave with evenly spaced supports would then yield an error of $`21`$%. This decreases to $`9`$% with six zones and $`5`$% using eight zones. We choose to set the minimum permitted Jeans length to six zones, admittedly somewhat arbitrarily.

In table 2 we list the local Jeans lengths for all model types and their resolution. Models of type $``$ start to be resolved at a resolution of $`N=256^3`$, whereas the ones of type $`𝒢`$ can already be regarded as resolved at $`N=128^3`$.

## 4 Results

In the previous section we considered under what circumstances numerical effects could allow unphysical gravitational collapse. In this section we now consider adequately resolved models in order to determine whether magnetized turbulence can prevent the collapse of regions that are not magnetostatically supported. We begin by demonstrating that supersonic turbulence does not cause a magnetostatically supported region to collapse, and then demonstrate that in the absence of magnetostatic support, MHD waves cannot completely prevent collapse, although they can retard it.

### 4.1 Magnetostatic Support

In a subcritical region with $`M<M_{cr}`$, the cloud is expected to collapse to a sheet, which in turn should be stable. Figure 3 shows the corresponding model $`𝒢i1s`$. These runs have been computed with a lower Mach number $`=5.0`$ in order to demonstrate the behavior of a magnetostatically supported cloud. The initially uniform magnetic field runs parallel to the $`z`$-axis. The field is strong enough to force significant anisotropy in the flow, although the dense sheets that form do not always lie perpendicular to the field lines as the driving can shift the sheets along the field lines without changing the mass-to-flux ratio. The sheets do not collapse further, because the shock waves cannot sweep gas across field lines and the cloud is initially supported magnetostatically.

Figure 4 demonstrates that this result is reasonably well resolved numerically. As in figure 1, we show peak densities and magnetic field strength for two models that differ only in their sound speeds, and thus by the number of zones in a Jeans wavelength $`\lambda _J`$. Whereas model $`i1s`$, with $`\lambda _J=3.2`$ zones does collapse, although it should be supported, model $`𝒢i1s`$, with $`\lambda _J=6.8`$ zones behaves as expected physically, just as our resolution criterion predicts. Note that the actual field strength in model $`𝒢i1s`$ always exceeds $`B_{sup}`$, the field strength necessary to support the region.

### 4.2 MHD Wave Support

A supercritical cloud with $`M>M_{cr}`$ is not magnetostatically supported, and could only be stabilized by MHD wave pressure, assuming ideal MHD. In this section we show that this appears to be insufficient to completely prevent gravitational collapse, although it can slow the process down.

#### 4.2.1 Morphology

In the upper panels of Figure 5 we compare the morphology of hydrodynamical ($`𝒟h1`$), weakly magnetized ($`h1w`$), and strongly magnetized ($`h1i`$) supercritical models at a resolution of $`256^3`$ zones. The figure presents two dimensional slices through the three dimensional simulation volume centered on the locations of the most massive clumps. In order to compare the models at similar stages of their evolution, we took the snapshots at a time when roughly 10% of the total mass had been accreted onto collapsing clumps. All three runs show well developed turbulence, rarefied regions, shocked regions, and at least one clump. However, model $`h1w`$, with $`M/M_{cr}=8.3`$, seems to contain more power on small scales than the pure hydro run, model $`𝒟h1`$ ($`M/M_{cr}=\mathrm{}`$). We discuss this more quantitatively below. In model $`h1i`$, with $`M/M_{cr}=1.8`$, the vertically oriented mean field (in the plane of the Figure) starts producing some anisotropy. This model represents a morphological transition from the pure hydrodynamical model $`𝒟h1`$ with completely randomly oriented motions to the magnetostatically supported model $`𝒢i1s`$, with its ordered strucures.

#### 4.2.2 Resolution Study

We must address the question of whether our magnetodynamic simulations are indeed well resolved. The parameters for models $``$ yield a global Jeans length of $`\lambda _J0.5`$, corresponding to a local, post-shock Jeans length of $`\lambda _J0.05`$ for isothermal shocks with Mach number $`10`$ (see Table 2). At $`N=128^3`$ zones, this results in a local Jeans length of only $`3.2`$ zones, but at $`256^3`$ the local Jeans length is $`6.4`$ zones, satisfying our resolution criterion. Instead of increasing the resolution, we increased the Jeans length in the models $`𝒢`$ discussed below in § 4.2.4 by increasing the sound speed. In these models we used a global Jeans length of $`\lambda _J1.1`$, corresponding to a local Jeans length of $`6.8`$ zones at $`N=128^3`$.

Figure 5 compares high resolution $`256^3`$ zone models with the corresponding lower resolution $`128^3`$ zone models in the same dynamical state, when 10% of the mass has been accreted onto cores. The lower resolution makes itself felt in broader shocks in all cases, so that the peak densities are lower than in the high-resolution runs. In the MHD models decreasing the resolution also leads to thicker collapsed sheets. Thus, unstable regions form at later times in both the MHD and hydrodynamical cases, as can be seen in Figure 6. There, we show the core mass accretion history for three weakly magnetized models varying only in their resolution $`l1w`$ ($`64^3`$), $`i1w`$ ($`128^3`$) and $`h1w`$ ($`256^3`$) (lower panel) and the corresponding hydrodynamical models ($`𝒟l1`$, $`𝒟i1`$, $`𝒟h1`$) (upper panel). Cores were determined using the modified CLUMPFIND algorithm of Williams et al. (1994) as described in Paper I.

Collapse occurs in both cases at all resolutions. However, increasing the resolution makes itself felt in different ways in hydrodynamical and MHD models. In the hydrodynamical case, higher resolution results in thinner shocks and thus higher peak densities. These higher density peaks form cores with deeper potential wells that accrete more mass and are more stable against disruption. If we increase the resolution in the MHD models, on the other hand, we can better follow short wavelength MHD waves, which appear to be able to delay collapse, although not to prevent it.

The turbulent formation of gravitationally condensed regions via shock interactions is a highly stochastic process. As in Paper I, we demonstrate this by choosing different random realizations of the driving velocity field with the same characteristic wavelengths. Figure 7 shows the core mass accretion history for a pure hydrodynamical model set and a MHD set. In the upper panel, we plotted the models $`𝒟l1`$, $`𝒟i1`$, $`𝒟h1`$, where the low-resolution model $`𝒟l1`$ has been repeated multiple times (solid thin lines). The dotted thick line denotes the average of these runs. We find that resolution effects are exceeded by statistical variations caused by random variations of the driving fields.

In the MHD case (lower panel) the thickness of the lines stands for the strength of the field, expressed in terms of the ratio $`M/M_{cr}`$. Dotted lines denote low-resolution runs computed with varying driving velocity fields, as in the upper panel. The high-resolution run with $`M/M_{cr}=1.8`$ was stopped at $`t=1.0t_{\mathrm{ff}}`$ because the Alfvén timestep became prohibitively small. Increasing the resolution makes itself felt for the runs with stronger fields in the same way as for the ones with weak fields: The higher the resolution, the better the small-scale MHD waves are resolved, and thus the slower the collapse. Collapse does always occur, however.

#### 4.2.3 Core Distribution

Although MHD waves cannot prevent local collapse entirely, the resulting collapse appears qualitatively different from collapse in the hydrodynamic case with corresponding global $`\lambda _J`$ and driving strength. In the hydrodynamical case driven at $`\lambda _D>\lambda _J`$, shocks are widely separated and sweep up substantial mass, producing isolated clusters of cores. In the presence of a weak supercritical field, the shock structure appears to have more small-scale structure, resulting in cores distributed more uniformly across the simulation volume, as shown in the middle panel of Figure 8. In fact, the weakly magnetized model driven on large scales with $`\lambda _D>\lambda _J`$ rather more resembles the hydro model driven on small scales with $`\lambda _D<\lambda _J`$ shown in Figure 11 of Paper I.

Figures 9 and 10 try to quantify this difference. Figure 9 shows the histogram of core distances for each panel of Figure 8, multiplied with the mean core mass. A clustered ensemble of high-mass cores should result in a peaked distribution at small distances, whereas a spread-out ensemble of low-mass cores should have a broader distribution. Figure 9 hints at such a behaviour, although we are well aware of the fact, that the statistics barely suffice. Nevertheless, the distribution for the magnetized runs is shifted to larger radii. Note that the total mass in the cores found is within $`10`$% the same for all three models. Figure 10 shows the weighted means of core distances, with their standard deviations as error bars. Again, we see a slight shift towards larger separations, suggesting a more uniform distribution of cores in the MHD cases, although larger simulations with greater core numbers will be needed to confirm this result.

#### 4.2.4 Energy Distribution

Further evidence for a qualitative difference between hydrodynamic and MHD collapse comes from Figure 11. Here, we show the time evolution of the ratio of kinetic to potential energy, decomposed into contributions from four spatial scales. (The time resolution is somewhat coarse as these models with $`256^3`$ zones were only dumped every $`t=0.3t_{\mathrm{ff}}`$.) Whereas the hydrodynamic model ($`𝒟h1`$) driven at $`k=12`$ collapses within less than $`0.5t_{\mathrm{ff}}`$, the weakly magnetized model $`h1w`$ ($`M/M_{crit}=8.3`$) is supported until $`t1t_{\mathrm{ff}}`$ before also collapsing. This is at least qualitatively similar to the behavior of model $`𝒟h3`$, which is driven at wavenumbers $`k=78`$ and thus has a denser network of shocks. On the other hand, Figure 11 suggests that the stronger-field model $`h1i`$ ($`M/M_{cr}=1.8`$) collapses even more thoroughly than its hydrodynamical counterpart due to the ordering influence of the strong mean field.

We use the model series $`𝒢`$ to follow the collapse to later times and with more frequent time sampling. These $`128^3`$ models still resolve the Jeans length in the densest clumps, as discussed in § 4.2.2. We reduced the Mach number to $`=5`$ in order to maintain an energy input comparable to models $``$. This decreases the rms post-shock density, so we actually expect the $`𝒢`$ series to form cores with somewhat more difficulty than the $``$ series. Figure 12 shows the ratio of kinetic to potential energy of the magnetized models in the $`𝒢`$ series. (The hydrodynamical model $`𝒢i1`$ collapses within half a free-fall time after gravity is turned on at $`t=0.0`$.) With increasing field strength (models $`𝒢i1w`$ and $`𝒢i1i`$), the collapse is delayed, but never prevented. This even applies for $`𝒢i1m`$, where the field strength is only marginally supercritical ($`M/M_{crit}=1.1`$). In this model, bound cores form, but are then destroyed by passing shocks, probably for the unphysical reason that they cannot continue collapsing to sizes smaller than a few zones (see discussion in Paper I). Increasing the field strength further leads to model $`𝒢i1s`$ ($`M/M_{crit}=0.8`$), where the field supports the cloud magnetostatically, and gravitationally bound cores do not form.

#### 4.2.5 Energy Spectra

We next examine Fourier spectra of the energy. Figure 13 presents the spectra of kinetic, potential and magnetic energies at times $`t=0.0`$ and $`t=1.5t_{\mathrm{ff}}`$ of the high-resolution models $`𝒟h1`$ (hydrodynamic), $`h1w`$ ($`M/M_{cr}=8.3`$) and $`h1i`$ ($`M/M_{cr}=1.8`$), all driven at wavenumber $`k=12`$. The spectra at time $`t=0.0`$ represent fully developed turbulence just before gravity is switched on. Here we find another reason for the fast collapse of the more strongly magnetized model $`h1i`$. The density enhancements due to shock interactions are larger for models $`𝒟h1`$ and $`h1i`$ than for $`h1w`$, as can be inferred from comparing the potential energies at $`t=0.0`$. Although still supercritical, the field in model $`h1i`$ is already strong enough to suppress motions perpendicular to the mean field, so that the field strength perpendicular to its initial mean direction is small, while strong shocks/waves can be formed parallel to the mean field, as observed as well by Smith et al. (2000). Thus, somewhat surprisingly, the density enhancements are larger for the stronger-field model $`h1i`$ than for the weak-field model $`h1w`$, leading to earlier collapse, as seen in Figure 11. In model $`h1w`$ the weaker, more turbulent magnetic field produces a more isotropic magnetic pressure that cushions and broadens the shocks, thus decreasing the density enhancements and delaying collapse.

We illustrate this effect in Figure 14, where we plot the $`x`$\- $`y`$\- and $`z`$-components of the magnetic energy against time for the magnetized models $`h1w`$ and $`h1i`$. In the weak-field case, the field is quickly tangled by the flow, so that it has no preferred direction by $`t=0.0`$, when gravity is switched on. The magnetic energy, and so the magnetic pressure, is isotropic. There is no secular increase of any component with time, supporting the picture of local collapse: In a globally collapsing environment, the magnetic field lines would follow the global gas flow and lead to a noticeable increase of magnetic energy. In the stronger-field case, on the other hand, the flow is dominated by the mean field oriented along the $`z`$-axis. The field allows matter to move more freely in the parallel than in the perpendicular direction. Matter thus collapses preferentially along field lines first, and then globally collapses.

After one free-fall time all the models have collapsed (Figure 13). Note that, as in Paper I, $`E_{pot}(k)>E_{kin}(k)`$ for all $`k`$ does not necessarily mean that the model becomes globally unstable. With increasing time, $`E_{pot}(k)`$ becomes constant for all $`k`$ just because this is the Fourier transform of a $`\delta `$-function, signifying that local collapse has produced point-like high-density cores. A similar argument applies for the kinetic energy: The flat spectrum stems from local concentrations of kinetic energy around the collapsing regions. Here, the spectral analysis no longer yields information on global stability.

#### 4.2.6 Conclusions

We conclude that the delay of local collapse seen in our magnetized simulations is caused mainly by weakly magnetized turbulence acting as a more or less isotropic pressure in the gas, decreasing density enhancements due to shock interactions. We feel justified in claiming that magnetic fields, as long as they do not provide magnetostatic support, can not prevent local collapse, even in the presence of supersonic turbulence.

We note that once bound cores form, we take this as evidence for local collapse, although subsequent shock interactions may destroy these cores again. In a real cloud, ambipolar diffusion would set in at the length scale of transient cores, so that any internal turbulence would be quickly dissipated, allowing further collapse, as discussed in Paper I.

A large-scale driver, such as interacting supernova remnants or galactic shear, together with magnetic fields, seems to act like a driver with a smaller effective scale in the sense that both yield a more uniform core distribution and a somewhat slower collapse rate. In weakly magnetized turbulence, a more or less isotropic magnetic pressure reduces the density enhancements behind shocks and thus slows down the process of isolated collapse. In strongly magnetized turbulence, however, the mean magnetic field dominates. The magnetic pressure is not isotropic any more, so the shocks perpendicular to the mean field direction cause high enough density enhancements for the regions to collapse within a free-fall time.

Thus, for small field strength, the effective additional pressure may be represented by a simple pressure term. However, in the regime of field strength interesting for molecular clouds, the field, although supercritical, is strong enough to result in an anisotropic magnetic pressure. Magnetic turbulence is an all-scale non-isotropic phenomenon, and the compression and perturbations on large scales make the cloud finally collapse.

## 5 Summary

In this paper, we investigated whether magnetized turbulence can prevent collapse of a Jeans-unstable region. From our high-resolution simulations we conclude that

1. In order to resolve self-gravitating MHD turbulence using a grid-based method like ZEUS-3D, the local Jeans length should not fall short of at least four grid zones for magnetostatic support and six grid zones for magnetodynamic support.
2. Local collapse cannot be prevented by magnetized turbulence in the absence of mean-field support. Strong local density enhancements due to shock interactions start collapsing at once.
3. However, the magnetic fields do delay local collapse by decreasing local density enhancements via magnetic pressure behind shocks.
4. Weakly magnetized turbulence appears qualitatively similar to hydrodynamic turbulence driven on a slightly smaller scale, while stronger fields close to but under the value for magnetostatic support tend to organize the flow into sheets and allow more clustered collapse.
5. The strength and wavelength of turbulent driving governs the behaviour of the cloud, overshadowing the effects of magnetic fields that do not provide magnetostatic support.
6. MHD turbulence can not prevent local collapse for much longer than a global free-fall time. Stars begin to form at a low rate as soon as local density enhancements contract. This result favours a dynamical picture of molecular clouds being a transient feature in the ISM (Ballesteros-Paredes et al. 1999, Elmegreen 2000) rather than living for many free-fall times.

We thank A. Burkert, H. Lesch, S. Phleps, E. Vázquez-Semadeni and E. Zweibel for valuable discussions. Thanks go to the referee, whose criticism helped clarify important points in the paper. M-MML is partially supported by an NSF CAREER fellowship, grant number AST 99-85392. Computations presented here were performed on SGI Origin 2000 machines of the Rechenzentrum Garching of the Max-Planck-Gesellschaft, the National Center for Supercomputing Applications (NCSA), and the Hayden Planetarium. ZEUS was used by courtesy of the Laboratory for Computational Astrophysics at the NCSA. This research has made use of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# Expansion-induced contribution to the precession of binary orbits

## Abstract

We point out the existence of new effects of global spacetime expansion on local binary systems. In addition to a possible change of orbital size, there is a contribution to the precession of elliptic orbits, to be added to the well-known general relativistic effect in static spacetimes, and the eccentricity can change. Our model calculations are done using geodesics in a McVittie metric, representing a localized system in an asymptotically Robertson-Walker spacetime; we give a few numerical estimates for that case, and briefly comment on ways in which the model should be improved.

The issue of whether the global cosmological expansion affects local gravitating systems, such as planetary systems and galaxies, has been studied for a long time (for recent work, and discussions of previous results, see e.g. Refs. ). The specific questions asked in the literature focus mostly on the extent to which local systems expand, and on the form and magnitude of the corrections to the effective forces felt by orbiting test bodies; while the approaches used vary, sometimes in conceptually important ways, the consensus is that there is an effect in principle, but in practice it is exceedingly, undetectably small. To a first approximation, a reasonable, physically motivated point of view, as expressed by Misner, Thorne, and Wheeler , for the case of galaxies, is that they are like rigid pennies attached to the expanding balloon representing the universe, and they do not themselves expand.

While we agree with this general conclusion, at least in terms of currently feasible observations, an increase in orbit size is not the only effect the cosmological expansion can have on local systems, nor the only cumulative one. In the search for such effects described in this paper, we use as a model spacetime an actual solution of Einstein’s equation which represents a gravitating body in an asymptotically expanding universe; restrict our attention to nearly Newtonian orbits in the weak-gravity, slow-expansion approximation; and consider the time variation of the parameters characterizing the corresponding Keplerian ellipses, their size, perihelion angle, and eccentricity. As we will show, this approach leads us to predict changes in all of those parameters as a consequence of global expansion.

Before we proceed, however, we will make a few comments on the systems involved. It is well known that Keplerian elliptical orbits in Newtonian gravity do not precess (the $`1/r`$ potential is one of two that yield closed orbits, together with the harmonic oscillator, as Bertrand’s theorem states) but the corresponding ones in general relativity do, while the size and eccentricity of the general relativistic orbits are constant. These results hold for the case of a single, spherically symmetric, isolated center of attraction, with a test body and no other matter around it. Corrections due to other orbiting bodies can be calculated, and they may affect all parameters of the orbits, but as far as the global expansion is concerned, under the above conditions there can be no local effects due to cosmological expansion in general relativity, regardless of the rate at which the rest of the universe expands, because Birkhoff’s theorem then concludes that locally the spacetime must be the static, Schwarzschild solution (Einstein and Straus  also showed this explicitly in one class of models, without using Birkhoff’s theorem). The cosmological expansion can be felt by a local system only if (i) The situation is not spherically symmetric, and/or (ii) There is matter (some gas, dust, or dark matter, say) within the local system, or a cosmological constant. This is one of the many ways in which one can see that general relativity obeys the spirit of Mach’s principle, but only up to a point (see, e.g., Ref. ); the global behavior of distant objects can affect the orbital dynamics, but the dynamical equations are local.

In a realistic model of a local system, both deviations from spherical symmetry and extra matter will have an effect; in the solar system, for example, the presence of other stars, solar multipoles such as oblateness, minor objects, interplanetary matter and/or the solar wind and radiation, are ways in which this happens. Studying the way the global expansion affects the local system then means: (i) Understanding how the dynamics of surrounding non-symmetric matter and/or the local internal environment of the system depends on the global behavior of the universe, and (ii) Calculating the effect of the local environment on the dynamics of the system of interest. Notice that, for a local system that is considered as part of a larger one, such as a planetary system inside a galaxy, the first part amounts to assuming that a similar problem has been solved one level up; for example, the spacetime metric for a planetary system should not be asymptotic to a Robertson-Walker metric expanding at the Hubble rate, but to a metric that is appropriate inside the galaxy, expanding, if at all, at a rate previously found for that system. The problem thus becomes somewhat involved, but it seems to us that this is the only way to avoid sneaking in an assumption about the local environment’s expansion as part of the setup, which has occasionally been done. On the other hand, once it has been formulated in this way, the problem is broken down into parts that can be separately investigated, and labelled by the structure scale under consideration (galactic, planetary, …), mechanism by which the effect may be transmitted, and type of effect.

Here, as a start, we will limit ourselves to considering a simplified model to show the effects that may arise. Spherical symmetry will be preserved, and the mechanism responsible for transmitting environment effects will be the presence of matter in the local system. In our numerical estimates we will use orbital parameters relevant for a couple of known planets and galaxies, and the Hubble constant to obtain the environment expansion rate; much more work will be needed on several aspects of this problem, including the type of matter and spacetime model used and the role of anisotropy, before we can claim to have a reasonable understanding of it along the lines described above.

While spherically symmetric, vacuum, asymptotically flat spacetimes and homogeneous, isotropic cosmological ones with fluid matter or a cosmological constant can be easily treated in general relativity and give rise, respectively, to the Schwarzschild solution and the Friedmann-Robertson-Walker or de Sitter spacetimes, solutions representing an isolated massive object embedded in an expanding universe are more difficult to get, and fully explicit forms are not usually given (see, for example, the general discussion in Ref. ). The known metrics representing isolated bodies in asymptotically FRW spacetimes are mainly of three types: McVittie , Einstein-Straus , and Ferraris-Francaviglia-Spallicci . In this paper, we will use McVittie metrics, although an explicit solution is given in Ref.  where a Schwarzschild and a FRW solution are matched smoothly across an extended region of spacetime represented by a transition metric, where the orbit of a test particle may also be analysed.

The (asymptotically) spatially flat McVittie model with $`k=0`$, in the isotropic coordinate form used by Hogan , has line element

$$\mathrm{d}s^2=\left(\frac{1M\mathrm{e}^{\beta /2}/2r}{1+M\mathrm{e}^{\beta /2}/2r}\right)^2\mathrm{d}t^2+\mathrm{e}^{\beta (t)}\left(1+\frac{M\mathrm{e}^{\beta /2}}{2r}\right)^4(\mathrm{d}r^2+r^2\mathrm{d}\mathrm{\Omega }^2).$$
(1)
Here, the constant $`M`$ is interpreted as the “mass at the singularity,” and $`\beta (t)`$ represents the (asymptotic) expansion rate of the universe. We will not give an explicit expression for $`\beta (t)`$, since it is related to the specific equation of state of the fluid, but will think of it instead in terms of an expansion $`\beta (t)=\beta _0+\dot{\beta }_0(tt_0)+𝒪((tt_0)^2)`$, around the present time $`t_0`$, say, where the first few coefficients can be fitted phenomenologically.

The above line element has the advantage that, as $`r\mathrm{}`$, the $`t`$ = const hypersurfaces become the surfaces of homogeneity of a Robertson-Walker model. For our purposes however, we find it more convenient to work with a radial coordinate with a geometrical meaning tied to the area of the corresponding 2-sphere, and use $`R:=r\mathrm{e}^{\beta /2}(1+GM\mathrm{e}^{\beta /2}/2r)^2`$ in the interval $`R(2M,\mathrm{})`$ . We can then rewrite the line element as

$$\mathrm{d}s^2=\left[f(R)+\left(\frac{R\dot{\beta }}{2c}\right)^2\right]c^2\mathrm{d}t^2+\frac{2}{\sqrt{f(R)}}\frac{R\dot{\beta }}{2c}c\mathrm{d}t\mathrm{d}R+\frac{\mathrm{d}R^2}{f(R)}+R^2\mathrm{d}\mathrm{\Omega }^2,$$
(2)
where $`f(R):=12GM/(c^2R)`$, and we have restored all $`c`$’s and $`G`$’s.

The use of the McVittie metric as a model for studying local systems in FRW universes has been the subject of criticism , but in a series of papers  Nolan made a good argument for the $`k=0`$ model by studying its global properties in detail. Nolan also showed that, in an appropriate sense, those solutions are the unique ones representing an isolated mass surrounded by a shear-free perfect fluid in a spatially flat FRW universe; the shear-free condition leads to special properties that we would not expect a black hole to have , but the metric does provide us with a viable explicit model for our purposes.

Because of the spherical symmetry, we consider as usual our orbits to lie in the $`\mathrm{sin}\theta =1`$ plane, and drop $`\theta `$ from the whole treatment from now on. For a particle of mass $`m`$ moving in the metric (2), one way of deriving the equations of motion is to use the Hamiltonian approach. The canonical momenta obtained from the action $`S[x]=(m/2)d\tau g_{\mu \nu }\dot{x}^\mu \dot{x}^\nu `$ are $`p_\mu =mg_{\mu \nu }\dot{x}^\nu `$, where the overdot on a dynamical variable, $`\dot{t}`$, $`\dot{R}`$, or $`\dot{\varphi }`$, denotes a derivative with respect to proper time $`\tau `$ along a trajectory (whereas $`\dot{\beta }`$ just indicates the derivative of the known function $`\beta (t)`$ with respect to its argument).

The test particle Hamiltonian $`=(2m)^1g^{\mu \nu }p_\mu p_\nu `$ in this case becomes

$$=\frac{1}{2m}\left(\frac{p_t^2}{f(R)c^2}+\frac{2}{\sqrt{f(R)}}\frac{R\dot{\beta }}{2c}\frac{p_tp_R}{c}+\left[f(R)\left(\frac{R\dot{\beta }}{2c}\right)^2\right]p_R^2+\frac{p_\varphi ^2}{R^2}\right),$$

supplemented by the constraint that $`g^{\mu \nu }p_\mu p_\nu =(mc)^2`$, or $`=\frac{1}{2}mc^2`$.

Angular momentum is conserved, $`p_\varphi =L`$, and as is common in celestial mechanics orbit problems, we reexpress $`R`$ in terms of $`u:=1/R`$ in the equations of motion. If we perturbatively expand those equations in powers of $`1/c`$ and $`\dot{\beta }`$, zeroth order terms correspond to a particle in Newtonian gravitation; neglecting all $`\dot{\beta }`$ terms corresponds to doing calculations in a static spacetime, and we recover the well-known Schwarzschild results; the expansion effects we are interested in arise when leading order terms in $`\dot{\beta }`$ are kept.

The radial equation can be recast into an orbit equation in terms of $`u(\varphi )`$,

$$u^{\prime \prime }+u=\frac{GMm^2}{L^2}+\frac{3GMu^2}{c^2}\frac{\dot{\beta }^2m^2}{4L^2u^3}\frac{GM\dot{\beta }m}{2c^2Lu}u^{}.$$

If we look for a solution of the form $`u(\varphi )=u_0+\xi (\varphi )`$, where $`u_0=GMm^2/L^2`$ corresponds to a Newtonian circular orbit and $`\xi (\varphi )`$ is a small perturbation (see, e.g., Refs. ), the linearized equation for the deviation $`\xi (\varphi )`$ is

$$\xi ^{\prime \prime }=\frac{3GMu_0^2}{c^2}\frac{\dot{\beta }^2m^2}{4L^2u_0^3}\left(1\frac{6GMu_0}{c^2}\frac{3\dot{\beta }^2m^2}{4L^2u_0^4}\right)\xi \frac{GM\dot{\beta }m}{2c^2Lu_0}\xi ^{}.$$
(3)
It should be noted that for a general McVittie metric, $`\dot{\beta }=\dot{\beta }(t(\varphi ))`$, so, to simplify the solution of the equation, we will now limit ourselves to considering the case $`\ddot{\beta }=0`$; in terms of comoving coordinates as in Eq. (1), this is the linear expansion case, although the metric for this model has a timelike Killing vector field, as can be seen from the fact that the line element (2) becomes time-independent (in fact, it becomes equivalent to a Schwarzschild-de Sitter solution ). Equation (3) is then of the form $`\xi ^{\prime \prime }=f_\mathrm{c}\omega ^2\xi +f_\mathrm{d}\xi ^{}`$, which can be viewed as that of an oscillator with a constant force $`f_\mathrm{c}`$, modified by a $`\dot{\beta }`$ term which affects the equilibrium displacement of the oscillator and thus the geodesic orbit size, but it does not lead to orbit expansion (when $`\ddot{\beta }=0`$); the frequency $`\omega `$ is also modified by a $`\dot{\beta }`$ term, which contributes to the orbital precession, and, interestingly, there is a damping term $`f_\mathrm{d}`$, which leads to a change in orbital eccentricity. The equation admits a general solution of the form $`\xi (\varphi )=𝒞+\mathrm{e}^\varphi \mathrm{cos}(𝒜\varphi \varphi _0)`$, depending on arbitrary constants $``$ and $`\varphi _0`$, and where

$$𝒜:=\sqrt{\omega ^2\frac{1}{4}f_\mathrm{d}^2}1\frac{3GMu_0}{c^2}\frac{3\dot{\beta }_0^2m^2}{8L^2u_0^4},$$

$$:=\frac{1}{2}f_\mathrm{d}=\frac{GM\dot{\beta }_0m}{4c^2Lu_0},$$

$$𝒞:=\frac{f_\mathrm{c}}{\omega ^2}\frac{3GMu_0^2}{c^2}\frac{\dot{\beta }_0^2m^2}{4L^2u_0^3},$$

from which, defining $`u_0^{}=u_0+𝒞`$ and $`\epsilon _0=/(u_0+𝒞)`$, we finally get for the inverse radial coordinate

$$u(\varphi )=u_0^{}\left[1+\epsilon _0\mathrm{e}^\varphi \mathrm{cos}(𝒜\varphi \varphi _0)\right].$$
(4)
The angle $`\sigma `$ by which the orbit (4) precesses during each revolution can be found from the fact that the change $`\mathrm{\Delta }\varphi `$ such that the argument of the cosine increases by $`2\pi `$ is $`2\pi +\sigma =2\pi /𝒜`$; i.e., for $`𝒜1`$, $`\sigma 2\pi (1𝒜)`$, or

$$\sigma \frac{6\pi GM}{a(1\epsilon ^2)c^2}+\frac{3\pi }{4}\frac{\dot{\beta }_0^2a^3(1\epsilon ^2)^3}{GM},$$
(5)
where $`\epsilon =\epsilon _0\mathrm{e}^\varphi `$ is the orbital eccentricity, and we have used the fact that for a Keplerian ellipse $`1/u_0=a(1\epsilon ^2)`$, with $`a`$ the semi-major axis. The first term in (5) is the Schwarzschild contribution, $`\sigma _0`$, as expected; the second term is the expansion-induced one, and has the same sign as the first term. The eccentricity also changes, by a fractional amount per revolution $`\mathrm{\Delta }\epsilon /\epsilon =\mathrm{e}^{2\pi }1`$, or

$$\frac{\mathrm{\Delta }\epsilon }{\epsilon }\frac{\pi GM\dot{\beta }_0m}{2c^2Lu_0}=\frac{\pi \dot{\beta }_0}{2c^2}\sqrt{GMa(1\epsilon ^2)}.$$
(6)
Equations (5) and (6) are our main results; they show that, even in a model with no orbit expansion (this may just be a feature of the $`\ddot{\beta }=0`$ McVittie models used), global expansion has an effect on periastron precession and eccentricity.

We can get a first crude estimate of the magnitude of these effects by evaluating $`\sigma `$ and $`\mathrm{\Delta }\epsilon /\epsilon `$ for a few known systems. For the solar planets Mercury and Pluto, even using a Hubble parameter value $`H_0`$ around 75 km/s/Mpc for $`\dot{\beta }_0`$, which is certainly a vast overestimate, one obtains for the relative size of the expansion-induced precession term in (5), values $`\sigma _\beta /\sigma _03\times 10^{21}`$ and $`7\times 10^{13}`$, respectively, and for the relative eccentricity change in (6), $`\mathrm{\Delta }\epsilon /\epsilon 1.5\times 10^{21}`$/rev and $`1.5\times 10^{20}`$/rev, respectively. From the $`M`$ and $`a`$ dependence of these quantities we see that none of the strong gravitational field situations in which other general relativistic effects are studied, such as the binary pulsar  or stellar motion around the galactic center  can give usably large values, and we conclude that all expansion effects are negligible for the astrophysics of stellar systems.

What about gravitationally bound pairs of galaxies? After all, as one might expect from cosmological effects (and in contrast to the Schwarzschild ones), both $`\sigma _\beta /\sigma _0`$ and $`\mathrm{\Delta }\epsilon /\epsilon `$ grow with orbit size $`a`$. We can get an idea if we apply equations (5) and (6) to the Large Magellanic Cloud, at $`a50`$ kpc from the Milky Way, with $`M10^{11}M_{}`$, which gives $`\sigma _\beta /\sigma _00.34`$ and $`\mathrm{\Delta }\epsilon /\epsilon 8\times 10^{11}`$/rev, both significantly larger than the previous values, by many orders of magnitude in the case of $`\sigma _\beta /\sigma _0`$, in a regime where the use of $`H_0`$ for $`\dot{\beta }`$, although probably still not appropriate, comes closer to being realistic.

So, if some of the predicted numbers are not necessarily small, are our expansion effects cosmologically relevant? The problem here lies with the time scales involved in galactic dynamics, which make any direct observation of time variations impossible. Still, the fact that the Hubble distance-redshift relation shows local deviations on small cosmological scales is an indication that on those scales local interactions and global expansion effects coexist. Similar ideas have already motivated, e.g., work on the effect of global expansion on the formation and evolution of clusters of galaxies , and more recently on the effect of inhomogeneities on the overall expansion  or of a cosmological constant on local dynamics , and one of the potential benefits of this line of research is the possibility of gaining a new handle on dark matter. The present approach can be seen as a model for the onset of this situation with simple, binary systems, and as we have seen it can provide predictions for such systems. The model can be improved to include more realistic sources, diffuse matter and anisotropy, as well as possible multipole moments of the orbiting mass . It may then be possible to identify statistical effects which can be measured by observing large numbers of galaxy pairs, or early universe effects which left and imprint on later evolution.

Acknowledgement

We would like to thank Kumar Bhatt for many stimulating conversations, without which this project might have never been started.


text:
# Why the 𝒇_𝟎(𝟗𝟖𝟎) is mostly 𝒔𝒔̄

## Abstract

We exploit the $`W`$-emission process to study the measured weak decay $`D_s^+f_0(980)\pi ^+`$ . We conclude that the scalar $`f_0(980)`$ meson is mostly $`s\overline{s}`$, which is supported by different model studies.

Of recent interest is the question of whether the scalar meson $`f_0(980)`$ consists mostly of nonstrange or of strange quarks .

In the naive $`q\overline{q}`$ picture, one can treat the isoscalar $`f_0(980)`$ as either mostly nonstrange and nearly degenerate with the isovector $`a_0`$(980) — like the vector pair $`\omega `$(782) and $`\rho `$(770) — or take the $`f_0(980)`$ as mostly $`s\overline{s}`$, in analogy to the close by, almost pure $`s\overline{s}`$ vector $`\varphi `$(1020). We prefer the latter interpretation, which is also obtained in the unitarized quark/meson model of Ref.  (see below). Another possibility is the four-quark scheme , in which both the $`f_0(980)`$ and $`a_0`$(980) have substantial nonstrange and strange components of comparable magnitude. We may study the prior question from the recently measured weak decay rate

$$\mathrm{\Gamma }\left(D_s^+f_0(980)\pi ^+\text{ }\right)=(2.39\pm 1.06)10^{14}GeV,$$
(1)
due to the branching ratio of $`(1.8\pm 0.8)`$ percent and the lifetime of $`(0.496\pm 0.01)10^{12}`$ $`s`$.

Given the rate of equation (1), we may express the decay amplitude magnitude as ($`p=732MeV`$ and $`m_{D_s^+}=1.9686GeV`$)

$$\left|\left(D_s^+f_0(980)\pi ^+\text{ }\right)\right|=m_{D_s^+}\sqrt{\frac{8\pi \mathrm{\Gamma }}{p}}=(1.78\pm 0.40)10^6GeV.$$
(2)
Assuming that the $`f_0(980)`$ is pure $`s\overline{s}`$, we may account for the experimental amplitude of equation (2) via the $`W^+`$-emission graph of figure (1).

This predicts, using the values given in  for the Fermi coupling constant, $`G_F`$, the cosine of the phase from the CKM weak mixing matrix, $`c_1`$, and the pion decay constant, $`f_\pi `$,

$$\left|\left(D_s^+f_0(980)\pi ^+\text{ }\right)\right|_{W^+\text{-emission}}=\frac{G_Fc_1^2}{\sqrt{2}}f_\pi \left(m_{D_s}^2m_{f_0}^2\right)\mathrm{\hspace{0.33em}2.13}10^6GeV.$$
(3)
We use here for the Fermi coupling constant, $`G_F`$, the value given in , as well as for $`c_1`$, the cosine of the Cabibbo angle, $`\mathrm{\Theta }_1`$ ($`\mathrm{\Theta }_112.8^{}`$, i.e., $`V_{ud}=0.975`$). We moreover assume $`V_{cs}=V_{ud}`$ at the accuracy of this paper. For the pion formfactor, $`f_\pi `$, we take the value $`93MeV`$.

Not only is equation (3) compatible with the data of equation (2), but the analogous $`\mathrm{\Delta }I=3/2`$ $`W^+`$-emission process $`K^+\pi ^+\pi ^0`$ is, too :

$$\left|\left(K^+\pi ^+\pi ^0\text{ }\right)\right|_{W^+\text{-emission}}=\frac{G_Fs_1c_1}{2\sqrt{2}}f_\pi \left(m_K^2m_\pi ^2\right)\mathrm{\hspace{0.33em}1.9}10^8GeV;$$
(4)
the data being  $`(1.834\pm 0.007)10^8GeV`$.

By comparing formula (2) with the prediction of formula (3), we could infer a small scalar mixing of $`n\overline{n}`$ and $`s\overline{s}`$, which is close to the theoretical estimates. Since $`\sigma |f_0=0`$, a mixing angle, $`\varphi _s`$, of about $`14^{}`$  to $`20^{}`$ , predicts an amplitude ratio

$$\frac{\left|\left(D_s^+f_0(980)\pi ^+\text{ }\right)\right|}{\frac{G_Fc_1^2}{\sqrt{2}}f_\pi \left(m_{D_s}^2m_{f_0}^2\right)\mathrm{cos}\left(\varphi _s\right)}\mathrm{\hspace{0.33em}0.9}\pm 0.2,$$
(5)
which certainly shows good agreement between the experimental decay amplitude ( 2) and the $`W^+`$-emission result ( 3).

As for earlier interpretations of $`f_0(980)`$ being mostly $`s\overline{s}`$, we refer to the 1982 and 1986 papers of the present authors,  and . In the latter case the nonstrange and strange low-lying scalar mesons are modelled jointly with the pseudoscalars and vectors including mesons which contain $`c`$ and $`b`$ quarks. The resulting four-parameter model, which fits very well low-energy $`S`$\- and $`P`$-wave meson-meson scattering cross sections as well as the bound state and resonance positions, predicts two complete scalar-meson nonets, one below and one above 1 GeV.

Also we note that the mostly $`s\overline{s}`$ structure of the $`f_0(980)`$ is stressed in the more recent work of Törnqvist  and of Törnqvist and Roos , the latter of which focuses on the nonstrange $`\sigma (4001200)`$. Let us recall here the DM2 data  measuring the nonstrange large $`\sigma `$ bump in the $`\pi \pi `$ mass distribution for $`J/\psi \omega \pi \pi `$, but only a small $`f_0(980)`$ pimple, which shows the smallness of the nonstrange content of the $`f_0(980)`$ meson.

Also the recent Fermilab E791 collaboration  suggests that the $`f_0(980)`$ is mostly $`s\overline{s}`$, but that the $`f_0(1370)`$ is mostly nonstrange — as hinted by the naive quark model — since $`D_s^+K^+K^{}\pi ^+`$ is not seen. It should be noted that the same is predicted in the unitarized meson model of Ref.  (see also Ref. ), contrary to e.g. Refs.  and , which interpret the $`f_0(1370)`$ as mostly $`s\overline{s}`$. This discrepancy between two unitarized models is quite striking and deserves a little more attention.

There are clearly too many scalar mesons below roughly 1.5 GeV than those needed to constitute one scalar nonet. Therefore, models that accept the light ($`<`$ 1GeV) scalars as genuine mesonic states with some kind of quark substructure — irrespective of the precise configuration — and not just as two-meson resonance effects due to strong $`t`$-channel exchanges and nearby thresholds (see e.g. Ref. ), normally consider the scalars above 1 GeV as excited states. However, in the unitarized models of Refs.  and  all scalar mesons below 1.5 GeV originate in one and only bare scalar $`q\overline{q}`$ nonet, somewhere in the mass region 1.3–1.5 GeV. Here, “bare” refers to the model situation where the coupling to two-meson states is switched off. Now, when this coupling is set to the model value which fits the data, the bare scalar spectrum gets deformed and even extra states show up, as for example the $`f_0(4001200)`$ or $`\sigma `$ meson, the $`a_0(980)`$, and, of course, the $`f_0(980)`$ (see Ref.  for more details).

In the model of Ref. , the extra states constitute a complete light nonet, including a $`K_0^{}(727)`$ or $`\kappa `$ meson, while the pattern of masses of the heavier scalars, which are nonetheless dominantly radial ground states, remains largely unaltered. For instance, the mostly strange $`f_0`$ stays around 1.5 GeV, while the mainly nonstrange $`f_0`$ settles close to 1.3 GeV, apparently quite compatible with the $`f_0(1370)`$ . On the other hand, Ref.  find neither a light $`\kappa `$, nor the $`f_0(1500)`$, and interpret the $`f_0(1370)`$ as mostly $`s\overline{s}`$. We do not wish to enter here into a detailed discussion on which interpretation is favored by experiment, and refer to Ref.  for an analysis of observed decay modes of the $`f_0(1370)`$ and $`f_0(1500)`$ that, we believe, support our view. In any case, this is evidently the most attractive scenario, in which there is a complete doubling of scalar states, but keeping the internal mass pattern of the nonets intact. To conclude the discussion of the $`f_0(1370)`$ , we would like to mention a very recent lattice calculation of scalar quarkonium masses , which is in agreement with our interpretation of the $`f_0(1370)`$ and $`f_0(1500)`$.

In conclusion and summarizing, by studying the $`W`$-emission process to describe the weak decay $`D_s^+f_0(980)\pi ^+`$ , we are able to reproduce the recently measured decay rate, provided that the $`f_0(980)`$ is assumed to be mostly $`s\overline{s}`$. A mixing angle of about $`14^{}`$, corresponding to a small $`n\overline{n}`$ admixture, is predicted. This interpretation of the $`f_0(980)`$ is in agreement with previous work of the authors , and also with other model studies . Finally, it is pointed out that a very recent experiment  suggests the $`f_0(1370)`$ is mostly $`n\overline{n}`$, as predicted by the unitarized meson model of Ref. .

Acknowledgement: We wish to thank Prof. R.L. Jaffe for useful comments on the first draft of our manuscript. This work is partly supported by the Fundação para a Ciência e a Tecnologia of the Ministério da Ciência e da Tecnologia of Portugal, under contract numbers PESO/P/PRO/15127/99, POCTI/35304/FIS/2000, and CERN/P/FIS/40119/2000.


text:
# One way to characterize the compact structures of lattice protein modelThis project was supported partly by Chinese Natural Science Foundation.

## 1 Introduction

The study of protein folding is fundamental on both theory and application. In order to tackle protein folding problem physically, it is important to pay much attention to concrete proteins and consider the details of interactions, such as for medical purpose. But there are also “global views” that should be noticed. For example, The possible configurations of folded proteins are enormous, while that can be observed in living form is rather limited. These protein structures generally can be described as belonging to a limit number of families. In each family, ignoring the details, the proteins possess similar overall conformations, and in many cases the structures show regular forms or approximate symmetry. Another example is that single domain proteins was observed only within a certain range of sequence length: the number of amino acid residues in single domain proteins seldom exceeds 200. Larger proteins usually fold into multi-domains native states.

With the accumulation of knowledge about the structures and functions of proteins, it was found that many proteins of similar structures pursue complete different functions, while proteins with different tertiary structures may perform similar functions. These suggested that to understand the protein folding problem physically, one should first get to know the properties of protein structures. Based on the the concepts from the physics of spin glass, study shows that to fold efficiently, proteins require a specially shaped energy landscape resembling a funnel. A heteropolymer with a completely random sequence generically possess a rugged energy landscape without a funnel. Goldstein et al have worked on optimizing energy functions for protein structure prediction. They found that some structures are more optimizable than others, i.e., there exist structures for which the funneled energy landscape can be obtained within a wide range of interaction parameters, while for some other structures the parameters for fast folding are much more restricted. The funneled landscape theory argued that the interactions in the folded structure must act in concert more effectively than expected in the most random cases. Accordingly, compared with most other structures, the superiority of highly optimizable structure should be that its geometric arrangement permit more sequences to reach the concert interaction states.

Other studies on the thermodynamic of lattice protein models also support the above idea. In the lattice HP models, a protein is represented by a self avoiding chain of beads placed on a discrete lattice with two types of beads: the Polar (P) and the Hydrophobic (H). A sequence is specified by a choice of monomer type at each position on the chain $`\{x_i\}`$. Where $`x_i`$ could be either H- or P-type, and $`i`$ is a monomer index. A structure is specified by a set of coordinates for all the monomers $`\{r_i\}`$. The energy has the form:

$$H=\underset{i<j}{}E_{x_ix_j}\mathrm{}(r_ir_j)$$

where $`\mathrm{}(r_ir_j)=1`$ when $`r_i`$ and $`r_j`$ are adjoining lattice sites while they are not adjacent along the sequence, and $`\mathrm{}(r_ir_j)=0`$ in other cases. Interaction parameter $`E_{x_ix_j}`$ differ according to the contact type HH, HP, or PP. Given the interaction parameters, it is possible to find out the ground state structure(s) of each sequence. Study shows that structures differ markedly in their tendency to be chosen by sequences as their unique ground states. The number of sequences which choice the structure as unique ground state is called the Designability of this structure. It was argued that only highly designable structures are thermodynamically stable and stable against mutation, and thus can be chosen by nature to fulfill the duty of life. Though interaction parameters used may differ strongly in different studies, the mostly designed structures do not depend strongly on the detail of interactions.

From above discussion we see that it should be essential to investigate the protein folding problem from structural point of view. To see the problem more clearly, we take square lattice HP model as an example. The total number of the most compact structures of 36 beads chain is 57337. Consider 36 beads homopolymer with interaction parameter $`E_{x_ix_j}=E_0<0`$. All the 57337 structures give the same energy when one such homopolymer fold onto each of them. Therefore the folded energy can not be used to distinguish the compact structures from each other. The essential here is that of discrimination, or characterization: give ways to tell how and why structures differ from each other. Nature’s way to break the symmetry is to replace homopolymer with heteropolymer. From this point of view, the success of lattice protein model is that it help to reveal this secret of nature.

Studies focusing on the properties of protein structures is still lack, in spite of some recent elaborations in this direction. In this article we present one way to break the symmetry, to distinguish the compact structures of lattice model without explicitly considering concrete interaction form. However, since only compact structures are considered here, an loose constraint is actually set on interactions: interactions under which compact structures are preferred as ground energy states. The method gives a number called partition number (Partnum) to each compact structure during a simple process. The Partnums of structures differ strongly, so giving one way to distinguish them from each other.

In the following section we will give the detail of the method, and compare the Partnum with designability. The statistical properties of Partnums are discussed in section II. The last section is for some remarks.

## 2 The definition and interpretation of Partnum

It is easy to find out all the compact structures of certain chain length with computer. Take 9 beads chain as an example. The search is self avoiding and restricted to the $`3\times 3`$ square lattice shown in Fig.1(A), and the resulting structures should not be related by rotation or reflection symmetry. As a result, there are only three starting points, $`(0,0)`$, $`(0,1)`$ and $`(1,1)`$, for the search of structures. To find the structures start at $`(0,0)`$, the first step is to go to $`(1,0)`$. This is the only choice, because $`(0,1)`$ is a symmetric point of $`(1,0)`$. We give all the structures following this step a number $`p_1=ln(1)`$. Now go to the next site. There are two possible choices: $`(2,0)`$ or $`(1,1)`$. Since the walk is self avoiding and restricted to the $`3\times 3`$ lattices, the walk following certain choice may fail to extend to 9 beads length. The choice that will reach to 9 beads length is called acceptable. Suppose that both $`(2,0)`$ and $`(1,1)`$ are acceptable. Then each compact structure which will be generated following $`(0,0)(1,0)(2,0)`$ or $`(0,0)(1,0)(1,1)`$ is given a number $`ln(1/2)`$. Generally speaking, restricting to $`3\times 3`$ lattice and beginning at a starting point, there are totally $`8`$ steps to finish a self avoiding walk. Each step is given one number according to the following rule: if the $`i`$-th step has totally $`C`$ acceptable choices not being symmetrically related, then the step is given a number called partnum of $`i`$-th step $`p_i=ln(1/C)`$. For 2D square lattice, the largest possible choice $`C_0`$ is 3.

Adding all the $`8`$ numbers and then dividing the sum by $`8`$, we get the Partnum $`P1`$. Here the structure is actually oriented . The consideration of oriented walk is reasonable in the case of protein structures, because the native protein structure would become unstable if the sequence is reversed, and also protein in life are produced successively from one end to another. However, if one consider the start and end reversal of the walk as a symmetric operation, then one oriented walk and its reverse together correspond to a structure that is not related with the direction. In the follows, oriented walk and non-oriented structures are used to distinguish the two different ways of viewing structure, and the Partnums corresponding to them are denoted as $`P1`$ and $`P2`$, respectively. However, when it is no need to distinguish them, simply structure is used and the Partnum is denoted as $`P`$. For the non-oriented structure, the Partnum can be define as: $`P2=P1(1)+P1(2)`$, where $`P1(1)`$ is the Partnum of one of two oriented walks and $`P1(2)`$ is that of its reverse.

The Partnums of structures of other chain length can be obtained similarly.

Since the original motivation of developing the Partnums of structures is to account for the difference of Designability of structures, in Fig.2 we give the plot of Designability against Partnum of orientd structures on $`5\times 5`$ lattice (the interaction parameters for calculating Designability is the same as used in Ref. ). There is not strict correspondence between Designability and Partnum. However the linear fit of the data revealed that Designability tends to increase with the increase of Partnum (see Fig.2). The same thing happens for other sequence length. In the case of $`6\times 6`$ lattice, the structure with highest Designability possess the second largest Partnum ($`P2`$).

According to Fig.1(B), an oriented walk corresponds to one path from the root to the top leave of the hierarchical tree. The value of Partnum of the structure is determined by the frequency of the path being disturbed by branches. If the path of a walk meet with fewer branches, the Partnum would be larger. This can be compared with the conclusion in Ref. . In Ref.  a simple version of HP model of protein is employed. A walk is reduced to a string of $`0`$s and $`1`$s, which represent the surface and core site respectively, as the backbone is traced. Each walk is therefore associated with a point in a high dimensional space. Sequences are represented by strings of their hydrophobicity and thus can be mapped into the same space. It was found that walks far away from other walks in the high dimensional space are highly designable and thermodynamically stable. For this reason, highly designable structures are called atypical in Ref. . Here the structures with large Partnum can also be called atypical (atypical average local freedom) since these structures correspond to paths on the hierarchical tree with fewer branches.

In an analog to the suggestion that nature selected out only highly designable structures, we assume that there exists a random process which selects out only the structures with the largest Partnum. It is interesting to see what this assumption will result in.

For concise we assume a critical Partnum $`P_c`$, so that only a small portion of oriented walks for which $`P1>P_c`$ can be selected out. Two oriented walks are called $`n`$-level similar if their first $`n1`$ steps are along the same path, and they branched at the $`n`$-th steps. Suppose $`s_1`$ is among the structures with the most highest Partnum satisfying $`P1(s_1)>P_c`$. This means that there are few branches along the path of $`s_1`$. As a result, it is difficult to find walks which show high level similarity to $`s_1`$. But if there do exist such walks, these walks should have high possibility to be selected out. For example, if $`s_2`$ is $`N1`$ level similar with $`s_1`$, $`N`$ being the chain length, then $`P1(s_2)=P1(s_1)>P_c`$. More generally, let $`n_{12}`$ being the similarity level between $`s_2`$ and $`s_1`$. We know that $`P1(s_2)=P1(s_1)(1\frac{n_{12}}{N1}ln(C_0))`$, $`C_0=3`$ being the maximal possible choices per step during the search of structures. According to this expression, the more similar $`s_2`$ is to $`s_1`$, the more possible it is to be selected out.

Assuming that $`s_3`$ is another walks with $`P1(s_3)>P_c`$, but it is dissimilar to $`s_1`$. From above discussion we know that there are two families, all the members of which are selected out. Within each family, the similarity level of two walks is much higher than $`n_{13}`$, while any two walks from different families are dissimilar from each other, and the similarity level is $`n_{13}`$. We thus come to the conclusion that the selected walks belong to separate families. walks within each family are similar, while walks belonging to different families are dissimilar.

For the non-oriented structures, there is no the convenience of the hierarchical tree to discuss their properties. But it is believable that the above result be kept once similarity between structures is properly defined. This is the case for the classification of real protein structures, where more or less arbitrary criteria are used to define the similarity between protein structures and to classify structure into families, superfamilies , folds, and so on.

## 3 The statistical properties of Partnums

Natural single domain proteins exist only within a limit range of sequence length. By both theoretical and numerical studies it is showed in Ref.  that the stability of folded sequences against mutation decrease with the increase of chain length. In that follows the dependence of the statistical properties of Partnum on chain length will be discussed. We will show how some structural properties are determined by general statistical principle.

The density distribution of $`P2`$ are shown in Fig.4. Things are similar for $`P1`$. In both cases, visually the distribution becomes more and more normal. Actually it will be shown that the distribution is Gauss distribution in the long chain limit. As the first step, however, let’s much generally, assume that the Partnums of chain length $`N`$ can be described by a density distribution function, $`F(P,v_1,v_2\mathrm{})`$, $`v_i`$ being the moment of $`i`$-th orders. It is easy to get the average $`v_1=<P>`$ and variance $`v_2=\mathrm{}P`$. The results of both oriented walk and non-oriented structures are shown in Fig.3.

Fig.3 shows that $`<P>`$ (both $`<P1>`$ and $`<P2>`$) decrease with the increase of chain length. However, from the definition of Partnum we know that $`<P1>`$ ($`<P2>`$) can not be smaller than $`ln3`$ ($`2ln3`$). So, for either oriented walks or non-oriented structures, there must exist $`\delta `$, so that

$$lim_N\mathrm{}<P>=\delta .$$

A similar argument applies to $`\mathrm{}P`$, where

$$lim_N\mathrm{}\mathrm{}P=ϵ,ϵ0.$$

It is known that the total number of compact structures $`M`$ increase exponentially with the increase of chain length $`N`$: $`M(N)(C_{av})^N`$, $`C_{av}<C_0=3`$ being the average number possible choices per step for the walks. This gives one way to estimate the value of $`C_{av}`$ using the knowledge of $`M(N)`$. Fig.5 show the fit of the data $`M(N)`$ to $`f(N)=(lnC_{av})N+b`$. The result is $`C_{av}=1.397`$. Viewing this value of $`C_{av}`$ as the value in long chain limit, we get that $`\delta =ln(1/C_{av})=0.3343`$ for oriented walks, a reasonable estimation (see Fig.3). It should be noticed that $`C_{av}`$ get this way is much larger that given by mean field consideration, where $`C_{av}=C_0/e=1.1`$. According to this $`C_{av}`$, $`\delta =0.099`$ for non-oriented walks. From Fig.3 we know that this is a value too large to be the long chain limit of $`<P1>`$. So it seems that the mean field treatment does not apply to the two dimensional protein model.

With the help of central limit theorem, we can argue that the density distribution is Gauss distribution in long chain limit, and $`ϵ=0`$. See follows.

In the space of compact structures, the Partnum $`P`$ of certain structure is the average of the partnums $`p_i`$ of all the $`T`$ steps. For oriented walks $`T`$ equals to the chain length subtracted by 1, and for non-oriented structures this value should be doubled further. Now divide the $`T`$ partnums into $`(T)/n`$ groups (suppose $`T/n`$ is an integer). In each group the $`n`$ members are chosen randomly within the total $`T`$ numbers. For each group we define a new random variable $`q_k=_ip_i/n`$, $`k`$ being the group index. Since the members in each group are chosen randomly form the total $`T`$ numbers, the $`T/n`$ newly defined random variable should have the same average and variance when $`n\mathrm{}`$. At the same time, since $`P=\frac{_kq_k}{(T1)/n}`$, applying the central limit theorem, we know that $`P`$ is a Gaussian random variable, and $`\delta P0`$ when $`T/n\mathrm{}`$.

From the above discussion we know that, according to Partnums, statistically all the compact structures become indistinguishable in long chain limit. Recalling the selection rule assumed above, we know that it becomes increasingly difficult to select out $`atypical`$ structures when chain length increases. These results show some connection to the work of Ejtihadi et al.. With a purely geometrical approach, they were able to reduce largely the candidates of structures that can be chosen as the ground states of sequences. They found that for the case of HP protein model the number of ground state candidates grows only as $`N^2`$, $`N`$ being the sequence length. While, as pointed out above, the total number of compact structures increase exponentially with the increase of $`N`$. So it becomes increasingly difficult to find the ground state candidates. This is in accordance with the statistical property of Partnum.

For fulfilling biology functions, proteins should possess some properties, for example fast folding, thermodynamically stable and stable against mutation. It was postulated that with the increase of sequence length, the folded structures become more and more difficult to possess these properties. Based on the study of Partnum, we propose that this property of proteins is determined by the statistical properties of protein structures, the detail of interaction having weak influence.

## 4 Conclusion Remarks

Protein structures seem to be a very special class among all the possible folded configurations of polypeptide chain. We now know something about how special it is, but little on why it be so. Ways of characterizing folded structures, from whatever point of view, will help to deepen our understanding about protein structures. In this paper, the study on Partnum itself is interesting, and more interesting when compared with the dynamic and thermodynamic study of proteins. The concept of Partnum is simple and can only be applied to lattice model. But the study on it reveals that it is possible to investigate protein structures with no consideration of interaction detail.

## ACKNOWLEDGMENTS

The authors would like to give thanks to Proff. Wei-Mou Zheng and Proff. Bai-Lin Hao for stimulating discussion. We also thank Mr. Guo-Yi Chen for many helps on computation.