text:
# Numerical evaluation of a two loop diagram in the cut-off regularization

## I Introduction

Dimensional regularization has become an efficient regularization scheme for theoretical calculations in high energy physics. It makes analytical results more available, and its gauge invariance almost decisively defeated the other schemes (say cutoff scheme) in QCD and other gauge theories.

However, things are not so simple. It is very natural to anticipate that DR might fail sooner or later somewhere as it is just a special regularization but never a true theory of the short distance physics. Such event is now known: the $`\overline{MS}`$ scheme does not provide a physical prediction for low energy nucleons scattering in case of long scattering lengths in nonperturbative context. To yield meaningful results, unconventional subtraction is devised and now known as PDS There are other schemes that could yield similar predictions as PDS did, see, e.g., Ref... Since this is obtained with hindsight, we should not discard the cutoff scheme without further exploiting its virtues. There have recently been some further investigations over the relations between the divergences in DR and cutoff schemes and the associated limitations. The conventional virtues of the DR scheme were shown to be due to an implicit subtraction. In nonperturbative contexts involving multiloop corrections the DR might even make renormalization infeasible.

Traditionally, it is very hard to obtain analytical results for multiloop amplitudes within the cutoff scheme. In this report we show that one could achieve this goal numerically at a pretty good precision. We will make use of the dynamical symmetry breaking of $`\lambda \varphi ^4`$ discussed in  to illustrate the numerical evaluation of the finite constant term in the two loop sunset diagram. Thus the report is organized as follows. First we quote some results of Ref. in Sec. II as preparation. In Sec. III we present our strategy and procedures for doing the numerical analysis of the problem and its possible prospective significance. The result is also given there. In Sec. IV we revisit the problem introduced in Sec. II to show the significance of the new result obtained in Sec. III. The summary will also be given there.

## II renormalization of the two loop effective potential

In, the dynamical symmetry breaking of massless $`\lambda \varphi ^4`$ model with $`Z_2`$ symmetry was studied in two loop effective potential. The algorithm in use is well known according to Jackiw,

$`V_{(2l)}\lambda \varphi ^4+{\displaystyle \frac{1}{2}}I_0(\mathrm{\Omega })+3\lambda I_1^2(\mathrm{\Omega })48\lambda ^2\varphi ^2I_2(\mathrm{\Omega }),`$ (1)
$`\mathrm{\Omega }\sqrt{12\lambda \varphi ^2};I_0(\mathrm{\Omega })={\displaystyle \frac{d^4k}{(2\pi )^4}\mathrm{ln}(1+\frac{\mathrm{\Omega }^2}{k^2})};I_1(\mathrm{\Omega })={\displaystyle \frac{d^4k}{(2\pi )^4}\frac{1}{k^2+\mathrm{\Omega }^2}};`$ (2)
$`I_2(\mathrm{\Omega })={\displaystyle \frac{d^4kd^4l}{(2\pi )^8}\frac{1}{(k^2+\mathrm{\Omega }^2)(l^2+\mathrm{\Omega }^2)((k+l)^2+\mathrm{\Omega }^2)}}.`$ (3)
Here we have Wick rotated all the loop integrals into Euclidean space. The loop integrals and the effective potential had been calculated in literature . For comparison, we list out below only the results for the sunset diagram defined by $`I_2(\mathrm{\Omega })`$ in Eq.( 3) , the only two loop diagram with overlapping divergences.

### A Loop amplitudes in DR and cutoff schemes

From  we find that

$`\mu ^{4ϵ}I_2^{(D)}(\mathrm{\Omega })={\displaystyle \frac{3\mathrm{\Omega }^2}{2(4\pi )^4ϵ^2}}\{1+(32\overline{L})ϵ+(2\overline{L}^26\overline{L}+7+6S{\displaystyle \frac{5}{3}}\zeta (2))ϵ^2\}`$ (4)
with $`S=_{n=0}^{\mathrm{}}\frac{1}{(2+3n)^2},\overline{L}=L+\gamma \mathrm{ln}4\pi `$ and $`L=\mathrm{ln}\frac{\mathrm{\Omega }^2}{\mu ^2}`$. While from ,

$`I_2^{(\mathrm{\Lambda })}(\mathrm{\Omega })={\displaystyle \frac{1}{(4\pi )^4}}\{2\mathrm{\Lambda }^2{\displaystyle \frac{3\mathrm{\Omega }^2}{2}}\mathrm{ln}^2{\displaystyle \frac{\mathrm{\Omega }^2}{\mathrm{\Lambda }^2}}+3\mathrm{\Omega }^2\mathrm{ln}{\displaystyle \frac{\mathrm{\Omega }^2}{\mathrm{\Lambda }^2}}+o(\mathrm{\Lambda }^2)\}.`$ (5)
Note that the finite double logarithmic term (leading term) disagree between the two regularization schemes, the disagreement is removed only after all subdivergences are removed.

We note in particular that the term proportional to $`\mathrm{\Omega }^2`$ in the two loop integral is not explicitly given in , and it is this finite term and its numerical estimation that is our main concern in this paper.

### B Prescription dependence

The renormalization were performed in Ref. in several prescriptions, $`\overline{MS}`$ for $`V_{(2l)}^{(D)}(\mathrm{\Omega })`$ , Jackiw, Coleman-Weinberg (CW) and $`\mu _\mathrm{\Lambda }^2`$ (defined in, Appendix B) for $`V_{(2l)}^{(\mathrm{\Lambda })}(\mathrm{\Omega })`$ with the renormalized potential taking the following form

$`V_{(2l)}(\mathrm{\Omega })=\mathrm{\Omega }^4\{{\displaystyle \frac{1}{144\lambda }}+{\displaystyle \frac{L1/2}{(8\pi )^2}}+{\displaystyle \frac{3\lambda }{(4\pi )^4}}[L^2+2(L1)^2+\alpha ]\},`$ (6)
where $`L=\overline{L}`$, $`\mathrm{ln}\frac{\mathrm{\Omega }^2}{\mu _\mathrm{\Lambda }^2}`$, $`\mathrm{ln}\frac{\mathrm{\Omega }^2}{12\lambda \mu _{Jackiw}^2}`$ and $`\mathrm{ln}\frac{\mathrm{\Omega }^2}{12\lambda \mu _{CW}^2}\frac{25}{6}`$ in respectively $`\overline{MS}`$, $`\mu _\mathrm{\Lambda }^2`$, Jackiw and Coleman-Weinberg schemes. In all the above formulas the scheme dependence of field strength and coupling constant are understood. The explicit intermediate renormalization prescriptions dependence of the effective potential expressed by $`\alpha `$ is listed in Table I.

The dynamical symmetry breaking solution in the two loop effective potential defined in Eq. (6) could be obtained from the first order condition $`\frac{V_{(2l)}(\mathrm{\Omega }(\varphi ))}{\varphi }=0`$ and the solutions read,

$`\varphi _+^2(\lambda ;[\mu ,\alpha ])={\displaystyle \frac{\mu ^2}{12\lambda }}\mathrm{exp}\left\{{\displaystyle \frac{1}{6}}\left[1{\displaystyle \frac{4\pi ^2}{3\lambda }}+\sqrt{{\displaystyle \frac{1}{3}}[436\alpha (1+{\displaystyle \frac{4\pi ^2}{\lambda }})^2]}\right]\right\},`$ (7)
whose existence requires that

$`\alpha <{\displaystyle \frac{1}{12}},\lambda \lambda _{cr}{\displaystyle \frac{4\pi ^2}{\sqrt{436\alpha }1}}.`$ (8)
That means not all schemes are consistent with symmetry breaking. In addition, for the symmetry breaking solutions to be stable, the coupling must be further constrained,

$`\lambda \widehat{\lambda }_{cr}{\displaystyle \frac{4\pi ^2}{\sqrt{436\alpha 27}1}}(>\lambda _{cr}).`$ (9)
Both Eq. (8) and (9) are summarized in Table II.

If the finite term proportional to $`\mathrm{\Omega }^2`$ in $`I_2(\mathrm{\Omega })`$ existed, the $`\alpha `$ in $`\mu _\mathrm{\Lambda }^2`$ and Jackiw prescriptions would consequently be replaced by $`2\frac{4C}{3}`$ and $`\frac{5}{4}\frac{4C}{3}`$ in the following parameterization of Eq. (5)§§§In Coleman-Weinberg prescription the definition of $`\alpha =\frac{49}{3}`$ is not altered due to its special renormalization condition. ,

$`I_2^{(\mathrm{\Lambda })}(\mathrm{\Omega })={\displaystyle \frac{1}{(4\pi )^4}}\{2\mathrm{\Lambda }^2{\displaystyle \frac{3\mathrm{\Omega }^2}{2}}\mathrm{ln}^2{\displaystyle \frac{\mathrm{\Omega }^2}{\mathrm{\Lambda }^2}}+3\mathrm{\Omega }^2\mathrm{ln}{\displaystyle \frac{\mathrm{\Omega }^2}{\mathrm{\Lambda }^2}}+C\mathrm{\Omega }^2+o(\mathrm{\Lambda }^2)\}.`$ (10)
It is not difficult to see that a different value of $`\alpha `$ means different symmetry breaking status, hence different physics, a nonperturbative scheme dependence problem as emphasized in ref.. It is therefore important to determine the value of $`C`$ in cutoff scheme in order to be sure that the Jackiw and $`\mu ^2`$ prescription are really consistent symmetry breaking in two loop effective potential. Thus the rest of this report will be devoted to estimate this constant term numerically.

## III Strategy for numerical estimation of $`c`$

First let us put the sunset diagram into the following form after integrating out the angular variables,

$`I_2(M)={\displaystyle \frac{1}{2(4\pi )^4}}{\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}{\displaystyle _0^\beta }{\displaystyle \frac{dy}{y+M}}\left\{x+y+M\sqrt{(x+y+M)^24xy}\right\}`$ (11)
with $`M\mathrm{\Omega }^2,xk^2,yl^2`$ and $`\beta \mathrm{\Lambda }^2`$. The asymptotic form of this integral as cutoff tends to infinity reads

$`I_2(M)`$ $`=`$ $`{\displaystyle \frac{1}{(4\pi )^4}}\left\{2\beta {\displaystyle \frac{3M}{2}}\mathrm{ln}^2{\displaystyle \frac{M}{\beta }}+3M\mathrm{ln}{\displaystyle \frac{M}{\beta }}+CM\right\}+o(M/\beta )`$ (12)
$``$ $`I_2^{(asy)}(M)+{\displaystyle \frac{CM}{(4\pi )^4}}+o(M/\beta ).`$ (13)
In order to determine the constant $`C`$ we compare $`I_2(M)`$ and $`I_2^{(asy)}(M)`$ after we evaluated the former numerically.

The feasibility of numerical treatment of an apparently divergent integral is guaranteed by the observation that the cutoff need not be very large. Eq.( 12) tells us that the analytical source of error comes from $`\frac{M}{\beta }`$. The magnitude of the constant $`C`$ should be of order 1 or 10 in the limit $`\mathrm{\Lambda }^2\mathrm{}`$. Thus taking $`\frac{M}{\beta }`$ to be $`10^4`$ or $`10^5`$ should be sufficient for our purpose. There is a main obstacle in the practice for numerical estimation of divergent integral in the cutoff scheme: the presence of power law pieces, which consume most of the computer’s capacity in order to get a finite number at a precision of $`10^4`$. The larger the cutoff is chosen, the less precision capacity is left over for finite numbers. We must get rid of such pieces. The remaining logarithmic terms save most of the computer capacity for the finite numbers, as even taking the cutoff as large as $`10^{10}M`$ will only yield a double log term at order $`100\mathrm{ln}10230`$, which consumes almost no precision capacity. This goal could be readily achieved at least in two ways: for the sunset diagram, the power law term $`2\beta `$ could be removed either by subtracting it explicitly in the integrand or by differentiating the whole integral with respect to $`M`$. In fact we will report both approaches’ result.

