text:
# Probability distribution of returns in the Heston model with stochastic volatility

## I Introduction

Stochastic dynamics of stock prices is commonly described by a geometric (multiplicative) Brownian motion, which gives a log-normal probability distribution function (PDF) for stock price changes (returns) Wilmott . However, numerous observations show that the tails of the PDF decay slower than the log-normal distribution predicts (the so-called “fat-tails” effect) Bouchaud-book ; Stanley-book ; Voit . Particularly, much attention was devoted to the power-law tails Mandelbrot ; Stanley-returns . The geometric Brownian motion model has two parameters: the drift $`\mu `$, which characterizes the average growth rate, and the volatility $`\sigma `$, which characterizes the noisiness of the process. There is empirical evidence and a set of stylized facts indicating that volatility, instead of being a constant parameter, is driven by a mean-reverting stochastic process Engle ; Papanicolaou . Various mathematical models with stochastic volatility have been discussed in literature HullWhite ; SteinStein ; Heston ; Baaquie ; Bakshi ; Duffie ; Pan .

In this paper, we study a particular stochastic volatility model, called the Heston model Heston , where the square of the stock-price volatility, called the variance $`v`$, follows a random process known in financial literature as the Cox-Ingersoll-Ross process and in mathematical statistics as the Feller process Papanicolaou ; Feller . Using the Fourier and Laplace transforms Feller ; Duffie , we solve the Fokker-Planck equation for this model exactly and find the joint PDF of returns and variance as a function of time, conditional on the initial value of variance. While returns are readily known from a financial time-series data, variance is not given directly, so it acts as a hidden stochastic variable. Thus, we integrate the joint PDF over variance and obtain the marginal probability distribution function of returns unconditional on variance. The latter PDF can be directly compared with financial data. We find an excellent agreement between our results and the Dow-Jones data for the 20-years period of 1982–2001. Using only four fitting parameters, our equations very well reproduce the PDF of returns for time lags between 1 and 250 trading days. In contrast, in ARCH, GARCH, EGARCH, TARCH, and similar models, the number of fitting parameters can easily go to a few tens McMillan .

Our result for the PDF of returns has the form of a one-dimensional Fourier integral, which is easily calculated numerically or, in certain asymptotical limits, analytically. For large returns, we find that the PDF is exponential in log-returns, which implies a power-law distribution for returns, and we calculate the time dependence of the corresponding exponents. In the limit of long times, the PDF exhibits scaling, i.e. it becomes a function of a single combination of return and time, with the scaling function expressed in terms of a Bessel function. The Dow-Jones data follow the predicted scaling function for seven orders of magnitude.

The original paper Heston  solved the problem of option pricing for the Heston model. Numerous subsequent studies Bakshi ; Duffie ; Pan ; options  compared option pricing derived from this model and its extensions with the empirical data on option pricing. They found that the Heston model describes the empirical option prices much better than the Black-Scholes theory, and modifications of the Heston model, such as adding discontinuous jumps, further improve the agreement. However, these papers did not address the fundamental question whether the stock market actually follows the Heston stochastic process or not. Obviously, if the answer is negative, then using the Heston model for option pricing would not make much sense. The stock-market time series was studied in Ref. Pan  jointly with option prices, but the focus was just on extracting the effective parameters of the Heston model. In contrast, we present a comprehensive comparison of the stock market returns distribution with the predictions of the Heston model. Using a single set of four parameters, we fit the whole family of PDF curves for a wide variety of time lags. In order to keep the model as simple as possible with the minimal number of fitting parameters, we use the original Heston model and do not include later modifications proposed in literature, such as jumps, multiple relaxation time, etc. Bakshi ; Duffie ; Pan . Interestingly, the parameters of the model that we find from our fits of the stock market data are of the same order of magnitude as the parameters extracted from the fits of option prices in Refs. Bakshi ; Duffie ; Pan .

## II The Model

We consider a stock, whose price $`S_t`$, as a function of time $`t`$, obeys the stochastic differential equation of a geometric (multiplicative) Brownian motion in the Itô form Wilmott ; Gardiner :

$$dS_t=\mu S_tdt+\sigma _tS_tdW_t^{(1)}.$$
(1)
Here the subscript $`t`$ indicates time dependence, $`\mu `$ is the drift parameter, $`W_t^{(1)}`$ is a standard random Wiener process111The infinitesimal increments of the Wiener process $`dW_t`$ are normally-distributed (Gaussian) random variables with zero mean and the variance equal to $`dt`$., and $`\sigma _t`$ is the time-dependent volatility.

Since any solution of (1) depends only on $`\sigma _t^2`$, it is convenient to introduce the new variable $`v_t=\sigma _t^2`$, which is called the variance. We assume that $`v_t`$ obeys the following mean-reverting stochastic differential equation:

$$dv_t=\gamma (v_t\theta )dt+\kappa \sqrt{v_t}dW_t^{(2)}.$$
(2)
Here $`\theta `$ is the long-time mean of $`v`$, $`\gamma `$ is the rate of relaxation to this mean, $`W_t^{(2)}`$ is a standard Wiener process, and $`\kappa `$ is a parameter that we call the variance noise. Eq. (2) is known in financial literature as the Cox-Ingersoll-Ross (CIR) process and in mathematical statistics as the Feller process Papanicolaou ; Feller . Alternative equations for $`v_t`$, with the last term in (2) replaced by $`\kappa dW_t^{(2)}`$ or $`\kappa v_tdW_t^{(2)}`$, have been also discussed in literature HullWhite . However, in our paper, we study only the case given by Eq. (2).

We take the Wiener process appearing in (2) to be correlated with the Wiener process in (1):

$$dW_t^{(2)}=\rho dW_t^{(1)}+\sqrt{1\rho ^2}dZ_t,$$
(3)
where $`Z_t`$ is a Wiener process independent of $`W_t^{(1)}`$, and $`\rho [1,1]`$ is the correlation coefficient. A negative correlation ($`\rho <0`$) between $`W_t^{(1)}`$ and $`W_t^{(2)}`$ is known as the leverage effect (Papanicolaou, , p. 41).

It is convenient to change the variable in (1) from price $`S_t`$ to log-return $`r_t=\mathrm{ln}(S_t/S_0)`$. Using Itô’s formula Gardiner , we obtain the equation satisfied by $`r_t`$:

$$dr_t=\left(\mu \frac{v_t}{2}\right)dt+\sqrt{v_t}dW_t^{(1)}.$$
(4)
The parameter $`\mu `$ can be eliminated from (4) by changing the variable to $`x_t=r_t\mu t`$, which measures log-returns relative to the growth rate $`\mu `$:

$$dx_t=\frac{v_t}{2}dt+\sqrt{v_t}dW_t^{(1)}.$$
(5)
Where it does not cause confusion with $`r_t`$, we use the term “log-return” also for the variable $`x_t`$.

Equations (5) and (2) define a two-dimensional stochastic process for the variables $`x_t`$ and $`v_t`$ Heston ; Duffie . This process is characterized by the transition probability $`P_t(x,v|v_i)`$ to have log-return $`x`$ and variance $`v`$ at time $`t`$ given the initial log-return $`x=0`$ and variance $`v_i`$ at $`t=0`$. Time evolution of $`P_t(x,v|v_i)`$ is governed by the Fokker-Planck (or forward Kolmogorov) equation Gardiner

$`{\displaystyle \frac{}{t}}P=\gamma {\displaystyle \frac{}{v}}\left[(v\theta )P\right]+{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{x}}(vP)`$ (6)
$`+\rho \kappa {\displaystyle \frac{^2}{xv}}(vP)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{x^2}}(vP)+{\displaystyle \frac{\kappa ^2}{2}}{\displaystyle \frac{^2}{v^2}}(vP).`$
The initial condition for (6) is a product of two delta functions

$$P_{t=0}(x,v|v_i)=\delta (x)\delta (vv_i).$$
(7)
The probability distribution of the variance itself, $`\mathrm{\Pi }_t(v)=𝑑xP_t(x,v)`$, satisfies the equation

$$\frac{}{t}\mathrm{\Pi }_t(v)=\frac{}{v}\left[\gamma (v\theta )\mathrm{\Pi }_t(v)\right]+\frac{\kappa ^2}{2}\frac{^2}{^2v}\left[v\mathrm{\Pi }_t(v)\right],$$
(8)
which is obtained from (6) by integration over $`x`$. Feller Feller  has shown that this equation is well-defined on the interval $`v[0,+\mathrm{})`$ as long as $`\theta >0`$. Eq. (8) has the stationary solution

$$\mathrm{\Pi }_{}(v)=\frac{\alpha ^\alpha }{\mathrm{\Gamma }(\alpha )}\frac{v^{\alpha 1}}{\theta ^\alpha }e^{\alpha v/\theta },\alpha =\frac{2\gamma \theta }{\kappa ^2},$$
(9)
which is the Gamma distribution. The parameter $`\alpha `$ in (9) is the ratio of the average variance $`\theta `$ to the characteristic fluctuation of variance $`\kappa ^2/2\gamma `$ during the relaxation time $`1/\gamma `$. When $`\alpha \mathrm{}`$, $`\mathrm{\Pi }_{}(v)\delta (v\theta )`$. The corresponding stationary PDF of volatility $`\sigma `$ is

$$\mathrm{\Pi }_{}^{(\sigma )}(\sigma )=\frac{2\alpha ^\alpha }{\mathrm{\Gamma }(\alpha )}\frac{\sigma ^{2\alpha 1}}{\theta ^\alpha }e^{\alpha \sigma ^2/\theta }.$$
(10)
Functions (9) and (10) are integrable as long as $`\alpha >0`$. The distributions $`\mathrm{\Pi }_{}(v)`$ and $`\mathrm{\Pi }_{}^{(\sigma )}(\sigma )`$ are shown in Fig. 1 for the value $`\alpha =1.3`$ deduced from the fit of the Dow-Jones time series and given in Table 1 in Sec. VIII.

## III Solution of the Fokker-Planck equation

Since $`x`$ appears in (6) only in the derivative operator $`/x`$, it is convenient to take the Fourier transform

$$P_t(x,v|v_i)=_{\mathrm{}}^+\mathrm{}\frac{dp_x}{2\pi }e^{ip_xx}\overline{P}_{t,p_x}(v|v_i).$$
(11)
Inserting (11) into (6), we find

$`{\displaystyle \frac{}{t}}\overline{P}=\gamma {\displaystyle \frac{}{v}}\left[(v\theta )\overline{P}\right]`$ (12)
$`\left[{\displaystyle \frac{p_x^2ip_x}{2}}vi\rho \kappa p_x{\displaystyle \frac{}{v}}v{\displaystyle \frac{\kappa ^2}{2}}{\displaystyle \frac{^2}{v^2}}v\right]\overline{P}.`$
Eq. (12) is simpler than (6), because the number of variables has been reduced to two, $`v`$ and $`t`$, whereas $`p_x`$ only plays the role of a parameter.

Since Eq. (12) is linear in $`v`$ and quadratic in $`/v`$, it can be simplified by taking the Laplace transform over $`v`$

$$\stackrel{~}{P}_{t,p_x}(p_v|v_i)=_0^+\mathrm{}𝑑ve^{p_vv}\overline{P}_{t,p_x}(v|v_i).$$
(13)
The partial differential equation satisfied by $`\stackrel{~}{P}_{t,p_x}(p_v|v_i)`$ is of the first order

$$\left[\frac{}{t}+\left(\mathrm{\Gamma }p_v+\frac{\kappa ^2}{2}p_v^2\frac{p_x^2ip_x}{2}\right)\frac{}{p_v}\right]\stackrel{~}{P}=\gamma \theta p_v\stackrel{~}{P},$$
(14)
where we introduced the notation

$$\mathrm{\Gamma }=\gamma +i\rho \kappa p_x.$$
(15)
Eq. (14) has to be solved with the initial condition

$$\stackrel{~}{P}_{t=0,p_x}(p_v|v_i)=\mathrm{exp}(p_vv_i).$$
(16)
The solution of (14) is given by the method of characteristics CourantHilbert :

$$\stackrel{~}{P}_{t,p_x}(p_v|v_i)=\mathrm{exp}\left(\stackrel{~}{p}_v(0)v_i\gamma \theta _0^t𝑑\tau \stackrel{~}{p}_v(\tau )\right),$$
(17)
where the function $`\stackrel{~}{p}_v(\tau )`$ is the solution of the characteristic (ordinary) differential equation

$$\frac{d\stackrel{~}{p}_v(\tau )}{d\tau }=\mathrm{\Gamma }\stackrel{~}{p}_v(\tau )+\frac{\kappa ^2}{2}\stackrel{~}{p}_v^2(\tau )\frac{p_x^2ip_x}{2}$$
(18)
with the boundary condition $`\stackrel{~}{p}_v(t)=p_v`$ specified at $`\tau =t`$. The differential equation (18) is of the Riccati type with constant coefficients Bender , and its solution is

$$\stackrel{~}{p}_v(\tau )=\frac{2\mathrm{\Omega }}{\kappa ^2}\frac{1}{\zeta e^{\mathrm{\Omega }(t\tau )}1}\frac{\mathrm{\Gamma }\mathrm{\Omega }}{\kappa ^2},$$
(19)
where we introduced the frequency

$$\mathrm{\Omega }=\sqrt{\mathrm{\Gamma }^2+\kappa ^2(p_x^2ip_x)}.$$
(20)
and the coefficient

$$\zeta =1+\frac{2\mathrm{\Omega }}{\kappa ^2p_v+(\mathrm{\Gamma }\mathrm{\Omega })}.$$
(21)
Substituting (19) into (17), we find

$`\stackrel{~}{P}_{t,p_x}(p_v|v_i)`$ (22)
$`=\mathrm{exp}\left\{\stackrel{~}{p}_v(0)v_i+{\displaystyle \frac{\gamma \theta (\mathrm{\Gamma }\mathrm{\Omega })t}{\kappa ^2}}{\displaystyle \frac{2\gamma \theta }{\kappa ^2}}\mathrm{ln}{\displaystyle \frac{\zeta e^{\mathrm{\Omega }t}}{\zeta 1}}\right\}.`$
## IV Averaging over variance

Normally we are interested only in log-returns $`x`$ and do not care about variance $`v`$. Moreover, whereas log-returns are directly known from financial data, variance is a hidden stochastic variable that has to be estimated. Inevitably, such an estimation is done with some degree of uncertainty, which precludes a clear-cut direct comparison between $`P_t(x,v|v_i)`$ and financial data. Thus we introduce the reduced probability distribution

