text:
# Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions

## Abstract

The nematic-to-isotropic orientational phase transition, or equivalently the $`RP^2`$ model, is considered in two dimensions and the question of the nature of the phase transition is addressed. Using powerful conformal techniques adapted to the investigation of critical properties of two-dimensional scale-invariant systems, we report strong evidences for a transition governed by topological defects analogous to the Berezinskii-Kosterlitz-Thouless transition in two-dimensional $`XY`$ model.

Liquid crystals may be seen as constituted of molecules essentially represented by long rigid rods. From maximization of entropy at high temperatures, all the molecule orientations are equally probable, independently of the neighbouring molecule directions and the system exists in an isotropic phase. At low temperatures a preferential orientation is more favourable in order to minimize interaction terms, and an ordered structure emerges. When order occurs along one space dimension only, the system is said to be nematic. Still at lower temperatures, other ordered phases can appear, e.g. smectic phases.

In a lattice model, each molecule may be represented by a unit vector $`𝝈_w`$ at site $`w`$ of an hyper-cubic lattice $`\mathrm{\Lambda }`$ of linear extent $`L`$. The $`𝝈`$’s live in a three-dimensional space attached to each lattice site. In the nematic phase, the preferential direction defines a unit vector, $`𝐧`$, called the director, and one can measure the deviation of molecule $`𝝈_w`$ with respect to the director by the scalar product $`𝝈_w𝐧=\mathrm{cos}\theta _w`$. Due to the local $`Z_2`$ symmetry (the rods are not oriented), one cannot distinguish between opposite directions $`\theta _w`$ and $`\theta _w+\pi `$, and $`\mathrm{cos}\theta _w`$ vanishes on average while $`\mathrm{cos}^2\theta _w`$ does not. In the disordered phase on the other hand, the angles are measured with respect to any arbitrary direction, and the thermal average of course leads to $`\mathrm{cos}\theta _w=0`$ and $`\mathrm{cos}^2\theta _w=\frac{1}{3}`$, so that $`\mathrm{cos}^2\theta _w\frac{1}{3}`$ represents a convenient order parameter. In the literature on liquid crystals, one usually defines the local order parameter by the second Legendre polynomial,

$$m(w)=P_2(𝝈_w𝐧)=P_2(\mathrm{cos}\theta _w).$$
(1)
This definition suggests to consider the following Hamiltonian to describe the nematic transition,

$$\frac{H}{k_BT}=\frac{J}{k_BT}\underset{w}{}\underset{\mu }{}P_2(𝝈_w𝝈_{w+\mu }),$$
(2)
where $`\mu `$ stands for the unit basis vectors of the lattice, $`𝝈_w𝝈_{w+\mu }=\mathrm{cos}(\theta _w\theta _{w+\mu })`$ is the scalar product between neighbouring vectors distant from one lattice spacing, and the interaction term $`JP_2(𝝈_w𝝈_{w+\mu })`$ is reminiscent from a dipole-dipole interaction. This Hamiltonian was introduced by Lebwohl and Lasher  as a lattice version of the mean field theory of Maier and Saupe , and its success came from its ability to reproduce the weak first order phase transition observed experimentally in the three dimensional nematic transition . In a more abstract context, this Hamiltonian is known as the $`RP^2`$ model, since at each lattice site is attached the manifold of directions in 3-dimensional space, also called the real projective space in 3 dimensions .

Like the non linear $`\sigma `$-model, this model possesses generically the symmetry group $`O(n)`$ which is non abelian for $`n3`$, and specifically the Lebwohl-Lasher or $`RP^2`$ model has a $`O(3)`$ symmetry. The question, as it was outstandingly formulated by Kunz and Zumbach , of the nature of the transition in the $`RP^2`$ model in two dimensions is still incompletely solved. The existence of a phase transition (at finite temperature) in two dimensional systems seems connected to the abelian nature of the underlying symmetry group, as both Ising and $`XY`$ models are famous examples, unlike the Heisenberg model. On the other hand, according to the Hohenberg-Mermin-Wagner theorem , models possessing a continuous symmetry group cannot exhibit any finite macroscopic magnetization with no magnetic field applied in dimensions 1 or 2 (we intentionally use the familiar terminology of magnetic systems). The two-dimensional $`XY`$ model is the most famous example and it exhibits a non conventional transition . In spite of the absence of long-range order (LRO) at low temperatures, the spin-spin correlation function decays algebraically with an exponent which increases monotonically with temperature  up to a temperature called after Berezinskii, Kosterlitz and Thouless . In this critical phase, macroscopic ordering is prevented by collective excitations, namely spin waves which nevertheless do not exclude a coherent orientation of spins at a smaller length scale. Together with the spin waves, localized excitations appear with increasing temperature. These are topological defects associated in pairs, like pairs of opposite charges in the low temperature phase of a two-dimensional Coulomb gas, and they perturb the spin field only locally. This is the usual meaning of the term quasi-long-range order (QLRO). The Berezinskii-Kosterlitz-Thouless transition is governed by unbinding of these defects which completely suppresses any type of long-range order at high temperature, hence the correlation functions decay exponentially like in an ordinary paramagnetic phase. This is the standard scenario of a topological transition. The puzzle becomes confused when one notices that $`XY`$ or Heisenberg models in three dimensions display conventional continuous transitions and the $`RP^2`$ model exhibits a first order transition while renormalization group treatment of non-linear $`\sigma `$ model predicts the absence of any transition at non-zero temperature in $`2d`$ (asymptotic freedom in the context of lattice gauge theories ) and a continuous one in $`3d`$ for all the three models. We also have to mention that the $`RP^1`$ model (the same as given in equation (2), but with two-component vectors $`𝝈_w`$) exactly coincides with the $`XY`$ model. The question of a possible topological transition in the two-dimensional (non abelian) Lebwohl-Lasher model is thus particularly attracting and was already addressed in previous studies . In their remarkable work, Kunz and Zumbach  concluded in 1992 in favor of such a topological transition scenario, essentially on the basis of qualitative arguments (pairing of topological defects at low temperature where they carry most of the energy in the system, sharp increase of the density of defects and apparent discontinuity of the rotational rigidity modulus at the transition, finite cusp in the specific heat, proliferation of unbinded defects at high temperature). Even though they performed a careful and sizeable study, they were unfortunately not able to decide conclusively between essential singularities or standard power laws - though their preference was for the first case - for the correlation length and the susceptibility when approaching the transition from the high temperature phase. This is essentially due to the limited possibilities of computers ten years ago, since the authors already took care about potential critical slowing down problems as they adapted the Wolff cluster algorithm to the $`RP^{n1}`$ model.

Ten years later, we want to address the same question of the nature of the phase transition of the two-dimensional Lebwohl-Lasher model. Since we believe that the conclusions of Kunz and Zumbach will hardly be improved, even with more powerful facilities, it is necessary to reconsider the problem from a different point of view. We will thus assume the existence of a critical phase at low temperatures and then follow the line of the behaviour of the $`XY`$ model to predict consequences of the above mentioned assumption, consequences which may be compared easily to numerical results. Of course, the test must discriminate between different scenarios, starting from the confirmation (and thus the characterization) of a topological transition, a conventional continuous transition, a first-order one, or no transition at all. The existence of a scale-invariant low temperature critical phase is characteristic from the first situation. Such a system must thus be conformally invariant at any temperature below the transition $`T_{KT}`$ (we will abusively keep the terminology adapted to the $`XY`$ model), so it becomes advantageous to deduce the functional expression of the correlation functions or density profiles in a restricted geometry adapted to numerical simulations from a conformal mapping $`w(z)`$:

$$𝝈_{w_1}𝝈_{w_2}=|w^{}(z_1)|^{x_\sigma }|w^{}(z_2)|^{x_\sigma }𝝈_{z_1}𝝈_{z_2}.$$
(3)
Here, $`w`$ labels the lattice sites in the transformed geometry (the one where the computations are really performed), $`z`$ are the corresponding points in the original one (usually the infinite plane where $`𝝈_{z_1}𝝈_{z_2}|z_1z_2|^{\eta _\sigma }`$ takes the standard power-law expression), and $`x_\sigma =\frac{1}{2}\eta _\sigma `$ is the scaling dimension associated to the scaling field $`\sigma `$. The interest of such an approach lies in the full inclusion in the functional expression of the changes due to shape effects. The most famous example is the exponential decay of the correlation functions at criticality along a strip of finite width, unlike the algebraic decay in the infinite plane. For simplicity reasons, it is even more convenient to work with density profiles in a finite system with symmetry breaking fields along some surfaces in order to induce a non vanishing local order parameter in the bulk. In the case of a square lattice $`\mathrm{\Lambda }`$ of size $`L\times L`$, with fixed boundary conditions along the four edges $`\mathrm{\Lambda }`$, one expects

$`m_{\mathrm{Fbc}}(w)`$ $`=`$ $`P_2(𝝈_w𝐡_{\mathrm{\Lambda }(w)})_{\mathrm{Fbc}}[\kappa (w)]^{\frac{1}{2}\eta _\sigma }`$ (4)
$`\kappa (w)`$ $`=`$ $`\mathrm{Im}\left[\mathrm{sn}{\displaystyle \frac{2\mathrm{K}w}{L}}\right]`$ (5)
$`\times `$ $`\left|\left(1\mathrm{sn}^2{\scriptscriptstyle \frac{2\mathrm{K}w}{L}}\right)\left(1k^2\mathrm{sn}^2{\scriptscriptstyle \frac{2\mathrm{K}w}{L}}\right)\right|^{\frac{1}{2}}`$
where $`\mathrm{Fbc}`$ specifies that the boundary conditions are fixed. This expression easily follows from the expression of the order parameter profile decaying in the upper half-plane from a distant surface of spins constantly fixed in a given direction, $`m(z)=P_2(𝝈_z𝐡_{\mathrm{\Lambda }(z)})_{\mathrm{uhp}}y^{x_\sigma }`$, and from the conformal transformation of the upper half-plane (uhp) $`z=x+iy`$ ($`0y<\mathrm{}`$) inside a square $`w=u+iv`$ of size $`L\times L`$ ($`L/2uL/2`$, $`0vL`$) with open boundary conditions along the four edges, realized by a Schwarz-Christoffel transformation

