text:
# Fluctuations of the Condensate in Ideal and Interacting Bose Gases

## I Introduction

The experimental achievement of Bose-Einstein condensation (BEC) in dilute alkali atoms , spin-polarized hydrogen  and recently in metastable helium  has enormously stimulated the theoretical research  on the ultracold bosons. In particular, fluctuations $`\delta ^2N_\mathrm{𝟎}`$ of the mean ground state occupation number $`N_\mathrm{𝟎}`$ have been recently thoroughly investigated in a series of papers. Apart from the intrinsic theoretical interest, it is foreseeable that such fluctuations will become experimentally testable in the near future .

It is well known that within grand canonical ensemble the fluctuations of the condensate are given by $`\delta ^2N_\mathrm{𝟎}=N_\mathrm{𝟎}\left(N_\mathrm{𝟎}+1\right)V^2`$, implying that $`\delta N_\mathrm{𝟎}`$ becomes of order $`N`$ when the temperature approaches zero. To avoid this sort of unphysically large condensate fluctuations, canonical (or microcanonical) ensemble has to be used to investigate the fluctuations of the condensate. Within microcanonical and canonical ensemble, the fluctuations of the condensate have been studied in a systematic way in the case of the ideal Bose gas . Recently, the question of how interatomic interactions affect the fluctuations of the condensate has been the object of several theoretical investigations . Giorgini et al.  found the anomalous behavior of the fluctuations in a weakly interacting Bose gas confined in a box within the traditional particle-number-nonconserving Bogoliubov approach. In  the fluctuations of the condensate follows the law $`\delta ^2N_\mathrm{𝟎}V^{4/3}`$. However, Idziaszek et al.  considered that the fluctuations are proportional to the volume. Recently, Kocharovsky et al.  supported and extended the results of the work of Giorgini et al.  using the particle-number-conserving operator formalism.

Although the correction to the ground state occupation number due to interatomic interaction has been clearly discussed within grand canonical ensemble  and canonical ensemble , the role of interaction on the condensate fluctuations of the weakly interacting Bose gas is still an open and unsolved problem. Different from the ground state occupation number, different models of describing the weakly interacting Bose gas will lead to vastly different prediction concerning the fluctuations of the condensate.

The purpose of this paper is to present a unified method to calculate the fluctuations of the condensate when different ways of approximation to the weakly interacting Bose gases are used. Within the canonical ensemble we give the distribution function of the ground state occupation number for the ideal and interacting Boson system in a box. We obtain the fluctuations of the condensate from the distribution function. In particular, we found that the distribution function is not Gaussian function in the case of interacting Boson system in a box. The paper is organized as follows. Sec.II is devoted to outline the canonical ensemble, which is developed to discuss the fluctuations of the condensate for the ideal Bose gas in a box. In Sec.III we investigate the fluctuations of the interacting Bose gas based on the lowest order perturbation theory. In Sec.IV the fluctuations are calculated based on the Bogoliubov theory. Finally, we give a discussion and summary of the results in Sec.V.

## II Mean Ground State Occupation Number and Fluctuations in the Ideal Bose Gases

Let us start our investigation on the fluctuations of the ideal Bose gases in the frame of canonical ensemble. According to the canonical ensemble the partition function of the $`N`$ non-interacting bosons in a box is given by

$$Z_{ideal}\left[N\right]=\underset{\mathrm{\Sigma }N_𝐧=N}{}\mathrm{exp}\left[\beta \left\{\mathrm{\Sigma }N_𝐧\epsilon _𝐧\right\}\right],$$
(1)
where $`N_𝐧`$ and $`\epsilon _𝐧`$ are the occupation numbers and energy level of the state $`𝐧=\{n_x,n_y,n_z\}`$ respectively. $`\beta =1/k_BT`$. In (1) the energy level of the system takes the form

$$\epsilon _𝐧=\frac{\pi ^2\mathrm{}^2\left(n_x^2+n_y^2+n_z^2\right)}{2mL^2}.$$
(2)
Separating out the ground state $`𝐧=\mathrm{𝟎}`$ from the state $`𝐧\mathrm{𝟎}`$, we have

$$Z_{ideal}\left[N\right]=\underset{N_\mathrm{𝟎}=0}{\overset{N}{}}\left\{\mathrm{exp}\left[\beta N_\mathrm{𝟎}\epsilon _\mathrm{𝟎}\right]Z_0\left(NN_\mathrm{𝟎}\right)\right\},$$
(3)
where $`Z_0\left(NN_\mathrm{𝟎}\right)`$ stands for the partition function of a fictitious system comprising $`NN_\mathrm{𝟎}`$ non-interacting bosons. Assuming $`A_0\left(NN_\mathrm{𝟎}\right)`$ stands for the free energy of the fictitious system,

$$A_0\left(NN_\mathrm{𝟎}\right)=k_BT\mathrm{ln}Z_0\left(NN_\mathrm{𝟎}\right)$$
(4)
From (3) and (4) the partition function $`Z_{ideal}\left[N\right]`$ becomes

$$Z_{ideal}\left[N\right]=\underset{N_\mathrm{𝟎}=0}{\overset{N}{}}\mathrm{exp}\left[q(N,N_\mathrm{𝟎})\right],$$
(5)
where $`q(N,N_\mathrm{𝟎})=\beta N_\mathrm{𝟎}\epsilon _\mathrm{𝟎}\beta A_0\left(NN_\mathrm{𝟎}\right)`$. Obviously $`\mathrm{exp}\left[q(N,N_\mathrm{𝟎})\right]/Z_{ideal}\left[N\right]`$ represents the probability to find $`N_\mathrm{𝟎}`$ atoms in the condensate. We will give the distribution function of the ground state occupation number in the following.

Let us first investigate the largest term in the sum of the partition function $`Z_{ideal}\left[N\right]`$. Assume the number of the condensed atoms is $`N_\mathrm{𝟎}^p`$ in the largest term of the partition function $`Z_{ideal}\left[N\right]`$. The largest term $`Z_0\left(NN_\mathrm{𝟎}^p\right)`$ is determined by the requirement that $`\frac{}{N_\mathrm{𝟎}}q(N,N_\mathrm{𝟎})|_{N_\mathrm{𝟎}=N_\mathrm{𝟎}^p}=0`$, ie.

$$\beta \epsilon _\mathrm{𝟎}\beta \frac{}{N_\mathrm{𝟎}^p}A_0\left(NN_\mathrm{𝟎}^p\right)=0$$
(6)
The calculations of the free energy $`A_0\left(NN_\mathrm{𝟎}^p\right)`$ is nontrivial because there is a requirement that the number of the particles is $`NN_\mathrm{𝟎}^p`$ in the summation of the partition function $`Z_0\left(NN_\mathrm{𝟎}^p\right)`$. Using the saddle-point method developed by Darwin and Fowler  it is straightforward to obtain the free energy $`A_0\left(NN_\mathrm{𝟎}^p\right)`$ of the fictitious system.

$$A_0\left(NN_\mathrm{𝟎}^p\right)=\left(NN_\mathrm{𝟎}^p\right)k_BT\mathrm{ln}z_\mathrm{𝟎}^pV\frac{k_BT}{\lambda ^3}g_{5/2}\left(z_\mathrm{𝟎}^p\right),$$
(7)
where $`\lambda =\sqrt{2\pi \beta \mathrm{}^2/m}`$ is the thermal wavelength. $`z_\mathrm{𝟎}^p`$ is the fugacity of the $`NN_\mathrm{𝟎}^p`$ non-interacting bosons and is determined by the equation

$$NN_\mathrm{𝟎}^p=\underset{_{𝐧\mathrm{𝟎}}}{}\frac{1}{\mathrm{exp}\left[\epsilon _𝐧/k_BT\right]\left(z_\mathrm{𝟎}^p\right)^11}=\frac{V}{\lambda ^3}g_{3/2}\left(z_\mathrm{𝟎}^p\right).$$
(8)
From (7) and (8) one finds

$$\beta \frac{}{N_\mathrm{𝟎}^p}A_0\left(NN_\mathrm{𝟎}^p\right)=\mathrm{ln}z_\mathrm{𝟎}^p.$$
(9)
Combining (6) and (9) one obtains $`\mathrm{ln}z_\mathrm{𝟎}^p=\beta \epsilon _\mathrm{𝟎}`$. Therefore the most probable value $`N_\mathrm{𝟎}^p`$ is determined by

$$N_\mathrm{𝟎}^p=N\underset{𝐧\mathrm{𝟎}}{}\frac{1}{\mathrm{exp}\left[\left(\epsilon _𝐧\epsilon _\mathrm{𝟎}\right)/k_BT\right]1}.$$
(10)
$`N_\mathrm{𝟎}^p`$ is exactly the mean occupation number of the condensate atoms in the frame of the grand canonical ensemble. For sufficiently large $`N`$, the sum $`_{N_\mathrm{𝟎}=0}^N`$ in (3) maybe replaced by the largest term, for the error omitted in doing so will be statistically negligible. In this case, (10) shows the equivalence between canonical ensemble and grand canonical ensemble for large $`N`$. From (10), below the critical temperature $`N_\mathrm{𝟎}^p`$ is determined by

$$N_\mathrm{𝟎}^p=N\frac{V}{\lambda ^3}\zeta \left(3/2\right)=N\left(1\left(\frac{T}{T_c^0}\right)^{3/2}\right),$$
(11)
where $`T_c^0=\frac{2\pi }{\left[\zeta \left(3/2\right)\right]^{2/3}}\frac{\mathrm{}^2}{mk_B}\left(\frac{N}{V}\right)^{2/3}`$ is the transition temperature of the ideal Bose gas.

Other terms in (5) will contribute to the fluctuations of the condensate, and lead to the deviation of the mean occupation number $`N_\mathrm{𝟎}`$ from the most probable value $`N_\mathrm{𝟎}^p`$. When $`N_\mathrm{𝟎}N_\mathrm{𝟎}^p`$, $`\frac{}{N_\mathrm{𝟎}}q(N,N_\mathrm{𝟎})0`$. Assuming

$$\frac{}{N_\mathrm{𝟎}}q(N,N_\mathrm{𝟎})=\alpha (N,N_\mathrm{𝟎}),$$
(12)
and repeating the saddle-point method, one obtains

$$N_\mathrm{𝟎}=N\underset{𝐧\mathrm{𝟎}}{}\frac{1}{\mathrm{exp}\left[\left(\epsilon _𝐧\epsilon _\mathrm{𝟎}\right)/k_BT\right]\mathrm{exp}\left[\alpha (N,N_\mathrm{𝟎})\right]1}.$$
(13)
From (11) and (13) one gets

$$\alpha (N,N_\mathrm{𝟎})=\frac{\lambda ^6}{V^2}\frac{\left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)^2}{4\pi }\theta \left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)$$
(14)
where $`\theta \left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)`$ is a Sign function. $`\theta \left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)=1`$ when $`N_\mathrm{𝟎}>N_\mathrm{𝟎}^p`$, and $`\theta \left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)=1`$ when $`N_\mathrm{𝟎}<N_\mathrm{𝟎}^p`$. To obtain (14) we have used the expansions $`g_{3/2}\left(1\delta \right)\zeta \left(3/2\right)2\sqrt{\pi \delta }`$  and the approximation $`\mathrm{exp}\left[\alpha (N,N_\mathrm{𝟎})\right]1\alpha (N,N_\mathrm{𝟎})`$. From (12) and (14) one obtains easily the following result for $`q(N,N_\mathrm{𝟎})`$.

