text:
# Equation of State of Supercooled Water from the Sedimentation Profile

## I Introduction

The physics of the liquid state has been the subject of intense research activity. Novel approaches and novel apparatus have made possible the study of liquids under extreme thermodynamic conditions and in wide windows of space and time. Interesting and unexpected new phenomena have emerged as a result of a combined effort involving experiments, theory and simulationsDebenedetti ; MarieClaire . One of these is the possibility of a liquid-liquid (LL) transition in one-component systems—in addition to the usual liquid-gas transition. Several liquids Katayama ; Ree ; Togaya  and several models StellHemmer ; GF ; Reza ; Swendsen ; Jagla ; Pablo ; Saika ; Lacks  have been studied in detail, and it appears that the class of materials where a LL transition can be observed is larger than the class of tetrahedral liquids that were originally considered Angell00  possible candidates for a LL transition.

Water is one of the liquids that might possess a LL transition. Indeed, the first conjecture of a LL transition was based on a numerical study of the ST2 potential Poole , a model designed to mimic the behavior of liquid water. In the case of ST2, the two coexisting phases differ in their local structure. The low density phase is formed by an open tetrahedrally-coordinated network of hydrogen bonds, while the high density phase has a more distorted network of hydrogen bonds. Recent theoretical work has shown that the interplay between local energy, entropy, and volume which may generate a LL transition can in principle be realized by spherically-symmetric potentials GF ; Jagla .

The evaluation of the $`P(\rho ,T)`$ equation of state (EOS) is key to test for phase coexistence (gas-liquid, liquid-crystal and LL), where $`P`$, $`\rho `$, and $`T`$ denote the pressure, density, and temperature. The numerical calculation of $`P(\rho ,T_0)`$ requires the study of the model for a variety of state points. The coexistence between two phases, in the appropriate temperature window, appears as a region of density values where $`P`$ is constant. In small-size numerical simulations is sometimes hard to observe phase transitions directly in one simulation box, in part because the free energy associated with creating an interface often stabilizes metastable phasesPanagiotopoulous . In these cases $`P(\rho ,T)`$ does not show any flat region negative .

A different approach for studying in one single numerical simulation an entire isotherm has been proposed in Ref. Barrat  (and later exploited in the experimental study of colloidal systems and in the study of crystallization profiles Ref. Piazza ; Lowen ; Allain ; Chaikin ). This approach simulates a semi-infinite tube in the presence of a very strong gravitational field and measures the density profile in equilibrium. A simple inversion of the density profile allows the model EOS to be constructed. This idea has also been applied in the experimental study of the EOS of colloidal particles, by inversion of the measured sedimentation equilibrium profile.

Here we apply the sedimentation profile method to ST2 water. We find that an interface separating two liquid states appears at low temperature, corresponding to the coexistence of two metastable liquid states of water within one simulation box, providing evidence for the presence of an LL transition.

## II Theory and Simulation Details

We study $`N=7680`$ rigid molecules of mass $`M`$ interacting through the ST2 water potential, a rigid, non-polarizable, 5-site potential Stillinger  that is able to reproduce qualitatively the thermodynamic anomalies of liquid water. “ST2 water” is characterized, on cooling, by isobaric density maxima, increasing compressibility, increasing constant $`P`$ specific heat, and evidence for a LL transition in the deeply supercooled regime (which is difficult to probe experimentally due to spontaneous crystallization Poole ).

To implement the sedimentation profile method, we use a column-shaped simulation box, semi-infinite along the $`z`$ axis and with periodic boundary conditions along the $`x`$ and $`y`$ axes. The top of the box is left open, while the bottom is assumed to be a repulsive soft-sphere surface, generating a short-range force proportional to $`z^{13}`$ acting on the molecule center of mass (Fig. 1). The box width in the $`x`$-axis and $`y`$-axis directions is 3 nm, corresponding to a bottom surface area of $`S_{xy}=9`$ nm2. A strong gravitational field $`g_s=2\times 10^{12}g`$ is applied downward, in the $`z`$-axis direction, where $`g=9.8`$ kg/m2 is the Earth’s gravitational field. The value of the field controls the range of $`P`$-values accessed in the simulations. The pressure at the bottom of the column is $`P=g_sNM/S_{xy}500`$MPa.

We perform the simulations using a multi-processor code on SGI Origin 2000, IBM SP, and IBM Regatta supercomputers. We choose a 1 fs time step and study different temperatures from $`T=300`$ K down to $`T=230`$ K. For T=$`230`$ K, we simulate four different systems to better estimate error and reproducibility of the results. A long ($`10ns`$) equilibration time precedes the actual calculation of the equilibrium density profile. To analyze equilibration we monitored the running average of the $`z`$-position of the center of mass. Production runs lasted at least 200 ps for the higher temperatures up to several ns per box when $`T=230`$ K.

We define a mass density field $`\rho (z)`$ by averaging the density over bins with height $`\mathrm{\Delta }z=1`$ nm to calculate the pressure field $`P(z)`$ from

$$P(z)=g_z^{\mathrm{}}\rho (z^{})𝑑z^{}.$$
(1)
A parametric plot of $`P(z)`$ vs. $`\rho (z)`$ provides the EOS.

## III Results and Discussion

Figure 2 shows the equilibrium density profile along the $`z`$-axis for three different temperatures. Each symbol represents one bin, where the height of each bin $`\mathrm{\Delta }z=1`$ nm. For $`T=300`$ K we see only one break in the density profile: the topmost points in the plots correspond to the gas-liquid interface. Above these points is a much less dense gas not shown on this density scale. For $`T=250`$ K, there appears an inflection. For $`T=230`$ K, we see a clear break in the density, associated with the interface between two different liquid states.

Figure 3 shows the corresponding $`P(\rho )`$ relations at the same three temperatures. We also show $`P(\rho )`$ as evaluated previously using standard MD for cubic boxes with periodic boundary conditions for systems with $`N=6^3=216`$ (Ref. peter ) and $`N=12^3=1768`$ (Ref.steve ). In Refs peter  and steve  $`P`$ was calculated with the standard virial relationnegative2 .

Consistent with the data shown in Fig. 2, we see a region of coexistence between two different phases at $`T=230K`$, in agreement with the estimated location of the LL transition in Ref. Sciortino . Note that when using the sedimentation equilibration method, we see no unphysical loops in the equation of state, unlike cubic box simulations where the boundary conditions may artificially stabilize metastable states.

To confirm that the two coexisting phases are both liquids, we calculate the mean square displacement (MSD) for different height values, both below and above the interface. We follow the evolution of each molecule for an average mean squared displacement smaller that $`\mathrm{\Delta }z`$, so that each height value can be unambiguously assigned to an average density value. Figure 4 shows that both phases are sufficiently diffusive. The low density phase has smaller diffusivity, in agreement with previous simulations of the density dependence of the dynamics Geiger .

## IV Conclusions

We studied the EOS of “ST2 water” with the sedimentation profile method, and presented evidence that the EOS at $`T=230K`$ shows a clear phase coexistence, between two phases which are both liquids. The observed LL coexistence in the ST2 potential phase diagram is consistent with estimate of the “critical point” being located around $`T=235`$ K, $`P=250`$ MPa, and $`\rho =1.05`$ g/cm3. Note that for most simple water models, the temperature and pressure scales are shifted relative to real values, thus they must be shifted to place them within an experimental contextcompare . The critical point estimated experimentally is around $`T_c=230`$ K and $`P_c=50`$ MPa with $`\rho =1.05`$ g/cm3Mishima .

## Acknowledgments

This work is supported by NSF grant CHE-0096892. MY is also supported by NSF grant GER-9452651. FS thanks INFM-Iniziative Parallel Computing for granting numerical resources.


format: text

text:
# On the origin of plankton patchiness

Center for Studies in Physics and Biology, The Rockefeller University, 1230 York Avenue, New York, New York 10021, USA

Department of Physics, FEN, Universitat Politécnica de Catalunya, Campus Nord, Mòdul B4, E-08034 Barcelona, Spain

Departament de Física Fonamental, Facultat de Física, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Spain

Keywords: complex ecosystems, plankton patchiness, prey-predator models, randomness, turbulent diffusion, scales.

PACS: 87.23.Cc, 05.40.-a, 87.10.+e, 05.45.-a

Abstract

Plankton is the productive base of aquatic ecosystems and plays a major role in the global control of atmospheric carbon dioxide. Nevertheless, after intensive study, the factors that drive its spatial distribution are still far from being clear. The models proposed so far show very limited agreement with actual data as many of their results are not consistent with field observations. Here we show that fluctuations and turbulent diffusion in standard prey-predator models are able to accurately and consistently explain plankton field observations at mesoscales (1-100 km). This includes not only the spatial pattern but also its temporal evolution. We explicitly elucidate the interplay between physical and biological factors, suggesting that the form in which small scale biotic fluctuations are transferred to larger scales may constitute one of the key elements in determining the spatial distribution of plankton in the sea.

Understanding how complex ecosystems work often relies on simplified models that disregard many details of the actual system while retaining the essential information . In the case of marine ecosystems, not only the simplest approaches failed to explain the spatial distribution of plankton populations but also more sophisticated models were unable to account consistently for the most remarkable features . Even now, it is not clearly understood why under apparent homogeneous conditions, such as temperature and nutrients, plankton is still patchilly distributed. This particular situation is of remarkable importance since in the absence of external sources of patchiness the pattern must arise as a mere consequence of the interactions between the individuals. The typical form of these patterns is illustrated in Fig. 1. The main trait is that zooplankton is more patchilly distributed than phytoplankton .

The most intuitive model that can be proposed to explain plankton dynamics  considers the population densities of prey (phytoplankton), $`N`$, and predators (zooplankton), $`P`$:

$$\frac{N}{t}=F_N(N,P)+D_N^2N,$$
(1)
$$\frac{P}{t}=F_P(N,P)+D_P^2P,$$
(2)
where $`D_N`$ and $`D_P`$ are diffusion coefficients; and $`F_N`$ and $`F_P`$ are functions that account for the interaction between both species. This class of models is the most frequently used in theories on pattern formation in ecology . In the case of plankton, they were able to display spatial heterogeneity under homogeneous conditions . Zooplankton, however, was less patchilly distributed than phytoplankton, in contradiction with the observed pattern .

There are two relevant features that are not taken into account by this kind of models. First, diffusion in the sea is not quantitatively well modeled by usual Fickian diffusion . Both types of diffusion processes will tend to spread and mix the populations, but the specific form in which this is achieved is different. Second, there is always an intrinsic stochasticity associated with the dynamics of the population . From birth to death, all processes share some degree of chance. The way randomness manifests in the dynamics of the individuals depends on the scale we are looking at ; deterministic equations are expected to be valid in the limit of high numbers of individuals . Therefore, a deterministic description may be a reasonable one for phytoplankton alone, but this does not need to be so for zooplankton which has much fewer individuals . More importantly, while phytoplankton interacts mainly with zooplankton, zooplankton interacts also with fish and whales which are far from being evenly distributed.

These two additional features have been incorporated in a prey-predator model:

$$\frac{N}{t}=F_N(N,P)\stackrel{}{v}\stackrel{}{}N,$$
(3)
$$\frac{P}{t}=F_P(N,P)\stackrel{}{v}\stackrel{}{}P+\xi (t),$$
(4)
where dispersal is given by advection with a velocity field $`\stackrel{}{v}`$ \[$`\stackrel{}{v}(\stackrel{}{r})`$\] that depends on the position $`\stackrel{}{r}`$, and where a noise term $`\xi (t)`$ has been included.

In general, the effects of the advective terms depend on the precise form of the velocity field. For some turbulent fields , the effect of advection can be simplified as follows: given a passive field $`f(\stackrel{}{r},t)`$ which evolves as

$$\frac{f}{t}=\stackrel{}{v}\stackrel{}{}f.$$
(5)
the spatial Fourier transform of $`f(\stackrel{}{r},t)`$ follows from

$$\frac{df_k}{dt}=D|k|^\beta f_k,$$
(6)
where $`k`$ is the wave number and $`D`$ a constant. In this case, advection can effectively be viewed as a diffusion process with a diffusion coefficient $`D_{eff}(k)=D|k|^{2\beta }`$ that depends on the scale. In contrast to usual Fickian diffusion, the variance of the field is not proportional to $`t`$ but is given by $`r^2t^{2/\beta }`$. This is the type of time dependence observed for the dispersion of tracers in the sea , from which one can obtain the explicit value of the parameter $`\beta `$.

