text:
# Effects of Differential Rotation on the Maximum Mass of Neutron Stars

## 1 Introduction

Neutron stars that are newly formed in the coalescence of binary neutron stars (or in supernova collapse) are likely to be differentially rotating (see, e.g., the fully dynamical simulations of Rasio & Shapiro 1992,1994, Shibata & Uryu 2000, Faber, Rasio & Manor 2001, as well as the review of Rasio & Shapiro 1999 and references therein). Differential rotation may play an important role for the stability of these remnants, since it can be very effective in increasing their maximum allowed mass.

Assuming that neutron stars in binaries have individual masses close to 1.4 $`M_{}`$ (Thorsett & Chakrabarty 1999), and assuming that the maximum allowed mass of a non-rotating neutron star is in the range of 1.8 – 2.3 $`M_{}`$ (Akmal, Pandharipande & Ravenhall 1998), one might conclude that the coalescence of binary neutron stars leads to immediate collapse to a black hole. However, both thermal pressure and rotation can increase the maximum allowed mass.

For coalescence from the innermost stable circular orbit, thermal pressure is believed to have a small effect. The maximum mass of uniformly rotating stars is limited by the spin rate at which the fluid at the equator moves on a geodesic; any further speed-up would lead to mass shedding (the “Kepler” limit). Uniform rotation can therefore increase the maximum allowed mass by at most about 20 % for very stiff equations of state (Cook, Shapiro & Teukolsky, 1992, hereafter CST1, and 1994, hereafter CST2), which is not sufficient to stabilize remnants of binary neutron star merger. Uniformly rotating equilibrium configurations with rest masses exceeding the maximum rest mass of nonrotating stars constructed with the same equation of state are referred to as “supramassive” stars (CST1, CST2).

Differential rotation, however, can be much more efficient in increasing the maximum allowed mass. In differentially rotating stars, the core may rotate faster than the envelope, so that the core can be supported by rapid rotation without the equator having to exceed the Kepler limit. This effect was demonstrated in Newtonian gravitation by Ostriker, Bodenheimer & Lynden-Bell (1966) for white dwarfs, and in general relativity by Baumgarte, Shapiro & Shibata (2000, hereafter BSS) for $`n=1`$ polytropes. BSS also showed by way of example that stars with about 60 % more mass than the maximum allowed mass of the corresponding non-rotating star can be dynamically stable against both radial and nonaxisymmetric modes. BSS refer to differentially rotating equilibrium configurations with rest masses exceeding the maximum rest mass of a uniformly rotating star as “hypermassive” stars.

In this paper we extend the findings of BSS for $`n=1`$ polytropes in two ways. We survey different polytropic indices, and study the effect of both differential rotation and stiffness of the equation of state on the maximum allowed mass. We also introduce a very simple model calculation that illustrates these effects qualitatively, and that moreover gives surprisingly accurate results for soft equations of state and moderate degrees of differential rotation.

This paper is organized as follows. In Section 2 we present qualitative considerations, leading to a simple estimate of the maximum mass of differentially rotating neutron stars. In Section 3 we construct numerical models of fully relativistic, differentially rotating neutron stars for different polytropic indices. We briefly discuss our findings in Section 4. We also include Appendix A with tables of our numerical results.

## 2 Qualitative Considerations

When magnetic fields and relativistic effects can be neglected, rotating stars in equilibrium satisfy the Newtonian virial theorem

$$W+2T+3\mathrm{\Pi }=W\left(12\frac{T}{|W|}\right)+3\mathrm{\Pi }=0$$
(1)
(see, e.g., Shapiro & Teukolsky 1983). Here the rotational kinetic energy $`T`$ scales as

$$T\frac{J^2}{MR_e^2},$$
(2)
the potential energy as

$$W\frac{M^2}{R_e}M^{5/3}\rho ^{1/3},$$
(3)
and the internal energy, computed from a volume integral of the pressure, as

$$\mathrm{\Pi }\rho ^{1/n}M.$$
(4)
where $`J`$ is the angular momentum, $`M`$ the mass, $`R_e`$ the equatorial radius and $`\rho `$ the mass density. We have also assumed a polytropic equation of state

$$P=K\rho ^{1+1/n},$$
(5)
where $`K`$ is a constant and $`n`$ the polytropic index. Here and throughout we set the gravitation constant $`G=1`$.

Inserting the above expressions into the virial equation (1), we find

$$\beta M^{5/3}\rho ^{1/3}\left(12\frac{T}{|W|}\right)+\gamma \rho ^{1/n}M=0,$$
(6)
where $`\beta `$ and $`\gamma `$ are the appropriate coefficients. In general, $`\beta `$ depends on the eccentricity of the star, but restricting our analysis to small values of $`T/|W|`$, and hence to nearly spherical stars, we assume that both $`\beta `$ and $`\gamma `$ are constant.

For non-rotating stars ($`T=0`$) we solve equation (6) to find

$$M=\left(\frac{\gamma }{\beta }\right)^{3/2}\rho ^{(3n)/(2n)},$$
(7)
while for rotating stars we obtain

$$M_{\mathrm{rot}}=M\left(12\frac{T}{|W|}\right)^{3/2}.$$
(8)
For small values of $`T/|W|`$ the right hand side can be expanded, and we then find for the mass increase $`\delta MM_{\mathrm{rot}}M`$

$$\frac{\delta M}{M}=3\frac{T}{|W|}$$
(9)
(cf. Shapiro & Teukolsky 1983, eqn 7.4.40). This expression determines the fractional mass increase as a function of $`T/|W|`$ for a constant value of the density $`\rho `$. We will use this result to estimate the increase in the maximum allowed mass of a neutron star, even though typically rotating stars assume their maximum masses at slightly different densities than the corresponding non-rotating stars.

To evaluate $`T`$ and $`|W|`$ we further simplify the problem by assuming that the star’s density profile is a step function with a constant density $`\rho _c`$ (equal to the original central density) inside a spherical core of radius $`R_c`$, and zero outside. From requiring that the mass of this model star,

$$M=\frac{4\pi }{3}\rho _cR_c^3,$$
(10)
be equal to the original mass

$$M=\frac{4\pi }{3}\overline{\rho }R_e^3,$$
(11)
where $`\overline{\rho }`$ is the average density, we find the following relation between the central condensation and the ratio of the radii

$$\frac{R_e}{R_c}=\left(\frac{\rho _c}{\overline{\rho }}\right)^{1/3}.$$
(12)
Further assuming that the core is uniformly rotating with the central angular velocity $`\mathrm{\Omega }_c`$, we find

$$T=\frac{1}{2}I\mathrm{\Omega }_c^2=\frac{MR_c^2}{5}\mathrm{\Omega }_c^2.$$
(13)
Inserting this relation together with the potential energy

$$W=\frac{3}{5}\frac{M^2}{R_c}.$$
(14)
into (9) yields

$$\frac{\delta M}{M}\frac{R_c^3\mathrm{\Omega }_c^2}{M}.$$
(15)
To find the increase in the maximum allowed mass, it is useful to assume that the star rotates at the mass-shedding limit

$$\mathrm{\Omega }_e^2=\frac{M}{R_e^3},$$
(16)
where $`\mathrm{\Omega }_e`$ is the equatorial angular velocity. This equation can now be used to eliminate the mass $`M`$ in equation (15), which, together with (12), yields

$$\frac{\delta M}{M}\frac{\overline{\rho }}{\rho _c}\left(\frac{\mathrm{\Omega }_c}{\mathrm{\Omega }_e}\right)^2.$$
(17)
Equation (17) provides a very simple estimate for the increase of the maximum allowed mass. It depends only on the central condensation of the non-rotating star, which is a function of the stiffness of the equation of state, and the ratio of the angular velocities at the center and equator, which is a function of the the degree of differential rotation. For uniformly rotating stars, the maximum mass increase is estimated to be simply the inverse of the central condensation. In Table 1 we compare this estimate with the numerical findings of CST2 and find remarkably good agreement for soft equations of state. Table 1 also illustrates an ambiguity; in Newtonian gravity, the central condensation is uniquely determined by the polytropic index, but in general relativity the central condensation of a star depends on the central density. We therefore compute a “relativistic” central condensation $`\rho _c/\overline{\rho }|_{\mathrm{TOV}}`$ from the central energy density $`\rho _c`$ and an average density defined as

$$\overline{\rho }=\frac{3M}{4\pi R^3}$$
(18)
of the non-rotating maximum mass model, where $`M`$ is the total mass-energy of the star, and $`R`$ is the circumferential radius. We find that this value yields better agreement with the numerical values of the maximum mass increase than adopting the Newtonian central condensation.

The ratio $`T/|W|`$ provides a useful criterion for the onset of secular ($`T/|W|0.14`$) or dynamical ($`T/|W|0.27`$) non-axisymmetric instabilities. Inserting these limits into equation (8) shows that mass increases of secularly stable stars are limited by $`\delta M/M1.63`$ while mass increases of dynamically stable stars are limited by $`\delta M/M3.2`$. These values agree quite well with the respective limits of 1.70 and 3.51 found by Shapiro & Teukolsky (1983; equations (7.4.41) and (7.4.42)), who also take stellar deformations into account.

## 3 Numerical Results

