text:
# Isospin dependence of the critical quark-deconfinement densities

## Abstract

We explore the dependence of the critical density, separating hadronic matter from a mixed phase of quarks and hadrons, on the ratio $`Z/A`$. We use both the MIT bag model and the Color Dielectric Model to describe the quark dynamics, while for the hadronic phase we employ various relativistic equations of state. We find that, if the parameters of quark models are fixed so that the existence of quark stars is allowed, then the critical density drops dramatically in the range $`Z/A`$ 0.3–0.4. Moreover, for $`Z/A`$ 0.3 the critical density is only slightly larger than the saturation density of symmetric nuclear-matter. This opens the possibility to verify the Witten-Bodmer hypothesis on absolute stability of quark matter using ground-based experiments in which neutron-rich nuclei are tested.

Hadronic matter is expected to undergo a phase transition into a deconfined phase of quarks and gluons at large densities and/or high temperatures. On very general grounds, the transition’s critical densities are expected to depend on the isospin of the system. Up to now, experimental data on the phase transition have been extracted from high-energy scattering of almost isospin-symmetric nuclei, having a proton fraction $`Z/A`$ 0.4–0.5. The analysis of observations of neutron stars, which are composed of $`\beta `$-stable matter for which $`Z/A`$ 0.1, can also provide hints on the structure of extremely asymmetric matter. No data on the quark deconfinement transition is at the moment available for intermediate values of $`Z/A`$. Recently it has been proposed by several groups to produce and study unstable neutron-rich nuclei (for a updated review of the status of these projects see e.g. Jonson (2002); Grunder (2002)). As we will show, these new experiments open the possibility to explore in laboratory the isospin dependence of the critical densities.

The information coming from experiments with ultra-relativistic heavy ions is that, for symmetric or nearly symmetric nuclear matter, the critical energy density must be considerably larger than nuclear matter energy density at saturation. Concerning non-symmetric matter, general arguments based on Pauli principle suggest that the critical density decreases with $`Z/A`$. We want to study in particular the range $`Z/A`$ 0.3–0.4, which can be explored in radioactive nuclear beam facilities. Here and in the following we are mainly interested in the critical density separating pure hadronic matter from a mixed phase of hadrons and quarks. The second critical density, separating the mixed phase from the pure quark matter phase, is in general reachable only in high energy experiments, while we are mainly interested in intermediate energy experiments.

A sistematic study of the isospin dependence of the critical densities has been performed up to now, to our knowledge, only by Mueller Mueller (1997). In that paper only one set of model parameters values is explored. The conclusion of Mueller (1997) is that, moving from symmetric nuclei to nuclei having $`Z/A0.3`$, the critical density is reduced by roughly 10$`\%`$. In this Letter we will explore in a more sistematic way the model parameters. In particular we will be interested in those parameters sets which would allow the existence of quarks stars Alcock et al. (1986); Haensel et al. (1986), i.e. parameters sets for which the so-called Witten-Bodmer hypothesis is satisfied Witten (1984); Bodmer (1971). According to that hypothesis, a state made of an approximately equal number of up, down and strange quarks can have an energy per baryon number $`E/A`$ smaller than that of iron ($`E_{Fe}930`$ MeV). To satisfy the Witten-Bodmer hypothesis, strong constraints on quark model parameters have to be imposed. For instance, using the MIT bag model, the so-called pressure-of-the-vacuum parameter $`B`$ must have a very small value, $`B^{1/4}`$ 140–150 MeV Farhi and Jaffe (1984). Assuming the Witten-Bodmer hypothesis to be true, ordinary nuclear matter would be metastable. In order not to contradict the obvious stability of normal nuclei, quark matter made of only two flavors must not be more stable than iron. As we shall see, slightly more strict boundaries on parameters’ value can be imposed by requiring not only iron, but also neutron rich nuclei like e.g. lead, to be stable. If the Witten-Bodmer hypothesis is satisfied, self-bound stars entirely composed of quark matter can exist Alcock et al. (1986); Haensel et al. (1986); Drago and Lavagno (2001). Recently several analysis of observational data have emphasized the possible existence of compact stars having very small radii, of the order of 9 kilometers or less Li et al. (1999); Dey et al. (1998); Pons et al. (2002); Drake et al. (2002). The most widely discussed possibility to explain the observed mass-radius relation is based on the existence of quark stars. It is therefore particularly interesting to envisage laboratory experiments testing the possible signatures of model parameters values that would allow the existence of these extremely compact stellar objects. What we are proposing in this Letter is to use radioactive nuclear beams to this purpose.

It is rather unlikely, at least in the near future, that neutron rich nuclei obtainable in radioactive nuclear beam facilities can be accelerated at very large energies, much larger than 1 GeV per nucleon. The scenario we would like to explore corresponds therefore to the situation realized in experiments at moderate energy, in which the temperature of the system is at maximum of the order of few tens MeV. In this situation, strange quarks cannot be produced and we need only to study the deconfinement transition from nucleonic matter into up and down quark matter.

In our analysis we have explored various hadronic and quark models. Concerning the hadronic phase, we have used the relativistic non-linear Walecka-type models of Glendenning-Moszkowski (GM1, GM2, GM3) Glendenning and Moszkowski (1991). We have also explored the possibility of enhancing the symmetry repulsion at high barion density introducing a coupling to a charged scalar $`\delta `$-meson. As remarked in Ref.Liu et al. (2002), this is fully in agreement with the spirit of the effective field theories, and of course with the phenomenology of the free nucleon-nucleon interaction. For the quark phase we have considered the MIT bag model at first order in the strong coupling constant $`\alpha _s`$ Farhi and Jaffe (1984) and the Color Dielectric Model (CDM) Pirner (1992); Birse (1990). In the latter, quarks develop a density dependent constituent mass through their interaction with a scalar field representing a multi-gluon state. After having chosen a model for the hadronic and for the quark EOS, the deconfinement phase transition is then described by imposing Gibbs equilibrium conditions Glendenning (1992). This technique can be considered as an effective way to describe the multi-quark correlations which generate the hadronization process in a more microscopic approach.

In Fig. 1 we show the critical density $`\rho _{cr}`$ separating nuclear matter from quark-nucleon mixed phase, as a function of the proton fraction $`Z/A`$. The figure has been obtained using GM2 (GM1) parametrization for the hadronic phase and the MIT bag model without gluons (with gluons and $`\alpha _s=0.3`$) in the upper and lower window, respectively. The most striking feature of the results shown in Fig.1 is the sharp decrease of $`\rho _{cr}`$ in the range $`Z/A`$ 0.3–0.4. The lower curves in each window correspond to parameters’ values satisfying the Witten-Bodmer hypothesis. In the latter case, and for $`Z/A0.3`$, the critical density is of the order of $`\rho _0`$. This opens the possibility to test the deconfinement transition in low energy experiments, such as the one performed in future RNB facilities.

The main features of Fig.1 can be easely understood if one recalls that we are investigating situations in which the minimum of pure quark matter EOS is at an energy just above or just below the minimum of the hadronic matter EOS. The first scenario is the one in which the absolute minimum of $`E/A`$, for a given value of $`Z/A`$, corresponds to the quark matter EOS (this situation corresponds to very small values of the parameter $`B`$, e.g. $`B^{1/4}=148`$ MeV in the top panel of Fig.1). In this case, the deconfinement transition starts at very small densities, smaller than nuclear matter saturation density. The numerical determination of these densities is rather delicate and we limit ourself to indicate with vertical arrows in Fig. 1 the behaviour of the critical density for a given value of $`B`$. If the value of $`B`$ is further reduced, the vertical arrow shifts towards larger values of $`Z/A`$ and therefore cannot correspond to a physically acceptable situation, since it would imply deconfinement into two flavor quark matter at low densities, even for almost symmetric nuclei. In particular we can exclude parameter’s values for which the dramatic drop in the critical density takes place for $`Z/A0.4`$. In this way we define a minimal value of, e.g., parameter $`B`$ of the MIT bag model. Notice that the limit on the value of the model parameters we get in this way is (slightly) more restrictive than the one generally adopted and based on the stability of Fe against decay into two flavor quark matter.

The second situation is the one in which the minimum of the quark EOS lies slightly above the hadronic minimum, as e.g. for $`B^{1/4}=150`$ MeV in the top panel of Fig.1. In this situation the deconfinement transition starts at a density slightly smaller than the one corresponding to the minimum of the quark EOS. The critical density cannot be further reduced, since at even smaller densities the energy $`E/A`$ in the quark phase rises dramatically, both in the MIT bag model and in the CDM, and therefore no mixing of hadronic matter with quark matter is possible at those densities. It is also important to notice that both in the MIT and in the CDM model, the position of the minimum of $`E/A`$ is near the value of nuclear matter saturation density when the energy of the minimum is near to the one obtained from the hadronic EOS Drago et al. (1995).

Finally, when the value of $`B`$ is further increased, the minima of the hadronic and of the quark EOS become more and more separated, the dependence of the critical density on the $`Z/A`$ fraction reduces progressively, and a situation similar to the one discussed in Ref.Mueller (1997) is reached.

In Fig.2 the same analysis is performed using the CDM for the quark EOS, obtaining results similar to those of Fig.1. In Fig.3 the effect of the exchange of the charged $`\delta `$ meson is also considered Liu et al. (2002). The $`\delta `$-exchange potential provides an extra isospin dependence of the EOS, and its effect shows up in a further reduction of the critical density in the region $`0.3Z/A0.4`$. We would like to remark that the NLH$`\rho `$ parametrisation gives results very similar to those obtained using the GM3 parameters set.

Let us now comment on the physical relevancy of the dramatic reduction of the critical density in neutron rich nuclei. Since in these nuclei neutrons presumably occupy an extended area around the core of the nucleus (neutron skin), the density of the very neutron-rich part can be considerably smaller than $`\rho _0`$. We cannot therefore expect to find a direct signal of deconfinement in the structure of these nuclei, but we can look for precursor signals. In particular, we can expect that the formation of clusters containing six or nine quarks will be enhanced due to the reduction of the critical deconfinement density. This enhancement can in turn be intepreted as a modification of single-nucleon properties due to the nuclear medium. The experimental search of effects like the one we are discussing here has a very long story, which includes the discovery of the EMC effect Aubert et al. (1983), that is the non trivial difference between free-nucleon and nuclear structure functions, for which many interpretations have been proposed (for a review see Arneodo (1994)). In particular, models for the EMC effect invoking the formation of multi-quark clusters have been Jaffe (1983); Carlson and Havens (1983) and are still quite popular Barshay and Kreyerhoff (2000). The decrease of the deconfinement critical density obtained in our analysis suggests a dependence of the EMC effect on the isospin, since the probability of forming virtual multiquark bags is enhanced in neutron rich nuclei. This dependence would add to the non-isoscalarity effect originated by the different structure functions of protons and neutrons, generally considered in the analyses Eskola et al. (1999). A full account of non-isoscalarity corrections arizing from genuine nuclear physics effects has been attempted in Barone et al. (2000b, a), and seems to be necessary in the light of the preliminary analysis of CHORUS data on deep inelastic scattering off a Pb target Oldeman (2000).

A more direct way to explore the reduction of the deconfinement critical density would be to test the EOS of matter via scattering of two neutron rich nuclei. This possibility is based on the analysis of intermediate-energy heavy-ion collisions, as discussed in Refs.Li et al. (1998); Li and Ko (1998); Li (2002a, b); Baran et al. (2002), and references therein.