To render our model analytically tractable, we consider the system around a stable state. Fluctuations in zooplankton, $`\xi `$, move the system away from equilibrium. If the fluctuations are not too large, we can perform a linear expansion of $`F_P`$ and $`F_N`$:

$`F_N(N,P)`$ $`=`$ $`c_Na_{11}Na_{12}P,`$
$`F_P(N,P)`$ $`=`$ $`c_P+a_{21}Na_{22}P.`$ (7)
Here $`c_N`$, $`c_P`$, $`a_{11}`$, $`a_{12}`$, $`a_{21}`$, and $`a_{22}`$ are positive constants. For the simplest form of the noise term , Gaussian white and uncorrelated in space \[$`\xi (\stackrel{}{r},t)=0`$ and $`\xi (\stackrel{}{r},t)\xi (\stackrel{}{r}^{},t^{})=2\sigma ^2\delta (\stackrel{}{r}^{}\stackrel{}{r})\delta (t^{}t)`$\], and for $`a_{22}0`$ the variance spectra are given by

$`S_N(k)={\displaystyle \frac{a_{12}^2\sigma ^2}{(\stackrel{~}{D}_N+\stackrel{~}{D}_P)\stackrel{~}{D}_N\stackrel{~}{D}_P}}\text{and}S_P(k)`$ $`=`$ $`{\displaystyle \frac{\sigma ^2}{\stackrel{~}{D}_P}}.`$ (8)
where $`\stackrel{~}{D}_ND|k|^\beta +a_{11}`$, $`\stackrel{~}{D}_PD|k|^\beta +a_{22}`$, and $`\sigma ^2`$ is the intensity of the noise source. The assumptions involved do not substantially constrain the applicability of the results. When $`a_{22}`$ is not negligible, the expressions become more involved but the qualitative behavior is still the same. In particular, the high wave-number limit remains unchanged. On the other hand, the type of noise we have considered is quite general and can arise, among others, from a random distribution of predators feeding on zooplankton or even from the birth process itself . Other types of noise with different properties, e.g. as those induced by turbulence , are certainly present but we assume that they are not relevant for the spectral properties of the pattern at the mesoscales.

The variance spectra obtained from previous equations display a power-law region with exponent $`3\beta `$ for the phytoplankton and $`\beta `$ for the zooplankton. The value of $`\beta 0.87`$ inferred from diffusion in the sea  leads to exponents $`2.6`$ and $`0.8`$, both in excellent agreement with field data . It is worth emphasizing that the power law behavior appears only for sufficiently high wavenumber (short scales); for low wavenumbers (long scales), the variance spectra is flat, as observed in most field data .

In the same way, one can compute the coherence between two patterns at different times , which provides information about the global dynamics. For the phytoplankton this quantity is given by

$$\frac{N(k,\mathrm{\Delta }t)}{N(k,0)}=\frac{\stackrel{~}{D}_Pe^{\stackrel{~}{D}_N\mathrm{\Delta }t}\stackrel{~}{D}_Ne^{\stackrel{~}{D}_P\mathrm{\Delta }t}}{\stackrel{~}{D}_P\stackrel{~}{D}_N},$$
(9)
where $`N(k,\mathrm{\Delta }t)=_0^{\mathrm{}}N(k,\omega )^2e^{i\omega \mathrm{\Delta }t}𝑑\omega `$. This result indicates that short scales lose their correlation faster than long ones and that eventually the whole pattern will be decorrelated, as observed in satellite measurements .

In Fig. 2 we plot the typical form of the variance spectra and the squared coherence for different time lags. Both of them are in excellent agreement with field data . Remarkably, the main properties of the pattern already appear in the linear regime. Therefore, nonlinear interactions that drive the system towards a stable state will lead to similar results. To study this aspect in more detail, we have performed numerical simulations for typical nonlinear interactions as explained in the caption of Fig. 3. The resulting two-dimensional spatial distribution, transects, and variance spectra (shown in Fig. 3) agree with both the linear model and field data. Other types of nonlinear interactions — e.g. different functional responses — as well as different types of noise — e.g. acting on zooplankton growth rate — also produce similar results (data not shown).

Field observations indicate that the power law region of the variance spectra and the value of the exponent of this power law are robust properties of the system; i.e. these properties are present under a wide variety of conditions. In our model, there is always a power law region whose exponent does not depend on biotic factors but is completely determined by the specific form in which turbulent diffusion acts on the system. This provides a straightforward explanation of the predominance of the observed exponents for the phytoplankton falling between $`3`$ and $`2`$. These are the values that arise for 2D ($`\beta =1`$) and 3D ($`\beta =2/3`$) isotropic turbulence, respectively . In the sea, the value of this exponent will depend on the particular situation, but it is reasonable to assume that it will be between those of 2D and 3D isotropic turbulence, as the available data shows . There are also non-robust properties, such as the region where the variance spectra turns flat. In the model, this depends on many factors: e.g. growth and death rate, and turbulence. Field data shows that, indeed, the position of this region exhibits great variability and that sometimes it is not even present in the range of scales observed.

Turbulent diffusion and noise are two obvious features that have already been considered in the context of marine ecosystems, but none of them by itself has been able to explain the mesoscale patterns. In particular, it is well known that noise generates variability, i.e. that noise can be a source of patchiness . For instance, reaction-diffusion prey-predator models with noise produce patterns that at a glance strongly resemble those observed in the sea . The exponents obtained ($`6`$ for the phytoplankton and $`2`$ for the zooplankton), however, are far from the observed ones. This quantitative, but not qualitative, disagreement is due to the dependence of the effective diffusion coefficient with the scale. Thus, reaction-diffusion models are unable to integrate correctly the scale dependence of the physical properties of the environment. When this is taken into account, noise not only generates patterns but is also able to produce the right ones.

Turbulence plays a somehow ambivalent role. It can act in the same way as diffusion does (transferring variance from smaller to larger scales) and also in the opposite way (from larger to smaller scales). These two types of processes are referred to as turbulent diffusion and turbulent stirring, respectively. The former is the one we have considered in our model. It was already considered in Ref.  together with the type of noise that turbulence induces but without the noise that can arise from biotic factors. The latter only plays an important role when some degree of environmental heterogeneity is present . Indeed, it has been shown that turbulent stirring can generate patterns that resemble the observed ones if spatial heterogeneity and time delays are introduced in the model . The type of time delays introduced, however, can lead non-realistic situations such as growing zooplankton in the absence of phytoplankton.

Finally, it is important to emphasize that in our case noise is the key element that allows moving from the individual to the population description. Our results suggest that zooplankton dynamics at lower scales affects the pattern at the mesoscale in the same way as noise does. Considering a more detailed description is not necessary to explain and to understand the main characteristics of the pattern. This does not mean that the actual dynamics of zooplankton is not important at all: its growth rate, its survival, and the intensity of noise itself depend, among other factors, on how zooplankton aggregates and on how it avoids its predators . It rather means that under a wide range of conditions all those intricate mechanisms will lead to patterns with properties as those induced by noise.

Acknowledgments
This work was supported in part by the DGICYT of the Spanish Government.

Figure Captions

* Transects obtained from field observations for (a) phytoplankton and (b) zooplankton (redrawn from Ref. ). Here Chlorophyll a is a measure of phytoplankton activity.
* (a) Variance spectra $`S(k)`$ for prey and predators from Eqs. (1), (2), and (7). (b) Square coherence $`[N(k,\mathrm{\Delta }t)/N(k,0)]^2`$ for prey \[Eq. (9)\] for time lags of $`1`$, $`6,`$ and $`7`$ days. The values of the parameters are $`D=12`$, $`a_{11}=0.3`$, $`a_{22}=0.05`$, and $`\beta =1`$. The length and time units are Km and days, respectively.
* Numerical simulations for a prey-predator model given by Eqs. (1) and (2), with $`F_N(N,P)rN\left(1N/K\right)cPf(N)`$ and $`F_P(N,P)=P\left(gf(N)ϵ\right)`$. The functional response is $`f(N)=N^2/(1+N^2)`$. Here, $`K`$, $`c`$, $`g`$, $`r`$ and $`ϵ`$ are positive constants. The velocity field $`\stackrel{}{v}\stackrel{}{v}(\stackrel{}{r})`$ consists of a series of vortices distributed as in Ref.  and the noise is assumed to be Gaussian with zero mean and correlation function $`\xi (\stackrel{}{r},t)\xi (\stackrel{}{r}^{},t^{})=2[\sigma P(\stackrel{}{r},t)]^2\delta (\stackrel{}{r}^{}\stackrel{}{r})\delta (t^{}t)`$. Typical transects for (a) prey and (b) predators. (c) Variance spectra for prey and predators. These results were obtained by discretizing the corresponding equations on a $`250\times 250`$ two-dimensional mesh , with periodic boundary conditions and then by using a standard method for integrating stochastic differential equations . The values of the parameters are $`r=0.3`$, $`K=4`$, $`c=2`$, $`g=0.1`$, $`\epsilon =0.05`$, and $`\sigma =3.5`$. The length and time units are Km and days, respectively. The size of each cell of the discretized mesh is $`0.25\times 0.25`$ Km2.

FIGURE 1

FIGURE 2

FIGURE 3


text:
# Thermal Model Description of Strangeness Enhancement at Mid-Rapidity in Pb–Pb collisions at 158 GeV A/c.

## I Introduction

The WA97 collaborationAntinori et al. (1999) at CERN has shown that strange particle yields per wounded nucleon reach a saturation level for the most central Pb-Pb collisions (the Pb-Pb data is presented in 4 centrality bins) and show a pronounced increase when compared to a p-Be system. A mixed canonical description of these data has been published by Hamieh et al.Hamieh et al. (2000); Redlich (2002) who showed that the density of strange particles derived from such a treatment depends on the size of the system consistent with the WA97 observationsAntinori et al. (1999). In addition Hamieh et al. Hamieh et al. (2000) have shown that if reasonable values are used for the thermal parameters of each system the behavior of the yield (in the Pb-Pb systems) per wounded nucleon normalized to the p-Be yield is in agreement with the values obtained by the WA97 collaboration. The present work uses also he mixed canonical ensemble Braun-Munzinger et al. (2002), and presents a more thorough analysis in that the thermal parameters for each system are obtained by minimizing $`\chi ^2`$ and a more explicit comparison between model and data is performed. The dependence of the parameters on the size of the system is presented, as well as the yields per wounded nucleon relative to p-Be.

## II Mixed Canonical Formalism

The mixed canonical partition function for a hadron gas conserving strangeness exactly is described by the following partition function Cleymans et al. (1997a):

$$Z_{S=0}^C=\frac{1}{2\pi }_\pi ^\pi 𝑑\varphi \mathrm{exp}\left(\underset{n=3}{\overset{3}{}}S_ne^{in\varphi }\right),$$
(1)
where $`S_n=V_kZ_k^1`$, $`V`$ is the volume, and the sum is over all particles and resonances carrying strangeness $`n`$. For a particle of mass $`m_k`$, with spin-isospin degeneracy factor $`g_k`$, carrying baryon number $`B_k`$ and charge $`Q_k`$ with baryon chemical potential $`\mu _B`$ and charge chemical potential $`\mu _Q`$, the one-particle partition function is expressed in the Boltzmann approximation as:

$$Z_k^1\frac{g_k}{2\pi ^2}m_k^2TK_2\left(\frac{m_k}{T}\right)\mathrm{exp}(B_k\mu _B+Q_k\mu _Q).$$
(2)
As described in Ref. Braun-Munzinger et al. (2002) the partition function for this ensemble may be rewritten in the following form which is well suited for numerical calculation:

$`Z_{S=0}^C=e^{S_0}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p=\mathrm{}}{\overset{\mathrm{}}{}}}a_3^pa_2^na_1^{2n3p}`$
$`\times I_n(x_2)I_p(x_3)I_{2n3p}(x_1),`$ (3)
where

$`a_i`$ $`=`$ $`\sqrt{S_i/S_i}`$ (4)
$`x_i`$ $`=`$ $`2\sqrt{S_iS_i}`$ (5)
and $`I_i`$ are modified Bessel functions.