We use a modified version of the numerical code of CST1 and CST2 to construct models of differentially rotating neutron stars. The code is based on similar algorithms developed by Hachisu (1986) and Komatsu, Eriguchi and Hachisu (1989), and we refer to CST1 for details. We adopt a polytropic equation of state

$$P=K\rho _0^{1+1/n},$$
(19)
where $`\rho _0`$ is the rest-mass density and where equation (19) reduces to equation (5) in the Newtonian limit. We take the polytropic constant $`K`$ to be unity without loss of generality. Since $`K^{n/2}`$ has units of length, all solutions scale according to $`\overline{M}=K^{n/2}M`$, $`\overline{\rho }_0=K^n\rho _0`$, etc., where the barred quantities are dimensionless quantities corresponding to $`K=1`$, and the unbarred quantities are physical quantities (compare CST1).

Constructing differentially rotating neutron star models requires choosing a rotating law $`F(\mathrm{\Omega })=u^tu_\varphi `$, where $`u^t`$ and $`u_\varphi `$ are components of the four-velocity $`u^\alpha `$ and $`\mathrm{\Omega }`$ is the angular velocity. We follow CST1 and assume the rotation law $`F(\mathrm{\Omega })=A^2(\mathrm{\Omega }_c\mathrm{\Omega })`$, where the parameter $`A`$ has units of length. Expressing $`u^t`$ and $`u_\varphi `$ in terms $`\mathrm{\Omega }`$ and metric potentials yields equation (42) in CST1, or, in the Newtonian limit,

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_c}{1+\widehat{A}^2\widehat{r}^2\mathrm{sin}^2\theta }.$$
(20)
Here we have rescaled $`A`$ and $`r`$ in terms of the equatorial radius $`R_e`$: $`\widehat{A}A/R_e`$ and $`\widehat{r}r/R_e`$. The parameter $`\widehat{A}`$ is a measure of the degree of differential rotation and determines the length scale over which $`\mathrm{\Omega }`$ changes. Since uniform rotation is recovered in the limit $`\widehat{A}\mathrm{}`$, it is convenient to parametrize sequences by $`\widehat{A}^1`$. In the Newtonian limit, the ratio $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ that appears in the estimate (17) is related to $`\widehat{A}^1`$ by $`\mathrm{\Omega }_c/\mathrm{\Omega }_e=1+\widehat{A}^2`$, but for relativistic configurations this relation holds only approximately.

We adopt this particular rotation law for convenience and for easy comparison with many other authors who have assumed the same law. We also compared with the remnants’ angular momentum distribution in the fully relativistic dynamical merger simulations of Shibata and Uryu (2000) and to the post-Newtonian simulations of Faber, Rasio and Manor (2001). We have found that their numerical results can be fit reasonably well by our adopted differential rotation law.

We modify the numerical algorithm of CST1 by fixing the maximum interior density instead of the central density for each model. This change allows us to construct higher mass models in some cases, since the central density does not always coincide with the maximum density and hence may not specify a model uniquely.

For a given a value of $`n`$ and $`\widehat{A}`$, we construct a sequence of models for each value of the maximum density by starting with a static, spherically symmetric star and then decreasing the ratio of the polar to equatorial radius, $`R_{pe}=R_p/R_e`$, in decrements of 0.025. This sequence ends when we reach mass shedding (for large values of $`\widehat{A}`$), or when the code fails to converge (indicating the termination of equilibrium solutions) or when $`R_{pe}=0`$ (beyond which the star would become a toroid). For each one of these sequences the maximum achieved mass is recorded. We repeat this procedure for different values of the maximum density, covering about a decade below the central density of the non-rotating maximum mass model, which yields the maximum allowed mass for the chosen values of $`n`$ and $`\widehat{A}`$. Our numerical results are tabulated in Appendix A. Our maximum mass increases are lower limits in the sense that even higher mass models may exist, but that we have not been able to construct them numerically.

## 4 Discussion and Summary

Our numerical results are tabulated in Tables 2 to 9 in Appendix A. We also compare the increases in the maximum allowed mass with the estimate (17), and find surprisingly good agreement for soft equations of state and moderate degrees of differential rotation.

In particular, we find that the fractional maximum rest mass increase $`\delta M/M`$ for uniformly rotating stars is well approximated by the inverse of the central concentration (see also Table 1). For moderate degrees of differential rotation, $`\delta M/M`$ increases approximately with the square of the ratio between the central and equatorial angular velocity $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$, in accord with equation (17).

For all equations of state we find that $`\delta M/M`$ increases with $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ only up to a moderate value of $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$, and starts to decrease again for larger values (at least with our code and algorithm we do not find monotonically increasing mass configurations). For stiff equations of state this turn-around occurs for smaller values of $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ than for soft equations of state. For moderate degrees of differential rotation $`\mathrm{\Omega }_c/\mathrm{\Omega }_e2`$, a given value of $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ will lead to a larger increase in the maximum allowed mass for a stiffer equation of state, as expected from the estimate (17).

We find the largest maximum mass increases for moderately stiff equations of state. Some of these configurations exceed the maximum allowed mass of the corresponding non-rotating star by more than a factor of two. These configurations typically have large values of $`T/|W|0.27`$, indicating that such stars may by dynamically unstable against bar formation (but see Shibata, Karino & Eriguchi 2002, who found mild bar mode instabilities at very small values of $`T/|W|`$ for extreme degress of differential rotation). They are also “toriodal”, i.e. assume their maximum density on a torus around the center of the star, which may indicate an $`m=1`$ instability at even smaller values of $`T/|W|`$ (Centrella, New, Lowe & Brown 2001). However, even restricting attention to those configurations that are not toriodal and have $`T/|W|<0.27`$, we find configurations with masses larger than the maximum mass of the corresponding non-rotating star by over 60 %. BSS demonstrated that at least some of these models are dynamically stable. Shibata & Uryu (2000) demonstrated that binary mergers may result in similarly stable hypermassive stars, when the progenitor masses are not too close to the maximum mass.

To summarize, we find that differential rotation is very effective in increasing the maximum allowed mass, especially for moderately stiff equations of state. The effect is probably large enough to stabilize the remnants of binary neutron star merger, which are likely to be differentially rotating. Binary neutron star coalescence may therefore lead to secularly stable, hypermassive neutron stars. As discussed in BSS (see also Shapiro 2000), magnetic braking is likely to bring such differentially rotating stars into uniform rotation, which reduces the maximum allowed mass and induces a delayed collapse to a Kerr black hole.

This work was supported in part by NSF Grants PHY-0090310 and PHY-0205155 and NASA Grant NAG5-10781 at the University of Illinois at Urbana-Champaign and NSF Grant PHY-0139907 at Bowdoin College. NDL gratefully acknowledges support through an Undergraduate Research Assistantship from the Department of Physics and Astronomy at Bowdoin College.

## Appendix A Numerical Results for Maximum Masses

We list below in Tables 2 to 9 values for the maximum rest mass increase for uniformly and differentially rotating polytropes. For each polytropic index $`n`$ we tabulate the differential rotation parameter $`\widehat{A}^1`$, the ratio of the central and equatorial angular velocity $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ (which reduces to (20) in the Newtonian limit, i.e. for soft equations of state), the numerically determined fractional rest mass increase $`(\delta M/M)_{\mathrm{num}}`$, the ratio of the (relativistic) rotational kinetic energy and the gravitational binding energy $`T/|W|`$, the ratio between polar and equatorial radius $`R_p/R_e`$, the maximum density $`\overline{\rho }_{\mathrm{max}}`$, and the estimate $`(\delta M/M)_{\mathrm{est}}`$ according to equation (17). In these estimates we used the numerically determined ratios $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ and the central condensations according to (18). For $`n1.25`$ some of the maximum mass configurations are toriodal, i.e. assume the maximum density on a toroid about the center. For these polytropic indices we also include the ratio $`\rho _c/\rho _{\mathrm{max}}`$.

All models are computed with the code of CST1 and CST2, using 64 zones both in the radial and angular direction, and truncating the Legendre polynomial expansion at $`\mathrm{}=16`$ (see CST1 for details of the numerical implementation). The accuracy of individual stellar models can be tested, for example, by computing a relativistic Virial theorem (Gourgoulhon & Bonazzola 1994; Cook, Shapiro & Teukolsky 1996; see also Nozawa et al. 1998 for a comparison of several different computational methods). In our analysis, however, the error in the maximum mass and related quantities is dominated by the finite step size in the sequences over $`R_p/R_e`$ and $`\overline{\rho }_{\mathrm{max}}`$, which result in errors typically in the order of a few percent. For soft equations of state the mass as a function of central density is a very slowly varying function, making it quite difficult to determine the central density of the maximum mass configuration very accurately. The error in the $`R_p/R_e`$ is determined by our stepsize of 0.025. We finally note that highly toroidal configurations depend very sensitively on the input parameters, so that those mass increases should only be taken as estimates.


format: text

text:
# Nonequilibrium transitions in complex networks: a model of social interaction

## I Introduction

