text:
# Molecular hydrogen formation in the interstellar medium

## 1 Introduction

Molecular hydrogen is the most abundant molecule in the universe and dominates the mass budget of gas in regions of star formation. H2 is also an important chemical intermediate in the formation of larger species and can be an important gas coolant in some conditions, particularly in the early universe. Yet, despite its importance, the H2 formation process is still not well understood. Observationally, it has been shown that H2 can be efficiently formed over a wide range of temperatures ( Jura 1974; Tielens & Hollenbach 1985a, 1985b; Hollenbach & McKee 1979). Theoretically, Gould and Salpeter (1963) showed the inefficiency of H2 formation in the gas phase and concluded that recombination of physisorbed H on “dirty” ice surfaces is efficient between 10 and 20K. Hollenbach and Salpeter (1970; 1971), recognizing that this small temperature range presents a problem, considered reactions involving H atoms bound to dislocations and impurities with energy exceeding normal physisorbed energies, and obtained a recombination efficiency $``$ 1 up to a critical temperature between 25 and 50K. Goodman (1978) calculated the quantum and thermal mobility of the atoms on graphite grains assuming that these atoms could only be physisorbed. Many studies focussed for various reasons on icy surfaces where H is physisorbed (Buch & Zhang 1991; Takahashi et al. 1999). However, most astrophysically relevant material (e.g., silicates, graphite) can bind H also in chemisorbed sites (Barlow & Silk, 1976, Aronowitz & Chang 1980; Leitch-Devlin & Williams 1984; Klose 1992; Fromherz et al. 1993). On these surfaces, binding can occur in a physisorption layer (E$``$500K) at a distance of some Z$``$3$`\AA `$ as well as in a chemisorption layer (E$``$10000K) deeper into the surface Z$``$0.5$`\AA `$ (Barlow & Silk, 1976; Zangwill, 1988). Recently, Katz et al. (1999) developed a model for H2 formation benchmarked by two sets of experiments. This model considers the atoms bound only in physisorbed sites and diffusing only thermally on the surface, colliding and recombining to form molecules. As for the ice models studied earlier, this allows molecule formation only below 15K for olivine grains and 20K for carbon grains which contradicts the ISM observations. Perusing these models, it is clear that the process of the H2 formation is governed by the binding of atomic H to the surface and the concomitant mobility of these atoms (Leitch-Devlin & Williams 1984; Tielens & Allamandola 1987). Surface diffusion can occur through quantum mechanical tunneling (at low temperatures) and through thermal hopping (at high temperatures). For a proper description of molecular hydrogen formation in the ISM both types of binding diffusion processes have to be taken into account (Cazaux & Tielens 2002).

## 2 Model

Recently, we have developed a model for H2 formation on grain surfaces based upon Langmuir kinetics where species accrete, migrate, react and the product species, possibly assisted by thermal energy, evaporate. This model is based on two main points: 1) The atomic mobilities are due to a combination of quantum tunneling and thermal diffusion and this is controlled by the temperature of the grain and the mass of the species. 2) Atoms can bind to the surface in two energetically different sites: a chemisorption site and a physisorption site. These interaction energies set the migration between the different sites and the reactions among the species. We have compared the results of our model with laboratory experiments on molecular hydrogen formation on olivine surfaces (Pirronello et al. 1999; Katz et al. 1999) to determine the relevant parameters. Here, we extrapolate this model to study H2 formation under astrophysically relevant conditions (e.g., low accretion rates and long timescales). We consider three different populations: physisorbed H, chemisorbed H and physisorbed H2. The surface concentrations of these species are described by three rate equations. Using the surface concentration of these species, the H2 desorption rate can be determined, which contains two contributions. First order desorption occurs when the grain temperature is high enough to allow evaporation of previously formed molecules. Second order desorption occurs when two atoms encounter and the newly formed molecule is directly released into the gas phase. We define $`H_P`$, $`H_C`$ and $`H_2`$ as the physisorbed H, chemisorbed H and molecular hydrogen concentrations, respectively; their evaporation rate are written as $`\beta _{HP}`$, $`\beta _{HC}`$, $`\beta _{H_2}`$. The different mobilities to go from a site $`i`$ to a site $`j`$ are given by $`\alpha _{ij}`$ where $`i`$ and $`j`$ are either a physisorbed site (P) or a chemisorbed site (C), and where the barrier between two sites is assumed to be square. However, the shape of the barrier is found out to be unimportant (see below). $`\mu `$ is the fraction of H2 which stays on the surface upon formation. The rate equation can then be written as (for details see Cazaux & Tielens 2002),

$$\dot{H}_P=F(1H_PH_2)\alpha _{pc}H_P2\alpha _{pp}H_{P}^{}{}_{}{}^{2}+\alpha _{cp}H_C(1H_P)\beta _{H_P}H_P$$
(1)
$$\dot{H}_C=\alpha _{pc}H_P(1H_C)\alpha _{pc}H_PH_C\alpha _{cp}H_C2\alpha _{cc}H_{C}^{}{}_{}{}^{2}\beta _{H_C}H_C$$
(2)
$$\dot{H}_2=\mu (2\alpha _{pp}H_{P}^{}{}_{}{}^{2}+\alpha _{pc}H_PH_C+\alpha _{cp}H_CH_P+\alpha _{cc}H_{C}^{}{}_{}{}^{2})\beta _{H_2}H_2$$
(3)
Where F is accretion rate in unit of monolayers per second (ML.$`s^1`$). The first order H2 desorption rate is given by:

$$k_1=\beta _{H_2}H_2$$
(4)
For second order desorption several processes contribute depending on whether physisorbed, chemisorbed atoms or both are involved:

$$k_2=(1\mu )(\alpha _{pc}H_PH_C+\alpha _{cp}H_CH_P+\alpha _pH_P^2+\alpha _{cc}H_C^2)$$
(5)
$`\beta _i`$, the desorption rate for a population $`i`$ is written $`\beta _i=\nu _i\mathrm{exp}\left[\frac{E_i}{kT}\right]`$, where $`\nu _i`$ is the frequency factor of species $`i`$, perpendicular to the surface, and $`E_i`$ the desorption energy of species $`i`$. The total H2 formation rate, R, is the sum of the first and second order terms, and the H2 recombination efficiency – the fraction of the accreted hydrogen which desorbs as H2– is defined by $`ϵ=\frac{2R}{F}`$ . The parameters in this model are summarized in Table 1. They have been discussed extensively elsewhere (see Cazaux & Tielens 2002). We have integrated this set of time dependent equations using a Runge-Kutta method with adaptative stepsize control, (Odeint, Subroutine package from Numerical Recipes 1992) until steady state was achieved. Under astrophysically relevant conditions, steady state is generally achieved. These steady state results are shown as a function of T in fig. 1 for three different H-fluxes. All curves are characterized by a very low efficiency at low temperatures, a steep rise around $``$5K – whose location depends on the actual H-flux –, an efficiency of unity until some 25K, followed by a gradual decline to $``$ 0.2 at $``$ 100K, and a final drop to very low efficiency around 300K, whose location again depends on the actual H accretion rate. In steady state, the system of equations 1-3 can be simplified by setting the time derivatives equal to zero. We can now discern two different temperatures regimes – low temperatures (T$``$100K) and high temperatures (T$``$100K) – where different processes dominate H2 formation. In these limits, the set of steady state equations can be further simplified and this yields much insight in the numerical results.

### 2.1 Physisorbed H

At low temperatures, by far the most efficient way to form H2 is when a physisorbed atom reacts with a chemisorbed atom. Other routes forward H2 as well as the evaporation of chemisorbed atoms are negligible. Realizing that $`H_P`$ is always small compare to 1 (physisorbed H quickly migrate to a chemisorbed well), the steady state equations can be simplified to,

$$F(1H_2)\alpha _{pc}H_P\beta _{H_P}H_P=0$$
(6)
$$\alpha _{pc}H_P(12H_C)=0$$
(7)
$$\mu (\alpha _{pc}H_PH_C)\beta _{H_2}H_2=0$$
(8)
The different populations can easily be determined in this system (e.g., $`H_C`$=$`\frac{1}{2}`$, $`H_P`$=$`(\frac{\mu \alpha _{pc}}{2\beta _{H_2}}+\frac{\alpha _{pc}}{F}+\frac{\beta _{H_2}}{F})^1`$ and $`H_2`$=$`\frac{\mu \alpha _{pc}}{2\beta _{H_2}}H_P`$), the H2 desorption rate is calculated and we deduce the recombination efficiency:

$$ϵ_{H_2}=\left(\frac{\mu F}{2\beta _{H_2}}+1+\frac{\beta _{H_P}}{\alpha _{pc}}\right)^1$$
(9)
We can recognize three different regimes in this temperature range. At really low temperatures, below some 10K, the high mobility of H atoms due to tunneling permits the recombination of incoming H atoms. The H2 formed stays on the surface, since the temperature is not high enough to allow evaporation and blocks further H-atom accretion. Then, the recombination efficiency is $`ϵ_{H_2}`$=$`\left[\frac{\mu F}{2\beta _{H_2}}\right]^1`$. At higher temperatures, (6-25K) the desorption rate depends only on the flux because all the incoming H atoms adsorb, recombine and desorb as H2. In this temperature regime, hydrogen accretes into a physisorbed site but quickly drops into a chemisorbed well. If this chemisorbed site is empty, the H atom will be trapped. If an H atom is already present, reaction will occur and the product will evaporate quickly either upon formation or thermally assisted; thus, $`H_C`$=$`\frac{1}{2}`$ and $`ϵ_{H_2}`$=1. Between $``$25K up to $``$100K, evaporation of physisorbed atoms competes with recombination, and the desorption rate decreases considerably.

### 2.2 Chemisorbed H

With increasing temperature, the physisorbed (at T$``$100K) and then the chemisorbed atoms(at T$``$300K) will start to evaporate. The most efficient reaction to form H2 is then the “collision” of two chemisorbed atoms. The system of equations reduces to,

$$F\alpha _{pc}H_P\beta _{H_P}H_P=0$$
(10)
$$\alpha _{pc}H_P2\alpha _{cc}H_C^2\beta _{H_C}H_C=0$$
(11)
$$\mu (\alpha _{cc}H_{C}^{}{}_{}{}^{2})\beta _{H_2}H_2=0$$
(12)
and the H2 recombination efficiency can be written as,

$$ϵ=\left(1+\frac{\beta _{H_P}}{\alpha _{pc}}\right)^1\xi $$
(13)
with $`\xi `$ the correction factor at high temperatures,

$$\xi =\left(1+\frac{\beta _{H_C}^2\beta _{H_P}}{2F\alpha _{pc}\alpha _{cc}}\right)^1,$$
(14)
which reflects the evaporation of chemisorbed H. This factor is flux dependent.

### 2.3 General expression for the recombination efficiency

The above expression for the H2 recombination efficiency can be combined into one general expression valid at any temperatures:

$$ϵ_{H_2}=\left(1+\frac{\mu F}{2\beta _{H_2}}+\frac{\beta _{H_P}}{\alpha _{pc}}\right)^1\xi $$
(15)
The mobilities, $`\alpha _{PC}`$ and $`\alpha _{CC}`$, are due to a combination of quantum mechanical tunneling and thermal hopping but considering the temperatures at which these parameters play a role, these mobilities are dominated by thermal hopping. Therefore, the H2 recombination efficiency is independent of the width of the considered barriers, and we can approximate $`\xi `$ and $`\frac{\beta _{H_P}}{\alpha _{pc}}`$ by the following expressions:

$$\xi =\left(1+\frac{\nu _{H_C}\mathrm{exp}\left[\frac{1.5E_{H_C}}{kT}\right]\left(1+\sqrt{\frac{E_{H_C}E_S}{E_{H_P}E_S}}\right)^2}{2F}\right)^1$$
(16)
$$\frac{\beta _{H_P}}{\alpha _{pc}}=\frac{1}{4}\left(1+\sqrt{\frac{E_{H_C}E_S}{E_{H_P}E_S}}\right)^2\mathrm{exp}\left[\frac{E_S}{kT}\right]$$
(17)
The parameters in these expressions have been determined from experimental data (Table 1; Cazaux & Tielens 2002). The H2 recombination efficiency for three different flux is reported figure 1. In astrophysical environments, the recombination rate is written:

$$R_{H_2}=\frac{1}{2}n_Hv_Hn_d\sigma _dϵ_{H_2}S_H(T)$$
(18)
where $`n_H`$ and $`v_H`$ are the number density and the thermal velocity of H atoms in the gas phase, $`n_d\sigma _d`$ is the total cross section of the interstellar grains and $`S_H(T)`$ is the sticking coefficient of the H atoms which can depend on temperature.

## 3 Discussion

Our study reveals the presence of two distinct regimes of H2 formation, which reflect directly the presence of two types of atomic H binding sites. At low temperature (T$``$100K), H2 formation involves migration of physisorbed H atoms. At higher temperatures (T$``$100K), H2 formation results from chemisorbed H recombination. The presence of these two types of binding sites allow H2 formation to proceed relatively efficiently even at elevated temperatures. The study of Hollenbach and Salpeter (1971) focused on icy surfaces on which H can only physisorb. As a result, H2 formation ceased at temperatures in excess of $``$20K. Recognizing this problem, Hollenbach and Salpeter involved the presence of enhanced binding sites on the ice with ill-determined parameters. These sites allowed H2 formation to proceed up to some 75K. Since their study, it has become abundantly clear that interstellar grains are not covered by ice in the diffuse interstellar medium (Whittet et al. 1983, 1988). Silicate and graphitic surfaces are now widely accepted as astrophysically relevant grain surfaces (Mathis 1991) and those surfaces intrinsically possess enhanced binding sites; e.g., chemisorbed sites. The parameters of these chemisorbed sites have not yet been well determined because experiments have focused on low temperature H2 formation (Pirronello et al. 1997; Katz et al. 1999). The values adopted in this study are however quite reasonable and illustrate the efficiency of H2 formation at elevated temperatures well. When future experiments determine the values of the parameters involved ($`E_{H_C}`$, $`E_{H_P}`$, $`\mu `$), the results can be directly adjusted. Inside dense clouds, interstellar grains are covered by ice. Of course, in such environments, almost all hydrogen is already in molecular form and H2 formation is perhaps only of academic interest. Nevertheless, we note that in such environments, molecular hydrogen formation may proceed mainly through H-abstraction from molecules such as H2S and N2H2 (Tielens and Hagen 1982). In a sense, these species act as “chemisorption” sites for hydrogen. Migrating atomic H may tunnel through the reaction barriers involved and form H2. Eventually, these ice covered grains are transported back into the diffuse ISM, when the molecular cloud is disrupted. Photodesorption and sputtering in strong shocks quickly remove their ice on a timescale of some $`10^6`$ yr (Tielens and Hagen 1982, Jones et al. 1994, Draine and Salpeter 1979, Barlow 1978). At that point, molecular hydrogen formation is again governed by the properties of bare grain surfaces. Similarly, any (thin) layer accreted in the diffuse interstellar medium will be quickly sputtered in even a modest velocity ($``$30km s-1) shock (Jones et al. 1994).
Finally, the formation efficiency of molecular hydrogen will also depend on the sticking coefficient of H atoms colliding with the grain. In our model and the formulae derived, the sticking coefficient is subsumed in the incident flux, F. Astrophysical studies of sticking of H on grain surfaces have concentrated on physisorbed interactions and the sticking coefficient is $``$1 at low temperatures and decreases with increasing temperature to about 0.3 at T=300K ( Hollenbach & Salpeter 1970; Hollenbach & McKee 1979; Burke & Hollenbach 1983; Leitch-Devlin & Williams 1985). However, if the interaction occurs through much stronger chemisorption, then the sticking coefficient might be large even at high temperatures (Tielens & Allamandola 1987; Duley & Williams 1984).

### 3.1 Summary and conclusions

Recently, we have modelled molecular hydrogen formation on grain surfaces. This model consider hydrogen atoms bound to the surface at two energy levels (i.e. chemisorption and physisorption). The H mobility from one site to another is a combination of tunelling effect and thermal diffusion. This model has been experimentally benchmarked (Pirronello et al. 1997a, 1997b, 1999) and the relevant surface characteristics have been determined. These characteristics allow us to extend our model for H2 formation under astrophysically relevant conditions. The results show efficient H2 formation from $``$6K to $``$300K. The different processes involved in H2 formation at different temperatures has been discussed. Until about 100K, H2 forms by recombination of a physisorbed H with a chemisorbed H and is highly efficient. At higher temperatures, when physisorbed atoms evaporate quickly, the recombination of two chemisorbed atoms is required to form H2. H2 formation is then less efficient, $`ϵ_{H_2}`$$``$0.2. The parameters involved in H chemisorption and H2 formation at high temperatures are presently not well known. The adopted values are very reasonable and the gross characteristics – H2 formation at high temperatures – are undoubtedly correct. Nevertheless, future experiments are very important to determine the maximum temperature to which H2 formation in the ISM can occur.


format: text

text:
# Active vs passive scalar turbulence

## Abstract

Active and passive scalars transported by an incompressible two-dimensional conductive fluid are investigated. It is shown that a passive scalar displays a direct cascade towards the small scales while the active magnetic potential builds up large-scale structures in an inverse cascade process. Correlations between scalar input and particle trajectories are found to be responsible for those dramatic differences as well as for the behavior of dissipative anomalies.

The transport of scalar fields by turbulent flows is a common physical phenomenon. The dynamics of the advecting velocity field often depends on the transported field, which is then dubbed active. That is, for instance, the case of temperature, affecting the velocity via buoyancy forces. Conversely, situations where the flow is not influenced by the scalar, e.g. for the concentration of a dilute tracer, are referred to as passive. Progress has been recently made for the passive case (see Ref. FGV01  and references therein). These studies put renewed emphasis on the known field-particle duality and long known concepts as intermittency, dissipative anomaly, direct and inverse cascades, usually defined in terms of field characteristics (the Eulerian description), can now be elegantly explained in terms of the statistical properties of particle trajectories (the Lagrangian description). In this Letter we address the problem of relating the Eulerian properties to the Lagrangian ones for active scalar transport.

The active field $`a(𝒙,t)`$ advected by an incompressible flow $`𝒗(𝒙,t)`$ is governed by the transport equation

$$_ta+𝒗\mathbf{}a=\kappa \mathrm{\Delta }a+f_a.$$
(1)
The passive field $`c(𝒙,t)`$ evolves in the same flow as

$$_tc+𝒗\mathbf{}c=\kappa \mathrm{\Delta }c+f_c.$$
(2)
The two scalars have the same diffusion coefficient $`\kappa `$ while $`f_a`$ and $`f_c`$ are independent realizations of the same random forcing, with characteristic lengthscale $`l_f`$. The difference between active and passive scalars resides in their relationship with the flow $`𝒗`$: the active field enters the velocity dynamics whereas the passive one does not. Here, the evolution equation for $`𝒗`$ is

$$_t𝒗+𝒗\mathbf{}𝒗=\mathbf{}p\mathrm{\Delta }a\mathbf{}a+\nu \mathrm{\Delta }𝒗.$$
(3)
Eqs. (1) and (3) are the two-dimensional magnetohydrodynamics equations (see, e.g., Ref. Bis ), with the following glossary: $`a`$ is the scalar magnetic potential; the magnetic field is $`𝒃=(_2a,_1a)`$; the coupling $`\mathrm{\Delta }a\mathbf{}a`$ is the Lorentz force $`(\mathbf{}\times 𝒃)\times 𝒃`$; the condition $`\mathbf{}𝒗=0`$ is ensured by the pressure term $`\mathbf{}p`$ and $`\nu `$ is the kinematic viscosity.

Besides its physical relevance (see e.g. Bis ), this example of active scalar is particularly interesting because of the conspicuous differences to its passive counterpart. Namely, whilst $`a`$ undergoes an inverse cascade process, i.e. it forms structures at increasingly larger scales Bis2 , $`c`$ cascades downscale (see Fig. 1). As a consequence, the dissipation of active scalar fluctuations by molecular diffusivity $`ϵ_a=\kappa |\mathbf{}a|^2`$ vanishes in the limit $`\kappa 0`$ : no dissipative anomaly for the field $`a`$. Since the scalar variance is injected at a constant rate $`F_0`$ by the source $`f_a`$, we have that $`e_a(t)=\frac{1}{2}a(𝒙,t)^2𝑑𝒙`$ increases linearly in time as $`\frac{1}{2}F_0t`$. On the contrary, passive scalar dissipation $`ϵ_c=\kappa |\mathbf{}c|^2`$ equals the input $`\frac{1}{2}F_0`$ and holds $`c`$ in a statistically stationary state (see Fig. 1).

The probability distribution function (pdf) of $`a`$ is Gaussian, with zero mean and variance $`F_0t`$ (see Fig. 2). This is a straightforward consequence of the vanishing of active scalar dissipation. Indeed, by averaging over the forcing statistics one obtains: $`_ta^{2n}=n(2n1)F_0a^{2n2}`$ (odd moments vanish by symmetry) whose solutions are the Gaussian moments: $`a^{2n}=(2n1)!!(F_0t)^n`$. This may be contrasted with the single-point pdf of $`c`$ which is stationary and supergaussian (see Fig. 2), as it generically happens for a passive field sustained by a Gaussian forcing in a rough flow (see, e.g., Ref. FGV01 ).

We now focus on the main purpose of this Letter, that is to investigate the statistical properties of active and passive scalar in the Lagrangian setting. The fundamental property which we will exploit is that the equations (1) and (2) can be formally solved in terms of particle propagators. Let $`𝑿(s)`$ be the trajectory of a fluid particle transported by the flow $`𝒗`$ and subject to a molecular diffusivity $`\kappa `$, landing at point $`𝒙`$ at time $`t`$. The particle moves according to the stochastic differential equation $`d𝑿(s)=𝒗(𝑿(s),s)ds+\sqrt{2\kappa }d𝜷(s)`$, where $`𝜷(s)`$ is a two-dimensional brownian motion. The probability density $`p(𝒚,s|𝒙,t)`$ of finding a particle at point $`𝒚`$ and time $`st`$ obeys the Kolmogorov equations (see, e.g., Risk ):

$$_sp(𝒚,s|𝒙,t)\mathbf{}_y[𝒗(𝒚,s)p(𝒚,s|𝒙,t)]=\kappa \mathrm{\Delta }_yp(𝒚,s|𝒙,t),$$
(4)
$$_tp(𝒚,s|𝒙,t)+\mathbf{}_x[𝒗(𝒙,t)p(𝒚,s|𝒙,t)]=\kappa \mathrm{\Delta }_xp(𝒚,s|𝒙,t).$$
(5)
The unusual minus signs in the l.h.s. of (4) are due to the fact that particles move backward in time. The solution of (1) and (2) can be written in terms of the propagator:

$`a(𝒙,t)`$ $`=`$ $`{\displaystyle _0^t}𝑑s{\displaystyle 𝑑𝒚f_a(𝒚,s)p(𝒚,s|𝒙,t)},`$ (6)
$`c(𝒙,t)`$ $`=`$ $`{\displaystyle _0^t}𝑑s{\displaystyle 𝑑𝒚f_c(𝒚,s)p(𝒚,s|𝒙,t)},`$ (7)
as it can be directly checked by inserting (6) and (7) in (1) and (2), respectively, and utilizing (5). It is important to notice that both the active and the passive scalar are expressed in terms of the same propagator $`p`$, since $`a`$ and $`c`$ are advected by the same velocity field. Eq. (6) can be written more compactly as $`a(𝒙,t)=_0^tf_a(𝑿(s),s)𝑑s_X`$, where $`\mathrm{}_X`$ denotes the average over particle trajectories, and similarly for eq. (7).

The interpretation of eqs. (6) and (7) is that the value of the scalar field at $`(𝒙,t)`$ is given by the superposition of the input along all trajectories eventually convening at that point and time. It is worth remarking that, even in the limit of vanishing diffusivity, the particle propagator does not collapse onto a single trajectory, as expected for incompressible rough flows FGV01 . The differences between active and passive fields are due to the correlations between the propagator $`p`$ and the forcing $`f_a`$, which are absent, by definition of passive transport, between $`p`$ and $`f_c`$. Indeed, via the coupling term in eq. (3), $`f_a`$ affects the velocity field and consequently the evolution of the propagator $`p`$. To further clarify the relationship between Lagrangian trajectories and active field forcing, it is necessary to investigate the evolution of the particle propagator. The problem is that $`p`$ evolves backward in time, according to (4) and the condition $`p(𝒚,t|𝒙,t)=\delta (𝒚𝒙)`$ is set at the final time $`t`$. Contrary to an usual forward integration, where $`𝒗`$ and $`p`$ can be advanced in parallel, a brute force backward integration would require the storage of velocity configurations for the whole lapse of integration which is quite unfeasible. To overcome this problem we devised a fast and low-memory demanding algorithm to integrate numerically eqs. (1), (3) and (4) alain .

Let us sketch the basic idea of the algorithm. We proceed inductively and first consider the case $`t=t_1=\mathrm{\Delta }t`$, i.e. a single step. The fields $`a(0)`$ and $`𝒗(0)`$ and the final condition $`p(t_1)`$ are given as input (we omit spatial dependences). The fields $`a(t_1)`$ and $`𝒗(t_1)`$ are obtained by a forward integration whose numerical cost is denoted $`c_f`$. The propagator is then evolved backward, at a cost denoted $`c_b`$, to obtain $`p(0)`$. The total cost of the procedure is $`C_1=c_f+c_b`$. Consider now the case $`t=t_N=(2^N1)\mathrm{\Delta }t`$. Our aim is to show that the algorithm for $`2^N1`$ steps can be split into two successive applications with $`2^{N1}1`$ steps and thus lends to a recursive procedure. Integrating forward $`a(0)`$, $`𝒗(0)`$ for $`2^{N1}`$ steps give $`a(t_{N1}+\mathrm{\Delta }t)`$ and $`𝒗(t_{N1}+\mathrm{\Delta }t)`$ at a cost $`2^{N1}c_f`$. Those fields, together with $`p(t_N)`$, are the input for the algorithm in the interval $`[t_{N1}+\mathrm{\Delta }t,t_N]`$, which gives $`p(t_{N1}+\mathrm{\Delta }t)`$ at a cost $`C_{N1}`$. A backward step is finally needed to obtain $`p(t_{N1})`$. With this field and $`a(0)`$, $`𝒗(0)`$ we can now apply the algorithm in the first half of the interval $`[0,t_{N1}]`$ to get $`p(0)`$. The computational cost is $`c_b+C_{N1}`$. Summing up all the costs we have the recursive relation $`C_N=2^{N1}c_f+2C_{N1}+c_b`$, with solution $`C_N=\frac{N}{2}2^Nc_f+(2^N1)c_b`$. Comparing $`C_N`$ with the cost of a forward integration of eqs. (1),(3) and (5), that is $`(c_f+c_b)(2^N1)`$, we have an increase by a factor $`\frac{c_f}{2(c_f+c_b)}N`$, i.e. logarithmic in the number of timesteps. As for memory requirements, only the $`N`$ couples of fields $`a(t_Nt_k),𝒗(t_Nt_k),k=1,\mathrm{},N`$ need to be stored – again a logarithmic factor.

A typical evolution of the propagator is shown in the central column of Fig. 3. We can now reconstruct the time sequence of the forcing contributions $`\varphi _{a,c}(s)=𝑑𝒚f_{a,c}(𝒚,s)p(𝒚,s|𝒙,t)`$ which, integrated over $`s`$, gives the amplitude of the scalar fields according to (6) and (7). As shown in Fig. 4, the time series of $`\varphi _a(s)`$ and $`\varphi _c(s)`$ are markedly different.

In the active scalar case, the sequence is strongly skewed towards positive values at all times. This signals that the trajectories preferentially select regions where $`f_a`$ has a positive sign, summing up forcing contributions to generate a typical variance of $`a`$ of the order $`F_0t`$. Conversely, $`f_c`$ can be positive or negative with equal probability on distant trajectories, and the time integral in eq. (7) averages out to zero for $`|st|>\tau `$ ($`\tau =l_f/v_{\mathrm{rms}}`$ is the eddy-turnover time). Therefore, typically $`c^2F_0\tau `$. The cumulative effect of the correlation between forcing and propagator is even more evident from the relation $`_0^s\varphi _a(s^{})𝑑s^{}=𝑑𝒚a(𝒚,s)p(𝒚,s|𝒙,t)`$ derived from (1) and (4). As shown in Fig. 4, the growth of the scalar variance is thus related to a strong spatial correlation between the propagator and the active field. That is further evidenced by comparing the first and the second column of Fig. 3: the distribution of particles closely follows the distribution of like-sign active scalar. This amounts to say that large-scale scalar structures are built out of smaller ones which “coalesce” together Bis2 . This situation has to be contrasted with the absence of large-scale correlations between the propagator and the passive scalar field (second and third column of Fig. 3).

As remarked previously, the absence of dissipative anomaly is closely related to the onset of an inverse cascade. Let us now discuss this issue from a Lagrangian viewpoint and consider the squared active field $`a^2`$. On one hand, it can be written as the square of (6). On the other hand, one can multiply (1) by $`2a`$ to obtain the equation $`_ta^2+𝒗a^2=\kappa \mathrm{\Delta }a^2+2af_a2ϵ_a`$ and solve it for $`ϵ_a=0`$ in terms of particle trajectories. The comparison of the two previous expressions yields

$`{\displaystyle _0^t}𝑑s{\displaystyle _0^t}𝑑s^{}{\displaystyle f_a(𝒚,s)f_a(𝒚^{},s^{})p(𝒚,s|𝒙,t)p(𝒚^{},s^{}|𝒙,t)}`$ $`=`$
$`{\displaystyle _0^t}𝑑s{\displaystyle _0^t}𝑑s^{}{\displaystyle f_a(𝒚,s)f_a(𝒚^{},s^{})p(𝒚^{},s^{};𝒚,s|𝒙,t)},`$ (8)
where $`p(𝒚^{},s^{};𝒚,s|𝒙,t)=p(𝒚,s|𝒙,t)p(𝒚^{},s^{}|𝒚,s)`$ denotes the probability that a trajectory ending in $`(𝒙,t)`$ were in $`(𝒚,s)`$ and $`(𝒚^{},s^{})`$. Integration over $`𝒚`$ and $`𝒚^{}`$ is implied. To obtain the r.h.s. of (8) consider the expression $`a^2(𝒙,t)=2_0^t𝑑s𝑑𝒚f_a(𝒚,s)a(𝒚,s)`$, insert $`a(𝒚,s)=_0^s𝑑s^{}𝑑𝒚^{}f_a(𝒚^{},s^{})`$, i.e. eq. (6) evaluated at time $`s`$, and exploit the symmetry under exchange of $`s`$ and $`s^{}`$ to switch from a time-ordered form to a time-symmetric one to get rid of the factor two.

Notice as a side remark that by integrating the r.h.s. of (8) over $`𝒙`$, exploiting the normalization $`𝑑𝒙p(𝒚,s|𝒙,t)=1`$ (as follows from(5)), and averaging over the forcing statistics, one obtains $`a^2=F(0)t`$. The result can be generalized to show that $`a^{2n}=(2n1)!!(F(0)t)^n`$, in agreement with the Eulerian argument for the Gaussianity of the single-point pdf of $`a`$ given above.

Eq. (8) can be recast as $`_0^tf_a(𝑿(s),s)𝑑s_X^2=[_0^tf_a(𝑿(s),s)𝑑s]^2_X`$, leading to the conclusion that $`_0^tf_a(𝑿(s),s)𝑑s`$ is a non-random variable over the ensemble of trajectories. This result has a simple interpretation: the absence of dissipative anomaly (i.e. $`ϵ_a=0`$) is equivalent to the property that along any of the infinite trajectories $`𝑿(s)`$ ending in $`(𝒙,t)`$ the quantity $`_0^tf_a(𝑿(s),s)𝑑s`$ is exactly the same, and equals $`a(𝒙,t)`$. Therefore, a single trajectory suffices to obtain the value of $`a(𝒙,t)`$ contrary to the passive case, where different trajectories contribute disparate values of $`_0^tf_c(𝑿(s),s)𝑑s`$, with a typical spread $`ϵ_ct`$. In this case, only the average over all trajectories yields the correct value of $`c(𝒙,t)`$.

The previous property extends to multiple points. Proceeding as for (8), one can show that $`m`$ trajectories are enough to calculate the product of arbitrary powers of $`a`$ at $`m`$ different points. In particular, structure functions $`\left(a(𝒓)a(\mathrm{𝟎})\right)^n`$ for any order $`n`$ can be shown to involve only two trajectories. This should be contrasted with the passive scalar, where the number of trajectories increases with $`n`$ and this is at the core of anomalous scaling of the field FGV01 . Indeed, we do not observe any anomalous scaling for the structure functions of the active field $`a`$ and the pdf’s of its increments are self-similar. Note that the same phenomenon and a lack of dissipative anomaly were also found for a passive scalar in a compressible flow GV . In that case, the ensemble of trajectories collapses onto a unique path, fulfilling in the simplest way the constraint (8) and its multi-point analogs. Here, even though the trajectories do not collapse, the constraints are satisfied due to the subtle correlation between forcing and trajectories peculiar to the active case.

We conclude with the following open issue. It is known that passive correlation functions are controlled by the asymptotically dominant terms in particle propagators (see, e.g., Ref. FGV01 ). Is this also true for active scalars ? The correlations between propagator and input in (6) suggest this be not the general rule. For some specific systems, those correlations might however be such that the asymptotics of the propagator still stands out at large times Itamar ; Matsu . Clarifying the dynamical conditions controlling this phenomenon will be the subject of future work.

We are grateful to A. Noullez and I. Procaccia for useful discussions. This work has been supported by the EU under the contract HPRN-CT-2000-00162, and by Indo-French Center for Promotion of Advanced Research (IFCPAR 2404-2). AM has been partially supported by Cofin2001 (prot.2001023842). Numerical simulations have been performed at IDRIS (project 021226) and at CINECA (INFM parallel computing initiative).


text:
# Cytoskeleton confinement of red blood cell membrane fluctuations

## Abstract

We analyze both the static and dynamic fluctuation spectrum of the red-blood cell in a unified manner, using a simple model of the composite membrane. In this model, the two-dimensional spectrin network that forms the cytoskeleton is treated as a rigid shell which is located at some constant average separation from the lipid bilayer. The cytoskeleton thereby confines both the static and dynamic fluctuations of the lipid bilayer. The predictions of the model account for the wavevector and frequency dependence of the experimental data. The observed amplitude of the thermal fluctuations is related to effects of ATP-driven fluctuations.