The expression for the particle density, $`n_i`$, may be obtained from the partition function Eq. (3) by following the standard method Cleymans et al. (1997a). For a particle $`i`$ having strangeness $`s`$ the result is

$`n_i={\displaystyle \frac{Z_i^1}{Z_{S=0}^C}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p=\mathrm{}}{\overset{\mathrm{}}{}}}a_3^pa_2^na_1^{2n3ps}`$
$`\times I_n(x_2)I_p(x_3)I_{2n3ps}(x_1).`$ (6)
## III Results

In this section we explore whether canonical strangeness suppression at small volumes (compared to grand canonical equilibrium particle yields) is able to explain the enhancement of strange particle yields per wounded nucleon from small to large systems, as measured by the WA97 collaboration at CERNAntinori et al. (1999). The thermal model is well suited for to $`4\pi `$ integrated data and to the central rapidity region if a boost invariant plateau exists around at mid-rapidityCleymans and Redlich (1999). We assume the validity of the latter in order to analyze the particle yields. A thorough comparison of $`4\pi `$ particle numbers with the thermal model has been made in Ref. Becattini et al. (2001). As has been the case in other calculations Braun-Munzinger and Stachel (1996, 1999), our findings at mid-rapidity are found to be consistent with full strangeness equilibration (i.e. $`\gamma _s=1`$). This is in contrast to $`4\pi `$ integrated particle yields which are incompatible with full strangeness chemical equilibriun and deviate from it by several standard deviations. Experimentally, ratios of particle to anti-particle yields111The reason for checking particle to anti-particle ratios is that the masses of the particles in the ratio are equal - this leads to these particles being equally affected by flow. This in turn minimizes the errors introduced by considering a limited kinematic region Yen and Gorenstein (1999) at mid-rapidity222In symmetric collisions the majority of new particles are produced at mid-rapidity. have been compared to $`4\pi `$ integrated yields and shown to be in agreement for S+S collisions at $`200`$ A GeV by the CERN NA35 collaboration Sollfrank et al. (1994). In addition, and perhaps of greater interest, the $`4\pi `$ integrated ratio $`\overline{\mathrm{\Xi }}/\mathrm{\Xi }`$ measured by CERN NA49 for Pb+Pb collisions at 158 A GeV has been shown to agree with the corresponding CERN WA97 mid-rapidity ratio Barton et al. (2001).

As the mixed canonical formalism has been derived here for the case of Boltzmann statistics, it is worth noting that the second term in a series expansion of the correct quantum statistical distribution function for kaons gives a correction to the kaon numbers of the order of $`3\%`$ at a temperature of 150 MeV. The corrections to the Boltzmann distribution functions due to quantum mechanics are expected to affect the kaons more than any other strange particle, as they are the lightest of the strange particles. With errors of $`3\%`$ and less, the use of Boltzmann statistics to describe the strange particles is justified. For all mixed canonical analyses the parameter $`\mu _Q`$ fit to zero, and has been subsequently removed as a free parameter. For a $`4\pi `$ integrated $`\mathrm{Pb}+\mathrm{Pb}`$ system, $`\mu _Q`$ is expected to be small and negative.

### III.1 Analysis A

The main results of analysis A are based on fitting the ratios of the strange particle yields to the yield of negatives for each case in Table 1. The negatives yield is used for each ratio due to the good statistics. In addition it provides a measure of the entropy of the collision. For negative pions, the quantum mechanical correction to the Boltzmann distribution function at 150 MeV is of the order of $`30\%`$ . It is unacceptable to ignore this correction when attempting to describe the data. Considering that the non-strange particle yields predicted by the model are independent of the exact conservation of strangeness, the full quantum mechanical grand canonical particle number expression is evaluated for all non-strange particles. The parameters obtained in this way are shown in Table 1 as analysis A333The table also includes the $`\chi ^2`$ and fit parameters for two other analyses described in detail below. The $`\chi ^2`$ for each system is good, and there are three (two) degrees of freedom for the p+Pb, and Pb+Pb systems (p+Be system). The large uncertainties in the radius parameter for the Pb+Pb systems shows that we have hit the grand canonical limit, and have no volume dependence, this is discussed briefly in a later section.

### III.2 Analysis B

Analysis B has only three data points ($`\frac{\overline{\mathrm{\Lambda }}}{\mathrm{\Lambda }},\frac{\overline{\mathrm{\Xi }}}{\mathrm{\Xi }}`$ and $`\frac{\overline{\mathrm{\Xi }}}{\overline{\mathrm{\Lambda }}}`$) for each fit. This analysis has been motivated by the reasons mentioned previously in the text. This analysis has no free parameters and shows reasonable agreement with analysis A. The large uncertainties of the fitted parameters are unavoidable. Interestingly, fitting these ratios proves impossible for the p+Pb system. The hypothesis of Hamieh et al. (2000) of a special ‘interaction volume’ for strange particles in this system cannot be tested, as the particle ratios would depend only on this interaction volume, as the regular volume cancels in the ratio of the particle multiplicities.

### III.3 Analysis C

Analysis C determines the fits to the data using the grand canonical ensemble with quantum statistics. In this case, the volume dependence canels out in the particle ratios and we have included $`\mu _S`$ and $`\mu _Q`$. This procedure fails for the p+Be system. As may be expected for a system of this size, a canonical treatment Cleymans et al. (1997b) is required. In this analysis, the parameters fitted for Bin 4 immediately catch the eye – with high uncertainty values even though $`\chi ^2<1`$. This is especially puzzling when one notices that this is not the case in analyses A or B. The $`T`$ and $`\mu _B`$ parameters agree with those obtained by Becattini et al. Becattini et al. (2001), for a grand canonical description of the $`4\pi `$-integrated data for central symmetric $`\mathrm{Pb}^{208}`$ collisions at 158 A GeV by the CERN NA49 collaboration.

The parameters in analyses A, B, and C agree within th each other one standard deviation.

### III.4 Thermal Model Parameters

In order to fully reproduce the data obtained by the WA97 collaboration , it is necessary to determine a relationship between the parameters $`\mu _B`$, the radius $`R`$ and the temperature $`T`$ obtained from the fits of analysis A and the number of wounded nucleons. The fits used are shown in Figs. 1 (temperature), and 2 ($`\mu _B`$), as well as the parameters obtained from a full Boltzmann treatment444Using Boltzmann statistics for strange, and non-strange particles. The volume paramaterisation is described by equation 7. The parameters obtained from the p+Pb data have not been considered because of their prediction of a volume smaller than that of the p+Be system.

#### III.4.1 Temperature

As may be seen in Fig. 1, it appears initially that the temperature increases with centrality. However, the last data point does not fit this hypothesis, and the variation of chemical freeze-out temperature with number of participants has been assumed to be constant ($`T=163`$ MeV). This is in agreement with the fit to the NA49 data by Becattini et al. Becattini et al. (2001) to central Pb+Pb data. This temperature value is slightly less than the one obtained by Hamieh et al. Hamieh et al. (2000).

#### III.4.2 Radius

In order to reproduce the WA97 data, knowledge of the variation of system size with average number of wounded nucleons is required. We have chosen a function of the form:

$$<N_{\mathrm{wound}}>=R^3+b$$
(7)
with $`b<0.5`$ adjusted to reproduce the data of the WA97 collaboration (Figs. 4 and 3). This differs from the dependence of the radius on $`A_{part}`$ assumed by Hamieh et al. Hamieh et al. (2000) where $`A_{part}1.31.7R^3`$.

The large uncertainty in the radius for the larger systems is not surprising, as the ratios being fitted are expected to show a small volume dependence, due only to canonical strangeness suppression. The effects of canonical strangeness suppression decrease with system size as one approaches the grand canonical limit where ratios of particle multiplicities have no volume dependence.

#### III.4.3 Baryon Chemical Potential

Figure 2 shows the variation of the baryon chemical potential with system size. This function increases rapidly before saturating. All the Pb+Pb bins are described by the saturation value of approximately 210 MeV. The obtained values of $`\mu _B`$ in the central Pb+Pb bin are in agreement with that of Becattini et al. Becattini et al. (2001). The $`\mu `$ used for the Pb+Pb bins by Hamieh et al. Hamieh et al. (2000) is just outside one standard deviation of the values obtained during this analysis. For the p+Be system the value of $`\mu _B`$ used in Ref. Hamieh et al. (2000) (150 MeV) is much larger than that obtained in this analysis ($`111\pm 7`$ MeV) .

The exact variation of $`\mu _B`$, $`T`$ and $`R`$ in the range between the p–Be and Pb–Pb may differ from what is shown here and new data from CERN NA57 in this range are eagerly awaited.

The figures mentioned above (1, and 2) also show the variation of the parameters $`T`$, and $`\mu _B`$ with centrality in the case where all particle multiplicities are assumed to be Boltzmann. These parameters are seen to be in agreement with the parameters from Analysis A. The value of the $`\chi ^2`$ parameter obtained for the purely Boltzmann fits were of order 1 for all systems considered.

### III.5 Predicted Evolution of Particle Yields

Figs. 3 and 4 show the ability of the model to reproduce the data using the functions in Figs 1 and 2, and equation 7 to predict the variation in chemical freeze-out temperature $`T`$, baryon chemical potential $`\mu _B`$, and radius $`R`$ with average number of wounded nucleons. As is clear, the agreement of the model with experimental data is good. The exact shape of the thermal model predictions shown in Figure 3 and 4 show a large dependence on the exact relationship between number of wounded nucleons and system size, and a weaker dependence on the functional fit to the other two parameters. This is rather unfortunate as the thermal model $`R`$ values obtained have uncertainties equal to their magnitude, so this relationship is not constrained by the model.

## IV Conclusions

A thermal model conserving baryon number and charge on average, and strangeness exactly has been used to describe particle yields from heavy ion collisions. The formalism includes all strange particles. The model has been applied to the CERN WA97 data, and is shown to be able to reproduce the data for all centrality classes (Figs. 3 and 4). The accuracy of the model parameters (temperature $`T`$, radius $`R`$, and baryon chemical potential $`\mu _B`$) is restricted by the lack of $`4\pi `$ integrated particle yields for multi-strange particles. In order to reproduce WA97 graphical presentation of their data accurately, some knowledge of the variation of $`T`$, $`R`$, and $`\mu _B`$ was required. A functional form of the variation of these parameters with the number of wounded nucleons has been determined and is presented in Figs 1 and 2, and equation 7.

The parameters obtained in this model for the central Pb+Pb bins are in agreement with the $`4\pi `$ thermal model application of Becattini et al. Becattini et al. (2001). This lends some weight to the accuracy of the model when fitting particle ratios in a limited kinematic region. This grand canonical model includes the factor $`\gamma _S1`$ to predict the yields of strange particles. To differentiate between canonical strangeness suppression, and suppression of strange particles by anomalous phase space occupancy555parameterised by $`\gamma _S1`$ in thermal models, the yield of the $`\varphi `$ particle could be used. Yields of this particle are not sensitive to canonical strangeness suppression, but as it contains an $`s`$ and an $`\overline{s}`$ quark, these yields will be sensitive to $`\gamma _S1`$.

The enhancement of strange particles measured by CERN WA97 has been considered a signal for QGP formation. At first glance, the ability of a full equilibrium thermal model to reproduce the data suggests the existence of a deconfined state, since equilibrium strange particle yields have been proposed as a possible signal for deconfinement Rafelski and Muller (1982). A deconfined phase is not, however, expected to be formed in the p+Be system, as it is only expected in large dense systems. The ability of the model to reproduce the p+Be data suggests that it is possible for strange particles to reach equilibrium yields by hadronic interactions alone.

###### Acknowledgements.

Sahal Yacoob wishes to thank the National Research Foundation (NRF) of South Africa for supporting this work.


text:
# Wave Packets can Factorize Numbers

## I Introduction