$$q(N,N_\mathrm{𝟎})=_{N_\mathrm{𝟎}^p}^{N_\mathrm{𝟎}}\alpha (N,N_\mathrm{𝟎})𝑑N_\mathrm{𝟎}+q(N,N_\mathrm{𝟎}^p)$$

$$=\frac{\lambda ^6}{12\pi V^2}\left|N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right|^3+q(N,N_\mathrm{𝟎}^p).$$
(15)
The partition function $`Z_{ideal}\left[N\right]`$ is thus

$$Z_{ideal}\left[N\right]=\underset{N_\mathrm{𝟎}=0}{\overset{N}{}}\left\{\mathrm{exp}\left[q(N,N_\mathrm{𝟎}^p)\right]G_{ideal}(N,N_\mathrm{𝟎})\right\},$$
(16)
where we have introduced a distribution function $`G_{ideal}(N,N_\mathrm{𝟎})`$,

$$G_{ideal}(N,N_\mathrm{𝟎})=\mathrm{exp}\left[\frac{\lambda ^6}{12\pi V^2}\left|N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right|^3\right].$$
(17)
Assuming $`P\left(N_\mathrm{𝟎}|N\right)`$ is the probability to find $`N_\mathrm{𝟎}`$ atoms in the condensate, the distribution function $`G_{ideal}(N,N_\mathrm{𝟎})`$ represents the ratio $`\frac{P\left(N_\mathrm{𝟎}|N\right)}{P\left(N_\mathrm{𝟎}^p|N\right)}`$, ie. the relative probability to find $`N_\mathrm{𝟎}`$ atoms in the condensate. From (16) and (17) one obtains the mean occupation number $`N_\mathrm{𝟎}`$ and fluctuations $`\delta ^2N_\mathrm{𝟎}`$ within the canonical ensemble,

$$N_\mathrm{𝟎}=\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}G_{ideal}(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_{ideal}(N,N_\mathrm{𝟎})},$$
(18)
$$\delta ^2N_\mathrm{𝟎}=N_\mathrm{𝟎}^2N_\mathrm{𝟎}^2=$$

$$\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}^2G_{ideal}(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_{ideal}(N,N_\mathrm{𝟎})}\left(\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}G_{ideal}(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_{ideal}(N,N_\mathrm{𝟎})}\right)^2.$$
(19)
From (11), (17) and (18), (19) it is easy to obtain $`N_\mathrm{𝟎}`$ and $`\delta ^2N_\mathrm{𝟎}`$ of the non-interacting Bose gases in a box. At the critical temperature $`T_c^0`$, $`N_\mathrm{𝟎}^p=0`$. Thus $`G_{ideal}\left(T=T_c^0\right)=\mathrm{exp}\left[\frac{\lambda _0^6}{12\pi V^2}N_\mathrm{𝟎}^3\right]`$, where $`\lambda _0`$ is the thermal wavelength at $`T_c^0`$. From (18) and (19) one obtains the analytical result for the condensate fluctuations at $`T_c^0`$.

$$\delta N_\mathrm{𝟎}^2_{T=T_c^0}=$$

$$\left[\frac{1}{3\mathrm{\Gamma }\left(4/3\right)}\left(\frac{\mathrm{\Gamma }\left(5/3\right)}{2\mathrm{\Gamma }\left(4/3\right)}\right)^2\right]\left(\frac{12\pi }{\lambda _0^6}\right)^{2/3}V^{4/3},$$
(20)
where $`\mathrm{\Gamma }\left(n\right)=_0^{\mathrm{}}e^tt^{n1}𝑑t`$ is Gamma function. $`\mathrm{\Gamma }\left(4/3\right)=0.893`$ and $`\mathrm{\Gamma }\left(5/3\right)=0.903`$. (20) clearly shows that there is anomalous behavior for the fluctuations of the condensate. When $`T0`$, from (17) one finds $`G_{ideal}(N,N_\mathrm{𝟎})0`$ when $`N_\mathrm{𝟎}N`$. Therefore when $`T0`$ one obtains $`N_\mathrm{𝟎}N`$ and $`<\delta ^2N_\mathrm{𝟎}>0`$.

In Fig.1 we plot $`N_\mathrm{𝟎}/N`$ as a function of temperature for the ideal Bose gases in a box. The solid line displays the mean ground state occupation number within the grand canonical ensemble (or $`N_\mathrm{𝟎}^p`$). When $`N>10^4`$, the mean ground state occupation number of the canonical ensemble agrees well with that of the grand canonical ensemble. Obviously, in the case of $`N\mathrm{}`$, the mean ground state occupation number of the canonical ensemble coincides with that of the grand canonical ensemble.

In Fig.2 we plot numerical result of $`\delta N_\mathrm{𝟎}=\sqrt{\delta ^2N_\mathrm{𝟎}}`$ (thick solid line) for the ideal Bose gas with $`N=1000`$ atoms confined in a three-dimensional box. The dashed line shows the result of , where the fluctuations is given by

$$\delta ^2N_\mathrm{𝟎}=2A\left(\frac{mk_BT}{\mathrm{}^2}\right)^2V^{4/3}.$$
(21)
The coefficient in (21) is $`A=2/\pi ^4\times \mathrm{\Sigma }_{𝐧\mathrm{𝟎}}1/𝐧^4=0.105`$. The dashed line is larger than our result because of the approximation in . In  $`\delta ^2N_\mathrm{𝟎}=_{𝐧\mathrm{𝟎}}f_𝐧^2`$, where $`f_𝐧=\left(\mathrm{exp}\left(\epsilon _𝐧/k_BT\right)1\right)^1`$. For the convenience of calculations, $`f_𝐧`$ is approximated as $`\epsilon _𝐧/k_BT`$ for low energy atoms. However, this approximation is used also for the atoms whose energy level is larger than $`k_BT`$. Obviously, this approximation will lead to the fluctuations become larger. On the other hand, (21) holds in the canonical ensemble except near and above $`T_c`$, while our analysis holds also for the temperature near $`T_c^0`$. In Fig.2 the dotted line shows the numerical result of Wilkens et al. .

In Fig.2 the arrow marks the temperature $`T_m`$ which corresponds to the maximum fluctuations $`\delta ^2N_\mathrm{𝟎}_{\mathrm{max}}`$. Below the temperature $`T_m`$, from (19), one obtains the analytical result for the fluctuations of the condensate.

$$\delta ^2N_\mathrm{𝟎}=A\left(\frac{mk_BT}{\mathrm{}^2}\right)^2V^{4/3}.$$
(22)
It is interested to find that the coefficient differs by a factor $`2`$, compared to (21). The thin solid line shows (22) in Fig.2.

We should note that our results are reliable although the disputable saddle-point method is used to investigate the fluctuations of the condensate. It is well known that the applicability of the saddle-point approximation for the condensed Bose gases has been the subject of a long debate . Recently, the analysis in  showed that the fluctuations are overestimated, and do not appear to vanish properly with temperature using the usual saddle-point method. Our discussions for fluctuations are reasonable because of two reasons. (i) As proved in , the free energy (7) of the non-interacting Bose gases is still correct, even when carefully deal with the failure of the standard saddle-point method below the critical temperature. (ii) In the usual statistical method $`N_\mathrm{𝟎}`$ and $`\delta ^2N_\mathrm{𝟎}`$ are obtained through the first and second partial derivative of the partition function respectively. When saddle-point approximation is used to calculate the partition function of the system, the error will be overestimated in the second partial derivative of the partition function so that we can not obtain correct fluctuations of the condensate in the usual method. However, in this paper what we used is the reliable result (7). The distribution function of the ground state occupation number is obtained directly from (7), without resorting to the second partial derivative of the partition function. $`N_\mathrm{𝟎}`$ and $`\delta ^2N_\mathrm{𝟎}`$ are obtained from the distribution function in this paper.

## III Fluctuations of the Condensate Based on the Lowest Order Perturbation Theory

In the case of interacting Bose gases, the role of interactions on the fluctuations of the condensate is still an open and unsolved problem. Giorgini et al.  predicted the anomalous behavior of the fluctuations in a weakly interacting Bose gas confined in a box, while Idziaszek et al.  considered that the fluctuations are normal. Researches show that different model of approximation to the interacting Bose gases will lead to different predictions concerning the fluctuations of the condensate. The method developed to obtain the fluctuations of the ideal Bose gas in this paper can be used straightforwardly to discuss the fluctuations of the interacting Bose gases when different models of approximation are adopted.

Let us first discuss the fluctuations of the condensate in the case of the lowest order perturbation theory, which is also discussed in . In terms of the lowest order perturbation theory, the partition function of the system within the canonical ensemble is given by

$$Z_{int}\left[N\right]=\underset{\mathrm{\Sigma }N_𝐧=N}{}\mathrm{exp}\left[\beta \left(\mathrm{\Sigma }N_𝐧\epsilon _𝐧+E_{int}\right)\right],$$
(23)
where the interaction energy of the system takes the form

$$E_{int}=\frac{4\pi a\mathrm{}^2}{mV}\left(N^2\frac{1}{2}N_\mathrm{𝟎}^2\right).$$
(24)
In (24) $`a`$ is the scattering length. Separating out the ground state $`𝐧=\mathrm{𝟎}`$ from the state $`𝐧\mathrm{𝟎}`$, one obtains the following form for the partition function

$$Z_{int}\left[N\right]=\underset{N_\mathrm{𝟎}=0}{\overset{N}{}}\left\{\mathrm{exp}\left[\beta N_\mathrm{𝟎}\epsilon _\mathrm{𝟎}\beta E_{int}\right]Z_0\left(NN_\mathrm{𝟎}\right)\right\},$$
(25)
where $`Z_0\left(NN_\mathrm{𝟎}\right)`$ stands for the partition function of a fictitious $`NN_\mathrm{𝟎}`$ non-interacting bosons. Using the free energy $`A_0\left(NN_\mathrm{𝟎}\right)`$ of the fictitious system, the partition function is thus

$$Z_{int}\left[N\right]=\underset{N_\mathrm{𝟎}=0}{\overset{N}{}}\mathrm{exp}\left[q(N,N_\mathrm{𝟎})\right],$$
(26)
where $`q(N,N_\mathrm{𝟎})`$ takes the form

$$q(N,N_\mathrm{𝟎})=\beta N_\mathrm{𝟎}\epsilon _\mathrm{𝟎}\beta E_{int}\beta A_0\left(NN_\mathrm{𝟎}\right).$$
(27)
Analogously to the case of ideal Bose gases, let us first investigate the largest term in the sum of $`Z_{int}\left[N\right]`$. The largest term is determined by the requirement $`\frac{}{N_\mathrm{𝟎}}q(N,N_\mathrm{𝟎})|_{N_\mathrm{𝟎}=N_\mathrm{𝟎}^p}=0`$. Therefore one obtains the most probable value $`N_\mathrm{𝟎}^p`$ of the interacting Bosons.