Lattice models are a powerful basic instrument in the study of phase transitions in equilibrium Statistical Mechanics, as well as in non-equilibrium systems Marro98 . Traditionally, equilibrium phase transitions have been studied in regular lattices, with the critical temperature being a non-universal quantity that depends on the particular lattice under consideration, while critical exponents and some amplitude ratios are universal quantities depending only on spatial dimension and some symmetries of the order parameter. The detailed structure of the regular network connections is, in most cases, irrelevant in the sense of the renormalization group. However, recent research in the structure and topology of complex networks rmp ; aip  has shown that social interactions and, more generally, biological and technological networks, are far from being regular as well as being also far from a random network or from a mean-field network linking all-to-all. This has triggered the study of standard models of Statistical Mechanics in these complex networks. In particular, recent studies of the Ising model in the so–called small-world Watts98  and the scale-free Barabasi99  networks have shown that the behavior of the model differs from that observed in a regular network.

A small world network Watts98  is generated by rewiring with a probability $`p`$ the links of a regular lattice by long-distance random links. The presence of a small fraction of “short cuts” connecting otherwise distant points, drastically reduces the average shortest distance between any pair of nodes in network, keeping the clustering high. The small-world networks generated by rewiring links have degree distributions with exponential tails. In contrast, scale-free networks Barabasi99  are characterized by a fat-tailed (power law) degree distribution. They are usually modelled by growing networks and preferential attachment of links.

The Ising model in small-world topologies shows a change of behavior from the regular case to the mean field characteristics. In Ref. Barrat00  it is shown analytically that for a small world lattice, obtained from rewiring with probability $`p`$ the links of a 1-d ring lattice with $`2k`$ nearest neighbors interactions, the crossover temperature to the mean field critical behavior varies for $`p<<1`$ as $`T_{co}(p)k(k+1)/\mathrm{ln}(p)`$ whereas the critical temperature scales as $`T_c(p)2k/\mathrm{ln}(p)`$, so that a ferromagnetic ordered phase exists for any finite value of $`p`$. The crossover to mean field behavior in small world ring lattices has been further discussed in Refs. Gitterman00 ; Hong02 , whereas numerical results in 2-d and 3-d lattices have been reported in Herrero02 . Interestingly, if directed links are considered, not only the critical temperature changes but the nature of the transition also switches from second to first order Sanchez02 .

A much different behavior is observed in scale-free networks. This can be related to the influence of the presence of so-called hubs, i.e.  units whose degree is much larger than average. This is well illustrated by the behavior of the Ising model in scale-free networks with degree distribution $`P(k)k^\gamma `$, $`\gamma >1`$. The results of Refs. Aleksiejuk02  and Goltsev02  show that equilibrium systems exhibit a phase diagram that is qualitatively different from the mean-field case. In particular, the Ising model in a random scale-free networks shows an infinite critical temperature in the thermodynamic limit of an infinite number of nodes. In fact it has been developed analytically a mean field theory connecting the exponent of the degree distribution and the critical behavior of the Ising model Dorogovtsev02 ; Leone02 ; Bianconi02 .

In this paper we address the question of the role played by the topology of complex networks in non-equilibrium transitions of models in which there is interaction between the variables associated with the nodes connected by links in the network. This is a natural next step beyond the analysis of equilibrium, Ising–type models in these complex networks. Simple non-equilibrium models closely related to percolation have been already considered Kuperman01 ; Pastor01 ; Eguiluz02 . Here, and given the social motivation and relevance of these complex networks, we have chosen to analyze the model proposed by Axelrod for the dissemination of culture Axelrod97 . The spreading process in this model cannot be reduced to a percolation process. The model rather describes a competition between dominance and spatial coexistence of different states in a non-equilibrium dynamics of coupled Potts-type models. The model was originally considered by Axelrod in a square lattice. The Statistical Mechanics analysis of the model in this regular two-dimensional network identifies a nonequilibrium order-disorder phase transition Castellano00 . However, it is interesting to notice that, in his original paper, Axelrod already discussed the relevance of the topology, speculating that “With random long-distance interactions, the heterogeneity sustained by local interaction cannot be sustained.” In particular we consider here this question.

In the next Section we introduce the original model proposed in Ref. Axelrod97  and summarize briefly the main results in regular 2-d networks. The model in small-world and scale-free networks is presented in Section III and IV respectively. The non-equilibrium transition is shown to disappear in the thermodynamic limit of a scale free network. We then consider in Section V a structured scale free network Klemm02  which shares characteristics of small-world and scale-free networks. A non-equilibrium transition is shown to persist for large systems in this network. Our conclusions are summarized in Section VI.

## II The model

The model we study is defined Axelrod97  by considering $`N`$ agents as the sites of a network. The state of agent $`i`$ is a vector of $`F`$ components (cultural features) $`(\sigma _{i1},\sigma _{i2},\mathrm{},\sigma _{iF})`$. Each $`\sigma _{if}`$ can take any of the $`q`$ integer values (cultural traits) $`1,\mathrm{},q`$, initially assigned independently and with equal probability $`1/q`$. The time-discrete dynamics is defined as iterating the following steps:

1. Select at random a pair of sites of the network connected by a bond $`(i,j)`$.
2. Calculate the overlap (number of shared features) $`l(i,j)=_{f=1}^F\delta _{\sigma _{if},\sigma _{jf}}`$.
3. If $`0<l(i,j)<F`$, the bond is said to be active and sites $`i`$ and $`j`$ interact with probability $`l(i,j)/F`$. In case of interaction, choose $`g`$ randomly such that $`\sigma _{ig}\sigma _{jg}`$ and set $`\sigma _{ig}=\sigma _{jg}`$.

In any finite network the dynamics settles into an absorbing state, characterized by the absence of active bonds. Obviously all the $`q^F`$ completely homogeneous configurations are absorbing. Homogeneous means here that all the sites have the same value of the cultural trait for each cultural feature. Inhomogeneous states consisting of two or more homogeneous domains interconnected by bonds with zero overlap are absorbing as well. A domain is a set of contiguous sites with identical cultural traits. We define an order parameter in this system Castellano00 ; Klemm02b  as the relative size of the largest homogeneous domain $`S_{max}/N`$, being $`N`$ the number of sites in the network.

Previous results have been obtained in square lattices with nearest neighbor interaction. A variation of the model with initial distribution of traits according to a Poisson rather than a uniform distribution shows a non-equilibrium order-disorder phase transition where the number of traits $`q`$ plays the role of a control parameter Castellano00 . The system reaches ordered absorbing states for $`q<q_c`$ ($`S_{max}=𝒪(N)`$) and disordered states for $`q>q_c`$ ($`S_{max}N`$). The same type of phase transition occurs in the original model with a uniform initial distribution of traits Klemm02b .

When comparing the effect of complex networks in this phase transition with the equilibrium Ising transition one should notice several conceptual differences. First, this is a sort of zero-temperature transition in which ordered or disordered states exists with no reference to thermal fluctuations. In fact, the effect of small noise in this system is essential, revealing the presence of metastable states and changing the phase diagram in a nontrivial way Klemm02b . A second related point is that the control parameter of the transition $`q`$, is here not a collective property of the system as temperature, but rather an ingredient of the definition of the system itself. In a way, the transition occurs going from one system to another as $`q`$ is changed. Finally, and in reference to critical properties and exponents, we note that the transition (except for F=2) is of first-order type: In figure 1 we plot the final values for the order parameter, obtained for $`32`$ different realizations of the dynamics. Notice that for $`q<50`$ all the systems end up in a homogenous state that basically fills up the entire lattice ($`S_{max}/N1`$) whereas for $`q>60`$ the maximum homogenous regions obtained are very small. This is the order–disorder phase transition discussed before. For $`50<q<60`$ we observe bistability in the sense that the system settles around any of two mean values for the order parameter. This bistability, which is usually associated with first order phase transitions, is clearly made explicit in the corresponding histogram shown in figure 2 where the two preferred values appear as maxima of the histogram. The transition point corresponds to $`q=q_c`$ for which these two values are equally probable.

## III Small world network

Social networks are far from being regular or completely random. However they also share some features with them. On the one hand, social networks are known to be small Milgram67 , i.e. any pair of nodes in the network can be connected following a number of links much smaller than the size of the network. This is also observed in random networks, where the average shortest distance between pair of nodes (the so called path length $`\mathrm{}`$) increases logarithmically with the size of the network $`\mathrm{}\mathrm{ln}N`$, while in regular lattice in $`d`$-dimensions $`\mathrm{}N^{1/d}`$. On the other hand, social networks are also known to form cliques Wasserman94 , i.e. groups of nodes highly connected between them. “Cliquishness” can be characterized by the so-called clustering coefficient $`C`$, which is defined as the relative number of closed triangles in the network. Regular lattices can show large clustering while in random networks $`CN^1`$. High clustering and short path length define a small world network.

The first model encompassing the small world effect was introduced in Ref. Watts98  proposing an algorithm that interpolates between a random and a regular lattice. First one generates a two-dimensional regular lattice with bonds between nearest neighbors and open boundary conditions. Then for each bond $`(ij)`$, with probability $`p`$ detach the bond from node $`j`$ and attach it to a node $`l`$ instead. Node $`l`$ is chosen at random with the restriction that duplicate and self-connections are excluded. The parameter $`p`$ interpolates between the original regular lattice ($`p=0`$, no rewiring) and a network very similar to a random graph ($`p=1`$). Thus in the limiting case $`p=0`$ we have a network with high clustering but also large path length; in the limit $`p=1`$ we have networks with small path length but also small clustering. For intermediate values of $`p`$ the algorithm generates networks with high clustering and small path length.

