format: text

text:
# A variable efficiency for thin disk black hole accretion

## 1. Introduction

In a wide class of black hole systems, including QSOs, Seyfert galaxies, and high-state Galactic Black Hole Candidates, the accretion disk is almost certainly geometrically thin and radiatively efficient. The radiation emitted from an annulus of such a disk does not simply equal the change in gravitational binding energy as gas flows across the annulus, but rather includes a contribution from energy transported into the annulus from elsewhere in the disk. As a consequence, the disk structure, along with related quantities such as the radiative efficiency (the fraction of rest mass energy which is radiated during the accretion process), depends upon the location and nature of the boundary conditions at the inner disk edge.

For a black hole accretion disk, a natural location to place the inner boundary condition is at the radius of marginal stability, $`r=r_{\mathrm{ms}}`$. This is the radius inside of which circular orbits are no longer stable, and is at $`6GM/c^2`$ for a non-rotating (Schwarzschild) black hole. In the standard model for black hole accretion (Novikov & Thorne 1973; Paczynski & Bisnovatyi-Kogan 1981; Abramowicz & Kato 1989; Paczynski 2000), one assumes that there is no angular momentum transport across $`r=r_{\mathrm{ms}}`$. The motivation for this zero torque boundary condition (ZTBC) is that the radial flow rapidly becomes supersonic once inside $`r=r_{\mathrm{ms}}`$, and so the material loses causal contact with the rest of the disk. Thereafter, material simply spirals ballistically into the black hole. For this reason, we will refer to the region $`r<r_{\mathrm{ms}}`$ as the plunging region.

The ZTBC limits the accretion efficiency to around 6% (for a non-rotating black hole) – already a substantial figure. However, it is possible that even higher efficiencies may arise. The suggestion is that magnetic fields, which are generated in the disk by the magnetorotational instability (Balbus & Hawley 1991) may become dynamically significant within the plunging region. These fields could then causally connect the plunging region to the rest of the accretion disk, allowing energy and angular momentum to be extracted from gas as it executes its final spiral into the black hole (Krolik 1999; Gammie 1999; Agol & Krolik 2000). Subsequent numerical simulations (none of which are General Relativistic) have confirmed the presence of magnetic torques at the marginally stable orbit, but have reached different conclusions as to their impact on the dynamics of the flow. Global MHD simulations of adiabatic accretion tori (Hawley 2000; Hawley & Krolik 2001) showed outward angular momentum transport continuing within $`r_{\mathrm{ms}}`$, while unstratified simulations suggested more modest effects (Armitage, Reynolds & Chiang 2001; Hawley 2001). The important parameters or numerical differences that cause this different behavior remain to be determined.

In this letter we explore the suggestion (Charles Gammie, private communication) that the dynamics of the flow within the plunging region may depend upon the strength of magnetic fields in the disk at $`r>r_{\mathrm{ms}}`$. We extend our earlier work by presenting ‘boosted’ simulations in which the strength of the magnetic field in the saturated turbulence is enhanced by starting with an initially large seed field. This is a numerical trick, though it may have a physical counterpart if the inner disk is threaded by a large-scale magnetic field. For the purposes of this letter, however, the aim is solely to allow us to study how the dynamics of the plunging region vary with field strength in a controlled fashion. We show that the dynamics of gas in the plunging region, and the validity of the ZTBC, do indeed depend upon the saturation field strength, and discuss the implications for observations of accreting black holes.

## 2. The simulations

We use the ZEUS code to solve the equations of ideal MHD (Stone & Norman 1992a, 1992b; Clarke, Norman & Fielder 1994; Norman 2000) in cylindrical co-ordinates, using a setup very similar to that described previously (Armitage, Reynolds & Chiang 2001). Stated briefly, the computational domain is the wedge given by $`r(r_{\mathrm{in}},r_{\mathrm{out}})`$, $`\varphi (0,\varphi _{\mathrm{max}})`$ and $`z(z_{\mathrm{max}},z_{\mathrm{max}})`$. The effects of General Relativity (in particular the existence of an innermost stable orbit) are modelled in this inherently non-relativistic code by using a pseudo-Newtonian potential (Paczynski & Wiita 1980), $`\mathrm{\Phi }=\frac{GM}{rr_{\mathrm{Sch}}}`$, where $`r_{\mathrm{Sch}}=2GM/c^2`$. The vertical component of gravity is neglected, so the disk has no vertical structure. There is therefore no time-averaged variation of density or magnetic field strength with $`z`$. The equation of state is isothermal, with a sound speed $`c_\mathrm{s}`$ which is assumed uniform and constant. At $`r_{\mathrm{ms}}`$, the ratio of the sound speed to the Keplerian velocity is $`c_s/v_\varphi =0.065`$, slightly cooler than our previous simulations.

The initial conditions for density and velocity were obtained by relaxing a Gaussian radial density profile, using ZEUS in its 1D hydrodynamic mode, to produce a density profile which is in accurate numerical equilibrium. The initial configurations of the 3D MHD simulations were then produced by adding a seed magnetic field. Table 1 shows the computational domain, resolution, initial value of $`\beta `$, and seed field configuration for the simulations. As usual, $`\beta `$ is defined as the ratio of thermal energy density to magnetic energy density. In all cases, we set $`r_{\mathrm{in}}=4GM/c^2`$ and $`r_{\mathrm{out}}=20GM/c^2`$. The radial boundary conditions were such as to permit outflow, and all other boundaries were made periodic. These simulations were then evolved until the magnetorotational instability (MRI) became non-linear and produced fully developed MHD turbulence (in the sense that the MRI has saturated and there is no further long-term growth of the magnetic energy). This turbulence drives accretion of material through the inner radial boundary, and the run is continued until an appreciable amount of material has left the computational domain..

## 3. The dynamics of the flow within the plunging region

### 3.1. Time averaged properties

We consider first the time-averaged properties of our runs. Figure 1 shows time-averaged radial profiles for $`\beta ^1`$ (i.e., the ratio of magnetic field energy density to thermal energy density) and the specific angular momentum $`l`$. As intended, the magnetic field in Run 2, which starts with a relatively strong vertical seed field ($`\beta _z=500`$), saturates with significantly higher relative magnetic field energy than those runs which start with either a weaker initial field or a comparable toroidal field. This is due to the continued boosting of the MRI from the conserved vertical flux that accompanies the initial vertical field (Hawley, Gammie & Balbus 1995). Increasing the size of the computational domain also increases the saturation value of the magnetic field, but by a smaller factor.

The dynamics of the flow – diagnosed using the specific angular momentum – correlates with the saturation field strength. Firstly, consider the unboosted runs (Runs 1,3 and 4; shown by the solid, dashed and dot-dashed lines respectively). As one approaches and enters the plunging region, the radial profile of $`l`$ flattens. If the ZTBC is strictly valid, $`l`$ should be constant within the plunging region. For these runs, we obtain the same result as before (Armitage, Reynolds & Chiang 2001), namely that the ZTBC yields a good but not perfect description of the inner disk dynamics. Indeed, examination of the $`r\varphi `$ component of the magnetic stress tensor shows that the Maxwell stresses at $`r=r_{\mathrm{ms}}`$ do not vanish but are less than the peak stress, which occurs at $`r1012GM/c^2`$, by a factor of 5–8.

On the other hand, angular momentum transport within the plunging region is much more significant in the boosted run (Run 2; dotted line in Fig. 1). The time averaged stress at $`r=r_{\mathrm{ms}}`$ is still less than the peak value (which occurs at $`r8GM/c^2`$), but only by a factor of 2. Significant stresses continue to operate within the region $`r<r_{\mathrm{ms}}`$, and there is evidently a clear violation of the ZTBC.

### 3.2. Variability within individual simulations

Turbulent accretion disks are inherently variable systems and so it is also interesting to examine variability during our simulations. For this purpose, we will compare and contrast Run 1 and Run 2. Both runs display order unity temporal and spatial fluctuations in $`\beta ^1`$. For the unboosted run these fluctuations have only a small effect on $`l`$. For the boosted run shown in Fig. 2, however, there are much larger fluctuations in $`l`$. In particularly, during one timeslice from Run 2, a large part of the plunging region becomes magnetically dominated ($`\beta <1`$). During this period (shown with a thick line in the bottom panels of Fig. 2) $`l`$ declines within the plunging region with the same slope as in the rest of the disk, with the Maxwell stresses staying fairly constant well into the plunging region. Events as dramatic as this are fairly rare (caught in 1 out of the 12 independent timeslices that we analyzed), but it is clear that the stress at $`r=r_{\mathrm{ms}}`$ is a rapidly varying function of time, with much more frequent occurrences of moderate stress.

Extraction of energy and angular momentum from within the plunging region can only occur if the material can remain causally connected to the rest of the disk. Analysis of the simulations shows that large stresses occur within the plunging region when the magnetic fields are strong enough for the peak radial MHD wave speed to exceed the radial inflow velocity throughout a large part of the plunging region. It is worth noting that even during the magnetically dominated event seen in Run 2, the azimuthally averaged MHD wave speeds are always slower that the inflow speeds. The angular momentum transport within this highly inhomogeneous plasma therefore appears to be mediated by low-$`\beta `$ (i.e., higher relative magnetic field) filaments.

## 4. Discussion

### 4.1. Radiative efficiency

Calculating the implied change to the radiative efficiency as a consequence of the magnetic torque at $`r_{\mathrm{ms}}`$ is not straightforward (see the discussion in Hawley & Krolik 2001). A simple approach is to note that the dissipation rate is

$$Q(r)=|\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}\mathrm{ln}r}|T^{r\varphi }𝑑z,$$
(1)
where $`T^{r\varphi }`$ is the local stress. We can then estimate the change in the radiative efficiency by comparing $`Q(r)`$, computed using the stress measured in the simulations, with the form expected in a steady disk if there were no stress at $`r_{\mathrm{ms}}`$.

In practice, the non-steady nature of the flow, and the limited radial extent of the simulations, creates difficulties. An indicative estimate, however, follows from noting that in run 2 the mean magnetic stress is roughly constant between $`r_{\mathrm{ms}}`$ and $`2r_{\mathrm{ms}}`$. Compared to a standard disc model, this implies a 50% increase in the dissipation in the region $`r_{\mathrm{ms}}<r<2r_{\mathrm{ms}}`$, and a 10% increase over $`r_{\mathrm{ms}}<r<10r_{\mathrm{ms}}`$. Significantly larger effects are seen in some individual timeslices.

### 4.2. Comparison with previous work

The simulations presented here suggest that the different dynamics within the plunging region, obtained by Hawley & Krolik (2001) and Armitage, Reynolds & Chiang (2001), are plausibly due to different values in the saturation levels of magnetic fields in the disk. The relatively strong magnetic fields in the former simulations, which obtained $`\beta 510`$, are more comparable to the boosted simulations in the current work, and lead to substantial angular momentum transport within the plunging region. Weaker disk fields, with $`\beta 20`$ and a Shakura-Sunyaev (1973) $`\alpha `$ parameter of $`\alpha 0.020.04`$, produce much less striking effects. Although there are undoubtably other significant differences between the simulations, all work to date is broadly consistent with the hypothesis that significant magnetic coupling to the plunging region requires disk fields with $`\alpha 0.1`$ (and $`\beta 10`$). This applies for disks that are moderately thin, with relative thickness $`h/r0.1`$, which are the only systems that have been simulated so far.

### 4.3. Scaling to thinner disks

