text:
# Systems of particles with interaction and the cluster formation in condensed matter

## 1  Introduction

Basically, statistical description of many-particle systems concerns homogeneous states. In monograph , a general method based on quantum field theory was successfully applied to homogeneous states of condensed media. An approach that took into account spatial nonhomogeneous states of a system of particles was first proposed in paper ; the approach was resting on the employment of the methodology of quantum field theory as well. The most general approach to the description of a system of particles, which included the interparticle interaction was offered in paper  in which the procedure of discount of the number of states providing for the nonhomogeneous particle distribution was proposed. However, a correct solution of nonlinear equations in the saddle-point is an impossible task for the most part: the general approach does not permit one obtains analytical solutions for a system of particles when the inverse operator of the interaction energy does not exist. In the interim we need to know a series of parameters which characterize possible clusters. The major parameters are the cluster size and the critical temperature of cluster stability. Because of that, one can arise the question about a more simplified approach, which, however, would allow the study of cluster formation in an initially homogeneous system of particles.

In the present work we study a system of interacting particles having regard to the concept proposed in papers . Such an approach reduces the number of variables, which determine the free energy and introduces a peculiar ”combined variable”  playing a role of a typical variable that describes the thermodynamics of the system with a spatial nonhomogeneous particle order.

For instance, as an analog of a certain combined variable one can consider the interaction of two waves on the water surface, which interfere in a limited range. However, after the interference the waves depart and hence in this case the combined variable does not bring into existence any instability of the system studied. On the other hand, the interaction of an incident electromagnetic wave with a polar crystal results into the polariton, a stable formation, which is indeed described by a combined (canonical) variable – a mixture of the electromagnetic field and optical phonon.

Can the combined variable give rise to the change in the homogeneous distribution of interacting particles? We analyze such an option considering concrete examples. Below we take a look at the specific systems with long-range repulsion/attraction and short-range attraction/repulsion. In particular, the approach provides a way of estimating of the cluster size, the number of particles in a cluster and the temperature of phase transition to the cluster state. The approach is examined in three examples: i) electrons on the liquid helium surface, ii) particles with the shielding Coulomb interaction, which are found in an applied elastic potential, and iii) gravitating masses with the Hubble expansion.

## 2  Statistical mechanics for model systems with interaction. The principle of states selection

Let us consider a model system of particles with arbitrary kinds of interactions. Nevertheless, the de Broglie thermal wavelength $`\lambda `$ of a particle is taken to be greater than the mean distant $`l`$ between particles. However on the other hand, the value of $`\lambda `$ must be smaller than the mean amplitude of scattering. The type of the statistics introduces significant peculiarities to the interparticle interaction, which, however, can be rendered by the classical methods allowing one does away with dynamic quantum correlations. In this case the macroscopic state of the system in question may be specified by filling numbers $`n_s`$, where the symbol $`s`$ describes the state of the system. In other words, particles can occupy only knots are still non-occupied. This means that the filling numbers $`n_s`$ can run only two meanings: 1 (the $`s`$th knot is occupied) or 0 (the $`s`$th knot is not occupied). Besides, the allocation of particles in the lattice invites the type of statistics, which set limits on the behaviour of particles, i.e. the statistics refers particles to either bosons or fermions, see below.

The Hamiltonian of the system of interacting particles can be written in the form

$$(n)=\underset{s}{}E_sn_s\frac{1}{2}\underset{ss^{}}{}v_{ss^{}}n_sn_s^{}+\frac{1}{2}\underset{ss^{}}{}u_{ss^{}}n_sn_s^{}$$
(1)
where $`E_s`$ is the additive part of the particle energy in the $`s`$th state; $`v_{ss^{}}`$ and $`u_{ss^{}}`$ are, respectively, the paired energy of attraction and repulsion between particles in the states $`s`$ and $`s^{}`$. Signs before the potentials in expression (1) reflect the proper signs of the attractive and repulsive paired energies and therefore both of the functions $`v_{ss^{}}`$ and $`u_{ss^{}}`$ are positive.

Such writing of the Hamiltonian corresponds to the model in which particles are found in sites of the three-dimensional lattice. As will be evident from the further consideration, a specific form of the lattice will be unneeded. However the medium will be supposed isotropic at the transition to the continual approximation.

Note that the Hamiltonian of the model of replacement of solid solutions  amounts to the Hamiltonian (1) as well. This means that the approach proposing can describe the behaviour of defects in solids, allows, liquids, etc. (in these cases one particle will characterize an unfilled place and another particle that of a filled one). The reduced solution Hamiltonian can be amounted to a general Ising model, which takes into account an arbitrary interaction allowing the spatial-nonhomogeneous distribution of particles.

The partition function of the system is

$$Z=\underset{\{n_s\}}{}\mathrm{exp}((n_s)/k_\mathrm{B}T)$$
(2)
where the summation provides for all possible values $`n_s`$, i.e., all states of the system. Such summation can be formally performed, if one passes to field variables according to Ref. . It is known from the theory of Gauss integrals that

$$\mathrm{exp}[\frac{\nu ^2}{2}\underset{ss^{}}{}w_{ss^{}}n_sn_s^{}]=\mathrm{Re}\underset{\mathrm{}}{\overset{\mathrm{}}{}}D\chi \mathrm{exp}[\nu \underset{s}{}n_s\chi _s\frac{1}{2}\underset{ss^{}}{}w_{ss^{}}^1\chi _s\chi _s^{}]$$
(3)
where $`D\chi \mathrm{\Pi }_s\sqrt{\mathrm{det}W_{ss^{}}}\sqrt{2\pi }d\chi _s`$ implies the functional integration with respect to the field $`\chi `$ and $`\nu ^2=\pm 1`$ in relation to the sign of interaction ($`+1`$ for attractive and $`1`$ for repulsion). Besides the obligatory condition

$$\underset{s^{\prime \prime }}{}w_{ss^{\prime \prime }}^1w_{s^{\prime \prime }s^{}}=\delta _{ss^{}}$$
(4)
is full-filled. In particular, condition (4) in principle could allow one to determine the inverse operator $`w_{ss^{}}^1`$ (in the case when it is exist) which in its turn would permit the construction of the Green function for interacting particles.

Let us introduce the dimensionless energies $`\stackrel{~}{w}_{ss^{}}w_{ss^{}}/k_\mathrm{B}T`$ and $`\stackrel{~}{E}_sE_s/k_\mathrm{B}T`$. Note that hereinafter the tilde over a symbol means the division of the symbol by $`k_\mathrm{B}T`$.

So the partition function (2) can be represented in the form

