text:
# An entanglement measure based on the capacity of dense coding

## Abstract

An asymptotic entanglement measure for any bipartite states is derived in the light of the dense coding capacity optimized with respect to local quantum operations and classical communications. General properties and some examples with explicit forms of this entanglement measure are investigated.

Quantum entanglement is one of the key ingredients in various types of quantum information processing. A notable example is dense coding , which doubles the capacity of transmission of classical information assisted by an maximally entangled pair of qubits shared between the sender (Alice) and receiver (Bob). Several authors have studied the capacity of dense coding in various situations . In this paper, the author derives an entanglement measure for any bipartite states in the light of the capacity of dense coding and investigates its properties systematically. First, the general scheme of dense coding with a mixed state on the Hilbert space $`C^dC^d`$ is described. Alice performs one of the local unitary transformations $`U_iU(d)`$ on her $`d`$-dimensional quantum system in order to put the initially shared entangled state $`\rho `$ in $`\rho _i=(U_i𝐈_d)\rho (U_i^{}𝐈_d)`$ with a priori probability $`p_i`$ $`(i=1,2,\mathrm{},i_{\mathrm{max}})`$, and then she sends her quantum system to Bob. Upon receiving this quantum system, Bob performs a suitable measurement on $`\rho _i`$ to extract the signal. The optimal amount of information that can be conveyed is known to be bounded from above by the Holevo quantity ,

$$\chi =S(\overline{\rho })\underset{i=1}{\overset{i_{\mathrm{max}}}{}}p_iS(\rho _i),$$
(1)
where $`S(\rho )=\mathrm{Tr}(\rho \mathrm{log}_2\rho )`$ denotes the von Neumann entropy and $`\overline{\rho }=_{i=1}^{i_{\mathrm{max}}}p_i\rho _i`$ is the average density matrix of the signal ensemble. Since the Holevo quantity is asymptotically achievable , Eq. (1) is used here as the definition of the capacity of dense coding. Capacity $`\chi `$ is maximized for signal states $`\rho _i`$ with mutually orthogonal unitary transformations, $`\mathrm{Tr}\left(U_i^{}U_j\right)=d\delta _{ij}`$ and equal probabilities $`p_i=d^2`$ ($`i_{\mathrm{max}}=d^2`$) . The optimal capacity is written as $`\chi ^{}(\rho )=\mathrm{log}_2d+I_B(\rho )`$, where $`I_B(\rho )=S(\rho ^B)S(\rho )`$ is the coherent information with $`\rho ^B=\mathrm{Tr}_A\rho `$. Since $`\mathrm{max}[S(\rho ^A)S(\rho ),S(\rho ^B)S(\rho )]E_R(\rho )`$ ,

$$I_B(\rho )E_R(\rho ),$$
(2)
and the capacity $`\chi ^{}(\rho )`$ satisfies $`\chi ^{}(\rho )\mathrm{log}_2d+E_R(\rho )`$ . Here, $`E_R(\rho )`$ is the relative entropy of entanglement  for states $`\rho `$ defined as $`E_R(\rho )=\mathrm{min}_{\sigma 𝒟}S(\rho ||\sigma )`$, where $`𝒟`$ the set of states with positive partial transpose (PPT states) and $`S(\rho ||\sigma )=\mathrm{Tr}\left[\rho (\mathrm{log}_2\rho \mathrm{log}_2\sigma )\right]`$ is the quantum relative entropy of $`\rho `$ with respect to $`\sigma `$.

Note that this capacity is optimal in the sense that Alice and Bob uses a given mixed state $`\rho `$ as a resource for dense coding without any changes. If they are allowed to perform local quantum operations and classical communications (LQCC) on the initially shared mixed state $`\rho `$ prior to the dense coding, the capacity could be enhanced further. The maximally achievable capacity thus obtained  is

$$\chi _{\mathrm{max}}^{}(\rho )=\mathrm{log}_2d+\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n}{sup}\frac{1}{n}I_B\left(\mathrm{\Lambda }_n(\rho ^n)\right)\mathrm{log}_2d+E_{dc}(\rho ).$$
(3)
Namely, $`E_{dc}(\rho )`$ is the asymptotic limit of the achievable (normalized) coherent information over the sequence of LQCC operations $`\left\{\mathrm{\Lambda }_n\right\}`$ or an LQCC protocol.

Hereafter the properties of $`E_{dc}(\rho )`$ defined in Eq. (3) are examined. $`E_{dc}(\rho )`$ is the maximal dense coding capacity subtracted by the classically achievable capacity $`\mathrm{log}_2d`$; it represents the maximal contribution of entanglement to the dense coding capacity. As shown in the following, $`E_{dc}(\rho )`$ is a measure of entanglement of $`\rho `$. Before proving this, the following inequalities must be proved.

$$E_D(\rho )E_{dc}(\rho )E_R^{\mathrm{}}(\rho ),$$
(4)
where $`E_D(\rho )`$ and $`E_R^{\mathrm{}}(\rho )`$ are, respectively, the distillable entanglement  and the asymptotic relative entropy of entanglement  of $`\rho `$, both of which are asymptotic entanglement measures. $`E_R^{\mathrm{}}(\rho )`$ is defined as the average relative entropy of entanglement per copy:

$$E_R^{\mathrm{}}(\rho )=\underset{n\mathrm{}}{lim}\frac{E_R(\rho ^n)}{n}.$$
(5)
Noting the subadditivity of the relative entropy of entanglement, i.e., $`E_R(\rho ^n)nE_R(\rho )`$ , a weaker version of Eq. (4) is obtained:

$$E_D(\rho )E_{dc}(\rho )E_R(\rho ).$$
(6)
Although the proof of the first part of Eq. (4) is essentially the same as that in , the proof is described here for completeness. It is always possible to consider that the distillation protocol is ended by $`UU^{}`$ twirling  so that the final state is an isotropic state of the form

$$\rho (F_n,d_n)=pP_+(C^{d_n})+(1p)\frac{1}{d_n^2}𝐈_{d_n},$$
(7)
where $`F_n=\mathrm{Tr}\left[\rho (F_n,d_n)P_+(C^{d_n})\right]`$ is the fidelity, $`𝐈_{d_n}`$ is the identity of dimensions $`d_n`$, and

$$P_+(C^d)=|\psi _+(C^d)\psi _+(C^d)|$$
(8)
is the maximally entangled state of a $`C^dC^d`$ system. In Eq. (8), $`|\psi _+(C^d)=\frac{1}{\sqrt{d}}_{i=1}^d|ii`$ with $`|i`$ are basis vectors in $`C^d`$. Because the protocol mentioned above is not necessarily optimal for $`E_{dc}`$, $`E_{dc}(\rho )lim_n\mathrm{}\frac{1}{n}I_B\left(\rho (F_n,d_n)\right)`$. The coherent information for $`\rho (F_n,d_n)`$ is easily calculated as

$$I_B\left(\rho (F_n,d_n)\right)=\mathrm{log}_2d_n+F_n\mathrm{log}_2F_n+(1F_n)\mathrm{log}_2\frac{1F_n}{d_n^21}.$$
(9)
By definition of the distillable entanglement , $`F_n1`$ and $`\frac{\mathrm{log}_2d_n}{n}E_D(\rho )`$ for $`n\mathrm{}`$. Therefore, $`E_{dc}(\rho )E_D(\rho )`$. The proof of the second part of Eq. (4) is as follows. Equation (2) and the weak monotonicity (see below) of $`E_R`$  give

$$E_{dc}(\rho )\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n}{sup}\frac{1}{n}E_R\left(\mathrm{\Lambda }_n(\rho ^n)\right)\underset{n\mathrm{}}{lim}\frac{1}{n}E_R(\rho ^n).$$
(10)
The right-hand side is, by definition, $`E_R^{\mathrm{}}(\rho )`$. Therefore, $`E_{dc}(\rho )E_R^{\mathrm{}}(\rho )`$.

The quantity thus defined is an entanglement measure; namely, it satisfies the following properties .

* $`E_{dc}(\rho )=0`$ for any separable state $`\rho `$.
* $`E_{dc}(\rho )0`$.
* For a pure state $`|\varphi \varphi |`$, $`E_{dc}`$ is the von Neumann entropy of the reduced density matrix, e.g., the entropy of entanglement;