Interference of waves instead of Newton’s light corpuscles — what a revolutionary idea of Thomas Young pronounced in his famous talk at the Royal Society in 1801 bib:young . The opposite idea — particles as waves — proposed by Louis de Broglie and cast into the appropriate mathematical form by Erwin Schrödinger more than 100 years later opened a new realm in atomic physics with more applications today than ever. “Can a quantum-mechanical description of physical reality be considered complete?” This famous question raised by Albert Einstein, Boris Podolsky and Nathan Rosen in their seminal paper bib:einstein  brought out most clearly the importance of entanglement bib:schroedinger  as the basic ingredient of quantum theory. Interference, wave mechanics and entanglement are the three ideas on which Peter Shor could build his factorization algorithm bib:shor1  with a polynomial rather than exponential speed. This discovery gave birth to the new field of quantum information bib:books . Guided by the same three central ideas we show in the present paper that the time evolution of elementary quantum systems carries intrinsically the tools to factorize numbers bib:varenna .

Stationary states have played an important role in the development of quantum mechanics. For this reason the time evolution of wave packets has for many years taken a back seat even in textbooks on quantum theory. In this context it is interesting to recall that already in 1927 Earle H. Kennard bib:kennard  studied the time evolution of wave packets in simple potentials bib:nieto . Recently, sophisticated techniques to create short laser pulses have opened a new era for a controlled creation and observation of wave packets in atoms bib:yeazell , molecules bib:vrakking  and cavities bib:eberly . Even wave packet dynamics of cold atoms bib:raithel  and Bose-Einstein condensates bib:leanhardt  has been observed. These examples illustrate convincingly that today’s technology allows to promote this gedankenexperiment of factorization to a real implementation in a quantum system.

The use of wave packet dynamics to factorize numbers has been recognized for special examples. Reference bib:dowling  uses a $`N`$-slit Young interference setup. Hence, the number to be factorized is encoded in the number of slits. The indicator of a factor is the interference pattern on a screen as a function of the wavelength of the light or matter wave. Whenever the intensity across the screen is constant, the appropriately scaled wavelength provides a factor. Similarly, Ref. bib:harter  considers the time evolution of a particle in a box. In contrast to these contributions we provide a general framework for factorization and moreover, exhaust the potential of interference by capitalizing on properties of quadratic phase factors.

Our paper is organized as follows. In Sec. II we briefly review the essential ingredients of wave packet dynamics bib:buch . In particular, we present an approximation for the autocorrelation function that illustrates its ability to factorize numbers. The richness of the autocorrelation function originates from quadratic phases which are also the origin of the Talbot effect bib:rohwedder , revivals and fractional revivals bib:leichtle1 , curlicues bib:curlicues  and quantum carpets bib:physworld . In Sec. III we turn to one of the most elementary quantum systems, namely the quantum rotor. In this case the autocorrelation function is a Gauß sum which has been discussed extensively in number theory bib:lang . We show that under appropriate conditions factors of a number manifest themselves in non-vanishing real or imaginary parts of the autocorrelation function. We explain this phenomenon using the properties of Gauß sums. Our conclusions and outlook presented in Sec. IV summarize the main results and provide an “idea for an idea” how to incorporate entanglement into our scheme.

## II Wave packets

In the present section we briefly review the essential features of wave packet dynamics. We start from the definition of the autocorrelation function and approximate the energy spectrum by one that is quadratic in the quantum number. Physical systems such as Rydberg atoms and diatomic molecules are well described by such a spectrum. We then show the possibility of factorizing a number.

### II.1 Autocorrelation function

In standard wave packet experiments many physical observables of an atom or molecule are accessible. However, one quantity is most appropriate to obtain information about the dynamics and is nowadays measured routinely: It is the so-called autocorrelation function

$$S(t)\psi (t=0)|\psi (t)=\underset{n=0}{\overset{\mathrm{}}{}}|\psi _n|^2e^{iE_nt/\mathrm{}},$$
(1)
which indicates the overlap between the time-evolved wave packet

$$|\psi (t)=\underset{n=0}{\overset{\mathrm{}}{}}\psi _ne^{iE_nt/\mathrm{}}|n$$
(2)
and the initial state. Here, the coefficients $`\psi _n`$ denote the probability amplitudes of the initial state in the basis of energy eigenstates $`|n`$ with corresponding energy eigenvalues $`E_n`$.

While so far the autocorrelation function has served as a tool for observing the dynamics, our goal is to use its intrinsic features to perform computational tasks. In particular, we show that autocorrelation functions can reveal prime factors of an appropriately encoded number.

For this purpose, we consider a wave packet which is localized in energy space around a quantum number $`\overline{n}`$, that is, the coefficients $`\psi _n`$ assume only significant values in a small vicinity of $`n=\overline{n}`$. This feature allows us to expand the energy spectrum

$$E_nE_{\overline{n}}+\frac{E_n}{n}|_{n=\overline{n}}(n\overline{n})+\frac{1}{2}\frac{^2E_n}{n^2}|_{n=\overline{n}}(n\overline{n})^2$$
(3)
around $`\overline{n}`$. Here we have assumed that the spectrum is smooth and dense. Moreover, we have omitted higher order terms.

The approximation, Eq. (3), is valid for smooth potentials and a narrow energy distribution of the wave packet. To be specific, we consider the sodium dimer $`\text{Na}_2`$. Here a proper laser pulse lifts the vibrational ground state in the potential of the electronic state $`\text{X}^{\mathrm{\hspace{0.17em}1}}\mathrm{\Sigma }_\text{g}^+`$ to the excited electronic state $`\text{A}^{\mathrm{\hspace{0.17em}1}}\mathrm{\Sigma }_\text{u}^+`$. In the potential curve of the excited electronic state the wave packet is no longer the ground state, but has a small spectral width around a mean value $`E_{\overline{n}}`$. Indeed, in this case third order terms can be neglected bib:leichtle1 .

To study the wave packet dynamics, Eq. (2), we use the approximation, Eq. (3), and obtain the phase

$$E_nt/\mathrm{}E_{\overline{n}}t/\mathrm{}+\frac{2\pi t}{T_{\mathrm{cl}}}(n\overline{n})+\frac{2\pi t}{T}(n\overline{n})^2.$$
(4)
Here we have introduced the classical period $`T_{\mathrm{cl}}`$ and revival time $`T`$.

Thus the autocorrelation function, Eq. (1), reads

$$S(t)e^{iE_{\overline{n}}t/\mathrm{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}W_m\mathrm{exp}\left[2\pi i\left(\frac{m}{T_{\mathrm{cl}}}+\frac{m^2}{T}\right)t\right]$$
(5)
with a shifted summation index $`mn\overline{n}`$ and the weight function $`W_m|\psi _{m+\overline{n}}|^2`$. Moreover, we have extended the lower bound of the sum to $`\mathrm{}`$, since in this domain the weight factors $`W_m`$ do not contribute significantly.

### II.2 An allusion to factorization

In order to demonstrate the ability of the autocorrelation function to factorize a number, we start from a Gaussian weight function $`W_m\mathrm{exp}\left[m^2/2\mathrm{\Delta }m^2\right]`$ and fix the ratio $`NT/T_{\mathrm{cl}}`$ to an integer. In this case, we have to deal with the sum

$$S_N(\tau )\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}W_m\mathrm{exp}\left[2\pi i\left(m+\frac{m^2}{N}\right)\tau \right]$$
(6)
for dimensionless time $`\tau t/T_{\mathrm{cl}}`$.

The sum Eq. (6) is already the key to factorization, although it might not be obvious at this moment. In order to bring this latent capability to light, we focus in Fig. 1 on the modulus squared $`\left|S_N(\tau )\right|^2`$ of the autocorrelation function for $`N=55`$. We show this quantity for times $`\tau `$ in the vicinity of prime numbers $`\mathrm{}`$. Only if $`\mathrm{}`$ is a factor of $`N`$ the function displays a clear symmetric maximum.

The phenomenon of fractional revivals bib:leichtle1  lies at the heart of this feature. It is well-known that wave packets can show revivals at multiples of the time $`T`$ and fractional revivals at times $`t=\frac{p}{q}T`$ with $`p`$ and $`q`$ being mutually prime integers. In our case, this fact implies that at integer multiples of $`T_{\mathrm{cl}}`$ — or in other words for integers $`\tau `$ — we find a fractional revival of the wave packet, that is, a symmetric sequence of Gaussian shaped humps. A detailed analysis of this structure bib:varenna  brings to light that in most cases neighboring Gaussians overlap and interfere. For this reason, at most times no distinguished revivals appear. Only for times $`\tau `$ equal to a factor of $`N`$ the peaks are sharp enough so that the central peak does not overlap with its neighbors and reveals a clear maximum.

## III Quantum rotor

The crucial point in our factorization scheme is the interference of quadratic phase factors. In the previous section, they were the consequence of an anharmonic potential. But quadratic phase factors even appear for a free particle with a continuous energy spectrum. However, in order to encode integers we need a discrete one. Such a spectrum can be achieved, for example, by periodic boundary conditions. A diatomic molecule rotating perpendicular to its molecular axis is one example of such a quantum rotor.

### III.1 Time evolution

In order to lay the foundations for the factorization scheme presented in this section we first briefly review the essential properties of the one-dimensional quantum rotor. Then we derive the corresponding autocorrelation function.

The time evolution of the quantum rotor is governed by the Hamiltonian

$$\widehat{H}\frac{\widehat{J}^2}{2I}$$
(7)
with the operator $`\widehat{J}`$ of the angular momentum and the moment of inertia $`I`$. Periodic boundary conditions enforce a discrete spectrum of $`\widehat{J}`$ with eigenstates $`|k`$ and eigenvalues $`\mathrm{}k`$ with $`k`$.

The resulting discrete energy spectrum $`E_k=\mathrm{}^2k^2/(2I)`$ gives rise to the time evolution

$$|\psi (t)=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}\psi _k\mathrm{exp}\left[i\mathrm{}\frac{k^2}{2I}t\right]|k$$
(8)
of the quantum rotor.

We model an initially well-localized wave packet using the expansion coefficients

$$\psi _k\{\begin{array}{cc}1/\sqrt{N}\hfill & \text{ for }0k<N\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array},$$
(9)
which leads to the explicit expression

$$S(t)S_N(t)=\frac{1}{N}\underset{k=0}{\overset{N1}{}}\mathrm{exp}\left[2\pi ik^2\frac{t}{T}\right]$$
(10)
for the autocorrelation function. Here we have introduced the revival time $`T4\pi I/\mathrm{}`$.

Again quadratic phases appear. However, the linear term present in the expression Eq. (6) is absent.

### III.2 Factorization

In order to bring out most clearly the number $`N`$ to be factorized, we cast the sum, Eq. (10), in the form

$$\stackrel{~}{S}_N(n)=\frac{1}{N}\underset{k=0}{\overset{N1}{}}\mathrm{exp}\left[2\pi ik^2\frac{n}{N}\right]$$
(11)
with $`nN\frac{t}{T}`$. We emphasize that only for integers $`n`$ the autocorrelation function $`\stackrel{~}{S}_N`$ is periodic with period $`N`$.

In Fig. 2(a) and (d) we show $`|\stackrel{~}{S}_N(n)|^2`$ for $`n=1`$ to $`N`$, with $`N=21`$ and $`N=15`$, respectively. We recognize dominant maxima for values of $`n`$ whenever $`n`$ and $`N`$ share a common factor. The appearance of these maxima already indicates that a sum of only quadratic phases allows us to find prime factors of $`N`$.

We can even go a step further and concentrate on the real and imaginary part of $`\stackrel{~}{S}_N`$, depicted in Fig. 2(b) and (c) for $`N=21`$, and (e) and (f) for $`N=15`$, respectively. It is remarkable that for $`N=21`$ only numbers which share a common factor with $`N`$ appear in the imaginary part, Fig. 2(c), whereas all the other ones vanish. This is not the case for $`N=15`$, where only one factor and its multiples appear in the real part, Fig. 2(e).

### III.3 Gauß sums

We recognize that the autocorrelation function, Eq. (11), is of the form of a Gauß sum

$$G(a,b)\frac{1}{b}\underset{m=0}{\overset{b1}{}}\mathrm{exp}\left[2\pi im^2\frac{a}{b}\right].$$
(12)
Therefore, the structures found in Fig. 2 result from the intrinsic properties of Gauß sums. In order to highlight this connection, we first briefly review the properties of Gauß sums and then apply the results to the autocorrelation function Eq. (11).

For odd $`b`$ and mutually prime $`a`$ and $`b`$ the sum can be evaluated bib:lang  and yields