$$P_t(x|v_i)=\underset{0}{\overset{+\mathrm{}}{}}𝑑vP_t(x,v|v_i)=\frac{dp_x}{2\pi }e^{ip_xx}\stackrel{~}{P}_{t,p_x}(0|v_i),$$
(23)
where the hidden variable $`v`$ is integrated out, so $`p_v=0`$. Substituting $`\zeta `$ from (21) with $`p_v=0`$ into (22), we find

$`P_t(x|v_i)={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dp_x}{2\pi }}e^{ip_xxv_i\frac{p_x^2ip_x}{\mathrm{\Gamma }+\mathrm{\Omega }\mathrm{coth}(\mathrm{\Omega }t/2)}}`$
$`\times e^{\frac{2\gamma \theta }{\kappa ^2}\mathrm{ln}\left(\mathrm{cosh}\frac{\mathrm{\Omega }t}{2}+\frac{\mathrm{\Gamma }}{\mathrm{\Omega }}\mathrm{sinh}\frac{\mathrm{\Omega }t}{2}\right)+\frac{\gamma \mathrm{\Gamma }\theta t}{\kappa ^2}}.`$ (24)
To check the validity of (IV), let us consider the limiting case $`\kappa =0`$. In this case, the stochastic term in (2) is absent, so the time evolution of variance is deterministic:

$$v_t=\theta +(v_i\theta )e^{\gamma t}.$$
(25)
Then process (5) gives a Gaussian distribution for $`x`$,

$$P_t^{(\kappa =0)}(x|v_i)=\frac{1}{\sqrt{2\pi t\overline{v}_t}}\mathrm{exp}\left(\frac{(x+\overline{v}_tt/2)^2}{2\overline{v}_tt}\right),$$
(26)
with the time-averaged variance $`\overline{v}_t=\frac{1}{t}_0^t𝑑\tau v_\tau `$. On the other hand, by taking the limit $`\kappa 0`$ and integrating over $`p_x`$ in (IV), we reproduce the same expression (26).

Eq. (IV) cannot be directly compared with financial time-series data, because it depends on the unknown initial variance $`v_i`$. In order to resolve this problem, we assume that $`v_i`$ has the stationary probability distribution $`\mathrm{\Pi }_{}(v_i)`$, which is given by (9). Thus we introduce the probability distribution function $`P_t(x)`$ by averaging (IV) over $`v_i`$ with the weight $`\mathrm{\Pi }_{}(v_i)`$:

$$P_t(x)=_0^{\mathrm{}}𝑑v_i\mathrm{\Pi }_{}(v_i)P_t(x|v_i).$$
(27)
The integral over $`v_i`$ is similar to the one of the Gamma function and can be taken explicitly. The final result is the Fourier integral

$$P_t(x)=\frac{1}{2\pi }_{\mathrm{}}^+\mathrm{}𝑑p_xe^{ip_xx+F_t(p_x)}$$
(28)
with

$`F_t(p_x)={\displaystyle \frac{\gamma \theta }{\kappa ^2}}\mathrm{\Gamma }t`$ (29)
$`{\displaystyle \frac{2\gamma \theta }{\kappa ^2}}\mathrm{ln}\left[\mathrm{cosh}{\displaystyle \frac{\mathrm{\Omega }t}{2}}+{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Gamma }^2+2\gamma \mathrm{\Gamma }}{2\gamma \mathrm{\Omega }}}\mathrm{sinh}{\displaystyle \frac{\mathrm{\Omega }t}{2}}\right].`$
The variable $`p_x`$ enters (29) via the variables $`\mathrm{\Gamma }`$ from (15) and $`\mathrm{\Omega }`$ from (20). It is easy to check that $`P_t(x)`$ is real, because Re$`F`$ is an even function of $`p_x`$ and Im$`F`$ is an odd one. One can also check that $`F_t(p_x=0)=0`$, which implies that $`P_t(x)`$ is correctly normalized at all times: $`𝑑xP_t(x)=1`$. The simplified version of Eq. (29) for the case $`\rho =0`$ is given in Appendix A.

Eqs. (28) and (29) for the probability distribution $`P_t(x)`$ of log-return $`x`$ at time $`t`$ are the central analytical result of the paper. The integral in (28) can be calculated numerically or, in certain regimes discussed in Secs. V, VI, and VII, analytically. In Fig. 2, the calculated function $`P_t(x)`$, shown by solid lines, is compared with the Dow-Jones data, shown by dots. (Technical details of the data analysis are discussed in Sec. VIII.) Fig. 2 demonstrates that, with a fixed set of the parameters $`\gamma `$, $`\theta `$, $`\kappa `$, $`\mu `$, and $`\rho `$, Eqs. (28) and (29) very well reproduce the distribution of log-returns $`x`$ of the Dow-Jones index for all times $`t`$. In the log-linear scale of Fig. 2, the tails of $`\mathrm{ln}P_t(x)`$ vs. $`x`$ are straight lines, which means that that tails of $`P_t(x)`$ are exponential in $`x`$. For short times $`t`$, the distribution is narrow, and the slopes of the tails are nearly vertical. As the time $`t`$ progresses, the distribution broadens and flattens.

## V Asymptotic behavior for long time $`t`$

Eq. (2) implies that variance reverts to the equilibrium value $`\theta `$ within the characteristic relaxation time $`1/\gamma `$. In this section, we consider the asymptotic limit where time $`t`$ is much longer than the relaxation time: $`\gamma t2`$. According to (15) and (20), this condition also implies that $`\mathrm{\Omega }t2`$. Then Eq. (29) reduces to

$$F_t(p_x)\frac{\gamma \theta t}{\kappa ^2}(\mathrm{\Gamma }\mathrm{\Omega }).$$
(30)
Let us change of the variable of integration in (28) to

$$p_x=\frac{\omega _0}{\kappa \sqrt{1\rho ^2}}\stackrel{~}{p}_x+ip_0,$$
(31)
where

$$p_0=\frac{\kappa 2\rho \gamma }{2\kappa (1\rho ^2)},\omega _0=\sqrt{\gamma ^2+\kappa ^2(1\rho ^2)p_0^2}.$$
(32)
Substituting (31) into (15), (20), and (30), we transform (28) to the following form

$$P_t(x)=\frac{\omega _0e^{p_0x+\mathrm{\Lambda }t}}{\pi \kappa \sqrt{1\rho ^2}}_0^{\mathrm{}}𝑑\stackrel{~}{p}_x\mathrm{cos}(A\stackrel{~}{p}_x)e^{B\sqrt{1+\stackrel{~}{p}_x^2}},$$
(33)
where

$$A=\frac{\omega _0}{\kappa \sqrt{1\rho ^2}}\left(x+\rho \frac{\gamma \theta t}{\kappa }\right),B=\frac{\gamma \theta \omega _0t}{\kappa ^2},$$
(34)
and

$$\mathrm{\Lambda }=\frac{\gamma \theta }{2\kappa ^2}\frac{2\gamma \rho \kappa }{1\rho ^2}.$$
(35)
According to formula 3.914 from GR , the integral in (33) is equal to $`BK_1(\sqrt{A^2+B^2})/\sqrt{A^2+B^2}`$, where $`K_1`$ is the first-order modified Bessel function.

Thus, Eq. (28) in the limit $`\gamma t2`$ can be represented in the scaling form

$$P_t(x)=N_te^{p_0x}P_{}(z),P_{}(z)=K_1(z)/z,$$
(36)
where the argument $`z=\sqrt{A^2+B^2}`$ is

$$z=\frac{\omega _0}{\kappa }\sqrt{\frac{(x+\rho \gamma \theta t/\kappa )^2}{1\rho ^2}+\left(\frac{\gamma \theta t}{\kappa }\right)^2},$$
(37)
and the time-dependent normalization factor $`N_t`$ is

$$N_t=\frac{\omega _0^2\gamma \theta t}{\pi \kappa ^3\sqrt{1\rho ^2}}e^{\mathrm{\Lambda }t},$$
(38)
Eq. (36) demonstrates that, up to the factors $`N_t`$ and $`e^{p_0x}`$, the dependence of $`P_t(x)`$ on the two arguments $`x`$ and $`t`$ is given by the function $`P_{}(z)`$ of the single scaling argument $`z`$ in (37). Thus, when plotted as a function of $`z`$, the data for different $`x`$ and $`t`$ should collapse on the single universal curve $`P_{}(z)`$. This is beautifully illustrated by Fig. 3, where the Dow-Jones data for different time lags $`t`$ follows the curve $`P_{}(z)`$ for seven orders of magnitude.

In the limit $`z1`$, we can use the asymptotic expression GR  $`K_1(z)e^z\sqrt{\pi /2z}`$ in (36) and take the logarithm of $`P`$. Keeping only the leading term proportional to $`z`$ and omitting the subleading term proportional to $`\mathrm{ln}z`$, we find that $`\mathrm{ln}P_t(x)`$ has the hyperbolic distribution (Bouchaud-book, , p. 14)

$$\mathrm{ln}\frac{P_t(x)}{N_t}p_0xz\mathrm{for}z1.$$
(39)
Let us examine Eq. (39) for large and small $`|x|`$.

In the first case $`|x|\gamma \theta t/\kappa `$, Eq. (37) gives $`z\omega _0|x|/\kappa \sqrt{1\rho ^2}`$, so Eq. (39) becomes

$$\mathrm{ln}\frac{P_t(x)}{N_t}p_0x\frac{\omega _0}{\kappa \sqrt{1\rho ^2}}|x|.$$
(40)
Thus, the PDF $`P_t(x)`$ has the exponential tails (40) for large log-returns $`|x|`$. Notice that, in the considered limit $`\gamma t2`$, the slopes $`d\mathrm{ln}P/dx`$ of the exponential tails (40) do not depend on time $`t`$. Because of $`p_0`$, the slopes (40) for positive and negative $`x`$ are not equal, thus the distribution $`P_t(x)`$ is not symmetric with respect to positive and negative price changes. According to (32), this asymmetry is enhanced by a negative correlation $`\rho <0`$ between stock price and variance.

In the second case $`|x+\rho \gamma \theta t/\kappa |\gamma \theta t/\kappa `$, by Taylor-expanding $`z`$ in (37) near its minimum in $`x`$ and substituting the result into (39), we get

$$\mathrm{ln}\frac{P_t(x)}{N_t^{}}p_0x\frac{\omega _0(x+\rho \gamma \theta t/\kappa )^2}{2(1\rho ^2)\gamma \theta t},$$
(41)
where $`N_t^{}=N_t\mathrm{exp}(\omega _0\gamma \theta t/\kappa ^2)`$. Thus, for small log-returns $`|x|`$, the PDF $`P_t(x)`$ is Gaussian with the width increasing linearly in time. The maximum of $`P_t(x)`$ in (41) is achieved at

$$x_m(t)=\frac{\gamma \theta t}{2\omega _0}\left(1+2\frac{\rho (\omega _0\gamma )}{\kappa }\right).$$
(42)
Eq. (42) gives the most probable log-return $`x_m(t)`$ at time $`t`$, and the coefficient in front of $`t`$ constitutes a correction to the average growth rate $`\mu `$, so that the actual growth rate is $`\overline{\mu }=\mu +dx_m/dt`$.

As Fig. 2 illustrates, $`\mathrm{ln}P_t(x)`$ is indeed linear in $`x`$ for large $`|x|`$ and quadratic for small $`|x|`$, in agreement with (40) and (41). As time progresses, the distribution, which has the scaling form (36) and (37), broadens. Thus, the fraction $`G_t`$ of the total probability contained in the parabolic (Gaussian) portion of the curve increases, as illustrated in Fig. 4. (The procedure of calculating $`G_t`$ is explained in Appendix B.) Fig. 4 shows that, at sufficiently long times, the total probability contained in the non-Gaussian tails becomes negligible, which is known in literature Bouchaud-book . The inset in Fig. 4 illustrates that the time dependence of the probability density at maximum, $`P_t(x_m)`$, is close to $`t^{1/2}`$, which is characteristic for a Gaussian evolution.

## VI Asymptotic behavior for large log-return $`x`$

In the complex plane of $`p_x`$, function $`F(p_x)`$ becomes singular at the points $`p_x`$ where the argument of the logarithm in (29) vanishes. These points are located on the imaginary axis of $`p_x`$ and are shown by dots in Fig. 5. The singularity closest to the real axis is located on the positive (negative) imaginary axis at the point $`p_1^+`$ $`(p_1^{})`$. Because the argument of the logarithm in (29) vanishes at these two points, we can approximate $`F(p_x)`$ by the dominant, singular term: $`F(p_x)(2\gamma \theta /\kappa ^2)\mathrm{ln}(p_xp_1^\pm )`$.

For large $`|x|`$, the integrand of (28) oscillates very fast as a function of $`p_x`$. Thus, we can evaluate the integral using the method of stationary phase Bender  by shifting the contour of integration so that is passes through a saddle point of the argument $`ip_xx+F(p_x)`$ of the exponent in (28). The saddle point position $`p_s`$, shown in Fig. 5 by the cross, is determined by the equation

