text:
# The Turbulent Star Formation Model. Outline and Tests

## 1. Introduction

Two key questions related to star formation are a) What is the origin and nature of star-forming (and non-star-forming, if they exist) cores in turbulent molecular clouds? b) The star formation efficiency (SFE), measured as the fraction of a molecular cloud’s mass that ends up in stars during its lifetime, is low, of order 5-10% (see, e.g., the review by Evans 1999). Why? Different answers to these questions are given in the main two competing models of star formation. In the so-called “standard” model (Mouschovias 1976; Shu 1977; see the review by Shu, Adams & Lizano 1987), low-mass star-forming cores begin their lives as larger, magnetostatic, magnetically subcritical, partially-ionized clumps of low density contrast relative to their parent cloud. These slowly lose magnetic support and contract as the neutrals “slip” through the ions in a process commonly referred to as “ambipolar diffusion”, with characteristic time scale $`\tau _{\mathrm{ad}}`$, until they finally become magnetically supercritical and proceed to gravitational collapse. The low SFE is explained in this case because $`\tau _{\mathrm{ad}}`$ is much larger than the free-fall time scale $`\tau _{\mathrm{ff}}`$, and thus constitutes the effective time scale for star formation, rather than $`\tau _{\mathrm{ff}}`$. High-mass star formation is assumed to occur through the collapse of clumps that agglomerate enough mass to become supercritical and thus collapse on a time scale $`\tau _{\mathrm{ff}}`$, although clumps this massive are assumed to be rare. However, this model suffers from a number of shortcomings (see, e.g., Mac Low & Klessen 2003; Hartmann, this conference). In particular, it does not address the formation of the clumps, which are simply taken as intial conditions, nor how can such magnetostatic structures survive in the turbulent environment of their parent clouds.

In the competing, more recent turbulent scenario (e.g., Elmegreen 1993; Padoan 1995; Ballesteros-Paredes, Vázquez-Semadeni & Scalo 1999a, Ballesteros-Paredes, Hartmann & Vázquez-Semadeni 1999b; Klessen, Heitsch & Mac Low 2000; Heitsch, Mac Low & Klessen 2001; Hartmann, Ballesteros-Paredes& Bergin 2001; Padoan & Nordlund 2002; Vázquez-Semadeni, Ballesteros-Paredes & Klessen 2003; see also the reviews by Vázquez-Semadeni et al. 2000; Mac Low & Klessen 2003) the answer to question (a) above is that the cores are turbulent density fluctuations in the molecular clouds (which have typical rms Mach numbers $`M_\mathrm{s}10`$) and that their statistical properties are therefore determined by the turbulent parameters in the parent cloud. The clouds are globally supported by turbulence, and in fact may be transient with short lifetimes (Ballesteros-Paredes et al. 1999b; Hartmann et al. 2001), but locally the density fluctuations induced by the turbulence may or may not collapse depending on their particular energy balance (Ballesteros-Paredes & Vázquez-Semadeni 1997; Vázquez-Semadeni, Shadmehri & Ballesteros-Paredes 2002). Then, the answer to question (b) above in the turbulent model is that the SFE is low because the cores contain only a fraction of the total cloud mass, only a fraction of them undergoes collapse, and, as is also the case in the standard model, only a fraction of the cores’ mass is involved in the collapse. That is, only a small fraction of the globally-turbulence-supported cloud is able to undergo collapse, through a kind of “turbulent colander”. Note that, while magnetic support is a key concept in the standard model of star formation, in the turbulent model it is not indispensible, as delay of the collapse is not necessary.

In this paper we first summarize the current status of the turbulent model of star formation, together with key work that has led to its formulation (§2), and then we discuss evidence that indeed turbulence can produce structures which reproduce some important observational properties of “quiescent” molecular cloud cores (§3), thus implying a high plausibility for the turbulent model. In this paper, we assume that the turbulence in molecular clouds is continuously driven, as otherwise it would decay in roughly one crossing time (Mac Low et al. 1998; Stone et al. 1998; Padoan & Nordlund 1999). This is in fact expected if molecular clouds are an intermediate-scale part of the turbulent cascade in the ISM, since energy is then continuously cascading through them from the large energy injection scales to the small dissipative ones.

## 2. The turbulent model of core and star formation

It was already forseen by von Weizsäcker (1951) and Sasao (1973) that compressible turbulence in gas clouds must have a dual role: turbulent modes of size scale $`\lambda `$ provide “pressure” towards structures with sizes $`>\lambda `$, while simultaneously they induce the formation of large-amplitude density fluctuations of sizes $`<\lambda `$. The latter process is generally called “turbulent fragmentation”. Note that, in order to be highly compressible, the turbulence must be supersonic, as is observed (see, e.g., Zuckerman & Palmer 1974), and therefore it overwhelms the thermal support at all but the smallest scales. This implies that turbulence-supported clouds may have masses much larger than their thermal Jeans mass.

Sasao (1973) already pointed out that the compressed regions may be driven to collapse, and Hunter & Fleck (1982) showed that within them the local Jeans length is strongly reduced. This implies that local collapsing cores generally may have masses smaller than even the thermal Jeans mass at the mean cloud density, and so their masses are in general much smaller than their parent cloud’s mass. Note additionally that the free fall time in the cores is given by $`\tau _{\mathrm{ff}}(G\rho _{\mathrm{core}})^{1/2}`$, and is therefore shorter than that of their parent cloud by a factor $`\sqrt{\rho _{\mathrm{cloud}}/\rho _{\mathrm{core}}}`$. This is the opposite to the case of linear (small-amplitude) perturbations, in which the fastest collapsing scale is the largest one (Larson 1985).

Moreover, not all the density fluctuations necessarily proceed to collapse; whether they do or do not depends on whether they acquire enough mass that their gravitational energy overwhelms all other forms of energy that provide support. If it does not, the fluctuations re-expand, as first shown numerically by Vázquez-Semadeni, Passot & Pouquet (1996) using non-isothermal, non-magnetic simulations subject to turbulent driving. In these, many generations of density fluctuations (cores) were seen to appear and disappear, until finally a (rather improbable) strong enough compressive event managed to produce a gravitationally unstable core, which then proceeded to collapse in roughly a local free-fall time. Moreover, the global density maximum in the simulations was seen to fluctuate chaotically in time (showing also the non-formation of any hydrostatic cores), until it suddenly began to increase at an accelerated pace, corresponding to the local, almost instantaneous collapse of a “core”.

Padoan (1995; see also Padoan & Nordlund 2002; Vázquez-Semadeni et al. 2003) proposed that the only fragments that collapse are those which are far enough down the turbulent cascade that their internal velocity dispersion is subsonic,111Note that the existence of a subsonic range is simply the consequence of a turbulent cascade, and does not necessarily imply that subsonic cores are at the scale of turbulence dissipation. and that in addition are gravitationally unstable. This scenario has the implication that, if one considers a succession of turbulent regimes in which subsonic velocity dispersions occur at progressively smaller scales, the fraction of the mass available for collapse, and therefore the SFE, should also be progressively smaller, because smaller regions contain smaller fractions of the total mass.

This suggestion was tested numerically by Vázquez-Semadeni et al. (2003). For a series of SPH simulations of turbulent clouds, they empirically measured the “sonic scale” $`\lambda _\mathrm{s}`$, i.e., the scale at which the turbulent velocity dispersion equals the sound speed. Then they also operationally defined and measured the SFE as the fraction of mass in collapsed sink particles after two turbulent crossing times. With these data, they showed that the SFE correlates well with $`\lambda _\mathrm{s}`$, and in fact scales as SFE $`\mathrm{exp}(\lambda _0/\lambda _\mathrm{s})`$, where $`\lambda _0`$ is a suitable reference scale. This suggests that indeed the sonic scale is one fundamental parameter in determining the SFE. Additionally, Vázquez-Semadeni et al. (2003) also showed that the SFE does not correlate well with the driving scale, and that previous results suggesting so (Klessen et al. 2000) were only an artifact of maintaining the same total kinetic energy in runs with different driving scales, as this caused higher energies per unit wavenumber interval in runs forced at smaller scales, decreasing the sonic scale.

However, it is also clear that not only subsonic cores can collapse (if they are super-Jeans), but also supersonic cores, if gravity anyway overwhelms all other forces. In this case, not contemplated by either Padoan (1995), Padoan & Nordlund (2002) nor Vázquez-Semadeni et al. (2003), turbulence will be able to still produce fragments, but not support the cloud as a whole. Given the shorter free-fall time of the smaller, nonlinear scales, these will collapse first. As soon as they do, they will collectively become a non-dissipative system, and thus will not collapse to a singularity, but instead form a bound cluster of collapsed objects (see Bate, this volume). Thus, one can tentatively associate the global collapse of an unsupported supersonic region with the “clustered” mode of star formation, while the individual collapse of subsonic regions may be associated with the “isolated” mode, in which each core may produce one or a few stars.

