format: text

text:
# Effect of pinning on the vortex-lattice melting line in type-II superconductors

## I Introduction

In three-dimensional high-$`T_c`$ superconductors with pinning, two phase transition lines are known to exist in the magnetic field $`H`$ \- temperature $`T`$ plane:1 ; 2 ; 3 ; 4 ; 5  The line $`H_m(T)`$ where a quasiordered Bragg glass6 ; 7  thermally melts into a flux-line liquid, and the order–disorder transition line $`H_{dis}(T)`$ separating the Bragg glass from an amorphous vortex state. The melting is caused by thermal vibrations of the lattice, while the order–disorder transition is induced by quenched disorder in the vortex system. These two lines merge at some point in the $`H`$-$`T`$ plane. Although both transitions are accompanied by a proliferation of dislocations in the vortex lattice, it was argued 8  that the dislocation density $`\rho `$ is essentially different in these cases: $`\rho a^2`$ for melting, and $`\rho R_a^2`$ for the order–disorder transition. Here $`a`$ is the spacing between flux lines, and $`R_a`$ is the so-called positional correlation length 9  within which the relative vortex displacements caused by the quenched disorder are of the order of $`a`$. In fact, an intersection of these two different phase transition lines occurs in this scenario, and the order–disorder line terminates at the intersection point while the melting line continues for some distance to higher $`H`$, see Fig. 1. Within this physical picture, the existence of the so-called slush phase10  can be naturally explained. Recent experiments 11 ; 12 ; 13 ; 14 ; 15  for YBaCuO seem to support this scenario. On the other hand, experimental data 16 ; 17 ; 18  for BSCCO strongly argue in favor of a different scenario which was implied, e.g., in the Refs. 2, ; 5, . In this second scenario, the dislocation densities for both lines coincide at the point where these lines merge, and in fact, one deals with only one phase transition line that describes both the order–disorder transition at low temperature and the melting near $`T_c`$, Fig. 1.

Phase diagrams of superconductors with pinning reflect the competition of three characteristic energies 3 : the elastic energy, the pinning energy, and the energy of thermal fluctuations. At melting, the cost in the elastic energy due to the proliferation of dislocations is mainly balanced by the entropy gain associated with thermal fluctuations, while the role of the pinning energy, as shown below, is determined by the parameter $`a/R_c`$. Here $`R_c`$ is the transverse collective pinning length 9 . On the other hand, at the order–disorder transition, the balance of pinning energy and elastic energy is most important, while the relative contribution of the entropy gain is negligible at low temperatures and, according to the scenario of Ref. 8, , is determined by the ratio $`a/R_a`$ near the intersection point. Thus, if the intersection of the melting and the order–disorder lines occurred sufficiently deep in the bundle pinning region (so that $`R_ca`$ at this point), the scenario of Ref. 8,  would lead to the conclusion that up to the intersection point, one can find the melting line by neglecting the pinning, and the order–disorder line by neglecting the entropy gain. Just this approximation was used in our paper4  for analyzing the phase diagrams of superconductors. However, our recent results19  point out that flux-line pinning can affect the melting line near the intersection point since at this point, the ratio $`R_c/a`$ has decreased to several units for any magnitude of the quenched disorder in the vortex lattice (even when the disorder is weak). As to the second scenario, the three energies are all the same order of magnitude in the temperature region where the order–disorder transition gradually transforms into melting. Thus, whatever scenario occurs in reality, a detailed investigation of the effect of pinning on the melting line is important to clarify the most intriguing part of the phase diagram.

The effect of pinning by point defects on the melting line was observed both in BSCCO20 and in YBaCuO 21 ; 22 ; 23  crystals. It was discovered that an increase of the quenched disorder in the vortex lattice leads to a noticeable shift of the intersection point to lower magnetic fields and simultaneously pushes the melting line in the $`T`$-$`H`$ plane slightly downwards, i.e., at a fixed temperature the appropriate magnetic field of the melting, $`H_m(T)`$, decreases. It is important that the shift of the intersection point is essentially more pronounced than the decrease of $`H_m(t)`$ itself.

Some theoretical results on this subject were obtained in Refs. 24, ; 2, ; 5, . Larkin and Vinokur 24  started from the assumption that for the vortex lattice with quenched disorder to melt, the temperature must match a characteristic barrier composed of the elastic energy, $`E_{\mathrm{el}}`$, and the pinning energy, $`E_{\mathrm{pin}}`$. So they estimated the effect of pinning on the melting line by considering the following balance of these three energies:

$$T=E_{\mathrm{pin}}+E_{\mathrm{el}},$$
(1)
where $`E_{\mathrm{pin}}`$ and $`E_{\mathrm{el}}`$ were calculated in the so-called cage model.1  But it follows from this equation that the melting line has to shift upwards when the quenched disorder increases. Another approach was used in Ref. 2, . To describe the unified phase transition line, Giamarchi and Le Doussal, 2  who implied the second scenario, put forward the following generalization of the Lindemann criterion:

$$u_{\mathrm{total}}^2=c_L^2a^2,$$
(2)
where $`c_L`$ is the Lindemann constant ($`c_L0.10.2`$), and $`u_{\mathrm{total}}[2u_T^2+u^2(a,0)]^{1/2}`$ is the rms displacement of neighboring flux lines caused both by the thermal fluctuations and by the quenched disorder in the lattice. Here $`u_T`$ is the magnitude of the thermal fluctuations, while $`u(a,0)`$ describes the mean relative displacement of neighboring flux lines caused by the disorder. This criterion leads to the usual condition for the order–disorder transition 1 ; 3 ; 4  at low temperatures when $`u_T`$ is negligible, and it goes over to the well known Lindemann criterion for pure melting when the disorder disappears. Equation (2) results in a qualitatively correct dependence of the melting line on pinning by point defects. However, this reasonable criterion is only an interpolation formula between the two limiting cases and has no serious justification. The same is true for the criterion of Ref. 5, . Radzyner et al.5  described the unified phase transition line using the criterion,

$$u_T^2+u^2(a,0)=c_L^2a^2,$$
(3)
which practically coincides with Eq. (2). The criterion (3) is equivalent to the following balance of energies:

$$T+E_{\mathrm{pin}}=E_{\mathrm{el}},$$
(4)
which evidently differs from Eq. (1). In calculations of the phase transition line, the energies $`E_{\mathrm{pin}}`$ and $`E_{\mathrm{el}}`$ were estimated in Ref. 5,  in the framework of the cage model. But it remained unclear why the pinning energy $`E_{\mathrm{pin}}`$ now enters into Eq. (4) with the opposite sign as compared to the energy balance (2) suggested by Larkin and Vinokur. 24

In the present paper, in order to find and to justify a criterion for the melting of the flux-line lattice with quenched disorder, we start with the analysis of melting in the ideal lattice and show that three different approaches lead to the same dependence $`H_m(T)`$. These approaches are: the Lindemann criterion, the energy balance, and the rigorous approach based on the Ginzburg-Landau (or on the London) Hamiltonian. In the case of the lattice with quenched disorder, we show that the pinning energy in the flux-line liquid is larger than the pinning energy in the Bragg glass. For this reason, the difference of these pinning energies, $`E_{\mathrm{pin}}`$, has the opposite sign as compared to Eq. (1), and we arrive at an equation similar to Eq. (4) but with an expression for $`E_{\mathrm{pin}}`$ that differs from the estimates previously published. Besides this, we estimate the effect of pinning on the entropy gain at melting and find the upper critical point of the melting line from the condition that this gain vanishes. We also show that the result for $`H_m(T)`$ based on this energy balance agrees with the result which can be derived from the Ginzburg-Landau Hamiltonian. Then, using the $`H_m(T)`$ obtained in the framework of the second and the third approaches, we find how the Lindemann criterion should be modified to give the same melting line. Interestingly, the presented energy balance clarifies the difference between the first and the second scenarios. Finally, we present the $`T`$ \- $`H`$ phase diagrams of superconductors with two types of flux-line pinning by point defects and compare these diagrams for the two scenarios.

In this paper we consider only magnetic fields exceeding considerably the lower critical field $`H_{c1}`$ and thus do not distinguish between the magnetic field $`H`$ and the magnetic induction $`B`$. Besides this, we deal only with uniaxial anisotropic three-dimensional superconductors, neglecting completely the decoupling of the superconducting layers. The anisotropy is characterized by the parameter $`ϵ=\lambda _{ab}/\lambda _c<1`$ where $`\lambda _{ab}`$ and $`\lambda _c`$ are the London penetration depth in the plane $`ab`$ perpendicular to the anisotropy axis and along this axis, respectively. The magnetic field is assumed to be directed along the anisotropy axis. As to the quenched disorder in the flux-line lattice, we assume that it is caused by point defects and is not too strong such that the melting line lies entirely in the bundle pinning region.

## II melting of the ideal vortex lattice

We begin with the analysis of melting in the ideal pinning-free vortex lattice and compare the results of various approaches.

### II.1 Lindemann criterion

According to the well known Lindemann criterion, the flux-line lattice melts when the magnitude of the thermal displacements of the lattice relative to its equilibrium position, $`u_T`$, reaches a certain fraction of the spacing between the flux lines, $`a`$:

$$u_T^2=c_L^2a^2,$$
(5)
where $`a=(\mathrm{\Phi }_0/H)^{1/2}`$, $`\mathrm{\Phi }_0`$ is the flux quantum, and $`c_L`$ is the Lindemann constant. The magnitude $`u_T`$ depends on the elastic moduli of the lattice 25  and was calculated in many papers; see, e.g., Refs. 9, , 25, ; 26, ; 27, ; 28, . It can be represented in the form:4

$$u_T^2\xi ^2t\left(\frac{Gi}{1t^2}\right)^{1/2}h^{1/2}f(h),$$
(6)
where $`\xi (t)`$ is the coherence length in the $`ab`$ plane, $`h=H/H_{c2}(t)`$, $`t=T/T_c`$, $`H_{c2}(t)=\mathrm{\Phi }_0/2\pi \xi ^2`$ is the upper critical field, $`Gi`$ is the Ginzburg number,

$$Gi=\frac{1}{2}\left(\frac{T_c}{H_c^2ϵ\xi _0^3}\right)^2,$$

which characterizes the width of the fluctuation region in zero magnetic field, $`\xi _0`$ and $`H_c`$ are the coherence length and thermodynamic critical magnetic field of the superconductor in the Ginzburg - Landau theory extrapolated to $`T=0`$. For definiteness, we implied in Eq. (6) and below that $`\xi ^2(t)=\xi ^2(0)/(1t^2)`$. \[Hence $`\xi _0=\xi (0)/\sqrt{2}`$\]. The complete expression for the function $`f(h)`$ was given in Ref. 28, , but for our further analysis it is sufficient to use a simplified form of this function 4  in which the contribution containing the compression modulus of the vortex lattice, $`c_{11}`$, is neglected:

$$f(h)=\frac{2\beta _A}{1h}\frac{[1+(1+\stackrel{~}{c})^2]^{1/2}1}{\stackrel{~}{c}(1+\stackrel{~}{c})},$$
(7)
with $`\stackrel{~}{c}=0.5[\beta _A(1h)]^{1/2}`$, and $`\beta _A=1.16`$. Note that this formula can be rewritten in the form:

$$f(h)=\frac{f_1(h)}{(1h)^{3/2}},$$
(8)
where the function $`f_1(h)`$ defined by this equality decreases monotonically with increasing $`h`$, and its variation in the interval $`0<h<1`$ is not large: $`f_1(0)2.34`$, $`f_1(1)1.78`$. Thus, to a first approximation, this function can be considered as a constant, $`f_1(h)f_12`$.

Combining formulas (5)-(8), we arrive at the equation for the normalized melting field $`h_m(t)=H_m(t)/H_{c2}(t)`$:

$$t\left(\frac{Gi}{1t^2}\right)^{1/2}h_m^{1/2}\frac{f_1(h_m)}{(1h_m)^{3/2}}=2\pi c_L^2.$$
(9)
This equation agrees with those obtained earlier9 ; 25 ; 26 ; 27 ; 28  and differs from them only in the form of the function $`f_1(h)`$ since different authors used slightly different approximations for the elastic moduli or took into account the contribution associated with the compression modulus of the vortex lattice.

When the normalized melting field $`h_m(t)`$ is small, $`h_m1`$, it follows from Eq. (9) that

$$\frac{H_m(t)}{H_{c2}(0)}=\left(\frac{2\pi c_L^2}{f_1(0)t}\right)^2\frac{(1t^2)^2}{Gi}.$$
(10)
This is a well-known result 9 ; 25 ; 26 ; 27 ; 28  that holds for temperatures near $`T_c`$ (but outside the fluctuation region in zero magnetic field), namely for

$$Gi1t^2(f_1(0)/2\pi c_L^2)^2Gi.$$
(11)
Note that for inequalities (11) to be fulfilled in a sufficiently wide temperature interval, the Ginzburg number should not be too small. In the opposite limiting case,