$$ix=\frac{dF(p_x)}{dp_x}|_{p_x=p_s}\frac{2\gamma \theta }{\kappa ^2}\times \{\begin{array}{cc}\frac{1}{p_sp_1^+},\hfill & x>0,\hfill \\ \frac{1}{p_sp_1^{}},\hfill & x<0.\hfill \end{array}$$
(43)
For a large $`|x|`$, such that $`|xp_1^\pm |2\gamma \theta /\kappa ^2`$, the saddle point $`p_s`$ is very close to the singularity point: $`p_sp_1^+`$ for $`x>0`$ and $`p_sp_1^{}`$ for $`x<0`$. Then the asymptotic expression for the probability distribution is

$$P_t(x)\{\begin{array}{cc}e^{xq_t^+},\hfill & x>0,\hfill \\ e^{xq_t^{}},\hfill & x<0,\hfill \end{array}$$
(44)
where $`q_t^\pm =ip_1^\pm (t)`$ are real and positive. Eq. (44) shows that, for all times $`t`$, the tails of the probability distribution $`P_t(x)`$ for large $`|x|`$ are exponential. The slopes of the exponential tails, $`q^\pm =d\mathrm{ln}P/dx`$, are determined by the positions $`p_1^\pm `$ of the singularities closest to the real axis.

These positions $`p_1^\pm (t)`$ and, thus, the slopes $`q_t^\pm `$ depend on time $`t`$. For times much shorter than the relaxation time ($`\gamma t2`$), the singularities lie far away from the real axis. As time increases, the singularities move along the imaginary axis toward the real axis. Finally, for times much longer than the relaxation time ($`\gamma t2`$), the singularities approach limiting points: $`p_1^\pm \pm iq_{}^\pm `$, which are shown in Fig. 5 by circled crosses. Thus, as illustrated in Fig. 6, the slopes $`q_t^\pm `$ monotonously decrease in time and saturate at long times:

$$q_t^\pm q_{}^\pm =\frac{\omega _0}{\kappa \sqrt{1\rho ^2}}\pm p_0\mathrm{for}\gamma t2.$$
(45)
The slopes (45) are in agreement with Eq. (40) valid for $`\gamma t2`$. The time dependence $`q_t^\pm `$ at short times can be also found analytically:

$$q_t^\pm \frac{2}{\kappa }\sqrt{\frac{\gamma }{t}}\mathrm{for}\gamma t2.$$
(46)
The dotted curve in Fig. 6 shows that Eq. (46) works very well for short times $`t`$, where the slope diverges at $`t0`$.

## VII Asymptotic behavior for short time $`t`$

For a short time $`t`$, we expand the equations of Sec. III to the first order in $`t`$ and set $`p_v=0`$. The last term in Eq. (22) cancels the penultimate term, and Eq. (18) gives $`\stackrel{~}{p}_v(0)=t(p_x^2ip_x)/2`$. Substituting this formula into (22) and taking the integral (23) over $`p_x`$, we find

$$P_t(x|v_i)=\frac{1}{\sqrt{2\pi v_it}}e^{\frac{(x+v_it/2)^2}{2v_it}}.$$
(47)
Eq. (47) shows that, for a short $`t`$, the probability distribution of $`x`$ evolves in a Gaussian manner with the initial variance $`v_i`$, because variance has no time to change.

Substituting (47) and (9) into (27), we find

$$P_t(x)=\frac{\alpha ^\alpha }{\mathrm{\Gamma }(\alpha )}\frac{e^{x/2}}{\sqrt{2\pi \theta t}}_0^{\mathrm{}}𝑑\stackrel{~}{v}_i\stackrel{~}{v}_i^{C1}e^{A\stackrel{~}{v}_iB/\stackrel{~}{v}_i}$$
(48)
where $`\stackrel{~}{v}_i=v_i/\theta `$, $`A=\alpha +\theta t/8`$, $`B=x^2/2\theta t`$, and $`C=\alpha 1/2`$. According to formula 3.471.9 from GR , the integral in (48) is $`2(B/A)^{C/2}K_C(2\sqrt{AB})`$ for Re$`A>0`$ and Re$`B>0`$, where $`K_C`$ is the modified Bessel function of the order $`C`$. Taking into account that $`A\alpha `$ (because $`t16\gamma /\kappa ^2`$ for short $`t`$), we obtain the final expression

$$P_t(x)=\frac{2^{1\alpha }e^{x/2}}{\mathrm{\Gamma }(\alpha )}\sqrt{\frac{\alpha }{\pi \theta t}}y^{\alpha 1/2}K_{\alpha 1/2}(y),$$
(49)
where we introduced the scaling variable

$$y=\sqrt{\frac{2\alpha x^2}{\theta t}}=\frac{2\sqrt{\gamma }}{\kappa }\frac{|x|}{\sqrt{t}}$$
(50)
In the limit $`y1`$, using the formula $`K_\nu (y)e^y\sqrt{\pi /2y}`$ in (49), we find

$$P_t(x)\frac{2^{1/2\alpha }}{\mathrm{\Gamma }(\alpha )}\sqrt{\frac{\alpha }{\theta t}}y^{\alpha 1}e^y.$$
(51)
Eqs. (50) and (51) show that the tails of the distribution are exponential in $`x`$, and the slopes $`d\mathrm{ln}P/dx`$ are in agreement with Eq. (46).

In the opposite limit $`y1`$, the small argument expansion of the Bessel function can be found from the following equations Stegun :

$$K_\nu (y)=\frac{\pi }{2}\frac{I_\nu (y)I_\nu (y)}{\mathrm{sin}(\nu \pi )},\frac{\pi }{\mathrm{sin}(\pi \nu )}=\mathrm{\Gamma }(\nu )\mathrm{\Gamma }(1\nu ),$$
(52)
and

$$I_\nu (y)\left(\frac{y}{2}\right)^\nu \underset{k=0}{\overset{\mathrm{}}{}}\frac{(y^2/4)^k}{k!\mathrm{\Gamma }(\nu +k+1)}.$$
(53)
Substituting (53) into (52), we find in the case $`1/2\alpha <3/2`$

$`K_{\alpha 1/2}(y)`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }(\alpha 1/2)}{2}}\left({\displaystyle \frac{y}{2}}\right)^{\alpha +1/2}`$ (54)
$`+{\displaystyle \frac{\mathrm{\Gamma }(\alpha +1/2)}{2}}\left({\displaystyle \frac{y}{2}}\right)^{\alpha 1/2}.`$
Substituting (54) into (49), we obtain

$$P_t(x)\frac{\mathrm{\Gamma }(\alpha 1/2)}{\mathrm{\Gamma }(\alpha )}\sqrt{\frac{\alpha }{2\pi \theta t}}\left[1\lambda \left(\frac{y}{2}\right)^{2\alpha 1}\right],$$
(55)
where we introduced the coefficient

$$\lambda =\frac{|\mathrm{\Gamma }(\alpha +1/2)|}{\mathrm{\Gamma }(\alpha 1/2)}.$$
(56)
Eq. (55) can be written in the form

$$\mathrm{ln}P_t(x)\mathrm{ln}P_t(0)\lambda \left(\frac{y}{2}\right)^{2\alpha 1}.$$
(57)
We see that $`\mathrm{ln}P_t(x)`$ approaches $`x=0`$ as a power of $`x`$ lower than 2 (for $`1/2\alpha <3/2`$). The slope $`d\mathrm{ln}P/dx`$ at $`x0`$ is zero for $`\alpha >1`$ and infinite for $`\alpha <1`$.

## VIII Comparison with the Dow-Jones time series

To test the model against financial data, we downloaded daily closing values of the Dow-Jones industrial index for the period of 20 years from 1 January 1982 to 31 December 2001 from the Web site of Yahoo Yahoo . The data set contains 5049 points, which form the time series $`\{S_\tau \}`$, where the integer time variable $`\tau `$ is the trading day number. We do not filter the data for short days, such as those before holidays.

Given $`\{S_\tau \}`$, we use the following procedure to extract the probability density $`P_t^{(DJ)}(r)`$ of log-return $`r`$ for a given time lag $`t`$. For the fixed $`t`$, we calculate the set of log-returns $`\{r_\tau =\mathrm{ln}S_{\tau +t}/S_\tau \}`$ for all possible times $`\tau `$. Then we partition the $`r`$-axis into equally spaced bins of the width $`\mathrm{\Delta }r`$ and count the number of log-returns $`r_\tau `$ belonging to each bin. In this process, we omit the bins with occupation numbers less than five, because we consider such a small statistics unreliable. Only less than 1% of the entire data set is omitted in this procedure. Dividing the occupation number of each bin by $`\mathrm{\Delta }r`$ and by the total occupation number of all bins, we obtain the probability density $`P_t^{(DJ)}(r)`$ for a given time lag $`t`$. To find $`P_t^{(DJ)}(x)`$, we replace $`rx+\mu t`$.

Assuming that the system is ergodic, so that ensemble averaging is equivalent to time averaging, we compare $`P_t^{(DJ)}(x)`$ extracted from the time-series data and $`P_t(x)`$ calculated in previous sections, which describes ensemble distribution. In the language of mathematical statistics, we compare our theoretically derived population distribution with the sample distribution extracted from the time series data. We determine parameters of the model by minimizing the mean-square deviation $`_{x,t}|\mathrm{ln}P_t^{(DJ)}(x)\mathrm{ln}P_t(x)|^2`$, where the sum is taken over all available $`x`$ and $`t=1`$, 5, 20, 40, and 250 days. These values of $`t`$ are selected because they represent different regimes: $`\gamma t1`$ for $`t=1`$ and 5 days, $`\gamma t1`$ for $`t=20`$ days, and $`\gamma t1`$ for $`t=40`$ and 250 days. As Figs. 2 and 3 illustrate, our expression (28) and (29) for the probability density $`P_t(x)`$ agrees with the data very well, not only for the selected five values of time $`t`$, but for the whole time interval from 1 to 250 trading days. However, we cannot extend this comparison to $`t`$ longer than 250 days, which is approximately 1/20 of the entire range of the data set, because we cannot reliably extract $`P_t^{(DJ)}(x)`$ from the data when $`t`$ is too long.

The values obtained for the four fitting parameters ($`\gamma `$, $`\theta `$, $`\kappa `$, $`\mu `$) are given in Table 1. We find that our fits are not very sensitive to the value of $`\rho `$, so we cannot reliably determine it. Thus, we use $`\rho =0`$ for simplicity, which gives a good fit of the data. On the other hand, a nonzero value of $`\rho `$ was found in Masoliver  by fitting the leverage correlation function introduced in Bouchaud-PRL  and in Bakshi ; Duffie ; Pan  by fitting the option prices.

All four parameters ($`\gamma `$, $`\theta `$, $`\kappa `$, $`\mu `$) shown in Table 1 have the dimensionality of 1/time. The first line of the Table gives their values in the units of 1/day, as originally determined in our fit. The second line shows the annualized values of the parameters in the units of 1/year, where we utilize the average number of 252.5 trading days per calendar year to make the conversion. The relaxation time of variance is equal to $`1/\gamma =22.2`$ trading days = 4.4 weeks $``$ 1 month, where we took into account that 1 week = 5 trading days. Thus, we find that variance has a rather long relaxation time, of the order of one month, which is in agreement with the conclusion of Ref. Masoliver .

The effective growth rate of stock prices is determined by the coordinate $`r_m(t)`$ where the probability density $`P_t(r_m)`$ is maximal. Using the relation $`r_m=x_m+\mu t`$ and Eq. (42), we find that the actual growth rate is $`\overline{\mu }=\mu \gamma \theta /2\omega _0\mu \theta /2=13`$% per year. \[Here we took into account that $`\omega _0\gamma `$, because $`\gamma \kappa /2`$ in Eq. (32).\] This number coincides with the average growth rate of the Dow-Jones index obtained by a simple fit of the time series $`\{S_\tau \}`$ with an exponential function of $`\tau `$, as shown in Fig. 7. The effective stock growth rate $`\overline{\mu }`$ is comparable with the average stock volatility after one year $`\sigma =\sqrt{\theta }=14.7\%`$. Moreover, as Fig. 1 shows, the distribution of variance is broad, and the variation of variance is comparable to its average value $`\theta `$. Thus, even though the average growth rate of stock index is positive, there is a substantial probability $`_{\mathrm{}}^0𝑑rP_t(r)=17.7`$% to have negative return for $`t=1`$ year.

According to (45), the asymmetry between the slopes of exponential tails for positive and negative $`x`$ is given by the parameter $`p_0`$, which is equal to 1/2 when $`\rho =0`$ (see also the discussion of Eq. (58) in Appendix A). The origin of this asymmetry can be traced back to the transformation from (1) to (4) using Itô’s formula. It produces the term $`0.5v_tdt`$ in the r.h.s. of (4), which then generates the second term in the r.h.s. of (6). The latter term is the only source of asymmetry in $`x`$ of $`P_t(x)`$ when $`\rho =0`$. However, in practice, the asymmetry of the slopes $`p_0=1/2`$ is quite small (about 2.7%) compared with the average slope $`q_{}^\pm \omega _0/\kappa 1/x_0=18.4`$.

By fitting the Dow-Jones data to our formula, we implicitly assumed that the parameters of the stochastic process ($`\gamma `$, $`\theta `$, $`\kappa `$, $`\mu `$) do not change in time. While this assumption may be reasonable for a limited time interval, the parameters generally could change in time. The time interval of our fit, 1982–2001, includes the crash of 1987, so one might expect that the parameters of the fit would change if we use a different interval. To verify this conjecture, in Figs. 2 and 3, we also compare the data points for the time interval 1990–2001 with the theoretical curves produced using the same values for the parameters as shown in Table 1. Although the empirical data points in the tails for long time lags decrease somewhat faster than the theory predicts, the overall agreement is quite reasonable. We find that changing the values of the fitting parameters does not visibly improve the agreement. Thus, we conclude that the parameters of the Heston stochastic process are essentially the same for 1980s and 1990s. Apparently, the crash of 1987 produced little effect on the probability distribution of returns, because the stock market quickly resumed its overall growth. On the other hand, the study Silva  indicates that the data for 2000s do not follow our theoretical curves with the same fitting parameters. The main difference appears to be in the average growth rate $`\mu `$, which became negative in 2000s, as opposed to +13% per year in 1980s and 1990s. Unfortunately, the statistics for 2000s is limited, because we have only few years. Nevertheless, it does seem to indicate that in 2000s the stock market switched to a different regime compared with 1980s and 1990s.

## IX Discussion and conclusions

We derived an analytical solution for the PDF $`P_t(x)`$ of log-returns $`x`$ as a function of time $`t`$ for the Heston model of a geometrical Brownian motion with stochastic variance. The final result has the form of a one-dimensional Fourier integral (28) and (29). (In the case $`\rho =0`$, the equations have the simpler form presented in Appendix A.) Our result agrees very well with the Dow-Jones data, as shown in Fig. 2. Comparing the theory and the data, we determine the four (non-zero) fitting parameters of the model, particularly the variance relaxation time $`1/\gamma =22.2`$ days. For time longer than $`1/\gamma `$, our theory predicts the scaling behavior (36) and (37), which the Dow-Jones data indeed exhibits over seven orders of magnitude, as shown in Fig. 3. The scaling function $`P_{}(z)=K_1(z)/z`$ is expressed in terms of the first-order modified Bessel function $`K_1`$. Previous estimates in literature of the relaxation time of volatility using various indirect indicators range from 1.5 days (Papanicolaou, , p. 80) to 73 days for the half-life of the Dow-Jones index Engle . Since we have a good fit of the entire family of PDFs for time lags from 1 to 250 trading days, we believe that our estimate, 22.2 days, is more reliable. A close value of 19.6 days was found in Ref. Masoliver .

An alternative point of view in literature is that the time evolution of volatility is not characterized by a single relaxation rate. As shown in Appendix C, the variance correlation function $`C_t^{(v)}`$ (72) in the Heston model has a simple exponential decay in time. However, the analysis of financial data (Bouchaud-book, , p. 70) indicates that the correlation function has a power-law dependence or superposition of two (or more) exponentials with the relaxation times of less than one day and more than few tens of days. (Large amount of noise in the data makes it difficult to give a precise statement.) Ref. Bacry  argues that volatility relaxation is multifractal and has no characteristic time. However, one should keep in mind that the total range (70) of variation of $`C_t^{(v)}`$ is only about 77% of its saturation value, not many orders of magnitude. As Figs. 2 of Refs. Bouchaud-EL  and Bacry  shows, the main drop of $`C_t^{(v)}`$ takes place within a reasonably well-defined and relatively short time, whereas residual relaxation is stretched over a very long time. In this situation, a simple exponential time dependence, while not exact, may account for the main part of relaxation and give a reasonable approximation for the purposes of our study. Alternatively, it is possible to generalize the Heston model by incorporating more than one relaxation times Duffie .

As Fig. 2 shows, the probability distribution $`P_t(x)`$ is exponential in $`x`$ for large $`|x|`$, where it is characterized by the time-dependent slopes $`d\mathrm{ln}P/dx`$. The theoretical analysis presented in Sec. VI shows that the slopes are determined by the singularities of the function $`F_t(p_x)`$ from (29) in the complex plane of $`p_x`$ that are closest to the real axis. The calculated time dependence of the slopes $`d\mathrm{ln}P/dx`$, shown in Fig. 6, agrees with the data very well, which further supports our statement that $`1/\gamma =22.2`$ days. Exponential tails in the probability distribution of stock log-returns have been noticed in literature before (Bouchaud-book, , p. 61), Miranda , however time dependence of the slopes has not been recognized and analyzed theoretically. As shown in Fig. 2, our equations give the parabolic dependence of $`\mathrm{ln}P_t(x)`$ on $`x`$ for small $`x`$ and linear dependence for large $`x`$, in agreement with the data. Qualitatively similar results were found in Ref. Serva  for a different model with stochastic volatility and in agreement with the NYSE index daily data. It suggests that the linear and parabolic behavior is a generic feature of the models with stochastic volatility. In Ref. Stanley-returns , the power-law dependence on $`x`$ of the tails of $`P_t(x)`$ was emphasized. However, the data for S&P 500 were analyzed in Ref. Stanley-returns  only for short time lags $`t`$, typically shorter than one day. On the other hand, our data analysis is performed for the time lags longer than one day, so the results cannot be directly compared.

Deriving $`P_t(x)`$ in Sec. IV, we assumed that variance $`v`$ has the stationary gamma-distribution $`\mathrm{\Pi }_{}(v)`$ (9). This assumption should be compared with the data. There were numerous attempts in literature to reconstruct the probability distribution of volatility from the time-series data Stanley-volatility ; Mantegna-volatility . Generally, these papers agree that the central part of the distribution is well described by a log-normal distribution, but opinions vary on the fitting of the tails. Particularly, Ref. Mantegna-volatility  performed a fit with an alternative probability distribution of volatility described in Ref. (Bouchaud-book, , p. 88). Unfortunately, none of these papers attempted to fit the data using Eq. (10), so we do not have a quantitative comparison. Taking into account that we only need the integral (27), the exact shape of $`\mathrm{\Pi }_{}(v)`$ may be not so important, and Eq. (9) may give a reasonably good approximation for our purposes, even if it does not fit the tails very precisely.

Although we tested our model for the Dow-Jones index, there is nothing specific in the model which indicates that it applies only to stock market data. It would be interesting to see how the model performs when applied to other time series, for example, the foreign exchange data Peinke , which also seem to exhibit exponential tails. The study Silva  indicates that our $`P_t(x)`$ also works very well for the S&P 500 and Nasdaq indices for 1980s and 1990s. However, in the 2000s the average growth rate $`\mu `$ of the stock market changed to a negative value, which complicates separation of fluctuations from the overall trend.

###### Acknowledgements.

We thank J.-P. Bouchaud and R. E. Prange for detailed discussions of the paper and numerous valuable comments, and A. T. Zheleznyak and A. C. Silva for careful reading of the manuscript and finding minor errors.

## Appendix A The case $`𝝆\mathbf{=}\mathrm{𝟎}`$

As explained in Sec. VIII, we fit the data using $`\rho =0`$ for simplicity. In this case, by shifting the variable of integration in (28) $`p_xp_x+i/2`$, we find

$$P_t(x)=e^{x/2}_{\mathrm{}}^+\mathrm{}\frac{dp_x}{2\pi }e^{ip_xx+F_t(p_x)},$$
(58)
where $`\alpha =2\gamma \theta /\kappa ^2`$,

$$F_t(p_x)=\frac{\alpha \gamma t}{2}\alpha \mathrm{ln}\left[\mathrm{cosh}\frac{\mathrm{\Omega }t}{2}+\frac{\mathrm{\Omega }^2+\gamma ^2}{2\gamma \mathrm{\Omega }}\mathrm{sinh}\frac{\mathrm{\Omega }t}{2}\right],$$
(59)
and

$$\mathrm{\Omega }=\sqrt{\gamma ^2+\kappa ^2(p_x^2+1/4)}\gamma \sqrt{1+p_x^2(\kappa ^2/\gamma ^2)}.$$
(60)
Now the function $`F_t(p_x)`$ is real and symmetric: $`F_t(p_x)=F_t(p_x)`$. Thus, the integral in (58) is a symmetric function of $`x`$, and the only source of asymmetry of $`P_t(x)`$ in $`x`$ is the exponential prefactor in (58), as discussed at the end of Sec. VIII.

In the second equation (60), we took into account that, according to values shown in Table 1, $`\kappa ^2/4\gamma ^21`$. Introducing the dimensionless variables

$$\stackrel{~}{t}=\gamma t,\stackrel{~}{x}=x/x_0,\stackrel{~}{p}_x=p_xx_0,x_0=\kappa /\gamma ,$$
(61)
Eqs. (58), (59), and (60) can be rewritten as follows:

$$P_t(x)=\frac{e^{x/2}}{x_0}_{\mathrm{}}^+\mathrm{}\frac{d\stackrel{~}{p}_x}{2\pi }e^{i\stackrel{~}{p}_x\stackrel{~}{x}+F_{\stackrel{~}{t}}(\stackrel{~}{p}_x)},$$
(62)
where $`\stackrel{~}{\mathrm{\Omega }}=\sqrt{1+\stackrel{~}{p}_x^2}`$ and

$$F_{\stackrel{~}{t}}(\stackrel{~}{p}_x)=\frac{\alpha \stackrel{~}{t}}{2}\alpha \mathrm{ln}\left[\mathrm{cosh}\frac{\stackrel{~}{\mathrm{\Omega }}\stackrel{~}{t}}{2}+\frac{\stackrel{~}{\mathrm{\Omega }}^2+1}{2\stackrel{~}{\mathrm{\Omega }}}\mathrm{sinh}\frac{\stackrel{~}{\mathrm{\Omega }}\stackrel{~}{t}}{2}\right].$$
(63)
It is clear from (61), (62), and (63) that the parameter $`\alpha `$ determines the shape of the function $`P_t(x)`$, whereas $`1/\gamma `$ and $`x_0`$ set the scales of $`t`$ and $`x`$.

In the limit $`\stackrel{~}{t}2`$, the scaling function (36) for $`\rho =0`$ can be written as

$$P_t(x)=N_te^{x/2}K_1(z)/z,z=\sqrt{\stackrel{~}{x}^2+\overline{t}^2},$$
(64)
where $`\overline{t}=\alpha \stackrel{~}{t}/2=t\theta /x_0^2`$ and $`N_t=\overline{t}e^{\overline{t}}/\pi x_0`$. Notice that Eq. (64) has only two fitting parameters, $`x_0`$ and $`\theta `$, whereas the general formula (62) and (63) has three fitting parameters. As follows from (39), for $`\stackrel{~}{x}\overline{t}`$ and $`\stackrel{~}{x}1`$, $`P_t(x)\mathrm{exp}(|x|/x_0)`$, so $`1/x_0`$ is the slope of the exponential tails in $`x`$.

## Appendix B Gaussian weight

Let us expand the integral in (58) for small $`x`$:

$$P_t(x)e^{x/2}\left(\mu _0\frac{1}{2}\mu _2x^2\right),$$
(65)
where the coefficients are the first and the second moments of $`\mathrm{exp}[F_t(p_x)]`$

$$\mu _0(t)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{dp_x}{2\pi }e^{F_t(p_x)},\mu _2(t)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{dp_x}{2\pi }p_x^2e^{F_t(p_x)}.$$
(66)
On the other hand, we know that $`P_t(x)`$ is Gaussian for small $`x`$. So, we can write

$$P_t(x)\mu _0e^{x/2}e^{\mu _2x^2/2\mu _0},$$
(67)
with the same coefficients as in (65). If we ignore the existence of fat tails and extrapolate (67) to $`x(\mathrm{},\mathrm{})`$, the total probability contained in such a Gaussian extrapolation will be

$$G_t=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑x\mu _0e^{x/2\mu _2x^2/2\mu _0}=\sqrt{\frac{2\pi \mu _0^3}{\mu _2}}e^{\mu _0/8\mu _2}.$$
(68)
Obviously, $`G_t<1`$, because the integral (68) does not take into account the probability contained in the fat tails. Thus, the difference $`1G_t`$ can be taken as a measure of how much the actual distribution $`P_t(x)`$ deviates from a Gaussian function.

We calculated the moments (66) numerically for the function $`F`$ given by (59), then determined the Gaussian weight $`G_t`$ from (68) and plotted it in Fig. 4 as a function of time. For $`t\mathrm{}`$, $`G_t1`$, i.e. $`P_t(x)`$ becomes Gaussian for very long time lags, which is known in literature Bouchaud-book . In the opposite limit $`t0`$, $`F_t(p_x)`$ becomes a very broad function of $`p_x`$, so we cannot calculate the moments $`\mu _0`$ and $`\mu _2`$ numerically. The singular limit $`t0`$ is studied analytically in Sec. VII.

## Appendix C Correlation function of variance

The correlation function of variance is defined as

$$C_t^{(v)}=v_{t+\tau }v_\tau =_0^{\mathrm{}}𝑑v_i_0^{\mathrm{}}𝑑vv\mathrm{\Pi }_t(v|v_i)v_i\mathrm{\Pi }_{}(v_i).$$
(69)
It depends only on the relative time $`t`$ and does not depend on the initial time $`\tau `$. The averaging $`\mathrm{}`$ is performed over the ensemble probability distribution, as written in (69), or over the initial time $`\tau `$ for time-series data. Eq. (69) has the same structure as in the influence-functional formalism of Feynman and Vernon Feynman , where $`\mathrm{\Pi }_t(v|v_i)`$ represents the conditional probability propagator from the initial value $`v_i`$ to the final value $`v`$ over the time $`t`$, and $`\mathrm{\Pi }_{}(v_i)`$ represents the stationary, equilibrium probability distribution of $`v_i`$.

Using Eqs. (69) and (9), it is easy to find the limiting values of $`C_t^{(v)}`$:

$`C_{\mathrm{}}^{(v)}=v^2=\theta ^2,`$ (70)
$`C_0^{(v)}=v^2=\theta ^2\left(1+{\displaystyle \frac{1}{\alpha }}\right)=\theta ^2(1+0.77),`$
where we used the numerical value from Table 1.

Differentiating Eq. (69) with respect to $`t`$ and using Eq. (8), we find that $`C_t^{(v)}`$ satisfies the following differential equation:

$$\frac{dC_t^{(v)}}{dt}=\gamma (C_t^{(v)}\theta ^2).$$
(71)
Thus, $`C_t^{(v)}`$ changes in time exponentially with the relaxation rate $`\gamma `$:

$$C_t^{(v)}=\theta ^2\left(1+\frac{e^{\gamma t}}{\alpha }\right).$$
(72)


format: text

text:
# Extraction of the CKM angle 𝛾 from the new ”mixed ” system of 𝐵⁺→𝜋⁺𝐾⁰ and 𝐵_𝑑⁰→𝜋⁰𝐾⁰ decays

## I  Introduction

As is well known, one of the main goals of the B-factories is to measure the CKM angles $`\alpha ,\beta `$ and $`\gamma `$. For the determination of the angle $`\gamma `$, $`B\pi K,\pi \pi `$ decay modes play a key role, and have been studied intensively in the literature . Up to now, many two-body charmless B meson decays have been observed by CLEO, BaBar and Belle collaborations. For the four $`BK\pi `$ decay modes considered here, the latest world average of the corresponding branching fractions are the following

$`(B\pi ^\pm K^{})`$ $`=`$ $`(17.3\pm 1.5)\times 10^6,`$ (1)
$`(B^\pm \pi ^0K^\pm )`$ $`=`$ $`(12.1\pm 1.7)\times 10^6,`$ (2)
$`(B^\pm \pi ^\pm K^0)`$ $`=`$ $`(17.4\pm 2.6)\times 10^6,`$ (3)
$`(B\pi ^0K^0)`$ $`=`$ $`(10.4\pm 1.7)\times 10^6,`$ (4)
The accuracy of the data is currently $`10\%`$ to $`20\%`$, and will be improved rapidly along with the progress of the experiments.

For the four $`BK\pi `$ decays, the isospin and $`SU(3)`$ flavor symmetries of strong interactions imply some important relations among their decay amplitudes . Based on these amplitude relations, three combinations of CP-averaged $`B\pi K`$ branching ratios and the corresponding ”pseudo-asymmetries” have been considered  to probe the angle $`\gamma `$:

$`\left(\begin{array}{c}R\hfill \\ A_0\hfill \end{array}\right)`$ $``$ $`{\displaystyle \frac{(B_d^0\pi ^{}K^+)\pm (\overline{B_d^0}\pi ^+K^{})}{(B^+\pi ^+K^0)+(B^{}\pi ^{}\overline{K^0})}}`$ (7)
$`\left(\begin{array}{c}R_c\hfill \\ A_0^c\hfill \end{array}\right)`$ $``$ $`2{\displaystyle \frac{(B^+\pi ^0K^+)\pm (B^{}\pi ^0K^{})}{(B^+\pi ^+K^0)+(B^{}\pi ^{}\overline{K}^0)}},`$ (10)
$`\left(\begin{array}{c}R_n\hfill \\ A_0^n\hfill \end{array}\right)`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{(B_d^0\pi ^{}K^+)\pm (\overline{B_d^0}\pi ^+K^{})}{(B_d^0\pi ^0K^0)+(\overline{B_d^0}\pi ^0\overline{K^0})}},`$ (13)
where the factors of $`2`$ and $`1/2`$ have been introduced to absorb the $`\sqrt{2}`$ factors originating from the wavefunctions of $`\pi ^0`$ meson. When CLEO firstly reported their observation of the decays $`B_d\pi ^\pm K^{}`$, $`B^\pm \pi ^\pm K`$, the measured ratio $`R=0.65\pm 0.40`$ lead to an interesting bound on angle $`\gamma `$. Since the measured $`R`$ is now very close to 1, however, the constraint on the angle $`\gamma `$ from this ”mixed” system is becoming weak now. For the possible constraints on $`\gamma `$ derived from the ”charged” and ”neutral” systems, one can see for example Refs. and references therein.