$`Z`$ $`=`$ $`\mathrm{Re}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}D\varphi {\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}D\psi {\displaystyle \underset{\{n_s\}}{}}\mathrm{exp}[{\displaystyle \underset{s}{}}\stackrel{~}{E}_sn_s+{\displaystyle \underset{s}{}}(\psi _s+i\varphi _s)n_s`$ (5)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ss^{}}{}}(\stackrel{~}{u}_{ss^{}}^1\varphi _s\varphi _s^{}+\stackrel{~}{v}_{ss^{}}^1\psi _s\psi _s^{})].`$
We do not include the dependence on momentum of particles to the functional integrals in the partition function (5); the sum takes into consideration only the spatial distribution of particles and their energy.

Now we can settle the quantity of particles in the system, $`_sn_s=N`$. The procedure is an equivalent to the consideration of the canonical ensemble. For this purpose one can use the formula known in the theory of complex variable

$$\frac{1}{2\pi i}𝑑zz^{N1_sn_s}=1.$$
(6)
It makes possible to introduce the sum of the canonical ensemble

$`Z_N`$ $`=`$ $`\mathrm{Re}{\displaystyle \frac{1}{2\pi i}}{\displaystyle }dz{\displaystyle }D\varphi {\displaystyle }D\psi \mathrm{exp}\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{ss^{}}{}}(\stackrel{~}{u}_{ss^{}}^1\varphi _s\varphi _s^{}+\stackrel{~}{v}_{ss^{}}^1\psi _s\psi _s^{})`$ (7)
$`+(N1)\mathrm{ln}z\}\times {\displaystyle }_{\{n_s\}}\mathrm{exp}\{{\displaystyle \underset{s}{}}n_s(\psi _s+i\varphi _s\stackrel{~}{E}_s\mathrm{ln}z)\}.`$
Summing over $`n_s`$ (note $`\{n_s\}=0,1`$) we get

$$Z=\mathrm{Re}\frac{1}{2\pi i}D\varphi D\psi 𝑑ze^{S(\varphi ,\psi ,z)}$$
(8)
where

$`S`$ $`=`$ $`{\displaystyle \underset{s}{}}\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{s^{}}{}}(\stackrel{~}{u}_{ss^{}}^1\varphi _s\varphi _s^{}+\stackrel{~}{v}_{ss^{}}^1\psi _s\psi _s^{})+\eta \mathrm{ln}\left|1+{\displaystyle \frac{\eta }{z}}e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{cos}\varphi _s\right|\}`$ (9)
$`+(N1)\mathrm{ln}z.`$
Here, the function $`\eta `$ is equal to $`\pm 1`$ (Fermi or Bose statistics), see e.g. Refs. \[7-9\]. This provides the means determining the necessary procedure of the states selection that introduces the most essential contribution to the partition function and defines the free energy of the system. Let us set $`z=\xi +i\zeta `$ and consider the action $`S`$ on the transit path passing through the saddle-point with a fixed imaginable variable $`\mathrm{Im}z=\zeta _0`$. In this case, it stands to reason that the action $`S`$, similarly to quantum field theory, is the variational functional that depends on three variables: the fields $`\varphi _s`$ and $`\psi _s`$, and the fugacity $`\xi `$ (here $`\xi =e^{\mu /k_\mathrm{B}T}`$ where $`\mu `$ is the chemical potential). The extremum of the functional must be realized at solutions of the equations $`\delta S/\delta \varphi _s=0`$, $`\delta S/\delta \psi _s=0`$, and $`\delta S/\delta \xi =0`$. These equations appear as the following

$$\underset{s^{}}{}\stackrel{~}{u}_{ss^{}}^1\varphi _s^{}=\frac{e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{sin}\varphi _s}{\xi +\eta e^{\stackrel{~}{E}}e^{\psi _s}\mathrm{cos}\varphi _s};$$
(10)
$$\underset{s^{}}{}\stackrel{~}{v}_{ss^{}}^1\psi _s^{}=\frac{e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{cos}\varphi _s}{\xi +\eta e^{\stackrel{~}{E}}e^{\psi _s}\mathrm{cos}\varphi _s};$$
(11)
$$\underset{s^{}}{}\frac{e^{\stackrel{~}{E}_s^{}}e^{\psi _s^{}}\mathrm{cos}\varphi _s^{}}{\xi +\eta e^{\stackrel{~}{E}_s^{}}e^{\psi _s^{}}\mathrm{cos}\varphi _s^{}}=N1.$$
(12)
Equations from (10) to (12) completely solve the problem of the statistical description of systems with any type of interaction. One or another of states of the system is realized in accordance with the solution of the nonlinear equations. Among the solutions there are solutions which correspond to the spatial nonhomogeneous distribution of particles. So the Bose condensation is realized in the real space because we can say that an aggregation of particles is the Bose condensation of a sort. One of the possibilities of cluster formation was demonstrated in Ref. ; the nonhomogeneous distribution of particles which were characterized by the shielded Coulomb interaction for both the attraction and the repulsion just yielded stable clusters. Such type of the interaction was chosen owing to the existence of the inverse operator of the interaction energy in an explicit form. Nevertheless, it is a very difficult problem to find the inverse operator among the most widespread types of interactions. This is why, there is the necessity to develop methods which are capable to describe the particle distribution at any interaction. In Refs.  the approach based on the selection of states in the form of a ”combined variable” has been proposed. The combined variable has given us the chance to find connections between different fields (such as $`\varphi _s`$ and $`\psi _s`$ in the present paper) along the extremal path that passes through the saddle-point.

## 3 Combined variable and the selection of states

Let us introduce the designation

$$\mathrm{\Gamma }_s=\frac{e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{cos}\varphi _s}{\xi +\eta e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{cos}\varphi _s}.$$
(13)
The value $`\mathrm{\Gamma }_s`$ includes the two field variables, $`\psi _s`$ and $`\varphi _s`$, and the fugacity $`\xi `$. Inserting $`\mathrm{\Gamma }_s`$ into the left hand side of Eq. (12) we will see that the total sum $`_s\mathrm{\Gamma }_s`$ is equal to the number of particles in the system studied

$$\underset{s}{}\mathrm{\Gamma }_s=N1.$$
(14)
This means that $`\mathrm{\Gamma }_s`$ directly specifies the quantity of particles in the $`s`$th state. Hence the value $`\mathrm{\Gamma }_s`$ might be considered as a characteristic of the number of particles contained in the $`s`$th cluster. So it is quite possible to consider $`\mathrm{\Gamma }_s`$ as the variable of particle number and, on the other hand, it may be called the combined variable owing to the fact that it combines the variables $`\psi _s`$, $`\varphi _s`$, and $`\xi `$. Probably the introduction of the combined variable is a rough approximation, but it admits to determine the clusterization conditions in cases when the inverse operator of the interaction potential can not be found.

Utilizing the variable (13) we can represent the action (9) as a function of only two variables, $`\mathrm{\Gamma }_s`$ and $`\xi `$, and it makes possible to do away with the inverse operators $`\stackrel{~}{v}_{ss^{}}^1`$ and $`\stackrel{~}{u}_{ss^{}}^1`$ in the final expression.

If we multiply two hand sides of Eq. (11) by the function $`\stackrel{~}{v}_{ss^{\prime \prime }}`$ and then summing Eq. (11) over $`s`$, we obtain

$$\psi _s^{}=\underset{s}{}\stackrel{~}{v}_{ss^{}}\mathrm{\Gamma }_s.$$
(15)
Now let us multiply the same equation (11) by $`\psi _s`$; then summing it over $`s`$ we acquire

$$\underset{ss^{}}{}\stackrel{~}{v}_{ss^{}}^1\psi _s^{}\psi _s=\underset{ss^{}}{}\stackrel{~}{v}_{ss^{}}\mathrm{\Gamma }_s^{}\mathrm{\Gamma }_s.$$
(16)
Multiplying Eq. (10) by $`\stackrel{~}{u}_{s^{\prime \prime }s}`$ and summing it over $`s`$ we have

$$\underset{ss^{}}{}\stackrel{~}{u}_{ss^{}}^1\stackrel{~}{u}_{s^{\prime \prime }s}\varphi _s^{}=\underset{s}{}\stackrel{~}{u}_{ss^{\prime \prime }}\frac{e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{sin}\varphi _s}{\xi +\eta e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{cos}\varphi _s}.$$
(17)
From relationship (13) one obtains

$$\frac{1}{\xi }e^{\stackrel{~}{E}_s}e^{\psi _s}\mathrm{cos}\varphi _s=\frac{\mathrm{\Gamma }_s}{1\eta \mathrm{\Gamma }_s}.$$
(18)
Let us substitute (18) into (17) taking into account Eq. (15) and using condition (4) for the left-hand side of Eq. (17). We gain

$$\varphi _s=\underset{s^{}}{}\stackrel{~}{u}_{s^{}s}\frac{1}{\xi }e^{\stackrel{~}{E}_s}e^{_{s^{\prime \prime }}\stackrel{~}{v}_{s^{\prime \prime }s}\mathrm{\Gamma }_{s^{\prime \prime }}}(1\eta \mathrm{\Gamma }_s^{})\mathrm{sin}\varphi _s^{}.$$
(19)
With the formula $`\mathrm{sin}\varphi _s^{}=\sqrt{1\mathrm{cos}^2\varphi _s^{}}`$, $`\mathrm{cos}\varphi _s`$ can be substituted from (18) into (19). As a result, instead of (19) we have

$$\varphi _s=\underset{s^{}}{}\stackrel{~}{u}_{s^{}s}\sqrt{\frac{1}{\xi ^2}e^{2\stackrel{~}{E}_s^{}}e^{2_{s^{\prime \prime }}\stackrel{~}{v}_{s^{\prime \prime }s}\mathrm{\Gamma }_{s^{\prime \prime }}}(1\eta \mathrm{\Gamma }_s^{})^2\mathrm{\Gamma }_s^{}^2}.$$
(20)
Consequently using Eqs. (10) and (13) we get

$$\underset{ss^{}}{}\stackrel{~}{u}_{ss^{}}^1\varphi _s^{}\varphi _s=\underset{s}{}\mathrm{\Gamma }_s\varphi _s\mathrm{tan}\varphi _s,$$

or in the explicit form

$`{\displaystyle \underset{ss^{}}{}}\stackrel{~}{u}_{ss^{}}^1\varphi _s^{}\varphi _s`$ $`=`$ $`{\displaystyle \underset{ss^{}}{}}\stackrel{~}{u}_{ss^{}}\sqrt{{\displaystyle \frac{1}{\xi ^2}}e^{2\stackrel{~}{E}_s^{}}e^{2_{s^{\prime \prime }}\stackrel{~}{v}_{s^{}s^{\prime \prime }}\mathrm{\Gamma }_{s^{\prime \prime }}}(1\eta \mathrm{\Gamma }_s^{})^2\mathrm{\Gamma }_s^{}^2}`$ (21)
$`\times \sqrt{{\displaystyle \frac{1}{\xi ^2}}e^{2\stackrel{~}{E}_s}e^{2_{s^{\prime \prime }}\stackrel{~}{v}_{ss^{\prime \prime }}\mathrm{\Gamma }_{s^{\prime \prime }}}(1\eta \mathrm{\Gamma }_s)^2\mathrm{\Gamma }_s^2};`$
here, allowance is made for the negative sign of the field $`\varphi _s`$ (see (19)).

Now using relationships (16) and (21) we can rewrite the action (9) in the point of extremum as follows:

$`S=`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ss^{}}{}}\stackrel{~}{v}_{ss^{}}\mathrm{\Gamma }_s^{}\mathrm{\Gamma }_s\eta {\displaystyle \underset{s}{}}\mathrm{ln}|1\eta \mathrm{\Gamma }_s|+(N1)\mathrm{ln}\xi {\displaystyle \frac{1}{2}}{\displaystyle \underset{ss^{}}{}}\stackrel{~}{u}_{ss^{}}`$ (22)
$`\times \sqrt{[{\displaystyle \frac{1}{\xi ^2}}e^{2\stackrel{~}{E}_s^{}}e^{2_{s^{\prime \prime }}\stackrel{~}{v}_{s^{}s^{\prime \prime }}\mathrm{\Gamma }_{s^{\prime \prime }}}(1\eta \mathrm{\Gamma }_s^{})^2\mathrm{\Gamma }_s^{}^2]}`$
$`\times \sqrt{[{\displaystyle \frac{1}{\xi ^2}}e^{2\stackrel{~}{E}_s}e^{2_{s^{\prime \prime }}\stackrel{~}{v}_{ss^{\prime \prime }}\mathrm{\Gamma }_{s^{\prime \prime }}}(1\eta \mathrm{\Gamma }_s)^2\mathrm{\Gamma }_s^2]}.`$
It has been shown in Ref.  that the minimum of the free energy based on expression (22) is lower than that obtained in the framework of the mean field approximation. Indeed, let us assume following Ref.  that $`N_u`$ and $`N_v`$ are numbers of particles getting into the influence of the repulsive and attractive force respectively. The number of particles $`\mathrm{\Gamma }_s`$ in a cluster can be put constant equals $`\mathrm{\Gamma }`$. Let interaction energies in the range of influence of these forces be change to their average values. In this case expression (22) becomes

$`S={\displaystyle \frac{1}{2}}(N1)\left\{N_v\mathrm{\Gamma }\overline{U}\mathrm{\Gamma }+N_v\mathrm{\Gamma }\overline{V}\mathrm{\Gamma }+2\eta \mathrm{ln}|1\eta \mathrm{\Gamma }|\right\}`$ (23)
$`{\displaystyle \frac{1}{2}}(N1)\left[N_u\overline{U}(1\eta \mathrm{\Gamma })e^{2\overline{E}}e^{2N_v\overline{V}\mathrm{\Gamma }^2}\xi ^22\mathrm{ln}\xi \right]`$
where $`\overline{U}`$ and $`\overline{V}`$ are average values of the repulsive and attractive potential respectively. Here, in the mean field approximation the expression in the curly brackets coincides with the free energy. The next term, i.e. the expression in the square brackets decreases the free energy at $`\xi <1`$. This inequality, $`\xi <1`$, corresponds to the coherent state of a system of interacting particles. Thus the approach that is developing can be considered as a general methodology of the mean field theory applying to the nonhomogeneous particle distribution.

Expression (22) is a function of only one variable $`\mathrm{\Gamma }`$ and the fugacity $`\xi `$ and it is applicable to any kind of the interaction, even though the inverse operator is unknown. The most advantage of expression (22) lies in the fact that it provides a way of finding of such characteristics of the system as its size, the number of particles in a cluster and the temperature of phase transition.

Now let us pass to continual variables into the action (22). Inasmuch as we are interested in the nonhomogeneous distribution of particles in an indeterminate volume, let a radius $`R`$ be the fitting parameter of the system studied. Assume that the density of particles is distinguished from zero only in the cluster volume. Then Eq. (14) can be written as

$$\mathrm{\Gamma }K=N1$$
(24)
where $`K`$ is the number of clusters in the system and $`\mathrm{\Gamma }`$ defined in expression (13) is the combined variable of the fields $`\varphi _s`$ and $`\psi _s`$, and the fugacity $`\xi `$ in a cluster. As is seen from expression (24) the variable $`\mathrm{\Gamma }`$ can be interpreted as the mean quantity of particles in a cluster. Then the passage to the continual presentation is realized by the substitution