$$G(a,b)=\frac{1}{\sqrt{b}}\left(\frac{2a}{b}\right)\mathrm{exp}\left[i\frac{\pi }{4}(b1)\right].$$
(13)
Here we have introduced the Jacobi symbol $`\left(\frac{2a}{b}\right)`$ which expresses the number-theoretical complexity of the quadratic sum and which can take on the values $`\pm 1`$.

A closer look at Fig. 2 reveals that the emergence of factors in either the real or imaginary part depends strongly on the representation of the factors as $`4s+1`$ or $`4s+3`$ with $`s`$ being a non-negative integer. Indeed, all odd numbers $`r`$ are either of the form $`4s+1`$ or $`4s+3`$, and therefore, they are members either of the set

$$_1\left\{r\right|r=4s+1\}$$
(14)
or

$$_3\left\{r\right|r=4s+3\}.$$
(15)
For the case $`b_1`$ the sum $`G(a,b)`$ is purely real. Indeed, for $`b=4s+1`$ we find

$$G\mathrm{exp}\left[i\frac{\pi }{4}\mathrm{\hspace{0.17em}4}s\right]=\mathrm{exp}\left[i\pi s\right]=(1)^s.$$
(16)
When $`b_3`$, that is $`b=4s+3`$, we arrive at

$$G\mathrm{exp}\left[i\frac{\pi }{4}(4s+2)\right]=\mathrm{exp}\left[i\pi \left(s+\frac{1}{2}\right)\right]=(1)^si$$
(17)
which is purely imaginary.

Moreover, Eq. (13) implies that the modulus

$$|G(a,b)|=\frac{1}{\sqrt{b}}$$
(18)
is independent of the value of $`a`$, provided that $`a`$ and $`b`$ are mutually prime.

### III.4 An explanation of the patterns

With the help of these properties of the Gauß sum we can now relate the structures in the autocorrelation function

$$\stackrel{~}{S}_N(n)=G(n,N)$$
(19)
apparent in Fig. 2 to the factors $`p`$ and $`q`$ of $`N`$. Here we assume that $`p`$ and $`q`$ are the only factors of the number $`N`$, that is $`N=pq`$. Moreover, we stipulate that $`p`$ and $`q`$ are both odd. For the case of even $`N`$ one can easily extract factors $`2`$.

According to Eq. (18) the modulus squared $`|\stackrel{~}{S}_N|^2`$ is given by $`1/N`$ if $`n`$ and $`N`$ do not share a common factor. Otherwise, the common factor cancels in Eq. (12) and we have $`b=1/p`$ or $`b=1/q`$ and thus with the help of Eq. (18) we have $`|\stackrel{~}{S}_N|^2=1/p`$ or $`|\stackrel{~}{S}_N|^2=1/q`$, respectively.

We are now in the position to explain the patterns in the real and imaginary parts using the classification of the prime factors into the sets $`_1`$ and $`_3`$. We have to distinguish three cases: (i) $`p,q_1`$, (ii) $`p,q_3`$ and (iii) $`p_1,q_3`$ or vice versa. The cases (i) and (ii) imply $`N_1`$, whereas (iii) results in $`N_3`$.

We start our discussion with an example for case (ii). For $`N=21=45+1=37`$ both factors $`p=3=40+3`$ and $`q=7=41+3`$ are elements of the set $`_3`$ and $`N_1`$. If $`n`$ and $`N`$ do not share a common factor then $`b=N`$ and the value of the Gauß sum is real. Otherwise, we can extract a factor and get $`b=q`$ or $`b=p`$, respectively. Consequently, the value of the sum is purely imaginary, since $`p,q_3`$. Since the Jacobi symbol mentioned in Eq. (13) can only take on $`\pm 1`$ the corresponding value of $`a=n/p`$ or $`a=n/q`$ does not change this behavior.

The same arguments hold true for the example $`N=15=43+3_3`$ illustrating case (iii). Indeed, the factors $`p=3=40+3_3`$ and $`q=5=41+1_1`$ are not members of the same set. Therefore, if $`n`$ is a multiple of $`p`$ we have $`b=5`$ and thus the sum is purely real. In all other cases, we have either $`b=N`$ or $`b=p`$ which are both members of $`_3`$ and yield purely imaginary results.

The last case (i), not shown in Fig. 2, discusses factors $`p`$, $`q`$ and $`N_1`$, resulting in purely real values of the sum for all $`n`$. Unfortunately, in this case we do not have vanishing values for integers that are not factors of $`N`$.

## IV Conclusions and outlook

The time evolution of a wave packet brings to light key aspects of number theory. This feature allows us to find prime factors $`p`$ and $`q`$ of an appropriately encoded number $`N=pq`$. These properties are intrinsic to the autocorrelation function, which can be easily measured. At the heart of this phenomenon lie the mathematical properties of Gauß sums.

Wave packets of the quantum rotor exhibit these number-theoretical properties in the cleanest form. The most remarkable result is, that in the case of $`p,q_3`$ we find all factors (and its multiples) of $`N=pq`$ in the imaginary part of the autocorrelation function. All other integers yield a vanishing imaginary part.

But how to make use of these unique features? Imagine that $`\text{Im}(\stackrel{~}{S}_N)`$ represents the statistics of an observable $`\widehat{O}`$. Thus, each single measurement of $`\widehat{O}`$ would provide one factor of $`N`$ in a non-probabilistic way. Similarly, in the mixed case $`p_1`$ and $`q_3`$ the measurement of the real part of $`\stackrel{~}{S}_N`$ would provide one factor with high probability.

These insights pave the way to an effective method of factorizing numbers. The next step is to find a composite system which (i) intrinsically includes quadratic phase factors (ii) allows parallel calculation of all possible factors and (iii) enables us to find a suitable observable, which reveals all the information about the factors in a deterministic way by use of entanglement.

###### Acknowledgements.

We thank I.Sh. Averbukh, M.V. Berry, H. Maier, I. Marzoli and F. Straub, for many fruitful discussions. Moreover, one of us (W.P.S.) is most grateful to M.S. Brandt for organizing a most stimulating Heraeus seminar. The work of H.M., M.F. and W.P.S. was supported by the Deutsche Forschungsgemeinschaft, by the European Commission through the IST network QUBITS and by the Landesstiftung Baden-Württemberg in the framework of the Quantum Information Highway A8.


text:
# Compton Heated Outflow from CDAFs

## 1 Introduction

Advection-dominated accretion flows (ADAF) nicely complement the classic thin disk accretion flows (Shakura & Sunyaev 1973; Ichimaru 1977; Rees et al. 1982; Narayan & Yi 1994, 1995; Abramowicz et al. 1995), and have been successfully applied to variety of cosmic objects, from galactic X-ray binaries to diffuse X-ray background (see Narayan, Mahadevan, & Quataert 1999 for review). However, analytic studies of ADAFs indicated the convectively unstable nature of ADAFs (Narayan & Yi 1994, 1995a; see Begelman & Meier 1982 for radiation-dominated ADAF), which has subsequently been proved in a series of numerical simulations (Igumenshchev, Chen, & Abramowicz 1996; Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle, & Begelman 1999; Igumenshchev, Abramowicz, & Narayan 2000). Especially, the numerical studies of ADAFs by Igumenshchev & Abramowicz (1999, 2000) show that an ADAF becomes a convection-dominated accretion flow (CDAF) whenever the viscosity parameter $`\alpha 0.1`$. Further analyses of self-similar CDAF solutions show their unique properties (Narayan, Igumenshchev & Abramowicz 2000, hereafter NIA; Quataert & Gruzinov 2000, hereafter QG): the density varies as $`\rho R^{1/2}`$ ($`R`$ is the radius), much flatter than the usual $`R^{3/2}`$ in an ADAF or in spherical accretion. Correspondingly, the mean radial velocity varies as $`vR^{3/2}`$, compared to $`R^{1/2}`$ in ADAF or in spherical flow. Energy generated at the inner part of the flow is transported to the outer part by convection. A CDAF has perhaps as much resemblance to the rotating stellar envelope of a massive star as to the usual accretion flow.

The other aspect of the two-dimensional nature of ADAF or CDAF solutions, often neglected, is the interaction between the outgoing radiation produced at smaller radii with the inflowing gas in the outer part of the flow. In optically thick stars this interaction plays a vital role in establishing the equilibrium states. The radiative interaction also plays a very important role in pure spherical accretion flows (Ostriker et al. 1976; Cowie, Ostriker, & Stark 1978; Wandel, Yahil, & Milgrom 1984; Park 1990a, 1990b; Nobili, Turolla, & Zampieri 1991; Zampieri, Miller, & Turolla 1996; Ciotti & Ostriker 1997, 2001). Park & Ostriker (1999, 2001) studied the same radiative interaction in the context of the ADAF solution and found that a polar outflow can be generated through Compton heating of electrons by high-energy photons emitted by the inner, hot part of the flow. The winds generated by the processes in the papers listed above are not momentum driven. Rather, they are caused by overheating of the gas in the slowly moving, low density polar regions. Outflows may also be generated from ADAF by other hydrodynamic processes (Narayan & Yi 1995; Xu & Chen 1997; Blandford & Begelman 1999).

In this work, we study the conditions for CDAFs to develop radiation driven outflows. We adopt the self-similar CDAF solution as the background flow structure (NIA; QG). The treatment in this paper is two dimensional, adopting the angular profile of NIA and QG, except that the radiation field is simplified and treated as spherically symmetric.

## 2 Flow Properties

### 2.1 Density and Temperature

Multi-dimensional numerical simulations and analytic analyses show that density and temperature profiles of CDAF follow a self-similar form in radius (NIA; QG; Ball, Narayan & Quataert 2001, hereafter BNQ):

$`\rho (r)`$ $``$ $`r^{1/2}`$ (1)
$`T_i(r)`$ $``$ $`r^1,`$ (2)
where $`r`$ is the radius in units of the Schwarzschild radius, $`R_S2GM/c^2=3.0\times 10^5m\mathrm{cm}`$, and $`m`$ is the black hole mass in solar units. This shallow radial dependence corresponds to marginal stability to convection. QG further showed that the density should follow a power law in $`\mathrm{sin}\vartheta `$ ($`\vartheta `$ is the angle from the pole) while the temperature is constant on spherical shells. So we shall assume that the density is given by

$`\rho (r,\vartheta )`$ $`=`$ $`\rho _0r^{1/2}\mathrm{sin}^2\vartheta `$ (3)
$`=`$ $`\rho _{out}(r/r_{out})^{1/2}\mathrm{sin}^2\vartheta .`$ (4)
The flow extends from $`r_{out}`$ down to $`r_{in}=1`$, with the outer radius (analogous to the stellar photosphere) fixed by energy balance considerations.

The electron temperature $`T_e`$ in an accretion flow is determined by the balance between various cooling and heating processes. In general, detailed studies of hot accretion flows (Narayan & Yi 1995b; Narayan, Barret, & McClintock 1997; BNQ) show that the electron temperature is equal to the ion temperature for $`T_e10^9\mathrm{K}`$, and then, at smaller radii, due to highly efficient relativistic bremsstrahlung and synchrotron emission, it flattens to somewhere between $`10^9\mathrm{K}`$ to $`10^{11}\mathrm{K}`$, depending on the amount of direct viscous heating to electrons. Hence, in this paper, we approximate the electron temperature as

$`T_e(r,\vartheta )`$ $`=`$ $`{\displaystyle \frac{T_0}{r}}\text{for}r>r_1`$ (5)
$`=`$ $`T_1\text{for}rr_1,`$ (6)
where $`T_010^{12}\mathrm{K}`$, $`T_1`$ is some constant between $`10^9\mathrm{K}`$ and $`10^{11}\mathrm{K}`$, and $`r_1T_0/T_1`$.

From equation (3), we find that the total Thomson optical depth along a given direction $`\vartheta `$ is

$$\tau _{es}\frac{1}{2}\dot{m}r_{out}^{1/2}\mathrm{sin}^2\vartheta ,$$
(7)
where $`\dot{m}\dot{M}/\dot{M_E}=\dot{M}c^2/L_E`$ is the dimensionless mass accretion rate. Our analysis, which is to be applied to low accretion rate solutions, will assume that the flow is optically thin to both scattering and redistribution processes.

### 2.2 Radiation Field and Radiation Temperature

