text:
# The Quantum-Classical Transition in Nonlinear Dynamical Systems

\[

## Abstract

Viewed as approximations to quantum mechanics, classical evolutions can violate the positive-semidefiniteness of the density matrix. The nature of this violation suggests a classification of dynamical systems based on classical-quantum correspondence; we show that this can be used to identify when environmental interaction (decoherence) will be unsuccessful in inducing the quantum-classical transition. In particular, the late-time Wigner function can become positive without any corresponding approach to classical dynamics. In the light of these results, we emphasize key issues relevant for experiments studying the quantum-classical transition.

\]

In recent years much effort has been expended, both theoretically and experimentally, to explore the transition from quantum to classical behavior in a controlled way. In this context, the interaction of trapped, cold atoms with optical potentials, both time-dependent and independent, has become a topic of considerable interest and activity . The experimental ability to systematically study dissipative quantum dynamics of nonlinear systems is an exciting new area where the frontier between classical and quantum mechanics may be carefully examined.

In this Letter we wish to explore some of the key qualitative features of the quantum-classical transition. We establish that with $`\mathrm{}`$ fixed at a finite value, and classical dynamical evolution equations for phase space distribution functions viewed as approximations to the underlying quantum equations, classical Liouville and Master equations violate the quantum constraint of positive-semidefiniteness of the density matrix: We refer to this property of the density matrix as ‘rho-positivity’. There is thus a global obstruction to the classical limit arising directly from quantum evolutions. We argue below that rho-positivity violation can (1) serve as a guide in classifying dynamical systems with regard to classical-quantum correspondence: weak violation as Type I, strong violation as Type II, and (2) explain robustness to decoherence in the sense of avoidance of the classical limit as exemplified by dynamical localization in the (open-system) quantum delta kicked rotor (QDKR). Our results impact directly on the interpretation and design of experiments to test various aspects of the quantum-classical transition.

As described in more detail below, certain nonclassical aspects of the dynamics of the QDKR turn out to be stable to decohering effects of external noise and decoherence due to spontaneous emission. A dynamical description in terms of the Wigner function leads to two alternatives to explain this stability: (1) diffusion in the quantum Master equation is simply not efficient at suppressing quantum interference terms present in the Wigner function, or (2) the much more intriguing possibility that the diffusion is successful at decohering the Wigner function, i.e., interference terms are suppressed and the Wigner function is (almost) everywhere positive, yet the late-time distribution is not the solution of the corresponding classical Fokker-Planck equation. We find that the second possibility is the one actually realized, and show how it arises as a consequence of the fact that the classical Fokker-Planck equation violates rho-positivity, while the quantum Master equation does not.

The singular nature of the classical limit $`\mathrm{}0`$ in quantum mechanics has been appreciated for a long time. However, what has not been stressed sufficiently is the reason for this singular behavior, that classical dynamics violates unitarity and rho-positivity, and thus $`\mathrm{}=0`$ cannot be connected smoothly to $`\mathrm{}0`$. A simple example suffices to make this point clear. Let us consider as initial condition a pure Gaussian state. Let us suppose that we evolve the corresponding (positive) Wigner function classically in some nonlinear potential (for linear systems classical and quantum dynamics are identical ), then the distribution becomes no longer Gaussian, but is still positive-definite. Three possibilities now present themselves: the evolved object can be interpreted as (1) a pure quantum state (unitarity is preserved), (2) a mixed quantum state (rho-positivity is preserved), and (3) cannot be interpreted as a quantum state (rho-positivity is violated). The first possibility can be dismissed using Hudson’s theorem: the only pure state with positive Wigner function is a Gaussian state with a (necessarily) Gaussian Wigner function . But our distribution is non-Gaussian. As to the second, we first note that the phase space integral of any function of the phase space distribution is preserved under a Liouville flow. In particular $`f^2(x,p)𝑑x𝑑p`$ remains constant. For Wigner functions this quantity is proportional to $`Tr\rho ^2`$ which is a direct measure of whether a state is mixed or not – since this measure cannot change, the evolved object is not interpretable as a mixed state. Thus we are forced to the third alternative, that the evolved object cannot be interpreted as a quantum state at all, i.e., the Weyl transform of the evolved classical distribution yields a ‘classical density matrix’ which is non-rho-positive .

The above analysis makes it clear that the classical Liouville equation can never arise as a formal limit of quantum theory. However, all real experiments deal with open systems, i.e., systems interacting with their environment, of which the particular case of a measuring apparatus (necessary to deduce classical behavior) is an important example. Quantum decoherence and conditioned evolution arising as a consequence of such system-environment couplings and the act of measurement and observation provide a natural pathway to the classical limit as has been quantitatively demonstrated, e.g., in Ref. . Thus, it is important to inquire into the role of rho-positivity violation in this context: When is it important, and when not?

Conditions have been previously derived that apply to the extraction of (noisy) classical trajectories via continuous quantum measurement . Once these (strict) conditions are satisfied (typically in the macroscopic limit $`\mathrm{}S`$, where $`S`$ is the system action), measurement induces classical behavior, and in this regime all systems are therefore Type I. However, when these conditions are not satisfied, as is the case in most current experiments, the differences are indeed important. The key point is that, under continuous measurement, Type I systems can violate the classicality conditions in the sense that individual classical trajectories cannot be extracted, yet expectation values are close to the classical results, whereas in Type II systems, violation of the classicality condition also implies a violation of correspondence at the level of expectation values. We will demonstrate this for the QDKR below.

A quantum Master equation representing a weakly coupled, high temperature environment often utilized in studies of decoherence is

$`{\displaystyle \frac{}{t}}f_W`$ $`=`$ $`L_{cl}f_W+L_qf_W+D{\displaystyle \frac{^2}{p^2}}f_W;`$ (1)
$`L_{cl}`$ $``$ $`{\displaystyle \frac{p}{m}}{\displaystyle \frac{}{x}}+{\displaystyle \frac{V}{x}}{\displaystyle \frac{}{p}},`$ (2)
$`L_q`$ $``$ $`{\displaystyle \underset{\lambda odd}{}}{\displaystyle \frac{1}{\lambda !}}\left({\displaystyle \frac{\mathrm{}}{2i}}\right)^{\lambda 1}{\displaystyle \frac{^\lambda V(x)}{x^\lambda }}{\displaystyle \frac{^\lambda }{p^\lambda }}.`$ (3)
An unraveling of the weakly-coupled, high temperature environment, Master equation (1) is provided by a continuous measurement of position. This process is described by a stochastic Master equation for the density matrix $`\rho (t)`$, conditioned on the measurement record $`X+\xi (t)`$ with $`\xi (t)(8\eta k)^{1/2}dW/dt`$,

$`\rho (t+dt)`$ $`=`$ $`\rho ({\displaystyle \frac{i}{\mathrm{}}}[H,\rho ]k[X,[X,\rho ]])dt`$ (5)
$`+\sqrt{2\eta k}([X,\rho ]_+2\rho Tr\rho X)dW,`$
where $`k`$ is a constant specifying the strength of the measurement, $`\eta `$ is the measurement efficiency and is a number between 0 and 1, and $`dW`$ is a Wiener process, satisfying $`(dW)^2=dt`$. When $`\eta =1`$, the evolution preserves the purity of the state and can be rewritten in a way which allows it to be understood as a series of diffuse projection measurements . Averages over the resulting Schrödinger trajectories reproduce expectation values computed using the reduced density matrix $`\rho `$ or the corresponding Wigner function $`f_W(x,p)`$ obtained from solving the Master equation (1). The strength of the measurement is related to the diffusion coefficient of Eq. (1) by $`D=\mathrm{}^2k`$.

When the diffusion constant $`D=0`$, Eq. (1) is just the quantum Liouville equation for the closed system. Note that the linearity of the quantum Liouville equation implies that in order for the evolution to be unitary, $`L_q`$ cannot be unitary since $`L_{cl}`$ is not (the sum $`L_{cl}+L_q`$ is unitary but not the operators separately). The familiar heuristic argument for obtaining the classical limit from the quantum Master equation is that the diffusion term smooths out the interference effects generated by $`L_q`$ in such a way that quantum corrections to the classical dynamics are much reduced. It has also been argued, that at finite $`\mathrm{}`$, the appropriate limiting case of the quantum Master equation is in fact the classical Fokker-Planck equation \[setting $`L_q=0`$ in Eqn. (1)\] rather than the classical Liouville equation . In any case, one immediately appreciates that if either of the classical equations are strongly rho-positivity-violating, i.e., are Type II then this implies the existence of compensatory ‘large’ quantum corrections in the quantum Master equation, and hence the above heuristic argument must fail: $`L_q`$ is responsible for more than just the generation of interference fringes in the Wigner evolution.

Previous work has already suggested the possibility that closed dynamical systems may be roughly divided into two types depending on the (dynamical) classical-quantum correspondence as follows: (1) Type I systems in which quantum expectation values and classical averages track each other relatively closely as a function of time , e.g., the driven Duffing oscillator with Hamiltonian,

$$H_{\text{duff}}=p^2/2m+Bx^4Ax^2+\mathrm{\Lambda }x\mathrm{cos}(\omega t),$$
(6)
and (2) Type II systems in which the quantum and classical averages diverge sharply after some finite time, e.g., dynamical localization in the QDKR . The QDKR Hamiltonian is

$$H_{\text{dkr}}=\frac{1}{2}p^2+\kappa \mathrm{cos}q\underset{n}{}\delta (tn).$$
(7)
We solved the classical and quantum Master equations corresponding to Eqs. (6) and (7) using a high-resolution spectral solver implemented on parallel supercomputers. The solver explicitly respects rho-positivity conservation.

We verified that in both the QDKR and the Duffing oscillator numerical examples discussed below the localization condition  necessary to obtain classical trajectories was violated. The relevant condition is $`8\eta k(_x^2F/F)\sqrt{_xF/2m}`$ where the force $`F`$ and its derivatives are evaluated at typical points in phase space. For both cases we have in fact, $`8\eta k(_x^2F/F)\sqrt{_xF/2m}`$ thus localization does not occur (direct numerical solution of the corresponding stochastic Schrödinger equation confirms this result) and, as discussed earlier, a meaningful distinction between Types I and II is possible. As the value of $`\mathrm{}`$ is reduced (with $`D`$ fixed and non-zero) one does expect an approach to the classical limit , though the trajectory in the space of $`D`$ and $`\mathrm{}`$ need not be simple .

Our numerical code returns us the classical distribution function, the quantum density matrix and the Wigner function as functions of time. We then numerically solve for the eigenvalues of the quantum density matrix and the eigenvalues of the Weyl transform of the classical phase space distribution (the ‘classical density matrix’). Results of one such computation are displayed in Fig. 1 for the QDKR (Type II) and Duffing system (Type I). For the QDKR, initial conditions are pure Gaussian Wigner functions characterized by the position width $`\mathrm{\Delta }x=2.5`$, momentum width $`\mathrm{\Delta }p=1`$, centered on the point $`(x,p)=(0,0)`$, and with $`\mathrm{}=5`$ and $`\kappa =10`$. The diffusion coefficient $`D=0.1`$, corresponding to $`k=0.004`$. The horizontal axis refers to the index $`i`$ corresponding to the eigenvalues $`\lambda _i`$, which are themselves plotted on the vertical axis. The solid line is a result from a numerical solution of the quantum Master equation, as expected all eigenvalues are positive (the pure initial state has one eigenvalue equaling unity, the rest being zero). The dashed line is the corresponding result from the classical Fokker-Planck equation, which is characterized by a strong contribution from negative eigenvalues. It is thus clear that the true quantum density matrix and that provided by the classical approximation are in fact quite different. In contrast, results from classical Duffing calculations show a very small contribution from negative eigenvalues \[Parameter values in the particular case shown in Fig. 1 were $`m=1`$, $`A=10`$, $`B=0.5`$, $`\mathrm{\Lambda }=10`$, $`\omega =6.07`$, $`\mathrm{\Delta }x=0.05`$, $`\mathrm{\Delta }p=1`$, $`(x,p)=(3,8)`$, $`\mathrm{}=0.1`$, $`D=0.02`$.\] These results show how rho-positivity violation may be used to distinguish the two types of dynamical systems. An important point to emphasize is that it is sufficient to only carry out the classical dynamical calculation in order to classify a dynamical system as being Type I or II. (The initial condition must of course be a Wigner function.) Also, it should be clear that non-violation of rho-positivity is a necessary but not sufficient condition for quantum-classical correspondence in terms of agreement of expectation values.