$$N_\mathrm{𝟎}^p=N\underset{𝐧\mathrm{𝟎}}{}\frac{1}{\mathrm{exp}\left[\beta \epsilon _𝐧\right]\left(z_\mathrm{𝟎}^p\right)^11},$$
(28)
where $`z_\mathrm{𝟎}^p`$ is determined by

$$lnz_\mathrm{𝟎}^p=\beta \epsilon _\mathrm{𝟎}+\beta \frac{}{N_\mathrm{𝟎}^p}E_{int}=\beta \epsilon _\mathrm{𝟎}\frac{2a\lambda ^2N_\mathrm{𝟎}^p}{V}$$
(29)
From (28) and (29) one obtains

$$N_\mathrm{𝟎}^pN\frac{V}{\lambda ^3}\left[\zeta \left(3/2\right)2\sqrt{\pi }\left(\frac{2a\lambda ^2N_\mathrm{𝟎}^p}{V}\right)^{1/2}\right]$$
(30)
Other terms in (26) will contribute to the fluctuations of the system. Assuming $`\frac{}{N_\mathrm{𝟎}}q(N,N_\mathrm{𝟎})=\alpha (N,N_\mathrm{𝟎})`$, one obtains the result for $`N_\mathrm{𝟎}`$

$$N_\mathrm{𝟎}=N\underset{𝐧\mathrm{𝟎}}{}\frac{1}{\mathrm{exp}\left[\beta \epsilon _𝐧\right]\left(z_\mathrm{𝟎}\right)^11},$$
(31)
where $`z_\mathrm{𝟎}`$ is determined by

$$\mathrm{ln}z_\mathrm{𝟎}=\beta \epsilon _\mathrm{𝟎}\frac{2a\lambda ^2N_\mathrm{𝟎}}{V}+\alpha (N,N_\mathrm{𝟎})$$
(32)
From (28), (29) and (31), (32) it is straightforward to obtain the distribution function $`G_{int}(N,N_\mathrm{𝟎})`$ of the interacting Bose gases.

$$G_{int}(N,N_\mathrm{𝟎})=G_{ideal}(N,N_\mathrm{𝟎})\times R_{int}(N,N_\mathrm{𝟎}),$$
(33)
where $`G_{ideal}(N,N_\mathrm{𝟎})`$ is the distribution function (17) of the ideal Bose gases, while $`R_{int}(N,N_\mathrm{𝟎})`$ takes the form

$$R_{int}(N,N_\mathrm{𝟎})=R_1(N,N_\mathrm{𝟎})\times R_2(N,N_\mathrm{𝟎})\times R_3(N,N_\mathrm{𝟎}).$$
(34)
In (34),

$$R_1(N,N_\mathrm{𝟎})=$$

$$\mathrm{exp}\left[\frac{\left(\zeta \left(3/2\right)\right)^{3/2}}{\sqrt{2\pi }}\left(\frac{a}{\lambda _0}\frac{N_\mathrm{𝟎}^p}{N}\right)^{1/2}\frac{\left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)^2}{Nt^2}\theta \left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)\right],$$
(35)
$$R_2(N,N_\mathrm{𝟎})=\mathrm{exp}\left[\frac{\zeta \left(3/2\right)a}{\lambda _0}\frac{N_\mathrm{𝟎}^2\left(N_\mathrm{𝟎}^p\right)^2}{Nt}\right],$$
(36)
$$R_3(N,N_\mathrm{𝟎})=\mathrm{exp}\left[\left|\frac{2\zeta \left(3/2\right)a}{\lambda _0}\frac{N_\mathrm{𝟎}^p\left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)}{Nt}\right|\right].$$
(37)
We should note that $`G_{int}(N,N_\mathrm{𝟎})`$ is not a Gaussian distribution function because of the non-Gaussian factors $`R_1(N,N_\mathrm{𝟎})`$ and $`R_2(N,N_\mathrm{𝟎})`$, while Idziaszek et al.  utilized the Gaussian distribution as an assumption to investigate the fluctuations of the interacting system. In $`R_{int}(N,N_\mathrm{𝟎})`$, $`R_1(N,N_\mathrm{𝟎})`$ comprises the factor $`\left(a/\lambda _0\right)^{1/2}`$ and represents the leading correction to the distribution function due to interatomic interaction, while $`R_2(N,N_\mathrm{𝟎})`$ and $`R_3(N,N_\mathrm{𝟎})`$ are high order correction to the distribution function. We should note that the leading contribution $`R_1(N,N_\mathrm{𝟎})`$ is not a Gaussian function.

From the distribution function $`G_{int}(N,N_\mathrm{𝟎})`$ the mean occupation number and fluctuations of the condensate are determined by

$$N_\mathrm{𝟎}=\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}G_{int}(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_{int}(N,N_\mathrm{𝟎})},$$
(38)
$$\delta ^2N_\mathrm{𝟎}=$$

$$\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}^2G_{int}(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_{int}(N,N_\mathrm{𝟎})}\left(\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}G_{int}(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_{int}(N,N_\mathrm{𝟎})}\right)^2.$$
(39)
From (30), (33) and (38), (39) we can obtain the fluctuations of the interacting Boson gases. In Fig.3 we give the numerical result for the repulsive interactions with $`a/\lambda _0=1\times 10^3`$, $`2\times 10^3`$, $`5\times 10^3`$, $`1\times 10^2`$. The crossover from interacting to ideal Bose gases (thick solid line) is clearly demonstrated in Fig.3, while in  the fluctuations (thin solid line) of the interacting Bose gases are irrelative to the scattering length. When $`a0`$ it is easy to recover the fluctuations of the ideal Bose gases. According to (39) the leading contributions to the fluctuations are anomalous, ie. proportional to $`N^{4/3}`$, while there are also normal contributions proportional to $`N`$ due to interatomic interaction. This conclusion contradicts with that of  which predicts normal behavior, where the lowest order perturbation theory is also used to discuss the fluctuations of the condensate.

## IV Fluctuations of the Condensate Based on Bogoliubov Theory

Let us investigate the fluctuations of the condensate in the framework of Bogoliubov theory of a uniform weakly interacting Bose gas confined in a box. According to Bogoliubov theory , the total number of particles out of the condensate is given by

$$N_T=\underset{𝐧\mathrm{𝟎}}{}N_𝐧=\underset{𝐧\mathrm{𝟎}}{}\left(u_𝐧^2+v_𝐧^2\right)f_𝐧,$$
(40)
where

$$u_𝐧^2+v_𝐧^2=\frac{\left(\left(\epsilon _𝐧^B\right)^2+g^2n_0^2\right)^{1/2}}{2\epsilon _𝐧^B},$$
(41)
$$u_𝐧v_𝐧=\frac{gn_0}{2\epsilon _𝐧^B},$$
(42)
and $`f_𝐧`$ is the number of quasi-particles present in the system at the thermal equilibrium.

$$f_𝐧=\frac{1}{\mathrm{exp}\left[\epsilon _𝐧^B/k_BT\right]1}$$
(43)
In addition, the energy of the quasi-particles entering (41) and (42) is given by the well known Bogoliubov spectrum

$$\epsilon _𝐧^B=\left(\left(\epsilon _𝐧+gn_0\right)^2g^2n_0^2\right)^{1/2},$$
(44)
where $`g=4\pi \mathrm{}^2a/m`$ is the coupling constant, and $`n_0=N_0/V`$ is the condensate density. At low $`\left|𝐧\right|=\sqrt{n_x^2+n_y^2+n_z^2}`$, one obtains $`u_𝐧^2v_𝐧^21/\left|𝐧\right|`$ and $`f_𝐧1/\left|𝐧\right|`$. This results in $`1/\left|𝐧\right|^2`$ divergence in (43) at low $`\left|𝐧\right|`$. Although this sort of divergence will not lead to large contribution to the number of low energy quasi-particles, it gives leading contribution to the fluctuations of the condensate, as pointed out in . We will investigate the fluctuations due to low energy quasi-particles in the following.

In (40) $`N_𝐧`$ can be regarded as the effective occupation number of the thermal atoms, while

$$N_𝐧^B=\frac{N_𝐧}{u_𝐧^2+v_𝐧^2}=f_𝐧$$
(45)
is the occupation number of the quasi-particles. From the form of $`f_𝐧`$, we can construct the partition function of the quasi-particles in the frame of canonical ensemble.

$$Z_B=\underset{\left\{𝐧\right\}}{}\mathrm{exp}\left[\beta \underset{𝐧}{}N_𝐧^B\epsilon _𝐧^B\right].$$
(46)
From (45) $`Z_B`$ becomes

$$Z_B=\underset{\left\{\mathrm{\Sigma }N_𝐧=N\right\}}{}\mathrm{exp}\left[\beta \underset{𝐧}{}N_𝐧\epsilon _𝐧^{eff}\right],$$
(47)
where $`\epsilon _𝐧^{eff}=\epsilon _𝐧^B/\left(u_𝐧^2+v_𝐧^2\right)`$ can be regarded as an effective energy level of the thermal atoms. In this case $`Z_B`$ is the partition function of a fictitious Boson system comprising $`N`$ non-interacting Bosons whose energy level is determined by $`\epsilon _𝐧^{eff}`$. From (47) the most probable value $`N_\mathrm{𝟎}^p`$ is given by

$$N_\mathrm{𝟎}^p=N\underset{𝐧\mathrm{𝟎}}{}\frac{1}{\mathrm{exp}\left[\left(\epsilon _𝐧^{eff}\epsilon _\mathrm{𝟎}^{eff}\right)/k_BT\right]1},$$
(48)
Obviously the occupation number of low $`\left|𝐧\right|`$ in (48) coincides with that of (40). Analogously, other $`N_\mathrm{𝟎}`$ is thus

$$N_\mathrm{𝟎}=N$$

$$\underset{𝐧\mathrm{𝟎}}{}\frac{1}{\mathrm{exp}\left[\left(\epsilon _𝐧^{eff}\epsilon _\mathrm{𝟎}^{eff}\right)/k_BT\right]\mathrm{exp}\left[\alpha (N,N_\mathrm{𝟎})\right]1}.$$
(49)
From (48) and (49) one gets

$$\alpha (N,N_\mathrm{𝟎})\frac{N_\mathrm{𝟎}N_\mathrm{𝟎}^p}{_{𝐧\mathrm{𝟎}}\left(u_𝐧^2+v_𝐧^2\right)^2f_𝐧^2},$$
(50)
where we have used the approximation $`f_𝐧k_BT/\epsilon _𝐧^B`$ for low energy quasi-particles. Therefore the Gaussian distribution function of the system is given by

$$G_B(N,N_\mathrm{𝟎})=\mathrm{exp}\left[\frac{\left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)^2}{2_{𝐧\mathrm{𝟎}}\left(u_𝐧^2+v_𝐧^2\right)^2f_𝐧^2}\right]$$