We now study the behavior of Axelrod’s model in dependence of $`p`$. A small world network is used from the beginning of each simulation run. Figure 3 shows the dependence of the order parameter on $`q`$, for three different values of $`p`$. For any fixed value $`p>0`$ we find a nonequilibrium phase transition which becomes sharp and well defined as the system size increases. There is a critical value $`q_c`$ of the control parameter which separates the ordered and the disordered state, just as in regular lattices. However, $`q_c`$ increases with the amount of spatial disorder. This is clearly shown in the $`(p,q)`$-phase diagram, Fig. 4. The filled area above the $`(p,q_c(p))`$-curve represents the disordered states, the area below the curve represents the ordered states. Consequently, for values $`q<q_c(p)`$ the outcome of the dynamics is always complete order, whereas for $`q>q_c(p)`$ only disordered frozen states are encountered. The density $`p`$ of rewired bonds determines the nature of these frozen states, but for $`q<q_c(p=1)`$ the system orders by increasing $`p`$, that is, the number of long distance links. We find a dependence $`q_c(p)q_c(p=0)p^\alpha `$ with $`\alpha =0.4`$ obtained from a best fit. This result is displayed in the inset of Fig. 4. Therefore, we find the same qualitative result that for the equilibrium Ising model, in the sense that the small world connectivity favors ordered states.

The robustness of the phase diagram is shown by performing a different dynamical scenario. First, a run of the dynamics in a regular lattice is performed. Only after an absorbing configuration has been reached the lattice is rewired according to the above rewiring procedure with the parameter $`p`$. After the rewiring, the configuration is not necessarily frozen because the rewiring can introduce active bonds connecting compatible cultures that have been disconnected before. Starting the dynamics again, the system may relax to a different absorbing configuration, which in general is more ordered than the configuration reached before the rewiring. After this second phase of relaxation, the order parameter is measured in the absorbing state. We find that the results of this alternative scenario (see Fig. 4) are in good agreement with the ones of the above original scenario, starting with a small world network in the initial condition.

## IV Scale-free networks

One important ingredient missing in the small world networks considered so far is that the degree distribution does not show a fat tail. Although it is not clear whether social networks present a power law distribution of degree, the evidence indicates that they are ubiquitous in biological and artificial networks Amaral00 ; Albert02 . Scale-free networks are characterized by a power law tail in the degree distribution of the form $`P(k)k^\gamma `$ where the exponent $`\gamma `$ lies in the range between 2 and 3. Two ingredients have been shown to be sufficient to generate such feature: growing number of nodes and preferential attachment of links.

The well established Barabási–Albert model is based in these two mechanisms Barabasi99 . At each time step a new node is added to the network and attaches $`m`$ links to an existing node with degree $`k`$ with probability $`\mathrm{\Pi }(k)k`$. This algorithm generates networks whose degree distribution follows $`P(k)=2m^2k^3`$, the path length increases as $`\mathrm{}\mathrm{ln}N`$, and the clustering decreases as $`C(\mathrm{ln}N)^2/N`$. We have studied the dynamics of Axelrod’s model for the diffusion of culture in scale-free networks following this algorithm. In Fig. 5 we show our results for the order parameter for different system sizes. For a given size $`N`$ we find a transition at $`q_c(N)`$. We can define the critical value $`q_c(N)`$ as the value where the standard deviation of the distribution of $`S_{\mathrm{max}}/N`$ reaches the maximum value. We find that $`q_c(N)N^{0.39}`$. Using this result we observe data collapse with a rescaling $`qN^\beta `$, see Fig. 6. The best result is obtained for $`\beta =0.39`$ in excellent agreement with the scaling obtained previously. This indicates that in the thermodynamic limit the transition disappears and the ordered monocultural state establishes in the system. This behavior is similar to the Ising model in regular and scale-free networks: While in a two-dimensional lattice the Ising model displays a phase transition at a finite critical temperature, in random scale-free networks an effective transition is observed for finite systems where the effective critical temperature diverges logarithmically with system size. This can be explained by the role of the hubs (nodes with a large number of links) in these networks. They are highly instrumental in establishing ferromagnetic order in the system. The same prominent role is played by the hubs in the case of the dissemination of culture. The hubs help the spreading of cultural traits as can be inferred from the observed dependence with system size. Note, however that the effective transition of Axelrod s model for a finite system in a scale free network displays the characteristics of a first order transition: We show in Figs. 5 and 7 the same type of behavior observed in Figs. 1 and 2 for the regular network. For a range of values of $`q`$ around $`q_c`$ a realization ends either in an ordered monocultural state or in a disordered frozen configuration, the two preferred values of the order parameter clearly seen in the histogram of Fig. 2.

## V Structured scale-free networks

The scale-free networks considered in the previous section underestimate the clustering observed in real networks Ravasz02 . A question that merits to be addressed is if this low clustering coefficient is responsible for the absence of the phase transition in the thermodynamic limit. In order to reproduce a high clustering along with a scale-free distribution of the degree, we employ the networks generated by the algorithm proposed in Ref. Klemm02 : Again at each time step we add a new node to the network. The node is attached to the $`m`$ active nodes in the network. The new node becomes active and one of the $`m+1`$ active nodes is deactivated with probability $`\mathrm{\Pi }(k)k^1`$. Starting from $`m`$ fully inter-connected active nodes, this algorithm generates scale-free networks with a clustering coefficient $`C5/6`$ independent of system size. In the following we call these networks structured scale-free networks because of the large clustering coefficient, the strong negative correlation between degrees of connected nodes Eguiluz02 , and the modular structure Ravasz02 . These properties are not found in the random scale-free networks of the previous section.

As shown in Fig. 9, in the structured scale-free networks the model displays a behavior different to what we observed in random scale-free networks in the previous section. For $`q10`$ the system settles into an ordered state. For increasing values of $`q`$, the order parameter undergoes a decay whose slope grows with system size. This suggests a phase transition at $`q_c10`$, in contrast to the absence of a transition point found for the randomly wired scale-free networks in the thermodynamic limit. As on large scales the structured scale-free networks have one-dimensional topology Klemm02a  it seems natural that this transition at $`q_cF`$ coincides with the behavior of the model found in one-dimensional regular lattices Klemm02c .

At difference with the regular lattices, in the structured scale-free networks for $`q>q_c`$ the order parameter does not tend to zero. It reaches a finite plateau value, indicating partial ordering of the system. Only for values $`qq_c`$ the order parameter drops below the plateau value and tends to zero. This behavior may be understood by relating the size $`S_{max}`$ of the largest cultural cluster with the largest degree $`k_{max}`$ present in the network, as shown in Fig. 10. In the intermediate range $`050<q<200`$, where the plateau of the order parameter is observed, we find $`S_{max}k_{max}`$ for almost all realizations. This suggests that the largest hub, the node with the largest degree, and its neighbors order such that they form the largest cluster in the absorbing state. As $`q`$ is reduced and its value approaches $`q_c`$ from above, the ordering goes beyond the largest hub and an increasing part of the network forms the largest cluster. On the other hand, for large values $`q>200`$, the neighborhood of the largest hub does no longer reach complete ordering and $`S_{max}<k_{max}`$. The $`q`$ value for the onset of the decay of $`S_{max}`$ below $`k_{max}`$ is expected to be dependent on system size: as Fig. 9 shows, for increasing system size, the plateau in the order parameter extends to larger values of $`q`$. With inrceasing system size, the value of the plateau is expected to decrease as $`k_{max}/N=N^{\beta 1}`$ with $`\beta =(\gamma 1)^1`$, where $`\gamma `$ is the exponent of the degree distribution. These results suggest that in the limit $`N\mathrm{}`$, the dynamics of the social interaction model in structured networks experiences a transition similar the one observed in a one-dimensional lattice.

## VI Conclusions

We have found that the nonequilibrium transition between order and disorder that exists in a regular d=2 network for Axelrod s model of cultural influence Axelrod97  is modified by underlying complex networks with similar qualitative features that an equilibrium thermal Ising-type transition. We have shown that the transition pertains also in the presence of random long-distance connections: with increasing density of long-distance connections in a small world network, the critical point $`q_c(p)`$ increases. Therefore the small world connectivity favors cultural globalization as described by the ordered state. The value of $`q_c`$ reaches a maximum for the random network obtained with a $`p=1`$ probability of rewiring in the small word network construction. A transition from disorder to order is obtained increasing $`p`$ for a fixed value of the control parameter $`q<q_c(p=1)`$. We have also found that, for a fixed finite system size, the scale free connectivity is more efficient than the limiting random connectivity of the small world network in promoting the ordered state of cultural globalization. However, there is a system size scaling in the transition observed for a free scale network, so that the transition disappears in the thermodynamic limit: In the presence of scale-free interactions the order state prevails due to the presence of hubs. The consideration of structured scale free-networks restores the order-disorder transitions in spite of the hubs, but the value of the order parameter for the disordered state reveals the existence of ordered clusters.


text:
# Orbital Migration and the Brown Dwarf Desert

## 1. Introduction