$$\underset{s}{}f_s=K\frac{1}{V}\underset{\mathrm{cluster}}{}f(\stackrel{}{r})𝑑\stackrel{}{r};$$
(25)
here the integration extends for the volume of a cluster and $`V=\frac{4\pi }{3}g^3`$ is the effective volume occupied by one particle where $`g`$ is the distance between particles, i.e., the lattice constant.

Thus in the continual presentation and with allowing for the assumption (14), we can transform the action (22) to the form

$`S`$ $`=`$ $`{\displaystyle \frac{1}{2V^2}}{\displaystyle 𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}\stackrel{~}{v}(\stackrel{}{r}\stackrel{}{r^{}})\mathrm{\Gamma }(\stackrel{}{r^{}})\mathrm{\Gamma }(\stackrel{}{r})}`$ (26)
$`{\displaystyle \frac{1}{2V^2}}{\displaystyle 𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}\stackrel{~}{u}(\stackrel{}{r}\stackrel{}{r^{}})\mathrm{\Gamma }(\stackrel{}{r^{}})\mathrm{\Gamma }(\stackrel{}{r})}`$
$`\times \sqrt{e^{2\stackrel{~}{E}(\stackrel{}{r})}e^{2{\scriptscriptstyle \stackrel{~}{v}(\stackrel{}{r}\stackrel{}{r^{}})\mathrm{\Gamma }(\stackrel{}{r^{\prime \prime }})𝑑\stackrel{}{r^{\prime \prime }}}}{\displaystyle \frac{1}{\xi ^2}}({\displaystyle \frac{1}{\mathrm{\Gamma }(\stackrel{}{r})}}\eta )^21}`$
$`\times \sqrt{e^{2\stackrel{~}{E}(\stackrel{}{r})}e^{2{\scriptscriptstyle \stackrel{~}{v}(\stackrel{}{r}\stackrel{}{r^{\prime \prime }})\mathrm{\Gamma }(\stackrel{}{r^{\prime \prime }})𝑑\stackrel{}{r^{\prime \prime }}}}{\displaystyle \frac{1}{\xi ^2}}({\displaystyle \frac{1}{\mathrm{\Gamma }(\stackrel{}{r})}}\eta )^21}`$
$`+\eta {\displaystyle \frac{1}{V}}{\displaystyle \mathrm{ln}|1\eta \mathrm{\Gamma }(\stackrel{}{r})|d\stackrel{}{r}}+(N1)\mathrm{ln}\xi `$
where the integrals are the same as in Eq. (25).

The integration is effected by the rule

$`{\displaystyle \frac{1}{V}}{\displaystyle \underset{\mathrm{cluster}}{}}f(\stackrel{}{r})𝑑\stackrel{}{r}`$ $`=`$ $`{\displaystyle \frac{1}{\frac{4\pi }{3}g^3}}{\displaystyle \underset{0}{\overset{2\pi }{}}}𝑑\phi {\displaystyle \underset{0}{\overset{\pi }{}}}\mathrm{sin}\theta d\theta {\displaystyle \underset{g}{\overset{R}{}}}f(r)r^2𝑑r`$ (27)
$`=`$ $`{\displaystyle \frac{1}{\frac{4\pi }{3}g^3}}4\pi {\displaystyle \underset{g}{\overset{R}{}}}f(r)r^2𝑑r;`$
We shift here the limits of integration from 0 and $`R`$ to $`g`$ and $`R`$ respectively as it allows one to eliminate the singularity of the integrand (the same procedure is made below in expression (29)). Thus, we can now normalize the integrals to the number of particles $`\mathrm{\Gamma }`$ in a cluster, that is,

$$\frac{1}{V}𝑑\stackrel{}{r}=\frac{1}{(4\pi /3)g^3}4\pi \underset{g}{\overset{R}{}}r^2𝑑r=\frac{R^3g^3}{g^3}=\mathrm{\Gamma }1\mathrm{\Gamma };$$
(28)
the last approximation is in fact correct as we assume that $`\mathrm{\Gamma }`$ satisfies the inequality $`\mathrm{\Gamma }1`$.

Relation (28) allow us to introduce the dimensionless variable $`x=r/g`$ in integral (27). Thereby the rule of transformation from summation to integration, i.e. (25), becomes

$$\frac{1}{K}\underset{s}{}f_s=3\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}f(gx)x^2𝑑x$$
(29)
(once again, in expressions (27) we shift the limits of integration from 0 and $`R`$ to $`g`$ and $`R`$ respectively as it allows one to eliminate the singularity of the integrand; the same shift is made in expression (29)).

Having integrated the action (26), we should exploit the following relationships

$`{\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}\stackrel{~}{v}(\stackrel{}{r}\stackrel{}{r^{}})\mathrm{\Gamma }(\stackrel{}{r})\mathrm{\Gamma }(\stackrel{}{r^{}})}={\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑\stackrel{}{r}\stackrel{~}{v}(\stackrel{}{r})\mathrm{\Gamma }(\stackrel{}{r})𝑑\stackrel{}{r^{}}\mathrm{\Gamma }(\stackrel{}{r^{}})}`$
$`={\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑\stackrel{}{r}\stackrel{~}{v}(\stackrel{}{r})\mathrm{\Gamma }\mathrm{\Gamma }𝑑\stackrel{}{r^{}}}={\displaystyle \frac{\mathrm{\Gamma }^2}{V}}{\displaystyle 𝑑\stackrel{}{r}\stackrel{~}{v}(\stackrel{}{r})}=3\mathrm{\Gamma }^2{\displaystyle \underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}}𝑑xx^2\stackrel{~}{v}(gx);`$ (30)
$$\frac{1}{V^2}𝑑\stackrel{}{r}𝑑\stackrel{}{r^{\prime \prime }}\stackrel{~}{u}(\stackrel{}{r}\stackrel{}{r^{\prime \prime }})\mathrm{\Gamma }(\stackrel{}{r})\mathrm{\Gamma }(\stackrel{}{r^{\prime \prime }})=\frac{1}{V^2}𝑑\stackrel{}{r}\stackrel{~}{u}(\stackrel{}{r})\mathrm{\Gamma }^2=3\mathrm{\Gamma }^2\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}𝑑xx^2\stackrel{~}{u}(gx);$$
(31)
With the transformations, we have used the step function: $`\mathrm{\Gamma }(r)=\mathrm{\Gamma }\vartheta (\mathrm{\Gamma }r^3/g^3)`$ where $`\vartheta (\mathrm{\Gamma }r^3/g^3)=0`$ if $`r^3/g^3\mathrm{\Gamma }`$ and $`\vartheta (\mathrm{\Gamma }r^3/g^3)=1`$ if $`r^3/g^3<\mathrm{\Gamma }`$. Besides,

$`{\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}\stackrel{~}{u}(\stackrel{}{r},\stackrel{}{r^{}})f(\stackrel{}{r})f(\stackrel{}{r^{}})}`$ $`=`$ $`{\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑\stackrel{}{r}\stackrel{~}{u}(\stackrel{}{r})f(\mathrm{\Gamma })𝑑\stackrel{}{r^{}}f(\mathrm{\Gamma })}`$ (32)
$`=`$ $`3f^2{\displaystyle \underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}}𝑑xx^2\stackrel{~}{u}(gx),`$
where $`f=\sqrt{\mathrm{exp}[2\stackrel{~}{E}+2\mathrm{\Gamma }𝑑\stackrel{}{r}\stackrel{~}{v}(\stackrel{}{r})]}`$.

Now if we introduce the designations

$$a=3\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}𝑑xx^2\stackrel{~}{u}(gx);b=3\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}𝑑xx^2\stackrel{~}{v}(gx),$$
(33)
we will represent the action (26) in the following final form

$$S=K\left\{\frac{1}{2}(ab)\mathrm{\Gamma }^2\frac{1}{2}\frac{1}{\xi ^2}\left(1\eta \mathrm{\Gamma }\right)^2e^{2\stackrel{~}{E}+2b\mathrm{\Gamma }}+\eta \mathrm{ln}|1\eta \mathrm{\Gamma }|\right\}+(N1)\mathrm{ln}\xi .$$
(34)
If we minimize the action (34) by the variable $`\mathrm{\Gamma }`$, we will be able to define the number of particles in a cluster, the cluster size, and the temperature of phase transition that triggers the nonhomogeneous distribution of particles. It must be emphasized that the definition of the parameters of nonhomogeneous formations requires only explicit forms of potentials of interparticle interactions. Then the system itself will select the realization which will provide for the minimum of the free energy. Nonetheless, we assume by first the availability of clusters and then the conditions and parameters of their existence are defined.

Studying the phase transition of clusters to the spatial homogeneous distribution of particles, we should take into account the behaviour of the fugacity $`\xi =e^{\mu /k_\mathrm{B}T}`$. Indeed, in the classical system the chemical potential

$$\mu =k_\mathrm{B}T\mathrm{ln}(\lambda _T^3n)$$
(35)
where $`n`$ is the particle concentration and the de Broglie thermal wavelength of a particle with a mass $`m`$

$$\lambda _T=h/\sqrt{3mk_\mathrm{B}T}.$$
(36)
Thus in the classical case the fugacity

$$\xi =\lambda _T^3n<<1.$$
(37)
In the case of a quantum system the strong inequality does not hold; in this case $`\xi <1`$ and $`\xi `$ is only a slight less than the unit. Thereby when we investigate the action (34) searching for clusters at the temperature $`T<T_c`$, we may neglect the last term as $`|\mathrm{ln}\xi |`$ is rather small in comparison with terms which include highest orders of $`\mathrm{\Gamma }`$. However, if we consider the action (34) scanning for the critical temperature $`T=T_c`$, we should retain the last term because in this case $`\xi `$ trends to a very small magnitude such that $`|\mathrm{ln}\xi |>>1`$ and yet other terms decrease owing to the diminution of the value of $`\mathrm{\Gamma }`$.

## 4 Clusters in condensed media

In this section we apply the methodology developed above to systems of particles with long-range repulsion (attraction) and short range attraction (repulsion).

### 4.1 Electrons on liquid helium surface

The energy of electrostatic repulsion between electrons can be written as $`u_{ss^{}}=1/4\pi ϵ_0Q^2/(r_sr_s^{})`$ where $`Q`$ is the effective charge of the electron in the helium film. The attraction between electrons caused by the deformation of the helium film can be taken as $`v_{ss^{}}=\frac{1}{2}\gamma _{\mathrm{film}}(r_sr_s^{})^2`$ where the coefficient $`\gamma _{\mathrm{film}}`$ characterizes the elasticity of the helium film in respect to the deformation caused by imbedded electrons (see, e.g. review ).

In this case the parameters $`a`$ and $`b`$ (33) of the action (34) are transformed to the following

$$a=3\frac{Q^2/(4\pi ϵ_0g)}{k_\mathrm{B}T}\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}\frac{1}{x}x^2𝑑x\frac{3}{8\pi ϵ_0}\frac{Q^2}{gk_\mathrm{B}T}\mathrm{\Gamma }^{2/3};$$
(38)
$$b=3\frac{\gamma _{\mathrm{film}}g^2}{2k_\mathrm{B}T}\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}x^2x^2𝑑x\frac{3}{10}\frac{\gamma _{\mathrm{film}}g^2}{k_\mathrm{B}T}\mathrm{\Gamma }^{5/3}$$
(39)
or