It is well-known that dynamical localization in the QDKR can be destroyed (in the sense that $`p^2(t)`$ no longer saturates at late times) by coupling to external noise or to dissipative channels (e.g., spontaneous emission) . However, what is important to note is that even in the presence of quite strong coupling to these decohering channels, the evolution does not go over to the classical one, and in this sense the QDKR is quite different from the Duffing system investigated in Ref. . In addition to numerically solving the Master equation (1) we have investigated in detail the effects of including amplitude and phase noise, timing jitter in the kicked system, and carried out more realistic simulations taking into account the effects of spontaneous emission. Stability to decoherence in the sense above was manifest in all of these cases. Since we have established that the QDKR is a Type II system (Fig. 1), this behavior is essentially forced: as long as the classical evolution strongly violates rho-positivity, it is impossible for the full evolution to ever become close to the classical limit as the quantum corrections must always be concomitantly large. The question remains whether the resulting Wigner function at least has a classical interpretation. In order to investigate this we computed as a function of time, the quantity $`\mathrm{\Gamma }=𝑑x𝑑p(|f_W|f_W)`$, which provides a global measure of negativity of the Wigner function. With $`D=0`$, one sees that $`\mathrm{\Gamma }`$ increases monotonically as the Wigner function develops the expected oscillatory structure as a consequence of quantum interference in phase space. When $`D0`$, diffusion in phase space wipes out the interference and produces an essentially positive distribution which one may interpret classically. However, because rho-positivity must be maintained, classical evolution cannot connect two such positive distributions. Thus, in Type II systems decoherence can be successful in rendering the Wigner distribution positive, but yet not lead to a classical limit. We note that in NMR systems there is an interesting question regarding when classical evolution of spins can reproduce quantum evolutions connecting spin states that have no entanglement and thus may be interpreted classically . We have shown that a similar situation can arise even in single-particle evolution where entanglement is not an issue.

Recent experiments have attempted to directly address the issue of environment-induced decoherence in the QDKR, in the context of cold atom optics . Despite some complications stemming from non-ideal realizations, the results indicate that classical and quantum evolutions agree only at inordinately large noise levels. Note that in these experiments, parametric noise or spontaneous emission was used as the decohering mechanism. (The non-selective Master equation for atomic motion in far-detuned laser light has a similar form to that of a particle subjected to continuous position measurement. However, arguments can be made that only the weak decoherence regime can be accessed in this manner.) The parameter values in our numerical work are close to those actually used in the experiments thus, as with our simulations, the experiments are not carried out in a classical regime in the sense of Ref. . And since we have demonstrated the strongly Type II nature of the QDKR, it follows immediately that to observe true classical behavior, either the current experiments have to switch to a Type I system or employ lower values of $`\mathrm{}`$. Simply increasing the measurement constant $`k`$, or equivalently $`D`$, while it produces localization, adds noise into the trajectory which must be kept small in order to achieve the classical limit. This final condition requires a reduction in $`\mathrm{}`$ as $`k`$ is increased .

The authors acknowledge helpful discussions with Tanmoy Bhattacharya, Doron Cohen, and Andrew Doherty. The work of BS was supported by the National Science Foundation and a grant from the City University of New York PSC-CUNY Research Award Program. Large-scale numerical simulations were carried out at the Advanced Computing Laboratory (ACL), LANL and at the National Energy Research Scientific Computing Center (NERSC), LBNL.


format: text

text:
# Ginsburg-Landau Expansion in a non-Fermi Superconductor

## 1 Introduction

Since the discovery of high temperature superconductivity (HTcS) the microscopic description of the HTcS was not a simple problem. It is generally accepted that the normal and superconducting phase are not describe by a Fermi theory. There are a lot of microscopical models describing the normal and superconducting phase. The nearly antiferromagnetic model proposed by Millis, Monien and Pines , non-Fermi Anderson model , marginal Fermi liquid proposed by C.M. Varma et.al.  are phenomenological models used to explain the normal state physical properties of high temperature superconductors. Using the BCS-like model, the Gorkov equations have been applied to describe the superconducting state in a non-Fermi liquid, proposed by the Anderson model . The superconducting state properties have been studied by different authors . The idea of this model is the existence of a state in bidimensional systems similar with the state described by the Luttinger liquid for one-dimensional systems. In such a system the spectral function $`A(𝐤,\omega )=ImG(𝐤,\omega )`$ satisfies the homogeneity relation $`A(\mathrm{\Lambda }𝐤,\mathrm{\Lambda }\omega )=\mathrm{\Lambda }^{1+\alpha }A(𝐤,\omega )`$ with an exponent $`\alpha `$ greater than zero . This spectral anomaly implies the break down of the Fermi liquid theory, the limit $`\alpha =0`$ is a special choice leading to the Fermi liquid model. Wen  has demonstrated that the exponent $`\alpha `$ is non-universal depending on coupling between the fermions.

In this paper we analyze the properties of the non-Fermi superconductor using the method developed in Ref.. We calculate the Ginsburg-Landau coefficients ( Sec. III ) and analyze the main physical properties of the non-Fermi superconductor (Sec. IV). Finally we compare the results with other theoretical models ( Sec. V ) for the cuprates superconductors.

## 2 Model

The non-Fermi behavior of the normal state of cuprates superconductors has been discussed by different authors  in order to explain the physical properties of the cuprate superconductors. In the following we consider the superconducting phase which is described by the BCS-like order parameter $`\mathrm{\Delta }_𝐤`$ which can be calculated from the Gorkov equations. In the superconducting phase the electrons are described by the normal and anomalous Green functions:

$$G(𝐤,\omega )=e^{i\pi (1\frac{\alpha }{2})sgn\left(\omega \right)}\frac{\left(\frac{g(\alpha )}{\omega _c^\alpha }\right)^1\left(\omega +\epsilon _𝐤i\delta \right)^{1\alpha }}{\left(\frac{g(\alpha )}{\omega _c^\alpha }\right)^2\left(\epsilon _𝐤^2\omega ^2i\delta \right)^{1\alpha }+\left|\mathrm{\Delta }_𝐤\right|^2},$$
(1)
$$F(𝐤,\omega )=\frac{\mathrm{\Delta }_𝐤}{\left(\frac{g(\alpha )}{\omega _c^\alpha }\right)^2\left(\epsilon _𝐤^2\omega ^2i\delta \right)^{1\alpha }+\left|\mathrm{\Delta }_𝐤\right|^2},$$
(2)
where $`\omega _c`$ is a cut-off energy, $`0<\alpha <1`$, g($`\alpha `$)=$`\pi \alpha /2\mathrm{sin}(\pi \alpha /2)`$. The equation for the transition temperature Tc can be obtained from the equation for the Cooper instability:

$$1\mathrm{\Pi }(\mathrm{𝟎},0)=0,$$
(3)
where the Cooper susceptibility $`\mathrm{\Pi }(0,0)`$ can be calculated following Ref.. As a result the equation for Tc is :

$$\frac{1}{V}=T_c\underset{n}{}\frac{d^2𝐤}{\left(2\pi \right)^2}\left[G(𝐤,i\omega _n)G(𝐤,i\omega _n)+F(𝐤,i\omega _n)F(𝐤,i\omega _n)\right],$$
(4)
where V is the attractive interaction between the electrons in the superconducting state. The critical temperature Tc was also calculated in Ref. and the final result is:

$$T_c=\omega _D\left(\frac{D\left(\alpha \right)}{C\left(\alpha \right)}\right)\mathrm{exp}\left[\frac{1}{V}\frac{1}{2\alpha A\left(\alpha \right)D\left(\alpha \right)}\left(\frac{\omega _c}{\omega _D}\right)^\alpha \right],$$
(5)
where $`A(\alpha )`$, $`C(\alpha )`$, $`D(\alpha )`$ are functions of $`\alpha `$ given by Ref.. The critical temperature Tc was calculated also in Ref.. The small difference from the result obtained in Ref. has the origin in the mathematical approximations made in the calculations. In deriving eq. (5) the sum over Matsubara frequencies is calculated first and than the integration over the momentum.

The result (5) gives in the limit $`\alpha 0`$ the well known result for Tc:

$$T_c=1.13\omega _D\mathrm{exp}\left[\frac{1}{N(0)V}\right],$$
(6)
which make us to believe that eq.(5) may be correct. The result given in Ref. also in the limit $`\alpha 0`$ gives the BCS result for the critical temperature.

## 3 Ginsburg-Landau Expansion

The difference between the superconducting and normal free energy densities is written in the usual Ginsburg-Landau expansion:

$$F_SF_N=a\left|\mathrm{\Delta }_𝐤\right|^2+\frac{b}{2}\left|\mathrm{\Delta }_𝐤\right|^4+c𝐤^2\left|\mathrm{\Delta }_𝐤\right|^2,$$
(7)
where $`a`$, $`b`$, $`c`$, are the Ginsburg-Landau coefficients. The calculation of the critical coefficients is done following Ref.. The coefficient $`a`$ is given by the relation:

$$a=\frac{1}{V}T\underset{n}{}\frac{d^2𝐤}{\left(2\pi \right)^2}\left[G(𝐤,i\omega _n)G(𝐤,i\omega _n)+F(𝐤,i\omega _n)F(𝐤,i\omega _n)\right].$$
(8)
Evaluation in eq. (8) was done using eq. (4) for the attractive potential between the fermions. The integration over $`𝐤`$ is transformed using the relation $`\frac{d^2𝐤}{\left(2\pi \right)^2}N(0)𝑑\epsilon `$ where $`N(0)`$ is the density of states. After a short calculation we get for $`a`$ the following expression:

$$a=2\pi ^{2\alpha 1}N\left(0\right)\left(\frac{g\left(\alpha \right)}{\omega _c}\right)^2\left(T^{2\alpha }T_c^{2\alpha }\right)\underset{0}{\overset{\mathrm{}}{}}\frac{dx}{\left(1+x^2\right)^{1\alpha }}\underset{n}{\overset{N}{}}\left(2n+1\right)^{2\alpha 1},$$
(9)
where $`x=\epsilon /\omega _n`$ is a new parameter and $`\omega _n=(2n+1)\pi T`$ are the fermionic Matsubara frequencies. $`N`$ is connected with the cutoff $`\omega _D`$ by the relation $`N=\omega _D/\pi T_c`$. After performing the integration over $`x`$ and the summation over $`n`$ we get the final result:

$$a=N\left(0\right)\frac{TT_c}{T_c}\left(\frac{\omega _D}{\omega _c}\right)^{2\alpha }g^2\left(\alpha \right)\frac{B(\frac{1}{2},\frac{1}{2}\alpha )}{\pi },$$
(10)
where Tc is the critical temperature given by eq.(5) and $`B(x,y)=\mathrm{\Gamma }\left(x\right)\mathrm{\Gamma }(y)/\mathrm{\Gamma }(x+y)`$, and $`\mathrm{\Gamma }\left(x\right)`$ are Euler functions. In the limit $`\alpha 0`$ we obtain for the coefficient $`a`$:

$$a_0=N\left(0\right)\frac{TT_{co}}{T_{co}},$$
(11)
which is the expression for the case of a two-dimensional isotropic superconductor. The coefficient $`b`$ can be written as :

$$b=T_c\underset{n}{}\frac{d^2𝐤}{\left(2\pi \right)^2}\left[G(𝐤,i\omega _n)G(𝐤,i\omega _n)+F(𝐤,i\omega _n)F(𝐤,i\omega _n)\right]^2.$$
(12)
The calculation for evaluating $`b`$ follows the same procedure described for the calculation of $`a`$, the only differences being the appearence of different exponents in the $`x`$ integration and summation over frequencies. After performing the calculations in eq.(12) we obtain for the $`b`$ the result:

$$b=N\left(0\right)\frac{7\zeta (3)}{8\pi ^2T_c^2}g^4\left(\alpha \right)\left(\frac{2\pi T_c}{\omega _c}\right)^{4\alpha }\frac{2B(\frac{1}{2},\frac{3}{2}2\alpha )}{\pi }\frac{2^{34\alpha }1}{7}\frac{\zeta \left(34\alpha \right)}{\zeta \left(3\right)},$$
(13)
where $`\zeta \left(x\right)`$ is Riemann function. In the limit $`\alpha 0`$ we obtain the well known result:

$$b_0=N\left(0\right)\frac{7\zeta (3)}{8\pi ^2T_{co}^2}.$$
(14)
In order to calculate $`c`$ we must perform a Taylor expansion in powers of $`𝐪`$ of the expression:

$$T_c\underset{n}{}\frac{d^2𝐤}{\left(2\pi \right)^2}\left[G(𝐤+\frac{𝐪}{2},i\omega _n)G(𝐤+\frac{𝐪}{2},i\omega _n)+F(𝐤+\frac{𝐪}{2},i\omega _n)F(𝐤+\frac{𝐪}{2},i\omega _n)\right]$$
(15)
The integration over $`𝐤`$ is now transformed using the relation $`\frac{d^2𝐤}{\left(2\pi \right)^2}N(0)𝑑\epsilon \underset{0}{\overset{2\pi }{}}\frac{d\theta }{2\pi }`$ because of the appearence of the $`\mathrm{cos}(\theta )`$ term where $`\theta `$ is the angle between $`𝐤`$ and $`𝐪`$. We have to consider the first and the second orders of Taylor expansion in order not to loose some important contributions. After performing the calculations in eq.(15) we get for $`c`$ the result:

$$c=N\left(0\right)\frac{v_F^27\zeta (3)}{32\pi ^2T_c^2}g^2\left(\alpha \right)\left(\frac{2\pi T_c}{\omega _c}\right)^{2\alpha }\frac{2^{32\alpha }1}{7}\frac{\zeta \left(32\alpha \right)}{\zeta \left(3\right)}$$

$$\times \frac{\left(3\alpha \right)\left(1\alpha \right)B(\frac{1}{2},\frac{5}{2}\alpha )\left(1\alpha ^2\right)B(\frac{3}{2},\frac{3}{2}\alpha )}{\pi }.$$
(16)
In the limit $`\alpha 0,`$ $`c`$ becomes:

$$c_0=N\left(0\right)\frac{v_F^27\zeta (3)}{32\pi ^2T_{co}^2}.$$
(17)
$`c_0`$ is the coefficient for the s-wave isotropic superconductor. From eqs.(11), (14), (17) we can see that we reobtain the results for the ordinary superconductors, for the Ginsburg-Landau coefficients in the limit $`\alpha =0`$.

## 4 Physical Characteristics of a non-Fermi Superconductor

In this section we present the behavior of the coherence length, the penetration depth of the magnetic field and the jump in the heat capacity at the transition point for a non-Fermi superconductor. The coherence length at a given temperature is given by the following expression between the Ginsburg-Landau coefficients:

$$\xi ^2\left(T\right)=\frac{c}{a}.$$
(18)
In the ordinary isotropic superconductor $`\xi _0\left(T\right)=\sqrt{c_0/a_0}=0.74\xi _0/\sqrt{1T/T_{co}}`$ where $`\xi _0=0.18v_F/T_{co}`$ . Using eqs. (10), (16) we obtain for the coherence length the expression:

$$\xi ^2\left(T\right)=\xi _0^2\left(T\right)\left(\frac{T_{co}}{T_c}\right)^2\left(\frac{1\frac{T}{T_{co}}}{1\frac{T}{T_c}}\right)\left(\frac{2\pi T_c}{\omega _D}\right)^{2\alpha }\frac{2^{32\alpha }1}{7}\frac{\zeta \left(32\alpha \right)}{\zeta \left(3\right)}$$

$$\frac{\left(3\alpha \right)\left(1\alpha \right)B(\frac{1}{2},\frac{5}{2}\alpha )\left(1\alpha ^2\right)B(\frac{3}{2},\frac{3}{2}\alpha )}{B(\frac{1}{2},\frac{1}{2}\alpha )}.$$
(19)
From eq. (19) we find the same temperature dependence of the coherence length $`\xi \left(T\right)=\xi /\sqrt{1T/T_c}`$. At a given temperature T the coherence length is a function of $`\alpha `$. The dependence of $`\xi \left(T\right)`$ as a function of $`\alpha `$ is presented in Fig.1 . The general expression for the penetration depth in terms of Ginsburg-Landau coefficients is given by:

$$\lambda ^2\left(T\right)=\frac{\overline{c}^2}{32\pi e^2}\frac{b}{ac}.$$
(20)
The general expression (20) can be simplified if we introduce the $`\lambda _0^2\left(T\right)`$ as the penetration depth in a ordinary superconductor  as:

$$\lambda _0\left(T\right)=\frac{1}{\sqrt{2}}\frac{\lambda _0}{\sqrt{1T/T_{co}}},$$
(21)
where $`\lambda _0^2=m\overline{c}^2/4\pi ne^2`$ is the penetration depth at T=0. Using eqs. (20) and (21) we have for $`\lambda ^2(T)`$:

$$\lambda ^2\left(T\right)=\lambda _0^2\left(T\right)\left(\frac{1\frac{T}{T_{co}}}{1\frac{T}{T_c}}\right)\left(\frac{2\pi T_c}{\omega _D}\right)^{2\alpha }\frac{2^{34\alpha }1}{2^{32\alpha }1}\frac{\zeta \left(34\alpha \right)}{\zeta \left(32\alpha \right)}$$

$$\times \frac{2\pi B(\frac{1}{2},\frac{3}{2}2\alpha )}{B(\frac{1}{2},\frac{1}{2}\alpha )\left[\left(3\alpha \right)\left(1\alpha \right)B(\frac{1}{2},\frac{5}{2}\alpha )\left(1\alpha ^2\right)B(\frac{3}{2},\frac{3}{2}\alpha )\right]}.$$
(22)
The temperature dependence of the penetration depth is the same as that in ordinary superconductors, but we have an $`\alpha `$ dependence of $`\lambda ^2\left(T\right)`$ as in Fig.2. Using eqs.(19) and (22) it is possible to calculate the Ginsburg-Landau parameter $`k=\lambda \left(T\right)/\xi \left(T\right)`$. The dependence of $`k`$ as function of $`\alpha `$ is the same with the $`\lambda \left(T\right)`$ dependence. We can also calculate the size of the discontinuity in the heat capacity at the transition point. The jump in the heat capacity is given by:

$$\frac{C_SC_N}{\mathrm{\Omega }}=\frac{T_c}{b}\left(\frac{a}{TT_c}\right)^2,$$
(23)
where $`C_S`$ and $`C_N`$ represent the superconducting and normal heat capacities, and $`\mathrm{\Omega }`$ is the volume of the system. If we consider the ratio between the jump in the heat capacity of a non-Fermi superconductor and the jump of a normal bidimensional superconductor we obtain:

$$\frac{\left(C_SC_N\right)_{T_c}}{\left(C_SC_N\right)_{T_{co}}}\frac{T_{co}}{T_c}=\left(\frac{T_c}{T_{co}}\right)^2\left(\frac{\sqrt{\omega _c\omega _D}}{2\pi T_c}\right)^{4\alpha }g^2\left(\alpha \right)\frac{B^2(\frac{1}{2},\frac{1}{2}\alpha )}{2\pi B(\frac{1}{2},\frac{3}{2}2\alpha )}\frac{7}{2^{34\alpha }1}\frac{\zeta \left(3\right)}{\zeta \left(34\alpha \right)}$$
(24)
The $`\alpha `$ dependence is presented in Fig.3. As we can see in Fig. 3 the jump in $`C_V`$ decreases with increasing $`\alpha `$.

## 5 Conclusions

We have calculated the critical temperature and the Ginsburg-Landau coefficients for a non-Fermi superconductor . We have found a similar behavior of the Ginsburg-Landau coefficients with those in the BCS case. The results are in agreement with previous calculations . A similar model was used by Sudbo  and Muthukumar  for the calculation of critical temperature. Using these results we calculated the coherence length, the penetration depth and the jump in the specific heat at the transition point. We have to mention that from a non-Fermi superconductor we can obtain the results for the BCS case by putting $`\alpha =0`$ in eqs.(19), (22), (24). Taking some reasonable value of the specific heat jump \[ see Fig. 3 \] we get that $`\alpha `$ should satisfy the condition $`0.2<\alpha <0.4`$.

In high temperature superconductivity the overdoped region is well described by the standard Fermi liquid theory. Our theory is applicable in the underdoped and slightly overdoped regions where the Fermi liquid theory breaks down. In these regions there is another problem which should be taken into consideration, the opening of a pseudogap in the normal state of these materials.

Finally we mention that these calculations give a qualitative idea on the superconductivity in the non-Fermi superconductor. A quantitative comparison with the experiment seems difficult because of the existence of many parameters in the theory.

## Acknowledgments

The author is grateful to Prof. M. Crisan for his comments on the physical aspects of the model.

Figure Captions

Fig. 1 The $`\alpha `$ dependence of the coherence length for T=0.1Tc and T=0.8Tc given by eq. (19)

Fig .2 The $`\alpha `$ dependence of the penetration depth for T=0.1Tc and T=0.8Tc given by eq. (22)

Fig. 3 The $`\alpha `$ dependence of the specific heat jump at the critical point given by eq. (24)


text:
# Small world effect in an epidemiological model

## I Introduction

How does the dynamics of an infectious disease depend on the *structure* of a population? A great amount of work has been done on the phenomenological description of particular epidemic situations . A classical mathematical approach to these problems deals with well mixed populations, where the subpopulations involved (typically susceptible, infected and removed) interact in proportion to their sizes. With these zero dimensional models it has been possible to study, among other epidemic features, the existence of threshold values for the spread of an infection , the asymptotic solution for the density of infected people , the effect of stochastic fluctuations on the modulation of an epidemic situation . A second classical approach describes spatially extended subpopulations. In this, the geographic spread of an epidemic can be analyzed as a reaction-diffusion process .