In this paper, we define and study a new ”mixed” system $`B^+\pi ^+K^0`$ and $`B_d^0\pi ^0K^0`$, to see if we can extract out or put some constraints on the angle $`\gamma `$ from the new observables, the ratio $`R_1`$ and the corresponding ”pseudo-asymmetry” $`A_1`$,

$`\left(\begin{array}{c}R_1\hfill \\ A_1\hfill \end{array}\right)`$ $``$ $`2{\displaystyle \frac{(B_d^0\pi ^0K^0)\pm (\overline{B}_d^0\pi ^0\overline{K}^0)}{(B^+\pi ^+K^0)+(B^{}\pi ^{}\overline{K}^0)}}`$ (16)
We will also make an update for the constraint on the angle $`\gamma `$ from the observables $`R`$ and $`A_0`$.

Using the CP-averaged branching ratios as given in Eq.(4) one finds that

$`R`$ $`=`$ $`0.99\pm 0.17,R_1=1.20\pm 0.36.`$ (17)
The central value of $`R`$ is very close to unit now.

This paper is organized as follows. In Sec.II, we present the general description of the $`B\pi K`$ decays, define the observables and make estimations about their magnitude. In Sec.III, we consider the new measured values of $`R`$ and $`A_0`$ to make an update for the bounds on $`\gamma `$ derived from the so-called ”mixed” system: $`B^+\pi ^+K^0`$ and $`B_d^0\pi ^+K^{}`$ decay modes. In Sec.IV, we study the new ”mixed” system, $`B^+\pi ^+K^0`$ and $`B_d^0\pi ^0K^0`$ decays, to find the possible bounds on $`\gamma `$ from this new combination. The conclusions are included in the final section.