Finally, there is the issue of the cores that do not collapse, but instead rebound after the compression is over, and merge back with their surronding medium. Vázquez-Semadeni et al. (2002) have performed a simple calculation to estimate the re-expansion time $`\tau _{\mathrm{exp}}`$, defined as the time it takes to double the initial radius. This is shown in fig. 1 as a function of the core radius at maximum compression, with $`\tau _{\mathrm{exp}}`$ given in units of the free-fall time, and the radius in units of the equilibrium radius at which internal pressure balances self-gravity. We see that the minimum re-expansion time is of the order of a few free-fall times, or $`1`$ Myr. This picture is consistent with observational estimates of core lifetimes, and, in particular, the proposed existence of failed cores is consistent with the known existence of a large number of starless cores (at least half of the total; e.g., Lee & Myers 1999). These cores have traditionally been considered to be in a pre-protostellar evolutionary stage, but in the turbulent scenario they may well never form stars (Taylor, Morata & Williams 1996; Vázquez-Semadeni et al. 2002).

We conclude this section by noting that it has all been based on the concept and numerical simulations of the non-magnetic case. Thus, within the context of the turbulent star formation scenario, the magnetic field is not an indispensible ingredient. The isolated and clustered modes of star formation, as well as the low SFE, can be well understood in terms of purely hydrodynamic, self-gravitating turbulence. We expect the presence of a magnetic field to introduce only quantitative modifications to this scenario.

## 3. Comparisons with observations

In order to assess the feasibility of the turbulent model of star formation, it is necessary to compare its results with the observations. A question that is frequently asked in relation to the turbulent model is whether it can explain the relatively quiescent nature of observed cores, where by “quiescent” it is normally meant subsonic, or transonic at most. However, as discussed in §2, the existence of subsonic velocity dispersions at small enough scales is a natural consequence of the turbulent cascade, and in fact, the “sonic scale” appears to play a key role in the control of the SFE, as described there.

Additionally, Barranco & Goodman (1998) have observationally found that the linewidth in a (small) sample of cores exhibits the interesting property of becoming independent of the offset from the core center at small offsets. This asymptotic value is slightly larger than the thermal value, and was interpreted by those authors as an indication that the non-thermal velocity dispersion becomes scale-independent in the cores, saturating at a slightly subsonic level. They referred to this phenomenon as “coherence” of the velocity field in the innermost regions of the cores, within scales $`<0.1`$ pc.

Finally, several recent observational works have fitted, with various degrees of success, Bonnor-Ebert (BE; Ebert 1955; Bonnor 1956) column density profiles to observed cores. BE spheres are hydrostatic, non-magnetic configurations, bounded by a hot, tenuous medium that contributes external pressure, but not to the gravitational potential of the system. The presence of such a hot confining medium is indispensible for the stability of the system, since otherwise the configurations need to smoothly extend to infinity, and in this case are generally unstable, as is, for example, the case of the singular isothermal sphere. Vázquez-Semadeni et al. (2002) have pointed out that such configurations cannot be realized in single-temperature media (as molecular clouds appear to be), since in this case the requisite of the hot, tenuous confining medium cannot be fulfilled. However, it is then puzzling that the observations seem to reasonably fit the column density profiles of BE spheres reasonably well.

In a couple of recent papers, we have addressed these two questions. Ballesteros-Paredes et al. (2003) considered a sample of 378 projections of the cores in the numerical simulations onto a hypothetical plane of the sky (POS). Column density profiles, averaged over the position angle on the POS, were produced, and the scatter around this average was recorded. Then, BE column density profiles were fitted to the average profiles, and accepted if the error between the fitted and average profile was smaller than the scatter observed in the profiles contributing to the average. Otherwise, the fits were rejected. This procedure yielded 65% of “acceptable” BE fits to the profiles, showing that the transient cores in the simulations have nevertheless angle-averaged profiles that often resemble BE ones. This can probably be attributed to the strong smoothing introduced by the angle-averaging.

The issue of a sub- or transonic nature of the cores has been addressed by Klessen et al. (2003), who investigated the velocity structure for the core sample studied by Ballesteros-Paredes et al. (2003). They found that, in simulations with gravity turned off, roughly half the cores have subsonic velocity dispersions. In simulations with gravity on, the number of subsonic cores decreases, but still the vast majority has velocity dispersions that do not exceed twice the sound speed, and roughly 10% have subsonic velocity dispersions. Moreover, the velocity dispersion was measured in this plot out to column densities of 1/5 times the maximum, which probably increases the measured velocity dispersion, in comparison to the standard procedure of measuring out to the FWHM only. Thus, the simulations are entirely consistent with the fact that molecular cloud cores are generally transonic and often even subsonic.

Klessen et al. (2003) also probed the cores in the simulations for “coherence”. Figure 2 shows the velocity dispersion vs. column density along all lines of sight through a few randomly selected cores, clearly showing a trend to acquire a constant value, often subsonic, towards high column densities. Thus, the cores in the simulations also exhibit “coherence” when projected on the POS. Since such independence of scale of the turbulent motions has no counterpart in normal turbulent flows, Klessen et al. (2003) suggested that the “coherence” may be an observational effect, probably due to the fact that near the projected core centers, the linewidth along a given line of sight is dominated by its length rather than by its typical scale on the POS (i.e., the contributing material lies within a cylinder that is longer than it is wide). If the column length does not vary much around the projected core maximum, neither will the linewidth, causing the apparent “coherence”.

## 4. Conclusions

In this paper we have given a global outline of the turbulent model of molecular clouds and star formation, and then summarized a number of tests of the simulations showing that the cores in the non-magnetic simulations studied exhibit similar density and velocity features to those observed in real cores, strongly supporting the validity of the turbulent model of molecular clouds. Provided that the turbulence in molecular clouds is continuously driven, the model explains in a natural way the origin and structure of the cores as transient turbulent density fluctuations in the clouds, and the low efficiency of star formation, since a) the cores contain a small fraction of the cloud’s mass, and b) not all cores are destined to collapse. Those that do not, rebound and merge back into their environment, and may generally correspond to the observed “starless” cores. This scenario does not necessarily rely on the magnetic field to operate, although the presence of such a field will most certainly introduce quantitative corrections into the model. Although further testing is still clearly necessary, and the role of the magnetic field needs to be quantified, a quantitative theory of star formation in turbulent molecular clouds appears to be under way.

### Acknowledgments.

The author gratefully acknowledges useful comments and criticisms from Susana Lizano and Javier Ballesteros-Paredes. This work has received financial support from CONACYT grant 36571-E.

## References

Ballesteros-Paredes, J. & Vázquez-Semadeni, E. 1997, in Star Formation, Near and Far, ed. S. Holt & L. Mundy (New York: AIP Press), 85

Ballesteros-Paredes, J. Vázquez-Semadeni, E. & Scalo, J. 1999, ApJ, 515, 286

Ballesteros-Paredes, J., Hartmann, L. & Vázquez-Semadeni, E. 1999, ApJ, 527, 285

Ballesteros-Paredes, J., Klessen, R. & Vázquez-Semadeni, E., 2003, ApJ, 592, 188

Barranco, J. A. & Goodman, A. A. 1998, ApJ 504, 207

Bonnor, W. B. 1956, MNRAS, 116, 351

Ebert, R. 1955, Zeitschrift für Astrophysik, 36, 222

Elmegreen, B. G. 1993, ApJL, 419, L29

Evans N. J. II 1991, in Frontiers of Stellar Evolution, ed. D. L. Lambert (San Francisco: Astronomical Society of the Pacific), 45

Hartmann, L., Ballesteros-Paredes, J. & Bergin, E. A. 2001, ApJ, 562, 852

Heitsch, F., Mac Low, M.-M. & Klessen, R. S. 2001, ApJ 547, 280

Hunter, J. H., Jr., & Fleck, R. C., Jr. 1982, ApJ, 256, 505

Klessen, R. S., Heitsch, F. & Mac Low, M.-M. 2000, ApJ 535, 887

Klessen, R. S., Ballesteros-Paredes, J., Vázquez-Semadeni, E. & Durán-Rojas, C. 2003, ApJ, submitted

Larson, R. B. 1985, MNRAS, 214, 379

Lee, C. W. & Myers, P. C. 1999, ApJS, 123, 233

Mac Low, M.-M., Klessen, R., Burkert, A., Smith, M. D. 1998, Phys., Rev. Lett., 80, 2754

Mac Low, M.-M. & Klessen, R. 2003, submitted to Rev. Mod. Phys., (astro-ph/0301093)

Mouschovias, T. Ch. 1976, ApJ, 207, 141

Padoan, P. 1995, MNRAS, 277, 377

Padoan, P. & Nordlund, Å. 1999, ApJ, 526, 279

Padoan, P. & Nordlund, Å. 2002, ApJ 576, 870

Sasao, T. 1973, PASJ 25, 1

Shu, F. H. 1977, ApJ 214, 488

Shu, F. H., Adams, F. C. & Lizano, S. 1987, ARAA, 25, 23

Stone, J. M., Ostriker, E. C., Gammie, C. F. 1998, ApJL, 508, L99