Real populations rarely fall into either of these categories, being neither well mixed nor lattices. Recently introduced by Watts and Strogatz , small world networks attempt to translate, into an abstract model, the complex topology of social interactions. Small worlds may play an important role in the study of the influence of the network structure upon the dynamics of many social processes, such as disease spreading, formation of public opinion, distribution of wealth, transmission of cultural traits, etc. . In relation to epidemiological models, it has been shown that small world networks present a much faster epidemic propagation than reaction-diffusion models or discrete models based on regular lattices of a social network .

In the original model of small worlds a single parameter $`p`$, running from 0 to 1, characterizes the degree of disorder of the network, respectively ranging from a regular lattice to a completely random graph. It has been shown that geometrical properties, as well as certain statistical mechanics properties, show a transition at $`p_c=0`$ in the limit of large systems, $`N\mathrm{}`$ . That is, any finite value of the disorder induces the small world behavior. In this article we show that a sharp transition in the behavior of an infection dynamics exists at a finite value of $`p`$.

## II Epidemic model

We analyze a simple model of the spread of an infectious disease. We want, mainly, to point to the role played by the network structure on the temporal dynamics of the epidemic. The disease has three stages: susceptible ($`S`$), infected ($`I`$), and refractory ($`R`$). An element of the population is described by a single dynamical variable adopting one of these three values. Susceptible elements can pass to the infected state through contagion by an infected one. Infected elements pass to the refractory state after an infection time $`\tau _I`$. Refractory elements return to the susceptible state after a recovery time $`\tau _R`$. This kind of system is usually called $`SIRS`$, for the cycle that a single element goes over. The contagion is possible only during the $`S`$ phase, and only by an $`I`$ element. During the $`R`$ phase, the elements are immune and do not infect.

The interactions between the elements of the population are described by a small world network. The links represent the contact between subjects, and infection can only proceed through them. As in Watts and Strogatz model, the small worlds we study are random networks built upon a topological ring with $`N`$ vertices and coordination number $`2K`$. Each link connecting a vertex to a neighbor in the clockwise sense is then rewired at random, with probability $`p`$, to any vertex of the system. With probability $`(1p)`$ the original link is preserved. Self-connections and multiple connections are prohibited. With this procedure, we have a regular lattice at $`p=0`$, and progressively random graphs for $`p>0`$. The long range links that appear at any $`p>0`$ trigger the small world phenomenon. At $`p=1`$ all the links have been rewired, and the result is similar to (though not exactly) a completely random network. This algorithm should be used with caution, since it can produce disconnected graphs. We have used only connected ones for our analysis.

Time proceeds by discrete steps. Each element is characterized by a time counter $`\tau _i(t)=0,1,\mathrm{},\tau _I+\tau _R\tau _0`$, describing its phase in the cycle of the disease. The epidemiological state $`\pi _i`$ of the element ($`S`$, $`I`$ or $`R`$) depends on this phase in the following way:

$$\begin{array}{cc}\pi _i(t)=S\hfill & \text{if}\tau _i(t)=0\hfill \\ \pi _i(t)=I\hfill & \text{if}\tau _i(t)(1,\tau _I)\hfill \\ \pi _i(t)=R\hfill & \text{if}\tau _i(t)(\tau _I+1,\tau _0)\hfill \end{array}$$
(1)
The state of an element in the next step depends on its current phase in the cycle, and the state of its neighbors in the network. The rules of evolution are the following:

$$\begin{array}{cc}\tau _i(t+1)=0\hfill & \text{if }\tau _i(t)=0\text{ and no infection occurs}\hfill \\ \tau _i(t+1)=1\hfill & \text{if }\tau _i(t)=0\text{ and }i\text{ becomes infected}\hfill \\ \tau _i(t+1)=\tau _i(t)+1\hfill & \text{if}\mathrm{\hspace{0.33em}1}\tau _i(t)<\tau _0\hfill \\ \tau _i(t+1)=0\hfill & \text{if}\tau _i(t)=\tau _R\hfill \end{array}$$
(2)
That is, a susceptible element stays as such, at $`\tau =0`$, until it becomes infected. Once infected, it goes (deterministically) over a cycle that lasts $`\tau _0`$ time steps. During the first $`\tau _I`$ time steps, it is infected and can potentially transmit the disease to a susceptible neighbor. During the last $`\tau _R`$ time steps of the cycle, it remains in state $`R`$, immune but not contagious. After the cycle is complete, it returns to the susceptible state.

The contagion of a susceptible element by an infected one, and the subsequent excitation of the disease cycle in the new infected, occur stochastically at a local level. Say that the element $`i`$ is susceptible, and that it has $`k_i`$ neighbors, of which $`k_{\text{inf}}`$ are infected. Then, $`i`$ will become infected with probability $`k_{\text{inf}}/k_i`$. Observe that $`i`$ will become infected with probability 1 if all its neighbors are infected. There are other reasonable choices for this mechanism. For example, if the susceptible had a probability $`q`$ of contagion with each infected neighbor, we would have a probability of infection $`[1(1q)^{k_{\text{inf}}}]`$. We have tested that both these criteria give qualitatively the same results for values of $`q0.2`$.

## III Numerical results

We have performed extensive numerical simulations of the described model. Networks with $`N=10^3`$ to $`10^6`$ vertices have been explored, with $`K=3`$ to $`10`$. A typical realization starts with the generation of the random network and the initialization of the state of the elements. An initial fraction of 0.1 infected, and the rest susceptible, was used in all the results shown here. Other initial conditions have been explored as well, and no changes have been observed in the behavior. After a transient a stationary state is achieved, and the computations are followed for several thousand time steps to perform statistical averages.

We show in Fig. 1 part of three time series displaying the fraction of infected elements in the system, $`n_{\text{inf}}(t)`$. The three curves correspond to systems with different values of the disorder parameter: $`p=0.01`$ (top), 0.2 (middle) and 0.9 (bottom). The three systems have $`N=10^4`$ and $`K=3`$, and infection cycles with $`\tau _I=4`$ and $`\tau _R=9`$. The initial state is random with $`n_{\text{inf}}(0)=0.1`$. The 400 time steps shown are representative of the stationary state. We can see clearly a transition from an endemic situation to an oscillatory one. At $`p=0.01`$ (top), where the network is nearly a regular lattice, the stationary state is a fixed point, with fluctuations. The situation corresponds to that of an endemic infection, with a low and persistent fraction of infected individuals. At high values of $`p`$—like the case with $`p=0.9`$ shown in the figure (bottom)—large amplitude, self-sustained oscillations develop. The situation is almost periodic, with a very well defined period and small fluctuations in amplitude. The period is slightly longer than $`\tau _0`$, since it includes the average time that a susceptible individual remains at state $`S`$, before being infected. Epidemiologically, the situation resembles the periodic epidemic patterns typical of large populations . A mean field model of the system, expected to resemble the behavior at $`p=1`$, can be easily shown to exhibit these oscillations. The transition between both behaviors is apparent—in this relatively small system—at the intermediate value of disorder $`p=0.2`$, shown in the middle curve. Here a low amplitude periodic pattern can be seen, appearing and disappearing again in a background of strong fluctuations. Moreover, the mean value of infection is seen to grow with $`p`$.

The formation of persistent oscillations corresponds to a spontaneous synchronization of a significant fraction of the elements in the system. Their phases $`\tau _i(t)`$ in the epidemic cycle become synchronized, and they go over the disease process together, becoming ill at the same time, and recovering at the same time. We can characterize this behavior with a synchronization parameter , defined as:

$$\sigma (t)=\left|\frac{1}{N}\underset{j=1}{\overset{N}{}}e^{i\varphi _j(t)}\right|,$$
(3)
where $`\varphi _j=2\pi (\tau _j1)/\tau _0`$ is a geometrical phase corresponding to $`\tau _j`$. We have chosen to let the states $`\tau =0`$ out of the sum in (3), and take into account only the deterministic part of the cycles.

When the system is not synchronized, the phases are widely spread in the cycle and the complex numbers $`e^{i\varphi }`$ are correspondingly spread in the unit circle. In this situation $`\sigma `$ is small. On the other hand, when a significant part of the elements are synchronized in the cycle, $`\sigma `$ is large. The synchronization would be strictly $`\sigma =1`$ if all the elements were at the same state at the same time. However, such a state would end up in $`N_{\text{inf}}=0`$ after $`\tau _0`$ time steps, and the epidemic would end since the system is closed and no spontaneous infection of susceptibles is being taken into account in the model.

In Fig. 2 we show the synchronization parameter $`\sigma `$, obtained as a time average of $`\sigma (t)`$ over 2000 time steps and a subsequent average over realizations of the system. Several curves are shown, corresponding to system sizes $`N=10^3,10^4,10^5`$ and $`10^6`$. A transition in the synchronization can be seen as $`p`$ runs from 0 to 1. The transition becomes sharper for large systems, at a value of the disorder parameter $`p_c0.4`$. It is worthwhile to remark that the transition to synchronization occurs as a function of the structure of the network, contrasting the phenomenon of synchronization as a function of the strength of the interaction, as in other systems of coupled oscillators . This behavior is observed for a wide range of values of $`\tau _I`$ and $`\tau _R`$. The amplest oscillations take place around $`\tau _I/\tau _R=1`$. They disappear when $`\tau _I`$ is significantly greater or smaller than $`\tau _R`$.

All the systems shown in figures 1 and 2 have $`K=3`$. We have explored higher values of $`K`$ as well. The picture is qualitatively the same, with a sharp transition from a quasi-fixed point to a quasi-limit cycle at a finite value of $`p`$. The critical value $`p_c`$ shifts toward lower values for growing $`K`$. This is reasonable, since higher values of $`K`$ approaches the system to a globally coupled one (at $`K=N`$ all the elements interact with every other one, even at $`p=0`$). So, the mean field behavior (the oscillations) can be expected to occur at lower values of $`p`$. In figure 3 this effect can be seen for a system with $`N=10^4`$ and growing values of $`K=3`$, 5 and 10.

## IV Discussion

Why does this transition to synchronization take place? Unfortunately we do not have analytical arguments to describe this yet. As mentioned before, a mean field model can be shown to predict the oscillations, but this can only be expected to describe the network at $`p=1`$, and it does not shed light on the nature of the transition at lower values of disorder. We can only advance some conjectures, based in our observation of the dynamical behavior of the system in real time on the computer screen.

An explanation involving $`L(p)`$, the typical distance between pairs of elements (defined as the size of the minimum path connecting two elements) is not plausible, since $`L`$ is known to behave critically with $`p_c=0`$, and we observe the transition at $`p>0`$. Complementarily to $`L`$, small world networks can be described by the degree of clusterization $`C(p)`$ . At low values of $`p`$, the networks are rather regular and highly clustered. As $`p`$ approaches 1, $`C`$ decreases. The crossover from high to low clusterization occurs at a higher value of $`p`$, compared to that observed in the decay of $`L`$. For systems with $`K=3`$, such as those mainly studied here, we have that this crossover is mostly concentrated between $`p=0.1`$ and $`p=0.5`$, precisely where the onset of oscillations occurs. Moreover, the change in the average clusterization $`C(p)`$ is accompanied by a corresponding one in the distribution of the clusterization at the element level, $`c_i(p)`$. Highly ordered networks, at low $`p`$, show not only a large value of the average clusterization $`C(p)`$ but also a small dispersion around it. On the other extreme, highly disordered networks, with $`p>0.5`$, exhibit a low average clusterization and also a low dispersion around it. There is an intermediate range of $`p`$, between 0.1 and 0.5, where the average clusterization shifts from high to low, but the distributions are wide. This indicates that, in this range, the system is a mixture of highly ordered, highly clustered regions and random, lowly clustered ones. If we consider that the clustering structure determines a partition of the whole system into smaller, interacting, subsystems, the global behavior of the system could be interpreted as a superposition of the subsystems’ dynamics. When $`p`$ is small the existence of large clusters (essentially one-dimensional in our model) inhibits the oscillatory behavior because, once the infection breaks into such a region, it remains restricted to it a long time until all the individuals have completed the cycle, orderly and deterministically. As $`p`$ increases, the number of big, ordered regions decreases. Elements within small regions with some long range links go through the infection cycle and become infected again before long. Some degree of synchronization can be seen here (Fig. 1, center, corresponds to $`p=0.2`$). When a critical value of $`p`$ is reached, and the system is essentially composed of enough small regions of low clusterization and a similar local dynamics, the synchronized, periodic global behavior establishes spontaneously.