$$a=\stackrel{~}{\alpha }\mathrm{\Gamma }^{2/3};b=\stackrel{~}{\beta }\mathrm{\Gamma }^{5/3}$$
(40)
where

$$\stackrel{~}{\alpha }=\frac{3}{8\pi ϵ_0}\frac{Q^2}{gk_\mathrm{B}T};\stackrel{~}{\beta }=\frac{3}{10}\frac{\gamma _{\mathrm{film}}g^2}{k_\mathrm{B}T}$$
(41)
and $`\alpha >>\beta `$ (hereafter the tilde over a symbol means the division by $`k_\mathrm{B}T`$). Since for electrons $`\eta =+1`$ and the inequality $`\stackrel{~}{E}E/k_\mathrm{B}T1`$ is legitimate (i.e., $`\mathrm{exp}(2E/k_\mathrm{B}T)1`$ where $`E`$ is the kinetic energy of an electron), the action (34) can be rewritten as

$$SK\left\{\frac{1}{2}\stackrel{~}{\alpha }\mathrm{\Gamma }^{8/3}\frac{1}{2}\stackrel{~}{\beta }\mathrm{\Gamma }^{11/3}+\mathrm{ln}\mathrm{\Gamma }+\mathrm{\Gamma }\mathrm{ln}\xi \right\}.$$
(42)
The extremum of the free energy is achieved at the solution of the equation $`S^{}(\mathrm{\Gamma })=0`$, or in the explicit form:

$$\frac{4}{3}\stackrel{~}{\alpha }\mathrm{\Gamma }^{5/3}\frac{11}{6}\stackrel{~}{\beta }\mathrm{\Gamma }^{8/3}+\frac{1}{\mathrm{\Gamma }}+\mathrm{ln}\xi =0;$$
(43)
Retaining leading terms in Eq. (43), i.e. highest powers to $`\mathrm{\Gamma }`$, the equation is reduced to

$$\frac{4}{3}\stackrel{~}{\alpha }\mathrm{\Gamma }^{5/3}\frac{11}{6}\stackrel{~}{\beta }\mathrm{\Gamma }^{8/3}0$$
(44)
and hence the solution to it is equal to

$$\mathrm{\Gamma }=\frac{8}{11}\frac{\alpha }{\beta }\frac{40}{33}\frac{nQ^2}{ϵ_0\gamma _{\mathrm{film}}}$$
(45)
where $`n=(\frac{4\pi }{3}g^3)^1`$ is the concentration of electrons in a cluster. The sufficiency is also satisfied, i.e., $`S^{\prime \prime }(\mathrm{\Gamma })>0`$ with the solution (45) and therefore it actually determines the minimum of the action (42).

The critical temperature $`T_c`$ is defined from the condition $`\mathrm{\Gamma }=2`$, i.e. when the temperature is so high that no more than two particles are able to interconnect. In this case Eq. (43) turns into

$`k_\mathrm{B}T_c`$ $`={\displaystyle \frac{4\alpha 2^{8/3}+\frac{11}{2}\beta 2^{11/3}}{3\mathrm{ln}\xi }}`$ (46)
$`{\displaystyle \frac{12}{\pi }}\left({\displaystyle \frac{2\pi }{3}}\right)^{1/3}{\displaystyle \frac{Q^2n^{1/3}}{ϵ_0\mathrm{ln}\xi }}.`$
In the first approximation we may allow that in the critical point the fugacity $`\xi e^{\mu /k_\mathrm{B}T}`$ reaches its classical meaning $`\lambda _T^3n`$ (see expressions (35) and (36)). Expanding function $`\mathrm{ln}\xi `$ in a Taylor series, we obtain $`\mathrm{ln}\lambda _T^3n3(\lambda _Tn^{1/3}1)`$. Substituting the expansion in Eq. (46) and taking into account expression (36) for $`\lambda _T`$, we arrive at the following equation for $`T_c`$

$$k_\mathrm{B}T_c\frac{4}{\pi }\left(\frac{2\pi }{3}\right)^{1/3}\frac{e^2n^{1/3}}{ϵ_0}.$$
(47)
Now, let assign the concrete numerical values to parameters. Let the charge $`Q`$ be the elementary charge $`e`$. Put $`n=2.4\times 10^{19}`$ m<sup>-3</sup> and $`\gamma _{\mathrm{film}}=8.4\times 10^{14}`$ N/m. Then setting these values into expression (45) we get the number of electrons in a cluster: $`\mathrm{\Gamma }10^8`$. For the radius of the cluster we have $`R=(3\mathrm{\Gamma }/4\pi n)^{1/3}10^4`$ m. The critical temperature $`T_c`$ estimated from expression (47) is of the order of 100 K that is far beyond the critical temperature (4 K) of the liquid helium film. Thus the thermodynamic condition of the liquid helium stability sets limits on the real temperature $`T_c`$ of the cluster existence.

The behaviour of electrons in the liquid helium film has been studied by many researchers . An experimental appraisal of the mean number of electrons in a near-surface bubblon was approximately equal to $`10^8`$. The radius of the bubblon was estimated as $`10^2`$ cm. In such a manner our qualitative evaluation of the values $`\mathrm{\Gamma }`$ and $`R`$ is in agreement with the experimental data.

### 4.2 Shielding Coulomb potential

Let particles be repulsed by the shielding Coulomb potential $`u_{ss^{}}=1/4\pi ϵ_0\times Q^2e^{\kappa |r_sr_s^{}|}\times (r_sr_s^{})^1`$, where $`\kappa `$ is the effective radius of screening of nucleons in a nucleus, and attracted by the potential $`v_{ss^{}}=\frac{1}{2}\gamma (r_sr_s^{})^2`$ that can be created by an outside reason (so $`\gamma `$ is the force constant of an off-site elastic field). In this case the parameters $`a`$ and $`b`$ (33) of the action (34) become

$`a`$ $`=`$ $`3{\displaystyle \frac{Q^2}{4\pi ϵ_0gk_\mathrm{B}T}}{\displaystyle \underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}}{\displaystyle \frac{1}{x}}e^{\kappa gx}x^2𝑑x`$ (48)
$`=`$ $`{\displaystyle \frac{3Q^2}{4\pi ϵ_0\kappa ^2g^3k_\mathrm{B}T}}\left[\left(\kappa g\mathrm{\Gamma }^{1/3}+1\right)e^{\kappa g\mathrm{\Gamma }^{1/3}}+\left(\kappa g+1\right)e^{\kappa g}\right].`$
We will consider weight nuclei ($`\mathrm{\Gamma }>>1`$) and therefore $`\mathrm{\Gamma }^{1/3}`$ should be several times as large as the unit. Therefore, at the first approximation we may neglect terms with $`e^{\kappa g\mathrm{\Gamma }^{1/3}}`$ and then expression (48) becomes

$$a\frac{3}{4\pi ϵ_0}\frac{Q^2}{\kappa ^2g^3k_\mathrm{B}T}(\kappa g+1)e^{\kappa g}.$$
(49)
Now

$$b=3\frac{\gamma g^2}{2k_\mathrm{B}T}\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}x^2x^2𝑑x\frac{3}{10}\frac{\gamma g^2}{k_\mathrm{B}T}\mathrm{\Gamma }^{5/3}.$$
(50)
Let us rewrite these expressions in the form

$$a=\stackrel{~}{\alpha };b=\stackrel{~}{\beta }\mathrm{\Gamma }^{5/3}$$
(51)
where

$$\stackrel{~}{\alpha }=\frac{3}{4\pi ϵ_0}\frac{Q^2}{\kappa ^2g^3k_\mathrm{B}T}(\kappa g+1)e^{\kappa g};\stackrel{~}{\beta }=\frac{3}{10}\frac{\gamma g^2}{k_\mathrm{B}T}$$
(52)
and $`\alpha >>\beta `$. Let us now substitute expressions (49) and (50) for $`a`$ and $`b`$ in the action (34). With regard for inequality $`\stackrel{~}{E}>>1`$ and with allowance for $`\eta =+1`$ we will get

$$SK\left\{\frac{1}{2}\stackrel{~}{\alpha }\mathrm{\Gamma }^2\frac{1}{2}\stackrel{~}{\beta }\mathrm{\Gamma }^{11/3}+\mathrm{ln}\mathrm{\Gamma }+\mathrm{\Gamma }\mathrm{ln}\xi \right\}.$$
(53)
The extremum of the action $`S`$ is reached with resolving the equation $`S^{}(\mathrm{\Gamma })=0`$, or explicitly

$$\stackrel{~}{\alpha }\mathrm{\Gamma }\frac{11}{3}\stackrel{~}{\beta }\mathrm{\Gamma }^{8/3}+\frac{1}{\mathrm{\Gamma }}+\mathrm{ln}\xi =0.$$
(54)
Retaining two highest order terms in Eq. (54), we obtain

$`\mathrm{\Gamma }`$ $``$ $`\left({\displaystyle \frac{6}{11}}{\displaystyle \frac{\alpha }{\beta }}\right)^{3/5}`$ (55)
$`=`$ $`{\displaystyle \frac{4}{3}}\pi ^{2/5}n\left\{{\displaystyle \frac{15}{11}}{\displaystyle \frac{Q^2}{ϵ_0\kappa ^2\gamma }}\left[\left({\displaystyle \frac{3}{4\pi }}\right)^{1/3}{\displaystyle \frac{\kappa }{n^{1/3}}}+1\right]\right\}^{3/5}\mathrm{exp}\left[{\displaystyle \frac{3}{5}}\left({\displaystyle \frac{3}{4\pi }}\right)^{1/3}{\displaystyle \frac{\kappa }{n^{1/3}}}\right].`$
The critical temperature $`T_c`$ of the cluster destruction is defined from Eq. (54) if one puts $`\mathrm{\Gamma }=2`$:

$`k_\mathrm{B}T_c{\displaystyle \frac{4\alpha }{\mathrm{ln}\xi }}.`$ (56)
However, having found the fugacity $`\xi `$, we now cannot follow the procedure described in the previous subsection. In a nucleus, the temperature and the thermal de Broglie wavelength are others. In other words, we should exploit the laws of high energy physics, i.e. that the total energy of a nucleon is linked with the momentum of the nucleon by the relation $`E=\sqrt{p^2c^2+m^2c^4}`$. Hence in this case

$$\frac{3}{2}k_\mathrm{B}T=c\sqrt{p^2+m^2c^2}$$
(57)
and therefore

$$p3k_\mathrm{B}T/2c.$$
(58)
Thus the de Broglie thermal wavelength of the nucleon in a nucleus is determined as

$$\lambda _T=\frac{2h}{3ck_\mathrm{B}T}.$$
(59)
Expanding function $`\mathrm{ln}\xi =\mathrm{ln}(\lambda _T^3n)`$ in Eq. (56) in a Taylor series we obtain $`\mathrm{ln}\xi 3`$ and, therefore, equation (53) becomes

