text:
# Comment on: “High-Pressure Elasticity of 𝛼-Quartz: Instability and Ferroelastic Transition”

Gregoryanz et al.  have analyzed the stability of $`\alpha `$-quartz under pressure in terms of the Born stability criterion (BSC) . The BSC requires a positive definite elastic matrix constant $`C_{ij}`$ for stable elastic solids. It is strictly valid only at zero external pressure. Important differences were discussed with theoretical predictions by Binggeli and Chelikowsky : The pressure regime in which $`\alpha `$-quartz could be expected to be mechanically stable if it was possible to suppress the quartz I - quartz II transition at $`p=22`$ GPa was found to be significantly larger than predicted theoretically. This conclusion is incorrect, because Gregoryanz et al. use an improper generalization of BSC to non-zero pressures. The correct criterion is to require that the matrix of Birch or stiffness coefficients  $`B_{ij}`$ is positive definite rather than $`C_{ij}`$ . Binggeli and Chelikowsky also state their BSC in terms of $`C_{ij}`$, but as it seems to us, they accidentally calculated $`B_{ij}`$, e.g., they might have inverted the matrix $`V\delta ϵ_i\delta ϵ_j/k_BT`$ with $`\delta ϵ_i`$ being the thermal fluctuation of the strain tensor around the equilibrium state at a given externally applied stress. This inversion gives the Birch coefficients rather than $`C_{ij}`$; the equality $`B_{ij}=C_{ij}`$ only holds at zero external stress. As we will show now by unpublished molecular dynamics simulations of the present authors, the predictions of Binggeli and Chelikowsky are correct. We use the potential energy surface from Beest et al. , which is known to be particularly well suited to describe (elastic) properties of $`\alpha `$-quartz. We focus on the BSC

$$B_3(C_{ij})=(C_{11}C_{12})C_{44}2C_{14}^2>0,$$
(1)
which is shown in Fig. 1.

In the inset of Fig. 1 it can be seen that the present simulations reproduce the experimentally measured pressure dependence of $`B_3(C_{ij})`$ (the central Fig. 3c of Ref. ) fairly well. We also want to note that we obtain the quartz I quartz II transition at the same pressure $`p\pm 0.5`$ GPa as in experiment. These agreements justify the use of the BKS potential for $`\alpha `$-quartz under pressure. If $`B_{\alpha \beta }`$ replace $`C_{\alpha \beta }`$ in Eq. 1, the behavior is similar to the one reported by Binggeli and Chelikowsky , e.g., a nearly vanishing slope of $`B_3(B_{ij})`$ at small pressures is found, and instability of $`\alpha `$ quartz can be expected to occur at $`p25`$ GPa. This is the correct generalization of BSC, which is equivalent to require that all eigenvalues of $`B_{ij}`$ are positive. The minimum eigenvalue of $`B_{ij}`$ as obtained in our simulations is shown in the main part of Fig. 1.

We further want to point out that the observation of $`C_{44}/p<0`$, which is discussed explicitly in Ref. , has no direct implication on the mechanical stability. The only relevant influence of $`C_{44}`$ on mechanical stability in the case of an isotropic pressure $`p`$ enters via the properly generalized Born stability criterion $`B_3(B_{ij})`$, see Eq. (1), with $`B_{44}=C_{44}p`$.

We would like to thank K. Binder for useful discussions.

M. H. Müser and P. Schöffel
Inst. für Physik, WA 331
Johannes Gutenberg-Universität
55099 Mainz, Germany


format: text

text:
# Nanowire Gold Chains: Formation Mechanisms and Conductance

## FIGURES

FIG. 1. (a) Energies ($`\mathrm{\Delta }E`$, in unit of eV plotted at the bottom, relative to the initial ds configuration at $`L=11.5`$ Å) and conductances ($`G`$, in unit of $`g_0`$, plotted at the top) of Au NWs versus the distance between the outermost layers of the opposing electrodes, $`L`$; the double-strand (ds, squares), bent-chain (bc, triangles), and linear-chain (lc, circles) intervals are marked at the top. The values for the rhombic (rh) isomer at $`L=13.3`$ Å and for the zig-zag (zz) isomers at $`L=14.8`$ Å and 15.8 Å are marked by crosses, and those for the end-dimer configurations at $`L=17.3`$ Å, 17.6 Å and 18 Å are depicted by stars. The elongation force between consecutive structures may be estimated by the corresponding slope $`\mathrm{\Delta }E/\mathrm{\Delta }L`$. The atomic configurations (with atom indices) shown as insets to (a) correspond to the ds (left, bottom), bc (left, top), and lc (right) structures at $`L=12.6`$ Å, 13.3 Å and 16.3 Å respectively. (b) atomic configurations for the rh isomer at $`L=13.3`$ Å ($`d(1,2)=2.76`$ Å, $`d(2,3)=2.59`$ Å, and $`\mathrm{}(1,2,4)=124^{}`$), and for the zz isomer at $`L=14.8`$ Å ($`d(1,2)=d(3,4)=2.55`$ Å, $`d(2,3)=2.52`$ Å and $`\mathrm{}(1,2,3)=119^{}`$).

FIG. 2 (color) (a) Iso-surface images of the magnitudes of the first ($`0.97g_0`$, top) and second ($`0.7g_0`$, bottom) transmitted eigenchannels for the highly-strained ds configuration (at $`L=12.6`$ Å). (b and c) HOMO (top) and HOMO(-1) (bottom) orbitals corresponding to the highly-strained ds NW (in b) and to the bc configuration at $`L=13.3`$ Å (in c). (d) HOMO (top) and HOMO(-2) (bottom) orbitals for the ed configuration at $`L=17.3`$ Å. In (a-d) Au atoms are depicted by red spheres. In (a) we display the entire atomic system including parts of the jellium-slab continuations of the electrodes (used in the conductance calculations), and in (b-d) we focus on the vicinity of the nanowires. The viewing orientation was chosen in each case to enhance visualization of the orbital characters on the wire atoms.


text:
# Untitled Document

##

Introduction. Fluid dynamical models are widely used to describe ultra-relativistic heavy ion collisions. Their advantage is that one can vary flexibly the Equation of State (EoS) of the matter and test its consequences on the reaction dynamics and outcome. In energetic collisions of large heavy ions, especially if Quark-Gluon Plasma (QGP) is formed in the collision, one-fluid dynamics is a valid and good description for the intermediate stages of the reaction. Here, interactions are strong and frequent, so that other models, (e.g. transport models, string models, etc., assuming binary collisions, with free propagation of constituents between collisions) have limited validity. On the other hand, the initial and final, Freeze-Out (FO), stages of the reaction are outside the domain of applicability of the fluid dynamical model.

In conclusion, the realistic, and detailed description of an energetic heavy ion reaction requires a Multi Module Model, where the different stages of the reaction are each described with suitable theoretical approaches. It is important that these Modules are coupled to each other correctly: on the interface, which is a 3 dimensional hyper-surface in space-time with normal $`d\sigma ^\mu `$, all conservation laws should be satisfied (e.g. $`[T^{\mu \nu }d\sigma _\nu ]=0`$), and entropy should not decrease, $`[S^\mu d\sigma _\mu ]0`$. These matching conditions were worked out and studied for the matching at FO in detail in refs. .

The final FO stages of the reaction, after hadronization, can be described well with kinetic models where the matter is already dilute.

The initial stages are more problematic. Frequently two or three fluid models are used to remedy the difficulties, and to model the process of QGP formation and thermalization.  Here, the problem is transferred to the determination of drag-, friction- and transfer- terms among the fluid components, and a new problem is introduced with the (unjustified) use of EoS in each component in a nonequilibrated situations, where EoS does not exist. Strictly speaking this approach can only be justified for mixtures of noninteracting ideal gas components. Similarly, the use of transport theoretical approaches assuming dilute gases with binary interactions is questionable, as due to the extreme Lorentz contraction, in the C.M. frame, enormous particle and energy densities, with the immediate formation of perturbative vacuum should be handled. Even in most parton cascade models these initial stages of the dynamics are just assumed in form of some initial condition, with little justification behind.

Our goal in the present work is to construct a model, based on the recent experiences gained in string Monte Carlo models and in parton cascades. One important conclusion of heavy ion research in the last decade is that standard ’hadronic’ string models fail to describe heavy ion experiments.

All string models had to introduce new, energetic objects: string ropes , quark clusters , fused strings , in order to describe the abundant formation of massive particles like strange antibaryons. Based on this, we describe the initial moments of the reaction in the framework of classical (or coherent) Yang-Mills theory, following ref.  assuming larger field strength (string tension) than in ordinary hadron-hadron collisions. In addition we now satisfy all conservation laws exactly, while in ref.  infinite projectile energy was assumed, and so, overall energy and momentum conservation was irrelevant. We do not solve simultaneously the kinetic problem leading to parton equilibration, but assume that the arising friction is such that the heavy ion system will be an overdamped oscillator, i.e. yo-yoing of the two heavy ions will not occur. This assumption is based on recent string and parton cascade results.

##