$$1t^2(f_1(1)/2\pi c_L^2)^2Gi,$$
(12)
the field $`h_m`$ is close to unity, and one obtains from Eq. (9):

$$1h_m\left(\frac{f_1(1)}{2\pi c_L^2}\right)^{2/3}t^{2/3}\left(\frac{Gi}{1t^2}\right)^{1/3},$$
(13)
or equivalently,

$$\frac{H_{c2}(t)H_m(t)}{H_{c2}(0)}\left(\frac{f_1(1)}{2\pi c_L^2}\right)^{2/3}t^{2/3}Gi^{1/3}(1t^2)^{2/3}.$$
(14)
Note that $`H_{c2}(0)Gi^{1/3}(1t^2)^{2/3}`$ is the width (along the $`H`$ axis) of the fluctuation region in not too small magnetic fields, 29 ; 30 ; 31  $`HGiH_{c2}(0)`$, and for inequality (12) to hold, the Ginzburg number should not be too large.

### II.2 Energy balance

At melting, proliferation of dislocations occurs in the vortex lattice. These dislocations create a network in the lattice, and we consider a mean unit cell of this network composed of edge and screw dislocations. The energies of these dislocations are of the order of

$`E_{\mathrm{edge}}`$ $``$ $`c_{66}a^2l_{},`$
$`E_{\mathrm{screw}}`$ $``$ $`(c_{44}c_{66})^{1/2}a^2l_{},`$ (15)
where $`c_{66}`$ and $`c_{44}`$ are the shear and tilt moduli of the flux-line lattice, while $`l_{}`$ and $`l_{}`$ are the dimensions of the cell in the longitudinal and transverse directions to $`H`$, respectively. Since at the melting these dimensions are of the order of $`a`$ (see below), we have omitted the logarithmic factors $`\mathrm{ln}(l_{}/a)`$, $`\mathrm{ln}(l_{}/a)`$ in the above formulas for the dislocation energies. The shear and tilt moduli may be expressed as,25 ; 9

$`c_{66}{\displaystyle \frac{\epsilon _0}{4a^2}}(1h)^2,`$
$`c_{44}{\displaystyle \frac{ϵ^2\epsilon _0}{a^2}}(1h),`$ (16)
where $`h=H/H_{c2}(t)`$, $`\epsilon _0=(\mathrm{\Phi }_0/4\pi \lambda _{ab})^2`$, and we have omitted the logarithmic factor of type $`\mathrm{ln}(a/\xi )`$ in the tilt modulus (this nonlocal modulus should be estimated at wave vectors $`k`$ of the order of $`l_{}^1`$, $`l_{}^1`$). The factors containing $`(1h)`$ take into account the softening of the vortex lattice near the $`H_{c2}(t)`$ line.32  Minimization of the elastic energy of the cell at its fixed volume leads to $`E_{\mathrm{edge}}E_{\mathrm{screw}}`$ and hence gives the relation between $`l_{}`$ and $`l_{}`$:

$$l_{}l_{}\left(\frac{c_{44}}{c_{66}}\right)^{1/2}.$$
(17)
Thus, up to a numerical factor, the elastic energy of the dislocation cell, $`E_{\mathrm{el}}`$, equals $`E_{\mathrm{screw}}`$.

At melting, the cost in the elastic energy due to the proliferation of the dislocation network is balanced by the entropy gain in the free energy of the flux-line lattice. We now give a simple estimate of this gain: The vortex-lattice degrees of freedom associated with shear undergo a change when the melting occurs. This change occurs for lattice modes with wave-lengths greater than $`l_{}`$, $`l_{}`$. There are $`l_{}^1l_{}^2`$ modes of this type in the unit volume of the lattice, and each of them contributes about $`T`$ to the entropy gain. Thus, up to a numerical factor, the gain per cell of the dislocation network is $`Tl_{}^1l_{}^2l_{}l_{}^2T`$.

We now can write down the change of the free energy per cell at melting:

$$\mathrm{\Delta }FC(c_{44}c_{66})^{1/2}a^2l_{}T,$$
(18)
where some constant $`C`$ is the ratio of the unknown numerical factors in $`E_{\mathrm{el}}`$ and in the entropy gain. Minimization of Eq. (18) with respect to the parameter $`l_{}`$ leads to the conclusion that this parameter should have the minimum possible value. It is clear that this value is of the order of $`a`$ in the lattice. Then, taking into account that $`\mathrm{\Delta }F=0`$ at the melting, we arrive at

$$C(c_{44}c_{66})^{1/2}a^3T=0,$$
(19)
where this constant $`C`$ may be slightly renormalized as compared to the $`C`$ of Eq. (18). Inserting Eqs. (II.2) for the elastic moduli and assuming $`\lambda _{ab}(t)/\xi (t)=`$const, one finds the equation for the melting line $`h_m(t)`$:

$$t\left(\frac{Gi}{1t^2}\right)^{1/2}\frac{h_m^{1/2}}{(1h_m)^{3/2}}=\frac{\sqrt{\pi }C}{4}.$$
(20)
As was noted earlier (see, e.g., Refs. 3, , 9, ), formula (10) for the melting line near $`T_c`$ can be obtained from the energy balance. Here we have taken into proper account the softening of the elastic moduli near $`H_{c2}(T)`$, and now equation (20) shows that not only formula (10) but also expression (13) can be derived by this method. Moreover, if one uses the approximation in which the function $`f_1(h)`$ is a constant,33  $`f_1(h)f_1`$, and put $`C=8\sqrt{\pi }c_L^2/f_1`$, Eqs. (9) and (20) completely coincide in the whole temperature interval.

### II.3 Some exact results

Within the mean-field theory, when one neglects fluctuations of the superconducting order parameter, the melting line $`H_m(T)`$ coincides with the $`H_{c2}(t)`$ line. It is the fluctuations that shift $`H_m(t)`$ downwards in the $`H`$-$`T`$ plane. As was mentioned above, see Eq. (14), the Lindemann criterion shows that at sufficiently strong magnetic fields, the distance between the melting line and the mean-field $`H_{c2}(t)`$ line is comparable with the width of the fluctuation region. But then the question arises about the applicability of this criterion (and of the energy balance) for determining $`H_m`$ in this region of the magnetic fields since expressions (II.2) for the elastic moduli were derived in the framework of the mean-field theory without accounting for the fluctuations. In this context, it should be noted that one cannot confine oneself to taking into account only the first fluctuation correction to the elastic moduli (in the fluctuation amplitude) since inside the fluctuation region the amplitude is large, and corrections of all orders are essential. In particular, the renormalized elastic moduli will vanish on a line which differs from the mean-field $`H_{c2}(t)`$. However, simple considerations 30 ; 34  show that at $`1h_m1`$, strong fluctuations can only renormalize the numerical factor in Eq. (13), but the dependences of $`H_m`$ on $`t`$ and on $`Gi`$ remain unchanged. We now briefly outline these considerations.

In dimensionless units the Ginzburg-Landau Hamiltonian depends on the three parameters: $`t`$, $`H/H_{c2}(0)`$ and the Ginzburg number $`Gi`$. As well known, the quadratic part (in the order parameter) of this Hamiltonian looks like the Hamiltonian of a particle with double electron charge in a magnetic field, and the energy spectrum of this particle is the so-called Landau levels. It is essential that in fields $`HGiH_{c2}(0)`$, the distance between these levels exceeds the width of the fluctuation region. Thus, if one expands the order parameter into the eigenfunctions of the particle, only the modes of the order parameter corresponding to the lowest Landau level strongly fluctuate near the melting line, and one may retain only these modes in the Hamiltonian to calculate the fluctuation part of the free energy of a superconductor. It turns out that the Hamiltonian thus obtained depends on a single combination of the parameters. In our notations, this combination, $`Q`$, can be represented in the form:

$$Q=\frac{(1h)(1t^2)^{1/3}}{t^{2/3}Gi^{1/3}h^{2/3}},$$
(21)
where $`h=H/H_{c2}(t)`$. Hence, in this region of the magnetic fields, the free energies of the vortex liquid, $`F_{\mathrm{liq}}`$, and the vortex lattice, $`F_{\mathrm{lat}}`$, are also determined only by this combination, and for the melting line $`h_m(t)`$ we arrive at the equation:

$$F_{\mathrm{liq}}(Q)=F_{\mathrm{lat}}(Q).$$
(22)
The solution of this equation has the form:

$$Q=\frac{(1h_m)(1t^2)^{1/3}}{t^{2/3}Gi^{1/3}h_m^{2/3}}=C_1,$$
(23)
where $`C_1`$ is some constant. Taking into account that $`h_m1`$, we find from equation (23):

$$1h_mC_1t^{2/3}\left(\frac{Gi}{1t^2}\right)^{1/3}.$$
(24)
It is seen that up to a numerical factor, this expression indeed coincides with Eq. (13).

The above consideration does not use any perturbation theory and therefore is exact, but it does not yield the constant $`C_1`$. Using some variants of perturbation theory (see, e.g., Ref. 35, ), approximate expressions for $`F_{\mathrm{liq}}`$ and $`F_{\mathrm{lat}}`$ were obtained in Refs. 36, , 37, . On this basis, an equation for the melting line, which agrees with Eq. (23), was derived in these papers together with the appropriate constant $`C_1`$. Hikami et al.36  estimated $`C_17`$, while Li and Rosenstein 37  who used a refined expression for $`F_{\mathrm{lat}}`$ found $`C_19.5`$.

A consideration similar to that presented above was also applied to the Hamiltonian of the vortex system in the London approximation, see Ref. 38, . In the region of the magnetic fields considered here, $`HH_{c1}`$, one has $`\lambda _{ab}a`$, and the Hamiltonian, as well as the free energy of the vortex system, are determined by a single combination of the physical parameters: $`ϵ\epsilon _0a/T`$. Thus, the melting line $`H_m(T)`$ is found from the equation:

$$\frac{ϵ\epsilon _0(T)a(H_m)}{T}=C_2,$$
(25)
with some constant $`C_2`$. In other words, the London approximation leads to $`H_m(t)[t\lambda _{ab}^2(t)]^2`$. Putting $`\lambda _{ab}(t)=\lambda _{ab}(0)(1t^2)^{1/2}`$ and taking into account the definition of $`Gi`$, we arrive at the formula:

$$\frac{H_m(t)}{H_{c2}(0)}=\frac{\pi }{(2C_2t)^2}\frac{(1t^2)^2}{Gi},$$
(26)
which agrees with Eq. (10).

To summarize, we have shown in this section that up to numerical factors, the three different approaches lead to the same dependences of the vortex-lattice melting field $`H_m`$ on the temperature $`T`$ and on the Ginzburg number $`Gi`$, Eqs. (10) and (14). Thus, after obtaining these dependences, e.g., from the energy balance, for the vortex lattice with pinning, one can guess the true form of the Lindemann criterion for this lattice. We shall use this procedure in the next section.

## III melting of the vortex lattice with quenched disorder

In this section we analyze the melting line of the vortex lattice with pinning assuming that the line is in the bundle pinning region (where the transverse collective pinning length $`R_c`$ is greater than $`a`$). As it follows from the figures of Sec. IV B, the melting line, as a rule, does entirely lie in this region, and only in the case of sufficiently strong $`\delta T_c`$ pinning can it enter the single vortex pinning region.

### III.1 Energy balance

We now consider the influence of pinning by point defects on the energy balance. The adjustment of the vortex system to the pinning potential decreases the total energy of the system, and the amount of this decrease is just of the order of the pinning energy. It will be shown in this section that at melting, the pinning energy in the flux-line liquid, $`E_{\mathrm{pin}}^{\mathrm{liq}}`$, is noticeably greater than the pinning energy in the flux-line lattice (Bragg glass), $`E_{\mathrm{pin}}^{\mathrm{lat}}`$. Thus, there is a gain in the free energy of the vortex system, $`E_{\mathrm{pin}}`$, associated with the pinning energy: $`E_{\mathrm{pin}}=E_{\mathrm{pin}}^{\mathrm{liq}}E_{\mathrm{pin}}^{\mathrm{lat}}E_{\mathrm{pin}}^{\mathrm{liq}}`$. This gain adds to the entropy gain, and equation (18) is modified as follows:

$$\mathrm{\Delta }FC(c_{44}c_{66})^{1/2}a^2l_{}T\mathrm{\Delta }SE_{\mathrm{pin}},$$
(27)
where the factor $`\mathrm{\Delta }S`$ takes into account the effect of pinning on the entropy gain per dislocation cell.

In order to estimate $`E_{\mathrm{pin}}^{\mathrm{lat}}`$, $`E_{\mathrm{pin}}^{\mathrm{liq}}`$ and $`\mathrm{\Delta }S`$, it is necessary in general to take into account the so-called thermal depinning9 ; 39  since the thermal displacement $`u_T`$ is sufficiently large at the melting, $`u_Tc_La`$. However, to explain the main ideas, we first carry out the estimates neglecting this depinning and then generalize the obtained results by taking it into account.