To check the practical feasibility of such an experiment, we have performed some simulations of the 132Sn + 132Sn collision (average $`Z/A`$=0.38) at various energies, for semicentral impact parameter, $`b`$=6 fm, just to optimize the neutron skin effect in order to get a large asymmetry in the interaction zone. We have used a Relativistic Transport Code Fuchs and Wolter (1995), adoptining the same effective interaction Liu et al. (2002) of the $`NLH(\rho +\delta )`$ EOS of Fig.3 to compute the critical deconfinement density (solid line). In Fig.3 the paths in the $`(\rho ,Z/A)`$ plane followed in the c.m. region during the collision are reported at energies of 300 A.MeV (crosses) and 1 A.GeV (circles). We see that already at 300 A.MeV we are reaching the border of the mixed phase, and we are well inside it at 1 A.GeV. In order to be sure that we are really testing the Nuclear Matter EOS, which is an equilibrium property, we have performed a check of the local thermalization in correspondence of the high baryon density regions reached during the collision. In Fig.4 we show that indeed, when the maximum density is reached ($`\rho 2.6\rho _0`$) the quadrupole momentum of the nucleon momentum distribution has dropped to $`10\%`$ of its initial value, a signal that the system has indeed thermalized.

The use of an harder hadronic EOS for symmetric matter at high density would correspond to a even more favourable situation than the one presented in Fig.3, as shown in the examples of Fig.1. Moreover the use of neutron-richer nuclei would allow to test the EOS at smaller values of $`Z/A`$. To this purpose, the most promising nuclei are the ones near the r-process path, in particular for neutron numbers near the magic values N=82 or 126. In these regions, the proton fraction is as low as 0.32–0.33 and these nuclei could be studied in future experiments with neutron-rich beams Rehm (1999).

In conclusion, our analysis supports the possibility of observing precursor signals of the transition to a mixed quark-hadron phase in the collision, central or semi-central, of exotic (radioactive) heavy ions in the energy range of a few hundred MeV per nucleon. A possible signature could be revealed through an earlier softening of the hadronic EOS for large isospin asymmetries, and it would be observed e.g. in the behaviour of the collective flows for particles having large transverse momentum.

###### Acknowledgements.

It is a pleasure to thank V. Barone for many valuable suggestions.


format: text

text:
# Microscopic theory of network glasses

## Abstract

A molecular theory of the glass transition of network forming liquids is developed using a combination of self-consistent phonon and liquid state approaches. Both the dynamical transition and the entropy crisis characteristic of random first order transitions are mapped out as a function of the degree of bonding and the density. Using a scaling relation for a soft-core model to crudely translate the densities into temperatures, the theory predicts that the ratio of the dynamical transition temperature to the laboratory transition temperature rises as the degree of bonding increases, while the Kauzmann temperature falls relative to the laboratory transition. These results indicate why highly coordinated liquids should be ”strong” while van der Waals liquids without coordination are ”fragile”.

While much of the phenomenology of the glass transition is universal, shared by all dense liquids, everyday glasses are formed from liquids that possess a network structure at the molecular level. The universal features of the liquid-glass transition such as the Kauzmann entropy crisis and the Vogel-Fulcher rate law have been explained using theories based on an underlying random first-order transitionP.G.Wolynes (1992); T.R.Kirkpatrick and D.Thirumalai (1995); Xia and P.G.Wolynes (2000), like that of mean field models with one step replica symmetry breakingD.J.Gross et al. (1985). Furthermore, the variations of the mean relaxation timeXia and P.G.Wolynes (2000) and degree of nonexponentialityXia and P.G.Wolynes (2001) of relaxation from substance to substance have been quantitatively predicted using this approach, provided the configurational heat capacity is known as input, thus quantitatively explaining the patterns of behavior verbally characterized using the expression ”liquid fragility”C.A.Angell (1985). Detailed microscopic calculations buttressing these developments are based on a model of a dense interacting atomic fluid in which the particles interact via van der Waals forces alone. Is this starting point appropriate for the network glasses? To address this issue, we explore the glass transition in a model of a network forming liquid. The calculations presented are similar in spirit to those made earlier for atomic fluids. Our model assumes cross-linking in a liquid is fluctuating, not quenched as in models of rubberGoldbart and Goldenfeld (1989), taking the mean degree of crosslinking as given from experiment or from a separate calculation explicitly accounting for the equilibria between the various different bonded species. These calculations suggest a connection between the glassy behavior of network systems as viewed by random first order transition theory and rigidity percolation. Rigidity percolation features prominently in Phillips’s discussion of the role of bonding constraints in glass formationJ.C.Phillips (1985), as well as in the extensive studies of ThorpeS.Feng and M.F.Thorpe (1985); P.M.Duxbury et al. (1999). A major result of the present calculations is that the density of non-bonded interactions needed to trigger energy landscape dominated glassy behavior decreases when bonding increases. It follows that, consistent with experiment, highly cross-linked network glasses should show activated dynamics even at very high temperatures. The theory still predicts a crossover to non-activated dynamics for sufficiently high temperature as has recently been inferred in computer simulations of silicaHorbach and Kob (2001); Coluzzi and Verrocchio (2002).

The self-consistent phonon method applied to the amorphous state provides our starting pointFixman (1969); J.P.Stoessel and P.G.Wolynes (1984). This scheme is closely connected to approaches based on density functionalsY.Singh et al. (1985). It has recently been justified by an elegant clone or replica approachFranz et al. (1999); Mezard and Parisi (1999). We actually use a hybrid of both density functional and self-consistent phonon calculations. For the hard sphere system, both approaches agree there is a critical density $`\rho _A`$, below which localized density distributions corresponding to amorphous packings are mechanically unstable to even small thermal fluctuations. $`\rho _A\sigma ^3`$ from the self-consistent theory is approximately 1.0 for hard spheres, where $`\sigma `$ is the hard sphere diameter. This is in agreement with the computer simulated dynamical transitionWoodcock and Angell (1981). This instability, from the vitrified side, coincides with the transition found in dynamical mode coupling approaches which start from the liquid sideP.G.Wolynes (1992); T.R.Kirkpatrick and D.Thirumalai (1995); Götze (1999); Bengtzilius et al. (1984). The thermal vibration allowed at $`\rho _A`$ in a frozen amorphous state gives mean square fluctuations $`<\delta r^2>=\alpha ^1`$, where $`\alpha `$ is an effective spring constant and is finite at the transition. $`\alpha `$ provides an analog of the Edwards-Anderson order parameter for spins and exhibits discontinuous replica symmetry breaking. This limit of stability for an amorphous system is comparable to the Lindemann criterion that describes melting of crystalline solids. At densities above $`\rho _A`$ many different mechanically stable amorphous packings or free energy minima can exist. At densities beyond $`\rho _A`$, however, the configurational entropy of these minima decreases until it vanishes at a density $`\rho _K`$. Analogously, for soft-sphere models at fixed densities one can find a temperature $`T_A`$, signaling the onset of glassy dynamics (“landscape determined”S.Sastry et al. (1998))and a still lower temperature $`T_K`$, the Kauzmann temperature where the configurational entropy would vanish. Between these densities or temperatures, the dynamics can be described using the concept of entropic dropletsP.G.Wolynes (1992); Xia and P.G.Wolynes (2000); T.R.Kirkpatrick and D.Thirumalai (1995); Xia and P.G.Wolynes (2001). The activation free energy diverges as $`T_K`$ is approached, consistent with the Vogel-Fulcher law $`\mathrm{\Delta }F^{}=DT_K/(TT_K)`$. A microscopic calculation of the energy cost of forming an entropic droplet quantitatively explains the experimentally observed variation of D with the heat capacity drop upon vitrification, $`\mathrm{\Delta }C_p`$Xia and P.G.Wolynes (2000). This correlation relies on the idea that the Lindemann ratio near $`T_K`$ should depend only weakly on the detailed force law. We will check this supposition for our network glass model.

The free energy is computed for a given amorphous packing as a function of $`\alpha `$. This free energy should then be averaged over packings. Each atom can be assigned to a cell and the density can be written as $`\rho (𝐫)=_i\rho _i=_i\left(\frac{\alpha _i}{\pi }\right)^{3/2}exp\left[\alpha _i(𝐫𝐑_i)^2\right]`$ where the locations $`\{𝐑_i\}`$, characterize the packing as a set of fiducial locations around which the atoms vibrate. The free energy as a function of $`\{\alpha _i\}`$, $`F(\{\alpha _i\})`$, is analogous to the inter replica potential given as a function of q for spin glassesS.Franz and G.Parisi (1997).

More generally, the thermal vibrations are described by a tensor $`𝜶_i`$ defined by

$`𝜶_i^1`$ $`=`$ $`(𝐫_i𝐑_i)(𝐫_i𝐑_i)_{W_i^{eff}}`$ (1)
where the averaging is with respect to the effective potential (dependent on $`\{\alpha _j\}`$)

$`\mathrm{exp}\left(\beta W_i^{eff}\right)\mathrm{\Pi }_{ij}\mathrm{exp}\left(\beta V^{eff}(|𝐫_i𝐑_j|)\right)=\mathrm{\Pi }_{ij}{\displaystyle \mathrm{𝐝𝐫}_j\rho _j(𝐫_j)\mathrm{exp}\left(\frac{\beta }{2}V(r_{ij})\right)}`$ (2)
where $`V(r)`$ is the potential energy. For an isotropic glass under the assumption that all $`\alpha _i`$ are the same $`V^{eff}`$ becomes essentially a function of a scalar distance R, $`V^{eff}(R)`$, and one findsJ.P.Stoessel and P.G.Wolynes (1984) $`\alpha =\frac{\rho }{6}𝑑Rg(R)^2\left(\beta V^{eff}(R)\right)`$and a free energy $`\frac{F(\alpha )}{Nk_BT}=f^{ideal}+\rho 𝑑Rg(R)V^{eff}(R)`$. Although a good approximation for large $`\alpha `$, the self-consistent phonon theory is not accurate for low $`\alpha `$. At low $`\alpha `$, an expansion around the uniform fluid can be used, but this depends on knowing the direct correlation function which is not available for network fluids. When no bonds are present, alternatively one can use the functional first suggested by TarazonaP.Tarazona (1984) and generalized by AshcroftA.Khein and N.W.Ashcroft (1999). This functional takes the form $`\frac{F}{Nk_bT}f=\frac{1}{N}𝑑r\rho (r)\left[\mathrm{log}(\rho (r))1\right]+\mathrm{\Psi }(\eta )`$ where $`\mathrm{\Psi }`$ is the interaction part of the free energy for a monatomic liquid or glass (which we take to be the Carnahan-Starling expression) and $`\eta `$ is the packing fraction ($`\pi \rho \sigma ^3/6`$). For the network forming liquids we will correct this free energy using a modification of the result of ErukhimovichErukhimovich et al. (2001), which provides for the entropy loss of forming bonds.