Formulation of model. Our basic idea is to generalize the model developed in , for collisions of two heavy ions and improve it by strictly satisfying conservation laws. First of all, we would create a grid in $`[x,y]`$ plane ($`z`$ – is the beam axes, $`[z,x]`$ – is reaction plane). We will describe the nucleus-nucleus collision in terms of steak-by-streak collisions, corresponding to the same transverse coordinates, $`\{x_i,y_j\}`$. We assume that baryon recoil for both target and projectile arise from the acceleration of partons in an effective field $`F^{\mu \nu }`$, produced in the interaction. Of course, the physical picture behind this model should be based on chromoelectric flux tube or string models, but for our purpose we consider $`F^{\mu \nu }`$ as an effective abelian field. Phenomenological parameters describing this field must be fixed from comparison with experimental data.

Let describe the streak-streak collision.

$$_\mu \underset{i}{}T_i^{\mu \nu }=\underset{i}{}F_i^{\nu \mu }n_{i\mu },$$
(1)
$$_\mu \underset{i}{}n_i^\mu =0,i=1,2,$$
(2)
$`n_i^\mu `$ is the baryon current of $`i`$th nucleus (we are working in the Center of Rapidity Frame (CRF), which is the same for all streaks. The concept of using target and projectile reference frames has no advantage any more). We will use the parameterization:

$$n_i^\mu =\rho _iu_i^\mu ,u_i^\mu =(\mathrm{cosh}y_i,\mathrm{sinh}y_i).$$
(3)
$`T^{\mu \nu }`$ is a energy-momentum flux tensor. It consists of five parts, corresponding to both nuclei and free field energy (also divided into two parts) and one defines the QGP perturbative vacuum.

$$T^{\mu \nu }=\underset{i}{}T_i^{\mu \nu }+T_{pert}^{\mu \nu }=\underset{i}{}\left[e_i\left(\left(1+c_0^2\right)u_i^\mu u_i^\nu c_0^2g^{\mu \nu }\right)+T_{F,i}^{\mu \nu }\right]+Bg^{\mu \nu },i=1,2.$$
(4)
$`B`$ – is the bag constant, the equation of state is $`P_i=c_0^2e_i`$, where $`e_i`$ and $`P_i`$ are energy density and pressure of QGP.

In complete analogy to electro-magnetic field

$$F_i^{\mu \nu }=^\nu A_i^\mu ^\mu A_i^\nu =\left(\begin{array}{cc}0& \sigma _i\\ \sigma _i& 0\end{array}\right),$$
(5)
$$\sigma _i=^3A_i^0^0A_i^3,$$
(6)
$`T_{F,i\mu \nu }=g_{\mu \nu }_{F,i}+{\displaystyle \underset{\beta }{}}{\displaystyle \frac{_{F,i}}{\left(^\mu A_i^\beta \right)}}_\nu A_i^\beta ,`$ (7)
$$_{F,i}=\frac{1}{4}F_{i\mu \nu }F_i^{\mu \nu }.$$
(8)
In our case the string tensions, $`\sigma _i`$, will have the same absolute value $`\sigma `$ and opposite sign (in complete analogy to the usual string with two ends moving in opposite directions), and $`\sigma _i`$ will be constant in the space-time region after string creation and before string decay.

To get the analytic solutions of the above equations, we use light cone variables

$$(z,t)(x^+,x^{}),x^\pm =t\pm z.$$
(9)
Following , we insist that $`e_1,y_1,\rho _1,A_1^\mu `$ are functions of $`x^{}`$ only and $`e_2,y_2,\rho _2,A_2^\mu `$ depend on $`x^+`$ only.

In terms of light cone variables:

$$n_i^\pm =n_{i,}=\rho _i(u_i^0\pm u_i^3)=\rho _ie^{\pm y_i},$$
(10)
$$\left(\begin{array}{cc}T_i^{++}& T_i^+\\ T_i^+& T_i^{}\end{array}\right)=\left(\begin{array}{cc}h_{i+}e^{2y_i}& h_i\\ h_i& h_{i+}e^{2y_i}\end{array}\right)+T_{F,i},$$
(11)
where

$$h_{i+}=(1+c_0^2)e_i,h_i=(1c_0^2)e_i.$$
(12)
$$\left(\begin{array}{cc}F_i^{++}& F_i^+\\ F_i^+& F_i^{}\end{array}\right)=\left(\begin{array}{cc}0& 2\sigma _i\\ 2\sigma _i& 0\end{array}\right).$$
(13)
$$T_{pert}=\left(\begin{array}{cc}0& 2B\\ 2B& 0\end{array}\right).$$
(14)
At the time of first touch of two streaks, $`t=0`$, there is no string tension. We assume that strings are created, i.e. the sting tension achieves the value $`\sigma `$ at time $`t=t_0`$, corresponding to complete penetration of streaks through each other.

##

Conservation laws — String rope creation. In light cone variables eq. (2) may be rewritten as

$$_{}n_1^{}+_+n_2^+=0.$$
(15)
So, we have a sum of two terms, depending on different independent variables, and the solution can be found in the following way.

$$\begin{array}{cc}_{}n_1^{}=a,\hfill & _+n_2^+=a,\hfill \\ n_1^{}=ax^{}+(n_1)_0,\hfill & n_2^+=ax^++(n_2)_0.\hfill \end{array}$$
(16)
Since both $`n_1^{}`$ and $`n_2^+`$ are positive (and also more or less symmetric) we can conclude that for our case $`a=0`$.

Finally

$$n_1^{}=\rho _1e^{y_1}=\rho _0e^{y_0},n_2^+=\rho _2e^{y_2}=\rho _0e^{y_0},$$
(17)
$$\rho _1=\rho _0e^{y_0+y_1},\rho _2=\rho _0e^{y_0y_2}.$$
(18)
Let us come back to the energy-momentum tensor $`T^{\mu \nu }`$. Based on eqs. (6, 7, 8) and taking into account the Lorentz gauge, $`^0A_i^0^3A_i^3=0`$, we can find

$$\left(\begin{array}{cc}T_{F,i}^{++}& T_{F,i}^+\\ T_{F,i}^+& T_{F,i}^{}\end{array}\right)=\left(\begin{array}{cc}\sigma _i^22\sigma _i\left(^0A_i^0\right)+2\sigma _i\left(^3A_i^0\right)& 0\\ 0& \sigma _i^2+2\sigma _i\left(^0A_i^0\right)+2\sigma _i\left(^3A_i^0\right)\end{array}\right).$$
(19)
As mentioned before, after string creation, i.e. $`t>t_0`$, and before string decay we choose the string tensions in the form:

$$\sigma _2=\sigma _1=\sigma >0.$$
(20)
To satisfy the above choice and the Lorentz gauge condition we take the vector potentials in the following form:

$$\begin{array}{cc}A_1^+=0,\hfill & A_1^{}=2\sigma x^{},\hfill \\ A_2^+=2\sigma x^+,\hfill & A_2^{}=0.\hfill \end{array}$$
(21)
In our calculations we used the parameterization:

$$\sigma =A\left(\frac{\epsilon _0}{m}\right)^2\rho _0\sqrt{l_1l_2},$$
(22)
where $`l_1,l_2`$ are the lengths of initial streaks. The typical values of $`A`$ are around $`0.050.06`$. Notice, that there is only one free parameter in parameterization (22). The typical values of $`\sigma `$ are $`815GeV/fm`$ for $`\epsilon _0=100GeV/nucl`$.

The problem with eq. (1) is that we do not know what the really conserved quantities are. Using the definition of $`F^{\mu \nu }`$, eq. (6), we can rewrite eq. (1) as

$$_\mu T^{\mu \nu }=\underset{i}{}F_i^{\mu \nu }n_{i,\mu }=\underset{i}{}\left(^\mu \left(A_i^\nu n_{i,\mu }\right)A^\nu ^\mu n_{i,\mu }^\nu \left(A_i^\mu n_{i,\mu }\right)+A^\mu ^\nu n_{i,\mu }\right).$$
(23)
The solution for $`n_1^{}`$ and $`n_2^+`$, eq. (17), shows as that second and fourth terms vanish. So, we can define new energy-momentum tensor $`\stackrel{~}{T}^{\mu \nu }`$, such that

$$_\mu \stackrel{~}{T}^{\mu \nu }=0,$$
(24)
$$\stackrel{~}{T}^{\mu \nu }=\underset{i}{}\stackrel{~}{T_i}^{\mu \nu }+T_{pert}^{\mu \nu }=\underset{i}{}\left(T_i^{\mu \nu }A_i^\nu n_i^\mu +g^{\mu \nu }A_i^\alpha n_{i\alpha }\right)+Bg^{\mu \nu }$$
(25)
Using the exact definition of $`A_i^\mu `$ – eqs. (21) – we obtain

$$\stackrel{~}{T}^{\mu \nu }=\left(\begin{array}{cc}h_{1+}e^{2y_1}+5\sigma ^2& h_1+4\sigma x^{}n_1^+\\ h_1+2\sigma x^{}n_1^+& h_{1+}e^{2y_1}+\sigma ^22\sigma x^{}n_1^{}\end{array}\right)+\left(\begin{array}{cc}h_{2+}e^{2y_2}+\sigma ^2+2\sigma x^+n_2^+& h_22\sigma x^+n_2^{}\\ h_24\sigma x^+n_2^{}& h_{2+}e^{2y_2}+5\sigma ^2\end{array}\right)$$

$$+\left(\begin{array}{cc}0& 2B\\ 2B& 0\end{array}\right).$$
(26)
Now the new conserved quantities are