Radial velocity surveys show that brown dwarfs are rarely found as close (within a few au) binary companions to Solar-type stars. This article makes the case that the existence of this brown dwarf desert may have nothing to do with the formation of brown dwarfs. We argue instead that a brown dwarf desert is an expected consequence of the orbital migration of brown dwarfs within an evolving protostellar disk (Armitage & Bonnell 2002). We make the basic assumption that brown dwarfs form in the same manner as stars – during gravitational collapse – and show that migration subsequently creates a desert by driving initially close star - brown dwarf binaries into mergers.

## 2. Why migration matters for brown dwarfs

A viscous accretion disk around a single star evolves such that the outward drift of a vanishingly small fraction of the mass absorbs all of the angular momentum, allowing the bulk of the disk to be accreted (e.g. Pringle 1981). The same reasoning applies to circumbinary disks (Pringle 1991). Given enough time, an arbitrarily small mass of gas, trapped outside the binary orbit, will soak up the binary orbital angular momentum, and drive the binary toward merger.

Protoplanetary disks, however, have limited lifetimes of around 1-10 Myr (Strom et al. 1989; Haisch, Lada & Lada 2001). For rapid migration to occur, we require (roughly speaking) that the mass of gas in the disk be at least comparable to the mass of the embedded substellar object111This condition can be easily understood. If the angular momentum of the substellar object greatly exceeds that of the disk, it will take many viscous times for the disk to transport the ‘extra’ angular momentum to large radii, and migration will be very slow.. This requirement is almost certainly met for brown dwarfs that are binary companions to Solar-type stars, since both observations (Osterloh & Beckwith 1995) and theoretical arguments (Lin & Pringle 1990) suggest that the disk mass at early times comfortably exceeds typical brown dwarf masses. We therefore expect migration to be significant for brown dwarfs orbiting Solar mass stars, though not for brown dwarfs embedded in the weaker disks surrounding very low mass stars.

## 3. Forming a brown dwarf desert

To make quantitative predictions for the impact of migration on brown dwarfs, we have used a one dimensional model that was developed to study planetary migration (Armitage et al. 2002; see also Trilling et al. 1998; Trilling, Lunine & Benz 2002). We solve for the coupled evolution of a viscous protoplanetary disk (using a variant of the disk model of Clarke, Gendrin & Sotomayor 2001) and an embedded brown dwarf. We assume that the brown dwarf orbit remains circular, not because we believe this to be true (Papaloizou, Nelson & Masset 2001; Goldreich & Sari 2002), but because of the difficulty of doing better. We also assume that the brown dwarf does not accrete further mass from the disk during migration. For circular orbits, this is likely to be a reasonable approximation for substellar objects with masses above a few Jupiter masses (Lubow, Seibert & Artymowicz 1999).

Figure 1 shows the extent of the resulting migration for $`0.04M_{}`$ brown dwarfs embedded in disks whose initial mass is $`0.1M_{}`$. Brown dwarfs with initial orbital radii smaller than about 5 au are all swept in to small radii and probably merge with the star (conceivably, as has been suggested for extrasolar planets, some survive at very small radii as the brown dwarf equivalent of hot Jupiters). At larger radii, where the disk is expanding, outward migration occurs and pushes brown dwarfs out to radii of $`10^2`$ au. Between these extremes, brown dwarfs formed within a small window of initial orbital radii spread out to populate all short period orbits. This dilution of $`\mathrm{d}N/\mathrm{d}\mathrm{log}a`$ at small radii leads to a brown dwarf desert within this model.

We do not expect the desert – if it results from migration – to be completely devoid of brown dwarfs. A handful of brown dwarfs should end up in short period orbits having migrated inward from initially larger disk radii. Depending upon the assumed initial (pre-migration) distribution of brown dwarf orbital radii, we have estimated that the frequency of brown dwarfs within 4 au should be depleted by a factor of 5-10 (Armitage & Bonnell 2002).

Why is there a brown dwarf desert at the same radii where there is an abundance of massive extrasolar planets (Marcy & Butler 2000)? The arguments given above suggest – correctly – that the lower masses of extrasolar planets should make them even more vulnerable to migration than brown dwarfs. Massive planets, however, probably do not form via the same mechanism as brown dwarfs (though for contrary opinions see Boss 1998; Armitage & Hansen 1999). Conventionally, they form instead via core accretion, a slow process that may take several Myr, especially in the outer parts of the disk. As shown in Figure 1, this allows planets to escape destruction via mergers by forming at a later epoch, when the protoplanetary disk is of lower mass and unable to force the planets into mergers.

## 4. Observational tests

We have argued that one should not conclude, based solely on the existence of the brown dwarf desert, that brown dwarfs cannot form as close binary companions to Solar-type stars. Migration at a later epoch can accomplish the same end result, and makes several additional predictions.

* Migration occurs if the binary mass ratio $`q=M_{\mathrm{BD}}/M_{}`$ is comparable to or less than the disk mass ratio $`M_{\mathrm{disk}}/M_{}`$. This is unlikely to be true for brown dwarfs around very low mass stars. Therefore, no brown dwarf desert is expected around M stars.
* The same argument applies to M star companions to stars of a few Solar masses – if the latter have disks that live for several viscous times. Migration could clear out a low mass stellar desert at small separations around stars modestly more massive than the Sun.
* Migration takes of the order of a Myr. No brown dwarf desert is expected around the youngest pre-main-sequence stars.

We emphasize that the first two predictions are not unique to this model. Accretion onto forming binary systems also abhors extreme mass ratios (Bate 2000), and provides an alternative theoretical explanation for the brown dwarf desert. The third prediction is more discriminatory. Detection of a higher frequency of brown dwarfs around young pre-main-sequence stars would provide strong support for the model presented here. Competing theories based on differential accretion (Bate 2000) and dynamical interactions (Reipurth & Clarke 2001) predict that the desert is established on a dynamical timescale, which is at least an order of magnitude shorter than the timescale on which migration occurs.

## References

Armitage, P. J., & Bonnell, I. A. 2002, MNRAS, 330, 11

Armitage, P. J., & Hansen, B. M. S. 1999, Nature, 402, 633

Armitage, P. J., Livio, M., Lubow, S. H., & Pringle, J. E. 2002, MNRAS, 334, 248

Bate, M. R. 2000, MNRAS, 314, 33

Boss, A. P. 1998, ApJ, 503, 923

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Goldreich, P., & Sari, R. 2002, ApJ, submitted

Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153

Lin, D. N. C., & Pringle, J. E. 1990, ApJ, 358, 515

Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001

Marcy, G. W., & Butler, R. P. 2000, PASP, 112, 137

Osterloh, M., & Beckwith, S. V. W. 1995, ApJ, 439, 288

Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263

Pringle, J. E. 1981, ARA&A, 19, 137

Pringle, J. E. 1991, MNRAS, 248, 754

Reipurth, B., & Clarke, C. J. 2001, AJ, 122, 432

Strom, K. M., Strom, S. E., Edwards, S., Cabrit, S., & Skrutskie, M. F. 1989, AJ, 97, 1451

Trilling, D. E., Benz, W., Guillot, T., Lunine, J. I., Hubbard, W. B., & Burrows, A. 1998, ApJ, 500, 428

Trilling, D. E., Lunine, J. I., & Benz, W. 2002, A&A, 394, 241


text:
# Temperature measurement at the end of a cantilever using oxygen paramagnetism in solid air

## Abstract

We demonstrate temperature measurement of a sample attached to the end of a cantilever using cantilever magnetometry of solid air “contamination” of the sample surface. In experiments like our Magnetic Resonance Force Microscopy (MRFM), the sample is mounted at the end of a thin cantilever with small thermal conductance. Thus, the sample can be at a significantly different temperature than the bulk of the instrument. Using cantilever magnetometry of the oxygen paramagnetism in solid air provides the temperature of the sample, without any modifications to our MRFM (Magnetic Resonance Force Microscopy) apparatus.

In many Magnetic Resonance Force Microscopy (MRFM)GaAsMRFM ; MRFM1 ; CaF ; MRFM2 ; MRFM3 ; MRFM4  and cantilever magnetometry experiments,magnetometry  the sample being measured is mounted at the end of a thin cantilever. The cantilever often has very poor thermal conductance because of its small cross section in comparison to its length. Even a small heat input into the sample can raise its temperature significantly from the temperature at the opposite end (base) of the cantilever. Obviously, the temperature of the sample can be critical to the results of an experiment because of many temperature dependent factors like spin polarization or time constants.MRFM1 ; CaF  Mounting a temperature sensor at the end of the cantilever would be difficult, requiring a sensor of $``$10 $`\mu `$m dimensions integrated with the cantilever. An indirect way to measure the temperature of a cantilever is to measure the thermal (Brownian) motion of the cantilever.brownian  Using simple Brownian motion requires that the spring constant of the loaded cantilever is known and that the temperature is uniform along the cantilever and equal to the sample temperature. In addition, the Brownian motion method relies on the absence of non-thermal sources of cantilever motion for correct temperature measurement. Sources of sample heating (such as RF magnetic fields and optical power) can also be sources of non-thermal motion.