For the network forming liquid and glass, we separate the Hamiltonian into bonding and non-bonding parts $`V=_{i,jbonded}V_b(r_{ij})+_{i,jnonbonded}V_{nb}(r_{ij})`$. Each pair i and j is considered to be either bonded on non-bonded and thus have an interaction described by either $`V_b`$ or $`V_{nb}`$. The molecular units which are considered to be bonded have a locally harmonic potential energy $`V_b(r_{ij})=V_0+\frac{1}{2}\gamma ^b(r_{ij}d_b)^2`$ with $`d_b`$ the equilibrium bond distance. $`V_{nb}`$, on the other hand, represents non-bonded interactions and therefore is taken to be a hard or soft sphere. Notice this potential does not explicitly describe, for any given configuration, whether a specific given pair should be considered bonded or not. This requires writing the potential using manybody switching functionsF.H.Stillinger and T.A.Weber (1985); L.R.Corrales and P.J.Rossky (1992). These manybody functions would reflect the saturation of chemical bonding tendencies (and would be needed for simulations.) For such a model the theory of associated liquids yields a (renormalized) equilibrium constant for the bonding reactions in terms of the underlying potentialL.W.Dahl and H.C.Andersen (1983). We will take the fraction of bonded pairs to be given and to be temperature and density independent. Of course, density and temperature dependent changes of the association equilibrium would be easy to include, but would distract from understanding how the glassy physics itself depends on temperature. Provided the association equilibrium constants are not paradoxically related to the density (i.e., decrease partially with $`\rho `$) the trends in any event will be as shown.

For the network system, the self-consistent phonon calculation must distinguish the two types of interactions. For specificity we assume any given atom can form 0-4 bonds, in a tetrahedral pattern. Other spatial patterns are easily treated as well. At large $`\alpha `$ the effective potential coming from a bonded pair is essentially $`V_b^{eff}(𝐑)=\frac{V_0}{2}+\frac{\gamma ^b}{4}\left(|𝐑|d_b\right)^2`$ where the factor of 1/2 comes from the self-consistent formulation. We set the bond distance $`d_b`$ equal to the hard sphere diameter for the non-bonded atoms and choose a zero of energy such that $`V_0=0`$. In numerical calculations we used $`\beta \gamma ^b=300`$, approximately the value appropriate for silica potentialsvan Beest et al. (1990). The probability that a given atom has n-bonds, $`p_n`$ is given by ErukhimovichErukhimovich et al. (2001) $`p_n=\left(\genfrac{}{}{0pt}{}{4}{n}\right)\mathrm{\Gamma }^n(1\mathrm{\Gamma })^{4n}`$ with $`\mathrm{\Gamma }=\frac{n_b}{4}`$ and $`n_b`$ is the average number of bonds per atom.

In our calculations the basic pair distribution for non-bonded interactions, $`g(R)`$ is taken to be the Percus-Yevick radial distribution function modified for high density in the standard way. To account for the presence of bonded interactions that block the approach of other partners, $`g(R)`$ for an atom with $`n`$ bonds (denoted $`g_n(R)`$) is scaled by $`\mathrm{max}(0,\frac{<nn>_\rho n}{<nn>_\rho })`$ for distances less than the first minimum of g(R). $`<nn>_\rho `$ is the number of nearest neighbors at total density $`\rho `$, obtained by integrating the first peak of g(R).

To each of the 5 types of atoms, those with 0, 1, 2, 3, or 4 bonds, we assign a mean $`\alpha _n`$. In this paper we use the uncoupled oscillator implementation of the self-consistent phonon theoryJ.P.Stoessel and P.G.Wolynes (1984) to determine the free energy. Essentially, this approximation treats the glass as a collection of Einstein oscillators. The derivation in Stoessel and Wolynes’s early paper is lengthy in detail but when applied to the present potential one obtains

$`F_{scp}/Nk_BT={\displaystyle \underset{n=0}{\overset{4}{}}}p_n\left[V_{b,n}^{eff}(d_b)+\rho {\displaystyle 𝑑Rg_n(R)V_{nb}^{eff}(R)}\right]`$
$`\mathrm{ln}{\displaystyle _v}𝑑r^Ne^{_{n=0}^4\frac{p_n}{2}r_i\left[\beta V_{b,n}^{eff}(R_i)+𝐈\alpha _n^{nb}\right]r_i}`$ (3)
where $`V_{b,n}^{eff}`$ is the sum of $`V_b^{eff}`$ for an atom with $`n`$ bonds over its bonds (and includes correlations between bonds through the use of a tetrahedral lattice), $`V_{nb}^{eff}`$ is defined by Eqn. 2 using only the non-bonding potential energy and $`g_n`$, and $`\alpha _n^{nb}`$ is the non-bonding contribution to $`\alpha _n`$. With a proper choice of local axes for each $`r_i`$, the argument of the exponential can be diagonalized, leading to 3 principal values of the vibrational tensor $`\alpha _{1,n}`$, $`\alpha _{2,n}`$, and $`\alpha _{3,n}`$. For the free energy to be consistent with Eqn. 1, we must have $`\alpha _n`$ satisfy $`3\alpha _n^1=\alpha _{1,n}^1+\alpha _{2,n}^1+\alpha _{3,n}^1`$ Note that for a glass transition to occur, all 3 components of $`\alpha _n`$ must be large. Given a bonding pattern, the $`\{\alpha \}`$ can be determined and, hence, the free energy.

$`F_{scp}`$ is the free energy of the large $`\alpha `$ glassy minimum. We compute the free energy of the network liquid with $`\alpha =0`$ using the Carnahan-Starling expression for a non-bonded liquid plus a modification of Erukhimovich’s result for equilibrium network materialsErukhimovich et al. (2001): $`F_{liq}/Nk_BT=\mathrm{\Psi }(\eta )+4\left[\mathrm{\Gamma }\mathrm{ln}\mathrm{\Gamma }+(1\mathrm{\Gamma })\mathrm{ln}(1\mathrm{\Gamma })\right]+\mathrm{ln}\left(\frac{\rho }{e}\right)\frac{n_b}{2}\mathrm{ln}\left(\frac{n_b\rho }{e}\right)+\frac{F_{scp}^b}{Nk_BT}`$ where $`F_{scp}^b`$ is obtained from Eqn. 3 with $`\alpha =0`$ and only the bonding terms. It corresponds to the entropy loss from bonding computed by ErukhimovichErukhimovich et al. (2001).

The self-consistent phonon theory alone allows us to find $`\rho _A`$, as the lowest density giving a non-zero $`\alpha `$, which thus corresponds to the onset of activated behavior while $`\rho _K`$, the ”Kauzmann” density at which the free energies of the liquid and glass match and the configurational entropy vanishes also relies on the bonding correction to $`F(\alpha =0)`$. The ratio $`\rho _A/\rho _K`$ dimensionlessly characterizes the thermodynamic aspects completely. To compare to the kinetic laboratory glass transition, we note the laboratory glass transition defined to occur when the viscosity reaches $`10^{14}`$ p. Previous work using random first order theory predicts this to occur when the configurational entropy is about 1.0 $`k_B`$. To make the translation of our thermodynamic results to the laboratory transition density in what follows we will therefore mean by $`\rho _G`$ the density where the liquid and glass free energies differ by 1.0 $`k_BT`$ per particle. The universality of configurational entropy at the laboratory transition is well confirmed, also by experiment.

$`\rho _A`$, $`\rho _K`$, and $`\rho _G`$ as functions of the average number of bonds $`n_b`$ are shown in Fig. 1. For $`n_b>3`$, $`\rho _A`$ vanishes, indicating that the rigidity percolation occurs when $`n_b=3`$ in our model. For comparison, Thorpe more precisely defining rigidity percolation as the absence of zero mode vibrational frequencies, finds $`n_b`$ = 2.4. For $`n_b>3`$, even at low densities the network is a rubber that can still vitrify at a higher densities. Fig. 2 shows the effective spring constants at these densities. For the hard sphere modelJ.P.Stoessel and P.G.Wolynes (1984) the uncoupled oscillator approximation leads to a somewhat lower $`\rho _A`$ than simulation and to relatively small, but nevertheless finite, values of $`\alpha _A`$. Quantitatively, the more complete coupled oscillator approximation (which uses a chain summation to treat the glass as a collection of Debye oscillatorsJ.P.Stoessel and P.G.Wolynes (1984)) leads to a larger $`\rho _A`$ in better agreement with computer studies. So we would expect quantitative modification by taking better account of the coupling but expect the same trend, namely that as the amount of bonding is increased, the transition densities decrease and therefore the transition temperatures increase. This theory then explains why the addition of a non-bonding impurity into a networked glass will lower the glass transition temperature, as known since ancient times for silicaPliny the Elder (xxxvi, 65). In physical terms, the non-bonding transition densities decrease because, as seen in Eqn. 3, as the contribution from bonding to the vibrational force increases, the remaining contribution needed from non-bonding interactions for self-consistency decreases.

It is clear that $`\rho _A`$ has a much sharper dependence on density than either $`\rho _K`$ or $`\rho _G`$. The degree of bonding has a very direct effect on the vibrational force constant, being proportional to the average Laplacian of the effective potential, but only an indirect effect on entropy, entering through the more slowly varying (in space) Boltzmann factor. The theory explains the change of fragility with bonding. To see this we note that real liquids have both repulsive cores of finite strength and attractions which conspire in setting the equation of state. Yet assuming the attractions provide merely a smooth background potential the density dependences of this model can be translated into temperature dependences. For a soft potential $`Vr^n`$ we have the scaling relation $`T_1/T_2=\left(\rho _2/\rho _1\right)^{3/n}`$. Using this relation for $`n=12`$ yields $`T_A/T_G`$ and $`T_K/T_G`$ as plotted in Fig. 3. Experimental measurements of $`T_K/T_G`$ give 0.82 for the fragile glass orthoterphenylAngell (1997), while for the strong glass silica the ratio has recently been determinedSipp et al. (2001) to be 0.36. The uncoupled oscillator approximation results in a variation from 0.82 for the non-networked glass to 0.6 at the approximate rigidity percolation limit, a variation somewhat smaller in magnitude. Doubtless, the crudeness of this result is due to our use of a scaling approximation appropriate to a simple, purely repulsive model and better agreement with experiment would be obtained by including attractive forces in our estimates. We see that $`T_A/T_G`$ increases dramatically with $`n_b`$, while $`T_K/T_G`$ decreases modestly. Thus highly bonded liquids have large temperature ranges of nearly Arrhenius activated behavior while liquids with few crosslinks have smaller ranges of strongly non-Arrhenius activated behavior.

For the uncoupled oscillator treatment $`\alpha _A`$ falls rapidly with density to very small values as crosslinking increases. The universality of the Lindemann ratio previously assumed by Wolynes and coworkersXia and P.G.Wolynes (2000) would seem in doubt. Yet this is not so. First, we know the coupled oscillator approximation will give a larger $`\rho _A`$ therefore decreasing the range of variation of $`\alpha _A`$. For $`n_b=0`$, the coupled oscillator gives a value of $`\alpha _A100`$, while the uncoupled theory gives $`\alpha _A20`$J.P.Stoessel and P.G.Wolynes (1984). More important, even at the uncoupled level $`\alpha _K`$ and $`\alpha _G`$ are roughly constant with increasing bonding. A careful examination of the earlier analysis shows that these are the appropriate Lindemann-like ratios to determine the activation barriers in the strongly glassy regime. The predicted universal relation between D and $`\mathrm{\Delta }C_p`$, relying on a universal $`\alpha _G`$, remains therefore sound.