In summary, we have observed a transition at a finite value of the disorder in a small world model. The dynamical behavior of an SIRS epidemiological model changes from an irregular, low-amplitude evolution at small $`p`$ to a spontaneous state of wide amplitude oscillations at large $`p`$. This may be related to observed patterns in real epidemics , where an effect of the social structure is observed in the dynamics of the disease.

The authors thank D. H. Zanette for fruitful discussions.


text:
# Cosmic 𝛾-Ray Bursts as a Probe of Star Formation History

## The Cosmic Star Formation Rate

The variation of the total cosmic star formation rate with redshift $`z`$ $``$ the SFR plot $``$ is conventionally determined by measurements of the H$`\alpha `$ luminosity of galaxies at $`z<0.2`$ and of UV luminosities at $`z>0.2`$ (Madau et al. 1996).

Two major complications in constructing the SFR plot at any $`z`$ are (i) the need to correct UV luminosities for dust extinction, and (ii) the need to assume a stellar IMF. However, both of these seem to have been addressed with considerable success: (i) by adopting the corrections of Calzetti (1997), which when applied to the CFRS data of Lilly et al. (1996) produce a SFR plot similar to that generated from infrared ISO observations (Flores et al. 1999), which are sensitive to the absorbed and re-radiated UV flux, and (ii) by using an IMF that flattens below 1 M e.g. that of Kroupa et al. 1993, which seems to be universal (Gilmore & Howell 1998). The SFR plot constructed with both of these assumptions, when integrated over cosmic time, reproduces the correct local stellar density of $`\mathrm{\Omega }_{}0.005`$ that we derive from a number of methods. We can therefore have considerable confidence in current determinations of the SFR plot at $`z<3`$ (Somerville et al. 2000).

At $`z>3`$, there are additional, potentially more serious, complications. The only types of high-$`z`$ galaxies whose contributions to the SFR plot have been unambiguously determined are the Lyman-break galaxies (LBGs; Steidel et al. 1999). Lyman-$`\alpha `$ (Ly$`\alpha `$) emitters may contribute too, but current indications are that their contribution is small (10-2 M yr-1 Mpc-3 compared to 10-1 M yr-1 Mpc-3 for the LBGs; Hu et al. 1998). However, these Ly$`\alpha `$-emitters and other galaxies which might exist but cannot be found by the selection techniques used to find either Ly-break or Ly$`\alpha `$ galaxies, could have their contributions systematically underestimated. Such galaxies might be lost due to $`(1+z)^4`$ surface-brightness dimming (Lanzetta et al. 2000), or the detectability of a small fraction of the members of a population might be enhanced by supernovae that happened to go off in those members causing us to miss the bulk of the sources representing that population (the supernovae are not dimmed by $`(1+z)^4`$ and can affect detectability if $`10\mathrm{log}(1+z)>\mathrm{mag}_{\mathrm{SN}}\mathrm{mag}_{\mathrm{gal}}`$). Therefore, indirect constraints are important at such high $`z`$. The constraint from $`\mathrm{\Omega }_{}`$ is weak since the contribution to the total integral of the SFR plot is small at high $`z`$ since $`\mathrm{d}t/\mathrm{d}z`$ is small there. Another constraint is that all the star formation at $`z>3`$ must produce enough metals to enrich the Ly$`\alpha `$ forest at $`z=3`$ (Cowie & Songaila 1998), but this is also of limited use here, since it is unclear what fraction of metals escape from the galaxies in which the stars form. A potentially more powerful probe is given by the $`\gamma `$-ray burst (GRB) formation rate plot, which must track the SFR plot closely if the majority of GRBs originate from the collapse of massive stars. This is what we investigate here.

## The $`\gamma `$-Ray Burst Rate

GRBs are detectable out to the farthest reaches of the observable Universe, and provide information about processes occurring at all cosmic epochs. Recent observations suggest that the long-duration GRBs and their afterglows are produced by highly relativistic jets emitted in core-collapse supernova explosions. Hence the redshift distributions of GRBs should track the cosmic star formation rate of massive stars accurately (Lamb & Reichart 2000; Blain & Natarajan 2000). At present, however, there are too few redshift measurements with which to estimate the global GRB formation rate. Nonetheless, these few measured redshifts can be used to calibrate properties of GRBs that might let them serve as standard candles. Fenimore & Ramirez-Ruiz (2000) have suggested that the spikiness of the burst time structure is correlated with luminosity, with smooth bursts being intrinsically less luminous. In principle, the measured spikiness combined with the observed flux can be used to obtain distances much like Cepheid observations give distance estimates from the pulsation period. Using a sample of 220 bright BATSE bursts for which high-resolution light curves were available, Fenimore & Ramirez-Ruiz (2000) estimated the evolution of the GRB formation rate from parameters measured solely at $`\gamma `$-ray energies.

## The SFR Plot at z $`<`$ 3.0

At low redshift, the GRB formation rate (BFR) plot tracks the SFR plot in shape (see Figure 1). This is the redshift regime where we have considerable confidence in the SFR plot (see $`\mathrm{\S }`$ 1). This concordance suggests: (i) the idea that GRBs do indeed track the SFR of the Universe. This would imply that the formation rate of massive stars that produce GRBs is proportional to the total SFR, as would follow if the IMF does not vary substantially over this $`z`$ range; (ii) that we can use GRBs as a high-$`z`$ probe. They are not affected by extinction or surface-brightness dimming. So one might therefore expect the high-$`z`$ BFR plot to track the high-$`z`$ SFR plot.

## The SFR plot at z $`>`$ 3.0

The high-$`z`$ BFR plot has a steep slope. This is different from the much flatter slope of the SFR plot as derived from LBG observations (Steidel et al. 1999). If the BFR plot does indeed track the SFR plot at high-$`z`$ (as it does at low-$`z`$), then this would suggest that much high-$`z`$ star formation is being missed in the optical surveys, even after corrections for dust extinction have been made.

Interestingly, Lanzetta et al. (2000) find evidence for a very steep SFR plot based on an analysis of high-$`z`$ galaxies in the Hubble Deep Fields in which they consider selection effects that can arise from surface-brightness dimming. An additional source of star-formation at high-$`z`$ could come from luminous SCUBA galaxies (Blain et al. 1999) which could be missing from optical surveys (e.g. Smail et al. 1999) because they are heavily extinguished by dust.

One important caveat is that the redshift range ($`z<3.0`$) in which the BFR plot tracks the SFR plot is similar to the redshift range in which GRBs were used to calculate the variability-luminosity relation (there is only one event with $`z>`$ 3.0). The application of this relation to faster variabilities and therefore higher luminosities (as is required to derive the BFR plot at $`z>3.0`$) depends on an extrapolation of this relation. If the extrapolation fails, this weakens our conclusion about the steep slope of the BFR at high-$`z`$. Also, if the fraction of massive stars that produce GRBs is a strong function of metallicity, then this could create a redshift-dependence of the normalization of the BFR plot relative to the SFR plot, which would also weaken the assertion that the SFR plot is steep at high-$`z`$.


text:
# Unzipping DNA : A hypothesis on correlation during replication *footnote **footnote *Proceedings of CMDAYS2K, held at Guru Ghasidas University, Bilaspur, India, Aug 29-31, 2K

## I Introduction

A double stranded DNA (dsDNA) needs to be opened up for its replication. Once done, the subsequent processes take place sequentially, like the binding of, say, helicases, single-strand-binding proteins etc, then the elongation of the new chain by polymerases etc.

Recently it has been proposed that the action of dnaA-type proteins could be thought of as a force acting at the initiation site (called “origin”). ( The function of DnaA is not to be confused with the “locomotive” or “sewing machine” action of the hypothetical replisome.) A simple model showed that there is a phase transition at a critical value of the force so that for forces less than the critical strength, DNA is a double-stranded object while for forces exceeding it, the chains get opened up.

In this paper we give a short description of this model, the principal results and the subsequent developments. We point out the features that seem to have led to “confusions”. We also show how this simple model can be extended to incorporate some other features.

The connection to in vivo replication and laboratory experiments are discussed in section III, where we propose that the correlation of unzipping is at the heart of the highly correlated collective process that goes by the name of replication.

## II Statistical Mechanical model

### A Model and results

The model proposed originally is that of a pair of flexible polymers bound together by a short-range interaction and pulled at one end by a force. We start with a homo-DNA, i.e. all base pairs identical. One end of the two strands are anchored or tied together and the other end is pulled by a force $`𝐠`$. The Hamiltonian is given by

$`{\displaystyle \frac{H}{k_\mathrm{B}T}}`$ $`=`$ $`{\displaystyle _0^N}𝑑s\left[{\displaystyle \frac{\epsilon _1}{2}}\left({\displaystyle \frac{𝐫_1(s)}{s}}\right)^^2+{\displaystyle \frac{\epsilon _2}{2}}\left({\displaystyle \frac{𝐫_2(s)}{s}}\right)^^2+V(𝐫(s))\right]+H_{\mathrm{force}},`$ (2)
$`V(𝐫)`$ $`=`$ $`v\delta (𝐫),`$ (3)
$`H_{\mathrm{force}}`$ $`=`$ $`𝐠(𝐫(N))=𝐠{\displaystyle _0^N}𝑑s{\displaystyle \frac{𝐫(s)}{s}},`$ (4)
where $`𝐫_i(s)`$ denotes the $`d`$-dimensional position of a monomer on chain $`i`$ (of elastic constant $`\epsilon _i`$) at a contour-length $`s`$ measured from the tied end ($`s=0`$), $`𝐫=𝐫_1𝐫_2`$ is the relative separation of the two chains at the same monomer $`s`$ and $`V(𝐫)`$ is a short-range potential simulating the interaction of the base pairs of two strands. Any realistic potential can be chosen here. But our interest is in the effect of the force on the bound double-strand case, and so for analytical simplicity we choose a delta function (or contact) potential, to be interpreted as a limit of a narrow square well. The strength of the well is chosen to be such that there is only one bound state in the problem.

As usual in statistical mechanics, there are two possible ensembles, namely a constant force ($`𝐠`$ constant) ensemble and a constant end-point separation($`𝐫(N)`$=constant) ensemble. For studies of phase transitions, it is useful to work in an ensemble of fixed intensive variable. We therefore choose the fixed force ensemble.

For the zero force case, there is a critical unbinding transition at $`v=v_c`$. The pulling force would like to align the strands in the direction of the force while the interaction would like to keep the strands together. There is a critical force below which the DNA remains in the dsDNA phase while for forces exceeding this critical value it gets unzipped.

Quantitatively, the partition function can be evaluated by a transfer matrix approach with the contour length as the transfer direction. The free energy per unit length is given by the largest eigenvalue of the transfer matrix and a phase transition takes place whenever there is a degeneracy of the largest eigenvalue.

