text:
# Data-driven derivation of the turbulent energy cascade generator

## I Introduction

The inertial-range dynamics of fully developed turbulent flows reveals an ubiquitous multiscale character, driven by large-scale forcing on one end and controlled by fluid viscosity and small-scale dissipation on the other. For the intermediate scale range, the underlying Navier-Stokes equations do not disclose any other distinguishing length scales, thus suggesting a picture of a scale-invariant inertial-range dynamics MON71 ; FRI95 ; comment . This scale-invariance should reflect itself in the scaling of structure functions, which are moments of velocity increments constructed from the measured velocity data. However, measured structure functions in shear flows seem to show only an approximate multiscaling character, even at the largest accessible Reynolds numbers, where the large-scale and the dissipation scale are separated by five orders of magnitude ARN96 ; SRE97 ; SREE98 . The reason for this distortion of scaling seems to be the sensitivity of structure functions to the mean shear that is inevitably present in most natural flows at high Reynolds numbers. While schemes have been proposed to account for these effects arad ; SREE98 , it has been demonstrated recently CLE03  that the lowest-order two-point correlation function of the surrogate energy dissipation of a high-Reynolds-number atmospheric boundary layer reveals a far more convincing power-law scaling over the entire inertial range, without having to resort to shear corrections of the sort needed for structure functions. This finding suggests that the phenomenological picture of the self-similar turbulent energy cascade may contain more truth than previously anticipated. Consequently, it is reasonable to reexamine simple empirical energy cascade models and extend to a higher level the comparison between the models and experimental data. This is the purpose of the paper.

Random multiplicative cascade processes (RMCP) present a particularly simple geometrical picture of the energy cascade FRI95  and are a natural archetype for multiscaling. In their binary version, for example MS , they introduce a hierarchy of length scales $`l_j=L/2^j`$ and a random multiplicative cascade generator is used to transport the energy flux from the integral scale $`L`$ through the inertial range scales $`Ll_j\eta `$ down to the dissipation scale $`\eta =L/2^J`$. Referring to the ultrametric branching structure of binary RMCPs, $`N`$-point correlation functions of arbitrary order have been calculated analytically with generating function techniques GRE95 ; GRE96 ; GRE98a ; GRE98b . In terms of the ultrametric two-point distance, which measures the number of cascade steps necessary to reach the last common branching, these models reveal rigorous multiscaling. However, in this form, the two-point RMCP statistics cannot be compared directly to quantities that can be extracted from the data because, from an experimentalist’s perspective, the two-point correlations are expressed as functions of the Euclidean two-point distance. In this respect, RMCPs are incomplete and have to be supplemented by model-dependent extensions, stating the conversion of the ultrametric to Euclidean statistics JOU02 ; EGG01 ; MEN90 . It may be expected that this unavoidable conversion will lead to some deviations from rigorous multiscaling of correlation functions for two-point distances within the inertial range. The challenge is to find the degree to which the theoretical, extended-RMCP, two-point correlation functions can match their experimental counterparts.

If this comparison is satisfactory, the two-point correlation statistics of the energy dissipation will allow for an unambiguous extraction of multiscaling exponents, these being related to moments of the RMCP generator. In addition, two-point cumulants of logarithmic energy dissipation EGG01  allow us to access properties of the cascade generator in another way. Taken together, reliable parametric estimations of the underlying cascade generator can then be given. This would be the best we can do to achieve a data-driven derivation of the turbulent energy cascade generator and to settle the question of whether the latter is bimodal SCH85 ; MEN87 , log-normal OBU62 ; KOL62 , log-stable KID91 ; SCH92 , log-infinite divisible SHE94 ; DUB94 ; SHE95  or another type of distribution. None of those prototype cascade generators has been directly derived from data.

Note that the breakup coefficients, sometimes called multipliers, have once been thought to represent a direct approach to derive the RMCP cascade generator from data CHAB92 ; SRE95 ; MOL95 ; PED96 ; NEL96 . In fact, generator-like scale-invariant distributions of breakup coefficients have been observed, but conditional distributions have been found to exhibit scale correlations. In a series of papers JOU99 ; JOU00 ; JOU02  these findings have been fully explained within ultrametric hierarchical RMCPs, once those are analyzed from an experimentalist’s perspective, including unavoidable small-scale resummation and restoration of spatial homogeneity via the ultrametric-Euclidean conversion. This work has demonstrated that cascade generators and distributions of breakup coefficients are not directly related.

The structure of this paper is as follows. Section II provides basic information about the two-point correlation of the lowest order for the surrogate energy dissipation in four fully developed turbulent flows—one atmospheric boundary layer and three wind-tunnel flows. This information serves to provide the needed guidance for the subsequent RMCP modelling of Sec. III, where analytic expressions for two-point correlators and cumulants are derived, taking into account the ultrametric-Euclidean conversion. More information on data analysis is provided in Sec. IV, leading in Sec. V to further comparisons between prototype distributions and experiments. As will also be demonstrated in this section, the extracted properties of the generator allow us to make parametric estimations of the RMCP cascade generator. A conclusion and outlook is given in Sec. VI.

## II Data analysis I: basic facts on two-point statistics

We analyze four different data sets, three of which have been recorded in a wind tunnel PEA02  and the fourth in an atmospheric boundary layer DHR00  about 35 m above the ground. We will refer to them as w1, w2, w3 and a, respectively. Characteristic quantities of all data sets are summarized in Table 1. The Reynolds number $`R_\lambda =\sqrt{u^2}\lambda /\nu `$ is based on the Taylor microscale $`\lambda =\sqrt{u^2/(u/x)^2}`$; $`\nu `$ is the kinematic viscosity and $`u`$ is the streamwise velocity component. Upon the application of the frozen flow hypothesis, the recorded time series were converted into one-dimensional spatial series. The energy spectra of all four records reveal more or less the typical $`5/3`$ slope in the inertial range. In contrast to wind-tunnel records, that from the atmospheric boundary layer reveals a white-noise behavior at very small scales; this noise, which comes from detailed electronic circuitry, has been removed by an appropriate Wiener filter. The energy dissipation was then calculated as the surrogate amplitude $`\epsilon =15\nu (\frac{du}{dx})^2`$. Various tests were made to ensure that the effect of Wiener filtering were not consequential.

Figure 1 illustrates the lowest-order two-point correlator $`r_{1,1}(d=|x_2x_1|)=\epsilon (x_1)\epsilon (x_2)/(\epsilon (x_1)\epsilon (x_2))`$ sampled from the four different experimental records. Well inside the inertial range $`\eta dL`$ the two-point correlators reveal a power-law behavior $`r_{1,1}(L/d)^{\tau _{1,1}}`$. Power-law fits are indicated by the shifted broken straight lines; see also the insets, where the local slopes are shown. The resulting scaling exponents are $`\tau _{1,1}=0.15`$ (w1), $`0.14`$ (w2), $`0.18`$ (w3) and $`0.20`$ (a). Note that there is a Reynolds number dependence of this exponent. This has been explored in greater detail elsewhere CLE03a .

For the records with the largest Reynolds number, there is a large scale range for which the two-point correlator exhibits a rigorous power-law scaling. However, with decreasing Reynolds number this scale range becomes smaller. As a rule of thumb, we observe that a good scaling range is confined between $`d2030\eta `$ and $`0.5L`$. If we could understand precisely the deviations from the power-law scaling beyond this intermediate inertial range, a more satisfactory extraction of multiscaling exponents would be possible, especially for turbulent flows with moderate Reynolds numbers. In the small-distance $`d2030\eta `$ region, however, the correlations observed exceed the power-law extrapolation; consult Fig. 1 again. As has been explained in Ref. CLE03 , this enhancement is a consequence of the unavoidable surrogacy of the experimentally measured energy dissipation. Without knowing the correct small-distance behavior based on the proper energy dissipation, a theoretical modelling of two-point correlations in the dissipative regime is not very meaningful. We are thus left to inspect only the remainder of the inertial range. Fortunately, this range is sufficiently large, especially for the atmospheric data.

At $`dL`$, the two-point correlator has not yet converged to unity, which is its asymptotic value as $`d\mathrm{}`$. For all four data records inspected, the decorrelation length appears to be around $`L_{\mathrm{dec}}4L`$ and matches the length observed in the autocorrelation function of the streamwise velocity component. These findings suggest that the two-point correlator can be described as

$$r_{1,1}(2030\eta <dL_{\mathrm{dec}})=a_{1,1}\left(\frac{L_{\mathrm{dec}}}{d}\right)^{\tau _{1,1}}f_{1,1}(d/L_{\mathrm{dec}}).$$
(1)
For two-point distances $`d`$ much smaller than the decorrelation length the finite-size function converges to unity; that is, $`f_{1,1}(dL_{\mathrm{dec}})1`$. One goal of this paper is to qualitatively and quantitatively reproduce this functional form with an extended modelling based on random multiplicative cascade processes. This approach, which is the subject of the next section, naturally suggests a physical interpretation of the decorrelation length $`L_{\mathrm{dec}}`$.

## III Two-point statistics of random multiplicative cascade processes

### III.1 Binary random multiplicative cascade process

In its simplest form an RMCP employs a binary hierarchy of length scales $`l_j=L_{\mathrm{casc}}/2^j`$. In the first cascade step the parent interval of length $`l_0=L_{\mathrm{casc}}`$ is split into left and right daughter intervals, both of length $`l_1`$. In subsequent cascade steps, each interval of generation $`0jJ1`$ is again split into a left and right subinterval of length $`l_{j+1}=l_j/2`$. Once the dissipation scale $`\eta =L_{\mathrm{casc}}/2^J`$ is reached, the interval splitting stops and has resulted into $`2^J`$ spatially ordered intervals of smallest size $`\eta `$. It is convenient to label them as well as their ancestors according to the binary notation $`\kappa ^{(j)}=\kappa _1,\kappa _2,\mathrm{},\kappa _j`$. The label refers to the hierarchical position of an interval of generation $`j`$, where $`\kappa _i=0`$ or $`1`$ stands for the left or right interval, respectively.

The binary interval splittings go together with a probabilistic evolution of the energy-flux field. From generation $`j`$ to $`j+1`$ the field amplitudes propagate locally as

$`\mathrm{\Pi }_{\kappa ^{(j)},0}`$ $`=`$ $`q_{\kappa ^{(j)},0}\mathrm{\Pi }_{\kappa ^{(j)}},`$
$`\mathrm{\Pi }_{\kappa ^{(j)},1}`$ $`=`$ $`q_{\kappa ^{(j)},1}\mathrm{\Pi }_{\kappa ^{(j)}}.`$ (2)
The two random multiplicative weights $`q_{\mathrm{left}}=q_{\kappa ^{(j)},0}`$ and $`q_{\mathrm{right}}=q_{\kappa ^{(j)},1}`$, with mean $`q_{\mathrm{left}}=q_{\mathrm{right}}=1`$, are drawn from a scale-independent bivariate probability density function $`p(q_{\mathrm{left}},q_{\mathrm{right}})`$, which is called the cascade generator. Initially, corresponding to $`j=0`$, the iteration (III.1) starts with a given large-scale energy flux $`\mathrm{\Pi }`$, which might itself be a random variable fluctuating around its normalized mean $`\mathrm{\Pi }=1`$. After the last iteration $`J1J`$, the energy-flux amplitude is interpreted as the amplitude $`\epsilon _{\kappa ^{(J)}}=\mathrm{\Pi }_{\kappa ^{(J)}}`$ of the energy dissipation supported at the interval of length $`\eta `$ at position $`\kappa ^{(J)}`$. As a result of (III.1), this amplitude is a multiplicative sum of the random weights, given by

$$\epsilon _{\kappa _1,\mathrm{},\kappa _J}=q_{\kappa _1,\mathrm{},\kappa _J}q_{\kappa _1,\mathrm{},\kappa _{J1}}\mathrm{}q_{\kappa _1}\mathrm{\Pi }.$$
(3)
In the following we assume the cascade generator to be of the factorized form $`p(q_{\mathrm{left}},q_{\mathrm{right}})=p(q_{\mathrm{left}})p(q_{\mathrm{right}})`$ with identical statistics for the left/right variable. Of course, the factorization is not the most general ansatz but, as already pointed out in Ref. JOU00 , it represents a reasonable approximation: the turbulent energy cascade takes place in three spatial dimensions and calls for a three-dimensional RMCP modelling, respecting energy conservation. Since the measured temporal records come from one-dimensional cuts, the three-dimensional RMCP has to be observed in unity subdimension. Because of this, the RMCP appears to be non-conservative and the two multiplicative weights appear to be almost decorrelated and independent of each other.