Very thin disks are difficult to simulate, and it is not obvious how the threshold $`\beta `$ seen in the current work scales with $`h/r`$. We note, however, that the Alfven speed in the disk scales as, $`v_A\left(\frac{h}{r}\right)\beta ^{1/2}`$. Unless $`\beta `$ in the disk decreases rapidly with decreasing $`h/r`$, the Alfven speed in the gas immediately outside the plunging region will decrease with decreasing disk thickness. Conversely, the inflow speed $`v_r`$ within the plunging region is a fixed function of $`r`$, and is independent of $`h/r`$. Hence, a larger increase in the Alfven speed is needed if thinner disks are to maintain causal contact into the plunging region. It is therefore probable that it becomes increasingly more difficult to extract energy from the plunging region as the disk becomes thinner, although further simulation work is required before drawing robust conclusions.

## 5. Conclusions

In this letter, we have used MHD simulations to investigate the effect of the magnetic field strength in the body of the accretion disk on the dynamics of the material within the plunging region. For a disk with a sound speed corresponding to an $`h/r0.1`$, we find that plausible saturation values of the magnetic field are within a factor of a few of a threshold between two regimes. In the runs which saturated with a fairly low magnetic field ($`\beta 20`$), there was only modest extraction of angular momentum from material within the plunging region, and the ZTBC was a reasonable approximation. In our ‘boosted’ runs, where the field saturated at a higher level ($`\beta 510`$), we found significantly stronger angular momentum transport within the plunging region. Furthermore, our boosted run showed sporadic intervals of very efficient angular momentum transport (with stresses that remained fairly constant well within the plunging region). We note that stratified simulations by other authors (Hawley 2000; Miller & Stone 2000; Hawley & Krolik 2001) suggest that magnetic fields could be stronger (i.e. lower $`\beta `$) above the disk midplane. Such fields, which are not included in our simulations, could contribute to more significant transport within the plunging region.

The extraction of energy and angular momentum from material in the plunging region can significantly increase the efficiency of accretion above the standard value, especially in the innermost regions. This raises the possibility of varying efficiencies, both across systems, and within a given system at different times.

In this paper, we have used a large scale flux threading the accretion disk as a numerical device for altering the saturation level of the MRI. Large scale fields could also, of course, be present in real systems, and are often invoked in models for the formation of jets and outflows (e.g. Ouyed, Pudritz & Stone 1997, and references therein). This raises the interesting possibility that periods of increased luminosity (due to a high radiative efficiency) may be related to periods during which strong outflows occur. Such associations have been reported for Galactic Black Hole Candidates such as GRS 1915+105 (Mirabel et al. 1998).

## Acknowledgments

We extend a special thanks to Charles Gammie for suggesting this problem. We are grateful to the developers of ZEUS and ZEUS-MP for making these codes available as community resources. CSR acknowledges support from Hubble Fellowship grant HF-01113.01-98A. PJA thanks JILA for hospitality.


text:
# Massive Star Evolution: Nucleosynthesis and Nuclear Reaction Rate Uncertainties

## 1 Introduction

Massive stars of more than 8 $`\mathrm{M}_{}`$ are the main source of oxygen and heavier elements in the universe. Availability of (theoretical and experimental) reaction rate data and computational resources make it now possible to follow the complete nucleosynthesis in massive stars from hydrogen burning till core collapse and through the supernova explosion (§2). However, significant uncertainties in several key nuclear reaction rates still exist. In §3 we discuss the influence of the $`{}_{}{}^{22}\mathrm{Ne}`$$`(\alpha ,\mathrm{n})`$ rate on the s-process in massive stars and in §4 we demonstrate the influence of the $`{}_{}{}^{12}\mathrm{C}`$$`(\alpha ,\gamma )`$ rate on the presupernova structure.

## 2 Complete nucleosynthesis study

We present the first calculations to follow the complete nucleosynthesis in massive stars from hydrogen ignition till onset of iron core collapse and through the supernova explosion. Figure 1 shows the average abundance of all ejecta, including stellar wind mass loss, relative to their solar values (production factor). The dashed line indicates the production factor for $`{}_{}{}^{16}\mathrm{O}`$, the dominant “metal” produced by massive stars, and the dotted lines indicate twice and half the production factor. The supernova explosion is simulated by a piston that resulted in $`1.7\times 10^{52}\mathrm{erg}`$ kinetic energy of the ejecta. The mass cut is then determined self-consistently from the hydrodynamical simulation. Note that Fig. 1 does not include a possible r-process contribution due to the neutrino wind from the nascent neutron star.

In the 25 $`\mathrm{M}_{}`$ star of Fig. 1, most isotopes from oxygen to the iron group are produced in about solar ratios relative to oxygen. The iron group itself is somewhat underproduced, due to fall-back during the explosion. Stars of $`15\mathrm{M}_{}`$ typically contribute more here and Type Ia supernovae have added to the solar abundance in the region. The s-process isotopes above the iron group till $`A=90`$ are slightly overproduced. Stars of lower mass and/or lower metallicity produce less here. Therefore these high yields are required to produce the solar abundances over the lifetime of the Galaxy. Above $`A100`$ many p-isotopes are produced by the $`\gamma `$-process in about solar abundances relative to $`{}_{}{}^{16}\mathrm{O}`$. For more details please refer to Rauscher et al. (2001).

## 3 The $`{}_{}{}^{22}\mathrm{Ne}`$$`(\alpha ,\mathrm{n})`$ rate

The $`{}_{}{}^{22}\mathrm{Ne}`$$`(\alpha ,\mathrm{n})`$ rate is the most important source of neutrons for the s-process in massive stars in the region $`60A90`$. Figure 1 shows the result for the lower limit given by Käppeler et al. (1994) (as used by Hoffman et al., 2000). In Fig. 2 we show the results for the same star, except that we use the low, high, and recommended values by Jaeger et al. (2001). The abundance ratio of the s-process does not change much by this; the abundances only shifting to a slightly higher absolute value. The lower panel shows the result when using the high limit given by Angulo et al. (1999) (and the rest of their rate set as well, which makes only a minor difference as compared to $`{}_{}{}^{22}\mathrm{Ne}`$$`(\alpha ,\mathrm{n})`$). The large overproduction of the s-process here cannot be balanced by galactochemical evolution. Even so, we do not find significant production of p-isotopes of $`\mathrm{Mo}`$ and $`\mathrm{Ru}`$ (cf. Costa, 2000). For more details please refer to Heger et al. (2001).

## 4 The $`{}_{}{}^{12}\mathrm{C}`$$`(\alpha ,\gamma )`$ rate

The uncertainty of the $`{}_{}{}^{12}\mathrm{C}`$$`(\alpha ,\gamma )`$ rate has significant influence on the late time evolution of massive stars (cf. Imbriani et al., 2001). It determines how much $`{}_{}{}^{12}\mathrm{C}`$ is left after core helium exhaustion (Fig. 3). Only for a sufficiently high value central carbon burning proceeds convectively while otherwise it burns radiatively. Similarly, the extent and duration of the carbon shell burning phases are affected. They set the stage for the later burning phases in that they produce carbon-free cores of different sizes, as a non-monotonous function of the carbon abundance. In Fig. 4 the iron core size varies by up to 30 %! Therefore the $`{}_{}{}^{12}\mathrm{C}`$$`(\alpha ,\gamma )`$ rate determines whether a neutron star or a black hole is formed. For more details please refer to Boyes et al. (2001).

## 5 Conclusions & outlook

The current extent of uncertainties in key nuclear reaction rates still has significant influence on stellar evolution and nucleosynthesis. Other major uncertainties comprise our understanding of mixing processes, rotation, and magnetic fields in stars. Some of these uncertainties are currently under re-investigation. The combined effort of refined rate determinations and stellar modeling allows both fields to profit from the synergy effects.

Acknowledgments: This research was supported, in part, by the DOE (W-7405-ENG-48), the National Science Foundation (AST 97-31569, INT-9726315), the Alexander von Humboldt Foundation (FLF-1065004), and the Swiss National Science Foundation (2000-061822.00). T.R. acknowledges support by a PROFIL professorship from the Swiss National Science Foundation (2124-055832.98).


text:
# Theoretical Uncertainties in Sparticle Mass Predictions

CERN-TH/2001-255

hep-ph/0110227

Snowmass-2001/P3-19

\]http://allanach.home.cern.ch/allanach/

## Abstract

We contrast the sparticle spectra obtained from three modern publicly available codes along model lines in minimal supersymmetric standard model (MSSM) parameter space. From this we gain an idea of the uncertainties involved with sparticle spectra calculations. The differences in predicted sparticle masses are typically at the several percent-level. In the focus-point scenario, there are differences of $`30\%`$ in the weak gaugino masses. These uncertainties need to be reduced in order to obtain accurate information about fundamental models of supersymmetry breaking.

preprint: P3-19preprint: CERN-TH/2001-255

MSSM phenomenology is so complicated that studies of the ability of future colliders to search for and measure supersymmetric parameters has focused on isolated ‘bench-mark’ model points Allanach et al. (2000a); Battaglia et al. (2001); Abramowicz et al. (1999). This approach, while being a start, is not ideal because one is not sure how many of the features used in the analyses will apply to other points of parameter space. In an attempt to cover more of the available parameter space, the Direct Investigations of SUSY Subgroup of SNOWMASS 2001 has proposed eight bench-mark model lines for study.

The lines have been defined to have the spectrum output from the ISASUGRA program (part of the ISAJET7.51 package Baer et al. (1999)) for $`m_t=175`$ GeV. ISASUGRA solves the MSSM normalization group equations subject to boundary conditions at low energy (measured Standard Model couplings and constraints from consistent electroweak symmetry breaking) and a higher energy scale (the theoretical boundary condition on the soft mass parameters). Knowledge of the uncertainties will be important when data is confronted with theory, i.e. when information upon a high-energy SUSY breaking sector is sought from low-energy data. Here, we intend to investigate the theoretical uncertainties in sparticle mass determination. To this end, we contrast the sparticle masses predicted by three modern publicly available supported codes: ISASUGRA, SOFTSUSY Allanach (2001) and SUSPECT Djouadi et al. (1998).

Each of the three packages calculates sparticle masses in a similar way, but with different approximations Moretti (2001). In each of the model line scenarios, we calculate the fractional difference for some sparticle $`s`$

$$f_s^{\text{CODE}}=\frac{m_s^{\text{ISASUGRA}}m_s^{\text{CODE}}}{m_s^{\text{ISASUGRA}}},$$
(1)
where CODE refers to SOFTSUSY1.2, or SUSPECT2.0. $`f_s^{\text{CODE}}`$ then gives the fractional difference of the mass of sparticle $`s`$ between the predictions of CODE and ISASUGRA. A positive value of $`f_s^{\text{CODE}}`$ then implies that $`s`$ is heavier in ISASUGRA than in CODE.

We focus upon model lines in scenarios which are currently supported by all three packages, i.e. supergravity mediated supersymmetry breaking (mSUGRA). At a high unification scale $`M_{GUT}1.9\times 10^{16}`$, the soft-breaking scalar masses are set to be all equal to $`m_0`$, the universal scalar trilinear coupling to $`A_0`$ and each gaugino mass $`M_{1,2,3}`$ is set. $`\mathrm{tan}\beta `$ is set at $`M_Z`$. The three choices of model lines are displayed in Table 1.

Model line A displays gaugino mass dominance, ameliorating the SUSY flavor problem. Model line B has non-universal gaugino masses and model line F corresponds to focus-point supersymmetry Feng et al. (2000), close to the electroweak symmetry breaking boundary.

We pick various sparticle masses that show a large difference in their prediction between the three calculations. For model line A, Fig. 1a shows $`f_{\tau _1,b_1,q_1,h^0,\chi _1^0,\chi _1^+,g}^{\text{SOFTSUSY}}`$ (the lightest stau, sbottom, squark, neutral Higgs, neutralino, chargino and gluino mass difference fractions respectively). Fig. 1b shows the equivalent results for the output of SUSPECT. Model line B differences are shown in Fig. 2.