## II General Description of $`B\pi K`$ decays

First of all, as illustrated in Fig.1, the Feynman diagrams contributing to the charmless $`B\pi K`$ decays can be classified as follows :

* a color-favored ”tree” amplitude $`T`$ and a color-suppressed ”tree” amplitude $`C`$;
* a QCD penguin amplitude $`P`$;
* an color-allowed electroweak (EW) penguin amplitude $`P_{EW}`$, and a color-suppressed EW penguin amplitude $`P_{EW}^C`$;
* an annihilation amplitude $`𝒜`$.

The possible rescattering diagrams are not shown in Fig.1, one can see Fig.12 of Ref. for some relevant rescattering diagrams.

Following Refs., the transition amplitudes for the four $`B\pi K`$ decays can be written as

$`A(B^+\pi ^+K^0)`$ $`=`$ $`P+c_dP_{EW}^C+𝒜,`$ (18)
$`\sqrt{2}A(B^+\pi ^0K^+)`$ $`=`$ $`\left[P+T+C+P_{EW}+c_uP_{EW}^C+𝒜\right],`$ (19)
$`A(B_d^0\pi ^{}K^+)`$ $`=`$ $`\left[P+T+c_uP_{EW}^C\right],`$ (20)
$`\sqrt{2}A(B_d^0\pi ^0K^0)`$ $`=`$ $`\left[PC\right]\left[P_{EW}c_dP_{EW}^C\right]`$ (21)
where $`c_u=2/3`$ and $`c_d=1/3`$ are the up- and down-type quark charges, respectively. Because of the small ratio $`|V_{us}V_{ub}^{}|/|V_{ts}V_{tb}^{}|0.02`$, the four $`B\pi K`$ decays are dominated by the QCD penguin $`P`$. Because of the large top quark mass, we have also to care about the EW penguins. The overall EW penguin amplitude should be $`𝒪(10\%)`$ that of the gluonic penguin $`P`$, module group-theoretic factors. The EW penguins contribute in the color-suppressed form to $`B^+\pi ^+K^0`$ and $`B_d^0\pi ^{}K^+`$ decays and are hence expected to play a minor role, whereas they contribute to $`B_d^0\pi ^0K^0`$ and $`B^+\pi ^0K^+`$ decays in the color-allowed form and may compete with tree-diagram-like topologies. Approximately, the contribution from $`P_{EW}^C`$ should be smaller than its color-allowed counterpart $`P_{EW}`$ by a factor of $`0.2`$ and is at most a $`5\%`$ effect in $`bs`$ transition relative to the dominant QCD penguin contribution.

The relative sizes of the diagrams corresponding to the $`\overline{b}\overline{u}u\overline{s}`$ and $`\overline{b}\overline{s}`$ transitions at the quark level have been estimated:

$`1`$ $`:`$ $`|P|,`$ (22)
$`𝒪(\lambda )`$ $`:`$ $`|T|,|P_{EW}|,`$ (23)
$`𝒪(\lambda ^2)`$ $`:`$ $`|C|,|P_{EW}^C|,`$ (24)
$`𝒪(\lambda ^3)`$ $`:`$ $`|A|,`$ (25)
where the parameter $`\lambda =0.22`$ is used as a measure of the approximate relative sizes of the various contributions. One can regard the above hierarchies as a simple estimation  since a modest enhancement or suppression due to hadronic matrix element for example can turn an effect of $`𝒪(\lambda ^n)`$ into an effect of $`𝒪(\lambda ^{n\pm 1})`$.

In Refs., the decay amplitudes of $`B\pi ^+K^0`$ and $`\pi ^{}K^+`$ have been parametrized as follows *** In the parameterization, the $`SU(2)`$ isospin symmetry of $`u`$ and $`d`$ quarks, the unitarity of the CKM matrix and the Wolfenstein parametrization of the CKM matrix have been applied.

$`A(B^+\pi ^+K^0)`$ $``$ $`𝒫=\left(1{\displaystyle \frac{\lambda ^2}{2}}\right)\lambda ^2A\left[1+\rho e^{i\theta }e^{i\gamma }\right]𝒫_{tc},`$ (26)
$`A(B_d^0\pi ^{}K^+)`$ $`=`$ $`\left[𝒫+T+P_{EW}^C\right],`$ (27)
with

$`𝒫_{tc}`$ $``$ $`|𝒫_{tc}|e^{i\delta _{tc}}=\left(P_tP_c\right)+\left(P_{EW}^{C(t)}P_{EW}^{C(c)}\right),`$ (28)
$`\rho e^{i\theta }`$ $`=`$ $`{\displaystyle \frac{\lambda ^2R_b}{1\lambda ^2/2}}\left[1\left({\displaystyle \frac{𝒫_{uc}+𝒜}{𝒫_{tc}}}\right)\right],`$ (29)
$`T`$ $``$ $`|T|e^{i\delta _T}e^{i\gamma },`$ (30)
$`P_{EW}^C`$ $``$ $`|P_{EW}^C|e^{i\delta _{ew}},`$ (31)
where $`P_q`$ and $`P_{EW}^{C(q)}`$ ($`q\{u,c,t\}`$) denote contributions from QCD penguin and color-suppressed electroweak penguin topologies with internal q quarks, respectively. $`𝒫_{uc}`$ is similar to $`𝒫_{tc}`$ in Eq.(28), $`\theta `$, $`\delta _{tc}`$, $`\delta _T`$ and $`\delta _{ew}`$ are CP-conserving strong phases, and

$`A`$ $`=`$ $`0.85\pm 0.04,\lambda =0.221\pm 0.002,R_b=0.38\pm 0.08`$ (32)
are the usual CKM factors.

The ratio $`R`$ and the corresponding ”pseudo-asymmetry” $`A_0`$ then take the form

$`R`$ $`=`$ $`1+r^2+ϵ^22rϵ\mathrm{cos}(\delta \mathrm{\Delta })\mathrm{cos}\gamma {\displaystyle \frac{2r}{\omega }}\left[\mathrm{cos}\delta \mathrm{cos}\gamma +\rho \mathrm{cos}(\delta \theta )\right]`$ (34)
$`+{\displaystyle \frac{2ϵ}{\omega }}\left[\mathrm{cos}\mathrm{\Delta }+\rho \mathrm{cos}(\mathrm{\Delta }\theta )\mathrm{cos}\gamma \right],`$
$`A_0`$ $`=`$ $`A_++{\displaystyle \frac{2r}{\omega }}\mathrm{sin}\delta \mathrm{sin}\gamma +{\displaystyle \frac{2ϵ}{\omega }}\rho \mathrm{sin}(\mathrm{\Delta }\theta )\mathrm{sin}\gamma +2rϵ\mathrm{sin}(\delta \mathrm{\Delta })\mathrm{sin}\gamma ,`$ (35)
where $`A_+`$ measures the direct CP violation in the decay $`B^+\pi ^+K^0`$

$`A_+`$ $``$ $`{\displaystyle \frac{BR(B^+\pi ^+K^0)BR(B^{}\pi ^{}0\overline{K^0})}{BR(B^+\pi ^+K^0)+BR(B^{}\pi ^{}\overline{K^0})}}`$ (36)
$`=`$ $`{\displaystyle \frac{2\rho }{\omega ^2}}\mathrm{sin}\theta \mathrm{sin}\gamma `$ (37)
with

$`\omega =\sqrt{1+2\rho \mathrm{cos}\theta \mathrm{cos}\gamma +\rho ^2}.`$ (38)
The parameters $`r`$ and $`ϵ`$, as well as the CP-conserving strong phases $`\delta `$ and $`\mathrm{\Delta }`$, have been defined as follows

$`r`$ $``$ $`{\displaystyle \frac{|T|}{\sqrt{<|𝒫|^2>}}},ϵ{\displaystyle \frac{|P_{EW}^C|}{\sqrt{<|𝒫|^2>}}},`$ (39)
$`\delta `$ $``$ $`\delta _T\delta _{tc},\mathrm{\Delta }\delta _{ew}\delta _{tc},`$ (40)
with

$`<|𝒫|^2>{\displaystyle \frac{1}{2}}\left(|𝒫|^2+|\overline{𝒫}|^2\right).`$ (41)
Here the $`\overline{𝒫}`$ is the CP-conjugate modes of $`𝒫`$ and obtained by performing the substitution $`\gamma \gamma `$.

For the decay $`B\pi ^0K^0`$, one parametrization presented in Ref. is

$`\sqrt{2}A(B_d^0\pi ^0K^0)`$ $``$ $`P_n=\left(1{\displaystyle \frac{\lambda ^2}{2}}\right)\lambda ^2A\left[1+\rho _ne^{i\theta _n}e^{i\gamma }\right]𝒫_{tc}^n,`$ (42)
where $`𝒫_{tc}^n`$ and $`\rho _ne^{i\theta _n}`$ take the form

$`𝒫_{tc}^n`$ $``$ $`|𝒫_{tc}^n|e^{i\delta _{tc}^n},`$ (43)
$`\rho _ne^{i\theta _n}`$ $`=`$ $`{\displaystyle \frac{\lambda ^2R_b}{1\lambda ^2}}\left[1\left({\displaystyle \frac{𝒫_{uc}𝒞}{𝒫_{tc}^n}}\right)\right]`$ (44)
By using the isospin symmetry of strong interactions for the $`u`$ and $`d`$ quarks and the decay amplitude relations as given in Eqs.(18,21), we here parametrize the $`B\pi ^0K^0`$ decay in a new way

$`\sqrt{2}A(B_d^0\pi ^0K^0)P_n=𝒫CP_{EW},`$ (45)
with

$`C|C|e^{i\delta _c}e^{i\gamma }P_{EW}|P_{EW}|e^{i\delta _{ew}^{}},`$ (46)
where the term $`C`$ denotes the contributions due to the color-suppressed ”tree” diagrams, the quantity $`P_{EW}`$ includes contributions from the color-allowed EW penguin topologies, and the $`\delta _c`$ and $`\delta _{ew}^{}`$ denote CP-conserving strong phases.

If we define the observables

$`r_1`$ $``$ $`{\displaystyle \frac{|C|}{\sqrt{<|𝒫|^2>}}},ϵ_1{\displaystyle \frac{|P_{EW}|}{\sqrt{<|𝒫|^2>}}},`$ (47)
$`\delta _1`$ $``$ $`\delta _c\delta _{tc},\mathrm{\Delta }_1\delta _{ew}^{}\delta _{tc}`$ (48)
we then find the expressions for $`R_1`$ and $`A_1`$

$`R_1`$ $`=`$ $`1+r_1^2+ϵ_1^2+{\displaystyle \frac{2r_1}{\omega }}\left[\mathrm{cos}\delta _1\mathrm{cos}\gamma +\rho \mathrm{cos}(\delta _1\theta )\right]`$ (50)
$`{\displaystyle \frac{2ϵ_1}{\omega }}\left[\mathrm{cos}\mathrm{\Delta }_1+\rho \mathrm{cos}(\mathrm{\Delta }_1\theta )\mathrm{cos}\gamma \right]2r_1ϵ_1\mathrm{cos}(\delta _1\mathrm{\Delta }_1)\mathrm{cos}\gamma `$
$`A_1`$ $`=`$ $`A_+{\displaystyle \frac{2r_1}{\omega }}\mathrm{sin}\delta _1\mathrm{sin}\gamma {\displaystyle \frac{2ϵ_1}{\omega }}\rho \mathrm{sin}(\mathrm{\Delta }_1\theta )\mathrm{sin}\gamma +2r_1ϵ_1\mathrm{sin}(\delta _1\mathrm{\Delta }_1)\mathrm{sin}\gamma ,`$ (51)
where the parameters $`\rho `$, $`\theta `$, $`\omega `$ and the CP asymmetry $`A_+`$ have been defined previously. In Eqs.(50,51), the electroweak penguin and rescattering effects are taken into account in a general way.

From the estimated hierarchy between the different diagrams as given in Eq.(25), we get to know that

$`rϵ_10.2,r_1ϵ0.04.`$ (52)
Evaluations based on the generalized factorization approach indicated that

$`r|_{fact}=0.16\pm 0.05,ϵ|_{fact}=0.010.03.`$ (53)
By direct calculations in the generalized factorization approach we find numerically that

$`r=0.140.20,ϵ=0.010.04,`$ (54)
and

$`r_1=0.0010.04,ϵ_1=0.070.15`$ (55)
in the case of neglected rescattering effects. Of course, above estimations may be affected severely by rescattering effects, which is unfortunately still unknown at present. A reliable theoretical evaluation of $`\rho `$ is indeed very difficult and requires insights into the dynamics of strong interactions. In Ref., Fleischer studied the rescattering processes of the kind