There is another subtle point that need careful consideration, the step length for the numerical integration. Taking the step too small will improve the precision analytically but cost a lot of time and might accumulate numerical errors, especially in personal computer. Usually a divergent integral has a smooth integrand except in a finite range where the integrand is steeply changing. The step must be small enough within such steep regions. The last concern is to spare both the computer time and numerical error we should perform the analytical calculation as far as we can. Thus we will first perform one integral analytically (this is feasible in most cases) while leave the remaining one to computer.

In the following we will describe one approach of our numerical treatment in detail, the one with differentiation with respect to $`M`$. That is, to extract $`C`$ we just consider the following double integral which at most possesses logarithmic divergence,

$`{\displaystyle \frac{I_2(M)}{M}}={\displaystyle \frac{1}{(4\pi )^4}}\left\{I_1I_2+{\displaystyle \frac{1}{2}}(I_3I_4)\right\}={\displaystyle \frac{I_2^{(asy)}(M)}{M}}+{\displaystyle \frac{C}{(4\pi )^4}}+o(M/\beta ),`$ (14)
$`I_1{\displaystyle _0^\beta }{\displaystyle \frac{dx}{(x+M)^2}}{\displaystyle _0^\beta }{\displaystyle \frac{dy}{y+M}}(x+y+M),`$ (15)
$`I_2{\displaystyle _0^\beta }{\displaystyle \frac{dx}{(x+M)^2}}{\displaystyle _0^\beta }{\displaystyle \frac{dy}{y+M}}\sqrt{(x+y+M)^24xy},`$ (16)
$`I_3{\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}{\displaystyle _0^\beta }{\displaystyle \frac{dy}{y+M}},`$ (17)
$`I_4{\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}{\displaystyle _0^\beta }{\displaystyle \frac{dy}{y+M}}{\displaystyle \frac{x+y+M}{\sqrt{(x+y+M)^24xy}}}.`$ (18)
It is not difficult to see that the $`2\beta `$ piece is removed in all the four integrals $`I_1,I_2,I_3`$ and $`I_4`$ there is at most double log terms. Then it is also easy to perform two folds integrals once. Moreover, $`I_1`$ and $`I_3`$ can be done analytically at all. Then the results read

$`I_1=\mathrm{ln}{\displaystyle \frac{\beta +M}{M}}\left\{\mathrm{ln}{\displaystyle \frac{\beta +M}{M}}+{\displaystyle \frac{M}{\beta +M}}1\right\}+{\displaystyle \frac{\beta ^2}{M(M+\beta )}},`$ (19)
$`I_2={\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}\left\{{\displaystyle \frac{\sqrt{R}}{Y}}+{\displaystyle \frac{x}{\sqrt{a}}}\mathrm{sinh}^1{\displaystyle \frac{axY}{2Y\sqrt{xM}}}+\mathrm{sinh}^1{\displaystyle \frac{Yx}{2\sqrt{xM}}}\right\}_{Y=M}^{Y=\beta +M},`$ (20)
$`I_3=\mathrm{ln}^2{\displaystyle \frac{\beta +M}{M}},`$ (21)
$`I_4={\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}\left\{{\displaystyle \frac{x}{\sqrt{a}}}\mathrm{sinh}^1{\displaystyle \frac{axY}{2Y\sqrt{xM}}}+\mathrm{sinh}^1{\displaystyle \frac{Yx}{2\sqrt{xM}}}\right\}_{Y=M}^{Y=\beta +M}`$ (22)
with $`Ra2xY+Y^2,ax^2+4xM`$. There is a seemingly singularity in the integrands ($`\frac{1}{\sqrt{xM}}`$ as $`x0`$) in both $`I_2`$ and $`I_4`$ which does not materialize after one completes the operations on the variable $`Y`$,

$`I_2={\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}\left\{{\displaystyle \frac{\sqrt{R}}{Y}}_{Y=M}^{Y=\beta +M}+{\displaystyle \frac{\sqrt{x}}{\sqrt{x+4M}}}\mathrm{ln}{\displaystyle \frac{u+\sqrt{u^2+1}}{d+\sqrt{d^2+1}}}+\mathrm{ln}{\displaystyle \frac{s+\sqrt{s^2+1}}{t+\sqrt{t^2+1}}}\right\},`$ (23)
$`I_4={\displaystyle _0^\beta }{\displaystyle \frac{dx}{x+M}}\left\{{\displaystyle \frac{\sqrt{x}}{\sqrt{x+4M}}}\mathrm{ln}{\displaystyle \frac{u+\sqrt{u^2+1}}{d+\sqrt{d^2+1}}}+\mathrm{ln}{\displaystyle \frac{s+\sqrt{s^2+1}}{t+\sqrt{t^2+1}}}\right\},`$ (24)
$`u{\displaystyle \frac{x(x+3M\beta )}{2(\beta +M)\sqrt{xM}}},v{\displaystyle \frac{x^2+3xM}{2M\sqrt{xM}}},s{\displaystyle \frac{\beta +Mx}{2\sqrt{xM}}},t{\displaystyle \frac{Mx}{2\sqrt{xM}}}.`$ (25)
The only necessary numerical works to be done are the two one fold integrals defined in Eq.( 2324). As the integrand become steeper as $`x`$ is closer to the origin, thus to reduce the error from numerical integration we separate the integral into several intervals in order to match between the step size and the steepness of the integrand. Taking $`M`$ to be 1 will make the numerical integration easy and the $`C`$ has been estimated with two cutoff scale, $`\beta /M\mathrm{\Lambda }^2/\mathrm{\Omega }^2=10^4,10^5`$. The match between the interval and step size is shown in Table III and IV respectively for $`\beta /M=10^4`$ and $`\beta /M=10^5`$.

Subtracting $`\frac{I_2^{(asy)}(M)}{M}`$ from the numerically obtained integral defined in Eq.( 14) we find the constant $`C`$: $`C=4.13048`$ for $`\beta /M=10^4`$ and $`C=4.15412`$ for $`\beta /M=10^5`$. The interval $`(0,10^6]`$ is not included as its contribution could at most be of order $`10^4`$. Noting that one cutoff is order greater than the other, the agreement of the two cases is striking. Moreover, if one try another approach, i.e., by directly subtracting $`I_2`$’s integrand a term that upon integration will yield the $`2\beta `$ piece, one could also obtain an estimate of $`C`$ with the results read: for $`\beta /M=10^4`$, $`C=4.14128`$ with a homogeneous step $`10^4`$, $`C=4.14274`$ with a homogeneous step $`10^3`$, $`C=4.156178`$ with a homogeneous step $`10^2`$; for $`\beta /M=10^5`$, $`C=4.160739`$ with a homogeneous step $`10^2`$, $`C=4.1633`$ with a homogeneous step $`0.005`$; for $`\beta /M=10^6`$, $`C=4.162993`$ with a homogeneous step $`10^2`$. Combining the two approaches’ results we could safely conclude that $`C`$ is not zero. To be more conservative we anticipate that the true value should lie in the interval $`[4,4.3]`$.

Before ending this section, we should point out and stress that in order to evaluate divergent multiloop integrals numerically in cutoff regularization: (i) one does not need a very large cutoff scale, (ii) the power law term can be easily removed in order to save capacity for more important task (like determination of the finite parts). The numerical workload could further be reduced if one varies the step size with steepness of the integrand.

## IV New predictions and summary

With the new constant we could reevaluate the critical couplings for symmetry breaking with the result summarized in Table V and VI. We take $`C`$ to be 4 when calculating the critical couplings. It is not difficult to see that the critical couplings in $`\mu _\mathrm{\Lambda }^2`$ and Jackiw prescriptions are dramatically lowered, or the symmetry breaking might take place at much lower values of the coupling. The interesting thing is, the $`\overline{MS}`$ critical coupling, which is smaller than the ones in the prescriptions based on old calculation of the sunset diagram , now becomes a larger one than the obtained with the new numerically determined constant $`C`$, and the larger the $`C`$ is, the smaller the critical couplings are. In this sense, the cutoff regularization is preferable to DR, just like the case in the EFT applications in nucleon interactions .

Generally one would expect the critical coupling to be close to 1 instead of being very large (like 10.7) for symmetry breaking to take place. Therefore it is a pleasing finding that constant $`C`$ is larger than 0 and is no less than 4. If $`C=20`$, a reasonable value, the critical coupling could be 1.5238 and 1.5542. Of course, as the potential is only calculated at two loop level, the precise value of the critical coupling could not be taken very serious. However, one should still be pleased to see that the critical couplings are dramatically lowered after more careful evaluation in the regularization scheme that is widely held as inferior to DR, even though the computation is not exact one.

In summary, we reevaluated the sunset diagram in $`\lambda \varphi ^4`$ in the cutoff regularization half analytically and half numerically and found that the finite local constant, which is usually taken to be zero, is not zero and this constant could dramatically lower the critical couplings for dynamical symmetry breaking in the two loop effective potential in a number of renormalization prescriptions. The important byproduct is that the numerical calculation could be efficiently done in cutoff regularization for divergent multiloop diagrams and we have illustrated that the cutoff scale need not be too large.

## Acknowledgement

Ji-Feng Yang thanks W. Zhu for his helpful discussions. This work is supported in part by the National Natural Science Foundation of China under Grant No.s 10075020 and 10205004.


format: text

text:
# Epidemics, disorder, and percolation

## 1 Introduction

Models for the spread of epidemics are an interesting application of non-equilibrium statistical physics. In this paper we discuss two aspects of this sort of study in cases when disorder is present. First we generalize the well-known  application of percolation theory to epidemics for which variation in the strength of coupling. Second, we consider the contact networks which carry epidemics. If there is large variation in the number of contacts between members of the population, then the network structure can lead to a percolation threshold of zero, for example if the distribution of contacts is a suitable power law. We show an example of a network which seems to have $`p_c0`$ for the same power-law distribution, and which might have some relevance to real infections.

## 2 SIR epidemics in an inhomogeneous population

The SIR (Susceptible, Infected, Recovered) model is a classic way to think about the spead of epidemic infections . It is based on the idea that in many epidemics an individual catches a disease, infects others at rate $`x`$ for a certain period, $`\tau `$, and then recovers or dies, and is thus removed from the infection network. A mean-field version of this model is given by the SIR equations for a perfectly mixed population:

$$dS/dt=xSI/N,dI/dt=xSI/NI/\tau $$
(1)
Here N is the size of the population, and we normalize $`x`$ in this way so that we can think about the infection rate per contact. Clearly, we have an outbreak if $`R_o=xS(0)\tau /N>1`$.

We can think of Eq. (1) as describing dynamics on a fully conjugated graph – every node sees every other equally. In this paper we will consider a model on a lattice, where each individual only sees its neighbors. Each site of the lattice is initially occupied by a susceptible ($`S`$) individual. One lattice site is initially infected, $`I`$. Any $`S`$ can be infected by $`I`$’s which are nearest neighbors. The probability of infection per unit time along the $`i^{th}`$ bond, $`x_i`$, is chosen from a distribution $`f(x)`$. For a case near the threshold for the epidemic (see below) the pattern of infection is shown in Figure 2. The similarity of this pattern to a percolation cluster is very clear.