$$\mathrm{exp}\left[\frac{\left(\zeta \left(3/2\right)\right)^{4/3}\left(N_\mathrm{𝟎}N_\mathrm{𝟎}^p\right)^2}{\left(2\pi \right)^2AN^{4/3}t^2}\right].$$
(51)
Obviously the mean occupation number $`N_\mathrm{𝟎}`$ and fluctuations $`\delta ^2N_\mathrm{𝟎}`$ is given by

$$N_\mathrm{𝟎}=\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}G_B(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_B(N,N_\mathrm{𝟎})},$$
(52)
$$\delta ^2N_\mathrm{𝟎}=$$

$$\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}^2G_B(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_B(N,N_\mathrm{𝟎})}\left(\frac{_{N_\mathrm{𝟎}=0}^NN_\mathrm{𝟎}G_B(N,N_\mathrm{𝟎})}{_{N_\mathrm{𝟎}=0}^NG_B(N,N_\mathrm{𝟎})}\right)^2.$$
(53)
From (51) and (52), (53) one obtains the fluctuations of the condensate based on Bogoliubov theory. At the critical temperature, $`G_B\left(T=T_c\right)=\mathrm{exp}\left[N_\mathrm{𝟎}^2/\theta \right]`$, where $`\theta =2_{𝐧\mathrm{𝟎}}\left(u_𝐧^2+v_𝐧^2\right)^2f_𝐧^2=\left(2\pi \right)^2AN^{4/3}/\left(\zeta \left(3/2\right)\right)^{4/3}`$. In this case, we obtains the analytical result of the condensate fluctuations.

$$\delta ^2N_\mathrm{𝟎}_{T=T_c}=\left(\frac{1}{2}\frac{1}{\pi }\right)\theta =\left(\frac{1}{2}\frac{1}{\pi }\right)\frac{\left(2\pi \right)^2A}{\left(\zeta \left(3/2\right)\right)^{4/3}}N^{4/3}.$$
(54)
(54) clearly shows that the anomalous behavior of the condensate fluctuations originates from the low energy quasi-particles, which gives the anomalous factor $`N^{4/3}`$ through $`\theta `$. In Fig.4 the solid line displays our results based on the Bogoliubov theory, while the dashed line shows the result of .

## V Discussion and Conclusion

In this paper we investigate the fluctuations of the condensate in a weakly interacting Bose gas confined in a box. Canonical ensemble is developed to calculate the fluctuations of the condensate when different models of interacting Bose gases are used. We found that both the lowest order perturbation theory and Bogoliubov theory give anomalous behavior of the fluctuations for the interacting Bose gases confined in a box.

Different from the usual method, the distribution function $`P\left(N_\mathrm{𝟎}|N\right)/P\left(N_\mathrm{𝟎}^p|N\right)`$ (ie. the ratio of the probability between $`N_\mathrm{𝟎}`$ and the most probable value $`N_\mathrm{𝟎}^p`$ ) of the ground state occupation number is obtained directly to calculate the fluctuations of the condensate. From some senses, we give a simple method to recover the applicability of the saddle-point approximation to discuss the condensate fluctuations, through the avoidance of the second derivative in the usual method.

For the present experiments of BEC, the harmonically trapped atoms are in a situation of almost complete isolation with the outer environment surrounding the trap, therefore canonical (or microcanonical) ensemble should be used to calculate the fluctuations of the condensate. On the other hand, one obtains more accurate mean ground state occupation number within the canonical ensemble, compared to the grand canonical ensemble. The present work may serve as another method to investigate the thermodynamic properties of the harmonically trapped interacting Bose gases such as critical temperature, condensate fraction, and fluctuations of the condensate.

The remain challenge is to extend the idea of this paper to the case of microcanonical ensemble where the energy of the system is also invariant. In addition, the role of interactions on the fluctuations of the condensate are expected to be much more dramatic in the case of attractive forces. We will investigate these problems in a subsequent work.

## Acknowledgments

This work was supported by the Science Foundation of Zhijiang College, Zhejiang University of Technology. This work was also supported by the National Natural Science Foundation of China.


format: text

text:
# Event anisotropy in 4.2A GeV/c C+C collisions

## Abstract

The directed and elliptic flow of protons and negative pions in 4.2A GeV/c C+C collisions is studied using the Fourier analysis of azimuthal distributions. It is found that the protons exhibit pronounced directed flow, while the flow of pions is either non existent or too weak to be detected experimentally. Also, it is found that in the entire rapidity interval the elliptic flow is very small if not zero. These results are confirmed by the Quark-Gluon-String Model (QGSM) and the relativistic transport model (ART 1.0), except that these models predict very weak antiflow of pions. The more detailed comparison with the QGSM suggests that the decay of resonances and rescattering of secondaries dominantly determine the proton and negative pion flow at this energy.

Event anisotropy, often called flow, has been observed in heavy-ion collisions at all incident energies . At Bevalac energies and below, the flow is usually studied in terms of the mean in-plane component of transverse momentum at a given rapidity, $`p^x(y)`$,  and additionally quantified in terms of derivative at the midrapidity $`F_y=dp^x/dy`$. At high energies, the Fourier expansion of the azimuthal distribution of particles is used . In this expansion the first harmonic, $`v_1`$, quantifies the directed flow while the second harmonic, $`v_2`$, quantifies the elliptic flow. Using the Fourier expansion, the anisotropic transverse flow was analysed for heavy symmetric systems at the AGS  and SPS  energies. It was found that this anisotropy, and particularly the elliptic flow, plays an important role for investigating properties of hadronic matter . However, it is still not clear whether the experimentally observed event anisotropy is of a dynamic origin or is due to the shadowing of spectator matter, passing time, etc.

In this paper the directed and elliptic flow of protons and negative pions in 4.2A GeV/c C+C collisions is studied using the Fourier analysis of azimuthal distributions. The analysis is performed using 9500 C+C semicentral and central collisions obtained with the 2-m propane bubble chamber, exposed at JINR, Dubna synchrophasotron. The data for semicentral and central collisions roughly correspond to the upper $`50\%`$ of the inelastic cross-section. Additionally, the same type of analysis is performed using the 400000 events generated by the QGSM  and 200000 events generated by the relativistic transport model (ART 1.0) . For these events the same centrality criterion is applied as in experiment, leading to the average impact parameter $``$ 2.6. In the 2-m propane bubble chamber practically all charged reaction products are detected. Negative particles, except identified electrons are considered to be $`\pi ^{}`$. Among them remains admixture of unidentified fast electrons ($`<`$ 5%). All positive particles with momenta less than 0.5 GeV/c are classified either as protons or $`\pi ^+`$ mesons according to their ionisation density and range. Positive particles above 0.5 GeV/c are taken to be protons, and because of this, the admixture of $`\pi ^+`$ of approximately 18% is subtracted statistically using the $`\pi ^+`$ and $`\pi ^{}`$ momentum distributions. From the resulting number of protons, the projectile spectators (protons with momenta $`p>`$ 3 GeV/c and emission angle $`\theta <4^\mathrm{o}`$) and target spectators (protons with momenta $`p<`$ 0.3 GeV/c) are further subtracted. The resulting number of participant protons still contains some 4% of deuterons (with $`p>`$ 0.48 GeV/c) which are statistically subtracted. The admixture of tritons (with $`p>`$ 0.65 GeV/c) is not considered. The experimental data are also corrected to the loss of particles emitted at small angles relative to the optical axes of chamber. The aim of this correction is to obtain isotropic distribution in azimuthal angle and smooth distribution in emission angle (both measured with respect to the direction of the incoming projectile).

The azimuthal distribution of particles may be represented with the first three terms of the corresponding Fourier expansion

$$\frac{dN}{d\varphi }\frac{1}{2\pi }\left[1+2v_1\mathrm{cos}(\varphi )+2v_2\mathrm{cos}(2\varphi )\right],$$
(1)
where the two coefficients, $`v_1`$ and $`v_2`$, quantify the directed and elliptic flow via $`v_1=cos(\varphi )`$ and $`v_2=cos(2\varphi )`$. In Eq. (1), $`\varphi =\varphi _{lab}\mathrm{\Phi }_{plane}`$ is the particle azimuthal angle determined with respect to the reaction plane, with $`\varphi _{lab}`$ denoting the azimuthal angle of particle in the laboratory frame and $`\mathrm{\Phi }_{plane}`$ denoting the azimuthal angle of the (true) reaction plane. Since both the projectile momentum and the impact parameter vectors are available in the QGSM simulation, they are used to determine the corresponding reaction plane. In the experiment the reaction plane is determined, for each event, using the projectile momentum vector and the vector $`𝐐`$ determined from

$$𝐐=\underset{i}{}𝐩_{Ti}(y>y_{cm}+\delta )\underset{j}{}𝐩_{Tj}(y<y_{cm}\delta ),$$
(2)
where $`𝐩_T`$ represents the transverse momentum of the proton emitted in the forward ($`y>y_{cm}+\delta `$), or backward ($`y<y_{cm}\delta `$), hemisphere. Here, $`y_{cm}`$ denotes the center of mass rapidity while the quantity $`\delta `$ (=0.2) removes the protons emitted around the $`y_{cm}`$ which are not contributing to the determination of the reaction plane. The reaction plane angle for a proton is determined using this expression only if this proton is not included in the above sum (i.e. if its rapidity lies in the interval from $`y_{cm}\delta `$ to $`y_{cm}+\delta `$). Otherwise, in order to avoid autocorrelation (which is an effect of the finite multiplicity), the $`𝐐`$ vector is constructed by the analogous expression in which the contribution of this proton is simply omitted . We found that the reaction plane angle distribution is essentially flat, thus confirming the absence of significant distorsions which could influence the magnitude of the extracted flow parameters.

The accuracy with which the reaction plane angle is determined, i.e. the reaction plane resolution, is evaluated by the subevent method . In this method, each event is divided randomly into two subevents, and then the corresponding two reaction planes are determined. Subsequently, the absolute value of the relative azimuthal angle, $`\mathrm{\Phi }_{12}`$, between these two estimated reaction planes is obtained. The width, $`\sigma `$, of the $`\mathrm{\Phi }_{12}`$ distribution determines the reaction plane resolution. For C+C collisions we find $`\sigma =50^0`$. The relative azimuthal angle distribution is the basis for the correction of the Fourier coefficients, $`v_n^{}`$, obtained with the estimated reaction plane. The relationship between the $`v_n^{}`$, and the Fourier coefficients $`v_n`$ obtained relative to the true reaction plane, is $`v_n^{}=v_ncos(n\mathrm{\Delta }\mathrm{\Phi })`$, where $`cos(n\mathrm{\Delta }\mathrm{\Phi })`$ is the correction factor determined from $`\mathrm{\Phi }_{12}`$ distribution following the prescription given in . We find $`cos(\mathrm{\Delta }\mathrm{\Phi })`$=0.56 and $`cos(2\mathrm{\Delta }\mathrm{\Phi })`$=0.24. The correctness of this procedure is checked using the QGSM. Using this model, the coefficients $`v_1`$ and $`v_2`$ are calculated with respect to the true reaction plane and also with respect to the estimated reaction plane. The results of the comparison will be discussed below.