$$w(z)=\frac{L}{2\mathrm{K}}\mathrm{F}(z,k),z=\mathrm{sn}\left(\frac{2\mathrm{K}w}{L}\right).$$
(6)
Here, $`F(z,k)`$ is the elliptic integral of the first kind, $`\mathrm{sn}(2\mathrm{K}w/L)`$ the Jacobian elliptic sine, $`\mathrm{K}=K(k)=F(1,k)`$ the complete elliptic integral of the first kind, and the modulus $`k`$ depends on the aspect ratio of $`\mathrm{\Lambda }`$.

Our strategy is now to fit numerical data of the order parameter profile against expression (4). Like in the previous study of Kunz and Zumbach  the resort to a cluster update algorithm is necessary in order to prevent the critical slowing down, all the spins of clusters (build through intermediate bond variables) being updated simultaneously. The algorithm becomes particularly efficient if the percolation threshold of the bond process occurs at the transition temperature of the spin model, which ensures the updating of clusters of all sizes in a single MC sweep. For $`O(n)`$ models, Ising variables are defined in the Wolff algorithm by the sign of the projection of the spin variables along some random direction. The bonds are then introduced through the Kasteleyn-Fortuin random graph representation . When one uses fixed boundary conditions, a difficulty occurs and the Wolff algorithm should become less efficient, since close to criticality the unique cluster will often reach the boundary and no update is made in this case. This is circumvented by the following trick: even when the cluster reaches the fixed boundaries, it is updated - and so are the boundary spins - and the order parameter profile is then measured with respect to the new direction of the boundary spins, $`m_{\mathrm{Fbc}}(w)=P_2(𝝈_w𝝈_\mathrm{\Lambda })_{\mathrm{Fbc}}`$.

Using this procedure, we studied systems of size $`48\times 48`$ up to $`200\times 200`$. For the measurement of the order parameter profile, we discarded $`10^6`$ Wolff sweeps for thermalization, and the measurements were performed on $`10^6`$ production sweeps. For reasons which are made obvious below, the energy density required a better statistics and the measurements were obtained over $`16.10^6`$ sweeps.

In order to underscore the existence of a line of marginal fixed points in the low temperature phase, we first check qualitatively the expression of the energy-energy correlations from the behaviour of the energy density profile. The energy density at site $`w`$ is for example defined as the average value of the energies of the four links reaching $`w`$:

$$\epsilon _w=\frac{1}{2d}\underset{\mu }{}\left[P_2(𝝈_{w\mu }𝝈_w)+P_2(𝝈_w𝝈_{w+\mu })\right].$$
(7)
The existence of a regular contribution in the energy density makes the calculation a bit more subtle than what presented in equation (4). This regular contribution $`\epsilon _0(T)`$ which depends on $`T`$ cancels after a suitable difference between profiles obtained with different conditions at the boundaries (free ($`\mathrm{f}\mathrm{bc}`$) and fixed ($`\mathrm{Fbc}`$) boundary conditions). Although it makes the numerical computation longer in order to reach some satisfactory accuracy, this makes possible to extract the singularity associated to the energy density:

$`\epsilon _z_{\mathrm{Fbc}}`$ $`=`$ $`\epsilon _0(T)+_{\mathrm{Fbc}}(T)y^{\eta _\epsilon (T)/2},`$ (8)
$`\epsilon _z_{\mathrm{f}\mathrm{bc}}`$ $`=`$ $`\epsilon _0(T)+_{\mathrm{f}\mathrm{bc}}(T)y^{\eta _\epsilon (T)/2}.`$ (9)
This is clearly illustrated in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 1 where convergence towards the same temperature-dependent constant $`\epsilon _0(T)`$ is shown, with amplitudes of the singular terms having opposite signs, therefore a simple difference of the quantities measured in the square geometry,

$$\mathrm{\Delta }\epsilon (w)=\epsilon _w_{\mathrm{Fbc}}\epsilon _w_{\mathrm{f}\mathrm{bc}}\mathrm{\Delta }\times [\kappa (w)]^{\frac{1}{2}\eta _\epsilon (T)}$$
(10)
leads to the value of the thermal scaling dimension $`\eta _\epsilon (T)`$. The right part in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 1 presents a log-log plot of the difference $`\mathrm{\Delta }\epsilon (w)`$ vs $`\kappa (w)`$ at two temperatures below $`T_{\mathrm{KT}}`$ and one above which shows that the functional expression used is no longer valid, as expected, in the paramagnetic phase. Due to the strong fluctuations, in the QLRO phase the data scatter around straight lines which represent the slopes $`[\kappa (w)]^2`$. This figure, though not definitely conclusive, confirms that the exponent of the decay of energy-energy correlations keeps a constant value $`\eta _\epsilon (T)=4`$ in the low-temperature phase of the $`RP^2`$ model, confirming that like in the case of the $`XY`$ model, the temperature is a marginal field, responsible for the existence of a critical line in the whole low-temperature phase. It thus implies a thermal scaling exponent $`x_\epsilon =dy_t=2`$ which ensures a vanishing RG eigenvalue $`y_t=0`$ (up to $`T_{\mathrm{KT}}`$ where it is consistent with an essential singularity of $`\xi `$ above the KT point, as suspected by Kunz and Zumbach ). The energy-energy correlation function in the plane should thus decay algebraically as

$$\epsilon _{z_1}\epsilon _{z_2}|z_1z_2|^{\eta _\epsilon },$$
(11)
with $`\eta _\epsilon (T)=2x_\epsilon =4`$ $`T<T_{\mathrm{KT}}`$.

Fig. 1: MC simulations of the $`2d`$ $`RP^2`$ model inside a square of $`100\times 100`$ spins ($`16.10^6`$ MC sweeps after cancellation of $`10^6`$ for thermalization). Several temperatures below the transition temperature are shown, and one above ($`T=0.6`$). Left: local energy density vs the rescaled variable $`\kappa \left(w\right)`$. Right: log-log plot of the difference $`\mathrm{\Delta }\epsilon \left(w\right)`$.

From the existence of a marginal line, one may suspect that the other scaling dimensions should continuously vary with the marginal field. The order parameter profile thus has to obey equation (4) in the whole low temperature phase, but with an exponent $`\eta _\sigma (T)`$ which depends on $`T`$. Equivalently, a log-log plot of $`m_{\mathrm{Fbc}}(w)`$ vs $`\kappa (w)`$ must display straight lines with different slopes below $`T_{\mathrm{KT}}`$. This is exactly what is observed in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 2. Again, the curves start to deviate from the straight line when the system enters the high temperature phase.

Fig. 2: MC simulations of the $`2d`$ $`RP^2`$ model inside a square of $`100\times 100`$ spins ($`4.10^6`$ MC sweeps after cancellation of $`10^6`$ for thermalization). The order parameter density is plotted against the rescaled distance at different temperatures below the KT transition temperature and slightly above. The numbers on the right give the value of $`k_BT/J`$, and 0.537 (the estimate of Kunz and Zumbach ) appears to be overestimated.

Fig. 3: Left: Correlation function decay exponent of the $`2d`$ $`RP^2`$ model as a function of the temperature. The vertical dashed line is the estimate of the transition temperature of Kunz and Zumbach . Right: value of the $`\chi ^2`$ per degree of freedom resulting from the power law fits of the data in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 2. The sharp increase is an indication of the location of the transition temperature (already below the estimate of Kunz and Zumbach ).

The slopes measured in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 2 are reported as a function of the temperature in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 3.

Fig. 4: Dependence of the rotational rigidity modulus with the $`\eta _\sigma `$ exponent.

Another signature of this mechanism is the behaviour of the rotational rigidity modulus. This latter quantity generalizes the helicity modulus $`\mathrm{{\rm Y}}`$ in the $`XY`$ model, which measures the quadratic response in the free energy of the system to a twist accross the sample. This is generalized to the $`RP^2`$ model by measuring the change in free energy when a rotation of angle $`\varphi `$ around some axis in spin space is applied to the system: $`F(\varphi )F(0)=\mathrm{\Sigma }\varphi ^2+O(\varphi ^4)`$. This expression defines the rotational rigidity modulus $`\mathrm{\Sigma }`$. In the $`XY`$ model, there exists a universal relation between the helicity modulus and the correlation function exponent, $`\eta _\sigma =k_BT/2\pi J\mathrm{{\rm Y}}`$ , which appears as a consequence of the Kosterlitz recursion relations. In the Lebwohl-Lasher model, the same type of behaviour is checked in figure   Evidence for a topological transition in nematic-to-isotropic phase transition in two dimensions 4 where one observes a linear dependence (the larger the system size, the better the linear behaviour) of $`\eta _\sigma `$ with $`T/\mathrm{\Sigma }`$, correctly fitted at low temperatures by $`\eta _\sigma =0.117k_BT/J\mathrm{\Sigma }`$.

To conclude, we mention that the most remarkable feature is a complete analogy with the two-dimensional $`XY`$ model  where the transition is mediated by the defects. The exponent $`\eta _\epsilon `$ keeps a constant value while $`\eta _\sigma `$, associated to the order parameter, starts from zero at $`T=0`$ and increases linearly with $`T`$ at low temperatures where a spin wave approximation should capture the essentials of the behaviour of the system. The influence of pairs of topological defects, which would appear in increasing number when the temperature increases, is probably responsible for the deviation from the spin wave approximation and of the sharper increase of $`\eta _\sigma `$, and the transition is presumably completely governed by unbinding of these defects, like in the Kosterlitz-Thouless scenario. The order parameter exponent at the transition $`k_BT_{\mathrm{KT}}/J0.52`$ takes a value $`\eta _\sigma (T)0.40(2)`$. The relation between rotational rigidity modulus and the exponent $`\eta _\sigma (T)`$ also seems completely coherent with what happens in the $`XY`$ model, giving one more evidence of the topological nature of the transition.

Acknowledgement: This work is supported by the french-venezuelian PCP program ‘Fluides pétroliers’. Support from the CINES under project c20020622309 is also gratefully acknowledged.


format: text

text:
# Formation of massive binaries

## 1. Introduction