$$Q_0=\stackrel{~}{T}^{00}𝑑V=\underset{i}{}_{\mathrm{\Omega }_i}\stackrel{~}{T_i}^{00}𝑑V,$$
(27)
$$Q_3=\stackrel{~}{T}^{03}𝑑V=\underset{i}{}_{\mathrm{\Omega }_i}\stackrel{~}{T_i}^{03}𝑑V.$$
(28)
Based on conservation of $`Q_0,Q_3`$ we can calculate rapidity, energy and baryon densities at the moment $`t=t_0`$, when the string with tension $`\sigma `$ is created. These new quantities are used as initial conditions for our differential eqs. (1, 2).

$$\left(\frac{\epsilon _1(t_0)}{m}\right)=\left(\frac{\epsilon _0}{m}\right)\frac{1}{1+c_0^2}\frac{\sigma ^2}{\rho _0\epsilon _0(1+c_0^2)}\frac{l_1+5l_2}{4l_1}\frac{\sigma e^{y_0}}{8\epsilon _0(1+c_0^2)}(l_1+2l_2)\frac{B}{\rho _0\epsilon _0(1+c_0^2)}\frac{l_1+l_2}{2l_1},$$
(29)
$$\left(\frac{\epsilon _2(t_0)}{m}\right)=\left(\frac{\epsilon _0}{m}\right)\frac{1}{1+c_0^2}\frac{\sigma ^2}{\rho _0\epsilon _0(1+c_0^2)}\frac{l_2+5l_1}{4l_2}\frac{\sigma e^{y_0}}{8\epsilon _0(1+c_0^2)}(l_2+2l_1)\frac{B}{\rho _0\epsilon _0(1+c_0^2)}\frac{l_1+l_2}{2l_2}.$$
(30)
Here the $`\epsilon _i`$ is energy per nucleon in CRF. Now the proper baryon density can be found, $`\rho _i(t_0)=\rho _0\frac{\epsilon _0}{\epsilon _i(t_0)}`$, $`\gamma _i=\frac{1}{\sqrt{1v_i^2}}=\frac{\epsilon _i(t_0)}{m}`$.

For $`x^\pm >x_0`$, where $`x_0=2t_0|z(0)|`$ defines the string creation surface $`t=t_0`$, for nucleon or cell element in the position $`z=z(0)`$ at the time $`t=0`$, we should solve eqs. (24), with boundary conditions

$$\begin{array}{cc}n_1^\pm (x^{}=x_0)=\rho _0e^{y_0}\hfill & n_2^\pm (x^+=x_0)=\rho _0e^{\pm y_0}\hfill \\ & \\ h_{1+}(x^{}=x_0)=e_1(t_0)(1+c_0^2)\hfill & h_{2+}(x^+=x_0)=e_1(t_0)(1+c_0^2)\hfill \\ & \\ y_1(x^{}=x_0)=y_1(t_0)\hfill & y_2(x^+=x_0)=y_2(t_0)\hfill \\ & \\ \sigma _1(x^{}=x_0)=\sigma \hfill & \sigma _2(x^+=x_0)=\sigma ,\hfill \end{array}$$
(31)
where $`e_i(t_0)=m\rho _i(t_0)`$ – energy density in the rest frame of $`ith`$ nuclei.

Let us present the complete analytical solution in the following form

$$e^{()^{i+1}2y_i}=\frac{d_i}{b_i}+\left(\frac{d_i}{b_i}+e^{()^{i+1}2y_i(t_0)}\right)\left(1\frac{x^ix_0}{\tau _i}\right)^{\frac{b_i}{\alpha a_j}},$$
(32)
$$h_{i+}=e^{()^{i+1}2y_i}e_i(t_0)(1+c_0^2)e^{()^{i+1}2y_i(t_0)}\left(1\frac{x^ix_0}{\tau _i}\right),$$
(33)
$$\rho _i=\rho _0e^{y_0}e^{()^{i+1}y_i},$$
(34)
where $`x^1=x^{},x^2=x^+`$, $`i,j=1,2,ij`$, and using the notations

$$b_i=\alpha a_j+2\sigma \rho _0e^{y_0},$$
(35)
$$d_i=c_i2\sigma \rho _0e^{y_0}e^{()^{i+1}2y_i(t_0)},$$
(36)
$$\tau _i=\frac{e_i(t_0)(1+c_0^2)}{e^{()^{i+1}2y_i(t_0)}a_j},$$
(37)
$$a_1=c_1+2\sigma \rho _0e^{y_0}2\sigma \rho _0e^{y_0}e^{2y_1(t_0)},a_2=c_2+2\sigma \rho _0e^{y_0}2\sigma \rho _0e^{y_0}e^{2y_2(t_0)},$$
(38)
$$c_i=\alpha ((1+c_0^2)e_i(t_0)e_0)/2t_0.$$
(39)
Then the trajectories of nucleons (or cell elements) for both nuclei are given by:

$$\begin{array}{c}x_1^+(x^{})=z(0)+_{t_0}^x^{}𝑑xe^{2y_1(x)}=\\ \\ z(0)\frac{d_1}{b_1}(x^{}x_0)+\left(\frac{d_1}{b_1}+e^{2y_1(t_0)}\right)\tau _1\frac{\alpha a_2}{2\sigma \rho _0e^{y_0}}\left[\left(1\frac{x^{}x_0}{\tau _1}\right)^{\frac{2\sigma \rho _0e^{y_0}}{\alpha a_2}}1\right],\end{array}$$

$$\begin{array}{c}x_2^{}(x^+)=z(0)+_{t_0}^{x^+}𝑑xe^{2y_2(x)}=\\ \\ z(0)\frac{d_2}{b_2}(x^+x_0)+\left(\frac{d_2}{b_2}+e^{2y_2(t_0)}\right)\tau _2\frac{\alpha a_1}{2\sigma \rho _0e^{y_0}}\left[\left(1\frac{x^+x_0}{\tau _2}\right)^{\frac{2\sigma \rho _0e^{y_0}}{\alpha a_1}}1\right],\end{array}$$

for nucleon or cell element in the position $`z=z(0)`$ at the time $`t=0`$.

##

Recreation of the matter. As we may see from the trajectories, eqs. (, ), nucleons (or cell domains) will keep going in the initial direction up to the time $`t=t_{i,turn}`$, then they will turn and go backwards until the two streaks again penetrate through each other and new oscillation will start. Such a motion is analogous to the ”Yo-Yo” motion in the string models. Of course, it is difficult to believe that such a process would really happen in heavy ion collisions, because of string decays, string-string interactions, interaction between streaks and other reasons, which are quite difficult to take into account. To be realistic we should stop the motion described by eqs. (, ) at some moment before the projectile and target cross again.

We assume that the final result of collisions of two streaks after stopping the string’s expansion and after its decay, is one streak with length $`l_f`$, homogeneous energy density distribution, $`e_f`$, and baryon charge distribution, $`\rho _f`$, moving like one object with rapidity $`y_f`$. We assume that this is due to string-string interactions and string decays. As it was mentioned above the typical values of the string tension, $`\sigma `$, are of the order of $`10GeV/fm`$, and these may be treated as several parallel strings. The string-string interaction will produce a kind of ”string rope” between our two streaks, which is responsible for final energy density and baryon charge homogeneous distributions. Now it is worth to mention that decay of our ”string rope” does not allow charges to remain at the ends of the final streak, as it would be if we assume full transparency.

The homogeneous distributions are the simplest assumptions, which may be modified based on experimental data. Its advantage is a simple expression for $`e_f,\rho _f,y_f`$.

The final energy density and rapidity, $`e_f`$ and $`y_f`$, may be determined from conservation laws.

$$\mathrm{cosh}^2y_f=\frac{(M^2(1+c_0^2)+2c_0^2)+\sqrt{(M^2(1+c_0^2)+2c_0^2)^2+4c_0^4(M^21)}}{2(1+c_0^2)(M^21)},$$
(40)
where we neglected $`Bl_f`$ next to $`Q_0`$ and introduced the notation $`M=(l_2+l_1)/(l_2l_1)`$,

$$e_f=\frac{\frac{Q_0}{xy}Bl_f}{((1+c_0^2)\mathrm{cosh}^2y_fc_0^2)l_f}.$$
(41)
The typical trajectory of the streak ends is presented in Fig. 1. From $`t_0`$ they move according to eqs. (, ) until they reach the rapidity $`y_i=y_f`$. Later the final string starts to move like one object with rapidity $`y_f`$.

The turning points can be found from the condition:

$$y_i=y_f,$$
(42)
which gives for $`i`$th nucleus ($`x_1=x^{},x_2=x^+`$)

$$x_{i,turn}=x_0+\tau _i\left[1\left(\frac{\frac{d_i}{b_i}+e^{()^{i+1}2y_i(t_0)}}{\frac{d_i}{b_i}+e^{()^{i+1}2y_f}}\right)^{\frac{\alpha a_j}{b_i}}\right].$$
(43)
##

Initial conditions for hydrodynamical calculations. In this section we present the results of our calculations. We are interested in the shape of QGP formed, when string expansions stop and their matter is locally equilibrated. This will be the initial state for further hydrodynamical calculations. We may see in Figs. 2, that QGP forms a tilted disk for $`b0`$. So, the direction of fastest expansion, the same as largest pressure gradient, will be in the reaction plane, but will deviate from both the beam axis and the usual transverse flow direction. So, the new flow component, called ”antiflow” or ”third flow component”, will appear in addition to the usual transverse flow component in the reaction plane. With increasing beam energy the usual transverse flow is getting weaker, while this new flow component is strengthened. The mutual effect of the usual directed transverse flow and this new ”antiflow” or ”third flow component” lead to an enhanced emission in the reaction plane. This was actually observed and widely studies earlier and referred to as ”elliptic flow”.