Taylor, S.D., Morata, O. & Williams, D.A. 1996, A&A 313, 269

Vázquez-Semadeni, E. Passot, T. & Pouquet, A. 1996, ApJ, 473, 881

Vázquez-Semadeni, E., Ostriker, E. C., Passot, T., Gammie, C. & Stone, J., 2000, in Protostars & Planets IV, ed. V. Mannings, A. Boss & S. Russell (Tucson: Univ. of Arizona Press), 3

Vázquez-Semadeni, E., Shadmehri, M. & Ballesteros-Paredes, J 2002, ApJ, submitted (astro-ph/0208245)

Vázquez-Semadeni, E. Ballesteros-Paredes, J. & Klessen, R. 2003, ApJL, 585, L131

von Weizsäcker, C. F. 1951, ApJ, 114, 165

Zuckerman, B. & Palmer, P. 1974, ARAA, 12, 279


format: text

text:
# Protein crystallization in confined geometries

## Abstract

We studied the crystallization of a globular protein, lysozyme, in the cubic phase of the lipid monoolein. The solubility of lysozyme in salt solution decreased by a factor of $`4`$ when confined in cubic phase. Calculations and Monte Carlo simulations show that this can be explained by the confinement of lysozyme molecules to the narrow water cells in the cubic phase.

We report a study of the crystallization of lysozyme, a globular protein, in the cubic phase of monoolein, an uncharged lipid. The protein solubility was found to decrease by a factor of $`4`$. We account for this observation by the confinement of lysozyme molecules within nanoscopic water ‘cells’ in the cubic phase. The resulting entropy-driven rise in their chemical potential can explain our observation of solubility suppression compared to the lipid-free case. The generic nature of this explanation means that our results are of potential interest to three rather different audiences. First, and most generally, scientists interested in the basic physics of confinement Gelb99  may find the lysozyme-monoolein system a useful model for detailed study. Secondly, there is growing interest in directing the self-assembly of synthetic or natural nanoparticles using various means of confinement (e.g. at an interface Dinsmore03  or in emulsion droplets Pine03 ). Our work suggests that lipid cubic phases can be useful in this regard. Thirdly, and most specifically, lipid cubic phases offer an environment for globular protein crystallization McPherson99 . (Their use for crystallizing membrane proteins Landau96 , which reside in the lipid bilayers, relies on different physics Caffrey00 ; grabe03 .)

Lysozyme crystallization in lipid cubic phase has been reported before Landau97 , although to date there is no complete phase diagram or physical explanation of the mechanism. We have studied the phase behavior of lysozyme solution as a function of protein and salt concentration in a monoolein cubic phase, and compared it with the case without lipids. The cubic phase was found to lower the solubility dramatically. This result was then explored using simple calculations and Monte Carlo simulations.

Lysozyme and monoolein (MO) were used as purchased from Sigma-Aldrich. Protein concentration in 0.05 M sodium acetate buffer was determined, after filtration through a 0.1 $`\mu `$m Millipore filter, by UV spectrophotometry (specific absorbance $`A_{280\mathrm{n}\mathrm{m}}=2.64`$ ml/mg.cm). NaCl solutions of known concentration were then added.

We first determined the bulk solubility boundary of lysozyme/salt without lipid from the dissolution of ‘seed’ crystals McPherson99 , which occurs at lower protein and salt concentrations than the bulk nucleation boundary (where crystals are first seen) pH4.5phase  due to the nucleation barrier McPherson99 . Next protein/salt solutions with concentrations below the bulk solubility boundary was added to and vigorously mixed with MO powder melted on a microscope slide. A highly viscous cubic phase was observed to form almost instantly. The excess solution was then removed and the cubic phase sealed under a cover slip (sample thickness $`100`$ $`\mu `$m) for observation in an optical microscope at 18C. All samples were transparent and optically isotropic under cross polarizers. Within two months, crystals were observed in many samples, e.g. Fig. 1.

Fig. 2 summarizes our observations in the form of a phase diagram. We give lysozyme concentrations, $`c`$, in mg/ml, and in terms of an effective volume fraction, $`\varphi `$, by modelling the globular protein as a sphere of diameter 3.4 nm (which reproduces the molecular volume). Note that for the data points recording behavior in cubo, the protein/salt concentrations given are those of the parent solution into which MO powder was added. The boundary above which we observed crystals in cubo (a kinetic nucleation boundary) lies entirely below the bulk solubility boundary. Since we expect the in cubo solubility boundary to lie at even lower protein/salt concentrations, our main experimental finding is that lipid cubic phase significantly depresses lysozyme solubility.

To understand this observation, we first note that the bicontinuous water in the cubic phase is made up of cells connected by narrow ‘necks’. The presence of lysozyme changes the MO cubic phase structure from Pn3m to Im3m Razumas96 , a schematic projection of which is shown in Fig. 3(a) with an inter-cell distance of 12.5 nm Razumas96  and bilayer thickness of 3.5 nm Caffrey00 . We have also shown a number of lysozyme molecules (hatched circles) modelled as spheres with diameter $`\sigma =3.4`$ nm. Since the bending rigidity of the bilayers is $`k_c9k_BT`$ at room temperature Chung94 , proteins cannot move between cells without incurring significant cost in distorting the connecting ‘necks’. Thus the lysozyme will largely be confined in nanoscopic water ‘cells’ with linear dimension $`10`$ nm, each of which can hold up to 7 lysozyme molecules Razumas96 ; structures . The center of a lysozyme molecule does not come closer than $`\sigma /2=1.7`$ nm to a bilayer. While such excluded volume is a small effect in the bulk, it is highly significant inside a water ‘cell’ in the cubic phase, leading to a large increase in the effective protein concentration. This explains why crystallization occurred at lower average protein concentration.

Enquiring where the crystals actually nucleate gives further understanding. The critical nucleus in our experiments is expected to consist of $`>10`$ lysozymes Galkin00 , more than an undistorted water ‘cell’ can hold. Moreover, a (tetragonal) unit cell contains 8 molecules Blake65 , and it is unlikely that crystals can nucleate in pores smaller than a few unit cells Mu96 . Experimentally, we often observed crystals next to solid inclusions or next to surfaces, Fig. 1. We therefore propose a role for defects.

The picture is as follows. The protein in bulk solution has a concentration-dependent chemical potential $`\mu (c)`$. Confinement in cubo increases the effective lysozyme concentration to $`c^{}>c`$, and its chemical potential to $`\mu (c^{})>\mu (c)`$. There are defects in the cubic phase of mesoscopic sizes ($``$ lysozyme diameters), due to grain boundaries, distortion by inclusions, or surfaces. In these defects, the chemical potential of lysozyme is initially at $`\mu (c)`$, but diffusion from the ‘loaded’ cubic phase rapidly increases this to close to $`\mu (c^{})`$. Crystals nucleate from this concentrated solution residing in the defects. Subsequently, as a crystal grows, it will break the surrounding matrix; the elasticity of the cubic phase may therefore affect the rate of crystal growth.

Support for the above picture comes from measuring the time-dependent protein concentration in lysozyme/salt solution in contact with the cubic phase formed by mixing it with MO. We observed (data not shown) a rapid increase in concentration that continued for many hours, showing that indeed $`\mu (c^{})>\mu (c)`$.

To explore quantitatively our proposed picture, we used Monte Carlo (MC) simulations to study attractive hard spheres confined within a ‘model cubic phase’. To do this, we first need a numerical estimate of the excluded volume effect. The center of a lysozyme molecule confined to a water ‘cell’ in the cubic phase cannot explore the whole of the cell volume, $`V_c`$, but only the accessible volume, $`V_a`$, given by taking off the volume of a layer of thickness $`\sigma /2`$ from $`V_c`$. The ratio $`\alpha =V_c/V_a`$ characterizes the excluded-volume effect.

Estimating $`\alpha `$ for the $`Im3m`$ structure involves non-trivial geometry and knowing the exact lattice parameter, $`a`$, pertaining to our experiments. The value of $`a`$ used in Fig. 3 is for fully-hydrated MO cubic phase equilibrated at 25C Razumas96 . These conditions almost certainly do not apply (after mixing at 25C we immediately removed excess solution and observed at 18C). Nevertheless, using $`a=12.5`$ nm Razumas96  as an estimate, we can approximate the accessible volume in each water ‘cell’ as a regular square bipyramid of side 8.8 nm, Fig. 3(b). The centres of lysozyme molecules can explore a smaller bipyramid of side $`8.8\sigma =5.4`$ nm, giving $`\alpha _{\mathrm{Im3m}}(8.8/5.4)^3=4.3`$ imagecharge . Literature data for MO under different conditions without Qui00  and with lysozyme Razumas96  suggest that we can expect variations of up to $`\pm 30\%`$ in the lattice parameter, so that $`\alpha 311`$.