Most massive stars, if not all, are (or were) members of binary and multiple systems. Ignoring the multiplicity of massive stars can introduce severe problems when comparing observations and stellar evolution theory. For example, the stellar luminosity in the range 20 – 60 solar masses scales with the third power of the mass, and the rate of production of Lyman continuum photons scales with the fourth power of the mass. This means that a 50 solar mass star has 4 times the luminosity of an (unresolved) pair of two 25 solar mass stars, and 8 times the amount of Lyman continuum photons per sec. This example, drastic as it may be, shows the importance of recognising how the mass of a system is distributed among its components.

The basic properties of stars are imprinted at the time of stellar formation. Although not the main theme of this IAU Symposium, any conference on the evolution of massive stars would be incomplete without considering the very early formative phases. Of particular interest for predicting the odyssey of an ensemble of massive stars are the initial mass function (IMF) and the initial binary star parameters. While the IMF of massive stars has been a hot topic for a long time (e.g. Lequeux 1979, Hunter 1995, Scalo 1998, Kroupa 2002), the binary properties of massive stars have been more of a side problem, despite several previous spectroscopic studies (see the brief review in Preibisch, Weigelt, and Zinnecker 2001)

However, an important clue to the origin of high-mass stars has come from the observation that they occur predominantly in binary and multiple systems, apparently more so than low-mass stars (Preibisch, Weigelt, and Zinnecker 2001) and often with very short orbital periods (Mermilliod and Garcia 2001). Therefore it seems that the understanding of the formation of high-mass stars is intimately tied to understanding the origin of their binary nature, be it via fragmentation of a massive accretion disk (Yorke and Sonnhalter 2002), via Bondi-Hoyle accretion onto a lower-mass binary (Maeder and Behrend 2002), via failed stellar or protostellar collisions (Zinnecker and Bate 2002) or via dynamical N-body interactions in dense stellar clusters (Bonnell and Bate 2002).

In this paper, we review the observational situation for unevolved massive binaries on the Main Sequence in young clusters and a number of theoretical ideas regarding their formation. Finally we discuss some implications and predictions, including a speculation of how the binary properties might change when we consider low metallicity environments.

## 2. Generalities and Motivation

### 2.1. The definition of massive stars and massive binaries

We believe there is general agreement about the lower mass limit of a massive star. It is 8 solar masses. Why is this a good definition? There are at least three reasons: The first is the fact that stars below 8 solar masses (B3V) do not emit any Lyman continuum photons. Thus in terms of spectral type we define massive stars to be earlier than B3V. (A B2.5V star emits 1044 Lyc photons/sec, an O6V star about 1049). Secondly, for initial stellar masses above 8 M there is no pre-Main Sequence phase (Palla and Stahler 1993), i.e. the star ignites H-burning before it finishes accumulating its final mass by accretion from the protostellar gas/dust envelope or circumstellar disk. Thirdly, 8 solar masses is the dividing line between the initial progenitor mass of white dwarfs and type II supernovae. Note that the onset of a significant stellar wind mass loss occurs at around 18 solar masses (B0V), which is the dividing between O and B stars. Thus only O-stars have stellar winds on the Main Sequence which can influence their evolution.

The definition of massive binaries is more tricky. The issue is whether both components must exceed the above 8 M limit, or just the primary of the two. I suggest “just the primary”. We will use this latter definition here, especially when it comes to the search for physical companions to OB-type stars. Thus a 50 M \+ 5 M pair or a 10 M \+ 1 M pair would both be considered as a massive binary according to our definition. Of course it will be difficult to detect such pairs as close spectroscopic binaries but we shall see that it is possible to resolve them as visual companions, e.g. in the Orion or M16 young cluster, using speckle interferometry or adaptive optics.

### 2.2. The theme of this contribution

The underlying theme of this contribution is to study the massive binary properties as a function of the “environment”, in the hope to understand the physics of their formation. For example, comparison of the properties of massive binaries, in young clusters vs. OB associations ought to reveal the role of dynamical gravitational interactions (N-body evolution). These effects are of great importance in dense star clusters, while they should play only a minor part in looser aggregates. Another aspect is the origin of runaway OB stars, such as $`\zeta `$ Puppis (O4If). Has it been ejected by an asymmetric supernova explosion or has it been released from a young cluster after suffering a tight binary-single star encounter? Finally, we are interested in the initial binary properties as a function of metallicity, particularly in view of the dynamics of the protoglobular clusters in the early universe and in the LMC/SMC. While we will not answer the questions, we want to at least raise them!

As we are concerned with the formation of massive binaries, we will concentrate on observations of young unevolved systems only. The initial or primordial binary properties, as they are sometimes called, consist of the binary frequency, the mass ratios, and the component separations. Measurements of eccentricities are also desirable, as they may hide clues to formation processes (e.g. tidal capture should initially result in large eccentricities for close binaries). One of the key questions is: what is the critical semi- major axis or critical orbital period for massive binaries to become interacting binaries later, i.e. for binary-related Roche lobe mass exchange to occur when the primary star inflates to become a supergiant? It depends on the mass ratio of the system, but for a mass ratio close to unity the critical stellar radius for Roche lobe overflow is around 0.4 of the semi-major axis. This implies that massive binaries, with primaries which can reach radii up to 6 AU during their post-Main Sequence evolution, can become semi-detached interacting binaries for systems with semi-major axes under 15 AU (i.e. periods under 10 years). What then is the fraction of all massive binaries that will turn into interacting ones? My guess is 25 – 50 %, considering the orbital period distribution and mass ratio distribution of bright O-stars found in the survey of Mason et al. (1998).

## 3. Review of observations

### 3.1. Definition of the companion star fraction (CSF)

Many massive stars are actually members of multiple systems, e.g. a hierarchical triple or a Trapezium type quadruple. In order to properly account for these higher order systems, we introduced the term “companion star fraction (CSF”) defined as (B + 2T + 3Q) / (S + B + T + Q), see Reipurth and Zinnecker (1993). Here B is the number of binaries (1 companion), T the number of triples (2 companions), and Q is the number of quadruple systems (3 companions). For example: 1 single, 1 binary, 1 triple, and 1 quadruple system in a cluster of 4 multiple stars like in the Trapezium system $`\theta ^1`$ Ori (D,C,B,A) in the center of the Orion Nebula cluster yields CSF = 1.5 (Preibisch et al. 1999).

This means that a massive star in Orion has on average more than 1 companion, which is significantly different from the case of low-mass pre-Main Sequence stars in Orion, where we typically find CSF = 0.5 (Prosser et al. 1994, Padgett et al. 1998). This discrepancy is the first hint that the formation process of low-mass and high-mass binary systems may be different.

### 3.2. Multiplicity of massive stars in young clusters

#### The Orion Nebula Cluster

The multiplicity of the 8 massive stars in the Orion Nebula Cluster has been studied by Preibisch et al. (1999), using bispectrum speckle interferometry in the near-infrared K-band. These include the 4 Trapezium stars $`\theta ^1`$ Ori A, B,C, D, $`\theta ^2`$ Ori A and B as well as LP Ori and NU Ori. $`\theta ^1`$ C, the exciting star of the Orion Nebula has a 33 mas visual companion, previously discovered by Weigelt et al. (1999). Note that $`\theta ^1`$ C is a peculiar star, showing strong optical and X-ray variability with a period of 15.4 days (Stahl et al. 1996, Babel & Montmerle 1997), probably related to a kG magnetic field misaligned with the rotation axis (Donati et al. 2002). It also exhibits radial velocity variations, which have been interpreted as a spectroscopic companion which could in fact be identical to the visual companion, if the orbit is sufficiently eccentric (Donati et al. 2002). $`\theta ^1`$ A and B have long been known to be spectroscopic binaries (Bossi et al. 1989, Abt et al. 1991), but in addition they also have one and two visual companions, respectively. Simon et al. (1999) even found another infrared companion to $`\theta ^1`$ B, so that this system is a Trapezium-type system within the Trapezium. $`\theta ^2`$ A is a triple system (spectroscopic binary with a 0.2” visual companion; Abt et al. 1991, Petr et al. 1998); the same is true for NU Ori. $`\theta ^1`$ D and $`\theta ^2`$ B are apparently single stars. Adding up the numbers, we find a total of 12 companions, implying CSF = 12/8 = 1.5. Interestingly, in all resolved systems the K-band flux ratios are quite small, i.e less than 1/3. This suggests that these visual companions are all significantly less massive than the primary OB stars.

#### The M16 cluster

Duchene et al. (2001) have searched for visual binaries among the high mass stars in the M16 cluster (NGC 6611), using adaptive optics. They list some 40 stars of spectral type earlier than B3V that are members of the cluster (age 2 Myr, distance 2 kpc). They find some 10 companions (25 %), with separations in the range 0.5” to 1.5” (1000 – 3000 AU), all but two with estimated mass ratios below q = 0.2. They conclude that these results are consistent with the assumption that companion masses are randomly drawn from the initial mass function. The lack of a strong dependence of the properties of the secondaries on the mass of the primaries lead them to suggest that the canonical picture of cloud fragmentation can explain the binary properties, and there is no need to invoke a different picture for high-mass stars and low-mass stars. This being said, we caution that the low-mass cluster members provide a serious source for false binaries (most companions are physical companions only at the 2 sigma confidence level).

As for spectroscopic binaries, Bosch, Morrell, and Niemela (1999) have studied the radial velocities of the 10 earliest type stars in the NGC 6611 cluster. One short-period (W412, 4 days) and two longer period SB2s (W175, W197) were found, but three more stars may be variable, among them the most massive object (W205, O4V).

#### NGC 6231 and other clusters

New statistical results about massive spectroscopic binaries in the young OB cluster were presented by Mermilliod and Garcia (2001) at the IAU-Symp. 200 (see also Garcia and Mermilliod 2001):

1) The spectroscopic binary frequency among 8 clusters rich in O-stars (i.e. more than 5 O-star members) varies from extremely high, of the order of 80 %, to rather low, of the order of 15 %, seemingly anti-correlated with stellar density in the clusters. NGC 6231, a loose cluster, contains 11 SB out of 14 stars, while Trumpler 14, the densest cluster in the sample, contains only 1 SB out of 7 O-stars. An intermediate case is Trumpler 16, with 7 SBs out of 20 O-star members (35 %).