A long-standing problem in the study of red blood cell (RBC) structure is the simultaneous softness of its membrane observed by thermal fluctuations sackman95 , and the strong shear elasticity found in static deformation experiments, such as micropipette aspiration pipette  and electrodeformation sackman84 . The membrane itself is a composite structure sackman94  with an outer, gel-like extracellular network of long sugar molecules (thought to be irrelevant to the structural strength), a mixed lipid/protein bilayer and an attached, intracellular network. Previous theoretical models of this membrane treated it as a single, polymerized network with the combined curvature bending modulus of the lipid bilayer $`\kappa `$ and the shear rigidity of the cytoskeleton $`\mu `$ peterson92 . Such models were successful in describing the response of the membrane in static deformation experiments, which give $`\mu 10^510^6`$J/m2 pipette ; sackman84 ; discher01 . However, comparing these models to the fluctuation data, leads to the conclusion that the membrane behaves as if the shear modulus vanishes $`\mu 0`$ brochard ; peterson92 . This surprising conclusion comes from the shape and amplitude analysis of the longest wavelength shape fluctuations sackman92 ; sackman95 . Various ideas have been raised in order to account for this observation, the main suggestion being that ATP-driven structural rearrangement in the spectrin network discher01 ; sackman94 ; sackman95  relaxes the shear-like deformations.

The previous studies were concerned with the shape fluctuations of largest wavelength. Here, we focus on the fluctuation spectrum, at the length-scales of $`10.01\mu `$m, where the effect of the cytoskeleton is observed sackman87 . The important question is to what extent the cytoskeleton effects are distinguishable from the fluctuations of a free, closed bilayer. We show that for a consistent description of both the static and the dynamic fluctuation spectrum we must include the confining effects of the cytoskeleton. A simple model where the only significant effect of the spectrin network is to confine the lipid bilayer membrane, consistently describes both the spatial and temporal spectra of the thermal fluctuations of the RBC membrane. We discuss how active processes that effectively increase the temperature of the membrane, can be accounted for within this description.

The curvature bending modulus $`\kappa `$ of the lipid bilayer seifert  is deduced from measurements of the amplitude of thermal fluctuations at the smallest measured wavelengths sackman87  $`\kappa 2\pm 1\times 10^{20}`$J. The RBC cytoskeleton is a two-dimensional, roughly triangular, network of spectrin proteins bennett , that is attached to the lipid bilayer at the nodes and at additional, random sites along the spectrin polymers. The cytoskeleton is well described as a network of entropic springs, of length $`l80`$nm bennett ; sackmangreenbook , with an effective spring constant: $`\mu 4\times 10^6`$J/m2, which is close to the measured static shear modulus: $`6\times 10^6`$J/m2 pipette . Compared with the lipid bilayer bending modulus $`\kappa `$, the curvature bending modulus of the cytoskeleton is negligible sackman95 ; landaumcdonald01discher98 : $`\kappa _{cyto}\mu w^210^{21}10^{22}`$J for a cytoskeleton thickness of $`w100500`$Å.

We now begin by analyzing the measured spatial correlations of the fluctuations sackman87 , and describe the effects of the cytoskeleton on the bilayer in terms of continuum mechanics. This is feasible since the cytoskeleton forms a rather open mesh that is attached to the bilayer at discrete points with a contact area that is small ($`1`$nm) compared to the inter-node distance ($`100`$nm). In a coarse grained picture, we can describe the thermal fluctuations of the bilayer using a continuum model of the free energy functional

$$F𝑑S\left[\frac{1}{2}\sigma \left(h\right)^2+\frac{1}{2}\kappa \left(^2h\right)^2+\frac{1}{2}\gamma h^2\right]$$
(1)
that includes the bending and effective surface tension energy of the bilayer (see for example safran ) in terms of the normal displacement $`h`$. The surface tension coefficient $`\sigma `$ arises from the constraint of constant surface area of the bilayer in the closed geometry of the RBC sackman97 ; seifert : $`\sigma _0=\alpha \kappa /R^21\times 10^9`$J/m2 (taking $`R4\mu `$m as the RBC radius, $`\alpha =1`$). The constant can vary $`\alpha 1100`$.

The last term in (1) is mathematically equivalent to a Lagrange multiplier that constrains the mean square amplitude of bilayer fluctuations to be equal to: $`d^2=k_BT/8(\gamma \kappa )^{1/2}`$, for an infinite bilayer. This term describes the effect of the cytoskeleton on the bilayer through a harmonic potential that maintains an average separation $`d`$ (of order $`w`$) between the lipid bilayer and the cytoskeleton safran , here treated as a separated, infinitely rigid shell, that does not participate in the thermal fluctuations. The discrete contacts that maintain the constant average separation are not specifically described in this continuum model; in a coarse-grained picture, these contacts are the physical origin of the constraint (potential) that determines the average membrane-spectrin network separation.

The attachment of the cytoskeleton to the bilayer also causes stretching and undulations of the bilayer sackman90 , partly due to steric repulsion between the spectrin and the bilayer around the point of attachment lipowsky . Balancing the spectrin stretching with a local curvature of the bilayer, results in a membrane with undulations of wavelength sackman90 ; steric  $`L\sqrt{\kappa /\mu }100200`$nm and amplitude $`10`$nm($`w`$) sackman95 ; landaumcdonald01discher98 . In our confining-shell model this length-scale is related to the potential-induced persistence length of the bilayer $`\xi _0=(\kappa /\gamma )^{1/4}L`$ safran , i.e. the wavelength below which the bilayer is freely fluctuating.

We now calculate the spatial correlations for a two-dimensional, flat bilayer, since for all but the largest wavelengths $`\lambda `$, the surface of the RBC is relatively flat: $`50nm<\lambda <1\mu m<R4\mu `$m. We determine the values of $`\sigma `$ and $`\gamma `$ by fitting to the experimental data. From Eq. (1) the equal-time (static) correlations of the normal deflections of the bilayer can be written safran ; sackman92

$$h_qh_q=\frac{k_BT}{\kappa _qq^4},\kappa _q=\kappa +\sigma q^2+\gamma q^4$$
(2)
In the inset of Fig.1 we plot the measured value of $`\kappa _q`$ sackman87  in the form $`(\kappa _q/\kappa 1)^1`$ as a function of the normalized wavevector $`(qd)^4`$ (where $`d`$ is determined by fitting the data to obtain $`\gamma d^4`$). From the linear slope in the limit of $`q0`$ we find the values of the parameter $`\gamma =7.5,1.0\times 10^7`$J/m4, for the two cells measured. These values correspond to mean amplitudes $`d200,350`$Å and $`\xi _0=(\kappa /\gamma )^{1/4}=130,220`$nm respectively. At larger values of $`q`$ there is a noticeable deviation from a straight line, which arises from the effective surface tension $`\sigma 7.7,2.8\times 10^7`$J/m2 for the two cells. Note that surface tension alone, without the confining effect of the cytoskeleton (i.e. $`\gamma =0`$), does not fit the data (dash-dot line, Fig.1 inset). In addition, there is a rather abrupt change at the crossover wavevector $`q_0=1/\xi _0`$ (indicated by the vertical dashed lines in Fig.1), above which the data are better described using $`\sigma \sigma _01.4\times 10^9`$J/m2 (solid lines in Fig.1).

The measured surface tension is consistent with: $`\sigma \kappa /\xi _0^2`$ (of order $`\mu /10`$). This expression gives the effect of bilayer shape constraint due to the static undulations of lateral size $`\xi _0\times \xi _0`$, described above. Indeed at length-scales shorter then $`\xi _0`$, there is no stretching of the cytoskeleton (Fig.1). Note that the effective surface tension of a closed bilayer is a very sensitive function of the excess area of the bilayer seifert ; sens , which is affected by the induced undulations. The spread in the measured parameters may be due to natural variations in the cytoskeleton network of normal RBC cells.

There is a qualitative difference in power-law of the wavevector dependence of $`\kappa _q`$ for RBC and empty vesicles sackman97 . The vesicles are well described (Fig.2) by equations (1,2) with $`\gamma =0`$, and an effective surface tension: $`\sigma _{vesicle}\kappa /R^22\times 10^{10}`$J/m2, where $`R27\mu `$m and $`\kappa 1.3\times 10^{19}`$J are the vesicle radius and bending modulus respectively. Both the RBC’s and vesicles data collapse when the wavelength is scaled by the r.m.s amplitude $`d`$ (Fig.2). For the vesicles of diameter $`50\mu `$m sackman97  the r.m.s. thermal amplitudes are $`d11.5\mu `$m (note that here $`d`$ is not related to confinement). The good scaling of the data indicates that there is indeed a single important length-scale in the problem, namely the persistence length $`\xi _0`$, that determines all the parameters appearing in the free energy ($`\gamma `$ and $`\sigma `$) and the r.m.s amplitude $`d`$.

We now use the same simple model of spectrin confinement of the bilayer to describe the temporal correlations of the membrane fluctuations. The shape fluctuations of the RBC membrane are driven by both thermal and metabolic energies. The active fluctuations have a frequency spectrum that is confined to the range 0.3-1Hz atprafi . For higher frequencies, our analysis shows that the active processes can be accounted for by an increase in the effective temperature of the fluctuations atprafi ; prost . The temporal height-height correlation function brochard ; anton  for a flat bilayer at a distance $`D_\omega `$ from a rigid wall, is

$$h_q(t)h_q(0)=\frac{k_BT}{\kappa _qq^4}e^{\omega (q)t}$$
(3)
where $`\kappa _q`$ is given in (2). The hydrodynamic interaction (Oseen interaction kernel landauhyd ) has a modified form: $`\mathrm{\Lambda }(q)=(1/4\eta q)\left[1(1+2qD_\omega )e^{2qD_\omega }\right]`$ ($`\mathrm{\Lambda }_f(q)=1/4\eta q`$ for a free bilayer), so that the relaxation frequency $`\omega (q)`$ is

$$\omega (q)=\left[\frac{1}{4\eta }\left(\kappa q^3+\sigma q+\frac{\gamma }{q}\right)\right]\left[1(1+2qD_\omega )e^{2qD_\omega }\right]$$
(4)
where $`\eta 3\eta _{water}`$ is some average viscosity of the cytoplasm and external solution. In the limit of short wavelengths ($`q\mathrm{}`$) we recover the free bilayer frequency: $`\omega (q)\kappa q^3/4\eta `$.

The mean square amplitude of the normal fluctuations, as a function of frequency $`\omega `$, is the Fourier transform