$`k_\mathrm{B}T_c`$ $`{\displaystyle \frac{4}{3}}\alpha `$ (60)
$`{\displaystyle \frac{4}{3}}{\displaystyle \frac{nQ^2}{ϵ_0\kappa ^2}}\left\{\left({\displaystyle \frac{3}{4\pi }}\right)^{1/3}{\displaystyle \frac{\kappa }{n^{1/3}}}+1\right\}\mathrm{exp}\left[\left({\displaystyle \frac{3}{4\pi }}\right)^{1/3}{\displaystyle \frac{\kappa }{n^{1/3}}}\right].`$
The results obtained in this section may account for the reasons of the atomic nucleus stability from the microscopic viewpoint. The shielding Coulomb potential $`u=1/4\pi ϵ_0\times Q^2e^{\kappa r}/r`$ is the typical nuclear (Yukawa) potential that provides for repulsion between protons. The potential $`v=\frac{1}{2}\gamma r^2`$, which is applied to nucleons, is ensured their mutual attraction. For example, setting $`\gamma =4\times 10^{17}`$ N/m we obtain from expression (55) $`\mathrm{\Gamma }30`$, which corresponds to the number of nucleons in a nucleus of zinc. The critical temperature of fission of a proton-neutron pair calculated from expression (60) results in $`T_c10^{11}`$ K. Note that this value is several times of magnitude greater than the typical temperature of a weight nucleus (the nucleus temperature is of the order of the Coulomb repulsion between protons in the nucleus).

### 4.3 Gravitating masses with Hubble expansion

Gravitational physics may also be assigned to the condensed matter . The gravitational attraction between particles, i.e. big masses, is beyond question. However the uniform character of the Hubble expansion is in doubt. We will use only the fact that an additional kinetic energy of particles $`E_{ss^{}}=\frac{1}{2}m(w_sw_s^{})^2`$ is associated with such an expansion. The relative velocity $`w_sw_s^{}`$ of particles which are found in points $`s`$ and $`s^{}`$, respectively, pertains to the relative distance between the particles, since

$$w_sw_s^{}=H(r_sr_s^{})$$

where $`H`$ is the Hubble constant. In this connection, we will consider a model system of identical gravitating masses, i.e. stars, with the attraction $`v_{ss^{}}=Gm^2/(r_sr_s^{})`$ and the effective repulsion $`u_{ss^{}}=\frac{1}{2}H^2m(r_sr_s^{})^2`$. For such kinds of interactions one has for the parameters $`a`$ and $`b`$ (33):

$`a=3{\displaystyle \frac{g^2H^2m}{2k_\mathrm{B}T}}{\displaystyle \underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}}x^2x^2𝑑x{\displaystyle \frac{3}{10}}{\displaystyle \frac{g^2H^2m}{k_\mathrm{B}T}}\mathrm{\Gamma }^{5/3};`$ (61)
$$b=3\frac{Gm^2}{gk_\mathrm{B}T}\underset{1}{\overset{\mathrm{\Gamma }^{1/3}}{}}\frac{1}{x}x^2𝑑x\frac{3}{2}\frac{Gm^2}{gk_\mathrm{B}T}\mathrm{\Gamma }^{2/3}$$
(62)
or

$$a=\stackrel{~}{\alpha }\mathrm{\Gamma }^{5/3};b=\stackrel{~}{\beta }\mathrm{\Gamma }^{2/3}$$
(63)
where

$$\stackrel{~}{\alpha }=\frac{3}{10}\frac{g^2H^2m}{k_\mathrm{B}T};\stackrel{~}{\beta }=\frac{3}{2}\frac{Gm^2}{gk_\mathrm{B}T}.$$
(64)
Inasmuch as gravitating masses are associated with the Bose statistics, i.e. $`\eta =1`$, the action (34) is rewritten as follows

$$SK\left\{\frac{1}{2}\stackrel{~}{\alpha }\mathrm{\Gamma }^{11/3}\frac{1}{2}\stackrel{~}{\beta }\mathrm{\Gamma }^{8/3}\mathrm{ln}\mathrm{\Gamma }+\mathrm{\Gamma }\mathrm{ln}\xi \right\}.$$
(65)
In expression (65) we omit the exponent term because of the assumption that the temperature $`T`$ of the universe is very close to the absolute zero (at small $`T`$ the exponent in expression (34) tends to zero). Then the equation $`S^{}(\mathrm{\Gamma })=0`$ becomes

$$\frac{11}{6}\stackrel{~}{\alpha }\mathrm{\Gamma }^{8/3}\frac{4}{3}\stackrel{~}{\beta }\mathrm{\Gamma }^{5/3}\frac{1}{\mathrm{\Gamma }}+\mathrm{ln}\xi =0.$$
(66)
Here, $`\alpha <<\beta `$. Retaining highest order terms in equation (66) we immediately obtain the solution

$$\mathrm{\Gamma }\frac{8}{11}\frac{\beta }{\alpha }=\frac{120\pi }{33}\frac{Gmn}{H^2}.$$
(67)
where $`n=(\frac{4\pi }{3}g^3)^1`$ is the concentration of particles, i.e. stars.

The equation for the critical temperature $`T_c`$ that determines the break-up of two coupled masses is obtained from Eq. (66) in which we set $`\mathrm{\Gamma }=2`$ (and the exponent term is neglected). With the inequality $`\alpha <<\beta `$, we get

$$k_\mathrm{B}T_c\frac{2^{2/3}}{3}\frac{\beta }{\mathrm{ln}\xi }\left(\frac{2\pi }{3}\right)^{1/3}\frac{Gm^2n^{1/3}}{\mathrm{ln}\xi }.$$
(68)
By definition, the fugacity $`\xi =e^{\mu /k_\mathrm{B}T}`$. At the transition from the Bose statistics to the classical one, the chemical potential gradually decreases with temperature starting from a value $`\mu 0`$ (but still $`\mu <0`$) to $`\mu <<0`$. So, the function $`\mathrm{ln}\xi `$ in Eq. (68) is strongly negative: $`\mathrm{ln}\xi =|\mu |/k_\mathrm{B}T`$. This signifies that the critical temperature in the system of gravitating masses is absent. In other words, in the universe, masses trends to the aggregation at any conditions.

Thus the results obtained demonstrate that the availability of the Hubble expansion is the deciding factor in the formation of galaxies in the universe (note that without $`H`$ masses in the universe will trend merely to the total mergence and any spatial distribution will not be observed). Actually, expression (67) shows that Hubble’s constant, $`H50`$ km/(s$``$Mpsc) = $`1.6\times 10^{18}`$ s<sup>-1</sup>, is one of the main parameters that determines the steady state of masses in a cluster, i.e., galaxy. So the greater the mean mass $`m`$, the greater the number $`\mathrm{\Gamma }`$ of stars in a galaxy.

## 5 Conclusion

In the present paper the statistical description of a system of interacting particles, which allows the study of their spatial nonhomogeneous distribution, has been proposed. One of the most important properties of our model is the presentation of the total potential of interacting particles in the form of the paired energy of attraction and repulsion. We have started from the Hamiltonian that describes interacting particles located in knots of the three-dimensional lattice and applied the quantum field theory methodology based on the functional integrals for the description of the partition function of systems studied.

We have shown that the solutions of equations characterizing the behaviour of the field variables yield the full resolving the initial statistical task. The solutions may correspond to both the homogeneous and nonhomogeneous distributions of particles and only the character of interparticle interaction and the temperature define the kind of the solution.

The free energy represented through field variables in the saddle-point has allowed us to investigate special conditions, which provide for the formation of clusters in a system with the initially homogeneous particle distribution. The conditions are answerable also for the size of clusters and the temperature of phase transition to the cluster state.

The approach proposed does not need fitting parameters: the cluster size and the temperature of phase transition to the cluster state are completely determined by the value and character of paired potentials in the systems of interacting particles. Besides, the introduction of a peculiar type of the combined variable has simplified the investigation procedure. Specifically, the procedure has provided a way for the obtaining the necessity criterion of the transition of a system of separate particles into the spatial nonhomogeneous state.

The authors are very thankful to the referees for the concrete constructive remarks that allowed a significant improvements of the contents of the paper.


format: text

text:
# Pattern of Reaction Diffusion Front in Laminar Flows

## Abstract

Autocatalytic reaction between reacted and unreacted species may propagate as solitary waves, namely at a constant front velocity and with a stationary concentration profile, resulting from a balance between molecular diffusion and chemical reaction. The effect of advective flow on the autocatalytic reaction between iodate and arsenous acid in cylindrical tubes and Hele-Shaw cells is analyzed experimentally and numerically using lattice BGK simulations. We do observe the existence of solitary waves with concentration profiles exhibiting a cusp and we delineate the eikonal and mixing regimes recently predicted.

The motion of interfaces and the propagation of fronts resulting from chemical reactions occur in a number of different areas scott , including population dynamics fisher ; kolmogorov  and flame propagation zeldo . It is known that autocatalytic reaction fronts between two reacting species propagate as solitary waves, namely at a constant front velocity and with a stationary concentration profile hanna ; martin . The important issue of the selection of the front velocity was addressed earlier on, but only a few cases are well understood, such as the pioneering works of Fisher fisher  and Kolmogorov-Petrovskii-Piskunov kolmogorov  on a reaction-diffusion equation with second-order kinetics scott ; zeldo ; saarloos . The effect of advective flow (inviscid and/or turbulent) on reacting systems was analyzed extensively in the propagation of flames in the context of combustion zeldo ; clavin . On the other hand, advective effects on the behavior of autocatalytic fronts have been only recently addressed berest ; abel ; edwards . B. F. Edwards edwards  studied theoretically the effect of a $`2D`$ laminar flow on an autocatalytic reaction front between two infinite planes separated by a gap $`b`$. In this geometry, the velocity profile is unidirectional in the direction $`z`$ of the flow and is given by Poiseuille’s equation, $`\stackrel{}{U}=U_M(1\zeta ^2)\stackrel{}{z}`$ where $`U_M=1.5\overline{U}`$ is the maximum velocity, $`\overline{U}`$ is the mean velocity, $`\zeta =2x/b`$ is the transverse normalized coordinate and $`\stackrel{}{z}`$ is the unit vector parallel to the flow, chosen as the direction of the front propagation in the absence of flow (see below). Consider the iodate-arsenous acid reaction described by a third-order autocatalytic reaction kinetics scott ; hanna ; martin :

$$\frac{C}{t}+\stackrel{}{U}\stackrel{}{}C=D_m\mathrm{}C+\alpha C^2(1C)$$
(1)
where $`C`$ is the concentration of the (autocatalytic) reactant iodide, normalized by the initial concentration of iodate, $`D_m`$ is the molecular diffusion coefficient, and $`\alpha `$ is the reaction rate kinetic coefficient. In the absence of hydrodynamics ($`\stackrel{}{U}=\stackrel{}{0}`$), Eq.1 admits a well-known solitary wave solution with front velocity $`V_0=\sqrt{\alpha D_m/2}`$ and front width $`L_0=D_m/V_0`$ hanna ; martin . The use of these two quantities to normalize velocities and lengths in Eq.1, leads to two independent parameters $`\eta =b/2L_0`$ and $`\epsilon =\overline{U}/V_0`$. Reference edwards  investigated numerically the solitary wave solution of Eq.1, and particularly its normalized front velocity, $`v=V_F/V_0`$, as a function of $`\epsilon `$, for different values of $`\eta `$. Of interest are the following asymptotic predictions:

In the narrow-gap regime ($`\eta 0`$ or $`\epsilon 0`$), it was found that $`v=1+\epsilon `$. Namely, when $`L_0>>b`$, mixing across the gap is significant, the concentration front is flat and advected by the mean flow, yielding: $`V_F=V_0+\overline{U}`$.

On the other hand, in the wide-gap regime ($`\eta 1`$), the front is thin and curved across the gap, and Eq.1 can be replaced by the eikonal equation:

$$\stackrel{}{V_F}.\stackrel{}{n}=V_0+\stackrel{}{U}.\stackrel{}{n}+D_m\kappa $$
(2)
where $`\stackrel{}{n}`$ is the unit vector normal to the thin front (oriented from reacted to unreacted species) and $`\kappa `$ the front curvature. In this regime, to leading order and neglecting the local front curvature, the front velocity is given by the simplified $`1D`$ eikonal equation:

$$V_F=V_0/\mathrm{cos}\theta +U(\zeta )$$
(3)
where $`\theta `$ is the angle between $`\stackrel{}{n}`$ and the flow direction and $`U(\zeta )`$ is the advection velocity. Under these conditions, reference edwards  predicted two behaviors depending on the flow direction: For a supporting flow ($`\epsilon >0`$), $`V_F=V_0+U_M`$, which means that the front is advected at the largest possible velocity. The front shape across the gap is then given by the solution of Eq.3. For an adverse flow ($`\epsilon <0`$), $`V_F=V_0`$, which also represents the maximum algebraic velocity one could have expected physically. The front is perpendicular to the walls ($`\theta =0`$ at $`\zeta =\pm 1`$), and presents a cusp in the middle of the gap (discontinuity of $`\theta `$ at $`\zeta =0`$). Here, the adverse flow elongates the front but does not slow it down. Note that curvature effects ($`D_m\kappa `$ in Eq.2) smooth the cusp, but do not otherwise alter these predictions. Note also that similar features would occur for other kinetics such as FKPP or Arrhenius ones berest .

The objective of the present letter is to experimentally test the above $`2D`$ predictions using two different devices, namely Hele-Shaw cells and cylindrical tubes. The case of the Hele-Shaw cell, consisting of two parallel plates separated by a gap $`b`$ small compared to the other dimensions, is supposed to be quantitatively addressed by edwards . Alternatively, the case of the cylindrical tube of inner radius $`a`$ (in which the flow field is also described by Poiseuille’s equation, with $`U_M=2\overline{U}`$ and $`\zeta =r/a`$) represents a genuine (axisymmetric) $`2D`$ situation. Experiments in Hele-Shaw cells are discussed with the help of lattice BGK simulations nous  of Eq.1 for a $`3D`$ flow.
In the experiment, the front is detected by using starch, at small concentrations, which reacts in the presence of iodine leading to a dark blue signature of the front passage martin ; nous . First, we consider the reaction in the absence of advection by the flow ($`\stackrel{}{U}=\stackrel{}{0}`$). As expected, we do observe solitary fronts propagating with flat shapes. In accordance with hanna ; martin , their velocity is $`V_00.02mm/s`$, from which one can estimate their front thickness $`L_0=D_m/V_00.1mm`$ ($`D_m2.10^9m^2/s`$). Because the reaction products have a lower density than the unreacted species, the hydrodynamically stable situation corresponds to descending fronts in vertical tubes. In the following, we focus on the interplay between advection and propagating fronts. To minimize the effect of density contrast, we studied the propagation of buoyantly stable fronts in small cells. We used vertical Hele-Shaw cells of size $`b\times W=0.1\times 1,\mathrm{\hspace{0.17em}\hspace{0.17em}0.2}\times 4,\mathrm{\hspace{0.17em}\hspace{0.17em}0.4}\times 8,\mathrm{\hspace{0.17em}\hspace{0.17em}1}\times 15mm^2`$ and circular capillary tubes of radius $`a=0.3,0.58,0.88`$ and $`1.9mm`$. A constant advecting flow, upwards or downwards, was fixed by a syringe. Note that these cells are small enough to prevent flattening of the front due to buoyancy, but large enough to enable a constant flow rate injection with our injection device. The average velocity of the imposed flow ranged between $`0`$ and $`60V_0`$.
In a $`3D`$ Hele-Shaw cell, the flow velocity profile is unidirectional and depends on the two transverse coordinates, $`x`$ and $`y`$ nicole3 . The profile across the gap is almost parabolic with a gap average value uniform over the width $`W`$, except in a boundary layer of order $`b`$, within which the velocity vanishes (see the gap-average profile on Fig.4). For the three aspect ratios studied, $`W/b=10,15,20`$, we can estimate from nicole3 , $`U_M/\overline{U}=1.60`$, $`1.57`$ and $`1.55`$, respectively.

We observed solitary waves in the whole range of flow rates investigated. Typical fronts are shown in Fig.1 in the plane of Hele-Shaw cells (top diagram) and in tubes (bottom diagram). For each geometry, two adverse flows (on the left) and two supportive flows (on the right) are displayed. The front shape always points toward the same direction as the underlying flow field, while its distortion increases with flow intensity (recall that the fronts are flat in the absence of flow, $`\overline{U}=0`$). The triangular shapes observed, in the case of adverse flows, are reminiscent of premixed flames mungal .

The two graphs in Fig.1 show the normalized front velocities $`v=V_F/V_0`$ versus $`\epsilon =\overline{U}/V_0`$, measured for different sizes of the Hele-Shaw cells and tubes. For each data set (given $`\eta `$), the front velocity increases linearly with the flow rate, but with a slope different for supportive and adverse flows. This change of slope is in accordance with edwards  as well as the observed evolution of the slopes with $`\eta `$. In addition, most of the data fall in the domain delimited by the asymptotic regimes described above. The exceptions correspond to supportive flows in the smaller tubes ($`1<\eta <6`$), which, unlike Edward’s $`2D`$ numerical data, fall below the mixing regime (predicted for $`\eta 1`$). This difference needs to be further analyzed, given that it is difficult to achieve experimentally a very low constant flow rate injection in the small tubes. At the same time, some authors abel  have suggested that the mixing straight line should be higher. Using the Peclet number ($`Pe=\overline{U}a/D_m=\epsilon \eta `$), which compares the relative importance of advection and diffusion, and the Damköhler number ($`Da=\alpha a/\overline{U}=2\eta /\epsilon `$), which is the ratio of advective to reactive time-scales, they predicted abel  that for $`Da1`$ (which is not attainable in our experiments), the front velocity should be the product of $`V_0`$ by the Taylor dispersion factor Taylor , which accounts for the coupling between advection and transverse diffusive mixing. This factor would then enhance the front velocities.
The front velocities measured in the Hele-Shaw cells are very close to the mixing regime. However, in the case of adverse flows, the measured values exhibit some departure toward the eikonal regime ($`V_F=V_0`$) when either $`\eta `$ or $`\epsilon `$ is increased. This trend is even more pronounced for the tubes, which present larger $`\eta `$ values than for the Hele-Shaw cells, in accordance with predictions edwards . On the other hand, all the values measured for supportive flows, even large, ($`\epsilon `$ up to $`50`$ for $`\eta =0.5`$ displayed in Fig.2) fall on the asymptotic mixing regime predicted by the strictly $`2D`$ gap analysis edwards  ($`V_F=V_0+\overline{U}`$ for $`\eta 0`$). This is all the more surprising since the flatness of the fronts, expected in the mixing regime, was not observed in the Hele-Shaw cell plane (top right of Fig.1). However, in our experiments, although the normalized gap $`\eta =b/2L_0`$ introduced in the $`2D`$ gap analysis edwards  is small ($`0.5,1,2,5`$), the normalized width $`W/2L_0`$ is large ($`5,20,40,75`$). Extrapolating Edwards’s $`2D`$ gap analysis to our $`3D`$ case would suggest that our experiments combine a mixing regime across the gap with an eikonal regime across the width. Under these assumptions, the shape and velocity of the front would obey an equation similar to Eq.2, namely

$$V_F=(V_0+D_m\kappa )/\mathrm{cos}\theta +U^{2D}(y)$$
(4)
where the effective advection velocity $`U^{2D}(y)`$ is the gap-averaged velocity defined in nicole3 , and where $`\kappa `$ represents now the curvature of the $`2D`$ front curve observed in the plane of the Hele-Shaw cell. The front velocity for $`\epsilon >0`$ would then be set by the maximum $`U_M^{2D}`$ of the profile $`U^{2D}(y)`$, found in the middle of the plane (for $`\theta (y=0)=0`$ and $`\kappa (y=0)=0`$). The so-obtained maximum front velocity, $`V_F=V_0+U_M^{2D}`$, expected in the asymptotic width-eikonal regime ($`W/2L_0>>1`$) can be compared to the velocity, $`V_F=V_0+\overline{U}`$, expected in the width-mixing regime ($`W/2L_0<<1`$). One finds that these two asymptotic velocities would be equal in the Hele-Shaw limit ($`W/b\mathrm{}`$), and are actually very similar in our experiments (as $`U_M^{2D}/\overline{U}=1.07`$, $`1.05`$ and $`1.03`$ for the three aspect ratios used). This could justify that the parameter $`\eta =b/2L_0`$ introduced in the $`2D`$ gap analysis edwards , actually controls the front velocity in $`3D`$ Hele-Shaw cells.

We tested the ability of the full description and the simplified one, given respectively by Eqs. 1 and 4, to account for both shape and velocity of the experimental fronts. As lattice BGK simulations have been used to obtain the solutions of Eq.1, we have first validated this numerical method, by reproducing Edwards’s results edwards  on the shape and velocity of the fronts propagating between two infinite planes ($`2D`$ simulations). Then, $`3D`$ lattice BGK simulations of Eq.1 and numerical integration of the $`2D`$ eikonal Eq.4 were performed, using respectively the analytical stationary $`3D`$ flow field given by nicole3 , and its gap-average $`U^{2D}(y)`$. The front obtained with the latter method was compared to the iso-concentration $`C=0.5`$ of the gap-averaged concentration map produced by the lattice BGK simulations. Fig.3 displays these fronts, for the same parameters ($`\epsilon `$, $`\eta `$, $`W/L_0`$) as in one typical experiment, in the cases of adverse ($`\epsilon =4.8`$) and supportive ($`\epsilon =+4.8`$) flows. Note that in these cases of interest, for which $`W/L_0`$ is finite, the integration of Eq.4 requires the value of the front velocity (thus fixing the value of the curvature $`\kappa `$ at the integration starting point). Hence, Eq.4 is not fully predictive, but links the shape of the front to its velocity. The shape and velocity predicted by lattice BGK simulations and the ones given by Eq.4 are found to compare fairly well with the experimental observations in the case of adverse flow. However, for the supportive flow case, although the front velocity is correctly predicted, the two numerical predictions, similar for the shape, fail to account for the experimental observations. We believe that this discrepancy might be due to an alteration of the flow velocity profile caused by a triangular meniscus which appears on the top of the solution in our supportive flow experiments. From its shape and its distance to the front (typically several tens of $`W`$), one can infer that the meniscus could introduce a few percents of excess fluid velocity in the middle of the cell plane. The resulting non-uniformity in $`U^{2D}(y)`$ may account for the rounded shape observed in the supportive flow experiments (see Fig.4).