$$E_{dc}\left(|\varphi \varphi |\right)=S\left(\mathrm{Tr}_A\left(|\varphi \varphi |\right)\right)=S\left(\mathrm{Tr}_B\left(|\varphi \varphi |\right)\right).$$

In particular, $`E_{dc}\left(P_+(C^d)\right)=\mathrm{log}_2d`$, where $`P_+(C^d)`$ is the maximally entangled state of a $`C^dC^d`$ system \[Eq. (8)\].
* Partial additivity: $`E_{dc}(\rho ^n)=nE_{dc}(\rho )`$.
* Weak monotonicity: $`E_{dc}\left(\mathrm{\Lambda }(\rho )\right)`$ $`E_{dc}(\rho )`$, where $`\mathrm{\Lambda }`$ is an LQCC operation. This is the most important property required of the entanglement measure.
* Convexity on pure state decomposition:

$$E_{dc}\left(\underset{i}{}p_i|\varphi _i\varphi _i|\right)\underset{i}{}p_iE_{dc}\left(|\varphi _i\varphi _i|\right),$$

with $`_ip_i=1`$ and $`p_i0`$
* Weak continuity: For any sequence of the pure state $`|\psi _n`$ and the mixed state $`\rho _n`$ of a system $`C^{d_n}C^{d_n}`$ such that $`\rho _n|\psi _n\psi _n|_10`$ and $`d_n\mathrm{}`$ for $`n\mathrm{}`$,

$$\underset{n\mathrm{}}{lim}\frac{E_{dc}(\rho _n)E_{dc}(|\psi _n\psi _n|)}{\mathrm{log}d_n}=0.$$

Here, $`A_1`$ denotes the trace norm of $`A`$; $`A_1=\mathrm{Tr}\sqrt{A^{}A}`$.

Properties (i)-(iii) are obvious from Eq. (6). The proof of property (iv) is as follows.

$`E_{dc}(\rho ^m)`$ $`=`$ $`\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_{nm}}{sup}{\displaystyle \frac{1}{n}}I_B\left(\mathrm{\Lambda }_{nm}(\rho ^{nm})\right)`$ (11)
$`=`$ $`m\underset{n^{}\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n^{}}{sup}{\displaystyle \frac{1}{n}}I_B\left(\mathrm{\Lambda }_n^{}(\rho ^n^{})\right)=mE_{dc}(\rho ).`$
Property (v) is obvious because $`E_{dc}(\rho )`$ is the optimized quantity with respect to LQCC protocols by definition and the tensor product of an LQCC operation is also an LQCC operation. Property (vi) follows from Eq. (6) and the fact that both $`E_{dc}(\rho )`$ and $`E_R(\rho )`$ coincide on pure states;

$`E_{dc}\left({\displaystyle \underset{i}{}}p_i|\varphi _i\varphi _i|\right)`$ $``$ $`E_R\left({\displaystyle \underset{i}{}}p_i|\varphi _i\varphi _i|\right)`$ (12)
$``$ $`{\displaystyle \underset{i}{}}p_iE_R\left(|\varphi _i\varphi _i|\right)={\displaystyle \underset{i}{}}p_iE_{dc}\left(|\varphi _i\varphi _i|\right).`$
The proof of property (vii) is given as follows. Noting the fact that $`I_B`$, $`E_{dc}`$, and $`E_R`$ coincide on pure states, the inequalities, $`I_B(\rho _n)E_{dc}(\rho _n)E_R(\rho _n)`$ give

$`I_B(\rho _n)I_B\left(|\psi _n\psi _n|\right)`$ $``$ $`E_{dc}(\rho _n)E_{dc}\left(|\psi _n\psi _n|\right)`$ (13)
$``$ $`E_R(\rho _n)E_R\left(|\psi _n\psi _n|\right).`$
Firstly,

$$\underset{n\mathrm{}}{lim}\frac{E_R(\rho _n)E_R(|\psi _n\psi _n|)}{\mathrm{log}d_n}=0,$$
(14)
because $`E_R`$ is continuous. Secondly, Fannes’ inequality ,

$$\left|S(\rho )S(\sigma )\right|\rho \sigma _1\mathrm{log}_2dim+\eta \left(\rho \sigma _1\right),$$
(15)
plays a key role. It holds for any two states $`\rho `$ and $`\sigma `$ acting on the Hilbert space $``$ provided that $`\rho \sigma _11/e`$. In Eq. (15), $`\eta (s)=s\mathrm{log}_2s`$. Noting the fact that the partial trace does not increase the trace norm and $`\eta (s)`$ is a monotonically increasing function for $`0s1/e`$, Fannes’ inequality \[Eq. (15)\] gives

$`\left|I_B(\rho _n)I_B\left(|\psi _n\psi _n|\right)\right|`$ $``$ $`3\rho _n|\psi _n\psi _n|_1\mathrm{log}_2d_n`$ (16)
$`+2\eta \left(\rho _n|\psi _n\psi _n|_1\right),`$
Therefore,

$$\underset{n\mathrm{}}{lim}\frac{\left|I_B(\rho _n)I_B\left(|\psi _n\psi _n|\right)\right|}{\mathrm{log}d_n}=0.$$
(17)
From Eqs. (13), (14), and (17), the following equation is obtained:

$$\underset{n\mathrm{}}{lim}\frac{E_{dc}(\rho _n)E_{dc}(|\psi _n\psi _n|)}{\mathrm{log}d_n}=0.$$
(18)
Namely, $`E_{dc}`$ is weakly continuous.

In addition to properties (i)-(vii), $`E_{dc}(\rho )`$ exhibits superadditivity. Namely, $`E_{dc}(\rho \sigma )E_{dc}(\rho )+E_{dc}(\sigma )`$. The proof is as follows. Because of the additivity of the coherent information,

$$E_{dc}(\rho )+E_{dc}(\sigma )=\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n^\rho \mathrm{\Lambda }_n^\sigma }{sup}\frac{1}{n}I_B\left(\left(\mathrm{\Lambda }_n^\rho \mathrm{\Lambda }_n^\sigma \right)(\rho ^n\sigma ^n)\right).$$
(19)
Here, even if the protocol $`\left\{\mathrm{\Lambda }_n^\rho \mathrm{\Lambda }_n^\sigma \right\}`$ is optimized, it is not necessarily the optimal one for $`\left\{(\rho \sigma )^n\right\}`$. Therefore,

$$E_{dc}(\rho )+E_{dc}(\sigma )\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n^{\rho \sigma }}{sup}\frac{1}{n}I_B\left(\mathrm{\Lambda }_n^{\rho \sigma }(\rho \sigma )^n\right)=E_{dc}(\rho \sigma ).$$
(20)
It is not clear at present if the equality (full additivity) holds. The convexity of the general form,

$$E_{dc}\left(\underset{i}{}p_i\rho _i\right)\underset{i}{}p_iE_{dc}(\rho _i),$$
(21)
is also doubtful. However, it should be noted that the breakdown of the full additivity and the general convexity is not a drawback; it is argued that these two requirements are too strong for asymptotic entanglement measures .

Although it is in general quite difficult to calculate $`E_{dc}(\rho )`$, there are some special mixed states in which $`E_{dc}(\rho )`$ is obtained explicitly.

Example 1—This is the example by Rains . It is called the maximally correlated state of a $`C^dC^d`$ system, and takes the form

$$\rho =\underset{i,j=1}{\overset{d}{}}\alpha _{ij}|iijj|.$$
(22)
The relative entropy of entanglement is calculated as

$$E_R(\rho )=I_B(\rho )=H(\alpha _{11},\alpha _{22},\mathrm{})S(\rho ),$$
(23)
where $`H(\alpha _{11},\alpha _{22},\mathrm{})`$ denotes the Shannon entropy of probability distribution $`\left\{\alpha _{ii}\right\}`$. From $`E_R(\rho )=I_B(\rho )`$, $`E_{dc}(\rho )=E_R^{\mathrm{}}(\rho )=E_R(\rho )`$, which is proved as follows:

$`E_R(\rho )E_{dc}(\rho )`$ $``$ $`\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n}{sup}{\displaystyle \frac{1}{n}}E_R\left(\mathrm{\Lambda }_n(\rho ^n)\right)`$ (24)
$``$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{1}{n}}E_R(\rho ^n)=E_R^{\mathrm{}}(\rho )E_R(\rho ).`$
The first inequality is obvious because $`I_B(\rho )E_{dc}(\rho )`$. The second inequality is a consequence of Eq. (2), and the third inequality follows from the weak monotonicity of $`E_R`$. The last inequality is a result of the subadditivity of $`E_R`$. The optimal LQCC operation for $`E_{dc}`$ is simply $`\mathrm{\Lambda }_n=I_{d^n}`$. It has been shown that $`E_R(\rho )`$ is exactly the PPT distillable entanglement (distillable entanglement with respect to positive partial transpose operations ). Since the set of LQCC operations is a subset of the set of PPT operations, $`E_D(\rho )`$ is the lower bound on the PPT distillable entanglement. Therefore, $`E_{dc}(\rho )=E_R^{\mathrm{}}(\rho )=E_R(\rho )=I_B(\rho )E_D(\rho )`$. When $`d=2`$, the maximally correlated state is equivalent to a mixture of two Bell states (a Bell diagonal state of rank two) if $`\mathrm{Tr}_A\rho =𝐈_2/2`$. For this state the hashing protocol of distillation yields the value of $`E_R(\rho )=1S(\rho )`$ for the distillable entanglement  so that $`E_{dc}(\rho )=E_D(\rho )`$.

Example 2—This is the example by Eisert et al. . Suppose that Alice and Bob share initially $`N=2J`$ ($`J=1,2,\mathrm{}`$) pair of qubits each in the same state $`|\varphi =\alpha |00+\beta |11`$. Hereafter in this example $`\alpha =\beta =1/\sqrt{2}`$ is assumed for simplicity. After the complete loss of the order of Bob’s particles, the initially shared pure state $`|\varphi ^N`$ becomes a mixed state of the form

$$\rho =\underset{j=0}{\overset{J}{}}\underset{\alpha _j,\beta _j=1}{\overset{d_j}{}}p_j|\psi _j(\alpha _j,\beta _j)\psi _j(\alpha _j,\beta _j)|,$$
(25)
where

$$|\psi _j(\alpha _j,\beta _j)=\frac{1}{\sqrt{2j+1}}\underset{m=j}{\overset{j}{}}|j,m,\alpha _j|j,m,\beta _j,$$
(26)
$`p_j=(2j+1)/(d_j2^{2J})`$, and $`d_j=\frac{2j+1}{2J+1}\left(\genfrac{}{}{0pt}{}{2J+1}{Jj}\right)`$ is the multiplicity of the $`j`$-representation in $`SU(2)^N`$. It is easy to calculate the coherent information;

$$I_B(\rho )=\underset{j=0}{\overset{J}{}}d_j^2p_j\left[\mathrm{log}_2(2j+1)\mathrm{log}_2d_j\right].$$
(27)
On the other hand, the relative entropy of entanglement and the distillable entanglement are calculated as

$$E_R(\rho )=E_D(\rho )=\underset{j=0}{\overset{J}{}}d_j^2p_j\mathrm{log}_2(2j+1)$$
(28)
so that $`E_{dc}(\rho )=E_R^{\mathrm{}}(\rho )=E_R(\rho )=E_D(\rho )I_B(\rho )`$. The first three equalities follow from Eqs. (4) and (6) and the subadditivity of $`E_R`$. The last equality holds only for $`J=1`$ ($`d_0=d_1=1`$). The optimal distillation is the optimal LQCC protocol for $`E_{dc}`$.

Two examples described above show that it is reasonable to conjecture that the optimal protocol for $`E_{dc}(\rho )`$ is either the identity \[$`E_{dc}(\rho )=I_B(\rho )E_D(\rho )]`$ or the optimal distillation protocols \[$`E_{dc}(\rho )=E_D(\rho )I_B(\rho )]`$. However, Horodecki et al. have conjectured that $`I_B(\rho )`$ does not exceed the one-way distillable entanglement (distillable entanglement with local operations plus one-way classical communications) for any state $`\rho `$ . If this conjecture (hashing inequality) is true,

$$E_D(\rho )I_B(\rho )$$
(29)
holds for any state $`\rho `$. This inequality implies that $`E_D`$ is weakly continuous, which is not proved yet. The proof of weak continuity follows from the same arguments of the proof of property (vii) of $`E_{dc}`$. Equation (29) also implies

$$E_{dc}(\rho )=E_D(\rho ).$$
(30)
The proof of Eq. (30) is simple and essentially the same as that in . The partial additivity and the weak monotonicity of $`E_D`$  give

$`E_D(\rho )=\underset{n\mathrm{}}{lim}{\displaystyle \frac{1}{n}}E_D(\rho ^n)`$ $``$ $`\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n}{sup}{\displaystyle \frac{1}{n}}E_D\left(\mathrm{\Lambda }_n(\rho ^n)\right)`$ (31)
$``$ $`\underset{n\mathrm{}}{lim}\underset{\mathrm{\Lambda }_n}{sup}{\displaystyle \frac{1}{n}}I_B\left(\mathrm{\Lambda }_n(\rho ^n)\right)=E_{dc}(\rho ).`$
From Eqs. (4) and (31), Eq. (30) is obtained. This is a satisfactory result. It strengthens the information-theoretic meaning of the distillable entanglement; namely, $`E_D`$ is the ultimate measure of resources for dense coding. Furthermore, the optimal entanglement distillation seems to be the best strategy to maximize the coherent information since it increases $`S\left(\mathrm{Tr}_A\left(\mathrm{\Lambda }_n(\rho ^n)\right)\right)`$ on one hand and decreases $`S\left(\mathrm{\Lambda }_n(\rho ^n)\right)`$ on the other hand while keeping the dimension of $`\mathrm{\Lambda }_n(\rho ^n)`$ as large as possible. According to the above reasonings, it is most likely that $`E_{dc}(\rho )=E_D(\rho )`$. Unfortunately, the assumed inequality $`E_D(\rho )I_B(\rho )`$, which is also a consequence of the equality $`E_{dc}(\rho )=E_D(\rho )`$ \[Eq. (30)\], is not proven yet. One of the possible counter-examples is Example 1. However, Rains has conjectured that for any maximally correlated state $`\rho `$ both the PPT distillable entanglement and the one-way distillable entanglement coincide , so $`E_D(\rho )=I_B(\rho )`$. It should be noted that this conjecture is also a consequence of the hypothetical hashing inequality.

In summary, in the light of the dense coding capacity optimized with respect to LQCC, an asymptotic entanglement measure $`E_{dc}`$ for any bipartite states was derived and its properties was investigated. Some examples of $`E_{dc}`$ with explicit forms were also given. Finally, it was argued that $`E_{dc}`$ coincides with the distillable entanglement. A possible counter-example for this conjecture was also given.

The author would like to thank Andreas Winter for helpful comments. This work was supported by CREST of Japan Science and Technology Corporation (JST).


format: text

text:
# Simple solutions of relativistic hydrodynamics for cylindrically symmetric systems

## 1 Introduction

The analytic resolution of 3 dimensional relativistic hydrodynamics is a difficult task because the equations are non-linear and the transverse and the longitudinal motions of the fluid are coupled in a way that is rather complicated to handle.

Historically, the first solutions of relativistic hydrodynamics were found for idealized, one dimensionally expanding sources: Landau’s one-dimensional analytical solution of relativistic hydrodynamics  and Bjorken’s boost-invariant solution , that correspond to ultrarelativistic bombarding energies. Both of these solutions are frequently utilized for various estimations of observables in ultra-relativistic nucleus-nucleus collisions.

The repertoir of relativistic hydrodynamics has been reviewed from the point of view of heavy ion physics in ref. . Recently, Biró has found self-similar exact solutions of relativistic hydrodynamics for cylindrically expanding systems . Bjorken’s one-dimensional solution can be naturally extended to a three-dimensional solution, assuming spherical symmetry instead of (1+1) dimensional longitudinal dynamics. Biró’s solutions generalize the (1+3) dimensional version of Bjorken’s solution to the case of cylindrical symmetry, interpolating between the one-dimensional and the three-dimensional Bjorken solution. However, Biró’s solutions are valid only when the pressure is independent of space and time, as happens at the softest point of the equation of state.

The numerical solution of the 1+3 dimensional relativistic hydrodynamics is possible with presently available computers and methods, for example the variational method of smoothed particle hydrodynamics has been developed to study event-by-event fluctuations in relativistic heavy-ion collisions .

Here we describe a new family of exact analytic solutions of relativistic hydrodynamics with a broad family of equations of state, assuming cylindrical symmetry of the particle emitting source. The physical motivation for this study is to consider the time evolution of central collisions in ultra-relativistic heavy-ion physics within the framework of an analytic approach.

Our solution was motivated by the analytic solution of non-relativistic hydrodynamics found by Zimányi, Bondorf and Garpman (ZBG) in 1978 for low-energy heavy-ion collisions with spherical symmetry . This solution has been generalized to the case of ellipsoidal symmetry by Zimányi and collaborators in ref. . In ref.  a Gaussian parameterization has been introduced to describe the mass dependence of the effective temperature and the radius parameters of the two-particle Bose-Einstein correlation functions in high-energy heavy-ion collisions. Later it has been realized that this phenomenological parameterization of data corresponds to an exact, Gaussian solution of non-relativistic hydrodynamics with spherical symmetry . The spherically symmetric self-similar solutions of non-relativistic hydrodynamics were generalized in ref.  in terms of an arbitrary scaling function for the temperature profile, and expressing the density distribution in terms of this scaling function. The spherically symmetric Gaussian solution has been generalized to ellipsoidally symmetric expansions , that provided simple analytic insight into the observables of non-central heavy-ion collisions . The family of ellipsoidally symmetric solutions of non-relativistic hydrodynamics has been expressed in terms of a general scaling function for the temperature profile in ref. .

Our approach corresponds to a generalization of these recently obtained analytic solutions  of non-relativistic fireball hydrodynamics to the case of relativistic longitudinal and transverse flows. In particular, an analytic approach, the Buda-Lund (BL) model has been developed to parameterize the single particle spectra and the two-particle Bose-Einstein correlations in high-energy heavy-ion physics in terms of hydrodynamically expanding, cylindrically symmetric sources . Here we attempt to find an exact family of solutions of relativistic hydrodynamics that may include the BL model as a particular limiting case.

We attempt here to solve the full (1+3) dimensional relativistic hydrodynamical problem, in trying to overcome two shortcomings of Bjorken’s well known solution. These two shortcomings of Bjorken’s solutions are that i) it describes ultra-relativistic collisions, it is independent of scale parameter in the longitudinal direction and the rapidity distribution is flat; ii) it contains no transverse flow.

We have described recently a new family of (1+1) dimensional solutions  that is able to describe an arbitrary rapidity distribution $`dn/dy`$ by introducing scale-dependent quantities in the longitudinal direction, overcoming shortcoming $`i)`$, but not addressing shortcoming $`\mathrm{𝑖𝑖})`$. Here we consider the opposite case, overcoming shortcoming ii), the development of a scale dependent transverse flow in relativistic hydrodynamics, but not addressing shortcoming i). We address both shortcomings simultaneously in a subsequent manuscript .

## 2 The equations of relativistic hydrodynamics

We solve the relativistic continuity and energy-momentum conservation equation:

$`_\mu (nu^\mu )`$ $`=`$ $`0,`$ (1)
$`_\nu T^{\mu \nu }`$ $`=`$ $`0.`$ (2)
Here the $`nn(t,𝐫)`$ is the number density, the four-velocity is denoted by $`u^\mu u^\mu (t,𝐫)=\gamma (1,𝐯)`$, normalized to $`u^\mu u_\mu =\gamma ^2(1𝐯^2)=1`$, and the energy-momentum tensor is denoted by $`T^{\mu \nu }`$. We assume perfect fluid,

$$T^{\mu \nu }=(ϵ+p)u^\mu u^\nu pg^{\mu \nu },$$
(3)
where $`ϵ`$ stands for the relativistic energy density and $`p`$ denotes the pressure.

We close this set of relativistic hydrodynamical equations with the equations of state. We assume a gas containing massive conserved quanta,

$`ϵ`$ $`=`$ $`mn+\kappa p,`$ (4)
$`p`$ $`=`$ $`nT.`$ (5)
The equations of state have two free parameters, $`m`$ and $`\kappa `$. Non-relativistic hydrodynamics of ideal gases corresponds to the limiting case of $`mT`$, $`𝐯^21`$ and $`\kappa =\frac{3}{2}`$.

The energy-momentum conservation equations can be projected into a component parallel to $`u^\mu `$ and components orthogonal to $`u^\mu `$, which are respectively the relativistic energy and Euler equations:

$`u^\mu _\mu ϵ+(ϵ+p)_\mu u^\mu `$ $`=`$ $`0,`$ (6)
$`u_\nu u^\mu _\mu p+(ϵ+p)u^\mu _\mu u_\nu _\nu p`$ $`=`$ $`0.`$ (7)
With the help of the equations of state and the continuity equation, the energy equation can be rewritten as an equation for the temperature,

$$u^\mu _\mu T+\frac{1}{\kappa }T_\mu u^\mu =0.$$
(8)
We adopt the following notational conventions: the coordinates are $`x^\mu =(t,𝐫)=(t,r_x,r_y,r_z)`$, $`x_\mu =(t,r_x,r_y,r_z)`$ and the metric tensor is $`g^{\mu \nu }=g_{\mu \nu }=\text{diag}(1,1,1,1)`$.

We solve 5 independent equations, the continuity, the (3 spatial components of) relativistic Euler, and the temperature equation, eqs. (1,7,8). The equations of state, eq. (4,5) close this system of equations in terms of 5 variables, $`n`$, $`T`$ and $`𝐯=(v_x,v_y,v_z)`$.

## 3 Self-similar, cylindrically symmetric solutions

As we are primarily interested in the effects of finite transverse size and the development of transverse flow, we assume that the longitudinal flow component is that of Bjorken’s type,

$$v_z(t,r_z)=\frac{r_z}{t}.$$
(9)
A similar assumption has been made when a new family of (1+1) dimensional solutions of relativistic hydrodynamics was obtained in ref. . However, in contrast to what has been done in , here we assume scale-invariance in the longitudinal direction.

We search for self-similar solutions, that are scale dependent in the transverse directions, and depend only on the transverse radius variable $`r_t=\sqrt{r_x^2+r_y^2}`$ and the longitudinal proper time $`\tau _z=\sqrt{t^2r_z^2}`$. Let us introduce the scaling variable $`x`$ as

$$x=\frac{r_x^2+r_y^2}{R^2},$$
(10)
and assume that, in the frame where $`v_z=0`$ (longitudinal proper frame), the transverse motion corresponds to a Hubble type of self-similar transverse expansion,

$`v_x^{}(\tau _z,r_z)`$ $`=`$ $`{\displaystyle \frac{\dot{R}(\tau _z)}{R(\tau _z)}}r_x,`$ (11)
$`v_y^{}(\tau _z,r_z)`$ $`=`$ $`{\displaystyle \frac{\dot{R}(\tau _z)}{R(\tau _z)}}r_y,`$ (12)
where $`\dot{R}=dR(\tau _z)/d\tau _z`$ and hereafter we will designate by starred symbols the variables in the longitudinal proper frame. We assume that the scale $`R`$ depends on time only trough the longitudinal proper time, $`\tau _z`$.

In a relativistic notation, the above form may be parametrized as

$`u^\mu `$ $`=`$ $`(\mathrm{cosh}\zeta \mathrm{cosh}\xi ,\mathrm{sinh}\xi {\displaystyle \frac{r_x}{r_t}},\mathrm{sinh}\xi {\displaystyle \frac{r_y}{r_t}},\mathrm{sinh}\zeta \mathrm{cosh}\xi ),`$ (13)
$`\mathrm{tanh}\xi `$ $`=`$ $`{\displaystyle \frac{\dot{R}(\tau _z)}{R(\tau _z)}}r_t=v_t^{}=\gamma _lv_t,`$ (14)
$`\mathrm{cosh}\xi `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1\dot{R}^2x}}}\gamma _t^{},`$ (15)
$`\mathrm{cosh}\zeta `$ $`=`$ $`{\displaystyle \frac{t}{\tau _z}}\gamma _l.`$ (16)
The space-time rapidity $`\eta `$ is defined as