The pinning energy of the flux-line lattice without dislocations, i.e., of the Bragg glass, in the volume equal to the cell of the dislocation network, $`E_{\mathrm{pin}}^{\mathrm{lat}}`$, can be estimated using the results of collective pinning theory:9

$$E_{\mathrm{pin}}^{\mathrm{lat}}(Wl_{}l_{}^2)^{1/2}u,$$
(28)
where $`W=f_{\mathrm{pin}}^2n\xi ^2/a^2`$, $`f_{\mathrm{pin}}`$ is the mean elementary pinning force exerted by one point defect, $`n`$ is the concentration of the defects, and the expression $`(Wl_{}l_{}^2)^{1/2}`$ is the mean pinning force per cell. The displacement $`uu(l_{},0)u(0,l_{})`$ is the rms relative shift of two line elements in the vortex lattice separated by a distance $`l_{}`$ transverse to the magnetic field, or by $`l_{}`$ along the field. This shift is caused by the random point defects. The magnitude of $`u`$ can be expressed in terms of the transverse collective pinning length, $`R_c`$, at which the relative displacements of points in the lattice are of the order of $`\xi `$. In particular, if the small bundle pinning regime occurs, i.e., $`R_c<\lambda _{ab}`$, one has:9

$$u^2\frac{\xi ^2}{1+\mathrm{ln}(R_c^2/a^2)},$$
(29)
where we have used that $`l_{}a`$. Note that according to Eq. (17), the longitudinal dimension $`l_{}`$ is of the order of $`ϵa/(1h)^{1/2}`$.

At melting, proliferation of dislocations occurs, and in the liquid vortices can adjust themselves to the pinning potential not only via their elastic deformations as in the Bragg glass but also via the plastic vortex-lattice displacements generated by dislocations. This additional adjustment mechanism increases the pinning energy in the liquid. To estimate the pinning energy per dislocation cell in the liquid, $`E_{\mathrm{pin}}^{\mathrm{liq}}`$, it is necessary to take into account that the displacements generated by the dislocation network are essentially larger than the displacements existing in the Bragg glass at the same temperature and magnetic field within the scales $`l_{}`$, $`l_{}`$. Indeed, in the lattice without dislocations, one has $`u(l_{},0)u(0,l_{})=u<\xi `$ \[see Eq. (29)\], while in the liquid the dislocations lead to displacements of about $`a>\xi `$ at the boundary of the dislocation cell. It is also essential that the displacements caused by the dislocations within a scale $`l<l_{}`$, $`l_{}`$ have a different dependence on $`l`$ than the elastic displacements in the Bragg glass. Then, we obtain the estimate:

$$E_{\mathrm{pin}}^{\mathrm{liq}}(Wl_{}l_{}^2)^{1/2}\xi .$$
(30)
Note that although the displacements generated by the dislocations are large and exceed $`\xi `$, we multiply the mean pinning force per cell only by $`\xi `$ in formula (30) since $`\xi `$ is the effective range of the elementary pinning force $`f_{\mathrm{pin}}`$. It follows from Eqs. (28)-(30) that $`E_{\mathrm{pin}}^{\mathrm{liq}}/E_{\mathrm{pin}}^{\mathrm{lat}}(\xi /u)[1+\mathrm{ln}(R_c^2/a^2)]^{1/2}`$, i.e., $`E_{\mathrm{pin}}^{\mathrm{liq}}`$ noticeably exceeds $`E_{\mathrm{pin}}^{\mathrm{lat}}`$ in the bundle pinning regime, and hence formula (30) gives the estimate of $`E_{\mathrm{pin}}`$ in Eq. (27).

It is convenient to rewrite expression (30) in terms of $`L_c`$, the single-vortex collective pinning length, using the relation, 9  $`f_{\mathrm{pin}}^2nϵ^4\epsilon _0^2/L_c^3`$. Then one arrives at the formula for $`E_{\mathrm{pin}}`$:

$$E_{\mathrm{pin}}ϵ\epsilon _0a[Dg_0(t)]^{3/2}h^{1/4}(1h)^{3/4},$$
(31)
where we have inserted the estimates for $`l_{}`$, $`l_{}`$, have taken into account that the quantity $`\epsilon _0\lambda _{ab}^2`$ should lead to an additional factor $`1h`$ when $`h=H/H_{c2}(t)`$ tends to unity,25  and have used the notation:4

$$\frac{ϵ\xi (t)}{L_c(t)}Dg_0(t)$$
(32)
with $`ϵ\xi (0)/L_c(0)=D`$. The function $`g_0(t)`$ is given4  by

$$g_0(t)=(1t^2)^{1/2}$$
(33)
for $`\delta l`$ pinning, and by

$$g_0(t)=(1t^2)^{1/6}$$
(34)
for $`\delta T_c`$ pinning. The parameter $`D`$ is a measure of the pinning strength and is estimated as 9  $`D(j_c/j_0)^{1/2}`$ where $`j_c`$ is the critical current density in the single vortex pinning regime and $`j_0`$ is the depairing current density, both taken at $`T=0`$.

Let us now estimate the entropy term in the energy balance (27). In Sec. II B we obtained the expressions for the melting line of the ideal vortex lattice, assuming that the main contribution to the entropy gain at melting is due to disappearance of the shear phonon modes in the liquid. The mechanism of this disappearance is the following: A shear stress in the liquid generates the so-called Peach-Köhler forces 40 ; 41  exerted on dislocations; the dislocations begin to move, and their shifts relax the shear stress in the liquid. In the vortex system with pinning a dislocation cannot move if the pinning force per dislocation cell, $`(Wl_{}l_{}^2)^{1/2}`$, exceeds the Peach-Köhler force $`ac_{66}u_{xy}l_{}`$ exerted on a dislocation segment of length $`l_{}`$ where $`u_{xy}`$ is the shear deformation of the vortex system in the plane normal to the magnetic field. Thus, for relaxation of the shear stress to occur, the deformation $`u_{xy}`$ must exceed the critical value

$$u_{xy}^{cr}=\frac{(Wl_{}l_{}^2)^{1/2}}{ac_{66}l_{}}[Dg_0(t)]^{3/2}h^{1/4}(1h)^{3/4}.$$
(35)
Here we again have inserted the relations used in deriving Eq. (31). In the volume with dimensions $`L_{}`$ and $`L_{}`$, thermal fluctuations generate shear deformations $`u_{xy}`$ that can be estimated from

$$c_{66}u_{xy}^2L_{}L_{}^2T.$$
(36)
Hence, for large $`L_{}`$ and $`L_{}`$, when

$$L_{}L_{}^2>[L_{}L_{}^2]^{cr}\frac{T}{c_{66}(u_{xy}^{cr})^2},$$

the deformations are less than the critical value; the relaxation of the shear stress in the liquid does not occur, and these modes do not contribute to the entropy gain (they do not differ essentially from the appropriate modes in the Bragg glass). Then, applying the considerations of Sec. II B, we obtain the factor $`\mathrm{\Delta }S`$ in Eq. (27): $`\mathrm{\Delta }S=1P(t,h)`$, where

$$P(t,h)=\frac{l_{}l_{}^2}{[L_{}L_{}^2]^{cr}}\frac{[Dg_0(t)]^3(1t^2)^{1/2}}{tGi^{1/2}h}.$$
(37)
As is well known,9 ; 39  thermal fluctuations of the flux-line lattice lead to a smoothing of the pinning potential and thereby affect the pinning. We take into account the effect of this thermal depinning on $`E_{\mathrm{pin}}`$ and $`P(t,h)`$ using the recipe of Ref. 9,  (see page 1214 in that paper). Then an additional factor $`\xi /r_p`$ appears in Eq. (30), while the right hand side of Eq. (35) is multiplied by $`(\xi /r_p)^2`$ where $`r_p=(\xi ^2+u_T^2)^{1/2}`$ is the new effective range of the pinning force when thermal fluctuations are allowed for. This modifies formulas (31) and (37) as follows:

$`E_{\mathrm{pin}}ϵ\epsilon _0a{\displaystyle \frac{[Dg_0(t)]^{3/2}h^{1/4}(1h)^{3/4}}{\left[1+t\left({\displaystyle \frac{Gi}{1t^2}}\right)^{1/2}{\displaystyle \frac{f(h)}{h^{1/2}}}\right]^{1/2}}},`$
$`P(t,h){\displaystyle \frac{[Dg_0(t)]^3(1t^2)^{1/2}}{tGi^{1/2}h\left[1+t\left({\displaystyle \frac{Gi}{1t^2}}\right)^{1/2}{\displaystyle \frac{f(h)}{h^{1/2}}}\right]^2}},`$ (38)
where we have used Eq. (6) for $`u_T^2`$.

Inserting formulas (III.1) in Eq. (27), expressing $`C`$ as $`8\sqrt{\pi }c_L^2/f_1`$ with a constant $`f_1`$, and taking into account that $`\mathrm{\Delta }F=0`$ at melting, we eventually find the equation for the melting line $`h_m(t)`$:

$`{\displaystyle \frac{F_T(t)h_m^{1/2}}{(1h_m)^{3/2}}}[1P(t,h_m)]`$ $`+`$ (39)
$`{\displaystyle \frac{A[Dg_0(t)]^{3/2}h_m^{1/4}}{\left[(1h_m)^{3/2}+F_T(t)h_m^{1/2}\right]^{1/2}}}`$ $`=`$ $`2\pi c_L^2,`$
where

$`P(t,h)={\displaystyle \frac{B[Dg_0(t)]^3(1t^2)^{1/2}}{tGi^{1/2}\left[h^{1/2}+F_T(t)(1h)^{3/2}\right]^2}},`$ (40)
$`A`$, $`B`$ are some numerical factors, and we have used the notation

$$F_T(t)=f_1t\left(\frac{Gi}{1t^2}\right)^{1/2},$$
(41)
putting $`f_1=2`$ below. Equation (39) generalizes Eq. (20) obtained in the case of the ideal lattice. Note that the factor $`1P=\mathrm{\Delta }S`$ in Eq. (39) naturally explains the existence of the upper critical point on the melting line. The temperature $`t_{\mathrm{up}}`$ of this point can be obtained from the condition $`\mathrm{\Delta }S=0`$, or

$$P(t_{\mathrm{up}},h_m(t_{\mathrm{up}}))=1,$$
(42)
i.e., the entropy jump in the vortex system at melting vanishes at this point.

Equation (39) is in agreement with the general conclusion made in Ref. 4, . Namely, it was stated4  that the phase diagrams of various three dimensional superconductors with point defects are determined only by the Ginzburg number $`Gi`$, the parameter $`D`$ characterizing the strength of pinning, and the function $`g_0(t)`$ which is defined by the type of pinning. Other physical quantities \[like $`T_c`$, $`H_{c2}(T)`$, $`ϵ`$, etc.\] either determine only scaling factors or do not appear explicitly in the appropriate equations at all. In particular, the anisotropy $`ϵ`$ is absorbed by the definitions of $`Gi`$ and $`D`$ and does not enter in Eq. (39).

### III.2 Analysis of the equation

Figures 2 and 3 present the melting lines $`H_m(t)`$ obtained numerically from Eq. (39) for different pinning strength $`D`$. Although the numerical factors $`A`$ and $`B`$ have remained unknown in the above derivation of Eq. (39), the figures show that it is possible to choose $`A_1=A/(2\pi )^{3/4}`$ and $`B_1=Bf_1/(2\pi )^{3/2}`$ so that the calculated melting lines have the properties observed in experiments (see Introduction). In particular, we find that the upper critical point noticeably shifts under the influence of pinning and tends to $`T_c`$ at reasonable values of $`D`$ when $`A_1`$ and $`B_1`$ are not too small ($`A_1`$, $`\sqrt{B}_11`$). With increasing $`D`$ the melting line shifts downwards if $`B_1<A_1^2`$, and this shift becomes small when $`B_1`$ does not differ too much from $`A_1^2`$. Of course, the choice of $`A_1`$ and $`B_1`$ on the basis of these requirements is not unique. Let us now analyze Eq. (39) in some limiting cases.

According to Sec. II, in the temperature region (12) one has $`1h_m1`$ for the melting line of the ideal vortex lattice. If the parameter $`D`$ is not too large, this property of the melting line remains true for the vortex lattice with pinning. But under the condition $`1h_m1`$ the lowest Landau level approximation is valid, and we can use the results of Appendix A where the effect of quenched disorder on the melting line has been analyzed in the framework of an approach similar to that of Sec. II C. On the other hand, putting $`h_m1`$ in Eq. (39), we arrive at

$`Q^{3/2}\left(1{\displaystyle \frac{Bf_1[G(t)]^2}{[1+Q^{3/2}]^2}}\right)`$ $`+`$
$`{\displaystyle \frac{AG(t)Q^{3/4}}{[1+Q^{3/2}]^{1/2}}}`$ $`=`$ $`2\pi c_L^2,`$ (43)
where $`Q^{3/2}=(1h_m)^{3/2}/F_T(t)`$ \[compare with formula (21)\], and