2) Almost all O-stars in the 14 young clusters with few O-stars (i.e. less than 3 O-star members) are spectroscopic binaries, often SB2, and even eclipsing. These massive binaries are usually members of hierarchical triple or quadruple systems, or of trapezia, and are often located at the cluster center. $`\theta ^1`$ Ori belongs to this category, although $`\theta ^1`$ C is not a massive close SB.

3) The orbital periods of the spectroscopic binaries in the O-star rich clusters are concentrated in the range 4 – 5 days, judging mainly from the NGC 6231 and Trumpler 16 clusters; while in the O-star poor clusters there is an accumulation of orbital periods around 3 $`\pm `$ 1 day. The orbital eccentricities are high, in the range 0.2 – 0.6, at least in NGC 6231.

These fascinating and surprising facts challenge our views of massive binary formation! In the next section, we take up the challenge, contrasting these clues to some formation models.

## 4. Review of theories

In essence, we can distinguish 4 types of formation theories:

### 4.1. Fragmentation of massive disks around massive stars

Although direct observational evidence for the presence of circumstellar disks around massive stars is missing (there are indirect hints such as massive outflows, e.g. Beuther et al. 2002, Henning & Stecklum 2002), let us assume that the disks in question exist, as is known to be the case around low-mass stars. Support for this scenario comes from the 2D rotating collapse calculations by Yorke and Sonnhalter (2002); see also Yorke (2002). These authors provided simplified numerical simulations of the accretion of a massive protostar for a series of initial cloud cores (30, 60, 120 M) with initially uniform rotation and no magnetic fields. Frequency dependent radiative transfer was included. Even with rotation, radiation pressure on dust is still a limiting factor for the final mass of the star-disk system. For example, the collapse of a 60 M cloud leads to a 33 M star-disk system, with a substantial fraction (about 1/3) still in the disk. Such a massive disk can be gravitationally unstable, provided it is cold enough (formally described by the Toomre criterion). 3D simulations would be needed to follow the evolution (fragmentation) of the disk, but it seems clear that the mass in the disk can only lead to binary systems with mass ratios less than 1/3, as most of the mass goes to the central object. Therefore, disk fragmentation cannot yield nearly equal mass spectroscopic binaries which are so prevalent among the observed massive binaries. Also, disk fragmentation produces only wide binaries, of order 10 – 100 AU, comparable to the radial disk extent (cf. the VLA observations of Shepherd et al. 2001). This mode may then explain the $`\theta ^1`$ C Ori companion, but not the SB2 short orbital periods, which correspond to sub-AU semi-major axes.

### 4.2. Bondi-Hoyle accretion onto a low-mass protobinary

Maeder and Behrend (2002) proposed a model for massive star formation based on a growing accretion rate of the central star, according to $`\dot{M}M^{1.5}`$. Realising that the model would have to accommodate the high incidence and short periods of massive binaries, they extended their model to include preferential binary formation among massive stars by Bondi-Hoyle accretion from the protocluster gas reservoir. Their idea is as follows: The Bondi-Hoyle accretion rate for an initial binary star (formed by fragmentation or otherwise) is higher than for a single star of the same mass, particularly if the Bondi-Hoyle radius of gravitational influence ($`b=2GM/v^2`$, $`M`$ is the binary mass and $`v`$ the center-of-mass velocity w.r.t. the gas) exceeds the semi-major axis of the initial binary system. Then these systems will “have more of a chance to move to the top of the mass scale and they will do so faster than the other single objects. Since mass is accreted on the protobinary system, the semi-major axis will decrease and the orbital periods will also decrease” (Maeder and Behrend 2002). The problem with this scenario is that the Bondi-Hoyle accretion rate under normal circumstances (gas density ca. 105cm-3 and relative velocity ca. 2 km/s, appropriate for a protocluster cloud) is far too small to supply the required growth of 10-5M /yr for solar mass protostars, unless the relative velocity goes to zero. This could happen if the orbital motion of the protobinary partly compensates the relative motion of the gas w.r.t. the centre-of-mass, and each binary component is effectively at rest w.r.t. the ambient gas, for part of the orbit. Clearly the situation is very complicated and detailed simulations of binary accretion rates are required (cf. Klessen 2001). Another issue here is whether the Bondi-Hoyle runaway accretion model can explain the evolution towards equal mass binary components, perhaps involving a circumbinary disk (cf. Bate and Bonnell 1997).

### 4.3. Encounters and tidal capture of intermediate-mass protostars

This scenario was originally invented to avoid the radiation pressure problem of massive star formation (Bonnell, Bate, and Zinnecker 1998). However, it was soon realised that the collisional build-up of massive stars could also help to form massive almost-touching eccentric binaries, if the kinetic energy of a parabolic encounter could be dissipated as tidal energy (cf. Fabian, Pringle, and Rees 1975) in a gravitationally focussed grazing fly-by (Zinnecker and Bate 2002). There are several problems with this scenario.

The first one is the extraordinary high stellar density in the center of a cluster that is needed, so that the collisions occur at frequent enough intervals to be shorter than the stellar evolution time of massive stars (few Myr). Stellar number densities of the order $`10^8stars/pc^3`$ are required, which can only be reached, if a protocluster of intermediate-mass stars undergoes significant contraction and revirialisation owing to accretion onto the stars (see Bonnell et al. 1998). A factor of 10 contraction boosts the stellar number density by a factor of 1000 and shortens the collision timescale by a factor of 1000. This can happen, but the question is: does it happen? No observational support is so far available, although some hot molecular cores (typical masses ca. 100 M and typical sizes ca. 0.1 pc) harbour hyper-compact HII regions (sizes ca. 0.01 pc), indicative of massive stars in a very density-peaked gas-rich environment (Kurtz & Franco 2002). No-one knows what is going on inside those very compact radio and mm continuum sources, as they are so heavily obscured by dust (e.g., IRc2-I in Orion OMC1, see Gezari et al. 1998). It is conceivable that collisional events and ’cannibalism’ are taking place inside, revealed by eruptive phenomena (e.g., the famous Orion H2 fingers), maser emission, and thermal infrared flares, as discussed by Bally (2002).

The second problem with this model is the difficulty to first promote and then avoid stellar mergers. Only in the last step of this runaway process do we require the stellar mergers to fail and to form a binary instead.

The third problem is that the formation of massive binaries seems to require an extreme fine tuning, especially in terms of the impact parameter. Indeed the cross-section for mergers is always larger than that for binary formation. One way out is disk formation during the tidal disruption of the lower-density star in a non head-on collision with a higher-density star (see the example, shown in Zinnecker and Bate 2002, of a 3 M pre-MS star disrupted by a 10 M MS-star in a grazing collision). Disk formation increases the cross-section for dissipative capture of the next intruder and may also be a general feature of stellar mergers due to the induced rapid stellar rotation (near break-up for realistic non-zero impact parameters). We conclude that tidal capture during near-miss mergers is a viable mechanism to form massive close binaries, but it requires exceptional circumstances.

Finally, there is the issue whether it is pre-Main Sequence stars (perhaps with circumstellar disks) or protostellar cores with extended accretion envelopes that are made to collide (Stahler, Palla, and Ho 2000). These protostellar envelopes offer a much higher interaction cross-section, and provide orbital drag for collisional binary formation (Silk 1978). This variant of the above failed merger picture must be kept in mind, as it works in less extreme and perhaps more realistic conditions.

### 4.4. N-body interactions in a young cluster

Bate and Bonnell (2002) recently simulated the gas accretion onto 1000 stars in a gas-rich protocluster using 3D SPH calculations. Accretion forces the cluster to contract, leading to a high-density core where stellar mergers can occur. Small-N groups of massive and intermediate mass stars form due to the self-gravity of the gas, and binary formation in these groups is common, occurring through dynamical three-body capture (different from tidal capture). To begin with, typically a massive star has a lower-mass wide companion. With time, an exchange reaction with a more massive third object occurs, thus forming an almost equal mass but still wide binary. This wide binary then shrinks (“hardens”) as it takes up most of the binding energy of the small group being kicked by the other lower-mass group members most of which tend to escape from the group. The typical final outcome is a tight equal mass massive binary, with a lower-mass wide companion (Bate, Bonnell, and Bromm 2002). This picture explains many aspects of massive binaries, including the average companion star fraction ($`CSF=1.5`$) and the prevalence of short-period spectroscopic systems. Also, the timescales for the N-body dynamics are short enough for the whole process to occur on or before the zero age Main Sequence. Therefore we favour this scenario as our best working hypothesis. There is one nagging doubt, however: one would expect the massive close binary frequency to be correlated with the overall stellar number density in young clusters, which is not what is observed. In fact, the contrary (an anti-correlation) is observed. We have no real explanation at present, how this environmental effect fits into the N-body picture, unless most of the most massive binaries in the densest young clusters have managed to merge into a single object (Zinnecker 1986).

## 5. Implications and Predictions

Massive binaries may be the main reason for the existence of many runaway OB stars and field OB stars, expelled from young clusters. The question whether all field O-stars (10 – 25 % of all O-stars) are of some runaway origin is an extremely crucial one; if not, then some O-stars must be born in isolation rather than in clusters, and many conclusion about massive star formation and indeed massive binary formation must be reconsidered. For example, the collision scenario is clearly in trouble if O-stars can form outside clusters. As it is, however, we know of at least two examples where B0 stars lie more or less equidistant and in opposite directions of a young cluster (S255, Zinnecker, McCaughrean, and Wilking 1993; MonR2, Carpenter 2000). The prediction (to be tested by future observations) is that at least one star in each pair of ejected stars must itself be a binary, to account for the required recoil energy and momentum; see the simulations of Kroupa (2000) who also modelled the fraction of massive stars escaping from a young cluster, thus mimicking an OB field star population (Fig. 1).