Fig. 1 (top) displays the experimentally determined $`v_1`$ coefficient vs. $`y`$ (with $`y`$ calculated in the center-of-mass frame), for protons and negative pions together with the $`v_1`$ calculated with QGSM relative to the true reaction plane and relative to the estimated reaction plane. For the proper comparison with the experiment, we excluded protons satisfying cuts for the proton spectators in the experiment. In the case of protons it is seen that the values of the two QGSM results for $`v_1`$ are quite close. The dependence of $`v_1`$ on rapidity is characterised by a curve with a positive slope and with the zero-crossing at $`y=0`$. The curve indicates a positive directed flow with magnitude $`v_10.17`$, at rapidities close to the beam rapidity ($`0.7<y<1.5`$). The QGSM reproduces satisfactorily the shape of $`v_1(y)`$ curve and the magnitude of the flow. Using the extracted values of $`v_1`$ and their relation to the mean transverse momentum projected onto the reaction plane, $`v_1=p_x/p_T`$, we can evaluate $`p_x`$ as a function of rapidity and determine the slope, $`F=dp_x/d(y/y_b)`$, with respect to rapidity normalised to beam rapidity in the lab frame ($`y_b`$=2.2). In the present analysis we find for the slope at the midrapidity $`F=`$ 144 MeV/c. After the normalisation to the mass number of the colliding system we obtain the so called scaled flow $`F_S=F/(A_1^{1/3}+A_2^{1/3})=`$ 31 MeV/c. This value is in agreement with the observed trend  that after reaching the maximum at beam energy around 0.7-2A GeV, the directed flow slowly decreases with increasing beam energy.

For negative pions the experimental values of $`v_1`$ indicate that the directed flow is non-existent. This result is confirmed by the model calculations of $`v_1`$ with respect to the estimated reaction plane. However, the model calculations of $`v_1`$ with respect to the true reaction plane show the existence of a very weak directed flow of pions with the sign of $`v_1`$ opposite to that of protons and with the maximum value of 0.02 around the target rapidity. This further suggests that in the collisions of light nuclei, like C+C, the very weak flow, if it exists, is not measurable because of the limited accuracy in determination of the reaction plane.

Fig. 1 (bottom) displays the experimentally determined $`v_2^{}`$ coefficient vs. $`y`$ for protons and negative pions. This coefficient is not corrected to the reaction plane resolution since the comparison of the model calculation of $`v_2^{}`$, obtained as in the experiment, and the model calculation of $`v_2`$, relative to the true reaction plane, indicates that the corresponding correction procedure for $`v_2`$, as outlined above, is not applicable. The reason for this is the lightness of the colliding nuclei and the smallness of the elliptic flow. The uncorrected values of $`v_2^{}`$ show that in the entire rapidity interval the elliptic flow is small ($`|v_2^{^{}}|0.02`$) if not zero, and this is consistent with the predictions of QGSM. The experimental values for $`v_2^{}`$ also show that, for both protons and pions, the elliptic flow depends on rapidity and that around the beam and target rapidities it is positive for protons and negative for pions. The positive sign for protons indicates an enhanced emission in the reaction plane, while the negative sign for pions indicates an enhanced emission perpendicular to the reaction plane. This behaviour points out to the shadowing by the nuclear matter as the origin of the elliptic flow.

Since the QGSM predictions are in fair agreement with the experimental results at 4.2A GeV/c, we use this model to clarify the question which of the processes are responsible for the flow effect. In this model, in collisions of light C+C nuclei, approximately $`40\%`$ of protons and $`70\%`$ of $`\pi ^{}`$ originate from decay of the lowest-lying resonances ($`\mathrm{\Delta }^{}`$s, $`\varrho ,\omega ,\eta `$ and $`\eta ^{^{}}`$). The rest originates from the ’non-resonant’ primary and secondary interactions of the type: $`NNNN\pi `$, $`\mathrm{\Delta }N\mathrm{\Delta }N`$, $`\pi N\pi N`$, $`\pi NNNN`$. The protons and pions from primary interaction escape the collision zone without further rescattering and comprise $`5\%`$ of the total. Therefore, according to QGSM, we separately evaluate the flow of protons and pions originating from the following sources:(i) decay of resonances, (ii ) primary non-resonant interactions (iii ) and secondary non-resonant interactions.

Figure. 2 (top) shows $`v_1`$ vs. rapidity for protons and negative pions originating from decay of resonances, and from primary and secondary non-resonant interactions, as well as the overall $`v_1`$ for protons and $`\pi ^{}`$. (In these model calculations, the experimental cuts for the proton spectators were not applied, and this leads to a small difference between the two curves for overall $`v_1`$ for protons in Figs. 1 and 2). The protons originating both from the decay of resonances and from the secondary interactions show the directed flow of similar intensity. The same applies to the antiflow of pions. The protons from the primary interactions show a relatively flat $`v_1(y)`$ distribution, while the pions from these interactions show a strong directed antiflow with magnitude $`0.13`$ around beam rapidity. The antiflow of pions can be also explained by the shadowing effect, but the shadowing matter is different for pions from primary and secondary interactions since these pions are produced at different collision stages. Initially (at times less than the passing time, $`t_p`$=4.2 fm/c) the pions are shadowed by the cold spectators. Later, after the spectator matter leaves the collision zone, the pions are shadowed by the participant nucleons. This may be the underlying mechanism that leads to the different behaviour of $`v_1`$ for pions. In heavy nuclei collisions, additionally generated by the QGSM, the protons from primary interactions, similarly to the case of pions from primary interactions, show strong antiflow caused by the shadowing from the cold nuclear matter. In the collisions of light nuclei this shadowing is small and there is no preferential emission of the protons.

Figure.2 (bottom) shows $`v_2`$ vs. rapidity for protons and negative pions originating from decay of resonances, and from primary and secondary non-resonant interactions, as well as the overall $`v_1`$ for protons and $`\pi ^{}`$. It is seen that the particles from secondary interactions and from decay of resonaces, exhibit similar behaviour. The particles from the primary interactions show a clear negative elliptic flow, and this out-of-plane emission can be attributed to the shadowing by the cold spectators.

In order to establish a less model dependent picture, the results of the experiment are also compared with the relativistic transport model, ART 1.0. These are shown in Fig. 3, where the calculations are performed both in the cascade and in the so-called ’mean field’ mode. The cascade mode underestimates the magnitude of the proton flow (some 20$`\%`$), and predicts a small directed antiflow of negative pions with magnitude $`v_10.04`$. In the mean field mode the model increases the magnitude of the proton flow. Also, the ART model predicts a very small proton and pion elliptic flow ($`|v_2|0.01`$).

In summary, the directed and elliptic flow of protons and negative pions in 4.2A GeV/c C+C collisions was examined using the Fourier analysis of azimuthal distributions of experimental events, and also by using the events generated by the QGSM and ART 1.0 model. It was found that the protons exhibit strong directed flow with magnitude $`v_10.17`$ at rapidities close to the beam rapidity. The QGSM reproduces satisfactorily the shape of the $`v_1(y)`$ curve and the magnitude of the flow. The ART model underestimates this magnitude in the cascade mode and increases this magnitude in the mean field mode. For negative pions the flow is either non existent, or too weak to be detected experimentally due to the limited accuracy in the determination of the reaction plane. The latter was suggested by the QGSM, where calculations with respect to the estimated reaction plane predicted non existant flow, while the calculations with respect to the true reaction plane predicted a small directed antiflow with magnitude $`v_10.02`$. The predictions of the ART model are similar. Also, it was found that in the entire rapidity interval the elliptic flow is small ($`|v_2^{^{}}|0.02`$) if not zero, and this is in agreement with the predictions of QGSM and ART 1.0. According to the QGSM, the two factors that dominantly determine the proton and negative pion flow, at this energy, are the decay of resonances and the rescattering of secondaries. The shadowing by the cold spectator matter affects only the flow of the particles produced at the early stage of the collision.

The authors are grateful to members of the JINR Dubna group that participated in data processing, G. Škoro for making the QGSM code available to us, and I. Mendaš for useful comments.


text:
# Transition from random to ordered fractals in fragmentation of particles in an open system

## Abstract

We consider the fragmentation process with mass loss and discuss self-similar properties of the arising structure both in time and space, focusing on dimensional analysis. This exhibits a spectrum of mass exponents $`\theta `$, whose exact numerical values are given for which $`x^\theta `$ or $`t^{\theta z}`$ has the dimension of particle size distribution function $`\psi (x,t)`$, where $`z`$ is the kinetic exponent. We obtained conditions for which the scaling and fragmentation process altogether break down and give explicit scaling solution for special case. Finally, we identify a new class of fractals ranging from random to non-random and show that the fractal dimension increases with increasing order and a transition to strictly self-similar pattern occurs when randomness completely ceases.

PACS number(s): 05.45.Df, 46.50.+a, 62.20.Mk

The kinetics of the irreversible and sequential breakup of particles occurs in a variety of physical processes and has important applications in science and technology. These include erosion , grinding and crushing of solids, polymer degradation and fiber length reduction , breakup of liquid droplets , etc. to name just a few. In recent years, there has been an increasing interest in studying fragmentation, allowing variations that increase the flexibility of the theory in matching such conditions of real phenomena such as extension to higher dimension , agglomerate erosion , mass loss , volume change  and fragmentation-annihilation . The kinetic equation approach of fragmentation is described by a linear integro-differential equation where mass or size is the only dynamical quantity. This theoretical approach is mean field in character since fluctuations are ignored altogether.

In this article, we study the kinetics of fragmentation with continuous mass loss, an interesting variant of the classical fragmentation equation introduced by Edwards et al. . This is relevant in all fragmentation processes mentioned earlier where mass loss might occur due to evaporation, oxidation, sublimation, dissolution, melting etc., or in Yule-Furry process of cosmic shower theory where energy loss occurs due to collision . Although there exist a series of paper devoted to this problem, the explicit scaling solution with exact numerical exponents, the geometric properties of the arising pattern and the interplay between fragmentation and mass loss remain unexplored. In addition to obtaining these, we give an alternative interpretation of the existing scaling theory by considering dimensional analysis, and we give an explicit bound to the exponent of the breakup rate for which the kinetic equation fails. A strong motivation for the present investigation came from the desire to know how fractal geometry of the resulting object changes with the degree of order, which is quantified by the global exponent and is typically known as fractal dimension.

The evolution of the particle size distribution function $`\psi (x,t)`$ for fragmentation with mass loss in one dimension is

$$\frac{\psi (x,t)}{t}=\psi (x,t)_0^xF(y,xy)𝑑y+2_x^{\mathrm{}}𝑑y\psi (y,t)F(x,yx)+\frac{}{x}(m(x)\psi (x,t)),$$