### III.2 Two-point correlators

Expressions for $`N`$-point moments $`\epsilon (x_1)^{n_1}\mathrm{}\epsilon (x_N)^{n_N}`$ are easily found. They can be calculated either by a straightforward approach or, more formally, by an iterative construction of the respective multivariate characteristic function GRE95 ; GRE96 ; a third and more elegant approach BIA99  makes use of the full analytic solution of the multivariate characteristic function for logarithmic cascade-field amplitudes GRE98a ; GRE98b . We simply state the results up to two-point correlations:

$`\epsilon (x_1)^{n_1}`$ $`=`$ $`q^{n_1}^J\mathrm{\Pi }^{n_1},`$ (4)
$`\epsilon (x_1)^{n_1}\epsilon (x_2)^{n_2}`$ $`=`$ $`q^{n_1+n_2}^{JD}q^{n_1}^Dq^{n_2}^D\mathrm{\Pi }^{n_1+n_2}.`$ (5)
Here, the binary notation $`\kappa ^{(J)}`$ has been transformed into a spatial bin label $`x=1+_{j=1}^J\kappa _j2^{Jj}`$, which runs over $`1x2^J`$ in units of $`\eta `$. Two bins $`x_1(\kappa _1,\mathrm{},\kappa _{JD},\kappa _{JD+1},\mathrm{},\kappa _J)`$ and $`x_2(\kappa _1,\mathrm{},\kappa _{JD},\kappa _{JD+1}^{},\mathrm{},\kappa _J^{})`$ are assigned an ultrametric distance $`D`$ once the first $`JD`$ $`\kappa `$’s are identical and $`\kappa _{JD+1}\kappa _{JD+1}^{}`$. In other words, after $`JD`$ common branches along the binary tree, the two bins separate into different branches.

For the extraction of scaling exponents

$$\tau _n=\mathrm{log}_2q^n,$$
(6)
it is enough to consider the two-point statistics (5). In normalized form, the two-point correlators are found to scale perfectly as

