text:
# RENORMALIZATION GROUP IN THE STATISTICAL THEORY OF TURBULENCE: TWO-LOOP APPROXIMATION

## I Introduction

One of the oldest open problems in theoretical physics is that of describing fully developed turbulence on the basis of a microscopic model. The latter is usually taken to be the stochastic Navier–Stokes (NS) equation subject to an external random force which mimics the energy input by the large-scale modes; see, e.g., . The aim of the theory is to verify the basic principles of the celebrated Kolmogorov–Obukhov phenomenological theory, study deviations from this theory, determine the dependence of various correlation functions on the times, distances, external (integral) and internal (viscous) turbulence scales, and derive the corresponding scaling dimensions. Most results of this kind were obtained within the framework of numerous semiphenomenological models which cannot be considered to be the basis for construction of a regular expansion in certain small (at least formal) parameter .

An important exception is provided by the renormalization group (RG) method that was earlier successfully applied in the theory of critical behaviour to explain the origin of critical scaling and to calculate universal quantities (critical dimensions and scaling functions) in the form of the $`\epsilon `$ expansions .

The RG approach to the stochastic NS equation, pioneered in \[4–7\], allows one to prove the existence of the infrared (IR) scale invariance with exactly known “Kolmogorov” dimensions and the independence of the correlation functions of the viscous scale (the second Kolmogorov hypothesis), and calculate a number of representative constants in a reasonable agreement with experiment. Detailed review of the RG theory of turbulence and more references can be found in .

In contrast to the standard $`\varphi ^4`$ model of critical behaviour , where the critical exponents are known up to the order $`\epsilon ^5`$ (five-loop approximation), all the calculations in the RG approach to the stochastic NS equation have been confined with the simplest one-loop approximation. The reason for this distinction is twofold. First, the multiloop calculations for this dynamical model are rather involved: one can say that the two-loop calculation for the stochastic NS equation is as cumbersome as the four-loop calculation for the conventional $`\varphi ^4`$ model. Second, the critical dimensions for the most important physical quantities (velocity and its powers, frequency, energy dissipation rate and so on) are given by the one-loop approximation exactly (the corresponding $`\epsilon `$ series terminate at first-order terms) and the higher-order calculations for them are not needed.

However, the $`\epsilon `$ series for other important quantities do not terminate and the calculation of the higher-order terms for them is of great interest. In this paper, we present the results of the two-loop calculation for a number of such quantities: the ultraviolet (UV) correction exponent $`\omega `$, the Kolmogorov constant $`C_K`$ and the inertial-range skewness factor $`𝒮`$.

## II The model, field theoretic formulation and renormalization

Detailed exposition of the RG theory of turbulence and the bibliography can be found in ; below we restrict ourselves to only the necessary information.

As the microscopic model of the fully developed, homogeneous, isotropic turbulence of an incompressible viscous fluid one usually takes the stochastic NS equation with a random driving force

$$_t\phi _i=\nu _0^2\phi _i_iP+F_i,_t_t+(\phi ).$$
(1)
Here $`\phi _i`$ is the transverse (due to the incompressibility) three-dimensional vector velocity field, $`P`$ and $`F_i`$ are the pressure and the transverse random force per unit mass (all these quantities depend on $`x\{t,𝐱\}`$), $`\nu _0`$ is the kinematical viscosity coefficient, $`^2`$ is the Laplace operator and $`_t`$ is the Lagrangian derivative. The problem (1) is studied on the entire $`t`$ axis and is augmented by the retardation condition and the condition that $`\phi _i`$ vanishes for $`t\mathrm{}`$. We assume for $`F`$ a Gaussian distribution with zero mean and correlator

$$F_i(x)F_j(x^{})=\delta (tt^{})(2\pi )^3𝑑𝐤P_{ij}(𝐤)d_F(k)\mathrm{exp}\left[\mathrm{i}𝐤\left(𝐱𝐱^{}\right)\right],$$
(2)
where $`P_{ij}(𝐤)=\delta _{ij}k_ik_j/k^2`$ is the transverse projector and $`d_F(k)`$ is some function of $`k|𝐤|`$ and model parameters.

The stochastic problem (1), (2) is equivalent to the field theoretic model of the doubled set of fields $`\mathrm{\Phi }\{\phi ,\phi ^{}\}`$ with action functional

$$S(\mathrm{\Phi })=\phi ^{}D_F\phi ^{}/2+\phi ^{}[_t\phi +\nu _0^2\phi (\phi )\phi ],$$
(3)
where $`D_F`$ is the random force correlator (2) and the required integrations over $`x=\{t,𝐱\}`$ and summations over the vector indices are understood.

The standard RG technique can be applied to the model (3) if the function $`d_F(k)`$ is chosen in the form

$$d_F(k)=D_0k^{12\epsilon }h(m/k),h(0)=1.$$
(4)
Here $`D_0`$ is the positive amplitude factor, $`m=1/L`$ is the reciprocal of the integral turbulence scale $`L`$, the function $`h(m/k)`$ provides the IR regularization and the exponent $`\epsilon >0`$ plays the part of the RG expansion parameter, similar to that played by $`\epsilon =4d`$ in Wilson’s theory of critical phenomena . The real (physical) value of this parameter is $`\epsilon =2`$: idealized energy injection by infinitely large eddies corresponds to $`d_F(k)\delta (𝐤)`$, and the function (4) for $`\epsilon 2`$ and the appropriate choice of the amplitude can be considered as a power-law model of the three-dimensional $`\delta `$ function.

The model (3) is logarithmic (the coupling constant $`g_0D_0/\nu _0^3`$ is dimensionless) at $`\epsilon =0`$, and the UV divergences have the form of the poles in $`\epsilon `$ in the correlation functions of the fields $`\phi `$ and $`\phi ^{}`$. Superficial UV divergences, whose removal requires counterterms, are present only in the 1-irreducible function $`\phi ^{}\phi `$, and the corresponding counterterm has the form $`\phi ^{}^2\phi `$. Thus for the complete elimination of the UV divergences it is sufficient to perform the multiplicative renormalization of the parameters $`\nu _0`$ and $`g_0=D_0/\nu _0^3`$ with the only independent renormalization constant $`Z_\nu `$:

$$\nu _0=\nu Z_\nu ,g_0=g\mu ^{2\epsilon }Z_g,Z_g=Z_\nu ^3(D_0=g_0\nu _0^3=g\mu ^{2\epsilon }\nu ^3).$$
(5)
Here $`\mu `$ is the reference mass in the minimal subtraction (MS) scheme, which we always use in what follows, $`g`$ and $`\nu `$ are renormalized analogues of the bare parameters $`g_0`$ and $`\nu _0`$, and $`Z=Z(g,\epsilon ,d)`$ are the renormalization constants. In the MS scheme they have the form “1 + only poles in $`\epsilon `$,” in particular,

$`Z_\nu =1+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}a_k(g)\epsilon ^k=1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}g^n{\displaystyle \underset{k=1}{\overset{n}{}}}a_{nk}\epsilon ^k,`$ (6)
with the one-loop coefficient $`a_{11}=1/40\pi ^2`$ .

## III Two-loop approximation for the RG functions, fixed point and the UV correction exponent

We have performed calculation of the constant $`Z_\nu `$ with the accuracy of $`O(g^2)`$ (two-loop approximation). The calculation is rather involved and will be presented elsewhere (some details can be found in ), and below we give only the results for the residues $`a_{22}`$ and $`a_{21}`$ at the second-order and first-order poles in $`\epsilon `$ in the representation (6):

$$a_{22}/a_{11}^2=1,a_{21}/a_{11}^21.65.$$
(7)
The knowledge of the renormalization constant $`Z_\nu `$ to order $`O(g^2)`$ allows for the calculation of the RG functions, the anomalous dimension $`\gamma _\nu `$ and the beta function $`\beta _g(\epsilon ,g)`$, with the following accuracy:

$`\beta (g,\epsilon )\stackrel{~}{𝒟}_\mu g=g\left(2\epsilon +3\gamma _\nu (g)\right),`$

$$\gamma _\nu (g)\stackrel{~}{𝒟}_\mu \mathrm{ln}Z_\nu =2g_ga_1(g)=2\left(a_{11}g+2a_{21}g^2\right)+O(g^3),$$
(8)
where $`\stackrel{~}{𝒟}_\mu `$ is the operation $`\mu /\mu `$ at fixed bare parameters. In the MS scheme only the residues at the first-order poles in $`\epsilon `$, that is, only the coefficients $`a_{k1}`$, contribute to the RG functions owing to the UV finiteness of the latter.

The coordinate of the fixed point is determined by the condition that $`\beta (g_{})=0`$. From (7) and (8) we thus obtain:

$$g_{}=(40\pi ^2\epsilon /3)(1+\lambda \epsilon )+O(\epsilon ^3),\lambda 2a_{21}/3a_{11}^21.10.$$
(9)
The correction exponent $`\omega `$ is determined by the slope of the beta function at the fixed point, $`\omega =\beta ^{}(g_{})`$. Thus from (8) and (9) we obtain the first and second terms of its $`\epsilon `$ expansion:

$$\omega =2\epsilon (1\lambda \epsilon )+O(\epsilon ^3).$$
(10)
## IV Two-loop calculation of the Kolmogorov constant

The Kolmogorov constant $`C_K`$ can be defined as the (dimensionless) coefficient in the inertial-range expression $`S_2(r)=C_K(\overline{ϵ}r)^{2/3}`$ for the second-order structure function, predicted by the Kolmogorov–Obukhov theory and confirmed by experiment . Here $`\overline{ϵ}`$ is the mean energy dissipation rate and the $`n`$-th order (longitudinal, equal-time) structure function is defined as

$$S_n(r)[\phi _r(t,𝐱+𝐫)\phi _r(t,𝐱)]^n,\phi _r(\phi _ir_i)/r,r|𝐫|.$$
(11)
Using the exact relation $`S_3(r)=4\overline{ϵ}r/5`$ that follows from the energy balance equation , the constant $`C_K`$ can be related to the inertial-range skewness factor:

$`𝒮S_3/S_2^{3/2}=(4/5)C_K^{3/2}.`$

Many studies have been devoted to the derivation of $`C_K`$ within the framework of the RG approach; see Refs. \[11–18\]. In order to obtain $`C_K`$, it is necessary to augment the solution of the RG equation for $`S_2`$ by some formula that relates the amplitude $`D_0`$ in the random force correlator (4) to the physical parameter $`\overline{ϵ}`$. In particular, in  the first-order term of the $`\epsilon `$ expansion for the pair correlator was combined with the so-called eddy damped quasi-Markovian approximation for the energy transfer function, taken at $`\epsilon =2`$. More elementary derivation, based on the exact relation between $`\overline{ϵ}`$ and the function $`d_F(k)`$ from (4) was given in ; see also . In spite of the reasonable agreement with the experiment, such derivations are not immaculate from the theoretical viewpoints. Their common flaw is that the relation between $`\overline{ϵ}`$ and $`D_0`$ is unambiguous only in the limit $`\epsilon 2`$:

$$\underset{\epsilon 2}{lim}\frac{D_0}{4\pi ^2(2\epsilon )}=\overline{ϵ},$$
(12)
so that the coefficients of the corresponding $`\epsilon `$ expansions appear in fact arbitrary; see the discussion in  and Sec. 2.10 of . The ambiguity is a consequence of the fact that the notion of the Kolmogorov constant has no definite extension to the nonphysical range $`\epsilon <2`$.

The experience on the RG theory of critical behaviour shows that unambiguous $`\epsilon `$ expansions can be written for universal quantities, such as critical exponents, normalized scaling functions and ratios of amplitudes in scaling laws . The constant $`C_K`$ extended to the range $`\epsilon <2`$ as in  involves a bare parameter, $`D_0`$, and hence is not universal.