$$\eta =0.5\mathrm{log}\left(\frac{t+r_z}{tr_z}\right).$$
(17)
For a scaling longitudinal Bjorken flow we obtain

$$\zeta =\eta .$$
(18)
Using the above ansatz for the flow velocity distribution, we find that the continuity equation is solved by the form

$$n(t,r_x,r_y,r_z)=n_0\left(\frac{\tau _{z0}R_0^2}{\tau _zR^2}\right)\frac{1}{\mathrm{cosh}\xi }𝒢(x),$$
(19)
where $`𝒢(x)`$ is an arbitrary non-negative function of the scaling variable $`x`$ and $`n_0`$, $`\tau _{z0}`$ and $`R_0`$ are normalization constants. We use the convention $`n_0=n(t_0,0,0,0)`$, $`\tau _{z0}=\tau _z(t_0,r_{z0})`$ and $`R_0=R(\tau _{z0})`$, where $`r_{z0}`$ is such that, together with $`t_0`$, satisfies eq. (9). This implies that $`𝒢(x=0)=1`$. The temperature equation, eq. (8), is solved by

$$T(t,r_x,r_y,r_z)=T_0\left(\frac{\tau _{z0}R_0^2}{\tau _zR^2}\frac{1}{\mathrm{cosh}\xi }\right)^{1/\kappa }(x).$$
(20)
The constants of normalization are $`T_0=T(t_0,0,0,0)`$ and $`(0)=1`$. We find that the solution is independent of the form of the function $`(x)`$, provided that $`(x)>0`$.