$`G(t)={\displaystyle \frac{[Dg_0(t)]^{3/2}}{[F_T(t)]^{1/2}}}={\displaystyle \frac{[Dg_0(t)]^{3/2}(1t^2)^{1/4}}{(f_1t)^{1/2}Gi^{1/4}}}.`$ (44)
Equation (III.2) shows that the quantity $`Q`$ on the melting line depends on $`t`$ only via the function $`G(t)`$. This conclusion completely agrees with the exact result of Appendix A.

When $`D=0`$ \[i.e., $`G(t)=0`$\], Eq. (III.2) gives the result for the ideal vortex lattice, $`Q^{3/2}=2\pi c_L^2`$, see Eq. (13). For nonzero $`D`$ the function $`G(t)`$ decreases with the temperature. In the case of $`\delta T_c`$ pinning one has $`G(t)t^{1/2}`$, while $`G(t)`$ is proportional to $`(1t^2)/\sqrt{t}`$ for $`\delta l`$ pinning. Hence in the equation for the melting line the role of terms associated with pinning increases with decreasing $`t`$. It follows from Eqs. (III.2), (42) that at the upper critical point the function $`G(t)`$ reaches a certain value, and condition (42) can be written in the form:

$`\sqrt{Bf_1}G(t_{\mathrm{up}})={\displaystyle \frac{1}{2}}+\left({\displaystyle \frac{1}{4}}+(2\pi c_L^2)^2{\displaystyle \frac{Bf_1}{A^2}}\right)^{1/2}.`$ (45)
Interestingly, up to a numerical factor (and the factor $`[1h_{sv}(t)]^{3/2}`$ which we do not take into account here 42 ) this condition reduces to Eq. (32) of Ref. 4,  for the intersection point of the melting line with the order–disorder line. In particular, it follows from Eqs. (44), (45) that $`t_{\mathrm{up}}`$, similarly to the temperature of the intersection point, depends on $`Gi`$ and $`D`$ only via the combination $`D^3/Gi^{1/2}`$ of these parameters.

The so-called depinning line9  along which $`u_T\xi `$, intersects the melting line of the ideal lattice, $`H_m^{\mathrm{id}}(t)`$, at $`1t^2(f_1^2/2\pi c_L^2)Gi`$, i.e., outside the temperature region defined by Eq. (12). In the region $`1t^2(f_1^2/2\pi c_L^2)Gi`$, the depinning is essential. In this depinning region, the terms $`T\mathrm{\Delta }S`$ and $`E_{\mathrm{pin}}`$ in the energy balance (27) become relatively small, and the melting line is close to that of the ideal lattice.

### III.3 Lindemann criterion

We now rewrite Eq. (39) in the form of the Lindemann criterion. As it follows from formulas (5) and (9), the factor $`F_T(t)h_m^{1/2}/(1h_m)^{3/2}`$ in the first term of Eq. (39) is simply $`2\pi u_T^2/a^2`$ where $`u_T`$ is described43  by Eq. (6). To rewrite the second term in Eq. (39), let us take into account the formula for the averaged relative shift of two line elements in the vortex lattice separated by a distance $`RR_c`$ transverse to the magnetic field:44

$$\frac{u^2(R,0)}{r_p^2}\left(\frac{ϵa}{L_c}\frac{\xi ^2}{r_p^2}\right)^3\frac{1+\mathrm{ln}(R^2/a^2)+ϵR/\lambda _{ab}}{(1h)^{3/2}}.$$
(46)
\[If the small bundle pinning regime occurs, $`R_c<\lambda _{ab}`$, the last term in this formula, $`ϵR/\lambda _{ab}`$, is small and may be omitted.\] At $`R=R_c`$, one has $`u(R_c,0)=r_p`$. Inserting $`r_p^2=\xi ^2+u_T^2`$ and the expression (6) for $`u_T`$ in the relationship thus obtained, and taking into account the definition (32), we arrive at

$`1+\mathrm{ln}\left({\displaystyle \frac{R_c^2}{a^2}}\right)+ϵ{\displaystyle \frac{R_c}{\lambda _{ab}}}`$ $`=`$ $`\left({\displaystyle \frac{h(1h)}{2\pi D^2g_0^2(t)}}\right)^{3/2}\times `$ (47)
$`[1`$ $`+`$ $`\left({\displaystyle \frac{Gi}{1t^2}}\right)^{1/2}{\displaystyle \frac{f_1t}{h^{1/2}(1h)^{3/2}}}]^3.`$
Simple manipulations using this formula show that the second term in Eq. (39), which tells the role of the pinning energy in the energy balance, is

$$2\pi A_1\frac{r_p^2}{a^2}\left[1+\mathrm{ln}\left(\frac{R_c^2}{a^2}\right)+\frac{ϵR_c}{\lambda _{ab}}\right]^{1/2},$$
(48)
where $`A_1=A/(2\pi )^{3/4}`$. Since the factor containing $`R_c`$ can be represented as $`u(a,0)/r_p`$ \[using Eq. (46) at $`R=a`$ and then at $`R=R_c`$\], we find one more form for the second term in Eq. (39):

$$2\pi A_1\frac{r_pu(a,0)}{a^2}.$$
(49)
As to the quantity $`P(t,h)`$ given by Eq. (40), it is expressed as $`B_1u^2(a,0)/u_T^2`$ where $`B_1=Bf_1/(2\pi )^{3/2}`$. Thus, equation (39) is equivalent to the following criterion:

$$u_T^2B_1u^2(a,0)+A_1r_pu(a,0)=c_L^2a^2,$$
(50)
which is true for the melting line in the bundle pinning region. Equation (50) can be also rewritten in the form:

$`u_T^2B_1u^2(a,0)+A_1u^2(a,0)\left[1+\mathrm{ln}\left({\displaystyle \frac{R_c^2}{a^2}}\right)+{\displaystyle \frac{ϵR_c}{\lambda _{ab}}}\right]^{1/2}`$
$`=c_L^2a^2.`$ (51)
Note that although criteria (2) and (3) are qualitatively close to Eq. (III.3), they underestimate pinning at large $`R_c`$.

Equation (50) \[or (III.3)\] is valid when $`u_T^2B_1u^2(a,0)`$ since condition (42) for the upper critical point now has the form:

$$u_T^2=B_1u^2(a,0).$$
(52)
Taking into account that $`B_11`$, this condition means that the upper critical point, as a rule, lies in the bundle pinning region. Indeed, if the melting line enters the single vortex pinning region before the upper critical point occurs, one finds $`u_T^2B_1r_p^2=B_1(\xi ^2+u_T^2)`$ at the boundary of this region. In other words, the melting line can intersect the boundary in the depinning region where $`u_T^2`$ is large. But in this region the melting line practically coincides with that of the ideal lattice, Sec. III B; the latter can cross the boundary only in the case of $`\delta T_c`$ pinning at sufficiently large $`D`$. 19

## IV Phase diagrams

### IV.1 The first and the second scenarios

To gain some insight into the character of the mergence of the melting line with the order–disorder line, let us consider more closely the dependence of the free energy difference between the liquid and the Bragg glass, $`\mathrm{\Delta }F`$, on the size of the dislocation cell, $`l_{}`$. We shall consider only such $`l_{}`$ that lie near the minimum possible size in the lattice $`l_{}^{\mathrm{min}}`$ ($`l_{}^{\mathrm{min}}a`$). In evaluating the entropy gain $`T\mathrm{\Delta }S`$, the elastic energy $`E_{\mathrm{el}}`$ and the pinning energy $`E_{\mathrm{pin}}`$, it is necessary to take into account that the modulus $`c_{44}`$ depends on the wave vector $`𝐤`$. According to Ref. 25, , one has $`c_{44}(k1/l_{})c_{44}(1/a)(l_{}/a)^2`$ at $`l_{}<\lambda _{ab}/ϵ`$, and hence $`l_{}l_{}^2`$, see formula (17). Then we can write $`\mathrm{\Delta }F`$ in the form (we still neglect logarithmic factors):

$$\mathrm{\Delta }FE_{\mathrm{el}}(a)\stackrel{~}{l}_{}^2E_{\mathrm{pin}}(a)\stackrel{~}{l}_{}^2+T[P(a)\stackrel{~}{l}_{}^41],$$
(53)
where $`\stackrel{~}{l}_{}=l_{}/l_{}^{\mathrm{min}}`$, $`P(a)`$ denotes $`P(t,h_m)`$ of Eq. (39), and $`E_{\mathrm{el}}(a)`$ and $`E_{\mathrm{pin}}(a)`$ are the elastic and the pinning energies per dislocation cell at $`l_{}=l_{}^{\mathrm{min}}`$, which have been inserted into formula (27) to derive Eq. (39). Equation (53) is valid when $`P(a)\stackrel{~}{l}_{}^41<0`$, otherwise the last term in this equation has to be omitted, and one has,

$$\mathrm{\Delta }FE_{\mathrm{el}}(a)\stackrel{~}{l}_{}^2E_{\mathrm{pin}}(a)\stackrel{~}{l}_{}^2,$$
(54)
for $`\stackrel{~}{l}_{}^4>1/P(a)`$.

As was mentioned above, in the equation (39) for $`H_m(t)`$ \[or equivalently $`T_m(H)`$\] the relative role of terms associated with pinning increases with decreasing $`T_m`$ for both types of pinning. When $`E_{\mathrm{pin}}(a)`$ and $`P(a)`$ are sufficiently small (i.e., for temperatures $`T_m`$ near $`T_c`$), the function $`\mathrm{\Delta }F(\stackrel{~}{l}_{})`$ is minimum at the lowest possible value $`\stackrel{~}{l}_{}=1`$ and then increases with increasing $`\stackrel{~}{l}_{}`$, Fig. 4a. This means that melting occurs at $`\stackrel{~}{l}_{}=1`$ as in the ideal lattice. This conclusion has been already used implicitly in deriving Eq. (39). When the temperature $`T_m`$ decreases, and hence $`P(a)`$ and the ratio $`E_{\mathrm{pin}}(a)/E_{\mathrm{el}}(a)`$ increase, the point $`\stackrel{~}{l}_{}=[1/P(a)]^{1/4}`$ shifts to smaller $`\stackrel{~}{l}_{}`$, and the slope of $`\mathrm{\Delta }F`$ in Eq. (54) on $`\stackrel{~}{l}_{}^2`$ decreases. Eventually, at the critical temperature $`t_{\mathrm{up}}`$, one arrives at the situation shown either in Fig. 4c or in Fig. 4d since the energy balance for the melting $`E_{\mathrm{el}}(a)E_{\mathrm{pin}}(a)=T_m[1P(a)]`$ leads to $`E_{\mathrm{el}}(a)E_{\mathrm{pin}}(a)=0`$ at this point. In the case of Fig. 4d when the dependence of $`\mathrm{\Delta }F`$ on $`\stackrel{~}{l}_{}^2`$ has a negative curvature at $`\stackrel{~}{l}_{}=1`$, we conclude that this case is necessarily preceded by the situation shown in Fig. 4b, which has to occur at some temperature $`t_i>t_{\mathrm{up}}`$ where $`P(a)<1`$ and $`E_{\mathrm{pin}}(a)<E_{\mathrm{el}}(a)`$. At this $`t_i`$ the size of the dislocation cell, $`l_{}`$, sharply increases because the absolute minimum of $`\mathrm{\Delta }F`$ now occurs at a larger $`\stackrel{~}{l}_{}>1`$. In other words, at this temperature we find an intersection of the melting line with the order–disorder line. The abrupt increase of $`l_{}`$ also means that the melting and the order–disorder transition cannot form a unified phase transition line. Thus, we arrive at the first scenario45  shown in Fig. 1. In the case of Fig. 4c two possibilities exist: Either at some temperature $`t_i>t_{\mathrm{up}}`$ the situation shown in Fig. 4b occurs, and we again arrive at the first scenario, or the function $`\mathrm{\Delta }F(\stackrel{~}{l}_{})\mathrm{\Delta }F(1)`$ remains positive for all $`\stackrel{~}{l}_{}`$ in the temperature interval $`t_{\mathrm{up}}t1`$. In the latter case we find that the mean size of the dislocation cell, $`l_{}`$, cannot gradually change in the process of the reduction of the entropy gain, and hence one has $`\stackrel{~}{l}_{}=1`$ down to $`t_{\mathrm{up}}`$. In this case, the melting line continuously transforms into the order–disorder line, and we arrive at the second scenario. Moreover, since $`t_{\mathrm{up}}`$ is not a specific temperature for $`E_{\mathrm{pin}}(a)`$ or $`E_{\mathrm{el}}(a)`$, it is quite probable that the result $`\stackrel{~}{l}_{}=1`$ remains true also at $`t<t_{\mathrm{up}}`$ (at least in the bundle pinning region). Then, at $`tt_{\mathrm{up}}`$ the order–disorder line (which is now a part of the unified phase transition line) is determined by the condition,

$$E_{\mathrm{el}}(a)=E_{\mathrm{pin}}(a).$$
(55)
As it follows from the results of Sec. III C, in the bundle pinning region condition (55) is equivalent to the Lindemann criterion:

$`A_1u^2(a,0)\left[1+\mathrm{ln}\left({\displaystyle \frac{R_c^2}{a^2}}\right)+{\displaystyle \frac{ϵR_c}{\lambda _{ab}}}\right]^{1/2}=c_L^2a^2.`$ (56)
If the order–disorder line enters the single vortex pinning region, where $`R_c=a`$, condition (55) yields:

$`A_1u^2(a,0)=c_L^2a^2.`$ (57)
(Note that the ratio $`ϵa/\lambda _{ab}`$ is negligible.) In Appendix B equations (56) and (57) are presented in explicit form.

To calculate the order–disorder line $`H_{dis}(t)`$ in the case of the first scenario, we used the criterion:4

$$u^2(a,0)=c_{LD}^2a^2,$$
(58)
where $`c_{LD}`$ is some new Lindemann constant describing the order–disorder transition. It should be emphasized that although Eq. (58) is similar to Eq. (57), the Lindemann constant $`c_{LD}`$ in Eq. (58) is independent of the constant $`c_L`$ defining the melting, while in the case of the second scenario Eq. (57) leads to a relationship between the Lindemann constants defining the melting and the order–disorder lines (the ratio of these constants is $`\sqrt{A}_1`$). Criterion (58) both in the single vortex and in the bundle pinning46  regions were analyzed in Ref. 4, , allowing for a smoothing of the pinning potential by thermal fluctuations. The appropriate equations in the explicit form are presented in Appendix B. For definiteness, in the analysis below we shall choose the constant $`c_{LD}`$ for the order–disorder line as equal to the constant $`c_L`$ for the melting line.

### IV.2 Numerical results

Although in this paper we have not analyzed the order–disorder line in detail, criteria (56), (57), and (58) enable one to evaluate the location of this line in the $`T`$ \- $`H`$ plane for the first and the second scenarios. In Figs. 5 \- 8 we compare the $`T`$-$`H`$ phase diagrams of type-II superconductors with different types of pinning for these two scenarios. In the construction of the figures, we use the values of the Ginzburg number $`Gi=0.01`$, $`0.0001`$ which are typical for high-$`T_c`$ superconductors, and we take into account the factor containing $`[1h_{sv}(t)]`$ that has been omitted in Sec. III.42  The complete set of the appropriate equations is given in Appendix B.

In Figs. 5 \- 8 we choose the constants $`A_1`$ and $`B_1`$ so that the following two requirements are satisfied: First, the properties of $`H_m(t)`$ observed in experiments are reproduced, namely, with increasing $`D`$ the upper critical point clearly shifts, while the downward shift of the melting line is small, see Sec. III B. Second, the upper critical point $`t_{\mathrm{up}}`$ of the melting line coincides with the intersection point $`t_i`$ of this line with the order–disorder line (calculated within the first scenario) at $`D=0.7`$ for the case of $`\delta l`$ pinning. The latter requirement is due to the following considerations: According to Ref. 12, , 15, , the coincidence of $`t_i`$ and $`t_{\mathrm{up}}`$ is observed in overdoped YBa<sub>2</sub>Cu<sub>3</sub>O<sub>y</sub> crystals ($`y>6.92`$) for which the upper critical point lies at sufficiently large magnetic fields, and so the case $`h_m1`$ appears to occur there. As it was mentioned in Sec. III B, in this case $`t_i`$ and $`t_{\mathrm{up}}`$ depend on $`Gi`$ and $`D`$ only via a single combination of these parameters, $`D^3/Gi^{1/2}`$. Therefore, one may expect that if the coincidence of the upper critical point with the intersection point occurs, it practically will not depend on the specific choice of $`D`$ or $`Gi`$. In other words, the coincidence will approximately occur for different $`D`$, $`Gi`$, and types of pinning as long as $`h_i1`$. Although $`h_i`$ is not too close to unity in Figs. 5 \- 7 (top), the data of these figures support this statement. As one might expect, the region of the coincidence is especially wide in $`D`$ for small $`Gi`$, Fig. 6. Thus, we have introduced the second requirement here in order to fit the phase diagrams calculated in the framework of the first scenario to the experimental situation observed in the overdoped YBa<sub>2</sub>Cu<sub>3</sub>O<sub>y</sub> crystals.

The presented figures show that when the intersection point is close to the upper critical point, both scenarios lead to qualitatively similar phase diagrams. However, when the strength of pinning $`D`$ increases, and the intersection point shifts toward $`T_c`$, the coincidence of the points fails especially for $`\delta T_c`$ pinning. In this case the first scenario leads to a noticeable extension of the melting line beyond the intersection point. Thus, if the first scenario really occurs, this result possibly explains the experimental findings 11 ; 12 ; 13 ; 14 ; 15  for optimally doped YBa<sub>2</sub>Cu<sub>3</sub>O<sub>y</sub> crystals ($`y6.92`$) in which the upper critical point lies at larger magnetic field than the intersection point. Note also that in the case of $`\delta T_c`$ pinning, the melting line enters the single vortex pinning region near $`T_c`$ at $`\nu (2\pi )^{3/2}D^3/Gi^{1/2}`$ of the order of several units.4 ; 19

In Figs. 5 \- 8 we also show the boundary of the single vortex pinning region, $`H_{sv}(t)`$. It should be noted that apart from the well-known lower single vortex pinning region, 9  an upper region exists where this pinning occurs. The upper region of single-vortex pinning was discussed by Larkin and Ovchinnikov 47  in the context of the origin of the peak effect in low-$`T_c`$ superconductors. Without account of thermal fluctuations, this region adjoins to the $`H_{c2}(t)`$ line and is caused by the softening of the vortex lattice near $`H_{c2}(t)`$. However, this softening also leads to an increase of $`u_T`$, which reduces the strength of pinning. As a result of these two opposite tendencies, the upper region does not extend to $`T_c`$ and has the shape of a “tongue”.4  Interestingly, when the strength of pinning increases, the upper and lower single-vortex pinning regions can merge at low temperatures. This merging occurs for $`D(4\pi )^{1/2}`$; this value is independent of $`Gi`$ and the type of pinning since the merging starts at $`T=0`$.

The authors of Ref. 5,  argued that only the second scenario can explain the decrease of $`H_{dis}`$ with the temperature that was observed in their experiments. However, the presented figures show that the order–disorder line found from criterion (58) (the first scenario) can decrease or increase with temperature and even can be nonmonotonic. This depends on whether the line is in the single vortex pinning region or in the bundle pinning region and also on the Ginzburg number. It should be noted that in contrast to Eq. (3) we take into account only the displacement caused by the quenched disorder and do not include $`u_T`$ explicitly in the equation for the order–disorder line, compare Eq. (58) with Eq. (3). However, the thermal depinning (which depends on $`u_T`$), the softening of the elastic moduli, and the possibility that the order–disorder line lies not only in the single vortex pinning region but also outside it, already produce the depicted variety of shapes of $`H_{dis}(t)`$. Hence, even in the framework of the first scenario the presented results can explain the fact 5  that the order–disorder lines observed in experiments have various shapes.

Finally, we briefly discuss the case of low-$`T_c`$ superconductors, which have a very small Ginzburg number. In this case one has $`D^3/Gi^{1/2}1`$ even for weak pinning strength $`D`$. Since the temperatures $`T_i`$ and $`T_{\mathrm{up}}`$ are mainly determined by this ratio of the parameters and increase when the ratio increases, we find that in low-$`T_c`$ superconductors these temperatures practically coincide with $`T_c`$, and only the order–disorder line can be observed.48  Interestingly, in this case the order–disorder lines for the first and for the second scenarios have the same functional dependences on $`t`$ and on the parameters $`D`$, $`Gi`$ \[this follows from equations (80) and (83)\], and moreover, they practically coincide with each other when $`A_1^22\pi c_L^2`$.

###### Acknowledgements.

This work was supported by the German Israeli Research Grant Agreement (GIF) No G-705-50.14/01.

## Appendix A Effect of disorder on the melting line: analysis beyond perturbation theory.

Assuming $`1t1`$, we shall consider the partition function of the vortex system as a functional integral with the Ginzburg-Landau Hamiltonian.30 ; 34  In dimensionless units this Hamiltonian has the form:

$`{\displaystyle \frac{H_{GL}}{T}}={\displaystyle }d𝐫[|_z\psi |^2+|(i+𝐀)\psi |^2+\tau |\psi |^2+`$
$`2^{3/2}\pi Gi^{1/2}t|\psi |^4],`$ (59)
where $`\psi `$ is the order parameter; $`\tau =t1`$; the coordinates $`x`$ and $`y`$ are measured in units of $`\xi _0`$, and $`z`$ in units of $`ϵ\xi _0`$; $`\xi _0`$ is the zero-temperature coherence length in the Ginzburg-Landau theory; the magnetic field is $`b=2\pi \xi _0^2H/\mathrm{\Phi }_0`$ and $`𝐀`$ is its vector potential; $`=(_x,_y)`$. Consider the melting line $`h_m(t)`$ in the temperature region described by Eq. (12) in which $`1h_m1`$, and hence the lowest Landau level approximation is valid (see Sec. II C). In this approximation the second term of Eq. (A) reduces to

$$b|\psi |^2.$$
(60)
Pinning is introduced into the dimensional Ginzburg -Landau Hamiltonian either via spatial disorder in the transition temperature $`T_c+\delta T_c(𝐫)`$ ($`\delta T_c`$ pinning) or by spatial variation of the effective mass $`m+\delta m(𝐫)`$ describing disorder in the mean free path $`l`$ of quasiparticles ($`\delta l`$ pinning).9  Here $`m`$ is the effective mass in the $`x`$-$`y`$ plane (since the magnetic field is along the $`z`$ axis). Thus, with quenched disorder in the vortex lattice, one should add to the dimensionless Hamiltonian (A) the term

$$\phi (𝐫)|\psi |^2,$$
(61)
where

$$\phi (𝐫)=\frac{\delta T_c(𝐫)}{T_c}$$
(62)
in the case of $`\delta T_c`$ pinning, and

$$\phi (𝐫)=\frac{\delta m(𝐫)}{m}b\frac{\delta m(𝐫)}{m}|\tau |$$
(63)
in the case of $`\delta l`$ pinning. In Eq. (63) we have used the relation

$$b=|\tau |h$$
(64)
that follows from the definitions of $`b`$ and $`h`$ and put $`h=1`$ since $`1h1`$ in the lowest Landau level approximation. For pinning by point defects, it is assumed in the collective pinning theory 9  that disorder in $`\delta T_c(𝐫)`$ and in $`\delta m(𝐫)`$ is short scale and described by a Gaussian distribution with zero mean value, $`\delta T_c(𝐫)=\delta m(𝐫)=0`$, and with the correlation function: A

$$\frac{\delta T_c(𝐫)\delta T_c(𝐫^{})}{T_c^2}=2\pi D^3\delta (𝐫𝐫^{})$$
(65)
for $`\delta T_c`$ pinning, and

$$\frac{\delta m(𝐫)\delta m(𝐫^{})}{m^2}=\frac{30\pi }{7}D^3\delta (𝐫𝐫^{})$$
(66)
for $`\delta l`$ pinning. Here $`\mathrm{}`$ means disorder averaging. Note that in agreement with Ref. 4,  and with Sec. III A of this paper, equations (A) - (66) again show that the phase diagrams of type-II superconductors with point defects depend only on $`D`$ and $`Gi`$.

To proceed further, let us rescale the coordinates and the order parameter $`\psi `$ similarly to Ref. 35, :

$`\stackrel{~}{x}=\sqrt{b}x,`$ (67)
$`\stackrel{~}{y}=\sqrt{b}y,`$ (68)
$`\stackrel{~}{z}=zGi^{1/6}(tb)^{1/3},`$ (69)
$`\stackrel{~}{\psi }^2=\psi ^2Gi^{1/6}t^{1/3}b^{2/3},`$ (70)
where $`\stackrel{~}{x}`$, $`\stackrel{~}{y}`$, $`\stackrel{~}{z}`$, $`\stackrel{~}{\psi }`$ are the new coordinates and order parameter. Beside this, to agree with the notation used in the main text of this paper, we put $`\tau =t^21`$ below. Then, the Hamiltonian (A) transforms into:

$`{\displaystyle \frac{H_{GL}}{T}}={\displaystyle 𝑑\stackrel{~}{𝐫}\left[|_{\stackrel{~}{z}}\stackrel{~}{\psi }|^2Q|\stackrel{~}{\psi }|^2+2^{3/2}\pi |\stackrel{~}{\psi }|^4\right]},`$ (71)
where $`Q`$ is defined by Eq. (21). In the case of the ideal lattice, Eq. (23) for $`h_m`$ follows from this Hamiltonian, see Sec. II C. In a lattice with quenched disorder, the additional term from (61) in the Hamiltonian has the form:

$$\stackrel{~}{\phi }(\stackrel{~}{𝐫})|\stackrel{~}{\psi }|^2,$$
(72)
with $`\stackrel{~}{\phi }(\stackrel{~}{𝐫})=0`$ and