The combination of self-consistent phonon and liquid state approaches provides a microscopic picture of network glasses that explains the dependence of the characteristic glass densities and, less directly, temperatures on degree of crosslinking and this results agrees with experiment. Doubtless, many of the detailed liquid state approximations made here can be improved but it is unlikely such modifications will change the qualitative story. In particular, the attractive forces that go beyond the soft-sphere model should probably be included to achieve quantitative agreement with experiment. This model combined with a theory of the association equilibria in forming the network provides a first principles treatment of glass properties starting from the underlying intermolecular forces.

RWH acknowledges NSF Grant. No. 9977124 for support.


text:
# The Black Hole–Bulge Relationship in QSOs

## 1 INTRODUCTION

The evolution of supermassive black holes and that of their host galaxies appear to be closely coupled. Every bulge system studied with high spatial resolution shows a central black hole, presumably a relic of AGN activity (Kormendy & Richstone 1995; Richstone et al. 1999; Kormendy & Gebhardt 2001). The black hole mass is roughly proportional to the luminosity of the bulge of the host galaxy, albeit with scatter of $`0.5`$ dex in $`M_{BH}`$ (Magorrian et al. 1998). Gebhardt et al. (2000a) and Ferrarese & Merritt (2000) found a tighter correlation involving $`M_{BH}`$ and $`\sigma _{}`$, where $`\sigma _{}`$ is the stellar line-of-sight velocity dispersion at radii outside the gravitational influence of the black hole. Tremaine et al. (2002) give this relationship as

$$M_{BH}=(10^{8.13}\mathrm{M}_{})(\sigma _{}/200\mathrm{km}\mathrm{s}^1)^{4.02}.$$
(1)
with an intrinsic scatter $`30\%`$. (The measure of scatter here is 1 $`\sigma `$ for $`M_{BH}`$ at fixed bulge luminosity or $`\sigma _{}`$, respectively.) Theoretical interpretations of this correlation (e.g., Silk & Rees 1998; Adams et al. 2001; Burkert & Silk 2001; Ostriker 2001; Balberg & Shapiro 2002; Haehnelt & Kauffmann 2000) variously would have the black hole form before, during, or after bulge formation. Given this uncertainty, measurements of $`M_{BH}`$ in galaxies with large look-back times would be valuable. Most promising would be a measurement of the $`M_{BH}\sigma _{}`$ relationship at high redshift, because the relatively small scatter of this relationship in the local universe should make it possible to discern even modest changes over time.

Unfortunately, measurements of $`M_{BH}`$ and $`\sigma _{}`$ at high redshift are challenging because spatially resolved kinematical studies are feasible only for nearby systems. However, recent work indicates that useful measures of $`M_{BH}`$ in AGN are possible (see Peterson 1997; Wandel, Peterson, & Malkan 1999). McLure & Dunlop (2001) and Wandel (2001) summarize the “reverberation” and “photoionization” methods to measure $`M_{BH}`$ in Seyfert galaxies. The results correlate with bulge luminosity in a way consistent with nearby galaxies, confirming the calibration of $`M_{BH}`$ for AGN. However, the host galaxy brightness is difficult to measure in the presence of a bright active nucleus, and a direct measurement of $`\sigma _{}`$ for quasar host galaxies using stellar absorption-lines is likewise difficult. Alternatively, the narrow emission-line widths can serve as a surrogate for $`\sigma _{}`$. Nelson & Whittle (1996) find agreement between narrow emission-line widths in AGNs and widths measured from stellar absorption-line kinematics. On this basis, Nelson (2000) proposes to use the width of the \[OIII\] $`\lambda \lambda 5007,4959`$ lines of AGN as a surrogate for $`\sigma _{}`$. Thus, one takes $`\sigma _{}\sigma _{[OIII]}`$, where $`\sigma _{[OIII]}\mathrm{FWHM}([\mathrm{O}\mathrm{III}])/2.35`$ and FWHM is the full width at half maximum. The divisor 2.35 relates $`\sigma `$ and FWHM for a Gaussian profile. Taking black hole masses from reverberation mapping, Nelson finds no systematic offset in the $`M_{BH}\sigma _{[OIII]}`$ relationship between low redshift AGN measured this way and measurements of galaxies based on observed stellar velocity dispersions.

In this paper, we use emission-line widths to study the $`M_{BH}\sigma _{[OIII]}`$ relationship for a sample of quasars at redshifts up to 3.3. Our goal is to assess the evolutionary history of black hole growth compared with that of the host galaxy. In Section 2, we describe the derivation of $`M_{BH}`$ and the use of the \[O III\] line width as a surrogate for $`\sigma _{}`$. Then we describe our adopted data set, which relies on published data covering a range of redshift together with our unpublished observations of low redshift QSOs. In Section 3, we present the results, and we examine the trend of $`M_{BH}`$ with $`\sigma _{[OIII]}`$ and the redshift dependence of this relationship. In Section 4, we discuss these results and the needed improvements. Throughout our discussion, we use a cosmology with $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1,\mathrm{\Omega }_\mathrm{M}=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. All values of luminosity used in this paper are corrected to our adopted cosmological parameters.

## 2 METHOD

### 2.1 Calculation of Black Hole Masses

The derivation of black hole masses from AGN broad line widths has been discussed by a number of authors, recently including Laor (1998), Wandel et al. (1999), Kaspi et al. (2000), Vestergaard (2002), and McLure & Dunlop (2001). The method relies on assuming that the line widths, at least for some lines, are dominated by orbital motion of the emitting gas in the gravitation potential of the black hole. This is supported by the generally symmetrical time variability of the wings, and the decrease in line width with increasing radius for different lines in a given object (see Wandel et al. 1999, and references therein). The black hole mass is then given by $`M_{BH}=v^2R/G`$, where $`v`$ and $`R`$ are an appropriate velocity and radius for the BLR. Deriving $`v`$ and $`R`$ from the observations is not entirely straightforward, as the emitting gas spans a range of radii and velocities. Pragmatically, the radius is derived from “reverberation mapping” or “echo mapping” studies that monitor the variation of the continuum and emission lines. The characteristic time lag between continuum variations and the response of a given line gives a measure of the radius of the region emitting that line (see review by Peterson 1993). Given this determination of $`R`$, one needs a measure of line width which, combined with this radius, gives the correct black hole mass. This choice can be parameterized as $`v=f\times \mathrm{FWHM}`$ for the line. Some authors use $`f=\sqrt{3}/2`$, appropriate for isotropic velocities. However, McLure & Dunlop (2001) argue that allowance for a flattened geometry of the BLR is preferable. Most work has employed the Balmer lines, in particular H$`\beta `$ . Black hole masses determined in this way for AGN with measured $`\sigma _{}`$ show overall agreement with the $`M_{BH}\sigma _{}`$ relationship (Gebhardt et al. 2000b; Ferrarese et al. 2001; McLure & Dunlop 2001).

The work of Wandel et al. (1999) and Kaspi et al. (2000) generated a set of BLR radii based on variability of the Balmer lines for 17 Seyfert 1 (Sy 1) galaxies and 17 PG QSOs. This work supports earlier indications that the BLR radius increases with luminosity, approximately as $`RL^{0.5}`$. Photoionization physics suggests $`RL^{0.5}`$, and this is also consistent with the idea that the BLR radius may be limited by the minimum radius at which dust grains survive (Netzer & Laor 1993). This relationship is enormously useful, for it opens the door to measuring black hole masses in large numbers of AGN from measurements of line width and luminosity at a single epoch. By determining the BLR radius in this way, one avoids the need for long series of observations of variability, which are in any case impractical for luminous, high redshift QSOs with long variability timescales. The use of such “photoionization masses”, calibrated in terms of echo results, has been discussed by several authors, including Wandel et al. (1999) and Vestergaard (2002). There is some controversy over the slope of the radius-luminosity relationship. Kaspi et al. (2000) find $`RL^{0.7}`$, based on echo radii. However, McLure & Jarvis (2002; see also Maoz 2002), find $`RL^{0.61}`$. On this basis, they fit the echo masses with $`M_{BH}=(10^{7.63}\mathrm{M}_{})v_{3000}^2L_{44}^{0.61},`$ where $`v_{3000}FWHM(\mathrm{H}\beta )/3000\mathrm{km}\mathrm{s}^1`$ and $`L_{44}\lambda L_\lambda (5100\mathrm{\AA })/(10^{44}\mathrm{erg}\mathrm{s}^1)`$, using the same cosmology that we have adopted. The difference in slope results from several factors, including a new echo mass for NGC 4051, different continuum luminosities, and the cosmological model. Here we take for our primary calibration the physically motivated slope $`RL^{0.5}`$. We adopt the $`L^{0.5}`$ fit shown in Figure 6 of Kaspi et al.,

$$M_{BH}=(10^{7.69}\mathrm{M}_{})v_{3000}^2L_{44}^{0.5},$$
(2)
where we have adjusted the coefficient to our cosmology. (This expression agrees with our own fitting of the Kaspi et al data and depends little on NGC 4051 because the slope is fixed.)

We use the H$`\beta `$ line width in this fashion to determine $`M_{BH}`$ , adopting the calibration of equation (2) for most of our discussion but illustrating the key results also for $`M_{BH}L^{0.61}`$. McLure & Jarvis and Vestergaard (2002) discuss the use of other lines. The use of the $`\lambda 5100`$ continuum follows Wandel et al. (1999) and Kaspi et al. (2000). McLure & Jarvis (2002) and Laor et al. (1997) argue that the $`\lambda 3000`$ Å continuum may work better, but the differences are not great. The $`\lambda 5100`$ continuum has the practical advantage that it can be measured in the same spectra as the H$`\beta `$ and \[O III\] lines, and it is less affected by dust extinction.

### 2.2 \[O III\] Lines in AGN

Nelson & Whittle (1996) have compared \[O III\] line widths and stellar velocity dispersions in AGN, finding generally good agreement. For the quantity log $`\sigma _{[OIII]}/\sigma _{}`$, they find a mean $`0.00\pm 0.01`$ and a dispersion $`\sigma =0.20`$, supporting the idea that the motions of the NLR gas are largely determined by the gravitational potential of the host galaxy. This is reinforced by the analysis by Nelson (2000), who shows essentially that $`\sigma _{[OIII]}`$ and $`M_{BH}`$ obey equation (1) for AGN with echo values of $`M_{BH}`$. These results support the use of $`\sigma _{[OIII]}`$ as a surrogate for $`\sigma _{}`$. A caution, however, is that \[O III\] profiles often have substantial asymmetry and non-Gaussian profiles, possibly resulting from outflow combined with extinction of the far side of the NLR (e.g., Wilson & Heckman 1985; Nelson & Whittle 1995).

Radio loud AGN tend to have stronger \[O III\] emission than radio quiet objects, as reflected in “Eigenvector 1” of Boroson & Green (1992). Radio jets may contribute to the motions of the NLR gas (Nelson & Whittle 1996). For this reason, we emphasize radio quiet objects in this paper, and we discuss radio loud objects separately.

### 2.3 New Observations

