text:
# Reversible stretching of homopolymers and random heteropolymers

## Abstract

We have analyzed the equilibrium response of chain molecules to stretching. For a homogeneous sequence of monomers, the induced transition from compact globule to extended coil below the $`\theta `$-temperature is predicted to be sharp. For random sequences, however, the transition may be smoothed by a prevalence of necklace-like structures, in which globular regions and coil regions coexist in a single chain. As we show in the context of a random copolymer, preferential solvation of one monomer type lends stability to such structures. The range of stretching forces over which necklaces are stable is sensitive to chain length as well as sequence statistics.

Experiments probing the mechanical response of individual protein molecules have demonstrated that certain domains can withstand significant stretching forces before unfolding . At strain rates much greater than unperturbed rates of unfolding, measured restoring forces imply that the native structures remain largely intact up to a threshold force, at which they unfold to an extensible state. Threshold forces for these domains are typically in the range $``$ 10-100 $`T/a`$, where $`T`$ is temperature (in units such that Boltzmann’s constant is unity) and $`a`$ is an average distance between neighboring, connected monomers in the chain. To be sure, these experimental systems are out of equilibrium, as highlighted by wide hysteresis upon releasing strain and broad distributions of threshold forces. Nevertheless, it is reasonable to expect that such proteins undergo sharp unfolding transitions under adiabatic stretching conditions as well.

The microscopic origin of this mechanical strength is not well understood. As Wolynes and coworkers have pointed out, small forces applied to the end-to-end distance of a protein may couple very weakly to the reaction coordinate for folding. Atomically detailed simulations of rapid stretching instead suggest an important role for certain backbone topologies that are stabilized by groups of hydrogen bonds . In this Letter we explore the equilibrium force-induced transition at a coarse-grained level, using simple estimates of the relevant free energetics. Specifically, we distinguish between aspects of the transition that are homopolymeric in nature, and those that arise from the quenched disorder characterizing random heteropolymers. We show that necklace-like structures, as depicted in Fig. 1 (c), occur with low probability in long homopolymers, but may be stabilized over a finite range of force and temperature by sequence heterogeneity. We describe the features of sequence statistics that affect this stability, and thereby determine mechanical strength. Our results may thus shed light on evolutionary design principles for proteins whose functions are mechanical in nature.

We consider first a chain of $`N`$ identical monomers in solution. In the absence of stretching force and below the $`\theta `$-temperature ($`\theta `$), the chain adopts a compact, globular conformation (Fig 1 (a)). Because this state exhibits only small fluctuations in monomer density, its free energy is well approximated using mean field theory:

$$F_\mathrm{g}BN+\gamma N^{2/3}TS_0.$$
(1)
Here, $`B`$ is the excess free energy per particle of a fluid of unconnected monomers at the same temperature and density as the globule. The second term in Eq. 1 includes surface energetics of the globule-solvent interface, as well as the conformational entropy $`S_{\mathrm{sph}}`$ lost upon confining the chain to a spherical volume with radius $`RN^{1/3}`$. $`S_0`$ is the entropy of an ideal, unconfined chain.

We imagine that the principal effect of a small stretching force, $`f`$, on the ends of a homopolymeric globule is to deform its spherical geometry (Fig. 1 (b)). The favorable energy of this deformation, $`f(R_{}R)`$, is offset by surface energetics as well as a loss of entropy, $`S_\mathrm{d}`$. We estimate this entropy loss through the statistics of a Gaussian chain confined to a deformed volume. In the long-chain limit, the free energy per monomer of such a chain is isomorphic with the ground state energy of a quantal particle confined to the same volume. Treating the asymmetric boundary condition as a perturbation, we sum an infinite class of terms in the ground-state expansion, obtaining

$$S_\mathrm{d}=S_{\mathrm{sph}}\left(\frac{2R_{}}{3R}+\frac{R^2}{3R_{}^2}1\right).$$
(2)
Because $`S_{\mathrm{sph}}N^{2/3}`$, the energy gained by a reasonable deformation of the globule ($`R_{}N^{1/3}`$), is comparable to $`S_\mathrm{d}`$ only for forces of magnitude $`N^{1/3}`$ or larger. For long chains, the deformation achieved by stretching forces $`T/a`$ is therefore negligible. We subsequently consider the globule to be undeformed, and the globular free energy to be unaffected by stretching.

For sufficiently large stretching forces, the dominant state of a homopolymer is an expanded coil (Fig. 1 (d)). In contrast to the globule, this state is characterized by extensive density fluctuations and is strongly susceptible to deformation. Considering short-ranged attractions between monomers to be unimportant in this case, we model the extended coil as a freely jointed chain. The free energy of such a chain, $`F_\mathrm{c}(N)=NT\mathrm{ln}[\mathrm{sinh}(fa/T)/(fa/T)]TS_0`$, is sensitive to the magnitude of $`f`$. At a given temperature $`T<\theta `$, a phase transition occurs at a force sufficient to lower the coil free energy below that of the globule. (See inset of Fig. 2.) For $`N^{2/3}1`$, the phase boundary is given by

$$B+\gamma N^{1/3}=T\mathrm{ln}[\mathrm{sinh}(fa/T)/(fa/T)].$$
(3)
The resulting phase diagram in the force-temperature plane is shown in Fig. 2 for various $`N`$. Here, we have taken $`B=T\theta `$ (accurate near $`T=\theta `$), which is unrealistic at low temperatures. It is reasonable, however, that as temperature decreases, stretching energy density of the strongly fluctuating coil grows more quickly than attractive energy density of the relatively placid globule. As a result, a reentrant coil phase appears at low temperature. The “coil” in this case is a nearly straight chain with small fluctuations in extension. Computer simulations of strained lattice heteropolymers appear to support our prediction of reentrance.

Although the globule-coil transition is second-order in the absence of stretching, the force-induced transition is here predicted to be first-order, since the average chain extension is discontinuous at the phase boundary. This result is certainly correct for $`T<\theta `$, where globule and coil phases are distinct. Near $`T=\theta `$ and $`f=0`$, however, our caricatures of these states are oversimplified: Density fluctuations are not negligible in the globule, and attractions are not negligible in the coil. In this region of the phase diagram, extension increases smoothly with force, and the stretching transition is second-order.

In constructing a phase diagram for homopolymer stretching, we did not consider chain structures that are necklace-like (i.e., coexisting globule and coil regions, as in Fig. 1 (c)). Neglecting surface effects, the free energy per monomer of a necklace lies between that of a globule and that of a coil. The entropy of a necklace is augmented by the freedom of globular regions to reside anywhere along the chain (provided they do not overlap). For a single globular region, however, this additional “translational” entropy scales as $`\mathrm{ln}N`$, and is insufficient to overcome the $`O(N)`$ deficit in free energy to either globule or coil, even for chains of modest length. In the case of many globular regions, the gained entropy, $`N\mathrm{ln}2`$, is considerable, but does not compensate for the cost of presenting an extensive surface to the solvent.

The stability and importance of necklace structures may be qualitatively different for heteropolymers. In particular, fluctuations in local sequence composition yield a preference for globule or coil that varies along the chain. We assess the strength of this effect and its influence on stretching behavior, for a two-letter random copolymer. In this model, each monomer has two possible identities, $`\sigma _i=\pm 1`$, perhaps denoting hydrophobic and hydrophilic moieties. For a given sequence $`\{\sigma _i\}`$, the energy of a chain conformation defined by monomer positions $`𝐫_i`$ is