Jagged curves in the figures are a result of numerical error in the ISASUGRA calulcation, and are at the 1$`\%`$ level for squarks, gluinos and the lightest neutralino. The stau, lightest Higgs and lightest chargino do not display any appreciable numerical error. We have checked that SOFTSUSY and SUSPECT do not have any numerical error that is detectable by eye for any of the sparticles.

Figs. 1,2 share some common features. In general, the largest discrepancies occur for low $`M_{1/2}`$, where the super-particle spectrum is lightest. The gluino and squark masses are consistently around 5$`\%`$ lower in ISASUGRA than the other two codes, which agree with each other to better than $`1\%`$. We note here that this uncertainty is not small, a $`5\%`$ error on the gluino mass at $`M_{1/2}=700`$ GeV in model A corresponds to an error on the predicted gluino mass of 75 GeV, for example. The lightest neutralino mass predicted by SUSPECT is in agreement with ISASUGRA to $`2\%`$, whereas SOFTSUSY predicts this neutralino to be lighter by another $`2\%`$ or so. The Higgs masses also agree to above 2$`\%`$ between each code and ISASUGRA, but SOFTSUSY and SUSPECT have opposite sign differences. This is probably due to the fact that SOFTSUSY uses a FeynhiggsFast calculation of the neutral Higgs masses with important two-loop effects added Heinemeyer et al. (1999), which predicts masses that tend to be higher than the one-loop calculation (as used in ISASUGRA or SUSPECT). The stau mass also looks fairly consistent between SUSPECT and SOFTSUSY (2% lighter than ISASUGRA), but this could conceivably be due to their one-loop running of scalar masses. The lightest neutralino mass shows large deviations between all codes up to 5$`\%`$, and this is not currently understood.

The focus-point scenario (model line F) is displayed in Fig. 3.

The overall view of spectral differences is similar to that in model lines A and B except for the masses of the lightest chargino and neutralino. They display large 10-30$`\%`$ differences. The large spikes come from numerical errors in the ISASUGRA calculation. In focus point supersymmetry, the bilinear Higgs mass parameter $`\mu `$ is close to zero and is very sensitive to threshold corrections to $`m_t`$ Allanach et al. (2000b). For small $`\mu <M_Z`$, the lightest chargino and neutralino masses become sensitive to its value. The predicted value of $`\mu (M_Z)`$ differs by 50$`\%`$-100$`\%`$ between all of the three code’s output. Only a few of the threshold corrections to $`m_t`$ are included in the ISASUGRA calculation, whereas SOFTSUSY, for example, includes all one-loop corrections with sparticles in the loop. SUSPECT also adds many of the sparticle loop corrections to $`m_t`$.

In order to discriminate high energy models of supersymmetry breaking, it will be necessary to have better than 1$`\%`$ accuracy in both the experimental and theoretical determination of superparticle masses Allanach et al. (work in progress). An alternative bottom-up approach Blair et al. (2001) is to evolve soft supersymmetry breaking parameters from the weak scale to a high scale once they are ‘measured’. The parameters of the high-scale theory are then inferred, and theoretical errors involved in the calculation will need to be minimised. We note that possible future linear colliders could determine some sparticle masses at the per-mille level Aguilar-Saavedra et al. (2001). An increase in accuracy of the theoretical predictions of sparticle masses by about a factor 10 will be necessary.

###### Acknowledgements.

I would like to warmly thank Werner Porod for his comments on the work.


text:
# The Consequences of Non-Normality

## 1 INTRODUCTION

Wilson-type lattice Dirac operators are non-normal, i.e.

$$[D,D^{}]0,$$
(1)
and this fact has several important consequences. The essential difference to normal operators is that they are not diagonalizable by a unitary transformation. Instead they are diagonalizable (if at all) by a similarity transformation of the form

$$\text{diag}(\lambda _1,\mathrm{},\lambda _N)=X^1DX,$$
(2)
with $`N=\text{dim}(D)`$. It follows that for a given eigenvalue $`\lambda _i`$ there are distinct left $`L_i|`$ and right $`|R_i`$ eigenvectors

$$D|R_i=\lambda _i|R_i,$$
(3)
$$L_i|D=\lambda _iL_i|,$$
(4)
in the sense that (in general) $`L_i|(|R_i)^{}`$.

This fact must be taken into account when performing the spectral decomposition

$$D=\underset{i=1}{\overset{N}{}}\lambda _i|R_iL_i|\underset{i=1}{\overset{N}{}}\lambda _i|R_iR_i|,$$
(5)
but also raises the question of an appropriate generalization of the other concepts familiar from the textbook (hermitian) quantum mechanics. For example, if an observable in the continuum is defined as $`\psi _i|A|\psi _i`$, where $`|\psi _i`$ is an eigenstate of the continuum operator $`D`$ — what is the proper generalization for the case where the left and the right eigenvectors are distinct?

In  we advocate the use of $`L_i|A|R_i`$ instead of $`R_i|A|R_i`$ as a kind of improved observable. Note that in lattice field theory this question is not of importance in the continuum limit, because the difference between the two definitions disappears. However, this question was also raised in the domain of solid state physics — see  for an interesting discussion.

## 2 TRUNCATED EIGENMODE EXPANSION

In this contribution we want to examine a different problem, which is also connected to the non-normality of the Wilson-type lattice Dirac operators. The idea to approximate the inverse operator $`D^1`$ by just a few eigenmodes with eigenvalues close to zero is physically appealing. One hopes that the lowest modes dominate physics and that the truncated eigenmode expansion

$$D_k^1=\underset{i=1}{\overset{kN}{}}\frac{1}{\lambda _i}|R_iL_i|$$
(6)
with the ordering

$$\frac{1}{|\lambda _1|}\frac{1}{|\lambda _2|}\mathrm{}\frac{1}{|\lambda _k|}\mathrm{}\frac{1}{|\lambda _N|}$$
(7)
will be a fair approximation of the proper inverse $`D^1=D_N^1`$.

However, in  it was found that such an expansion is unstable and that it does not uniformly saturate by increasing $`k`$. This is illustrated in Fig. 1 where the values of the correlation function for different temporal distances $`\mathrm{\Delta }t`$ are plotted versus $`k`$. (The illustration is just for two dimensional QED on a $`16^2`$ lattice for $`\beta =2`$, but it should be stressed that for the discussion which follows this is irrelevant — the arguments are pure linear algebra and they will be valid for any physical system and parameter set.)

By a closer examination of pseudospectra (for more details about pseudospectra see ) of $`D`$, we identified that the spikes in Fig. 1 correspond essentially to the instabilities due to the non-normality of $`D`$. For non-normal operators the truncated eigenmode expansion in the original biorthogonal basis seems to be problematic. However, one can do better than that.

## 3 SINGULAR VALUE DECOMPOSITION

For any matrix A there exists a so-called singular value decomposition (SVD)

$$A=U\mathrm{\Sigma }V^{}=\underset{i}{}\sigma _i|u_iv_i|,$$
(8)
where $`U`$ and $`V`$ are unitary and $`\mathrm{\Sigma }`$ is the diagonal matrix with positive or zero entries on the diagonal — the “singular values” $`\sigma _i`$. We assume the following ordering:

$$\sigma _1\sigma _2\mathrm{}\sigma _N0.$$
(9)
Among other interesting properties, one can prove the following theorem (Th. 2.5.3. in ):

If $`k<r=\text{rank}(A)`$ and

$$A_k=\underset{i=1}{\overset{k}{}}\sigma _i|u_iv_i|$$
(10)
then

$$\underset{\text{rank}(B)=k}{\mathrm{min}}AB_2=AA_k_2=\sigma _{k+1},$$
(11)
i.e., $`A_k`$ is the closest matrix to $`A`$ that has rank $`k`$. In other words, the truncation of SVD is the optimal approximation to the original matrix for a given $`k`$.

Actually, the drawback is that for a general matrix $`A`$ one needs two sets of vectors: $`|u`$-s and $`|v`$-s. However, the situation is better for the Wilson-type lattice Dirac operators (and their inverses). Using the $`\gamma _5`$-hermiticity property

$$D^{}=\gamma _5D\gamma _5$$
(12)
we are going to prove that the truncated SVD of $`D^1`$ is identical to

$$D_k^1=\gamma _5\underset{i=1}{\overset{k}{}}\frac{1}{\mu _i}|w_iw_i|=\gamma _5H_k^1,$$
(13)
with the ordering according to the absolute values of inverse $`\mu `$-s:

$$\frac{1}{|\mu _1|}\frac{1}{|\mu _2|}\mathrm{}\frac{1}{|\mu _k|}\mathrm{}\frac{1}{|\mu _N|}.$$
(14)
The $`\mu `$-s and $`|w`$-s are the eigenvalues and eigenvectors of the hermitian Wilson-Dirac operator

$$H=\gamma _5D=WMW^{}=\underset{i=1}{\overset{N}{}}\mu _i|w_iw_i|.$$
(15)
The hermiticity of $`H`$ follows from the $`\gamma _5`$-hermiticity property (12). Thus $`H`$ is diagonalizable by a unitary transformation $`W`$ and has real eigenvalues $`\mu _i`$, ordered on the diagonal of the diagonal matrix $`M`$.

Per construction, $`D^{}D`$ is a hermitian matrix and it is easy to see that

$$D^{}D=WMW^{}WMW^{}=WM^2W^{}.$$
(16)
However, by using the singular value decomposition (8) of $`D`$ it follows that

$$D^{}D=V\mathrm{\Sigma }U^{}U\mathrm{\Sigma }V^{}=V\mathrm{\Sigma }^2V^{}.$$
(17)
As $`V`$ is unitary and $`\mathrm{\Sigma }^2`$ diagonal, one must have (with appropriate ordering)

$$\sigma _i^2=\mu _i^2\sigma _i=|\mu _i|,$$
(18)
i.e., the singular values of $`D`$ are just the absolute values of the eigenvalues of the hermitian Wilson-Dirac operator $`H`$. Also, the $`|v`$-s should correspond to the $`|w`$-s (there is an arbitrary phase factor, but it cancels in (13)). Finally, we can rewrite:

$$D=U\mathrm{\Sigma }V^{}=U\text{sign}(M)MW^{}=\gamma _5WMW^{}$$
(19)
which makes it obvious that

$$U=\gamma _5W\text{sign}(M)|u_i=\text{sign}(\mu _i)\gamma _5|w_i.$$
(20)
Hence, the truncated singular value expansion of $`D`$ is