$$\stackrel{~}{\phi }(\stackrel{~}{𝐫})\stackrel{~}{\phi }(\stackrel{~}{𝐫}^{})=\frac{2\pi D^3}{Gi^{1/2}t}\delta (\stackrel{~}{𝐫}\stackrel{~}{𝐫}^{})$$
(73)
in the case of $`\delta T_c`$ pinning, and with

$$\stackrel{~}{\phi }(\stackrel{~}{𝐫})\stackrel{~}{\phi }(\stackrel{~}{𝐫}^{})=\frac{30\pi D^3(1t^2)^2}{7Gi^{1/2}t}\delta (\stackrel{~}{𝐫}\stackrel{~}{𝐫}^{})$$
(74)
in the case of $`\delta l`$ pinning.

Up to numerical factors, the right hand sides of Eqs. (73) and (74) coincide with the function $`[G(t)]^2`$ defined by Eq. (44). Then, the free energies of the vortex liquid, $`F_{\mathrm{liq}}`$, and the vortex lattice, $`F_{\mathrm{lat}}`$, are determined by $`Q`$ and $`G(t)`$, and along the melting line $`h_m(t)`$, one has:

$$F_{\mathrm{liq}}(Q,G(t))=F_{\mathrm{lat}}(Q,G(t)).$$
(75)
Thus, on the melting line, the quantity $`Q`$ is some function of $`G(t)`$. Note that the functional form of Eq. (III.2) agrees with this conclusion, see Sec. III B.

## Appendix B Equations for calculation of the phase diagrams

In section III where we considered the melting line in the bundle pinning region, the normalization factor containing $`[1h_{sv}(t)]`$ was omitted.42  The origin of this factor is due to the additional power of $`1h`$ in the shear modulus $`c_{66}`$ as compared with the tilt modulus $`c_{44}`$, see Eqs. (II.2). In the single vortex pinning region where the shear modulus does not play any role, the additional power should not manifest itself in physical properties. Therefore, if a physical quantity in the bundle pinning region contains a factor $`(1h)^n`$ caused by this additional power, it is necessary to introduce a normalization factor $`[1h_{sv}(t)]^n`$ in this quantity to provide its continuity at the boundary of the single vortex pinning region $`h_{sv}(t)`$. (Here $`n`$ is some power.) Although the effect of the factor $`[1h_{sv}(t)]^n`$ on the melting line $`H_m(t)`$ and on the order–disorder line $`H_{dis}(t)`$ in the bundle pinning regime is small, it essentially influences the boundary of the single vortex pinning region, $`H_{sv}(t)`$, and the line $`H_{dis}(t)`$ inside this region. In this Appendix B, taking into account this factor,B  we compile the complete set of equations used for the construction of Figs. 5 \- 8.

The boundary of the single vortex pinning region $`h_{sv}(t)=H_{sv}(t)/H_{c2}(t)`$ is described by equation (19) of Ref. 4, :

$`h_{sv}^{1/2}(t)+{\displaystyle \frac{F_T(t)}{[1h_{sv}(t)]^{3/2}}}=(2\pi )^{1/2}Dg_0(t),`$ (76)
where we have used the notation (41). The boundary of the upper single vortex pinning region, discussed in Sec. IV B, can be obtained from the equation: 4

$`{\displaystyle \frac{1h}{1h_{sv}(t)}}\left[h^{1/2}+{\displaystyle \frac{F_T(t)}{(1h)^{3/2}}}\right]^2=2\pi \left(Dg_0(t)\right)^2,`$ (77)
which in addition reproduces the root of Eq. (76), i.e., it yields the entire boundary of the single vortex pinning regions. Inside this upper region the vortex lattice is in a state where $`R_c=a`$ and $`u(a,0)=r_p`$, i.e., a borderline state between the single vortex pinning and bundle pinning regimes occurs there.

Equation (39) for the melting line in the bundle pinning region is rewritten as follows:

$`{\displaystyle \frac{F_T(t)h_m^{1/2}}{(1h_m)^{3/2}}}\left[1P(t,h_m)\right]`$ $`+`$
$`{\displaystyle \frac{A[Dg_0(t)]^{3/2}h_m^{1/4}[1h_{sv}(t)]^{3/4}}{\left[(1h_m)^{3/2}+F_T(t)h_m^{1/2}\right]^{1/2}}}`$ $`=`$ $`2\pi c_L^2,`$ (78)
with

$$P(t,h_m)=\frac{B[Dg_0(t)]^3(1t^2)^{1/2}[1h_{sv}(t)]}{tGi^{1/2}\left[h_m^{1/2}+{\displaystyle \frac{F_T(t)}{(1h_m)^{3/2}}}\right]^2}.$$
(79)
Equation (B) is valid in the interval $`t_{\mathrm{up}}t1`$ where the temperature $`t_{\mathrm{up}}`$ defines the position of the upper critical point of the melting line and is given by the condition (42). Note that in this paper we consider situations when the upper critical point lies in the bundle pinning region.

In the framework of the second scenario the equation for the order–disorder transition line (which is the continuation of the melting line to $`t<t_{\mathrm{up}}`$) in the bundle pinning region follows from criterion (56) as:

$`{\displaystyle \frac{A[Dg_0(t)]^{3/2}h_m^{1/4}[1h_{sv}(t)]^{3/4}}{\left[(1h_m)^{3/2}+F_T(t)h_m^{1/2}\right]^{1/2}}}=2\pi c_L^2,`$ (80)
i.e., the first term in Eq. (B) disappears. If this order–disorder line enters the single vortex pinning regions, it is described by equations that result from criterion (57). This criterion reads in explicit form

$`A_1h_{sv}^{3/5}(t)h_m^{2/5}\left[1+{\displaystyle \frac{F_T(t)}{h_{sv}^{1/2}(1h_m)^{3/2}}}\right]=2\pi c_L^2`$ (81)
for the order–disorder line in the lower single vortex pinning region, while in the upper single vortex pinning region one finds:

$`A_1F_T(t){\displaystyle \frac{h_m^{1/2}}{(1h_m)^{3/2}}}+A_1h_m=2\pi c_L^2.`$ (82)
Here we have used that $`u(a,0)=r_p`$ in the upper region.

In the framework of the first scenario, equations for the order–disorder line follow from condition (58).4  With the use of our approximation $`f_1=2`$, equation (24) of Ref. 4,  for the order–disorder line $`H_{dis}(t)`$ in the bundle pinning region can be written in the form:

$$h_{dis}^{1/2}(1h_{dis})^{3/2}=F_T(t)K_\pm (t),$$
(83)
where $`h_{dis}=H_{dis}(t)/H_{c2}(t)`$,

$$K_\pm (t)=G_1^21\pm \left[(G_1^21)^21\right]^{1/2},$$
(84)
and

$$G_1=\frac{\pi ^{1/4}G(t)[1h_{sv}(t)]^{3/4}}{2^{1/4}c_L}$$
(85)
with $`G(t)`$ from Eq. (44). The order–disorder line $`H_{dis}(t)`$ in the single vortex pinning region is given by Eq. (22) of Ref. 4, :

$`h_{dis}\left[1+{\displaystyle \frac{F_T(t)}{(1h_{dis})^{3/2}[h_{sv}(t)]^{1/2}}}\right]^{5/2}=`$
$`2\pi c_L^2\left({\displaystyle \frac{2\pi c_L^2}{h_{sv}(t)}}\right)^{3/2}.`$ (86)
Note that apart from the factor $`A_1`$ Eq. (B) is equivalent to Eq. (81).


text:
# Is Galaxy Distribution Non-extensive and Non-Gaussian?

## 1 Introduction

The Long-range unshielded force gravity forms unique structures in the Universe called self-gravitating systems (SGS). Almost all the relevant structures such as stars, galaxies, clusters and super-clusters of galaxies, belong to this category.

These SGS generally (a) show non-extensive properties, and (b) have no absolute equilibrium state; the systems never stop their evolution toward singularity.

The former property (a) is apparent in principle. Actually, the mass density $`\rho `$ of iso-thermal SGS of size $`r`$ is given by

$`\rho ={\displaystyle \frac{M}{4\pi r^3/3}}={\displaystyle \frac{<v^2>}{4\pi Gr^2/3}}r^2`$ (1)
where we used the virial theorem: $`K`$=$`V/2`$, where $`K`$ and $`V`$ are respectively the kinetic energy and the potential energy of the system. Since the general extensive structures have definite basic elements and therefore have constant density, the above scale-dependent density clearly shows the non-extensivity of SGS. This property reflects the long distance singularity of gravity. In other words, SGS cannot be divided into uniform independent ingredients (non-additive) in general.

The latter property (b) is also apparent in principle. Actually, the partition function $`Z`$ of SGS in the temperature $`T(=1/\beta )`$, again using the virial theorem, is given by

$`Z`$ $`=`$ $`{\displaystyle e^{\frac{\left(K+V\right)}{kT}}𝑑x_1𝑑p_1\mathrm{}𝑑x_N𝑑p_N}`$
$``$ $`{\displaystyle e^{\beta V/2}𝑑x_1\mathrm{}𝑑x_N}`$
$`=`$ $`{\displaystyle \mathrm{exp}[\beta G/2\underset{i<j}{}\left|x_ix_j\right|^1]𝑑x_1\mathrm{}𝑑x_N}`$
which has an essential singularity at the short distance limit. Therefore the absolute equilibrium state does not exist. This property reflects the short distance singularity of gravity.

Since the ordinary Boltzmann statistical mechanics is based on the extensivity and the existence of the absolute equilibrium state, the above properties of gravity make it inapplicable in principle.

Despite the property (b), SGS frequently show quasi-equilibrium states which are often characterized by scaling properties; especially (c) the long-tail in various distribution functions are marvelous. This fact may suggest the existence of any underling fundamental statistical mechanics for SGS.

There are many examples of the above properties (a) and (c) in SGS: The actual mean mass density of a structure is proportional to $`r^{1.7}`$ (de Vaucouleurs 1970)(Eq.(1) is too simplified). Moreover in the extreme gravitating case of Black Hole, the ‘entropy’ is not proportional to the mass $`M`$ but to the area $`M^2`$ . The property (c) is based on the absence of the characteristic scale in the Hamiltonian of SGS. The well known Holtzmark distribution of force $`F`$ acting on a single star in a uniform star cluster is exactly the stable Levy distribution with index 3/2 and behaves as $`F^{5/2}`$ in large force limit (Binney et al. 1987). The correlation function of galaxies and the clusters show power law behavior $`\xi r^{1.86}`$ (Börner 2002).

As explained in the above, the ordinary Boltzmann statistical mechanics cannot be applied to SGS in the naive form, and therefore any generalization or new statistical mechanics would be necessary in order to explain the quasi-equilibrium states of SGS.

The necessary next step will be as follows, on which we will concentrate in this paper:

1. To choose astronomical data which reflect the above properties (a), (b), and (c) in the fundamental level. The data should be investigated not simply by the apparent form of the distribution functions but by internal structure of the theory that distinguish various statistical mechanics.
2. To introduce a fair criterion which can specify the correct statistical mechanics among various proposed theories which have different number of free parameters.

For the purpose (A), we focus on the data of CfA II South (Huchra J., et al., 1999) galaxy distribution survey 111 We do not claim that this is the best available data, though it is uniform and large to some extent. . Especially, we use count-in-cell method, in which the probability $`f_V\left(N\right)`$ of finding a fixed number of galaxies in a fixed volume plays a central role. For the purpose (B), we focus on the Akaike Information Criteria (AIC) (Akaike 1973) 222 We do not claim that this is the best available method; there are many extensions and generalizations of AIC. These are unnecessarily complicated for the present purpose and therefore we choose the simplest version of AIC. . This method will be most appropriate to compare different theories with different number of free fitting parameters.

We have chosen the following four theories/models of statistical mechanics, classified according to the properties (a) non-extensivity and (c) long-tails. (Table 1.)

* (1) The Boltzmann statistical mechanics which is of course extensive from its construction. The associated Boltzmann distribution function does not possess long tail. Since this theory alone never explain the CfAII data (see below), we introduce an extra parameter $`b`$ as in the reference by Saslaw W. C. & Hamilton A. J. S. (1984). This parameter measures the possible deviation from the complete virial equilibrium.
* (2)The (mono-)fractal-space model as a simple example of non-extensive theory with short tail distribution.
* (3)The Rényi statistical mechanics as an appropriate example of extensive statistical mechanics with long-tail distribution.
* (4)The Tsallis statistical mechanics (Nakamichi A., Joichi I., Iguchi O., Morikawa M., 2001), as a typical non-extensive theory. This has the same long-tail distribution as Rényi. It should be noted that we use proper Tsallis statistical mechanics with Escote averaging, which should be distinguished from the older version of the normal averaging. Both Tsallis and Rényi statistical mechanics posses an extra parameter $`q`$ which measures the deviation from the ordinary Boltzmann statistical mechanics; they reduce to the latter in the limit $`q1`$ .