We have also analyzed the dynamics of the shape formation in the case of adverse flows. Fig.4 displays the time development of the iso-concentration $`C=0.5`$, initially flat, toward the stationary triangular shape. The sequence shows that an early determination of both the final front velocity and the final angle $`\theta `$ is achieved as soon as the profile is altered over a typical distance $`b`$ from the side walls. This supports the contention that $`W`$ plays no role in the determination of both shape ($`\theta `$) and velocity in the regimes under consideration. This was confirmed by simulations in wider lattices which produced the same values of the velocity and $`\theta `$. Thus $`\eta =b/2L_0`$ is effectively the relevant parameter in Hele-Shaw cells (for which $`W/b1`$).
In conclusion, we have performed experiments and lattice BGK simulations of autocatalytic reaction fronts in laminar advective flow fields in Hele-Shaw cells and circular tubes. Solitary waves were observed in the entire range of flow rates. For flows adverse to the chemical front propagation, we observed cusp-like fronts in tubes and triangular fronts in the plane of Hele-Shaw cells. Our measurements of the front velocity agree with the $`2D`$ asymptotic predictions edwards , in the limiting cases where either diffusion overcomes reaction ($`\eta 1`$) or it is negligible ($`\eta 1`$). It would be interesting to extend the range of the cell sizes. Larger cells could be used to study the buoyancy stabilizing effect and smaller cells within the scope of microfluidics.

This paper benefited from discussions with Professor K. Showalter. The work was partly supported by IDRIS (project 024052), CNES No 793/CNES/00/8368, ESA (No AO-99-083). MRT grant (M. L.). All these sources of support are gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

The problem of how a charm-quark fragments into a charmed hadron is challenging both from the theoretical and the experimental point of view. Indeed, perturbative QCD is not applicable at energies comparable with the charm-quark mass. Therefore, parameterizations have to be used and the parameters are determined experimentally. In the experiments the direct identification of the charmed hadron in the final state is only possible through the visual observation of the hadron decay and the measurement of the kinematical variables, practically feasible only with nuclear emulsions.

In the following, we focus on the so-called charm-production fractions ($`f_h`$’s); i.e. the probability that a charm-quark fragments into a charmed hadron $`h`$ ($`=D^0,D^+,D_s,\mathrm{\Lambda }_c`$).

In this paper we review all existing data on charm-production fractions as measured by $`e^+e^{}`$, $`\pi N`$ and $`\gamma N`$ experiments and predict $`f_h`$ in neutrino induced charm-production. Indeed, data on charmed fractions in interactions induced by neutrinos are rather scarce. Only one experiment, E531 , measured $`f_h`$ with a statistics of 122 events with an identified charmed hadron in the final state.

After a thorough overview of all available data, we present a method to extract $`f_h`$ for neutrino experiments. The major difference of this method with respect to other presented in the past, i.e. see Ref. , is that$`f_h`$’s from $`e^+e^{}`$ are not extrapolated to neutrino experiments straightway, but some peculiarities of $`\nu `$-induced interactions are accounted for. Indeed, neutrinos may undergo to quasi-elastic and diffractive processes with production of charmed hadrons in the final state. These results are also used to estimate the semi-muonic branching ration of charmed hadrons. Finally, we discuss the statistical and systematic uncertainties associated with our predictions and compare them with the E531 results.

## 2 Charm-production fraction in DIS interactions

### 2.1 Measurement of the $`D^+/D^0`$ ratio

In this Section we give an overview of the available data on $`D^+/D^0`$ measurements in $`e^+e^{}`$, $`\pi N`$, $`pN`$ and $`\gamma N`$ experiments. Furthermore, we also study the correlated variable $`F_V`$, which is defined as $`V/(V+P)`$ , where $`V`$ and $`P`$ signify vector and pseudo-scalar charmed mesons, respectively. When needed, we recomputed $`F_V`$ by using the latest charmed-hadron branching-ratios (BR). Finally, we extract for the first time $`F_V`$ from neutrino induced charm-production data.

In the following we assume that $`D^0`$ and $`D^+`$ production cross-sections are equal at the parton level, as well as the direct production of $`D^0`$ and $`D^+`$.

From here on, we indicate with $`D`$ the sum of both prompt and decayed ($`D^{}D`$) $`D`$ meson production.

Under the previous assumptions, the measurement of $`D^+/D^0`$ and $`D^{}/D`$ ratios allows us to extract $`F_V`$. By using the formulae given in Ref. , $`F_V`$ can be extracted from the relations

$$F_V\times B_{}=\frac{1R_1}{1+R_1}$$
(1)
$$F_V\times B_{}=\frac{R_2}{1R_2}$$
(2)
where $`R_1D^+/D^0`$, $`R_2(D^0\text{from}D^+)/D^0`$ and $`B_{}=BR(D^+D^0\pi ^+)=(67.7\pm 0.5)\%`$ .

#### 2.1.1 $`e^+e^{}`$ experiments

The basic principle to reconstruct $`e^+e^{}`$ events with charmed hadrons in the final state is common to all the experiments. Being a $`c\overline{c}`$ pair produced in the annihilation, one charmed hadron is used to tag the events, while the other one to study the decay properties. The decay modes used to tag the event are 111Here and in the following by $`D^+(D^+)`$ we implicitly indicate also the charge conjugated decay modes.

$$D^+D^0\pi ^+(K^{}\pi ^+)\pi ^+(0.026\pm 0.006)$$

$$D^0K^{}\pi ^+(0.0383\pm 0.0009)$$

$$D^0K^{}\pi ^+\pi ^+\pi ^{}(0.0749\pm 0.0031)$$

$$D^+K^{}\pi ^+\pi ^+(0.090\pm 0.006)$$

where the corresponding branching ratios are also given in brackets .

Results on the total cross-sections for inclusive production of the charmed particles $`D^0`$, $`D^+`$, $`D^0`$ and $`D^+`$ at various $`\sqrt{s}`$ are shown in Table 1.

A complete review of the probabilities ($`f(cC)`$) that a c-quark fragments into a $`D^{}`$, $`D^0`$, $`D^+`$ and other charmed hadrons as measured in $`Z^0`$ decays is given in Ref.  and reported in Table 2.

From Tables 1 and 2 we can extract both $`R_1`$ and $`F_V`$, the latter being estimated by using Eqs. (1) and (2). The results are given in Table 3 and show that within the experimental errors, both $`R_1`$ and $`F_V`$ are independent of the energy.

#### 2.1.2 $`\pi N`$ and $`\gamma N`$ experiments

Several experiments have studied charm-production and extracted $`R_1`$ and $`F_V`$ by using $`\pi `$ beams of different energies impinging onto different targets. A non-exhaustive list of all available data in $`\pi N`$ experiments is given in Table 4.

The NA14/2 photo-production experiment  measured both $`R_1=0.37\pm 0.10`$ and $`R_2=0.26\pm 0.04`$ from which we can extract the weighted average $`F_V=0.57\pm 0.09`$.

From these data we can conclude that, within the experimental errors, $`R_1`$ is both process- and energy-independent.

#### 2.1.3 $`\nu N`$ experiments

Recently two measurements which allowed us to extract for the first time $`F_V`$ from neutrino experiments became available. In Ref.  the CHORUS Collaboration presented a measurement of the production rate of $`D^0`$ based on a sample of about $`26000`$ $`\nu _\mu `$ charged-current events interactions located and analyzed so far in the target emulsions. After reconstruction of the event topology in the vertex region, $`283D^0`$ decays were observed with an estimated background of 9.2 events from $`K^0`$ and $`\mathrm{\Lambda }`$ decays. The CHORUS Collaboration measured the $`D^0`$ production cross-section times $`BR(D^0V2)+BR(D^0V4)]`$ . The total cross-section has been extracted by accounting for the $`D^0`$ decays into all neutrals. The value we used is $`BR(D^0\text{all neutral})=(25\pm 5)\%`$ .

Therefore, the $`D^0`$ production cross-section normalized to $`\nu _\mu `$ charged-current (CC) interactions is

$$\frac{\sigma (D^0)}{\sigma _{CC}}=(2.65\pm 0.18\pm 0.24\pm 0.50)\times 10^2$$
(3)
at $`27GeV`$ average $`\nu _\mu `$ energy. Notice that this measurement includes both $`D^0`$ prompt and $`D^0`$ from the decay of $`D^{}`$ mesons.

The $`D^+`$ production in $`\nu _\mu `$ charged-current interactions has been measured, with a similar $`\nu _\mu `$ beam, by the BEBC  and NOMAD  experiments to be $`(1.22\pm 0.25)\times 10^2`$ and $`(0.79\pm 0.20)\times 10^2`$, respectively. The weighted average $`D^+`$ production rate normalized to $`\nu _\mu `$ charged-current interactions is

$$\frac{\sigma (D^+)}{\sigma _{CC}}=(0.96\pm 0.16)\times 10^2$$
(4)
From the measured ratios (3) and (4), and by knowing $`B_{}`$, we can compute $`R_2=0.25\pm 0.06`$. From the latter value and from Eq. (2) we can extract

$$F_V=0.50\pm 0.12$$

It is worth to notice that $`R_1`$ and $`F_V`$ extracted from neutrino experiments can be compared straightway to $`e^+e^{}`$, $`\pi N`$ and $`\gamma N`$ results, being $`D^+`$ and $`D^0`$ either produced promptly or from the decay of prompt $`D^+`$ and $`D^0`$. Namely, processes peculiar of $`\nu `$ interactions do not affect $`R_1`$ and $`F_V`$ 222 In neutrino interactions $`D^{()+}`$ may also be produced diffractively but, due to the $`V_{cd}`$ suppression, its rate is expected to be about $`(1.6\pm 0.3)\times 10^4`$ with respect to CC interactions and therefore negligible..

#### 2.1.4 Summary and discussion of all available data on $`R_1`$ and $`F_V`$

From results reported in the previous Sections, we can argue that within the experimental errors $`R_1`$ is constant over a wide range of energies ($`\sqrt{s}4÷90`$ GeV) and independent of the process. The constant behavior of $`R_1`$ down to $`\sqrt{s}4`$ GeV can be derived with simple arguments:

* the masses of $`D^+(1869.3)`$, $`D^0(1864.5)`$, $`D^+(2010.0)`$ and $`D^0(2006.7)`$ are very similar. Therefore, the threshold suppression of $`D^{}`$ mesons which tends to enhance $`D^+`$ contribution, is very little;
* whatever $`D^{}/D`$ meson is produced a pion should be always created. Therefore, all charmed mesons have the same threshold behavior, which cancels out in the ratio.