To circumvent this difficulty, we propose below an alternative derivation that relates $`C_K`$ to an universal quantity and thus does not involve any relation between $`D_0`$ and $`\overline{ϵ}`$, and calculate $`C_K`$ to second order of the expansion in $`\epsilon `$ (previous attempts have been confined with the first order). Consider the ratio

$$Q(\epsilon )r_rS_2(r)/|S_3(r)|^{2/3}=r_rS_2(r)/(S_3(r))^{2/3}.$$
(13)
The operation $`r_rr/r`$ kills the constant contribution $`\phi ^2`$ in $`S_2`$ that diverges as $`\mathrm{\Lambda }\mathrm{}`$ for $`\epsilon <3/2`$; in $`S_3`$ such constant contributions are absent.

Solving the RG equations for the quantities in $`Q(\epsilon )`$ in the IR range ($`\mathrm{\Lambda }r1`$) for general $`0<\epsilon 2`$ gives

$$S_3(r)=D_0r^{3\mathrm{\Delta }_\phi }f_3(\epsilon ),r_rS_2(r)=D_0^{2/3}r^{2\mathrm{\Delta }_\phi /3}f_2(\epsilon ),\mathrm{\Delta }_\phi =12\epsilon /3;$$
(14)
see e.g. . Thus the quantity $`Q(\epsilon )=f_2/(f_3)^{2/3}`$ in (13) does not depend on $`D_0`$ and can be calculated in the form of a regular $`\epsilon `$ expansion. We calculated the scaling functions $`f_{2,3}`$ to the second order of the $`\epsilon `$ expansion, which corresponds to the two-loop approximation in (9), (10), and obtained:

$$Q(\epsilon )=(1/3)(20\epsilon )^{1/3}\left[1+0.525\epsilon +O(\epsilon ^2)\right].$$
(15)
It is worth noting that the $`\epsilon `$ expansion for $`f_3`$ can be obtained not only from the direct perturbative calculation, but also from the exact expression

$$S_3(r)=\frac{6\mathrm{\Gamma }(2\epsilon )}{2^{2\epsilon }\pi ^{3/2}\mathrm{\Gamma }(3/2+\epsilon )}D_0r^{3\mathrm{\Delta }_\phi }$$
(16)
that follows from the energy balance equation and in the limit $`\epsilon 2`$, along with the formula (12), reproduces the correct coefficient $`4/5`$ (see above).

The value of $`Q(\epsilon )`$ at $`\epsilon =2`$ determines the Kolmogorov constant and skewness factor through the exact relations $`C_K=610^{2/3}Q(2)`$, $`𝒮=[1.5Q(2)]^{3/2}`$, which follow from the definitions and the identity $`r_rr^\lambda =\lambda r^\lambda `$ for any $`\lambda `$. Substituting the value of (15) we obtain $`C_K=3.02`$, $`𝒮=0.15`$. If we retained only the first-order term in (15) we would have obtained $`C_K=1.47`$ and $`𝒮=0.45`$. We also recall the experimental estimates recommended in : $`C_K1.9`$ and $`𝒮0.28`$.

## V Conclusion

We have accomplished the complete two-loop calculation of the renormalization constant and RG functions for the stochastic problem (1)–(4) and derived the coordinate of the fixed point, the UV correction exponent $`\omega `$, the Kolmogorov constant $`C_K`$ and the inertial-range skewness factor $`𝒮`$ to the second order of the corresponding $`\epsilon `$ expansions. The new point is not only the inclusion of the second-order correction, but also the derivation of $`C_K`$ through an universal (in the sense of the theory of critical behaviour) quantity.

Of course, one should have not expected that the second-order terms of the $`\epsilon `$ expansions would be small in comparison to the first-order terms. The experience on the RG theory of critical behaviour shows that such corrections are not small for dynamical models (in contrast to static ones) and for amplitudes (in contrast to exponents); see . It is thus rather surprising that in our case the account of the two-loop contributions leads to reasonable changes in the results.

Although the $`\epsilon ^2`$ correction to $`\omega `$ in (10) is rather large, it does not change its sign and hence does not destroy the IR stability of the fixed point.

The first-order approximation $`C_K=1.47`$ underestimates, and the second-order approximation $`C_K=3.02`$ overestimates the conventional experimental value of the Kolmogorov constant $`C_K1.9`$ . Thus the experimental value of $`C_K`$ (and hence for $`𝒮`$) lies in between of the two consecutive approximations. A similar situation is encountered for the well-known Heisenberg model , where the analogue of the Kolmogorov constant is known exactly and lies between the first-order and second-order approximations given by the corresponding $`\epsilon `$ expansion . If we assume, by the analogy with the Heisenberg model, that the (unknown) exact predictions for $`C_K`$ and $`𝒮`$ lie between the first two approximations, we may conclude that our calculation has given a very good estimate for these quantities.

The work was supported by the Nordic Grant for Network Cooperation with the Baltic Countries and Northwest Russia No. FIN-18/2001 and the GRACENAS Grant No. E00-3-24.


format: text

text:
# 1 Introduction

## 1 Introduction

The investigation of deep inelastic scattering ($`DIS`$) of leptons off the nucleon is an important tool to get fundamental information on the structure of quark distributions in the nucleon. At large values of the Bjorken variable $`x`$ the ratio $`d/u`$, or equivalently the ratio $`F_2^n/F_2^p`$, is known to have a high degree of theoretical significance . However, our present knowledge of the large-$`x`$ $`n/p`$ ratio is quite poor, mainly because its extraction from inclusive deuteron $`DIS`$ data is inherently model dependent. This fact has led to the suggestion of new strategies. One of them is to try to exploit the mirror symmetry of $`A=3`$ nuclei; in other words, thanks to nuclear charge symmetry, one expects that the magnitude of the $`EMC`$ effect

$`_{EMC}^A{\displaystyle \frac{F_2^A}{F_2^D}}{\displaystyle \frac{F_2^p+F_2^n}{ZF_2^p+NF_2^n}}`$ (1)
is very similar in $`{}_{}{}^{3}He`$ and $`{}_{}{}^{3}H`$ and hence the so called super-ratio

$`𝒮_{EMC}{\displaystyle \frac{_{EMC}^{{}_{}{}^{3}He}}{_{EMC}^{{}_{}{}^{3}H}}}={\displaystyle \frac{F_2^{{}_{}{}^{3}He}}{F_2^{{}_{}{}^{3}H}}}{\displaystyle \frac{2F_2^n+F_2^p}{2F_2^p+F_2^n}},`$ (2)
should be very close to unity regardless of the size of the $`EMC`$ ratios . In this case the $`n/p`$ ratio could be extracted directly from the ratio $`F_2^{{}_{}{}^{3}He}/F_2^{{}_{}{}^{3}H}`$ without significant nuclear modifications. However, we observe that, even if nuclear charge symmetry were exact, the motion of protons and neutrons in a non isosinglet nucleus (say $`{}_{}{}^{3}He`$) is somewhat different due to the spin-flavor dependence of the nuclear force.

We have therefore performed  a detailed study of inclusive $`DIS`$ from mirror $`A=3`$ nuclei at large $`x`$ with the aim of estimating the theoretical uncertainties on the extraction of the $`n/p`$ ratio from such measurements. To this end we have considered a variety of $`EMC`$ models both with and without modifications of the nucleon structure function in the medium.

## 2 $`EMC`$ models with no modification of bound nucleon

We have explored  in greater details the $`EMC`$ model used in , which is based on the Virtual Nucleon Convolution ($`VNC`$) approach with no modification of the nucleon structure function in the medium.

We have analyzed the effects of: i) charge symmetry breaking terms in the nucleon-nucleon ($`NN`$) interaction; ii) finite $`Q^2`$ effects in the impulse approximation; iii) the role of different prescriptions for the nucleon Spectral Function normalization providing baryon number conservation; and iv) the role of different parton distribution function sets. Additionally we have also compared the predictions of the $`VNC`$ model with the ones obtained within the light-cone ($`LC`$) formalism.

Within the $`VNC`$ model, in which no modification of the bound nucleon is considered, the deviation of the super-ratio (2) from unity is found to stay within $`1\%`$ only for $`x\text{ }<0.75`$ in agreement with  in which the effects of (i)-(iv) have not been considered.

## 3 $`EMC`$ models with $`F_2^N`$ modifications in the medium

However the above estimations are not complete, since $`VNC`$ and $`LC`$ models are just two of the many models of the $`EMC`$ effect and, moreover, they underestimate significantly the $`EMC`$ data at large $`x`$ for a variety of nuclei. In particular, the convolution approach is not able to reproduce the minimum of the $`EMC`$ ratio around $`x0.7`$ as well as the subsequent sharp rise at larger $`x`$ (see ). Thus in order to draw final conclusions about the size of the deviation of the super-ratio $`𝒮_{EMC}`$ from unity one should investigate effects beyond those predicted by the convolution approach.

Furthermore, it is very important to asses any isospin dependence of the $`EMC`$ effect in order to extract $`F_2^n`$ from $`{}_{}{}^{3}He`$ and $`{}_{}{}^{3}H`$ data in a reliable way. An isospin dependence of the $`EMC`$ effect is naturally expected from the differences in the relative motion of $`pn`$ and $`nn`$ ($`pp`$) pairs in $`{}_{}{}^{3}H`$ ($`{}_{}{}^{3}He`$). Since the interaction of a $`pn`$ pair is more attractive than the one of a $`nn`$ pair, the overlapping probability to find a $`NN`$ pair with $`r_{NN}1fm`$, may be $`40\%`$ larger for a $`pn`$ pair than for a $`nn`$ pair (see ). This is a very important isospin effect in mirror $`A=3`$ nuclei, which can lead to deviations of the super-ratio (2) from unity depending on the size of the $`EMC`$ effect itself. It is worth noting that the nuclear charge symmetry will not limit such deviations.

Hence we have carried out a detailed analysis of the super-ratio within a broad range of models of the $`EMC`$ effect, which take into account possible modifications of the bound nucleons in nuclei, like: i) a change in the quark confinement size (including swelling); ii) the possible presence of clusters of six quarks; and iii) the suppression of point-like configurations due to color screening. Our main result  is that one cannot exclude the possibility that the cancellation of the nuclear effects in the super-ratio may occur only at a level of $`3\%`$, resulting in a significant uncertainty (up to $`12\%`$ for $`x0.7÷0.8`$) in the extraction of the free $`n/p`$ ratio from the ratio of the measurements of the $`{}_{}{}^{3}He`$ and $`{}_{}{}^{3}H`$ $`DIS`$ structure functions (see Fig. 1).

In  it was suggested that, once the ratio $`F_2^{{}_{}{}^{3}He}/F_2^{{}_{}{}^{3}H}`$ is measured, one can employ an iterative procedure to extract the $`n/p`$ ratio which can almost eliminate the effects of the dependence of the super-ratio $`𝒮_{EMC}`$ on the large-$`x`$ behavior of the specific structure function input. Namely, after extracting the $`n/p`$ ratio assuming a particular calculation of $`𝒮_{EMC}`$, one can use the extracted $`F_2^n`$ to get a new estimate of $`𝒮_{EMC}`$, which can then be employed for a further extraction of the $`n/p`$ ratio. Such a procedure can be iterated until convergence is achieved and self-consistent solutions for the extracted $`F_2^n/F_2^p`$ and the super-ratio $`𝒮_{EMC}`$ are obtained. In  a good convergence was achieved for $`x`$ up to $`0.8`$. However this result depends on the assumed validity of the convolution approach at large $`x`$. Within the $`EMC`$ models considered above we have found  that: i) the consistency between the $`n/p`$ ratio used as an input and the extracted one is not guaranteed; ii) the iteration diverges already at values of $`x`$ ($`0.7`$) smaller than the ones obtained in ; moreover, below $`x0.7`$ the iteration procedure converges to a value of the $`n/p`$ ratio which differs from the input one exactly by the amount of the $`EMC`$ effect implemented in the calculation of $`F_2^{{}_{}{}^{3}He}/F_2^{{}_{}{}^{3}H}`$.