We end with a speculation on the frequency and orbital separation of massive binaries in lower metallicity star forming regions (e.g. LMC/SMC). We suspect that at lower metallicity the angular momentum problem of star formation is even more severe than for molecular clouds with solar composition, implying that fragmentation must lead to more and/or wider binaries; perhaps even to faster rotation of the binary components. The chain of thought goes as follows: lower metallicity implies fewer heavy elements which are the main donors of free electrons in quasi-neutral clouds. Fewer electrons imply a lesser fractional degree of ionisation. The smaller the degree of cloud ionisation, the shorter is the timescale of retaining a magnetic field during cloud contraction and collapse. Thus the effect of magnetic braking to help solve the angular momentum problem (e.g. Mouschovias 1977) is reduced, and the final step of star formation has to get rid of the extra angular momentum by forming wider binary systems in which the remaining excess angular momentum is stored. If this is true, it will have many consequences for low-metallicity populations, but this beyond the scope of this review.

## Questions/Answers:

Peredur Williams:
What is the influence of selection effects on the incidence of companions as a function of type?

Hans Zinnecker:
At present we don’t know if there are close binaries with extreme mass ratios (not detectable as spectroscopic binaries) and we don’t know at all which sort of companions exist in the range of intermediate separations, too close to be seen as visual binaries in high-spatial resolution observations and too wide to show up as spectroscopic binaries. This gap in our knowledge may be bridged to some extent by future interferometric studies reaching a spatial resolution of the order of 1 mas (1 AU at 1 kpc distance) and brightness ratios 1 : 100. In summary, I believe that the census of massive stars with lower mass companions is still quite incomplete.

Cedric Foellmi:
There is an impressive number of scenarios of interacting massive binaries. We have undertaken in Montreal a systematic search for Wolf-Rayet binaries in the Magellanic Cloud. We found: ca. 40 % binary frequency in the SMC, 15 – 25 % in the LMC, no evidence for accretion in the secondaries, and intrinsic hydrogen in the wind of most of the (single and binary) WR stars. This argues against binary evolution for WR stars where it is theoretically the most expected (i.e. in a low Z environment).

Hans Zinnecker:
Your new data are very important constraints for stellar evolution and the close binary frequency of massive stars. However, they do not constrain the frequency of wide binaries. In my theoretical speculation, a lower Z environment tends to favour wider separations of binaries and likely faster rotation of the components, i.e. in my view a lower Z environment produces a higher angular momentum population of massive stars. As an aside, I never claimed that WR stars originate through the binary channel. In fact, this old concept probably was invented at a time when we had no idea about the importance of winds from massive stars. (Even at SMC/LMC metallicities stellar winds can do the trick to form WR stars, without a binary channel).

Claus Leitherer:
I do not see why I would make a large error in the interpretation of data of unresolved starburst populations if I assumed all stars were single. Unless stellar evolution changes (which may indeed be the case) very little should change.

Hans Zinnecker:
It all depends on the frequency of interacting close massive binaries which can lead to stellar mass exchange and even stellar mergers for the closest systems. Consequently some rejuvenation of a fraction of the massive stars can occur, as pointed out by Van Bever and Vanbeveren (1998), implying that binary evolution can make a starburst look younger than it really is, say 5 Myr instead of 10 Myr. A high binary frequency can make a difference here.

Anatol Cherepashchuk:
Ignoring binaries in starbursts will make a big difference for the presence of compact objects in post-starbursts, and the associated accretion effects (X-ray binaries).

Hans Zinnecker:
Indeed. Thank you, Mr. chairman. HZ shows the Chandra X-ray image of the Galactic Centre region with some 1000 point sources (Wang et al. 2002, Nature 415, 148)

## Acknowledgement

I would like to thank Karel van der Hucht for his invitation to speak at IAU-Symp. 212 in Lanzarote and for his patience with the manuscript. I acknowledge helpful information from Anthony Brown and Pavel Kroupa. I also owe a great deal of thanks to my massive stars co-investigators M. Bate, I. Bonnell, Th. Preibisch, G. Weigelt, and also my colleagues A. Maeder, F. Palla, and H. W. Yorke for many stimulating discussions which helped to shape my ideas on massive star formation. Finally, I thank the Brazilian participants of the conference to share with me the experience of the 2002 soccer world-cup final on June 30 (Brazil - Germany 2 : 0) in the hotel TV room in Lanzarote.

## References

Abt, H. A., Wang, R., & Cardona, O. 1991, ApJ 367, 155

Babel, J. & Montmerle, T. 1997, ApJ 485, L29

Bally, J. 2002, in ASP Conf. Ser. Vol. 267, The earliest stages of massive star birth, ed. P. Crowther (San Francisco: ASP), 219

Bate, M. R. & Bonnell, I. A. 1997, MNRAS 285, 33

Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS 336, 705

Beuther, H., Schilke, P., Sridharan, T. K. et al. 2002, A&A 383, 892

Bonnell, I. A., & Bate, M. R. 2002, MNRAS 336, 659

Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS 298, 93

Bosch, G., Morrell, N., & Niemela, V. 1999, Rev Mx A&A 35, 85

Bossi, M., Gaspani, A., Scardia, M., & Tadini, M. 1989, A&A 222, 117

Carpenter, J. M. 2000, AJ 120, 3139

Donati, J.-F., Babel, J., Harries, T. J. et al. 2002, MNRAS 333, 55

Duchene, G., Simon, T., Eisl ffel, J., & Bouvier, J. 2001, A&A 379, 147

Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS 172, 15

Garcia, B. & Mermilliod, J. C. 2001, A&A 368, 122

Gezari, D. Y., Backman, D. E., & Werner, M. W. 1998, ApJ 509, 283

Henning, Th. & Stecklum, B. 2002, in ASP Conf. Ser., Modes of star formation and the origin of field populations, ed. E. K. Grebel & W. Brandner, in press

Hunter, D. A. 1995, Rev Mx A&A 3, 1

Klessen, R. S. 2001, ApJ 550, L77

Kroupa, P. 2000, in Massive Stellar Clusters, ed. A. Lancon & C. Boily,
ASP Conf. Ser. 211, 233

Kroupa, P. 2002, Science 295, 82

Kurtz, S. & Franco, J. 2002, Rev Mx A&A 12, 16

Lequeux, J. 1979, A&A 71, L1

Maeder, A. & Behrend, R. 2002, in ASP Conf. Ser. Vol. 267, The earliest stages of massive star birth, ed. P. A. Crowther (San Francisco: ASP), 179

Mason, B. D., Gies, D. R., Hartkopf, W. I. et al. 1998, AJ 115, 821

Mermilliod, J. C. & Garcia, B. 2001, in IAU-Symp. 200, The Formation of Binary Stars, ed. H. Zinnecker & R. D. Mathieu (Dordrecht: Reidel), 191

Mouschovias, T. Ch. 1977 ApJ 211, 147

Padgett, D. L., Strom, S, E., & Ghez, A. 1997, ApJ 477, 705

Palla, F. & Stahler, S. W. 1993, ApJ 418, 414

Petr, M. G., Coude Du Foresto, V., Beckwith, St. V. W. et al. 1998, ApJ 500, 825

Preibisch, Th., Balega, Y., Hofmann, K.-H. et al. 1999, NewA 4, 531

Preibisch, Th., Weigelt, G., & Zinnecker, H. 2001, in IAU-Symp. 200, The Formation of Binary Stars, ed. H. Zinnecker & R. D. Mathieu (Dordrecht: Reidel), 69

Prosser Ch. F., Stauffer, J. R., Hartmann, L. et al. 1994, ApJ 421, 517

Reipurth, B. & Zinnecker, H. 1993, A&A, 278, 81

Scalo, J. 1998, in ASP Conf. Ser. Vol. 142, The Stellar Initial Mass Function,
ed. G. Gilmore & D. Howell (San Francisco: ASP), 201

Shepherd, D. S., Claussen, M. J., & Kurtz, S. E. 2001, Science 292, 1513

Silk, J. 1978, in Protostars and Planets I, ed. T. Gehrels (Tucson: Univ. of Arizona Press), 187

Simon, M., Close, L., & Beck, T. 1999, AJ 117, 1375

Stahl, O., Kaufer, A., Rivinius, T. et al. 1996, A&A 312, 539

Van Bever, J. & Vanbeveren, D. 1998, A&A 334, 21

Weigelt, G., Balega, Y., Preibisch, Th. et al. 1999, A&A 347, L15

Yorke, H. W. 2002, in ASP Conf. Ser. Vol. 267, The earliest stages of massive star birth, ed. P. Crowther (San Francisco: ASP), 165

Yorke, H. W. & Sonnhalter, C. 2002, ApJ 569, 846

Zinnecker, H. 1986, in IAU Symp. 116, Luminous stars and associations in Galaxies, ed. C. W. H. de Loore, A. J. Willis, & P. Laskarides (Dordrecht: Reidel), 271

Zinnecker, H., McCaughrean, M. J., & Wilking, B. M. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine (Tucson: Univ. of Arizona Press), 429

Zinnecker, H. & Bate, M. R. 2002, in ASP Conf. Ser. Vol. 267, The earliest stages of massive star birth, ed. P. Crowther (San Francisco: ASP), 209


text:
# The Dynamical Cluster Approximation (DCA) versus the Cellular Dynamical Mean Field Theory (CDMFT) in strongly correlated electrons systems

## Abstract

We are commenting on the article Phys. Rev. B 65, 155112 (2002) by G. Biroli and G. Kotliar in which they make a comparison between two cluster techniques, the Cellular Dynamical Mean Field Theory (CDMFT) and the Dynamical Cluster Approximation (DCA). Based upon an incorrect implementation of the DCA technique in their work, they conclude that the CDMFT is a faster converging technique than the DCA. We present the correct DCA prescription for the particular model Hamiltonian studied in their article and conclude that the DCA, once implemented correctly, is a faster converging technique for the quantities averaged over the cluster. We also refer to their latest response to our comment where they argue that instead of averaging over the cluster, local observables should be calculated in the bulk of the cluster which indeed makes them converge much faster in the CDMFT than in the DCA. We however show that in their original work, the authors themselves use the cluster averaged quantities to draw their conclusions in favor of using the CDMFT over the DCA.

In their article,biroli  G. Biroli and G. Kotliar compare two cluster methods, the Dynamical Cluster Approximation (DCA) and the Cellular Dynamical Mean Field Theory (CDMFT) for the simplified one-dimensional large-N modelaffleck  and conclude that the CDMFT converges faster with cluster size than the DCA.