The inner part of the accretion flow generally produces higher energy photons, while the outer part produces lower energy photons. Hence the radiation field at a given position is determined by the contribution from inner radiating shells plus the outer shells. Although in each shell, the density, and therefore, the emissivity (i.e. cooling function) is a function of poloidal angle $`\vartheta `$, the energy density is more uniform over angle than is the emissivity. So, we will simplify the calculation of radiative transfer by assuming each radiating shell is homogeneous, using an angle average over the sphere.

Therefore, the radiation energy density at given $`r`$, $`E_X(r)`$, is given by

$$cE_X(r)=R_s\left[\frac{1}{4\pi r^2}_{r_{in}}^r\mathrm{\Lambda }(r^{})4\pi r^2𝑑r^{}+_r^{r_{out}}\mathrm{\Lambda }(r^{})\frac{r^{}}{r}\mathrm{ln}\sqrt{\frac{r^{}+r}{r^{}r}}dr^{}\right],$$
(8)
where the first term is the contribution from inner radiating shells and the second term from outer radiating shells (see Appendix A). The emissivity per unit volume is denoted as the cooling function $`\mathrm{\Lambda }`$, which will be described in §3.1.

For relativistic un-Comptonized bremsstrahlung, with $`\mathrm{\Lambda }\rho ^2T_e`$ and CDAF scaling $`\rho r^{1/2}`$ and $`T_er^1`$, the incremental luminosity $`dL_X(r)/dr=\text{constant}`$, and the luminosity at a given radius $`r`$ has the largest contribution from the outermost shell. Nevertheless, equation (8) gives a peak energy density in the inner region (cf. dashed line in Fig. 1) and a significant radiation pressure gradient. In addition, when there is strong Comptonization, the contribution from the inner region can be significantly enhanced.

Moreover, the amount of Compton heating is determined by the luminosity times the photon energy, and is largely dominated by the inner hot region. The radiation temperature, defined as $`kT_X(r)`$ being the energy-weighted mean photon energy, at a given radius $`r`$ is determined by the equation

$$cT_X(r)E_X(r)=R_S\left[\frac{1}{4\pi r^2}_{r_{in}}^rT_X(r^{})\mathrm{\Lambda }(r^{})4\pi r^2𝑑r^{}+_r^{r_{out}}T_X(r^{})\mathrm{\Lambda }(r^{})\frac{r^{}}{r}\mathrm{ln}\sqrt{\frac{r^{}+r}{r^{}r}}dr^{}\right].$$
(9)
This quantity is more concentrated than $`E_X(r)`$, especially in the presence of strong Comptonization \[see dotted curve in Fig. 1 for $`T_X(r)`$\].

In presence of an overheated wind, some region along the pole in CDAF will not be able to accrete matter, while the remaining equatorial region normally accretes (see §3.3). In most of the relevant parameter space, the overheated ‘funnel’ is very narrow, and we assume (with regard to energy generation) that the self-similar CDAF flow is filling the whole space including the polar region, omitting the small correction due to the solid angle $`\mathrm{\Omega }_W`$ ($`0.2`$ sr) occupied by the outflowing wind/jet.

The resulting profiles of radiation moments and the radiation temperature show the aforementioned characteristics of CDAFs. The luminosity,

$$L(r)=R_S^3_{r_{in}}^r\mathrm{\Lambda }(r^{})4\pi r^2𝑑r^{},$$
(10)
increases at large radius (solid curve in Fig. 1), while the radiation temperature is roughly flat or decreasing slowly outward after peaking at some intermediate radius (dotted curve in Fig. 1).

### 2.3 Outer Boundary

The mathematical self-similar CDAF has zero mass accretion rate (NIA), and this, of course, corresponds to zero luminosity. However, the accretion flow must lose energy by radiative cooling, and this energy loss is very likely provided by the convective energy transport (NIA, BNQ). We follow BNQ to assume that the total luminosity of ADAF is proportional to the mass accretion rate:

$$L(r_{out})=\eta _cϵ_c\dot{M}c^2,$$
(11)
where the convective efficiency $`ϵ_c`$, determined by numerical simulations to be in the range $`10^210^3`$ and $`\eta _c`$ is the fraction of convected energy to be radiated away by electrons. In subsequent examples, we nominally adopt $`ϵ_c\eta _c=10^2`$ (BNQ). We also tried smaller values of $`ϵ_c\eta _c`$ or a scaling of $`ϵ_c\eta _c`$ with $`\dot{m}`$ (cf. §3.3), and the results were qualitatively the same.

We determine the outer boundary for a given mass accretion rate $`\dot{m}`$ by finding the radius $`r_{out}`$ at which equation (11) is satisfied. This gives a somewhat different relation between the total luminosity and the mass accretion rate from that in BNQ because we have used Comptonized relativistic bremsstrahlung (and synchrotron) while BNQ employ only un-Comptonized non-relativistic bremsstrahlung. However, the general characteristic of the relation is the same: higher mass accretion flows have smaller outer boundary radii. This is due to the special radial profile of density in CDAF solutions. Because of the weak density gradient with radius, most of the luminosity contribution is from the outer part of the flow, which makes CDAF flows quite distinct from other accretion flows including ADAFs. The dimensionless total luminosity, $`L(r_{out})/L_E`$, is directly proportional to the dimensionless mass accretion rate $`\dot{m}\dot{M}/\dot{M}_E`$ by equation (11) while the luminosity accumulated from bremsstrahlung is very roughly proportional to $`\dot{m}^2r_{out}`$ (with a higher power in $`\dot{m}`$ if Comptonization is important). Since these two luminosities must be equal, a higher $`\dot{m}`$ flow implies a smaller $`r_{out}`$ (see solid lines in Fig. 2). Therefore lower mass accretion rate CDAF solutions have lower luminosity, yet are more extended, if we treat $`ϵ_c\eta _c`$ as fixed, independent of $`\dot{m}`$.

## 3 Outflow

### 3.1 Cooling and Heating

The main cooling mechanism we consider in this work is Comptonized bremsstrahlung and Comptonized synchrotron. We first focus on the Comptonized bremsstrahlung. This process is not only the source of electron cooling but also the source of the Compton heating radiation field.

The Comptonized bremsstrahlung cooling rate per unit volume can be reasonably approximated by (Svensson 1982; Stepney & Guilbert 1983)

$$\mathrm{\Lambda }_{br}=\sigma _Tc\alpha _fm_ec^2n^2[F_{ei}(T_e)+F_{ee}(T_e)],$$
(12)
where

$`F_{ei}`$ $`=`$ $`4\left({\displaystyle \frac{2}{\pi ^3}}\right)^{1/2}\theta _e^{1/2}(1+1.781\theta _e^{1.34})\text{for}\theta _e<1`$
$`=`$ $`{\displaystyle \frac{9}{2\pi }}\theta _e\left[\mathrm{ln}(1.123\theta _e+0.48)+1.5\right]\text{for}\theta _e>1`$
$`F_{ee}`$ $`=`$ $`{\displaystyle \frac{5}{6\pi ^{3/2}}}(443\pi ^2)\theta _e^{3/2}(1+1.1\theta _e+\theta _e^21.25\theta _e^{5/2})\text{for}\theta _e<1`$
$`=`$ $`{\displaystyle \frac{9}{\pi }}\theta _e\left[\mathrm{ln}(1.123\theta _e)+1.2746\right]\text{for}\theta _e>1,`$
$`\theta _ekT_e/m_ec^2`$, $`n`$ is the electron (ion) number density, $`\sigma _T`$ the Thomson cross section, $`\alpha _f`$ the fine structure constant, and $`m_e`$ the electron mass. In this work, we assume pure hydrogen gas.

The emitted bremsstrahlung photons are upscattered by inverse Compton scattering. The amplification of photon energy by single scattering is $`A1+4\theta _e+16\theta _e^2`$. The probability of single scattering is $`P=1\mathrm{exp}(\tau _{es})`$ (see e.g., Dermer et al. 1991). Therefore the mean amplification factor by Compton scattering would be

$$\eta _0=1P+PA=1+P(A1).$$
(15)
The fraction $`1P`$ of photons remain unscattered while the fraction $`P`$ of photons are upscattered to $`A`$ times their initial energy. This prescription is valid only for a single scattering Comptonization. Since most CDAF flow considered here has $`\tau _{es}<1`$, we use this amplification factor for Comptonized bremsstrahlung. Although Dermer et al. (1991) provides a handy formula for $`\eta `$, applicable to diverse regimes, it can give incorrect values for $`A1`$ and $`\tau _{es}1`$, the main parameter regime of CDAF. Although $`\tau _{es}`$ is a hard-to-define quantity in a complex flow, we simply adopt $`\tau _{es}=n_e\sigma _Tr`$ in this work for the purpose of estimating the Comptonization.

When $`PA1`$, Comptonization becomes saturated, and all photons that are not absorbed are upscattered to $`3T_e`$, that is $`\eta \eta _{sat}=3\theta _e/x`$ where $`xh\nu /m_ec^2`$ (Dermer et al. 1991). Then the fully saturated Comptonized bremsstrahlung emission is

$`\mathrm{\Lambda }_{Cbr}^{sat}`$ $`=`$ $`{\displaystyle _{x_{abs}}^{3\theta _e}}\eta _{sat}ϵ_{br}(x)𝑑x+{\displaystyle _{3\theta _e}^{\mathrm{}}}ϵ_{br}(x)𝑑x`$ (16)
$`=`$ $`\mathrm{\Lambda }_{br}\left[{\displaystyle _3^{\mathrm{}}}e^t𝑑t+{\displaystyle _{x_{abs}/\theta _e}^3}3t^1e^t𝑑t\right]`$ (17)
where $`ϵ_{br}(x)`$ is the bremsstrahlung spectrum

$$ϵ_{br}(x)=\mathrm{\Lambda }_{br}\mathrm{exp}\left(\frac{x}{\theta _e}\right)dx$$
(18)
and $`x_{abs}`$ is the absorption frequency $`\nu _{abs}`$ in unit of $`m_ec^2/h`$. Since $`\mathrm{\Lambda }_{br}`$ is smaller by the factor $`\mathrm{exp}(x_{abs}/\theta _e)`$ in presence of absorption, we take the maximum amplification factor for bremsstrahlung as

$$\eta _{br}^{sat}=e^{x_{abs}/\theta _e}\times \left\{e^3+3\left[E_1\left(\frac{x_{abs}}{\theta _e}\right)E_1(3)\right]\right\},$$
(19)
where $`E_n_1^{\mathrm{}}t^n\mathrm{exp}(xt)𝑑t`$ is the exponential integral. The final Comptonized bremsstrahlung emission is then

$$\mathrm{\Lambda }_{Cbr}=\mathrm{min}(\eta _0,\eta _{br}^{sat})\mathrm{\Lambda }_{br}.$$
(20)
The energy-weighted mean photon energy for unsaturated Comptonized bremsstrahlung is (see e.g. PO2)

$$\frac{4kT_X}{m_ec^2}=\frac{_0^{\mathrm{}}\eta _{br}^2xϵ_{br}(x)𝑑x}{_0^{\mathrm{}}\eta _{br}ϵ_{br}(x)𝑑x}=\eta _{br}\theta _e,$$
(21)
and that for saturated Comptonized bremsstrahlung is simply

$$T_X=T_e,$$
(22)
since the spectrum approaches Wien spectrum which has $`T_X=T_e`$. Therefore, the radiation temperature of locally radiated bremsstrahlung emission is

$$T_X^{Cbr}=\frac{1}{4}T_e\mathrm{min}(\eta _{br},4).$$
(23)
The frequency $`\nu _{abs}`$ is chosen to be the frequency at which the free-free absorption optical depth is equal to 1,

$$ra_{ff}(\nu _{abs})=1,$$
(24)
where $`a_{ff}`$ is the absorption coefficient given by Dermer et al. (1991)

$$a_{ff}(x)=\sqrt{8\pi }\frac{\alpha _f^2\sigma _Tr_e^3}{x^2\theta _e^{3/2}\left[1+(8/\pi )^{1/2}\theta _e^{3/2}\right]}n_i^2\overline{g}$$
(25)
with