Thus we conclude that in spite of the presented limitations the measurements with mirror $`A=3`$ nuclei will significantly improve our knowledge of the $`n/p`$ ratio at large $`x`$. However, it will be very important to complement them with the measurements of semi-inclusive processes off the deuteron, in which the momentum of the struck nucleon is tagged by detecting the recoiling one. Imposing the kinematical conditions that the detected momentum is low ($`\text{ }<150MeV/c`$), which means that the nucleons in the deuteron are initially far apart , it is possible to minimize significantly the nuclear effects. Furthermore, all the unwanted nuclear effects can be isolated by using the same reaction for the extraction of the proton structure function by detecting slow recoiling neutrons and comparing the results with existing hydrogen data, as well as by performing tighter cuts on the momentum of the spectator proton and then extrapolating to the neutron pole.


text:
# Continuum in the Excitation Spectrum of the 𝑆=1 Compound CsNiCl3

## Acknowledgments

We would like to thank A.A. Nersesyan, P. Pieri and G. Sierra for useful discussions.

—————————————–


text:
# Classification of scale free networks

## Abstract

While the emergence of a power law degree distribution in complex networks is intriguing, the degree exponent is not universal. Here we show that the betweenness centrality displays a power-law distribution with an exponent $`\eta `$ which is robust and use it to classify the scale-free networks. We have observed two universality classes with $`\eta 2.2(1)`$ and $`2.0`$, respectively. Real world networks for the former are the protein interaction networks, the metabolic networks for eukaryotes and bacteria, and the co-authorship network, and those for the latter one are the Internet, the world-wide web, and the metabolic networks for archaea. Distinct features of the mass-distance relation, generic topology of geodesics and resilience under attack of the two classes are identified. Various model networks also belong to either of the two classes while their degree exponents are tunable.

Emergence of a power law in the degree distribution $`P_D(k)k^\gamma `$ in complex networks is an interesting self-organized phenomenon in complex systems strogatz ; review1 ; review2 . Here, the degree $`k`$ means the number of edges incident upon a given vertex. Such a network is called scale-free (SF) physica . Real world networks which are SF include the author collaboration network newman  in social systems, the protein interaction network (PIN) pin  and the metabolic network metabolic  in biological systems, and the Internet internet  and the world-wide web (WWW) www ; ha  in communication systems. The power-law behavior means that most vertices are sparsely connected, while a few vertices are intensively connected to many others and play an important role in functionality. While the emergence of such a SF behavior in degree distribution itself is surprising, the degree exponent $`\gamma `$ is not universal and depends on the detail of network structure. As listed in Table 1, numerical values of the exponent $`\gamma `$ for various systems are diverse but most of them are in the range $`2<\gamma 3`$. From the viewpoint of theoretical physics, it would be interesting to search a universal quantity associated with SF networks.

Recently a physical quantity called “load” was introduced as a candidate for the universal quantity in SF networks. It quantifies the load of a vertex in the transport of data packet along the shortest pathways in SF networks load . It was shown that the load distribution exhibits a power law, $`P_L(\mathrm{})\mathrm{}^\delta `$, and the exponent $`\delta `$ is robust as $`\delta 2.2`$ for diverse SF networks with various degree exponents in the range $`2<\gamma 3`$. Since the universal behavior of the load exponent was obtained empirically, fundamental questions such as how the load exponent is robust in association with network topology or the possibility of any other universal classes existing, have not been explored yet. In this paper, we address those issues in detail.

While the load is a dynamic quantity, it is closely related to a static quantity, the “betweenness centrality (BC)”, commonly used in sociology to quantify how much a given person is influential in a society betweenness . To be specific, BC is defined as follows. Let us consider the set of the shortest pathways, or geodesics, between a pair of vertices $`(i,j)`$ and denote their number by $`C(i,j)`$. Among them, the number of the shortest pathways running through a vertex $`k`$ is denoted by $`C_k(i,j)`$. The fraction $`g_k(i,j)=C_k(i,j)/C(i,j)`$ may be interpreted as the amount of the role played by the vertex $`k`$ in social relation between two persons $`i`$ and $`j`$. Then the BC of the vertex $`k`$ is defined as the accumulated sum of $`g_k(i,j)`$ over all ordered pairs for which a geodesic exists, i.e.,

$$g_k=\underset{ij}{}g_k(i,j)=\underset{ij}{}\frac{C_k(i,j)}{C(i,j)}.$$
(1)
Because of only slight difference between load and BC, both quantities behave very closely. In fact, the BC $`g_k`$ of each vertex is exactly the same as the load for tree graphs. In general, distributions of the two are indistinguishable within available resolutions. The BC distribution follows a power law,

$$P_B(g)g^\eta ,$$
(2)
where $`g`$ means BC and the exponent $`\eta `$ is the same as the load exponent $`\delta `$. Since the topological feature of a network can be grasped more transparently using BC, we deal with BC in this work.

Based on numerical measurements of the BC exponent for a variety of SF networks, we find that SF networks can be classified into only two classes, say, class I and II. For the class I, the BC exponent is $`\eta 2.2(1)`$ and for the class II, it is $`\eta 2.0(1)`$. We conjecture the BC exponent for the class II to be exactly $`\eta =2`$ since it can be derived analytically for simple models. We show that such different universal behaviors in the BC distribution originate from different generic topological features of networks. Moreover, we study a physical problem, the resilience of networks under an attack, showing different behaviors for each class, as a result of such difference in generic topologies. It is found that the networks of class II are much more vulnerable to the attack than those for class I.

To obtain our results, we use available data for real world networks or existing algorithms for model networks. Once a SF network is constructed, we select a pair of vertices $`(i,j)`$ on the network and identify the shortest pathways between them. Next, BC is measured on each vertex along the shortest pathways using the modified version of the breath-first search algorithm introduced by Newman newmanalgorithm . We measure the BC distribution and the exponent $`\eta `$ for a variety of networks both in real world and in silico.

Real World and Artificial Networks Investigated
The networks that we find to belong to the class I with $`\eta =2.2(1)`$ include: (i) The co-authorship network in the field of the neuroscience, published in the period 1991-1998 coauthor , where vertices represent scientists and they are connected if they wrote a paper together. (ii) The protein interaction network of the yeast S. cerevisiae compiled by Jeong et al. pin  (PIN1), where vertices represent proteins and the two proteins are connected if they interact111The network is composed of disconnected clusters of different sizes, viz., small isolated clusters as well as a giant cluster. For both (ii) and (xi), the degree distribution is likely to follow a power law but there needs an exponential cutoff to describe its tail behavior for finite system. However, it converges to a clean power law for (xi) as system size increases, but the converging rate is rather slow jeenu . Despite this abnormal behavior in the degree distribution for finite system, the BC distribution follows a pure power law with the exponent $`\eta 2.2(1)`$ in (ii) and (xi). . (iii) The core of protein interaction network of the yeast S. cerevisiae obtained by Ito et al.222In contrast to (ii), the degree distribution obeys a power law. (PIN2) ito . (iv) The metabolic networks for 5 species of eukaryotes and 32 species of bacteria in Ref.metabolic , where vertices represent substrates and they are connected if a reaction occurs between two substrates via enzymes. The reaction normally occurs in one direction, so that the network is directed. (v) The Barabási-Albert (BA) model ba  when the number of incident edges of an incoming vertex $`m2`$. (vi) The geometric growth model by Huberman and Adamic ha . (vii) The copying model copying , whose degree exponent is in the range of $`2<\gamma 3`$. (viii) The undirected or the directed static model load , whose degree exponent is in the range of $`2<\gamma 3`$ or $`2<(\gamma _{\mathrm{in}},\gamma _{\mathrm{out}})3`$, respectively. (ix) The accelerated growth model proposed by Dorogovtsev et al. accel . (x) The fitness model fitness  with a flat fitness distribution. (xi) The stochastic model for the protein interaction networks introduced by Solé et al. sole . All those networks (i)-(xi) exhibit a power-law behavior in the BC distribution with the exponent $`\eta 2.2(1)`$. Detailed properties of each network are listed in Table 1. The representative BC distributions for real world networks (i), (iii), and (iv) are shown in Fig.1a.

The networks that we find to belong to the class II with $`\eta =2.0`$ include: (xii) The Internet at the autonomous systems (AS) level as of October, 2001 internet2 . (xiii) The metabolic networks for 6 species of archaea in Ref.metabolic . (xiv) The WWW of nd.edu www . (xv) The BA model with $`m=1`$ ba . (xvi) The deterministic model by Jung et al. jung . In particular, the networks (xv) and (xvi) are of tree structure, where the edge BC distribution can be solved analytically. The detailed properties of each network are listed in Table 1. The BC distributions for real world networks (xii) and (xiv) are shown in Fig.1b. Since the BC exponents of each class are very close numerically, one may wonder if there exist really two different universal classes apart from error bar. To make this point clear, we plot the BC distributions for the BA model with $`m=1,2`$ and $`3`$ in Fig.2, obtained from large system size, $`N=3\times 10^5`$. We can see clearly different behaviors between the two BC distributions for the cases of $`m=1`$ (class II) and of $`m=2`$ and $`3`$ (class I).

Topology of the Shortest Pathways
To understand the generic topological features of the networks in each class, we particularly focus on the topology of the shortest pathways between two vertices separated by a distance $`d`$. Along the shortest pathways, we count the total number of vertices $`(d)`$ lying on these roads, averaged over all pairs of vertices separated by the same distance $`d`$. Adopting from the fractal theory, $`(d)`$ is called the “mass-distance” relation. We find that it behaves in different ways for each class; For the class I, $`(d)`$ behaves nonlinearly (Figs.3a-b), while for the class II, it is roughly linear (Fig.3c-d).

For the networks belonging to the class I such as the PIN2 (iii) and the metabolic network for eukaryotes (iv), $`(d)`$ exhibits a non-monotonic behavior (Fig.3a-b), viz., it exhibits a hump at $`d_h10`$ for (iii) or $`d_h14`$ for (iv). To understand why such a hump arises, we visualize the topology of the shortest pathways between a pair of vertices, taken from the metabolic network of a eukaryote organism, Emericella nidulans (EN), as a prototypical example for the class I. Fig.4a shows such a graph with linear size 26 edges ($`d=26`$), where an edge between a substrate and an enzyme is taken as the unit of length. From Fig.4a, one can see that there exists a blob structure inside which vertices are multiply connected, while vertices outside are singly connected. What is characteristic for the class I is that the blob is localized in a small region. To see this, we measure the mass density $`m(r;d)`$, the average number of substrates or enzymes located at position $`r`$ $`\left(_{r=1}^dm(r;d)=(d)\right)`$. The average is taken over all possible pairs of vertices ($`56`$ pairs), separated by the same distance $`d=26`$. Note that the metabolic network is directed, so that the position $`r`$ is uniquely defined. As shown in Fig.5, we find that $`m(r;d)`$ is sharply peaked at $`r=3`$, and is larger than 1 only at $`r=2,4,`$ and 6 for substrates. Thus the blob structure is present, even after taking averages, and is localized in a small region of size $`d_b4`$$``$$`5`$, centered at almost the same position $`r3`$$``$$`4`$ for different pairs of vertices. The blob size $`d_b`$ can be measured in another way. In a given graph of the shortest pathways, we delete singly connected substrates successively until none is left and measure the linear size of the remaining structure. When averaged over all pairs of vertices with separations $`d>10`$, it comes out to be $`d_b4.5`$, well consistent with the value obtained previously for $`d=26`$ only. Due to this blob structure, the mass-distance relation increases abruptly across $`d=4`$ as shown in Fig.3b.