This is a surprising result and in contradiction with exact results we have previously published in Ref. DCA\_maier  where we have shown that for large linear cluster sizes $`L_c`$, the DCA converges with corrections of $`𝒪(1/L_c^2)`$ while the CDMFT converges with corrections of $`𝒪(1/L_c)`$. The objective of this comment is to resolve this controversy.

The difference in scaling behaviors of the DCA and the CDMFT is a direct consequence of their different boundary conditions of the cluster. The DCA DCA\_hettler1 ; DCA\_hettler2  and the CDMFT kotliar  are both approximative methods incorporating nonlocal correlations in correlated electron systems relinquished by fully local techniques such as the Dynamical Mean Field Approximation (DMFA) metzner-vollhardt ; muller-hartmann ; review1 ; review2  and the Coherent Potential Approximation (CPA). taylor  The DCA is formulated in momentum space by dividing the first Brillouin zone into a number of subcells. It reduces the complexity of the problem by coarse-graining over all the momenta within a subcell. As a result, the real lattice is mapped onto a cluster with periodic boundary conditions coupled to a mean field. The CDMFT on the other hand is formulated in real space and maps the lattice onto a cluster with open boundary conditions.

We realized that the implementation of the DCA algorithm in Ref. biroli  is incorrect as they do not coarse-grain the nonlocal interaction and their analysis lacks an accurate scaling comparison in terms of the cluster size for these two methods. In addition, instead of presenting physically relevant lattice quantities they plot the corresponding cluster quantities which in the DCA in general are not identical to their lattice counterparts. We repeat the calculations of Ref. biroli  by implementing the DCA correctly and conclude that for any cluster size $`L_c`$, the DCA results are closer to the exact solution of the studied toy model than the CDMFT results, consistent with our previous findings. DCA\_maier

The exact mean field treatment of the simplified one-dimensional large-N model Hamiltonian (Eq. 4 in Ref. biroli ) yields

$$E_k=2(t+\chi )\mathrm{cos}k+\mu ,\chi =\frac{1}{L}\underset{k}{}f(\beta E_k)\mathrm{cos}k,$$
(1)
as self-consistent equations for the dispersion $`E_k`$, with the number of sites $`L`$, the chemical potential $`\mu `$, the inverse temperature $`\beta `$ and the Fermi function $`f(\beta E_k)`$. In the dispersion $`E_k`$, one can derive the self-energy as

$$\mathrm{\Sigma }(k)=2\chi \mathrm{cos}k,$$
(2)
where the $`\mathrm{cos}k`$ factor comes from the interaction between nearest neighbors.

In the DCA, the first Brillouin zone is divided into $`L_c`$ subcells with a linear size of $`\mathrm{\Delta }k=2\pi /L_c`$. The irreducible quantities are constructed from the coarse-grained propagators. Both the Green function, and the interaction (when it is non-local) must be coarse grainedDCA\_hettler2 ; aryanpour . Thus, the self-energy part in Eq. 2 changes to its DCA counterpart

$$\mathrm{\Sigma }_{DCA}(K_c)2\chi _{cl}^{DCA}\overline{\mathrm{cos}(K_c)},$$
(3)
where

$$\overline{\mathrm{cos}(K_c)}=\frac{L_c}{L}\underset{\stackrel{~}{k}}{}\mathrm{cos}(K_c+\stackrel{~}{k}),$$
(4)
comes from coarse-graining the nearest neighbor interaction over all the momenta $`\stackrel{~}{k}`$ within each subcell around the cluster momenta $`K_c`$ representing the subcells. G. Biroli and G. Kotliar in Ref. biroli  (Eq. 6 for the dispersion) do not coarse-grain $`\mathrm{cos}k`$ in the nonlocal interaction and instead, they just substitute it by $`\mathrm{cos}K_c`$. This implementation of the DCA is incorrect and overestimates the effects of the interaction. The correct self-consistent DCA equations for the dispersion are

$`\chi _{cl}^{DCA}={\displaystyle \frac{1}{L_c}}{\displaystyle \underset{𝐊_c}{}}\mathrm{cos}(K_c)\times \overline{f(\beta E(K_c))},`$ (5)
$`\overline{f(\beta E(K_c))}={\displaystyle \frac{L_c}{L}}{\displaystyle \underset{\stackrel{~}{k}}{}}f(\beta E_{K_c+\stackrel{~}{k}}),`$
$`E_{K_c+\stackrel{~}{k}}2t\mathrm{cos}(K_c+\stackrel{~}{k})2\chi _{cl}^{DCA}\overline{\mathrm{cos}(K_c)}+\mu .`$
After self-consistency is achieved, one can use the cluster quantity $`\chi _{cl}^{DCA}`$ to compute its lattice counterpart

$$\chi _{latt}^{DCA}=\frac{1}{L}\underset{k}{}f(\beta E_{K_c+\stackrel{~}{k}})\mathrm{cos}k.$$
(6)
G. Biroli and G. Kotliar plot the cluster quantity $`\chi _{cl}^{DCA}`$ instead of the physically relevant lattice quantity $`\chi _{latt}^{DCA}`$ (Fig. 1 in Ref. biroli ). In the DCA, non-local cluster quantities like $`\chi _{cl}^{DCA}`$ are not identical to their corresponding lattice quantities and thus have no physical meaning. In the CDMFT, since there exists one to one correspondence between the lattice quantity $`\chi _{latt}^{CDMFT}`$ and the cluster quantity $`\chi _{cl}^{CDMFT}`$, they are identical.

By numerically solving the self-consistent Eqs. 1 and 5 and the one for the CDMFT (Eq. 8 in Ref. biroli ), we can make a comparison between the exact solution (cf. Eq. 1 in this article) and the predictions by the DCA and the CDMFT.

In Fig. 1, we plot $`\chi `$ versus $`\beta `$ for $`t=1`$ and fixed filling of $`n=0.37`$ ($`\mu 1`$) using the correct implementation of the DCA (compared to the results of Fig. 1 in Ref. biroli ) and CDMFT for $`L_c=2,4,8`$ and $`16`$. It is more physical to do the calculations at fixed filling $`n`$ instead of fixed chemical potential $`\mu `$ as it was done in Ref. biroli . As illustrated, it is clear that the DCA results converge faster than the CDMFT ones to the exact solution even for small $`L_c`$. The CDMFT results converge consistently whereas for small $`L_c`$ ($`L_c=8`$ in particular) the DCA results fluctuate around the exact solution. Nevertheless, at large enough $`L_c`$, the DCA behaves consistently and converges from below. We have avoided odd cluster sizes in these calculations because for the DCA, the periodic boundary conditions on the cluster suppress $`Q=\pi `$ fluctuations by frustrating the cluster.

In Fig. 2, we plot $`\chi `$ as a function of the inverse cluster size $`1/L_c`$ computed by the DCA and the CDMFT at fixed inverse temperature $`\beta =5`$ for $`t=1`$ and $`n=0.37`$ ($`\mu 1`$) and compare the results with the exact solution (solid line). The DCA begins to converge quadratically and the CDMFT linearly in $`1/L_c`$ to the exact solution beyond $`L_c=32`$ and $`L_c=10`$ respectively (the fits in the insets). The scaling behavior of both the DCA and CDMFT at large enough $`L_c`$ are consistent with our previous results. DCA\_maier

It is noteworthy that for this Hamiltonian in particular, the values of $`L_c`$ beyond which the scaling behaviors of the DCA and the CDMFT are observed are relatively high. E.g., in the one dimensional Falicov-Kimball model at half filling, these scaling behaviors are manifest beyond $`L_c=2`$. DCA\_maier

The results in both Figs. 1 and 2 indicate the superiority of the DCA over the CDMFT in terms of faster convergence with the cluster size $`L_c`$. The DCA also significantly reduces the computational CPU time compared to the CDMFT. In the CDMFT, because of the open boundary conditions on the cluster, the Green function in both real and momentum spaces is not diagonal. Thus, expensive diagonalizations of large matrices are required for large cluster sizes.

In their latest response to this comment, biroli.resp  G. Biroli and G. Kotliar argue that in the CDMFT, whenever the correlation length is finite, local observables such as $`\chi `$ in this problem are best estimated in the bulk than near the boundaries of the cluster. In fact, their values in the bulk converge exponentially in the CDMFT while in the DCA they still converge as $`1/L_c^2`$. Because the CDMFT breaks the translational invariance inside the cluster, local observables calculated in the bulk of the cluster differ from those on the boundary. Provided that the system is far from a transition, the sites in the center of the CDMFT cluster couple to the mean-field only through propagators which fall exponentially with distance. Local observables, when extracted from the center of the cluster instead of from their flat average, hence converge much faster, i.e. exponentially as opposed to $`1/L_c`$. As a matter of fact, this point was made previously by some of us maier:dca4 ; however, for most models, it is difficult to formulate a causal self energy from quantities measured on the center of the cluster. Nevertheless, in their original work in Fig.1 of Ref. biroli , G. Biroli and G. Kotliar themselves perform a flat averaging of the local observable $`\chi `$ and by comparing its results with their incorrect DCA results, they conclude that the CDMFT converges better. To prove that the results in Fig.1 in Ref. biroli  are indeed the flat averaged $`\chi `$ and not its value in the bulk, we plot in Fig. 3 both the flat averaged and bulk $`\chi `$ versus $`\beta `$ for $`L_c=2,4`$. We have also extracted the CDMFT curves in Fig.1 of Ref. biroli  for $`L_c=2,4`$ (filled circles and squares respectively) using Engauge digitizing package. Their perfect match with the flat averaged curves and their deviation from the bulk result for $`L_c=4`$ best indicate that $`\chi `$ in Fig.1 of Ref. biroli  was flat averaged and not calculated in the bulk.

Another significant difference between these two techniques appears in the calculation of the self-energy. The DCA approximates the lattice self-energy by a constant within a DCA subcell in momentum space and as a result, the self-energy becomes a step function in $`k`$-space. Thus, one may use an approach such as the Akima spline which is the smoothest curve interpolating through these step functions to generate an analytic form for the self-energy. Using the Akima spline is consistent with the assumption that the self-energy is a smoothly varying function in momentum space. In the CDMFT, since the lattice self-energy is approximated by the cluster self-energy in real space, its Fourier transform to momentum space yields a smooth self-energy curve. Due to the open boundary conditions of the cluster, the translational invariance is violated and the self-energy in momentum space is not diagonal. Nevertheless, by using the Random Phase Approximation (RPA), senechal  the translational invariance can be restored and the self-energy becomes a diagonal function in momentum space.