In our simulations, we used small spheres to model lysozyme, while 125 larger spheres fixed on a primitive cubic lattice modelled the geometric confinement effect of the MO cubic phase. Such a model has been used before to study fluids in porous media Page96 . The small sphere diameter and large-sphere simple cubic lattice constant were taken to be 3.4 nm and 12.5 nm respectively, to agree with the known lysozyme molecular volume and the cubic phase with lysozyme lattice constant at 25C Razumas96 . We simulated a model with a particular $`\alpha `$ by using a value for the large sphere diameter such that the small-sphere centers can access $`\alpha ^1`$ of the total space between large spheres. Periodic boundary conditions apply.

The small spheres (diameter $`\sigma `$) interact via a hard core and a square-well attraction:

$$u(r)=\{\begin{array}{cc}+\mathrm{},\hfill & r<\sigma \hfill \\ ϵ\hfill & \sigma r<(1+\lambda )\sigma \hfill \\ 0\hfill & r(1+\lambda )\sigma \hfill \end{array},$$
(1)
where $`\lambda `$ and $`ϵ`$ are the range and depth of the square well in units of $`\sigma `$ and $`k_BT`$ respectively. We took $`\lambda =0.08`$ Kulkarni03 . Increasing the salt concentration experimentally corresponds to increasing $`ϵ`$; in our range of salinity, we expect $`0<ϵ<3`$ Kulkarni03 . The small spheres behave as hard spheres to the fixed large spheres.

We investigated the coexistence of a fluid of such attractive small hard spheres confined by a fixed array of large hard spheres with an unconfined crystal of the attractive hard spheres. This is the thermodynamic situation of crystals nucleating in mesoscopic defects filled with concentrated protein solution ‘fed’ from the cubic phase. Rather than simulating the free energies, we used an approximation to calculate coexistence by simply equating chemical potentials and neglecting the equal-pressure condition (see Sear99  for a careful justification).

First, in the spirit of the cell model Buehler51  we can approximate the volume accessible to the centre of a particle in the crystal to be $`\pi \lambda ^3\sigma ^3/6`$, so that the chemical potential of the crystal can be written as,

$$\mu _s=\mu _0n_s\left(\frac{ϵ}{2}\right)3k_BT\mathrm{ln}\lambda ,$$
(2)
where $`n_s`$ is the number of particles whose centers lie in the region $`\sigma r<(1+\lambda )\sigma `$ from a given particle. Here, $`\mu _0=k_BT\mathrm{ln}(6\mathrm{\Lambda }^3/\pi \sigma ^3)`$ where $`\mathrm{\Lambda }`$ is the de Broglie (thermal) wavelength. Note that with this choice of $`\mu _0`$, the ideal-solution chemical potential is given by $`\mu _i=\mu _0+k_BT\mathrm{ln}\varphi `$.

Next, consider the fluid of confined small spheres. The leading order correction to the ideal limit is

$$\mu _f=k_BT\mathrm{ln}(\rho _a\mathrm{\Lambda }^3)=\mu _i+k_BT\mathrm{ln}\alpha ,$$
(3)
where $`\rho _a`$ is the number density of small spheres in the volume accessible to them, and $`\alpha ^1`$ is the ratio of this accessible volume to the total volume between large spheres. In our simulations, $`\alpha =2.0,4.0,6.0`$ or 8.0. We used the ‘Widom insertion method’ widom  to obtain $`\mathrm{\Delta }\mu =\mu _f\mu _i`$ at different values of $`ϵ`$ with and without confinement, Fig. 4. Note that the confined results indeed go to $`\mathrm{\Delta }\mu /k_BT=\mathrm{ln}\alpha `$ in all cases when $`\varphi 0`$.

Now we calculate fluid-crystal coexistence at each value of $`ϵ`$ by equating the crystal chemical potential, Equation 2 with $`n_s=12`$, to the simulated fluid chemical potential Sear99 . The low-density branches with and without confinement are shown in Fig. 5. Increasingly pronounced confinement indeed systematically decreases the solubility in a way very reminiscent of the data shown in Fig. 2.

Quantitatively, we show in the inset that multiplying the volume fraction $`\varphi `$ by $`\alpha `$ leads to overlap of the phase boundaries for $`\alpha \varphi 0.2`$. Note that the same result is obtained (data not shown) if we map $`ϵ`$ to the salt concentration via measured second virial coefficients Kulkarni00 ; Kulkarni03 .

This result suggests that we should attempt the same procedure on the experimental data. This is problematic. The solubility boundary in cubo cannot be measured: Fig. 2 shows a nucleation boundary. If we nevertheless proceed, we find that multiplying the lysozyme concentration by a factor of 4 shifts the confined nucleation boundary to the bulk solubility boundary (inset, Fig. 2). The in cubo solubility boundary will undoubtedly occur at lower concentrations than the nucleation boundary. We therefore expect that a ‘shift factor’ of $`>4`$ will be needed to translate the (unavailable) in cubo solubility boundary to overlap with the bulk solubility boundary caveat . Given that we expect $`\alpha 311`$ in our cubic phase, this expectation is not inconsistent with our simulation results.

In reporting our experimental findings, we have so far not mentioned any time-dependent features. Two aspects are worthy of note. First, we often observed amorphous precipitates before the appearance of crystals. Unlike crystals, these precipitates would occur rather uniformly throughout a sample, and often appeared as ‘stars’ in the microscope. Sometimes, growing crystals would ‘consume’ surrounding amorphous precipitates. The origin and nature of these precipitates are unclear. Secondly, although it is difficult to be quantitative, it is our impression that the in cubo growth of crystals is a slower process compared to that in bulk solution. This may, of course, be a reflection of the fact that ‘feeding’ lysozyme into the growth medium residing in defects is a slow process.

To conclude, we have studied the phase behavior of lysozyme/salt in MO cubic phase. The latter lowers the solubility of the protein. This is predominantly a confinement effect: experiments and simulations suggest that a single ‘excluded volume scaling factor’ can shift the solubility boundary back to its bulk position. The altered kinetics of crystallization in cubo, relying as we suggest on defects and other properties of the cubic phase, may offer a new means of controlling protein crystal nucleation and growth. The generic, entropic, nature of our proposed explanation suggests that lipid cubic phases may be useful for manipulating a range of nanoparticles, and for studying the basic statistical physics of confinement.

ST held a JSPS Overseas Research Fellowship. We thank Mike Cates for discussions relating to Sear99  and critical reading of the manuscript.


text:
# Untitled Document

Another Proof of a Famous Inequality

1. The mean-value theorem of differential calculus, in its simplest form, states that if a function $`\varphi `$ has a derivative at all points between $`c`$ and $`d`$, and including those points, then for some $`e`$ strictly between $`c`$ and $`d`$

$$\varphi (d)\varphi (c)=(dc)\varphi ^{}(e).$$

The word “between” was used as the order of $`c`$ and $`d`$ is unimportant. However in the following we will assume that $`c<d`$. 2. The basic inequality of the title is: if $`ab`$ are positive and if $`u,v>0`$ then

$$\frac{ua+vb}{u+v}>\left(a^ub^v\right)^{1/(u+v)}.$$
$`(GA)`$
This is the simplest case of the inequality between the arithmetic and geometric means, from which the general case follows by a simple induction; see \[2, p.81\].

Two things should be remarked for our later discussion: (i) as both sides of (GA) are symmetric in $`a`$ and $`b`$ there is no loss in generality in assuming $`a<b`$;

(ii) as (G) is unchanged if $`u,v`$ are replaced by $`\lambda u,\lambda v`$ respectively, for any $`\lambda >0`$, we can assume, without loss in generality, that $`u,v1`$ 3. Consider then the distinct power functions $`\varphi _1(x)=x^u,\varphi _2(x)=x^v`$; where $`uv`$, $`u,v1`$ and $`x>0`$ . Applying the mean value theorem to both of these functions we get, on cancelling the common factor $`dc`$,

$$\frac{d^uc^u}{d^vc^v}=\frac{ue_1^{u1}}{ve_2^{v1}},$$
$`(1)`$
for some $`e_1,e_2`$, with $`c<e_1,e_2<d`$.

Using the various conditions (1) leads to

$$\frac{d^uc^u}{d^vc^v}>\frac{uc^{u1}}{vd^{v1}}>\frac{uc^u}{vd^v};$$
$`(2)`$
or on multiplying out

$$vd^v(d^uc^u)>uc^u(d^vc^v)$$
$`(3)`$
To get (3) we have assumed that $`uv`$ but it is easy to check that (3) remains valid when $`u=v`$.

Rewriting (3) gives

$$uc^{u+v}+vd^{u+v}>(u+v)c^ud^v.$$
$`(4)`$
Inequality (4) is the just (GA) as can be seen by putting $`a=c^{u+v},b=d^{u+v}`$. 3. It is possible to have $`e_1=e_2`$ in (1). However then we must use the more sophisticated Cauchy mean-value theorem: if $`\varphi _1,\varphi _2`$ satisfy the condition that $`\varphi `$ satisfies section 1 above, and if in addition $`\varphi _2^{}`$ is never zero then for some $`e`$ strictly between $`c`$ and $`d`$,

$$\frac{\varphi _1(d)\varphi _1(c)}{\varphi _2(d)\varphi _2(c)}=\frac{\varphi _1^{}(e)}{\varphi _2^{}(e)}.$$

Now choosing the $`\varphi _1,\varphi _2`$ as in section 2 above we get (1) with $`e_1=e_2`$, that is

$$\frac{d^uc^u}{d^vc^v}=\frac{ue^{u1}}{ve^{v1}}=\frac{ue^u}{ve^v}.$$

The only advantage of this approach is that we now need not assume $`u,v>1`$, $`u,v>0`$ will suffice. 4. This proof of (GA) is based on an idea in ; see also \[2,p. 78\]. There $`u,v`$ are integers bigger than $`1`$ and to get (2) the identity

$$\frac{d^nc^n}{dc}=c^{n1}+c^{n2}d+\mathrm{}+cd^{n2}+d^{n1},$$

obtained by summing a geometric series, is used. For then in the case $`c<d`$ and $`u,v2`$, integers,

$$\frac{d^uc^u}{d^vc^v}=\frac{c^{u1}+\mathrm{}+d^{u1}}{c^{v1}+\mathrm{}+d^{v1}}>\frac{uc^{u1}}{vd^{v1}},$$

which is (2).

Bibliography

1. T.Arnold Brown. Elementary inequalities,  Math. Gazette, 263 (1941), 2–11.

2. P.S.Bullen. Handbook of Means and Their Inequalities, Kluwer Academic Publishers, Dordrecht/Boston/London. 2003.


text:
# Monte Carlo simulation of an experiment looking for radiative solar neutrino decays

## 1 Introduction

In the last few years it has become clear that neutrinos have non-vanishing masses, and that the neutrino flavor eigenstates ($`\nu _e`$, $`\nu _\mu `$ and $`\nu _\tau `$) are superpositions of mass eigenstates ($`\nu _1`$, $`\nu _2`$ and $`\nu _3`$). For a recent review, see . In this context, neutrinos could undergo radiative decays, e.g. $`\nu _2\nu _1+\gamma `$, as initially suggested in . The present status of decaying theory is sumarized in . Such decays request that the involved neutrinos have a non-vanishing electric dipole moment; the very stringent existing experimental limits refer to the flavor neutrino eigenstates and they are not directly applicable to possible dipole moments of mass neutrino eigenstates.

In a pioneering experiment performed during the Total Solar Eclipse (TSE) of October 24, 1995 a search was made for visible photons emitted through possible radiative decays of solar neutrinos during their flight between the Moon and the Earth . In the analysis of the data, the authors assumed that all neutrinos have masses of the order of few eV, $`\mathrm{\Delta }m_{12}^2=m_2^2m_1^210^5`$ eV2, and an average neutrino energy of 860 keV; furthermore they assumed that all decays would lead to visible photons, which would travel nearly in the same direction as the parent neutrinos, thus leading to a narrow spot of light coming from the direction of the center of the dark disk of the Moon.

Subsequently, Frère and Monderen made more accurate calculations on the shape of the expected signal , considering the Sun as an extended source of electron neutrinos with typical energies of the order of 1 MeV, a Lagrangian formalism for the radiative decay, different neutrino masses (1 eV or 0.5 eV), and mass squared differences $`\mathrm{\Delta }m^2`$ of $`10^5`$, 0.25 and 1 eV2; they have shown that the expected signal could have an extended angular pattern.

Some of the authors of this paper were involved in two experiments along the line of , during the total solar eclipses of August 11, 1999 (in Romania) , and of June 21, 2001 (in Zambia) . In 1999 the bad weather conditions did not allow the planned observations, but we could use a videotape filmed by a local television (Râmnicu Vâlcea). The analysis of the data was performed in the hypothesis of a possible

$$\nu _2\nu _1+\gamma $$
(1)
decay, with $`m_2>m_1`$. For the analysis of the 1999 data we have chosen the two $`\mathrm{\Delta }m^2`$ values suggested by the MSW SMA (Small Mixing Angle) and LMA (Large Mixing Angle) solutions of the Solar Neutrino Problem (SNP), allowed by the then available experimental data from solar neutrino experiments. We developed a Monte Carlo (MC) simulation of the radiative solar neutrino decay, considering the solar neutrino energy spectrum predicted by the Standard Solar Model (SSM)  and the mass of the $`\nu _1`$ in the range of 1 - 10 eV, as it was expected at that time. Since the angular resolution of the data was not very good, we considered the Sun as a pointlike neutrino surce. The simulation has shown that the expected signal should be a narrow spot of light in the direction of the center of the Sun, and allowed an evaluation of the fraction of decays yielding visible photons as function of the chosen neutrino $`\nu _1`$ mass $`m_1`$ and $`\mathrm{\Delta }m^2`$ values.

The 2001 experiment lead to better quality data, so the real spatial distribution of the solar neutrino yield had to be considered. Furthermore, the recent SNO results  favour the LMA solution and could indicate also the presence of a $`\nu _\tau `$ component in the solar neutrino flux at the Earth level.

The WMAP (Wilkinson Microwave Anisotropy Probe) results after the first year of flight  limit the sum of the masses of the three neutrino species to 0.23 eV ( 95% Confidence Level).

In this paper we present a more complete simulation. We assume that $`m_1<m_2<m_3`$ where $`m_1`$, $`m_2`$, $`m_3`$ are the masses of the $`\nu _1`$, $`\nu _2`$ and $`\nu _3`$ mass eigenstates, respectively; but we restrict our analysis to a two generation mixing scenario, assuming the present mass differences obtained from solar neutrino experiments, the LMA solution with $`\mathrm{\Delta }m_{12}^2=6\times 10^5`$ eV2. Since SNO suggests also the presence of $`\nu _3`$ in the solar neutrino flux, we considered also the mass difference measured by atmospheric neutrino experiments : $`\mathrm{\Delta }m_{13}^2\mathrm{\Delta }m_{23}^2=2.5\times 10^3`$ eV2. We included in our MC all the details of neutrino production in the Sun, as given by the “BP2000” SSM .

The aim of our simulation is to give information on the shape of a possible decay signal of solar neutrinos (angular and energy distributions of the emitted photons) and to estimate the probabilities for the decay and the geometrical detection efficiency. This information should allow to obtain limits on the neutrino lifetimes from the experimental observations. As it is shown in the following sections, some of the input parameters of the code refer to the characteristics of a specific experiment. The results presented in this paper are obtained in the conditions of the observations made in 2001 in Zambia  with a digital videocamera.

## 2 The Monte Carlo simulation

The assumed geometry for the simulation of the radiative decay of solar neutrinos during a TSE is shown in Fig. 1. The notations in the figure will be used in the following subsections.

### 2.1 Neutrino production inside the Sun

In simulating the solar neutrino production, we used the “BP2000” SSM  in its numerical form directly available from . The first step consists in chosing a specific reaction/decay yielding neutrinos (both from the p-p and the CNO cycles) according to its predicted contribution to the solar neutrino flux at the Earth. The neutrino energy and the radius $`R`$ of its production point are then randomly generated according to the SSM. In Fig. 2 we present the neutrino solar energy spectrum obtained from $`3\times 10^7`$ generated neutrinos, and the distribution of the distance $`R`$ from the center of the Sun of the production points.

In order to ensure a uniform distribution of the birth points of the neutrinos ($`P_1`$ in Fig. 1) on each shell $`\mathrm{\Delta }R`$ at radial distance $`R`$, we determine for each neutrino the angular spherical coordinates generating uniformely the cosinus of the zenithal angle ($`\theta _S`$ in Fig. 1, assuming the “z” axis oriented from the centre of the Sun towards the centre of the Earth) and the polar angle ($`\varphi _S`$ in Fig. 1.).

### 2.2 Propagation of solar neutrinos and of the decay photons

The solar neutrinos are produced as $`\nu _e`$’s, thus as a superposition of neutrino mass states. After leaving the Sun, during their flight in the interplanetary space, their weak flavor is irrelevant.

In our simulation we consider the $`\nu _2`$ (or $`\nu _3`$) component of the neutrino flux, that may decay into the lower mass state $`\nu _1`$ and a photon. As we are interested in the possibility to observe such photons during a total solar eclipse, we impose that the decay processes take place in the space between the Moon and the Earth, inside the shadow cone of the Moon (otherwise the separation of the decay photons from the solar light backround would be impossible). Furthermore, the decay photons must reach the detector (a CCD-equipped telescope, a digital camera or some other observation systems of the same kind).

The next step in the simulation is to define the incidence direction of the decay photon on the detector, situated on the Earth surface at the location $`P_3`$ (see Fig. 1). We generate uniformly the cosinus value of the “local” zenith angle $`\theta _E`$ inside the shadow cone (or in the cone subtended by the telescope or by the analysis apperture if smaller) and then the value of the “local” azimuthal angle $`\varphi _E`$, as defined in Fig. 1.

The probability of neutrino radiative decays (thus of their lifetime) depends on their electric dipole momentum, which, for the weak flavor eigenstates, is limited by the existing experimental data . Those limits are not directly applicable for the mass eigenstates. We may assume that the lifetime of $`\nu _2`$ is much larger than the time of flight from the Sun to the Earth (otherwise the experiment would be impossible). This implies that the decay points of the massive solar neutrinos are uniformely distributed along their path from the Moon to the Earth.

As we already defined the direction of the “detected” photon, we cannot choose the decay point of the neutrino using the above assumption; instead we take advantage of the nearly cylindrical symmetry of our problem to choose the decay point (point $`P_3`$ in Fig. 1) randomly along the photon path and then determine the neutrino path as the segment starting in $`P_1`$ and ending in $`P_3`$. This does not affect the expected uniformity of the distribution of the decay points along the neutrino path from the Moon to the Earth, as shown in Fig. 3a. We also checked the isotropy on the incidence of solar neutrinos on the Moon surface (see Fig. 3b).

### 2.3 Neutrino decays

For every simulated event we know (from the steps described in the above subsections) the 4-vector of the “heavy” neutrino $`\nu _2`$, its birth and decay points, and the direction of the photon emission; from the 4-momentum conservation we obtain the energy (in the laboratory reference frame, LRF) of the emitted $`\gamma `$

$$E_\gamma =\frac{\mathrm{\Delta }m^2}{2}\frac{1}{E_\nu p_\nu \mathrm{cos}\theta },$$
(2)
where $`E_\nu `$ and $`p_\nu `$ are the energy and the momentum of the parent $`\nu _2`$ neutrino, respectively, and $`\theta `$ is the emission angle of the photon relative to the direction of the $`\nu _2`$ momentum in the LRF (see Fig. 1).

In the center of mass (CM) of the decaying $`\nu _2`$ neutrino, the probability density of the emission of a photon at the zenithal angle $`\theta ^{}`$ is given by:

$$\frac{d\mathrm{\Gamma }}{d(\mathrm{cos}\theta ^{})}=K(1+\alpha \mathrm{cos}\theta ^{})$$
(3)
where the $`\alpha `$ parameter depends on the polarization state of the initial $`\nu _2`$ flux: $`\alpha =0`$ for unpolarized (Majorana) neutrinos, and $`\alpha =1`$ for left and right handed (Dirac) $`\nu _2`$’s, respectively. The constant $`K`$ comes from the general description of the decay of a fermion into a fermion and a boson

$$K=\frac{\alpha _e^2}{\pi ^2}\frac{m_2}{(\mathrm{\Delta }m^2)^3}(m_1^2+m_2^2+m_1m_2),$$
(4)
where $`\alpha _e`$ is the fine structure constant.

In order to estimate the probability of each simulated event, one would have to integrate Eq. 3 inside the solid angle under which the emitted photon “sees” the detector. This is unpractical since the physical dimension of any usable detector is too small compared to the distance scales involved in the simulation. An equivalent approach is based on the observation that any imaging system (such as a CCD in the focal plane of some optical system) has a limited angular resolution, as the images of all point sources inside the solid angle covered by a single pixel will be superimposed. We can then numerically integrate Eq. 3 on the area around the decay point $`P_2`$ (see Fig. 4 for a sketch; the integration area is represented as the dark square centered in $`P_2`$) that corresponds to the angular aperture of a single pixel. The probabilities obtained in this way are used as weights for the MC generated events.

## 3 Results and discussions

The results presented in this Section are obtained assuming an angular resolution of each pixel of 10”, and an angular aperture of the analysis (the maximum value of $`\theta _E`$ in Fig. 1) of 480”, as in the case of the digital video-camera used during the 2001 TSE . The mass $`m_1`$ of the $`\nu _1`$ mass eigenstate was considered in a range between $`10^3`$ eV and 0.3 eV; MC simulations were made for $`\mathrm{\Delta }m^2=6\times 10^5`$ eV2 and also for $`\mathrm{\Delta }m^2=2.5\times 10^3`$ eV2. For the polarization parameter $`\alpha `$ in Eq. 3, we considered three possible values: -1, 0 and 1.

For each combination of neutrino mass and $`\mathrm{\Delta }m^2`$ we requested $`3\times 10^4`$ unweighted events yielding photons in the visible range. Assuming $`\mathrm{\Delta }m^2=6\times 10^5`$ eV2 the total number of iterations was about $`1.6\times 10^9`$, while for $`\mathrm{\Delta }m^2=2.5\times 10^3`$ eV2, about $`3.3\times 10^7`$ generated events were needed. The simulated events were then weighted according to the integral of Eq. 3, as discussed in the previous Section. In weighting the events the factor $`K`$ in Eq. 3 was set to 1, as it is a constant for each $`(m_2,\mathrm{\Delta }m^2)`$ hypothesis. For this analysis we also imposed the initialization of the random number generator to be the same for each run.

Fig. 5 shows the fraction of visible photons, versus the mass of the lighter neutrino $`\nu _1`$. Fig. 5a corresponds to $`\mathrm{\Delta }m^2=6\times 10^5`$ eV2, while the results for $`\mathrm{\Delta }m^2=2.5\times 10^3`$ eV2 are presented in Fig. 5b. The light triangles are obtained assuming $`\alpha =1`$, while the dark and light circles correspond to $`\alpha =0`$ and $`\alpha =+1`$, respectively.

As discussed in the previous Section, for each simulated neutrino decay event we numericaly integrated Eq. 3, thus obtaining a probability that includes the contributions from the kinematics of the decay itself as well as from the a priori request included in the simulation that the emitted photon reaches the detector. In Fig. 6 we present these probabilities, averaged over all Monte Carlo events yielding visible photons, versus the mass of the $`\nu _1`$ neutrino. Fig. 6a refers to $`\mathrm{\Delta }m^2=6\times 10^5`$ eV2, and Fig. 6b to $`\mathrm{\Delta }m^2=2.5\times 10^3`$ eV2. The symbols used for different polarization states are the same as in Fig. 5.

The probabilities shown in Fig. 6 are very small, but one should consider that they apply to all the solar massive neutrinos that cross the “acceptance cone” of the detector between the Earth and the Moon (see Fig. 1).

Both Figs. 5 and 6 show a strongly non-linear behaviour, as the condition imposed to the photons to be in the visible spectrum selects different regions of the solar neutrino spectrum for different mass or polarization hypothesis. This effect is illustrated in Fig. 7, for three values of the neutrino mass $`m_1`$, 0.001, 0.01 and 0.1 eV: the solid, dashed and dotted histograms respectively. Fig. 7a corresponds to $`\mathrm{\Delta }m^2=6\times 10^5`$ eV2, while Fig. 7b to $`\mathrm{\Delta }m^2=2.5\times 10^3`$ eV2. In all cases the polarization parameter $`\alpha `$ was assumed -1.

Fig. 7 suggests that the high energy solar $`{}_{}{}^{8}B`$ neutrino tail does not yield visible photons through radiative decays. Such process could instead happen for low energy $`pp`$ neutrinos, with some contributions from the $`{}_{}{}^{13}N`$, $`{}_{}{}^{15}O`$ neutrinos and from the $`{}_{}{}^{7}Be`$ and $`pep`$ lines; notice that these last contributions are more noticeable in Fig. 7a than in Fig. 7b.

In conducting an experiment searching for visible photons from a hypothetical radiative solar neutrino decay during TSE’s and in choosing the appropriate data analysis methodology, the simulation of the expected signal is very important. Fig. 8 presents such simulations, in the same conditions and with the same conventions as in Fig. 7. The shape of the signals corresponding to different neutrino masses and mass differences (assuming left handed neutrinos, thus $`\alpha =1`$) is presented as the radial distribution of the “detected” luminosity. This is equivalent to the distribution of the weighted number of visible photons, averaged on circular corona around the same value of the incidence zenith angle $`\theta _E`$ (see Fig. 1).

The histograms in Fig. 8a correspond to the simulated $`\nu _2\nu _1+\gamma `$ decays. For all neutrino masses, the expected signal is concentrated at small $`\theta _E`$ angles (about 50 arcsec). The widths and shapes of the signals are sensitive to the mass assumed: the larger the mass, the narrower the signal band. The two peaks seen at about 300 and 400 arcsec. for $`m_1=10^3eV`$ could be correlated with the contribution of $`{}_{}{}^{7}Be`$ or $`pep`$ neutrinos; they are about a factor 100 lower than the central maxima, so their contribution is not measurable.

In the case of $`\nu _3\nu _1+\gamma `$ simulated decays (Fig. 8b), the signal is broader (about 250 arcsec) and is less sensitive to the mass choice. The peaks observed at about 250” could have a similar origin as those in Fig. 8a, but could also be statistical fluctuations.

If an experiment on solar neutrinos has also a good energy resolution, then further information could be obtained from the analysis of the spectra of the observed signal. Such spectra are shown in Fig. 9, assuming only left-handed neutrinos, and considering the same examples as in Figs. 7 and 8.

As for the shape of the signal, its spectral decomposition seems to be more sensitive to the neutrino mass values for the $`\nu _2\nu _1+\gamma `$ decays (Fig. 9a). In this case most of the visible photons are “detected” in the red part of the spectrum, while the spectra in Fig. 9b suggest a dominant signal in the green, due to the larger mass difference between $`\nu _3`$ and $`\nu _1`$.

Let us assume that an experiment as those simulated in this work would measure, during a TSE, an excess of visible photons $`\mathrm{\Phi }_{obs.}`$ from the direction of the center of the sun, or, if the search yields a negative result, $`\mathrm{\Phi }_{obs.}`$ is the photon flux of the observed fluctuations. In the first case, estimates of the neutrino lifetime could be done; otherwise, a lower experimental limit could be deduced. Assuming that solar (electron) neutrinos are superpositions of only two mass eigenstates,

$$|\nu _e>=|\nu _1>\mathrm{cos}\theta +|\nu _2>\mathrm{sin}\theta ,$$
(5)
where $`m_2>m_1`$ and $`\theta `$ is the mixing angle, the average lifetime (or its lower limit) $`\tau `$ of the $`\nu _2`$ neutrino could be computed from

$$N_\gamma =P\mathrm{\Phi }_2S_Mt_{obs}\left(1e^{\frac{<t_{ME}>}{\tau }}\right)e^{\frac{t_{SM}}{\tau }},$$
(6)
where $`N_\gamma `$ is the number of decay visible photons observed, $`P`$ are the mass - dependent probabilities shown in Fig. 6 a, $`\mathrm{\Phi }_2=\mathrm{\Phi }_\nu \mathrm{sin}^2\theta `$,(where $`\mathrm{\Phi }_\nu `$ is the flux of solar neutrinos at the Earth (or Moon) level), $`S_M`$ is the area of the Moon surface covered by the analysis (the base of the cone of angle $`\theta _E`$ in Fig. 1) and $`t_{obs}`$ is the time of observation. $`<t_{ME}>`$ is the average time spent by solar neutrinos inside the observation cone (about one third of the flight time from the Moon to the Earth), and $`t_{SM}`$ is the time of flight of the neutrinos from the Sun to the Moon. The low numerical values of the probabilities $`P`$ are compensated by the large solar neutrino flux, combined with the large area $`S_M`$, so an experiment as the one simulated here could yield at least significant upper limits on the $`\nu _2`$ lifetime $`\tau `$. In the conditions of a 3.5 minutes long TSE (as that of 2001) observed with an instrument with characteristics similar to those considered in this simulation, one would expect a $`\nu _2`$ lifetime (in the proper reference frame) sensitivity ranging from few seconds to about $`10^4`$ seconds, assuming $`\nu _2`$ neutrino masses of few $`10^2`$ eV.

## 4 Conclusions

In this paper we presented a Monte Carlo simulation for an experiment looking for visible photons emitted by a possible solar neutrino radiative decay, during a total solar eclipse. It was shown that for neutrino masses smaller than 0.1 eV and assuming squared mass differences in agreement with the Large Mixing Angle Solution (LMA) of the solar neutrino oscillations , such an experiment could give also an estimate of the neutrino mass.

The analysis of the experimental data collected by some of the authors during the 2001 TSE in Zambia is being completed using the simulation results obtained in this paper.

## 5 Acknowledgements

We would like to aknowledge many colleagues for useful comments and discussions.

This work was founded by NATO Grant PST.CLG.977691 and partially supported by the Italian Space Agency (ASI) and the Romanian Space Agency (ROSA).


text:
# The PLANET microlensing campaign: Implications for planets around galactic disk and bulge stars

## References

Dominik, M., et al. 2002, P&SS, 50, 299

Gaudi, B. S., et al. 2002, ApJ, 566, 463

Horne, K. 2003, MNRAS, submitted


text:
# A Fallback Disk Model for Ultraluminous X-ray Sources

## 1 Introduction

Gamma-ray bursts (GRBs) are the most brilliant of all astronomical explosions in our universe. Although we lack direct evidence on the nature of the central engine, it is widely accepted that GRBs may involve black holes accreting rapidly through a disk or torus. In this process some fraction of the accreted mass energy and the rotational energy of the black holes is extracted and goes into powering jets that induce GRBs.

There has been accumulating evidence supporting that long-duration ($`>2`$ s) GRBs are associated with the deaths of massive stars (Galama et al., 1998; Fruchter et al., 1999; Bloom et al., 1999; Garnivich et al., 2001). The temporal coincidence of SN2003dh with GRB030329 (Hjorth et al., 2003) strongly suggests that core collapse events can give rise to GRBs, thereby favoring the ”collapsar” model (Woosley, 1993; MacFadyen & Woosley, 1999).

The study of collapsars or ”failed supernovae” was initiaed by Bodenheimer & Woosley (1983) and the model has been explored as a GRB progenitor (Woosley, 1993). During the core collapse of stars with initial mass higher than $`2025M_{}`$, the black hole may form either promptly, since a successful outgoing shock fails to be launched by the collapsed iron core, or in a mild explosion by fallback. In both cases a significant fraction of the core with mass of $`0.15M_{}`$ may fall back onto the collapsed remnants over a period of minutes to hours (Woosley, 1993). Such a fallback would be induced by a deep gravitational potential of relatively massive stars or by a reverse shock propagating inward from the outer composition interface (Chevalier, 1989). If most of the infalling gas possesses sufficient angular momentum, the formation of a gaseous disk is highly feasible.

Although the formation of an accretion disk around a newborn black hole is a natural outcome in most popular models of GRBs, there has not been direct observational evidence for the existence of such a disk. The influence of the disk accretion on the supernova light curve (Balberg, Zampieri, & Shapiro, 2000), and possible gravitational radiation (Mineshige et al., 2002) have been discussed in the literature.

Theoretically, evolution of the fallback disk around an isolated compact star has been considered by e.g., Meyer & Meyer-Hofmeister (1989), Cannizzo, Lee, & Goodman (1990), Mineshige, Nomoto, & Shigeyama (1993), and Mineshige et al. (1997). The accretion rates were found to generally have the power-law dependence on time as $`\dot{M}t^a`$ with $`1<a<1.5`$. This is expected by the self-similar solutions for a viscous disk under the assumption of constant angular momentum. If disk accretion occurs, it is natural to expect that the collapsar remnants would eventually become luminous X-rays sources after the initial, rapidly accreting phase, for sufficient long time with the disk luminosity around the Eddington luminosity (see below). In this Letter we argue that some of the ultraluminous X-ray sources (ULXs) discovered in nearby galaxies are likely to evolve from collapsars.

## 2 ULXs as black holes accreting from a fallback disk

ULXs are the most luminous point-like extra-nuclear X-ray sources found in nearby galaxies, with observed (isotropic) X-ray luminosities in excess of $`10^{39}`$ ergs-1, roughly corresponding to the Eddington limit accretion luminosity for a 10 $`M_{}`$ star (Fabbiano, 1989). The majority of ULXs are believed to be accreting black holes, according to their spectral and timing characteristics (Makishima et al., 2000; Fabbiano et al., 2001). Currently there are two categories of models proposed for ULXs. (1) They belong to a new type of $`10^210^5M_{}`$ intermediate-mass black holes accreting at sub-Eddington rates (Colbert & Mushotzky, 1999; Miller & Hamilton, 2002). (2) They are traditional stellar-mass black holes in binary systems, but with either truly super-Eddington X-ray emission (Begelman, 2002), or anisotropic (King, 2001)/relativistically beamed emission (Körding, Falcke, & Markoff, 2002). Though some individual ULXs could be intermediate-mass black holes, the latter model seems to be more favored by the association of ULXs with star formation regions (Zezas et al., 2002).

Here our interest is focused on the ULXs that were found to be associated with supernova remnants (SNRs). IC 342 X-1 is one of the nearest and most comprehensively studied ULXs with an X-ray luminosity $`>10^{40}`$ ergs-1 in $`0.510`$ keV (Makishima et al., 2000). Recent optical observations (Pakull & Mirioni, 2002; Roberts et al., 2003) show a large diameter ($`110`$ pc) SNR encircling the position of the ULX, suggesting that the ULX may be physically related to the SNR. Further evidence for the association comes from the X-ray ionized nebula inside the SNR (Roberts et al., 2003).

Since the age of the SNR is $`10^5`$ years (Roberts et al., 2003), the presence of the ULX implies that it was possibly formed very recently. If it is in a binary system, mass transfer would occur in terms of stellar wind capture and/or Roche lobe overflow from the donor star. However, to produce a luminous X-ray source only $`10^5`$ years after the black hole was formed requires that both the primary (the black hole progenitor) and the secondary (donor) stars should have extremely fine-tuned initial masses (see Bhattacharya & van den Heuvel (1991) for a discussion on the formation and evolution of X-ray binaries). A much more plausible explanation is that the black hole is accreting from the disk originating from the fallback supernova debris, which implies that this source had probably experienced a GRB. The inferred initial energy input to the SNR is indeed at least 2$``$3 times greater than the canonical energy for an ordinary SNR (Roberts et al., 2003). The presence of the ULX within an unusually energetic SNR strengthens the suggestion that it may be the aftermath of a GRB.

Another example is the ULX in the SNR MF 16 in NGC 6946 (Roberts & Colbert, 2003), which is possibly $`3500`$ years old (Blair & Fesen, 1994). If the association is correct, the accreting material is almost definitely from the fallback disk.

We note that nebulae have been reported to be present at or around the positions of several ULXs by Pakull & Mirioni (2002). Some of them could be SNRs related to the birth of ULXs. These nebulae have kinematic ages of $`\stackrel{<}{}`$million years and appear to be directly linked to the highly energetic formation process of the ULXs or being inflated by ongoing stellar wind/jet activity. Fallback disk accretion seems to be a viable mechanism for the X-ray emission from these sources.

Radio emission was recently discovered from a ULX in the dwarf irregular galaxy NGC 5408 (Kaaret et al., 2003). The X-ray, radio and optical fluxes as well as the X-ray spectral shape are consistent with beamed relativistic jet emission from an accreting stellar-mass black hole, which is a strong indicator of disk accretion. The ULX is found to be close to the star formation regions of NGC 5408, with a projected displacement of 280 pc, and the ULX could traverse this distance in $`\stackrel{<}{}3`$ Myr with a speed of 100 kms-1 (Kaaret et al., 2003). This suggests that the source may belong to young stellar population, i.e., high-mass X-ray binaries (HMXBs). However, the high X-ray to optical flux ratio $`f_\mathrm{x}/f_\mathrm{v}>380`$ observed is in sharp conflict with this argument since Galactic HMXBs have very luminous optical counterpart with $`f_\mathrm{x}/f_\mathrm{v}\stackrel{<}{}1`$ (Bhattacharya & van den Heuvel, 1991). Again, this puzzle may be removed if the black hole is isolated but accreting from the surrounding disk of the supernova fallback matter. The optical emission in this case mainly comes from the outer part of the disk, with a flux much smaller than that in X-ray, as in low-mass X-ray binaries. A simple estimate from the standard $`\alpha `$-disk model (Shakura & Sunyaev, 1973) gives the optical luminosity of the disk to be

$$L_{\mathrm{opt}}6\times 10^{35}(\frac{M}{3M_{}})^{2/3}(\frac{\dot{M}}{10^{19}\mathrm{gs}^1})^{2/3}\mathrm{ergs}^1,$$
(1)
where $`M`$ and $`\dot{M}`$ are the mass of the black hole and the mass accretion rate, respectively. This value is smaller than the X-ray luminosity by a factor of $`\stackrel{>}{}1000`$ with $`\dot{M}10^{19}\mathrm{gs}^1`$. Actually the optical luminosity may be somewhat higher than the above value because a part of the hard radiation flux of the inner disk is reradiated by the X-ray heated, outer disk. But in this case the influence is mainly limited to infrared radiation, since X-ray irradiation becomes important at $`R>10^{12}`$ cm (Sansuichi, Yamada, & Fukue, 1993).

All of the above evidence points to the fact that some of the ULXs are likely to be accreting black holes in isolation, or in binary systems but with no or little mass transfer between the two stars. The source of the accreting material results from the fallback matter during the (failed) supernovae that produce the black holes. From the work of Mineshige et al. (1997), we can roughly estimate the mass accretion rate onto the black hole as,

$$\dot{M}2.7\times 10^{19}(\frac{M}{3M_{}})(\frac{\mathrm{\Delta }M}{0.5M_{}})(\frac{\alpha }{0.01})^{1.35}(\frac{t}{10^5\mathrm{yr}})^{1.35}\mathrm{gs}^1,$$
(2)
where $`\mathrm{\Delta }M`$ is the amount of fallback material and $`\alpha `$ the viscosity parameter. It is clearly seen that the black holes can be ultraluminous, super-Eddington, X-ray sources with an X-ray lifetime up to $`10^5`$ yr for typical values of the parameters, if their X-ray emission is mildly or relativistically beamed. The birth rate of hypernovae (Paczynski, 1998), as the observational consequences of the collapsar model, could be $`10^3`$ yr-1 in the Galaxy (Hansen, 1999), compared to a few $`10^5`$ yr-1 for HMXBs. The ultraluminous X-ray lifetime of the latter may last for several $`10^7`$ yr (Podsiadlowski, Rappaport, & Han, 2003). From this one may expect that a considerable ($`\stackrel{>}{}10\%`$) fraction of the ”young” ULXs may produce X-ray emission with the fallback disk accretion. This provides observational support to the prediction of the collapsar model.

## 3 Discussion and conclusions

We have argued that fallback disk accretion may play an important role in producing the X-ray emission of some ULXs that possibly associated with SNRs or star formation region. This association rules out mass transfer in binaries as the main energy source, since there is little time for the donor stars (if they exist) to evolve and transfer enough mass to the black holes.

The discovery of the optical counterpart to ULX NGC 3031 X-11 is not in conflict with our argument, which was identified as a 23 $`M_{}`$ O8V star (Liu, Bregman, & Seitzer, 2002). For this ULX system, since the companion is a main-sequence star, its stellar wind should not be so strong to power the X-ray emission, suggesting that the black hole may be accreting gas via Roche lobe overflow. But both observations and binary evolution calculations have shown that, it is much more likely that the secondary is evolved during mass transfer in luminous HMXBs (Bhattacharya & van den Heuvel, 1991; Podsiadlowski, Rappaport, & Han, 2003). In this case accretion from the fallback disk could be more important than mass transfer between the two stars.

If some of the ULXs are really related to collapsars (and hence GRBs), then the post-GRB objects should possess a fallback disk with a high mass accretion rate, and emit intense high-energy radiation. However, no such emission has been observed so far from these objects. The reason is generally attributed to photon trapping during hypercritical accretion (e.g. Mineshige et al., 1997). But more recent slim disk calculations with photon trapping effect taken into account (Ohsuga et al., 2002) show that the disk luminosity could be kept around the Eddington luminosity, even if the mass accretion rate is extremely high. Perhaps not all collapsars can produce GRBs. According to Narayan, Piran, & Kumar (2001), the accretion flow in a putative GRB central engine can be either a convection-dominated accretion flow (CDAF) if the accreting mass is introduced at a large outer radius $`R_{\mathrm{out}}`$, or a neutrino-dominated accretion flow (NDAF) if $`R_{\mathrm{out}}`$ is relatively small, i.e., $`R_{\mathrm{out}}\stackrel{<}{}10100R_\mathrm{g}`$, where $`R_\mathrm{g}=2GM/c^2`$ is the Schwarzschild radius (see also Kohri & Mineshige, 2002). Both types of accretion flows are radiatively inefficient, but the latter is much more plausible for a GRB engine since the mass loss during accretion is much smaller than in the former case. Kohri & Mineshige (2002) have studied the properties of an NDAF with mass of $`1M_{}`$ and size of $`5\times 10^6`$ cm. From their work we can estimate the disk luminosity to be as low as $`10^{31}10^{32}\mathrm{erg}\mathrm{s}^1`$. On the other hand, when the disk evolves to extend to larger radius (or the disk size is initially relatively large), the accretion may occur as a CDAF (Narayan, Piran, & Kumar, 2001). It is not yet clear how bright the CDAF is. But its radiation luminosity should also be low because the viscous heating rate is very small in a CDAF since there is no shear stress (Narayan, Igumenshchev, & Abramowicz, 2000). So in both cases we may expect that the collapsar remnants are dim sources at their early stage. Full multidimensional MHD simulations coupled with radiation hydrodynamics are necessary to settle this issue.

ULXs were also found in elliptical galaxies (Sarazin, Irwin, & Bregman, 2001; Colbert et al., 2003). It has been suggested that these ULXs may be soft X-ray transients, i.e., wide low-mass X-ray binaries, at their high (outburst) state (King, 2002). Recent analysis by Colbert et al. (2003) shows that more than 20% of all ULXs found in spirals originate from the older (pop II) stellar populations, indicating that many of the ULXs that have been found in spirals galaxies are in fact pop II ULXs, like those in elliptical galaxies. This would require a fairly large birth rate of ULXs (see King (2001) for a discussion), which is not easily accounted for by wide low-mass X-ray binaries with an estimated birth rate of only $`10^910^8`$ yr-1 in the Galaxy (Kulkarni & Narayan, 1988). We instead suggest that part of these ULXs may be results of the mergers of compact star binaries (e.g., neutron star-neutron star binaries, black hole-white dwarf binaries, black hole-neutron star binaries). Here the accretion disks originate from the disrupted, less massive stars. But the accretion processes are likely to be similar as in the collapsar model - after the initial hypercritical accretion (probably also giving rise to GRBs), the accretion rate decreases to be around the Eddington accretion rate, and the accreting compact stars (most probably black holes) eventually appear as ULXs. The event rate of these mergers may be high enough (Narayan, Paczynski, & Piran, 1992; Fryer, 1999) to account for pop II ULXs.

I am grateful to the anonymous referee for insightful comments. This work was supported by NSFC and NKBRSF.