Next, we measure the average mass of blob, that is, the number of vertices inside a blob for a given graph of the shortest pathways with separation $`d`$, averaged over all pairs of vertices with the same separation. We find that the average blob mass is broadly distributed in the range $`3<m_b<23`$. In particular, relatively heavy blob masses, $`m_b=15`$$``$$`23`$, mainly comes from the graphs whose linear size is $`d=8`$$``$$`14`$. Due to those blobs with heavy mass, the mass-distance relation exhibits a hump, and decreases at around $`d=14`$$``$$`16`$, beyond which, the mass $`(d)`$ increases linearly by the presence of singly connected vertices. In short, the anomalous behavior in the mass-distance relation is due to the presence of a compact and localized blob structure in the topology of the shortest pathways between a pair of vertices for the metabolic network of eukaryotes. We have checked the mass-distance relations and the graphs of the shortest pathways for other networks belonging to the class I such as the PIN2 and metabolic networks for other organisms, and found that such topological features are generic, generating the anomalous behavior in the mass-distance relation. It still remains a challenge to derive the BC exponent $`\eta 2.2`$ analytically from such structures.

For the class II, the mass depends on distance linearly, $`(d)Ad`$ for large $`d`$ (Fig.3c-d). Despite the linear dependence, the shortest pathway topology for the case of $`A>1`$ is more complicated than that of the simple tree structure where $`A=1`$. Therefore, the SF networks in the class II are subdivided into two types, called the class IIa and IIb, respectively. For the class IIa, $`A>1`$ and the topology of the shortest pathways includes multiply connected vertices (Figs.4b and 4c), while for the class IIb, $`A1`$ and the shortest pathway is almost singly connected (Fig.4d). Examples in real world networks in the class IIa are the Internet at the AS level ($`A4.5`$) and the metabolic network for archaea ($`A2.0`$), while that in the class IIb is the WWW ($`A1.0`$).

Let us examine the topological features of the shortest pathways for the networks in the class IIa and IIb more closely. First, for the class IIa, we visualize in Fig.4b a shortest pathway in the Internet system between a pair of vertices separated by 10 edges, the farthest separation. It contains a blob structure, but the blob is rather extended as $`d_b=5`$, comparable to the maximum separation $`d=10`$. We obtain $`d_b=5`$ for $`d=11`$ for another system. For comparison, $`d_b4.5`$ for $`d=26`$ in the class I. Moreover, the feature-less mass-position dependence $`m(r;d)`$ we have found implies that while most blobs are located almost in the middle of the shortest pathways, which seems to be caused by the geometric effect, there are a finite number of blobs located at the verge of the shortest pathways. Note that $`m(r;d)=m(dr;d)`$ since the Internet is undirected. Owing to the extended structure and the scattered location of the blob, the mass-distance relation exhibits the linear behavior, $`(d)Ad`$ with $`A4.5`$. The extended blob structure is also observed in the metabolic network for archaea (Fig.4c). Since the network in this case is directed, the symmetry in $`m(r;d)`$ does not hold. However, the blob structure extends to almost one half of maximum separation, and the shortest pathways are very diverse, so that their topological property such as the mass-distance relation $`(d)`$ is similar to that of the Internet.

The WWW is an example belonging to the class IIb. For this network, the mass-distance relation exhibits $`(d)1.0d`$, suggesting that the topology of the shortest pathway is almost singly connected, which is confirmed in Fig.4d. When a SF network is of tree structure, one can solve the distribution of BC running through each edge analytically, and obtain the BC exponent to be $`\eta =2`$. A derivation of this exact result is presented in the Appendix.

Comparison of the Resilience under Attack
So far we have investigated the topological features of the shortest pathways of SF networks of each class. Then what would be distinct physical phenomena originated from such different topological features? Associated with this question, we investigate a problem of the resilience of network under a malicious attack. It is known that SF networks are extremely vulnerable to the intentional attack to a few vertices with high degree, while it is very robust to random failures attack ; havlin . To compare how vulnerable a network in each different class is under such attacks, we first construct a directed network whose numbers of vertices and edges, and the degree distribution are identical to those of the WWW (xiv), but whose BC exponent is $`2.2`$. It can be generated, for example, by following the stochastic rule introduced in the directed static model load . For both the WWW in real world and the artificial model network, we remove vertices in the descending order of BC successively. As vertices are removed, both the mean distance $`d`$ between two vertices, known as the diameter, and the relative size of the giant cluster $`S`$ are measured as a function of the fraction of removed vertices $`f`$. As can be seen in Fig.6a, the diameter of the WWW with $`\eta =2.0`$ (class IIb) increases more rapidly than that with $`\eta =2.2`$ (class I) and shows discrete jumps while vertices are removed. Also the relative size of the largest cluster decreases more rapidly for $`\eta =2.0`$ than for $`\eta 2.2`$ (Fig.6b). This behavior arises from the fact that the shortest pathway consists of mainly singly connected vertices for the class IIb, so that there is no alternative pathways with the same distance when a single vertex lying on the shortest pathway is removed. For the Internet in real world with $`\eta 2.0`$ in the class IIa and an artificial network with $`\eta 2.2`$ with the same numbers of vertices and edges and the identical degree distribution, the differences in the diameter $`d`$ and in the relative size $`S`$ of the largest cluster appear to be rather small (Figs.6c-d), in comparison to the case of the WWW (Figs.6a-b). This is because the shortest pathways are multiply connected for the class IIa.

Conclusions
In conclusion, we have found that the betweenness centrality can determine the universal behavior of SF networks. By examining a variety of real world and artificial SF networks, we have observed two distinct universality classes whose BC exponents are $`\eta 2.2(1)`$ (class I) and $`2.0`$ (class II), respectively. The mass-distance relation is introduced to characterize the topological features of the shortest pathways. It shows a hump for the class I networks due to compact and localized blobs in the shortest pathway topologies, while it is roughly linear for the class II ones which are more or less tree-like. The class II networks can further be divided into two types depending on whether the shortest pathway topology contains diversified pathways (class IIa) or mostly singly connected ones (class IIb). Distinct features of the resilience under attack arising from the different topologies of the shortest pathways are also identified. Since SF networks show the small world property, the topology of the shortest pathways should be of relevance for characterizing the network geometry. Indeed the mass-distance relations for different universality classes show different behaviors. Such a relation between the universality class and the topological features of the shortest pathways may be understood from the perspective of the fact that the geometric fractal structure of the magnetic domains in equilibrium spin systems at criticality can classify the universality classes. Further characterizations in static and dynamic properties and possible evolutionary origin of the universality classes are interesting questions left for future study.

Acknowledgments
This work was supported by the Korean Research Foundation (Grant No. 01-041-D00061) and BK21 program of MOE, Korea.

APPENDIX
Here we present the analytic derivation of the BC distribution for a tree structure, however, the derivation is carried out for the edge BC rather than the vertex version333See Szabó et al. (cond-mat/0203278, 2nd version) for a mean field treatment of the vertex BC problem on trees.. The edge betweenness centrality is defined on edges as in Eq.(1), with the subscript $`k`$ now denoting a bond. Without any rigorous proof, we assume that the distributions of vertex BC and edge BC behave in the same manner particularly on tree structures, which is confirmed by numerical simulations. We have also checked the identity between the vertex BC and the edge BC distributions for a deterministic model of scale-free tree introduced by Jung et al. jung , which will be published elsewhere.

We consider a growing tree network such as the BA type model with $`m=1`$, where a newly introduced vertex attaches an edge to an already existing vertex $`j`$ with the probability proportional to its degree as $`(k_j+a)/_{\mathrm{}}(k_{\mathrm{}}+a)`$. Then the network consists of $`N(t)=t+1`$ vertices and $`L(t)=t`$ edges at time $`t`$. The stationary degree distribution is of a power law with $`\gamma =3+a`$ bu ; porto . Each edge of a tree divides the vertices into two groups attached to either sides of the edge. Let $`P_s(m,t)`$ be the probability that the edge born at time $`s`$ bridges a cluster with $`m`$ vertices on the descendant side and another with remaining $`t+1m`$ vertices on the ancestor side. Due to the tree structure, the BC running through that edge born at $`s`$ is given as $`g=2m(t+1m)`$, independent of the birth time $`s`$. The probability $`P_s(m,t)`$ evolves as a new vertex attaches to one of the two clusters. The rate equation for this process is written as

$$P_s(m,t+1)=r_1(m,t)P_s(m,t)+r_2(m1,t)P_s(m1,t),$$
(A1)
where $`r_1(m,t)`$ is the probability that a new vertex attaches to the cluster with $`(t+1m)`$ vertices on the ancestor side, and $`r_2(m1,t)`$ with $`(m1)`$ vertices on the descendant side. They are given explicitly as

$$r_1(m,t)=1r_2(m,t)=\frac{(2+a)(tm)+at+1}{2t+a(t+1)}.$$
(A2)
Since the amount of the BC on the edge $`s`$ is independent of the birth time, we introduce $`P(m,t)`$,

$$P(m,t)=\frac{1}{t}\underset{s=1}{\overset{t}{}}P_s(m,t),$$
(A3)
which is the probability for a certain edge to locate between two clusters with $`m`$ and $`t+1m`$ vertices averaged over its birth time. The BC on that edge is still given by $`2m(t+1m)`$. In terms of $`P(m,t)`$, Eq.(A1) can be written as

$`(t+1)P(m,t+1)`$ $`=`$ $`r_1(m,t)tP(m,t)`$ (A5)
$`+r_2(m1,t)tP(m1,t).`$
In the limit of $`t\mathrm{}`$, one may rewrite $`P(m,t)`$ in a scaling form, $`P(m,t)=𝒫(m/t)`$ and then Eq.(A5) is rewritten as

$$(t+1)𝒫(x)t𝒫(x)x\frac{d𝒫(x)}{dx}𝒫(x)$$

where $`x=m/t`$ and the approximation $`𝒫(x1/t)𝒫(x)(1/t)d𝒫(x)/dx`$ has been used. From this we obtain that

$$𝒫(x)\frac{1}{x^2},$$
(A6)
independent of the tuning parameter $`a`$. Using $`g=2(t+1m)m2t^2x`$ for large $`t`$ and finite $`m`$, Eq.(A6) becomes

$$P_B(g)\frac{1}{g^2}.$$
(A7)
Thus $`\eta =2`$ is obtained for the tree structure, independent of $`\gamma >2`$. General finite size scaling relations for $`P_B(g)`$ are discussed in Ref.kolkata .

Figure Legends

Fig. 1 The BC distributions of real world networks.

(a) Networks belonging to the class I: Co-authorship network ($`\times `$), Core of PIN of yeast by Ito et al. (+), and metabolic network of EN ($``$). The solid line is a fitted line with a slope $`2.2`$. (b) Networks belonging to the class II: WWW of nd.edu ($``$) and Internet AS as of October, 2001 ($`\mathrm{}`$). The solid line has a slope $`2.0`$.

Fig. 2 Comparison of the BC distributions for the two classes.

BA model with $`m=1`$, 2, and 3 are simulated for large system size, $`N=3\times 10^5`$, averaged over 10 configurations. The dotted line has a slope $`2.0`$ and the dashed one, $`2.2`$.

Fig. 3 The mass-distance relation $``$$`(d)`$.