###### Acknowledgements.

We would like to acknowledge M. H. Hettler and H. R. Krishnamurthy for the review of this comment and useful suggestions. This work was supported by NSF grants DMR-0073308 and DMR-9704021.


text:
# Triggers and alerts with GLAST

## 1. Description of the instruments

GLAST is a next generation high-energy gamma ray observatory, onboard a satellite scheduled for launch in 2006. It consists of a Large Area Telescope (LAT, see Brez, S. et al. 2001) and a GLAST Burst Monitor (GBM, see Kippen et al. 2001), and is designed for making observations of gamma-ray sources in the 10 keV to 300 GeV energy range.

The LAT is a pair conversion telescope, operating from 10 MeV to more than 300 GeV. It is composed of three subsystems: the Anti-Coincidence Detector (ACD) is responsible for vetoing the enormous charged cosmic-ray background and the Earth albedo secondary electrons and nuclei. The tracker (TKR) consists of a four-by-four array of tower modules, each of which being composed of 19 pairs of interleaved planes of silicon strip detectors and tungsten converter sheets. The silicon strip detectors track the electron positron pair created by the incident gamma. The calorimeter (CAL) is a segmented arrangement of CsI(Tl) bars, designed to give both longitudinal and transverse information on the energy deposition pattern.

The tracker will have the ability to determine the location of an object in the sky to within 0.5 to 5 arc minutes. The pair conversion signature is also used to help reject the background of charged cosmic rays. The GBM provides overlapping energy coverage in the range 10 keV to 25 MeV for bright transients such as bursts and solar flares. It is composed of 4 triplets of NaI(Tl) scintillators(low energy band) and 2 BGO scintillators (high energy band).

Due to its angular resolution (5 times better than EGRET) and its improved afterglow sensitivity, GLAST will extend the knowledge about GRB’s. Its large field of view (about 2 sr) and its improved effective area (10 times EGRET’s) will provide better burst statistics (up to 100 bursts/year with many more photons detected.)

The GBM will enhance GRB science with GLAST by providing low energy context measurements with high time resolution. This will improve GLAST’s wide band spectral sensitivity, allow for comparison of low energy versus high energy temporal variability, and provide a continuity with current GRB knowledge base (e.g GRO-BATSE.) Moreover, the GBM provides fast time/location triggers: within 2 seconds, it will provide approximate locations for ground or space follow up observations. In case of specifically interesting events, it will allow, within ten minutes, for repointing of GLAST for afterglow observations.

## 2. GLAST and GRB Simulations

The LAT team has set up a complete simulation chain, including generation of the incoming flux, full simulation of the detector resonse, reconstruction chain and analysis of the final trigger and alert. It is to be noted that the flux simulations include background (from albedo, cosmic and diffuse gamma ray events) and two possible ways of signal generation: the first based on a complete physical model, the second based on a phenomenological extrapolation from BATSE data. Orbit, tilting and correct exposure are also taken into account.

The physical model has been validated in the following way: the BATSE observed distribution of duration has been fitted with a double gaussian, for short and long bursts (see Figure 1, left panel). Then, we used this bimodal distribution as input parameter for our model, which explores the parameter space and returns a burst with a given duration. The resulting simulated catalog is depicted in Figure 1, right panel. Finally, the simulated fluences are compared with the observed ones, for the 4 BATSE channels (Figure 2). In this way we have a calibration method to correlate the emission at high energy with the low energy observations. Scaling laws are also reproduced.

## 3. Trigger and Alerts

After calibrating the model, we investigated the efficiency of GLAST in triggering and reconstructing photons coming from GRB events. In this preliminary study we have tested various efficiencies for a sample of gamma ray bursts without background. A set of GRB’s with different fluences above 100 MeV has been generated, and the resulting signal has been processed through the LAT instrument simulator. Figure 3 (left) shows the number of incident, triggered and reconstructed gammas as a function of fluence above 100 MeV. The general trend is that for a typical burst of $`10^8`$ erg/cm2 we expect approximately 100 incident photons, and a photon trigger efficiency of 30%. More than 90% of the triggered photons will then be reconstructed.

For the purpose of this study, we use a ’strawman’ real time unbinned trigger algorithm: the sky is scanned with a sliding 20 event window, looking for clusters of events in a 35 degree circle; then a joint spatial-temporal likelihood is formed. The threshold for this trial situation was set adopting an arbitrarily low background event rate of 3 Hz, to produce less than one false trigger per week.

The preliminary results of this trigger simulation are presented in Figure 3. For the quoted background event rate, efficiency reached about 85% for the BASTE-like bursts, with around 78% of those bursts detected by GLAST within one second of GRB onset.

## 4. Conclusion

We have investigated the performances of GLAST in detecting GRB events, using a complete simulation chain. From our preliminary results, it appears that GLAST should be able to substantially improve on the EGRET detection rate of about one GRB per year. In addition to burst alerts and localizations provided by the GBM, the LAT should also be able to provide real time triggers and GRB localizations.

Refinements to the trigger simulations will include simulations at higher background event rates, spatial dependence of diffuse gamma flux, and temporal dependence of the on board background rates.

### Acknowledgments.

J.D Scargle acknowledges support from NASA’s Applied Information Systems Research Program.

## References

Brez, S. et al. 2001, in Proceedings of the 7th International Conference on Advanced Technology and Particle Physiscs, (World Scientific)

Kippen, R. M. et al. 2001, in AIP Conf. Proc., 587, Gamma-Ray Astrophysics 2001, eds. S. Ritz, N. Gehrels, & C. R. Shrader, 801


text:
# Effects of parametric noise on a nonlinear oscillator

## 1 Introduction

A simple model of a nonlinear stochastic system is obtained by including a non-quadratic confining potential in the classical Langevin equation . In this work, we address the problem of a nonlinear oscillator with random linear frequency. Due to the continuous injection of energy by the random force, observables such as the oscillator’s mechanical energy, position and velocity, grow algebraically with time in the absence of dissipation. In Section 2, we calculate the associated growth exponents in the case of Gaussian white noise, derive the long time asymptotic behavior of the probability distribution function (P.D.F.) in phase space and match it with the (non-Gibbsean) stationary P.D.F. in presence of small dissipation. In Section 3, we prove that if the noise has a non-vanishing correlation time, the dynamics of the nonlinear oscillator is drastically modified: colored noise yields anomalous diffusion exponents equal to half the values found for white noise.

## 2 White multiplicative noise

We consider a nonlinear oscillator of amplitude $`x(t)`$, submitted to parametric noise and confined by a potential $`𝒰(x)`$. For $`|x|\mathrm{}`$, a suitable rescaling permits us to write $`𝒰\frac{x^{2n}}{2n}`$ with $`n1.`$ Thus, for large amplitudes, the stochastic differential equation that governs the system becomes

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+\gamma \frac{\text{d}}{\text{d}t}x(t)+x(t)^{2n1}=\xi (t)x,\text{ with }\xi (t)\xi (t^{})=𝒟\delta (tt^{}),$$
(1)
where $`\gamma `$ is the damping rate and $`\xi (t)`$ is a Gaussian white noise of zero mean-value and of autocorrelation $`𝒟`$. This equation is interpreted according to the rules of Stratonovich Calculus . For a quadratic potential (i.e. $`n=1`$) and for small $`\gamma `$, the energy of the oscillator grows exponentially with time . The effect of nonlinear restoring forces must be taken into account to avoid this exponential amplification.

We shall analyse the motion of the nonlinear stochastic oscillator following the method explained in . Defining the energy and the angle variables as

$$E=\frac{1}{2}\dot{x}^2+\frac{1}{2n}x^{2n},\text{ and }\varphi =\frac{\sqrt{n}}{(2n)^{1/2n}}_0^{x/E^{1/2n}}\frac{\text{d}u}{\sqrt{1\frac{u^{2n}}{2n}}},$$
(2)
we transform the coordinates in phase space from position and velocity to energy and angle. Introducing an auxiliary variable $`\mathrm{\Omega }`$, with $`\mathrm{\Omega }=(2n)^{\frac{n+1}{2n}}E^{\frac{n1}{2n}},`$ we rewrite Eq. (1) as

$`\dot{\mathrm{\Omega }}`$ $`=`$ $`\gamma {\displaystyle \frac{n1}{(2n)^{\frac{1}{n}}}}𝒮_n^{}(\varphi )^2\mathrm{\Omega }+(n1)𝒮_n(\varphi )𝒮_n^{}(\varphi )\xi (t),`$ (3)
$`\dot{\varphi }`$ $`=`$ $`\gamma {\displaystyle \frac{𝒮_n(\varphi )𝒮_n^{}(\varphi )}{(2n)^{\frac{1}{n}}}}+{\displaystyle \frac{\mathrm{\Omega }}{(2n)^{\frac{1}{n}}}}{\displaystyle \frac{𝒮_n(\varphi )^2}{\mathrm{\Omega }}}\xi (t),`$ (4)
where the hyperelliptic function $`𝒮_n`$ satisfies $`𝒮_n(\varphi )=x/E^{1/(2n)}`$. These equations are rigorous and have been derived without any hypothesis on the noise term. In the small damping limit $`\gamma 0`$, we deduce from Eqs. (3,4) that $`\mathrm{\Omega }t`$ and $`\varphi t^2`$, i.e. $`\varphi `$ is a fast variable as compared to $`\mathrm{\Omega }`$. Averaging the dynamics over $`\varphi `$ yields the statistical equipartition identities :

$$E=\frac{n+1}{2n}\dot{x}^2\mathrm{and}x^{2n}=\dot{x}^2.$$
(5)
This averaging procedure  allows us to derive a closed equation for the stochastic evolution of the slow variable $`\mathrm{\Omega }.`$ We start with the Fokker-Planck equation for the P.D.F. $`P_t(\mathrm{\Omega },\varphi )`$ associated with the system (3, 4), and average it under the hypothesis that the probability density becomes uniform in $`\varphi `$ when $`t\mathrm{}`$. Hence, the reduced probability density $`\stackrel{~}{P}_t(\mathrm{\Omega })`$ satisfies