##

Conclusions. Based on earlier Coherent Yang-Mills field theoretical models, and introducing effective parameters based on Monte-Carlo string cascade and parton cascade model results, a simplified model is introduced to describe the pre fluid dynamical stages of heavy ion collisions at the highest SPS energies and above. The model predicts limited transparency for massive heavy ions.

Contrary to earlier expectations, — based on standard string tensions of 1 GeV/fm which lead to the Bjorken model type of initial state, — effective string tensions are introduced for collisions of massive heavy ions, as a consequence of collective effects related to QGP formation. These collective effects in central and semi central collisions lead to an effective string tension of the order of 10 GeV/fm and consequently cause much less transparency than earlier estimates. The resulting initial locally equilibrated state of matter in semi central collisions takes a rather unusual form, which can be then identified by the asymmetry of the caused collective flow. Our prediction is that this special initial state may be the cause of the recently predicted ”antiflow” or ”third flow component”.

Detailed fluid dynamical calculations as well as flow experiments at semi central impact parameters for massive heavy ions are needed at SPS and RHIC energies to connect the predicted special initial state with observables.


text:
# References

Self-organized critical topology
of stock markets
N.Vandewalle1, F.Brisbois1 and X.Tordoir2

1 GRASP, Institut de Physique B5, Université de Liège, B-4000 Liège, Belgium.

2 IPNE, B15, Université de Liège, B-4000 Liège, Belgium.

We have analyzed the cross-correlations of daily fluctuations for $`N=\mathrm{6\hspace{0.17em}358}`$ US stock prices during the year 1999. From those $`N(N1)/2`$ correlations coefficients, the Minimum Spanning Tree (MST) has been built. We have investigated the topology exhibited by the MST. Eventhough the average topological number is $`n2`$, the variance $`\sigma `$ of the topological distribution $`f(n)`$ diverges! More precisely, we have found that $`f(n)n^{2.2}`$ holding over two decades. We have studied the topological correlations for neighbouring nodes: an extremely broad set of local configurations exists, confirming the divergence of $`\sigma `$.

In a series of works , Mantegna and coworkers have studied the Minimum Spanning Tree (MST) of financial data, mainly the cross-correlations between $`N`$ different stocks. Building a MST is a common method in spin glasses for putting into evidence the primary structure of a complex system with a non-trivial dynamics . The MST puts into evidence the cooperative behaviours in the stock market.

Herein, we report some analysis of the topology exhibited by the MST. Due to Initial Public Offerings (IPOs) and mergers, it is difficult to perform some MST-analysis over long periods. The number of companies is indeed always evolving. We focussed our work on $`N=\mathrm{6\hspace{0.17em}358}`$ US companies during the year 1999. Those stocks are all traded on the Nasdaq, NYSE and AMEX places. Cross-correlations $`\rho _{ij}`$ for pairs $`ij`$ of stocks are computed as

$$\rho _{ij}=\frac{<Y_iY_j><Y_i><Y_j>}{\sqrt{(<Y_i^2><Y_i>^2)(<Y_j^2><Y_j>^2)}}$$
(1)
where the $`Y`$’s are the daily fluctuations of the log of the stock prices. The MST has been built following the Kruskal’s algorithm  in order to find out the $`N1`$ most important correlated pairs of stocks among the $`N(N1)/2=\mathrm{20\hspace{0.17em}208\hspace{0.17em}903}`$ possible pairs! Figure 1 presents three parts of the MST: AMZN (Amazon.com Inc.), YHOO (Yahoo! Inc.), KO (Coca-Cola Company) and their nearest neighbours. AMZN is connected to EBAY (eBay Inc.), CMGI (CMGI Inc.) and SOFN (SoftNet Systems Inc.). YHOO is connected to ARBA (Ariba Inc.) and also connected to CMGI such that YHOO is the next nearest neighbour of AMZN. KO is connected to PEP (PepsiCo, Inc.), CCE (Coca-Cola Enterprises, Inc.) and AVP (Avon Products Inc.). In addition to KO, AVP is connected to 42 others stocks!

One can observe a clustering of companies being in the same sector such as energy, technology, internet, transportation, food and cosmetics, etc. The clustering in various sectors, each being a branch of the MST, has been previously observed in . In addition to the sectorization, we have observed a trend for a global clustering of stocks traded on different market places. Most Nasdaq stocks are located around the QQQ ticker (Nasdaq Composite ticker), and largest stocks of NYSE are located around DIA (Dow Jones ticker).

One could distinguish three types of topological configurations for companies: (i) important nodes, (ii) links, and (iii) dangling ends. The important nodes (like AVP) should correspond to companies controlling or mediating the fluctuations of their neighbourhood. The links (like YHOO) are mediating the information (fluctuation) along the branches. The dangling ends (like PEP) are supposed to be the less influent companies, or to be the less influenced by the market fluctuations (like indices).

From the mathematical point of view, we can associate a “topological number” $`1n(N1)`$ at each node of the tree. The topological number $`n`$ is defined as the number of connections starting from the considered site, i.e. the coordination number of that site. The statistics of the tree is thus determined by the distribution of topological numbers $`f(n)`$. The mean topological number is given by

$$n=\underset{n=1}{\overset{\mathrm{}}{}}nf(n).$$
(2)
If one multiplies both sides of this equation by $`N`$, the summation becomes equivalent to count twice all the $`N1`$ connections along the tree. As a consequence, one has

$$n=2\frac{N1}{N}$$
(3)
which is approximately equal to 2 when the number of sites $`N`$ becomes large, i.e. in our case. The above mathematical property is quite general and holds for all types of trees. In the following, we will consider that $`n=2`$. Thus, the first moment $`n`$ of $`f`$ cannot give any information about the topology of trees since it does not reflect the shape of the tree. However, the second moment $`\sigma `$ can be quite useful. Let us assume a topological correlation in between neighbouring nodes. The Aboav law, a well known law in the Physics of foams , stipulates that topological correlations can be captured by the following relationship

$$m_n=A+\frac{B}{n}$$
(4)
where $`m_n`$ is the mean topological number of the nearest neighbours of a node characterized by a topological number $`n`$. $`A`$ and $`B`$ are positive constants to be determined. The form of the Aboav law implies that important nodes with large $`n`$ are more favorably connected to dangling ends. Indeed, when the coordination number becomes large and reaches $`nN1`$, which correspond to a star structure, one should obtain $`m_n=1`$. This leads to the constraint 111This constraint has a counterpart $`A=5`$ in the case of foams . $`A=1`$. For a random structure, one expects that the statistical average of the product $`nm_n`$ is similar to $`n^2`$, one has

$$n+B=nm_nn^2=\sigma ^2+4$$
(5)
and thus $`B=\sigma ^2+2`$ and the Aboav law becomes

$$m_n=1+\frac{\sigma ^2+2}{n}$$
(6)
for any kind of tree with a large number $`N`$ of nodes. The above equation is valid in the hypothesis of random trees.

Here, we have defined some mathematical tools ($`f(n)`$, $`\sigma `$ and the Aboav law) for describing the topology of trees in full generality. Let us investigate the specific topology of financial MST. The log-log plot of the topological distribution is shown in Figure 2. The distribution $`f(n)`$ is broad and looks like a power law

$$f(n)n^\alpha $$
(7)
holding over two decades with an exponent $`\alpha =2.2\pm 0.1`$. This result means that the variance of the $`f(n)`$ distribution diverges even the mean $`n`$ remains finite (see below)! In other words, nodes with high topological numbers are not so rare. For random trees, one expects an exponential decay of $`f(n)`$ and quite different statistics. The market is naturally self-organized in a coordination invariant structure, i.e. a self-organized critical structure ! One should also remark that our results have an important graphical consequence: large $`n`$ values for nodes are not so rare such that it is extremely hard to draw the entire MST.

The power law (7) and the divergence of $`\sigma `$ are not numerical artefacts. Indeed, we have checked the Aboav-like law for which the variance of the topological distribution is relevant (see above). Figure 3 presents the $`(m_n,\frac{1}{n})`$ diagram in which the Aboav law should be a line with a slope $`\sigma ^2+2`$. However, the $`N`$ values of the measured $`m_n`$ on the MST are sparsely distributed in the diagram such that no relevant linear fit can be performed. Eq.(6) is not valid expressing that statistical averages cannot be taken as for usual random topological structures. We argue that the extremely broad dispersion of local configurations in $`(m_n,\frac{1}{n})`$ expresses the divergence of $`\sigma `$.

We put into evidence the emergence of some order out of the apparent disorder of the stock markets! Our results are robust. We did the same analysis for subsets of the data and we obtained the same behaviour (a power law for $`f(n)`$ and a divergence of $`\sigma `$). Also, we did the same for the “anti-tree”, i.e. the tree for the strongest anticorrelations, and we have found similar results with the same exponent $`\alpha 2.2`$. As found in , the connections between pairs of stocks are mostly unchanged by taking other time periods. Of course, IPOs ($`NN+1`$) and mergers ($`NN1`$) would influence a little the MST but we stress that the topological properties of the MST would not dramatically change.

Acknowledgements

NV thanks the FNRS for financial support. Valuable discussion with R.N.Mantegna is gratefully acknowledged.

Any correspondance should be sent to Dr. Nicolas Vandewalle.