$`\overline{g}`$ $`=`$ $`(1+2\theta _e+2\theta _e^2)\mathrm{ln}\left[{\displaystyle \frac{4\eta _E(1+3.42\theta _e)\theta _e}{x}}\right]`$
$`+({\displaystyle \frac{3\sqrt{2}}{5}}+2\theta _e)\theta _e\mathrm{ln}\left[{\displaystyle \frac{4\eta _E(11.2+10.4\theta _e^2)\theta _e}{x}}\right],`$
$`\alpha _f`$ is the fine structure constant, and $`r_e=e^2/m_ec^2`$ the electron radius (Svensson 1984). This expression is valid for $`x\theta _e`$ and we used $`(\pi /2)^{1/2}\theta _e^{1/2}[1+(8/\pi )^{1/2}\theta _e^{3/2}]`$ to approximate $`\mathrm{exp}(1/\theta _e)\text{K}_2(1/\theta _e)`$.

Flow at a given position is heated (or cooled) by the inverse Comptonization off electrons. We use the heating rate (Levich & Sunyaev 1971)

$$\mathrm{\Gamma }_C=4\sigma _Tc[\theta _X(r)\theta _e(r)]E_X(r)n_e(r,\vartheta ),$$
(27)
where $`E_X(r)`$ is the radiation energy density from equation (8) and $`\theta _XkT_X/m_ec^2`$ the radiation temperature from equation (9).

### 3.2 Equilibrium Temperature and Overheating

In original CDAF, the temperature of the gas is determined by the balance between the viscous heating plus the convective energy transport versus the radiative cooling. However, if the radiative heating is dominant in some region of the flow (the condition under which this assumption is valid is discussed in §3.4.), the temperature in that region will change to reach a new equilibrium. The new equilibrium temperature of the flow will be determined by the balance between radiative heating and radiative cooling.

In CDAF considered in this work, Compton heating versus Comptonized bremsstrahlung cooling are the main radiative processes. The thermal equilibrium temperature $`T_{eq}`$, then, satisfies

$$\mathrm{\Gamma }_C(T_{eq})=\mathrm{\Lambda }_{Cbr}(T_{eq})$$
(28)
at given position $`(r,\vartheta )`$. From equations (20) and (27), $`\theta _{eq}kT_{eq}/m_ec^2`$ is determined by

$$4cE_X(r)[\theta _X(r)\theta _{eq}]=\alpha _fm_ec^3n_e(r,\vartheta )\eta _{Cbr}(n_e,r,\theta _{eq})[F_{ei}(\theta _{eq})+F_{ee}(\theta _{eq})].$$
(29)
Since the number density $`n_e(r,\vartheta )`$ decrease as $`\vartheta 0`$ (toward the pole), the derived electron temperature increases toward the pole as long as $`\theta _X>\theta _e`$.

For smaller enough $`\vartheta `$, the thermal equilibrium temperature can be higher than the virial temperature,

$$T_{eq}(r,\vartheta )>T_{vir}(r).$$
(30)
The virial temperature is defined as $`(5/2)kT_{vir}=m_pGM/r`$, where $`m_p`$ is the proton mass. Once electrons are heated above the virial temperature, and therefore above the ion temperature, electrons are likely (via collisions and instabilities) to heat ions to above the virial temperature, thereby, creating winds, especially because the dynamical time of CDAF is much longer than the free-fall flow. Therefore, we adopt equation (30) as the condition for overheating and producing a wind.

What equation (30) means is that some region of CDAF will be radiatively heated from the background temperature $`T_e`$ to $`T_{eq}`$ to achieve the thermal equilibiurm. However, as the temperature of the flow approaches $`T_{eq}`$, it will become unbound when the condition (30) is met, and thereby producing the outflow.

There is also a trivial, yet additional constraint for overheating: the radiation temperature must be higher than the virial temperature

$$T_X(r)>T_{vir}(r),$$
(31)
otherwise the radiation field will cool the flow.

### 3.3 Outflow

Outflows will extend from the polar axis to the angle $`\vartheta _c`$ at which the equilibrium temperature is equal to the virial temperature,

$$T_{eq}(r,\vartheta _c)=T_{vir}(r).$$
(32)
Regions with $`\vartheta <\vartheta _c`$ are overheated to above the virial temperature. Hence, the shape of the outflow is determined by the angle $`\vartheta _c(r)`$ as a function of $`r`$.

One example of the outflow is shown in Figure 3 for the mass accretion rate of $`\dot{m}=10^3`$, the electron temperature of inner region $`T_1=10^{11}\mathrm{K}`$. The solid curve shows $`\vartheta _c(r)`$, and the outer dotted circle is the outer boundary of the flow. The outflow starts at $`r40`$ and the opening angle $`\vartheta _c`$ reaches maximum of $`17^{}`$ around $`r700`$ and ends with $`8.3^{}`$ at the outer boundary. Figure 4 shows the same for $`\dot{m}=2.0\times 10^3`$ and $`T_1=3.0\times 10^{10}\mathrm{K}`$ with $`\vartheta _c=6.6^{}`$ at the boundary. The inner isothermal region and the immediate surroundings have $`T_X<T_{vir}`$, and the overheating does not occur. However, at larger radius Compton heating is strong enough to overheat the polar region of the gas and produce outflow. The fraction of the sphere covered by the outflowing gas is $`2\pi \vartheta ^2/4\pi 0.01`$ for $`\vartheta =8.3^{}`$.

We have considered combinations of $`\dot{m}`$ and $`T_1`$. Outflow solutions exist within a limited range of these parameters. Too high $`\dot{m}`$ makes the outer boundary so small that the whole region is isothermal, i.e., $`T_X=T_e`$, and the flow is not heated. The outer boundary is smaller for higher $`\dot{m}`$ because radiative luminosity is proportional to $`\dot{m}^2`$ while CDAF flows have a fixed (assumed) radiation efficiency, which makes the total luminosity roughly proportional to $`\dot{m}`$. The only way to reconcile this is $`r_{out}`$ being smaller for higher $`\dot{m}`$ (BNQ). In a too low $`\dot{m}`$ flow, Comptonization is not strong enough to keep the radiation temperature high. Lower energy photons from larger radii dilute the radiation field and flows are not heated. In Figure 5, we show region of space in ($`\dot{m}`$,$`T_1`$) where outflow solutions are successfully produced. Circles denote ($`\dot{m}`$,$`T_1`$) for which outflows are found. The size of circle represents the opening angle of the outflow at the outer boundary, the majority being $`10^{}`$. Crosses denote ($`\dot{m}`$,$`T_1`$) values for which outflows do not exist.

As expected, flows with high electron temperature are more likely to develop outflows. Outflows are not expected for $`T_1<m_ec^2/k6\times 10^9\mathrm{K}`$: Comptonization is less efficient when $`kT_e/m_ec^2<1`$. Since the adopted radiation efficiency is $`ϵ_c\eta _c10^2`$, we expect CDAFs with outflows to exist in the luminosity range $`8\times 10^7L/L_E4\times 10^5`$ for $`T_1=10^{11}\mathrm{K}`$ and $`4\times 10^5L/L_E10^4`$ for $`T_1=10^{10}\mathrm{K}`$. These results do not depend on the mass of the black holes.

So far, we have assumed that the total radiation efficiency $`ϵ_c\eta _c`$ is constant. In fact, current theoretical understanding of the CDAF solution is not yet secure on this point. We now relax that assumption and study an alternative case when $`ϵ_c\eta _c`$ is assumed to be proportional to $`\dot{m}`$. We adopt simply $`ϵ_c\eta _c=10^2(\dot{m}/10^2)`$ so as to make $`ϵ_c\eta _c=10^2`$ for $`\dot{m}=10^2`$. Now for small $`\dot{m}`$ the total luminosity determined from equation (11) and the total bremsstrahlung emission from equation (10) are both proportional to $`\dot{m}^2`$, and the outer boundary is now at a roughly constant radius (see dotted lines in Fig. 2). At high $`\dot{m}`$, Comptonization adds additional power of $`\dot{m}`$ to the total emission, and $`r_{out}`$ depends on $`\dot{m}`$ (see dotted lines in Fig. 2). The outflow solutions for this choice of radiation efficiency are shown in Figure 6. Qualitatively, they are not much different from the solutions with constant $`ϵ_c\eta _c`$. Nonetheless, the flow is much less extended due to lower radiation efficiency for lower $`\dot{m}`$ and the shape of the outflow funnel in this case is almost straight. The whole flow structure is quite self-similar except for the magnitude of the opening angle. Most of the solutions still have the opening angle at the outer boundary $`20^{}`$.

### 3.4 Convection, Heating, and Cooling Timescales

In the original self-similar CDAF, the viscous heating is balanced by the convective energy transport. However, radiative heating can be dominant in some region of the flow. The importance of each heating and cooling process is measured by the corresponding timescale. The radiative heating timescale is given by the ratio between the internal energy of the electron, $`\epsilon _e`$, and the Compton heating rate,

$$t_H\frac{\epsilon _e}{\mathrm{\Gamma }_C}.$$
(33)
The radiative cooling timescale is similarly given by

$$t_C\frac{\epsilon _e}{\mathrm{\Lambda }_{Cbr}}.$$
(34)
The timescale for convective energy transport at a certain two-dimensional position is not trivial to evaluate because the two-dimensional velocity profile including the convective motion is needed. Although two-dimensional numerical simulations (Stone, Pringle, & Begelman 1999; Igumenshchev, Abramowicz, & Narayan 2000; Igumenshchev & Abramowicz 2000) suggest very limited flow motion near the pole as in self-similar two-dimensional ADAF (Narayan & Yi 1995), the two-dimensional flow motion has yet to be expressed in simple analytic form. Here, we adopt the height-averaged convective energy flux given by NIA,

$$F_c=\alpha _c\frac{c_s^2}{\mathrm{\Omega }_K}\rho T\frac{ds}{dr},$$
(35)
where $`\alpha _c`$ is the convection coefficient analogous to the usual Shakura & Sunyaev $`\alpha `$, $`c_s`$ the isothermal sound speed, $`\mathrm{\Omega }_K`$ the Keplerian angular velocity, and $`s`$ the entropy of the flow. In a self-similar solution that is marginally stable to convection, $`\alpha _c=3\alpha `$ is expected when $`\alpha 0.05`$ (NIA). The cooling rate per volume due to this convective energy flux is then

$$Q_{cv}\frac{F_c}{r}.$$
(36)
Corresponding timescale $`t_{cv}`$ is now

$$t_{cv}=\frac{\epsilon _i}{Q_{cv}}.$$
(37)
Since the actual convective motion near the pole is generally much smaller than the average convective motion near the disk midplane as is assumed in equation (35), the actual timescale for the convective cooling could be much longer than $`t_{cv}`$ evaluated by equation (37).

For given dimensionless mass accretion rate $`\dot{m}`$ and the inner electron temperature $`T_1`$, small $`\alpha `$ (of Shakura-Sunyaev), or equally small $`\alpha _c`$, means slower radial motion, and therefore, slower convective motion. So we searched for the critical value of $`\alpha `$ below which $`t_{cv}>t_H`$ is satisfied within the overheated region. For higher electron temperature of $`T_1=10^{11}\mathrm{K}`$. the value of the critical $`\alpha `$ ranges from $`2.0\times 10^3`$ for $`\dot{m}=2\times 10^3`$ to $`1.2\times 10^5`$ for $`\dot{m}=8\times 10^5`$, all at a constant radiation efficiency of $`10^2`$. For lower electron temperature of $`T_1=3\times 10^{10}\mathrm{K}`$. the value of the critical $`\alpha `$ ranges from $`2.7\times 10^2`$ for $`\dot{m}=8\times 10^3`$ to $`3.0\times 10^4`$ for $`\dot{m}=8\times 10^4`$. For even lower electron temperature of $`T_1=10^{10}\mathrm{K}`$ at which overheated region is very limited in radius, the value of the critical $`\alpha `$ ranges from $`1.6\times 10^2`$ for $`\dot{m}=10^2`$ to $`2.8\times 10^3`$ for $`\dot{m}=4\times 10^3`$. Expected values of $`\alpha `$ for which CDAF exists are in the range $`\alpha <0.05`$ (NIA). We, therefore, conclude that the Compton heated outflow is expected to be produced in CDAF if the viscosity parameter $`\alpha `$ is small enough, exact value being dependent on the mass accretion rate. However, as already mentioned, the real convective motion near the pole is expected to be much smaller than the average convective motion in the disk, the Compton heated outflow near the pole may develop even for much larger value of $`\alpha `$.