The existence of a large epidemic depends on the parameters just as in the mean-field case. For fixed $`x_i`$, if the recovery time, $`\tau `$, is too small the epidemic will die out, and if $`\tau `$ is large enough it will persist. A critical value, $`\tau _c`$ marks the threshold. However, even for large $`\tau `$ there is a finite probability that the epidemic dies. We define a spanning probability, $`M(\tau )`$, the probability that an infection started from a single site in the middle of a large lattice reaches the edge. $`M`$ grows rapidly for $`\tau >\tau _c`$. In Figure 2 we show $`M`$ for a 200 x 200 lattice as a function of $`\tau `$ for two different choices of $`f(x)`$, one rectangular, and one strongly asymmetric.

This model is not obviously related to percolation theory. For a very narrow distribution, $`f(x)=\delta (xx_o)`$ Grassberger  showed that we do have percolation as follows: the probability to infect a neighbor is $`p(\tau )=1(1x_o)^\tau `$. Thus $`p`$ is the fraction of bonds completed in a given epidemic, and the infection process is mapped onto bond percolation with $`p(\tau )`$ playing the role of the percolation probability. For a square lattice in two dimensions, $`p_c=1/2`$, so we expect that when $`p`$ is greater than $`1/2`$ we will have an epidemic. This gives a critical recovery time: $`\tau _c=\mathrm{ln}2/\mathrm{ln}(1x_o)`$. For any other network, $`p_c=1(1x_o)^{\tau _c}`$ determines the threshold for epidemic spread. For arbitrary $`f`$’s, we can write:

$$p(\tau )=1(1x)^\tau f(x)𝑑x=p_c.$$
(2)
Thus, for any $`f`$ we have mapped our infection problem to percolation. The solution of Eq. (2) gives $`\tau _c`$. The mapping of a dynamic problem onto static percolation is not trivial: different runs for infection spreading in a system with a given realization of $`x_i`$ give different, but not independent realizations of a percolation problem, i.e. of a set of a completed bonds. In previous work  we have verified Eq. (2) by doing a data collapse for various different choices of $`f(x)`$. After data collapse, the function $`M`$, expressed as a function of $`p`$ using Eq. (2) is identical to the mass of the infinite cluster in ordinary percolation theory. All of the results of this work may be applied easily for any network.

In order to understand the effect of the distribution of bond strengths we studied $`f`$’s which are either weakly or strongly asymmetric. An example of a symmetric $`f`$ is:

$$f_w(x)=1/x_{max},0xx_{max}.$$
(3)
In cases like this where $`f(x)`$ is concentrated near its mean value, it is easy to show the spread of $`x_i`$ does influence the infection propagation very much.

However, the dependence $`p(\tau )`$, Eq.(2), is strongly nonlinear and is dominated by the behavior of $`f(x)`$ in vicinity of $`x=0`$: the spread of the infection is controlled by weakest bonds. Thus strongly asymmetric distributions concentrated near zero should be very different from the case discussed above. Consider, for example, power-laws $`f(x)=(\alpha 1)x^\alpha `$ with $`0<\alpha <1`$. In this case $`(1x)^\tau f(x)𝑑x=(\alpha 1)B(\tau +1,1\alpha )`$, where $`B(x,y)`$ is the beta-function. Strongly asymmetric distributions correspond to values of $`\alpha `$ in vicinity of $`1`$. Using the fact that for such distributions $`\tau 1`$, we can use the Stirling formula to get in leading order: $`\tau _c\left[\mathrm{\Gamma }(2\alpha )(1p_c)\right]^{1/(1\alpha )}`$ which diverges for $`\alpha 1`$.

An even more extreme case is that of functions $`f(x)`$ with $`\alpha 1`$. These are not acceptable probability densities since they are not normalizable on . However, one can truncate the $`f`$ in the vicinity of zero to allow normalization. For example:

$$f_s(x)=C/xx_{min}x1;C=1/|\mathrm{ln}(x_{min})|.$$
(4)
The two distributions considered in Figure 2 are given by Eq. (4) with $`x_{min}=e^{15}`$, an extremely asymmetric case, and Eq.(3), the symmetric case.

In Figure 2 we see that the effect of strong asymmetry is to spread out the critical region, i.e. the region where $`M`$ is not close to either 1 or 0. This is the parameter region where the epidemic pattern is patchy, as in Figure 2. For a strongly asymmetric $`f`$ the patchy pattern corresponds to a large range of the parameter $`\tau `$. We recall from percolation theory that the regime where $`M`$ is neither very small nor very close to 1 is where the fractal nature of the percolation cluster at $`p_c`$ persists over large length scales. This could have implications for practical epidemiology: if the actual distribution of $`x_i`$ is very asymmetric, we would be very likely to observe patchiness. For such distributions mean-field theory is very inaccurate  and sampling of populations must be done with care because of the spatial correlations.

## 3 The network of contacts: power-law distributions

In the previous section we considered networks for which the number of contacts between individuals does not vary much, though we supposed that the strength of such contacts varied a good deal. However, there is data to suggest  that some infection networks have a wide distribution of the number of contacts, so that the bond number, $`k`$, is drawn from a distribution of the form $`P(k)1/k^\alpha `$ with $`2<\alpha <3`$. In this case the mean value of the number of contacts, $`<k>`$, exists, but the fluctuation, $`<k^2>`$ diverges.

We have seen that it is sufficient to consider percolation on such networks to understand SIR epidemics. This problem has been discussed a good deal . There is general agreement that $`p_c`$ is zero for many cases. The graphs considered by these authors are essentially random graphs with a given degree distribution. Thus the model has a kind of random mixing.

The random graphs dominate the the ensemble of all graphs with a given $`P(k)`$. The study of this ensemble was introduced by Molloy and Reed  and a criterion for the formation of a giant component was derived. Callaway, et. al  and Cohen, et al.  express a related result for percolation: $`p_c(<k^2>/<k>^22)=1`$ Thus if $`<k^2>`$ diverges we expect $`p_c=0`$. We can see this by starting at a random node on the graph so that there are $`k_0`$ bonds of which $`pk_0`$ are active. If neighbors are chosen at random, the probability for the number of bonds on the next neighbor is $`p^2k_1P(k_1)/<k>`$. The factor of $`k_1`$ comes from the fact that highly bonded sites are more likely to have a bond with the first site. Now the number of infected increases by a factor $`<k_1>/<k_0>=p<k^2>/<k>^2`$ at each such step. If $`<k^2>`$ is finite, then this can be made less than unity for $`p`$ small. and the infection will die out. If $`<k^2>`$ is infinite, then $`p`$ must be 0.

Now consider a model which has the same $`P(k)`$ but has *geography* . We assume that individuals infect their neighbors, but that the number of infections is drawn from $`P`$. We construct this by taking each point on a $`d`$-dimensional lattice with probability $`p`$ and connecting them all other points within a $`d`$-sphere of volume $`k`$, centered on the point. In 2 dimensions we surround each point with a disk of radius of roughly $`R=[k/\pi ]^{1/2}`$. The probability distribution of $`R`$ is given by $`P(R)=P(k)dk/dRR^{(d\alpha d+1)}`$. Our model is a variant of an old two-dimensional circle model of continuum percolation  with variable radius of circles. See Figure 3. This is a small part of the Molloy-Reed ensemble, but we suggest that it may be relevant for real-world epidemics: sometimes infection is primarily local.

Note that a site can have the bonds it generates itself (‘proper’ bonds) and those from other disks. We show that the total number of bonds per node has the same power-law behavior as the proper bonds if $`\alpha >2`$ by calculating the mean number of bonds connecting a node $`i`$ to nodes $`j`$ whose center is a distance $`r>R_i`$ away. The probability to find $`R_j`$ larger than $`r`$ is $`_r^{\mathrm{}}P(r^{})𝑑r^{}`$. Thus, the mean number of nodes connected to $`i`$ from outside is $`k_{out}_{R_i}^{\mathrm{}}𝑑rr^{d1}_r^{\mathrm{}}P(R_j)𝑑R_j.`$ Thus $`k_{out}R_i^{d(\alpha 2)}`$ which tends to zero for large $`R`$, as long as $`\alpha >2`$. This means that the probability distribution of the number of the ‘proper’ bonds and the actual number of the bonds of a node is the same for $`k`$ large, and that bonds from outside are rare.

For $`d=1`$ it is easy to see that the system never percolates, because for any $`p<1`$ there will eventually be a gap of length $`\mathrm{\Delta }`$ in the chain of connected disks. Consider disks with radii between $`k`$ and $`k+dk`$; the concentration of their centers is $`pP(k)dk`$. The probability that the gap is not covered by one of these disks is equal to the probability that none of their centers are found in the interval of length $`2k+\mathrm{\Delta }`$ (in units of the lattice constant) centered at the gap, $`1pP(k)(2k+\mathrm{\Delta })dk\mathrm{exp}[pP(k)(2k+\mathrm{\Delta })dk]`$ The overall probability that no disk overlaps the gap is the product over $`k`$:

$$\mathrm{exp}\left(_0^{\mathrm{}}pP(k)(2k+\mathrm{\Delta })𝑑k\right)=\mathrm{exp}\left(2pkp\mathrm{\Delta }\right),$$
(5)
which is finite as long as $`k`$ exists.

In 2 dimensions we did simulations to find $`p_c`$. Disks were randomly added to an $`L`$ x $`L`$ lattice at different sites until a cluster spanned the lattice. Figure 4 shows the average percolation threshold $`p_c(L)`$ as a function of $`1/L`$. It appears from the plot that there is a finite percolation threshold for not just $`\alpha >3`$ but $`2<\alpha <3`$ as well, the region of interest for real world epidemics. For $`\alpha <2`$, the results are consistent with $`p_c=0`$ and seem to scale according as $`L^{2(\alpha 2)}`$ for sufficiently large $`L`$. We have shown  that this scaling is related to the presence of giant disks that span the whole sample. These occur with finite probability for $`\alpha <2`$.

We can try to understand the results in Figure 4 by calculating $`<k_1>`$ as above. Most of the bonds emanating from the center of disk point within it. These have degrees drawn from the distribution $`P(k_1)`$, not $`k_1P(k_1)`$. Thus, for the typical behavior we $`<k_1>=<k>`$ for $`\alpha >2`$. However, there are rare events which correspond to the overlap of large disks from outside with the center of the starting disk. These give a large contribution. In fact, in 2 dimensions for a disk of size $`R_0`$ we add up the contributions of outside disks of radius $`R_1`$:

$$<k_1>_{R_0}^{\mathrm{}}r𝑑r_r^{\mathrm{}}R_1^2P(R_1)𝑑R_1=_{R_0}^{\mathrm{}}𝑑rr^{(2\alpha 5)}$$
(6)
This diverges for $`\alpha <3`$. Note that this is an overestimate since have not attempted to exclude already infected sites.

We have a confusing situation: the typical behavior is that the cluster of disks is finite in extent, but the average mass added diverges. There is a rigorous theorem in the study of continuum percolation which treats exactly this situation . It states that if $`<k^{21/d}>`$ is finite, but $`<k^2>`$ diverges ( in 2 dimensions $`2.5<\alpha <3`$), then there is a $`p_c`$ such that for $`p<p_c`$ *with probability unity all clusters are finite, but the expected number of points of the largest cluster is infinite.* (The probability to span a large lattice is not addressed by this theorem.) We interpret this odd situation, at least for $`\alpha >2.5`$ as a form of intermittancy. Typical epidemics are finite for small $`p`$, but huge, rare events occur. It is interesting to ask the physical question of what we would expect to observe in this case. We claim (perhaps too hopefully) that the situation is best represented by the computer simulations. Figure 4 shows that at $`\alpha =2.5`$ there is quite solid evidence for a non-zero $`p_c`$, and this seems to extend down to $`\alpha =2`$. (Note that the theorem does not exclude this possibility).