We include here results from an unpublished set of spectra of QSOs obtained at McDonald Observatory with the Large Cassegrain Spectrograph (LCS) on the 2.7-meter telescope. These objects correspond to the X-ray sample studied by Laor et al. (1997), which is in turn taken from the Bright Quasar Survey based on the Palomar Green (PG) survey (Boroson & Green 1992; Schmidt & Green 1983). We take values of FWHM($`\mathrm{H}\beta `$) from Boroson & Green. The instrumental resolution was 150 to 180 $`\mathrm{km}\mathrm{s}^1`$ FWHM, depending on redshift. Emission from Fe II blends was removed from the spectra with the aid of the Boroson & Green (1992) template. Table 1 gives the \[O III\] widths for 2 RL and 14 RQ objects in the redshift range 0.09 to 0.33 derived from our spectra. We use here the results of directly measuring the half-maximum point of the observed line profile rather than by fitting any kind of curve. The width of $`\lambda 4959`$ was noted for corroboration. (See below for a discussion of methods of measuring the \[O III\] line width.) Our own measurements of FWHM for H$`\beta `$ ranged from 0.95 to 1.35 times the BG92 values, with a mean ratio of 1.09. This may give some indication of systematic uncertainties arising from Fe II subtraction and fitting procedure, although real temporal variations may affect individual objects. For our PG sample, continuum luminosities were taken from Laor et al. (1997). Table 2 lists the redshift, adopted line widths, continuum luminosity, and black hole mass derived from equation (2).

### 2.4 Data from the literature

We have drawn observations of H$`\beta `$ and \[O III\] from several published sources. These are listed below, roughly in order of decreasing redshift. The results are given in Table 2. The quoted line widths are intrinsic values after subtraction of the instrumental width.

(1) Dietrich et al. (2002, “D02”) present infrared spectra for six QSOs with $`z3.4`$, placing H$`\beta `$ and \[O III\] in the infrared $`K`$-band wavelength region. Of these, Q0256-0000 and Q0302-0019 have \[O III\] lines of adequate strength to measure the line width. We have measured the FWHM of \[O III\] and H$`\beta `$ from the original data, using a direct measurement of the FWHM. We subtract in quadrature the instrumental width of 400 $`\mathrm{km}\mathrm{s}^1`$. We find $`\mathrm{FWHM}([\mathrm{O}\mathrm{III}])=838\pm 29\mathrm{km}\mathrm{s}^1`$ for Q0256-0000 and $`743\pm 19\mathrm{km}\mathrm{s}^1`$ for Q0302-0019, where the quoted errors reflect only the noise in the data. (For a single Gaussian fit to the \[O III\] profiles, we find 886 and 846 $`\mathrm{km}\mathrm{s}^1`$, respectively. These are 0.02, 0.05 dex larger than the direct measurements. The sense of the difference is typical, but the magnitude is not significant for our purposes. The appropriateness of Gaussian profiles is discussed below.) We corrected for the narrow component of $`\mathrm{H}\beta `$ by assuming a typical ratio of $`\lambda 5007`$ to narrow H$`\beta `$ of 10 to 1, as discussed in Brotherton (1996b, “B96b”) and references therein. The quoted noise errors and the likely error resulting from subtraction of the narrow H$`\beta `$ component are small compared to scatter among objects described below. Continuum luminosities are discussed below.

(2) McIntosh et al. (1999, “M99”) give $`H`$-band spectra of 32 luminous QSOs at $`2.0z2.5`$, covering H$`\beta `$ and \[O III\]. The instrumental resolution is about 500 $`\mathrm{km}\mathrm{s}^1`$. The signal-to-noise ratio (S/N) of the spectra is poor in some cases. Therefore, we present in Figures 1 and 2 two different data selections. For the “full” sample, we include all objects for which M99 give FWHM for both H$`\beta `$ and \[O III\]. For the “select” sample, we include those QSOs for which the greatest of the positive and negative errors in FWHM for both H$`\beta `$ and \[O III\], as a fraction of FWHM, is 0.33 or less. This results in a list of 4 radio loud (RL) and 4 radio quiet (RQ) objects. The “full” sample offers a larger number of objects for averaging purposes, and it avoids the question of bias in selecting the better measurements. However, it contains some rather uncertain measurements and may exaggerate the true dispersion of the data. The “select” sample is consistent with a reasonable judgment of the most reliable measurements, based on the S/N of the spectra. (The spectra are given in M99.) M99 quote a FWHM for the “total H$`\beta `$” profile and for the “broad H$`\beta `$” component. The “broad $`\mathrm{H}\beta `$” width exceeds the “total” width by more than can be accounted for by the removal of a narrow H$`\beta `$ component with a typical intensity ratio to \[O III\] (see above). We have carried out our own fits to the M99 data for some representative objects, including a broad H$`\beta `$ component plus a narrow H$`\beta `$ component with the same width as the \[O III\] lines. We find that the width of the broad H$`\beta `$ component is fairly close to M99’s “total H$`\beta `$ ” width and much narrower than their quoted “broad $`\mathrm{H}\beta `$” width. Accordingly, we use FWHM for $`\lambda 5007`$ and H$`\beta _{total}`$ from Table 4 of M99 and a RL/RQ type from their Table 2. See Vestergaard (2002) for a discussion of the appropriate measure of H$`\beta `$ width for determinations of black hole mass. Average errors in $`FWHM([OIII])`$ are $`\pm 0.12,0.16`$ dex for the RL, RQ objects in the “full” sample and 0.07, 0.10 dex for the “select” sample, respectively. Errors are in most cases smaller for H$`\beta `$ than \[O III\] width. Error bars are omitted from our figures for clarity, and the reader should bear these errors in mind.

(3) Brotherton (1996a, “B96a”) gives results of infrared spectroscopy of 18 RL and 14 RQ QSOs ranging in redshift from 0.7 to 2.5. Of these, 11 RQ and 7 RL objects were observed with the Cryogenic Spectrometer (CRSP) on the 2.1 meter telescope at Kitt Peak National Observatory (KPNO). The remaining objects were observed with other setups giving inadequate spectral resolution, or are included in B96b. All our data sources give line widths corrected for instrumental resolution except B96a. For B96a, we subtract in quadrature the instrumental resolution of 470, 600, and 460 $`\mathrm{km}\mathrm{s}^1`$ for I, J, and H band observations, respectively, as estimated by Brotherton (2002). Our “full” data set includes all 11 CRSP radio quiet objects. For the “select” sample, we use all objects for which the fractional error in $`FWHM([OIII])`$ is $`0.10`$ (Brotherton 2002). This is based on the error in the corrected line widths resulting from the formal error in the original Gaussian fits by B96a. These errors are not directly comparable with those of M99, which involve template fitting and a Monte Carlo error analysis. However, the chosen cutoff agrees with a reasonable judgment of the most reliable measurements, based on the S/N of the spectra. Among the radio quiet objects, the selected ones comprise four of the six CRSP objects for which B96a gives a quality flag of “A”. The select RL sample consists of 4 of the 5 CRSP “A” objects. (The fifth is PKS 0424-131, which is also in the M99 select sample. We include the B96a measurement in Table 2 for comparison.) The \[O III\] widths tabulated by B96a involve single Gaussian fits. From the spectra, we have estimated $`FWHM([OIII])`$ by direct measurement and find typical agreement within about 10% with the adopted Gaussian fits for the RL and RQ select objects. B96a also made no correction for the NLR contribution to the H$`\beta `$ profile. None of the spectra shows a prominent narrow line component of H$`\beta `$ . For a typical ratio of $`\lambda 5007`$ to narrow H$`\beta `$ of 10 to 1, the correction would on average increase FWHM(H$`\beta `$) by 11% for the RL select objects and 7% for the RQ select objects. These differences are not significant for our purposes. In addition to FWHM for H$`\beta `$ and \[O III\] we took from Chapter 6 of B96a the apparent $`V`$ magnitudes and RL/RQ classifications of the subject QSOs. We note that 1120+01 is a lens candidate, although it may be a true binary QSO (Kochanek 2002). Michalitsianos et al. (1997) suggest that the amplification factor could be as large as 100. Correction for lensing would lower $`M_{BH}`$ as $`L^{0.5}`$ from the value in Table 2 and Figure 1. Perhaps the $`M_{BH}\sigma _{[OIII]}`$ relation has potential as an indicator of lensing.

(4) B96b studied profiles of H$`\beta `$ and \[O III\] in 60 radio loud quasars with $`0.05z0.93`$. The instrumental resolution was 150 to 300 $`\mathrm{km}\mathrm{s}^1`$. We used the FWHM for H$`\beta `$ and \[O III\] from his Tables 3 and 4. Brotherton gives a quality flag A or B for 37 objects. We judged 6 of these to be unsuitable by visual inspection of the spectra, typically because the broad H$`\beta `$ line was weak or had an unusual profile. (The spectra are given in in B96b.) We also eliminated those objects for which we could not obtain a flux density at 5100 Å rest wavelength from available spectra, as described below. This left a sample of 23 QSOs that are included in Table 2.

(5) Grupe et al. (1999, “G99”) discuss optical emission-line properties of 76 bright soft X-ray selected AGN, based on spectra with a resolution of $`5`$ Å FWHM ($`300`$ $`\mathrm{km}\mathrm{s}^1`$). Redshifts and line widths are taken from Table 1 of G99. The \[O III\] line width was measured using a single Gaussian fit. We selected the 33 objects for which the fractional uncertainty in FWHM was better than 0.10 both for H$`\beta `$ and \[O III\]. Four of these were eliminated as unsuitable on the basis of inspection of the spectra as reproduced in Grupe (1996), and one for radio confusion. Two of the remaining objects were found from NED <sup>1</sup><sup>1</sup>1The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.  to be radio loud, leaving a sample of 26 radio quiet objects.

### 2.5 Continuum Luminosity

The luminosities for our samples come from heterogeneous sources. An assumed power law $`F_\nu \nu ^{0.5}`$ was used when necessary to scale the measurements to 5100 Å rest wavelength. Magnitudes were converted to flux densities $`F_\nu `$ as prescribed by Allen (1973). When absolute magnitudes or specific luminosities were given, these were adjusted to our adopted cosmology.

For D02, the continuum flux density at 5100 Å rest wavelength was measured from the published spectra. For M99 we used $`L_\nu (\mathrm{V})`$ from their Table 2, and for B96a we used the apparent V magnitude given in his Chapter 6. For B96b we used flux densities taken from spectrophotometry described by Netzer et al. (1995) and Wills et al. (1995) and unpublished spectrophotometry (Wills and Brotherton 2002) when available and otherwise took them from Figure 1 of B96b. Objects were omitted when neither of these was available, or when the result differed by more than 0.5 dex from the flux density at 5100 Å estimated from the absolute visual magnitude $`M_{abs}`$ quoted by B96b (taken from Véron-Cetty & Véron 1989). For G99, we measured the optical continuum flux density from the published spectra.

We have not attempted a correction for a contribution to the continuum from the host galaxy. This is unlikely to have a large effect on our conclusions, because $`M_{BH}`$ scales only as $`L_\nu ^{1/2}`$. We have estimated the host galaxy luminosity from the value of $`\sigma _{[OIII]}`$ for our individual AGN, using the Faber-Jackson relation (Forbes & Ponman 1999; Bernardi et al. 2001), $`M_V=20.572.5[4\mathrm{log}(\sigma _{[\mathrm{O}\mathrm{III}]}/220\mathrm{km}\mathrm{s}^1)]`$. For the adopted RL and RQ objects in the D02, M99, B96, B96b, and PG samples, the average value of $`L_{\mathrm{bulge}}/L_{\mathrm{observed}}`$ is only 2 or 3 percent, and 7 percent for G99. Thus, we have assumed that the measured flux is dominated by the AGN continuum for the adopted spectra of all objects.

### 2.6 Radio Loudness