The transfer matrix can conveniently be written in the form of a quantum Hamiltonian, albeit non-hermitian, for a particle of co-ordinate $`𝐫`$ (CM behaviour is like a free chain or a free particle)

$$H_q(𝐠)=\frac{1}{2}(𝐩+i𝐠)^2+V(𝐫),$$
(5)
in units of $`\mathrm{}(k_BT)=1`$ and mass $`=1`$, with $`𝐩`$ as momentum. For long chains $`N\mathrm{}`$ the free energy is the ground-state energy of this non-hermitian Hamiltonian. A phase transition takes place whenever the ground state is degenerate. The analysis done in Ref.  shows that if the ground state energy ( i.e. the binding energy of the dsDNA per unit length or per base pair ) is $`E_0`$, then the critical force is given by

$$g_c=2\sqrt{E_0}|vv_c|^{1/|2d|}.$$
(6)
where the $`v`$-dependences of $`E_0`$ close to $`v_c`$, for general $`d`$, is used. In absence of any force, DNA can be denatured either by changing temperature or by changing pH of the solution. In this respect the effective interaction parameter $`v`$ is a better variable than the temperature itself. The generic phase diagram is shown in Fig. 1.

The degeneracy of the eigenvalues leads to a diverging length scale that could be associated. for $`gg_c`$. as the characteristic length of the opened Y-fork (the replication fork). This length scale measured along the strand of the DNA gives a measure of the number of monomers opened by the force. Denoting this number by $`m`$, we get

$$mg_cg^{\nu _m}\mathrm{with}\nu _m=1.$$
(7)
($`\nu _m`$ was denoted by $`\nu _\tau `$ in Ref .)

### B Questions on the model and the results

It is known for centuries that a certain minimum pulling force is needed to move an object on a surface with friction, or for adhesion, and so, is the occurrence of a critical force an obvious result? An understanding of this threshold phenomena beyond phenomenological level is still in its infancy. Our results are to be seen in this context, and the model is different from the stick-slip models used in friction or adhesion problems. In fact, previous force measurements for unzipping of DNA have been interpreted in terms of equilibrium statistical mechanics, and not as a friction problem.

A classical diatomic molecule or a bond (or a particle in a potential well) under a force does not show any phase transition. At a quantum mechanical level, a bond with a stretching force is still described by a hermitian hamiltonian and therefore will not show the phase transition of the previous subsection. It is the polymeric property or connectivity that is crucial for the unzipping transition.

The result of Eq. 7 is not to be confused with the pulling of a polymer discussed in text books. In the case of a polymer pulled at the ends by a stretching force, there would be an extension of the chain and for a small force, the overall dimension of the polymer would scale linearly with the stretching force. In contrast, we see that, in the dsDNA situation, a small force will have no major effect and the law in Eq. 7 is a dependence on the deviation from the critical force.

Is continuum model a good approach? The answer to this question lies in the recognition of various length scales. The Hamiltonian of Eq. 2 can be discretized and studied with appropriate bond lengths and realistic $`V(𝐫)`$. We adopted a coarse-grained approach to focus on the phase transition behaviour. In other words, a continuum approach by itself is not a necessity to study the phase transition and other details of unzipping.

Previous studies seem to indicate that DNA is better represented by worm-like chains but again a renormalization group argument shows that in the long length limit such a Hamiltonian will generate the elastic term and therefore Eq. 2 is a representative of the universality class. In the same spirit the delta-function potential can be replaced by any realistic potential. But for the binding-unbinding ( melting) transition, it is only the integral of the potential that matters and so a square well or delta function is a valid starting point.

Is the restriction to homo-DNA any good? It is worth remembering that the genetic code was deciphered from homoDNA’s. For a controlled experiment to detect the unzipping transition, large homo-DNA’s should be studied. The model is easily generalizable (see below) to consider specific base-sequence.

### C Subsequent developments

After the original proposal of model for unzipping in Ref. , several studies have been made. The dynamics of pulling a polymer from a potential well, i.e., the dynamics of unzipping in one dimension has been studied at a mean field level by Sebastian. In addition to recovering the results of previous section, Ref  shows that the unzipping transition is also characterized by a diverging time scale. Denoting the characteristic time scale by $`𝗍`$, one finds

$$𝗍g_cg^{\nu _t}\mathrm{with}\nu _t=1.$$
(8)
A divergence of time scales indicates a long-range correlation in time in dynamics of unzipping.

DNA generally consists of inhomogeneous sequence of base pairs. Such a case in this model can be considered by taking the interaction energy to be monomer position dependent. The interaction term in Eq. 2 can be written as $`_0^N𝑑sb(s)v\delta (𝐫)`$, with $`b(s)`$ a variable depending only on monomer (base) position $`s`$ describing the specific details of the sequence. A case of random $`b(s)`$, a random interaction model (RANI model) has been considered earlier for the zero force case where the randomness has been found to be marginal. In the case of unzipping transition, the randomness does not destroy the sharpness of the transition but changes the exponent $`\nu _m`$ defined in Eq. 7 to $`\nu _m=2`$. The dynamics is yet to be studied.

The nature of the transition has also been considered. Maritan, Orlandini and Seno showed by Monte Carlo simulations and analytical methods that the noncrossing constraint on the chains can lead to a first order transition, and more interestingly to a “cold denaturation”. It has recently been argued from a direct evaluation of the partition function that the transition for the Hamiltonian of Eq. 2 for $`d=1`$ could be first order. It is known that depending on the reunion exponent, the phase transition for polymers could be first order, and in the quantum context, the order of transition (Eq. 5) is determined by the normalizability of the “critical” wavefunction. Several exact results for the unzipping transition in various dimensions have been obtained in Ref..

### D Extensions

The simple model does not take into account the self avoidance of the chains. Effects of self-avoidance has been studied in Ref. . In a laboratory experiment, one can take Eq. 2 as a model for DNA in a theta-solvent. In a good solvent, the self and mutual avoidance of each chain can be introduced by an imaginary random potential and averaging the partition function. Adding a term to Eq. 2

$`{\displaystyle \frac{H_{\mathrm{imag}}}{k_\mathrm{B}T}}`$ $`=`$ $`{\displaystyle _0^N}𝑑s\left[i𝒱(𝐫_1(s))+i𝒱(𝐫_2(s))\right],`$ (9)
with $`𝒱(𝐫)`$ as an annealed Gaussian random potential with zero mean and variance $`𝒱(𝐫)𝒱(𝐫^{})`$ = $`u\delta (𝐫𝐫^{})`$, the averaged partition function $`\overline{Z}=\overline{𝒟\mathrm{exp}(H)}`$, with overline denoting averaging over $`𝒱`$, leads to a hamiltonian of the form of Eq. 2 with additional terms, $`H_{\mathrm{int}}`$

$$H_{\mathrm{int}}=𝑑s𝑑s^{}\underset{i,j=1,2}{}u\delta (𝐫_i(s)𝐫_j(s^{})).$$
(10)
This represents excluded volume interaction for the chains. This imaginary potential of Eq. 9 would lead to a complex scalar potential in the quantum hamiltonian of Eq. 5.

We have ignored the double helical nature of DNA. In actual replication, topoisomerases act when there is super-coiling. Once the chains are opened up in a region by DnaA forming a Y-fork, the remaining chain is brought back to the native sate by the topoisomerase. This gives a justification for the absence of the helical configuration in the hamiltonian in Eq. 2. The topological linking number of the two chains can be described by the Gauss integral which is known to act like a real vector potential. It would be appropriate to take into account the linking as an independent quantity only if there is a change in the coiling in the process.

A toy model is a quantum Hamiltonian corresponding to a polymer with winding as

$$H_\mathrm{q}=\frac{1}{2}(p+𝒜)^2+V(𝐫)+i𝒱(r),$$
(11)
where $`𝒜=A_0(𝐫)+ig`$, with the real vector potential describing windings and the imaginary part, as before, the exerted force. Studies of the most general quantum Hamiltonian with both complex scalar and vector potential are needed.

## III Implication for biological replication

We explain the hypothesis on the role of unzipping in replication.

The replication process starts with a DnaA protein attaching itself, in interaction with the membrane, at the “origin” to start the Y-fork ( or “eye”). The next step is the binding of the various enzymes/proteins like single-strand binding (SSB) protein, helicase, topoisomerase, polymerases etc. Unlike the latter proteins/enzymes, the functionality of DnaA protein is not well understood. It is apparent that the whole process requires a strong correlation in space and time, whose source or origin is not obvious. In order to explain such correlations in subsequent processes, a structural unit “replisome” has been postulated which has never been isolated.

We make the hypothesis that

* The function of DnaA is actually to exert a pulling force (let us take this time $`t=0`$) with a force which is close to, but slightly below, the critical force for unzipping. This force tends to form a Y-fork or an eye-type bubble (depending on the location of the origin).
* The large length and time scales for the unzipping process, leads to a space time correlation in the initial nonequilibrium dynamics of unzipping. We postulate that this correlation controls the subsequent processes especially the dynamics of bindings of the subsequent enzymes/proteins.

The correlation in the unzipping dynamics can be characterized quantitatively by the early time behaviour of, say, $`m(t)m(t^{}0)`$ where $`m(t)`$ is the number of monomers unzipped at time $`t`$, or a more microscopic correlation $`r_i(s,t)r_j(s^{}N,t^{}0)`$, $`i,j`$ being the spatial components.

The dynamics of unzipping is described by a Langevin equation

$$\mu \frac{𝐫(s,t)}{t}=\frac{\delta H}{\delta 𝐫(s,t)}+\eta (s,t),$$
(12)
where $`t`$ represents time, $`H`$ is given by Eq. 2, $`\eta `$ is the thermal noise related to friction coefficient $`\mu `$ by $`\eta _i(s,t)\eta _j(s^{},t^{})=2\delta _{ij}\mu k_BT\delta (ss^{})\delta (tt^{})`$. The binding of the next set of proteins in the opened region may be represented symbolically by a kinetic equation of the type

$$\frac{c}{t}=𝒞(s,t)𝖿\{c\},$$
(13)
where $`c`$ denotes the concentration of the free enzymes in the bath, $`𝖿\{c\}`$ describes the dynamics in absence of any correlations, and $`𝒞(s,t)`$ stands for the initial or early time behaviour of the unzipping correlations. This is equivalent to making the rate constant dependent on the initial correlation. This way of viewing the replication process as coupled system driven by the long-range space-time correlated unzipping dynamics is different from the prevalent contact-based approach.

It is necessary to (a) study long homo-DNA’s under well-controlled pulling force to map out the phase diagram , and (b) probe jointly the correlation and binding kinetics in presence of a pulling force.

A few appealing features of such a scenario are to be noted. Since the critical force is dependent on the structural details of the DNA, which is different for different species, our hypothesis naturally requires species variations of the DnaA-type protein for proper functionality. This protein should be such that it exerts a force slightly below the critical value for that particular DNA. Finally, in the proposed scenario, the dynamic correlation plays the most important role, and therefore a composite structure is not a necessity. In other words, the so-far-elusive hypothesized replisome need not exist as a structural unit.


text:
# References

Optimal shapes of compact strings

Amos Maritan, Cristian Micheletti, Antonio Trovato & Jayanth R. Banavar

International School for Advanced Studies (S.I.S.S.A.), Via Beirut 2-4, 34014 Trieste, Istituto Nazionale di Fisica della Materia and the Abdus Salam International Center for Theoretical Physics, Trieste, Italy