Figure Captions

Figure 1 — Three typical configurations of the nearest neighborhood of three companies on the MST: (a) AMZN (Amazon.com Inc.) connected to EBAY (eBay Inc.), CMGI (CMGI Inc.) and SOFN (SoftNet Systems Inc.); (b) YHOO (Yahoo! Inc.) is connected to ARBA (Ariba Inc.) and CMGI; (c) KO (Coca-Cola Company) is connected to PEP (PepsiCo, Inc.), CCE (Coca-Cola Enterprises, Inc.) and AVP (Avon Products Inc.). The topological indices $`n`$ are given.

Figure 2 — Log-log plot of the topological distribution $`f(n)`$ of the MST. The continuous line is a fit with the power law (7).

Figure 3 — The $`(m_n,\frac{1}{n})`$ diagram for the $`N=\mathrm{6\hspace{0.17em}358}`$ companies. No Aboav (linear) behaviour emerges.


text:
# A Fitting Formula for the Non-Linear Evolution of the Bispectrum

## 1 Introduction

The large-scale structure of the universe observed in galaxy surveys is thought to arise from the growth of primordial fluctuations due to gravitational instability. Making accurate predictions of gravitational clustering is essential if we are to compare theories with observations and put constraints on cosmological parameters.

Gravitational clustering is a non-linear process: when fluctuations are sufficiently small the equations of motion can be solved perturbatively, but at small enough scales one has resort to numerical simulations to accurately follow the dynamics. Non-linearity of the equations of motion implies that the mass distribution becomes non-Gaussian, even if the initial conditions are Gaussian (Peebles 1980, Fry 1984).

The power spectrum, or two-point correlation function, has become the most common statistical measure of gravitational clustering. Its importance rests on the fact that for a Gaussian random field, the power spectrum completely characterizes its statistical properties. Hamilton et al. (1991) suggested a useful way of thinking about the non-linear evolution of the power spectrum, summarized by a simple fitting formula obtained from numerical simulations. The accuracy of this procedure has been improved (e.g. Jain, Mo & White 1995; Peacock & Dodds 1996) and extensively used in the literature to predict statistical properties related to the density power spectrum.

A complete description of the non-Gaussian distribution that results from gravitational instability, however, requires knowledge of higher-order correlation functions. A number of useful statistics have emerged in the last decade to become standard probes of non-Gaussianity in large-scale structure. These include the higher-order moments of the smoothed density field (e.g. skewness and kurtosis; Peebles 1980; Juszkiewicz, Colombi & Bouchet 1993; Bernardeau 1994) and the three-point correlation function in Fourier space, the bispectrum (Fry 1984; Scoccimarro et al. 1998, hereafter SCFFHM) or in real space (Jing & Börner 1998; Frieman & Gaztañaga 1999).

A fitting formula for the non-linear evolution of one-point higher-order moments is given by extended perturbation theory (EPT, Colombi et al. 1997). The EPT prescription uses the perturbative expressions for skewness, kurtosis, and so on, but with an “effective” spectral index that becomes a function of the variance $`\sigma ^2`$ of the smoothed density field. This function is obtained numerically by fitting to measurements of the density field moments in N-body simulations. The EPT prescription is able to fit reasonably well the non-linear evolution of higher-order moments and thus the one-point probability distribution function with this single function obtained from numerical simulations.

Probing non-Gaussianity beyond one-point statistics gives precious additional information that can be used to learn about cosmology. The bispectrum is the lowest-order statistic that probes the shape of large-scale structures generated by gravitational clustering. Measurements of the bispectrum in galaxy redshift catalogs have been proposed as a way to constrain galaxy biasing (Fry 1994; Hivon et al. 1995; Matarrese, Verde & Heavens 1997) and thus break the degeneracies present in power spectrum measurements between the matter density of the universe $`\mathrm{\Omega }_m`$ and galaxy bias. Such measurements have recently been done in IRAS redshift surveys (Scoccimarro et al. 2000; Feldman et al. 2000), providing constraints on galaxy bias, $`\mathrm{\Omega }_m`$ and non-Gaussian primordial fluctuations. In addition, upcoming large redshift surveys such as 2dF and SDSS will determine the bispectrum with unprecedented accuracy.

In this work we study the description of the non-linear evolution of the bispectrum. A fitting formula for the bispectrum was proposed by Scoccimarro & Frieman (1999), hereafter SF, that interpolates between perturbative and the highly non-linear regime, assuming scale-free initial conditions. Since current models for structure formation are not scale-free, these results could not be used to make realistic predictions to be compared with observations. In this work, we extend the SF results to provide a fitting formula for the non-linear evolution of the bispectrum in CDM models.

## 2 The Bispectrum

We denote by $`\stackrel{~}{\delta }`$ the Fourier transform of the density contrast $`\delta `$ obtained from the density field $`\rho =\overline{\rho }(1+\delta )`$, where $`\overline{\rho }`$ is the mean density. The power spectrum $`P(k)`$ and the bispectrum $`B_{123}=B(k_1,k_2,k_3)`$ are defined as

$`\stackrel{~}{\delta }(𝐤_1)\stackrel{~}{\delta }(𝐤_2)`$ $`=`$ $`\delta _D(𝐤_1+𝐤_2)P(k)`$ (1)
$`\stackrel{~}{\delta }(𝐤_1)\stackrel{~}{\delta }(𝐤_2)\stackrel{~}{\delta }(𝐤_3)`$ $`=`$ $`\delta _D(𝐤_1+𝐤_2+𝐤_3)B_{123}.`$ (2)
In this paper, we assume Gaussian initial conditions. In this case, the bispectrum generated by gravitational instability at large scales is given by tree-level (second-order) perturbation theory, and reads (Fry 1984)

$$B_{123}=2F_2(𝐤_1,𝐤_2)P_1P_2+\mathrm{cyc}.,$$
(3)
where $`P_iP(k_i)`$, and

$$F_2(𝐤_1,𝐤_2)=\frac{5}{7}+\frac{1}{2}\frac{𝐤_1𝐤_2}{k_1k_2}\left(\frac{k_1}{k_2}+\frac{k_2}{k_1}\right)+\frac{2}{7}\left(\frac{𝐤_1𝐤_2}{k_1k_2}\right)^2.$$
(4)
Note that Eq.(2) implies that the bispectrum is defined only for closed triangles formed by the wave vectors $`𝐤_1,𝐤_2,𝐤_3`$. It is convenient to define the reduced bispectrum $`Q_{123}`$ as

$$Q_{123}=\frac{B_{123}}{P_1P_2+P_2P_3+P_3P_1},$$
(5)
which takes away most of the dependence on scale and cosmology.

SF gave a fitting formula for the non-linear evolution of the bispectrum from scale-free initial conditions. They replaced the kernel in Eq. (4) by the effective kernel