$`B^+`$ $`\{F_c^{(s)}\}\pi ^+K^0,`$ (56)
$`B^+`$ $`\{F_u^{(s)}\}\pi ^+K^0,`$ (57)
where $`F_c^{(s)}\{\overline{D^0}D_s^+,\overline{D^0}D_s^+,\mathrm{}\}`$ and $`F_u^{(s)}\{\pi ^0K^+,\pi ^0K^+,\mathrm{}\}`$, and found that (a) $`\rho 0`$ if rescattering processes of type (56) played the dominant role in $`B^+\pi ^+K^0`$ decay; (b) $`\rho =𝒪(10\%)`$ if (57) is dominant; and (c) $`\rho =𝒪(\lambda ^2R_b)0.04`$ if both (56) and (57) were similarly important.

## III Bounds on $`\gamma `$ from $`R`$ and $`A_0`$: an update

In Refs., the strategies to extract the CKM angle $`\gamma `$ from the ratio $`R`$ have been studied. Because of the refinement of the measured $`R`$ and the first measurement of $`A_+`$, we here make an update for this approach.

Very recently, CLEO, BaBar and Belle collaboration reported their first measurement about the CP violating asymmetries of $`B\pi ^+K^{},\pi ^+K^0`$ and $`\pi ^0K^+`$ decays , as listed in Table I. Although the measured CP-violating asymmetries of three $`B\pi K`$ decay modes have large uncertainty and therefore are still consistent with zero, we believe that they will be measured with a good accuracy within one or two years. From the measured $`R`$ and $`𝒜_{CP}(B\pm \pi ^\pm K^{})`$, we find that

$`A_0𝒜_{CP}(B^\pm \pi ^\pm K^{})R=0.12\pm 0.13.`$ (58)
A recent theoretical calculation based on PQCD approach predicted that $`𝒜_{CP}(B_d^0\pi ^\pm K^{})0.19`$ for $`\gamma 60^{}`$ , which is consistent with the experimental measurements. It is reasonable for us to assume that $`|A_0|0.2`$.

The ratio $`R`$ and the asymmetry $`A_0`$ as given in Eqs.(34,35) depend on seven parameters: $`r`$, $`ϵ`$, $`\gamma `$, $`\rho `$ and CP-conserving strong phases $`\theta `$, $`\delta `$ and $`\mathrm{\Delta }`$. Although parameters $`r`$ and $`ϵ`$ can be fixed through theoretical arguments, the parameter $`\rho `$ is most possibly smaller than $`0.15`$, other four parameters still remain unknown.

In Fig.2 we show the dependence of $`R`$ on the angle $`\gamma `$ for $`\rho =0.1`$, $`ϵ=0.04`$, $`r=0.2`$, while assuming $`(\theta ,\delta ,\mathrm{\Delta })=0^{}`$ (curve 1), $`90^{}`$ (curve 2) or $`180^{}`$ (curve 3). The solid curve corresponds to the standard model prediction obtained by employing the generalized factorization approach and using the input parameters as specified in Ref.. The band between two horizontal dots lines shows the experimental measurement: $`R^{exp}=0.99\pm 0.17`$. From this figure, one can see that

* The ratio $`R`$ in the generalized factorization approach has a similar dependence on $`\gamma `$ with the ratio $`R`$ as given in Eq.(34) in the case of $`(\theta ,\delta ,\mathrm{\Delta })=0^{}`$.
* The ranges of $`\gamma <45^{}`$ and $`\gamma >120^{}`$ can be excluded for extreme values ($`0^{}`$ or $`180^{}`$) of three strong phases. But such constraint on $`\gamma `$ will disappear for $`\theta =\delta =\mathrm{\Delta }=90^{}`$. In other words, no bounds on $`\gamma `$ can be extracted directly from the $`R\gamma `$ plane due to our ignorance of the strong phases.

As shown in Ref., however, the observable $`A_0`$ allow us to eliminate the strong phase $`\delta `$ in the expression of $`R`$. By assuming both the parameter $`r`$ and the strong phase $`\delta `$ in the expression of $`R`$ as ”free” parameters, one found the minimal value of $`R`$ as follows

$`R_{min}=\kappa \mathrm{sin}^2\gamma +{\displaystyle \frac{1}{\kappa }}\left({\displaystyle \frac{A_0}{2\mathrm{sin}\gamma }}\right)^2,`$ (59)
with

$`\kappa ={\displaystyle \frac{1}{\omega }}\left[1+2(ϵ\omega )\mathrm{cos}\mathrm{\Delta }+(ϵ\omega )^2\right].`$ (60)
where $`\omega `$ has been given in Eq.(38). Now $`R_{min}`$ is independent of both $`r`$ and $`\delta `$, the effects of electroweak penguin and rescattering are included through the parameters $`ϵ`$, $`\rho `$ and $`\theta `$ appeared in $`\kappa `$ and $`\omega `$. By comparing the plots of $`R_{min}\gamma `$ with the measured $`R`$ as given in Eq.(17), one may draw the bounds on the CKM angle $`\gamma `$.

As a first approximation, we neglect the rescattering and electroweak penguin effects (i.e. setting $`\rho =ϵ=0`$ ). The value of $`R_{min}`$ therefore depends on $`A_0`$ and $`\gamma `$ only. The $`\gamma `$dependence of $`R_{min}`$ for $`A_0=0,0.1,0.2,0.3`$ and $`0.4`$ is shown in Fig.3, which is identical with the Fig.1 of Ref.. From Fig.3 we get to know that if $`R`$ is found to be smaller than $`1`$, the values of $`\gamma `$ implying $`R_{min}>R`$ would be excluded. The current measured value of $`R`$ is unfortunately very close to unit, it is also unlikely to become smaller than $`1`$ when more B decay events are available. The bound on $`\gamma `$ from measurement of $`R`$ is therefore not as promising as three years ago.

For given $`A_0`$, the dependence of $`r`$ on the angle $`\gamma `$ is of the form

$`r=\sqrt{\left(R+\mathrm{cos}^2\gamma \right)\pm 2\sqrt{\mathrm{cos}^2\gamma \left(R{\displaystyle \frac{A_0^2}{4\mathrm{sin}^2\gamma }}\right)}}`$ (61)
in the case of $`\rho =ϵ=0`$. Fig.4a show such dependence for $`R=0.99`$, and $`A_0=0`$ (solid curve), $`0.1`$ (dots curve), $`0.2`$ (short-dashed curve). Fig.4b shows the same dependence of $`r`$ on the angle $`\gamma `$ but for $`R=0.65`$ as being used in Ref.. The contours as shown in Figs.(4a,4b) are rather different. For $`R=0.65`$, the value of $`r`$ can not be smaller than $`0.2`$. For $`R=0.99`$, $`r0.15`$ as indicated by theoretical calculations based on ”factorization” is natural. For given $`R`$ and $`A_0`$, the allowed ranges of $`\gamma `$ could be read off from the figures if the parameter $`r`$ can be fixed by using an additional input. For $`|A_0|=0.2`$, $`R=0.99`$ and $`r=0.2`$, for example, the regions

$`0^{}\gamma 32^{},82^{}\gamma 98^{},and150^{}\gamma 180^{}`$ (62)
should be excluded. The above constraints are consistent with the limit on $`\gamma `$ obtained from the global fit : $`43^{}<\gamma <87^{}`$ at the 95% C.L.

For more details of $`\rho `$ and $`ϵ`$ dependence of $`R_{min}`$, one can see the original papers. The new measurement of $`R`$ and $`A_0`$ do not affect previous discussions.

## IV Bounds on $`\gamma `$ from $`R_1`$ and $`A_1`$

Analogous to the cases of $`R`$ and $`A_0`$, the observables $`R_1`$ and $`A_1`$ may also lead to interesting bounds on $`\gamma `$. By comparing the expressions of $`R_1`$ and $`A_1`$ with those of $`R`$ and $`A_0`$, we find three special features

* Between the observables $`(R_1,A_1)`$ and $`(R,A_0)`$, there is a direct transformation relation: $`rr_1`$ and $`ϵϵ_1`$.
* The parameter $`r_1`$ which describes the contributions of ”color-suppressed” tree diagram is small in size: its ”factorized” value is $`(r_1)_{fact}=0.0010.04`$ and can be neglected.
* The parameter $`ϵ_1`$ which describes the contributions of ”color-allowed” electroweak penguins, however, may be large in size as given in Eq.(55), and usually can not be neglected.

Like the ratio $`R`$ and the asymmetry $`A_0`$, $`R_1`$ and $`A_1`$ also depend on seven parameters: $`r_1`$, $`ϵ_1`$, $`\rho `$, CKM angle $`\gamma `$ and the strong phases $`\theta `$, $`\delta _1`$ and $`\mathrm{\Delta }_1`$ as defined in Eqs.(29,48). The parameters $`r_1`$ and $`ϵ_1`$ can be fixed through theoretical arguments. If one can neglect or treat $`\rho `$ parameter and three strong phases properly, one may determine or put constraint on the angle $`\gamma `$ from the measured $`R_1`$.

In Fig.5 we show the general dependence of the ratio $`R_1`$ in Eq.(50) on the angle $`\gamma `$ for $`\rho =0.1`$, $`r_1=0.04`$, $`ϵ_1=0.1`$, while assuming $`(\theta ,\delta _1,\mathrm{\Delta }_1)=0^{}`$ (curve 1), $`90^{}`$ (curve 2) or $`180^{}`$ (curve 3). The solid curve corresponds to the standard model prediction of $`R_1`$ obtained by employing the generalized factorization approach and using the input parameters as specified in Ref.. The band between two horizontal dots lines shows the data: $`R_1^{exp}=1.20\pm 0.36`$. Obviously no constraint on the angle $`\gamma `$ can be obtained by comparing the measured $`R_1`$ with the theory directly. It seems that the current data prefer large strong phases (curve 3).

In case of neglected rescattering and the color-suppressed tree diagrams, the expression of $`R_1`$ in Eq.(50) can be greatly reduced into the form

$`R_1=1+ϵ_1^22ϵ_1\mathrm{cos}\mathrm{\Delta }_1.`$ (63)
Now it depends on two parameters only. If one can fix the value of $`ϵ_1`$ from theoretical arguments, the measured $`R_1`$ will imply limits on the strong phase $`\mathrm{\Delta }_1`$. In Fig.6 we show the dependence of $`R_1`$ on the strong phase $`\mathrm{\Delta }_1`$ for various values of $`ϵ_1`$. For given $`ϵ_1=0.2`$, the lower limit on $`\mathrm{\Delta }_1`$ can be read off directly from this figure

$`\mathrm{\Delta }_1115^{},`$ (64)
for $`R_1=1.2`$. In other words, the measured $`R_1`$ prefers $`\mathrm{\Delta }_190^{}`$, i.e. $`\mathrm{cos}\mathrm{\Delta }_1<0`$.

Following the same procedure of Ref., one can eliminate the strong phase $`\delta _1`$ in Eq.(50) and find the minimal value of the ratio $`R_1`$. For this purpose, we rewrite Eq.(50) and Eq.(51) as

$`R_1`$ $`=`$ $`R_0+2r_1\left(h\mathrm{cos}\delta _\mathit{1}+k\mathrm{sin}\delta _\mathit{1}\right)+r_1^2,`$ (65)
$`A`$ $`=`$ $`\left(B\mathrm{sin}\delta C\mathrm{cos}\delta \right)r_1,`$ (66)
where the quantities

$`R_0`$ $`=`$ $`12{\displaystyle \frac{ϵ_1}{\omega }}\left[\mathrm{cos}\mathrm{\Delta }+\rho \mathrm{cos}(\mathrm{\Delta }_1\theta )\mathrm{cos}\gamma \right]+ϵ_1^2,`$ (67)
$`h`$ $`=`$ $`{\displaystyle \frac{1}{\omega }}\left(\mathrm{cos}\gamma +\rho \mathrm{cos}\theta \right)ϵ_1\mathrm{cos}\mathrm{\Delta }_1\mathrm{cos}\gamma ,`$ (68)
$`k`$ $`=`$ $`{\displaystyle \frac{\rho }{\omega }}\mathrm{sin}\theta ϵ_1\mathrm{sin}\mathrm{\Delta }_1\mathrm{cos}\gamma ,`$ (69)
$`A`$ $`=`$ $`{\displaystyle \frac{A_1A_+}{2\mathrm{sin}\gamma }}+{\displaystyle \frac{ϵ_1\rho }{\omega }}\mathrm{sin}(\mathrm{\Delta }_1\theta ),`$ (70)
$`B`$ $`=`$ $`\left({\displaystyle \frac{1}{\omega }}ϵ_1\mathrm{cos}\mathrm{\Delta }_1\right),`$ (71)
$`C`$ $`=`$ $`ϵ_1\mathrm{sin}\mathrm{\Delta }_1,`$ (72)
are independent of $`r_1`$. From Eq.(66), we get

$`\mathrm{sin}\delta _1`$ $`=`$ $`{\displaystyle \frac{AB\pm C\sqrt{(B^2+C^2)r_1^2A^2}}{(B^2+C^2)r_1}},`$ (73)
$`\mathrm{cos}\delta _1`$ $`=`$ $`{\displaystyle \frac{AC\pm B\sqrt{(B^2+C^2)r_1^2A^2}}{(B^2+C^2)r_1}},`$ (74)
and then eliminate the strong phase $`\delta _1`$ in Eq.(65):

$`R_1=R_0AD\pm E\sqrt{(B^2+C^2)r_1^2A^2}+r_1^2,`$ (75)
with

$`D=2\left({\displaystyle \frac{hCkB}{B^2+C^2}}\right),E=2\left({\displaystyle \frac{hB+kC}{B^2+C^2}}\right)`$ (76)
Treating now $`r_1`$ in Eq.(75) as a free variable, we find the minimal value of $`R_1`$

$`(R_1)_{min}=t\mathrm{sin}^2\gamma +{\displaystyle \frac{1}{t}}\left({\displaystyle \frac{A_1}{2\mathrm{sin}\gamma }}\right)^2`$ (77)
with

$`t={\displaystyle \frac{1}{\omega ^2}}\left[12ϵ_1\omega \mathrm{cos}\mathrm{\Delta }_1+ϵ_1^2\omega ^2\right].`$ (78)
This is an exact formulae derived without any approximation. The effects of electroweak penguin and rescattering processes are included through the parameter $`ϵ_1`$, $`\rho `$ and $`\theta `$, respectively.

For $`B\pi ^+K^0`$ and $`\pi ^+K^{}`$ decays the EW penguin contributes in the ”color-suppressed” form only and therefore play a minor role. For $`B\pi ^0K^0`$ decay, however, the ”color allowed” electroweak penguin is important and should be taken into account. This is the main difference between two sets of observables $`(R,A_0)`$ and $`(R_1,A_1)`$. For $`\rho =0`$ and $`ϵ_10`$, the minimal value of $`R_1`$ will depend on the angle $`\gamma `$, the parameter $`ϵ_1`$, the pseudo-asymmetry $`A_1`$ and the strong phase $`\mathrm{\Delta }_1`$, as illustrated in Figs.(7,8,9).