(a) Core of PIN of yeast obtained by Ito et al.. (b) Metabolic networks of eukaryotes. Data are averaged over all 5 organisms in Ref.metabolic . Note that in this case we count only substrates for $`(d)`$. (c) Internet AS as of October, 2001. (d) WWW of nd.edu.

Fig. 4 Topology of the shortest pathways.

(a) The metabolic network of EN (eukaryote) of length 26. (b) The Internet AS of length 10. (c) The metabolic network of Methanococcus jannaschii (archae) of length 20. (d) WWW of nd.edu of length 20. In (a) and (c), circles denote substrates and rectangles denote intermediate states.

Fig. 5 The mass density.

$`m(r;d)`$ for EN with $`d=26`$. Circles denote substrates and rectangles intermediate states.

Fig. 6 Attack vulnerability of the scale-free networks.

The WWW ($`\eta =2.0`$) ($`\mathrm{}`$) and the artificial directed SF network with $`\eta =2.2`$ ($`\mathrm{}`$), the Internet ($`\eta =2.0`$) ($``$) and the artificial undirected SF network with $`\eta =2.2`$ ($``$): Changes in network diameter (a, c) and the relative size of the largest cluster (b, d) are shown as a function of $`f`$, the fraction of removed vertices measured in percent (%).

Table Legend

Table 1 Natures of diverse SF networks.

Tabulated for each network are the size $`N`$, the mean degree $`k`$, the degree exponent $`\gamma `$, and the betweenness centrality exponent $`\eta `$.


text:
# 1 Introduction

## 1 Introduction

The dynamic questions in solving some soliton problems are relegated to the background, at least, at the initial stage, and the problem of determining the static stresses in the system becomes the key point in the investigation. This problem directly reduces to finding quasi-static soliton configurations described using a set of non-Goldstone collective coordinates. In this context, the problem of determining the intersoliton potential in various models can be mentioned. The non-Goldstone collective coordinates in these problems are the separation parameter and also the variables describing the relative orientation of the soliton pair in the internal space in multidimensional cases. Furthermore, the non-Goldstone coordinates, which in essence are nonlinear amplitudes, describe the dynamics of internal soliton modes. The problem of determining the quasi-static configurations in these problems has been repeatedly discussed, and the ways to solve it were related to various model constructions. In this connection, we mention the investigations of the skyrmion interaction begun as early as 1962 . Moreover, the soliton interaction was also studied for some other models , where some important properties of intersoliton potentials were established. At, the same time, it is clear that an approach based on model semiphenomenological notions does not ensure a comprehensive investigation of the problem, for which a consistent theoretical scheme is required.

Touching upon this problem, we indicate some difficulties related to introducing non-Goldstone collective coordinates. For example, the notion of the distance between extended objects undergoing mutual deformation acquires its unambiguous visual meaning only if these objects are at sufficiently large distances from each other. Therefore, a consistent theory must not only describe the process of deformation of solitons as they approach each other but also answer the question of interpreting the ambiguity in determining the intersoliton distance at the stage where the solitons lose their identity as a result of this deformation. The difficulties in constructing a theoretical scheme suitable for describing intersoliton interaction that are due to the above-mentioned factor have been extensively mentioned in the literature and were most distinctly stated in . Some simple considerations were later presented in  in favor of the universal variational approach to the method of collective variables, which permits stating various problems of soliton dynamics in the framework of the Hamiltonian procedure. In , the general equation for quasi-static configurations determined by an arbitrary set of collective coordinates was also derived. It is essential that the resulting equation is invariant under reparameterizing the collective variables, which ensures that the physical quantities are independent of the way in which the intersoliton distance is defined.

In Sec. 2, we present a new covariant method for deriving an equation for the quasi-static ansatz using basic notions related to Riemannian manifolds. The following two sections are devoted to elaborating solution methods for this equation in some specific cases. In Sec. 3, the problem of determining the configurations and static stresses of the $`\mathrm{\Phi }^4`$-kink whose deformation is produced by the excitation of the internal soliton mode is considered. This problem can be solved exactly. Similar problems were considered in many papers. The corresponding references can be found in . It is unlikely that there is a model with a nontrivial potential that admits an exact solution of the problem of determining the intersoliton interaction. A solution of this problem that allows calculating an asymptotic expression for the intersoliton interaction can most likely be found for an arbitrary model. An example of a calculation of this kind is given in Sec. 4, where the equation for the quasi-static ansatz is used to determine the kink interaction for the nonlinear Klein-Gordon equation in the domain of large intersoliton distances. We stress that, apart from independent interest, knowing the intersoliton potential in the asymptotic domain permits introducing a physically acceptable gauge in the numerical analysis of the equation throughout the range of intersoliton distances. In Sec. 5, we discuss the investigation results and indicate the main elements of the suggested theoretical scheme.

## 2 Equation for the quasi-static ansatz

We consider the question of finding the quasi-static ansatz for a soliton system described by an arbitrary set of collective variables. Let

$$L=g_{ik}\dot{q}^i\dot{q}^kH(q,0)$$

be a nondegenerate Lagrangian quadratic with respect to velocities that admits soliton solutions, let $`q^i`$ be the field coordinates, and let $`i\{i,x\}`$ be the set of discrete indices and spatial coordinates. It is convenient to interpret the basic configuration space as a Riemannian manifold with a metric $`g_{ik}`$ generated by the kinetic term. In this case, if $`Q`$ are collective variables, then the quasi-static ansatz $`q_c^i(Q)`$ describing the deformed soliton state can be interpreted as a submanifold parameterized by the variables $`Q`$. We assume that this ansatz is determined by the stationarity condition for the Hamilton function defined on $`q_c^i(Q)`$ for the system in question relative to arbitrary infinitesimal variations $`\delta q_c^i(Q)`$ orthogonal to the manifold surface, i.e.,

$$\frac{H(q_c,0)}{q_c^i}\delta q_c^i=0,$$
(1)
where the variations $`\delta q_c^i(Q)`$ satisfy the conditions

$$g_{ik}(q_c)\frac{q_c^i(Q)}{Q}\delta q_c^k(Q)=0.$$

To cancel the variations in Eq. (1), we introduce the operator $`𝒫`$ of projection onto the subspace tangent to the submanifold $`q_c^i(Q)`$ at the point $`Q`$

$$𝒫_k^i=g_{kl}(q_c)\frac{q_c^l}{Q}g^{QQ}\frac{q_c^i}{Q},𝒫^2=𝒫,$$
(2)
where

$$g_{QQ}=\frac{q_c^i}{Q}g_{ik}(q_c)\frac{q_c^k}{Q}$$

is the metric on the submanifold $`q_c^i(Q)`$ and $`g^{QQ}`$ is the matrix inverse to $`g_{QQ}`$. Because the variations $`\delta q_c^i(Q)`$ in the orthogonal directions to $`q_c^i(Q)`$ are arbitrary, Eq. (1) can then be brought to the form

$$(1𝒫)\frac{H(q_c,0)}{q_c}=0.$$

Substituting expression (2) for the operator $`𝒫`$ in this formula, we obtain

$$\frac{H(q_c,0)}{q_c^i}=g_{ik}(q_c)\frac{q_c^k}{Q}g^{QQ}\frac{E(Q)}{Q},$$
(3)
where $`E(Q)H(q_c(Q),0)`$ is the energy of the quasi-static configuration. By the definition of the tensor $`g^{QQ}`$, the right-hand side of Eq. (3), which takes the existence of static stresses in the collective subsystem into account, is invariant under an arbitrary nondegenerate change of variables $`QQ^{}(Q)`$. This indicates that the variational procedure under consideration is covariant. In this connection, we note that in relation to the problem of determining two-soliton configurations, for example, Eq. (3) exactly describes these configurations in the domain of distances where the deformed solitons lose their identity. The choice of the gauge only determines the specific form of the parameterization describing the configurations. For convenience, it can be required that the separation parameter coincide with the natural intersoliton distance if this distance notably exceeds the soliton size. We note that if the collective coordinates are Goldstone variables, i.e., if the coordinates $`Q`$ coincide with the degeneration parameters of vacuum solutions, then the energy $`E`$ does not depend on $`Q`$, and Eq. (3) becomes the equation $`H(q_c,0)/q_c^i=0`$ for static solutions.

## 3 Static stresses of the deformed $`\mathrm{\Phi }^4`$-kink.

We use Eq. (3) to determine the configuration of the deformed kink in the $`\mathrm{\Phi }^4`$ model. The Lagrangian of the model and the corresponding static kink solution $`u_c(xX)`$ defined to within the spatial coordinate $`X`$ have the forms

$$L=(1/2)𝑑x\left(\dot{\mathrm{\Phi }}^2(x)\mathrm{\Phi }^2(x)(\mathrm{\Phi }^2(x)1)^2\right),u_c(x)=\mathrm{tanh}x.$$
(4)
The eigenfunctions of the operator $``$ describing linear fluctuations on the background of the kink solution include the mode $`\sigma (x)`$ belonging to the discrete-spectrum states,

$$\sigma (x)=3\sigma (x),\sigma (x)=\sqrt{3/2}\frac{\mathrm{tanh}x}{\mathrm{cosh}x},=\frac{^2}{x^2}+46\mathrm{cosh}^2x.$$
(5)
In describing nonstationary states, it is usually assumed that the $`\mathrm{\Phi }^4`$-kink possesses an internal degree of freedom whose small oscillations can be identified with the mode $`\sigma (x)`$. Therefore, apart from the translational coordinate $`X`$, the set of collective coordinates describing the kink state must also include the variable $`\tau `$ related to the excitation of the internal mode of the $`\mathrm{\Phi }^4`$-kink. For this set of collective coordinates, we can expect that the directing vectors $`\mathrm{\Phi }_c^{}\mathrm{\Phi }_c/x`$ and $`\mathrm{\Phi }_{c\tau }\mathrm{\Phi }_c/\tau `$ constructed on the solutions of Eq. (3) are orthogonal in the tangent space,

$$𝑑x\mathrm{\Phi }_c^{}(x)\mathrm{\Phi }_{c\tau }(x)<\mathrm{\Phi }_c^{}\mathrm{\Phi }_{c\tau }>=0.$$
(6)
Also, taking into account that the energy is independent of the cyclic variable, we can represent Eq. (3) in the case under consideration as

$$\mathrm{\Phi }_c^{\prime \prime }(x)+2\mathrm{\Phi }_c(x)(\mathrm{\Phi }_c^2(x)1)=\frac{\mathrm{\Phi }_c(x)}{\tau }<\mathrm{\Phi }_{c\tau }^2>^1\frac{E}{\tau }.$$
(7)
The right-hand side of Eq. (7) is invariant under the reparameterization $`\tau \tau ^{}(\tau )`$. If the collective variable is fixed using the gauge

$$<\mathrm{\Phi }_{c\tau }^2>^1E/\tau =1,$$
(8)
then Eq. (7) becomes the nonlinear diffusion equation

$$\frac{\mathrm{\Phi }_c(x)}{\tau }\mathrm{\Phi }_c^{\prime \prime }(x)+2\mathrm{\Phi }_c(x)(\mathrm{\Phi }_c^2(x)1)=0.$$
(9)
Equation (9) was investigated in , where its exact solutions were found. In particular, the solutions of the type of a deformed kink that tend to $`+1`$ as $`x+\mathrm{}`$ can be written in the form