A simple alternative is to use cantilever magnetometry to measure the paramagnetism of oxygen in frozen air on the sample. This technique does not require any modification of the experimental probe used for MRFM or cantilever magnetometry, and provides a relative measure of the sample temperature up to the sublimation temperature of the solid air. In essence, the temperature is measured by using the temperature dependence of the magnetic susceptibility of the air “contamination” of the sample surface. The paramagnetic electron spins of the oxygen molecules in the solid air have a large, easily measureable temperature dependence of magnetic susceptibility. This technique could also be done just using the susceptibility of the sample itself if there is enough temperature dependence of the susceptibility.

The experimental setup is the same as used for our sample on cantilever MRFM experiments.GaAsMRFM  The sample for MRFM (in this case GaAs) is mounted with thermally conductive silver epoxy on a gold-coated SiNx cantilever with spring constant, k$``$0.05 N/m. The SiNx cantilever has very poor thermal conductance because of the low thermal conductivity of SiNx and thin cross section (the cantilever is 0.6 $`\mu `$m thick, 22 $`\mu `$m wide, and 320 $`\mu `$m long).cantilever  In order to increase the thermal conductivity of the cantilever, it was coated on each side with either 300 Å/700 Å or 300 Å/1700 Å of Ti/Au. Coating on both sides is required to reduce bending of the cantilever from the differential contraction of the gold and SiNx on cooling. The position of the cantilever is measured with an optical interferometer.CaF  Below the sample is the source of the magnetic field gradient, a 250 $`\mu `$m diameter iron cylinder. The magnetic field gradient ($`\delta B_z/\delta z=6000`$ to 9000 T/m) applies a static force,

$$F_z=\mu _z\frac{\delta B_z}{\delta z}$$
(1)
to the net magnetic moment of the oxygen electrons, $`\mu `$, which we measure by the deflection of the cantilever. Raw data from the interferometer is shown in figure 1. As the magnetic field is increased, the oxygen electrons become more polarized and the force between the oxygen and the magnetic field gradient becomes larger. The cantilever deflects towards the iron magnet. Raising the temperature decreases the susceptibility and thus the amount of deflection (figure 1(c)).

This technique of temperature measurement could be done with the sample itself, if the sample has sufficient temperature dependence of its magnetic susceptibility for easy measurement. However, “contamination” of the sample and cantilever surface with solid air can be used with an arbitrary sample. In addition, air is a resource available in almost every lab. Coating the sample with a small amount of frozen air can easily be done by not fully pumping out the probe before cool down. Because of the sensitivity of the cantilever, very little air is required. Just 1011 oxygen molecules ($``$monolayer) will result in several interferometer fringes of deflection (hundreds of nm) during a magnetic field sweep to 9 T. Deliberately allowing air into the vacuum system with a non-regulating valve typically resulted in too much air contamination and the magnetic force would stick the cantilever to the iron magnet.

This provides a cautionary warning for cantilever magnetometry or MRFM experiments. Oxygen contamination must be avoided for sensitive magnetometry or MRFM experiments because of the large magnetic susceptibility of oxygen. For this reason, the temperature measurements were done first and then the air contamination removed before the MRFM experiments. The air contamination could be removed by heating the probe head to about 20 K while pumping. The best technique for avoiding air contamination was to use the heater to keep the probe head warm, while the rest of the probe was cooled down.

We should also note that solid air is a better choice for large and simple temperature dependence of the magnetic susceptibility than pure oxygen. Pure solid oxygen has antiferromagnetic interactions between the molecules and undergoes two phase transitions in the solid state at low temperatures (43.8 and 23.9 K).EPR ; adsorbed  By contrast, the oxygen electrons in solid air have been studied recently with EPR at 5 K and shown to behave as $`S=1`$ paramagnetsEPR  following the Hamiltonian

$$H=g\mu _BBS_z+D\left(S_z^2\frac{S(S+1)}{3}\right)$$
(2)
where $`g=2.0`$ and the zero field splitting $`D=5.1`$ K. The Hamiltonian is written for the case where the axis of the oxygen molecule is along the magnetic field. In the frozen air, the oxygen molecules have a distribution of their molecular axes with respect to the magnetic field.EPR  Apparently, the oxygen molecules are dispersed enough among the nitrogen molecules in the solid air to prevent antiferromagnetic ordering. Similarly, dilute adsorbed oxygen results in roughly Curie-Weiss law behavior of the susceptibility, $`1/(T\theta )`$, with $`\theta `$ small and negative.adsorbed  Figure 2 shows the deflection of the cantilever up to high field. The saturation of the susceptibility at high field is consistent with spin 1 and $`g2`$ electronic moments. To provide a calibration for this temperature measurement method, we heated the probe head, which results in a nearly uniform temperature of the probe head, cantilever, and sample. The probe head itself has its temperature measured with a resistive sensor.tsensor  As shown in figure 3, the temperature dependence of the magnetic susceptibility does indeed roughly follow a Curie-Weiss law with the fitted $`\theta =2`$ K.

Then, we applied this technique to measure the sample temperature under conditions when it is different from the temperature of the rest of the probe head. The magnetic susceptibility of the sample with air “contamination” is measured by the deflection of the cantilever by a change in magnetic field. The temperature is derived by comparison with the calibration done by heating the entire probe head. The calibration step could be avoided by assuming the Curie-Weiss law temperature dependence of the oxygen susceptibility and the dominance of the oxygen susceptibility over any other magnetic susceptibility of the sample or cantilever. In these MRFM experiments,GaAsMRFM  the GaAs sample can be heated from two major sources. First, we are applying RF magnetic field which can heat the sample and gold-coated cantilever directly by eddy current heating and indirectly through gas conduction. Second, and more significantly, we are shining light on the sample for optical pumping.OP  Another smaller heat source is the optical power used in the interferometer measurement. We measured the susceptibility and found, for instance, that roughly 4 Gauss of RF (51.5 MHz) magnetic field raised the sample temperature about 1 K. Even more dramatic is the effect of optical pumping (Fig. 4). The $``$80 $`\mu `$W of optical power on the GaAs sample at the end of the cantilever can heat it by 3 K, while the main bulk of the probe has warmed by less than 0.1 K. This local heating is consistent with the estimated thermal conductance of the cantilever, $``$25 $`\mu `$W/K, which is dominated by the 700 Å gold coating.

We have presented a technique for measuring the temperature of a sample at the end of a cantilever. In general, oxygen contamination must be avoided for cantilever magnetometry and MRFM experiments. But here, by measuring the susceptibility of the oxygen in solid air contamination, we can determine the temperature without any contact to the sample. This technique can be extended by using the susceptibility of some other “contaminant” or just the sample itself. In addition, these measurements are done without any modification of the apparatus which we use primarily for MRFM measurement of the nuclear polarization of the sample.

This work was supported by the DARPA Defense Science Office Spins in Semiconductors program. The cantilevers were coated by Monica Taysing-Lara. The authors would like to thank John A. Marohn, John Sidles, and Dan Rugar for many helpful discussions.


text:
# 1 Introduction

## 1 Introduction

Hadronic reactions at low momentum transfer and high energies (for charge-even exchange) are described in QCD in terms of QCD-Pomeron described by the BFKL equation . The charge-odd exchange is less well understood although the corresponding BKP equations  have attracted much attention recently , thus reviving the relevance of phenomenological studies of the Odderon exchange pointed out years ago in Ref. . Unfortunately, the recent studies of specific channels where the QCD Odderon contribution is expected to be singled out have turned out to be very disappointing. In particular, recent experimental studies at HERA of exclusive $`\pi ^0`$ photoproduction  indicate a very small cross section for this process which stays in contradiction with theoretical predictions based on the stochastic vacuum model .

The general feature of all meson production processes is that scattering amplitude describing Odderon exchange enters quadratically in the cross section. This observation lead to the suggestion in Ref. , that the study of observables where Odderon effects are present at the amplitude level - and not at the squared amplitude level - is mandatory to get a convenient sensitivity to a rather small normalization of this contribution. This may be achieved by means of charge asymmetries, as for instance in open charm production . Since the final state quark-antiquark pair has no definite charge parity both Pomeron and Odderon exchanges contribute to this process. Another example  is the charge asymmetry in soft photoproduction of two pions. Bearing in mind perturbative QCD (pQCD) description, we calculated the ”hard” analogs of these asymmetries , and supplemented them by the single spin asymmetries, which may be studied at HERA with polarized lepton beam .

## 2 Method and Results

We consider the process $`e^{}(p_e)N(p_N)e^{}(p_e^{})\pi ^+(p_+)\pi ^{}(p_{})N^{}(p_N^{})`$. The application of pQCD for the calculation of a part of this process is justified by the presence of a hard scale: the squared mass $`Q^2=(p_ep_e^{})^2`$ of the virtual photon, $`Q^2`$ being of the order of a few GeV2. The amplitude of this process (Fig.1)

includes the convolution of a perturbatively calculable hard subprocess with two non-perturbative inputs, the 2-pion generalized distribution amplitude (GDA) and the Pomeron-Odderon (P/O) proton impact factors. Since the $`\pi ^+\pi ^{}`$ system is not a charge parity eigenstate, the GDA includes two charge parity components and allows for a study of the corresponding interference term. The relevant GDA is just given by the light cone wave function of the two pion system .

### 2.1 Charge asymmetry

We define the forward-backward or charge asymmetry $`A(Q^2,t,m_{2\pi }^2,y,\alpha )`$ by