where $`F(x,y)`$ is the breakup kernel describing the rules and the rate at which a particle of size $`(x+y)`$ breaks into sizes $`x`$ and $`y`$. Equation $`(1)`$ describes a process whereby cuts are equivalent to seeds being sown on the fragmenting objects, thus producing two new segments. This immediately creates two more new ends belonging to the two different, newly created fragments; in doing so, fragments start losing their masses immediately from either end until they encounter another seed or become dustlike, thereby stopping the loss of their masses. Equation $`(1)`$ also describes the sequential deposition of point-size particles that grow once deposited successfully and stop growing upon collision with another point or growing particle. In other words, the present model describes nucleation and growth when $`\psi (x,t)`$ describes the gap size distribution of size $`x`$ at time $`t`$ or how space is covered by growing objects. This occurs in a number of natural phenomena, including phase separation, wetting, droplet growth, and growth of breath figures. Recently, a variant of the present model was considered in  in which the deposition of growing rods instead of growing seeds was addressed.

We consider the breakup kernel to be $`F(x,y)=(xy)^\beta (x+y)^{\sigma 1}`$, for which the breakup rate $`a(x)=_0^xF(y,xy)𝑑y=px^{\lambda (\beta ,\sigma )}`$, where $`p=\frac{(\mathrm{\Gamma }(\beta +1))^2}{\mathrm{\Gamma }(2\beta +2)}`$, the homogeneity index $`\lambda (\beta ,\sigma )=2\beta +\sigma `$, and $`\mathrm{}(\beta )>1`$. To comply with the present choice of breakup rate $`a(x)`$, it is essential to consider a similar power law form for $`m(x)`$, hence we choose $`m(x)mx^\gamma `$, with $`m`$ a positive real constant. However, $`a(x)`$ being the quantity describing the rate at which particles are fragmenting, $`x^{\lambda (\beta ,\sigma )}\tau (x)`$ must bear the dimension of time, and this put a strong constraint on the exponent $`\gamma `$. That is, the dimensional consistency requires $`\gamma =\lambda (\beta ,\sigma )+1`$. This dimensional consistency has been ignored in all previous studies ; instead $`\gamma <\lambda (\beta ,\sigma )+1`$ was identified as the recession regime and $`\gamma >\lambda (\beta ,\sigma )+1`$ as the fragmentation regime. A closer look at the rate equation reveals that $`x`$ and $`t`$ are inextricably intertwined via the dimensional consistency, and hence the system becomes stochastic in nature. So, it is obvious that either of the two can be taken to be an independent parameter when the other one is expressible in terms of this. For example, if $`x`$ is chosen to be the independent parameter, then using the fact that the dimension of a physical quantity is always expressed as power monomial, we can define the dimensionless quantity $`\xi =\frac{t}{x^\alpha }`$. It is obvious that the dimension of the governed quantity $`\psi (x,t)`$ as well can be expressed in terms of $`x`$ alone, and hence we can define another dimensionless quantity $`\mathrm{\Pi }=\frac{\psi (x,t)}{x^\theta }x^\theta \psi (x,x^\alpha \xi )F(x,\xi )`$. Since $`\xi `$ and $`\mathrm{\Pi }`$ are dimensionless quantities upon transition from one system of units of measurement to another inside a given class, their numerical values must remain unchanged meaning $`\frac{F}{x}=0`$. This implies that $`\mathrm{\Pi }`$ is independent of $`x`$ and can be completely expressed in terms of $`\xi `$ alone. Thus we can define $`\mathrm{\Pi }=\mathrm{\Phi }(\xi )`$ to enable us to write the spatial scaling ansatz $`\psi (x,t)x^\theta \mathrm{\Phi }(\xi )`$. Had we chosen time $`t`$ to be the independent quantity, a similar argument which would lead us to write the temporal scaling ansatz $`\psi (x,t)t^{\theta z}\varphi (\eta )`$ with $`\eta =\frac{x}{t^z}`$. We know that $`x^{\lambda (\beta ,\sigma )}`$ has the dimension of time, therefore $`t^{\frac{1}{\lambda (\beta ,\sigma )}}\delta (t)`$ must have the dimension of $`x`$. This gives us $`z=\frac{1}{\lambda (\beta ,\sigma )}`$, which is known as the kinetic exponent since $`\delta (t)`$ describes the mean or typical cluster size and $`\alpha =\lambda (\beta ,\sigma )`$. Inserting the temporal scaling ansatz into Eq. (1) and requiring scaling to exist would give the same result for kinetic exponent $`z`$. However, dimensional consideration proved to be very instructive and time saving, yet rich in physics.

Note that the exponent $`\theta `$ takes the value for which $`x^\theta `$ and $`t^{\theta z}`$ have a dimension of $`\psi (x,t)`$, and hence it is called the mass exponent. The existence of scaling or a self-similar solution actually means that we can choose self-similar coordinates $`\frac{\psi }{t^{\theta z}}`$ (or $`\frac{\psi }{x^\theta }`$) and $`x/\delta (t)`$ (or $`t/\tau (x)`$) such that their plots for any initial condition collapse into one single curve. It is very instructive to note that at $`\lambda (\beta ,\sigma )=0`$, the governing parameters $`x`$, $`t`$ and the governed parameter $`\psi `$ all lose their dimensional characters. As a result, the system loses its stochastic nature. This means that one can no longer define self-similar coordinates at $`\lambda (\beta ,\sigma )=0`$, which is conceptually very important and crucial for scaling to exist, and hence scaling at $`\lambda (\beta ,\sigma )0`$ breaks down. Note the inherent properties of the fragmentation process, which are the typical or mean cluster size $`\delta (t)`$, must decrease as time proceeds, and the number density must be an increasing function of time. In order to be so, the system must be governed by some conservation laws. However, this is not true as the homogeneity index $`\lambda (\beta ,\sigma )<0`$, in which case $`\delta (t)`$ becomes an increasing function of time. This simply goes against the principle of kinetics of fragmentation. Combining all this, we argue that the system fails to show scaling not only at $`\lambda (\beta ,\sigma )<0`$, as was reported by Cheng and Redner , but also at $`\lambda (\beta ,\sigma )=0`$. Many authors noticed further anomalous behaviour in this regime, e.g., lack of conservation of mass  and an absence of self-averaging . Due to this anomalous behaviour, this regime was termed a transition to shattering . These heurestic arguments actually imply a possible failure to describe physically meaningful fragmentation process by the breakup kernel, for which the homogeneity index $`\lambda (\beta ,\sigma )0`$ and the shattering is actually an articulate term for this regime. The above discussion is true for Eq. $`(1)`$ even when the mass loss term is absent. Therefore, we shall restrict the rest of the discussion to the regime $`\lambda (\beta ,\sigma )>0`$ with $`\mathrm{}(\beta )>1`$.

We now turn to find the mass exponent $`\theta `$, which can only be found if the system follows some conservation laws. For example, for pure fragmentation ($`m=0`$), the mass or size of the system is a conserved quantity and gives $`\theta =2`$. Defining the $`n^{th}`$ moment $`M_n(t)=_0^{\mathrm{}}x^n\psi (x,t)𝑑x`$ and combining it with the rate equation $`(1)`$ for the present choice of $`F(x,y)`$ and $`m(x)`$ yields

$$\frac{dM_n(t)}{dt}=[\frac{(\mathrm{\Gamma }(\beta +1))^2}{\mathrm{\Gamma }(2\beta +2)}\frac{2\mathrm{\Gamma }(\beta +1)\mathrm{\Gamma }(n+\beta +1)}{\mathrm{\Gamma }(n+2\beta +2)}+mn]M_{n+2\beta +\sigma }(t).$$
(1)
Note that the number density $`M_0(t)`$ evolves in the same fashion as it would in the absence of the mass loss term. This means that particles keep losing their mass in a continuous manner by some mechanisms that do not alter the number density. The interesting feature of the above equation is that for $`m>0`$, there are infinitely many $`n=D_f(\beta ,m)`$ values for which $`M_{D_f(\beta ,m)}(t)`$s are conserved quantities. However, for $`m=0`$, there is only one conserved quantity $`M_1(t)`$, i.e., size or mass of the system, and this does not depend on $`\beta `$. The dynamics of the system is governed by conservation laws and, as a consequence, the system shows scale invariance. These conserved quantities in fact are the intrinsic agent responsible for tuning the numerical value of the mass exponent $`\theta `$, and introduce universality to the process. However, at $`\lambda 0`$ there is no evidence that there exists any conservation law for which the mean cluster size can be a decreasing function of time while the mean number density can behave in the opposite way. We can find the numerical value of $`D_f(\beta ,m)`$ by searching for the positive and real root of the equation obtained by setting the term in the bracket of the Eq. (2) equal to zero, which is polynomial in $`n`$ of a degree determined by the $`\beta `$ value. Substituting the temporal scaling ansatz into the definition of $`M_n(t)`$ gives $`M_n(t)t^{(n(\theta 1))z}_0^{\mathrm{}}\eta ^n\varphi (\eta )𝑑\eta `$, and demanding $`M_{D_f}(\beta ,m)`$ be a conserved quantity immediately gives $`\theta =(1+D_f(\beta ,m))`$, which clearly depends on $`\beta `$ and $`m`$ if and only if $`m>0`$. Owing to the random nature of the process and due to the presence of mass loss term, it is clear that when the process continues ad infinitum, it creates a distribution of points (dust) along a line at an extreme late stage that is different from any known set . To measure the size of the set created in the long-time limit, we define a line segment $`\delta (t)=\frac{M_1(t)}{M_0(t)}t^{\frac{1}{\lambda (\beta ,\sigma )}}`$, which is typical cluster size. We can count the number of such segments needed to cover the set, and in the limit $`\delta 0`$ (i.e., $`t\mathrm{}`$), the number $`N(\delta )`$ will simply measure the set and appear to scale as $`N(\delta )\delta ^{D_f(\beta ,m)}`$. The exponent $`D_f(\beta ,m)`$ is known as the Hausdorff-Besicovitch  dimension or the fractal dimension of the arising pattern.

To get a physical picture of the role played by $`m`$, we set $`\beta =0`$ for the time being, for which the polynomial equation becomes quadratic in $`n`$ and the real positive root is $`D_f(m)=\frac{1}{2}(1+1/m)+\frac{1}{2}\sqrt{(1+1/m)^2+4/m}`$ when the second root is $`D=(D_f(m)+1+1/m)`$. Therefore, the exponent $`\theta `$ is also function of $`m`$. The expression for $`D_f(m)`$ reveals that as value of $`m`$ increases, the fractal dimension decreases very sharply and in the limit $`m\mathrm{}`$, $`D_f(m)0`$. This means that as $`m`$ increases, the size of the corresponding arising set decreases sharply due to fast disappearance of its member, whereas, as $`m0`$, $`D_f(m)1`$, that is, we recover the full set (pure fragmentation) that describes a line. On the other hand, had we kept $`m`$ fixed and let $`p`$ decrease, the effect would have been the same as we observed for increasing $`m`$ with $`p=1`$ (i.e. $`\beta =0`$). Thus, it is the ratio between $`m`$ and $`p`$ that matters rather rather than their individual increases or decreases. To give further physical picture of what these results mean, we define mass length relation for the object as $`M_0\delta ^{D_f(m)}`$ and $`M_e\delta ^d`$ for the space where the object is being embedded, where $`d`$ describes the Euclidean space. The density of the property of the object $`\rho `$ then scales as $`\rho \delta ^{D_f(m)d}`$. It is thus clear that for a given class of set created by a specific rule, when $`D_f(m)`$ decreases it means that it is increasingly moving away from $`d`$ and hence more and more members are removed from the full set. This in turn creates increasingly ramified or stringy objects, since $`D_f(m)=d`$ describes the compact object with uniform density. Therefore, this shows that increasing $`m/p`$ ratio means that the mass loss process becomes stronger than the fragmentation process and vice versa.