Using the ansatz for the flow profile and the solution for the density and the temperature profile, the relativistic Euler equation reduces to a complicated non-linear equation that contains $`R`$, $`\dot{R}`$ and $`\ddot{R}`$ and the variable $`x`$. Taking this equation at $`x=0`$ we express $`\ddot{R}`$ as a function of $`R`$ and $`\dot{R}`$. Substituting this back to the Euler equation we obtain a transcendental equation for $`\dot{R}^2`$, and $`x`$. This equation has a particular solution if

$$\dot{R}=\dot{R}_0.$$
(21)
In this case, the acceleration of the radius parameter vanishes, $`\ddot{R}=0`$, and the solution is $`R=R_0+\dot{R}_0(\tau _z\tau _{z0})`$. Then the relativistic Euler equation reduces to

$$\left(1+\frac{1}{\kappa }\right)\left(\frac{R\dot{R}}{\tau _z}+3\dot{R}^2\right)=2(1x\dot{R}^2)\left[\mathrm{log}𝒢(x)(x)\right]^{},$$
(22)
where the lhs depends only on $`\tau _z`$ while the rhs is only a function of the variable $`x`$, hence both sides are constant. This implies that

$$\frac{R}{\tau _z}=\dot{R}_0,$$
(23)
thus $`R_0=\dot{R}_0\tau _{z0}`$. Thus the origin of the time axis (fixed by the assumption of the scaling longitudinal Bjorken flow profile) coincides with the vanishing value of the transverse radius parameters.

The solutions can be casted in a relatively simple form by introducing the proper time $`\tau `$,

$`\tau `$ $`=`$ $`\sqrt{\tau _z^2r_t^2}=\sqrt{t^2r_x^2r_y^2r_z^2}.`$ (24)
Using this natural variable we find that

$`𝐯`$ $`=`$ $`{\displaystyle \frac{𝐫}{t}},`$ (26)
$`u^\mu `$ $`=`$ $`{\displaystyle \frac{x^\mu }{\tau }}.`$ (27)
Thus the velocity field of our solution corresponds to the flow field of the spherically symmetric scaling solution. However, in the scaling solution the temperature and the pressure distributions are dependent only on the proper time $`\tau `$, while in our case both the density and the temperature distributions are generally dependent on the scale variable $`x`$ in the transverse direction.

As the solution is relativistic, and it is defined in the positive light-cone, given by $`\tau 0`$, we obtain a constraint for the transverse coordinate, $`r_t\tau _z`$. This together with eq. (23), the solution for the scale $`R`$, implies that the scaling variable has to satisfy the constraint $`x\dot{R}_0^21`$, which corresponds to the limitation that the velocity of the fluid can not exceed the speed of light.

By replacing eqs. (21,23) into the Euler equation, eq.(22), one obtains

$$\frac{d}{dx}\mathrm{log}\left[(1x\dot{R}_0^2)^{2(1+1/\kappa )}𝒢(x)(x)\right]=0,$$
(28)
which gives, together with the condition $`𝒢(0)(0)=1`$,

$$𝒢(x)(x)=(1\dot{R}_0^2x)^{2(1+1/\kappa )}.$$
(29)
In this family of solutions, the scaling functions for the temperature and the density distribution are thus not independent. However, a constraint is given for their product, hence one of them can be chosen as an arbitrary positive function.

For clarity, let us introduce new forms of the scaling functions as

$`𝒯(x)`$ $`=`$ $`(x)(1\dot{R}_0^2x)^{2/\kappa },`$ (30)
$`𝒱(x)`$ $`=`$ $`𝒢(x)(1\dot{R}_0^2x)^2.`$ (31)
Then the constraint can be casted to the simple form of

$$𝒱(x)𝒯(x)=1.$$
(32)
This construction for the scaling functions of the transverse density and temperature profiles coincides with the method, that we developed for the solution of the relativisitic hydrodynamical equations in the (1+1) dimensional problem, but here the transverse flow has a two-dimensional distribution, so the exponents and the scaling variables had to be re-defined accordingly.

Let us summarize our new family of solutions of the 1+3 dimensional relativistic hydrodynamics for cylindrically symmetric systems by substituting the results to the density, temperature and pressure profiles.

We obtain

$`𝐯`$ $`=`$ $`{\displaystyle \frac{𝐫}{t}},\text{for}|𝐫|t,`$ (33)
$`x`$ $`=`$ $`{\displaystyle \frac{r_t^2}{\dot{R}_0^2\tau _z^2}},\text{for}r_t\tau _z,`$ (34)
$`n(t,𝐫)`$ $`=`$ $`n_0\left({\displaystyle \frac{\tau _{z0}}{\tau }}\right)_0^3𝒱(x),`$ (35)
$`p(t,𝐫)`$ $`=`$ $`p_0\left({\displaystyle \frac{\tau _{z0}}{\tau }}\right)^{3+3/\kappa },`$ (36)
$`T(t,𝐫)`$ $`=`$ $`T_0\left({\displaystyle \frac{\tau _{z0}}{\tau }}\right)^{3/\kappa }{\displaystyle \frac{1}{𝒱(x)}},`$ (37)
where $`p_0=n_0T_0`$. Note that the scaling variable $`x`$ is invariant for boosts in the longitudinal direction, and it is rotation-invariant in the transverse direction, but $`x`$ is not boost-invariant in the transverse directions. Hence we have generated cylindrically symmetric, longitudinally boost invariant solutions of relativistic hydrodynamics. In the longitudinal direction, these solutions are homogeneous, boost-invariant and also scale-invariant. Due to this reason, the observable rapidity distribution is

$$\frac{dn}{dy}=\text{const},$$
(38)
a flat distribution, corresponding to the ultra-relativistic nature of the solution in the longitudinal direction (where $`y=0.5\mathrm{log}[(E+k_z)/(Ek_z)]`$ is the rapidity of a particle with four-momentum $`(E,𝐤)`$ and $`dn/dy`$ is the rapidity distribution of particle density).

A new hydrodynamical solution is assigned to each non-negative function $`𝒱(x)`$, similarly to the cases of the non-relativistic solutions of ref.  and the 1+1 dimensional relativistic solution of ref. . Note that the solutions are valid also for massive particles, the form of the solution is independent of the value of the mass $`m`$. The form of solutions depends parameterically on $`\kappa `$, that characterizes the equation of state.

We have obtained new solutions of the (1+3) dimensional relativistic hydrodynamical equations which describe a self-similar, streaming flow. In the case of $`\dot{R}=1`$ and $`𝒱(x)=1`$ we recover the spherically symmetric scale-invariant solution. This means that, in this limiting case, the pressure, the density and the temperature profiles depend only on the proper time $`\tau `$.