Objects in the sample were divided into radio-quiet and radio-loud using the radio-to-optical ratio, $`R_{ro}=F_\nu (5\mathrm{GHz})/F_\nu (4400\mathrm{\AA })`$ as defined by Kellermann et al. (1989), where $`R_{ro}>10`$ designates a radio-loud QSO. For our PG sample, we used values for $`R_{ro}`$ from Boroson & Green (1992); for B96a we used values from that reference; for D02 we used values from Hooper et al. (1995); and for M99 we used values from that paper.

For G99 we obtained flux densities at 5 GHz, 4.85 GHz, or 1.4 GHz from photometric points in NED; otherwise, we measured the flux density from NRAO VLA Sky Survey (NVSS) radio flux density maps at 1.4 GHz (Condon et al. 1998). We assumed unresolved single, or very occasionally, double unresolved components. The angular resolution of NVSS of 45<sup>′′</sup> (Condon et al. 1998) includes the extended structure for nearly all QSOs. When no detection was listed or evident on the maps, we assumed that the object was radio quiet. For radio flux densities at frequencies other than 5 GHz we extrapolate to $`F_\nu (5\mathrm{GHz})`$ assuming $`F_\nu \nu ^{0.7}`$. Subtleties in the definition of $`R_{ro}`$ and continuum slope have little effect for our purposes, because of wide range in $`R_{ro}`$ values between RQ and RL objects. For purposes of computing $`R_{ro}`$, optical continuum flux densities were found as described above, or taken from NED, and corrected to rest wavelength $`\lambda 4400`$ using $`L_\nu \nu ^{0.5}`$.

## 3 RESULTS

### 3.1 The $`𝐌_{\mathrm{𝐁𝐇}}\sigma _{[𝐎\mathrm{𝐈𝐈𝐈}]}`$ relationship

Figures 1 and 2 show results for the RQ objects, based on equation (2). Figure 1 includes results for the full samples described above for the M99 and B96a data sets, and Figure 2 shows the select samples. (The data for the other sources is the same.) On average, the QSO results show good agreement with equation (1), shown as the straight line in Figures 1 and 2. This agreement supports the use of photoionization masses for $`M_{BH}`$ and the use of the \[O III\] line width as a surrogate for $`\sigma _{}`$ to masses over $`10^9\mathrm{M}_{}`$. However, we emphasize that the close alignment of the massive, high redshift objects in the upper right part of Figure 2 is better than expected from the errors (see discussion above) and must be to some degree fortuitous. The scatter in the PG sample may give a better indication of the scatter resulting from measurement errors together with true variations among objects, as discussed below. We assign significance only to the average agreement of the high redshift points as a group with the $`M_{BH}\sigma _{}`$ trend.

We also show in Figures 1 and 2 the objects listed in Table 1 of Nelson (2000). The Nelson points agree in the mean with our results but show less scatter around the Tremaine et al. (2002) line. The echo masses used by Nelson are presumably more accurate than the photoionization masses used here. Also, Nelson restricted his study to data with spectral resolution R $`>`$ 1500. (Included in the Nelson points in Figures 1 and 2 are two objects in common with our PG sample, namely 0953+414 and 1411+442.)

If the high redshift and Nelson (2000) points are removed from Figure 1, the scatter in the remaining points is too large to define accurately the slope of the $`M_{BH}\sigma _{[OIII]}`$ trend. However, our philosophy here is that the $`M_{BH}\sigma _{}`$ relationship is well established, and we are interested in examining the agreement of the AGN data with that relationship. The lower redshift points agree in the mean with equation 1, and this supports our working hypothesis that photoionization masses and $`\sigma _{[OIII]}`$ can be used as a substitute for direct measurements of $`M_{BH}`$ and $`\sigma _{}`$ in AGN. This is the basis for our discussion of the redshift dependence below.

Nelson & Whittle (1996) found that AGN with powerful, linear radio sources sometimes have $`FWHM([OIII])`$ larger than expected for the value of $`\sigma _{}`$. Figure 3 shows the radio loud objects from our data set, using the select sample for B96a and M99. Objects from most of the data sets scatter around equation (1). However, the measurements of B96b generally stand above and to the left of the trend. One issue may be the procedure for measuring the FWHM of \[O III\]. Some workers fit a Gaussian and quote the FWHM of this fit. B96b instead gives a direct measurement of the width at the half-maximum level of the observed line profile, which is feasible only when the data have good spectral resolution and signal-to-noise ratio. For four objects in common between the single Gaussian fit of B96a and the direct measurement of B96b, the B96b width on average is narrower by $`0.2`$ dex. A single Gaussian fit may be affected by broad wings or asymmetry of the line profile and by a tendency to smooth the line peak. This may explain some of the displacement of the B96b points in Figure 2. On the other hand, we also made a direct measurement of $`FWHM([OIII])`$ for our PG sample, and in Figures 1 and 2 our results do not share the displacement of the B96b sample. Consistent with this, we find that Gaussian fits to \[O III\] for our PG sample are only 0.04 dex wider on average than the direct measurements used here. As a further illustration, we measured $`FWHM([OIII])`$ using a Gaussian fit for many of the M99 objects, using the original data. The Gaussian widths are always wider than the M99 values, which should approximate a direct measurement because of the multiple template fit used by M99. For 3 RQ objects in our select sample for M99, the Gaussians are wider by an average of 0.08 dex. However, a few objects, typically with weak or strongly asymmetrical \[O III\], show differences in FWHM approaching a factor 2. Accurate results using \[O III\] as a surrogate for $`\sigma _{}`$ will require careful attention to the shape of the line profile.

### 3.2 Redshift Dependence

The $`M_{BH}\sigma _{[OIII]}`$ relationship shown in Figures 1 and 2 provides a basis for assessing the evolution of the relationship with redshift. The B96a, D02, and M99 objects fall in the redshift range 1 to 3.3. These objects do not systematically depart from the overall $`M_{BH}\sigma _{[OIII]}`$ relationship. In order to display the trend as a function of redshift, we use a fictitious “\[O III\] mass” $`M_{[OIII]}`$ defined as the mass given by equation (1) with $`\sigma _{[OIII]}`$ in place of $`\sigma _{}`$ . The quantity $`\mathrm{\Delta }\mathrm{log}\mathrm{M}\mathrm{log}\mathrm{M}_{\mathrm{BH}}\mathrm{log}\mathrm{M}_{[\mathrm{O}\mathrm{III}]}`$ is then a measure of the vertical deviation of a given object from the Tremaine et al. (2002) mean trend in Figures 1 - 3. In Figures 4 and 5, we plot this quantity as a function of redshift for our full and select RQ samples, respectively. In the mean, the high redshift objects show little deviation from the $`M_{BH}\sigma _{}`$ trend. Because the most massive objects are the high redshift objects, we cannot independently assess the slope of the $`M_{BH}\sigma _{[OIII]}`$ relationship and the redshift dependence. Our point is that the data are consistent with a $`M_{BH}\sigma _{[OIII]}`$ relationship that extends to high masses with a slope close to that of Tremaine et al. (2002) and that does not evolve strongly with time.

A quantitative measure of the adherence of the high redshift objects to the low redshift trend is difficult, given the small number of high redshift objects in the select sample and the various systematic errors. We may estimate the intrinsic scatter of the high redshift samples on the basis of the PG sample, which has the higher luminosity of the low redshift samples. The RQ objects in our PG sample show a mean $`\mathrm{\Delta }\mathrm{log}\mathrm{M}`$ of +0.10 and a dispersion 0.5 ($`1\sigma `$). The mean of $`\mathrm{\Delta }\mathrm{log}\mathrm{M}`$ is -0.11 for the 9 high redshift RQ objects in the select sample. The standard deviation of the mean for these 9 objects, if drawn from a population with similar dispersion to the PG sample, would be about 0.2 dex. For the full RQ sample, there are 23 objects with a mean $`\mathrm{\Delta }\mathrm{log}\mathrm{M}`$ of -0.15, which suggests that there is no significant bias in our select sample. These results indicate that, at times corresponding to redshifts $`z2`$ to 3, the average $`M_{BH}`$ at a given $`\sigma _{[OIII]}`$ was within a factor 2 or 3 of the present day value.

One systematic uncertainty is the slope of the radius–luminosity relationship for the BLR. The high redshift RQ select objects are on average 1.8 dex more luminous than the PG sample. Use of the McLure and Jarvis (2002) relationship, with $`RL^{0.61}`$, elevates the high redshift points by $`0.2`$ dex in Figures 1 - 5 while having little effect on the low redshift samples. This brings the high redshift and PG samples into even closer agreement, with a mean $`\mathrm{\Delta }\mathrm{log}\mathrm{M}`$ of +0.15 for the RQ PG sample and +0.14 for the RQ high redshift select sample. Figure 6 shows $`\mathrm{\Delta }\mathrm{log}\mathrm{M}`$ versus redshift for the McLure and Jarvis calibration (given before equation 2 above). Another uncertainty is the slope of the $`M_{BH}\sigma _{}`$ relationship, discussed by Tremaine et al. (2002). The high redshift objects in Figure 1 have $`\sigma _{}500\mathrm{km}\mathrm{s}^1`$. If the slope were $`M_{BH}\sigma _{}^{3.5}`$ rather than $`\sigma _{}^{4.0}`$, for example, the predicted mass for the high redshift objects would be lowered by about 0.2 dex. This would raise these objects by this amount in Figures 4 - 6.

## 4 DISCUSSION

There are two main points that come from our analysis. First, Figures 1 and 2 show that we can place AGN on the $`M_{BH}\sigma _{}`$ correlation using the \[OIII\] emission-line width as a surrogate for $`\sigma _{}`$. This bypasses the need to obtain high signal-to-noise spectra in order to obtain the host galaxy absorption line kinematics. Nelson (2000) was the first to point this out. We confirm it with our independent data set including luminous QSOs and using photoionization masses rather than reverberation masses.

Second, given that one can measure \[OIII\] line profiles at fairly modest signal-to-noise levels in AGN, one can now study the relationship between black holes and their host galaxies at high redshift. Although a larger sample of high redshift objects is needed, Figures 4 - 6 suggests that the $`M_{BH}\sigma _{[OIII]}`$ relationship is not a strong function of redshift. For the adopted cosmology, the age of the universe was 2.0 Gyr at $`z=3.3`$ and 3.3 Gyr at $`z=2.0`$. The objects in this redshift range have black hole masses up to $`10^{10}\mathrm{M}_{}`$ and implied host galaxy masses up to $`10^{13}\mathrm{M}_{}`$, if $`M_{\mathrm{bulge}}10^{2.8}M_{BH}`$ (Kormendy & Gebhardt 2001). Our results suggest that black holes typically grow contemporaneously with their host galaxy bulges, or else that both were well formed by $`z3`$.

Recent results on black hole demographics suggest that high mass black holes acquire most of their mass from luminous accretion during episodes of QSO activity (Yu & Tremaine 2002, and references therein). Yu & Tremaine find that half the black hole mass is accreted before a redshift $`z1.8`$, and only 10% before $`z=3`$. Thus, our results suggest that the black hole-galaxy bulge relationship is roughly obeyed, at least by very massive holes, at a time when much of the growth of present-day black holes lay in the future.

Previous workers have noted a tendency for \[O III\] widths to increase with luminosity for AGN (e.g., B96b; M99; Véron-Cetty, Véron, & Gonçalves 2001). In the present context, this is a natural consequence of the tendency to have larger $`M_{BH}`$ in more luminous objects, together with the $`M_{BH}\sigma `$ relationship. B96b suggested that the mass of the host galaxy might be involved in the increase of \[O III\] width with luminosity.