$`d(\omega )^2`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle q𝑑qh_q(t)h_q(0)e^{i\omega t}𝑑t}`$ (5)
$`={\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle h_q(0)h_q(0)\frac{\omega (q)}{\omega (q)^2+\omega ^2}q𝑑q}`$
For a free bilayer this expression (5) gives an anomalous frequency dependence anton : $`d(\omega )^2\omega ^{5/3}`$. We integrate the expression (5) numerically in the range $`\pi /R<q<\pi /a`$ ($`a50`$Å), and compare with the experimental data atprafi . In the inset of Fig.3 we plot $`d(\omega )`$ using the parameters of the two cells of Fig.1, and both the high and low values of the effective surface tension ($`\sigma `$ and $`\sigma _0`$ respectively). We find a reasonable agreement between the calculation and the measurements, when taking $`D_\omega 0.4d`$. The similar magnitude of the bilayer-rigid shell separation from both static and dynamic experiments, shows the overall consistency of our confinement model.

In the limit of high frequencies, the earlier result of Brochard et.al. brochard , gave $`d(\omega )^2\omega ^{4/3}`$. Since $`qD_\omega 1`$ in the measured range, we find in this limit (4): $`\omega (q)\kappa q^5D_\omega ^2/4\eta `$, leading to $`d(\omega )^2\omega ^{7/5}`$. It is difficult to distinguish between these two values using the newer data atprafi . Our calculation has the advantage of consistently describing both the static (spatial) and dynamic (temporal) fluctuation data. Note that the case of a pure bilayer with large effective surface tension $`\sigma `$, but without the effect of the rigid wall, is in complete disagreement with the data (dotted line, Fig.3 inset).

In Fig.3 we show that the normal RBC, ATP depleted RBC and RBC ghost, are all well described by the same expression (5) (using the softer cell from Fig.1, i.e. the smaller $`\gamma `$ and $`\sigma `$), differing by an effective temperature factor of up to $`3`$. This is similar to the amplitude enhancement factor of $`2.5`$ found in a previous study prost .

The largest effect of the rigid shell is to increase the effective viscosity of the water near the bilayer, by constraining its flow. Defining an effective viscosity by: $`\omega =\kappa q^3/4\eta _{eff}`$, we get from (4): $`\frac{\eta _{eff}}{\eta }=\frac{\kappa q^2}{4d^2\left(\gamma +\sigma q^2\right)}`$. At the crossover wavevector $`q_0`$ this function has its peak: $`\eta _{eff}/\eta 4530`$, depending on the value of $`\sigma `$. These values are in close agreement with the value $`\eta _{eff}/\eta 50`$ found from the relaxation times of an electrodeformed RBC peterson92 ; sackman84 . In these experiments, the cytoplasm flows through the cytoskeleton mesh, setting up a flow field at the crossover wavevector $`q_0`$. Thus, a rigid cytoskeletal wall, separated by a fixed distance from the bilayer, accounts for the larger effective viscosity required to fit these dynamical experiments.

While this model accounts for the wavevector dependence of the statics and the frequency dependence of the dynamics, the absolute amplitude of the fluctuations and the different values observed in active and ATP-depleted cells, must still be explained. One possibility is that ATP driven fluctuations completely determine the amplitude of the largest wavelength fluctuations atprafi  through the process of spectrin-actin disconnections and reconnections atprafi ; sackman95 ; bennett . These ATP-driven conformational changes give rise to defects in the triangular spectrin network, resulting in nodes with more or less than 6 spectrin bonds. The local curvature of the cytoskeleton may change at the site of a defect, from being locally flat (6 bonds) to having a $`80`$nm deviation out of the plane of the flat cytoskeleton (5-fold node). The effect of this random buckling is to increase the mean bilayer-rigid shell separation by a factor of $`4`$. According to our model, this will increase the amplitude of the $`q0`$ modes by a factor of $`4^4`$, as measured peterson92 .

###### Acknowledgements.

We thank R. Korenstein, E. Sackmann and H. Strey for useful discussions. This work was supported by the ISF grant for Center on Self-Assembly. The authors are grateful to the donors of the Petroleum Research Fund administered by the American Chemical Society and to the Schmidt Minerva Center for their support. N.G.’s research is being supported by the Louis L. and Anita M. Perlman Postdoctoral Fellowship.


text:
# 1 Introduction

## 1 Introduction

Thermal neutrons are ideally suited to probing the excitations in condensed-matter systems, as their wavelengths $`0.5<\lambda <10`$ Å are comparable to the interatomic distances, and their energies are of the same magnitude as the excitations. In particular, the small-angle neutron scattering (SANS) gives valuable information on thermal fluctuations in liquids.

The only previous SANS experiment on liquid helium has been performed by Egelstaff and London long ago, in 1953 . Unfortunately, there are only a few experimental points below the lambda point, and it is difficult to come to a definite conclusion about the small-angle scattering behaviour in a wide temperature region. Modern high-flux reactors and advanced instruments have increased the detection sensitivity by orders of magnitude. Additionally, previous detail experiments on the static structure factor of helium-4 have been performed only for the momentum transferred above 1 Å-1 .

This was a motivation for us to remeasure SANS from liquid helium in a wide temperature range.

## 2 Experiment

The experiments were carried out at the instrument D22 of the Institut Laue-Langevin (ILL), Grenoble, France. We used a standard ILL Orange Cryostat modified for SANS with quartz windows. We filled the inner part of the cryostat, which is delimited by transparent quartz windows at the neutron entrance and exit sides, with liquid helium by condensing it into the sample volume. We then lowered the helium temperature first by opening the cold valve, and when we came close to the phase transition to He-II, by directly pumping on the helium.

The experiments were performed with 4.6 Å neutrons. The sample-to-detector distance was 2.8 m (just avoiding shadowing on the detector that was centred with respect to the beam). The ”collimation” distance, i.e. the distance between the end of the neutron guide of 40 mm (wide) by 55 mm (high) and the sample, was 5.6 m. The sample aperture had a diameter of 10 mm. The scattering from the empty cryostat at low temperature, and the neutron and electronic noise were measured in order to correct the data for backgrounds. The data were radially averaged and the container and noise contributions were subtracted by standard programs; they were put on an absolute scale by normalizing with the scattering from 1 mm of water. We also checked by measurements with a detector setting of 14 m that there was no central increase in the helium scattering pattern.

## 3 Results and Discussion

Figure 1 shows radially averaged scattered neutron intensities versus momentum transfer for some helium temperatures and water normalized scattering patterns with sample background subtracted.

The momentum transferred is calculated on the assumption that the neutron scattering is elastic and is related to the full scattering angle $`\theta `$ through the expression

$$|Q|=\frac{4\pi }{\lambda }\mathrm{sin}(\theta /2),$$
$`(1)`$
where $`\lambda `$ is the wavelength of the neutron.

Although the background from the cryostat is rather small in comparison with the intensity scattered from helium, we see an increase of the background intensity at low momentum transfer ($`Q<0.03`$Å-1) and, as a consequence, a decrease of an intensuty of scattered from helium neutrons as seen in Fig.1. A reasonable explanation is that the small-$`Q`$ scattering from the cryostat is a result of the existence of dust, cracks, scratches, that are small in size, on the surfaces of the quartz windows. Special experiments are needed to clarify this situation, and that is why we analyse only experimental data in the $`Q`$-range $`(0.03÷0.3)`$ Å-1.

The essential quantity in the description of a liquid is its structure factor $`S(Q)`$ that is a measure of the correlation between positions of the atoms in fluid. The macroscopic cross section $`d\mathrm{\Sigma }/d\mathrm{\Omega }`$ of neutron scattering is related to the static structure factor $`S(Q)`$ by

$$\frac{d\mathrm{\Sigma }}{d\mathrm{\Omega }}(Q)=\frac{\sigma _b}{4\pi }\rho _{\mathrm{at}}S(Q),$$
$`(2)`$
where $`\sigma _b`$ is the bound atom scattering cross section of helium (= 1.172 barn) and $`\rho _{\mathrm{at}}`$ is the atomic density of helium.

According to Goldstein , at finite temperature the structure factor is given by

$$S(Q)=n_0k_\mathrm{B}T\chi _T+\underset{n=1}{\overset{\mathrm{}}{}}(1)^nQ^{2n}r_g^{2n}[(2n+1)!]^1,$$
$`(3)`$
where $`n_0`$ is the number density, $`k_\mathrm{B}`$ the Boltzmann constant, $`\chi _T`$ the isothermal compressibility, and $`r_g^{2n}`$ is the moment of the pair correlation function $`g(r)`$ defined by

$$r_g^{2n}=n_0r^{2n}[1g(r)]d^3r.$$
$`(4)`$
In particular, to order $`Q^2`$ Eq.(3) is just

$$S(Q)n_0k_BT\chi _Tr_g^2\frac{Q^2}{6},S(0)=n_0k_BT\chi _T.$$
$`(5)`$
This approximation looks like the Guinier approximation of the decay of the intensity near $`Q=0`$

$$I(Q)=I(0)\mathrm{exp}(Q^2R_g^2/3)I(0)(1Q^2R_g^2/3),$$
$`(6)`$
where $`R_g`$ is the radius of gyration.

On the other hand, as it shown by Ornstein and Zernike, in the range of critical point the intensity of scatterred neutrons is described as

$$S(Q)=\frac{S(0)}{1+r_c^2Q^2}S(0)(1r_c^2Q^2).$$
$`(7)`$
We see that there is a direct relation between $`r_c`$, $`R_g`$ and $`r_g`$. The critical temperature of helium equals 5.2 K, and thus we can definitely say that at temperature around 4 K the value of $`R_g`$, as well as $`r_g`$, reflects correlation between helium atoms and we can consider them as correlation radius. As to smaller temperatures, the physical meaning of these quantities not so apparent.

The parabolic behaviour of the intensity of scattered neutrons at small $`Q`$ is clearly seen in the experiment (Fig.2).

It is necessary here to note that we have not take into account the scattering from zero-point fluctuations that leads to the linear function of $`Q`$ that has been first obtained by Feynman , but this process takes place at lower temperatures.

We calculated the zero-angle cross section by extrapolation of the scattered intensity to $`Q`$=0, and the results are shown in Fig.3. It is seen that the zero angle cross section is monotonically decreasing with decreasing temperature and does not show any peculiarity at the lambda-point. In the same figure we show the theoretical zero angle cross section $`S(0)`$ calculated using the known data for the temperature dependence of the helium density and the isothermal compressibility . At temperatures below the lambda-point both density and isothermal compressibility are practically independent of the temperature, and thus the zero-angle cross section is proportional to the temperature that allows it to fit linearly at low temperatures. As seen, there is a temperature dependent background, but we don’t know its origin.

We have measured SANS with a temperature interval 0.1-0.2 K that has provided the first comprehensive data on the correlation radius in the normal and superfluid state of helium. The results are shown in Fig. 4. The presented in Fig.1 data were corrected for constant background that corresponds to $`d\mathrm{\Sigma }/d\mathrm{\Omega }(0)`$ at zero temperature. It is seen that the temperature dependence of the correlation radius behaves like for a normal liquid, i.e. it decreases with decreasing temperature. However, at the superfluid transition temperature (2.2 K) we see a pronounced change of its slope.

At a qualitative level the observed loss of spatial correlations in liquid helium is due to a rearrangement of atoms to create the necessary spaces for delocalisation to occur. Nevertheless, this intuitive explanation has not theoretical basis. For instance, using Monte Carlo method Mayers  argued that the link between spatial correlations and Bose-Einstein condensate fraction is a geometrical consequence of the hard-core repulsion between atoms. The excitation of rotons has also been proposed by Masserini et al.  as an alternative explanation for the increase in spatial correlation as the temperature is rased in the superfluid state.

In this connection, we have to emphasize that in our SANS experiment the distance scale probed is greater than interatomic distances, and the correlation length obtained characterizes long-range order in liquid helium whereas in all previous experimental and theoretical works the behaviour of the first peak in the $`S(Q)`$ was considered, i.e. short-range correlations at $`Q`$-value about 2 angstrom-1. In our SANS experiment we deal with $`Q<0.3`$Å-1. As opposed to previous experiments we have studied long-range correlation between helium atoms.

## 4 Acknoledgements

The authors are deeply grateful A.S.Stepanov and S.M.Apenko for fruitful discussions of this work.


text:
# Chiral critical behavior in two dimensions from five-loop renormalization-group expansions

## I Introduction

The critical behavior of two-dimensional frustrated spin models with noncollinear or canted order has been the object of intensive theoretical and experimental studies, being of prime interest for those investigating layered magnetic systems with special structures and superconducting Josephson-junction arrays in an external magnetic field.

In physical magnets, frustration (leading to noncollinear order) may arise either because of the special geometry of the lattice or from the competition of different kinds of interactions. An example of the first kind is the two-dimensional triangular antiferromagnet, whereas for the latter case typical models are the fully frustrated $`XY`$ model on a square lattice and models with different nearest-neighbor and next-nearest-neighbor interactions. In these systems the Hamiltonian is minimized by the noncollinear configurations with a 120o spin structure. As a consequence, at criticality, there is a breakdown of O($`N`$) symmetry in the high-temperature phase to O($`N2`$) symmetry in the ordered phase.

Field-theoretical (FT) studies of systems with noncollinear order are based on the O($`N`$)$`\times `$O($`M`$) symmetric Landau-Ginzburg-Wilson (LGW) Hamiltonian

$`={\displaystyle d^dx}`$ $`\left\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{a}{}}\left[(_\mu \varphi _a)^2+r\varphi _a^2\right]+{\displaystyle \frac{1}{4!}}u_0\left({\displaystyle \underset{a}{}}\varphi _a^2\right)^2+{\displaystyle \frac{1}{4!}}v_0{\displaystyle \underset{a,b}{}}\left[(\varphi _a\varphi _b)^2\varphi _a^2\varphi _b^2\right]\right\},`$ (1)
where $`\varphi _a`$ ($`1aM`$) are $`M`$ sets of $`N`$-component vectors. We will consider the case $`M=2`$, that, for $`v_0>0`$, describes frustrated systems with noncollinear ordering. Negative values of $`v_0`$ correspond to simple ferromagnetic or antiferromagnetic ordering, and to magnets with sinusoidal spin structures .

The physically relevant cases are described by the Hamiltonian (1) with $`N=2`$ (frustrated $`XY`$ model) and $`N=3`$ (frustrated Heisenberg model). Despite the very intensive theoretical and experimental efforts to fully understand the critical behavior of these models, the situation is still controversial. For the frustrated $`XY`$ model the strongly debated point is either the single critical temperature $`T_c`$ exists at which both the $`SO(2)`$ and the $`Z_2`$ symmetries are simultaneously broken or there are two successive phase transitions at different critical temperatures. In the latter case, the order at which two phase transitions occur and the numerical values of the critical exponents are doubtful too. Obviously, studying the Hamiltonian (1), we cannot clarify the point just mentioned, since this Hamiltonian describes only the chiral phase transition. In more detail, these and other relevant issues are reviewed in Ref. , where one can find a complete list of references; we mention here only the very recent studies  (not quoted in Ref. ).

The critical behavior of the frustrated Heisenberg antiferromagnet in two dimensions is clearer: it undergoes one phase transition mediated by $`Z_2`$ topological defects. Again, a complete list of references can be found in Ref. . Here the debated point concerns whether the LGW Hamiltonian is able to keep the non trivial topological excitations present in this model, to the contrary of other approaches like the nonlinear $`\sigma `$ model ($`NL\sigma `$). In Ref.  it is claimed that it may be the case and we will assume it in the following.

The LGW Hamiltonian (1) has been extensively studied in the framework of $`ϵ=4d`$ , and by means of the $`1/N`$  expansion. The existence and the stability properties of the fixed points were found to depend on $`N`$ and on the spatial dimensionality . Within the $`ϵ`$ expansion, for sufficiently large values of $`N`$ the renormalization group equations possess four fixed points: the Gaussian fixed point ($`v=u=0`$), the $`O(2N)`$ ($`v=0`$) one and two anisotropic fixed points located in the region $`u,v>0`$ and usually named chiral and antichiral. The chiral fixed point is the only stable one. There is a critical dimensionality $`N_c`$ where the chiral and antichiral fixed points coalesce and disappear for $`N<N_c`$. In the last case, under the absence of any stable fixed point the system undergoes a weak first-order phase transition, since the associated RG flows run away from the region of stability of the fourth-order form in the free energy expansion. The three-loop estimate of $`N_c`$ is

$$N_c=21.823.4ϵ+7.1ϵ^2+O(ϵ^3),$$
(2)
that, after an appropriate resummation, results in $`N_c>3`$ in three dimensions. This inequality leads to the conclusion that for the physical models with $`N=2,3`$ the three-dimensional chiral transition is first order, as corroborated also by some other RG studies . To the contrary, both for $`N=2`$ and $`N=3`$ the highest-order six-loop calculations in three dimensions  reveal a strong evidence for a stable chiral fixed point that, however, is not related to its counterpart found within the $`ϵ`$ expansion. The four-loop two-dimensional analysis shows an equivalent topology of the fixed points location in the renormalized coupling constants plane . Furthermore, in a recent work  we claimed that both in two and three dimensions the controversial situation may reflect the quite unusual mode of critical behavior of the $`N`$-vector chiral model under the physical values of $`N`$. It was shown that this critical behavior is governed by a stable fixed point which is a focus, attracting RG trajectories in a spiral-like manner. Approaching the fixed point in a nonmonotonic way, looking somewhat irregular, may result in a large variety of the crossover and near-critical regimes.

In this paper, we report the perturbative five-loop two-dimensional RG expansions and carry out their analysis. Some of the results for $`N=2`$ and $`N=3`$ was already anticipated in . We point out the idea that the critical behavior for large $`N`$ is the one predicted by the $`ϵ`$ expansion. With decreasing $`N`$, the stable chiral fixed point becomes focuslike at some marginal value of $`N`$, $`N^{}`$ that differs, in principle, from $`N_c`$ where both chiral and antichiral fixed points disappear. With the precision of our calculations, we are not able to fully clarify if there is a region of values of $`N`$, where the chiral transition is first order. It is worthy to note that the parallel work on the three-dimensional models  reveals a full analogy between chiral critical behaviors in two and three dimensions.

The paper is organized as follows. In Sec. II we derive the perturbative series for the renormalization-group functions up to five-loop order. The resummations methods and the results of the analysis are presented in Sec. III. In Sec. IV we draw our conclusions.

## II PERTURBATIVE EXPANSIONS IN TWO DIMENSIONS

### A Renormalization of the theory

The fixed-dimension field-theoretical approach  represents an effective procedure in the study of the critical properties of systems that undergo second-order phase transitions . In the case under consideration, the expansion is performed in two quartic coupling constants of the Hamiltonian (1). The theory is renormalized by a set of zero-momentum conditions for the one-particle irreducible two-point correlation function, four-point correlation function, and two-point correlation function with an insertion of the operator $`\frac{1}{2}\varphi ^2`$:

$$\mathrm{\Gamma }_{ai,bj}^{(2)}(p)=\delta _{a,b}\delta _{i,j}Z_\varphi ^1\left[m^2+p^2+O(p^4)\right],$$
(3)
$$\mathrm{\Gamma }_{abcd}^{(4)}(0)=Z_\varphi ^2m\left[\frac{u}{3}\left(\delta _{ab}\delta _{cd}+\delta _{ac}\delta _{bd}+\delta _{ad}\delta _{bc}\right)+vC_{ai,bj,ck,dl}\right],$$
(4)
$$\mathrm{\Gamma }_{ab}^{(1,2)}(0)=\delta _{ab}Z_t^1,$$
(5)
where $`C_{ai,bj,ck,dl}`$ is an appropriate combinatorial factor .

The perturbative knowledge of the functions $`\mathrm{\Gamma }^{(2)}`$, $`\mathrm{\Gamma }^{(4)}`$, and $`\mathrm{\Gamma }^{(1,2)}`$ allows one to relate the renormalized physical parameters ($`u,v,m`$) to the bare ones ($`u_0,v_0,r`$ ).

The fixed points of the model are defined by the common zeros of the $`\beta `$ functions,

$`\beta _u(u,v)=m{\displaystyle \frac{u}{m}}|_{u_0,v_0},\beta _v(u,v)=m{\displaystyle \frac{v}{m}}|_{u_0,v_0}.`$ (6)
The stability properties of these points are determined by the eigenvalues $`\omega _i`$ of the matrix:

$$\mathrm{\Omega }=\left(\begin{array}{cc}\frac{\beta _u(u,v)}{u}& \frac{\beta _u(u,v)}{v}\\ \frac{\beta _v(u,v)}{u}& \frac{\beta _v(u,v)}{v}\end{array}\right).$$
(7)
A fixed point is stable if both the eigenvalues have a positive real part. If the eigenvalues possess nonvanishing imaginary parts, the fixed point is called a focus and the corresponding RG trajectories are spirals. The eigenvalues $`\omega _i`$ are connected to the leading scaling corrections, which go as $`\xi ^{\omega _i}|t|^{\mathrm{\Delta }_i}=|t|^{\nu \omega _i}`$, with few exceptions as the two-dimensional Ising model (see for a discussion Ref. ).

The values of the critical exponents $`\eta `$, $`\nu `$ and $`\gamma `$ are related to the RG functions $`\eta _\varphi `$ and $`\eta _t`$ evaluated at the stable fixed point:

$`\eta =\eta _\varphi (u^{},v^{}),`$ (8)
$`\nu =\left[2\eta _\varphi (u^{},v^{})+\eta _t(u^{},v^{})\right]^1,`$ (9)
where

$`\eta _\varphi (u,v)`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}Z_\varphi }{\mathrm{ln}m}}|_{u_0,v_0}=\beta _u{\displaystyle \frac{\mathrm{ln}Z_\varphi }{u}}+\beta _v{\displaystyle \frac{\mathrm{ln}Z_\varphi }{v}},`$ (10)
$`\eta _t(u,v)`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}Z_t}{\mathrm{ln}m}}|_{u_0,v_0}=\beta _u{\displaystyle \frac{\mathrm{ln}Z_t}{u}}+\beta _v{\displaystyle \frac{\mathrm{ln}Z_t}{v}}.`$ (11)
### B Five-loop series

In this section we extend the four-loop perturbative series  for the RG functions $`\beta _u`$, $`\beta _v`$, $`\eta _\varphi `$, and $`\eta _t`$ up to five-loop order using the numerical values of two-dimensional five-loop integrals evaluated in Ref. .

In order to obtain finite fixed point values in the limit of infinite components of the order parameter ($`N\mathrm{}`$) we use the rescaled couplings

$$u\frac{8\pi }{3}R_{2N}\overline{u},v\frac{8\pi }{3}R_{2N}\overline{v},$$
(12)
where $`R_N=9/(8+N)`$.

The resulting RG series are :

$`\overline{\beta }_{\overline{u}}=`$ $`\overline{u}+\overline{u}^2+{\displaystyle \frac{1N}{(4+N)}}\overline{u}\overline{v}{\displaystyle \frac{1N}{(8+2N)}}\overline{v}^2+{\displaystyle \underset{i+j3}{}}b_{ij}^{(u)}\overline{u}^i\overline{v}^j,`$ (13)
$`\overline{\beta }_{\overline{v}}=`$ $`\overline{v}{\displaystyle \frac{1}{2}}{\displaystyle \frac{6N}{(4+N)}}\overline{v}^2+{\displaystyle \frac{6}{(4+N)}}\overline{u}\overline{v}+{\displaystyle \underset{i+j3}{}}b_{ij}^{(v)}\overline{u}^i\overline{v}^j,`$ (14)
$`\eta _\varphi =`$ $`{\displaystyle \frac{1.83417(1+N)}{(8+2N)^2}}\overline{u}^2+{\displaystyle \frac{1.83417(1N)}{(8+2N)^2}}\overline{u}\overline{v}{\displaystyle \frac{1.37563(1N)}{(8+2N)^2}}\overline{v}^2+{\displaystyle \underset{i+j3}{}}e_{ij}^{(\varphi )}\overline{u}^i\overline{v}^j,`$ (15)
$`\eta _t=`$ $`{\displaystyle \frac{2(1+N)}{(4+N)}}\overline{u}{\displaystyle \frac{(1N)}{(4+N)}}\overline{v}+{\displaystyle \frac{13.5025(1+N)}{(8+2N)^2}}\overline{u}^2+{\displaystyle \frac{13.5025(1N)}{(8+2N)^2}}\overline{u}\overline{v}`$ (17)
$`{\displaystyle \frac{10.1269(1N)}{(8+2N)^2}}\overline{v}^2+{\displaystyle \underset{i+j3}{}}e_{ij}^{(t)}\overline{u}^i\overline{v}^j,`$
where

$$\overline{\beta }_{\overline{u}}=\frac{3}{16\pi }R_{2N}^1\beta _u,\overline{\beta }_{\overline{v}}=\frac{3}{16\pi }R_{2N}^1\beta _v.$$
(18)
Only the five-loop coefficients $`b_{ij}^{(u)}`$, $`b_{ij}^{(v)}`$, $`e_{ij}^{(\varphi )}`$ and $`e_{ij}^{(t)}`$ are reported in Table I, since up to the fourth order they are reported in Ref.. Note that, due to the last rescaling, the element of the $`\mathrm{\Omega }`$ matrix (7) are two times the derivative of $`\overline{\beta }`$ with respect to $`\overline{u}`$ and $`\overline{v}`$.

We verify the exactness of our series by checking them in particular limits ($`N=1`$, $`\overline{v}=0`$, $`N=M=2`$, $`N\mathrm{}`$) as in Ref. .

## III RESUMMATION AND ANALYSIS OF THE FIVE-LOOP SERIES

The field theoretical perturbative expansions are known to be divergent asymptotic series. Quantitative informations may be extracted from these series exploiting the properties of Borel summability for $`\varphi ^4`$ theories in two and three dimensions and resumming them by a Borel transformation combined with a method for the analytical extension of the Borel transform. This last procedure can be realized by using Padé approximants or via a conformal mapping, which maps the domain of analyticity of the Borel transform (cut at the instanton singularity) $`\overline{u}_b`$ onto a circle .

Let us consider a perturbative series in $`\overline{u}`$ and $`\overline{v}`$

$$R(\overline{u},\overline{v})=\underset{k=0}{\overset{l}{}}\underset{h=0}{\overset{lk}{}}R_{hk}\overline{u}^h\overline{v}^k,$$
(19)
which we want to resum (i.e., one of the $`\eta `$ or $`\beta `$ functions). In the Padé-Borel method we use the trick of the resolvent series considering the Borel transformed function,

$$\stackrel{~}{R}(\lambda ,\overline{u},\overline{v},b)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{\lambda ^n}{\mathrm{\Gamma }(n+b+1)}\underset{l=0}{\overset{n}{}}R_{l,nl}\overline{u}^l\overline{v}^{nl}.$$
(20)
Calling $`P_M^N(\lambda ,\overline{u},\overline{v},b)`$ the $`[N/M]`$ Padé approximant of the series $`\stackrel{~}{R}(\lambda ,\overline{u},\overline{v},b)`$ in the variable $`\lambda `$, one obtains the estimate of the desired quantity as

$$\stackrel{~}{P}_M^N(b,\overline{u},\overline{v})=_0^{\mathrm{}}e^tt^bP_M^N(t,\overline{u},\overline{v},b)𝑑t.$$
(21)
In this way we produce several approximants of the function $`R(\overline{u},\overline{v})`$ with varying the parameters $`b`$, $`N`$, and $`M`$.

In the conformal mapping method, the idea is to exploit, in the course of the resummation, the knowledge of the large order behavior of the series $`R(\overline{u},\overline{v})`$ at fixed $`z=\overline{v}/\overline{u}`$ . This large-order behavior is related to the singularity of the Borel transform closest to the origin . For the Hamiltonian (1), the value of the Borel transform singularity closest to the origin $`\overline{u}_b`$ has been computed in Ref , obtaining:

$`{\displaystyle \frac{1}{\overline{u}_b}}=`$ $`aR_{2N}`$ $`\mathrm{for}0<z<4,`$ (22)
$`{\displaystyle \frac{1}{\overline{u}_b}}=`$ $`aR_{2N}\left(1{\displaystyle \frac{1}{2}}z\right)`$ $`\mathrm{for}z<0,z>4,`$ (23)
where $`a=0.238659217\mathrm{}`$. For $`z>2`$ there is a singularity on the real positive axis which, however, is not the closest one to the origin for $`z<4`$. Thus, for $`z>2`$ the series are not Borel summable. We resum the RG functions also for the values of $`z`$ for which the series are not Borel summable since the method should provide a reasonable estimate if $`z<4`$, because we take into account the leading large-order behavior.

With the knowledge of $`\overline{u}_b`$, one can perform the mapping

$$y(\overline{u};z)=\frac{\sqrt{1\overline{u}/\overline{u}_b(z)}1}{\sqrt{1\overline{u}/\overline{u}_b(z)}+1},$$
(24)
in order to extend $`R(\overline{u},\overline{u}z)`$ to all positive values of $`\overline{u}`$. Consequently one obtains the set of approximants

$$E(R)_p(\alpha ,b;\overline{u},\overline{v})=\underset{k=0}{\overset{p}{}}B_k(\alpha ,b;\overline{v}/\overline{u})\times _0^{\mathrm{}}𝑑tt^be^t\frac{y(\overline{u}t;\overline{v}/\overline{u})^k}{[1y(\overline{u}t;\overline{v}/\overline{u})]^\alpha },$$
(25)
where the coefficient $`B_k`$ is determined by the condition that the expansion of $`E(R)_p(\alpha ,b;\overline{u},\overline{v})`$ in powers of $`\overline{u}`$ and $`\overline{v}`$ gives $`R(\overline{u},\overline{v})`$ to order $`p`$.

An important issue in the fixed dimension approach to critical phenomena concerns the analytic properties of the $`\beta `$ functions. As shown in Ref.  for the $`O(N)`$ model, the presence of confluent singularities in the zero of the perturbative $`\beta `$ function causes a slow convergence of the resummation of the perturbative series to the correct fixed-point value. The apparent stability of the results when analyzing a finite number of terms of the perturbative expansion may not provide a reliable indication of the uncertainty of the overall estimates. The $`O(N)`$ two-dimensional field-theory valuations of physical quantities  are less accurate than the three-dimensional counterparts, due to the bigger critical value of the quartic coupling constant and stronger nonanalyticities at the fixed point . In Ref.  it is explicitly shown that the nonanalytic terms cause a large deviation in the estimate of the right correction to the scaling, i.e., the exponent $`\omega `$. At the same time, the perturbative results for the fixed-point location turn out to be rather good approximations for the exact ones. We think that a similar situation holds for frustrated models.

### A Fixed points, stability, and critical exponents

In order to calculate the fixed points of the theory, we resum the perturbative expression for each $`\beta `$ function in the whole plane $`(\overline{u},\overline{v})`$ following the procedure explained above. With the conformal mapping method, we choose the approximants that stabilize the series for small values of the coupling constants with varying the considered perturbative order, i.e., $`\alpha =0,1,2`$ and $`b=5,7,\mathrm{},15`$. In this manner we have $`18`$ approximants for each $`\overline{\beta }`$ function. To obtain reasonable values and error bars for the fixed-point coordinates, we divide the domain $`0\overline{u}4`$, $`0\overline{v}6`$ in $`40^2`$ rectangles, and we mark all the sites in which at least two approximants for $`\overline{\beta }_{\overline{u}}`$ and $`\overline{\beta }_{\overline{v}}`$ vanish. This procedure is applied for the four-loop and five-loop series and the results for the zeros of the $`\beta `$ functions for $`N=2`$ and $`N=3`$ are displayed in Figs. 1 and 2, respectively. We remind the reader that the four-loop results were previously reported in Ref. , where similar figures were shown also for the three-loop case. It was also shown that no upper branch of zeros of $`\overline{\beta }_{\overline{u}}`$ appears in the three-loop approximation for all considered values of $`N`$.

As one can see from Figs. 1 and 2, for these physical systems the resummed RG expansions yield similar pictures for the zeros of the $`\overline{\beta }`$ functions. The existence of three fixed points – the gaussian, the $`O(2N)`$, and the chiral – is clear in both cases, while the presence of an unstable antichiral fixed point looks doubtful. In fact, if this last point exists, it should be located in the domain where the resummation procedure fails to give reliable quantitative results ($`z>4`$).

The location of the chiral fixed point for $`N=2`$ is

$`[u^{},v^{}]=`$ $`[2.3(2),5.0(5)]`$ $`(\text{4-loop CM}),`$ (26)
$`[2.25(25),4.80(45)]`$ $`(\text{5-loop CM}),`$ (27)
$`[2.35(30),4.35(25)]`$ $`(\text{Padé}),`$ (28)
and for $`N=3`$

$`[u^{},v^{}]=`$ $`[2.3(3),3.9(5)]`$ $`(\text{4-loop CM}),`$ (29)
$`[2.2(2),3.9(3)]`$ $`(\text{5-loop CM}),`$ (30)
$`[2.25(20),3.60(15)]`$ $`(\text{Padé}).`$ (31)
The stability properties of the chiral fixed point depend on the signs of the real parts of the eigenvalues $`\omega _i`$ of the $`\mathrm{\Omega }`$ matrix. We take the above 18 approximants for each $`\overline{\beta }`$ function and compute the numerical derivatives of each couple of approximants of the two $`\overline{\beta }`$ functions at their common zero, obtaining 324 possible combinations. To have reliable numerical estimates for $`\omega _i`$, we limit ourselves by the approximants that yield the chiral fixed point coordinates compatible with their final, properly weighted values.

Both for $`N=2`$ and $`N=3`$, we find that the chiral fixed point is a stable focus since it possesses stability eigenvalues with nonvanishing imaginary parts and positive real parts. At the five-loop (four-loop) approximations for $`N=2`$ the complex eigenvalues are produced by 87% (87%) of the working approximants within the conformal mapping method, while for $`N=3`$ by 89% (74%) of the approximants employed. When estimating the averaged eigenvalues of the $`\mathrm{\Omega }`$ matrix at the chiral fixed point, we disregard the approximants leading to purely real eigenvalues. For $`N=2`$ the final results of our calculations are:

$`\omega _\pm =`$ $`1.50(25)\pm i\mathrm{\hspace{0.17em}1.00}(45),`$ (4-loop CM) (32)
$`\omega _\pm =`$ $`2.05(35)\pm i\mathrm{\hspace{0.17em}0.80}(55),`$ (5-loop CM) (33)
and for $`N=3`$

$`\omega _\pm =`$ $`1.30(25)\pm i\mathrm{\hspace{0.17em}0.50}(35),`$ (4-loop CM) (34)
$`\omega _\pm =`$ $`1.55(25)\pm i\mathrm{\hspace{0.17em}0.55}(35),`$ (5-loop CM) (35)
The presence of imaginary parts comparable in magnitude with the real ones leads to the conclusion that the critical behaviors of frustrated two-dimensional Heisenberg and $`XY`$ systems are driven by focus fixed points. To strengthen this new and important issue, we repeat the analysis of the stability properties of the chiral fixed point by a Padé study, which confirms the previous results, although the statistics in this case is much less significant due to a large number of defective Padé approximants. All the results obtained are summarized in Table II. Note that in the past such a behavior was found in three dimensions within the three-loop RG approximation but only for unphysical values of $`N`$ ($`N=5,6,7`$.

In order to fully characterize the two physical systems discussed, we calculate the critical exponents as well. For the exponent $`\eta `$ we find values that agree with the results of the four-loop analysis . At the same time, we can not give a reasonable estimate for the critical exponent $`\nu `$ because of the strong oscillations observed with varying the working approximants.

Finally, we report in Table II the five-loop estimate of the crossover exponent $`\omega _{O(2N)}`$, finding that the Heisenberg fixed point is unstable for both $`N=2,3`$, in agreement with the previous four-loop calculation .

We apply the same kind of analysis also to the models with $`N4`$. All these systems are free of topological defects  and are well described by appropriate $`NL\sigma `$ models  which provide the existence of a stable chiral fixed point with critical exponents and corrections to the scaling equal to the ones of a system with an infinite number of components of the order parameter.

The zeros of the two $`\beta `$s for $`N=32,16,8,6,4`$ are shown in Fig. 3, 4, 5, 6, 7. From these figures it is clear that for large enough values of $`N`$ four fixed points exist: the gaussian, the $`O(2N)`$, the chiral, and the antichiral, with a topological structure similar to the one found in $`ϵ`$ expansion for $`N>N_c`$. The location of the chiral fixed point for various $`N`$ is reported in Tab. III.

From the $`NL\sigma `$ model it is expected $`\omega _1=\omega _2=2`$ . For this reason $`N^{}`$ should be less or at least equal to $`4`$. From our analysis we find that the chiral fixed point is stable for all the considered $`N`$. However, the estimates for $`\omega _i`$ are slighty different from $`2`$, also for large $`N`$ (e.g., for $`N=16`$ we find $`\omega _1=1.59(5)`$ and $`\omega _2=2.01(11)`$). Such a disagreement was already observed in the four-loop approximation, and it may be attributed to the nonanalyticities discussed in the previous subsection. We observe that for $`N8`$ all the stability eigenvalues are purely real, but for $`N=6`$ about the $`50\%`$ of the five-loop approximants possess a nonvanishing imaginary part, suggesting $`N^{}6`$. This result seems to disagree with the $`NL\sigma `$ prediction. Again this might be due to nonanalyticities.

To further approve the above conclusions, we apply, along with the conformal-mapping based, the resummation by means of the Padé-Borel technique. In the region where the fixed point is expected, many Padé approximants turn out to be defective, i.e., possessing dangerous poles. We find that the approximants $`[4/1]`$, $`[3/1]`$, and $`[2/2]`$ are nondefective in the majority of the cases. Considering these three Padé approximants with $`b=0,1,2`$ (being the more stable under variation of the number of loops), we finally have nine approximants for each $`\beta `$ function. The location of the chiral fixed point is reported in Table III, with all the defective cases discarded.

We also compute the exponents $`\eta `$ and $`\eta _t`$ with both the resummation techniques. We consider a lot of approximants and finally we choose the more stable ones with varying the number of loops. The found values are close to those of the $`O(N)`$ model for $`N3`$, but differ from the exact ones coming from $`NL\sigma `$ model $`\eta =0`$ and $`\eta _t=2`$ (e.g., we find $`\eta _t=1.95(20)`$ for $`N=16`$ and $`\eta =0.16(2)`$ for $`N=6`$). We think these discrepancies are due to nonanalyticities, as in the case of $`O(N)`$ models. To conclude the analysis of the fixed points, we study the stability of the Heisenberg fixed point by calculating the crossover exponent $`\omega _{O(2N)}`$. The five-loop results show that this fixed point is unstable, in agreement with the conclusion drawn earlier from the four-loop series .

### B Renormalization-group flow and crossover

In this subsection we demonstrate the structure of the RG flow for physical values of $`N`$. As it was already shown, a stable focus governs the critical behavior of these systems. In Fig. 8 ($`N=2`$) one can see the spiral-like approach of the stable fixed point which is peculiar for the fixed points with stability eigenvalues possessing nonvanishing imaginary parts. The shaded area is the region where the resummation procedure is expected to fail (i.e., $`z>4`$), since we do not take into account the singularity of the Borel transform closest to the origin that is on the real positive axis (cf. Eq (22)). Consequently, the RG trajectories given by the resummed perturbative series are quantitatively correct only within the unshaded areas. We report, nevertheless, the RG flows in other parts of the coupling constant plane in order to present a complete qualitative picture. Note, that the domain usually referred to as the sector of the first-order transitions, where the positiveness of the quartic form in the free-energy expansion breaks down, is given by the inequality $`z>2`$. In fact, however, because of the presence of the higher-order terms in this expansion keeping the system globally stable at any temperature, the true domains of the first-order transitions should be more narrow.

From Fig. 8 we see that for some initial values of $`(\overline{u},\overline{v})`$ the RG trajectories have the coordinate $`\overline{v}`$ that grows very fast at the beginning and seems to reach the region of the first-order phase transitions, but just before arriving there these trajectories drastically curve moving toward the stable chiral fixed point. These bizarre trajectories can imitate both the discontinuous phase transitions and strongly scattered effective critical exponents observed in numerous physical and computer experiments. To substantiate the former scenario, we calculate the numerical values which the effective exponent $`\eta _{eff}`$ takes along the RG trajectories running from the Gaussian fixed point to the chiral one. Corresponding curves are shown in Fig. 9. As is seen from this figure, the exponent oscillates within a large range, even near the stable fixed point.

We demonstrate also the RG flows for $`N=3`$ (Fig. 10) and for $`N=8>N^{}`$ (Fig. 11). In the former (physical case) the stable fixed point is a focus and phase trajectories approach it in a spiral-like manner. In the latter (unphysical case) the stable fixed point turns out to be a node and the trajectories are certainly not spirals. Also for $`N=3`$ we estimate the effective exponent $`\eta _{eff}`$ along the RG trajectories running from the Gaussian fixed point to the chiral one. The corresponding curves (cf. Fig. 12) show big oscillations and sometimes reach negative values. We note that these negative values are obtained only in the region where the resummation is not sure, i.e., $`z>4`$. Nevertheless, this strange crossover, if reproduced from the three-dimensional series , might explain some doubtful experimental results, where a negative $`\eta `$ was found.

To conclude, it is worthy to note that earlier the focuslike stable fixed points were found on the RG flow diagrams of the model describing critical behavior of liquid crystals  and of the $`O(n)`$-symmetric systems undergoing first-order phase transition close to the tricritical point . In those cases, however, the independent coupling constants had different scaling dimensionality and played essentially different roles in critical thermodynamics. Moreover, recently complex subleading exponents were found from the high-temperature expansion of Dyson’s hierarchical Ising model .

## IV CONCLUSIONS

We studied the critical thermodynamics of the two-dimensional $`N`$-vector chiral model in the five-loop RG approximation. Using the advanced resummation technique based upon Borel transformation combined with conformal mapping and use of Padé approximants, it was shown that the fixed point governing the chiral critical behavior is a stable node only for large enough $`N`$, $`N>N^{}`$. For smaller $`N`$ a chiral fixed point turns out be a stable focus and the RG trajectories approach this point in a spiral-like manner. Of particular importance is the fact that a focus-driven critical behavior was found for systems with physical values of $`N`$, $`N=2,3`$. Evaluating the critical exponent $`\eta `$ along some RG trajectories, we demonstrate that spiral-like approach of the chiral fixed point results in unusual crossover and near-critical regimes that may imitate varying critical exponents seen in physical and computer experiments.

## ACKNOWLEDGMENTS

We would like to thank Ettore Vicari for many useful discussions. The financial support of the Russian Foundation for Basic Research under Grant No. 01-02-17048 (E.V.O., A.I.S.) and the Ministry of Education of Russian Federation under Grant No. E00-3.2-132 (E.V.O., A.I.S.) is gratefully acknowledged. A.I.S. has benefited from the warm hospitality of Scuola Normale Superiore and Dipartimento di Fisica dell’Università di Pisa, where the major part of this research was done.


text:
# 1 Data Analysis

## 1 Data Analysis

The Cyg X-1 data presented here were obtained by the PCA onboard RXTE, mostly gained through our monitoring campaign 1998–2002 (for a description of the data see ). We split each light curve into segments of 1 s length and determine the mean flux of each segment. The segments are binned into 41 equally segmented flux bins, for each of which we calculate the power spectral density (PSD) of all contained lightcurve segments via the DFT (see, e.g., ).

From all periodograms of each flux bin the mean PSD is determined using standard methods. For better statistics we chose to take into account only flux bins containing at least 20 periodograms. Integrating the Poisson noise corrected PSD over the range $`\nu =1`$$`32`$ Hz, we arrive at the squared fractional rms variability. In a next step we multiply the fractional rms variability by the mean flux of the bin, to obtain the absolute rms variability $`\sigma `$.

## 2 rms-flux-Relation

For all bins of each observation we plot $`\sigma `$ over flux $`F`$ and fit a linear function in two different representations. By fitting $`\sigma =kF+a`$, two characterizing values are determined: slope $`k`$ and intercept $`a`$ on the $`\sigma `$ axis. The gradient of the $`\sigma `$-$`F`$-trend $`k`$ is equivalent to the fractional rms variability of the light curve.

Alternatively to intercept $`a`$, it was proposed to determine the physically meaningful intercept $`C`$ on the $`F`$ axis by fitting $`\sigma =k(FC)`$ . If $`k`$ is interpreted as the variable component of the lightcurve, then $`C`$ represents a second component of it which does not follow the linear $`\sigma `$-$`F`$-trend. In deriving results, it should be kept in mind that $`C`$, which is identical to $`a/k`$, and $`k`$ are not independent values for the characterization of the rms-flux-relation. Generally, there is a good linear relationship between $`\sigma `$ and $`F`$ in the hard state, which is reflected by the stable Pearson correlation coefficient in Fig. 1. We also notice coincidence of the change of the general long term behavior of Cyg X-1 from a “quiet hard state” to a “flaring hard state” in 1998 May  with a shift in the values of slope $`k`$ and the intercepts.

## 3 In the Soft State and during “Failed State Transitions”

The shape of the rms-flux-relation in the soft states of 1996 and 2001/2002 changes continuously between waviness and approximate linearity (see Fig. 2).

During “failed state transitions”(FST), i.e., during times when the source reached its intermediate state , the linear rms-flux relation breaks down. We examined four FST: 1998 July 15, 1999 Dec 05, 2000 Nov 03, and 2001 Jan 29 (indicated in Fig. 1 by dotted lines). Comparing these events with the neighboring hard state observations, we notice that the rms-flux-relations of 1999 December 05 and 2000 November 03 – being strictly linear before and after – change to an arch-like appearance during the FST itself (see Fig. 3). The other two examined FST also deviate from the linear relation, but these deviations are not as prominent.

## 4 Discussion and Conclusions

An explanation of the rms-flux-relation seems to be the modification of a theory which explains the PSD of GBHC by the superposition of accretion disk (AD) instabilities occuring at several radii in the AD that are propagated towards the innermost region of the disk . Local $`\dot{M}`$ variations, which are caused by short term changes of the disk viscosity and whose time scale is dependent on radius, are superimposed on long term $`\dot{M}`$ modulations. If the fractional amplitude of the $`\dot{M}`$ variations is independent of the long term $`\dot{M}`$, a linear flux-rms relation will be observed . The rms-flux-relation will break down if the corona is severely disturbed. This seems to be the case during the intermediate state of Cyg X-1, in agreement with earlier results for the behavior of X-ray time lags . Here, the observed X-ray time lag is much larger than during the normal hard state, which could indicate changes in the geometry of the AD corona, possibly related to the observed radio emitting outflow.