We now attempt to find the spatial scaling solution for $`\mathrm{\Phi }(\xi )`$. Note that the dimension of the arising pattern is independent of $`\sigma `$ and consequently independent of how fast or slow the system performs the process. Therefore, we can set $`\sigma =1`$ without risk of missing any physics, but it certainly simplifies our calculation. Substituting the spatial scaling ansatz into the rate equation (1) for $`F(x,y)=1`$ and $`m(x)=mx^2`$ and differentiating it with respect to $`\xi `$, transforms the partial integro-differential equation into an ordinary differential equation,

$$\xi (1m\xi )\mathrm{\Phi }^{\prime \prime }(\xi )+\{(1\theta )\xi [2m(2\theta )1]\}\mathrm{\Phi }^{}(\xi )[m(2\theta )(1\theta )(3\theta )]\mathrm{\Phi }(\xi )=0.$$
(2)
For $`m=1`$, this is hypergeometric differential equation  whose only physically acceptable linearly independent solutions are $`{}_{2}{}^{}F_{1}^{}(1,(1+2D_f);D_f;\xi )`$ and $`\xi ^{(1+D_f)}{}_{2}{}^{}F_{1}^{}(2+D_f,D_f;2+D_f;\xi )`$, where $`D_f=0.414213`$. From these exact solutions for spatial scaling function, we can obtain the asymptotic temporal scaling function $`\varphi (\xi )e^{D_f\xi }`$ which satisfies the condition $`\varphi (\xi )0`$ as $`\xi \mathrm{}`$.

We now attempt to see the role of $`\beta `$ on the system. To judge its role, it is clear from the previous discussion that we ought to give equal weight to all the terms in Eq. (1) so that each of them can compete on an equal footing. This can be done if we set $`m=p=\frac{(\mathrm{\Gamma }(\beta +1))^2}{\mathrm{\Gamma }(2\beta +2)}`$ so that the relative strength between fragmentation and the mass loss process stays the same as the value of $`\beta `$ increases. This is a very crucial point to be emphasized. We can obtain the fractal dimension for different values of $`\beta `$, which is simply the real positive root of the polynomial equation in $`n`$ of degree $`\beta `$. A detailed survey reveals that the fractal dimension increases monotonically with increasing $`\beta `$. To find the fractal dimension in the limit $`\beta \mathrm{}`$, we can use the Stirling’s approximation in the polynomial equation to obtain $`\mathrm{ln}[n+1]=(1n)\mathrm{ln}[2]`$ when $`n=0.4569997`$ solves this equation. In order to give a physical picture of the role of $`\beta `$ in the limit $`\beta \mathrm{}`$, we consider the following model: $`F(x,y)=(x+y)^\gamma \delta (xy)`$. This model describes that cuts are only allowed to be in the middle in order to produce two fragments of equal size at each time event. This makes $`a(x)=\frac{1}{2}x^\gamma `$, so we need to choose $`m(x)=\frac{1}{2}x^{\gamma +1}`$, ($`m=\frac{1}{2}`$ gives the same weight as for the fragmentation process). Then the rate equation for $`M_n(t)`$ becomes $`M_n^{}(t)=[(n+1)/22^n]M_{n+\gamma }(t)`$. As before, we set the numerical factor on the right-hand side of this equation equal to zero and then take the natural logarithm on both sides to obtain the $`n`$ value for which $`M_n(t)`$ is time independent. In doing so, we arrive at the same functional equation for $`n`$ as we found for $`\beta \mathrm{}`$. This shows that the kernel $`F(x,y)=(xy)^\beta (x+y)^{\sigma 1}`$ behaves exactly in the same fashion as for $`F(x,y)=(x+y)^\gamma \delta (xy)`$. We thus find that in the limit $`\beta \mathrm{}`$, the resulting distribution of points is a set with fractal dimension $`D_f=0.4569997`$, which is a strictly self-similar fractal as randomness is seized by dividing fragments into equal pieces. We are now in a position to give a physical picture of the role played by $`\beta `$. First of all, the process with $`\beta =0`$ that describes the frequency curve of placing cuts about the size of the fragmenting particles is Poissonian in nature. Consequently, the system enjoys the maximum randomness and the corresponding fractal dimension is $`D_f=0.414213`$. For $`\beta >0`$, the frequency curve of placing cuts about the size of the fragmenting particles is Gaussian in nature, meaning that as the value of $`\beta `$ increases, the particles are increasingly more likely to break in the middle than on either end. That is, as $`\beta `$ increases, the variance decreases in such a manner that in the limit $`\beta \mathrm{}`$, the variance of the frequency curve becomes infinitely narrow, meaning a $`\delta `$-function distribution for which the fragments are broken into two equal pieces. Therefore, there is a spectrum of fractal dimensions between $`\beta 0`$ when $`D_f=0.414213`$ and $`\beta \mathrm{}`$ when $`D_f=0.4569997`$. A detailed numerical survey that we do not present here, confirms that the fractal dimension increases monotonically with $`\beta `$ and reaches a constant value when $`\beta \mathrm{}`$. The previous density-dimension relation implies that increasing $`\beta `$ vis-a-vis increasing order also means that the system loses less and less mass, and this happens despite the fact that now the $`\frac{m}{p}`$ ratio stays the same. This shows that there exists an interplay between fragmentation and the mass loss process that can be tuned either by changing the ratio of $`m`$ and $`p`$, which is obvious of course, or by changing the degree of order alone, which is indeed a nontrivial result.

In summary, we have identified a new set with a wide range of subsets produced by tuning the degree of randomness only. The process starts with an initiator of unit interval $`[0.1]`$, and the generator divides the interval into two pieces, and deleting some parts from either side of both pieces at each time step. The amount of parts to be deleted is determined by the parameter that controls the intensity of randomness. We quantified the size of the resulting set obtained in this way by fractal dimension and showed that the fractal dimension increases with increasing order and reaches its maximum value when the pattern described by the set is perfectly ordered, which is contrary to some recently found results . We also discussed the scaling theory of the process, emphasizing dimensional analysis, and we showed that the shattering is in fact is an articulate term whereby the equation fails to describe physically meaningful fragmentation process. We have also shown that the interplay between fragmentation and mass-loss arises not only from the ratio of their strengths determined by their respective numerical coefficients, but also by the degree of order. We give exact numerical value of the mass exponent, which has never been reported, and we obtained the explicit scaling function for special case of interest. Finally, we argue on the basis of our findings that fractal dimension, degree of order and the extent of ramifications of the arising pattern are interconnected.

The author is grateful to R. M. Ziff for sending valuable comments and acknowledges inspiring correspondence with P. L. Krapivsky. MKH acknowledges the Alexander von Humboldt Foundation for granting the fellowship


text:
# Star Formation Modes in Low-Mass Disk Galaxies

## 1. Introduction

Low-mass spiral galaxies exist within the luminosity range 0.1$`L^{}`$-0.01$`L^{}`$, between the standard giant spirals and the typical dwarf irregulars. These objects have morphological types of Sc-Sdm, are reasonably organized, sometimes have thin stellar disks, display modest spiral structure, and often contain large amounts of HI gas (e.g., Matthews & Gallagher 1997). Under normal conditions their star formation rates tend to be low, but their substantial gas contents make them excellent hosts for starbursts.

In this paper we first consider quiescent star formation in low-mass spiral galaxies as observed in nearby, dynamically cool ‘superthin’ disk galaxies. We then turn to properties of high surface brightness, blue starbursts, which frequently occur in small disk systems, and may be an important contributor to the large populations of faint blue galaxies found at moderate redshifts.

## 2. Star Formation in Slowly Evolving Galactic Disks

### 2.1. Basic Characteristics

The low mass Sd spirals include excellent examples of pure disk galaxies; galaxies which show no evidence for a stellar spheroid. Observationally we most readily identify such galaxies when they are seen near edge-on, allowing us to judge the thickness of the stellar disk. Some low-mass spirals have thin stellar disks–i.e, they are superthin galaxies with little or no stellar halo, and large ratios of radial to vertical disk scale lengths (Gallagher & Hudson 1976, Goad & Roberts 1981, Matthews et al. 1999, Dalcanton & Bernstein 2000, Matthews 2000). These objects are also rich in HI gas, and may have low metal abundances (Bergvall & Rönnback 1995, Giovanelli et al. 1997, Matthews & van Driel 2000). Figure 1 compares WIYN Telescope R-band and H$`\alpha `$ images of UGC 7321, an outstanding example of a superthin disk.

This class of galaxy is of special interest since they would appear to have been assembled from objects containing only gas, and therefore experienced only minimal star formation before their their disks were fully formed. Such galaxies are difficult to make in standard cold dark matter (CDM) scenarios (Gnedin et al. 2000, Wyse 2000). The presence of a thin stellar disk further indicates that little dynamical heating occurred over most of the lifetime of the disk, while the presence of large gas content demonstrates that in a global sense star formation has been inefficient (Matthews 2000). Low-mass Sd spirals are thus the least evolved examples of well-organized disk galaxies in the nearby universe.

### 2.2. Star Formation Patterns in UGC 7321

Figure 1 shows the overall pattern of star formation in UGC 7321 as traced by its HII regions. These extend over almost the entire optically visible disk, but are most frequent and brightest in the central regions of the galaxy. Further insight into products of recent star formation and the state of the ISM come from high angular resolution HST WFPC2 images. As seen in Figure 2, the central region of UGC 7321 looks like a normal galactic disk. While there is no continuous dust lane, dusty dark clouds are obvious and suggestive of the presence of molecular gas, a view that is supported by the detection of molecular CO 1-0 emission by Matthews & Gao (2001). The HST images also reveal distinct clusters and associations of luminous stars in the inner parts of this galaxy.

The outer parts of the UGC 7321 exhibit somewhat different characteristics (Figure 2). Dark nebulae are larger and have lower contrast, while the distribution of luminous stars becomes more diffuse. Possibly two modes of star formation exist within the same galaxy. At radii of $`<`$2 kpc, star formation appears to occur in a normal cycle of well-defined clouds embedded within a stellar disk. In this region we could be seeing self-regulated star formation, where energy supplied by massive stars limits the overall density of the gaseous disk. Further out, where the scale height of the gaseous component becomes comparable to that of the stars, we are probably observing an internally supported gas disk (e.g., by turbulence; Elmegreen & Parravano 1994), where star formation could be playing less of a role in the vertical support of the ISM. The slow evolution of the outer parts of this galaxy then is a natural result of inefficient star formation in a diffuse gaseous disk (see Gallagher et al. 2001 for details).