The conclusion that galaxy growth is contemporaneous with black hole growth is consistent with chemical abundances in QSOs. Heavy element abundances in luminous QSOs are solar or several times solar (Hamann et al. 2002). The highest abundances are typically found in the centers of large galaxies (e.g., Garnett et al. 1997).

We have crudely estimated bolometric luminosities for the AGN in our sample by taking $`L_{bol}=9\nu L_\nu (5100\mathrm{\AA })`$ (Kaspi et al. 2000). On average the various data sets have $`\mathrm{log}L_{bol}/L_{\mathrm{Ed}}0.4`$, where the Eddington limit $`L_{\mathrm{Ed}}`$ is calculated from $`M_{\mathrm{BH}}(H\beta )`$. This ratio shows little systematic difference between the samples at higher and lower luminosity or redshift. Some individual objects exceed the Eddington limit on this basis, but only by amounts that may be consistent with uncertainties in $`M_{BH}`$ and $`L_{\mathrm{bol}}`$.

The use of $`\sigma _{[OIII]}`$ provides a new way to estimate black hole masses in AGN. Once the systematics of the \[O III\] profile are understood in this context, black hole masses from $`\sigma _{[OIII]}`$ may be as reliable as “photoionization masses” based on H$`\beta `$ width and continuum luminosity. The availability of two independent estimates of $`M_{BH}`$ will allow workers to compare each for consistency with other approaches, such as accretion disk fitting of the continuum energy distribution (e.g., Mathur et al. 2001). Both measures of $`M_{BH}`$ may likewise be examined for systematic correlations with various properties of AGN, such as the correlation of $`L/L_{\mathrm{Ed}}`$ with “eigenvector 1” found by Boroson (2002).

The $`M_{BH}\sigma _{}`$ relationship offers, in principle, a new standard candle for cosmology. The value of $`M_{BH}`$ derived from the broad H$`\beta `$ line width varies as $`L^a`$, with $`a0.5`$ to 0.6 (see discussion preceding equation 2). If one assumes that either $`\sigma _{}`$ or $`\sigma _{[OIII]}`$ correctly predicts $`M_{BH}`$ through equation (1), then equation (2) can be solved for $`L`$ in terms of $`\sigma _{}`$ and FWHM(H$`\beta `$). This gives

$$L_{44}=10^{0.88}\sigma _{200}^{8.02}v_{3000}^4.$$
(3)
However, the use of equation (3) for cosmological measures will be challenging, given the large exponents of $`\sigma _{200}`$ and $`v_{3000}`$ and the uncertainties in the calibration of the $`M_{BH}\sigma _{[OIII]}`$ and $`M_{\mathrm{BH}}Ł`$ relationships.

Future work on this topic should include a systematic examination of the best way of characterizing the width of the \[O III\] lines for use in predicting $`\sigma _{}`$ and calibration of the chosen measure. Use of shorter wavelength narrow lines such as \[O II\] would allow study of higher redshift objects at a given observed wavelength. An obvious need is for more and better observations of the H$`\beta `$ and \[O III\] region in high redshift QSOs. In order to minimize the uncertainty associated with the scaling of BLR radius with luminosity, observations of high and low redshift QSOs should be made as nearly as possible at similar intrinsic luminosity. High resolution imaging studies may provide direct measures of the host galaxy luminosity as a check on $`\sigma _{[OIII]}`$.

G.A.S. gratefully acknowledges the hospitality of Lick Observatory during summer 2001 and 2002. This material is based in part upon work supported by the Texas Advanced Research Program under Grant No. 003658-0177-2001. B.J.W. acknowledges financial support through NASA Long Term Space Astrophysics grant NAG5-3431. M.D. acknowledges support through NASA Long Term Space Astrophysics grant NAG5-3234 and NSF research grant AST99-84040. We are grateful to G. Blumenthal, T. Boroson, A. Gonçalves, D. Grupe, E. Hooper, J. Kormendy, A. Laor, J. Miller, and D. Osterbrock for helpful discussions.


text:
# NTUTH-02-505F August 2002 A computational system for lattice QCD with overlap Dirac quarksThis work was supported in part by National Science Council R.O.C. under the grant number NSC90-2112-M002-021

## 1 Introduction

It is well known that extracting physics from lattice QCD requires computing power exceeding that of any desktop personal computer currently available in the market. Therefore, for one without supercomputer resources, building a computational system seems to be inevitable if one really wishes to pursue a meaningful number of any physical quantity from lattice QCD. However, the feasibility of such a project depends not only on the funding, but also on the theoretical advancement of the subject, namely, the realization of exact chirally symmetry on the lattice . Now, if we also take into account of the current price/performance of PC hardware components, it seems to be the right timing to rejuvenate the project  with a new goal - to build a computational system for lattice QCD with exact chiral symmetry. In this paper, we outline the essential features of a Linux PC cluster which is now being developed at National Taiwan University, and discuss how to optimize its hardware and software for lattice QCD with overlap Dirac quarks. More detailed descriptions have been given in Ref. .

First, we start from quenched QCD, and compute quark propagators in the gluon field background, for a sequence of configurations generated stochastically with weight $`\mathrm{exp}(𝒜_g)`$. Then the hardronic observables such as meson and baryon correlation functions can be constructed, and from which the hadron masses and decay constants can be extracted. In general, one requires that any quark propagator coupling to physical hadrons must be of the form

$`(D_c+m_q)^1,`$ (1)
where $`m_q`$ is the bare quark mass, and $`D_c`$ is a chirally symmetric and anti-hermitian Dirac operator111Here we assume that $`D_c`$ is doubler-free, and has correct continuum behavior, and $`D=D_c(1+raD_c)^1`$ is exponentially local for a range of $`r>0`$.. For any massless lattice Dirac operator $`D`$ satisfying the Ginsparg-Wilson relation

$`D\gamma _5+\gamma _5D=2raD\gamma _5D,`$ (2)
it can be written as $`D=D_c(1+raD_c)^1`$ , and the bare quark mass is naturally added to the $`D_c`$ in the numerator ,

$`D(m_q)=(D_c+m_q)(1+raD_c)^1.`$
Then the quenched quark propagator becomes

$`(D_c+m_q)^1=(1rm_qa)^1[D(m_q)^1ra].`$
If we fix one of the end points at $`(\stackrel{}{0},0)`$ and use the Hermitcity $`D^{}=\gamma _5D\gamma _5`$, then only 12 (3 colors times 4 Dirac indices) columns of

$`D(m_q)^1=D^{}(m_q)\{D(m_q)D^{}(m_q)\}^1`$ (3)
are needed for computing the time correlation functions of hadrons. Now our problem is how to optimize a PC cluster to compute $`D(m_q)^1`$ for a set of bare quark masses.

For overlap Dirac quarks, we need to solve the following linear system by (multi-mass) conjugate gradient (CG),

$`DD^{}Y=\left(\alpha +\eta P_\pm H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}\right)Y=\text{1I}`$ (4)
where $`P_\pm =(\gamma _5\pm 1)/2`$, $`\alpha =2m_0^2+m_q^2/2`$, $`\eta =2m_0^2m_q^2/2`$, and $`m_0`$ is fixed to $`1.3`$ in our computations. Then the quark propagators can be obtained via (3). With Zolotarev optimal rational approximation  to $`(H_w^2)^{1/2}`$, the multiplication ($`h_wH_w/\lambda _{min}`$)

$`H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}Yh_w(h_w^2+c_{2n}){\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \frac{b_l}{h_w^2+c_{2l1}}}Y`$ (5)
can be evaluated by invoking another (multi-mass) conjugate gradient process to the linear systems

$`(h_w^2+c_{2l1})Z_l=Y,l=1,\mathrm{},n.`$ (6)
where the coefficients $`c_l`$, $`b_l`$ and $`d_0`$ have been given explicitly in Ref. .

In order to improve the accuracy of the rational approximation as well as to reduce the number of iterations in the inner CG loop, it is crucial to project out the largest and some low-lying eigenmodes of $`H_w^2`$. Denoting these eigenmodes by

$`H_wu_j=\lambda _ju_j,j=1,\mathrm{},k,`$ (7)
then we project the linear systems (6) to the complement of the vector space spanned by these eigenmodes

$`(h_w^2+c_{2l1})\overline{Z_l}=\overline{Y}(1{\displaystyle \underset{j=1}{\overset{k}{}}}u_ju_j^{})Y.`$ (8)
Thus the r.h.s. of (5) can be rewritten as

$`h_w(h_w^2+c_{2n}){\displaystyle \underset{l=1}{\overset{n}{}}}b_l\overline{Z}_l+{\displaystyle \underset{j=1}{\overset{k}{}}}{\displaystyle \frac{\lambda _j}{\sqrt{\lambda _j^2}}}u_ju_j^{}YS`$ (9)
Then the breaking of exact chiral symmetry (2) can be measured in terms of

$`\sigma ={\displaystyle \frac{|S^{}SY^{}Y|}{Y^{}Y}}.`$ (10)
In practice, one has no difficulties to attain $`\sigma <10^{12}`$ for most gauge configurations on a finite lattice (Table 1).

Now the computation of overlap Dirac quark propagator involves two nested conjugate gradient loops: the so-called inner CG loop (8), and the outer CG loop (4).

## 2 Optimization

Next we address the question how to configure the hardware and software of a PC cluster such that it can attain the optimal price/performance ratio for the execution of the nested CG loops.

The vital observation is that not all column vectors are used simultaneously at any step of the nested CG loops, and also the computationally intense part is at the inner CG loop (8). Thus we can use the hard disk as virtual memory for the storage of the intermediate solution vectors and their conjugate gradient vectors at each iteration of the outer CG loop, while the CPU is working on the inner CG loop. Then the amount of required memory at each node can be greatly reduced. It is easy to figure out  the minimum memory required for accommodating the link variables and all relevant vectors for the inner CG loop,

$`N_{cg}^{min}=384\times N_s(n+3)\text{byte},`$ (11)
where $`n`$ is the degree of the Zolotarev rational polynomial, and $`N_s`$ is the total number of lattice sites. For computations of quark propagators on the $`16^3\times 32`$ lattice with $`n=16`$, (11) gives $`N_{cg}^{min}=0.912\text{Gbyte}`$, which can be implemented in a single node with one Gbyte of memory (i.e., four stripes of 256 Mbyte memory modules). Then the maximum speed of the cluster can be attained since there is no communication overheads. Moreover, the time for disk I/O (at the interface of inner and outer CG loops) only constitutes a few percent of the total time for the entire nested CG loops (Table 1). Further, to take advantage of the vector unit of Pentium 4, we rewrite the computationally intense part ($`H_w`$ times a vector $`Y`$) of our program in SSE2 codes, which yields a speed-up by a factor of $`1.8`$.

## 3 Conclusions

The speed of our system of 30 nodes is higher than 30 Gflops, and the total cost of the hardware is less than US$30,000. This amounts to price/performance ratio better than $1.0/Mflops for 64-bit (double precision) computations with overlap Dirac quarks. The basic idea of optimization is to let each node compute one of the 12 columns of the quark propagators (for a set of bare quark masses), and also use the hard disk as virtual memory for the vectors in the outer CG loop, while the CPU is working on the inner CG loop.