$`F_2^{\mathrm{eff}}(𝐤_1,𝐤_2)`$ $`=`$ $`{\displaystyle \frac{5}{7}}a(n,k_1)a(n,k_2)`$ (6)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{𝐤_1𝐤_2}{k_1k_2}}\left({\displaystyle \frac{k_1}{k_2}}+{\displaystyle \frac{k_2}{k_1}}\right)b(n,k_1)b(n,k_2)`$
$`+`$ $`{\displaystyle \frac{2}{7}}\left({\displaystyle \frac{𝐤_1𝐤_2}{k_1k_2}}\right)^2c(n,k_1)c(n,k_2),`$
where the functions $`a(n,k)`$, $`b(n,k)`$, and $`c(n,k)`$ are chosen to interpolate between these two regimes according to the one-loop PT and N-body results (SCFFHM). This yields ($`2n0`$)

$`a(n,k)`$ $`=`$ $`{\displaystyle \frac{1+[0.7Q_3(n)]^{1/2}(kR_0)^{n+6}}{1+(kR_0)^{n+6}}},`$ (7)
$`b(n,k)`$ $`=`$ $`{\displaystyle \frac{1+0.2(n+3)(kR_0)^{n+3}}{1+(kR_0)^{n+3.5}}},`$ (8)
$`c(n,k)`$ $`=`$ $`{\displaystyle \frac{1+4.5/[1.5+(n+3)^4](kR_0)^{n+3}}{1+(kR_0)^{n+3.5}}},`$ (9)
where $`R_0`$ is a measure of the non-linear scale, defined from the linear variance $`\sigma _L^2(R_0)1`$ with a Gaussian window function. The function $`Q_3(n)`$ is given by

$$Q_3(n)=\frac{(42^n)}{(1+2^{n+1})}.$$
(10)
Note that at large scales, when the functions $`a=b=c=1`$, we recover the tree-level PT expression for the bispectrum; on the other hand, at small scales $`a^2=(7/10)Q_3`$ and $`b=c=0`$, so Eqs. (3) and (6) give $`Q_{123}=Q_3(n)`$ independent of triangle configuration: this is the prediction of hyperextended perturbation theory (HEPT), which was shown to describe well the bispectrum (and higher order moments) in the strongly non-linear regime (SF).

Summarizing, the non-linear bispectrum for scale-free intial conditions is given by Eq.(3) and Eqs.(6-9), together with the non-linear power spectrum obtained, for example, from the fitting formula of Peacock & Dodds (1996). In the next section we provide a generalization of the bispectrum fitting formula for CDM models.

## 3 A Fitting Formula for CDM Models

The numerical simulations used in this work correspond to cluster-normalized cold dark matter (CDM) models run by the Virgo collaboration (see e.g., Jenkins et al. 1998 for more details). There are four different models, SCDM, $`\tau `$CDM, OCDM, and $`\mathrm{\Lambda }`$CDM, see table 1 for a description of the parameters. These simulations each contain $`256^3`$ particles in a box $`240h^1`$ Mpc on a side. They were run with an adaptive P3M code (Couchman, Thomas & Pearce, 1995; Pearce & Couchman, 1997). There is only one realization for each cosmological model, and all these have the same initial phases. We have also made measurements in a set of smaller simulations run by the Hydra Consortium (Couchman, Thomas & Pearce, 1995), with very similar results. The bispectrum was measured in this set of simulations as described in SCFFHM.

Before generalizing the scale-free fitting formula (SF) reviewed in the previous section, a number of issues related to the CDM spectrum must be addressed. First, since the initial spectrum is no longer scale-free, there is no unique spectral index to be used in Eqs. (7-9). Here we adopt the local spectral index obtained from the linear power spectrum. Another possibility would be to use at scale $`k`$ the spectral index at scale $`k_L`$ obtained from the HKLM mapping, $`k=k_L[1+4\pi k^3P(k)]^{1/3}`$ (Hamilton et al 1991), which maps a non-linear scale to a linear scale with the same number of enclosed pairs. Both approaches are reasonable but can only be simple approximations to a more complex situation; in practice the transfer of power from large to small scales depends not only on the local spectral index but on the full shape of the power spectrum (e.g. Jain & Bertschinger 1994) and non-linear power at a given wavenumber is not just a function of the linear power at another scale. In the transition to the non-linear regime, one-loop calculations (SCFFHM) and numerical simulations with initial power spectra cutoff beyond some characteristic scale (Fry et al. 1993) suggest that the local index approach is a better choice. On the other hand, at small scales one would expect the HKLM choice to be clearly better, since for models such as CDM where the initial power spectrum decreases strongly with $`k`$, most of the power is generated by cascading from larger scales. Despite this, approximations for one-point higher order moments in the non-linear regime based on the local spectral index such as EPT (Colombi et al. 1997) and HEPT (SF) show very good agreement even for CDM spectra. Furthermore, attempts to use the HKLM choice gave poor results (Colombi et al. 1997).

Another issue is the dependence of the fitting formula on cosmological parameters such as $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$. Perturbative results show that the reduced bispectrum $`Q`$ is very weakly dependent on cosmology at large scales, in particular, $`F_2(𝐤_1,𝐤_2)`$ in Eq. (4) is very well approximated by its Einstein-DeSitter value (Bouchet et al. 1995). In addition, SCFFHM showed that as long as $`\mathrm{\Omega }_m/f^2(\mathrm{\Omega }_m)1`$, where $`f(\mathrm{\Omega }_m)\mathrm{\Omega }_m^{0.6}`$, all higher-order moments become independent of cosmology at all orders in perturbation theory for a fixed $`\sigma _8`$. For these reasons, we shall assume that the fitting formula does not depend on cosmology; this is certainly a good approximation within the accuracy of our measurements in simulations. Thus, all the cosmology dependence of the bispectrum $`B_{123}`$ is encoded in the cosmology dependence of the non-linear power spectrum.

We should also note that the asymptotic behavior of the bispectrum at small scales is not yet well tested against N-body simulations. Here we shall follow the standard assumption of the hierarchical ansatz, in which $`Q`$ becomes constant in the deeply non-linear regime, as implicitly assumed in Eq. (7). This behavior is consistent with our measurements in simulations in the limited range of scales that we can test.

To obtain a fitting formula for the bispectrum in CDM models we proceeded as follows. First, we used the fitting formula for the scale-free case by replacing $`kR_0`$ with $`q=k/k_{nl}`$, where $`4\pi k_{nl}^3P_L(k_{nl})=1`$, and $`P_L(k)`$ is the linear power spectrum at the corresponding redshift. This provides a reasonable, but not very accurate, fit to the simulation results, with an average deviation of $`d_{\mathrm{avg}}=40\%`$ from the N-body measurements. We then perturbed away from this fit by using two methods. We introduced six free parameters, $`a_i`$, such that $`qa_1q`$ in Eq. (7), $`n+6n+a_2`$ in Eq. (7), $`0.20.2a_3`$ in Eq. (8), $`4.54.5a_4`$ in Eq. (9), $`qa_5q`$ in Eq. (9), and finally $`[0.7Q_3(n)]^{1/2}[0.7Q_3(n)]^{1/2}\sigma _8^{a_6}(z)`$, which quantifies deviations from HEPT that depend on redshift. We then obtained the best fit parameters by doing a $`\chi ^2`$ minimization between the predicted and measured bispectrum for all triangles, using a large grid of models (containing about $`10^6`$ of them). We also looked for the best fit using the Levenberg & Marquardt Method (Press et al. 1992), the solution was consistent with the grid approach but somewhat sensitive to the initial guess for $`a_i`$. As a result of this procedure we get a best fit model with parameters $`a_1=0.25`$, $`a_2=3.5`$, $`a_3=2`$, $`a_4=1`$, $`a_5=2`$, $`a_6=0.2`$, that is

$`a(n,k)`$ $`=`$ $`{\displaystyle \frac{1+\sigma _8^{0.2}(z)\left[0.7Q_3(n)\right]^{1/2}(q/4)^{n+3.5}}{1+(q/4)^{n+3.5}}}`$ (11)
$`b(n,k)`$ $`=`$ $`{\displaystyle \frac{1+0.4(n+3)q^{n+3}}{1+q^{n+3.5}}}`$ (12)
$`c(n,k)`$ $`=`$ $`{\displaystyle \frac{1+4.5/\left[1.5+(n+3)^4\right](2q)^{n+3}}{1+(2q)^{n+3.5}}},`$ (13)
and $`q=k/k_{nl}(z)`$. The average deviation of this fitting formula from the measurements in the simulations is about $`d_{\mathrm{avg}}=15\%`$. We also found similarly good fits with $`d_{\mathrm{avg}}=16\%`$ for $`a_1=0.25`$, $`a_2=4`$, $`a_3=1.5`$, $`a_4=1`$, $`a_5=1.5`$ and $`a_6=0`$, so the small deviations from HEPT implied by Eq.(11) do not seem statistically significant.

Figures 1-5 illustrate the results. Although we have measured, and used in the fitting procedure, the bispectrum for many triangle shapes, we show for simplicity only triangles where $`k_2=2k_1`$ as a function of angle $`\theta `$ between $`𝐤_1`$ and $`𝐤_2`$, and equilateral triangles as a function of scale $`k`$. Since we only have a single realization for each cosmological model, we cannot show rigorous error bars, which in any case would only be part of the story, since different triangles in a single panel are correlated: highly so in the non-linear regime. However, we can obtain a reasonable estimate of the errors in a single realization by scaling the variance obtained from multiple realizations of $`\mathrm{\Lambda }`$CDM lower resolution runs (Scoccimarro et al 2000b): this gives an error in $`Q`$ of about $`\pm 0.25`$ for scales $`k\begin{array}{c}>\hfill \\ \hfill \end{array}0.3`$ h/Mpc approximately independent of scale.

Figures 1-4 show a comparison of the fitting formula to the N-body measurements for triangles of different scales (as denoted in each panel) as a function of triangle shape. The reduced bispectrum $`Q`$ shows the familiar configuration dependence at large scales (top-left panel), and then flattens out at small scales (bottom-right panel). This behavior is well reproduced by the fitting formula, and it is also well understood physically (SCFFHM): at large scales the bispectrum probes the shapes of large-scale structures, which leads to an anisotropic $`Q`$, whereas at small scales all configurations become equally probable due to virialization. Note that the flattening is more evident in higher $`\sigma _8`$ models, as expected, since at a given scale these models have a stronger level of non-linearity.

Comparison between corresponding panels in Fig. 3 and Fig. 4 and also between the top panels in Fig. 5 illustrate the lack of sensitivity of $`Q`$ (for a fixed level of non-linearity, or $`\sigma _8`$) to cosmological parameters such as $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, even in the non-linear regime. Note that in some cases the agreement at the smallest scales between the fitting formula and the simulations is not as good as at larger scales, especially in the $`\tau `$CDM case; however, this does not seem very serious given the estimated error of $`\pm 0.25`$ and that different triangles are strongly correlated. More numerical work would be necessary to assess this issue.

## 4 Conclusions

We provided a fitting formula for the non-linear evolution of the bispectrum in CDM models, which fits the measurements in N-body simulations to an accuracy of $`15\%`$. This will be a useful tool for making predictions about the non-linear behavior of non-Gaussian statistics related to the bispectrum of the density field, and is analagous to the role played by fitting formulas for the non-linear evolution of the power spectrum (Hamilton et al. 1991; Jain, Mo & White 1995; Peacock & Dodds 1996). Indeed, Van Waerbeke et al. (2000) apply the results obtained here to the calculation of the skewness of the convergence field in weak gravitational lensing surveys.

There are some issues in this work that merit further investigation. One obvious aspect is improving the accuracy of the fitting function; in view of the upcoming large galaxy and weak lensing surveys it would be useful to have a more accurate prediction for the non-linear evolution of the bispectrum. The main need in this regard is the availability of multiple independent numerical realizations of each cosmological model, as the bispectrum is rather sensitive to the presence of massive clusters in the simulation volume (Colombi et al. 1996). Another issue is the extension towards smaller scales, and the validity of the hierarchical ansatz. Analytical models based on dark matter halos predict that at small scales, $`k\begin{array}{c}>\hfill \\ \hfill \end{array}3`$ h/Mpc, $`Q`$ should increase with $`k`$ (Ma & Fry 2000, Scoccimarro et al. 2000b). Testing this prediction reliably will require the use of large-volume, higher-resolution simulations.

## 5 acknowledgments

R.S. is supported by endowment funds from the Institute for Advanced Study and NSF grant PHY-0070928 at IAS. H.M.P.C. Acknowledges support from NSERC. The numerical simulations analyzed in this paper were carried out by the Virgo Supercomputing Consortium (http://star-www.dur.ac.uk/ frazerp/virgo/virgo.html) using computers based at the Max Plank Institut fur Astrophysik, Garching and the Edinburgh Parallel Computing Centre.


text:
# A Model for Bone Strength and Osteoporotic Fractures

## Abstract

Inner porous regions play a critical role in the load bearing capability of large bones. We show that an extension of disordered elastic networks \[Chung et. al., Phys. Rev. B, 54, 15094 (1996)\] exhibits analogs of several known mechanical features of bone. The “stress-backbones” and histograms of stress distributions for healthy and weak networks are shown to be qualitatively different. A hereto untested relationship between bone density and bone strength is presented.

Osteoporosis is a multi-faceted metabolic disease that reduces bone strength and leads to a significant number of fractures occurring in older adults . Unfortunately, therapeutic agents available for prevention and treatment of osteoporosis often induce adverse effects in patients . Thus, non-invasive diagnostic tools to determine the need for therapeutic intervention are essential for effective management of osteoporosis. Simple models can form a useful complement to traditional studies of osteoporosis. In this Letter we introduce such a model, and describe some of its features.

Large bones such as thigh bones and vertebrae consist of an outer cylindrical shaft (cortex) and an inner porous region (trabecular architecture) . The cortex is made of compact bone and has a thickness of several millimeters. The structure of the trabecular architecture (TA) is that of a disordered cubic network of “trabeculae” whose axial and cross sectional dimensions are of the order 1mm and 0.1mm respectively  (see Figure 1).

Routine activities (e.g., climbing stairs) inflict micro-damage on bone. Material in the neighborhood of these fractures is resorbed by a class of cells known as “osteoclasts.” Their presence attracts a second group of cells, “osteoblasts,” which help regenerate lost bone . This sequence of events, referred to as bone remodeling, reduces the accumulation of micro-damage shallower than about 0.1mm . The full restoration of bone strength typically takes a period of 2-3 months.

In the cortex, deeper fractures, created during occasional trauma, are not repaired through remodeling. The resulting micro-fracture accumulation leads to lower bone quality and fracture toughness , reducing the load bearing capacity of the cortex with aging.

Since the thickness of trabecalae is $``$0.1 mm, fractures that do not sever them can be remodeled. Clinical studies also indicate that perforated trabeculae are seldom regenerated  and that cross sections of those surviving changes little . Thus, even though the connectivity of the TA reduces with aging, the surviving trabeculae can (mostly) be expected to retain their quality. Indirect evidence for this conjecture is provided by mechanical studies which have shown that a TA from a given skeletal location fractures at a fixed level of strain independent of the age of the bone . (In contrast, the breaking stresses of such samples reduce with age, though much less than the corresponding degradation of a cortex.) Due to efficacy of its remodeling, the TA becomes the principal load carrier in bones of older adults.

Reductions in estrogen and testosterones lead to increased bone turnover and an imbalance in remodeling. The accompanying weakening of a TA increases the susceptibility of bone to fracture . The management of osteoporosis can be greatly aided by the availability of characteristics that can identify this weakening.

Bone Mineral Density (BMD), or the effective bone density is the principal clinical measurement used as a surrogate of bone strength . It is estimated using x-ray imaging via an evaluation of the opacity of bone. Large clinical studies have shown an exponential increase in fracture incidence with the loss of BMD . Among other features known to be relevant for bone strength is the architecture of the porous bone. Structural properties of a TA such as the average width of the trabeculae, mean trabecular spacing, connectivity of the network and its fractal dimension have been proposed as additional surrogates of the ultimate (or breaking) strength of bone . Alternative approaches to estimate bone strength involve the use of finite element computations  and cellular models . All previous analyses have been carried out on TAs from bone samples. This makes it nearly impossible to isolate effects of individual factors on bone decay.

We introduce a simple mechanical model of a TA to complement these analyses. Since model parameters can be varied independently, we expect the task of elucidating the essential differences between healthy and osteoporotic bone to be simplified. Possible diagnostic tools for osteoporosis can then be identified using their ability to quantify these differences.

Images (e.g., Figure 1) suggest that disordered elastic networks  may be used to model TAs. Our initial studies are conducted on a square network of linear springs. An unstressed configuration is generated by displacing the vertices of the square grid (of side $`D`$) randomly by an amount less than $`D\mathrm{\Delta }`$. Externally imposed deviations increase the potential energy of the network via a combination of elastic ($`\frac{1}{2}k(\mathrm{\Delta }l)^2`$) and bond bending ($`\frac{1}{2}\kappa (\mathrm{\Delta }\theta )^2`$) contributions. Here $`\mathrm{\Delta }l`$ and $`\mathrm{\Delta }\theta `$ are changes in the length of a spring and the bond angle between adjacent springs respectively . The differences between trabeculae are modeled by assigning random values for $`k`$ and $`\kappa `$; specifically $`k[k_0(1\eta _e),k_0(1+\eta _e)]`$ and $`\kappa [\kappa _0(1\eta _b),\kappa _0(1+\eta _b)]`$, where $`\eta _e`$ and $`\eta _b`$ are predetermined. Since the network is constrained to lie on a plane, torsional effects are ignored. The following conditions are included to model known features of bone.

* As discussed above, the fracture criterion for trabeculae should be based on the level of strain . Following mechanical studies of fracture, we assume that fracture strains $`\gamma `$ are distributed on a (two-parameter) Weibull distribution , with cumulative probability

$$C(\gamma )=1exp\left[(\gamma /\gamma _e)^m\right].$$
(1)
An elastic element which is strained beyond its fracture strain is removed from the network along with all bond-bending forces it contributes to. The bonds are assumed to fracture when changed beyond an amount distributed on a second Weibull distribution (with parameters $`\gamma _b`$ and $`m`$); the bond and its side with the smaller spring constant are removed from the network.
* Osteoporosis is modeled by a random removal of springs from the network. The probability of removal, $`\nu `$, is used to quantify the “level of osteoporosis;” links eliminated at a given value of $`\nu `$ are not re-introduced later. The surviving springs are assumed to retain their strength.
* Clinical studies show a dramatic increase of bone strength (disproportionate to the increase of BMD) following therapeutic intervention . It is explained by assuming that bone remodeling is preferentially carried out in locations of high stress (Wolff’s Law ). In the model, elasticity of the spring under maximum strain is increased by a small fraction $`p`$; those strained by a fraction $`f>f_0`$ of the maximum are strengthened by $`p\left(\frac{ff_0}{1f_0}\right)^\beta `$. This algorithm is implemented over repeated application of a given external strain.

Below, we discuss properties of networks subjected to uniform strain $`\zeta `$, reached through a sequence of small, adiabatic increments . When elastic elements are removed from the network due to fracture, equilibrium is recalculated prior to increasing $`\zeta `$. The sides of the network are constrained to move vertically, as is the case for a TA because it is connected to the cortex. The force $`T(\zeta )`$ required to sustain a given strain is estimated by adding the vertical forces on the upper surface. Results analogous to those presented are observed under four-point bending which is a leading mechanism of fracture of long bones .

The stress-strain relationship $`T(\zeta )`$, which is initially linear becomes nonlinear beyond the “yield point” and reaches a maximum (ultimate stress) $`T_{max}`$ prior to failure of the network. The yield point coincides with the first fracture of springs. Indirect evidence has been presented to suggest that yielding of a bone coincides with cracking of trabeculae . Failure of a network is accompanied by a crack propagating across the entire network.

Figure 2(a) shows the stress distributions on a “healthy” network strained below yield. Springs experiencing large stresses (forming the “stress backbone”) are distributed evenly throughout the network. The histogram of stresses for this configuration contains a broad peak, similar to histograms for elastic networks . When $`\zeta `$ is increased, the stress backbone occupies a smaller subset following fracture of a group of springs. Similar changes have been observed in elastic networks  and in finite element computations on digitized images of bone .

The nature of the stress backbone is very different in weaker networks; i.e., those with larger values of $`\nu `$. As shown in Figure 2(b), the stress backbone for such networks consists of a few coherent pathways. Beyond $`\nu >\nu _020\%`$ there is no peak in the stress-histogram. We expect these differences to prove useful in identifying new diagnostic tools for osteoporosis.

In Figure 2, the X’s denote long horizontal fractures; specifically, locations where four or more consecutive vertical bonds are absent. They prevent the propagation of stress in significant vertical slices, even when bond-bending forces are included. When $`\nu `$ is large, the stress backbone can be seen to avoid these regions, and consequently contains only a small fraction of available bonds.

The last observation can be used to estimate the decay of bone strength with increasing $`\nu `$. Since an externally applied stress passes through every horizontal layer, we need to consider stress propagation on each one dimensional slice. This “chain-of-bundles” model  further assumes that the stress on each edge of a fracture containing $`k`$ sites is enhanced by a factor $`(1+\frac{k}{2})`$; i.e., the edges share the load assigned to bonds in the fracture.

Consider an elastic network of size $`N\times N`$. The probability of forming a fracture of length $`k`$ (with links on either side) is $`\nu ^k(1\nu )^2`$. Hence the size of the largest fracture $`k_m`$ (which occurs with a probability of $`1`$ ) can be estimated by

$$N^2\nu ^{k_m}(1\nu )^21.$$
(2)
For sufficiently large fractures $`k_m2\mathrm{ln}N/\mathrm{ln}\nu `$. In the chain-of-bundles model, the stress on each side of this fracture is $`(1+\frac{k_m}{2})\frac{T}{N}`$, where $`T`$ is the externally applied stress. When this stress is the typical breaking stress of an elastic element, the fracture will propagate. Thus $`T_{max}`$ is related to $`\nu `$ through

$$T_{max}\frac{1}{k_m}\mathrm{ln}\nu .$$
(3)
Figure 3 shows that the decay of $`T_{max}`$ in model networks is consistent with this expression.

A similar analysis of “penny-shaped” fractures in three dimensional networks gives

$$T_{max}\sqrt{\mathrm{ln}\nu }.$$
(4)
Results from mechanical studies of bone samples (with an assumed power law relationship) show that $`T_{max}(BMD)^\alpha `$, where $`\alpha 2.6`$ . There is no theoretical basis for this power law, and Eqn. (4) also provides a satisfactory fit to the data. More sensitive experiments need to be conducted to discriminate between these possibilities.

We finally discuss “therapeutic regeneration” (introduced above) of a network. There has not been quantitative studies on this issue, but it has been suggested that therapeutic regeneration will be more effective in increasing the strength of healthy bone . This suggestion is based on an assumed absence of hysteresis with changes in BMD. The validity of this proposal can be tested in the model system . Figure 4 shows the fractional increases of $`T_{max}`$ experienced by two networks that are formed by degrading the same lattice by $`\nu =10\%`$ and $`\nu =20\%`$. As seen in clinical studies, there is a dramatic enhancement of bone strength ($``$ factor 2), disproportionate to the increase of BMD ($`2\%`$). For the example shown, the weaker network is seen to strengthen by a larger factor under increase of BMD. There are large fluctuations of the enhancement between distinct networks, and in occasional examples the growth of the stronger network is larger.

In this Letter, we have argued that an extension of elastic networks can be used to model the inner porous regions of bone. Since the latter is the principal load carrier in bones of older adults, the model can prove useful in identifying characteristics of osteoporotic bone. Eqn. (4) is compatible with published data, but more sensitive tests are needed to discriminate it from a possible power law decay. Examples from the model provide evidence that bone regeneration under Wolff’s law is not compatible with the absence of hysteresis in bone strength with bone density. In fact, generally, regeneration is more efficient in strengthening networks that are weaker.

We propose to use the model identify new diagnostic tools for osteoporosis, using characteristics that can differentiate between stress backbones (Figure 2). Preliminary studies indicate that the ratio of responses of a network to static and periodic strain is a suitable surrogate of bone strength .

The authors would like to thank M. P. Marder for discussions on fracture and S. R. Nagel for suggesting the possible use of response measurements to quantify variations in stress backbones. We would also like to thank the referees of earlier drafts of the manuscript for their comments. This research is partially funded by the Office of Naval Research (GHG), National Science Foundation (GHG and KEB), Texas Higher Education Coordinate Board (GHG, KEB, KKM and SJW), Department of Energy (KKM), and National Institutes of Health (SJW).


text:
# High Speed Interconnects and Parallel Software Libraries: Enabling Technologies for the NVO11footnote 1This work is sponsored by DARPA, under Air Force Contract F19628-00-C-0002. Opinions, interpretations, conclusions and recommendations are those of the author and are not necessarily endorsed by the United States Air Force.

## 1. Introduction

The National Virtual Observatory (NVO) will directly or indirectly touch upon all steps in the process of transforming raw observational data into “meaningful” results. These steps include:

* Acquisition and storage of raw data.
* Data reduction (i.e. translating raw data into source detections).
* Acquisition and storage of detected sources.
* Multi-sensor/multi-temporal data mining of the products of steps (1), (2) and (3).

The highly distributed nature of the NVO places new twists on all of these steps. Future NVO research is likely to focus on developing the software tools necessary for Step (4) as well as the methods for “federating” data from Steps (1) and (3). However, past experience with individual surveys indicates that Step (2) has dominated computer software and hardware costs and may have a large impact on the NVO.

Federation of data sets from multiple institutions, which is a primary NVO goal, will be made significantly easier if improvement of the data reduction pipeline software is also undertaken. Addressing the challenges of Step (2) for the NVO can be accomplished by significantly improving the software environment for data reduction pipelines. Although the NVO can and should influence this effort it may be outside of the NVOs core activities.

The rest of this paper presents a further analysis of the computing and networking requirements of the NVO and provides a discussion of some of the challenges and solutions for addressing data reduction for these massive NVO data sets.

## 2. Large Survey Requirements

Large area imaging surveys are generating data at an exponentially increasing rate. Reducing this data is a significant hardware and software challenge that is likely to push the limits of computing for some time to come. The computing requirements of these surveys are best understood in terms of the number of operations required to process a single pixel or a single detection. By looking at the data reduction pipeline in these terms it is easy to determine the overall computing needs of a given data stream by multiplying the data rate by the per pixel processing requirement.

### 2.1. Computing Requirements

Figure 1 shows the steps that are common to most data reduction pipelines. The primary driver is the matched filtering step which involves 2D convolutions of the entire image. Such convolutions can easily result in as many as $`10^4`$ operations per pixel and requires a sizable amount of computing to keep up with the data rate. Detection processing also involves 2D convolutions but of a sparser nature. Although the real-time processing requirements of a survey may be readily satisfied by a few tens of computers, it is quite common to re-run the data processing pipeline on several years worth of data after the pipeline software has been upgraded. In such cases, hundreds or even thousands of computers will be needed to re-process the data in a reasonable time.

For example, to support a camera with a real time data rate of 100 million pixels per second (e.g. the LSST) requires a computer system that can deliver 100 Gigaflops, which is approximately 1000 of todays state of the art workstations. To re-process such a data stream could require millions of such computers.

### 2.2. Software Requirements

The primary difficulty in developing high performance pipeline software is the large number of systems architecture issues (number of processors, total memory, network bandwidth, disk bandwidth, …) that need to be considered in order to keep up with the data rate (see Figure 2). In other words, the software pipeline becomes highly tuned to the system, which in turn increases the size of the code and the expertise necessary to maintain the code. In addition, upgrading to a new system can require a significant re-write. All of these issues mean that there is very little code re-use across different data pipelines.

Unfortunately, the need for highly tuned software will only grow in the future as next generation computers incorporate more complex features (e.g., multi-processor chips, on chip vector units, multi-threading, …). While it is now possible to get 1% or 10% of a computers peak performance with little or no optimization, in the future an un-optimized program can expected to give less than 0.1% of peak performance. In other words, the price of not being optimized to the hardware will increase.

## 3. Enabling Technologies

The primary goals of a computer system architect are to manage complexity and to create systems that avoid the need for heroic programming efforts to meet the program milestones. In the context of a real-time data reduction pipeline this means developing simplified abstractions of the hardware and software so as to isolate the specifics of the application from the specifics of the hardware.

In the previous section, the argument for complex, highly parallel computing systems was presented. Without such systems it will require months to analyze a single nights worth of observations, and significantly longer to re-process a significant portion of a survey.

### 3.1. High Speed Interconnects

One way to create a less complex parallel computing systems is to invest in extremely capable networks. Minimizing network usage is the single most common optimization step performed in writing parallel software. High performance networks greatly alleviate the pressure on the programmer to implement to the specific hardware. In other words, it is incumbent the system architect to assemble a “balanced” system (i.e. computation vs. communication). Fortunately, the rapid rise of the Internet and cluster computing has driven the need for ever more capable interconnects. As a result a variety of technologies will be available for producing systems that are well balanced for data processing pipelines (see Figure 3).

### 3.2. High Performance Parallel Software Libraries

A fast method for implementing high performance data pipelines is to re-use already optimized code. The best way is to leverage existing libraries (e.g. Lapack, ScaLapack, FFTW, VSIPL, …) developed by other communities (see Figure 4). These software packages remove the majority of the effort required to achieve optimal performance on a given computer. In addition, it is important for the community to increase the capability of the data pipeline applications it has developed (e.g. IRAF, IDL, …). Currently, these tools provide a variety of application specific functions. Unfortunately, they are not designed for real-time parallel data pipelines. It would be highly beneficial to upgrade (and add to) these tools so that they can exploit the hardware technology that will be required to effectively process large surveys.

## 4. Summary

The NVOs core data mining and archive federation activities are heavily dependent on the underlying data pipeline software necessary to translate the raw data into scientifically relevant source detections. The data pipeline software dictates: the raw data storage and retrieval mechanisms, the meaning and format of the fields in the source catalogs, and the ability of the NVO users to re-analyze raw data for their own purposes. Increasing the performance of the core data pipeline software so that it can address the needs of current and future high data rate surveys is an important activity that should be addressed in concert with the development of the NVO.