$`_t\stackrel{~}{P}=\gamma {\displaystyle \frac{n1}{n+1}}_\mathrm{\Omega }\left(\mathrm{\Omega }\stackrel{~}{P}\right)+{\displaystyle \frac{\stackrel{~}{𝒟}}{2}}\left(_\mathrm{\Omega }^2\stackrel{~}{P}{\displaystyle \frac{2}{n1}}_\mathrm{\Omega }{\displaystyle \frac{\stackrel{~}{P}}{\mathrm{\Omega }}}\right).`$ (6)
This averaged Fokker-Planck equation is associated with the following effective Langevin dynamics for the variable $`\mathrm{\Omega }`$,

$$\dot{\mathrm{\Omega }}=\frac{\stackrel{~}{𝒟}}{n1}\frac{1}{\mathrm{\Omega }}\gamma \frac{n1}{n+1}\mathrm{\Omega }+\stackrel{~}{\xi }(t)$$
(7)
where the effective Gaussian white noise $`\stackrel{~}{\xi }`$ has an amplitude $`\stackrel{~}{𝒟}`$ given by

$$\stackrel{~}{𝒟}=𝒟(2n)^{\frac{2}{n}}\frac{(n1)^2}{n+1}\frac{\mathrm{\Gamma }\left(\frac{3}{2n}\right)\mathrm{\Gamma }\left(\frac{3n+1}{2n}\right)}{\mathrm{\Gamma }\left(\frac{1}{2n}\right)\mathrm{\Gamma }\left(\frac{3n+3}{2n}\right)}.$$
(8)
If the dissipation rate $`\gamma `$ is taken to be zero, the solution of Eq. (6) is :

$$\stackrel{~}{P}_t(E)=\frac{1}{\mathrm{\Gamma }\left(\frac{n+1}{2(n1)}\right)}\frac{n1}{nE}\left(\frac{(2n)^{\frac{n+1}{n}}E^{\frac{n1}{n}}}{2\stackrel{~}{𝒟}t}\right)^{\frac{n+1}{2(n1)}}\mathrm{exp}\left\{\frac{(2n)^{\frac{n+1}{n}}E^{\frac{n1}{n}}}{2\stackrel{~}{𝒟}t}\right\}.$$
(9)
From this P.D.F., we derive analytical expressions for statistical averages of physical observables . The leading scaling behaviour is given by

$$E(𝒟t)^{\frac{n}{n1}},x(𝒟t)^{\frac{1}{2(n1)}},\dot{x}(𝒟t)^{\frac{n}{2(n1)}}.$$
(10)
These formulae are in excellent agreement with numerical simulations . We show in Fig. 1 that the averaged energy scales like $`t^{\frac{n}{n1}}`$ as predicted. The amplitude $`x`$ undergoes an anomalous diffusion with an exponent that diverges as $`n1`$, which is consistent with the exponential growth of the linear oscillator.

In the presence of dissipation, the stationary solution of the Fokker-Planck equation (6) is

$$\stackrel{~}{P}_{\mathrm{st}}(E)=\frac{1}{\mathrm{\Gamma }\left(\frac{n+1}{2(n1)}\right)}\frac{n1}{nE}\left(\frac{(2n)^{\frac{n+1}{n}}E^{\frac{n1}{n}}}{\frac{n+1}{\gamma (n1)}\stackrel{~}{𝒟}}\right)^{\frac{n+1}{2(n1)}}\mathrm{exp}\left\{\frac{(2n)^{\frac{n+1}{n}}E^{\frac{n1}{n}}}{\frac{n+1}{\gamma (n1)}\stackrel{~}{𝒟}}\right\}$$
(11)
This measure is not  the canonical Boltzmann-Gibbs distribution. The crossover time $`t_c`$ from the asymptotic distribution function for the undamped oscillator (9) to the stationary measure in presence of damping (11) is $`t_c=\frac{1}{2\gamma }\frac{n+1}{n1}.`$ Thus, Eq. (9), although derived for a non-dissipative system, is physically relevant: it provides an intermediate time asymptotic behavior  for the P.D.F. before the stationary state is reached; the energy of the system is distributed according to the time-dependent P.D.F. (9) for times $`t`$ such that $`1tt_c`$ and then follows the stationary law (11) when $`tt_c`$.

## 3 Colored multiplicative noise

We now consider a system without damping and submitted to a multiplicative noise with non-zero correlation time $`\tau `$. In , we infered from dimensional analysis that the white noise scalings are modified. Here, we shall prove this conjecture by calculating the anomalous exponents from a self-consistent approach.

We study the equation:

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+x(t)^{2n1}=\eta (t)x(t),$$
(12)
where $`\eta (t)`$ is a colored, Ornstein-Uhlenbeck, Gaussian noise, obeying

$$\frac{\text{d}\eta (t)}{\text{d}t}=\frac{1}{\tau }\eta (t)+\frac{1}{\tau }\xi (t).$$
(13)
When $`t,t^{}\tau `$, we find:

$$\eta (t)=0\mathrm{and}\eta (t)\eta (t^{})=\frac{𝒟}{2\tau }\mathrm{e}^{|tt^{}|/\tau }.$$
(14)
As stated earlier, Eqs. (3) and (4) remain valid if $`\xi (t)`$ is replaced by $`\eta (t)`$. We assume that, in the long time limit, $`\mathrm{\Omega }`$ grows algebraically with time with a scaling exponent $`\alpha `$, i.e., $`\mathrm{\Omega }t^\alpha .`$ We deduce from Eq. (4) that $`\varphi t^\nu `$ with $`\nu =\alpha +1`$. Substituting this scaling of $`\varphi `$ in Eq. (3), we obtain (leaving aside all proportionality constants):

$$\mathrm{\Omega }_0^t\mathrm{dz}𝒮_\mathrm{n}(\mathrm{z}^\nu )𝒮_\mathrm{n}^{}(\mathrm{z}^\nu )\eta (\mathrm{z}).$$
(15)
The asymptotic behavior of this expression can be estimated as in  by discretizing time. This simplification does not alter the critical exponents but modifies only the prefactors. We thus replace the colored Gaussian noise $`\eta `$ by a discrete random variable $`ϵ_k`$ which takes the values $`\pm \sqrt{𝒟/(2\tau )}`$ randomly during the time interval $`[k\tau ,(k+1)\tau ].`$ From Eq. (15), we deduce that

$$\mathrm{\Omega }^2\underset{k=0}{\overset{t/\tau }{}}\left(_{k\tau }^{(k+1)\tau }\mathrm{dz}𝒮_\mathrm{n}(\mathrm{z}^\nu )𝒮_\mathrm{n}^{}(\mathrm{z}^\nu )\right)^2\underset{\mathrm{k}=1}{\overset{\mathrm{t}/\tau }{}}\frac{1}{\left(\mathrm{k}\tau \right)^{2\nu 2}}\mathrm{t}^{32\nu },$$
(16)
where the integral is evaluated by integrating by parts and retaining only the leading terms . Assuming that the variable $`\mathrm{\Omega }`$ is not multifractal, we also have $`\mathrm{\Omega }^2t^{2\alpha }.`$ We then deduce that $`32\nu =2\alpha `$. Since $`\nu =\alpha +1`$, we obtain $`\alpha =1/4`$ and $`\nu =5/4`$. Finally, the scaling exponents for colored noise are

$$Et^{\frac{n}{2(n1)}},xt^{\frac{1}{4(n1)}},\dot{x}t^{\frac{n}{4(n1)}}.$$
(17)
These predictions are confirmed by numerical simulations (see in Fig. 2 the scaling behavior of the averaged energy).

We observe that the exponents are halved when the noise is colored. Further, we expect that the colored noise scaling will be observed for large enough time as soon as the correlation time $`\tau `$ is nonzero. This result may seem paradoxical at first sight because it is often wrongly thought that “in the long time limit, colored noise appears white”. In fact the period $`T`$ of the underlying deterministic oscillator must be compared with the correlation time $`\tau `$ of the noise. If $`\tau T,`$ the noise is uncorrelated over a period of the deterministic oscillator and acts as if it were white: the scalings found in Eq. (10) are satisfied. But, for $`T\tau `$, the noise remains coherent over a period and its effect is perceptible only when it is considered over a large number of periods and hence the diffusion slows down: the scalings of Eq. (17) now apply. The crossover between the two regimes is observed when $`T\tau `$. The deterministic period $`T`$ decreases with the energy as $`TE^{\frac{n1}{2n}}`$ . From Eq.(10), we have $`E^{\frac{n1}{2n}}(𝒟t)^{1/2}`$. Therefore, the crossover time scales as $`t_c(𝒟\tau ^2)^1`$. As shown in Figure 3, numerical simulations confirm this scenario. Such a crucial relevance of the correlation time of the noise was also found in the study of the bifurcation threshold of the Duffing oscillator with multiplicative noise .

## 4 Conclusion

A particle trapped in a confining potential and subject to multiplicative noise undergoes anomalous diffusion before the dissipative time scale is reached. The critical exponents associated with this behaviour are calculated exactly in the case of Gaussian white noise and an analytical expression for the intermediate time asymptotics of the P.D.F. in phase space is derived.

Our approach relies on the integrability of the underlying deterministic system: we derive exact stochastic equations in energy-angle variables and, after averaging out the phase variations, a projected dynamics for the energy variable is obtained. Our analytical predictions agree perfectly with numerical simulations. For colored multiplicative noise, the anomalous diffusion exponents have been calculated in a self-consistent manner. We find that these exponents are halved in presence of time correlations, however small the correlation time $`\tau `$ may be.

We observe a crossover from the white noise regime at short times to the colored noise scalings at long times. The associated crossover time scales like $`\tau ^2`$ and corresponds to the regime when the period of the associated deterministic oscillator is of the order of $`\tau `$. In the colored noise case, however, the averaging technique does not lead to conclusive results because the noise itself is averaged out to the leading order. A precise calculation in the case of colored Gaussian noise still remains to be done.