With Zolotarev optimal rational approximation to $`(H_w^2)^{1/2}`$, projection of high and low-lying eigenmodes of $`H_w^2`$, the multi-mass CG algorithm, the SSE2 acceleration, and the memory management, we are able to compute overlap Dirac quark propagators of $`16`$ bare quark masses on the $`16^3\times 32`$ lattice, with the precision of quark propagators up to $`10^{11}`$ and the precision of exact chiral symmetry up to $`10^{12}`$, at the rate of two gauge configurations ( $`\beta =6.0`$ ) per two days, with our present system of 30 nodes (Table 1). This demonstrates that an optimized Linux PC cluster can be a viable computational system to extract physical quantities from lattice QCD with overlap Dirac quarks .


text:
# Logarithmic rate dependence in deforming granular materials

(Oct. 8, 2002)

## Abstract

Rate-independence for stresses within a granular material is a basic tenet of many models for slow dense granular flowswood90 ; nedderman92 ; jaeger96 ; herrmann97 ; behringer97 . By contrast, logarithmic rate dependence of stresses is found in solid-on-solid frictionovarlez01 ; nasuno98 ; losert00 , in geological settingsdieterich79 ; ruina83 , and elsewhereberthoud99 ; baumberger99 ; carlson96 ; gnecco00 . In this work, we show that logarithmic rate-dependence occurs in granular materials for plastic (irreversible) deformations that occur during shearing but not for elastic (reversible) deformations, such as those that occur under moderate repetitive compression. Increasing the shearing rate, $`\mathrm{\Omega }`$, leads to an increase in the stress and the stress fluctuations that at least qualitatively resemble what occurs due to an increase in the density. Increases in $`\mathrm{\Omega }`$ also lead to qualitative changes in the distributions of stress build-up and relaxation events. If shearing is stopped at $`t=0`$ , stress relaxations occur with $`\sigma (t)/\sigma (t=0)A\mathrm{log}(t/t_o)`$ . This collective relaxation of the stress network over logarithmically long times provides a mechanism for rate-dependent strengthening.

Slow granular flows, the subject of this letter, are typically described in the context of Mohr-Coulomb friction modelsnedderman92  that resemble those used for describing friction between two solid bodiesheslot94 . In the well-known solid friction scenariopersson98 ; meyer98 , an object on a frictional surface will resist a force and remain at rest provided the magnitude of the tangential force is less than the product of a static friction coefficient and the normal force. This picture was translated into the granular context (for dense granular systems characterized by networks of force chainshowell99a ) by Coulombcoulomb76  and more recent authorswood90 ; nedderman92 , where the normal and tangential forces are replaced by corresponding normal and shear stresses, and the surface of interaction is replaced by a plane within the material. For large enough tangential force (shear stress) relative to the normal force (normal stress) sliding friction (failure in a granular material) occurs. In these pictures, sliding friction (deformation following failure) is independent of the speed of sliding (the shear rate).

In reality, experiments in diverse contexts have shown that solid friction exhibits a logarithmic dependence on rate and that static frictional contacts strengthen logarithmically with age. These experiments span a vast range of lengths, and include studies at the atomicgnecco00 , labnasuno98 ; losert00 ; berthoud99 ; baumberger99 ; carlson96 , and geological scaledieterich79 ; ruina83 . Recent experiments by Ovarlez *et al.*ovarlez01  using granular materials sliding against the interior wall of a piston showed clear rate dependence that was associated with aging effects of individual solid-friction contacts and with the force network. Additionally, experiments by Nasuno *et al.*nasuno98 ; losert00  in which a solid surface was pushed across a granular bed showed a slow strengthening (or aging) with time of the (quasi-static) force. In the present experiments, which zoom in uniquely on the grains, we show that there is a logarithmic rate-dependence in slowly sheared granular materials that is associated with irreversible rearrangements of the grain contacts within the material itself. This is manifest as a strengthening with rate, and it is unique to granular systems and possibly other jammed systemsliu98 ; i.e., it does not depend *per se* on the logarithmic strengthening of frictional contacts.

The experiments described here were carried out with a 2D realisation of a granular system. The particles were made of a photoelastic material and were either relatively thin disks or flat particles with a pentagonal cross section. By using a photoelastic material, we could determine forces at the grain scale. More detailed descriptions of the experimental apparatus and methodsheywood52  have been given elsewhere for related experimentshowell99a ; howell99b ; geng01 , and here we provide only essential information for two qualitatively different experiments. For each experiment, the humidity varied by no more than $`\pm 5\%`$. In the first of these two experiments, the photoelastic particles were sheared in an annular geometry (inset, Fig. 2) that was bounded on the inside by a rough wheel, the source of shearing, and on the outside by an equally rough static ring. The particles sat on a horizontal, flat powder-lubricated Perspex sheet, and were contained from above by a similar sheet. The whole was placed in a polariscope consisting of a pair of right and left circular polarizers, a light source and a video camera. Photoelastic images were obtained at rates up to 15 Hz, captured by a framegrabber, and processed on-the-fly to determine the forces on the particles within the field of view of the camera ($`100`$ to $`200`$ particles). Specifically, we computed the force on each particle in each frame and summed over all particles in a frame. We stored this integrated force for each image to produce a time series over very long sampling times, between $`2^{14}`$ and $`2^{16}`$ points depending on the rate of shearing. To a reasonable approximation ($`10`$%), this quantity is proportional to the pressure, and for simplicity, we will refer to it as the stress.

Key control parameters are the rate of shearing, $`\mathrm{\Omega }=2\pi /T`$, where $`T`$ is the period for one rotation of the shearing wheel, and the density of the system, which we give in terms of the area fraction, $`\gamma `$, occupied by the particles. As shown previouslyhowell99a ; howell99b , there is a critical value of $`\gamma _c`$, below which the shearing ceases, and we typically reference $`\gamma `$ to this value. ($`\gamma _c`$ depends on properties such as the particle shape and possibly the geometry of the container.) It is also useful to give the rotation rate, $`f=1/T=\mathrm{\Omega }/2\pi `$. In these experiments, $`0.03\text{ mHz}f60\text{ mHz}`$.

In Fig. 1, we show several examples of stress time series for $`\mathrm{\Omega }`$’s ranging over about three orders of magnitude in the shearing rate. The three images of this figure give an indication of the effect of shearing rate: the time series for the slowest rate is clearly of lower amplitude and more intermittent than that for the highest rate.

The rate dependence is seen clearly in the mean (time-averaged) stress vs. $`\mathrm{\Omega }`$, Fig. 2. These data are consistent with a logarithmic rate dependence for all $`\gamma `$’s. Similar effects have been seen in solid-on-solid friction experiments but with a significant difference: in the present experiments, there is a logarithmic strengthening of the mean stress with rate, while other experimentsovarlez01 ; losert00 ; berthoud99  show a logarithmic strengthening with waiting time and/or humidity. If we fit these data to the form $`\overline{\sigma }=A\mathrm{log}(\mathrm{\Omega }/\mathrm{\Omega }_o)`$, then the amplitude, $`A`$, is sensitive to the density, but the characteristic frequency, $`\mathrm{\Omega }_0`$, depends only weakly on $`\gamma `$, and corresponds to a typical time of $`2\pi /\mathrm{\Omega }_o30\text{ days}`$. However, there is no compelling reason to believe that this relation holds in the limit $`\mathrm{\Omega }0`$. A more realistic description may be $`\overline{\sigma }=A\mathrm{log}[(\mathrm{\Omega }+\mathrm{\Omega }_1)/\mathrm{\Omega }_o]`$, where $`\mathrm{\Omega }_1>\mathrm{\Omega }_o`$.

If a sample is sheared under steady state conditions and the shearing wheel is abruptly stopped, the stress relaxes over very long time scales. In Fig. 3, we show the relaxation of stress vs. time on a semi-log plot. Note that there are two time scales: a quick initial decay lasting $`20\text{ s}`$ and a much longer process such that the stress network may still be relaxing after 20 hrs. Fig. 3 shows fits to the stress decrease with time, of the form $`\sigma (t=0)/\sigma (t)A\mathrm{log}(t/t_0)`$. Although the relaxation is not completely uniform in time, a logarithmic functional form is not unreasonable for many cases (after an initial rapid relaxation). This was not always the case, and an interesting example of a major relaxation event is seen in one of these runs about 350 s after the stopping of shearing; thereafter, continued slow relaxation occurred. A key insight is that, while the micro-contacts between particles may be strengthening in timelosert00 ; berthoud99 ; heslot94 , the stress network as a whole is relaxing over very long time scales.

We gain additional insight into the effect of rate changes by determining the distribution of build-up events (monotonic increases in stress) and release events (”avalanches” – monotonic decreases in stress). Here, there are two aspects that are of interest: the size of the build ups and avalanches, $`\mathrm{\Delta }\sigma `$, and the angle, $`\theta `$, through which the shearing wheel turns during the course of a build-up or an avalanche. In the top part of Fig. 4, we show an example of an ”avalanche” event, and in the bottom, we show distributions for $`\mathrm{\Delta }\sigma `$ and $`\theta `$ for avalanches. Note that for near-critical $`\gamma `$’s and for the slowest $`\mathrm{\Omega }`$’s, the avalanche distributions are consistent with power laws. At higher $`\gamma `$ and $`\mathrm{\Omega }`$ the distributions for both build-up (not shown due to space constraints) and avalanche events are roughly exponentials. Interestingly, both increases in $`\mathrm{\Omega }`$ and in $`\gamma `$ tend to have the same qualitative effect on the statistics of avalanches and build-ups.

In order to further explore the origin of this rate dependence, we have carried out a second set of experiments, using the same photoelastic techniques, in which the particles were confined to a box with three fixed sides. The sample was then compressed and released repeatedly on the fourth side by means of an oscillating piston, as in the sketch of the inset of Fig. 5. The driving velocity of the piston was chosen to match the corresponding velocities of the shear wheel in the first experiment. The key difference between the two sets of experiments was that in the first set there was continuous plastic deformation while in the latter set there was none; i.e., the particles maintained their relative position. As shown in the body of the figure, which gives $`\overline{\sigma }`$ vs. $`\mathrm{\Omega }`$, there was no rate dependence of this process, within experimental error.

To conclude, we have demonstrated that sheared granular materials exhibit logarithmic rate dependence that is tied in an essential way to particle rearrangements (plastic deformation). The nature of this rate dependence is further illuminated by the slow relaxation that occurs when shearing is halted. The latter experiments in particular suggest that there are slow collective rearrangements of the particles that can occur over large time scales. This slow relaxation of stresses over long times is consistent with increases in with rate: force chains generated by shearing cannot completely relax on the time scales over which new chains are formed, an effect that is exacerbated by increasing $`\mathrm{\Omega }`$. It seems likely that this effect is intimately related to the long time scales reported for compaction in granular materialsnowak98 . For compaction, collective rearrangements become progressively more difficult over time because each new rearrangement requires the involvement of larger collections of particles. Similarly, in the present experiments, relaxation events involving collections of particles become less probable over time due to geometric effects and to the fact that the available elastic energy is gradually reduced with each rearrangement. An interesting question is whether similar properties are seen in other jammed systemsliu98  such as colloids or foams.

###### Acknowledgements.

Acknowlegements: This work was supported by NSF Division of Materials Research, NSF Division of Materials Science & by NASA. We appreciate helpful conversations with Prof. Sue Coppersmith.