For these reasons we assume that $`R_1`$ measured in $`e^+e^{}`$ (see Table 3) can be used in neutrino induced charm-production, too.

As an important by product of our study we have also extracted $`F_V`$ from different processes and at several energies (see Table 5). The simplest model to predict $`F_V`$ is based on the spin-counting. Namely, vector mesons are spin-one states ($`{}_{}{}^{3}S_{1}^{}`$), while pseudo-scalar mesons are spin-zero states ($`{}_{}{}^{1}S_{0}^{}`$), therefore $`F_V=0.75`$. The discussion of more refined models (UCLA, JETSET, HERWIG and others) is beyond the purposes of this paper. For details we refer to .

From Table 5 we can see that the measured $`F_V`$ is independent of the processes and of the energy. This means that the probability for a c-quark to fragment into a $`D`$ or a $`D^{}`$ meson is universal and does not depend neither on the process nor on the energy. Notice that the UCLA model is the best in describing available $`e^+e^{}`$ data.

### 2.2 Measurement of the $`D_s`$ to $`D^0`$ and $`\mathrm{\Lambda }_c`$ to $`D^0`$ ratios

The ratio $`D_s/D^0`$ has been measured in $`e^+e^{}`$, $`\pi N`$ and $`\gamma N`$ experiments. A summary of the available data is given in Table 6.

The decay mode used to tag the event is

$$D_s^+\varphi \pi ^+(K^{}K^+)\pi ^+$$

whose BR, as reported by the Particle Data Group, are :

$$BR(D_s^+\varphi \pi ^+)=0.036\pm 0.009$$

$$BR(\varphi K^{}K^+)=0.492\pm 0.007$$

From Table 6 we can see that the $`D_s/D^0`$ ratio is, within the errors, independent of the energy and of the process.

Data on the ratio $`\mathrm{\Lambda }_c/D^0`$ are very poor. Indeed, it has been measured only in $`e^+e^{}`$ experiments (see Table 6).

The decay mode used to tag the event is

$$\mathrm{\Lambda }_cpK^{}\pi ^+$$

whose BR, as reported by the Particle Data Group, is :

$$BR(\mathrm{\Lambda }_cpK^{}\pi ^+)=0.050\pm 0.013$$

From Table 6 we can see that, although with a smaller statistical accuracy, both $`R_s`$ and $`R_c`$ are, within the errors, independent of the energy.

#### 2.2.1 Summary and discussion of all available data on $`R_s`$ and $`R_c`$

Although from Table 6 it seems that $`R_s`$ is constant over a wide range of energies ($`\sqrt{s}4÷90`$ GeV) and independent of the process, some comments are in order.

Given the quark composition of the $`D_s^+(c\overline{s})`$ meson, it has to be created always together with at least one $`K`$ meson. Therefore, being $`m_K500`$ MeV, we expect that the threshold effect for $`D_s`$ production is more pronounced than for $`D`$ mesons, when at least one $`\pi `$ has to be produced ($`m_\pi 100`$ MeV). To account for the different threshold effect at low energies, in the following we do not use the $`R_s`$ values measured at the $`Z^0`$ peak. Furthermore, under the assumption that $`\sqrt{s}`$ in collider experiments can be replaced with $`W`$ in fixed target experiments, $`W`$ being the final state hadronic mass, and noting that neutrino-induced charm events at present experiments are characterized by values of $`W`$ in the range $`4÷10`$ GeV, we can argue a $`R_s`$ value of $`R_s=0.171\pm 0.029`$ which corresponds to the weighted average of the first three results in Table 6.

The available measurements on $`R_c`$ are very poor. Nevertheless, as it will be discussed in Section 3, we do not use the $`R_c`$ value to predict the charmed fractions in events induced by neutrinos.

## 3 Charm-production fractions in neutrino experiments

### 3.1 The method

From the previous Sections we can argue that once a charm-quark has been produced in deep-inelastic interactions (i.e. the energy of the process is higher than the threshold), it has a probability to produce a charmed hadron $`C_h`$ which is, within the experimental errors, independent of the process and of the energy. Therefore, as far as the deep-inelastic scattering (DIS) is concerned, we can write the charm production rate as

$`{\displaystyle \frac{\sigma _c(DIS)}{\sigma _{CC}}}`$ $`=`$ $`{\displaystyle \frac{\sigma _{D^0}}{\sigma _{CC}}}+{\displaystyle \frac{\sigma _{D^+}}{\sigma _{CC}}}+{\displaystyle \frac{\sigma _{D_s}}{\sigma _{CC}}}+{\displaystyle \frac{\sigma _{\mathrm{\Lambda }_c}}{\sigma _{CC}}}`$ (5)
$`=`$ $`{\displaystyle \frac{\sigma _{D^0}}{\sigma _{CC}}}\times \left(1+R_1+R_s+R_c\right).`$
In the case of neutrino induced charm-production, Eq. (5) is not correct. Indeed, in this case we also have to account for diffractive and quasi-elastic charm-production. Therefore, the inclusive charm-production rate in neutrino interactions can be written as

$$\frac{\sigma _c(\nu )}{\sigma _{CC}}=\frac{\sigma _c(DIS)}{\sigma _{CC}}+\frac{\sigma _{D_s}}{\sigma _{CC}}_{diff}+\frac{\sigma _{\mathrm{\Lambda }_c}}{\sigma _{CC}}_{QE}.$$
(6)
If one accepts that $`R_1`$, $`R_s`$ and $`R_c`$ from $`e^+e^{}`$ (with $`\sqrt{s}=4.1÷90`$ GeV) and other experiments can be used to described the fragmentation of charm-quarks produced in DIS neutrino interactions with average final state hadronic mass $`W10`$ GeV, then

$`{\displaystyle \frac{\sigma _c(\nu )}{\sigma _{CC}}}`$ $`=`$ $`{\displaystyle \frac{\sigma _{D^0}}{\sigma _{CC}}}\times \left(1+R_1+R_s+R_c\right)+{\displaystyle \frac{\sigma _{D_s}}{\sigma _{CC}}}_{diff}+{\displaystyle \frac{\sigma _c}{\sigma _{CC}}}_{QE}`$ (7)
From Eq. (7) charm-production fractions in neutrino interactions can be written as

$$f_{D^0}(E_\nu )=\frac{\sigma _{D^0}}{\sigma _{CC}}(E_\nu )\times \frac{1}{\frac{\sigma _c(\nu )}{\sigma _{CC}}(E_\nu )}$$

$$f_{D^+}(E_\nu )=R_1\times \frac{\sigma _{D^0}}{\sigma _{CC}}(E_\nu )\times \frac{1}{\frac{\sigma _c(\nu )}{\sigma _{CC}}(E_\nu )}$$

$$f_{D_s^+}(E_\nu )=\left(R_s\times \frac{\sigma _{D^0}}{\sigma _{CC}}(E_\nu )+\frac{\sigma _{D_s}}{\sigma _{CC}}_{diff}(E_\nu )\right)\times \frac{1}{\frac{\sigma _c(\nu )}{\sigma _{CC}}(E_\nu )}$$

Notice that, given the poor knowledge on $`R_c`$ and on the quasi-elastic charm-production cross-section, we derive $`f_{\mathrm{\Lambda }_c}`$ by using the normalization constrain $`f_{D^0}+f_{D^+}+f_{D_s^+}+f_{\mathrm{\Lambda }_c^+}=1`$. In order to estimate the charmed fractions and their energy dependence, we use the following inputs

* the inclusive charm-production rate derived in Ref. ;
* the energy dependence of the $`D^0`$ production rate reported in Ref.  properly scaled to account for the effect discussed in Section 2.1.3;
* the energy dependence of the diffractive $`D_s`$ production given in , properly scaled in order to reproduce the average diffractive charm-production cross-section as measured by BEBC and NuTeV  ($`\frac{\sigma _{D_s}}{\sigma _{CC}}_{diff}=(0.31\pm 0.05)\times 10^2`$).

### 3.2 Results on $`f_h`$ and $`B_\mu `$ and comparison with the data

By using the method described in the previous section, we derived the charm-production fractions, as a function of the neutrino energy, as reported in Table 7 and shown in Fig. 1. Our results are in good agreement with the charm-production fractions extracted from the E531 data, see Fig. 1. It is worth noticing that $`f_{\mathrm{\Lambda }_c^+}`$ shows a dependence on the energy, higher values at low neutrino energies, consistent with the expectations. Indeed, quasi-elastic charm-production, which yields only $`\mathrm{\Lambda }_c`$, is expected to largely contribute to $`\sigma _c(\nu )`$ for $`E_\nu <25`$ GeV . In Table 7 the expected charm-production fractions in the CHORUS experiment are also given.

Having determined the $`f_h`$’s, we can also estimate the semi-muonic branching ratio $`B_\mu `$ of the charmed hadrons as a function of the neutrino energy. $`B_\mu `$ is a very important quantity, being the input variable needed to extract from the dimuon data the element of the CKM matrix $`V_{cd}`$.

Recently, a direct measurement of $`\overline{B}_\mu `$ has been performed by the CHORUS Collaboration by using a statistics of about 1000 charm events reconstructed in the nuclear emulsions. Out of these, ($`88\pm 10\pm 8`$) dimuon events have been reconstructed, which correspond to

$$\overline{B}_\mu =(9.3\pm 1.3)\%.$$

This measurement has to be compared with our prediction obtained by convoluting the charm-production fractions with the CHORUS neutrino flux

$$\overline{B}_\mu =(8.8\pm 1.0)\%.$$

In Table 8 we derived, by using the $`V_{cd}`$ value obtained by imposing the unitarity constraint to the CKM matrix and the measurements of $`B_\mu V_{cd}^2`$ from various experiments, $`B_\mu `$. Given the fact that the different experiments exploit different neutrino energy spectra, we can probe the sensitivity of $`B_\mu `$ to the neutrino energy. As expected, the higher the neutrino beam energy the larger the value of $`B_\mu `$.

## 4 Conclusions

We have presented a method to extract the charmed fractions in neutrino induced events. The method relies on the fact that, apart from processes peculiar to neutrinos such as quasi-elastic and diffractive charm-production, the charm-production and fragmentation mechanism is believed to be process-independent. We have verified this natural assumption going through a complete review of the available data from different experiments. Moreover, the $`D^+`$ over $`D^0`$ ratio is constant over the large energy range spanned by the collider and fixed target experiments reviewed in this paper. By using recent data from neutrino experiments, we have assessed the consistency of this ratio with the one predicted by other experiments. On the other side, the energy-independent behavior of the ratio itself is clear from the review of the experiments.

Threshold effects for the $`D_s`$ production are seen and accounted for. By introducing the diffractive charm-production and using the unity constraint, we have predicted the charm fragmentation as a function of the neutrino energy in the range useful to present neutrino experiments. In particular, a prediction for the CHORUS experiment has been made. The determination of the fragmentation function also allows the prediction of the semi-muonic branching ratio of charmed hadrons. The prediction given is in good agreement with a recent measurement made by the CHORUS experiment.