$``$ $`=`$ $`_0+\mathrm{\Gamma }{\displaystyle \underset{i\mathrm{exposed}}{}}\sigma _i𝐟(𝐫_N𝐫_1),`$ (4)
$`_0`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{N}{}}}\delta (𝐫_i𝐫_j)(B_0+\chi \sigma _i\sigma _j),`$ (5)
where $`B_0`$ is a mean attractive energy density stabilizing the globule at low temperatures. We consider $`\chi <0`$, so that attractions are strongest between monomers of the same type. The second sum in Eq. 4 includes only monomers that are exposed to solvent. For $`\mathrm{\Gamma }>0`$, monomers of type $`\sigma _i=1`$ are favorably solvated, while exposure of $`\sigma _i=1`$ monomers is unfavorable. In Ref. , the phase diagram was determined for the model defined by $`_0`$ at the mean field level, using the replica trick to average over random sequences of monomer types. The one-step replica symmetry breaking demonstrated in that work is consistent with a suitably chosen random energy model. In other words, the average thermodynamics are reproduced by drawing $`(a^3/v)^N`$ energy levels at random from a distribution $`P(E)=(\pi N\mathrm{\Delta }^2)^{1/2}\mathrm{exp}[(E\overline{E})^2/N\mathrm{\Delta }^2],`$ with variance $`\mathrm{\Delta }=|\chi |\mu ^2\overline{\rho }`$ and mean $`\overline{E}=B_0N\overline{\rho }`$. Here, $`v`$ is the volume occupied by a monomer, $`\overline{\rho }v^3`$ is the mean density in the globule core, and $`\mu `$ is the variance of the monomer distribution. Here we describe the effect of solvation and stretching energetics on the effective distribution of energies.

For a compact, spherical globule with surface area $`A`$, the solvation contribution in Eq. 4 broadens the energy distribution, increasing $`\mathrm{\Delta }`$ by $`(\mathrm{\Gamma }^2\overline{\rho }/4|\chi |)A/N`$. This result is obtained following the analysis of Ref. , with the additional assumption that the spatial pattern of solvent-exposed monomers is independent of compact chain conformation. As in the case of the homopolymer, we neglect the energetic contribution of globule deformation. The resulting distribution of energies is dominated by states in the interval $`\overline{E}N^{1/2}\mathrm{\Delta }<E<\overline{E}+N^{1/2}\mathrm{\Delta }`$. At energies just below a critical value, $`E^{}=\overline{E}N\mathrm{\Delta }(\mathrm{ln}a^3/v)^{1/2}`$, the number of states is $`O(1)`$, while just above $`E^{}`$ the number is exponentially large. The ground states of particular random sequences are distributed narrowly about $`E^{}`$. Introducing solvation energetics thus lowers the average ground state energy by an amount $`(\mathrm{\Gamma }^2\overline{\rho }/4|\chi |)(\mathrm{ln}a^3/v)^{1/2}A`$. If $`\chi `$ and $`\mathrm{\Gamma }`$ are comparable in magnitude, this shift is a significant fraction of the energy gained by exposing only monomers with $`\sigma _i=1`$ to the solvent. The density of states around $`E^{}`$ in this case is sufficiently large that a ground state with favorable solvation energetics may always be selected.

In an extended coil, nearly all monomers are exposed to solvent. For sequences with fixed total composition, $`_{i=1}^N\sigma _i=0`$, thermodynamics of the coil state are unaffected by solvation energetics. We consider the heteropolymer coil to have essentially the same physics as the homopolymer coil described above.

For necklace structures of a random heteropolymer, free energy depends on the positions of globular regions along the chain. In effect, these globules move in a random potential generated by sequence disorder. (See Fig. 3.) The scale and correlation length of this potential are determined by the size of the globule, and by statistics of the sequence. Although the total composition of the chain is fixed, a globular region with $`n<N`$ monomers, situated at a given point on the chain, has in general an excess of one monomer type:

$$q\frac{1}{n}\underset{i\mathrm{globule}}{}\sigma _i0.$$
(6)
Because the local composition $`q`$ is a sum of many independent random variables, its distribution is Gaussian, with variance $`\mu n^{1/2}`$. This distribution of sequence compositions leads to localization of globules along the chain.

At a particular location, a globular region has an apparent distribution of monomer types that is modified according to the local value of $`q`$, $`p(\sigma _i;q)\mathrm{exp}[(\sigma _iq)^2/2\mu ^2(1q^2)]`$. The ground state energy of the globule at this location, determined from the random energy model described above, is

$$E^{}=(B_0|\chi |q^2)\overline{\rho }n\mathrm{\Gamma }q\overline{\rho }n^{2/3}\left(\mathrm{ln}\frac{a^3}{v}\right)^{1/2}\left[2|\chi |\mu ^2(1q^2)\overline{\rho }n+\frac{\mathrm{\Gamma }^2\overline{\rho }}{4|\chi |}n^{2/3}\right].$$
(7)
Since local composition is not fixed for a heteropolymer, $`E^{}`$ is a random variable for different globule locations. But its fluctuations are weak: $`(\delta E^{})^2^{1/2}=\mathrm{\Gamma }\overline{\rho }\mu n^{1/6}+O(n^0)`$, where angled brackets denote an average over the distribution of $`q`$.

For a single globular region incorporating $`n`$ monomers, the remaining $`Nn`$ monomers of the necklace belong to coil regions. When the globule resides at a given chain location with composition $`q`$, the complementary composition in coil regions is $`q`$. Consequently, the solvation energy of these fully exposed regions, $`\mathrm{\Gamma }nq`$, is also a random variable, with fluctuations of magnitude $`\mathrm{\Gamma }\mu n^{1/2}`$. These relatively strong fluctuations in coil solvation energy establish the scale of the random potential for globule motion along the chain. Below a critical temperature $`T_\mathrm{c}`$, the globule will become localized in the deepest minimum of this potential. Following Derrida’s analysis of randomly distributed energies, we find that $`T_\mathrm{c}=\mathrm{\Gamma }\mu n^{1/2}/\sqrt{2s}`$, and that the free energy due to globule motion is

$$F_{\mathrm{rand}}=\{\begin{array}{cc}sT[1+(\frac{T_\mathrm{c}}{T})^2],\hfill & T>T_\mathrm{c}\hfill \\ 2sT_\mathrm{c},\hfill & TT_\mathrm{c}.\hfill \end{array}$$
(8)
In Eq. 8, $`s`$ is the logarithm of the number of possible globule locations. The number of statistically independent locations is approximately a factor of $`n`$ smaller, so that $`s`$ is appropriately renormalized, $`s\mathrm{ln}(N/n)1`$. For $`TT_\mathrm{c}`$, randomness of the potential is irrelevant, and the homopolymer result, $`T\mathrm{ln}N`$ is recovered.

The analysis for a single globule is readily generalized for a necklace with several globular regions. In this case, $`e^sn^M(Nn+1)(N2n+1)\mathrm{}(NMn+1)/M!`$, where $`M`$ is the number of globules. For small globules, $`n=O(1)`$, the minimum value of $`F_{\mathrm{rand}}`$ is obtained for $`M=O(N)`$, giving $`\mathrm{min}F_{\mathrm{rand}}=O(N^{1/2})`$. For large globules, $`n=O(N)`$, $`M`$ is necessarily $`O(1)`$, and again $`\mathrm{min}F_{\mathrm{rand}}=O(N^{1/2})`$. Globule sizes are thus expected to be distributed broadly, with a preference for a small number of large globules dictated by surface tension. Due to the $`O(N^{1/2})`$ stabilization provided by globule localization, a phase rich in necklace structures covers an appreciable range of force and temperature for $`N=10^2`$ (Fig. 4 (a)). Indeed, simulated stretching of a short, nearly random heteropolymer involves partially extended structures. For $`N=10^3`$ (Fig. 4 (b)), the necklace phase covers a much smaller region of the phase diagram, due to the extensive free energy cost of mixing globule and coil structures. Although necklaces become unstable as $`N\mathrm{}`$, our results suggest that they may be prevalent for macromolecules of finite size.

From our results for uncorrelated sequence statistics, we may deduce the basic effects of introducing correlations. For “blocky” sequences, in which monomers within a correlation length $`\xi `$ are likely to be of the same type, fluctuations in local composition are large. Specifically, $`(\delta q)^2=O(n^0)`$ for $`n<\xi `$. If $`\xi `$ scales with chain length (or if sequence correlations decay algebraically with distance along the chain), the resulting free energy contribution due to localization is $`O(N)`$, and necklace structures persist as $`N\mathrm{}`$. By contrast, sequences that are anticorrelated on a scale $`\xi `$ exhibit small fluctuations in local composition. If $`(\delta q)^2=O(n^{1/2})`$, fluctuations in solvation energy are $`O(1)`$, and the stretching behavior will be homopolymeric.

The equilibrium behavior of strained polymers described in this Letter is not directly related to the intrinsically dynamical, nonequilibrium response measured in experiments. Nonetheless, it is reasonable to associate the equilibrium stability of various structures with their kinetic accessibility. In particular, when necklaces are stable over a wide range of stretching forces, the kinetic transition from compact globule to extended coil is likely not sharp. Instead, we expect that the chain visits an ensemble of necklace structures as it passes through the transition region. The breadth we have predicted for this ensemble in the case of random heteropolymers suggests that the restoring force to applied strain should exhibit large fluctuations as the chain is stretched. Such a scenario has in fact been observed for the protein barnase. The observation that certain protein domains do undergo sharp stretching transitions is thus an indication of evolutionary design for mechanical strength, as is reflected by the roles of such proteins in cell adhesion and muscle elasticity.


format: text

text:
# Comparison of the fits to data on polarised structure functions and spin asymmetries

## Abstract

In order to obtain polarised parton densities we have made next to leading order QCD fit using experimental data on deep inelastic structure functions on nucleons. This fit is compared with the updated fit to corresponding spin asymmetries. We get very similar results for all fits also for different data samples. It seems that only polarised parton densities for non-strange quarks $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ are relatively well determined from the present polarised deep inelastic experiments. Integrated gluon contribution at $`Q^2=1\mathrm{GeV}^2`$ is, as in our previous fits, very small.

In order to get polarised parton (i.e. quark and gluon) densities one uses data on deep inelastic scattering of polarised lepton on polarised nucleon targets. Quite a lot of data exist for such scattering on different nucleon targets. The data come from experiments made at SLAC \[1-10\], CERN \[11-16\] and DESY . Recently the data from SLAC from E155  experiment on protons has been published. The experimental groups present data for spin asymmetries as well as for polarised structure functions.

The analysis of the EMC group results  started an interest in studying such data. So called ”spin crisis” was connected with the fact that only very little of spin of nucleon was carried by quarks. The suggestion from ref.  was that polarised gluons may be responsible for this phenomena. Many next to leading (NLO) order QCD analysis \[20-27\] were performed and polarised parton distributions were determined. The main purpose of this paper is to use data for polarised structure functions $`g_1(x,Q^2)`$ on proton, neutron and deuteron targets in order to determine polarised parton distributions. This fit will be compared with our updated (in which we take into account recently published data on protons from E155  experiment in SLAC) fits to spin asymmetries. It was advocated by us  and  that using spin asymmetries for determination of polarised parton densities one avoids the problem with higher twist contributions. On the other hand it is the polarised structure fuctions and polarised parton distributions that we want to determine.

We compare these two ways of making fits with similar technical assumptions and the same parton functions used for fitting. We will see that both methods give very similar results for parton distributions. Our previous  important conclusion that the integrated gluon contribution is negligible at $`Q^2=1\mathrm{GeV}^2`$ does not change. Most of the groups used experimental data for spin asymmetries to determine polarised parton distributions. In addition to spin asymmetries one has also experimental data for polarised structure functions calculated from the spin asymmetries in a specific way chosen by an experimental group. There were also several fits to the data on polarised structure functions . As in  we will make fits to two samples of the data. In the first group we will have data for the same $`x`$ (strictly speaking for the near values) and different $`Q^2`$ and in the second the ”averaged” data where one averages over $`Q^2`$ (the errors are smaller and $`Q^2`$ dependence is smeared out). In most of the fits to experimental data only second sample (namely with averaged $`Q^2`$ dependence) was used. Our fits use the both sets of data. The data for polarised structure functions are usually given only for averaged sample of data (the exception is E155 experiment for deuteron and proton data).

Experiments on unpolarised targets provide information on the unpolarised quark densities $`q(x,Q^2)`$ and $`G(x,Q^2)`$ inside the nucleon. These densities can be expressed in term of $`q^\pm (x,Q^2)`$ and $`G^\pm (x,Q^2)`$, i.e. densities of quarks and gluons with helicity along or opposite to the helicity of the parent nucleon.

$$q=q^++q^{},G=G^++G^{}.$$
(1)
$`q`$ stands for quark and antiquark contributions.

The polarised parton densities, i.e. the differences of $`q^+`$, $`q^{}`$ and $`G^+`$, $`G^{}`$ are given by:

$$\mathrm{\Delta }q=q^+q^{},\mathrm{\Delta }G=G^+G^{}.$$
(2)
We will try to determine $`q^\pm (x,Q^2)`$ and $`G^\pm (x,Q^2)`$, in other words, we will try to connect unpolarised and polarised data.

In our fits we will use functions for the polarised parton densities that are suggested by the fit to unpolarised data. We risk that asymptotic behaviour of our parton distributions is not correct but it seems to us not so important when we limit ourselves to the measured region of x. Maybe it is also not bad to use parametrisation completely different from other groups to check how the results depend on that. Our starting point are the formulas for unpolarised quark and gluon distributions gotten (at $`Q^2=1\mathrm{GeV}^2`$) from the fit performed by Martin, Roberts, Stirling and Thorne  (they use $`\mathrm{\Lambda }_{\overline{MS}}^{n_f=4}=0.3`$ GeV and $`\alpha _s(M_Z^2)=0.120`$)

We will not use small and high x behaviour of unpolarised parton distributions as fitted parameters as some other groups do. We will split $`q`$ and $`G`$, as was already discussed in ref. , into two parts in such a manner that the distributions $`q^\pm (x,Q^2)`$ and $`G^\pm (x,Q^2)`$ remain positive. Our polarised densities for quarks and gluons are parametrised as follows:

$`\mathrm{\Delta }u_v(x)=x^{0.5911}(1x)^{3.395}(a_1+a_2\sqrt{x}+a_4x),`$
$`\mathrm{\Delta }d_v(x)=x^{0.7118}(1x)^{3.874}(b_1+b_2\sqrt{x}+b_3x),`$
$`2\mathrm{\Delta }\overline{u}(x)=0.4\mathrm{\Delta }M(x)\mathrm{\Delta }\delta (x),`$ (3)
$`2\mathrm{\Delta }\overline{d}(x)=0.4\mathrm{\Delta }M(x)+\mathrm{\Delta }\delta (x),`$
$`2\mathrm{\Delta }\overline{s}(x)=0.2\mathrm{\Delta }M_s(x),`$
$`\mathrm{\Delta }G(x)=x^{0.0829}(1x)^{6.587}(d_1+d_2\sqrt{x}+d_3x).`$
whereas for the antiquarks ( and sea quarks ):

$`\mathrm{\Delta }M(x)=x^{0.7712}(1x)^{7.808}(c_1+c_2\sqrt{x}),`$
$`\mathrm{\Delta }M_s=x^{0.7712}(1x)^{7.808}(c_{1s}+c_{2s}\sqrt{x}),`$ (4)
$`\mathrm{\Delta }\delta (x)=x^{0.183}(1x)^{9.808}c_3(1+9.987x33.34x^2),`$
We also have

$`\mathrm{\Delta }u=\mathrm{\Delta }u_v+2\mathrm{\Delta }\overline{u},`$
$`\mathrm{\Delta }d=\mathrm{\Delta }d_v+2\mathrm{\Delta }\overline{d},`$ (5)
$`\mathrm{\Delta }s=2\mathrm{\Delta }\overline{s},`$
We use as before  additional independent parameters for the strange sea contribution with the same as for non-strange sea functional dependence. Maybe not all parameters are important in the fit and it could happen that some of the coefficients in eq.(3,4) taken as free parameters in the fit are small or in some sense superfluous. Putting them to zero (or eliminating them) increase $`\chi ^2`$ only a little but makes the quantity $`\chi ^2/N_{DF}`$ smaller. We will see that that is the case with some parameters introduced in eq. (3,4).

In order to get the unknown parameters in the expressions for polarised quark and gluon distributions (eqs.(3,4)) we calculate the spin asymmetries (starting from initial $`Q^2`$ = 1 $`\text{GeV}^2`$) for measured values of $`Q^2`$ and make a fit to the experimental data on spin asymmetries for proton, neutron and deuteron targets. The spin asymmetry $`A_1(x,Q^2)`$ can be expressed via the polarised structure function $`g_1(x,Q^2)`$ as

$$A_1(x,Q^2)\frac{(1+\gamma ^2)g_1(x,Q^2)}{F_1(x,Q^2)}=\frac{g_1(x,Q^2)}{F_2(x,Q^2)}[2x(1+R(x,Q^2))],$$
(6)
where $`R=[F_2(1+\gamma ^2)2xF_1]/2xF_1`$ whereas $`F_1`$ and $`F_2`$ are the unpolarised structure functions and $`\gamma =2Mx/Q`$ ($`M`$ stands for proton mass). We will take the new determined value of $`R`$ from the . The factor $`(1+\gamma ^2)`$ plays non negligible role for $`x`$ and $`Q^2`$ values measured in SLAC experiments. In calculating $`g_1(x,Q^2)`$ and $`F_2(x,Q^2)`$ in the next to leading order we use procedure described in  following the method described in  performing calculations with Mellin transforms and then calculating Mellin inverse. Having calculated the asymmetries according to equation (6) for the value of $`Q^2`$ obtained in experiments we can make a fit to asymmetries on proton, neutron and deuteron targets. The other possibility is to use directly the data for the polarised structure function $`g_1(x,Q^2)`$ (the values of $`g_1`$ were given for the averaged values of $`Q^2`$) on proton, neutron and deuteron targets to determine unknown coefficients in expressions for polarised parton distributions. The problem is that different experimental groups used their own specific methods to obtain the values of polarised structure functions $`g_1(x,Q^2)`$ from the measured asymmetries. At the end we want to know polarised structure functions and polarised parton distributions. In order to calculate them from spin asymmetries we have to choose what shall we take for $`F_1(x,Q^2)`$ or $`F_2(x,Q^2)`$ and $`R(x,Q^2)`$. As was already mentioned we calculated $`F_2(x,Q^2)`$ in NLO for actual values of $`x`$ and $`Q^2`$ using quark and gluon contributions for $`Q^2=1\mathrm{GeV}^2`$ given by MRST . The values of R were in the earlier fits taken from Whitlow  and later from E143 group .We have treated all experiments in the same way. There is also a problem of higher twist corrections (power low corrections to R were included). We will not include the higher twist corrections because of still big experimental errors. The spread of the results could be a measure of uncertainties in both methods. We will compare fits using determination of parameters from polarised structure functions and from spin asymmetries. As will be seen later the results obtained by two methods are very similar.

It was already seen from our previous paper  that $`g_1(x,Q^2)`$ calculated from spin asymmetries fits not bad the data points for $`g_1`$. Data points for polarised structure functions were given for averaged data set so it is natural to compare fit to $`g_1`$ (called by us fit $`g`$) with the fit to spin asymmetries (fit $`A_1`$) with the same number of points. We will take into account in this case 197 data points (We will take E155 proton and deuteron data without averaging). These fits will be compared with the fit (fit $`A_2`$) to spin asymmetries for non averaged data where we take into account 465 experimental data points. At the beginning we will not put any constrains from hyperon decays. Later we will also not fix $`a_8=\mathrm{\Delta }u+\mathrm{\Delta }d2\mathrm{\Delta }s`$ value but we will add experimental point $`a_8=0.58\pm 0.1`$ with enhanced (to 3$`\sigma `$) error. That means we will simply add to $`\chi ^2`$ corresponding to experimental points for spin asymmetries the term connected with experimental point from hyperon decays. We will discuss how this additional experimental point influences our results.

Not all parameters are important in the fits. It seems that some of the parameters of the most singular terms are superfluous and we can eliminate them. We will put $`d_1=0`$ (such assumption gives that $`\delta G/Gx^{1/2}`$ for small $`x`$), $`b_1=0`$ (the most singular term in $`\mathrm{\Delta }d_v`$) and assume $`c_{1s}=c_1`$ (i.e. the most singular terms for strange and non-strange sea contributions are equal). Fixing these four parameters in the fit practically does not change the value of $`\chi ^2`$ but improves $`\chi ^2/N_{DF}`$. We also have to make some remarks about parameters $`c_3`$ and $`d_2`$. In the fit for $`g_1`$ these parameters are important i.e. when we eliminate them $`\chi ^2`$ per degree of freedom increases. That is not the case in fits for spin asymmetries. In this case these parameters are superficial (that means for example that splitting $`\overline{u}\overline{d}`$ is not well defined by data to spin asymmetries. We will leave this parameters for comparison with fit to $`g_1`$ (do not eliminate them) but they could be not well determined and cause some artificial shifts in other parameters.

Table 1. The parameters of our three fits calculated at $`Q^2=1\mathrm{GeV}^2`$ together with $`\chi ^2`$ per degree of freedom

| fit | $`a_1`$ | $`a_2`$ | $`a_4`$ | $`b_2`$ | $`b_3`$ | $`c_1`$ |
| --- | --- | --- | --- | --- | --- | --- |
| fit $`g`$ | $`0.61`$ | $`7.05`$ | $`17.1`$ | $`2.02`$ | $`0.34`$ | $`0.34`$ |
|  | $`\pm 0.0`$ | $`\pm 0.23`$ | $`\pm 0.23`$ | $`\pm 0.0`$ | $`\pm 0.24`$ | $`\pm 0.03`$ |
| fit $`A_1`$ | $`0.56`$ | $`5.50`$ | $`14.7`$ | $`1.67`$ | $`0.173`$ | $`0.338`$ |
|  | $`\pm 0.13`$ | $`\pm 1.22`$ | $`\pm 1.66`$ | $`\pm 0.03`$ | $`\pm 0.26`$ | $`\pm 0.10`$ |
| fit $`A_2`$ | $`0.49`$ | $`5.34`$ | $`14.66`$ | $`2.02`$ | $`0.35`$ | $`0.32`$ |
|  | $`\pm 0.13`$ | $`\pm 1.20`$ | $`\pm 1.63`$ | $`\pm 0.0`$ | $`\pm 0.25`$ | $`\pm 0.10`$ |

| fit | $`c_2`$ | $`c_{2s}`$ | $`c_3`$ | $`d_2`$ | $`d_3`$ | $`\chi ^2/N_{DF}`$ |
| --- | --- | --- | --- | --- | --- | --- |
| fit $`g`$ | $`4.15`$ | $`4.15`$ | $`1.05`$ | $`29.0`$ | $`87.1`$ | $`0.87`$ |
|  | $`\pm 0.0`$ | $`\pm 0.0`$ | $`\pm 0.56`$ | $`\pm 8.6`$ | $`\pm 36.1`$ |  |
| fit $`A_1`$ | $`3.23`$ | $`4.15`$ | $`0.617`$ | $`15.4`$ | $`42.2`$ | $`0.81`$ |
|  | $`\pm 0.80`$ | $`\pm 0.22`$ | $`\pm 0.58`$ | $`\pm 0.22`$ | $`\pm 15.0`$ |  |
| fit $`A_2`$ | $`3.69`$ | $`4.15`$ | $`0.39`$ | $`14.0`$ | $`27.0`$ | $`0.84`$ |
|  | $`\pm 0.79`$ | $`\pm 0.20`$ | $`\pm 0.48`$ | $`\pm 0.04`$ | $`\pm 11.4`$ |  |

In the Table 1 we present the values of parameter from the fit to the data on polarised structure functions and spin asymmetries for averaged and non averaged data together with $`\chi ^2/N_{DF}`$ values.

For example in the case of first row in table 1 corresponding to the fit to polarised structure functions the obtained quark and gluon distributions lead for ($`Q^2`$ =1 $`GeV^2`$) to the following integrated (over $`x`$) quantities: $`\mathrm{\Delta }u=0.72,\mathrm{\Delta }d=0.64,\mathrm{\Delta }s=0.05,\mathrm{\Delta }u_v=0.54,\mathrm{\Delta }d_v=0.65,2\mathrm{\Delta }\overline{u}=0.18,2\mathrm{\Delta }\overline{d}=0.01.`$

We have positively polarised sea for up and down quarks and positively polarised sea for strange quarks. The gluon polarisation is small. The value of $`a_3=1.36`$ was not assumed as an input in the fit (as is the case in nearly all fits ) and comes out slightly higher than the experimental value. The value of $`a_8=0.01`$ is completely different from the experimental figure. Taking into account that fits to polarised structure functions and spin asymmetries use different methods to calculate $`F_2(x,Q^2)`$ and $`R(x,Q^2)`$ there is no reason to expect that they give exactly the same results. The obtained values of parameters are very close and practically within experimental errors. The parameters calculated in fits to spin asymmetries (line 2 and 3) are closer in comparison with fit to $`g_1`$ (line 1) but are not identical. The spread of parameters measures small differences in the $`Q^2`$ evolution, differences in experimental errors and influence of our specific functions used in fits.

As was already mentioned in  the asymptotic behaviour at small $`x`$ of our polarised quark distributions is determined by the unpolarised ones and hence do not have the expected theoretically Regge type behaviour. Some of the quantities specially integrated sea contributions and also some valence contributions in our fit change rapidly for $`x0.003`$. that is not something that we expect from Regge behaviour with small exponent.

Hence, we will present quantities integrated over the region from $`x`$=0.003 to $`x`$=1 (it is practically integration over the region which is covered by the experimental data, except of non controversial extrapolation for highest $`x`$). The values of integrated quantities in the measured region we consider as more reliable then those in the whole region.

The corresponding quantities for three our fits together with $`\chi ^2/N_{DF}`$ are presented in table 1a.

Table 1a. The values of quark and gluon polarisations at $`Q^2=1\mathrm{GeV}^2`$ for our three fits

| fit | $`\mathrm{\Delta }u`$ | $`\mathrm{\Delta }d`$ | $`\mathrm{\Delta }u_v`$ | $`\mathrm{\Delta }d_v`$ | $`2\mathrm{\Delta }\overline{u}`$ | $`2\mathrm{\Delta }\overline{d}`$ | $`2\mathrm{\Delta }\overline{s}`$ | $`\mathrm{\Delta }G`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| fit $`g`$ | $`0.74`$ | $`0.49`$ | $`0.44`$ | $`0.63`$ | $`0.30`$ | $`0.14`$ | $`0.11`$ | $`0.15`$ |
| fit $`A_1`$ | $`0.75`$ | $`0.48`$ | $`0.57`$ | $`0.57`$ | $`0.18`$ | $`0.09`$ | $`0.11`$ | $`0.01`$ |
| fit $`A_2`$ | $`0.76`$ | $`0.47`$ | $`0.54`$ | $`0.63`$ | $`0.22`$ | $`0.16`$ | $`0.12`$ | $`0.19`$ |

From the table we see that there are changes in valence and sea contributions in different fits but the values of $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ practically do not differ. We use the parametrisation where the most singular term in sea contribution is very similar to valence quark terms and that maybe is the reason why splitting into valence and sea contribution is fragile and changes for different fits but in the case of $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ do not differ much. In the first fit at $`Q^2=1\mathrm{GeV}^2`$ we get $`\mathrm{\Gamma }_1^p=0.124`$, $`\mathrm{\Gamma }_1^n=0.052`$, $`a_3=1.23`$ and $`\mathrm{\Delta }\mathrm{\Sigma }=0.36`$ comparing with the second fit to averaged spin asymmetries where we have $`\mathrm{\Gamma }_1^p=0.125`$, $`\mathrm{\Gamma }_1^n=0.051`$, $`a_3=1.24`$ and $`\mathrm{\Delta }\mathrm{\Sigma }=0.38`$. These results are very close. The value of $`a_3`$ in the measured region without any assumption comes out close to the value measured in hyperon decays.

For illustration we present in Fig.1 the distributions $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ for our sets of parameters calculated from polarised structure functions and the sample of averaged spin asymmetries data. The corresponding values for $`\mathrm{\Delta }u_v`$ and $`\mathrm{\Delta }d_v`$ differ much stronger. In our previous paper  we already presented how the values of polarised structure functions for proton, deuteron and neutron calculated from the fits to spin asymmetries compare with the experimental data. To see what is the difference in the values fitted directly to the polarised structure functions and the values calculated from the fit to spin asymmetries we present in Fig.2 the corresponding curves in comparison with experimental points for $`g_1^p`$, $`g_1^d`$ and $`g_1^n`$ at the values of $`Q^2`$ in corresponding experiment. As an example we show comparison with experimental points for polarised structure functions $`g_1^p`$, $`g_1^d`$ from SMC from CERN and E143 from SLAC and $`g_1^n`$ from E142 and E154 from SLAC. There are some differences but they are not big in comparison with experimental errors. The comparisons for other experimental sets look very similar. It is of course not astonishing that the fitted curves with the parameters from the first fit are closer to experimental values.

As we already pointed out before we have not made any assumptions about $`a_8`$. We got from the fits that the value of $`a_8`$ is near zero very far from the experimental value and we got positive values for $`\mathrm{\Delta }s`$. The value of $`a_3`$ that also was not constrained in the fit is close to experimental value (in the measured region of $`x`$). In order to make more direct comparison with other fits as before we will also not fix $`a_8=\mathrm{\Delta }u+\mathrm{\Delta }d2\mathrm{\Delta }s`$ value but we will add experimental point $`a_8=0.58\pm 0.1`$ with enhanced (to 3$`\sigma `$) error. That means we will simply add to $`\chi ^2`$ corresponding to experimental points for spin asymmetries the term connected with experimental point from hyperon decays. The parameters of our three new fits (called $`g^{}`$, $`A_1^{}`$ and $`A_2^{}`$ ) are now presented in Table 2 and results in the measured region in $`0.003x1`$ in the Table 2a.

Table 2. The parameters of three new fits calculated at $`Q^2=1\mathrm{GeV}^2`$

| fit | $`a_1`$ | $`a_2`$ | $`a_4`$ | $`b_2`$ | $`b_3`$ | $`c_1`$ |
| --- | --- | --- | --- | --- | --- | --- |
| fit $`g^{}`$ | $`0.61`$ | $`6.84`$ | $`16.83`$ | $`1.86`$ | $`0.13`$ | $`0.31`$ |
| fit $`A_1^{}`$ | $`0.56`$ | $`5.51`$ | $`14.73`$ | $`1.65`$ | $`0.21`$ | $`0.34`$ |
| fit $`A_2^{}`$ | $`0.50`$ | $`5.39`$ | $`14.72`$ | $`1.98`$ | $`0.29`$ | $`0.33`$ |

| fit | $`c_2`$ | $`c_{2s}`$ | $`c_3`$ | $`d_2`$ | $`d_3`$ | $`\chi ^2/N_{DF}`$ |
| --- | --- | --- | --- | --- | --- | --- |
| fit $`g^{}`$ | $`4.15`$ | $`0.54`$ | $`1.04`$ | $`34.5`$ | $`102.2`$ | $`0.88`$ |
| fit $`A_1^{}`$ | $`3.67`$ | $`0.63`$ | $`0.64`$ | $`15.8`$ | $`43.5`$ | $`0.80`$ |
| fit $`A_2^{}`$ | $`4.15`$ | $`0.56`$ | $`0.44`$ | $`14.2`$ | $`27.6`$ | $`0.84`$ |

There are some small changes in the parameters in comparison to the fits g, $`A_1`$, $`A_2`$ and the biggest change is in $`c_{2s}`$ the parameter responsible for the strange sea. Strange quark contribution is not well determined by the polarised deep inelastic data alone and it is easy by additional experimental point on $`a_8`$ from hyperon decays to shift the value of $`a_8`$ from nearly zero to correct experimental value with only small changes in non strange parton parameters. Comparing Table 1a and Table 2a we can see what is the influence of this additional experimental point $`a_8`$ on integrated parton densities for our three fits. This additional experimental point causes shifts of integrated parton values. As is seen from Table 2a the valence non-strange quark distributions nearly cancel and the value of $`a_8=\mathrm{\Delta }u+\mathrm{\Delta }d2\mathrm{\Delta }s=0.58`$ is built up from relatively high sea contributions. The values of $`\mathrm{\Gamma }_1^p`$, $`\mathrm{\Gamma }_1^n`$ and $`a_3`$ do not change in comparison to previous fits. One can easily calculate from the Table 2a that $`\mathrm{\Delta }\mathrm{\Sigma }=0.24`$ in the fits $`g^{}`$ and $`A_1^{}`$.

Table 2a. The values of quark and gluon polarisations at $`Q^2=1\mathrm{GeV}^2`$ for three new fits (where one includes experimental point from hyperon decays).

| fit | $`\mathrm{\Delta }u`$ | $`\mathrm{\Delta }d`$ | $`\mathrm{\Delta }u_v`$ | $`\mathrm{\Delta }d_v`$ | $`2\mathrm{\Delta }\overline{u}`$ | $`2\mathrm{\Delta }\overline{d}`$ | $`2\mathrm{\Delta }\overline{s}`$ | $`\mathrm{\Delta }G`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| fit $`g^{}`$ | $`0.79`$ | $`0.44`$ | $`0.47`$ | $`0.60`$ | $`0.32`$ | $`0.16`$ | $`0.11`$ | $`0.16`$ |
| fit $`A_1^{}`$ | $`0.80`$ | $`0.44`$ | $`0.57`$ | $`0.57`$ | $`0.23`$ | $`0.13`$ | $`0.12`$ | $`0.01`$ |
| fit $`A_2^{}`$ | $`0.80`$ | $`0.43`$ | $`0.54`$ | $`0.62`$ | $`0.26`$ | $`0.19`$ | $`0.11`$ | $`0.19`$ |

However it is not clear whether the general conclusion that the sea contributions for quarks ( both non-strange and strange) are very big in the measured $`x`$ region is correct. It is specific for our model that the leading singularity for polarised quark valence and sea contributions are comparable. That means that splitting into valence and sea contributions could be not well determined in our fits. That of course could be connected with the functional form of polarised parton densities used by us but what is important is to reproduce the results of the experiment in the measured region. That means that only $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ values can be determined using our parametrisation from the polarised deep inelastic data and not valence and sea contributions separately. The value of $`\mathrm{\Delta }s`$ is determined only with additional $`a_8`$ value from hyperon decays. In our model $`\mathrm{\Delta }G`$ is small. The integrated values $`\mathrm{\Delta }u=0.80`$, $`\mathrm{\Delta }d=0.44`$ are actually not unexpected. Similar values follow from the other models. In other fits  with completely different assumptions for example when we assume the values of $`a_3`$ and $`a_8`$ by fixing the parameters of the fitted parton distributions (normalisation constants) and with the assumption of $`SU(3)`$ symmetry for quark sea we get (using completely different parametrisation from that we use for polarised parton densities):

$`\mathrm{\Delta }u_v\mathrm{\Delta }d_v`$ $`=`$ $`1.26,`$
$`\mathrm{\Delta }u_v+\mathrm{\Delta }d_v`$ $`=`$ $`0.58`$ (7)
If follows that $`\mathrm{\Delta }u_v=0.92`$, $`\mathrm{\Delta }d_v=0.34`$ and in order to get $`\mathrm{\Delta }\mathrm{\Sigma }=0.20`$  we get $`2\mathrm{\Delta }\overline{u}=2\mathrm{\Delta }\overline{d}=2\mathrm{\Delta }\overline{s}=0.13`$ and that means $`\mathrm{\Delta }u=0.79`$, $`\delta d=0.47`$ ( the values not very different from our values). In such models sea contribution is relatively big and negative contrary to our model where we have at least in the measured region big and positive sea contribution. This type of splitting into big valence and relatively big negative sea contribution is caused by the assumptions of the model. We have not made such assumptions taking $`a_8`$ as additional experimental point. Our solution with relatively small valence contribution and relatively large positive sea contributions is different but the values of $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ in both models are very close. It seems that our assumptions are less restrictive. The fact that we can get completely different splitting into valence and sea contributions using the same experimental data shows that this splitting is not well determined by experimental data what is more reliable are distributions of $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ and their integrated values.

We have made fits to data on polarised structure functions on proton, neutron and deuteron targets and we have determined polarised parton distributions. These fits are compared with our previous fits for corresponding asymmetries (improved by usage of recent data from E155 proton experiment in SLAC). As a check we also have made fits to spin asymmetries for all (non averaged in $`Q^2`$) data on spin asymmetries. These fits lead to very similar results with small integrated gluon contribution. The fits were made without inclusion of information on $`a_3`$ and $`a_8`$ from hyperon decays and then repeated with additional experimental point on $`a_8`$. In the first case $`a_8`$ is close to zero and $`\mathrm{\Delta }s`$ is positive. In the second case additional experimental point on $`a_8`$ changes practically only parameters of strange quark and causes small shifts in other parameters. The value of $`a_3`$ at least in the measured region of $`x`$ without any assumptions comes out very close to experimental value. It seems that with the parametrisation used by us only $`\mathrm{\Delta }u(x)`$ and $`\mathrm{\Delta }d(x)`$ are well determined (not the splitting into valence and sea parts). Polarised strange quark distributions, gluon distributions and also $`\overline{u}\overline{d}`$ splitting are not well determined by polarised deep inelastic experimental data.


text:
# Structure functions and energy dissipation dependence on Reynolds number

## I Introduction

There is considerable experimental  and numerical  evidence that the longitudinal velocity difference structure functions

$$S_q(\mathrm{})\delta u(\mathrm{})^q=C_q\epsilon ^{q/3}\mathrm{}^{q/3}\left(\frac{\mathrm{}}{L}\right)^{\zeta _qq/3}$$
(1)
are affected by intermittency corrections in the scaling exponents $`\zeta _q`$ which deviate from the Kolmogorov self-similar prediction $`\zeta _q=q/3`$  ($`C_q`$ is a constant possibly depending on the Reynolds number, $`\epsilon `$ is the turbulent kinetic energy dissipation, $`\mathrm{}`$ is a generic length scale and $`L`$ is the integral length scale). The estimate of $`\zeta _q`$ is now available for a wide range of Reynolds numbers and different flow configurations. The analysis of experimental data show that structure functions display the scaling behavior (1) for sufficiently high Reynolds number, $`Re`$, and that in such conditions the longitudinal exponents $`\zeta _p`$ become independent on $`Re`$ . There is also experimental support for the scaling exponents to be $`Re`$-independent also at small $`Re`$ (when the scaling behavior in (1) is even not observable) if one make use of the so-called extended self similarity analysis (ESS) . However, the precise relation between this empirical result and the classical scaling behavior (1), expected from dimensional arguments, is not well understood. We will not further discuss this point here.

The basic fundamental property of fully developed turbulence is that the average energy dissipation $`\epsilon `$ entering in (1) is asymptotically $`Re`$-independent, when adimensionalized with large scale variables . On the other hand, it is well known that the statistics of the local energy dissipation, which assuming isotropy can be defined in terms of its $`1D`$ surrogate

$$\epsilon (x)=15\nu \left(\frac{u}{x}\right)^2,$$
(2)
strongly depends on $`Re`$, becoming more and more intermittent with increasing $`Re`$ . This effect is also reflected in the tails of the pdf of velocity increments at very small separations . A statistical description of the behavior of the local energy dissipation (2) can be obtained in terms of quantities defined in the inertial range, such as (1), within the multifractal approach. This approach, originally introduced as a phenomenological model for the inertial range statistics , has been extended to the prediction of dissipative scale statistics . In particular, this approach predicts a scaling behavior of the moments of (2) with $`Re`$ , which has been recently measured in a simplified model of turbulence .

In the present paper we investigate the Reynolds dependence of the statistics of (2) in a experimental water jet. We have examined a series of data obtained from LDA measurements which cover about one decade of Reynolds numbers. The quality of the statistics allows us to compute high order structure functions with reasonable accuracy and to partially resolve the dissipative scales. Our main result is that the dependence of the statistics of dissipation on the Reynolds number is found to be consistent with the multifractal prediction obtained by assuming a fluctuating dissipative scale.

The remaining of the paper is organized as follows. In Section II we introduce the theoretical models for the statistics of dissipation. In Section III, the experimental set-up and the data analysis procedure. Section IV is devoted to the presentation of the results, whereas concluding remarks are given in Section V.

## II Multifractal description of energy dissipation

A phenomenological description of intermittency is given by the multifractal model . This model introduces a continuous set of scaling exponents $`h`$ which locally relates the velocity fluctuations at scale $`\mathrm{}`$ entering in (1) with a large scale velocity fluctuation $`u^{}`$:

$$\delta u(\mathrm{})u^{}\left(\frac{\mathrm{}}{L}\right)^h.$$
(3)
The local exponent $`h`$ is realized with a probability which scales with $`(\mathrm{}/L)^{Z(h)}`$ where $`Z(h)`$ is the codimension of the fractal set on which the $`h`$-scaling holds. The scaling exponent of structure functions (1) are obtained by a steepest descent argument over $`h`$:

$$\zeta _q=\underset{h}{inf}\left[qh+Z(h)\right].$$
(4)
The scaling region of (3) is bounded from below by the Kolmogorov dissipative scale $`\eta `$ at which dissipation starts to dominate, i.e. at which the local Reynolds number is of order $`1`$:

$$\frac{\eta \delta u(\eta )}{\nu }1.$$
(5)
From (3) and (5) one obtains that in the multifractal description of intermittency the dissipative scale is a fluctuating quantity, i.e. depends on the local scaling exponent $`h`$ according to

$$\eta L\left(\frac{u^{}L}{\nu }\right)^{\frac{1}{1+h}}LRe_\lambda ^{\frac{2}{1+h}}$$
(6)
where $`Re_\lambda =u^{}\lambda /\nu `$ is the Reynolds number based on the Taylor microscale $`\lambda =\sqrt{15\nu u^2/\epsilon }`$.

Below the dissipative scale the flow can be assumed smooth and one can replace the derivative in (2) with

$$\epsilon =15\nu \left(\frac{\delta u(\eta )}{\eta }\right)^2.$$
(7)
Assuming that the multifractal model can be pushed down to the dissipative scale, one can evaluate the statistics of (2) by inserting (3) and (6) into (7). One ends with the expression

$$\epsilon ^p\epsilon ^p𝑑\mu (h)Re_\lambda ^{2[3php+Z(h)]/(1+h)}\epsilon ^pRe_\lambda ^{2\theta _p}$$
(8)
where again the integral has been evaluated by a steepest descent argument as

$$\theta _p=\underset{h}{inf}\left[\frac{3php+Z(h)}{1+h}\right]$$
(9)
The standard inequality in the multifractal model (following from the exact result $`\zeta _3=1`$), $`Z(h)13h`$ , implies for (9) $`\theta _1=0`$ which is nothing but the request of finite non-vanishing dissipation in the limit $`Re\mathrm{}`$. For $`p>1`$, $`\theta _p<0`$, i.e. the tail of the distribution of $`\epsilon `$ becomes wider with the Reynolds number.

Let us remark that the prediction (9) is based on the assumption of a fluctuating dissipative scale $`\eta `$ according to (6). If one assume, on the contrary, that dissipation enters in (7) as an average quantity one ends up with a different predictions for the exponents $`\theta _p`$ . Numerical simulations with a simplified model of turbulent cascade have shown that the exponents (6) are indeed observed and the alternative prediction is ruled out . In the next Section, we will see that also our experimental data are in agreement with prediction (9).

In the following we will consider data analysis of one component of the velocity in water jets at different Reynolds numbers. Because of the discretization of the acquisition, we are forced to replace spatial derivatives with velocity differences at small scales. The key quantity for our discussion is a generalization of (2) over a finite scale $`\mathrm{}`$:

$$E(\mathrm{})15\nu \left(\frac{\delta u(\mathrm{})}{\mathrm{}}\right)^2.$$
(10)
$`E(\mathrm{})`$ is a convenient definition of surrogate energy dissipation if the scale $`\mathrm{}`$ is sufficiently small. The average dissipative scale dependence on Reynolds number can be obtained from (7) and (1) as

$$\overline{\eta }LRe_\lambda ^{\frac{2}{2\zeta _2}}.$$
(11)
In the analysis of experimental data it is convenient to normalize separations with $`\overline{\eta }`$ and we will consider (10) at fixed $`\mathrm{}^{}=\mathrm{}/\overline{\eta }`$.

In the limit of very small separations (i.e. $`\mathrm{}^{}1`$) (10) recovers the finite difference representation of the 1D-energy dissipation (7) and thus, from (8)

$$\frac{E(\mathrm{})^p}{E(\mathrm{})^p}=\frac{\epsilon ^p}{\epsilon ^p}Re_\lambda ^{2\theta _p}$$
(12)
with the exponents $`\theta _p`$ given by (9).

In the case of separations in the inertial range, from (1) we have $`E(\mathrm{})^p/E(\mathrm{})^p\mathrm{}^{\zeta _{2p}p\zeta _2}`$. Because $`\mathrm{}=\mathrm{}^{}\overline{\eta }`$, using (11) we end with the prediction

$$\frac{E(\mathrm{})^p}{E(\mathrm{})^p}Re_\lambda ^{\frac{2\zeta _{2p}2p\zeta _2}{2\zeta _2}}$$
(13)
Let us remark that the exponents in (12) and (13) are expected to be not very different (they are both zero for non-intermittent turbulence) and thus a discrimination between the two prediction require good accuracy.

## III Experimental set-up and data analysis

The experimental setup consists of a water jet in a closed circuit facility as shown in Figure 1. A centrifugal pump moves water from a primary tank into a settling chamber which is equipped with a valve to damp oscillations due to the pump. A series of contractions leads to a 1.5 m long pipe, with an inner diameter of 14 cm which is followed by the final contraction (1:50 in area) to the jet. The jet (diameter $`d=20`$ mm) exits into a large water-filled tank (height $`30d`$, width $`30d`$, length $`60d`$) from which the water returns to the primary tank. The pipe, the contraction and large tank are made of perspex to allow optical access to the flow. At the jet exit, the flow is axisymmetric and has no swirl; preliminary measurements also confirmed that it is unaffected by any external forcing due to the pump. The jet has a top-hat velocity profile at the nozzle exit with a boundary layer shape parameter (defined as the ratio between the displacement and momentum thickness at the outlet) equal to $`3.29`$ and a turbulent intensity equal to about $`0.021`$ on the jet axis.

Velocity measurements are performed by means of a forward-scatter Laser Doppler Anemometer (LDA) equipped with two Bragg cells. The fringe spacing is $`3.416\mu m`$ and the measurement volume size is about $`0.1mm`$, $`0.1mm`$, $`0.8mm`$ along the $`x`$, $`y`$ and $`z`$ axes respectively. The LDA data are randomly distributed in time; therefore, they are resampled by using a linear interpolation to obtain evenly spaced samples and also to provide unbiased statistics . The value of $`Re_\lambda `$ is sufficiently large to expect the LDA noise not to affect the behavior of structure functions in the IR . The LDA data on longitudinal and transverse structure functions and on scaling exponents were compared to data obtained by Hot Wire Anemometry in a similar jet; the agreement was good up to $`8^{th}`$ order .

Measurements were made at $`x/d40`$ ($`x`$ is measured from the nozzle exit plane), where the flow field may be considered to be approximately self-preserving and isotropic . The jet exit velocity $`U_0`$ was selected so that the exit Reynolds number $`Re_dU_0d/\nu `$ changes from about $`2\times 10^4`$ to $`2\times 10^5`$. As a consequence, $`Re_\lambda `$ changes from about $`200`$ to almost $`2000`$ at the measurement location. The number of collected samples is about $`10^6`$ in each run. As usually done in practice, we resort to Taylor hypothesis to transform time differences into space differences by means of the local mean velocity $`U`$. The number of independent samples, given by $`UT_S/2L`$ (where $`T_S`$ is the total record duration), is about $`10^4`$. Probability density functions of the longitudinal velocity increments have been calculated at different values of $`\mathrm{}^{}`$. The distributions indicate that the number of samples is adequate for achieving a closure of the integrand associated with the structure functions $`S_q(\mathrm{})`$ at least up to $`q=8`$. In Table I we summarize some parameters of the experiments.

Because the acquisition rate is constant for the different runs ($`\delta t=5\times 10^4s`$) the smallest resolved scale $`\delta x=U\delta t`$ varies with $`Re_\lambda `$. By increasing $`Re_\lambda `$ we have both an increasing of $`\delta x`$ (as $`Re_\lambda ^2`$) and a decreasing of the dissipative scale $`\overline{\eta }`$ according to (11). Despite these limitations, we are confident that are able to resolve, at least partially, the dissipative range in all the reported runs. In the extreme case $`Re_\lambda =1840`$ we are able to resolve down to $`\mathrm{}^{}10`$ which is just at the border of the dissipative range . In Figure 2, the second and fourth-order structure functions obtained for $`Re_\lambda =230`$ and for $`Re_\lambda =1840`$ compensated with the scaling behavior (1) are plotted (Kolmogorov length and velocity scales are used to adimensionalise the horizontal and vertical axes). The large scale asymptotic values are also functions of the Reynolds number and, for the second-order structure functions, increase almost linearly as $`Re_\lambda `$; this result is in agreement with that simply derived by dimensional arguments on rms and Kolmogorov velocities, i.e. that $`u^2/\delta u(\eta )^2Re_\lambda `$. Scaling exponents can be derived from structure functions of Figure 2 within the inertial range (approximately from $`\mathrm{}80`$ to $`\mathrm{}500`$, using a criterion based on $`3\%`$ difference from the maximum in the third-order structure function divided by $`\mathrm{}`$). The results ($`\zeta _2=0.708\pm 0.020`$ and $`\zeta _4=1.265\pm 0.035`$ for $`Re_\lambda =230`$) agree with those obtained numerically and experimentally by other authors .

## IV Results and comparison with theoretical predictions

In Figure 3 we plot the average surrogate energy dissipation $`E(\mathrm{})`$, computed from (10) and adimensionalized with large scale quantity $`u^3/L`$, as a function of $`Re_\lambda `$. A constant value of energy dissipation is expected from energy balance arguments . Thus, the observed independence of $`E(\mathrm{})`$, computed at fixed $`\mathrm{}^{}=10`$, on $`Re_\lambda `$ is a confirmation that the dissipative scales are resolved. We can give a different estimation of the dissipation on the basis of the energy spectrum. The advantage in this case is that we can reduce the noise by interpolation of the spectrum down to the dissipative scales . The result for our data, also plotted in Figure 3, is consistent with the previous one thus confirming that the evaluation of the energy dissipation statistics from the present data is reliable.

In order to verify the prediction for $`\theta _p`$ of Section II, the scaling laws (12) of the moments of $`E(\mathrm{})`$ with $`Re_\lambda `$ must be computed. In Figure 4, an example of this scaling, computed at $`\mathrm{}^{}=10`$, is given for different values of $`p`$. We are confident that at least up to $`p=3`$ the scaling exponents can be evaluated with an error on the fit less than $`5\%`$.

The resulting exponents $`\theta _p`$ are given in Figure 5. They have been computed at two different values of the separation $`\mathrm{}^{}=10,100`$. As expected, for separations moving from the dissipative range into the inertial range, the curve of the exponents $`\theta (p)`$ becomes flatter and approaches the prediction (13). In Figure 5 we also plot the two predictions (9) and (13). These are been obtained by computing the structure function scaling exponents $`\zeta _q`$ (1) as discussed above, at the largest $`Re_\lambda `$ available. Given the exponents $`\zeta _q`$, the codimension $`Z(h)`$ has been obtained by numerical inversion of (4). Finally, (9) is used to predict the values of $`\theta _p`$.

The agreement of our data with the two predictions (12) and (13) is remarkable. The value of $`\theta _p`$ obtained from $`E(\mathrm{})`$ at $`\mathrm{}^{}=10`$ are very close to prediction (9-12). Some deviations are observed, starting from moment $`p2.5`$. These are probably due to the relatively large value of $`\mathrm{}^{}=10`$, which is at the border of the dissipative range and to the statistics. Indeed, according to (10), the computation of $`\theta _3`$ corresponds roughly to the computation of the $`6`$-th order structure function.

Let us remark again that the difference between predictions (9) and (13) is based on the assumption of a dissipative scale fluctuating with the local energy dissipation. Our experimental data demonstrate that this is indeed the case.

## V Remarks and conclusions

In this paper, the scaling of velocity increments at very small scales in a turbulent flow is investigated with special focus on Reynolds number dependence. Within the multifractal framework it is possible to derive a relation between the scaling of the energy dissipation statistics with the Reynolds number and the scaling exponents of the velocity structure functions at fixed Reynolds number. This relation is found to fit very well experimental data taken into a water jet over about a decade of $`Re_\lambda `$. The one-dimensional surrogate of energy dissipation is estimated from velocity differences at very small scale, at the border of the dissipative range. It is shown that our statistics changes from prediction (12) to prediction (13) as the separation $`\mathrm{}`$ moves into the inertial range. The difference between the two predictions is the demonstration of the fluctuation of the dissipative scale.

A Reynolds number dependence of the moments of the velocity increments also in the limit of very high values ($`Re_\lambda >1000`$) has been recently reported by different authors . Due to the fact that the Reynolds number is expected to affect mainly the dissipative scales, the previous findings confirm experimentally the connection between small and inertial range scalings. This could support ESS (for example in respect to the enhanced range of scales involved), but simultaneously invalidate it (in respect of the independence on the Reynolds number). If these results are dependent on the still finite value of the Reynolds number or on the anisotropy of the flow field is left as a subject for future investigations.


text:
# Untitled Document

Paper is withdrawn due to errors (Proposition 11.4 is false). Section 12 is false, and the main statement is true only for bundles $`B`$ with $`c_1(B)`$ $`SU(2)`$-invariant. The rest is OK, but it’s presented in a better way in math.AG/0604303


text:
# Rotational Brownian Motion of a Massive Binary

## 1 Introduction

Two sorts of Brownian motion are defined by statistical physicists. Translational Brownian motion is the irregular motion exhibited by a massive particle as it collides with molecules in a fluid. Collisions also affect the spin and orientation of particles; theories of rotational Brownian motion are concerned with the evolution of these quantities due to encounters. Rotational Brownian motion was first discussed by Debye (1913, 1929) in the context of dielectric theory. The polarization of a dielectric material is a competition between torques due to the imposed field, which tend to align polar molecules, and collisions, which tend to destroy the alignment.

The likely existence of black holes with masses much greater than those of stars has led to a renewed interest in Brownian motion among astrophysicists (Merritt, 2001; Milosavljevic & Merritt, 2001; Dorband, Hemsendorf & Merritt, 2001). Translational Brownian motion of black holes in galactic nuclei may be directly observable (e.g. Reid et al. (1999); Backer & Sramek (1999)), allowing constraints to be placed on black hole masses. Translational Brownian motion may also affect the rate of tidal disruption of stars by a nuclear black hole (e.g. Young (1977)) or the efficiency with which a binary black hole can interact with stars as it wanders through a nucleus (Quinlan & Hernquist, 1997; Milosavljevic & Merritt, 2001). In Paper I (Merritt, 2001), the translational Brownian motion of a massive binary was discussed; it was shown that the amplitude of the wandering can be increased by a modest factor compared with that of a single point mass due to exchange of energy between field stars and the binary.

Here we consider the second kind of Brownian motion experienced by a binary, the irregular variation of the orientation of the binary’s spin axis due to encounters. A single field star that passes within a distance $`a`$ of the binary, with $`a`$ the binary’s semi-major axis, will exchange orbital angular momentum with the binary, leading both to a change in the binary’s orbital eccentricity as well as a change in the orientation of the binary’s spin axis. Repeated encounters will cause the binary’s orbital eccentricity to evolve (usually in the direction of increasing eccentricity) and will also cause the orientation of the binary to undergo a random walk. The first process has been discussed extensively (e.g. Mikkola & Valtonen (1992); Quinlan (1996)); the second, discussed here for the first time, will be called rotational Brownian motion, by analogy with the similar process that occurs in a fluid of polar molecules (e.g. McConnell (1980)).

The orbital orientation of a black-hole binary influences the direction of the spin axis of the single black hole that forms after coalescence of the binary due to emission of gravitational radiation (e.g. Flanagan & Hughes (1998)). The spin axis, in turn, is thought to determine the plane of the inner accretion disk which forms around the black hole (Bardeen & Petterson, 1975), the direction of jets launched from the accretion disk (Rees, 1978), etc. Thus changes in the orientation of a black-hole binary as it coalesces may have observable consequences.

In §2 the conservation equations for interaction of a field star with a massive binary are presented; most of these have been seen before and were used in earlier studies to derive the rates of change of the binary’s separation and eccentricity. The rotational diffusion equation is presented in §3 and the order-of-magnitude of the rotational diffusion coefficient $`\mathrm{\Delta }\vartheta ^2`$ is derived. §4 describes the numerical derivation of $`\mathrm{\Delta }\vartheta ^2`$ from scattering experiments. $`\mathrm{\S }5`$ presents solutions of the diffusion equation and derives a relation between the degree of hardening of a binary and the change in its orientation.

## 2 Encounter Kinematics

Consider an encounter of a single particle (the “field star”) with a binary. Assume that the binary remains bound during the encounter, appropriate if the mass of the field star is much less than the mass of each of the components of the binary or if the binary is sufficiently hard. At early times, the field star has velocity $`𝐯_0`$ with respect to the center of mass of the field star-binary system and its impact parameter is $`p`$. A long time after the encounter, the velocity $`𝐯`$ of the field star attains a constant value. Conservation of linear momentum implies that the change $`\delta 𝐕`$ in the velocity of the binary’s center of mass is given by

$$\delta 𝐕=\frac{m_f}{M_{12}}\delta 𝐯$$
(1)
where $`m_f`$ is the mass of the field star and $`M_{12}=M_1+M_2`$ is the mass of the binary. This relation was used in Paper I in conjunction with scattering experiments to investigate the translational Brownian motion of a binary due to encounters.

The energy of the field star-binary system, expressed in terms of pre-encounter quantities, is

$`E_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_fv_0^2+{\displaystyle \frac{1}{2}}M_{12}V_0^2{\displaystyle \frac{GM_1M_2}{2a_0}}`$ (2a)
$`=`$ $`{\displaystyle \frac{1}{2}}m_f\left(1+{\displaystyle \frac{m_f}{M_{12}}}\right)v_0^2{\displaystyle \frac{GM_1M_2}{2a_0}}`$ (2b)
with $`𝐕_0=(m_f/M_{12})𝐯_0`$ the initial velocity of the binary’s center of mass and $`a_0`$ the binary’s initial semi-major axis. After the encounter,

$$E=\frac{1}{2}m_f\left(1+\frac{m_f}{M_{12}}\right)v^2\frac{GM_1M_2}{2a}$$
(3)
and $`E=E_0`$, so that

$`\delta \left({\displaystyle \frac{1}{a}}\right)`$ $`=`$ $`{\displaystyle \frac{m_f(v^2v_0^2)}{GM_1M_2}}\left(1+{\displaystyle \frac{m_f}{M_{12}}}\right)`$ (4a)
$``$ $`{\displaystyle \frac{m_f(v^2v_0^2)}{GM_1M_2}}`$ (4b)
where the latter relation assumes $`M_{12}m_f`$. This result was used by Hills (1983, 1992), Mikkola & Valtonen (1992), and Quinlan (1996) to compute the hardening rate, $`(d/dt)(1/a)`$, of a massive binary due to encounters.

The angular momentum of the field star-binary system about its center of mass, expressed in terms of pre-encounter quantities, is

$$_0=m_f\left(1+\frac{m_f}{M_{12}}\right)\mathrm{}_0+\mu _{12}\mathrm{}_{b0}$$
(5)
where $`\mu _{12}M_1M_2/M_{12},\mathrm{}_0p𝐯_0`$ and $`\mathrm{}_{b0}_{b0}/\mu _{12}`$ with $`_b`$ the binary’s spin angular momentum. Conservation of angular momentum during the encounter gives

$`\delta \mathrm{}_b`$ $`=`$ $`{\displaystyle \frac{m_f}{\mu _{12}}}\left(1+{\displaystyle \frac{m_f}{M_{12}}}\right)\delta \mathrm{}`$ (6a)
$``$ $`{\displaystyle \frac{m_f}{\mu _{12}}}\delta \mathrm{}.`$ (6b)
Changes in $`|\mathrm{}_b|`$ correspond to changes in the binary’s orbital eccentricity $`e`$ via the relation $`e^2=1\mathrm{}_b^2/GM_{12}\mu ^2a`$ (Mikkola & Valtonen, 1992; Quinlan, 1996). Changes in the direction of $`\mathrm{}_b`$ correspond to changes in the orientation of the binary, leading to rotational diffusion.

## 3 Rotational Diffusion Equation

Let $`F(\theta ,\varphi ,t)d\mathrm{\Omega }`$ be the probability that the spin axis of the binary is oriented within the solid angle $`d\mathrm{\Omega }`$ at time $`t`$. We seek an equation describing how $`F`$ evolves with time due to encounters of field stars with the binary, if each encounter is assumed to produce a negligibly small change in $`\theta `$ and $`\varphi `$. If we imagine that $`_b`$ is initially oriented parallel to the $`\theta =0`$ axis, and that the encounters leading to changes in the orientation of the binary are isotropic in velocity and direction, then $`F`$ will evolve in such a way as to remain a function of $`\theta `$ alone, $`F=F(\theta ,t)`$. These are the same assumptions made by Debye (1929) in his theory of the rotational Brownian motion of spherical molecules. Debye showed that the evolution equation for $`F`$, in the absence of an external torque, is

$$\frac{F}{t}=\frac{1}{\mathrm{sin}\theta }\frac{}{\theta }\left(\mathrm{sin}\theta \frac{\mathrm{\Delta }\vartheta ^2}{4}\frac{F}{\theta }\right)$$
(7)
with $`\mathrm{\Delta }\vartheta ^2`$ the rotational diffusion coefficient,

$$\mathrm{\Delta }\vartheta ^2=\mathrm{\Psi }(d\mathrm{\Omega },d\mathrm{\Omega }^{})\vartheta ^2𝑑\mathrm{\Omega }^{};$$
(8)
$`\mathrm{\Psi }(d\mathrm{\Omega },d\mathrm{\Omega }^{})d\mathrm{\Omega }^{}`$ is the probability that, during a unit interval of time, a binary whose angular momentum $`_b`$ is directed toward $`d\mathrm{\Omega }`$ will reorient itself such that $`_b`$ lies within $`d\mathrm{\Omega }^{}`$, and $`\vartheta `$ is the angular separation between $`d\mathrm{\Omega }`$ and $`d\mathrm{\Omega }^{}`$.

From its definition, $`\mathrm{\Delta }\vartheta ^2`$ is the sum, over a unit interval of time, of $`(\delta \vartheta )^2`$ due to encounters with field stars. It can be expressed in terms of mean changes in the components of the binary’s angular momentum as follows. Let the binary be initially oriented with its spin vector along the $`z`$ axis, $`\mathrm{}_{b0}=\mathrm{}_{b0}𝐞_z`$. After an encounter, the angle $`\vartheta `$ between initial and final spins is

$$\mathrm{cos}\vartheta =\frac{\mathrm{}_{b,z}}{\mathrm{}_b}=\frac{\mathrm{}_{b,z}}{\sqrt{\mathrm{}_{b,x}^2+\mathrm{}_{b,y}^2+\mathrm{}_{b,z}^2}}.$$
(9)
Expanding equation (9) to second order in $`\vartheta `$ and taking averages, we find

$$\mathrm{\Delta }\vartheta ^2=\frac{\mathrm{\Delta }\mathrm{}_{b,x}^2+\mathrm{\Delta }\mathrm{}_{b,y}^2}{\mathrm{}_{b0}^2}$$
(10)
which can further be expressed in terms of changes in the angular momentum components of the field star via equation (6b) as

$`\mathrm{\Delta }\vartheta ^2`$ $`=`$ $`{\displaystyle \frac{m_f^2}{\mu _{12}^2}}{\displaystyle \frac{\mathrm{\Delta }\mathrm{}_x^2+\mathrm{\Delta }\mathrm{}_y^2}{\mathrm{}_{b0}^2}}`$ (11a)
$`=`$ $`{\displaystyle \frac{m_f^2}{GM_1M_2\mu _{12}a(1e^2)}}\left[\mathrm{\Delta }\mathrm{}_x^2+\mathrm{\Delta }\mathrm{}_y^2\right].`$ (11b)
The last relations assume $`M_{12}m_f`$. Specializing to the case of an equal-mass, circular-orbit binary, equation (11b) becomes

$$\mathrm{\Delta }\vartheta ^2=16\left(\frac{m_f}{M_{12}}\right)^2\frac{\mathrm{\Delta }\mathrm{}_x^2+\mathrm{\Delta }\mathrm{}_y^2}{GM_{12}a}.$$
(12)
Equation (12) is a prescription for computing the rotational diffusion coefficient from scattering experiments. This prescription will be applied in the next section.

We can estimate the order of magnitude of the rotational diffusion coefficient by noting that $`\sqrt{GM_{12}/a}=V_{bin}`$, the relative velocity of $`M_1`$ and $`M_2`$, and that a single close encounter produces a change in $`\mathrm{}`$ of the field star of order $`V_{bin}\times a`$. The rate of close encounters is $`2\pi GM_{12}n_fa/\sigma _f`$ with $`n_f`$ and $`\sigma _f`$ the number density and 1d velocity dispersion of the field stars (Paper I). Thus $`\mathrm{\Delta }\mathrm{}^2(aV_{bin})^2\times 2\pi GM_{12}n_fa/\sigma _f`$ and

$$\mathrm{\Delta }\vartheta ^232\pi \frac{m_f}{M_{12}}\frac{G\rho _fa}{\sigma _f}$$
(13)
where $`\rho _fm_fn_f`$.

The same encounters that induce rotational diffusion will also cause the binary to harden, according to equation (4b). The hardening rate may be written

$$\frac{d}{dt}\left(\frac{1}{a}\right)=\frac{G\rho _fH}{\sigma _f}$$
(14)
with $`H15`$ for a hard binary (e.g. Quinlan (1996)). Thus

$$\mathrm{\Delta }\vartheta ^2\frac{32\pi }{H}\frac{m_f}{M_{12}}\left[a\frac{d}{dt}\left(\frac{1}{a}\right)\right]\frac{32\pi }{H}\frac{m_f}{M_{12}}\frac{1}{t_{hard}}$$
(15)
with $`t_{hard}^1a(d/dt)(1/a)`$. Equation (15) suggests that rotational diffusion occurs on a time scale that is of order $`M_{12}/m_f`$ times the hardening time, or that the binary will rotate by an angle of order $`\sqrt{m_f/M_{12}}`$ in the time it takes to harden significantly. This result will be assessed more quantitatively in §5.

## 4 Numerical Computation of the Diffusion Coefficient

The same set of scattering experiments described in Paper I were used here to evaluate $`\mathrm{\Delta }\mathrm{}_x^2`$ and $`\mathrm{\Delta }\mathrm{}_y^2`$. As described more fully in that paper, field stars were treated as massless particles moving in the potential of the binary whose center of mass and orbital parameters remained fixed. Orbits were integrated using the routine D0P853 of Hairer, Nørsett & Wanner (1991), an 8(6)th order embedded Runge-Kutta integrator; each field star was assumed to begin at a position $`(x,y,z)=(\mathrm{},p,0)`$ and was advanced from $`r=\mathrm{}`$ to $`r=50a`$ along a Keplerian orbit about a point mass $`M_{12}`$. The integrations were terminated when the star had moved a distance from the binary that was at least $`100`$ times its initial distance with positive energy. The binary’s orbit was assumed to be circular; the orientation of the binary’s orbital plane with respect to the initial velocity vector of the field star, and the initial phase of the binary, were chosen randomly for each integration. Unless otherwise indicated, distances and velocities are given below in program units of $`a`$ and $`(GM_{12}/a)^{1/2}`$ respectively. The orbit integrations, and the expressions given above relating changes in the field-star’s angular momentum to the orientation of the binary, were checked by carrying out a limited set of experiments using a fully general, three-body integrator and various values for $`m_f/M_{12}`$. The expressions derived above for $`\delta \vartheta `$ etc. were found to be accurately reproduced in the limit of small $`m_f/M_{12}`$.

Consider first the rotational diffusion coefficient corresponding to interactions with stars of a single velocity $`V`$ at infinity. Multiplying $`(\delta \vartheta )^2`$ in a single encounter by $`2\pi pn_fV`$, with $`p`$ the impact parameter and $`n_f`$ the density of field stars, and integrating over $`p`$ gives

$`\delta \vartheta ^2(V,p_{max})`$ $`=`$ $`2\pi n_fV{\displaystyle _0^{p_{max}}}𝑑pp\overline{\delta \vartheta ^2}`$ (16a)
$`=`$ $`2\pi n_fV{\displaystyle _0^{p_{max}}}𝑑pp\left({\displaystyle \frac{\overline{\delta \mathrm{}_{b,x}^2}+\overline{\delta \mathrm{}_{b,y}^2}}{\mathrm{}_b^2}}\right)`$ (16b)
$`=`$ $`{\displaystyle \frac{32\pi n_fV}{GM_{12}a}}\left({\displaystyle \frac{m_f}{M_{12}}}\right)^2{\displaystyle _0^{p_{max}}}𝑑pp\left(\overline{\delta \mathrm{}_x^2}+\overline{\delta \mathrm{}_y^2}\right)`$ (16c)
(cf. equation 12), where $`\overline{\delta \mathrm{}_x^2},\overline{\delta \mathrm{}_y^2}`$ are defined as mean square changes at fixed $`(V,p_{max})`$, averaged over many encounters in which the orientation of $`𝐕`$ with respect to the binary’s spin axis is taken to be random. Figure 1 shows the distribution of field-star velocity changes as a function of $`P`$ for $`V=0.5V_{\mathrm{bin}}`$.

We define a dimensionless, velocity-dependent diffusion coefficient $`D_R`$ as

$`\delta \vartheta ^2(V,p_{max})`$ $`=`$ $`\pi a^2n_fV_{\mathrm{bin}}\left({\displaystyle \frac{m_f}{M_{12}}}\right)^2D_R(V,p_{max}),`$ (17a)
$`D_R(V,p_{max})`$ $``$ $`16\left({\displaystyle \frac{V}{V_{bin}}}\right){\displaystyle _0^{p_{max}/a}}d\left({\displaystyle \frac{p}{a}}\right)^2{\displaystyle \frac{\left(\overline{\delta \mathrm{}_x^2}+\overline{\delta \mathrm{}_y^2}\right)}{a^2V_{bin}^2}}`$ (17b)
with $`V_{\mathrm{bin}}=\sqrt{GM_{12}/a}`$. Figure 2 plots $`D_R`$ as a function of $`V`$ and $`p_{max}`$. The most effective encounters are those with small $`V`$ and $`p_{max}`$; encounters with impact parameters $`p\begin{array}{c}>\hfill \\ \hfill \end{array}\sqrt{2GM_{12}a}/V`$ have a closest-approach distance to the binary of $`r_p\begin{array}{c}>\hfill \\ \hfill \end{array}a`$ and produce almost no change in the binary’s orientation. Hence $`D_R`$ reaches a maximum value at $`p\sqrt{2GM_{12}a}/V`$ before levelling off at larger $`p`$.

Multiplying $`\delta \vartheta ^2`$ by $`f_f(v_f)`$, the distribution of field-star velocities, and integrating $`\mathrm{𝐝𝐯}_f`$ gives $`\mathrm{\Delta }\vartheta ^2`$. We assume a Maxwellian distribution of field star velocities,

$$f_f(v_f)=\frac{1}{(2\pi \sigma _f^2)^{3/2}}e^{v_f^2/2\sigma _f^2},$$
(18)
so that

$$\mathrm{\Delta }\vartheta ^2=\frac{4\pi }{\left(2\pi \sigma ^2\right)^{3/2}}_0^{\mathrm{}}𝑑VV^2e^{V^2/2\sigma _f^2}\delta \vartheta ^2.$$
(19)
This can be written in terms of a dimensionless coefficient $`L`$,

$`\mathrm{\Delta }\vartheta ^2`$ $`=`$ $`{\displaystyle \frac{m_f}{M_{12}}}{\displaystyle \frac{G\rho _fa}{\sigma _f}}L,`$ (20a)
$`L`$ $``$ $`\sqrt{2\pi }S^1{\displaystyle _0^{\mathrm{}}}𝑑zz^2e^{z^2/2}D_R(zV_{\mathrm{bin}}/S,Rp_f).`$ (20b)
$`L`$ is defined in analogy with the dimensionless coefficients $`(H,J,K)`$ that describe respectively the rates of change of the binary’s energy, the rate of mass ejection by the binary, and the binary’s eccentricity growth rate (cf. Quinlan (1996)). Note that $`L`$ is a function of two parameters, $`L=L(R,S)`$, where

$$R\frac{p_{max}}{p_f}=\frac{p_{max}\sigma _f^2}{GM_{12}},S\frac{V_{bin}}{\sigma _f}.$$
(21)
$`R`$ is the maximum impact parameter in units of the radius of gravitational influence of the binary and $`S`$ is the dimensionless hardness of the binary. Figure 3 plots $`L(R,S)`$. At a given hardness $`S`$, $`L`$ reaches its limiting value $`L_{\mathrm{}}L(\mathrm{},S)`$ by $`R3/S`$. For all but the softest binaries, this limiting value is reached before $`R=1`$, which is roughly the value of $`R`$ we expect for a massive binary at the center of a steeply-falling density cusp (Paper I). Hence we ignore the dependence of $`L`$ on $`p_{max}`$ in what follows and write $`L(S)=L(\mathrm{},S)`$. $`L(S)`$ varies from $`25`$ for the softest binaries ($`S1`$) to $`60`$ for $`S10`$.

## 5 Solutions of the Evolution Equation

Equation (7) describes the evolution of $`F(\theta ,t)`$, the probability distribution for the binary’s orientation, given the rotational diffusion coefficient $`\mathrm{\Delta }\vartheta ^2`$. The latter (equation 20) is a function of the binary separation $`a`$, directly via the factor $`G\rho _fa/\sigma _f`$, and indirectly via the weak dependence of $`L`$ on binary hardness (Figure 3). Furthermore the timescale for the binary to reorient itself is shorter than the hardening time $`\left|(1/a)(da/dt)\right|^1`$ (cf. equation 15), hence the change in $`a`$ can not be neglected when solving for the evolution of the orientation.

We nevertheless begin by considering an idealized model in which $`\mathrm{\Delta }\vartheta ^2`$, $`a`$, and the parameters $`(\rho _f,\sigma _f)`$ that describe the background are assumed constant in time. The evolution equation (7) becomes

$$\frac{F}{\tau }=\frac{1}{2}\frac{}{\mu }\left[\left(1\mu ^2\right)\frac{F}{\mu }\right],$$
(22)
where $`\mu \mathrm{cos}\theta `$, $`\tau t/t_0`$, and

$$t_02\left(\frac{M_{12}}{m_f}\right)\frac{\sigma _f}{G\rho _faL}.$$
(23)
Consider the evolution of $`\overline{\mu }(t)=2\pi _1^1F(\mu ^{},t)\mu ^{}𝑑\mu ^{}`$, the expectation value of $`\mu `$. We have

$`{\displaystyle \frac{d\overline{\mu }}{d\tau }}`$ $`=`$ $`2\pi {\displaystyle _1^1}{\displaystyle \frac{F}{\tau }}\mu ^{}𝑑\mu ^{}`$ (24a)
$`=`$ $`\pi {\displaystyle _1^1}\mu ^{}{\displaystyle \frac{}{\mu ^{}}}\left[\left(1\mu ^2\right){\displaystyle \frac{F}{\mu ^{}}}\right]`$ (24b)
$`=`$ $`2\pi {\displaystyle _1^1}\mu ^{}F𝑑\mu ^{}=\overline{\mu }`$ (24c)
or

$$\overline{\mu }(\tau )=\overline{\mu }_0e^\tau .$$
(25)
Hence $`t_0`$ is the time constant for relaxation to a uniform distribution of orientations, $`\overline{\mu }=0`$, given a fixed binary separation and unchanging stellar background. It may be expressed in physical units as

$$t_0=6.1\times 10^9\mathrm{yr}\left(\frac{M_{12}}{10^6m_f}\right)\left(\frac{L}{50}\right)^1\left(\frac{\rho _f}{10^6M_{}\mathrm{Mpc}^3}\right)^1\left(\frac{\sigma _f}{200\mathrm{km}\mathrm{s}^1}\right)\left(\frac{a}{1\mathrm{pc}}\right)^1.$$
(26)
Next consider the more interesting case where $`a`$, as well as the parameters $`\rho _f`$ and $`\sigma _f`$ that describe the field stars, may be changing with time. This time dependence is not easily specified but we can make progress by changing evolution variables from $`t`$ to $`x\mathrm{log}(a/a_0)`$. First combining equations (7) and (20a), we find

$$\frac{F}{t}=\frac{L}{4}\frac{m_f}{M_{12}}\frac{G\rho _fa}{\sigma _f}\frac{}{\mu }\left[\left(1\mu ^2\right)\frac{F}{\mu }\right].$$
(27)
Now changing variables,

$$\frac{F}{t}=\frac{F}{x}\frac{dx}{dt}=\frac{F}{x}\frac{1}{a}\frac{da}{dt}=\frac{F}{x}\frac{G\rho _fHa}{\sigma _f}$$
(28)
(cf. equation (14)) and so

$$\frac{F}{x}=\frac{1}{4}\frac{L}{H}\frac{m_f}{M_{12}}\frac{}{\mu }\left[\left(1\mu ^2\right)\frac{F}{\mu }\right].$$
(29)
This equation gives the evolution of the binary’s orientation in terms of changes in its semi-major axis $`a`$, with no explicit dependence on the parameters that describe the stellar background. Ignoring the weak dependence of $`H`$ and $`L`$ on $`a`$, appropriate for a hard binary, we find as before

$$\overline{\mu }(a)=\overline{\mu }_0\left(\frac{a}{a_0}\right)^{\frac{1}{2}\frac{L}{H}\frac{m_f}{M_{12}}}.$$
(30)
The exponent in this expression is of order $`m_f/M_{12}1`$. Hence we can write

$$\overline{\mu }(a)1+\frac{1}{2}\frac{L}{H}\frac{m_f}{M_{12}}\mathrm{log}\left(\frac{a}{a_0}\right)$$
(31)
where $`\overline{\mu }_0`$ has been set to unity, corresponding to an initial orientation parallel to the $`\theta =0`$ axis.

For the case of binary supermassive black holes, if we define $`a_0`$ as the separation when the binary first forms a bound pair, we expect gravitational radiation coalescence to occur when $`a/a_010^2`$ (Merritt, 2000). Since $`L/2H2`$ for a hard binary, we find for the expectation value of $`\mu `$ at coalescence:

$$\overline{\mu }110\frac{m_f}{M_{12}}.$$
(32)
Writing $`\delta \theta \sqrt{2(1\overline{\mu })}`$, the rms change in the angle defined by the binary’s spin axis, this becomes

$$\delta \theta \sqrt{\frac{20m_f}{M_{12}}},m_fM_{12}.$$
(33)
For $`m_f/M_{12}=10^6`$, e.g. a binary of mass $`10^6M_{}`$ surrounded by $`1M_{}`$ stars, this predicts $`\delta \theta 0.25^{}`$. For a binary of so-called “intermediate-mass” black holes with $`M_{12}10^3M_{}`$, $`\delta \theta 8^{}`$.

Milosavljevic & Merritt (2001) give Information about the rotational Brownian motion of a massive binary in a set of $`N`$-body simulations. They show (in their Figure 11) the evolution of the orientation of an equal-mass, circular-orbit binary in three simulations with $`m_f/M_{12}=\{164,328,655\}`$. The simulations all follow the evolution of the binary until $`a/a_010^1`$; hence equation (32) predicts $`\delta \theta \{14^{},10^{},7^{}\}`$. These values are in excellent agreement with the angular deflections shown in Figure 11 of Milosavljevic & Merritt (2001).

Under what circumstances could the roeorientation angle of a supermassive black hole binary be much greater than the modest values predicted here? Since $`\delta \theta m_f^{1/2}`$, one possibility is for the perturbing objects to be much more massive than $`1M_{}`$; for instance, they could be supergiant stars ($`m_f10^2M_{}`$), star clusters ($`m_f10^3M_{}`$), or giant molecular clouds ($`m_f10^4M_{}`$). Such physically large objects would presumably produce most of their perturbing torques when the binary separation was greater than or comparable to their own size.

An even more extreme possibility is that the perturbing objects are themselves supermassive black holes. Consider for example a binary supermassive black hole at the center of a giant elliptical galaxy. The stellar densities at the centers of such galaxies are low and one might expect the coalescence of the binary to stall due to a shortage of interacting stars (e.g. Milosavljevic & Merritt (2001)). In this case, further decay of the binary might await the infall of a third and subsequent black holes from later galactic mergers (e.g. Valtonen (1996)). Mass ratios $`m_f/M_{12}`$ would then be of order $`10^210^0`$ and signficant reorientation of the binary could occur before coalescence.

Orientations of radio jets in Seyfert galaxies are known to be almost random with respect to the plane of the stellar disk (e.g. Kinney et al. (2000)). If jets are launched parallel to the spin axes of rotating black holes, this implies that the spins of nuclear black holes bear almost no relation to disk spins. While the modest reorientation rates derived here are probably insufficient to explain this phenomenon, a related mechanism may suffice. The spin orientation of a supermassive black hole that forms via repeated mergers with smaller black holes would undergo a random walk due to the random orbital inclinations of the infalling black holes. If nuclear black holes in Seyfert galaxies grew through multiple mergers before the formation of the stellar disks, one might expect their spin axes to be uncorrelated with those of the stars and gas.

This work was supported by NSF grants AST 96-17088 and 00-71099 and by NASA grants NAG5-6037 and NAG5-9046.


text:
# DARK MATTER IS BARYONS

## 1 Cosmoids

The existence of a meteoroid population in highly eccentric orbits was discovered in the results of the three ”dust” experiments carried aboard the PIONEER 10 and 11 missions to the outer planets. A re-examination of the results from the Asteroid/Meteoroid Detector (nicknamed Sisyphus) established beyond statistical doubt that a new population termed cosmoids, a contraction of cosmic meteoroids, had been measured (Dubin & Soberman 1991). The three-way conflict among the conclusions of the three meteoric experiments (Soberman et al. 1976) was resolved by the recognition that all had measured cosmoids (Dubin & Soberman 1991).

The seminal velocity measurements of Grün et al. (1992) established what was previously suspected, that cosmoids were extra-solar in origin. That was confirmed when velocities in excess of 100 km/s were measured for members of a “new population” of terrestrial radio meteors (Taylor et al. 1996) first discovered utilizing VLF (2 MHz) radar (Olsson-Steel & Elford 1987).

The properties that allowed cosmoids to escape detection through more than a half century of intensive search for their predicted existence includes their extremely low comet-nucleus-like albedo (about 2 - 3

## 2 Are Non-baryons Permitted?

An appreciable fraction of non-baryonic dark matter cannot be excluded by the discovery of cosmoids. However, only “Big Bang” theory predicts sizable masses in alternate forms. Apart from detected MACHOS (MAssive Compact Halo Objects) which are likely large gravity compressed cosmoids, the nature of other hypothesized candidates might have been drawn from “Alice In Wonderland” (Cowen 2000). Occam’s razor argues that they should not be assumed if unneeded.

The validity of the “Big Bang” has been long questioned by a small distinguished group of astrophysicists beginning with Fritz Zwicky. He coined its name to deride the concept. For three and one-half decades Halton Arp (1966) has presented hundreds of examples where the red shift conflicts with apparently associated objects. Only Penzias’ and Wilson’s (1965) discovery of the approximate 3 K Cosmic Microwave Background (CMB) shifted most to believe in an expanding universe and fostered a major cosmology industry. Neglected or dismissed is the need for “inflation,” faster-than-light galactic aggregation and most recently the resurrection of Einstein’s cosmological constant investing free space with anti-gravity power.

An alternate explanation for the red shift that ties it generally with distance is the interaction of photons with cosmoids (Soberman & Dubin 2001). In this model the CMB results from a thermal (size-independent) eigenvalue plateau for cosmoids being drawn into the Galaxy. The fractional kelvin discrepancy may be removed by laboratory determination. This contrasts with the decade discrepancy between the measured CMB and prediction (Alpher et al. 1948).