Department of Physics and Center for Materials Physics, 104 Davey Laboratory, The Pennsylvania State University, University Park, Pennsylvania 16802, USA

Optimal geometrical arrangements, such as the stacking of atoms, are of relevance in diverse disciplines . A classic problem is the determination of the optimal arrangement of spheres in three dimensions in order to achieve the highest packing fraction; only recently has it been proved  that the answer for infinite systems is a face-centred-cubic lattice. This simply stated problem has had a profound impact in many areas , ranging from the crystallization and melting of atomic systems, to optimal packing of objects and subdivision of space. Here we study an analogous problem–that of determining the optimal shapes of closely packed compact strings. This problem is a mathematical idealization of situations commonly encountered in biology, chemistry and physics, involving the optimal structure of folded polymeric chains. We find that, in cases where boundary effects  are not dominant, helices with a particular pitch-radius ratio are selected. Interestingly, the same geometry is observed in helices in naturally-occurring proteins.

The problem of placing spheres in three dimensions in order to attain the highest density was first posed by Kepler and has attracted much interest culminating in its recent rigorous mathematical solution . The close-packed hard sphere problem may be restated in an alternative manner, more convenient for numerical implementation, as the determination of the arrangement of a set of points in a given volume that results in the minimum distance between any pair of points being as large as possible . It is notable that the resulting ‘bulk’ optimal arrangement exhibits translational invariance in that, far from the boundaries, the local environment is the same for all points.

In this letter, we introduce a new problem pertaining to the optimal shapes of compact strings. Consider a string (an open curve) in three dimensions. We will utilize a geometric measure  of the curve, the ‘rope-length’, defined as the arc length measured in units of the thickness, which has proved to be valuable in applications of knot theory . The thickness $`\mathrm{\Delta }`$ denotes the maximum radius of a uniform tube with the string passing through its axis, beyond which the tube either ceases to be smooth, owing to tight local bends, or it self-intersects. Our focus is on finding the optimal shape of a string of fixed arc length, subject to constraints of compactness, which would maximize its thickness, or equivalently minimize its rope length.

Following the approach of Gonzalez and Maddocks , who studied knotted strings, we define a global radius of curvature as follows. The global radius of curvature of the string at a given point is computed as the minimum radius of the circles going through that point and all other pairs of points of the string. It generalizes the concept of the local radius of curvature (the radius of the circle which locally best approximates the string) by taking into account both local (bending of the string) and non-local (proximity of another part of the string) effects. For discretized strings the local radius of curvature at a point is simply the radius of the circle going through the point and its two adjoining points. The minimum of all the global radii then defines the thickness, i.e. the minimum radius of the circles going through any triplet of discrete points. This coincides with the previous definition in the continuum limit, obtained on increasing the number of discretized points (assumed to be equally spaced) on the string keeping the string length fixed .

We used several different boundary conditions to enforce the confinement of the string. The simplest ones discussed here are the confinement of a string of length $`l`$ within a cube of side $`L`$ or constraining it to have a radius of gyration (which is the root-mean-square distance of the discretized points from their centre of mass) that is less than a pre-assigned value $`R`$. Even though different boundary conditions influence the optimal string shape, the overall features are found to be robust. Examples of optimal shapes, obtained from numerical simulations, for different ratios of $`l/L`$ and $`l/R`$ are shown in Fig. 1. In both cases, two distinct families of strings, helices and saddles, appear. The two families are close competitors for optimality and different boundary conditions may stabilize one over the other. For example, if optimal strings of fixed length are constrained to have a radius of gyration less than $`R`$, then upon decreasing $`R`$, the string goes from a regime where the trivial linear string is curled into an arc, then into a portion of helix and finally into a saddle. When the string is constrained to lie within a cube of size $`L`$, as $`L`$ decreases first saddles are observed and then helices.

We have also been able to find bulk-like solutions which are not influenced by boundary effects. Such solutions can be obtained by imposing uniform local constraints along the string. On imposing a minimum local density on successive segments of the string (for example, constraining each set of six consecutive beads to have a radius of gyration that is less than a preassigned value $`R`$), we obtained perfectly helical strings, as in Fig. 2, confirming that this is the optimal arrangement. Note that, in close analogy with the sphere-packing problem, the optimal shape displays translational invariance along the chain. In all cases, the geometry of the chosen helix is such that there is an equality of the local radius of curvature (determined by the local bending of the string) and the radius associated with a suitable triplet of non-consecutive points lying in two successive turns of the helix. This is a feature that is observed only for a special ratio $`c^{}`$ of the pitch, $`p`$, and the radius, $`r`$, of the circle projected by the helix on a plane perpendicular to its axis. When $`p/r>c^{}=2.512`$ the global radius of curvature is equal to the local radius with the helix thickness given by $`\mathrm{\Delta }=r(1+p^2/(2\pi r)^2)`$. If $`p/r<c^{}`$, the global radius of curvature is strictly lower than the local radius, and the helix thickness is determined basically by the distance between two consecutive helix turns: $`\mathrm{\Delta }p/2`$ if $`p/r1`$. Optimal packing selects the very special helices corresponding to the transition between the two regimes described above. A visual example is provided by the optimal helix of Fig. 2.

For discrete strings, the critical ratio $`p/r`$ depends on the discretization level. A more robust quantity is the ratio $`f`$, averaged over all the points of the string, of the minimum radius of the circles going through each point and any two non-adjacent points and the local radius. For discretized strings, $`f=1`$ just at the transition described above, whereas $`f>1`$ in the ‘local’ regime and $`f<1`$ in the ‘non-local’ regime. In our computer-generated optimal strings, $`f`$ differed from unity by less than a part in a thousand.

It is interesting to note that, in nature, there are many instances of the appearance of helices. For example, many biopolymers such, as proteins and enzymes, have backbones which frequently form helical motifs. (Rose and Seltzer  have used the local radii of curvature of the backbone as input in an algorithm for finding the peptide chain turns in a globular protein.) It has been shown  that the emergence of such motifs in proteins (unlike in random heteropolymers which, in the melt, have structures conforming to gaussian statistics) is the result of the evolutionary pressure exerted by nature in the selection of native state structures that are able to house sequences of amino acids which fold reproducibly and rapidly  and are characterized by a high degree of thermodynamic stability . Furthermore, because of the interaction of the amino acids with the solvent, globular proteins attain compact shapes in their folded states.

It is then natural to measure the shape of these helices and assess if they are optimal in the sense described here. The measure of $`f`$ in $`\alpha `$-helices found in naturally-occurring proteins yields an average value for $`f`$ of $`1.03\pm 0.01`$, hinting that, despite the complex atomic chemistry associated with the hydrogen bond and the covalent bonds along the backbone, helices in proteins satisfy optimal packing constraints (for the measure of $`f`$ we considered $`\alpha `$-helices extracted from the unrelated proteins 1erv, 1beo and 2end in the Protein Data Bank). This result implies that the backbone sites in protein helices have an associated free volume distributed more uniformly than in any other conformation with the same density. This is consistent with the observation  that secondary structures in natural proteins have a much larger configurational entropy than other compact conformations. This uniformity in the free volume distribution seems to be an essential feature because the requirement of a maximum packing of backbone sites by itself does not lead to secondary structure formation . Furthermore, the same result also holds for the helices appearing in the collagen native state structure, which have a rather different geometry (in terms of local turn angles, residues per turn and pitch ) from average $`\alpha `$-helices. In spite of these differences, we again obtained an average $`f=1.01\pm 0.03`$ (Fig. 3), very close to the optimal situation.

Acknowledgements This work was supported by INFN, NASA and The Donors of the Petroleum Research Fund administered by the American Chemical Society. A.T. thanks the Physics Department of Università degli Studi di Padova, Padova, Italy, for its hospitality.


text:
# Dependence of Gas Phase Abundances in the ISM on Column Density

## 1 Introduction

The study of elemental abundances in the interstellar medium has a long history, both in the optical and with space-based UV spectrographs, such as Copernicus and the Hubble Space Telescope (HST). These studies show that most elements have apparent abundances relative to hydrogen substantially below those in the Sun. This is interpreted as due to depletion of ions onto dust grains.

Savage & Sembach (1996a) summarized high-quality abundance studies based on HST data. They concluded that the depletion of an element depends on the nature of the interstellar medium along the sightline. In cold, dense disk gas, S has near-solar abundance, while Si, Mg, Mn, Cr, Fe and Ni are progressively more depleted. In warm-disk gas, the gas-phase abundances appear higher than in the cool gas, while in halo gas they are even larger.

Jenkins, Savage & Spitzer (1986) used Copernicus to measure accurate column densities for Mg II, P II, Cl I, Cl II, Mn II and Fe II, for stars with distances, $`d`$, between 0.1 and 2 kpc. They measured N(H I) from the damping wings of Ly$`\alpha `$ and found a relation between the ionic abundance relative to H I, $`A`$, and the average density, $`n_\mathrm{H}`$ ($``$ N(H I)/$`d`$). An explanation was provided by Spitzer (1985), who suggested that the observed depletion in a line of sight depends on the fraction of warm and cold (dense) gas encountered, where each component has a specific gas phase abundance of each ion.

Crinklaw, Federman, & Joseph (1994) observed Ca II and Ti II toward 12 stars with distances between 150 and 650 pc. The slope of the relation between log\[$`n_\mathrm{H}`$\] and log\[$`A`$(ion)\] is similar for Ca II and Ti II, but it differs from those for Fe II, Mn II, Mg II and P II. They proposed that Ca II and Ti II mostly occur in the warm intercloud medium, and that the apparent depletion is very sensitive to the inclusion of high-density regions (with high depletion) in the sightline.

These studies were limited to sightlines with N(H I) in the range 1020 to 7$`\times `$1021 cm-2. To extend this range, we used the high- and intermediate-velocity gas (HVCs and IVCs, see Kuntz & Danly 1996, Wakker & van Woerden 1997, Wakker 2000). In 21-cm emission these stand out in velocity from gas in the Milky Way disk, so that weak components (down to 1018 cm-2) can be measured. At low velocities, parcels of gas with such low column densities can not be separated from other high column density concentrations in the line of sight.

Wakker (2000) analyzed all the published absorption-line data pertaining to the high- and intermediate-velocity gas. Improved H I column densities (and thus improved ion abundances) were determined for about 250 sightlines. Ca II was measured for 97 components, Ti II for 6, Fe II for 13, Mn II for 5, Mg II for 10, and Na I for 34 (low-velocity Na I was taken from Ferlet et al. 1985). Most IVCs have intrinsically near-solar abundance (Wakker 2000). However, HVC complex C has $`A`$$``$0.1 solar (Wakker et al. 1999), and the Magellanic Stream has $`A`$$``$0.25 solar (Gibson et al. 2000). Including these HVC complexes could lower some of the apparent abundances (in 8 cases for Ca II, 4 for Fe II and 10 for Mg II), but no clear effect is visible. \[An anomalously low value for Mg II, in the VHVC probed by Mrk 205, was excluded as this cloud may have low intrinsic abundance, and N(H I) is uncertain.\]

No distances are known for many of the HVCs/IVCs, and they fill only a small fraction of the pathlength. Thus, we looked at the relation between $`A`$ and N(H I), rather than $`n_\mathrm{H}`$. Unexpectedly, we found tight anticorrelations of N(H I) with the abundances of Mg II, Ca II, Ti II, Mn II and Fe II. A minimal correlation is expected because both quantities are integrated along lines of sight with varying physical conditions. Variations include ($`a`$) the ionization of H, ($`b`$) widely differing histories of grain destruction by interstellar shocks, and ($`c`$) considerable differences of the interstellar radiation field impinging on the gas. Ca II and Na I are not even the dominant stage of ionization in neutral regions, so the radiation field is of critical importance for their ionization fraction. In §2 we show the correlations between column density and abundances, while in §3 we discuss some implications.