Fig.7 shows the dependence of $`(R_1)_{min}`$ on the angle $`\gamma `$ for $`\rho =0`$, $`ϵ_1=0.2`$, $`\mathrm{\Delta }_1=180^{}`$ (maximal effect) and $`|A_1|=0`$, $`0.2`$ and $`0.4`$. From Fig.7, we find that

* Due to the contribution from the ”color allowed” electroweak penguin, the value of $`(R_1)_{min}`$ can be larger than unit now:

$`R_1t\mathrm{sin}^2\gamma `$ (79)
for $`A_1=0`$. Here the function $`t`$ is larger than unit if $`\mathrm{cos}\mathrm{\Delta }_1<0`$. The values of $`\gamma `$ implying $`(R_1)_{min}>R_1`$ would be excluded. Numerically, the ranges around $`\gamma =90^{}`$, i.e.

$`58^{}\gamma 122^{},`$ (80)
$`65^{}\gamma 115^{},`$ (81)
$`80^{}\gamma 100^{},`$ (82)
would be excluded for $`R_1=1,1.2`$ and $`1.4`$, respectively. According to current experimental measurements, $`R_11.2`$ is indeed natural, the corresponding bounds on the angle $`\gamma `$ are thus practical and interesting. Such bounds are also complimentary to the limits from global fit.
* Around $`\gamma =90^{}`$, the bound on $`\gamma `$ is approximately independent of $`A_1`$. The excluded regions around $`\gamma =0^{}`$ and $`180^{}`$, however, depend on the values of $`A_1`$.

For the minimal value of $`R_1`$, the EW penguin contribution is included through $`t=12ϵ_1\mathrm{cos}\mathrm{\Delta }_1+ϵ_1^2`$ in the case of $`\rho =0`$. Fig.8 shows the dependence of $`(R_1)_{min}`$ on the angle $`\gamma `$ for $`\rho =0`$, $`|A_1|=0.2`$, $`\mathrm{\Delta }_1=180^{}`$, $`ϵ_1=0.05,0.10,0.15`$ and $`0.20`$. It is easy to see that $`(R_1)_{min}`$ has a moderate dependence on the value of $`ϵ_1`$. Using $`ϵ_1=0.1`$ and $`R_1=1.20`$, the range of $`83^{}\gamma 97^{}`$ could be excluded.

Fig.9 shows the dependence of $`(R_1)_{min}`$ on the angle $`\gamma `$ for $`\rho =0`$, $`|A_1|=0.2`$, $`ϵ_1=0.2`$, $`\mathrm{\Delta }_1=90^{},120^{}`$, $`150^{}`$ and $`180^{}`$. Obviously, $`(R_1)_{min}`$ and thus the bound on $`\gamma `$ has a strong dependence on phase $`\mathrm{\Delta }_1`$. The bound given in Eq.(64) is the first limit on $`\mathrm{\Delta }_1`$ from the measured $`R_1`$, but its uncertainty is also large. To get a reliable bound on $`\gamma `$ from this strategy, one has to determine the value of $`\mathrm{\Delta }_1`$ with a reasonable accuracy.

Now we check the effects of the rescattering. For the minimal value of $`R_1`$, the rescattering effects are included through $`\omega =\sqrt{1+2\rho \mathrm{cos}\theta \mathrm{cos}\gamma +\rho ^2}`$, and maximum for $`\theta =0^{}`$ or $`180^{}`$. In Fig.10 we show the dependence of $`(R_1)_{min}`$ on the angle $`\gamma `$ for $`\rho =0`$, $`0.1`$, and $`0.2`$, while assuming $`ϵ_1=0.2`$, $`\mathrm{\Delta }_1=180^{}`$, $`|A_1|=0.2`$, and $`\theta =0^{}`$ or $`180^{}`$. As shown in Fig.10, $`(R_1)_{min}`$ has a weak dependence on $`\rho `$ only: its maximum around $`\gamma =90^{}`$ is almost independent of $`\rho `$. The uncertainty of the bound on the angle $`\gamma `$ for a given $`R_1`$ is at most $`10^{}`$ in the range of $`\rho =00.2`$.

According to the estimated hierarchy and direct calculation in generalized factorization approach, the parameter $`r_1`$ should be very small: $`0r_10.04`$. By treating the $`\delta _1`$ in Eq.(65) as a free parameter, on the other hand, one can also put the lower and upper bounds on $`r_1`$ from the measured $`R_1`$

$`(r_1)_{min}^{max}=\left|\sqrt{R_0t\mathrm{sin}^2\gamma }\pm \sqrt{R_1t\mathrm{sin}^2\gamma }\right|.`$ (83)
where $`R_0`$ and $`t`$ have been given in Eqs.(67,78). Fig.11 shows the allowed regions of $`r_1`$ for $`\rho =0`$, $`ϵ_1=0.2`$, $`\mathrm{\Delta }_1=180^{}`$ and for various values of $`R_1`$ corresponding to its currently allowed experimental range $`R_1=1.20\pm 0.36`$. From this figure, one can see that

* Small values of $`R_1`$ requires large values of $`r_1`$. For $`R_1=0.84`$, for example, the minimal value of $`r_1`$ is $`0.28`$, which is much larger than the theoretical estimations: $`r_10.04`$. For $`R_11.2`$, however, the minimal values of $`r_1`$ become compatible with the theoretical estimations.
* If $`r_1`$ could be fixed by using an additional input, the bounds on $`\gamma `$ can be read off directly from Fig.11. For $`R_1=1.40`$ and $`r_1=0.05`$, for example, the range of

$`78^{}\gamma 102^{}`$ (84)
could be excluded.

## V Conclusions

In this paper we defined and studied a new ”mixed” system, $`B^+\pi ^+K^0`$ and $`B_d^0\pi ^0K^0`$ decays, to extract or constrain the CKM angle $`\gamma `$ from the measured ratio $`R_1`$ and the corresponding ”pseudo-asymmetry” $`A_1`$. We also made an update for the constraints on the angle $`\gamma `$ from the observables $`R`$ and $`A_0`$.

In the parameterization, the $`SU(2)`$ isospin symmetry of strong interactions for $`u`$ and $`d`$ quarks, the unitarity of the CKM matrix and the Wolfenstein parametrization of the CKM matrix have been applied. From the theoretical calculations and currently available experimental measurements, we found the following results:

* The measured value of $`R`$ is now very close to unit, the bound on the angle $`\gamma `$ from measurement of $`R`$ is therefore not as promising as three years ago. If $`A_0`$ is measured with good accuracy and $`r`$ can be fixed by theoretical arguments, however, some bounds on $`\gamma `$ are still possible. For $`\rho =ϵ=0`$, $`R=0.99`$, $`|A_0|=0.2`$ and $`r=0.2`$, for example, the allowed regions of angle $`\gamma `$ can be read off from Fig.4

$`30^{}\gamma 82^{},and98^{}\gamma 150^{}.`$ (85)
* For $`B_d^0\pi ^0K^0`$ decay, the color-allowed EW penguin play an important role. By direct calculations in the generalized factorization approach we find that

$`r_1=0.0010.04,ϵ_1=0.070.15`$ (86)
in the case of neglected rescattering effects, which agrees well with the estimated hierarchy of different Feynman diagrams.
* As shown in Fig.6, the measured $`R_1`$ implies a limit on the strong phase $`\mathrm{\Delta }_1`$ in case of neglected rescattering and the color-suppressed tree diagrams. For given $`ϵ_1=0.2`$ and $`R_1=1.2`$, the lower limit on $`\mathrm{\Delta }_1`$ is $`\mathrm{\Delta }_1115^{}`$. This is the first limit on $`\mathrm{\Delta }_1`$ from the measured $`R_1`$, but its uncertainty is still large.
* Due to the contribution from the ”color allowed” electroweak penguin, the minimal value of $`R_1`$ can be larger than unit, as illustrated in Fig.7. Using $`ϵ_1=0.2`$ and $`R_1=1.20`$, the range of $`65^{}\gamma 115^{}`$ could be excluded, which is also approximately independent of $`A_1`$. These bounds on the angle $`\gamma `$ are interesting and complimentary to the limits from global fit.
* The bound on $`\gamma `$ has a moderate dependence on the value of $`ϵ_1`$. Using $`R_1=1.20`$ and $`ϵ_1=0.1`$, the range of $`83^{}\gamma 97^{}`$ could be excluded.
* The bound on $`\gamma `$ has a strong dependence on the strong phase $`\mathrm{\Delta }_1`$. For $`\mathrm{\Delta }_1<90^{}`$, for example, the bound may disappear. To get a reliable bound on $`\gamma `$ from this strategy, one has to determine the value of $`\mathrm{\Delta }_1`$ with a reasonable accuracy.
* The bound on $`\gamma `$ has a weak dependence on the rescattering effects only. The uncertainty of the bound on the angle $`\gamma `$ for a given $`R_1`$ is at most $`10^{}`$ in the range of $`\rho =00.2`$.
* The bounds on $`\gamma `$ can be read off directly from $`r_1\gamma `$ plane if $`r_1`$ could be fixed by using an additional input.

## ACKNOWLEDGMENTS

Z.J. Xiao acknowledges the support by the National Natural Science Foundation of China under Grant No.10075013, and by the Research Foundation of Nanjing Normal University under Grant No.2001WLXXGQA916.


text:
# Counting peaks at height 𝑘 in a Dyck path

## Abstract

A Dyck path is a lattice path in the plane integer lattice $`\times `$ consisting of steps $`(1,1)`$ and $`(1,1)`$, which never passes below the $`x`$-axis. A peak at height $`k`$ on a Dyck path is a point on the path with coordinate $`y=k`$ that is immediately preceded by a $`(1,1)`$ step and immediately followed by a $`(1,1)`$ step. In this paper we find an explicit expression for the generating function for the number of Dyck paths starting at $`(0,0)`$ and ending at $`(2n,0)`$ with exactly $`r`$ peaks at height $`k`$. This allows us to express this function via Chebyshev polynomials of the second kind and the generating function for the Catalan numbers.

Keywords: Dyck paths, Catalan numbers, Chebyshev polynomials.

## 1. Introduction and main results

The Catalan sequence is the sequence

$$\{C_n\}_{n0}=\{1,1,2,5,14,132,429,1430,\mathrm{}\},$$

where $`C_n=\frac{1}{n+1}\left(\genfrac{}{}{0pt}{}{2n}{n}\right)`$ is called the $`n`$th Catalan number. The generating function for the Catalan numbers is denoted by $`C(x)=\frac{1\sqrt{14x}}{2x}`$. The Catalan numbers provide a complete answer to the problem of counting certain properties of more than $`66`$ different combinatorial structures (see Stanley \[S, Page 219 and Exercise 6.19\]). The structure of use to us in the present paper is Dyck paths.

Chebyshev polynomials of the second kind are defined by

$$U_r(\mathrm{cos}\theta )=\frac{\mathrm{sin}(r+1)\theta }{\mathrm{sin}\theta }$$

for $`r0`$. Evidently, $`U_r(x)`$ is a polynomial of degree $`r`$ in $`x`$ with integer coefficients. Chebyshev polynomials were invented for the needs of approximation theory, but are also widely used in various other branches of mathematics, including algebra, combinatorics, number theory, and lattice paths (see \[K, Ri\]). For $`k0`$ we define $`R_k(x)`$ by

$$R_k(x)=\frac{U_{k1}\left(\frac{1}{2\sqrt{x}}\right)}{\sqrt{x}U_k\left(\frac{1}{2\sqrt{x}}\right)}.$$

For example, $`R_0(x)=0`$, $`R_1(x)=1`$, and $`R_2(x)=1/(1x)`$. It is easy to see that for any $`k`$, $`R_k(x)`$ is a rational function in $`x`$.

A Dyck path is a lattice path in the plane integer lattice $`\times `$ consisting of up-steps $`(1,1)`$ and down-steps $`(1,1)`$, which never passes below the $`x`$-axis (see Figure 1). Let $`P`$ be a Dyck path; we define the weight of $`P`$ to be the product of the weights of all its steps, where the weight of every step (up-step or down-step) is $`\sqrt{x}`$. For example, Figure 1 presents two Dyck paths, each of length $`12`$ and weight $`x^6`$.

A point on the Dyck path is called a peak at height $`k`$ if it is a point with coordinate $`y=k`$ that is immediately preceded by a up-step and immediately followed by a down-step. For example, Figure 1 presents two Dyck paths; the path on the left has two peaks at height $`2`$ and two peaks at height $`3`$; and the path on the right has one peak at height $`1`$, one peak at height $`2`$, and one peak at height $`3`$. A point on the Dyck path is called a valley at height $`k`$ if it is a point with coordinate $`y=k`$ that is immediately preceded by down-step and immediately followed by up-step. For example, in Figure 1, the path on the left has two valleys at height $`1`$ and one valley at height $`2`$, and the path on the right has only two valleys at height $`0`$. The number of all Dyck paths starting at $`(0,0)`$ and ending at $`(2n,0)`$ with exactly $`r`$ peaks (resp. valleys) at height $`k`$ we denote by $`\mathrm{peak}_k^r(n)`$ (resp. $`\mathrm{valley}_k^r(n)`$). The corresponding generating function is denoted by $`\mathrm{Peak}_k^r(x)`$ (resp. $`\mathrm{Valley}_k^r(x)`$).

Deutsch \[D\] found the number of Dyck paths of length $`2n`$ starting and ending on the $`x`$-axis with no peaks at height $`1`$ is given by the $`n`$th Fine number: $`1`$, $`0`$, $`1`$, $`2`$, $`6`$, $`18`$, $`57`$, $`\mathrm{}`$ (see \[D, DS, F\] and \[SP, Sequence M1624\]). Recently, Peart and Woan \[PW\] gave a complete answer for the number of Dyck paths of length $`2n`$ starting and ending on the $`x`$-axis with no peaks at height $`k`$. This result can be formulated as follows.

###### Theorem 1.1.

(see \[PW, Section 2\]) The generating function for the number of Dyck paths of length $`2n`$ starting and ending on the $`x`$-axis with no peaks at height $`k`$ is given by

$$\frac{1}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1{\displaystyle \frac{\mathrm{}}{1{\displaystyle \frac{x}{1x^2C^2(x)}}}}}}}}},$$

where the continued fraction contains exactly $`k`$ levels.

Theorem 1.1 is in fact a simple consequence of Theorem 1.2 (as we are going to show in Section 3).