$`{\displaystyle \frac{\underset{\lambda =+,}{}\mathrm{cos}\theta d\sigma (s,Q^2,t,m_{2\pi }^2,y,\alpha ,\theta ,\lambda )}{\underset{\lambda =+,}{}𝑑\sigma (s,Q^2,t,m_{2\pi }^2,y,\alpha ,\theta ,\lambda )}}`$ (1)
as a weighted integral over polar angle $`\theta `$ of the relative momentum of two pions. Although this asymmetry depends on the full set of the kinematical variables, different dependencies, due to factorization, come from different sources. The most clean one are the dependencies on $`Q^2`$, coming from the hard subprocess, and on dipion mass $`m_{2\pi }=\sqrt{(p_++p_{})^2}`$, coming from GDA. The specific form of $`m_{2\pi }`$ dependence is explained by the fact, that the phase of GDA should add to the phase shift between Pomeron and Odderon. The $`Q^2`$ dependence, which is a subject to corrections from BFKL, BKP and ERBL evolution, is due to the different $`Q^2`$ dependence of Pomeron and Odderon coefficient functions (perturbative impactfactors), the latter having the extra propagators, leading to the decrease of asymmetry with $`Q^2`$. The typical scale of $`\alpha _s`$ is determined by the gluons transverse momenta $`\stackrel{}{k_i}`$ (see Fig.1) and we checked that the result is not changed substantially if ”hard” $`\alpha _S(Q^2)`$ is used for coefficient function and ”soft” $`\alpha _S0.5`$ is used for non-perturbative proton impactfactors.

The latter objects, together with the coefficient functions, define the dependence of asymmetry on $`t=(p_Np_N^{})^2`$. As to the dependence on the standard leptonic variables $`y`$ (the lepton energy fraction carried by the virtual photon) and $`\alpha `$ (the angle between lepton and photon scattering planes), it is determined by the relative size of the virtual photon helicity amplitudes and allows to estimate the promising kinematical region .

### 2.2 Spin Asymmetry

The single spin asymmetry $`A_S(Q^2,t,m_{2\pi }^2,y,\alpha )`$ is defined by

$`{\displaystyle \frac{\underset{\lambda =+,}{}\lambda \mathrm{cos}\theta d\sigma (s,Q^2,t,m_{2\pi }^2,y,\alpha ,\theta ,\lambda )}{\underset{\lambda =+,}{}𝑑\sigma (s,Q^2,t,m_{2\pi }^2,y,\alpha ,\theta ,\lambda )}}`$ (2)
and requires to fix the lepton beam polarization $`\lambda `$. Contrary to charge asymmetry, this effect is proportional to the imaginary, rather than to real part of the interference term. As the Pomeron amplitude is imaginary and the Odderon one is real the relative phase between them is the maximal one for the emergence of single spin asymmetries . The effect should be therefore maximal for zero relative phase between isoscalar and isovector GDA’s, providing a complementary probe. This complementarity can be seen from the dependence of spin asymmetry on $`m_{2\pi }`$ and $`Q^2`$(Fig.2)

The smaller numerical value of spin asymmetry is to a large extent due to the kinematical factor $`\sqrt{t}`$. This asymmetry may be therefore important in the region of large $`t`$, which would become a relevant QCD scale and allow to have smaller $`Q^2`$.

## 3 Conclusions

We found that a sizable charge asymmetry may be a useful tool to look for QCD Odderon contribution at HERA. The spin asymmetry is smaller, but it can be important at larger $`t`$ and smaller $`Q^2`$. Note finally, that the numerical value of our predictions depend on the adopted model for proton impact factors, so the observation of the predicted effects with the different magnitude and/or $`t`$dependence might be considered as an indirect experimental determination of these important non-perturbative objects.


text:
# Photoproduction of the 𝜔 meson on the proton at large momentum transfer

## Abstract

The differential cross section, $`d\sigma /dt`$ for $`\omega `$ meson exclusive photoproduction on the proton above the resonance region ($`2.6<W<2.9`$ GeV) was measured up to a momentum transfer $`t=5`$ GeV2 using the CLAS detector at Jefferson Laboratory. The $`\omega `$ channel was identified by detecting a proton and $`\pi ^+`$ in the final state and using the missing mass technique. While the low momentum transfer region shows the typical diffractive pattern expected from Pomeron and Reggeon exchange, at large $`t`$ the differential cross section has a flat behavior. This feature can be explained by introducing quark interchange processes in addition to the QCD-inspired two-gluon exchange.

\[

\]

In this paper we report results of the first comprehensive measurement of the cross section for $`\omega `$ meson photoproduction on protons for $`E_\gamma `$ between 3.19 and 3.91 GeV over the $`t`$ range $`0.15.0`$ GeV2. Previous studies at DESY, SLAC, and NINA  are sparse and cover a limited kinematic range of $`t<1`$ GeV2  and $`tt_{max}`$ (4-5 GeV2) . The low momentum transfer data ($`t<1`$ GeV2) shows a diffractive behavior that can be interpreted in the framework of the Vector Meson Dominance (VMD) model  as the elastic scattering of vector mesons off the proton target. In a more recent approach, this process is also described by the $`t`$-channel exchange of the Pomeron and the dominating $`\pi `$ Regge trajectory . Other approaches  based on effective Lagrangians and inclusion of nucleon resonances as predicted by quark model calculations, are able to reproduce the data at lower photon energies. At high $`t`$, where the cross section is sensitive to the microscopic details of the interaction, the underlying physics can be described using parton degrees of freedom. The onset of this regime can be tested by a combined analysis of different flavor channels. The recent JLab measurements of $`\varphi `$  and $`\rho `$  photoproduction cross sections at large momentum transfer show a behavior consistent with a QCD-inspired framework . At large $`t`$, the small impact parameter ($`1/\sqrt{t}`$) prevents the constituent gluons (quarks) of the exchange from interacting and forming a Pomeron (Reggeon). Because of the dominant $`s\overline{s}`$ component of the $`\varphi `$, quark exchange is strongly suppressed in this channel by the OZI rule and the two-gluon mechanism dominates (Fig. 1-a-b) . In contrast, the light quark composition of the $`\rho `$ allows valence quarks to be exchanged between the baryon and the meson states (Fig. 1-c) . The same quark exchange mechanism is predicted to dominate the $`\omega `$ sector. Complete and detailed measurements of the $`\omega `$ differential cross section are therefore a stringent test of this conjecture.

The measurement was performed at Jefferson Lab with a bremsstrahlung photon beam produced by a continuous electron beam of $`E_0`$ = 4.1 GeV hitting a gold foil of $`10^4`$ radiation lengths. A bremsstrahlung tagging system , with a photon energy resolution of 0.1$`\%`$ $`E_0`$, was used to tag photons in the energy range from $`34`$ GeV. The target cell, a mylar cylinder of 6 cm in diameter and 18 cm long, was filled with liquid hydrogen at 20.4 K. The high-intensity photon flux ($`410^6\gamma `$/s) was continuously monitored during data taking by an $`e^+e^{}`$ pair spectrometer located downstream of the target. The efficiency of this device was determined during dedicated low intensity ($`10^5\gamma `$/s ) runs by comparison with a 100% efficient lead-glass total absorption counter. The systematic uncertainty of the photon flux has been estimated to be 5$`\%`$.

The hadrons were detected in CLAS  (CEBAF Large Acceptance Spectrometer), a spectrometer with nearly 4$`\pi `$ coverage with a toroidal magnetic field ($`1T`$) generated by six superconducting coils. The field was set to bend the positive particles away from the beam into the acceptance of the detector. Three drift chamber regions allowed tracking of charged particles , and time-of-flight scintillators (TOF) were used for hadron identification . The momentum resolution was of the order of a few percent, while the detector geometric acceptance was about 70% for positive hadrons. Low energy negative particles, however, were mainly lost at forward angles because they were bent out of the acceptance. Coincidences between the photon tagger and the CLAS detector (TOFs) triggered the recording of hadronic interactions. From a total of 70M triggers, 100k events were identified as $`p\omega `$ candidates.

For this analysis we chose the most sizeable $`\omega `$ decay mode ($`\omega \pi ^+\pi ^{}\pi ^0`$, b.r. 88.8$`\%`$), requiring detection of both the proton and the $`\pi ^+`$ in CLAS. The data analysis consisted of two main steps: two-pion background rejection and $`\omega `$ yield extraction from the multi-meson background. Due to the different dynamics governing the low and the high $`t`$ domains, we divided the data set into two samples corresponding to low ($`t<1`$ GeV2) and high ($`t>1`$ GeV2) momentum transfer. The analysis procedure was then performed and optimized independently for the two samples.