Individual models/theories are explained in detail below.

## 2  Various models/theories of statistical mechanics

### 2.1 Boltzmann statistical mechanics

We consider the galaxy distribution which is supposed to obey the Boltzmann statistical mechanics with grand canonical ensemble. In the ordinary Boltzmann statistical mechanics, the tail of the distribution function exponentially reduces. We phenomenologically generalize the theory and introduce the virial parameter $`b`$, which measures the deviation from the complete dynamical-equilibrium:

$$pV=NT(1b),$$
(3)
where $`p`$ denotes the pressure, $`V`$ the volume of the system, $`N`$ the number of galaxies contained in this volume, and $`T`$ the temperature defined by the velocity dispersion.

Usually, collision-less SGS attain the dynamical equilibrium before the thermal equilibrium: $`\tau _{\mathrm{thermal}}\frac{N}{\mathrm{ln}N}\tau _{\mathrm{dynamical}}\tau _{\mathrm{dynamical}}`$ (see J. Binney et al. 1987 for detail). Therefore the parameter $`b`$ phenomenologically represents some class of quasi-equilibrium .

The distribution function $`f(0)`$ is defined to be the probability of finding no galaxy in any part of the volume $`V`$. The explicit form is given by

$`f\left(0\right)=e^{\overline{N}\left(1b\right)}=e^{n\left(1b\right)4\pi r^3/3}`$ (4)
where $`n`$ is the galaxy number density and $`\overline{N}=nV`$ (Saslaw 1984). This is the generating functional of the general function $`f(N)`$, which is defined to be the probability to find the $`N`$ galaxies in the fixed volume V. The general expression of $`f(N)`$ is given in the Appendix.

From the observational data of CfA II South galaxy redshift survey, we calculate the probability $`f(N)`$ using the count-in-cell method. The best fit parameter $`b`$ is calculated to minimize the Akaike Information Criteria (AIC) for $`f(0)`$. The comparison of the theory and observation is shown in Fig. 1. The diamonds with error bars are CfA II South Observations. The void probability $`f(0)`$ calculated from the Boltzmann statistical mechanics is plotted by the broken line with the best fit parameter $`b=0.523`$. Quite a large deviation between the theory and the observation implies the inapplicability of the Boltzmann statistical mechanics even in the modified version.

We fix the value of the parameter and with this value, general probabilities $`f(N)`$ are calculated and are compared with the CfA II South data in Fig. 2.

### 2.2 Boltzmann statistical mechanics with fractal matter distribution

Various observational data suggest the idea that the matter distribution in the Universe is (multi-)fractal at least in some limited scale region (Kurokawa et al. 1999; 2000). In this section, we investigate a simple fractal model with the ordinary Boltzmann statistical mechanics. In this model, the system is non-extensive in the sense that the simple addition of the two identical fractal distribution in the three dimensional space generally destroys the scaling property of the original fractal system. On the other hand, the tail of the distribution function exponentially reduces since the model is based on the ordinary Boltzmann statistical mechanics.

We consider a mono-fractal matter distribution with the fractal dimension $`\alpha `$, and suppose the system obeys Boltzmann statistical mechanics with grand canonical ensemble, as in the previous section.

The void probability is simply given by $`f\left(0\right)=e^{\overline{N}}=e^{nV}`$ , but the number $`\overline{N}`$ non-trivially depends on the scale $`r`$ in our fractal model. Thus $`f(0)`$ is given by

$`f\left(0\right)=e^{\overline{N}}=\mathrm{exp}\left[n{\displaystyle \frac{\pi ^{\alpha /2}r^\alpha }{\mathrm{\Gamma }\left(1+\alpha /2\right)}}\right]`$ (5)
This expression is shown in Fig.3 with the observational data of CfAII redshift survey. By using the same method as in the previous section, we found the best fit value for the parameter $`\alpha =2.8821`$.

With this best fit parameter value, general probabilities $`f(N)`$ are calculated using the method in the Appendix, and they are compared with the CfA II South data in Fig. 4.

### 2.3 Rényi statistical mechanics

In this section, we study Rényi statistical mechanics, which is extensive but the associated distribution function has long-tail. Rényi statistical mechanics is a generalization of the ordinary Boltzmann theory by introducing a new form of entropy called Rényi entropy. This entropy is very similar to the information measure often used in multi-fractal models. However the use of this entropy as the starting point of the solid statistical mechanics is not clear. Therefore we concentrate on the application of this statistical mechanics in this paper without discussion of its physical foundation.

The Rényi entropy is defined by

$$S[p]=\frac{\mathrm{ln}(_ip_{i}^{}{}_{}{}^{q})}{1q},$$
(6)
where the parameter $`q`$ measures the deviation from the ordinary Boltzmann entropy. Actually in the limit$`q1`$ , this reduces to the Boltzmann entropy.

When we compose two independent systems $`A`$ and $`B`$, the distribution function is given by $`p_{i,j}[A,B]=p_i[A]p_j[B]`$, and the composed entropy becomes

$$S_{A+B}=S_A+S_B,$$
(7)
which clearly shows extensivity.

Let us suppose that the galaxy distribution obeys the Rényi statistical mechanics with the grand canonical ensemble. Then the distribution function is given by

$$p_{N,E}=[1\frac{1q}{T}\{E\overline{E}\mu (N\overline{N})\}]^{\frac{1}{1q}}$$
(8)
which maximizes the above entropy.

As is clearly seen from the above expression, the distribution function has a long-tail with power-law shape; a significant characteristic which distinguish this formalism and the Boltzmann formalism.

Since the total entropy for $`N`$ galaxies, $`S_N`$, is simply the sum of $`N`$ entropies for a single galaxy: $`S_N=sN`$, the Euler relation becomes the ordinary one:

$$S=\frac{E+pV\mu N}{T}.$$
(9)
Using the above expression, the void probability becomes

$$f(0)=\{1+(1q)Ns\}^1.$$
(10)
Because the change of the parameter $`s`$ can be renormalized into the change of the parameter $`q`$, the probability is independent of $`s`$. This fact naturally reflects the extensiveness of the present entropy.

As in the previous section, using the count-in-cell method, we compare the theoretical model and the observational data of CfAII South.

The best fit value of the parameter $`q`$ is $`q=0.042`$. By fixing this value, we further obtain general probabilities $`f(N)`$ and they are compared with the CfA II South data.

### 2.4 Tsallis statistical mechanics

Tsallis C. (1988) proposed the non-extensive entropy

$$S[P]=\frac{(_iP_{i}^{}{}_{}{}^{1/q})^q1}{1q},$$
(11)
where the distribution function $`P_i`$ is the Escort distribution and is related with the bare distribution function $`p_i`$ by:

$$P_i=p_{i}^{}{}_{}{}^{q}/\underset{j=1}{\overset{W}{}}p_{j}^{}{}_{}{}^{q}p_{i}^{}{}_{}{}^{q}/C.$$
(12)
The Escort distribution is identified to be the physical distribution and is used to obtain observable averaging. The distribution function is derived by maximizing the above entropy functional with appropriate constraints.

When we compose two independent systems $`A`$ and $`B`$, and for the distribution $`p_{i,j}[A,B]=p_i[A]p_j[B]`$, composed Tsallis entropy satisfies the non-extensive relation:

$$S_{A+B}=S_A+S_B+(1q)S_AS_B,$$
(13)
as is easily calculated by using the above entropy form.

Let us suppose the galaxy distribution obeys the Tsallis statistical mechanics with grand canonical ensemble. As in the case of fractal model and Rényi statistical mechanics, we do not need the virial parameter $`b`$ in Tsallis statistics.

The distribution function

$$p_{N,E}=\frac{1}{\mathrm{\Xi }_q}[1\frac{1q}{\stackrel{~}{T}}\{E\overline{E}\mu (N\overline{N})\}]^{\frac{1}{1q}}$$
(14)
maximizes Tsallis entropy. Here $`\stackrel{~}{T}CT`$ should be identified as the ‘temperature’ of the system. The factor $`C`$ appears everywhere in this theory.

Note that the parameter $`q`$ measures the extent of deviation from the ordinary Boltzmann statistical mechanics; the distribution function reduces to the Boltzmann distribution for$`q1`$.

Tsallis distribution is identical to the Rényi distribution except the ‘temperature’ and the normalization factor. Therefore the Tsallis distribution also has the long-tail.

Since the system is non-extensive, the total entropy of the $`N`$ galaxies should be calculated by using the composition law of the entropy given above. The expression is given by (Nakamichi A., et al., 2001)

$$S_N=\frac{\{1+(1q)s\}^N1}{1q},$$
(15)
where the parameter $`s`$ is identified to be the entropy for a single galaxy.

Then we obtain the generalized Euler relation à la Tsallis:

$$\frac{\{1+(1q)S\}\mathrm{log}[1+(1q)S]}{1q}=\frac{E+pV\mu N}{T}.$$
(16)
As in the previous two sections, using the count-in-cell method, probability of void becomes

$$f(0)=\{1+(1q)s\}^{\frac{N}{1q}}[1+N\mathrm{log}\{1+(1q)s\}]^{\frac{q}{1q}}.$$
(17)
This prediction is compared with the observational data of CfAII South. We found the best fit parameters $`q=5.7`$ and $`s=0.16`$. Almost complete fitting is remarkable.

Further, we can calculate the general probability of finding $`N`$ galaxies in a given volume $`V`$, and compare with the CfA II South data in the figures 7,8.

We observe that for the higher order probability functions, the prediction by the Tsallis statistical mechanics is worse than the previous Rényi statistical mechanics, but is better than the Boltzmann statistical mechanics.

We have used the Escort averaging and not used the normal (old type) averaging in Tsallis statistical mechanics because the latter averaging is not a consistent theory. Tsallis originally proposed the non-extensive entropy using the normal averaging $`p_i`$ (see for example, Tsallis C., Mendes R. S., Plastino A. R., 1998)

$$S[p]=\frac{(_ip_{i}^{}{}_{}{}^{q})1}{1q},$$
(18)
where distribution function $`p_i`$ in normal averaging is determined so as to maximize the entropy functional.

If we suppose the galaxy distribution obeys the above normal-averaging Tsallis statistical mechanics with grand canonical ensemble, the distribution function

$$p_{N,E}=\frac{1}{\mathrm{\Xi }_q}[1\frac{1q}{\stackrel{~}{T}}\{E\overline{E}\mu (N\overline{N})\}]^{\frac{1}{1q}}$$
(19)
maximizes the Tsallis entropy. Here $`\stackrel{~}{T}CT`$ is the ‘temperature’ of the system. Also in this case, the distribution function approaches the ordinary Boltzmann distribution function for the parameter $`q`$ approaches to 1.

The void probability is obtained in the same manner as the case of Escote averaging, but the expression has a slightly different form:

$$f(0)=\{1+(1q)s\}^{\frac{N}{1q}}[1+N\mathrm{log}\{1+(1q)s\}]^{\frac{1}{1q}}.$$
(20)
This time we have found the best fit values $`q=1000`$ and $`s=10^{200}`$. Even in this best fit case, there exists still significant difference between the theory and observational data. The extremely large values above imply the absurdity of the theory.

## 3 Comparison of statistical mechanics using Akaike Information Criterion

In order to calculate $`f(N)`$, we have studied several kinds of statistical mechanics, in which the number of free parameters are different. In general, a model with much more parameters can fit the observational data better. However good physical theories must be simple and include the minimum number of free parameters. We need a fair measure to choose the correct statistical model to describe the galaxy distributions. One such measure will be the Akaike Information Criterion (AIC), (Akaike H. 1973; see for example, K. P. Burnham and D. R. Anderson 2002), which imposes a penalty on larger number of free parameters.

AIC measure is given by

$`AIC`$ $`=`$ $`2\times (\mathrm{maximum}\mathrm{log}\mathrm{liklihood})`$ (21)
$`+`$ $`2\times (\mathrm{number}\mathrm{of}\mathrm{parameters})`$
which is the unbiased estimator for expected mean log-liklihood. The first term on the right hand side tends to decrease as more parameters are introduced for better fitting, while the second term then becomes larger. This is the tradeoff between the number of parameters and the variance. Better model has smaller AIC measure.

Here the “parameters” denotes the free parameters to be estimated by the maximum likelihood analysis. For the present argument, they are the variance estimated by the maximum likelihood analysis, and the other free parameters in each theoretical model.

Here the “maximum log liklihood” is calculated as follows. When we compare observational data sets $`\{x_i,y_i\}_{1iN}`$ and the theoretical model $`y=P(x)`$, we assume that the errors (distance between the data and the theoretical model) distribute according to the Normal distribution: $`f(x_i,y_i)=\frac{1}{\sqrt{2\pi \sigma ^2}}e^{\frac{1}{2\sigma ^2}\left(y_iP(x_i)\right)^2}.`$ Then the log-likelihood function becomes