$$\mathrm{\Phi }_c(x,\tau )=\frac{e^{2x}1}{e^{2x}+1\pm e^{x3\tau }}.$$
(10)
Integrating, we can easily verify that these solutions satisfy orthogonality condition (6) and gauge condition (8). We now note that the physical domain of small deformations, where solutions (10) only slightly differ from static kink solution (4), corresponds to $`\tau \mathrm{}`$. In other words, the collective variable $`\tau `$ defined by gauge condition (8) is inconvenient from the physical standpoint. It plays only an intermediate role as it permits representing Eq. (7) in the form of a nonlinear diffusion equation. It is therefore advisable to combine the two solutions in (10) to form one solution,

$$\mathrm{\Phi }_c(x,q)=\frac{e^{2x}1}{e^{2x}+1\sqrt{6}qe^x}u_c(x)+q\sigma (x)+O(q^2),$$
(11)
where the new collective variable $`q`$, which takes both positive and negative values, has the meaning of the amplitude of the internal mode under small deformations . The configuration energy $`E(q)H(\mathrm{\Phi }_c,0)`$ determining the static stresses of the deformed kink has the form

$$E(q)=\frac{4}{3}+\frac{\overline{q}^2}{1\overline{q}^2}\frac{\overline{q}^3}{(1\overline{q}^2)^{3/2}}\mathrm{arccos}\overline{q}\frac{4}{3}+\frac{3}{2}q^2+O(q^3),\overline{q}=\sqrt{3/2}q.$$

Here, the first term is the kink mass, and the quadratic term with respect to $`q`$ describes the potential of a harmonic oscillator with the internal mode frequency $`\omega ^2=3`$. The range of the collective variable where the energy remains real is confined to the interval $`\sqrt{2/3}<q<+\sqrt{2/3}`$. Solution (11) has no singularities with respect to the spatial coordinate in this range of $`q`$. We also stress that the choice of the physical collective variable tending to the amplitude of the internal mode in the limit of small deformations is ambiguous. Namely, it is determined for convenience and is connected with explicit form (10) of the solution for the deformed kink in the situation in question.

## 4 Kink interaction.

An important role in studying the properties of soliton states is played by the solutions of Eq. (9) that are defined by the asymptotics $`\mathrm{\Phi }_c(x)1`$ as $`x\pm \mathrm{}`$. These solutions describe the kink-antikink configurations and permit establishing the properties of the intersoliton interaction. In this case, the collective variable relates to the separation parameter $`r`$, $`2r`$ is the intersoliton distance. The excitations of internal modes are not taken into account. Unfortunately, it seems that exact solutions of Eq. (9) with given boundary conditions have not yet been found . Approximate solutions in the domain of large distances can be obtained in the general form if the corresponding equation

$$\mathrm{\Phi }_c^{\prime \prime }(x)+U^{}(\mathrm{\Phi }_c(x))=\frac{\mathrm{\Phi }_c(x)}{r}<\mathrm{\Phi }_{cr}^2>^1\frac{E}{r},$$
(12)
for the nonlinear Klein-Gordon model is considered instead of (7). Here, $`U(\mathrm{\Phi })`$ is the potential of the model. If the origin is placed at the center of the bisoliton pair and its symmetry $`\mathrm{\Phi }_c(x,r)=\pm \mathrm{\Phi }_c(x,r)`$ is taken into account, then the entire configuration can be described by considering the domain $`x0`$. The even functions describe the kink-antikink states and the odd ones describe the kink-kink states (in the case of a periodic potential $`U(\mathrm{\Phi })`$). If the intersoliton distance is large, then to the left of the point $`x=0`$, the ansatz $`\mathrm{\Phi }_c(x,r)`$ only slightly differs from the one-kink solution $`u_c(x+r)`$ of the static Klein-Gordon equation

$$u_c^{\prime \prime }(x)+U^{}(u_c(x))=0.$$

Therefore, representing the ansatz $`\mathrm{\Phi }_c(x,r)`$ in the form

$$\mathrm{\Phi }_c(x,r)=u_c(x+r)+\eta (x+r,r),x0,$$
(13)
where $`\eta (x+r,r)`$ tends to zero as $`r\mathrm{}`$, and taking into consideration that in the domain of large values of $`r`$,

$$<\mathrm{\Phi }_{cr}^2>2\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑xu_c^2(x)2M,$$
(14)
where $`M`$ is the kink mass, we derive an equation for $`\eta (x)`$ from (12),

$$\eta (x)=\frac{1}{2M}\frac{E}{r}u_c^{}(x),xr,\frac{^2}{x^2}+U^{\prime \prime }(u_c(x)).$$
(15)
Here, $`E/r`$ is a small parameter in the domain of large values of $`r`$. Multiplying Eq. (15) by the zero mode $`u_c^{}(x)`$ of the operator $``$ and integrating by parts from $`\mathrm{}`$ to $`x`$ , we can reduce the order of the equation. As a result, we obtain the linear first-order equation for $`\eta (x)`$

$$u_c^{}(x)\eta ^{}(x)+u_c^{\prime \prime }(x)\eta (x)=\frac{M(x)}{2M}\frac{E}{r},M(x)=\underset{\mathrm{}}{\overset{x}{}}𝑑x^{}u_c^2(x^{}).$$
(16)
The solution of Eq. (16) can be written as

$$\eta (x)=u_c^{}(x)\left(\frac{\eta (r)}{u_c^{}(r)}+\frac{1}{2M}\frac{E}{r}\underset{x}{\overset{r}{}}𝑑x^{}\frac{M(x^{})}{u_c^2(x^{})}\right).$$
(17)
The function $`\eta (r)`$ playing the role of the integration constant in (17) is found from the boundary conditions. For the kink-antikink and kink-kink pairs, these conditions have the respective forms $`\mathrm{\Phi }_c^{}(0,r)=0`$ and $`\mathrm{\Phi }_c(0,r)=0`$. Using these conditions and Eq. (16), we obtain the following expressions for the boundary values $`\eta (r)`$ and $`\eta ^{}(r)`$ in these two cases:

$$\eta ^{}(r)=u_c^{}(r),\eta (r)=\frac{1}{u_c^{\prime \prime }(r)}\left(\frac{M(r)}{2M}\frac{E}{r}u_c^2(r)\right)$$
(18)
for the kink-antikink pair and

$$\eta (r)=u_c(r),\eta ^{}(r)=\frac{1}{u_c^{}(r)}\left(\frac{M(r)}{2M}\frac{E}{r}+u_c^{\prime \prime }(r)u_c(r)\right)$$
(19)
for the kink-kink pair. Solution (17) describes the distortion of the kink under the effect of the second soliton when the distance between the kinks substantially exceeds their sizes and depends on the small parameter $`E/r`$, which determines the interaction of the solitons. We note that formulas (17-19) include terms of higher order than the leading asymptotic terms with respect to $`r`$. However, because the function $`\eta (x)`$ contains competing spatial exponentials, it is inadvisable to select these terms in considering the field variable $`\eta (x)`$. Formulas (17-19) should therefore be used to correctly determine the leading approximation term for the integrated characteristic, the interkink potential $`E(r)`$. Because $`(1/2)u_c^2=U(u_c)`$, substituting (13) in the formula for the energy $`E(r)=H(\mathrm{\Phi }_c,0)`$, expanding with respect to $`\eta (x)`$, and integrating by parts with regard to Eq. (16) result in

$$\begin{array}{c}E(r)2M=2\underset{r}{\overset{\mathrm{}}{}}𝑑xu_c^2(x)+\frac{1}{2M}\frac{E}{r}\underset{\mathrm{}}{\overset{r}{}}𝑑x\eta (x)u_c^{}(x)\hfill \\ \hfill +\eta (r)(2u_c^{}(r)+\eta ^{}(r)).\end{array}$$
(20)
Substituting expression (17) for $`\eta (x)`$ in (20) and using boundary values (18) and (19) for the kink-antikink and kink-kink pairs, we see that the linear terms with respect to $`E/r`$ cancel. In this case, relation (20) becomes

$$E(r)2M=A(r)+B(r)\left(\frac{1}{2M}\frac{E}{r}\right)^2,$$
(21)
where the functions $`A(r)`$ and $`B(r)`$ are given by the formulas

$`A(r)`$ $`=`$ $`{\displaystyle \frac{u_c^3(r)}{u_c^{\prime \prime }(r)}}2{\displaystyle \underset{r}{\overset{\mathrm{}}{}}}𝑑xu_c^2(x),`$
$`B(r)`$ $`=`$ $`{\displaystyle \frac{M^2(r)}{u_c^{}(r)u_c^{\prime \prime }(r)}}+{\displaystyle \underset{\mathrm{}}{\overset{r}{}}}𝑑x{\displaystyle \frac{M^2(x)}{u_c^2(x)}}`$
for the kink-antikink pair and

$`A(r)`$ $`=`$ $`u_c^2(r){\displaystyle \frac{u_c^{\prime \prime }(r)}{u_c^{}(r)}}+2{\displaystyle \underset{r}{\overset{\mathrm{}}{}}}𝑑xu_c(x)u_c^{\prime \prime }(x),`$
$`B(r)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{r}{}}}𝑑x{\displaystyle \frac{M^2(x)}{u_c^2(x)}}`$
for the kink-kink pair. The asymptotic expression for the kink solution has the form

$$u_c(x)|_{x\pm \mathrm{}}u_\pm a\mathrm{exp}(mx),$$

where $`m=\sqrt{U^{\prime \prime }(u_\pm )}`$ is the pion mass. In a case of an antisymmetric kink-kink configuration the vacuum constant $`u_+`$ is equal to zero. Taking into consideration that $`M(\mathrm{})=M`$ and $`M(\mathrm{})=0`$ and calculating the leading approximation terms with respect to $`e^{mr}`$ for the coefficients $`A(r)`$ and $`B(r)`$, we see that the functions $`A(r)`$ vanish and $`B(r)`$ coincide to within the sign. Therefore the leading approximation for Eq. (21) has the form

$$E(r)2M=\frac{1}{8m^3a^2}e^{2mr}\left(\frac{E}{r}\right)^2.$$

whence we find the interaction energy of the bisoliton pair at large distances,

$$E(R)=2M2ma^2e^{2mr}.$$
(22)
This formula describes attraction (repulsion) for the kink-antikink (kink-kink) configuration. The expression for the energy in (22) in the case of the $`\mathrm{\Phi }^4`$ and sine-Gordon systems was found in  using an approximate equation of type (12) whose right-hand side was imitated with a phenomenological delta-shaped source .

To determine the soliton configuration and interaction throughout the range of distances, the gauge should be fixed, and Eq. (12) should be solved numerically. We note that the representation of the solution in form (13) partly fixes the gauge in the domain of large intersoliton distances. Formulas (14) and (22) obtained using expansion (13) can be used to introduce the natural gauge

$$<\mathrm{\Phi }_{cr}^2>^1\frac{E}{r}=\pm \frac{2m^2a^2}{M}e^{2mr},$$
(23)
ensuring the coincidence of the doubled parameter $`r`$ with the true intersoliton distance in the asymptotic domain $`r\mathrm{}`$. After gauge (23) is fixed, Eq. (12) takes the form of a nonlinear diffusion equation,

$$\mathrm{\Phi }_c^{\prime \prime }(x)+U^{}(\mathrm{\Phi }_c(x))=\pm \frac{2m^2a^2}{M}e^{2mr}\frac{\mathrm{\Phi }_c(x)}{r}.$$
(24)
Equation (24) can be used to find the corrections to potential (22) in the asymptotic domain for the given gauge, and its numerical solution permits reproducing the entire set of physically admissible two-soliton configurations. We stress that solving this equation in the case of the kink-antikink configuration corresponds to describing diffusion in the direction of positive ”time” $`r`$. This process is stable, and the mutual approach in the kink-antikink pair terminates with its ”annihilation” under which $`\mathrm{\Phi }_c(x)`$ goes into the vacuum constant . Actually, the calculation of the evolution of the kink-antikink pair performed in  in relation to the investigation of the nonlinear diffusion equation is an example of rigorous determination of the entire set of two-soliton configurations in problems related to studying intersoliton interaction on the basis of Eq. (3).