###### Theorem 1.2.

(see \[RV, Proposition 1\]) For given a Dyck path $`P`$ we give every up-step the weight $`1`$, every down-step from height $`k`$ to height $`k1`$ not following a peak the weight $`\lambda _k`$, and every down-step following a peak of height $`k`$ the weight $`\mu _k`$. The weight of $`w(P)`$ of the path $`P`$ is the product of the weights of its steps. Then the generating function $`_Pw(P)`$, where the sum over all the Dyck paths, is given by

$$\frac{1}{1(\mu _1\lambda _1){\displaystyle \frac{\lambda _1}{1(\mu _2\lambda _2){\displaystyle \frac{\lambda _2}{1(\mu _3\lambda _3)\mathrm{}}}}}}$$

In this paper we find an explicit formulas for the generating functions $`\mathrm{Peak}_k^r(x)`$ and $`\mathrm{Valley}_k^r(x)`$ for any $`k,r0`$. This allows us to express these functions via Chebyshev polynomials of the second kind $`U_k(x)`$ and generating function for the Catalan numbers $`C(x)`$. The main result of this paper can be formulated as follows:

###### Main Theorem 1.1.

* For all $`k2`$,

$$\mathrm{Peak}_k^r(x)=\mathrm{Valley}_{k2}^r(x);$$
* For all $`k,r0`$,

$$\mathrm{Valley}_k^r(x)=\delta _{r,0}R_{k+1}(x)+\frac{x^rC^{r+1}(x)}{U_{k+1}^2\left(\frac{1}{2\sqrt{x}}\right)\left(1x(R_{k+1}(x)1)C(x)\right)^{r+1}};$$
* For all $`r0`$,

$$\mathrm{Peak}_1^r(x)=\delta _{r,0}+\frac{x^{3r+2}C^{2r+2}(x)}{(1x^2C^2(x))^{r+1}}.$$

We give two proofs of this result. The first proof, given in Section 2, uses a decomposition of the paths under consideration, while the second proof, given in Section 3, uses the continued fraction theorem due to Roblet and Viennot (see Theorem 1.2) as the starting point.

###### Remark 1.3.

By the first part and the second part of the Main Theorem, we obtain an explicit expression for the generating function for the number of Dyck paths starting at $`(0,0)`$ and ending on the $`x`$-axis with no peaks at height $`k2`$, namely

$$\mathrm{Peak}_k^0(x)=R_{k1}(x)+\frac{x^rC^{r+1}(x)}{U_{k1}^2\left(\frac{1}{2\sqrt{x}}\right)\left(1x(R_{k1}(x)1)C(x)\right)^{r+1}}.$$

We also provide a combinatorial explanation for certain facts in Main Theorem. For example, we provide a combinatorial proof for the fact $`(ii)`$ in the Main Theorem for $`r=k=0`$.

Acknowledgments. The author expresses his appreciation to the referees for their careful reading of the manuscript and helpful suggestions.

## 2. Proofs: directly from definitions

In this section we present a proof for the Main Theorem which is based on the definitions of the Dyck paths.

Proof of the Main Theorem$`(i)`$. We start by proving the first part of the Main Theorem by introducing a bijection $`\mathrm{\Psi }`$ between the set of Dyck paths of length $`2n`$ with $`r`$ peaks at height $`k`$ and the set of Dyck paths of length $`2n`$ with $`r`$ valleys at height $`k2`$.

###### Theorem 2.1.

$`\mathrm{Peak}_k^r(x)=\mathrm{Valley}_{k2}^r(x)`$ for all $`k2`$.

###### Proof.

Let $`P=P_1,P_2,\mathrm{},P_{2n}`$ be a Dyck path of length $`2n`$ with exactly $`r`$ peaks at height $`k2`$ where $`P_j`$ are the points of the path $`P`$. For any point $`P_j`$ we define another point $`\mathrm{\Psi }(P_j)=Q_j`$ as follows. If $`P_j`$ appears as a point of a valley at height $`k2`$ then we define $`Q_j=P_j+(0,2)`$. If $`P_j`$ appear as a point of a peak at height $`k`$ then we define $`Q_j=P_j(0,2)`$ (this is possible since $`k2`$). Otherwise, we define $`Q_j=P_j`$. Therefore, we obtain a new path $`Q=Q_1,Q_2,\mathrm{},Q_{2n}`$, and by definition of $`Q`$ it is easy to see that $`Q`$ is a Dyck path of length $`2n`$ with exactly $`r`$ valleys at height $`k2`$ (see Figure 2).

In fact, it is easily verified that the map which maps $`P`$ to $`Q`$ is a bijection. This establishes the theorem. ∎

Formula for $`\mathrm{Valley}_0^0(x)`$. Let $`P`$ be a Dyck path with no valleys at height $`0`$. It is easy to see that $`P`$ has no valleys at height $`0`$ if and only if there exists a Dyck path $`P^{}`$ of length $`2n2`$ such that

$$P=\mathrm{up}\text{-}\mathrm{step},P^{},\mathrm{down}\text{-}\mathrm{step}.$$

Let $`P^{\prime \prime }`$ be the path that results by shifting $`P^{}`$ by $`(1,1)`$. Then the map $`\mathrm{\Theta }`$ which sends $`PP^{\prime \prime }`$ is a bijection between the set of all Dyck paths starting at $`(0,0)`$ and ending at $`(2n,0)`$ with no valleys, and the set of all Dyck paths starting at point $`(0,0)`$ and ending at $`(2n2,0)`$. Hence

$$\mathrm{Valley}_0^0(x)=1+xC(x),$$

where we count $`1`$ for the empty path, $`x`$ for the up-step and the down-step, and $`C(x)`$ for all Dyck paths $`P^{\prime \prime }`$.

Proof of the Main Theorem$`(ii)`$. First of all, let us present two facts. The first fact concerns the generating function for the number of Dyck paths from the southwest corner of a rectangle to the northeast corner.

###### Fact 2.2.

(see \[K, Theorem A2 with Fact A3\]) Let $`k0`$. The generating function for the number of Dyck paths which lie between the lines $`y=k`$ and $`y=0`$, starting at $`(0,0)`$ and ending at $`(n,k)`$ is given by

$$F_k(x):=\frac{1}{\sqrt{x}U_{k+1}\left(\frac{1}{2\sqrt{x}}\right)}.$$

The second fact concerns the generating function for the number of Dyck paths starting at $`(a,k+1)`$ and ending at $`(a+n,k+1)`$ with no valleys at height $`k`$.

###### Fact 2.3.

The generating function for the number of Dyck paths starting at $`(a,k+1)`$ and ending at $`(a+n,k+1)`$ with no valleys at height $`k`$ is given by

$$\frac{C(x)}{1x(R_{k+1}(x)1)C(x)}.$$

###### Proof.

Let $`P`$ be a Dyck path starting at $`(k+1,0)`$ and ending at $`(k+1,n)`$ with no valleys at height $`k`$. It is easy to see that $`P`$ has a unique decomposition of the form

$$P=W_1,\mathrm{down}\text{-}\mathrm{step},V_1,\mathrm{up}\text{-}\mathrm{step},W_2,\mathrm{down}\text{-}\mathrm{step},V_2,\mathrm{},\mathrm{up}\text{-}\mathrm{step},W_m,$$

where the following conditions holds for all $`j`$:

* $`W_j`$ is a path consisting of up-steps and down-steps starting and ending at height $`k+1`$ and never passes below the height $`k+1`$;
* $`V_j`$ is a path consisting of up-step and down-steps starting and ending at height $`k`$ and never passes over the height $`k`$ (see Figure 3).

Using \[K, Theorem 2\] we get that the generating function for the number of paths of type $`V_j`$ (shift for a Dyck path) is given by $`R_{k+1}(x)1`$. Using the fact that $`W_j`$ is a shift for a Dyck paths starting and ending on the $`x`$-axis we obtain the generating function for the number of Dyck paths of type $`W_j`$ is given by $`C(x)`$. If we sum over all the possibilities of $`m`$ then we have

$$C(x)\underset{m0}{}\left(xC(x)(R_{k+1}(x)1)\right)^m=\frac{C(x)}{1x(R_{k+1}(x)1)C(x)}.$$

Now we are ready to prove the second part of the Main Theorem.

###### Theorem 2.4.

The generating function $`\mathrm{Valley}_k^r(x)`$ is given by

$$\delta _{r,0}R_{k+1}(x)+\frac{x^rC^{r+1}(x)}{U_{k+1}^2\left(\frac{1}{2\sqrt{x}}\right)\left(1x(R_{k+1}(x)1)C(x)\right)^{r+1}}.$$

###### Proof.

Let $`P`$ be a Dyck path starting $`(0,0)`$ and ending at $`(2n,0)`$ with exactly $`r`$ valleys at height $`k`$. It is easy to see that $`P`$ has a unique decomposition of the form

$$\begin{array}{c}P=E_1,\mathrm{up}\text{-}\mathrm{step},D_0,\mathrm{down}\text{-}\mathrm{step},\mathrm{up}\text{-}\mathrm{step},D_1,\mathrm{down}\text{-}\mathrm{step},\mathrm{up}\text{-}\mathrm{step},\mathrm{},D_r,\mathrm{down}\text{-}\mathrm{step},E_2,\hfill \end{array}$$

where the following conditions holds:

* $`E_1`$ is a Dyck path that lies between the lines $`y=k`$ and $`y=0`$, starting at $`(0,0)`$, and ending at point on height $`k`$;
* $`D_j`$ is a Dyck path starting and ending at points on height $`k+1`$ without valleys at height $`k`$, for all $`j`$;
* $`E_2`$ is a Dyck path that lies between the lines $`y=k`$ and $`y=0`$, starting at point on height $`k`$, and ending at $`(2n,0)`$.

Using Fact 2.2 and Fact 2.3 we get the the desired result for all $`r1`$. Now, if we assume that $`r=0`$, then we must consider another possibility which is that all the Dyck paths lie between the lines $`y=k`$ and $`y=0`$, starting at $`(0,0)`$, and ending on the $`x`$-axis. Hence, using \[K, Theorem 2\] we get that the generating function for the number of these paths is given by $`R_{k+1}(x)`$. ∎

As a corollary of the Main Theorem$`(ii)`$ for $`k=0`$ (using \[M, Example 1.18\]) we get

###### Theorem 2.5.

For all $`r0`$,

$$\mathrm{Valley}_0^r(x)=\delta _{r,0}+x^{r+1}C^{r+1}(x).$$

In other words, the number of Dyck paths starting at $`(0,0)`$ and ending at $`(2n,0)`$ with exactly $`r`$ valleys at height $`0`$ is given by

$$\frac{r+1}{n}\left(\genfrac{}{}{0pt}{}{2nr1}{n+1}\right).$$

Proof of the Main Theorem$`(iii)`$. If we merge the first two parts of Main Theorem, then we get an explicit formula for $`\mathrm{Peak}_k^r(x)`$ for all $`r0`$ and $`k2`$. Besides, by definition there are no peaks at height $`0`$. Thus, it is left to find $`\mathrm{Peak}_1^r(x)`$ for all $`r0`$.

###### Theorem 2.6.

For all $`r0`$,

$$\mathrm{Peak}_1^r(x)=\delta _{r,0}+\frac{x^{3r+2}C^{2r+2}(x)}{(1x^2C^2(x))^{r+1}}.$$

###### Proof.

Let $`P`$ be a Dyck path starting at $`(0,0)`$ and ending at $`(2n,0)`$ with exactly $`r`$ peaks at height $`1`$. It is easy to see that $`P`$ has a unique decomposition of the form

$$P=D_0,\mathrm{up}\text{-}\mathrm{step},\mathrm{down}\text{-}\mathrm{step},D_1,\mathrm{up}\text{-}\mathrm{step},\mathrm{down}\text{-}\mathrm{step},\mathrm{},\mathrm{up}\text{-}\mathrm{step},\mathrm{down}\text{-}\mathrm{step},D_r,$$

where $`D_j`$ is a nonempty Dyck path starting and ending at point on the $`x`$-axis with no peaks at height $`1`$. Hence, the rest is easy to obtain by using \[D\]. ∎

For example, for $`r=0`$ the above theorem yields the main result of \[D\].

## 3. Proofs: Directly from Theorem 1.2

In this section we present another proof for the Main Theorem which is based on Roblet and Viennot \[RV, Proposition 1\] (see Theorem 1.2).

Let $`\lambda _j=x`$ for all $`j`$, $`\mu _j=x`$ for all $`jk`$, and $`\mu _k=z`$. Theorem 1.2 yields

$`(1)`$
$$\underset{r0}{}\mathrm{Peak}_k^r(x)z^r=\frac{1}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1{\displaystyle \frac{\mathrm{}}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1(zx){\displaystyle \frac{x}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1\mathrm{}}}}}}}}}}}}}}}}}},$$

where $`z`$ appears in the $`k`$th level. On the other hand, $`xC^2(x)=C(x)1`$, we have that

$`(2)`$
$$C(x)=\frac{1}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1\mathrm{}}}}}}.$$

Using the identities (1) and (2) with $`xC^2(x)=C(x)1`$ we get

###### Theorem 3.1.

The generating function $`_{r0}\mathrm{Peak}_k^r(x)z^r`$ is given by

$$\frac{1}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1{\displaystyle \frac{\mathrm{}}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1zx^2C^2(x)}}}}}}}}}}},$$

where the continued fraction contains exactly $`k`$ levels.

For example, Theorem 3.1 yields for $`z=0`$ the generating function $`\mathrm{Peak}_k^0(x)`$ as in the statement of Theorem 1.1. More generally, Theorem 3.1 yields an explicit expression for $`\mathrm{Peak}_k^r(x)`$ for any $`r1`$ by using the following lemma.

###### Lemma 3.2.

For all $`k1`$,

$$\frac{1}{1{\displaystyle \frac{x}{1{\displaystyle \frac{x}{1{\displaystyle \frac{\mathrm{}}{1{\displaystyle \frac{x}{1zxA}}}}}}}}}=R_k(x)\frac{1zR_{k1}(x)xAR_{k1}(x)}{1zR_k(x)xAR_k(x)}.$$

where the continued fraction contains exactly $`k`$ levels.

###### Proof.

Immediately, by using the identity $`R_{m+1}(x)=1/(1xR_m(x))`$ and induction on $`k`$. ∎

Therefore, using Theorem 3.1, the above lemma, and the identity $`R_{m+1}(x)=1/(1xR_m(x))`$, together with definitions of $`R_k(x)`$, we get the explicit expression for the generating function $`\mathrm{Peak}_k^r(x)`$ for any $`r1`$ (see the Main Theorem).