In the general case, our solution depends on the scale $`R`$ (or we can say it contains some characteristic scale $`R`$) and also on an arbitrary scaling function $`𝒱(x)`$.

## 4 Summary

We have found a new family of solutions of 1+3 dimensional relativistic hydrodynamics with cyllindrical symmetry. This family solves the continuity equation and the conservation of the energy - momentum tensor of a perfect fluid, assuming a simple equation of state, similarly to the case of the recently found 1+1 dimensional solutions. The mass of the particles $`m`$ and the constant of proportionality between the kinetic energy density and the pressure, $`\kappa `$, are free parameters of the solution.

As compared to the well-known case of the scale-invariant solution, we have solved one more equation, the continuity equation. We have considered equations of state that have two free parameters, the mass $`m`$ and $`\kappa =ϵ/p`$, while the scale-invariant solution is obtained in the $`m=0`$ approximation. Interestingly, our generalizations resulted in additional freedom in the solution. The new solutions, similarly to the scale-invariant case, prescribe scaling 3-dimensional flow and pressure distribution. However, in our case, the pressure is a product of the local number density and the local temperature, hence one of these can be chosen in an arbitrary manner.

The essential result of our paper is that we found a rich family of exact analytic solutions of relativistic hydrodynamics that contain both a longitudinal Bjorken flow (that is frequently utilized in estimations of observables in high energy heavy ion collisions) and a relativistic transverse flow (whose existence is evident from the analysis of the single particle spectra at RHIC and SPS energies).

Acknowledgments: T. Cs. would like to thank L.P. Csernai, B. Lukács and J. Zimányi for inspiring discussions during the initial phase of this work, and to Y. Hama and G. Krein for kind hospitality during his stay at USP and IFT, Sao Paulo, Brazil. This work has been supported by a NATO Science Fellowship (Cs.T.), by the OTKA grants T026435, T029158 and T034269 of Hungary, the NWO - OTKA grant N 25487 of The Netherlands and Hungary, and the grants FAPESP 00/04422-7, 99/09113-3 of São Paulo, Brazil.


text:
# Critical phenomena in networks

## Abstract

We develop a phenomenological theory of critical phenomena in networks with an arbitrary distribution of connections $`P(k)`$. The theory shows that the critical behavior depends in a crucial way on the form of $`P(k)`$ and differs strongly from the standard mean-field behavior. The critical behavior observed in various networks is analyzed and found to be in agreement with the theory.

05.10.-a, 05.40.-a, 05.50.+q, 87.18.Sn

Many complex, interacting systems such as the brain, Internet, social systems, etc. are recognized as networks. They demonstrate a spectrum of unique effects , which distinguish them from all other known structures of condensed matter. One could expect that such systems have peculiar cooperative phenomena. Only the first studies of particular cooperative models on particular networks were performed recently . They were focused on traditional models in statistical physics, such as the Ising and $`XY`$ models , percolation , epidemic spreading, et al. in scale-free networks . It was revealed that their critical behavior is more rich and extremely far from that expected from the standard mean-field theory. A deviation from the usual mean-field critical behavior appears when, depending on the model, the fourth, $`k^4`$, or third, $`k^3`$, moments of the degree distribution $`P(k)`$ diverge. Here, degree is the number of connections of a vertex. When $`k^2`$ diverges, all models undergo an unusual phase transition of infinite order. This case is of primary importance because many real networks such as the Internet and biological nets are described by $`P(k)`$ with infinite $`k^2`$.

Why do critical phenomena in networks differ so much from those in usual substrates and what is their common origin? Why do all investigated models demonstrate universal behavior when $`k^2`$ diverges? In order to answer the questions raised above and analyze results of simulations and experiments from a general point of view, it is necessary to have a general theory which is not restricted by specific properties of any model.

The advantage of the phenomenological approach is that the origin of interactions and nature of interacting objects are irrelevant. These can be spins, percolating clusters, biological objects and etc. In this Letter, we develop a phenomenological theory of cooperative phenomena in networks. Our approach is based on concepts of the Landau theory of continuous phase transitions. It is shown that the critical behavior in networks has a universal character and is determined by two general properties: (i) the structure of a network and (ii) the symmetry underlying a model. We find that in networks described by a degree distribution with divergent moments, the thermodynamic potential of an interacting system with cooperative effects is a singular function of an order parameter. The theory is in complete agreement with exact results for the Ising and Potts models and percolation on ‘equilibrium’ uncorrelated random networks. On the basis of the theory we discuss results of theoretical studies and simulations of the Ising and $`XY`$ models and epidemic spreading on evolving and small-world networks.

Let us consider a system of interacting objects. Interactions or links between these objects form a net. We assume that some kind of ‘order’ can emerge. This ‘ordered’ phase may be characterized by some quantitative characteristic $`x`$ while it will vanish in a ‘disordered’ phase above a critical point. In order to study the critical behavior we assume that the thermodynamic potential $`\mathrm{\Phi }`$ of the system is not only a function of the order parameter $`x`$ but also depends on the degree distribution:

$$\mathrm{\Phi }(x,H)=Hx+\underset{k}{\overset{\mathrm{}}{}}P(k)\varphi (x,kx).$$
(1)
Here $`H`$ is a field conjugated with $`x`$. It should be emphasized that Eq. (1) is not obvious a priori. The function $`\varphi (x,kx)`$ may be considered as a contribution of vertices with $`k`$ connections. There are arguments in favor of this assumption. Let us consider the interaction of an arbitrary vertex $`0`$ with $`k`$ neighboring vertices. In the framework of a mean-field approach, $`k`$ neighboring sites having a spontaneous ‘moment’ $`x`$ produce an effective field $`kx`$ acting on the site $`0`$. This indicates that the expansion is over not only $`x`$ but also $`kx`$.

We assume that $`\varphi (x,y)`$ is a smooth function of $`x`$ and $`y`$ and can be represented as a series in powers of both $`x`$ and $`y`$:

$$\varphi (x,y)=\underset{m=0}{\overset{\mathrm{}}{}}\underset{l=0}{\overset{\mathrm{}}{}}\varphi _{ml}x^my^l,$$
(2)
where $`\varphi _{ml}`$ are functions of ‘temperature’ $`T`$ and $`H`$. We use the term ‘temperature’ for convenience. If a change of another parameter leads to the emergence of ‘order’, then this control parameter can be used as $`T`$.

There are general restrictions on the coefficients $`\varphi _{ml}`$. For an arbitrary $`P(k)`$, at zero field $`H=0`$, the expansion of $`\mathrm{\Phi }(x,H=0)`$ over $`x`$ must contain only the second and higher powers of $`x`$. Therefore, $`\varphi _{00}=\varphi _{01}=\varphi _{10}=0`$.

$`\mathrm{\Phi }`$ must be finite if $`k`$ is finite. This condition is satisfied if at $`y1`$ and an arbitrary $`x`$ the function $`\varphi (x,y)`$ increases slower than $`y`$:

$$\varphi (x,y)g(x)y,$$
(3)
where $`g(x)`$ is a function of $`x`$.

Eqs. (1)–(3) are the basis of our theory.

In the framework of the Landau theory, it is assumed that the thermodynamic potential $`\mathrm{\Phi }`$ near a phase transition point can be represented as a series in powers of the order parameter $`x`$:

$$\mathrm{\Phi }(x,H)=Hx+f_2x^2+f_3x^3+f_4x^4+\mathrm{}.$$
(4)
Higher terms in the expansion (4) are irrelevant for the critical behavior. The coefficients $`f_n`$ are related to the coefficients $`\varphi _{ml}`$ in the expansion (2). Let us suppose first that the critical temperature $`T_c`$ is finite. (The case of an ‘infinite’ critical temperature will be considered separately.) Following Landau, we assume that near the critical temperature the coefficient $`f_2`$ can be written as $`a(TT_c)`$ where $`a>0`$ due to the condition of stability of the disordered phase at $`T>T_c`$. $`f_2=0`$ at the transition point. This gives the following equation for $`T_c`$:

$$\varphi _{20}(T_c)+\varphi _{11}(T_c)k+\varphi _{02}(T_c)k^2=0.$$
(5)
According to this equation, in the case $`\varphi _{20}=0`$ the critical temperature $`T_c`$ depends only on the ratio $`k^2/k`$. Such a dependence was revealed in all studied models on networks . In the ordered phase the order parameter $`x`$ is determined by the condition that $`\mathrm{\Phi }(x,H)`$ is minimum:

$$d\mathrm{\Phi }(x,H)/dx=0.$$
(6)
Solving Eq. (6), we find $`x(T,H)`$ and all other physical parameters: the response function $`\chi _x=dx/dH`$, the specific heat $`C=Td^2\mathrm{\Phi }/dT^2`$, et. al.

The condition of the stability of the ordered phase near the continuous phase transition demands that either $`f_3>0`$ or if $`f_3=0,`$ then $`f_4>0`$. In the general case, $`f_{2n+1}0`$. If for symmetry reasons $`\mathrm{\Phi }(x,H)=\mathrm{\Phi }(x,H)`$, then $`f_{2n+1}=0`$. This condition takes place if $`\varphi _{ml}=0`$ at $`m+l=2n+1`$.

Let us consider analytical properties of $`\mathrm{\Phi }`$ in the general case. The $`n`$-th derivative of $`\mathrm{\Phi }(x,H)`$ at $`x=0`$ is equal to

$$\mathrm{\Phi }^{(n)}(0)d^n\mathrm{\Phi }(x,H)/dx^n|_{x=0}=n!\underset{l=0}{\overset{n}{}}\varphi _{nl,l}k^l.$$
(7)
If all the moments $`k^l=_kP(k)k^l`$ converge, then all the derivatives $`\mathrm{\Phi }^{(n)}(0)`$ are finite. Therefore, $`\mathrm{\Phi }(x,H)`$ is a smooth function of $`x`$, at least at small $`x`$. However, in most interesting real networks, $`P(k)`$ has divergent moments. This leads to a very important consequence. If moments $`k^l`$ with $`l<p`$ converge and moments $`k^l`$ with $`l`$ $`p`$ diverge, then $`\mathrm{\Phi }^{(l)}(0)`$ with $`l<p`$ are finite while $`\mathrm{\Phi }^{(l)}(0)\mathrm{}`$ for $`lp`$. In this case $`\mathrm{\Phi }(x,H)`$ is a singular function of $`x`$ and can be written as follows:

$$\mathrm{\Phi }(x,H)=Hx+\underset{n=2}{\overset{p1}{}}f_n(k)x^n+x^ps(x),$$
(8)
where $`s(x)`$ is a singular function: $`s(0),`$ $`s^{(n)}(0)\mathrm{}`$, but $`xs(x)0`$ at $`x0`$. This result is of primary importance. It is the singular function that can lead to a deviation from the standard mean-field behavior. A method for determining $`s(x)`$ is proposed below. For convenience we use a power-law degree distribution $`P(k)k^\gamma `$. Then, $`k^4`$ diverges for $`\gamma 5`$, $`k^3`$ diverges for $`\gamma 4`$, and $`k^2`$ is divergent for $`\gamma 3`$.

We begin with the case $`f_{2n+1}=0`$. Let us consider the critical behavior for $`f_4>0`$ and different $`\gamma `$.

(a) $`\gamma >5`$. Here, $`k^4`$ converges. The coefficient $`f_4`$ is finite and $`\mathrm{\Phi }`$ has a usual form (4). At $`H=0`$ Eq. (6) leads to the standard critical behavior:

$$x\tau ^{1/2}\text{}\mathrm{\Delta }C0\text{}\chi _x(H=0)\tau ^1.\text{ }$$
(9)
Here, $`\tau =1T/T_c`$, $`\mathrm{\Delta }C1/k^4`$ is the jump of the specific heat at $`T_c`$. It tends to zero when $`\gamma 5`$.

(b) $`3<\gamma 5`$. Here, $`k^4`$ diverges. According to Eq. (8), a singular term $`x^4s(x)`$ is expected to appear in $`\mathrm{\Phi }`$. The function $`\varphi (x,y)`$ in Eq. (2) may be divided into two parts:

$$\varphi (x,y)=\underset{m=0}{\overset{\mathrm{}}{}}\underset{1l<\gamma 1}{}\varphi _{ml}x^my^l+G(x,y).$$
(10)
The first term determines a non-singular contribution to $`\mathrm{\Phi }`$ in Eq. (1). The second term, $`G(x,y)`$, gives the singular contribution. For $`4<\gamma 5`$ the function $`G(x,y)`$ has the following properties: $`G(0,y)y^4`$ at small $`y`$, $`d^nG(0,y)/dy^n0`$ for $`n4`$ at $`y=0`$. Owing to the condition (3), the function $`G(x,y)`$ increases slower than $`y^3`$ at $`y1`$.

It is convenient to use a continuous approximation for a degree distribution $`P(k)=Ak^\gamma `$, where $`A`$ is a normalization factor. Then Eq. (1) may be rewritten in the following form:

$$\mathrm{\Phi }(x,0)=f_2x^2+\phi (x)x^4+A_m^{\mathrm{}}\frac{dk}{k^\gamma }G(x,xk),$$
(11)
where $`m`$ is the smallest degree in $`P(k)`$, $`\phi (x)`$ is a smooth function which is determined by the convergent moments $`k^l`$ with $`l3`$, $`\phi (0)0`$. In order to find $`\mathrm{\Phi }`$ at small $`x`$, one can put $`G(x,y)G(0,y)`$ into Eq. (11). Considering the integral over a variable $`y=xk`$, one can show that the region $`mxyb`$ gives a leading contribution. $`b`$ is a model parameter. We obtain

$`\mathrm{\Phi }(x,H)`$ $`=`$ $`Hx+a(TT_c)x^2+Ax^4\mathrm{ln}(b/x),`$ (12)
$`\mathrm{\Phi }(x,H)`$ $`=`$ $`Hx+a(TT_c)x^2+Bx^{\gamma 1},`$ (13)
for $`\gamma =5`$ and $`3<\gamma <5`$, respectively. $`A`$ and $`B`$ are constants at $`T=T_c`$. The critical behavior for the case $`f_3=0`$, $`f_4>0`$ at $`\gamma >3`$ is summarized in Table I. It was observed in the Ising model on ‘equilibrium’ uncorrelated random networks .

Note that the divergence of $`k^4`$ does not change the critical behavior of $`\chi _x`$.

(c) $`2<\gamma 3`$. Now $`k^2`$ diverges. Equation (5) shows that in this situation there is no phase transition at any finite temperature. According to Eq. (8), a singular term $`x^2s(x)`$ is expected to appear in $`\mathrm{\Phi }`$. Calculations of $`\mathrm{\Phi }`$ may be carried out in a similar way as above, using the function $`G(x,y)=\varphi (x,y)_{m=1}\varphi _{m1}x^my`$, Eq. (10). We find

$`\mathrm{\Phi }(x,H)`$ $`=`$ $`Hx+Cx^2Dx^2\mathrm{ln}(b/x),`$ (14)
$`\mathrm{\Phi }(x,H)`$ $`=`$ $`Hx+C^{}x^2D^{}x^{\gamma 1}`$ (15)
for $`\gamma =3`$ and $`2<\gamma <3`$, respectively. The coefficients $`C,`$ $`D`$ , $`C^{}`$ and $`D^{}`$ are functions of $`T`$ and $`H`$. If $`C,`$ $`D`$ , $`C^{}`$ and $`D^{}`$ are positive, then at arbitrary finite temperature $`T1`$ in zero field $`H=0`$, Eq. (6) has stable nontrivial solutions:

$`x`$ $`=`$ $`b\mathrm{exp}((2C+D)/(2D)),\text{ }\chi _x=1/(2D),`$ (16)
$`x`$ $`=`$ $`\left({\displaystyle \frac{2C^{}}{(\gamma 1)D^{}}}\right)^{1/(3\gamma )},\text{ }\chi _x={\displaystyle \frac{1}{2(3\gamma )C^{}}}`$ (17)
for $`\gamma =3`$ and $`2<\gamma <3`$, respectively. As in the Landau theory, the results (16) and (17) are obtained in terms of the coefficients of the thermodynamic potential. Note that in the situation where a phase transition is absent at any finite temperature, the phenomenological theory can not determine the temperature behavior of the coefficients $`C,`$ $`D`$, $`C^{}`$ and $`D^{}`$. In this situation, the temperature dependences of the coefficients can be found only by a microscopic theory (see below).