$`r_{n_1,n_2}(D)`$ $`=`$ $`{\displaystyle \frac{\epsilon (x_1)^{n_1}\epsilon (x_2)^{n_2}}{\epsilon (x_1)^{n_1}\epsilon (x_2)^{n_2}}}`$ (7)
$`=`$ $`{\displaystyle \frac{\mathrm{\Pi }^{n_1+n_2}}{\mathrm{\Pi }^{n_1}\mathrm{\Pi }^{n_2}}}\left({\displaystyle \frac{q^{n_1+n_2}}{q^{n_1}q^{n_2}}}\right)^{JD}`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Pi }^{n_1+n_2}}{\mathrm{\Pi }^{n_1}\mathrm{\Pi }^{n_2}}}\left({\displaystyle \frac{L_{\mathrm{casc}}}{l_D}}\right)^{\tau _{n_1,n_2}},`$
where $`l_D=2^{D1}\eta `$ represents the characteristic two-bin distance corresponding to the ultrametric distance $`D>0`$ and

$$\tau _{n_1,n_2}=\tau _{n_1+n_2}\tau _{n_1}\tau _{n_2}.$$
(8)
For an experimentalist, the expression (7) does not present an observable result. Different pairs of bins, all having an identical Euclidean distance $`\eta d<L_{\mathrm{casc}}`$, do not have an unequivocal ultrametric distance. Depending on their position within the binary ultrametric cascade tree, the two bins might share a cascade history that is long (small $`D`$) or short (large $`D`$). Consequently, as an experimentalist analyzes the two-point statistics in terms of $`d`$, the ultrametric expression (7) has to be averaged over all $`D`$ that contribute to the same value of $`d`$. In order to perform this conversion from an ultrametric to an Euclidean distance and, by this means, to restore spatial homogeneity, we introduce the discrete conditional probability distribution

$$p(D|d)=\{\begin{array}{ccc}0\hfill & & (1D<A=\mathrm{log}_2d)\hfill \\ 1d2^D\hfill & & (D=A)\hfill \\ d2^D\hfill & & (A<DJ)\hfill \\ 0\hfill & & (J<D<\mathrm{})\hfill \\ d2^J\hfill & & (D=\mathrm{})\hfill \end{array}$$
(9)
of finding the ultrametric distance $`D`$ for a given Euclidean distance $`d`$ in units of $`\eta `$ EGG01 . This expression has been derived by employing the chain picture of independent cascade configurations; consult Fig. 2. It roughly goes as $`p(D|d)2^{\mathrm{log}_2dD}`$. The sum $`_{D=0}^Jp(D|d)=1d2^J`$ does not add up to unity, since $`p(\mathrm{}|d)=d2^J`$ represents the probability that the two bins belong to different $`L_{\mathrm{casc}}`$-domains.

Since the one-point statistics $`\epsilon (x)^n=\epsilon (x+d)^n`$ do not depend on the spatial index $`x`$, the ultrametric-Euclidean conversion of the normalized two-point correlator (7) leads to

$`r_{n_1,n_2}(d0)`$ $`=`$ $`{\displaystyle \underset{D=1}{\overset{J}{}}}p(D|d)r_{n_1,n_2}(D)+p(\mathrm{}|d)`$ (10)
$`=`$ $`{\displaystyle \frac{\mathrm{\Pi }^{n_1+n_2}}{\mathrm{\Pi }^{n_1}\mathrm{\Pi }^{n_2}}}\left(1{\displaystyle \frac{d}{2^A}}+{\displaystyle \frac{d}{2^A}}\left(2{\displaystyle \frac{q^{n_1+n_2}}{q^{n_1}q^{n_2}}}1\right)^1\right)\left({\displaystyle \frac{q^{n_1+n_2}}{q^{n_1}q^{n_2}}}\right)^{JA}`$
$`+\left(1{\displaystyle \frac{\mathrm{\Pi }^{n_1+n_2}}{\mathrm{\Pi }^{n_1}\mathrm{\Pi }^{n_2}}}\left(2{\displaystyle \frac{q^{n_1+n_2}}{q^{n_1}q^{n_2}}}1\right)^1\right){\displaystyle \frac{d}{2^J}}.`$
This expression holds for every $`\eta dL_{\mathrm{casc}}`$. For $`d=0`$ and $`d>L_{\mathrm{casc}}`$ the normalized two-point density simply becomes $`r_{n_1,n_2}(d=0)=\mathrm{\Pi }^{n_1+n_2}q^{n_1+n_2}^J/(\mathrm{\Pi }^{n_1}q^{n_1}^J\mathrm{\Pi }^{n_2}q^{n_2}^J)`$ and $`r_{n_1,n_2}(d>L_{\mathrm{casc}})=1`$, respectively. The two-point density (10) does not reveal perfect scaling anymore. Usually the second term, scaling as $`d/L_{\mathrm{casc}}`$, is small when compared to the first term, except for $`dL_{\mathrm{casc}}`$. The modulations, observed for the first term, are an artifact of the discrete scale invariance SOR98  of the binary random multiplicative cascade model implementation. In the following we will discard these modulations by first considering only dyadic distances $`d_m=L_{\mathrm{casc}}/2^m`$ with integer $`0m<J`$, and then switching again to continuous $`d`$ by interpolating between the discrete $`d_m`$. The expression (10) then simplifies to

$$r_{n_1,n_2}(d)=a_{n_1,n_2}\left(\frac{L_{\mathrm{casc}}}{d}\right)^{\tau _{n_1,n_2}}f_{n_1,n_2}(d/L_{\mathrm{casc}}),$$
(11)
with

$$a_{n_1,n_2}=\frac{\mathrm{\Pi }^{n_1+n_2}}{\mathrm{\Pi }^{n_1}\mathrm{\Pi }^{n_2}}\left(2\frac{q^{n_1+n_2}}{q^{n_1}q^{n_2}}1\right)^1$$
(12)
and the finite-size scaling function

$$f_{n_1,n_2}(d/L_{\mathrm{casc}})=1+\left(a_{n_1,n_2}^11\right)\left(\frac{d}{L_{\mathrm{casc}}}\right)^{1+\tau _{n_1,n_2}}.$$
(13)
Figure 3 compares the expressions (10) and (11) for the order $`n_1=n_2=1`$.

The finite-size scaling function has the property $`f_{n_1,n_2}(dL_{\mathrm{casc}})=1`$ as long as the condition $`1+\tau _{n_1,n_2}0`$ or, equivalently, $`q^{n_1+n_2}/(q^{n_1}q^{n_2})>1/2`$ is fulfilled. This is the case for all positive combinations $`n_10,n_20`$. However, combinations with negative orders do exist, for which the second term on the right hand side of (13) then dominates over the first term in the limit $`d/L_{\mathrm{casc}}0`$. This implies that the normalized two-point density (11) asymptotically scales as $`r_{n_1,n_2}(d)(L_{\mathrm{casc}}/d)^1`$, giving rise to the effective scaling exponents $`\tau _{n_1,n_2}^{\mathrm{eff}}=sup\{1,\tau _{n_1,n_2}\}`$. This scaling transition is again a pure consequence of the ultrametric-Euclidean conversion. More discussions on this subject can be found in Refs. MEN90 ; ONE93 ; BEN97 .

Upon studying the expression (13) more closely, we realize that two effects, ultrametric-Euclidean conversion and large-scale fluctuations, contribute to the finite-size scaling function. They have a tendency to cancel each other. Once we have

$$\frac{\mathrm{\Pi }^{n_1+n_2}}{\mathrm{\Pi }^{n_1}\mathrm{\Pi }^{n_2}}=2\frac{q^{n_1+n_2}}{q^{n_1}q^{n_2}}1,$$
(14)
the finite-size scaling function becomes $`f_{n_1,n_2}(d/L_{\mathrm{casc}})=1`$ exactly, showing no $`d`$-dependence.

We also wish to point out an interesting mathematical observation following from the specific expressions (11)-(13). Since the finite-size scaling function (13) reveals the simple scaling behavior

$$\left(f_{n_1,n_2}(d/L_{\mathrm{casc}})1\right)\left(\frac{L_{\mathrm{casc}}}{d}\right)^{\tau _{n_1,n_2}}\left(\frac{d}{L_{\mathrm{casc}}}\right)^\mu $$
(15)
with $`\mu =1`$, we find

$$r_{n_1,n_2}(d)\frac{1}{\xi ^\mu }r_{n_1,n_2}(\xi d)=a_{n_1,n_2}\left(1\frac{1}{\xi ^{\mu +\tau _{n_1,n_2}}}\right)\left(\frac{L_{\mathrm{casc}}}{d}\right)^{\tau _{n_1,n_2}},$$
(16)
where the normalized two-point correlator with the rescaled two-point distance $`\xi d`$ has been subtracted from itself. As a function of the two-point distance $`d`$ this quantity exhibits a rigorous power-law behavior with scaling exponents $`\tau _{n_1,n_2}`$, which is independent of the chosen rescaling parameter $`\xi `$ and is free of large-scale effects.

### III.3 Two-point cumulants

Because experimental data yield limited statistics, the two-point correlation densities (11) will be restricted to lowest orders $`1n_1+n_23`$ or $`4`$. This limits us to indirect information on the cascade generator, namely the scaling exponents $`\tau _1,\tau _2,\tau _3`$ and, perhaps, $`\tau _4`$ of (6). In order to do better, we need to accumulate additional and complementary information. In fact, as proposed already in Ref. EGG01 , this can be achieved by switching to the logarithmic amplitude $`\epsilon (x)\mathrm{ln}\epsilon (x)`$, and from two-point correlation densities to two-point cumulants

$`C_{1,1}(x_2x_1)`$ $`=`$ $`\mathrm{ln}\epsilon (x_1)\mathrm{ln}\epsilon (x_2)\mathrm{ln}\epsilon (x)^2,`$
$`C_{2,1}(x_2x_1)`$ $`=`$ $`\mathrm{ln}^2\epsilon (x_1)\mathrm{ln}\epsilon (x_2)2\mathrm{ln}\epsilon (x_1)\mathrm{ln}\epsilon (x_2)\mathrm{ln}\epsilon (x)\mathrm{ln}^2\epsilon (x)\mathrm{ln}\epsilon (x)+2\mathrm{ln}\epsilon (x)^3,\mathrm{},`$
$`C_{n_1,n_2}(x_2x_1)`$ $`=`$ $`{\displaystyle \frac{^{n_1+n_2}}{\lambda ^{n_1}\lambda ^{n_2}}}\mathrm{ln}\epsilon (x_1)^{\lambda _1}\epsilon (x_2)^{\lambda _2}|_{\lambda _1=\lambda _2=0}.`$ (17)
Explicit RMCP expressions have already been derived in Ref. EGG01  within the ultrametric view as well as the converted ultrametric-Euclidean view. Here, we summarize only the latter lowest-order results, which hold for $`\eta dL_{\mathrm{casc}}`$:

$$C_{n1,1}(d)=G_1(J,d)c_n+G_0(J,d)\mathrm{ln}^n\mathrm{\Pi }_c.$$
(18)
The geometric functions $`G_n(J,d)=(1/2^J)_{D=1}^J(JD)^np(D|d)`$ are related to moments of the conditional probability distribution (9) and are fingerprints of the hierarchical RMCP tree structure. They are given by the expressions

$`G_0(J,d)`$ $`=`$ $`1{\displaystyle \frac{d}{L_{\mathrm{casc}}}},`$
$`G_1(J,d)`$ $`=`$ $`(JA)2{\displaystyle \frac{d}{\eta }}\left(2^A2^J\right)\mathrm{log}_2\left({\displaystyle \frac{L_{\mathrm{casc}}}{d}}\right)2+2{\displaystyle \frac{d}{L_{\mathrm{casc}}}},`$ (19)
with the last step neglecting small log-oscillations. The cumulants of the logarithmic multiplicative weight

$$c_n=\mathrm{ln}^nq_c=\frac{^nQ(\lambda )}{\lambda ^n}|_{\lambda =0}$$
(20)
are generated by the logarithm of the Mellin transform of the cascade generator, i.e.,

$$Q(\lambda )=\mathrm{ln}\left(𝑑qp(q)q^\lambda \right).$$
(21)
The cumulants $`\mathrm{ln}^n\mathrm{\Pi }_c`$ of the initial large-scale fluctuation are defined analogous to $`\mathrm{ln}^nq_c`$.

### III.4 Multifractal sum rules

The cumulants $`c_n`$ of (20) and the scaling exponents $`\tau _n`$ of (6) are not independent of each other. Combining Eqs. (6), (20) and (21), we arrive at

$$Q(n)=\mathrm{ln}2\tau _n=\underset{k=1}{\overset{\mathrm{}}{}}c_k\frac{n^k}{k!}.$$
(22)
In the lowest order, this translates to

$$\mathrm{ln}2\tau _1=0=c_1+\frac{c_2}{2}+\frac{c_3}{6}+\mathrm{}.$$
(23)
These multifractal sum rules can be used, for example, to estimate $`c_1`$, which can not be extracted from the two-point statistics (18).

Another approach to estimate the value of $`c_1`$ is given by the well-known replica trick:

$$c_1=\mathrm{ln}q=\frac{Q(\lambda )}{\lambda }|_{\lambda =0}=\frac{q^\lambda }{\lambda }|_{\lambda =0}=\underset{\lambda 0}{lim}\frac{q^\lambda 1}{\lambda }=\underset{\lambda 0}{lim}\frac{2^{\tau _\lambda }1}{\lambda }.$$
(24)
Another form of sum rules follows from (22) and states that

$$\mathrm{ln}2\frac{\tau _\lambda }{\lambda }|_{\lambda =n}=\underset{k=0}{\overset{\mathrm{}}{}}c_{k+1}\frac{n^k}{k!}=c_1+nc_2+\frac{n^2}{2}c_3+\mathrm{}.$$
(25)
It can be seen as a generalization of (24).

## IV Data analysis II: more on two-point statistics

In this section, we provide further aspects of data analysis that are relevant. The goal is three-fold: first and already pointed out in Sect. II, to test the RMCP expression (11) with the proposed finite-size scaling for two-point correlators, second, to test the expression (18) for two-point cumulants derived from the RMCP theory, and, third, to extract reliable values for the scaling exponents $`\tau _n`$ and cumulants $`c_n`$ from various turbulent records discussed in Sec. II.

The expressions (11) and (18) come with parameters $`L_{\mathrm{casc}}`$, $`\tau _{n_1,n_2}`$, $`a_{n_1,n_2}`$, $`c_n`$ and $`\mathrm{ln}^n\mathrm{\Pi }_c`$. The parameter $`L_{\mathrm{casc}}`$ depends on neither the order $`n_1,n_2`$ nor the choice of the two-point statistics, i.e., whether we use the correlator or the cumulant. The fits of (11) and (18) to their counterparts from experimental records have to respect this independence. In addition to the common parameter $`L_{\mathrm{casc}}`$, the fit of each order has two more parameters: these are $`\tau _{n_1,n_2}`$, $`a_{n_1,n_2}`$ for the two-point correlator and $`c_n`$, $`\mathrm{ln}^n\mathrm{\Pi }_c`$ for the two-point cumulant. We demonstrate the quality of fits by choosing the data record w2. For this data set, two-point correlators with orders from $`n_1,n_2=1,1`$ to $`2,2`$ are illustrated in Fig. 4, while the two-point cumulants with orders $`n_1,n_2=1,1`$ and $`2,1`$ are illustrated in Fig. 5. Except for very small two-point distances $`\eta d2030\eta `$, where, as already noted, the surrogacy effect of the energy dissipation corrupts the experimental two-point statistics CLE03 , the agreement between the experimental two-point correlators and cumulants and the best-fit expressions (11) and (18) is remarkable. Table II lists the best-fit parameters $`L_{\mathrm{casc}}`$, $`\tau _n`$ and $`c_n`$. Note that, due to (8) and $`\tau _1=0`$, the scaling exponents $`\tau _{n_1,n_2}`$ have been converted into $`\tau _2=\tau _{1,1}`$, $`\tau _3=\tau _{2,1}+\tau _2`$ and $`\tau _4=\tau _{3,1}+\tau _3=\tau _{2,2}+2\tau _2`$. For w2 the two values $`\tau _4=0.79`$ and $`0.77`$, the first value extracted from $`r_{3,1}(d)`$ and the second from $`r_{2,2}(d)`$, are consistent with each other, although the statistical convergence of the two two-point correlators of order $`n_1+n_2=4`$ is already beyond the limit of acceptability. For the record w3 a similar statement can be made, but the other records w1 and a are definitely confined to $`n_1+n_23`$. Their best-fit parameters $`L_{\mathrm{casc}}`$, $`\tau _n`$ and $`c_n`$ are also listed in Table II.

The cumulant $`c_1`$ cannot be directly extracted from two-point cumulants. However, an indirect extraction is possible via the relations (23) and (24). Using the numerical values of $`c_2`$ and $`c_3`$ determined already, the truncated multifractal sum rule (23) leads to the $`c_1`$ values listed in the second-last column of Table 2. Input for the replica-trick formula (24) are the scaling exponents $`\tau _1`$, $`\tau _2`$ and $`\tau _3`$ determined already. The application of a cubic spline for the extrapolation results in the last column of Table 2. If a linear spline is used instead, the resulting $`c_1`$ values of the records w1, w2, w3 shift closer to the values of the second-last column, but would do worse for the record a. For an order-of-magnitude estimate, however, it is safe to say that $`c_10.05`$ for the records w1, w2, w3 and $`c_10.07`$ for the record a.

This comparison of the RMCP theory with data demonstrates that the parameter $`L_{\mathrm{casc}}`$ is a meaningful quantity and deserves some fundamental consideration. For two-point distances $`2030\eta dL_{\mathrm{casc}}`$, the fitted expressions (11) and (18) are in good qualitative agreement with their experimental counterparts. At $`d=L_{\mathrm{casc}}`$ all two-point correlations decorrelate and become identical to their asymptotic values corresponding to $`d\mathrm{}`$; $`L_{\mathrm{casc}}>L`$ is somewhat larger than the operationally defined integral length $`L`$ of Table I, which is calculated as the correlation length of the velocity autocorrelation function. The physical interpretation of the extracted parameter $`L_{\mathrm{casc}}`$ is that of a turbulent cascade length which, according to Fig. 2, describes the spatial extension of a hierarchical RMCP domain.

The finite-size scaling $`f_{n_1,n_2}(d/L_{\mathrm{casc}})`$ of two-point correlators, predicted by RMCP, now allows for an unambiguous derivation of scaling exponents, even for fully developed turbulent flows with a moderate Reynolds number. Consequently, as a closer inspection of Table II shows, reliable statements that the intermittency exponent $`\tau _2`$ might show a weak dependence on the Reynolds number appear to be within reach. Of course, to make this statement solid, the analysis of many more records is needed; this effort will be reported in CLE03a .

## V Parametric estimation of the RMCP generator

If the scaling exponents $`\tau _n`$, or the cumulants $`c_n`$, exist and are known for all orders $`1n<\mathrm{}`$, the binary RMCP generator could in principle be reconstructed via the inverse transform of (21). Unfortunately, as we have seen in the previous section, reliable information is limited to the lowest orders. Hence, the best we can do is to use sophisticated parametric estimates. Section V.A lists some of the most popular parametrizations and compares their performance with the results listed in Table II. Section V.B introduces the so-called log-normal inverse Gaussian distribution BAR98 ; BAR01 , which represents a broader and more flexible parametrization class, with the purpose of finding a suitable approximation to the true cascade generator.

### V.1 Dictionary of prototype cascade generators

Here we list a number of popular generators $`p(q)`$ for binary random multiplicative cascade processes. They all have the property that the expectation value $`q=1`$.

The log-normal distribution

$$p_{\mathrm{lognormal}}(q)=\frac{1}{\sqrt{2\pi }\sigma q}\mathrm{exp}\left[\frac{1}{2\sigma ^2}\left(\mathrm{ln}q+\frac{\sigma ^2}{2}\right)^2\right]$$
(26)
is classic OBU62 ; KOL62 . Its log-stable generalization has also been considered KID91 ; SCH92 , but does not qualify for our purposes, since the cumulants $`c_n`$ do not exist for this distribution beyond some order. For comparison, we will also employ the rescaled gamma distribution

$$p_{\mathrm{gamma}}(q)=\frac{\gamma ^\gamma }{\mathrm{\Gamma }(\gamma )}q^{\gamma 1}e^{\gamma q}$$
(27)
and the asymmetric beta distribution JOU00

$$p_{\mathrm{beta}}(q)=\frac{\mathrm{\Gamma }(8\beta )}{\mathrm{\Gamma }(\beta )\mathrm{\Gamma }(7\beta )}\mathrm{\hspace{0.33em}8}^{18\beta }q^{\beta 1}(8q)^{7\beta 1}.$$
(28)
The bimodal distribution

$$p_{\mathrm{alpha}}(q)=\frac{\alpha _2}{\alpha _1+\alpha _2}\delta \left(q(1\alpha _1)\right)+\frac{\alpha _1}{\alpha _1+\alpha _2}\delta \left(q(1+\alpha _2)\right),$$
(29)
although discrete, has also been used extensively SCH85 ; MEN87 . Another popular discrete representative is the log-Poisson distribution

$$p_{\mathrm{logPoisson}}(q)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{2^{\nu _1}(\nu _1\mathrm{ln}2)^k}{k!}\delta \left(q2^{\nu _1(1\nu _2)}\nu _2^k\right),$$
(30)
which was originally derived with $`\nu _1=2`$ and $`\nu _2=2/3`$ from some plausible reasoning on the structure of the most singular objects in fully developed turbulent flows SHE94 ; DUB94 ; SHE95 .

For all parametrizations (26)–(30) it is straightforward to determine analytic expressions for the scaling exponents and cumulants via (20)–(22). The free parameter of the one-parametric distributions (26)–(28) is then fixed to reproduce the observed intermittency exponent $`\tau _2`$, listed in Table II. The two-parametric distributions (29) and (30) need also conform to $`\tau _3`$ in addition to $`\tau _2`$. No further freedom is left for the scaling exponents of higher order and cumulants of all orders. Table III summarizes their predicted values.

It is difficult to rate the prototype cascade generators because of ambiguity inherent in the data. Within the one-parametric distributions (26)-(28) the log-normal distribution performs better: for all the records, the predicted values for $`c_1`$ and $`c_2`$ are close to the observed cumulants. However, the log-normal distribution without skewness is unable to reproduce the observed positive values for $`c_3`$ and for the record w2 it also overestimates the scaling exponent $`\tau _4`$. Furthermore, the difficulties of the log-normal distribution for high-order moments is now well known FRI95 . Compared to the log-normal distribution, the rescaled gamma distribution and the asymmetric beta distribution have the tendency to overestimate the first two cumulants. Furthermore, $`c_3`$ is predicted with an opposite sign. Rather surprisingly, the simple two-parametric bimodal distribution (29) shows the closest agreement for all records. The scaling exponent $`\tau _4`$, if observed, as well as the cumulants $`c_1`$ and $`c_2`$ almost match their observed counterparts. Moreover, $`c_3`$ comes with the correct sign, although it is about a factor $`2`$ too low for the records w1, w2, w3 and roughly a factor $`2`$ too large for the atmospheric boundary layer record. Like the distributions (27) and (28), the two-parametric log-Poisson distribution overestimates the second cumulant $`c_2`$ and, except for record w3, predicts $`c_3`$ with the wrong sign. It is interesting to note that the parameter-free log-Poisson distribution SHE94  with $`\nu _1=2`$ and $`\nu _2=2/3`$ matches well the scaling exponents $`\tau _3`$ and $`\tau _4`$ of record a with the largest Reynolds number, but disagrees with all cumulants.

### V.2 Log-normal inverse Gaussian distribution

A broader and more flexible parametrization class is the so-called normal inverse Gaussian distribution BAR98 ; BAR01

$$p(x;\alpha ,\beta ,\mu ,\delta )=a(\alpha ,\beta ,\mu ,\delta )s\left(\frac{x\mu }{\delta }\right)^1K_1\left\{\delta \alpha s\left(\frac{x\mu }{\delta }\right)\right\}e^{\beta x},$$
(31)
with $`s(x)=\sqrt{(1+x^2)}`$, $`a(\alpha ,\beta ,\mu ,\delta )=\pi ^1\alpha \mathrm{exp}(\delta \sqrt{\alpha ^2\beta ^2}\beta \mu )`$ and $`\mathrm{}<x<\mathrm{}`$. $`K_1`$ is the modified Bessel function of the third kind and index 1. The domain of variation of the four parameters is given by $`\mu 𝐑`$, $`\delta 𝐑_+`$ and $`0|\beta |<\alpha `$. The distribution is denoted by $`\mathrm{NIG}(\alpha ,\beta ,\mu ,\delta )`$, and its cumulant generating function $`\mathrm{Q}(\lambda ;\alpha ,\beta ,\mu ,\delta )=\mathrm{ln}e^{\lambda x}`$ has the simple form

$$\mathrm{Q}(\lambda ;\alpha ,\beta ,\mu ,\delta )=\delta \left(\sqrt{\alpha ^2\beta ^2}\sqrt{\alpha ^2(\beta +\lambda )^2}\right)+\mu \lambda .$$
(32)
If $`x_1,\mathrm{},x_m`$ are independent normal inverse Gaussian random variables with common parameters $`\alpha `$ and $`\beta `$ but individual location-scale parameters $`\mu _i`$ and $`\delta _i`$ $`(i=1,\mathrm{},m)`$, then $`x_+=x_1+\mathrm{}+x_m`$ is again distributed according to a normal inverse Gaussian law with parameters $`\alpha `$, $`\beta `$, $`\mu _+`$ and $`\delta _+`$. Furthermore, we note that the $`\mathrm{NIG}`$ distribution (31) has semi-stretched tails

$$p(x;\alpha ,\beta ,\mu ,\delta )\left|x\right|^{3/2}\mathrm{exp}\left(\alpha \left|x\right|+\beta x\right)$$
(33)
as $`x\pm \mathrm{}`$. This result follows from the asymptotic relation $`K_\nu (x\mathrm{})\sqrt{\pi /2}x^{1/2}e^x`$.

For our purposes we assume the random multiplicative weight to be distributed according to

$$\mathrm{ln}q\mathrm{NIG}(\alpha ,\beta ,\delta ,\mu ),$$
(34)
which turns normal inverse Gaussian statistics into log-normal inverse Gaussian statistics. With (20), (22) and (32), the scaling exponents and cumulants yield

$$\tau _n\mathrm{ln}2=\mathrm{Q}(n;\alpha ,\beta ,\mu ,\delta )$$
(35)
and

$$c_1=\mu +\frac{\delta \rho }{\sqrt{1\rho ^2}},c_2=\frac{\delta }{\alpha (1\rho ^2)^{3/2}},c_3=\frac{3\delta \rho }{\alpha ^2(1\rho ^2)^{5/2}},\mathrm{},$$
(36)
where $`\rho =\beta /\alpha `$.

For each of the records w1, w2, w3 and a, the four NIG parameters $`\alpha `$, $`\beta `$, $`\delta `$ and $`\mu `$ are determined so as to reproduce $`\tau _1=0`$ and the observed values for $`\tau _2`$, $`\tau _3`$ and $`c_2`$ listed in Table 2. Since the respective expressions (35) and (36) are nonlinear, real solutions for the parameters are not guaranteed. Where complex-valued solutions resulted in the first attempt, which occurred for w2, w3 and a, the values for $`c_2`$, $`\tau _3`$ and $`\tau _2`$ are relaxed, in this order and to some small extent, until real-valued parameter solutions are obtained. The outcome is listed in Table 4. The results are very close to the log-normal values listed in 3. The log-normal inverse Gaussian distribution has the tendency to overestimate the fourth-order scaling exponent $`\tau _4`$. The magnitude of the third cumulant $`c_3`$ is strongly underestimated, so that its predicted sign shows only random scatter. Figure 6 compares the extracted log-normal inverse Gaussian distributions of Table 4 with the extracted log-normal distributions of Table 3. For all four records the extracted distributions are very similar.

As a summary of this section, we reiterate that the bimodal distribution (29) produces the best overall agreement with the observed scaling exponents and cumulants. However, the true cascade generator will not be discrete. From the set of continuous generator representatives tested, the log-normal and log-normal inverse Gaussian distributions perform best and about equally well.

## VI Conclusion

Random multiplicative cascade processes are able to describe the observed two-point correlation structure of the surrogate energy dissipation of fully developed turbulent flows beyond simple power-law scaling. Keeping in mind the need for a satisfactory comparison between modelling and experimental data, a useful transformation has been introduced: this transformation converts model-dependent and unobservable ultrametric two-point distances to Euclidean two-point distances, the latter reflecting the horizontally sampled $`n`$-point statistics of the experimental records. The predictions of RMCP for finite-size scaling of two-point correlation functions are confirmed by experimental data from three (wind tunnel) shear flows and one atmospheric boundary layer; a new physical length scale characterizing the upper end of the inertial range, called here the turbulent cascade length, has been shown to appear naturally. Furthermore, the quantitative classification of the deviations from a rigorous scaling of two-point correlators allows for an unambiguous extraction of multiscaling exponents, even for flows with moderate Reynolds numbers. When complemented with additional information extracted from two-point cumulants of the logarithmic energy dissipation, reliable parametric estimates of the RMCP generator have become feasible.

RMCPs produce a consistent geometrical modelling of the self-similar turbulent energy cascade. This is further supported by recent investigations on the scaling part of three-point statistics SCH03  and previous investigations on scale correlations JOU99 ; JOU00 ; JOU02 ; CLE00 . Thus, the self-similar and RMCP-like cascade process appears to be universal. However, as shown by variations among the four records considered, the strength of the cascade generator appears to depend on the Reynolds number and perhaps also the flow geometry. This needs further clarification by means of a separate and extended effort CLE03a ; needless to say, there is a strong need for the experimentalists to produce clean and longer records with converged statistics.

###### Acknowledgements.

The authors acknowledge fruitful discussions with Hans C. Eggers and Markus Abel.


format: text

text:
# Extrasolar Planet Orbits and Eccentricities

## 1 1. Mass and period distributions

The mass distribution of extrasolar planets is sharply cut off above $`10M_J`$, where $`M_J`$ is the Jupiter mass Zucker & Mazeh (2001); Jorissen et al. (2001). The absence of companions with masses in the range $`10M_J<M<100M_J`$ (the “brown dwarf desert”) provides the strongest reason to believe that the extrasolar planets are formed by a different mechanism than low-mass companion stars, and thus the criterion $`M<10M_J`$ is probably the least bad definition of an extrasolar planet. The observed mass distribution for $`M<10M_J`$ can be modeled as a power law $`dnM^\alpha dM`$ with $`\alpha =1.1\pm 0.1`$ Zucker & Mazeh (2001); Heacox (1999); Stepinski & Black (2001); Tabachnik & Tremaine (2002); in other words the distribution is approximately flat in $`\mathrm{log}M`$, at least over the range $`M_J\stackrel{<}{}M\stackrel{<}{}10M_J`$. Weak evidence that this distribution can be extrapolated to smaller masses comes from the solar system: fitting the distribution of the masses of the 9 planets to a power law over the interval from $`M_{\mathrm{Pluto}}`$ to $`M_J`$ yields $`\alpha =1.0\pm 0.1`$, with a high Kolmogorov-Smirnov (KS) confidence level (over 99.9%).

The period distribution can also be described by a power law to a first approximation: $`dnP^\beta dP`$, where $`\beta =0.73\pm 0.06`$ Tabachnik & Tremaine (2002). Some analyses Heacox (1999); Stepinski & Black (2001); Mazeh & Zucker (2002) find a steeper slope, $`\beta 1`$, probably because they do not correct for selection effects (the reflex motion of the star due to a planet of a given mass declines with increasing period). Stepinski & Black (2001) and Mazeh & Zucker (2002) stress that the distribution of periods (and eccentricities) is almost the same for extrasolar planets and spectroscopic binaries—only the mass distributions are different—and this striking coincidence demands explanation (see §2.7). Recall that the standard minimum solar nebula has surface density $`\mathrm{\Sigma }(r)r^{1.5}`$ Hayashi (1981). If the formation process assigned this mass to equal-mass planets without migration we would have $`\beta =0.67`$, consistent with the exponent found above for extrasolar planets.

In Figure 1 we plot the masses and periods of extrasolar planets.111The data for figures are taken from http://www.obspm.fr/planets as of October 2003. Although power laws are good first approximations to the mass and period distributions for $`M<10M_J`$, there have been a number of suggestions of additional structure. (i) Zucker & Mazeh (2002) and Udry et al. (2003) suggest that there are too few massive planets ($`M\stackrel{>}{}3M_J`$) on short-period orbits ($`P\stackrel{<}{}50`$ days), although (ii) Zucker & Mazeh argue that this shortfall is not present if the host star is a member of a wide binary system (denoted by star symbols in Figure 1). In addition, Udry et al. (2003) suggest that there is (iii) a deficit of planets in the range $`10\text{ d}\stackrel{<}{}P\stackrel{<}{}100\text{ d}`$; and (iv) a deficit of low-mass planets ($`M\stackrel{<}{}0.8M_J`$) on long-period orbits ($`P\stackrel{>}{}100\text{ d}`$). We do not find the latter two features persuasive, since their significance has not yet been verified by statistical tests, and of course low-mass planets on long-period orbits produce the smallest reflex velocity and so are most difficult to detect.

Whether or not these features exist, in our view the most impressive feature of the period distribution is that it is described so well by a power law. In particular, the mechanisms that might halt planetary migration often operate most effectively close to the star Terquem (2003), so it is natural to expect a “pile-up” of planets at orbital periods of a few days; however, there is only weak evidence for an excess of this nature compared to a power-law distribution.

## 2 2. Eccentricity distribution

One of the remarkable features of extrasolar planets is their large eccentricities: the median eccentricity of 0.28 is larger than the maximum eccentricity of any planet in the solar system. The eccentricity distribution of extrasolar planets differs from that of solar system planets at the 99.9% KS confidence level.

Observational selection effects in radial-velocity surveys could either favor or disfavor detection of high-eccentricity planets: as the eccentricity grows at fixed period the periastron velocity grows, but the time spent near periastron declines. The limited simulations that have been done Fischer & Marcy (1992); Mazeh et al. (1996) suggest that these two effects largely cancel, so the eccentricity-dependent selection effects are small. Thus, the observed distribution should accurately reflect the eccentricities of massive, short-period extrasolar planets.

Planets that form from a quiescent disk are expected to have nearly circular orbits, and thus the high eccentricities of extrasolar planets demand explanation. A wide variety of eccentricity excitation mechanisms has been suggested so far.

### 2.1 2.1. Interactions with the protoplanetary gas disk

Gravitational interactions between a planet and the surrounding protoplanetary disk can excite or damp the planet’s eccentricity. These interactions are concentrated at discrete resonances, the most important of which are given by the relation Goldreich & Tremaine (1978, 1980, 1981); Ward (1986); Artymowicz (1993)

$$\mathrm{\Omega }+\frac{ϵ}{m}(\mathrm{\Omega }\dot{\varpi })=\mathrm{\Omega }_p+\frac{ϵ_p}{m}(\mathrm{\Omega }_p\dot{\varpi }_p);$$
(1)
here $`\mathrm{\Omega }`$ is the mean angular speed, $`\dot{\varpi }`$ is the apsidal precession rate, the integer $`m>0`$ is the azimuthal wavenumber, $`ϵ=0,\pm 1`$, and the unsubscripted variables refer to the disk while those with the subscript “p” refer to the planet. Resonances with $`ϵ=0`$ are called corotation resonances and those with $`ϵ=+1`$ and $`1`$ are called outer and inner Lindblad resonances, respectively. Resonances with $`ϵ=ϵ_p`$ are called “coorbital” resonances since the resonant condition (1) is satisfied when the gas orbits at the same angular speed as the planet; coorbital resonances are absent if the planet opens a gap in the disk. Resonances that are not coorbital are called external resonances. Apsidal or secular resonances have $`m=1`$, $`ϵ=ϵ_p=1`$ so the resonant condition is simply $`\dot{\varpi }=\dot{\varpi }_p`$; in this case the collective effects in the disk (mostly pressure) only have to compete with differential precession rather than with differential rotation, so the resonance is much broader.

Interactions at external Lindblad and corotation resonances excite and damp the planet’s eccentricity, respectively Goldreich & Tremaine (1980). In the absence of other effects, damping exceeds driving by a small margin ($`5\%`$). However, corotation resonances are easier to saturate, by trapping the orbital angular speed into libration around the pattern speed; if the corotation resonances are saturated and the Lindblad resonances are not, then the eccentricity can grow.

Interactions at coorbital resonances damp the planet’s eccentricity Ward (1986); Artymowicz (1993), but these only operate in a gap-free disk. Apsidal resonances can also damp eccentricity Ward & Hahn (1998, 2000), but at a rate that depends sensitively on the disk properties Goldreich & Sari (2003).

Given these complexities, about all we can conclude is that either eccentricity growth or damping may occur, depending on the properties of the planet and the disk. Numerical simulations have so far provided only limited insight: two-dimensional simulations Nelson et al. (2000); Papaloizou et al. (2001) show eccentricity growth if and only if the planet mass exceeds 10–$`20M_J`$, but the relation of these results to the resonant behavior described above remains unclear.

### 2.2 2.2. Close encounters between planets

In the final stages of planet formation dynamical instability can develop, either as the masses of the planets increase due to accretion or as their orbital separation decreases due to differential migration. The instability usually leads one or more planets to be ejected from the system or to collide with one another or with the central star (this last outcome is relatively rare). The surviving planets are left on eccentric orbits Rasio & Ford (1996); Weidenschilling & Marzari (1996); Lin & Ida (1997); Ford et al. (2001).

A great attraction of this mechanism is that it makes calculable predictions; unfortunately, the predictions have some difficulty matching the observations: (i) It has been suggested that this mechanism could produce the ‘hot Jupiters’—planets such as 51 Peg B that are found on low-eccentricity, very short-period orbits (a few days) Rasio & Ford (1996)—by close encounters of planets at distances of a few AU, which throw one planet onto a highly eccentric orbit that is later circularized by tidal dissipation. However, the frequency of hot Jupiters is far larger than this mechanism could produce Ford et al. (2001). (ii) The median eccentricity of the surviving planet following an ejection is about 0.6, compared to 0.3 in the observed systems Ford et al. (2001); on the other hand, so far the simulations have focused on equal-mass planets and the eccentricity is likely to be lower when a more massive planet ejects a less massive one. (iii) Collisions between planets lead to a population of collision remnants on nearly circular orbits, which is not observed Goldreich & Sari (2003); Ford et al. (2001). To avoid collisions would require that the planets are much more compact than Jupiter (i.e. the escape speed from the planet must be much larger than the orbital speed).

Another potential problem is the predicted mass-eccentricity relation: in this scattering process, low-mass planets are naturally expected to be excited to higher eccentricities. Unfortunately, so far no simulations have looked for mass-eccentricity correlations in a large sample of simulated planetary systems, so we do not know how strong the correlation should be. The eccentricities and masses of the known extrasolar planets show a weak correlation (at about 90% confidence level) in the opposite sense: more massive planets seem to have higher eccentricities (Figure 2).

A related possibility is that the eccentricity is excited by interactions with massive planetesimals Murray et al. (1998, 2002). Normally, close encounters with small bodies tend to damp the eccentricity and resonant interactions excite eccentricity. Because this process relies on interactions with solid bodies, rather than gas giant planets, it requires a rather massive planetesimal disk, comparable to the mass of the planet, which in turn requires a gaseous protoplanetary disk that is even more massive, $`\stackrel{>}{}0.1M_{}`$.

### 2.3 2.3. Resonant interactions between planets

Many satellites in the solar system, as well as Neptune and Pluto, are locked in orbital resonances. The formation and evolution of these resonances has been thoroughly studied (see for example Peale (1999) and references therein).

Typically, resonance capture in the solar system occurs because of convergent outward migration, in which an inner satellite migrates outward faster than an outer satellite, thereby entering a mean-motion resonance. Continued migration after resonance capture excites the eccentricities of the resonant bodies; for example, this process is believed to be responsible for Pluto’s eccentricity of 0.25 Malhotra (1993).

In contrast, the known extrasolar planets have presumably migrated inward. In this case, migration after resonance capture excites eccentricities even more efficiently. For example, a test particle that is captured into the 2:1 resonance with an exterior planet will be driven to an eccentricity of unity when the planet has migrated inward by a factor of four following capture Yu & Tremaine (2001). Lee & Peale (2002) have analyzed the dynamics of the planets in GJ 876, which are locked in a 2:1 resonance; they argue that to avoid exciting the eccentricities above the observed values requires either (i) extremely strong eccentricity damping, or (ii) fine-tuning the resonance capture to occur just before migration stops.

Chiang et al. (2002) point out that eccentricities can also be excited by divergent migration, in which the ratio of the semimajor axes of two planets increases; the planets traverse a series of resonances, each of which excites additional eccentricity.

There are at least two concerns with models based on resonant interactions: (i) this mechanism requires at least two planets, whereas only a single planet has been discovered so far in most of the known extrasolar planetary systems; (ii) like close encounters, resonant interactions are expected to excite low-mass planets to higher eccentricities, but the observed mass-eccentricity correlation is small, and if anything in the opposite sense (Figure 2).

### 2.4 2.4. Secular interactions with a distant companion star

A planet’s eccentricity can be excited by secular interactions with a distant companion that does not lie in the planet’s orbital plane—this is sometimes called the Kozai mechanism Kozai (1962); Holman et al. (1997); Innanen et al. (1997); Mazeh et al. (1997); Krymolowski & Mazeh (1999); Ford et al. (2000). This mechanism has the interesting property that the separation and mass of the companion affect the period but not the amplitude of the eccentricity oscillation; thus even a weak tidal force from a distant stellar or planetary companion can excite a large eccentricity. However, the precession rate due to the companion must exceed the precession rate due to all other effects: other planets in coplanar orbits, any residual planetesimal disk, general relativity, etc.

There are several straightforward predictions for this mechanism: (i) There should be no correlation between the planet’s mass and eccentricity, which is marginally consistent with the weak correlation seen in Figure 2. (ii) High-eccentricity planets should be found in binary star systems. This does not appear to be the case (Figure 2), although two of the highest eccentricities, 16 Cyg B and HD 80686, are found in binary systems. Perhaps some of the other high eccentricities are excited by unseen companions—planets or brown dwarfs. (iii) Coplanar multi-planet systems should have low eccentricities, since their mutually induced precession is far larger than the precession due to a distant companion, thereby suppressing the Kozai mechanism. This prediction is not supported by the dozen or so multi-planet systems, which have substantial eccentricities, although in most cases we have no direct evidence that these systems are coplanar (see §3).

### 2.5 2.5. Phenomenological diffusion

Given the limited success of the models we have described so far, it is useful to ask what constraints we can place on the eccentricity excitation mechanism from simple phenomenological models. Here we explore the possibility that the excitation can be modeled as a diffusion process in phase space, with an eccentricity-dependent diffusion coefficient $`D(e)`$.

We shall work with the eccentricity vector $`𝐞=(e\mathrm{cos}\varpi ,e\mathrm{sin}\varpi )`$, since its components are approximately canonical variables for small eccentricities. Let $`n(𝐞)d𝐞`$ be the number of planets with eccentricity vector in the range $`d𝐞`$. Then the diffusion equation in the eccentricity plane can be written as

$$\frac{n(e,t)}{t}=(Dn)=\frac{1}{e}\frac{}{e}\left[eD(e)\frac{n(e,t)}{e}\right].$$
(2)
We solve this equation for two cases: (i) a single initial burst of formation at $`t=0`$ ($`n(e,0)\delta (e)/e`$); (ii) continuous formation at a constant rate ($`n(0,t)=`$const). Since planets are ejected from the system when they reach $`e=1`$, there is also a boundary condition $`n(1,t)=0`$. We assume that the diffusion coefficient is a power law, $`D(e)e^p`$.

The best-fit ccentricity distributions that emerge from the diffusion models are plotted in Figure 3. They are shown at late times when the planets in the initial-burst model have diffused away from the center of the eccentricity plane ($`e=0`$) and the continuous-formation model has reached its steady state $`n(e)e^{1p}e`$. Both the initial-burst model and the continuous-formation model produce satisfactory fits (with KS confidence level about 60%), so we cannot distinguish between the two types of models by comparing to the current observations; however, the predicted behavior is very different near $`e=0`$ so the models should be distinguishable with more data.

In both cases, the best-fit values of $`p`$ ($`p=1.9`$ for the initial-burst model and $`p=1.1`$ for the continuous-formation model) suggest that eccentricity diffusion is produced by a mechanism that is ineffective for circular orbits.

### 2.6 2.6. Propagation of eccentricity disturbances from other planets

Stars in the solar neighborhood approach passing stars to within about 400 AU during a lifetime of $`10^{10}`$ y. Such encounters cannot excite the eccentricity of a single planet on a orbit like those of the known extrasolar or solar system planets. However, disks of gas and dust around young stars are observed to extend to radii of hundreds of AU, and it is plausible to assume that these disks form planetesimals or planets that survive to the present at radii $`10^2`$ AU. If so, then passing stars can efficiently excite the eccentricities of the outer bodies in the system, and the eccentricity disturbance can then propagate inward via secular interactions between the planets or planetesimals, much like a wave. This process is described more fully in the article “Propagation of eccentricity disturbances in planetary systems” Zakamska & Tremaine (2003) elsewhere in these proceedings. The most important disk properties that determine the efficiency of this excitation mechanism are the radial extent of the disk and the steepness of its smeared-out surface-density profile; the total mass of the disk is not a factor. If the surface density is a power law, $`\mathrm{\Sigma }(r)r^q`$, efficient eccentricity excitation at small radii requires a rather flat profile, $`q1`$, compared to the minimum solar nebula profile $`q1.5`$.

### 2.7 2.7. Formation from a collapsing protostellar cloud

It is conventional to assume that the components of binary star systems form simultaneously from condensations in a collapsing protostellar cloud (and hence have high-eccentricity orbits), whereas planets form later from the protoplanetary disk (and hence initially have low-eccentricity orbits). The two populations of companions seem to be clearly separated by the brown-dwarf desert.

Thus it is remarkable that the eccentricity and period distributions of extrasolar planets and low-mass secondaries of spectroscopic binaries are almost indistinguishable Zucker & Mazeh (2001); Stepinski & Black (2001). This observation suggests that perhaps the extrasolar planets, like binary stars, formed in the collapsing protostellar cloud Black (1997); Papaloizou & Terquem (2001); Terquem & Papaloizou (2002).

Simulations by Papaloizou & Terquem (2001) suggest that this process preferentially forms either “hot Jupiters” or planets with high eccentricities and large semimajor axes; although Papaloizou & Terquem did not carry out a detailed comparison, it seems likely that, with some tuning, their simulations could reproduce the distributions of periods and eccentricities of the known extrasolar planets. A possible concern is that the lower limit for opacity-limited fragmentation is a few Jupiter masses, and many planets are known to have smaller masses. Terquem & Papaloizou (2002) postulate that the low-mass planets formed in a disk and their eccentricities were later excited by encounters with high-mass planets formed by fragmentation.

## 3 3. Inclination Distribution

So far there is little direct evidence that the planets in multi-planet systems have small mutual inclinations, as we might expect if they formed from a disk. Chiang et al. (2001) argue that if the apsidal alignment between $`\upsilon `$ And C and D ($`\mathrm{\Delta }\omega =5^{}\pm 5^{}`$) is not a coincidence, then these two planets must have mutual inclination $`\stackrel{<}{}20^{}`$. Laughlin et al. (2002) conclude that stability considerations restrict the mutual inclination between the two planets in 47 UMa to $`\stackrel{<}{}40^{}`$.

In multi-planet systems with large mutual inclinations the Kozai mechanism can lead to large inclination and eccentricity oscillations, which promote instability by enabling close encounters between planets or with the central star. Thus, even if planetary systems are formed from collapsing clouds, as suggested in §2.7, multi-planet systems may be restricted to a disk-like configuration, containing planets on both prograde and retrograde orbits.

## 4 4. Conclusions

Planetary eccentricities can be excited by a variety of mechanisms, and it is likely that more than one of these plays a role in determining the eccentricity distribution of extrasolar planets. As the planet search programs continue to find planets of lower and lower masses, and longer and longer periods, which resemble more closely the giant planets in the solar system, the issue of why our planets have much lower eccentricities than their extrasolar analogs becomes more acute.


text:
# Anisotropic simple-cubic Ising lattice: extended phenomenological renormalization-group treatment

## 1 Introduction

The phenomenological renormalization-group (RG) method in which the transfer-matrix technique and finite-size scaling (FSS) ideas are combined is a powerful tool for investigation of critical properties in different two-dimensional systems . Unfortunately, its application in three and more dimensions is sharply retarded due to huge sizes of transfer matrices which arise in approximations of $`d`$-dimensional lattices by $`L^{d1}\times \mathrm{}`$ subsystems.

Indeed, even in the simplest case of systems with only two states of a site (a spin-$`\frac{1}{2}`$ Ising model) the order of transfer matrix in three dimensions ($`d=3`$) increases according to the law $`2^{L^2}`$ (instead of the essentially more sparing law $`2^L`$ in two dimensions). Hence, for the cluster $`3\times 3\times \mathrm{}`$ it is needed to solve the eigenproblem of the transfer matrix $`512\times 512`$, for the $`4\times 4\times \mathrm{}`$ subsystem — $`\mathrm{65\hspace{0.17em}536}\times \mathrm{65\hspace{0.17em}536}`$, and for the $`5\times 5\times \mathrm{}`$ cluster it is required to find the eigenvalues and eigenvectors of dense matrices with huge sizes of $`\mathrm{33\hspace{0.17em}554\hspace{0.17em}432}`$ by $`\mathrm{33\hspace{0.17em}554\hspace{0.17em}432}`$.

One can now solve the full eigenproblem for the transfer matrices of Ising parallelepipeds $`L\times L\times \mathrm{}`$ with the side length $`L4`$. Our aim in this paper is to use such solutions with the largest effect and extract as much as possible accurate information about physical properties of the bulk system.

The ordinary phenomenological RG is based on the FSS equations for correlation lengths . However, it is known  that the phenomenological RG can be built up using other quantities with a power divergence at the phase-transition point. It is remarkable that the such modified renormalizations can provide more precise results by the same sizes of subsystems .

In this article we give improved estimates for the critical temperature of the anisotropic three-dimensional (3D) Ising model by use of the extended phenomenological RG schemes found in . Achieved estimates are taken to calculate the values of different invariants of the 3D Ising universality class. Some of them are also given in the paper.

## 2 Basic equations

Start from the ordinary FSS equations  for the inverse correlation length $`\kappa _L(t,h)`$ and singular part of the dimensionless free-energy density $`f_L^s(t,h)`$, but we write them out for the derivatives with respect to the reduced temperature $`t=(TT_c)/T_c`$ and external field $`h`$:

$$\kappa _L^{(m,n)}(t,h)=b^{my_t+ny_h1}\kappa _{L/b}^{(m,n)}(t^{},h^{})$$
(1)
and

$$f_L^{s(m,n)}(t,h)=b^{my_t+ny_hd}f_{L/b}^{s(m,n)}(t^{},h^{}).$$
(2)
Here $`\kappa _L^{(m,n)}(t,h)=^{m+n}\kappa _L/t^mh^n`$ and the same for $`f_L^{s(m,n)}`$; $`y_t`$ and $`y_h`$ are, respectively, the thermal and magnetic critical exponents of the system; $`b=L/L^{}`$ is the rescaling factor. Deriving Eqs. (1) and (2) we used a linearized form of RG equations $`t^{}b^{y_t}t`$ and $`h^{}b^{y_h}h`$.

In the traditional phenomenological RG theory , Eq. (1) with $`m=n=0`$ is considered as an RG mapping $`(t,h)(t^{},h^{})`$ for a cluster pair $`(L,L^{})`$. By this, the critical temperature $`T_c`$ is estimated from the equation

$$L\kappa _L(T_c)=L^{}\kappa _L^{}(T_c).$$
(3)
Phenomenological renormalization $`(t,h)(t^{},h^{})`$ can be also realized by using any of the relations (1) and (2) or their combination. It has been shown by the author  that some of such extended renormalizations lead to more rapid convergence in $`L`$ than the standard phenomenological RG transformation. In particular, test examples on the fully isotropic systems  exhibited that the relations

$$\frac{\kappa _L^{^{\prime \prime }}}{L^{d1}\chi _L}|_{T_c}=\frac{\kappa _L^{}^{^{\prime \prime }}}{(L^{})^{d1}\chi _L^{}}|_{T_c}$$
(4)
and

$$\frac{\chi _L^{(4)}}{L^d\chi _L^2}|_{T_c}=\frac{\chi _L^{}^{(4)}}{(L^{})^d\chi _L^{}^2}|_{T_c}$$
(5)
locate $`T_c`$ more accurately in comparison with the ordinary RG equation (3). In the relations (4) and (5), the derivative $`\kappa _L^{^{\prime \prime }}=^2\kappa _L/h^2`$, the zero-field susceptibility $`\chi _L=f_L^{s(0,2)}`$, and the nonlinear susceptibility $`\chi _L^{(4)}=f_L^{s(0,4)}`$ can be evaluated by standard formulas via the eigenvalues and eigenvectors of transfer matrices (see, e. g., ).

To get the thermal critical exponent $`y_t`$ we applied two approaches. Firstly, we used again the standard finite-size expression

$$y_t=\frac{\mathrm{ln}[L\dot{\kappa }_L/(L^{}\dot{\kappa }_L^{})]}{\mathrm{ln}(L/L^{})}$$
(6)
which follows from Eq. (1) by $`m=1`$, $`n=0`$; $`\dot{\kappa }_L=\kappa _L/t`$. Secondly, we took the formula

$$y_t=\frac{\kappa _L^{}\dot{\kappa }_L\kappa _L\dot{\kappa }_L^{}}{(\kappa _L\kappa _L^{}\dot{\kappa }_L\dot{\kappa }_L^{})^{1/2}\mathrm{ln}(L/L^{})}.$$
(7)
This expression is a direct sequence of the well-known Roomany-Wyld approximant to the Callan-Symanzik $`\beta `$-function .

To calculate the magnetic critical exponent $`y_h`$ we also used two ways:

$$y_h=\frac{d}{2}+\frac{\mathrm{ln}(\chi _L/\chi _L^{})}{2\mathrm{ln}(L/L^{})}$$
(8)
and

$$y_h=\frac{1}{2}+\frac{\mathrm{ln}(\kappa _L^{^{\prime \prime }}/\kappa _L^{}^{^{\prime \prime }})}{2\mathrm{ln}(L/L^{})}$$
(9)
\[these finite-size relations follow from Eqs. (1) and (2)\].

In addition, we calculated the universal ratios of critical FSS amplitudes. Such a kind of the ratios can be identified from the Privman-Fisher functional expressions  which for partly anisotropic systems read

$$\kappa _L(t,h)=L^1G_0𝒦(C_1tL^{y_t},C_2hL^{y_h})$$
(10)
and

$$f_L^s(t,h)=L^dG_0(C_1tL^{y_t},C_2hL^{y_h}).$$
(11)
Scaling functions $`𝒦(x_1,x_2)`$ and $`(x_1,x_2)`$ are the same within the limits of a given universality class but they may depend on the boundary conditions and the subsystem shape (a cube, infinitely long parallelepipeds, etc.). All nonuniversality including the lattice anisotropy is absorbed in the geometry prefactor $`G_0`$ and metric coefficients $`C_1`$ and $`C_2`$. The critical amplitude ratios from which the parameters $`G_0`$, $`C_1`$, and $`C_2`$ drop out should be universal. In particular, the amplitude combinations

$$U=\frac{A_{\chi ^{(4)}}A_\kappa }{A_\chi ^2}=\frac{\kappa _L\chi _L^{(4)}}{L^{d1}\chi _L^2}$$
(12)
(Binder-like ratio for the spatially anisotropic systems),

$$Y_1=\frac{A_{\kappa ^{^{\prime \prime }}}}{A_\chi }=\frac{\kappa _L^{^{\prime \prime }}}{L^{d1}\chi _L},$$
(13)
and

$$Y_2=\frac{A_{\kappa ^{(4)}}}{A_{\chi ^{(4)}}}=\frac{\kappa _L^{(4)}}{L^{d1}\chi _L^{(4)}}$$
(14)
are expected to be not depend on the lattice anisotropy parameter $`\mathrm{\Delta }=J^{}/J`$.

## 3 Results and discussion

In the present paper we carried out calculations for the subsystems $`L\times L\times \mathrm{}`$ with $`L=3`$ and 4. To avoid undesirable surface effects the periodic boundary conditions have been imposed in both transverse directions of parallelepipeds $`L\times L\times \mathrm{}`$. Thus, the transfer matrices for which the eigenproblems was solved were dense matrices of sizes up to $`\mathrm{65\hspace{0.17em}536}\times \mathrm{65\hspace{0.17em}536}`$. To solve the eigenproblem we took into account the internal and lattice symmetries of subsystems and used the block-diagonalization method (see, e. g., ). Calculations were performed on an 800 MHz Pentium III PC running the FreeBSD operating system.

Note that the matrix quasidiagonalization procedure is rather tedious in realization. Another possible way for treating the full eigenproblem of such large-scale matrices is to use supercomputers (say, Ref.  where the direct diagonalizations have been carried out for the Hamiltonian matrices of sizes $`2^{15}`$ by $`2^{15}`$).

### 3.1 Critical temperature

The critical temperature estimates coming from solutions of the transcendental equations (4) and (5) are collected in Table 1.

In the purely isotropic case $`(J^{}=J)`$ there are high precision numerical estimates for the critical point of the 3D Ising model. The most precise values for it have been obtained by Monte Carlo simulations ; $`K_c=\mathrm{0.221\hspace{0.17em}654\hspace{0.17em}59}(10)`$, i. e. $`k_BT_c/J=1/K_c=\mathrm{4.511\hspace{0.17em}5240}(21)`$.

Inspecting Table 1 one can see that the estimates for $`J^{}=J`$ which follow from Eq. (4) and (5) are the lower and upper bounds respectively. By this, their mean value has the accuracy of $`0.01\%`$. Note also that our mean estimate is better than the value $`k_BT_c/J=\mathrm{4.533\hspace{0.17em}71}`$ obtained in Ref.  (see also ) for the fully isotropic lattice by using the ordinary phenomenological renormalization of the bars with $`L=4`$ and 5.

Discuss now the anisotropic case. Here there is well-known exact asymptotic formula for the critical temperature

$$(k_BT_c/J)_{asym}=2/[\mathrm{ln}(J/2J^{})\mathrm{ln}\mathrm{ln}(J/2J^{})+O(1)]$$
(15)
as $`J^{}/J0`$. It is a direct consequence of the molecular-field approximation in which the linear Ising chain is taken as a cluster.

Unfortunately, simple formula (15) yields considerable errors in the region $`10^3J^{}/J1`$. Its modifications in spirit of Ref. , $`k_BT_c/J2/[\mathrm{ln}(J/J^{})\mathrm{ln}\mathrm{ln}(J/J^{})]`$, lead to a loss of monotonous convergence when $`J^{}/J`$ varies from unity to zero.

We choose infinitely long clusters $`L\times L\times \mathrm{}`$ stretched in a lattice direction with the dominant interaction $`J`$. Such a cluster geometry reflects the physical situation in the system. Therefore one may expect more precise results for the critical temperature as the anisotropy of the quasi-one-dimensional lattice increases. We may also expect the monotonous convergence for the estimates from Eq. (4) and (5) because there are must be physical reasons (finite length of clusters in the longitudinal direction, etc.) for the non-monotonous or oscillatory character of behavior; they are absent in our approximations. That is, if Eq. (4) yields the lower bound in the most unfavorable case $`J^{}=J`$ then it should preserve such behavior for all $`J^{}<J`$. Similar arguments are valid for the estimates following from Eq. (5); these estimates are upper.

Note that the mean values from Table 1 are better not only than the estimates of $`k_BT_c/J`$ calculated with the $`(3,4)`$ cluster pair by the standard phenomenological RG method, but than their improvements found by means of three-point extrapolations from the sizes $`L=2`$, 3, and 4 to the bulk limit .

In the range $`10^2J^{}/J1`$, there are also the data for the critical temperature of a simple-cubic Ising lattice which were extracted from the Padé-approximant analysis of the high-temperature series . For $`J^{}=J`$ according to these data, $`k_BT_c/J=4.5106`$ that is lower by $`0.014\%`$ in comparison with the results of Ref. . For $`J^{}/J=0.1`$ the authors of Ref.  found the value $`k_BT_c/J=1.343`$. This quantity overestimates somewhat the mean value from Table 1. At last, for $`J^{}/J=0.01`$ the series method  yields $`k_BT_c/J=0.65`$ that goes out of our lower bound. This is not surprising because the calculations based on the high-temperature series rapidly deteriorate owing to the very limited number ($`11`$) of terms available in such series for the anisotropic lattices.

So, we may treat the values found from Eqs. (4) and (5) as lower and upper bounds on the real critical temperature. Their mean value for each $`J^{}/J`$ yields the best estimate which we achieve in this paper for the reduced critical temperature $`k_BT_c/J`$ (the last column in Table 1). By this, its absolute error is not larger in any case than the half difference of the corresponding upper and lower bounds. Using data from Table 1 we establish that the relative errors for $`k_BT_c/J`$ monotonically decrease from $`0.72\%`$ to $`0.14\%`$ as $`J^{}/J`$ goes from 1 to $`10^3`$.

### 3.2 Invariants of the 3D Ising universality class

Taking the improved estimates for the critical temperature of anisotropic simple-cubic lattice we calculate now some invariants of the three-dimensional Ising model universality class.

#### 3.2.1 Critical exponents

According to the RG theory, critical exponents are determined entirely by a fixed point and do not depend on the lattice anisotropy. For a three-dimensional Ising model the universality of critical exponents has been confirmed for $`\mathrm{\Delta }[0.2,5]`$ by the high-temperature series calculations .

At present, the most precise estimates of critical exponents are provided by the high-temperature expansions for models with improved potentials characterized by suppressed leading scaling corrections. For the fully isotropic 3D Ising lattice these methods yield $`\nu =\mathrm{0.630\hspace{0.17em}12}(16)`$ and $`\gamma =1.2373(2)`$ . Hence, $`y_t=1/\nu =1.5870(4)`$ and $`y_h=(d+\gamma /\nu )/2=\mathrm{2.481\hspace{0.17em}80}(18)`$.

In Table 2 we report our estimates for the critical exponents $`y_t`$ and $`y_h`$. It follows from those data that when the lattice anisotropy parameter $`\mathrm{\Delta }`$ varies in three orders (from unity to $`10^3`$), the estimates of critical exponents are changed only on a few per cent or less. In particular, calculations by Eqs. (6) and (7) with the cluster pair $`(3,4)`$ yield respectively $`y_t=1.47(6)`$ and $`y_t=1.60(7)`$. Their variations are over the range $`44.4\%`$. Similar calculations of the magnetic critical exponent carried out by use of Eqs. (8) and (9) also with the pair $`(3,4)`$ leads to $`y_h=2.586(5)`$ and $`y_h=2.579(5)`$, correspondingly. Relative dispersions of these estimates are about $`0.2\%`$.

Thus, our calculations confirm the universality of both critical exponents in a remarkably wider range of $`\mathrm{\Delta }`$ than it was done in earlier investigations. Systematic errors of the achieved estimates arise due to small sizes, $`L`$, of subsystems used.

#### 3.2.2 Critical FSS amplitude ratios

Critical amplitudes are determined by scaling functions. As a result, their “universal ratios” like $`A_{\kappa ^{(4)}}/A_{\chi ^{(4)}}=𝒦^{(0,4)}(0,0)/^{(0,4)}(0,0)`$ depend, generally speaking, upon the lattice anisotropy because it can change the shape of subsystems. But in the case of parallelepipeds $`L^{d1}\times \mathrm{}`$ with unchanged (between themselves) transverse coupling constants the shape of a sample (all its aspect ratios) will be independent of the interaction in longitudinal direction. Such a kind of the universality is studied here.

Table 3 contains results of our calculations for the critical FSS amplitude ratios $`U=A_{\chi ^{(4)}}A_\kappa /A_\chi ^2`$, $`Y_1=A_{\kappa ^{^{\prime \prime }}}/A_\chi `$, and $`Y_2=A_{\kappa ^{(4)}}/A_{\chi ^{(4)}}`$. Calculations have been carried out for $`\mathrm{\Delta }[10^3,1]`$ by use of a cyclic cluster $`4\times 4\times \mathrm{}`$.

In accord with data of Table 3, the average ratio $`U=4.900(3)`$. Hence, when the anisotropy parameter $`\mathrm{\Delta }`$ varies in three orders, this quantity changes only on $`0.06\%`$. With such accuracy we may consider the given ratio as a constant. In the case of fully isotropic lattice, $`A_\kappa =1.26(5)`$ and $`A_{\chi ^{(4)}}/A_\chi ^2=3.9(2)`$  and therefore $`A_{\chi ^{(4)}}A_\kappa /A_\chi ^2=4.9(5)`$. Our values of $`U`$ from Table 3 are in good agreement with this estimate.

It is follow from Table 3 that $`Y_1=A_{\kappa ^{^{\prime \prime }}}/A_\chi =1.759(2)`$. Hence, the constancy of this universal amplitude ratio is estimated at least as a few times $`10^3`$. Our average value for $`Y_1`$ agrees well with the estimate for the isotropic lattice, $`A_{\kappa ^{^{\prime \prime }}}/A_\chi =1.749(6)`$ .

According to data of Table 3 the amplitude ratio $`Y_2=A_{\kappa ^{(4)}}/A_{\chi ^{(4)}}=2.0133(6)`$. Thus, this quantity is most stable out of all invariants of the 3D Ising universality class which were investigated in the paper. Note that we are not aware of any quantitative estimates for the $`A_{\kappa ^{(4)}}/A_{\chi ^{(4)}}`$.

## 4 Conclusions

In this paper the large-scale transfer-matrix computations of the phase-transition temperatures, as well as the critical exponents and critical FSS amplitudes have been carried out. Application of the extended phenomenological RG schemes allowed us to find the tight limits on the critical temperature in the anisotropic simple-cubic Ising lattice and improve the available estimates for it.

We also calculated the thermal and magnetic critical exponents. Our results confirm the universality of $`y_t`$ within $`44.4\%`$ and of $`y_h`$ within $`0.2\%`$ over a wide range of $`\mathrm{\Delta }`$, $`10^3\mathrm{\Delta }1`$.

Finally, the presented results give an obvious evidence that the critical FSS amplitude ratios $`U=A_{\chi ^{(4)}}A_\kappa /A_\chi ^2`$, $`Y_1=A_{\kappa ^{^{\prime \prime }}}/A_\chi `$, and $`Y_2=A_{\kappa ^{(4)}}/A_{\chi ^{(4)}}`$ do not depend on the lattice anisotropy parameter $`\mathrm{\Delta }=J^{}/J`$ with accuracies at least $`0.1\%`$. We give likely for the first time in the literature an estimate for the universal quantity $`Y_2`$.

## Acknowledgment

This work was supported by the Russian Foundation for Basic Research under project 03-02-16909.


text:
# Charmed-strange Mesons Experimental Results

## 1 Introduction

In a simplified picture, the charmed-strange meson $`c\overline{s}`$ (generically denoted as $`D_{sJ}`$ in this paper) is an atom of a massive charm quark and a light anti-strange quark. The mass splitting of different states is the result of interaction of the spin angular momenta of the two quarks, $`\stackrel{}{s_c}`$ and $`\stackrel{}{s_s}`$, and the orbital angular momentum $`\stackrel{}{L}`$ between them. According to HQET De Rujula et al. (1976); Isgur and Wise (1991), in the limit that the charm quark is infinitively heavy, its spin is totally decoupled from the light degree of freedom. Then the spin of charm quark $`\stackrel{}{s_c}`$ and $`\stackrel{}{j}=\stackrel{}{L}+\stackrel{}{s_s}`$ are conserved separately by strong interactions. This is the so-called heavy quark symmetry (HQS).

The charm quark, however, is not infinitively heavy, but it is heavier than the QCD scale $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. Thus taking $`\stackrel{}{J}=\stackrel{}{L}+\stackrel{}{s_s}+\stackrel{}{s_c}`$ as a good quantum number, the two ground states ($`L=0,J^P=0^{},1^{}`$) can be considered as $`j=1/2`$ doublets and the four first orbital excited states ($`L=1`$) can be treated as $`j=1/2`$ doublets ($`J^P=0^+,1^+`$) and $`j=3/2`$ doublets ($`J^P=1^+,2^+`$Isgur and Wise (1991); Godfrey and Kokoski (1991).

Before this year only four of these six states had been observed. All the observed ones are narrow. The $`0^{}`$ state, $`D_s^+`$, is the lightest $`D_{sJ}`$ meson and thus can decay only weakly Chen et al. (1983). The $`1^{}`$ state, $`D_s^+`$, was discovered in the electromagnetic radiative mode $`D_s^+D_s^+\gamma `$ Albrecht et al. (1984). The kinematically allowed strong transition $`D_s^+D_s^+\pi ^0`$ is isospin suppressed, and has branching fraction of only $``$ 6% Gronberg et al. (1995). The two observed $`L=1`$ states are $`D_{s1}(2536)^+D^{}K`$, and $`D_{sJ}(2573)^+DK`$ Albrecht et al. (1989); Kubota et al. (1994). Being members of $`j=3/2`$ doublets, they decay in D-wave not S-wave, explaining their relatively narrow widths.

The two missing $`L=1`$ states ($`0^+`$ and $`1^+`$) were predicted by most potential models Godfrey and Isgur (1985); Zeng et al. (1995); Gupta and Johnson (1995); Ebert et al. (1998); Di Pierro and Eichten (2001) to be massive enough that they would decay to $`DK`$ and $`D^{}K`$, respectively, in a S-wave. The widths were thus expected to be very broad, $``$200-300 MeV. There were, however, a few predictions that these states would have masses below $`D^{()}K`$ threshold that evidently were not paid much attention Bardeen and Hill (1994); Fayyazuddin and Riazuddin (1993); Deandrea et al. (2001). Effectively “everyone” thought that $`D^{()}K`$ were the modes to look for these two states and they were difficult to find due to the large width. The recent discoveries reveal a different picture Aubert et al. (2003a); Stone and Urheim (2003); Besson et al. (2003); Abe et al. (2003); Krokovny et al. (2003); Aubert et al. (2003b).

## 2 Discovery of $`D_{sJ}^{}(2317)^+`$

The BaBar collaboration observed a $`D_s^+\pi ^0`$ structure in their $`e^+e^{}`$ continuum event sample Aubert et al. (2003a). The center of peak is $`2317.3\pm 0.4\pm 0.8`$ MeV as shown in Figure 1. The width of the peak is $`8.6\pm 0.4`$ MeV, consistent with their detector resolution. The structure is observed in different $`D_s^+`$ decay modes. It does not appear in their generic Monte Carlo simulated sample, and thus it is not a reflection of a previously known decay.

Since the decay products of this new state must contain a charm and an anti-strange quark, it is natural to think that this is one of the $`L=1`$ $`D_{sJ}`$ mesons that are still missing. Thus it is named as $`D_{sJ}^{}(2317)^+`$. Furthermore, the $`1^+`$ meson is forbidden to decay into $`0^{}0^{}`$, whereas the $`0^+`$ meson is allowed in S-wave. The decay angular distribution is flat after reconstruction efficiency correction, which means either $`D_{sJ}^{}(2317)^+`$ is generated unpolarized or it is a spin-0 state. So this new state is probably the $`0^+`$ $`D_{sJ}`$ meson, though higher spin is not ruled out.

The mass of $`D_{sJ}^{}(2317)^+`$, however, is much lighter than the $`0^+`$ $`D_{sJ}`$ meson predicted by most potential models. For example, the model in reference Di Pierro and Eichten (2001) worked quite well with known $`D`$ and $`D_{sJ}`$ mesons at the time it was created, and successfully predicted the mass of $`0^+`$ and $`1^+D`$ mesons that were later discovered. It predicted the mass of $`0^+`$ $`D_{sJ}`$ meson to be 2487 MeV. The newly observed $`D_{sJ}^{}(2317)^+`$ is 170 MeV lower than the expectation, it is even $`40`$ MeV below the $`DK`$ threshold. And the width is much narrower ($`<10`$ MeV) than the prediction of $``$ 200-300 MeV.

## 3 Discovery of $`D_{sJ}(2460)^+`$

The CLEO collaboration confirms the $`D_s^+\pi ^0`$ resonance observed by BaBar Stone and Urheim (2003); Besson et al. (2003). They find that the measured width of the peak is $`8.0_{1.2}^{+1.3}`$ MeV, somewhat broader than their detector resolution of $`6.0\pm 0.3`$ MeV. More interestingly they also observe another state, $`D_{sJ}(2460)^+`$, at 2463 MeV that decays into $`D_s^+\pi ^0`$ (Figure 2).

Figure 2.a shows the invariant mass difference, $`\mathrm{\Delta }M=M(D_s^+\gamma \pi ^0)M(D_s^+\gamma )`$. Requiring $`D_s^+\gamma `$ consistent with $`D_s^+`$, they find $`55\pm 10`$ events in the peak. The center of peak is measured to be $`349.8\pm 1.3`$ MeV, similar to that of $`D_{sJ}^{}(2317)^+`$ that CLEO finds at $`349.4\pm 1.0`$ MeV in the $`\mathrm{\Delta }M=M(D_s^+\pi ^0)M(D_s^+)`$ spectrum. The width of peak is $`6.1\pm 1.0`$ MeV, close to the detector resolution of $`6.6\pm 0.5`$ MeV. The BaBar data also shows excess in $`D_s^+\gamma \pi ^0`$ invariant mass spectrum Aubert et al. (2003a), although the conclusion reached in the publication was that further study is needed due to the complexity of the reflection from the $`D_{sJ}^{}(2317)^+`$.

The $`\mathrm{\Delta }M`$ values are very close for $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$. When the $`D_s^+`$ from a $`D_{sJ}^{}(2317)^+`$ decay picks up a random photon, the invariant mass of the two can fall in the selection window of $`D_s^+`$. Because of the equality of the mass difference, when the $`\pi ^0`$ of the same $`D_{sJ}^{}(2317)^+`$ decay is added, the total invariant mass is consistent with $`D_{sJ}(2460)^+`$. Thus $`D_{sJ}^{}(2317)^+`$ could reflect into $`D_{sJ}(2460)^+`$ peak, but simulation shows that this peak has width of $`15`$ MeV, much broader than the real $`D_{sJ}(2460)^+`$ signal peak. Checking the event sample from $`D_s^+`$ sidebands (Figure 2.b), CLEO find that the reflection of $`D_{sJ}^{}(2317)^+`$ could only account for 1/5 to 1/4 of events in $`D_{sJ}(2460)^+`$ peak.

The reflection also exists in the opposite direction, when the single photon from $`D_{sJ}(2460)^+`$ decay is “ignored” and a fake $`D_{sJ}^{}(2317)^+`$ peak is created. With the MC simulation event sample CLEO estimates the cross reflection efficiencies, and then extract the true number of reconstructed $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ signals. In both peaks, about 20% events are due to reflection. The number of $`D_{sJ}(2460)^+`$ signal is $`41\pm 12`$, consistent with estimation using $`D_s^+`$ sidebands.

The Belle collaboration confirms both $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ states in continuum event sample as well as in $`B`$ decays that will be discussed in next section Abe et al. (2003); Krokovny et al. (2003). They also observed $`D_{sJ}(2460)^+`$ in $`D_s^+\gamma `$ and $`D_s^+\pi ^+\pi ^{}`$ modes. After careful study of cross reflection the BaBar collaboration also confirms the $`D_{sJ}(2460)^+`$ meson Aubert et al. (2003b). So there is no doubt about the existence of the $`D_{sJ}(2460)^+`$ state. As $`D_{sJ}(2460)^+`$ decays to $`1^{}0^{}`$, it is most probably the missing $`J^P=1^+`$ state decays in a S-wave. It can not be a $`0^+`$ state, though other possibilities are not ruled out. Further investigation is needed.

## 4 Observation of $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ in $`B`$ decays

Cross reflection of the two new $`D_{sJ}`$ states in continuum data complicates the investigation. The cross reflection, however, is eliminated in $`B`$ decays as extra constraints are applied. Belle searches for $`B\overline{D}D_{sJ}^+`$ decays of both charged and neutral $`B`$ Krokovny et al. (2003). For events whose mass and beam energy constraints are consistent with the $`B\overline{D}D_{sJ}^+`$ decay, the invariant mass spectrum of $`D_s^+\pi ^0`$, $`D_s^+\pi ^0`$ and $`D_s^+\gamma `$ are shown in Figure 3. Belle observes both the $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ in $`B`$ decays. The peak in Figure 3.c is the first observation of $`D_{sJ}(2460)^+`$ $``$ $`D_s^+\gamma `$ mode. The ratio of partial width of this mode to that of $`D_{sJ}(2460)^+D_s^+\pi ^0`$ is measured to be $`0.38\pm 0.11\pm 0.04`$, consistent with $`0.55\pm 0.13\pm 0.8`$ measured in continuum data by Belle. The branching fractions are measured to be:

$`\begin{array}{cc}\hfill (B\overline{D}D_{sJ}^{}(2317)^+)& \times (D_{sJ}^{}(2317)^+D_s^+\pi ^0)\hfill \\ & =(8.5_{1.9}^{+2.6}\pm 2.6)\times 10^4,\hfill \\ \hfill (B\overline{D}D_{sJ}(2460)^+)& \times (D_{sJ}(2460)^+D_s^+\pi ^0)\hfill \\ & =(17.8_{3.9}^{+4.5}\pm 5.3)\times 10^4,\hfill \\ \hfill (B\overline{D}D_{sJ}(2460)^+)& \times (D_{sJ}(2460)^+D_s^+\gamma )\hfill \\ & =(6.7_{1.2}^{+1.3}\pm 2.0)\times 10^4.\hfill \end{array}`$
The $`B`$ decay provides a much better laboratory to study the spin parity of the new $`D_{sJ}`$ states. In $`B\overline{D}D_{sJ}^+`$ decay, the $`D_{sJ}`$ is totally longitudinally polarized as both $`B`$ and $`D`$ are spin-0 particles. Belle measures the helicity angular distribution of $`D_{sJ}(2460)^+`$ in $`D_s^+\gamma `$ mode shown in Figure 4. The measurement strongly supports the $`1^+`$ assignment.

## 5 Possible Explanation and Search of Other Decay Modes

The world averaged mass difference are $`349.1\pm 0.6`$ MeV and $`346.7\pm 0.8`$ MeV for $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ respectively. Adding the PDG value of M($`D_s^+`$) = $`1968.5\pm 0.6`$ MeV and M($`D_s^+`$) = $`2112.4\pm 0.7`$ MeV, the masses are $`2317.6\pm 0.8`$ MeV and $`2459.1\pm 1.0`$ MeV. The upper limits of width at 90% CL are 4.6 and 5.5 MeV respectively set by Belle Krokovny et al. (2003).

Since the discovery of $`D_{sJ}^{}(2317)^+`$ state several possible explanations appeared. Cahn and Jackson use non-relativistic vector and scalar exchange forces and recalculate within potential model to explain the mass Cahn and Jackson (2003). Van Beveran and Rupp use a unitarized meson model to explain the low mass as a threshold effect van Beveren and Rupp (2003). Bardeen et al explains that it is a normal $`c\overline{s}`$ state Bardeen and Hill (1994); Bardeen et al. (2003). Barnes et al suggest that it is a $`DK`$ molecule Barnes et al. (2003). Several others propose different multi-quark models Szczepaniak (2003); Cheng and Hou (2003); Terasaki (2003); Nussinov (2003); Datta and O’Donnell (2003).

Due to the low mass and narrow width, $`D_{sJ}^{}(2317)^+`$ has difficulty fit in the potential models, nor does $`D_{sJ}(2460)^+`$. They could be $`DK`$ and $`D^{}K`$ molecules as they are about just 40 MeV below the thresholds. The mass difference between $`D`$ and $`D^{}`$ is $``$140 MeV, explaining the mass difference between $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ of $``$142 MeV. Inside the molecule, $`D^{()}`$ and $`K`$ are pre-formed. As the direct decay mode $`D^{()}K`$ is closed, quark antiquark pairs of the two have to be broken to form a $`D_s^+`$ and a $`\pi ^0`$, thus the decay is weak.

The molecule picture suggests the existence of $`D_s^{()+}\pi ^\pm `$ resonances. Observation of these resonances would strongly support molecule hypothesis as they are not conventional $`q\overline{q}`$ meson due to their quark content. The CDF collaboration studies $`D_s^+\pi ^\pm `$ modes and find no narrow structure. The CLEO collaboration has searched for $`D_s^{()+}\pi ^\pm `$ structures as shown in Figure 5. No narrow structure is found. The productions of narrow $`D_s^{()+}\pi ^\pm `$ states are at least a factor of ten lower than the $`D_s^{()+}\pi ^0`$ modes. This proves that $`D_{sJ}^{}(2317)^+`$ and $`D_{sJ}(2460)^+`$ are iso-scalers. It, however, does not totally rule out the molecule scenario as an iso-vector molecule is expected to be broad, although there is no indication of the existence of such structure in $`B`$ decay sample.

The new $`D_{sJ}`$ states fit in well the quark model as normal $`0^+`$ and $`1^+`$ $`c\overline{s}`$ mesons except for maybe the low masses. Bardeen et al couple chiral perturbation theory with a quark model representing HQET, and in fact predicted the masses of the $`0^+`$ and $`1^+`$ $`c\overline{s}`$ mesons below $`D^{()}K`$ thresholds. The narrow widths are due to isospin violation in the decays. They infer that $`D_{sJ}^{}(2317)^+`$ is indeed the $`0^+`$ $`c\overline{s}`$ meson. It has an $`1^+`$ partner with mass splitting identical to that between $`0^{}`$ and $`1^{}`$ $`c\overline{s}`$ mesons, which is backed up by the measurements. They also calculate partial width of other decay modes as shown in Table. 1. The measured ratios and limits (at 90% C.L.) from CLEO and Belle are also listed. The predictions are consistent with the measurements, and thus this explanation is favored.

Factorization implies that the branching fractions of $`B\overline{D}D_{sJ}^+`$ for the new $`D_{sJ}`$ states be similar to that of $`D_s^+`$ and $`D_s^+`$, which are $``$ 1%. The measurements are about a factor of ten lower. This casts a shadow on the favored conventional $`c\overline{s}`$ explanation. Four-quark or molecule states, however, would have branching fraction consistent with the measurements Chen and Li (2003); Cheng and Hou (2003); Datta and O’Donnell (2003). Browder et al propose that these states are mixtures of $`c\overline{s}`$ and four-quark states Browder et al. (2003). More experimental measurements and theoretical ideas are needed to reveal the true identity of these two new states.

The author would like to thank all those, especially Prof Sheldon Stone for helping prepare this note.