$`D_k`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{k}{}}}\sigma _i|u_iv_i|={\displaystyle \underset{i=1}{\overset{k}{}}}\sigma _i\text{sign}(\mu _i)\gamma _5|w_iw_i|`$ (21)
$`=`$ $`\gamma _5{\displaystyle \underset{i=1}{\overset{k}{}}}\mu _i|w_iw_i|=\gamma _5H_k.`$
As $`D^1`$ also satisfies the $`\gamma _5`$-hermiticity property (12), the same reasoning applies and (13) follows.

## 4 CONCLUSION

Fig. 2 illustrates that, as expected, by using the truncated singular value expansion, non-normal artefacts disappear and the saturation as a function of $`k`$ is superior to the expansion in the eigenmodes of $`D`$. This approach has already been successfully used in a realistic QCD setting .

Acknowledgments: I. H. thanks the DOE’s Institute of Nuclear Theory at the University of Washington for its hospitality and the DOE for partial support during the completion of this work. W. S. is supported by the DFG Graduiertenkolleg “Feldtheoretische Methoden in der Elementarteilchentheorie und Statistischen Physik”.


text:
# SUBSTRUCTURE IN CDM HALOS AND THE HEATING OF STELLAR DISKS

## 1 Introduction

One recent highlight of cosmological N-body simulations has been the discovery that during the merger and accretion events that characterize the assembly of dark matter halos the central regions of accreted halos may survive for several orbital times as dynamically distinct, self-bound entities in the parent halo (Klypin et al 1999, Moore et al 1999, hereafter K99 and M99, respectively). The population of surviving halo cores, or “subhalos”, typically contributes less than $`10\%`$ of the total mass of the system, with the bulk of the mass in a smooth monolithic structure, as envisioned in the pioneering analytical work of White & Rees (1978). Despite the small fraction of the total mass they make up, at any given time a large number of subhalos are expected within the virialized region of a cold dark matter halo. For example, M99 find that up to $`500`$ satellites with circular velocities exceeding $`10`$ km/s may have survived within $`300`$ kpc from the center of the Galaxy. Comparing this with the dozen or so known Milky Way satellites of comparable velocity dispersion implies that most subhalos must have failed to form a significant number of stars.

Even if luminous galaxies fail to “turn on” in most subhalos a potential difficulty has been cited by M99: the fluctuating gravitational potential induced by the clumpy structure of the halo may act to heat and thicken fragile stellar disks beyond observational constraints. I present here the results of an attempt to quantify the effects of substructure on the dynamical evolution of thin stellar disks embedded in dark matter halos.

## 2 Substructure in CDM halos

The mass function of substructure halos, their internal structure, and the parameters of their orbits are the main properties of the subhalo population that determine the tidal effects of substructure on stellar disks. The short-dashed curves in Figure 1 (left panel) show the subhalo velocity function corresponding to galaxy-sized dark matter halos formed in the $`\mathrm{\Omega }=1`$ CDM and in the $`\mathrm{\Lambda }`$CDM cosmogonies. Circular velocity (instead of mass) is used to characterize subhalos because of its weaker dependence on the exact way in which substructure is identified. Still, circular velocities do change as a function of radius from the center of a subhalo, and there is no unique way of defining subhalo circular velocities. Figure 1 reports results for $`V_{\mathrm{peak}}`$ and $`V_{\mathrm{outer}}`$, which correspond to the maximum circular speed within the subhalo, and to its value at the outermost bound radius, respectively. Scaled to the virial velocity of the halo, the substructure velocity function is roughly independent of the mass of the parent halo (M99) and of the cosmological parameters (K99). The results of K99 are shown by the dot-dashed line in Figure 1 (left) and are in reasonable agreement with our determination.

How does the subhalo velocity function compare with the observed number of satellites in the vicinity of the Milky Way? The solid squares in Figure 1 (joined by a solid line) illustrate the cumulative number of known Milky Way satellites as a function of the circular velocity of their halos, as plotted by M99. Here, circular velocities are derived for the halos of dwarf spheroidals assuming that stars in these systems are on isotropic orbits in isothermal potentials. This is a plausible, but nevertheless questionable, assumption. Dark halos differ significantly from simple isothermal potentials, and numerical simulations indicate that circular velocities decrease substantially near the center. If stars populate the innermost regions of subhalos their velocity dispersions may substantially underestimate the subhalo peak circular velocities. This has been discussed by White (2000), who finds, using the mass model proposed by Navarro, Frenk & White (1997, NFW), that dwarf spheroidals may plausibly inhabit potential wells with circular velocities up to a factor of $`3`$ times larger than inferred under the isothermal assumption. Such correction may reconcile, at the high mass end, the Milky Way satellite velocity function with the subhalo function, as shown by the dotted line in Figure 1. Thus, the possibility remains that the number of massive satellites expected in the CDM scenario might not be in gross conflict with observation.

On the other hand, this does not erase the large number of low mass “dark” satellites that should inhabit the main halo. Can thin stellar disks survive unscathed in this clumpy environment? The tidal heating rate by substructure can be shown to scale as $`dE/dtn(m_s)m_s^2𝑑m_s`$, where $`m_s`$ is the subhalo mass (White 2000). Since, according to numerical simulations, $`n(m_s)m_s^{1.8}`$ (Ghigna et al 2000, Springel et al 2000), tidal heating by substructure is dominated by the few most massive subhalos. One consequence of this result is that stochastic effects are expected to dominate the heating rate and a statistically significant sample of systems should be simulated in order to obtain conclusive results. The results reported here (which are based on a couple of realizations) should then be viewed as preliminary until confirmed by further studies.

## 3 The Evolution of Stellar Disks in a $`\mathrm{\Lambda }`$CDM halo

In order to assess the dynamical influence of subhalos on stellar disks we have carried out a series of numerical simulations that follow the evolution of a disk of particles within a dark matter halo with substructure similar to that of the halos shown in Figure 1 (left). We choose parameters so that the model reproduces many of the observational characteristics of the Milky Way galaxy, as in Velázquez & White (1999). The galaxy model is first evolved without substructure halos in an attempt to assess the deviations from equilibrium induced by noise in the particle distribution. N-body disks are notoriously unstable, and substantial numbers of particles are needed to ensure stability over several rotation periods. The simulations we report here use $`40,000`$ particles in the disk, $`13,333`$ in the bulge, and $`2.2\times 10^6`$ in the halo, all of equal mass, $`m_p=1.4\times 10^6M_{}`$. Within this system we insert substructure halos, using a procedure which ensures that the subhalos have approximately the same masses, densities and orbits as those in the CDM halos.

We evolve the disk model first for about $`3.5`$ Gyrs without including substructure in order to quantify the heating rate due to noise in the particle distribution. The ($`R`$, $`z`$, $`\varphi `$) velocity dispersions of disk particles at the solar circle (i.e., $`8.5`$ kpc from the center) grow from ($`31`$, $`27`$, $`26`$) km s-1 to ($`43`$, $`30`$, $`33`$) km s-1 over the same period (see solid lines in Figure 4). Quantifying this heating rate by the usual expression $`\sigma _{tot}^2=\sigma _0^2+Dt`$, we find $`D200`$ km2 s-2 Gyr-1, about a factor of two less than inferred for stars in the solar neighborhood from the age-velocity dispersion (Fuchs et al 2000).

The two-body heating rate in our equilibrium stellar disk model is thus low compared with the actual heating experienced by stars in the disk of the Galaxy. This implies that our disk model is stable enough to verify numerically whether substructure in the halo leads to heating rates inconsistent with observational constraints. This is shown by the dashed lines in Figure 4, which show the evolution of the disk velocity dispersion when the substructure halos are added to the system. Clearly, the heating rate is approximately the same with and without substructure, a result that may be traced to the fact that there are no satellites more massive than $`10^{11}M_{}`$ and that their orbits seldom take them near the disk. If substructure in the halos we consider is representative of galaxy-sized CDM halos (and we have no reason to suspect it is not), this would imply that tidal heating rates of thin stellar disks by substructure halos may be consistent with the observational evidence.

## 4 Discussion

These results suggest that concerns regarding excessive tidal heating of thin stellar disks by substructure in CDM halos may be less serious than previously thought. We conclude that the substructure observed in virialized CDM halos is not clearly inconsistent with the existence of thin stellar disks such as that of the Milky Way. These conclusions are subject to a number of caveats. The most obvious one is that our study explores only two numerical realizations of a disk galaxy within clumpy dark matter halos, and it is always hazardous to extrapolate from such a small number of cases. Our study does show, however, that it is possible at least in some cases to maintain a stellar disk in spite of substructure. A second caveat is that we have explored a model motivated by the present-day structure of the Milky Way and by the $`z=0`$ substructure of a CDM halo. Models that take into account the ongoing formation of the disk and a more realistic treatment of the evolution of substructure are clearly desirable in order to refine the conclusions presented here.

Finally, as discussed by Navarro & Steinmetz (2000), it is quite difficult to account simultaneously for the masses, luminosities, rotation speeds, and angular momenta of galaxy disks in cosmogonies such as CDM, where much of the mass of a virialized halo tends to be assembled through a sequence of mergers. Until these issues are fully resolved it would be premature to extend a clean bill of health to the CDM paradigm regarding the formation and evolution of spiral galaxies like our own Milky Way.

## Acknowledgments

I thank my collaborators, Andreea Font, Joachim Stadel and Tom Quinn for allowing me to present the results of our collaboration. This research has been supported by CIAR, CFI, NSERC and the Alfred P. Sloan Foundation.

## References


text:
# Growing Scale-Free Networks with Tunable Clustering

## Abstract

We extend the standard scale-free network model to include a “triad formation step”. We analyze the geometric properties of networks generated by this algorithm both analytically and by numerical calculations, and find that our model possesses the same characteristics as the standard scale-free networks like the power-law degree distribution and the small average geodesic length, but with the high-clustering at the same time. In our model, the clustering coefficient is also shown to be tunable simply by changing a control parameter—the average number of triad formation trials per time step.

A great number of systems in many branches of science can be modeled as large sparse graphs, sharing many geometrical properties reviews . For example: social networks, computer networks, and metabolic networks of certain organisms all have a logarithmically growing average geodesic (shortest path) length $`\mathrm{}`$ and an approximately algebraically decaying distribution of vertex degree. In addition to this, social networks typically show a high clustering, or local transitivity: If person $`A`$ knows $`B`$ and $`C`$, then $`B`$ and $`C`$ are likely to know each other.

Works on the geometry of social networks, which is the main focus of the present paper, have originated from Rapoport’s studies of disease spreading rapoport , and have been further developed in Refs. post\_rapoport ; FRO . General mathematical models for random graphs with a structural bias are called the Markov graphs and were studied in Ref. markov . In the physics literature, networks with high clustering are commonly modeled by the small-world network model of Watts and Strogatz (WS) WS , while networks with the power-law degree distribution by the scale-free network model of Barabási and Albert (BA) SF . Although both models have a logarithmically increasing $`\mathrm{}`$ with the network size, each model lacks the property of the other model: the WS model shows a high clustering but without the power-law degree distribution, while the BA model with the scale-free nature does not possess the high clustering. In this work, we propose a network model which has both the perfect power-law degree distribution and the high clustering. Furthermore, in our model, the degree of the clustering, measured by the clustering coefficient (see below), is shown to be tunable and thus controllable by adjusting a parameter of the model.

We start from the definition of a network as a graph $`𝒢=(𝒱,)`$, where $`𝒱`$ is the set of vertices and $``$ is the set of edges somegraphtheorybook . An edge connects pairs of vertices in $`𝒱`$ and not more than one edge may connect a specific pair of vertices. To quantify the clustering, Watts and Strogatz introduced the clustering coefficient $`\gamma \gamma _v`$ with the average $`\mathrm{}`$ for all vertices in $`𝒱`$. The local clustering coefficient $`\gamma _v`$ for the vertex $`v`$ is defined as follows: Suppose that the vertex $`v`$ has $`k_v`$ neighbors ($`k_v`$ is called the degree of the vertex $`v`$, a neighbor is a vertex separated by exactly one edge). For those $`k_v`$ neighbors, there can exist at most $`\left({\displaystyle \genfrac{}{}{0pt}{}{k_v}{2}}\right)=k_v(k_v1)/2`$ edges connecting two of $`k_v`$ vertices. If one defines $`|(\mathrm{\Gamma }_v)|`$ as the number of actual edges existing in the network connecting those neighbors, the local clustering coefficient is written as WS

$$\gamma _v\frac{|(\mathrm{\Gamma }_v)|}{\left({\displaystyle \genfrac{}{}{0pt}{}{k_v}{2}}\right)}.$$
(1)
¿From the above definition, it is clear that $`\gamma `$ is a measure of the relative number of triads (fully connected subgraphs of three vertices). Note also that $`\gamma `$ is strictly in the interval $`[0,1]`$ with the upper limit attained only for a fully connected graph. In a social acquaintance network, for example, $`\gamma =1`$ if everyone in the network knows each other. It should be noted that even though the BA model successfully explains the scale-free nature of many networks, it has $`\gamma 0`$ and thus fails to describe correctly networks with the high clustering, such as social networks.

We below review briefly the BA model of the scale-free network and present our model for the scale-free network with the high clustering. The BA model SF  is defined as follows:

* Initial condition: To start with, the network consists of $`m_0`$ vertices and no edges.
* Growth: One vertex $`v`$ with $`m`$ edges is added at every time step. Time $`t`$ is identified as the number of time steps.
* Preferential attachment (PA): Each edge of $`v`$ is then attached to an existing vertex with the probability proportional to its degree, i.e. the probability for a vertex $`w`$ to be attached to $`v`$ is 111In practice one may use $`P_w=(k_w+1)/_{v𝒱}(k_v+1)`$, to make it possible for disconnected vertices to be connected. This does not change the resultant network geometry significantly.

$$P_w=\frac{k_w}{_{v𝒱}k_v}.$$
(2)

In the BA model, the growth step is then iterated $`N=|𝒱|`$ times, and for each growth step the PA step is iterated $`m`$ times for $`m`$ edges of the newly added vertex $`v`$.

In order to incorporate the high clustering we modify the above BA algorithm by adding an additional step:

* Triad formation (TF): If an edge between $`v`$ and $`w`$ was added in the previous PA step, then add one more edge from $`v`$ to a randomly chosen neighbor of $`w`$. If there remains no pair to connect, i.e., if all neighbors of $`w`$ were already connected to $`v`$, do a PA step instead.

When a vertex $`v`$ with $`m`$ edges is added to the existing network, we first perform one PA step, and then perform a TF step with the probability $`P_t`$ or a PA step with the probability $`1P_t`$. The average number $`m_t`$ of the TF trials per added vertex is then given by $`m_t=(m1)P_t`$, which we take as the control parameter in our model (see Fig. 1). It should be noted that our model reduces to the original BA model when $`m_t=0`$.

The standard scale free network model not only generates networks with certain geometrical properties, it suggests a mechanism for the emergence of power-law degree distributions in evolving networks: New actors (vertices) in a social context prefers to attach to more connected (“well known”) actors. The sociological interpretation for the triad formation step is that after being acquainted with (linked to) $`w`$ an actor $`v`$ is likely to be acquainted to $`w`$’s acquaintances as well. This mechanism of the emergence of clustering is well-known, and was discussed under the name “sibling bias” already in Ref. FRO . Recently, Ref. coauthorship\_networks  provided empirical evidence for both the mechanisms of triad formation and preferential attachment used in our construction algorithm.

The clustered scale-free network algorithm defined above gives the same degree distribution as the standard scale free network, at least if every TF step follows a PA step. To see this, first observe that in a PA step an arbitrary vertex $`v`$ increases its degree with the rate

$$\frac{\mathrm{\Delta }k_v}{\mathrm{\Delta }t}=A\frac{k_v}{_{w𝒱}k_w}\text{ for a PA step},$$
(3)
where the normalization factor $`A`$ for one edge is determined to be unity following Ref. SF . For a TF step the average increase of $`k_v`$ is proportional to the probability that a vertex in the neighborhood $`w`$ is linked in the PA step before, times the inverse of that vertex’s degree (the probability that $`v`$ is linked from $`w`$):

$$\frac{\mathrm{\Delta }k_v}{\mathrm{\Delta }t}=\frac{_{w\mathrm{\Gamma }_v}k_w(1/k_w)}{_{w𝒱}k_w}=\frac{k_v}{_{w𝒱}k_w}\text{ for a TF step},$$
(4)
where we have used the same normalization as in Eq. (3) and $`\mathrm{\Gamma }_v`$ is the neighborhood of $`v`$ (we use that the number of vertices in $`\mathrm{\Gamma }_v`$ is $`k_v`$). From Eqs. (3) and (4) the total rate for one time step, composed of $`m_t`$ TF steps and $`mm_t`$ PA steps, is expressed as

$$\frac{\mathrm{\Delta }k_v}{\mathrm{\Delta }t}=m_t\left(\frac{k_v}{_{w𝒱}k_w}\right)+(mm_t)\left(\frac{k_v}{_{w𝒱}k_w}\right)=\frac{k_v}{2t},$$
(5)
which has the same form as the original BA model and thus results in

$$k_vt^{1/2}.$$
(6)
Consequently, the degree of an arbitrary vertex increases as the square root of the time, which then yield the power-law degree distribution: $`P(k)k^3`$ SF .

In the above discussion we have assumed that a TF step always follows a PA step. If a TF step would be proceeded by another TF step the factor $`k_w(1/k_w)`$ in Eq. (4) would be replaced by $`k_w(1/(k_w1))`$ which is a small correction when $`k_w`$ is large (which it is likely to be by the definition of the PA step). And thus the resulting degree distribution would not differ much from a power-law. In Fig. 2, the degree distributions $`P(k)`$ at various values of $`m_t`$ are displayed and we find that at any value of $`m_t`$, the distribution is well described by the power law with the exponent $`a3`$ in $`P(k)k^a`$, as is expected from the above analytic consideration.

The parameter $`m_t`$ in our model introduces the clustering effect into the system by allowing the formation of triads. We only focus on the case of $`m=3`$ with expectation that other values of $`m`$ should give qualitatively the same behavior. One expects then that for any $`m_0`$ a finite $`m_t`$ gives a finite clustering coefficient $`\gamma `$ in the thermodynamic limit of $`N\mathrm{}`$, whereas for $`m_t=0`$ (the BA scale-free network model) $`\gamma `$ goes to zero as $`N`$ becomes larger. In Fig. 3(a), $`\gamma `$ at various values of $`m_t`$ is shown as a function of system size $`N`$. As expected, we find that $`\gamma `$ approaches to a finite nonzero value as $`N`$ is increased at nonzero $`m_t`$, whereas the BA model, which corresponds to the limiting case of $`m_t=0`$ in our model, is confirmed to to have $`\gamma =0`$. Furthermore, we also observe that the relation between $`m_t`$ and $`\gamma `$ is almost linear, as depicted in Fig. 3(b).

¿From the above observations, we conclude that our model exhibits both the scale-free nature and the high-clustering at the same time, while the WS model (the BA model) lacks the former (the latter) property. We note that in many real networks, both properties usually coexist, and thus believe that our model is more realistic. The triad formation step in our model, which inevitably gives a high clustering coefficient, is expected to make the average geodesic length smaller than the BA network, since the edge for the triad could have been used to connect two vertices separated by a large distance if only the preferential attachment step was allowed. However, the characteristic path length, defined as the average of the geodesic length, $`\mathrm{}`$, is found to behave logarithmically with the size $`N`$, the same behavior as the WS model and the BA model. In Fig. 4, we present $`\mathrm{}`$ versus $`N`$ at various values of $`m_t`$. It is shown that $`\mathrm{}`$ becomes larger as $`m_t`$ is increased, as expected. Furthermore Fig. 4 shows that the increase of $`\mathrm{}`$ is logarithmic for all $`m_t`$.

By mimicking principles in network formation, a generation algorithm can construct graphs with certain topological statistics, such as a degree distribution, clustering coefficient, and so on. However, it should be emphasized that these kinds of algorithms cannot claim to uniformly sample the ensemble of networks with specific statistical properties. This drawback exists even in more general classes of random graphs where structural biases, such as clustering, are imposed. markov ; Str

Recently, Klemm and Equíluz klemm  have proposed a network model based on a finite memory of vertices, i.e., vertices become inactive and do not get new edges after a finite number of time steps, and have shown that their growth and deactivation model exhibits both the high clustering and the scale-free nature. Our model provides an alternative possibility to achieve the same feature, the clustered scale-free nature, based on our frequent everyday experience on how we are acquainted by newcomers: $`B`$ becomes $`A`$’s new friend since $`B`$ is introduced by one of $`A`$’s friends. Even in the network of scientific citations, it is likely that authors of paper $`A`$ refer paper $`B`$ since they have found $`B`$ when they read a famous review paper $`C`$ citation\_networks . This then has close resemblance to our model, the TF step accompanied by the PA step. In Ref. puniyani , a model with both the high clustering and the scale-free distribution has also been suggested. However, the power-law degree distribution was assigned to the network to start with, and the next following steps were devised not to change the degree at each vertex. In other words, the power-law distribution in Ref. puniyani  was not an emerging property in the model, which is different from the BA model as well as our model in this work. Very recently, we have learned about the work by Davidsen et aldavidsen , which is based on the same observation of triad formation as ours and has been shown to possess similar network properties, i.e. the high clustering, small average geodesic length, and a scale-free distribution. We believe, however, that our model has some advantage in describing networks which grow in time, whereas the network model in Ref. davidsen  has fixed network size.

In conclusion, we have proposed an algorithm for generation of growing networks with power-law degree distribution, a logarithmic increase of the average geodesic length, and a finite clustering. The last two properties make the generated graphs qualify as a small-world network in the Watts and Strogatz sense, in addition to their scale-freeness. The simple relation between the coefficient $`m_t`$ and $`\gamma `$ further increases the usefulness of the suggested algorithm, making it possible to tune the clustering coefficient in a systematic way.

###### Acknowledgements.

This work was supported in part by the Swedish Natural Research Council through Contracts Nos. F 5102-659/2001 and E 5106-1643/1999.


text:
# Violations of local equilibrium and linear response in classical lattice theories

## I Introduction

Non-equilibrium dynamics appears in almost all areas of physics, from physics of the early universe, heavy ion collisions to transport in matter. In many of these situations, we would like to find out how physical quantities behave in various quantum field theories under non-equilibrium conditions. In an ideal world, we would like to analyze their dynamics from first principles, with no approximations. Let us first discuss what this entails: We would use no dynamical assumptions, the computations will be fully non-perturbative, transport properties will not be restricted only to linear response and the non-equilibrium properties of the system arise dynamically if some non-equilibrium conditions are imposed at the boundaries. These are some of the requirements that need to be satisfied.

In practice, this is too ambitious. Usually, we resort to computations of the Green–Kubo formulas for the transport coefficients or use the Boltzmann equation, which is truncated from the full BBGKY hierarchy. The equations are solved essentially in equilibrium so that the region of validity of the results are unclear, not to mention that the results can not be applied to situations far from equilibrium. Even then, further approximations need to be made. For instance, the transport coefficients of $`\varphi ^4`$ theory to leading order in the coupling was completed only recently . The reason for this is clear; the computations are quite involved and results can not be obtained in closed form, except for its asymptotic behaviors.

In this work, we study the non-equilibrium dynamics of two types of theories, $`\varphi ^4`$ theory and the FPU $`\beta `$ model, in $`D=1,2,3`$ spatial dimensions. We shall work with classical theories on the lattice and from first principles. In particular, we shall make no  dynamical assumptions and our computations will be fully non-perturbative. While the $`\varphi ^4`$ theory has a bulk limit, the FPU $`\beta `$ model does not, and the temperature dependence of the transport coefficients qualitatively differs. It hardly needs to be mentioned, but the $`\varphi ^4`$ theory is a classic prototypical field theory and the FPU $`\beta `$ model is another classic model which has been studied widely. These distinct theories in various dimensions should allow us to gauge how applicable the results are to more general theories. Putting these systems under weak and strong thermal gradients, we shall study how we can measure deviations from local equilibrium and analyze its behavior quantitatively. We also examine the validity of the linear response law and analyze its violations quantitatively. Further, we study the relations of these violations to the thermodynamic properties of the systems, such as the equations of state. We then discuss how to extend these results to quantum field theories.

The main reason we can analyze the systems from first principles is because we are working with classical systems on a finite lattice (of varying sizes). In the classical theory, the problem essentially reduces to computing the solutions to the equations of motion for the degrees of freedom of the system. We can perform this task without any approximations such as perturbation theory by numerically integrating the equations of motion. While the behavior of a classical system is different from that of a quantum system, the latter often reduces to the former under certain conditions, such as for high temperatures. Also, the classical theory can play a major role in obtaining the results in the full quantum theory, as can be seen from the results in the electroweak theory at finite temperatures . Furthermore, we believe that one should understand the behavior of classical systems first, before we can hope to understand the full quantum theory. We work with lattice theories which we can think of as a real lattice as in solids. We can also think of the lattice as a discrete approximation of the continuum theory, which is quite natural and corresponds to a theory with a cutoff. Of course, it hardly needs to be mentioned that the behavior of classical lattice systems is of interest on its own right!

## II The systems

The Hamiltonians of the models we study are

$$H=\underset{𝐫}{}\left[\frac{1}{2}\pi _𝐫^2+\frac{1}{2}(\varphi )_𝐫^2+V(\varphi _𝐫)\right],V_{\varphi ^4}(\varphi _𝐫)=\frac{\varphi _𝐫^4}{4},V_{FPU}(\varphi _𝐫)=\frac{(\varphi )_𝐫^4}{4}$$
(1)
in $`1`$$`3`$ spatial dimensions. The sum is taken over all lattice points $`𝐫`$ and $`\varphi _𝐫`$ is the lattice gradient. As should be clear, the $`\varphi ^4`$ theory can be thought of as a discretized version of the standard continuum $`\varphi ^4`$ theory. The mass has been set to zero for simplicity here, although we have also studied cases where the mass is non-zero.

In this work, we are interested in the steady state properties of the systems under thermal gradients. As such, we apply thermostats at the boundaries $`x=0,L`$ and impose periodic boundary conditions for the other directions for $`D=2,3`$. It is important to note that the thermostats are applied only at the ends so that in the interior, $`0<x<L`$, the degrees of freedom are solely  those of the $`\varphi ^4`$ theory or the FPU model. The thermostats we use are generalized variants of the thermostats of Nosé and Hoover. Explicit forms of equations of motion including the thermostat degrees of freedom can be found in . The thermostats are deterministic, rather than stochastic, and the thermostats effectively impose thermal distributions at the thermostatted sites when time averaged. A physical observable, $`𝒪`$, such as correlation functions or currents, is computed by sampling them along the classical trajectory of the whole system in the phase space and taking its time average, $`𝒪`$. We studied lattice sizes in the $`x`$ direction up to $`L=10000`$ for $`D=1`$ and up to $`L=2000`$ and transverse lattice sizes of $`L_{}=320`$ when $`D=2,3`$.

Since we will be interested in thermal transport, we shall need an expression for the energy flow,

$$𝒯_{\varphi ^4}^{0x}=\pi _𝐫_x\varphi _𝐫,𝒯_{FPU}^{0x}=\left(\pi _𝐫_x\varphi _𝐫\right)\left[1+(\varphi )_𝐫^2\right]$$
(2)
which satisfies the usual continuity equation.

## III Near equilibrium behavior

When the two end point temperatures, are equal, $`T_1=T_2`$, the interior of the system thermalizes at the same temperature. This can be explicitly confirmed by checking that the momentum distribution is Maxwellian. The local temperature is defined by $`T(𝐫)=\pi _𝐫^2`$, which is called the ideal gas temperature. This is well defined and is a robust local measure, since the Hamiltonians in (1) are quadratic in $`\pi _𝐫`$ and $`\pi _𝐫`$ interact with the neighboring sites only indirectly through $`\varphi _𝐫`$.

When we make the two boundary temperatures slightly different, we recover a straight temperature profile which can be well described by the Fourier’s law,

$$𝒯^{0x}=\kappa (T)_xT$$
(3)
as can be seen in Fig. 1 (left).While it will not be discussed here, in general, temperature jumps will arise at the boundaries which can be understood. By choosing various small temperature differences around the same central temperature, we directly extract the thermal conductivity $`\kappa `$ for a given temperature in each system at each lattice size from measurements of the energy flow and the temperature gradient, using (3). In the $`\varphi ^4`$ theory, we find a bulk limit exists for $`\kappa `$, namely the limit $`\kappa (L\mathrm{})`$ exists. In the FPU model, the bulk limit does not  exist. The results can be summarized for the $`\varphi ^4`$ theory as (see Fig. 2 (left))

$$\kappa (T)=\frac{A}{T^\gamma },\begin{array}{ccc}& A& \gamma \\ D=1& 2.8(1)& 1.32(2)\\ D=2& 5.2(4)& 1.42(6)\\ D=3& 8.9(5)& 1.75(5)\end{array}$$
(4)
For the FPU model, the behavior is more complicated as can be seen in Fig. 2 (right); $`\kappa `$ depends on L and $`T`$ dependence is no longer a simple power. The $`L`$ dependence is compatible with a simple power law and the asymptotic behavior

$$\kappa \{\begin{array}{cc}1.2L^\delta T^1\hfill & (T0.1)\hfill \\ 2L^\delta T^{1/4}\hfill & (T50)\hfill \end{array},\delta =0.37(3).$$
(5)
can be understood from scaling arguments.

We can also compute the thermal conductivity independently of the Fourier’s law (3) through the integral of the autocorrelation of the energy current using the Green–Kubo formula. The results agree with the direct computation from Fourier’s law. For $`D=1`$ $`\varphi ^4`$ theory, the Green–Kubo results are shown in Fig. 2 (left) along with those from direct computations. We note in passing that the applicability of the Green–Kubo formalism is regarded as being most subtle for lower dimensional systems due to the existence of persistent correlations.

As we make the temperature gradients larger, the temperature profiles become visibly non-linear, as in Fig. 1 (right). This does not necessarily signal the breakdown of linear response.§§§By “linear response”, we refer strictly to the validity of (3). In particular, we do not presume that the existence of a bulk limit.Since $`\kappa `$ depends on the temperature (cf. Fig. 2), as the temperature changes, so does the gradient, and it can do so within the same system. This reasoning leads to the following formula for the temperature profile by applying the linear response law (3) locally.

$$T(x)=T_1\left[1\frac{x}{L}+\left(\frac{T_2}{T_1}\right)^{1\gamma }\frac{x}{L}\right]^{\frac{1}{1\gamma }}$$
(6)
Here, we assumed that the temperature dependence of the thermal conductivity can be described by $`\kappa (T)=AT^\gamma `$, which is true for the $`\varphi ^4`$ theory and for most temperature regions in the FPU model. It should be noted that the profile has the scaling behavior that it depends on $`x`$ only through $`x/L`$. This formula for the profile works quite well when the we are not too far from equilibrium. For example in Fig. 1 (right), the theoretical curve (6) along with the profile is presented and we see that it provides an excellent description of the curved temperature profile for not so large gradients. We shall turn to the question of quantifying how far we are from equilibrium and how good the description is, in the next section.

## IV Breakdown of local equilibrium and linear response

When the system moves further away from equilibrium, the formula (6) is no longer adequate to describe the temperature profile, as in Fig. 1 (right). One reaction might be to employ some kind of ‘non-linear’ response theory. However, we would first like to understand the overall behavior of the system quantitatively and judge whether such a description is valid. To this end, we would like to also assess if local equilibrium is satisfied so that for instance, the temperature is well defined. For this, we need an explicit measure of how well the local equilibrium holds, which we find in the cumulants of $`\pi _𝐫`$. Since the Hamiltonian of the system (1) is quadratic in $`\pi _𝐫`$, the thermal distribution should be Maxwellian if local equilibrium holds. The cumulants provide an unambiguous local measure of local equilibrium breakdown. We shall also need a local measure of how “far” we are from equilibrium and for this, we use $`_xT/T`$. We note that it is natural to use the rescaled gradient here. This is quite evident in Fig. 1 (right); the linear response prediction is less applicable in the small $`T`$ region, wherein $`_xT`$ is smaller but $`_xT/T`$ is larger.

Before proceeding further, we would like to discuss the general picture: We find that a physical observable, $`𝒜`$, generically deviates from its local equilibrium value in a manner

$$\delta _𝒜\frac{\delta 𝒜}{𝒜}=C_𝒜\left(\frac{_xT}{T}\right)^2+C_𝒜^{}\left(\frac{_xT}{T}\right)^4+\mathrm{}$$
(7)
While this seems quite natural since the intrinsic behavior should not depend on the direction of the gradient, it has not been established analytically. Furthermore, the behavior of the coefficients is more subtle than what one would expect, as we shall see below.

Let us analyze first the violations of local equilibrium and linear response. For deviations from local equilibrium, we use the observable, $`\pi ^4/T^2`$ where the fourth cumulant is $`\pi ^4\pi ^43\pi ^2^2`$. Likewise, for linear response, we use the observable $`(𝒯^{0x}\kappa _xT)/\kappa _xT`$. We show some of the results in Fig. 3. These computations need to be performed with care since we are looking at differences.

From these analyses, we extract the coefficients $`C_{LE},C_{LR}`$ in (7) for local equilibrium and linear response, respectively. $`C_{LE},C_{LR}`$ can in general depend on $`T`$ and $`L`$. We have computed these coefficients at various $`T`$ and $`L`$ for both the $`\varphi ^4`$ theory and the FPU model in various dimensions. While we do not have yet a global picture of how these coefficients behave, a few comments are in order. For the $`\varphi ^4`$ theory, which has a bulk limit, we would naively expect that these coefficients will be independent of $`L`$. However, this turns out to be not the case. For instance, at $`T=1`$, $`C_{LE}=aL`$ with $`a=3,2,4`$ for $`D=1,2,3`$. There is thus a strong dependence on the size of the system, so that a simple local understanding of (7) does not seem possible. In general, both for the $`\varphi ^4`$ theory and for the FPU model in $`D=1,2,3`$ and for various $`L`$, we find that $`\delta _{LE}\delta _{LR}`$. This means, in particular, that a simple non-linear response theory is not applicable to these systems since the breakdown of local equilibrium also needs to be considered simultaneously. We note that a priori, it could have been such that $`\delta _{LE}`$ was much smaller that $`\delta _{LR}`$, in which case a non-linear response theory would have seemed quite appropriate. Generically, we find the dependence of $`C_{LE},C_{LR}`$ on $`L`$ to be close to linear, and the temperature dependence to be weak, which supports the rescaled gradient $`_xT/T`$ as a natural measure of how far we are from equilibrium. We have further computed the deviations in the equations of state for pressure and local energy density from local equilibrium relations. We find that the deviations can again be understood using the description in (7).

## V Towards quantum field theory in non-equilibrium

As we have seen in the previous sections, in classical lattice systems, we can ask practically any question regarding the physics behavior of non-equilibrium systems and obtain an answer, at least under steady state conditions. We are not restricted to weak coupling, nor do we have to assume linear response theory and we can be arbitrarily far from equilibrium. Clearly, we want to extend this situation to quantum field theories. Ideally, we would like to solve quantum dynamics in a similar fashion, leading to the solution of quantum field theory in non-equilibrium, in the continuum limit. This is, of course, too ambitious at this point; even for moderately large systems (say $`100`$ degrees of freedom), it is computationally prohibitive to carry out such a program. Our goal here will be more modest; we would first like to apply the results of the classical theory to quantum field theory and find out how they fit in. We shall see that even this task is quite non-trivial. Typically, classical results are applicable to quantum theories in some region or in some limit . At first sight, the prospect seems quite promising for the lattice $`\varphi ^4`$ theory; amongst other features, the theory has a bulk limit with the behavior appearing already for system sizes of order hundred or less for moderate temperatures and the temperature profile (6) has a smooth continuum limit.

However, in practice, we immediately encounter an obvious problem: when one tries to take the naive continuum limit by taking the lattice spacing to zero, quantities such as the energy density, thermal conductivity diverge. They diverge for the trivial reason that the system will have an infinite number of oscillators in the unit volume.This is the “Rayleigh–Jeans” problem, known from the early days of the 20th century which afflicts all classical theories in the continuum. While it might seem that this will make it impossible to obtain finite results in any classical continuum theory at finite temperatures, it has been demonstrated that one can obtain finite results for physical observables in the classical continuum theory that matches with the quantum results in the appropriate region. This is done by renormalizing the theory classically using the appropriate matching conditions.

Next point to be resolved before we can start applying our results to quantum field theory is the continuum limit. How far are we from the continuum limit and how do we know? While we do not have a definitive answer to this difficult question, we can study the physical mass in lattice units, $`m_{ph}a`$, which effectively measures how big our lattice mesh is in physical units. (The lattice size, $`a`$, has been set to one here.) To avoid lattice artifacts, we would like this quantity to be small. The physical mass can be measured from the correlation function as

$$\varphi (0)\underset{𝐫_{}}{}\varphi (x,𝐫_{})e^{m_{ph}|x|}$$
(8)
where $`𝐫_{}`$ denotes the coordinates transverse to $`x`$. In Fig. 4 (left), we study the correlation function averaged in the transverse direction. We see that the correlation function indeed does decay exponentially, which allows us to extract the physical mass. We find, as in Fig. 4 (right), that the physical mass has a temperature dependence which we find to be

$$m_{ph}=\{\begin{array}{cc}0.97(2)T^{0.30(1)}\hfill & D=1\hfill \\ 0.69(5)T^{0.44(3)}\hfill & D=3\hfill \end{array}$$
(9)
To put things in perspective, in lattice QCD, the state of the art simulations are typically run with $`m_{ph}a=0.10.5`$. So, for $`T1`$, we should be able to overcome lattice artifacts.

One more aspect needs to be addressed: In the continuum limit, when the lattice no longer exists and the theory becomes relativistic, we have to decide relative to what reference frame we have transport in the bulk. Here, a natural situation is provided if we consider a system with more conserved charges. For instance, it is quite straightforward to extend our analysis to a complex scalar field theory with U$`(1)`$ symmetry. In this case, we can address energy transport in the system with no net charge flow, which would allow us to extract the thermal conductivity in the usual sense. Since our methods can be used in practically any classical lattice system, we should be able to address the problem of relativistic transport, with judicious choice of models.

While more work is definitely needed to apply our results to the appropriate regimes of quantum field theory, none of the obstacles seems insurmountable. The lattice artifacts should be reasonably small for $`T1`$, the divergences of the classical continuum theory can be overcome and the relativistic limit is quite compatible.

## VI Discussions

In summary, we have been able to obtain the physical behavior of $`\varphi ^4`$ theory and the FPU $`\beta `$ model under thermal gradients, near and far from equilibrium. Close to equilibrium, we see that linear response theory works well and local equilibrium description is applicable. As we move away from equilibrium, linear response breaks down, but so does local equilibrium. We have measured quantitatively how much they deviate from equilibrium values and they are rather similar in all the systems we studied. The deviation of these quantities from the equilibrium behavior is quadratic in the thermal gradient, as in (7). This behavior is rather generic in the systems we studied, but no analytic explanation of this behavior is known. It is interesting to note that in the non-linear dependence of the viscosity on the shear rate, it was initially argued analytically and further found numerically that the behavior is not  quadratic but rather has a non-analytic behavior with a power of $`3/2`$. Research in this area is still in progress and the question whether the dependence is analytic still seems quite open .

We have further discussed how to apply our results to quantum field theory. Obviously, our results do not incorporate essentially quantum behavior so that they can be applied only to certain regimes in the quantum theory. Even though the task is non-trivial, we found that the program is quite feasible.

Much more needs to be done: While we can answer practically any question thrown at us regarding the steady state behavior of classical lattice systems under thermal gradients, the results by themselves lack analytic understanding. We have tried to fill in the missing links by analyzing these results, but we feel that a deeper understanding of the phenomena and their relation to the intrinsic non-equilibrium dynamics is highly desirable. For instance, as we mentioned above, even the quadratic behavior of the deviations from equilibrium is something that has not been established analytically in any model. While this behavior might be intuitive and understandable, the fact that we have the breakdown of local behavior in these deviations is quite subtle and needs to be explained. In applying the classical results to the quantum theory, we still need to complete the program. Furthermore, if one wanted to perform even moderately large scale quantum simulations from first principles, how to approach it is still unclear. Also, even in classical systems, we have not dealt with transient phenomena, which is an interesting and important avenue of study. We analyzed the steady state physics as a first step, but many non-equilibrium phenomena, such as the theory of the early universe or heavy ion collisions, are intrinsically time dependent. We feel that there is still much more interesting physics to be uncovered in the area of non-equilibrium systems, even classically.

Acknowledgements: We would like to thank Gert Aarts, Jan Smit and Larry Yaffe for stimulating discussions regarding the material in the last section.


text:
# On the gravitational instability of a dissipative medium.

## I Introduction

Structure formation in fluid systems is closely related to the existence of growing modes while considering density fluctuations in the hydrodynamical regime. Relevant growing modes at very small wave numbers have long ago been identified in gravitational systems, leading to key concepts in present astrophysics such as the Jeans wave number and the Jeans mass . These parameters correspond to the conditions that a density fluctuation must satisfy in order to grow up to the point of forming a structure.

Entropy production plays an important role in realistic descriptions of thermodynamical processes, and some discussion of its physical sources may yield interesting results in the analysis of the time evolution of statistical fluctuations in astrophysical systems. This work is mainly devoted to the study of the dispersion relation that arises from the fluctuation analysis of the thermodynamical variables of a simple dissipative fluid under the influence of a gravitational field. The analysis is reviewed in the case of non-vanishing viscosity and different features of the fluid are established. In particular, this work shows the possibility of viscosity-affected growing modes for some astrophysical fluids, which may be of interest for both astrophysicists and cosmologists.

This paper is divided as follows: section two is dedicated to the hydrodynamical model analysis and the establishment of the dispersion relation in which the main results are based. Section three is devoted to the mathematical treatment needed for its solution. The Jeans instability is recovered in the presence of viscosity for typical situations present in protogalaxies. Finally, a discussion of the results obtained and an outline of future work are included in section four of the paper.

## II Hydrodynamic model

In order to establish the basic equations that lie at the core of this work, we remind the reader that the so-called Navier-Stokes-Fourier equations of hydrodynamics for a simple fluid arise from the structure of what is now called Linear Irreversible Thermodynamics (LIT) by supplementing the two conservation equations for mass and momentum, respectively and the balance equation for the internal energy, with additional information. Indeed, these equations read .

$$\frac{\rho }{t}+\frac{}{x^i}(\rho u^i)=0$$
(1)
$$\rho \frac{Du^i}{Dt}+\frac{\mathrm{\Xi }^{ij}}{x^j}=F^i$$
(2)
$$\rho \frac{D\epsilon }{Dt}+\frac{Q^j}{x^j}=\mathrm{\Xi }^{ij}\frac{u_i}{x^j}$$
(3)
where $`\frac{D}{Dt}=\frac{}{t}+u^i\frac{}{x^i}.`$

Here, $`\rho (x^j,t),u^i(x^j,t)`$ and $`\epsilon (x^j,t)`$ are the local density, velocity and internal energy, respectively, $`\mathrm{\Xi }^{ij}`$ the momentum current (or stress tensor) and $`Q^j`$ the heat flux. All indices run from 1 to 3. In general, for isotropic fluids, $`\mathrm{\Xi }^{ij}=p\delta ^{ij}+\tau ^{ij}`$, where $`p`$ is the local hydrostatic pressure and $`\tau ^{ij}`$ the viscous tensor ($`\delta ^{ij}`$ is the unitary dyadic). Notice that Eqs. (1-3) contain fifteen unknowns (including the pressure) and there are only five equations, so the system is not well determined. If we arbitrarily choose to describe the states of the fluid through the set of variables $`\rho (x^j,t),u^i(x^j,t),`$ $`T(x^j,t),`$ where $`T`$ is the local temperature, we need nine dynamic equations of state (or constitutive equations) relating $`\mathrm{\Xi }^{ij}`$ and $`Q^j`$ to the state variables plus two local equations of state $`p=p(\rho ,T)`$ and $`\epsilon =\epsilon (\rho ,T).`$ If according to the tenets of LIT we choose the so-called linear constitutive laws, namely

$$\tau ^{ij}=2\eta \sigma ^{ij}\zeta \theta \delta ^{ij}$$
(4)
$$Q^j=\kappa \delta ^{ji}\frac{T}{x^i}$$
(5)
with $`\eta `$ and $`\zeta `$, the shear and bulk viscosities, respectively, $`\sigma ^{ij}`$ the symmetrical traceless part of the velocity gradient, $`\theta =\frac{u^i}{x^i}`$ and $`\kappa `$ being the thermal conductivity.

Eqs. (4) and (5) are the well known constitutive equations of Navier-Stokes and Fourier, respectively. Substitution of these equations into Eqs. (1-3) yields a set of second order in space, first order in time non-linear coupled set of partial differential equations for the chosen variables $`\rho ,u^i`$ and $`T`$. This set, which the reader may seek in the literature , , is the so called Navier-Stokes-Fourier system of hydrodynamic equations.The non-linearities appearing in such equations have two sources, the inertial terms $`u^i\frac{}{x^i}`$ arising from the hydrodynamic time derivatives, plus quadratic terms in the gradients of velocity arising from Eqs. (4) and (5). Moreover, it should be mentioned that this set of equations is consistent with the second law of thermodynamics, the Clausius uncompensated heat , or entropy production, is strictly positive definite.

Nevertheless, for the purpose of this paper, this set of equations is too complicated. In order to deal with fluctuations around the equilibrium state, one assures that for any of two state variables, call them $`X(x^j,t)`$ one can write that,

$$X(x^j,t)=X_o+\delta X(x^j,t)$$
(6)
where $`X_o`$ is the equilibrium value of $`X`$ and $`\delta X`$ the corresponding fluctuation. Neglecting all terms of order $`(\delta X(x^j,t))^2`$ and higher in the NSF non-linear set one finds the linearized NSF equations of hydrodynamics -,

$$\frac{}{t}\left(\delta \rho \right)+\rho \theta =0$$
(7)
$$\rho _o\frac{u_k}{t}=\frac{1}{\rho _o\kappa _T}\frac{}{x^k}\left(\delta \rho \right)\frac{\beta }{\kappa _T}\frac{}{x^k}\left(\delta T\right)+2\eta \delta ^{ij}\frac{u_k}{x^ix^j}(\frac{2}{3}\eta \varsigma )\frac{}{x^k}(\theta )+F_k$$
(8)
$$\frac{}{t}\left(\delta T\right)=D_T\delta ^{ij}\frac{\left(\delta T\right)}{x^ix^j}\frac{\beta T_o}{\rho _oC_p\kappa _T}\theta $$
(9)
since $`u_{ko}=0`$, $`u_k=\delta u_k`$ and $`\frac{\beta }{\kappa _T}=(\frac{p_0}{T_o})_{\rho _o}`$, $`\kappa _T=\rho _o(\frac{p_0}{\rho _o})_{T_o}`$ and $`D_T=\frac{k}{\rho _oC_v}`$ is the thermal diffusivity. $`C_p`$ and $`C_v`$ are the specific heats at constant pressure and constant volume, respectively.

Note now that from Eq. (8), $`u_{kTrans}Rot\left(u_k\right)`$ uncouples from the hydrodynamic modes, whence, one formally arrives at the desired set, namely:

$$\frac{}{t}\left(\delta \rho \right)+\rho _o\theta =0$$
(10)
$$\rho _o\frac{\theta }{t}=\frac{1}{\rho _o\kappa _T}\delta ^{ij}\frac{\left(\delta \rho \right)}{x^ix^j}\frac{\beta }{\kappa _T}\delta ^{ij}\frac{\left(\delta T\right)}{x^ix^j}+D_v\delta ^{ij}\frac{\left(\theta \right)}{x^ix^j}\delta ^{ij}\frac{\left(\phi \right)}{x^ix^j}$$
(11)
$$\frac{}{t}\left(\delta T\right)=D_T\delta ^{ij}\frac{\left(\delta T\right)}{x^ix^j}\frac{\beta T_o}{\rho _oC_p\kappa _T}\theta $$
(12)
assuming that $`F_k=\frac{\phi }{x^k}`$ . Also, $`D_v=(\frac{4}{3}\eta +\varsigma )\frac{1}{\rho _o}`$. Now, we notice that $`\frac{}{t}\left(\delta \rho \right)=\rho _o\theta `$ so we can reduce this set to only two equations, Eq. (11) and Eq. ( 12). We now perform a few minor transformations introducing the speed of sound of the fluid $`C_o^2`$ $`=\left(\frac{p}{\rho }\right)_s`$ through the relationship

$$k_T=\gamma k_s=\frac{C_p}{C_v}k_s$$
(13)
where $`k_s`$ is the adiabatic compressibility $`\left(\frac{\rho }{p}\right)_s`$ whence, $`k_T=\frac{\gamma }{C_o^2}`$. This leaves us finally with the set:

$$\frac{^2(\delta \rho )}{t^2}+\frac{C_o^2}{\gamma }\delta ^{ij}\frac{\left(\delta \rho \right)}{x^ix^j}+\frac{C_o^2\beta }{\gamma }\delta ^{ij}\frac{\left(\delta T\right)}{x^ix^j}+D_v\delta ^{ij}\frac{}{x^ix^j}(\frac{}{t}\left(\delta \rho \right))+\rho _o\frac{\left(\phi \right)}{x^ix^j}=0$$
(14)
$$\frac{}{t}\left(\delta T\right)=D_T\delta ^{ij}\frac{\left(\delta T\right)}{x^ix^j}\frac{\gamma 1}{\beta }\frac{}{t}\left(\delta \rho \right)=0$$
(15)
where we have used that $`C_pC_v=\frac{\beta ^2T_o}{\rho _o^2\kappa _T}`$.

Eqs.(14-15) form a set of coupled equations for the density and temperature fluctuations in the fluid under the action of a conservative force whose nature need not to be specified for the time being. They are the basis for studying the properties of the time correlation functions of thermodynamic fluctuations. Those of the density will be of particular interest here.

## III Solution to the hydrodynamic equations

The results derived in the previous section are far from being new. Aside from the term $`^2\varphi `$ which arises from the presence of an external conservative force, they are identical to the ones that have been widely discussed in the literature. The question here is if the gravitational potential introduces any substantial modification in the correlation functions for the thermodynamic fluctuations. To examine this possibility we recall that, if we are considering only the fluctuations, the gravitational potential satisfies the Poisson equation:

$$\delta ^{ij}\frac{\left(\phi \right)}{x^ix^j}=4\pi G\delta \rho $$
(16)
Now, equation (15) reads

$$\frac{^2(\delta \rho )}{t^2}+\frac{C_o^2}{\gamma }\delta ^{ij}\frac{\left(\delta \rho \right)}{x^ix^j}+\frac{C_o^2\beta }{\gamma }\delta ^{ij}\frac{\left(\delta T\right)}{x^ix^j}+D_v\delta ^{ij}\frac{}{x^ix^j}(\frac{}{t}\left(\delta \rho \right))4\pi G\rho _o\delta \rho =0$$
(17)
The solution to Eqs.(15) and (17) proceeds in the standard fashion. We reduce them to a set of algebraic equations by taking their Laplace-Fourier transform, choose to set the static temperature fluctuations equal to zero, eliminate the temperature leading to an equation for $`\delta \rho (\stackrel{}{k},s)`$, which is the ratio of two polynomials in $`s`$ . To compute $`\delta \rho (\stackrel{}{k},t)`$ one must take the inverse Laplace transform of the former quantity, which demands the knowledge of the roots of the denominator, which is a cubic equation in $`s^3`$ (dispersion equation). Next, we proceed with intuition, we assume that the Rayleigh peak is nos substantially modified by the external potential, so that according to the standard theory, one of the roots of the cubic equation is $`s=D_Tk^2,`$ where $`D_T`$ is the thermal diffusivity. This allows us to simplify the denominator to the form,

$$(s+D_Tk^2)(s^2+D_vk^2s+C_o^2k^24\pi G\rho _o)=0$$
(18)
the two roots of the quadratic equation are:

$$s_{1,2}=\frac{D_vk^2}{2}\pm i\left[\left(C_o^2k^24\pi G\rho _o\right)\frac{D_v^2k^4}{4}\right]^{1/2}$$
(19)
Iff $`4\pi G\rho _o=0`$ and viscosity dominates over the term in $`k^2`$, this result reduces to the one giving rise to the standard Brillouin peaks which correspond to density fluctuations of the type -,

$$\delta \rho (\stackrel{}{k},t)=\delta \rho (\stackrel{}{k},0)\frac{1}{\gamma }e^{D_vk^2t}Cos\left[C_okt\right]$$
(20)
which are the acoustic modes damped by the Stokes-Kirchhoff factor $`D_v.`$

On the other hand, if $`4\pi G\rho _o0`$, the threshold value for $`k`$ distinguishing between damped oscillations and growing modes is given by

$$\left(C_o^2k^24\pi G\rho _o\right)\frac{D_v^2k^4}{4}=0$$
(21)
or,

$$k^2=\frac{2C_o^2}{D_v^2}\left(1\pm \sqrt{1\frac{4\pi G\rho _oD_v^2}{C_o^4}}\right)$$
(22)
Eq. (22) is a generalization of Jeans wave number when dissipative effects due to viscosity are non-negligible, and is the main result of the paper. We can note that, if $`\frac{4\pi G\rho _oD_v^2}{C_o^4}1,`$ and taking the $``$ sign for the square root, we have:

$$k^2\frac{2C_o^2}{D_v^2}\left[11+\frac{1}{2}\left(\frac{4\pi G\rho _oD_v^2}{C_o^4}\right)\right]$$
(23)
or

$$k^2\frac{4\pi G\rho _o}{C_o^2}=K_J^2$$
(24)
This is the square value of the Jeans wave number  . It was derived by Jeans in 1902 and rederived in many other waves by several authors . Here we simply show that it is almost a trivial consequence of LIT.

Let us expand the square root in Eq. (22) up to second order in order to have a better grasp of the effects of viscosity on the Jeans wave number. Now we have that,

$$k^2\frac{2C_o^2}{D_v^2}\frac{2C_o^2}{D_v^2}\left[1\frac{2\pi G\rho _oD_v^2}{C_o^4}\frac{1}{8}\left(\frac{4\pi G\rho _oD_v^2}{C_o^4}\right)^2\right]$$
(25)
or, using Eq. (24):

$$k^2K_J^2+\frac{4D_v^2}{C_o^6}\left(\pi G\rho _o\right)^2=K_J^2(1+\frac{2D_v^2}{C_o^4}\pi G\rho _o)$$
(26)
so that the dimensionless number

$$\epsilon =\frac{2D_v^2}{C_o^4}\pi G\rho _o=\frac{2(\frac{4}{3}\eta +\varsigma )^2}{\rho _oC_o^4}\pi G$$
(27)
is a good measure of the second order corrections to the Jeans wave number due to viscous effects.

For the corrected Jeans mass we have :

$$M_c=\frac{4}{3}\pi \left(\frac{\pi }{k}\right)^3=\frac{4}{3}\pi \left(\frac{\pi }{K_J^{}(1+\frac{2D_v^2}{C_o^4}\pi G\rho _o)^{\frac{1}{2}}}\right)^3\rho _o=M_J(1+\frac{2D_v^2}{C_o^4}\pi G\rho _o)^{\frac{3}{2}}$$
(28)
or

$$M_c=M_J(1+\epsilon )^{\frac{3}{2}}$$
(29)
where the standard Jeans mass $`M_J`$ is given by:

$$M_J=\frac{4}{3}\pi \left(\frac{\pi }{K_J^{}}\right)^3\rho _o=\frac{4}{3}\pi \left(\frac{\pi }{\sqrt{4\pi G\rho _o}}C_o\right)^3\rho _o=\frac{\left(\pi \right)^{\frac{5}{2}}}{6\left(\left(\rho _o\right)^{\frac{1}{2}}\left(G\right)^{\frac{3}{2}}\right)}C_o^3$$
(30)
One should compare these results with Weinberg  and Peebles p. 108 . The question here is if in the primeval Universe, the plasma had a substantial viscosity.

Eq. (30) deserves more attention. As shown by Weinberg  a similar result follows from a more intuitive calculation addressing the behavior of the dispersion relation in sound waves. Using numerical values for the density of barions (protons) and the velocity of sound in a background of photons and protons, he is able to estimate the temperature range for which the protogalactic mass can grow to $`10^{11}`$ solar masses. Similar reasoning may be applied to Eq. (29) and, of course, predict a new $`M_J`$. The question here is if the viscous corrections have any real meaning at all or they just turn out to be an academic curiosity. For instance, in accretion disks  $`\rho _o4\times 10^4\frac{kg}{m^3}`$, $`\eta `$ can be estimated as high as $`10^{10}Pas`$and $`C_o5\times 10^3\frac{m}{s}`$ . This would yield a value of $`D_v10^{14}\frac{m^2}{s},`$ so that

$$\epsilon =\frac{2D_v^2}{C_o^4}\pi G\rho _o10^1$$
(31)
the estimated correction is at the order of 10%.

## IV Final Remarks

Dissipation is always present in real processes. Whenever the thermodynamical functions of a self-gravitating viscous system fluctuate around an equilibrium value, dissipation will affect the general equilibrium condition for the stability of those fluctuations. We have been able to establish a simple expression, Eq. (27) that measures how important this correction will be, assuming that the Rayleigh peak remains essentially unaffected by the the gravitational field. If this were not the case, then a full treatment of the dispersion relation would be needed, involving interesting effects of the gravitational field in the Rayleigh-Brillouin spectrum, for the case of damped systems. This modifications may have some importance in astrophysics, since natural astrophysical masers are already well-known  and its scattering in fairly dense systems will provide information about the transport properties of the scatterer.

We thank Ana Laura García-Perciante, and Roy Maartens for valuable comments. AS-V was supported by a CONACyT postdoctoral grant, and thanks the Relativity and Cosmology Group at Portsmouth for hosting him while this work was done.