In summary we have a counterexample to the prevailing notion that $`p_c`$ must be zero if $`<k^2>`$ is infinite, certainly in one dimension. In 2 dimensions the situation is complex, and will take more work to fully sort out.

We have had many useful conversations with Mark Newman and Carl Simon.


text:
# The PDF of fluid particle acceleration in turbulent flow with underlying normal distribution of velocity fluctuations

## I Introduction

In the context of Tsallis formalism Tsallis  and generalized statistics approach Johal ; Aringazin ; Beck3  to developed turbulence, C. Beck Beck4  used the underlying log-normal distribution, $`f(\beta )`$, to describe fluctuation of the parameter $`\beta `$ entering some Langevin equation for the acceleration $`a`$ of fluid particle in the Lagrangian frame; see also Beck2 ; Wilk . For the case of a linear drift force, the resulting marginal probability density function of the acceleration is Beck4

$$P(a)=\frac{1}{2\pi s}_0^{\mathrm{}}𝑑\beta \beta ^{1/2}\mathrm{exp}\left[\frac{(\mathrm{ln}\frac{\beta }{m})^2}{2s^2}\right]e^{\frac{1}{2}\beta a^2},$$
(1)
where $`m=\mathrm{exp}[s^2]`$ provides a unit variance and $`s`$ is a fitting parameter, $`s^2=3.0`$. This distribution was found to be in a very good agreement with the recent Lagrangian experimental data by Porta, Voth, Crawford, Alexander, and Bodenschatz Bodenschatz , the new data by Crawford, Mordant, and Bodenschatz (the Taylor microscale Reynolds number is $`R_\lambda =690`$, the normalized acceleration range is $`[60,60]a`$, the Kolmogorov timescale $`\tau _\eta `$ is resolved) Bodenschatz2 , Mordant, Delour, Leveque, Arneodo, and Pinton ($`R_\lambda =740`$, $`a[20,20]`$, $`\tau _\eta `$ is not resolved) Mordant , and direct numerical simulations of the Navier-Stokes equations by Kraichnan and Gotoh ($`R_\lambda =380`$, $`a[150,150]`$Gotoh .

The same approach with the underlying $`\chi ^2`$ distribution, $`f(\beta )`$, of fluctuations resulted in an analytically explicit distribution Aringazin3 ,

$$P(a)=\frac{C\mathrm{exp}[a^2/a_c^2]}{(1+\frac{1}{2}V_0(q1)a^2)^{1/(q1)}},$$
(2)
where $`C`$ is a normalization constant, $`V_0=4`$, and $`q=3/2`$ (Tsallis entropic index) are due to the theory, and $`a_c`$ is a free parameter used for a fitting, $`a_c=36.0`$. This distribution was found to be in a good agreement with the experiments Bodenschatz2 . The result (2) is based on the consideration of Ref. Beck .

In general, the marginal distribution is defined as

$$P(a)=_0^{\mathrm{}}𝑑\beta P(a|\beta )f(\beta ),$$
(3)
where $`P(a|\beta )`$ is a probability density function associated to a surrogate dynamical equation, the Langevin equation for the acceleration $`a`$ Beck ,

$$_ta=\gamma F(a)+\sigma L(t).$$
(4)
Here, $`F(a)`$ is a ”drift force”, $`L(t)`$ is delta-correlated Langevin source (Gaussian white noise, $`L(t)=0`$, $`L(t)L(t^{})=2\delta (tt^{})`$), $`\beta =\gamma /\sigma ^2`$ is assumed to be a fluctuating real positive parameter. For constant parameters $`\gamma `$ and $`\sigma `$, this model assumes that the stochastic process (4) is Markovian, and $`P(a|\beta )`$ is found as a stationary solution of the associated Fokker-Planck equation. For a linear drift force, $`F(a)=a`$, the stationary conditional distribution, $`P(a|\beta )`$, is found to be

$$P(a|\beta )=C(\beta )\mathrm{exp}[\beta a^2/2],$$
(5)
where $`C(\beta )`$ is a normalization constant, $`a[\mathrm{},\mathrm{}]`$.

The function $`f(\beta )`$ entering Eq. (3) is probability density function arising from the assumption that the parameter $`\beta `$ is a stochastic variable. For a constant $`\beta `$, the Gaussian solution (5) meets the large scale (large time increment) statistics but it fails to describe observed Reynolds number dependent stretched exponential tails of the experimental small scale $`P(a)`$ Bodenschatz  that correspond to high probability to find extremely big values of the fluid particle acceleration in the developed turbulent flow.

The interest in studying Langevin type equations to describe developed turbulence is motivated by the recent high precision Lagrangian experiments Bodenschatz ; Bodenschatz2 ; Mordant , which give an important dynamical information and new look to the intermittency in fluid turbulence. Time response characteristics of the tracer polystyrene particle and the precision in the experiments Bodenschatz ; Bodenschatz2  allow to resolve about 1/20 of the Kolmogorov time and 1/40 of the Kolmogorov length ($`R_\lambda =970`$) so that the acceleration is really resolved and particle follows rare violent events (within 7% of the ideal value of acceleration even at the high Reynolds numbers studied there), with the collected statistical data sufficient to establish convergence of the fourth moment of acceleration.

In Sec. 2, we construct a simple model that allows one to derive the marginal probability density function, $`P(a)`$, of the acceleration of fluid particle in turbulent flow, with underlying normally distributed velocity fluctuations. In Sec. 3, we reproduce the result of Ref. Beck4  as a particular case, and propose a generalization of the log-normal model. In Sec. 4, we briefly summarize the results and make a few remarks.

## II The marginal distribution with underlying normally distributed velocity fluctuations

While the conditional distribution $`P(a|\beta )`$ can be found starting from the Langevin equation (4), we do not have sufficiently strong theoretical requirements to determine a unique form of the distribution $`f(\beta )`$ except for that it should obey general conditions, such as that the integral in Eq. (3) should converge and the resulting marginal distribution $`P(a)`$ should be normalizable. To select a particular form of the distribution $`f(\beta )`$, one can use ad hoc statistical distributions or make assumptions stemming from considerations of the developed turbulence.

As a first step, one can associate the parameter $`\beta `$ with the mean energy dissipation rate $`ϵ`$, and use the relationship of $`ϵ`$ with the fluid velocity $`v_i`$, within the framework of Kolmogorov scaling theory of developed turbulence.

Second, as shown by Renner, Peinke, and Friedrich Friedrich , the small scale intermittency can be traced back to a stochastic nature of the averaged energy dissipation rate Landau , and the Markovian condition for the stochastic process is fulfilled at scales larger than the Taylor microscale $`\lambda `$. Thus, a stochastic nature of $`\beta `$ can be related to the stochastic energy dissipation rate (measured in Eulerian experiments), which we denote by $`\stackrel{~}{ϵ}`$.

For example, the assumptions $`\beta \stackrel{~}{ϵ}`$ and $`\stackrel{~}{ϵ}u_i^2`$, discussed at length by Beck Beck  (we do not repeat it here for brevity), for the fluctuating, averaged over a ball of radius $`r`$, energy dissipation $`\stackrel{~}{ϵ}`$ and velocity fluctuations $`u_i`$ at the Kolmogorov scale were used to propose $`\beta u_i^2`$ and to select the $`\chi ^2`$ distribution of $`\beta `$. Three-dimensionality of the space was used explicitly to determine the value of the free parameter (Tsallis entropic index) in this model. A comparison of the formalism based on Eq. (4) with the Sawford model Sawford  implies $`\beta =2C_0^2a_0\nu ^{1/2}ϵ^{3/2}`$ Beck4  so that with the replacement of $`ϵ`$ by the fluctuating $`\stackrel{~}{ϵ}`$ and using again $`\stackrel{~}{ϵ}u_i^2`$, one can obtain a different relation, $`\beta u_i^3`$. The Kolmogorov’s 1941 approach implies $`ϵ\overline{u}^3`$, where $`\overline{u}`$ stands for the rms velocity, which can also be used to try to relate $`\stackrel{~}{ϵ}`$ to the velocity fluctuations. Although successful in capturing main features of the experimental data, different models use different powers of $`ϵ`$ to represent $`\beta `$ that makes a theoretical problem when selecting an appropriate distribution of $`\beta `$.

It is however common to the above representations that velocity fluctuations, $`u_i`$, are assumed to be normally distributed with zero mean, and that statistics of $`\beta `$, $`\stackrel{~}{ϵ}`$, and $`u_i`$ are related to each other due to some functional dependencies between their characteristics holding in the inertial range. A direct dependence of statistical properties of the acceleration $`a`$ on velocity fluctuations was established in the wellknown Heisenberg-Yaglom theory. This gives a different look to the problem since the Lagrangian description refers to individual trajectories of fluid particles as compared with the well established Eulerian framework, in which the intermittency is understood in terms of anomalous scaling of the moments of velocity increments in space related to the stochastic nature of the energy dissipation rate (non-dynamical description). From a simplified dynamical point of view, one can think of that the tracer particle trajectory extends to a large region (few integral length scales, in the experiments) crossing during the course subregions characterized by different local amplitudes of velocity fluctuations which are randomly distributed in space and varies with time (Brownian like motion of the particle driven by the stochastic delta-correlated force with stochastically slow varying intensity). The Lagrangian velocity autocorrelation function is known to decay very slowly, to vanish at times bigger than the integral time scale. This view is valid for sufficiently high Reynolds numbers and meets that provided by the stochastic energy dissipation rate argued earlier Beck . This is also in an agreement with the remark made in Ref. Bodenschatz  that extremely high accelerations seem to be associated with coherent structures (geometrically these may be thought of as very thin tornadoes) which persist for many Kolmogorov times, substantially longer than the correlation time of the acceleration components. At lower Reynolds numbers, $`R_\lambda <500`$, the acceleration is increasingly coupled to large scales of the flow in a dynamical way so that a different approach may be required to describe such a situation to a good accuracy.

In the Lagrangian frame, it is then natural to relate parameters of the stochastic dynamical equation (4), which describes the acceleration in the context of a generalized Brownian motion, to velocity fluctuations due to the Heisenberg-Yaglom picture of developed turbulence. Ultimately of course we deal with a statistical description of the acceleration and the velocity fluctuations within the framework of Fokker-Planck equation associated to Eq. (4). Guided by this observation, one can develop a general consideration along this line of reasoning.

In the present paper, we formulate some general requirements for $`f(\beta )`$ using the assumption that the velocity fluctuations $`u_i`$ are Gaussian distributed random variables with zero means, $`u_i=0`$, and variances $`s_i`$,

$$g(u_i)=C(s_i)\mathrm{exp}\left[u_i^2/(2s_i^2)\right],$$
(6)
where $`C(s_i)`$ is normalization constant; $`u_i[\mathrm{},\mathrm{}]`$, $`i=1,2,3`$. This assumption meets the Lagrangian experimental data Bodenschatz , for each component of the velocity fluctuations.

It should be stressed that our treatment is made regardless to a particular functional dependence of $`\beta `$ on $`\stackrel{~}{ϵ}`$, or a particular form of $`f(\beta )`$, and provides a simple way to formalize the general model in accord to the Heisenberg-Yaglom picture, which relates some statistical properties of the acceleration to that of the velocity fluctuations.

We assume that $`\beta `$ depends on the velocity $`v_i`$ the dynamics of which is decoupled from that of the acceleration so that the well known fluctuating character of $`v_i`$ for the turbulent flow determines fluctuations of $`\beta `$. Below we drop the index $`i`$ to simplify notation. Also, we treat the developed turbulent flow to be statistically isotropic and translation invariant thus discarding skewness effects; the statistical anisotropy is known to become very small with the increase of the Reynolds number, e.g., at $`R_\lambda >900`$.

In the mean field approximation, we can represent $`v=v+u`$, where the fluctuation $`u`$ is characterized by zero mean, $`u=0`$. Hence, we can write

$$\beta =\beta _c(v)+B(v,u),$$
(7)
where we formally separated out $`\beta _c`$ which depends only on $`v`$.

We suppose that the function $`B(u)B(v,u)`$ satisfies the following natural requirements: (i) It is a sufficiently smooth invertible function with respect to $`u`$ so that the inverse function $`u(\beta \beta _c)`$ can be found; (ii) it maps the interval $`[\mathrm{},\mathrm{}]`$ to $`[\beta _c,\mathrm{}]`$, to provide positiveness of $`\beta `$; also one can restrict consideration by the requirement that (iii) there is one-to-one correspondence between $`\beta `$ and $`u`$.

In general, for $`\beta `$ to be a stochastic variable the function $`B(u)`$ should be Borel function of the stochastic variable $`u`$. We assume that $`u`$ has an absolutely continuous distribution function and the associated probability density function, in particular, that given by Eq. (6). Not any Borel function $`B(u)`$ provides existence of a probability density function for $`\beta `$. However, a probability density function for $`\beta `$ a priori exists if $`B(u)`$ is a monotonic function.

Using the general relation

$$g(u)du=g(u(\beta \beta _c))\left|\frac{du}{d\beta }\right|d\beta $$
(8)
and Eq. (6), we can make the identification $`f(\beta )=g(u(\beta \beta _c))|du/d\beta |`$, i.e.,

$$f(\beta )=C(s)\mathrm{exp}\left[\frac{[u(\beta \beta _c)]^2}{2s^2}\right]\left|\frac{du}{d\beta }\right|,$$
(9)
where $`C(s)`$ is a normalization constant. The inverse function $`u(\beta \beta _c)`$ should provide positive definiteness and normalizability of the above function $`f(\beta )`$. From Eq. (3) we then finally have

$$P(a)=C(s)_0^{\mathrm{}}𝑑\beta P(a|\beta )\mathrm{exp}[\frac{[u(\beta \beta _c)]^2}{2s^2}]\left|\frac{du}{d\beta }\right|.$$
(10)
This equation allows one to calculate the distribution $`P(a)`$ for given stationary solution $`P(a|\beta )`$ associated to the Langevin equation (4) and the invertible (monotonic) function $`B(u)`$. By construction, the structure of this model is such that the velocity fluctuation $`u`$, underlying the dynamics, is normally distributed with zero mean and variance $`s`$. Note that only the amplitude of $`u`$ contributes $`P(a)`$. Eq. (10) can be viewed as a specific case of the general equation (3) which captures one of the well known features of the turbulent flow, and thus enables to rule out some ad hoc types of distribution of $`\beta `$.

Note that in Eq. (10) we did not change the variable over which the integral is performed as compared to Eq. (3). The main idea was to introduce the variable $`u`$ which follows Gaussian distribution and encode various underlying distributions in the functional dependence $`\beta =B(u)`$. This significantly reduces the class of admissible probability density functions $`f(\beta )`$ that we consider as one of the main results of our approach.

In the next Section, we study a specific choice of $`B(u)`$ which is relevant both from the mathematical and physical points of view.

## III The log-normal model and beyond

As mentioned by Beck Beck4 , selection of the log-normal distribution $`f(\beta )`$ provides that for any power law, $`\beta \stackrel{~}{ϵ}^\kappa `$, the function $`\mathrm{ln}\beta =\kappa \mathrm{ln}\stackrel{~}{ϵ}`$ guarantees that the functional form of the log-normal distribution does not change, which is viewed as a hint for a physical relevance of a log-normally distributed $`\beta `$. Also, it is known that the probability density function of the averaged energy dissipation rate is log-normal in agreement with Kolmogorov’s assumption in K62 theory.

The exponential dependence,

$$B(u)=e^{u/u_0},$$
(11)
where $`u_0`$ is constant, and $`\beta _c=0`$, provides a relevant example for which we get

$$u=u_0\mathrm{ln}\beta ,\left|\frac{du}{d\beta }\right|=\frac{|u_0|}{\beta }.$$
(12)
The distribution (9) then becomes the log-normal distribution, and Eq. (10) for $`P(a|\beta )`$ given by Eq. (5) reproduces the marginal distribution (1), with $`u_0=1`$.

Also, supposing $`\beta \stackrel{~}{ϵ}^\kappa `$, we obtain the relation

$$\stackrel{~}{ϵ}e^{u/\kappa },$$
(13)
where we put $`u_0=1`$ for simplicity. Again, one can see that the only measurable effect of the parameter $`\kappa `$ is that it scales the variance, $`ss/\kappa `$, in Eq. (10).

One can make the following refinement of the dynamical equation (4) using the relation (11), which corresponds to the log-normal model (1). The defining feature of the considered model is the assumption of slow fluctuating character of the composite parameter $`\beta =\gamma /\sigma ^2`$. It can be easily shown that a noisy character of the drift coefficient $`\gamma `$, with $`F(a)=a`$ and delta-correlated Gaussian white noise $`\gamma `$ with non-zero mean, can be used to describe intermittency effects (power law tails) because it is related to wellknown random multiplicative processes extensively studied in the literature. Below, we focus on a stochastic nature of the additive noise amplitude $`\sigma `$ and study its contribution to the intermittency separately. A joint effect of both the noisy $`\gamma `$ and random intensity of the additive noise $`\sigma `$ is out of scope of the present paper and can be studied elsewhere.

From Eq. (11) it follows $`\sigma ^2=e^{u/u_0}`$ so that the Langevin equation (4) becomes

$$_ta=\gamma F(a)+\mathrm{exp}\left[u/(2u_0)\right]L(t),$$
(14)
where the velocity fluctuation $`u`$ is treated statistically and follows Gaussian distribution with zero mean, $`F(a)=a`$, and $`\gamma `$ is taken to be constant to simplify consideration of the contribution of the additive noise intensity. The Langevin model (14) meets the Lagrangian experiments to a very good accuracy as it implies the probability density function of the form (1).

Equation (14) assumes that the dynamics of Lagrangian acceleration $`a`$ (the acceleration in the comoving frame, $`a_i_tv_i+v_k_kv_i`$), which comes mainly from small scales, and that of the fluctuating energy dissipation rate $`\stackrel{~}{ϵ}`$, or the related amplitude of velocity fluctuations, $`u_i=v_iv_i`$, which is associated mainly to large scales, are weakly related to each other at high Reynolds numbers, $`R_\lambda >500`$, so that in the lowest approximation it is natural to take them independent. At small time scale, the velocity fluctuation is characterized by approximately constant amplitude while its directional part changes wildly (cf. Bodenschatz ). However, in the statistical context, which reflects Lagrangian dynamics at all the time scales, the acceleration is related to the amplitude of velocity fluctuation as it follows, for example, from the Heisenberg-Yaglom scaling, $`a^2a_0\overline{u}^{9/2}\nu ^{1/2}L^{3/2}`$, where $`\overline{u}`$ is the rms velocity and $`L`$ is the integral scale length. This longstanding universal $`\overline{u}^{9/2}`$ scaling was confirmed by the recent Lagrangian experiments Bodenschatz  to a very high accuracy, for $`R_\lambda >500`$. The large scale origin of the additive noise stochastic intensity, slowly varying in time, is precisely the reason of the stretching of exponential tails of the acceleration probability density function $`P(a)`$. The observed Reynolds number dependence of the stretching can be thus qualitatively understood as the result of increasing coupling of the acceleration dynamics to large scale dynamics at lower Reynolds numbers.

In summary, the dynamical equation (14) reflects the above point of view by incorporating exponential of the dynamically independent $`u`$, with the prescribed external statistics, as the intensity of the additive noise. A full dynamical treatment of Eq. (14) with $`u`$ viewed as a random function of time can be made elsewhere. Here, we note only that Langevin type equation with the additive noise of the form $`e^{\omega (t)}L(t)`$, with $`\omega (t)`$ a longtime correlated stochastic process, was recently proposed Mordant  to describe the Lagrangian intermittency in the context of new multifractal random walk model by Bacry and Muzy, a continuous extension of discrete cascade models. Our approach (14) meets this result. The difference is that we approximate the random process $`\omega `$ to be trivially stationary and it obeys normal statistics.

In accord to Eqs. (5) and (14), the associated stationary distribution for $`F(a)=a`$ is given by

$$P(a|u)=C(u)\mathrm{exp}[e^{u/u_0}a^2/2],$$
(15)
where $`C(u)`$ is a normalization constant. From this point of view, small deviations of the theoretical $`P(a)`$ given by Eq. (1) from the experimental data Bodenschatz2  (see Aringazin3  for details) could be attributed to small deviations from Gaussian distribution of the velocity fluctuation $`u`$; the flatness was measured to be approximately 2.8 for the axial and 3.2 for the longitudinal velocity component Bodenschatz  compared with the flatness value 3 for an exact Gaussian distribution. Eq. (15) can be used to account for such deviations in a phenomenological way,

$$P(a)=_{\mathrm{}}^{\mathrm{}}𝑑uC(u)\mathrm{exp}[e^{u/u_0}a^2/2]g_1(u),$$
(16)
where $`g_1(u)`$ is the distribution which is approximately Gaussian one. However, it should be noted that small departure from the experimental data can also arise from shortcomings of the considered one-dimensional Langevin model of developed turbulence.

While it is evident that the three-dimensional Navier-Stokes equation with a Gaussian random forcing belongs to a class of non-linear stochastic equations for which one can associate some generalized Fokker-Planck equations, it is a theoretical challenge to make a link between the Navier-Stokes equation and such a surrogate one-dimensional Langevin model, which is, of course, far from being a full model of the Lagrangian dynamics of fluid in the turbulent regime. This problem is addressed in the forthcoming paper Aringazin5 . Here, we note only that the Lagrangian description simplifies the problem of finding the stationary probability measure since the advection term, which dominates in the developed turbulent flow, is incorporated to the definition of Lagrangian acceleration. Also, large scale dynamics and an interaction between the large scale and small scale dynamics within the inertial range are important to provide understanding of the origin of noises (multiplicative and additive ones) used in such Langevin models and a weak dependence of the noise statistics on the small scale dynamics at high Reynolds numbers.

## IV Conclusions

(i) The presented approach implies a specific form of the probability density function $`f(\beta )`$, given by Eq. (9), which stems from the assumption that velocity fluctuations are normally distributed with zero mean, Eq. (5), as it is confirmed by the experimental data Bodenschatz ; Bodenschatz2  for high-Reynolds-number turbulent flows, to a good accuracy.

(ii) Given the exponential dependence (11) the presented formalism leads to a consideration of the log-normal distribution $`f(\beta )`$, which was proved to be relevant from the experimental point of view Beck4 . Also, such a distribution is in agreement with the log-normal distribution of mean energy dissipation in Kolmogorov 1962 picture. An exponential relation of the type (11) requires a physical treatment in the context of Langevin, or Fokker-Planck, equation which can be made elsewhere. Here, we note that the exponential dependence on the velocity fluctuation, $`\mathrm{exp}[u]`$, indicates a specific variance of the absolute value of velocity increment during the time-evolution, and can be thought of as a kind of Lyapunov instability in the velocity space.

(iii) One may also be interested in using some other distributions, instead of the normal one (6), to derive $`f(\beta )`$ and the associated theoretical distribution $`P(a)`$.

(iv) Finally, one can try other appropriate functions $`B(u)`$, for instance,

$$B(u)=e^{u/u_0}c_ku^k,$$
(17)
where $`c_k`$ are constants, instead of (11), to reproduce the experimental data on the acceleration statistics of fluid particle in turbulent flow. In particular, for $`\beta _c=0`$ and the function $`B(u)=u^2`$ the map doubly covers the interval $`[0,\mathrm{}]\beta `$, with each covering being a monotonic map, and we obtain $`u=\pm \sqrt{\beta }`$ and $`|du/d\beta |=\beta ^{1/2}`$. The basic equation (9) yields $`f(\beta )\beta ^{1/2}\mathrm{exp}[\beta /(2s^2)]`$, which is $`\chi ^2`$ probability density function of order one.

Acknowledgement. The authors are grateful to referee for valuable comments, which allowed to improve physical argumentation in the revised version of the paper.


text:
# Double-well potentials in current qubits

## I Introduction

The rf-SQUID has attracted great attention in connection with the study of quantum two-level systems Leggett . In particular it has offered prospects for macroscopic quantum phenomena, such as quantum tunneling and quantum coherence between two quantum states Caldeira . A frequently appearing structure is the well-known double-well potential with a smooth barrier in some parameter regime. Surprisingly, however, the inspection into the rf-SQUID with precise boundary conditions has not been achieved.

A periodic boundary condition can be obtained by a gauge invariance and the phase evolution by the wave vector of Cooper pairs which is identical with the fluxoid quantization condition. Using this boundary condition, we can obtain the effective potential of an rf-SQUID and three-Josephson junction loop which consist of the energy of moving Cooper pairs and the Josephson coupling energy. The external flux induces a persistent current of Cooper pairs in superconducting loop. When the self inductance of the superconducting loop and the Josephson coupling energy are large, the total flux, sum of the external flux and induced flux, become quantized, while the fluxoid quantization condition is always satisfied.

Considering the energy of Cooper pairs, the effective potential of the system shows a new periodic structure consisting of cusp and smooth barriers. In one extreme, the tunneling through smooth-barrier potential is dominant, which is the usual case Leggett ; Caldeira ; Kurk ; Friedman ; Mooij ; Wal , while, in the other extreme, the cusp barrier potential is more appropriate in tunneling phenomena. Both may be used as a qubit in quantum computing Nielsen ; Makhlin . We compare the tunneling rates for the usual smooth-barrier potential and the cusp-barrier potential and show that the latter may offer more effective qubit.

## II rf-SQUID

We first consider an rf-SQUID shown in Fig. 1 (a). We derive the effective potential of the system showing double-well shapes with both smooth and cusp barriers. The current through the Josephson junction in Fig. 1 (a) is written as

$`I=I_c\mathrm{sin}\varphi +{\displaystyle \frac{V}{R}}+C\dot{V},`$ (1)
where $`\varphi `$ is the superconducting phase difference across the Josephson junction and $`I_c`$, $`R`$, and $`C`$ are critical current, resistance, and capacitance, respectively. Using the relation, $`I=(\mathrm{\Phi }_t\mathrm{\Phi }_{\mathrm{ext}})/L_s`$, and the Josephson voltage-phase relation, $`V=(\mathrm{\Phi }_0/2\pi )\dot{\varphi }`$, Eq. (1) becomes

$$\mathrm{\Phi }_t=\mathrm{\Phi }_{\mathrm{ext}}L_sI_c\mathrm{sin}\varphi L_s\frac{\mathrm{\Phi }_0}{2\pi }\left(\frac{\dot{\varphi }}{R}+C\ddot{\varphi }\right).$$
(2)
Here $`\mathrm{\Phi }_t`$ is the total flux, $`\mathrm{\Phi }_{\mathrm{ext}}`$ the external flux, $`\mathrm{\Phi }_{\mathrm{ind}}=L_sI`$ the induced flux, $`L_s`$ the self inductance of the superconducting loop, and $`\mathrm{\Phi }_0=h/2e`$ the unit flux quantum. When $`L_sI_c/\mathrm{\Phi }_0<<1`$, one can use the approximation $`\varphi 2\pi \mathrm{\Phi }_t/\mathrm{\Phi }_0`$ in Eq. (2), which results in the usual smooth double-well potential Kurk .

The fluxoid quantization condition for a thin superconducting loop Tinkham ,

$$\mathrm{\Phi }_t+(m_c/q_c)\stackrel{}{v}_c𝑑\stackrel{}{l}=n\mathrm{\Phi }_0$$
(3)
with $`q_c=2e`$, $`m_c=2m_e`$, and the average velocity of Cooper pairs $`\stackrel{}{v}_c`$, can be represented as the periodic boundary condition,

$$k_nl=2\pi n+2\pi f_t,$$
(4)
determining the wave vector at a quantum state $`n`$, where $`k_n`$ is the wave vector of the Cooper pairs, $`l`$ the circumference of the loop, and $`f_t\mathrm{\Phi }_t/\mathrm{\Phi }_0`$. Here, $`k_nl`$ is the phase accumulated along the circumference of the loop Matveev . For a rf-SQUID in Fig. 1 (a) the boundary condition with the phase difference across the Josephson junction $`\varphi `$ becomes

$`k_nl=2\pi n+2\pi f_t\varphi .`$ (5)
The current, on the other hand, is given by $`I=(n_cAq_c/m_c)\mathrm{}k_n`$, where $`n_c`$ is the Cooper pair density, $`A`$ the cross section, and the induced flux $`\mathrm{\Phi }_{\mathrm{ind}}=L_sI`$ becomes

$`{\displaystyle \frac{\mathrm{\Phi }_{\mathrm{ind}}}{\mathrm{\Phi }_0}}={\displaystyle \frac{\gamma l}{2\pi }}k_n,`$ (6)
where $`\gamma (q_c^2n_cA/lm_c)L_s`$. Combining the expression for the total flux,

$`f_t=f+{\displaystyle \frac{L_sI}{\mathrm{\Phi }_0}},`$ (7)
with Eq. (5), the wave vector is written as

$`(1+\gamma )k_nl=2\pi n+2\pi f\varphi ,`$ (8)
where $`f\mathrm{\Phi }_{\mathrm{ext}}/\mathrm{\Phi }_0`$.

Neglecting the dissipative term originating from the quasi-particle tunneling through the junction, Eq. (1) becomes

$`(n_cAq_c/m_c)\mathrm{}k_n=I_c\mathrm{sin}\varphi C(\mathrm{\Phi }_0/2\pi )\ddot{\varphi }.`$ (9)
With the representation for $`k_n`$ in Eq. (8) we get

$`\left({\displaystyle \frac{\mathrm{\Phi }_0}{2\pi }}\right)^2C\ddot{\varphi }={\displaystyle \frac{ϵ_0}{(1+\gamma )\pi }}\left(n+f{\displaystyle \frac{\varphi }{2\pi }}\right)E_J\mathrm{sin}\varphi ,`$ (10)
where $`E_J=I_c\mathrm{\Phi }_0/2\pi `$ is the Josephson coupling energy and $`ϵ_02\pi ^2\mathrm{}^2n_cA/m_cl`$. Eq. (10) describes the motion of a particle with kinetic energy $`E_C=Q^2/2C`$ with $`Q=C(\mathrm{\Phi }_0/2\pi )\dot{\varphi }`$ in an effective potential $`U_{\mathrm{eff}}(\varphi )`$,

$`U_{\mathrm{eff}}(\varphi )=E_n+E_J(1\mathrm{cos}\varphi ),`$ (11)
where

$`E_n`$ $`=`$ $`{\displaystyle \frac{ϵ_0}{1+\gamma }}\left(n+f{\displaystyle \frac{\varphi }{2\pi }}\right)^2.`$ (12)
The meaning of $`E_n`$ can be understood more clearly when we represent $`E_n`$ as $`E_n=(1+\gamma )E_{\mathrm{kin}}`$ using Eq. (8), where $`E_{\mathrm{kin}}=(\mathrm{}^2k_n^2/2m_c)n_cAl`$ is the kinetic energy of Cooper pairs. Since the induced energy of the loop can be represented as $`E_{\mathrm{ind}}=(1/2)L_sI^2=\gamma E_{\mathrm{kin}}`$, we can see that $`E_n`$ is the sum of the kinetic energy and the induced energy Imry ; Bloch , i.e., $`E_n=E_{\mathrm{kin}}+E_{\mathrm{ind}}`$. If we introduce the kinetic inductance $`L_km_cl/An_cq_c^2`$, the kinetic energy becomes $`E_{\mathrm{kin}}=(1/2)L_kI^2`$ and the parameter $`ϵ_0`$ in Eq. (10) is represented by $`ϵ_0=\mathrm{\Phi }_0^2/2L_k`$. If we introduce the energy scale $`E_K\mathrm{\Phi }_0^2/2L_k=ϵ_0`$ and $`E_L\mathrm{\Phi }_0^2/2L_s`$, the parameter $`\gamma `$ becomes

$`\gamma ={\displaystyle \frac{L_s}{L_k}}={\displaystyle \frac{E_K}{E_L}}`$ (13)
and the energy $`E_n`$ in Eq. (12) can be rewritten as

$`E_n`$ $`=`$ $`{\displaystyle \frac{E_K}{1+\gamma }}\left(n+f{\displaystyle \frac{\varphi }{2\pi }}\right)^2.`$ (14)
The similar expression was obtained for a superconducting loop without Josephson junction Majer .

The number of excess Cooper pair charges on Josephson junction $`\widehat{N_c}\widehat{Q}/q_c`$ is conjugate to the phase difference $`\widehat{\varphi }`$ and the commutation relation $`[\widehat{\varphi },\widehat{N_c}]=i`$ gives the quantum phase fluctuations of the junction. The Hamiltonian $`\widehat{H}=\widehat{P}^2/2M+U_{\mathrm{eff}}(\widehat{\varphi })`$ describes the dynamics of the phase difference, where $`\widehat{P}\widehat{N_c}\mathrm{}=i\mathrm{}/\widehat{\varphi }`$ and $`M=(\mathrm{\Phi }_0/2\pi )^2C`$. The charging energy $`E_C`$ causes quantum tunneling between two wells in the potential $`U_{\mathrm{eff}}(\widehat{\varphi })`$ of the rf-SQUID and the following three-Josephson junction loop.

The double-well structures shown in Fig. 2 are given by taking the lowest energy for the quantum states $`n=1,0,`$ and $`1`$ for the effective potential $`U_{\mathrm{eff}}(\varphi )`$. The shape of the effective potential is determined by the parameter $`\eta `$ defined by

$`\eta {\displaystyle \frac{E_K}{(1+\gamma )E_J}}.`$ (15)
For large $`\gamma `$, the value of $`\eta `$ can be determined by two parameters $`L_s`$ and $`E_J`$. We calculate the value of the parameter $`\gamma `$ to obtain $`\gamma 600`$ for the parameters of the rf-SQUID of Ref. . At much low temperature T $``$ O(10mK), we roughly estimate the Cooper pair density in a pure metal $`n_cn/2`$ Tinkham  with the conduction electron density $`n`$. Even if we consider the effect reducing the Cooper pair density such as band structure of the materials, finite temperature, and so on, we expect that the value of $`\gamma `$ is still large such that $`\gamma >>1`$. In this case the parameter $`\eta `$ is approximately represented as follows,

$`\eta {\displaystyle \frac{E_L}{E_J}}={\displaystyle \frac{\mathrm{\Phi }_0^2}{2L_sE_J}}={\displaystyle \frac{\pi \mathrm{\Phi }_0}{L_sI_c}},`$ (16)
where $`E_J=I_c\mathrm{\Phi }_0/2\pi `$ with critical current of the Josephson junction $`I_c`$. We obtain that $`\eta 8.5`$ for the rf-SQUID Friedman .

For the value of $`\eta =10`$, we can see the double-well structure in Fig. 2 (a) with a smooth barrier centered at $`\varphi /2\pi =0.5`$ for $`\mathrm{\Phi }_{\mathrm{ext}}/\mathrm{\Phi }_0`$=0.49. The usual double-well structure with a smooth barrier corresponds to this one. The left (right) well corresponds to the diamagnetic (paramagnetic) current state. Thus the total flux, $`\mathrm{\Phi }_t=\mathrm{\Phi }_{\mathrm{ext}}+\mathrm{\Phi }_{\mathrm{ind}}`$, for the state at left (right) well approaches $`0(\mathrm{\Phi }_0)`$ as the self inductance becomes large. For a smaller $`\eta `$, however, a new double-well potential emerges around the cusp as shown in Fig. 2 (b). The high cusp barrier shrinks while the low smooth barrier grows up as $`\eta `$ decreases. Thus we get two wells separated by a cusp barrier. The state of current at right (left) well is diamagnetic (paramagnetic).

At the local minima of $`U_{\mathrm{eff}}(\varphi )`$, the condition $`U_{\mathrm{eff}}(\varphi )/\varphi =0`$ gives $`(n_cA\mathrm{}^2/m_c)k_n+E_J\mathrm{sin}\varphi =0`$ and the current relation becomes

$`{\displaystyle \frac{n_cAq_c}{m_c}}\mathrm{}k_n+I_c\mathrm{sin}\varphi =0.`$ (17)
Actually Eq. (7) is also a current relation, i.e., $`I=(\mathrm{\Phi }_t\mathrm{\Phi }_{\mathrm{ext}})/L_s`$. If one consider both $`\mathrm{\Phi }_t`$ and $`I`$ as independent dynamic variables while $`\mathrm{\Phi }_{\mathrm{ext}}`$ is a control variable and minimize the energy of the system with respect to the variable $`\mathrm{\Phi }_t`$, the current relation, Eq. (7), can be obtained Imry ; Bloch . Now the remaining dynamic variables are $`\varphi `$ and $`I`$ or equivalently $`k_n`$. The current relation, Eq. (17), is obtained by minimizing the total energy with respect to the variable $`\varphi `$. This current relation warrants the validity of our effective potential.

## III Three-Josephson junction loop

Recently the superconducting loop with three-Josephson junctions of Fig. 1 (b) has been proposed as an efficient qubit Mooij . The periodic boundary condition similar to Eq. (8) is written as

$`(1+\gamma )k_nl=2\pi n+2\pi f\varphi _1\varphi _2\varphi _3.`$ (18)
The effective potential $`U_{\mathrm{eff}}(\varphi _1,\varphi _2,\varphi _3)`$ for this qubit is given by

$`U_{\mathrm{eff}}(\varphi _1,\varphi _2,\varphi _3)`$ $`=`$ $`E_n+E_J(1\mathrm{cos}\varphi _1)+E_J(1\mathrm{cos}\varphi _2)`$ (19)
$`+`$ $`E_{J3}(1\mathrm{cos}\varphi _3),`$
where

$`E_n`$ $`=`$ $`{\displaystyle \frac{E_K}{1+\gamma }}\left(n+f{\displaystyle \frac{\varphi _1+\varphi _2+\varphi _3}{2\pi }}\right)^2.`$ (20)
Here we introduced the rotated coordinates $`(\varphi _a,\varphi _b,\varphi _c)`$ such that $`\varphi _a(1/2)(\varphi _1+\varphi _2),\varphi _b(1/2)(\varphi _1\varphi _2),`$ and $`\varphi _c(1/3)(\varphi _1+\varphi _2+\varphi _3)`$. Then the effective potential of this qubit is written as,

$`U_{\mathrm{eff}}={\displaystyle \frac{E_K}{1+\gamma }}\left(n+f{\displaystyle \frac{3}{2\pi }}\varphi _c\right)^2`$
$`+2E_J(1\mathrm{cos}\varphi _a\mathrm{cos}\varphi _b)+E_{J3}\left[1\mathrm{cos}(3\varphi _c2\varphi _a)\right].`$
Fig. 3 shows the effective potential at $`\varphi _c/2\pi =f/3`$, where we can see that the value of $`\varphi _b`$ is strongly localized around $`\varphi _b=n\pi `$ since the barrier between $`Q`$ and $`R`$ is much higher than that between $`P`$ and $`Q`$. This means that the phase differences at the potential minima in identical Josephson junctions are $`2\pi n`$. For $`\varphi _b=0`$ the effective potential is shown in Fig. 4. The barriers between $`A`$ and $`C`$ as well as between $`B`$ and $`D`$ are smooth, as can be seen in Fig. 4 (c). The two current states at $`A`$ and $`C`$ correspond to the qubit states of Ref. , where $`\eta `$ is so large that the cusp barrier is high enough for tunneling across the barrier to be neglected. The cusp barrier between $`A`$ and $`D`$ shrinks while the smooth barrier grows as $`\eta `$ decreases. In Fig. 4 (d), thus, we can see the cusp barrier between $`A`$ and $`D`$. The tunneling through this cusp barrier becomes dominant as $`\eta `$ becomes small.

The current in the loop is given by $`I=I_c\mathrm{sin}\varphi _a`$. Thus the current states at $`A`$ and $`B`$ ($`C`$ and $`D`$) are in a paramagnetic (diamagnetic) current state. The values of $`\varphi _a`$ and $`\varphi _c`$ at potential minima can be obtained through the current relations, $`U_{\mathrm{eff}}/\varphi _i=0`$ which give

$`{\displaystyle \frac{E_K}{(1+\gamma )\pi }}\left(n+f{\displaystyle \frac{\varphi _1+\varphi _2+\varphi _3}{2\pi }}\right)=E_{Ji}\mathrm{sin}\varphi _i,`$ (22)
where $`E_{J1}=E_{J2}=E_J`$. In the rotated coordinate we obtain the relations

$`{\displaystyle \frac{\eta }{\pi }}\left[n+f{\displaystyle \frac{3}{2\pi }}\varphi _c\right]`$ $`=`$ $`\mathrm{sin}\varphi _a`$ (23)
$`=`$ $`\lambda \mathrm{sin}(3\varphi _c2\varphi _a),`$
where $`\lambda E_{J3}/E_J`$ and we set $`\varphi _b=0`$. From Eqs. (18) and (23), we calculate $`\stackrel{~}{f}(1+\gamma )k_nl/2\pi `$. We calculate the values of the parameter $`\gamma `$ to obtain $`\gamma 200`$ and $`\eta 3300`$ for the three-Josephson junction qubit of Ref. . For $`\gamma >>1`$, the induced flux in Eq. (6) can be represented such that $`\mathrm{\Phi }_{\mathrm{ind}}/\mathrm{\Phi }_0=\gamma k_nl/2\pi \stackrel{~}{f}`$.

In Fig. 5 the induced flux is shown when $`f=0.45`$ and $`\lambda =0.8`$. The ground states in $`C`$ and $`D`$ induce a diamagnetic current while the ground states in $`A`$ and $`B`$ induce a paramagnetic current. In the inset of Fig. 5 we can see that $`\stackrel{~}{f}`$ is very small for the qubit of Ref.  so that the induced flux is negligible and the boundary condition, Eq. (18), becomes $`2\pi n+2\pi f=(\varphi _1+\varphi _2+\varphi _3)`$ which is used in Ref. . On the contrary, when the value of $`\eta `$ is much smaller, the induced flux become larger. The value of total flux approaches 0 ($`\mathrm{\Phi }_0`$) for diamagnetic (paramagnetic) current, which clearly shows the flux quantization condition. Thus we can see that the total flux becomes quantized as $`L_sE_J`$ increases, while the fluxoid quantization condition, Eq. (3), is always satisfied. For superconducting ring without Josephson junction, the total flux is quantized when $`L_s`$ is sufficiently large Imry .

## IV Tunneling frequency

The quantum states in a cusp potential may be used to make a more efficient qubit because of high tunneling frequency across the barrier. This allows more single qubit operations in a short dephasing time. The geometric configuration of the cusp potential, such as its higher curvature and smaller width, causes higher tunneling frequency. We compare the tunneling rate calculated by the WKB approximation, the instanton method Weiss , and the Fourier Grid Hamiltonian Method(FGHM) FGHM  in Fig. 6. The instanton method calculates the tunneling rate $`t`$ as follows,

$`t=\mathrm{}\omega _0C_0\sqrt{{\displaystyle \frac{S_{\mathrm{inst}}}{2\pi \mathrm{}}}}\mathrm{exp}\left[{\displaystyle \frac{S_{\mathrm{inst}}}{\mathrm{}}}\right],`$ (24)
where $`S_{\mathrm{inst}}=M_{\mathrm{}}^{\mathrm{}}𝑑\tau \dot{\varphi }^2(\tau )`$, $`M=(\mathrm{\Phi }_0/2\pi )^2C=\mathrm{}^2/8E_c`$, $`\omega _0=\sqrt{k/M}`$ is the frequency of small oscillation in a well with $`k=^2U_{\mathrm{eff}}(\varphi )/\varphi ^2|_{\varphi _{\mathrm{min}}}`$, and $`C_0`$ is a factor which depends on the shape of the double-well potential. The instanton action is evaluated in such a way that $`S_{\mathrm{inst}}=\sqrt{2M}_{\varphi _{\mathrm{min}}}^{\varphi _{\mathrm{min}}}\sqrt{U_{\mathrm{eff}}(\varphi )}𝑑\varphi `$, where $`\varphi _{\mathrm{min}}`$ and $`\varphi _{\mathrm{min}}`$ represent two potential minima. The instanton satisfies the equation of motion $`(1/2)M\dot{\varphi }^2(\tau )=U_{\mathrm{eff}}[\varphi (\tau )]`$. Integrating the equation of motion to obtain the instanton trajectory and comparing it with the asymptotic behavior, $`C_0`$ is given by

$`C_0=\varphi _{\mathrm{min}}\sqrt{{\displaystyle \frac{2M\omega _0}{S_{\mathrm{inst}}}}}\mathrm{exp}[\omega _0I],`$ (25)
where $`I`$ is a constant obtained numerically.

The values of $`E_C`$ are chosen through the parameter $`\alpha E_J/E_C`$ which is roughly proportional to $`1/\eta `$ in experiments Mooij ; Wal ; Friedman . In Fig. 6 we see that the frequency of tunneling through the cusp barrier for low $`\eta `$ is much higher than that for the smooth barrier in both rf-SQUID and the three-Josephson junction qubit. The tunneling through the potential barriers in the effective potential, $`U_{\mathrm{eff}}`$, introduces the tunneling term in the single qubit Hamiltonian,

$`H_{\mathrm{qubit}}=F(\mathrm{\Phi }_{\mathrm{ext}})\sigma ^zt\sigma ^x,`$ (26)
where $`F(\mathrm{\Phi }_{\mathrm{ext}})`$ is the energy level separation between two current states and $`\sigma ^{z(x)}`$ are the Pauli matrices. The first (second) term in Eq. (26) causes the single qubit rotation, $`R_z(R_x)`$. The operation time $`\tau _{\mathrm{op}}`$ for the $`R_x`$ gate is inversely proportional to the tunneling frequency. Since the quality factor is $`Q_\varphi =\tau _\varphi /\tau _{\mathrm{op}}`$ with the dephasing time $`\tau _\varphi `$, the present results show that the cusped barrier allows the quality factor that is an order of magnitude larger than occur with the smooth barrier.

Lower $`\eta `$ requires higher $`E_J`$ as shown in Eq. (16). Thus, a higher inductive flux in Fig. 5 could be detrimental to the phase coherence of the Cooper pair wave function. However, recent experiments for a three-Josephson junction loop with a Cooper-pair island Vion  and a current-biased Josephson junction qubit Martinis  show that using a large Josephson junction the high quality factor can be achieved. A current-biased qubit of dc SQUID type with large Josephson junctions can be described by the present theory for low value of $`\eta `$.

The quantum phase slip may be relevant in a thin superconducting loop when the energy levels are degenerate Matveev . The phase slip process causes the change of tunneling frequency and, thus, induces the decoherence. For the lower value of $`\eta `$ in Eq. (16) the loop inductance is larger, which means the loop is thicker. Therefore the qubit using the cusp barrier with the lower value of $`\eta `$ as shown in Fig. 2 may have a longer dephasing time than that using the smooth one.

## V Summary

Using the boundary condition for the phase evolution of the wave function of the Cooper pairs, we obtain a new double-well potential with cusp barriers for both rf-SQUID and three-Josephson junction loop. At the minimum of the effective potential we show that the current relation is satisfied. This picture allow us to calculate the wave vector of Cooper pairs, $`k_n`$, as a function of $`\eta `$ and, thus, obtain the induced flux to show that the flux quantization condition is satisfied when the value of $`\eta `$ becomes small. Moreover we calculate the tunneling frequency as a function of $`\eta `$ which is the characteristic parameter responsible for the shape of the potential. The tunneling frequency in the cusp potential is much higher than that in the smooth potential, thus, the cusp potential offers a more efficient qubit.

ACKNOWLEDGMENTS

This work is supported by the Korean Science and Engineering Foundation through the Center for Strongly Correlated Materials Research (2002) at Seoul National University and the Brain Korea 21 Project in 2002.


text:
# 1 Abstract

## 1 Abstract

A quantum statistical system with energy dissipation is studied. Its statisitics is governed by random complex-valued non-Hermitean Hamiltonians belonging to complex Ginibre ensemble. The eigenenergies are shown to form stable structure in thermodynamical limit (large matrix dimension limit). Analogy of Wigner and Dyson with system of electrical charges is drawn.

## 2 Summary

A complex quantum system with energy dissipation is considered. The quantum Hamiltonians $`H`$ belong the complex Ginibre ensemble. The complex-valued eigenenergies $`Z_i`$ are random variables. The second differences $`\mathrm{\Delta }^1Z_i`$ are also complex-valued random variables. The second differences have their real and imaginary parts and also radii (moduli) and main arguments (angles). For $`N`$=3 dimensional Ginibre ensemble the distributions of above random variables are provided whereas for generic $`N`$\- dimensional Ginibre ensemble second difference distribution is analytically calculated. The law of homogenization of eigenergies is formulated. The analogy of Wigner and Dyson of Coulomb gas of electric charges is studied.

## 3 Introduction

We study generic quantum statistical systems with energy dissipation. The quantum Hamiltonian operator $`H`$ is in given basis of Hilbert’s space a matrix with random elements $`H_{ij}`$ . The Hamiltonian $`H`$ is not hermitean operator, thus its eigenenergies $`Z_i`$ are complex-valued random variables. We assume that distribution of $`H_{ij}`$ is governed by Ginibre ensemble . $`H`$ belongs to general linear Lie group GL($`N`$, C), where $`N`$ is dimension and C is complex numbers field. Since $`H`$ is not hermitean, therefore quantum system is dissipative system. Ginibre ensemble of random matrices is one of many Gaussian Random Matrix ensembles GRME. The above approach is an example of Random Matrix theory RMT . The other RMT ensembles are for example Gaussian orthogonal ensemble GOE, unitary GUE, symplectic GSE, as well as circular ensembles: orthogonal COE, unitary CUE, and symplectic CSE. The distributions of the eigenenergies $`Z_1,\mathrm{},Z_N`$ for $`N\times N`$ Hamiltonian matrices is given by Jean Ginibre’s formula :

$`P(z_1,\mathrm{},z_N)=`$ (1)
$`={\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{\pi j!}}{\displaystyle \underset{i<j}{\overset{N}{}}}|z_iz_j|^2\mathrm{exp}({\displaystyle \underset{j=1}{\overset{N}{}}}|z_j|^2),`$
where $`z_i`$ are complex-valued sample points ($`z_i𝐂`$). For Ginibre ensemble we define complex-valued spacings $`\mathrm{\Delta }^1Z_i`$ and second differences $`\mathrm{\Delta }^2Z_i`$:

$$\mathrm{\Delta }^1Z_i=Z_{i+1}Z_i,i=1,\mathrm{},(N1),$$
(2)
$$\mathrm{\Delta }^2Z_i=Z_{i+2}2Z_{i+1}+Z_i,i=1,\mathrm{},(N2).$$
(3)
The $`\mathrm{\Delta }^2Z_i`$ are extensions of real-valued second differences

$$\mathrm{\Delta }^2E_i=E_{i+2}2E_{i+1}+E_i,i=1,\mathrm{},(N2),$$
(4)
of adjacent ordered increasingly real-valued energies $`E_i`$ defined for GOE, GUE, GSE, and Poisson ensemble PE (where Poisson ensemble is composed of uncorrelated randomly distributed eigenenergies) .

There is an analogy of Coulomb gas of unit electric charges pointed out by Eugene Wigner and Freeman Dyson. A Coulomb gas of $`N`$ unit charges moving on complex plane (Gauss’s plane) C is considered. The vectors of positions of charges are $`z_i`$ and potential energy of the system is:

$$U(z_1,\mathrm{},z_N)=\underset{i<j}{}\mathrm{ln}|z_iz_j|+\frac{1}{2}\underset{i}{}|z_i^2|.$$
(5)
If gas is in thermodynamical equilibrium at temperature $`T=\frac{1}{2k_B}`$ ($`\beta =\frac{1}{k_BT}=2`$, $`k_B`$ is Boltzmann’s constant), then probability density function of vectors of positions is $`P(z_1,\mathrm{},z_N)`$ Eq. (1). Complex eigenenergies $`Z_i`$ of quantum system are analogous to vectors of positions of charges of Coulomb gas. Moreover, complex-valued spacings $`\mathrm{\Delta }^1Z_i`$ are analogous to vectors of relative positions of electric charges. Finally, complex-valued second differences $`\mathrm{\Delta }^2Z_i`$ are analogous to vectors of relative positions of vectors of relative positions of electric charges.

The $`\mathrm{\Delta }^2Z_i`$ have their real parts $`\mathrm{Re}\mathrm{\Delta }^2Z_i`$, and imaginary parts $`\mathrm{Im}\mathrm{\Delta }^2Z_i`$, as well as radii (moduli) $`|\mathrm{\Delta }^2Z_i|`$, and main arguments (angles) $`\mathrm{Arg}\mathrm{\Delta }^2Z_i`$.

## 4 Second Difference Distributions

We define following random variables for $`N`$=3 dimensional Ginibre ensemble:

$$Y_1=\mathrm{\Delta }^2Z_1,A_1=\mathrm{Re}Y_1,B_1=\mathrm{Im}Y_1,$$
(6)
$$R_1=|Y_1|,\mathrm{\Phi }_1=\mathrm{Arg}Y_1,$$
(7)
and for the generic $`N`$-dimensional Ginibre ensemble : $`W_1=\mathrm{\Delta }^2Z_1`$.

Their distributions for are given by following formulae :

$`f_{Y_1}(y_1)=f_{(A_1,B_1)}(a_1,b_1)=`$ (8)
$`={\displaystyle \frac{1}{576\pi }}[(a_1^2+b_1^2)^2+24]\mathrm{exp}({\displaystyle \frac{1}{6}}(a_1^2+a_2^2)).`$
$$f_{A_1}(a_1)=\frac{\sqrt{6}}{576\sqrt{\pi }}(a_1^4+6a_1^2+51)\mathrm{exp}(\frac{1}{6}a_1^2),$$
(9)
$$f_{B_1}(b_1)=\frac{\sqrt{6}}{576\sqrt{\pi }}(b_1^4+6b_1^2+51)\mathrm{exp}(\frac{1}{6}b_1^2),$$
(10)
$`f_{R_1}(r_1)=`$ (11)
$`\mathrm{\Theta }(r_1){\displaystyle \frac{1}{288}}r_1(r_1^4+24)\mathrm{exp}({\displaystyle \frac{1}{6}}r_1^2),`$
$`f_{\mathrm{\Phi }_1}(\varphi _1)={\displaystyle \frac{1}{2\pi }},\varphi _1[0,2\pi ].`$
$`P_3(w_1)=`$ (12)
$`=\pi ^3{\displaystyle \underset{j_1=0}{\overset{N1}{}}}{\displaystyle \underset{j_2=0}{\overset{N1}{}}}{\displaystyle \underset{j_3=0}{\overset{N1}{}}}{\displaystyle \frac{1}{j_1!j_2!j_3!}}I_{j_1j_2j_3}(w_1),`$
$`I_{j_1j_2j_3}(w_1)=`$ (13)
$`=2^{2j_2}{\displaystyle \frac{^{j_1+j_2+j_3}}{^{j_1}\lambda _1^{j_2}\lambda _2^{j_3}\lambda _3}}F(w_1,\lambda _1,\lambda _2,\lambda _3)|_{\lambda _i=0},`$
$`F(w_1,\lambda _1,\lambda _2,\lambda _3)=`$ (14)
$`=A(\lambda _1,\lambda _2,\lambda _3)\mathrm{exp}[B(\lambda _1,\lambda _2,\lambda _3)|w_1|^2],`$
$`A(\lambda _1,\lambda _2,\lambda _3)=`$ (15)
$`={\displaystyle \frac{(2\pi )^2}{(\lambda _1+\lambda _2\frac{5}{4})(\lambda _1+\lambda _3\frac{5}{4})(\lambda _11)^2}},`$
$`B(\lambda _1,\lambda _2,\lambda _3)=`$ (16)
$`=(\lambda _11){\displaystyle \frac{2\lambda _1\lambda _2\lambda _3+\frac{1}{2}}{2\lambda _1+\lambda _2+\lambda _3\frac{9}{2}}}.`$
## 5 Conclusions

We compare second difference distributions for different ensembles by defining following dimensionless second differences:

$$C_\beta =\frac{\mathrm{\Delta }^2E_1}{<S_\beta >},$$
(17)
$$X_1=\frac{A_1}{<R_1>},$$
(18)
where $`<S_\beta >`$ are the mean values of spacings for GOE(3) ($`\beta =1`$), for GUE(3) ($`\beta =2`$), for GSE(3) ($`\beta =4`$), for PE ($`\beta =0`$) , and $`<R_1>`$ is mean value of radius $`R_1`$ for $`N`$=3 dimensional Ginibre ensemble .

On the basis of comparison of results for Gaussian ensembles, Poisson ensemble, and Ginibre ensemble we formulate homogenization law : Eigenenergies for Gaussian ensembles, for Poisson ensemble, and for Ginibre ensemble tend to be homogeneously distributed. The second differences’ distributions assume global maxima at origin for above ensembles. For Coulomb gas the vectors of relative positions of vectors of relative positions of charges statistically vanish. It can be called stabilisation of structure of system of electric charges.

## 6 Acknowledgements

It is my pleasure to most deeply thank Professor Jakub Zakrzewski for formulating the problem.