Constructing the kink-kink configuration reduces to studying diffusion in the direction of negative ”time” $`r`$. In the domain of large intersoliton distances, the presence of nonlinearity in Eq. (24) stabilizes the instability characteristic of this procedure in the case of the linear theory. The stabilization weakens during the mutual approach, and the limiting field configuration is determined by the point where the stability of the process terminates.

## 5 Conclusion

We have covariantly derived an equation for a quasi-static ansatz on the basis of the variational approach. A characteristic feature of this approach is the invariance under reparameterization with a change of collective variables. The cyclic variables or the variables related to the degeneration of vacuum states are introduced in the standard way. The arbitrariness in the choice of the non-Goldstone variables can be eliminated based on some physical considerations in the domain where these variables can be given an immediate physical interpretation. Otherwise they are chosen for convenience.

As an example, we have considered the static deformations of the $`\mathrm{\Phi }^4`$-kink related to excitation of the internal mode. This problem admits an exact solution. The configurations of the deformed kink, its energy, and the physically admissible range of the amplitude of nonlinear oscillations have been determined. Determining the intersoliton potential for the nonlinear Klein-Gordon equation has also been touched upon. Linearizing the equation for the quasi-static ansatz in the domain of large distances permits establishing the deformation of the kinks and the asymptotic character of their interaction, which allows introducing the natural gauge and representing the equation for the quasi-static ansatz as a nonlinear diffusion equation suitable for numerical determination of the entire set of physically admissible two-soliton configurations

In relation to determining the intersoliton interaction, we stress that although two-soliton configurations have been comprehensively studied using various qualitative approaches for the majority of actual models, investigations on the basis Eq. (3) permit including the exact effect of soliton deformation at small distances, i.e., in the domain most sensitive to the defects characteristic of artificial field constructions. This is particularly important, for instance, in the Skyrme model because the nucleon-nucleon potential in this domain has an insignificant minimum forming the nuclear coupling energy. Although Eq.(3) does not take the spin effects into account, the study of skyrmion-skyrmion interaction on the basis of this equation is a useful stage in the procedure of determining the nucleon-nucleon potential.
The author is grateful to P. 0. Mchedlov-Petrosyan and V.V. Gann for a number helpful discussions.


text:
# Fully-Renormalized QRPA fulfills Ikeda sum rule exactly

## Abstract

The renormalized quasiparticle-RPA is reformulated for even-even nuclei using restrictions imposed by the commutativity of the phonon creation operator with the total particle number operator. This new version, Fully-Renormalized QRPA (FR-QRPA), is free from the spurious low-energy solutions. Analytical proof is given that the Ikeda sum rule is fulfilled within the FR-QRPA.

The well-known Random Phase Approximation (RPA)  has found many applications in different fields of physics, from chemistry and condensed matter  to relativistic field theory . In nuclear physics RPA has been extensively exploited to model properties of the excited nuclear states that allows to calculate intensities of various nuclear reactions, including the double beta decay (see reviews ).

The neutrinoless double beta decay ($`0\nu \beta \beta `$-decay), which violates the total lepton number by two units, is a sensitive low-energy probe for new physics beyond the Standard Model . The observation of the $`0\nu \beta \beta `$-decay would give unambiguous evidence that at least one of the neutrinos is a Majorana particle with non-zero mass . The current experimental upper limits on the $`0\nu \beta \beta `$-decay half-life impose stringent constraints, e.g., on the parameters of Grand Unification and supersymmetric extensions of the Standard Model.

It is important that the same nuclear structure methods can be applied to calculate intensities of $`2\nu \beta \beta `$-decay, which is a second order process allowed within the Standard Model. In this case the quality of the RPA calculations can be directly checked by comparison with the corresponding experimental data available at the moment for a number of nuclei . Most of the nuclei are open-shell ones and one uses the quasiparticle version of RPA (QRPA) which takes into account nucleon pairing correlations.

The QRPA is able to reproduce the experimental data , but with the strength of a particle-particle interaction which is very close to the point where QRPA collapses. This well-known drawback of the RPA is due to overestimation of ground state correlations when one uses the quasiboson approximation (QBA) for the bifermionic operators and the Pauli exclusion principle (PEP) is violated. The proximity to the collapse makes the entire QRPA scheme unstable and calculation results unreliable. To cure this problem, the so-called renormalized QRPA (RQRPA) was invented , which takes into account PEP in an approximate way. Namely, the expectation values of the quasiparticle number operators are consistently treated to be non-vanishing in the RQRPA vacuum. The RQRPA does not collapse for physical values of the particle-particle interaction strength. It has been shown within schematic models that by including PEP in the QRPA good agreement with the exact solution of the many-body problem can be achieved even beyond the critical point of the standard QRPA (see, e.g.  and references therein). The RQRPA has been used in previous studies of the double beta decay .

Nevertheless, it has been shown within the modern versions of RQRPA that the model-independent Ikeda sum rule (ISR) is violated . The Ikeda sum rule states that the difference between the total Gamow-Teller strengths in the $`\beta ^{}`$ and $`\beta ^+`$ channels is $`3(NZ)`$ . This drawback of the model is quite serious, because the derivation of the ISR is based just upon two principles: closure relation and average conservation of the particle number. Violation of the ISR immediately means violation of one (or both) of these important principles. The drawback has not been cured neither in the self-consistent QRPA  or in second-QRPA . The restoration of the ISR claimed in the Ref.  actually works only for a special one-level model and can not be generalized for the realistic case as shown in Ref. . At the moment, the second-QRPA is able to reduce the violation of the ISR to the level of few percents only when contributions from the three (and more) boson states are neglected in the boson expansion of the bifermionic operators .

It can be conjectured that the closure relation is not fulfilled within the standard formulation of the RQRPA. The vacuum of the RQRPA contains quasiparticles. The RQRPA takes into account the fact that the quasiparticles block creation of the bifermionic QRPA bosons but does not consider the possibility, that the quasiparticles in the ground state take part in the transitions. This leads to violation of the unitarity and underestimates the ISR which should be fulfilled exactly.

In view of these remarks, modification of the phonon operator in order to include so-called scattering terms (describing transitions of the quasiparticles) is unavoidable if one wants to restore the ISR within RQRPA. This has been understood for some time and attempts have been undertaken to treat the terms as additional independent constituents of the phonon operator , in spirit of the thermal extensions of RPA  for description of the particle-hole excitations in heated nuclei. That leads to increasing of the dimension of the QRPA equation system and the appearance of low-lying energy roots originating from the transitions of the quasiparticles. Although such a consideration has allowed to restore the ISR for finite temperatures , the low-lying spectrum has no physical interpretation when the even-even nuclei are considered. For instance, the model would predict a number of $`2^+`$ excitations below the first quadrupole state, that have never been observed in experiments (see also the discussion in Ref. ). Authors of a recent paper  have succeeded in removing these spurious states by introducing new quasiparticles but they have again used the QBA in terms of the new quasiparticles and formulated the QRPA with respect to the HFB vacuum.

The above-mentioned examples show that a correct formulation of the RQRPA is quite a delicate task. The goal of this Letter is to introduce a version of the RQRPA for even-even nuclei which is able to overcome mentioned difficulties, i.e. it fulfills the ISR and at the same time is free from the spurious solutions.

Let us start from the general basis of the RPA, which assumes that an excited state of the nucleus in question, with the angular momentum $`J`$ and the projection $`M`$, is created by applying the phonon-operator $`Q_{JM}^{}`$ to the vacuum state $`|0_{RPA}^+`$ of the initial, even-even, nucleus with the same nucleon number $`A`$:

$$|JM=Q_{JM}^{}|0_{RPA}^+\text{with}Q_{JM}|0_{RPA}^+=0.$$
(1)
At this point one usually introduces the following ansatz for the phonon-operator $`Q_{JM}^{}`$:

$$Q_{JM}^{}=\underset{\tau \tau ^{}}{}\left[X_{(\tau \tau ^{},J)}A^{}(\tau \tau ^{},JM)Y_{(\tau \tau ^{},J)}\stackrel{~}{A}(\tau \tau ^{},JM)\right].$$
(2)
where $`X_{(\tau \tau ^{},J)}`$, $`Y_{(\tau \tau ^{},J)}`$ denotes free variational amplitudes, which are calculated by solving the RQRPA equations; $`A^{}(\tau \tau ^{},JM)`$ and $`A(\tau \tau ^{},JM)`$ ($`\tau `$ and $`\tau ^{}`$ denote single-particle quantum numbers including isospin projection) are the two quasi-particle creation and annihilation operators:

$`A^{}(\tau \tau ^{},JM)=\left[a_\tau ^{}a_\tau ^{}^{}\right]_{JM}`$ $``$ $`{\displaystyle \underset{m_\tau ,m_\tau ^{}}{}}C_{j_\tau m_\tau j_\tau ^{}m_\tau ^{}}^{JM}a_{\tau m_\tau }^{}a_{\tau ^{}m_\tau ^{}}^{}.`$ (3)
The quasiparticle creation and annihilation operators $`a_{\tau m_\tau }^+`$ and $`a_{\tau m_\tau }`$ are defined by the Bogoliubov transformation

$$\left(\begin{array}{c}a_{\tau m_\tau }^+\\ \stackrel{~}{a}_{\tau m_\tau }\end{array}\right)=\left(\begin{array}{cc}u_\tau & v_\tau \\ v_\tau & u_\tau \end{array}\right)\left(\begin{array}{c}c_{\tau m_\tau }^+\\ \stackrel{~}{c}_{\tau m_\tau }\end{array}\right),$$
(4)
where $`c_{\tau m_\tau }^+`$ ($`c_{\tau m_\tau }`$) denotes the particle creation (annihilation) operator. The tilde symbol indicates the time-reversal operation, e.g. $`\stackrel{~}{a}_{\tau m_\tau }`$ = $`(1)^{j_\tau m_\tau }a_{\tau m_\tau }`$.

However it turns out, that the most appropriate way is to analyze the phonon structure in terms of the particle creation and annihilation operators (instead of the quasiparticle ones). Let us first mention two quite general conditions which should be kept in mind while constructing the phonon operator. Namely, because the phonon operator is assumed to create a nuclear state of $`A`$ nucleons acting on the $`0^+`$ ground state of the initial nucleus (also of $`A`$ nucleons), the operator has to

* have good angular momentum $`J`$ and the projection $`M`$ (it has been taken into account in Eqs. (2),(3));
* commute with the total particle number operator $`\widehat{A}=\widehat{N}+\widehat{Z}=\underset{\tau ,m_\tau }{}c_{\tau m_\tau }^{}c_{\tau m_\tau }`$. Neutral excitations can be further distinguished from the charge-exchange ones using their commutativity with $`\widehat{N}\widehat{Z}`$.

The first rule tells us that the simplest building blocks for the phonon operator are $`\left[c_\tau ^{}\stackrel{~}{c}_\tau ^{}\right]_{JM}`$, $`\left[c_\tau ^{}c_\tau ^{}^{}\right]_{JM}`$ together with their hermitian conjugates. Applying the second rule leads to a further restriction of the phonon operator structure:

$$Q_{JM}^{}=\underset{\tau \tau ^{}}{}\left[x_{(\tau \tau ^{},J)}C^{}(\tau \tau ^{},JM)y_{(\tau \tau ^{},J)}\stackrel{~}{C}(\tau \tau ^{},JM)\right].$$
(5)
with $`C^{}(\tau \tau ^{},JM)=\left[c_\tau ^{}\stackrel{~}{c}_\tau ^{}\right]_{JM}`$ and $`\stackrel{~}{C}(\tau \tau ^{},JM)=()^{JM}C(\tau \tau ^{},JM)`$. For the neutral excitations $`\tau `$ and $`\tau ^{}`$ have the same isospin projections, whereas for the charge-exchange ones - opposite. Therefore, in terms of the original particle creation and annihilation operators $`Q`$ consists of scattering terms only. Going now into the quasiparticle representation, one gets:

$`C^{}(\tau \tau ^{},JM)`$ $`=`$ $`u_\tau v_\tau ^{}A^{}(\tau \tau ^{},JM)+u_\tau ^{}v_\tau \stackrel{~}{A}(\tau \tau ^{},JM)`$ (6)
$`+`$ $`u_\tau u_\tau ^{}B^{}(\tau \tau ^{},JM)v_\tau v_\tau ^{}\stackrel{~}{B}(\tau \tau ^{},JM)`$ (7)
$`=`$ $`u_\tau v_\tau ^{}\overline{A}^{}(\tau \tau ^{},JM)+v_\tau u_\tau ^{}\stackrel{~}{\overline{A}}(\tau \tau ^{},JM).`$ (8)
with the following definitions: $`B^{}(\tau \tau ^{},JM)=\left[a_\tau ^{}\stackrel{~}{a}_\tau ^{}\right]_{JM}`$ for the quasiparticle scattering term and $`\overline{A}^{}=A^{}+\left(u_\tau ^{}v_\tau ^{}B^{}u_\tau v_\tau \stackrel{~}{B}\right)/(v_\tau ^{}^2v_\tau ^2)`$ for the bifermionic operator $`\overline{A}^{}`$.

Now we are able to rewrite (5) in the form which is similar to (2):

$$Q_{JM}^{}=\underset{\tau \tau ^{}}{}\left[X_{(\tau \tau ^{},J)}\overline{A}^{}(\tau \tau ^{},JM)Y_{(\tau \tau ^{},J)}\stackrel{~}{\overline{A}}(\tau \tau ^{},JM)\right],$$
(9)
where $`X=u_\tau v_\tau ^{}xv_\tau u_\tau ^{}y,Y=u_\tau v_\tau ^{}yv_\tau u_\tau ^{}x`$. Thus, instead of $`A^{},A`$, bifermionic operators $`\overline{A}^{},\overline{A}`$ are now the basic terms in which the RQRPA should be formulated. The key point is that the latter automatically contain the quasiparticle scattering terms which, however, are not associated with any additional degrees of freedom. That means that there are no spurious low-lying solutions in the present theoretical scheme.

From this point we can follow the usual way to formulate the RQRPA  substituting everywhere $`A`$ by $`\overline{A}`$. The forward- and backward- going free variational amplitudes X and Y satisfy the equation:

$$\left(\begin{array}{cc}𝒜& \\ & 𝒜\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)=_{QRPA}\left(\begin{array}{cc}𝒰& 0\\ 0& 𝒰\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right),$$
(10)
where

$`𝒜`$ $`=`$ $`0_{RPA}^+|[\overline{A},[H_F,\overline{A}^{}]]|0_{RPA}^+,`$ (11)
$``$ $`=`$ $`0_{RPA}^+|[\overline{A},[H_F,\overline{A}]]|0_{RPA}^+,`$ (12)
$`H_F`$ is the nuclear Hamiltonian and the renormalization matrix $`𝒰_{\tau \tau ^{}}`$ is

$`𝒰_{\tau \tau ^{}}`$ $`=`$ $`0_{RPA}^+|[\overline{A}(\tau \tau ^{},JM),\overline{A}^{}(\sigma \sigma ^{},JM)]|0_{RPA}^+=\delta _{\tau \sigma }\delta _{\tau ^{}\sigma ^{}}𝒟_{\tau \tau ^{}},`$ (13)
$`𝒟_{\tau \tau ^{}}`$ $`=`$ $`1+\left((u_\tau ^{}^2v_\tau ^{}^2)𝒩_\tau ^{}(u_\tau ^2v_\tau ^2)𝒩_\tau \right)/(v_\tau ^{}^2v_\tau ^2)`$ (14)
Here, the exact (fermionic) expressions of the commutators are taken into account and

$$𝒩_\tau =\widehat{j}_\tau ^20_{RPA}^+|\underset{m_\tau }{}a_{\tau m_\tau }^{}a_{\tau m_\tau }|0_{RPA}^+=\widehat{j}_\tau ^10_{RPA}^+|B(\tau \tau ,00)|0_{RPA}^+$$
(15)
gives the relative quasiparticle occupation number for the level $`\tau `$ in the RQRPA vacuum ($`\widehat{j}\sqrt{2j+1}`$). The method to calculate these occupation numbers can be found elsewhere . It is useful to introduce the notation:

$$\overline{X}=𝒰^{1/2}X,\overline{Y}=𝒰^{1/2}Y,$$
(16)
$$\overline{𝒜}=𝒰^{1/2}𝒜𝒰^{1/2},\overline{}=𝒰^{1/2}𝒰^{1/2}.$$
(17)
Then the amplitudes $`\overline{X}`$ and $`\overline{Y}`$ satisfy the equation of usual QRPA:

$$\left(\begin{array}{cc}\overline{𝒜}& \overline{}\\ \overline{}& \overline{𝒜}\end{array}\right)\left(\begin{array}{c}\overline{X}\\ \overline{Y}\end{array}\right)=_{QRPA}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\overline{X}\\ \overline{Y}\end{array}\right),$$
(18)
Solving the FR-QRPA equations, one gets the fully renormalized amplitudes $`\overline{X}`$, $`\overline{Y}`$ with the usual normalization and closure relation, which will be used below to prove the fulfillment of the ISR:

$`{\displaystyle \underset{\tau \tau ^{}}{}}\overline{X}_{(\tau \tau ^{},J)}^m\overline{X}_{(\tau \tau ^{},J)}^k\overline{Y}_{(\tau \tau ^{},J)}^m\overline{Y}_{(\tau \tau ^{},J)}^k`$ $`=`$ $`\delta _{km},`$ (19)
$`{\displaystyle \underset{m}{}}\overline{X}_{(\tau \tau ^{},J)}^m\overline{X}_{(\tau _1^{}\tau _2^{},J)}^m\overline{Y}_{(\tau \tau ^{},J)}^m\overline{Y}_{(\tau _1^{}\tau _2^{},J)}^m`$ $`=`$ $`\delta _{\tau \tau _1^{}}\delta _{\tau ^{}\tau _2^{}},`$ (20)
$`{\displaystyle \underset{m}{}}\overline{X}_{(\tau \tau ^{},J)}^m\overline{Y}_{(\tau _1^{}\tau _2^{},J)}^m\overline{Y}_{(\tau \tau ^{},J)}^m\overline{X}_{(\tau _1^{}\tau _2^{},J)}^m`$ $`=`$ $`0.`$ (21)
Here, $`m`$ and $`k`$ mark different roots of the QRPA equations for a given $`J^\pi `$.

Now we are ready to prove the fulfillment of the ISR within the FR-QRPA. Let restrict ourselves to the charge-exchange pn-FR-QRPA and write down the definition for the ISR:

$$ISR=\underset{m}{}\left|JM,m|\beta _{JM}^{}|0_{RPA}^+\right|^2\underset{m^{}}{}\left|JM,m^{}|\beta _{JM}|0_{RPA}^+\right|^2.$$
(22)
Here, we would like to consider the Fermi and Gamow-Teller transitions simultaneously by defining:

$$\beta _{JM}^{}=\underset{pn,m_pm_n}{}pm_p|q_{JM}|nm_nc_p^{}c_n=\widehat{J}^1\underset{pn}{}pq_JnC^{}(pn,JM),$$
(23)
where $`q_{00}=1`$ (Fermi transitions) and $`q_{1\mu }=\sigma _\mu `$ (Gamow-Teller transitions) and one has for both cases $`ISR=NZ`$ (in the latter case - just because of the only fixed angular momentum projection). One can use Eqs. (8),(9) to get:

$`JM,m|\beta _{JM}^{}|0_{RPA}^+`$ $`=`$ $`\widehat{J}^1{\displaystyle \underset{pn}{}}pq_Jn\left(u_pv_n\overline{X}_{(pn,J)}^m+v_pu_n\overline{Y}_{(pn,J)}^m\right)\sqrt{𝒟_{pn}},`$ (24)
$`JM,m|\beta _{JM}|0_{RPA}^+`$ $`=`$ $`()^{JM}\widehat{J}^1{\displaystyle \underset{pn}{}}nq_Jp\left(u_pv_n\overline{Y}_{(pn,J)}^m+v_pu_n\overline{X}_{(pn,J)}^m\right)\sqrt{𝒟_{pn}}.`$ (25)
Substituting (25) into (22) and making use of the closure conditions (21), one ends up with

$$ISR=\widehat{J}^2\underset{pn}{}\left|pq_Jn\right|^2(v_n^2v_p^2)𝒟_{pn}.$$
(26)
Within the usual QBA one has $`𝒟_{pn}=1`$ and

$$ISR=\widehat{J}^2\left(\underset{np}{}\left|pq_Jn\right|^2v_n^2\underset{np}{}\left|pq_Jn\right|^2v_p^2\right)=NZ,$$
(27)
making use of $`\underset{p}{}\left|pq_Jn\right|^2=\widehat{J}^2\widehat{j}_n^2,\underset{n}{}\left|pq_Jn\right|^2=\widehat{J}^2\widehat{j}_p^2`$ and of the equations for the BCS chemical potentials $`\lambda _n,\lambda _p`$

$`0_{HFB}^+|\widehat{N}|0_{HFB}^+={\displaystyle \underset{n}{}}\widehat{j}_n^2v_n^2=N,0_{HFB}^+|\widehat{Z}|0_{HFB}^+={\displaystyle \underset{p}{}}\widehat{j}_p^2v_p^2=Z.`$ (28)
This corresponds to the fulfillment of the ISR for the usual QRPA.

In the case of the FR-QRPA with $`𝒟_{pn}`$ determined by Eq. (14) one has

$`ISR`$ $`=`$ $`{\displaystyle \underset{n}{}}\widehat{j}_n^2\left(v_n^2+(u_n^2v_n^2)𝒩_n\right)`$ (29)
$``$ $`{\displaystyle \underset{p}{}}\widehat{j}_p^2\left(v_p^2+(u_p^2v_p^2)𝒩_p\right)`$ (30)
$`=`$ $`NZ`$ (31)
because of the modified FR-QRPA equations for the chemical potentials:

$`0_{RPA}^+|\widehat{N}|0_{RPA}^+`$ $`=`$ $`{\displaystyle \underset{n}{}}\widehat{j}_n^2\left(v_n^2+(u_n^2v_n^2)𝒩_n\right)=N,`$ (32)
$`0_{RPA}^+|\widehat{Z}|0_{RPA}^+`$ $`=`$ $`{\displaystyle \underset{p}{}}\widehat{j}_p^2\left(v_p^2+(u_p^2v_p^2)𝒩_p\right)=Z.`$ (33)
In conclusions, we have reformulated the RQRPA into Fully-Renormalized QRPA (FR-QRPA) for even-even nuclei using restrictions imposed by the commutativity of the phonon creation operator with the total particle number operator. The FR-QRPA is free from the spurious low-energy solutions. We have also shown analytically that the Ikeda sum rule is fulfilled within the FR-QRPA.

###### Acknowledgements.

V.R. would like to thank the Graduiertenkolleg ”Hadronen im Vakuum, in Kernen und Sternen” GRK683 for supporting his stay in Tübingen.