Another timescales we want to check are $`t_H`$ versus $`t_C`$ at the temperature of the background flow $`T_e(r,\vartheta )`$ given by equation (5). Only when radiative heating is greater than cooling, i.e., $`t_H|_{T_e}<t_C|_{T_e}`$, the flow is heated from $`T_e`$ to overheated temperature $`T_{eq}`$ of equation (28). The ratio $`(t_H/t_C)_{T_e}`$ can be expressed as

$$\frac{t_C}{t_H}=\frac{\mathrm{\Gamma }_C(T_e)}{\mathrm{\Lambda }_{Cbr}(T_e)}=\frac{\mathrm{\Gamma }_C(T_{eq})}{\mathrm{\Lambda }_{Cbr}(T_{eq})}\frac{\mathrm{\Gamma }_C(T_e)}{\mathrm{\Gamma }_C(T_{eq})}\frac{\mathrm{\Lambda }_{Cbr}(T_{eq})}{\mathrm{\Lambda }_{Cbr}(T_e)}=\frac{\mathrm{\Gamma }_C(T_e)}{\mathrm{\Gamma }_C(T_{eq})}\frac{\mathrm{\Lambda }_{Cbr}(T_{eq})}{\mathrm{\Lambda }_{Cbr}(T_e)},$$
(38)
which is always less than 1 as long as $`T_{eq}>T_e`$ because $`\mathrm{\Gamma }_C(T_e)`$ is a monotonically decreasing function of $`T_e`$ and $`\mathrm{\Lambda }_{Cbr}(T_e)`$ a monotonically increasing function of $`T_e`$. The ratio $`\mathrm{\Gamma }_C(T_{eq})/\mathrm{\Lambda }_{Cbr}(T_{eq})=1`$ by the definition of $`T_{eq}`$ in equation (28).

### 3.5 Synchrotron Emission

The synchrotron emission of hot electrons can produce copious soft photons some of which are subsequently Compton upscattered. This will lead to an increase in the total cooling rate as well as a decrease in the radiation temperature of the emitted radiation. Increased cooling rate will reduce the outer boundary of the CDAF whereas decreased radiation temperature will reduce the Compton heating, or even cool the flow if the radiation temperature falls below the electron temperature.

So we recalculated the condition for overheating in presence of magnetic field for the flow parameters considered in Fig. 5. The treatment of Comptonized synchrotron is described in Appendix B. We fix the black hole mass to be $`10^8M_{}`$. We find that the synchrotron emission can significantly affect the high temperature flow whereas the low temperature flow is less affected. The overheating and subsequent production of the outflow is still possible for flow with $`T_1=10^{10}\mathrm{K}`$ as long as the magnetic field is less than 10% of the equipartition (the gas pressure being equal to the magnetic pressure) field. The outer boundary of the flow and the opening angle of the outflow are hardly changed if this condition is met. For $`T_1=3\times 10^{10}\mathrm{K}`$, we find that the magnetic field has to be smaller than 3% of the equipartition field to produce the outflow. For the highest temperature $`T_1=10^{11}\mathrm{K}`$ flow, the production of the outflow is possible only when the magnetic field is less than 1% of the equipartition field. So we conclude that the Compton preheated outflow is possible in the magnetic CDAFs as long as the magnetic field is less than from 10% to 1% of the equipartition field.

## 4 Summary and Discussion

Hot accretion flows like ADAFs have a number of physical characteristics that compliment the classic low-temperature disk flows. In previous work, we have explored the consequences of the high temperature and the two-dimensional density structure of these flows, and have found that ADAFs may be able to produce radiatively driven outflows. We subsequently have noted that self-similar CDAFs have even more suitable properties for producing outflows as compared to ADAFs: steeper poloidal density gradients and higher radiation efficiencies. In this paper, we have studied the conditions for self-similar two-dimensional CDAFs (NIA; QG) to develop radiatively heated polar outflows.

1. We have found that CDAFs produce enough luminosity and photon energy to drive polar outflows via Compton heating for a reasonable range of mass accretion rate, or, equally, luminosity, as long as the magnetic field is less than from 10% to 1% of the equipartition field. When the electron temperature saturates around $`10^{11}\mathrm{K}`$ at the inner region, polar outflows are possible for $`8\times 10^7L/L_E4\times 10^5`$ for radiation efficiency of $`10^2`$, where $`L_E`$ is the Eddington luminosity. The luminosity range for which outflow exists is narrower for lower electron temperature flows and disappears completely for electron temperature $`6\times 10^9\mathrm{K}`$

2. In most cases, outflows are well collimated along the rotation axis, with an opening angle typically in the range $`10^{}`$.

3. If we, instead of taking efficiency as constant, assume that it is proportional to the mass accretion rate $`\dot{m}`$, the solutions are qualitatively the same but are more self similar, i.e., opening angle and outer boundary (in Schwarzschild units) depend less strongly on $`\dot{m}`$.

4. Outflow is more probable for small viscosity parameter $`\alpha `$.

The treatment in this work is not completely satisfactory in the sense that the dynamics and the temperature profiles of CDAFs are not self-consistently solved. However, it was not our intention to solve fully three-dimensional, self-consistent gloabl CDAFs with proper consideration for all gas, radiative, and magnetic processes, which will be eventually needed to fully understand CDAFs. Rather, our goal has been to show that even in the framework of simple self-similar solutions, radiatively heated polar outflows appear as natural consequences of the physical characteristics of CDAFs.

This work is the result of research activities (Astrophysical Research Center for the Structure and Evolution of the Cosmos) supported by Korea Science & Engineering Foundation.

## Appendix A Radiation Field Inside An Optically Thin Spherical Shell

Let’s consider an optically thin, uniform radiating shell with radius $`r^{}`$ and the thickness $`\mathrm{\Delta }r^{}`$. The specific intensity at radius $`r`$ ($`r<r^{}`$) in the direction of $`\vartheta `$ (see Figure 7) is given by integration of radiative transfer equation along that direction. If we denote the emissivity (per unit volume per unit solid angle) of the shell as $`ϵ_\nu /4\pi `$, then the specific intensity is simply the emissivity times the path length,

$$I_\nu =\frac{ϵ_\nu }{4\pi }\mathrm{\Delta }s=\frac{ϵ_\nu }{4\pi }\frac{\mathrm{\Delta }r^{}}{\mathrm{cos}\varphi },$$
(A1)
where $`\varphi `$ is the angle between the ray and the normal of the shell. From the law of sines, $`\mathrm{sin}\varphi =(r/R)\mathrm{sin}\vartheta `$. The radiation energy density, $`E_\nu (r)`$, is given by the integral of $`I_\nu `$ over the solid angle $`d\mathrm{\Omega }`$,

$`E_\nu (r)`$ $`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \frac{ϵ_\nu }{4\pi }\frac{\mathrm{\Delta }r^{}}{\mathrm{cos}\varphi }𝑑\mathrm{\Omega }}`$ (A2)
$`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \frac{ϵ_\nu }{4\pi }}\mathrm{\Delta }r^{}{\displaystyle \frac{\mathrm{sin}\vartheta d\vartheta }{\sqrt{1\frac{r^2}{r^2}\mathrm{sin}^2\vartheta }}𝑑\phi }`$ (A3)
$`=`$ $`{\displaystyle \frac{1}{c}}ϵ_\nu \mathrm{\Delta }r^{}{\displaystyle \frac{r^{}}{r}}\mathrm{ln}\sqrt{{\displaystyle \frac{r^{}+r}{r^{}r}}}.`$ (A4)
Equation (A2) formally diverges when $`rr^{}`$ whereas the correct value saturates. In the limit when $`r^{}r\mathrm{\Delta }r^{}`$, the effect of shell’s curvature has to be incorporated in $`\mathrm{\Delta }s`$. However, since the radiation energy density is calculated by the sum of contributions from discrete shells in $`r`$ (eq. ), the divergence is automatically avoided.

## Appendix B Comptonized Synchrotron Emission

The angle-averaged synchrotron emission by relativistic Maxwellian electrons is given by (Pacholczyk 1970)

$$ϵ_{syn}(\nu )d\nu =\frac{2\pi }{\sqrt{3}}\frac{e^2}{c}\frac{n_e\nu }{\theta _e^2}I^{}(x_M)d\nu $$
(B1)
where $`x_M2\nu /(3\nu _0\theta _e^2)`$, $`\nu _0eB/(2\pi m_ec)`$, and

$$I^{}(x_M)=\frac{4.0505}{x_M^{1/6}}\left(1+\frac{0.40}{x_M^{1/4}}+\frac{0.5316}{x_M^{1/2}}\right)\mathrm{exp}(1.8899x_M^{1/3})$$
(B2)
is a fitting formula (Mahadevan, Narayan, & Yi 1996). When absorption is not important, the cooling rate due to the optically thin synchrotron emission is obtained by integrating the equation (B1)

$$\mathrm{\Lambda }_{syn}^0=213.6\frac{e^2}{c}n_e\nu _0^2\theta _e^2.$$
(B3)
However, a large fraction of the low energy synchrotron photons are generally absorbed by synchrotron self-absorption. The synchrotron emission in the presence of absorption can be approximated as

$$\mathrm{\Lambda }_{syn}=f_{syn}\mathrm{\Lambda }_{syn}^0,$$
(B4)
where

$$f_{syn}_{x_M^{abs}}^{\mathrm{}}x_MI^{}(x_M)𝑑x_M/_0^{\mathrm{}}x_MI^{}(x_M)𝑑x_M$$
(B5)
is the fraction of synchrotron emission above the synchrotron self-absorption frequency $`\nu _{abs}`$ which satisfies

$$\tau _{syn}=\frac{1}{4\sqrt{3}}\frac{e^2cn_e(r)r}{\nu _{abs}kT_e\theta _e^2}I^{}(\nu _{abs})=1.$$
(B6)
Thus, we only consider the optically thin part of synchrotron emission as the cooling function for the gas and as a contribution to the preheating radiation field.

Locally emitted synchrotron photons are upscattered by inverse Comptonization off hot electrons. Here, we adopt a simple estimate of Comptonized synchrotron, which is reasonable in physical conditions considered in this paper,

$$\mathrm{\Lambda }_{CS}=\eta _{syn}f_{syn}\mathrm{\Lambda }_{syn}^0,$$
(B7)
where

$$\eta _{syn}=\mathrm{min}(\eta _0,\eta _{syn}^{sat})$$
(B8)
is the Comptonized energy enhancement factor (see e.g., Dermer et al. 1991) for synchrotron emission. The enhancement factor for fully saturated Comptonized synchrotron in the presence of absorption is

$$\eta _{syn}^{sat}=\frac{_{x_M^{abs}}^{\mathrm{}}\frac{3kT_e}{h\nu }x_MI^{}(x_M)𝑑x_M}{_{x_M^{abs}}^{\mathrm{}}x_MI^{}(x_M)𝑑x_M},$$
(B9)
(see §3.1 for bremsstrahlung case). The synchrotron absorption frequency is compared with that for free-free absorption, and the larger of the two is chosen.

The radiation temperature of locally emitted Comptonized synchrotron is similarly

$$T_X^{CS}=\mathrm{min}(\eta _0\frac{3}{8}\frac{h\nu _0}{k}\theta _e^2<x_M>_{ϵ_\nu },T_e),$$
(B10)
where

$$<x_M>_{ϵ_\nu }\frac{_{x_M^{abs}}^{\mathrm{}}x_M^2I^{}(x_M)𝑑x_M}{_{x_M^{abs}}^{\mathrm{}}x_MI^{}(x_M)𝑑x_M}$$
(B11)
is the energy-weighted mean photon energy in $`x_M`$ unit. This treatment of Comptonized synchrotron radiation is almost the same as that adopted in Park & Ostriker (2001), except that the energy enhancement factor and the treatment of saturated Comptonization are slightly improved.

In the presence of both bremsstrahlung and synchrotron emission, $`\mathrm{\Lambda }`$ is replaced by $`\mathrm{\Lambda }_{Cbr}+\mathrm{\Lambda }_{CS}`$ in equation (8) and $`T_X\mathrm{\Lambda }`$ by $`T_X^{Cbr}\mathrm{\Lambda }_{Cbr}+T_X^{CS}\mathrm{\Lambda }_{CS}`$ in equation (9).