## 2 Correlations between abundances and N(H I)

### 2.1 The data

Figure 1 shows scatter plots of log\[N(H I)\] vs log\[$`A`$(ion)\], the gas-phase abundance relative to H I, for Mg II, Ca II, Ti II, Mn II, Fe II and Na I. The least-squares fit and dispersion (labeled “rms”) are given in each panel. Symbols indicate the origin of the datapoints, as detailed in the figure caption. The Local ISM points (open squares) were not used in the fits (see §2.5 below). The high-quality data obtained with the “Goddard High Resolution Spectrograph” (GHRS) are presented separately, in order to show the difference between halo and disk sightlines more clearly (§2.5).

The currently preferred oscillator strengths of Mn II-1197, 1201 and Mg II-1239, 1240 are 0.20 and 0.67 dex higher than those used by Jenkins et al. (1986) (see Savage & Sembach 1996a). We therefore corrected their Mn II and Mg II column densities downward by these amounts.

As the highest abundances are associated with high-velocity gas, we might be seeing a manifestation of the Routly & Spitzer (1952) effect (this is the increase in the ratio N(Ca II)/N(Na I) with LSR velocity). This effect is interpreted as showing that Ca is less depleted at higher peculiar velocities. However, it is based on nearby ($`<`$100 pc), low-velocity ($`<`$20 km s-1) gas (e.g. Vallerga et al. 1993). For the HVCs/IVCs the LSR velocity is not a good measure of the peculiar velocity relative to their surroundings. Further, the IVCs for which both Na I and Ca II have been measured do not show the effect (Wakker 2000).

We now discuss some implications of the $`A`$ vs N(H I) relations.

### 2.2 Result 1: small dispersions

For Mg II, Mn II and Fe II the standard deviations of log\[$`A`$(ion)\] from the least-squares fit at a given N(H I) are only 0.27 in the log. That is less than a factor 2 either way, even though the correlation extends over 4 orders of magnitude in N(H I) and over 2 orders of magnitude in $`A`$(Fe II). The standard deviations are slightly larger for Ca II and Ti II (0.4 in the log, or a factor 2.6 either way). The largest scatter is seen for Na I, where it is 0.52 in the log or a factor 3.5 either way.

The larger scatter for Ca II and especially Na I may be related to the fact that, unlike the other ions, these are not the dominant ionization stage in the diffuse ISM, so that ionization effects may play a larger role. However, Ti II is the dominant ionization stage, and it shows a scatter that is comparable to that of Ca II.

Ferlet et al. (1985) proposed a slope near zero for the relation between Na I abundance and N(H I) for the low-velocity gas. This has been used extensively to estimate N(H I) from N(Na I). The lower column density IVC data tend to show slightly higher abundances. However, the spread in $`A`$(Na I) at any given value of N(H I) is rather large. For a given N(Na I), the derived N(H I) will be a factor $`>`$3.5 off either way in 33% of the cases. The most deviant points in the Na I diagram differ from the mean relation by a factor $`>`$20. Thus, although on average N(Na I)/N(H I) is independent of N(H I), N(Na I) is a rather poor predictor of N(H I).

### 2.3 Result 2: differing slopes

The slopes of the correlations differ substantially: $``$0.78$`\pm `$0.04 for Ca II, $``$0.69$`\pm `$0.08 for Ti II, $``$0.59$`\pm `$0.04 for Fe II, $``$0.39$`\pm `$0.04 for Mn II, $``$0.24$`\pm `$0.04 for Mg II, and $``$0.16$`\pm `$0.06 for Na I. We find slopes of $``$0.63$`\pm `$0.12 and $``$0.66$`\pm `$0.23 for low- and high-velocity Ca II; for Fe II we find $``$0.54$`\pm `$0.05 and $``$0.70$`\pm `$0.07, respectively. For other elements there are insufficient datapoints to make separate fits. The errors are larger because of the reduced range in N(H I), but the resulting slopes are the same to within the formal error. This shows that on average the high column density, low-velocity disk gas behaves in the same manner as the low column density high-velocity halo gas.

### 2.4 Result 3: no obvious ionization effects

The ionization potentials (I.P.) of Fe II, Mn II, and Mg II are 16.2, 15.6, and 15.0 eV, respectively (compared to 13.6 eV for H I), so these ions can co-exist with H+. The Fe III/Fe II ratio should not depend on N(H I), but only on the surrounding radiation field (which only depends on location). Thus, if the interstellar radiation field has similar values in the sightlines that were observed, and if H+ were only present in shells around H I cores, then the fraction of H+ should increase with decreasing N(H I), resulting in an overestimate of the abundance, when calculated as N(ion)/N(H I). This would produce an upturn at low N(H I) in the $`A`$ vs N(H I) relation. However, such an effect is not obvious. Ca II (I.P. 11.9 eV) should ionize more easily than H, leading to a downturn at low N(H I). This is also not seen.

### 2.5 Result 4: environmental dependence

For Mg II, Mn II and Fe II the GHRS disk points (open circles in Fig. 1) scatter around those of Jenkins et al. (1986). GHRS halo sightlines (filled circles) tend to lie above the average correlation. The small scatter in $`A`$(Fe II) remarked upon by Savage & Sembach (1996a) is consistent with the fact that their halo points span only a factor 10 in N(H I). However, some more recent halo points (filled triangles) have higher Fe II/H I ratios, as do some HVCs.

The low N(H I) Local ISM (LISM; open squares) shows abundances far below halo points at the same N(H I), although they are near the high end for disk gas and the observed depletion is consistent with the average Local ISM density of 0.1 cm-3. We surmise that both the Sun and the LISM probes lie within a cloud, so that the LISM paths do not pass through all of it, unlike the case for clouds in the sightlines to more distant stars. Thus, the Local ISM points suggest that density is the relevant physical parameter determining depletion. Over sufficiently long pathlengths column density can substitute because it is dominated by the densest regions.

The largest depletions are seen toward $`\zeta `$ Oph and $`\xi `$ Per, which were used by Sembach & Savage (1996) as the archetypical cold disk sightlines. These have high N(H I) and $`A`$(Mg II) and $`A`$(Fe II) lie below the average relation. Including N(H2) would move these points 0.2 and 0.4 dex to the right and down, respectively, shifting them $``$0.1 and 0.2 dex further below the average.

We thus see the environmental dependence found by Sembach & Savage (1996): disk sightlines tend to have low gas-phase abundances, while halo sightlines show high gas-phase abundances. We can understand this in large part as a column density dependence, as disk sightlines tend to have higher column density than halo sightlines.

## 3 Discussion

The general trend that the gas phase abundances of refractory ions decreases with N(H I) is, of course, well known; the ions are removed from the gas phase by sticking to grains. However, the low dispersions of the depletions at a particular value of N(H I), continuing to low values to N(H I), are remarkable. The depletion patterns provide a diagnostic of the density fluctuation within the ISM, since the depletions require encounters between the grains and ions. The low dispersions of $`A`$ at a given N(H I) suggest additional properties of the ISM:

1. Large values of N(H I) arise from one or a few regions of relatively large local density rather than from a superposition of many regions of lower density. In this way, there is unexpected physical meaning to the integrated column density: a good correlation with the mean local density of the typical H atom along the sightline. If large values of N(H I) were built up from many regions of low density along the line of sight, the averaged $`A`$ would correspond to the lower densities. The association of large N(H I) with a large local density is strengthened by the fact that the anticorrelation of gas phase abundances with $`n_\mathrm{H}`$ is better than with N(H I) (Savage & Bohlin 1979).

2. The low dispersions in log\[$`A`$\] suggest that grains are not destroyed in strong shocks, as some models have suggested. If such destruction occurred, a few points with large gas phase abundances would be expected at large N(H I). Rather, the detachment of a refractory element from grains seems to be a continuous function of the mean density, and is common enough to produce a low dispersion about the mean. Grain destruction may occur by means of many small shocks with a broad distribution of velocities.

3. Both the rates of ions sticking to grains and becoming detached from them (probably in atomic form) must be rapid enough to achieve rather good statistical equilibrium. These balancing rates lead to a particular value of $`A`$ for a given density; otherwise, the past history of the H I along the line of sight would strongly influence $`A`$.

4. In view of (3) above, we think that the best possibility is that attachment takes place within regions of far greater than average density distributed along each line of sight. Encounters between grains and ions are very slow if the mean density of the diffuse ISM is assumed (see below). Possibly magnetohydrodynamic waves create non-thermal motions of the ions relative to the grains, since the charge per mass of the grains is vastly less that for the ions. As a result, ions encounter more dust grains and the net attachment rate increases.

Under the assumption that relative motions were thermal, Weingartner & Draine (1999) calculated average gas phase abundances of a highly depleted species (Ti II) for large N(H I). In their model, dense ($`n_\mathrm{H}`$$``$103 cm-3) molecular clouds contain half the mass of the ISM, 30% is in cool, neutral clouds ($`n_\mathrm{H}`$=30 cm-3) and 20% is in warm neutral clouds ($`n_\mathrm{H}`$=0.4 cm-3). The attachment of the Ti onto grains occurs in the molecular and cool clouds. This model shows that the observed $`A`$(Ti II) for large N(H I) can only be achieved by circulation into very dense regions. To achieve rapid attachment within the diffuse ISM, such circulation seems necessary, operating at both low and high column densities. An additional condition for rapid attachment of ions is that the gas be cold as well as dense, in order to provide efficient focusing by the Coulomb force onto negatively charged grains.

We appreciate comments and clarifications by an anonymous referee, Alex Lazarian, Chris Howk, and Bruce Draine. B.P.W. is partially supported by NASA grant NAG5-9179.

Figure 1. Correlation plots of the abundances of Na I, Mg II, Ca II, Ti II, Mn II and Fe II vs N(H I). The solid lines show the least-squares fits, whose coefficients are given by the label in the lower left corner. “rms” is the rms of the residual, also indicated by the two parallel dotted lines. “$`\rho `$” is the correlation coefficient. Filled circles are for the HVC and IVC data from Wakker (2000). If measurement errors were given for N(ion), error bars are shown; usually no errors are available. Crosses show data for low-velocity gas from Ferlet et al. (1985; Na I), Jenkins et al. (1986; Mg II, Mn II and Fe II) and Crinklaw et al. (1994; Ca II, Ti II). Errors are usually smaller than the size of the dot. The right panels for Mg II, Mn II and Fe II separately show measurements with relatively low systematic errors obtained with the GHRS by Savage et al. (1993), Sofia et al. (1994), Cardelli et al. (1995), Spitzer & Fitzpatrick (1995), Savage & Sembach (1996b) and Sembach & Savage (1996). For these data, halo sightlines are shown by filled circles, disk sightlines by open circles. Filled/open triangles show halo/disk sightlines from Fitzpatrick & Spitzer (1997; after combining their fit components into the 5 H I components visible in the spectrum) and Howk et al. (1999). Open squares show Local ISM sightlines from Linsky et al. (1995), Hébrard et al. (1999) and Howk (priv. comm.); these were not used in the fits, see §2.4. For all low-velocity data the error bars are smaller than the symbol sizes. The horizontal lines show the solar reference abundance, from Anders & Grevesse (1989).

Figure 1