$`\mathrm{}\left(\theta \right)`$ $``$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}f(x_i,y_i)`$
$`=`$ $`{\displaystyle \frac{N}{2}}\mathrm{ln}\left(2\pi \sigma ^2\right){\displaystyle \frac{1}{2\sigma ^2}}{\displaystyle \left(y_iP(x_i)\right)^2},`$
where $`\theta `$ represents all the parameters in the model. The maximum likelihood estimate equations are

$`{\displaystyle \frac{\mathrm{}\left(\theta \right)}{\theta }}`$ $`=`$ $`\mathrm{\hspace{0.33em}0},`$ (23)
$`{\displaystyle \frac{\mathrm{}\left(\theta \right)}{\sigma ^2}}`$ $`=`$ $`{\displaystyle \frac{N}{2\sigma ^2}}+{\displaystyle \frac{1}{2\sigma ^4}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(y_iP(x_i)\right)^2=\mathrm{\hspace{0.33em}0}.`$
The first equation determines the parameters $`\widehat{\theta }`$. The second equation determines the variance $`\sigma ^2`$ of the error to be estimated value $`\widehat{\sigma }^2`$ as

$`\widehat{\sigma }^2={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(y_iP(x_i)\right)^2.`$ (24)
Then the maximum log-likelihood function becomes $`\mathrm{}\left(\widehat{\theta }\right)=\frac{N}{2}\frac{N}{2}\mathrm{ln}2\pi \widehat{\sigma }^2.`$ Therefore AIC measure reduces to

$`AIC`$ $`=`$ $`(\mathrm{\#}\mathrm{data}\mathrm{points})\times (1+\mathrm{log}_{10}2\pi (\mathrm{variance}))`$ (25)
$`+`$ $`2\times (1+\mathrm{number}\mathrm{of}\mathrm{free}\mathrm{parameters})`$
We have compared the AIC values for the probability of finding N galaxies in each theoretical model. In the figure 9, such AIC values are plotted.

For the void probability ($`N=0`$), we observe that the Tsallis statistical mechanics with the Escote averaging has a remarkably small AIC. On the other hand for higher order probabilities (($`1N`$)), the Rényi statistical mechanics has much smaller AIC than the Tsallis with the Escote averaging. Therefore only from these data, it is not possible to select the correct model among Tsallis and Rényi statistical mechanics. Therefore we cannot determin the (non)-extensivity of SGS. However, it is clear that the other models have large AIC values for all the probabilities. Therefore we may conclude that the long tail of the distribution function would be essential for describing the distribution of galaxies.

## 4 Discussions

We have studied two specific properties of SGS, non-extensivity of the entropy $`S`$, and the long-tails in the distribution functions, by explicitly introducing several models of statistical mechanics. By using the data of CfA II South redshift survey and the count-in-cell method, we have studied four statistical mechanics: (1) Boltzmann, (2) Fractal, (3) Rényi, and (4) Tsallis. These statistical mechanics have been evaluated by the Akaike information criteria (AIC) for the fair comparison.

In our study, it has been essential that the expression of the probabilities $`f(N)`$ depends not only the distribution function but also the explicit form of the entropy. The long-tail property of distribution function affects the former dependence and the (non-)extensivity of entropy affects the latter dependence. Therefore we could distinguish the Rényi and Tsallis theories despite that they have almost the same form of distribution functions $`P_{N,E}`$ (see for example, Arimitsu T. & Arimitsu N. preprint 2002 ).

First we have seen that in Boltzmann statistical mechanics, it is unable to fit the CfAII South data of galaxy distribution even if we introduce the free parameter $`b`$ which measures the deviation from the complete dynamical-equilibrium, or the free parameter $`\alpha `$ which measures the fractal property of the space.

Then we have investigated two statistical mechanics which have long tail distributions. Both the Rényi extensive statistical mechanics and the Tsallis non-extensive statistical mechanics are far better than the above two models based on Boltzmann statistical mechanics. Therefore the long-tail in the distribution function will be essential for the correct statistical mechanics to describe SGS.

On the other hand for the (non-) extensivity, we cannot have clear conclusions. This is because the Tsallis non-extensive theory can fit $`f(0)`$ better than Rényi extensive theory, and the latter theory can fit $`f(N)`$ ($`N1`$ ) better than the former theory. This incomplete conclusion may reflect that the size of the CfA survey data may still be small and much comprehensive survey would derive the complete answer to our question. We believe our method can in principle select the correct theory of statistical mechanics to describe SGS out of many candidate theories.

Acknowledgement One of the authors (MM) would like to thank Prof. Hiroshi Hasegawa for very useful discussions and suggestions.

## Appendix A Probabilities $`P\left(N\right)`$ and the generating functional

When we calculate the provability $`P\left(N\right)`$ of finding N galaxies in a given volume V, we use the expression

$$P\left(N\right)\frac{(n)^N}{N!}\frac{d^N}{dn^N}P(0)$$
(26)
where $`n`$ is the number density of galaxies. (White S. D. M. 1979)

The proof of this expression is originally given by S. White in the above reference with cumbersome arguments. In this Appendix, we extract the essence of the proof and try to clarify the logic as much as possible.

For a general statistical variable $`\varphi \left(x\right)`$, which is a field on three-dimensional space, the partition function is given by

$`Z\left[J\right]`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle _V}J\left(x\right)\varphi \left(x\right)d^3x\right]`$
$`=`$ $`\mathrm{exp}{\displaystyle _V}J\left(x\right)\varphi \left(x\right)d^3x_c`$
$``$ $`\mathrm{exp}W\left[J\right]`$
where the brackets represent the functional integral of the field $`\varphi \left(x\right)`$ (or the sum over possible functions $`\varphi \left(x\right))`$, and $`J\left(x\right)`$ is a source field. The whole part of the cumulants or the connected correlation functions, denoted as the brackets with a suffix c, are defined by the above equation.

We now specify the field $`\varphi \left(x\right)`$ as (discontinuous, non-averaged) number density of galaxies:

$$\varphi \left(x\right)\underset{i=1}{\overset{\mathrm{}}{}}\delta ^3\left(x_ix\right)$$
(28)
and accordingly the functional integration is reduced to the following form

$$\mathrm{}_V\mathrm{}\frac{d^3x_1}{V}\frac{d^3x_2}{V}\frac{d^3x_3}{V}\mathrm{}$$
(29)
since a distribution of galaxies at $`\{x_1,x_2,x_3,\mathrm{}x_n,\mathrm{}\}`$ determines a field $`\varphi \left(x\right)`$. It is obvious that the probability of finding one galaxy within a small volume $`d^3x_1`$ around the space position $`x_1`$ is $`P\left(x_1\right)d^3x_1=\varphi \left(x_1\right)d^3x_1`$, and the joint probability of finding one galaxy within a small volume $`d^3x_1`$ around the space position $`x_1`$ and the other within $`d^3x_2`$ around $`x_2`$ is $`P(x_1,x_2)d^3x_1d^3x_2=\varphi \left(x_1\right)\varphi \left(x_2\right)d^3x_1d^3x_2`$, and so on. Much more useful quantities are the probability of finding no galaxy $`P\left(0\right)`$ within a fixed volume $`V`$. Suppose the volume $`V`$ is divided into small pieces $`\{v_1,\mathrm{},v_M\}`$, each of which is of order $`v=V/M`$. Then this void probability is

$`P\left(0\right)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}\left(1\varphi \left(x_m\right)v_m\right)`$
$`=`$ $`\mathrm{exp}\left[{\displaystyle \underset{m=1}{\overset{M}{}}}\varphi \left(x_m\right)v_m\right]+O\left(Mv^2\right)`$
Taking the limit $`M\mathrm{}`$ with fixed $`V`$, and therefore $`O\left(Mv^2\right)=O\left(V^2/M\right)`$, the void probability simply becomes

$`P\left(0\right)`$ $`=`$ $`{\displaystyle \underset{xV}{}}\left(1\varphi \left(x\right)d^3x\right)`$
$`=`$ $`\mathrm{exp}\left[{\displaystyle _V}\varphi \left(x\right)d^3x\right]=Z\left[1\right]`$
Similarly, the probability of finding one galaxy within a small volume $`d^3x_1`$ around the space position $`x_1`$ and finding no other galaxies is given by

$`P(x_1;\mathrm{\hspace{0.33em}0})d^3x_1`$
$`=`$ $`\underset{M\mathrm{}}{lim}\varphi \left(x_1\right)v_1{\displaystyle \underset{m=2}{\overset{M}{}}}\left(1\varphi \left(x_m\right)v_m\right)`$
$`=`$ $`\underset{M\mathrm{}}{lim}\varphi \left(x_1\right)v_1\mathrm{exp}\left[{\displaystyle \underset{m=1}{\overset{M}{}}}\varphi \left(x_m\right)v_m\right]+O\left(Mv^2\right)`$
$`=`$ $`\varphi \left(x_1\right)\mathrm{exp}\left[{\displaystyle _V}\varphi \left(x\right)d^3x\right]d^3x_1`$
General probability for N galaxies is given by

$`P(x_1,\mathrm{},x_N;\mathrm{\hspace{0.33em}0})d^3x_1\mathrm{}d^3x_N`$ (33)
$`=`$ $`\begin{array}{c}\varphi \left(x_1\right)\mathrm{}\varphi \left(x_N\right)d^3x_1\mathrm{}d^3x_N\times \hfill \\ \underset{xV\{x_1,x_2,\mathrm{},x_N\}}{}\left(1\varphi \left(x\right)d^3x\right)\hfill \end{array}`$
$`=`$ $`\varphi \left(x_1\right)\mathrm{}\varphi \left(x_N\right)\mathrm{exp}\left[{\displaystyle _V}\varphi \left(x\right)d^3x\right]d^3x_1\mathrm{}d^3x_N`$
On the other hand, the partition function can be expanded as

$`Z\left[J\right]`$ $`=`$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left(1\right)^l}{l!}}{\displaystyle \begin{array}{c}\varphi \left(x_1\right)\varphi \left(x_2\right),\mathrm{}\varphi \left(x_l\right)\hfill \\ \times J\left(x_1\right)J\left(x_2\right),\mathrm{}J\left(x_l\right)\hfill \\ \times dV_1dV_2\mathrm{}dV_l\hfill \end{array}}`$ (40)
$`=`$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left(1\right)^l}{l!}}{\displaystyle \begin{array}{c}P(x_1x_2,\mathrm{},x_l)\hfill \\ \times J\left(x_1\right)J\left(x_2\right),\mathrm{}J\left(x_l\right)\hfill \\ \times dV_1dV_2\mathrm{}dV_l\hfill \end{array}}`$ (44)
In general,

$`P(x_1,\mathrm{},x_N)`$ $`=`$ $`\varphi \left(x_1\right)\mathrm{}\varphi \left(x_N\right)`$
$`=`$ $`\left(1\right)^N{\displaystyle \frac{^NZ\left[J\right]}{J\left(x_1\right)\mathrm{}J\left(x_N\right)}}|_{J0}`$
and

$`P(x_1,\mathrm{}x_N;\mathrm{\hspace{0.33em}0})`$ $`=`$ $`\varphi \left(x_1\right)\mathrm{}\varphi \left(x_N\right)e^{_V\varphi }`$
$`=`$ $`\left(1\right)^N{\displaystyle \frac{^NZ\left[J\right]}{J\left(x_1\right)\mathrm{}J\left(x_N\right)}}|_{J1}`$
The probability of finding exactly $`N`$ galaxies with the volume $`V`$, using $`P(x_{1,}x_2,\mathrm{}x_N;\mathrm{\hspace{0.33em}0})`$ but not $`P(x_{1,}x_2,\mathrm{}x_N)`$, is given by

$`P\left(N\right)`$ $`=`$ $`{\displaystyle \frac{1}{N!}}{\displaystyle _{V_1}}\mathrm{}{\displaystyle _{V_N}}P(x_{1,}x_2,\mathrm{}x_N;\mathrm{\hspace{0.33em}0})`$
$`=`$ $`{\displaystyle \frac{\left(1\right)^N}{N!}}{\displaystyle _{V_1}}\mathrm{}{\displaystyle _{V_N}}{\displaystyle \frac{^NZ\left[J\right]}{J\left(x_1\right)\mathrm{}J\left(x_N\right)}}|_{J1}`$
If we factor out the mean number density $`n`$ from the density field $`\varphi \left(x\right)`$ as $`\varphi \left(x\right)=n\left(1+\delta \left(x\right)\right)`$, where $`\delta \left(x\right)`$ denotes deviation from the average, then we have

$$_V\frac{Z\left[J\right]}{J\left(x\right)}|_{J1}=n\frac{Z\left[1\right]}{n}$$
(48)
and similarly

$$_{V_1}\mathrm{}_{V_N}\frac{^NZ\left[J\right]}{J\left(x_1\right)\mathrm{}J\left(x_N\right)}|_{J1}=n^N\frac{^NZ\left[1\right]}{n^N}$$
(49)
Using this form and Eq.(A) for Eq.(A), we finally have

$$P\left(N\right)=\frac{\left(n\right)^N}{N!}\frac{^N}{n^N}P\left(0\right)$$
(50)
as desired.