In the general case the symmetry of the model admits non-zero $`f_3`$. In the case $`f_3>0`$ the analysis of the analytical properties of $`\mathrm{\Phi }`$ can be performed as above. For $`\gamma >4`$, $`\mathrm{\Phi }`$ is a smooth function of $`x`$ and leads to the standard critical behavior. At $`\gamma 4`$, $`\mathrm{\Phi }`$ contains a singular term: (a) $`x^3\mathrm{ln}x`$ at $`\gamma =4;`$ (b) $`x^{\gamma 1}`$ in the ranges $`4>\gamma >3`$ and $`3>\gamma >2`$; (c) $`x^2\mathrm{ln}x`$ at $`\gamma =3.`$ The corresponding critical behavior for $`\gamma >3`$ is represented in Table I. This behavior was observed for percolation on uncorrelated random networks . In the range $`2<\gamma 3`$ when $`k^2`$ diverges, $`\mathrm{\Phi }`$ has the universal form, Eqs. (14) and (15).

In the case $`f_3<0`$, or $`f_4<0`$ if $`f_3=0`$, when $`k^2`$ converges, a first-order phase transition occurs. In agreement with this prediction we found such a transition in the $`q`$state Potts model with $`q3`$ on ‘equilibrium’ uncorrelated random networks by use of the approach of Ref. . In the limit $`\gamma 3`$ the jump of the order parameter at the transition tends to zero, and the transition transforms into the infinite order phase transition. A detailed study of the transition will be given elsewhere.

In order to complete Table I, let us discuss the temperature behavior in the case $`2<\gamma 3`$ within the microscopic theory of the Ising model and percolation on ‘equilibrium’ uncorrelated random networks . For this purpose we use the more general ferromagnetic $`q`$state Potts model which at $`q=1`$ and 2 is equivalent to percolation and the Ising model, respectively (see, for example Ref. ). Using the approach of Ref.  we obtain that in the Potts model a continuous phase transition occurs at the exact critical temperature

$$T_c=2/\mathrm{ln}\frac{k^2+(q2)k}{k^22k}.$$
(18)
Hereafter, we set the energy of ferromagnetic interaction between nearest neighbors $`J=1`$. The parameter $`p_c=1\mathrm{exp}(2/T_c)`$ determines the percolation threshold  and establishes the relation between the Potts model and percolation.

Let us consider the case $`2<\gamma 3.`$ Here, $`k^2`$ diverges. $`T_c`$ is infinite for the infinite networks. In any finite network, $`k^2<\mathrm{}`$, and $`T_c`$ is finite, although it may be very high, $`T_c2k^2/(kq)`$. In the temperature region $`T1`$, where $`x1`$, the Potts model demonstrates universal behavior at all $`q1`$:

$`x`$ $``$ $`(q/k)e^{qT/k},\text{ }\chi _xT^1,`$ (19)
$`x`$ $``$ $`T^{1/(3\gamma )},\text{ }\chi _xT^2`$ (20)
for $`\gamma =3`$ and $`2<\gamma <3`$, respectively. Without the continuum approximation, we have instead of $`k`$ in the exponential, a constant which is determined by the complete form of $`P(k)`$. In Ehrenfest’s classification, this transition is of infinite order, as all temperature derivatives of $`\mathrm{\Phi }`$ are finite. The results (16) and (17) agree with Eqs. (19) and (20) if we put $`C,C^{}T^2`$, $`D,D^{}T`$. The exponential behavior (19) has been revealed in epidemic spreading within scale-free networks with $`\gamma =3`$  and in percolation on these networks .

At small $`x`$ there is the following relationship between the response function $`\chi _x`$ and the susceptibility $`\chi =dM/dH`$: $`\chi 2/qT+k\chi _x/q`$. At $`\gamma >3`$ the critical behavior of $`\chi `$ is determined by $`\chi _x`$: $`\chi k\chi _x/q`$. At $`2<\gamma 3`$ we have $`\chi 1/qT,`$ as the paramagnetic contribution $`2/qT`$ is of the order of $`\chi _x`$ or larger.

Let us discuss the results of theoretical studies and simulations of critical phenomena in different networks on the basis of the phenomenological theory.

The phenomenological theory as well as the Landau theory assumes that the contribution of fluctuations to the thermodynamic potential is small. Above we have shown that the theory gives the exact critical behavior of the Ising model, percolation and Potts model on ‘equilibrium’ uncorrelated random networks . The reason for this is that these networks have a local tree-like structure and vertices are statistically equivalent. Due to these properties, vertices in the networks can be regarded as forming a Bethe-lattice structure for which a mean-field approach is exact . It means that the fluctuation contribution is negligibly small. Note that in a graph with a Cayley tree-like structure, vertices are statistically inequivalent, and a mean-field approach is valid only for vertices deep within such a graph .

The origin of a deviation from the standard mean-field critical behavior is different for regular lattices and networks. In a regular lattice the deviation is caused by strong critical fluctuations which depend crucially on the space dimension. In networks the deviation is brought about by the most connected vertices which induce strong correlations in their close neighborhood. With decreasing exponent $`\gamma `$ the relative number of highly connected vertices increases and their role turns out to be more important.

Recent simulations  of the Ising model on the Barabási-Albert scale-free network with the degree distribution $`P(k)=Ak^3`$ revealed a temperature behavior described by Eq. (19) as well as in the Potts model on ‘equilibrium’ uncorrelated random networks. Unlike the latter networks, the Barabási-Albert net is correlated. Nevertheless, results of the simulations agree with the universal temperature behavior predicted by the phenomenological theory for networks with the divergent moment $`k^2`$.

The exact results for percolation on small-world networks were obtained in Ref. . The Ising and $`XY`$ models on small-world networks were studied analytically and by use of Monte Carlo simulations. These networks, introduced by Watts and Strogatz , have a Poisson-like degree distribution and large clustering coefficients. It was found that the critical behavior is characterized by the standard mean-field critical exponents. This result is in agreement with the prediction of the phenomenological theory that the standard mean-field critical behavior should occur in networks described by a degree distribution with convergent moments.

Thus, the analysis of critical behavior of the Ising, Potts and $`XY`$ models, percolation and epidemic spreading on uncorrelated random, scale-free and small-word networks shows that the critical behavior in networks depends crucially on the form of the degree distribution and the symmetry underlying a model in a complete agreement with the phenomenological theory. The studied networks differ in clustering coefficients, degree correlations, etc. However, the situation looks like these characteristics are not relevant for critical behavior. Their role in the formation of the critical exponents is still unclear. Further investigations in this topic are necessary. It would be interesting to find a network where critical behavior differs from a mean-field one due to strong fluctuations. However, even in this case it is expected that when $`T`$ tends to $`T_c`$, a temperature region of the mean-field critical behavior precedes a region of strong fluctuations.

Real systems like the Internet, WWW or biological networks have a network structure with a power-law degree distribution with exponent $`\gamma `$ below 3, see Refs. . The phenomenological theory predicts for this case the power-law critical behavior (20), see also the last line in Table I. This behavior agrees with an empirical observation  for the nd.edu domain of the WWW, where the variations of the size of the giant component under random damage were studied. This size and the number of damaged vertices play a role of the order parameter $`x`$ and the control parameter, respectively.

In conclusion, the phenomenological theory shows that the deviation of the critical behavior of interacting systems with a network structure from the standard mean-field behavior emerges when a degree distribution has divergent moments. The theory predicts different classes of critical behavior in agreement with microscopic studies of the Ising and Potts models, and percolation on ‘equilibrium’ uncorrelated random networks. It also agrees with results previously obtained for various models on small-world and evolving networks. The theory can easily be generalized for models with a multicomponent order parameter and can give useful insight into collective effects in different networks discussed in connection with the Internet, biological networks, etc. Using this approach, one can also study the critical relaxation in models on networks.

S.N.D and J.F.F.M. were partially supported by the project POCTI/99/FIS/33141. A.G. acknowledges the support of the NATO program OUTREACH. We also thank A.N. Samukhin for many useful discussions.

E-mail address: goltsev@pop.ioffe.rssi.ru
E-mail address: sdorogov@fc.up.pt
E-mail address: jfmendes@fc.up.pt