Gallagher et al. (2001) also estimate a star formation rate from the number of dark clouds to be $`\psi =`$10$`{}_{}{}^{1.5\pm 0.5}M_{}^{}\mathrm{yr}^1`$. For a stellar mass of 2 $`\times `$10$`{}_{}{}^{9}M_{}^{}`$ the Roberts time to form the existing mass of stars is roughly a Hubble time. UGC 7321 is most likely an example of a slowly evolving galaxy rather than a young galaxy, a model that also is supported by the complex vertical structure of its stellar disk (Matthews 2000).

## 3. Interactions and Starbursts in Small Disks

### 3.1. The Basics

In the transition region between giant spirals and star-forming dwarfs we find galaxies with a variety of structures, including the Magellanic irregulars with their thick stellar disks and frequently asymmetric structures, as well as dwarf S0 systems where star disk star formation appears to be winding down. Yet in other ways, such as peak rotation velocities, these galaxies appear to be similar types of objects. A natural suggestion then is that interactions join with initial conditions in determining the structures of the current populations of intermediate luminosity disk (and possibly dwarf E; Mayer et al. 2001) galaxies.

To test this hypothesis we consider observational signatures of interactions and of any possibly associated starbursts:

$`\mathrm{}`$ Interactions between galaxies directly perturb the outer regions of the systems which collide, and may yield wide-spread changes in cases where mergers occur. After-effects may include an increase in disk vertical scale height, shell or ripple structures, or tidal tails (Schweizer & Seitzer 1988, Howard et al. 1993, Velázquez & White 1999). Even a distant collision can suffice to drive gas inwards thereby triggering a starburst (e.g., Mihos & Hernquist 1994, Barton et al. 2000).

$`\mathrm{}`$ Starbursts occur when there is a major (factor of $``$10) and unsustainable increase in a galaxy’s star formation rate. Ongoing starbursts are obvious from the signatures of their huge populations of high mass stars, including strong emission lines, large scale gas outflows, and high power outputs per unit baryonic mass. Products of starbursts include compact super star clusters which may be long-lived.

### 3.2. Starbursts in Small Disk Galaxies

In the nearby Universe, small galaxies host high surface brightness, compact starbursts. That these starbursts reside in disk galaxies can be established from the kinematics and presence of spirals arms and other kinematic features associated with dynamically cool stellar populations (e.g., NGC 3310: Mulder et al. 1995; NGC 7673: Homeier & Gallagher 1999). This class of starburst also includes the nearby archetype M82, where due to its highly inclined orientation the disky structure of the system is obvious while much of the starburst is hidden from optical view.

Preliminary results from a small survey of nearby compact starburst galaxies include evidence for starburst triggering through interactions; e.g., structural or kinematic features that suggest a minor merger or recent interaction or the presence of obviously disturbed companions. In these systems it often appears that a relatively minor perturbation has led to large scale starbursts. Furthermore, the starbursts sometimes persist well after the interaction is past; e.g., in NGC 7673.

A complete picture of the conditions leading to compact starbursts is not yet in hand. However, it has been known for some time that these systems frequently contain large amounts of HI (Smoker et al. 2000); we suspect their gas-rich nature is a key factor in their susceptibility to starbursts. Our model is that compact starbursts take place in gassy disks that were on the margins of stability before they were perturbed (e.g., galaxies like UGC 7321). The external perturbation then drives gas inwards, either directly or via the production of bars or spiral arms, thereby fueling the starburst (e.g., Noguchi 1988).

The presence of large amounts of gas also provides a basis for understanding the large scale clumpy patterns of star formation which are frequently seen in this class of starburst. As emphasized by Elmegreen & Efremov (1997) and Noguchi (1999), a disturbed, gas-rich disk will be unstable and tend to break up into large bound clumps, which we may see as huge regions of star formation. Our high angular resolution observations with HST of NGC 7673, for example, reveal that the clumps often are associations of dense star clusters embedded in backgrounds consisting of luminous massive stars (see also de Grijs et al. 2001). Evidently in these circumstances the star forming hierarchy has gone one level beyond that found in OB associations with dense star clusters paralleling the role of the most massive stars under more normal circumstances.

### 3.3. The Initial Mass Function

In most astrophysical environments where the stellar mass function can be constrained, it appears to be relatively normal. A similar initial mass function (IMF) appears to hold from low density dwarf spheroidal galaxies to dense star clusters, such as R136a. However, as the locations and conditions of star formation become more extreme, might the IMF change?

One well studied case is that of M82, where some data suggested a deficiency of low mass stars within the starburst, while later work implied a relatively normal IMF. In the latest round of this saga, Smith and Gallagher (2001) find the dense super star cluster M82-F appears to be $``$7 times more luminous than predicted for a model of its age with a Salpter IMF. Evidently this object currently lacks the usual population of low mass stars, a situation that also seems to occur in the Galactic center Arches star cluster.

Even if IMF peculiarities in starbursts were limited to only some super star clusters, they could have important observational consequences. For example, the clusters would then be over-luminous for their mass and boost the intensity of a starburst for a fixed astration rate. This would also imply short lifetimes for the super star clusters as bound objects. As a result, much of the evidence of a starburst could dissolve in only a few Gyr.

### 3.4. Starbursts, Tully-Fisher Relation, & Global Properties

M82 provides a useful example of the relationship between starbursts and global properties of galaxies. Figure 3 shows an observed Tully-Fisher (TF) diagram, including the approximate location of M82. We have made no correction for internal absorption and so the K-band luminosity of M82 is a lower limit. As expected, M82 lies well above the standard Tully-Fisher correlation, but as it is involved in a short-lived starburst, galaxies in this state should be rare.

It is then interesting to ask how M82 will evolve on this diagram? Note that if M82 were on the correlation before its starburst and its velocity width does not change duringthe starburst, then the postburst galaxy, that must be more luminous than the pre-burst system, will lie above the TF relationship mean. This presents a problem; post-burst galaxies should be neither obviously peculiar (and thus excluded from normal galaxy TF samples) nor particularly rare. Yet there is little evidence for an asymmetric upwards scatter in the TF diagram that would result from postburst versions of M82.

How can this be? The problem mainly occurs if the velocity width for a galaxy like M82 does not change during a starburst cycle. For example, if starbursts increase both the luminosity and rotation velocity, then galaxies would move up roughly along the mean TF relationship after a starburst. Large deviations would occur only when the stellar luminosity was high, near the peak of a burst, as in M82. This process, shown schematically in Figure 3, is similar to the requirements set by theoretical models of galaxy formation (cf. Steinmetz & Navarro 1999, van den Bosch 2000). We should be able to test this type of model in the rich populations of small disk galaxies in the Local Supercluster, which include objects covering the full range of starburst evolutionary cycles.

## 4. Conclusions

Low-mass spiral galaxies with similar global properties, such as rotation speeds, sizes, or HI masses, can be found with star formation rates extending from $`<`$0.1 to $`>`$1 M$`_{}`$ yr<sup>-1</sup>. The range in levels of star formation activity within these galaxies covers an even wider span. We observe examples of minimal star formation rates in the outer disks of some superthin Sd systems, where the gas disks probably are largely self-supported, and normal, feedback-dominated star formation in their central regions.

At the other extreme, low-mass disks tend to be gas-rich and thus are prime targets for starbursts. When in this runaway mode of making new stars, large star forming ‘clumps’ are common features, possibly resulting from the presence of huge, gravitationally bound gas complexes within their disks. The star-forming clumps are prime locations for the formation of super star clusters, which in turn could have unusual properties, such as top-heavy IMFs.

Most starbursts in low-mass disks are associated with some form of collisional perturbation, which is likely to modify the mass distribution within the starbursting system. Some form of dissipation to concentrate mass within the optical galaxy seems to be essential to preserve the Tully-Fisher relationship for postburst galaxies. Interactions may also play a role in converting low-mass disk galaxies from slowly evolving systems with superthin disks into more ragged and more active forms of late-type galaxies, such as the Magellanic irregulars.

#### Acknowledgments.

We thank NASA for support of much of this work through funding associated with various Hubble Space Telescope projects. JSG also acknowledges partial funding for this research by the National Science Foundation grant AST-9803018 to the University of Wisconsin and by the Vilas Trustees through the University of Wisconsin Graduate School. LDM thanks the National Radio Astronomy Observatory for her Jansky Postdoctoral Fellowship.

References

Barton, E. J., Geller, M. J., & Kenyon, S. J. 2000, ApJ, 530, 660

Bergvall, N. & Rönnback, J. 1995, MNRAS, 273, 603

Dalcanton, J. J. & Bernstein, R. A. 2000, AJ, 120,203

Elmegreen, B. G. & Parravano, A. 1994, ApJ, 435, L121

Elmegreen, B. G. & Efremov, Y. N. 1997, ApJ, 480, 235

Gallagher, J. S. & Hudson, H. S. 1976, ApJ, 209, 389

Gallagher, J. S. et al. 2001, in preparation

Giovanelli, R., Avera, E., & Karaschentsev, I. D. 1997, AJ, 114, 122

Gnedin, N., Norman, M. L., & Ostriker, J. P. 2000, ApJ, 540, 32

Goad, J. W. & Roberts, M. S. 1981, ApJ, 250, 79

de Grijs, R., O’Connell, R. W., & Gallagher, J. S. 2001, AJ, 121, 768

Homeier, N. L. & Gallagher, J. S. 1999, ApJ, 522, 199

Howard, S., Keel, W. C., Byrd, G., & Burkey, J. 1993, ApJ, 417, 512

Ichikawak, T., Yanagisawa, K., Itoh, N., Tarusawa, K., van Driel, W., & Ueno, M. 1995, AJ, 109, 2038

Matthews, L. D. 2000, AJ, 120, 1764

Matthews, L. D. & Gallagher, J. S. 1997, AJ, 114, 1899

Matthews, L. D., Gallagher, J. S., & van driel, W. 2000, AJ, 118, 2751

Matthews, L. D. & van Driel, W. 2000, A&AS, 143, 421

Matthews, L. D. & Gao, Y. 2001, ApJ, 549, L191

Mayer, L. et al. 2001, ApJ, 547, L123

Mihos, J. C. & Hernquist, L. 1994, ApJ, 425, L13

Mulder, P. S., van Driel, W., & Braine, J. 1995, A&A, 300, 687

Noguchi, M. 1988, A&A, 203, 259

Noguchi, M. 1999, ApJ, 514, 77

Schweizer, F. & Seitzer, P. 1988, ApJ, 328, 88

Smith, L. J. & Gallagher, J. S. 2001, MNRAS, in press

Smoker, J. V., Davies, R. D., Axon, D. J., & Hummel, E. 2000, A&A, 361, 19

Steinmetz, M. & Navarro, J. F. 1999, ApJ, 513, 555

van den Bosch, F. C. 1999, ApJ, 530, 177

Velázquez, H. & White, S. D. M. 1999, MNRAS, 304, 254

Verheijen, M. A. W. 1997, The Ursae Major Cluster of Galaxies, Ph.D. thesis, Groningen University

Wyse, R. F. G. 2000, in “The Galactic Halo: from Globular Clusters to Field Stars”, in press.