The two-pion background is dominated by the $`\gamma pp\rho ^0`$ channel since its cross section is five times larger than that for $`\gamma pp\omega `$ for $`E_\gamma `$3-4 GeV, and the mass of the $`\rho `$ meson (770 MeV) is very close to the $`\omega `$ mass (783 MeV). Even though the $`\rho `$ has a larger width ($``$150 MeV FWHM) compared to the $`\omega `$ ($``$8 MeV FWHM enlarged to $``$55 MeV FWHM by the experimental resolution), the missing mass for the reaction $`\gamma ppX`$ alone does not allow separation of the two channels. The two-pion background was rejected by requiring that the missing mass for the reaction $`\gamma pp\pi ^+X`$ be larger than 0.3 GeV. We estimated that the $`\omega `$’s surviving this cut were around 99$`\%`$. Figure 2 shows the ($`p\pi ^+`$) missing mass squared spectrum: the missing $`\pi ^{}`$ peak was easily removed (the hatched area corresponds to the retained events). The small contamination surviving the cut (estimated to be around 5$`\%`$ by the simulations) is spread over a wide proton missing mass interval, and it was reduced to a negligible level in the second step of the analysis. The $`\omega `$ yield extraction from the multi-meson background was performed on the proton missing mass spectrum by using two different procedures: a Gaussian fit to the $`\omega `$ peak and a side-band subtraction.

Both of them rely on the hypothesis of a smooth and continuous background variation from one side-band region to the other. The two methods were not totally independent, but the comparison of their results allowed estimation of the systematic error related to the $`\omega `$ identification. The proton missing mass in each $`t`$ bin was fitted to a Gaussian curve (the $`\omega `$ peak) plus a fourth order polynomial (the multi-meson background). The $`\omega `$ yields in each $`t`$ bin were the area under the Gaussian. Figure 3 shows the fitted spectra in a low and high $`t`$ bin. The side-band subtraction procedure allowed extraction of a localized signal over an extended background subtracting the regions on either side of the peak (side bands) after a proper normalization. The middle region was fixed at 6$`\sigma `$ of the Gaussian curve describing the $`\omega `$ peak ($`\sigma 24`$ MeV) while the side bands had a width of 3$`\sigma `$ each. The $`\omega `$ yield was obtained as average of the two procedures while the maximum difference, 8$`\%`$, was used as an estimate of the systematic error.

The CLAS acceptance and reconstruction efficiency were evaluated with Monte Carlo simulations using the event generator of Ref. . This code included the main contributions to the $`p\pi ^+\pi ^{}`$ ($`\gamma pp\rho ^0`$, $`\gamma p\mathrm{\Delta }^{++}\pi ^{}`$, and $`\gamma pp\pi ^+\pi ^{}`$ in $`s`$-wave) and $`p\pi ^+\pi ^{}\pi ^0`$ final states ($`\gamma pp\omega `$ and $`\gamma pp\pi ^+\pi ^{}\pi ^0`$ phase space), along with background reactions with four or more pions. The generated events were processed by a GEANT-based code simulating the CLAS detector, and reconstructed using the same analysis procedure that was applied to the raw data. The acceptance was derived as a function of $`E_\gamma `$ and the momentum transfer $`t`$, integrating over the remaining independent variables. To minimize the model dependence in the acceptance calculation, the $`\gamma pp\omega `$ differential cross section was iteratively determined from the data and implemented in the Monte Carlo code. The final state ($`p`$ and $`\pi ^+`$ detected) did not allow us to measure the $`\omega `$ decay, therefore the available experimental data about the decay matrix elements , as well as the general decay property of vector mesons , was implemented in the event generator. The systematic error associated with the efficiency calculation was estimated by comparing the results obtained after generating events with slightly different distributions both in production and decay. The resulting systematic uncertainty was estimated to be $`10\%`$. The average acceptance of CLAS for detected $`p\pi ^+`$ ranged from $`8\%`$ to $`10\%`$. For the very forward angles ($`t<0.1`$ GeV2) and the very backward angles ($`tt_{max}`$) the CLAS detector had no acceptance for this reaction.

The $`\omega `$ photoproduction cross section as a function of $`t`$ was extracted in four energy bins in the range 3.19-3.91 GeV. Data are shown in Fig. 4 and 5: vertical error bars include both the statistical uncertainties (ranging from 2$`\%`$ to 25$`\%`$) and the overall systematic error (14$`\%`$) summed in quadrature, while the horizontal bars reflect the bin sizes. In the low momentum transfer region, $`0.1<t<0.5`$ GeV2, good agreement with the previous measurement of Ref.  in a similar energy range is evident. At higher $`t`$ the CLAS data lie between the two data sets taken respectively at smaller and larger energy. Assuming an exponential $`Ae^{Bt}`$ behavior in the range $`0.1<t<0.5`$ GeV2, the coefficient resulting from this experiment, $`B=5.4\pm 0.6`$ GeV-2, is consistent with the values $`B=5.1\pm 1.4`$ GeV-2 and $`B=7.1\pm 1.7`$ GeV-2 obtained by fitting respectively, the $`E_\gamma =2.8`$ GeV and $`E_\gamma =4.7`$ GeV data sets reported in Ref. . Good agreement is also found with existing data at the largest momentum transfer taken at NINA  with a bremsstrahlung photon beam and a single arm spectrometer.

Predictions of the QCD-inspired model of Refs.  are also shown in Fig. 4. Here the Pomeron exchange has been replaced by the exchange of two non-perturbatively dressed gluons (dotted line). The low momentum transfer region is dominated by the pion exchange that, added to the two-gluon and $`f_2(1270)`$ trajectory exchanges, gives good agreement up to $`t`$ 1 GeV2. The $`\pi `$ exchange gives a strong contribution because of the large coupling constant $`g_{\omega \pi \gamma }`$ (0.334). Close to the upper kinematic limit ($`tt_{max}`$) the cross section is well reproduced by the exchange of the nucleon Regge trajectory in the $`u`$ channel . At intermediate momentum transfer, the two gluon exchange contribution underestimates (by an order of magnitude) the experimental cross section. The calculation uses the same expression as in our phi-photoproduction work , where only the relevant mass and radiative decay width have been changed. In contrast to the $`\varphi `$ meson, quark interchanges (Fig. 1-c) are not forbidden in $`\omega `$ production. As explained in Refs.  these hard-scattering mechanisms can be incorporated in an effective way by using the so called “saturated” trajectory that is independent of $`t`$ at large momentum transfer . Regge trajectories are usually assumed to be linear in $`t`$, but there are both phenomenological and theoretical arguments supporting the idea of non-linear trajectories . Saturated trajectories lead to the asymptotic quark counting rules  that, model independently, determine the energy behavior of the cross section at large $`t`$. This approach was successfully adopted to explain the large momentum transfer $`hadronhadron`$ interactions , as well as several photon-induced reactions . The pion saturating trajectory ($`\alpha _\pi ^{sat}(t)`$ -1 when $`t`$ \- $`\mathrm{}`$) is in a form that reproduces the $`\gamma pn\pi ^+`$ reaction around $`\theta _\omega ^{}`$=90 . The solid line in Fig. 4 shows the full calculation, including such a saturating trajectory, while the dot-dashed line corresponds to the same calculation with linear trajectories. Quark exchange increases the cross section at large $`t`$ by more than one order of magnitude.

The measured $`d\sigma /dt`$ in the other three photon energy bins are shown in Fig. 5. From the four data sets, the cross section at $`\theta _\omega ^{}`$=90 was extracted as a function of energy. A power law fit $`s^C`$ to $`d\sigma /dt`$ at $`\theta _\omega ^{}`$=90 was performed also using the only other datum available in the literature (SLAC datum at $`s`$ = 6.13 GeV2 ). The experimental points include both statistical and systematic errors summed in quadrature. The fit yields $`C=7.2\pm 0.7`$ ($`\chi ^2=0.5`$). It is the first time that such a power law behavior, seen for other exclusive reactions , has been observed in the $`\omega `$ channel. The quark exchange diagrams of Fig. 1-c-left (point-like interaction) and 1-c-right (hadronic component of the photon) have a $`s^7`$ and $`s^8`$ power-law behavior, respectively, both in dimensional counting  and in recent models . Note that the saturated $`\pi `$ Regge trajectory behaves like $`s^8`$, too. Besides the differential cross section at fixed energy, the $`s`$ dependence is a strong hint of the presence of quark interchange hard mechanisms in addition to the two-gluon exchange process.

In conclusion, elastic photoproduction of the $`\omega `$ mesons from the proton was measured for the first time with nearly complete kinematic coverage. The energy power law behavior of the differential cross section at $`\theta _\omega ^{}`$=90 was observed. The comparison with a QCD-inspired model, able to reproduce the $`\varphi `$ and the $`\rho ^0`$ photoproduction data, provides further evidence for the presence of hard processes. Adopting a QCD language in this energy region, the two-gluon exchange mechanism (that fully describe the $`\varphi `$ photoproduction data) badly miss the cross section at large momentum transfer and its energy dependence. Good agreement is achieved when quark interchange processes, suppressed in the $`\varphi `$ channel and weakly contributing in the $`\rho `$ case, are included in an effective way in the calculation of the $`\omega `$ cross section.

We would like to acknowledge the outstanding efforts of the staff of the Accelerator and the Physics Divisions at JLab that made this experiment possible. This work was supported in part by the Italian Istituto Nazionale di Fisica Nucleare, the French Centre National de la Recherche Scientifique and the Commissariat à l’Energie Atomique, the U.S. Department of Energy and the National Science Foundation, and the Korea Science and Engineering Foundation. The Southeastern Universities Research Association (SURA) operates the Thomas Jefferson National Accelerator Facility for the United States Department of Energy under contract DE-AC05-84ER40150.