format: text

text:
# Even-visiting random walks: exact and asymptotic results in one dimension

## 1 Introduction

Random walk is one of the most ubiquitous concepts of statistical physics. The large-time behavior of random walks is rather insensitive to details such as the presence of an underlying lattice or of an elementary time step, so that a vast range of them belongs to the universality class of Brownian motion, their continuum limit being described by the diffusion equation.

Recently, in the context of surface growth, Noh, Park, and den Nijs  have introduced a model of random walks with the constraint that every visited site should be visited an even number of times. Such walks have been named even-visiting walks . It has been argued that the non-local constraint of being even-visiting changes the universality class of random walks, with the typical extension of a walk of $`n`$ steps scaling as $`n^{1/(d+2)}`$ in dimension $`d`$, instead of the usual $`n^{1/2}`$.

Derrida’s argument  for this behavior, to be recalled just below, demonstrates a deep analogy between the even-visiting walk problem, the trapping problem  of a diffusive particle which is absorbed by traps located at random positions, and the Lifshitz tails  of the density of states of electron or phonon spectra of disordered solids near band edges.

Consider for definiteness a random walk of $`n`$ steps on the hypercubic lattice in dimension $`d`$, with unit lattice spacing. The probability that this walk is confined within a sphere with radius $`Rn^{1/2}`$ scales as

$$p(R)\mathrm{exp}\left(\frac{j^2n}{2dR^2}\right),$$

where $`j`$ is the first positive zero of the Bessel function $`J_{(d2)/2}`$, so that $`j^2/R^2`$ is the lowest eigenvalue of the Laplace operator in the sphere, with Dirichlet boundary conditions. Then, if $`nR^d`$, most sites in the sphere must have been visited many times. Hence, to a good approximation, the probability that any given site has been visited an even number of times is $`1/2`$, independently of other sites. Thus, among all walks that remain within the sphere, the fraction of even-visiting ones is about

$$f(R)\mathrm{exp}\left(\mathrm{\Omega }\mathrm{ln}2R^d\right),$$

where $`\mathrm{\Omega }=\pi ^{d/2}/\mathrm{\Gamma }(d/2+1)`$ is the volume of the unit ball. Consequently, for large $`n`$, the main contribution to the number of even-visiting walks comes from ‘Lifshitz spheres’, i.e., optimal spheres whose radius $`R_\mathrm{L}`$ maximizes the product $`p(R)f(R)`$, hence

$$R_\mathrm{L}\left(\frac{j^2n}{d^2\mathrm{\Omega }\mathrm{ln}2}\right)^{1/(d+2)}.$$
(1.1)
Both inequalities needed to justify the argument are satisfied by this solution, in any dimension $`d`$. Eq. (1.1) gives an estimate for the typical maximum extent of an $`n`$-step even-visiting walk from its origin. The total number of such walks scales as $`(2d)^np(R_\mathrm{L})f(R_\mathrm{L})`$, hence

$$𝒩(n)(2d)^n\mathrm{exp}\left(\frac{d+2}{2}(\mathrm{\Omega }\mathrm{ln}2)^{2/(d+2)}\left(\frac{j^2n}{d^2}\right)^{d/(d+2)}\right).$$
(1.2)
In this paper we reconsider the one-dimensional case  in detail. We map the even-visiting walk problem onto a model with quenched disorder, analogous to a non-Hermitian Anderson model. Many concepts and techniques of one-dimensional disordered systems are then available . We develop efficient numerical tools to enumerate even-visiting walks. The parameter in the generating function of the walks governs the strength of disorder. The disordered system comes with an invariant measure, whose support is bounded for small disorder, and becomes unbounded at some critical point, at which the generating function of even-visiting walks develops a singularity. We first compute the escape probability (4.9), which is similar to the integrated density of states of the Anderson model. This allows to take advantage of the existence of a simpler problem, the random-mirror problem, for which many quantities can be computed analytically, and then translated in terms of even-visiting walks. Numerical checks show that the analogy between both problems indeed works at a quantitative level. We thus obtain in particular the estimate (6.11) for the number of closed even-visiting walks, and eqs. (7.3), (7.4) for the cumulants of the maximum height reached by such a walk.

## 2 Definitions and mapping to a disordered system

An $`n`$-step walk on a graph is a sequence $`v_0v_1\mathrm{}v_n`$ of vertices of the graph such that $`v_0v_1,\mathrm{},v_{n1}v_n`$ are edges of the graph. An even-visiting walk is such that for each vertex $`vv_0,v_n`$ of the graph, the number of indices $`i`$ such that $`v=v_i`$ is even. Note that this is a non-local constraint on the walk, so that it is not unlikely that statistical properties of even-visiting random walks are rather different from those of standard random walks.

In this paper we shall concentrate on one-dimensional graphs, namely the segments with vertices $`0,1,\mathrm{},N`$ or the half-line with vertices $`0,1,2,\mathrm{}`$ In each case, the edges connect nearby integers.

We shall count even-visiting closed walks starting and ending at the origin (vertex $`0`$). This problem can be reformulated as a one-dimensional disordered system . Consider walks on a half-line with vertices numbered $`0,1,2,\mathrm{}`$ Each walk comes with a weight, which is computed as follows: each traversed edge contributes a factor $`\sqrt{t}`$, and each arrival at site $`i`$ contributes a factor $`w_i`$ for $`i=0,1,2,\mathrm{}`$ We consider for the time being that $`t`$ and $`w_0,w_1,\mathrm{}`$ are commuting indeterminates.

Let $`F(t,w_0,w_1,\mathrm{})`$ be the generating function for walks $`P`$ starting and ending at the origin, each walk being counted with its weight. Formally

$$F(t,w_0,w_1,\mathrm{})=\underset{P}{}\sqrt{t}^{S(P)}w_0^{V_0(P)}w_1^{V_1(P)}\mathrm{},$$

where $`S(P)`$ is the number of steps, and $`V_i(P)`$ the number of arrivals to site $`i`$. The only walk with $`0`$ steps has weight $`1`$. Any other walk can be decomposed uniquely as a succession of elementary blocks of the following type: a step from site $`0`$ to site $`1`$, giving weight $`\sqrt{t}w_1`$, a closed walk on the half line made of the sites $`1,2,\mathrm{}`$, contributing to $`F(t,w_1,w_2,\mathrm{})`$, and a step from site $`1`$ to site $`0`$, giving weight $`\sqrt{t}w_0`$. So

$`F(t,w_0,w_1,\mathrm{})=1`$ $`+`$ $`\sqrt{t}w_1F(t,w_1,w_2,\mathrm{})\sqrt{t}w_0`$
$`+`$ $`\left(\sqrt{t}w_1F(t,w_1,w_2,\mathrm{})\sqrt{t}w_0\right)^2+\mathrm{}`$
A resummation of the geometric series gives

$$F(t,w_0,w_1,\mathrm{})=\frac{1}{1tw_0w_1F(t,w_1,w_2,\mathrm{})},$$
(2.1)
and leads to a continued-fraction expansion

$$F(t,w_0,w_1,\mathrm{})=\frac{1}{1{\displaystyle \frac{tw_0w_1}{1{\displaystyle \frac{tw_1w_2}{1\mathrm{}}}}}}.$$
(2.2)
These formulæ are the starting point of all further considerations. Note that, in the expansion of $`F(t,w_0,w_1,\mathrm{})`$, the even-visiting random walks correspond exactly to monomials such that $`V_1(P),V_2(P),\mathrm{}`$ are all even. Then $`V_0(P)`$ as well, because $`V_0(P)+V_1(P)+\mathrm{}`$ is the number of steps, which is even for any closed walk. So, if the weights $`w_i`$ are turned into independent random variables with vanishing odd moments and all even moments equal to $`1`$, the average of $`F(t,w_0,w_1,\mathrm{})`$, denoted by $`\overline{F(t,w_0,w_1,\mathrm{})}`$, is the generating function for even-visiting walks on the half-line. In the same way, for finite segments, $`w_0,\mathrm{},w_N`$ are chosen as above, and $`w_{N+1}=0`$.

Now, a random variable with vanishing odd moments and all even moments equal to $`1`$ is simply a random sign: it takes values $`\pm 1`$ with probability $`1/2`$. The above formulation can therefore be further simplified: if $`w_0,w_1,\mathrm{}w_N`$ are independent random signs, then so are $`\epsilon _N=w_0w_1`$, $`\epsilon _{N1}=w_1w_2`$, $`\mathrm{}`$, $`\epsilon _1=w_{N1}w_N`$.

The generating series of even-visiting walks on the segment $`0,\mathrm{},N`$ reads therefore

$$F_N(t,\epsilon _1,\mathrm{},\epsilon _N)=\frac{1}{1{\displaystyle \frac{t\epsilon _N}{1{\displaystyle \frac{t\epsilon _{N1}}{1{\displaystyle \frac{\mathrm{}}{1t\epsilon _1}}}}}}}.$$
(2.3)
It is clear from this expression that averages over the quenched disorder represented by the $`\epsilon _i`$’s will be even functions of $`t`$. As each power of $`t`$ counts for two steps, the length (number of steps) of closed even-visiting walks is a multiple of four.

Eq. (2.3) implies the recursion formula

$$F_N(t,\epsilon _1,\mathrm{},\epsilon _N)=\frac{1}{1t\epsilon _NF_{N1}(t,\epsilon _1,\mathrm{},\epsilon _{N1})}.$$
(2.4)
If we write $`F_N=\mathrm{\Psi }_N/\mathrm{\Psi }_{N+1}`$, with $`\mathrm{\Psi }_0=\mathrm{\Psi }_1=1`$, then we have

$$\mathrm{\Psi }_{N+1}\mathrm{\Psi }_N+t\epsilon _N\mathrm{\Psi }_{N1}=0.$$
(2.5)
This three-term linear recursion relation defines a non-Hermitian one-dimensional Anderson model with off-diagonal binary disorder. In this context, $`t`$ can be interpreted as the strength of the disorder. The variables $`F_N`$, which obey the recursion relation (2.4), are the associated Riccati variables .

It is clear that the $`t`$-expansions of $`F_N(t,\epsilon _1,\mathrm{},\epsilon _N)`$ and $`F_{N+1}(t,\epsilon _0,\epsilon _1,\mathrm{},\epsilon _N)`$ coincide up to order $`N`$ (a closed walk of $`2N`$ steps cannot visit vertices higher than $`N`$), so all formulæ have a well-defined $`N\mathrm{}`$ limit.

The $`t`$-expansion of $`F_N`$ starts with $`1`$, so its power $`F_N^\alpha `$ is well-defined as a formal power series in $`t`$, for $`\alpha `$ an arbitrary complex number, or even an indeterminate. Surely, only when $`\alpha =1`$ does this quantity have a simple interpretation as a counting problem. Nevertheless, considering arbitrary $`\alpha `$ also has its own interest (see section 8). Anyway, we set

$$f_N(t,\alpha )=\overline{F_N^\alpha (t)},f(t,\alpha )=\underset{N\mathrm{}}{lim}f_N(t,\alpha ).$$

## 3 Weak-disorder expansion

If we formally expand $`F_N^\alpha (t)`$ from eq. (2.4) in powers of $`\epsilon _N`$, we find

$$F_N^\alpha (t)=1+\underset{m1}{}\epsilon _N^mF_{N1}^mt^m\frac{\mathrm{\Gamma }(\alpha +m)}{\mathrm{\Gamma }(\alpha )m!},$$

so that after averaging (and with $`m=2n`$)

$$f_N(t,\alpha )=1+\underset{n1}{}f_{N1}(t,2n)t^{2n}\frac{\mathrm{\Gamma }(\alpha +2n)}{\mathrm{\Gamma }(\alpha )(2n)!},$$

and for $`N\mathrm{}`$

$$f(t,\alpha )=1+\underset{n1}{}f(t,2n)t^{2n}\frac{\mathrm{\Gamma }(\alpha +2n)}{\mathrm{\Gamma }(\alpha )(2n)!}.$$
(3.1)
Let $`𝒩_k`$ (resp. $`𝒩_{N,k}`$) denote the number of even-visiting closed walks of $`4k`$ steps starting at the origin on the half-line (resp. on the segment $`[0,N]`$), and let us introduce the $`t`$-expansions

$$f_N(t,\alpha )=\underset{n0}{}f_{N,2n}(\alpha )t^{2n},f(t,\alpha )=\underset{n0}{}f_{2n}(\alpha )t^{2n}.$$
(3.2)
Eq. (3.1) yields the following recursion relation:

$$f_{N,2m}(\alpha )=\underset{n=1}{\overset{m}{}}f_{N1,2m2n}(2n)\frac{\mathrm{\Gamma }(\alpha +2n)}{\mathrm{\Gamma }(\alpha )(2n)!},$$

with initial conditions $`f_{N,0}(\alpha )=1`$ for $`N0`$. The above formula is recursive both in $`N`$ and $`m`$, and it holds for $`N1`$, $`m1`$. All the above quantities are non-decreasing functions of $`N`$. Letting $`N\mathrm{}`$, we obtain

$$f_{2m}(\alpha )=\underset{n=1}{\overset{m}{}}f_{2m2n}(2n)\frac{\mathrm{\Gamma }(\alpha +2n)}{\mathrm{\Gamma }(\alpha )(2n)!}$$
(3.3)
for $`m1`$, with the initial condition $`f_0(\alpha )=1`$.

This formalism can be used to compute the number of even-visiting walks in a very efficient way. To do so, it is better to trade the variable $`\alpha `$ for integers, to compute recursively an array $`f_{2m}(2\mathrm{})`$ for $`0\mathrm{}mk1`$, using

$$f_{2m}(2\mathrm{})=\underset{n=1}{\overset{m}{}}f_{2m2n}(2n)\left(\genfrac{}{}{0pt}{}{2n+2\mathrm{}1}{2n}\right),$$

and finally to obtain the number of closed even-visiting walks with $`n=4k`$ steps on the half-line, starting at $`0`$, as

$$𝒩_k=f_{2k}(1)=\underset{\mathrm{}=1}{\overset{k}{}}f_{2k2\mathrm{}}(2\mathrm{}).$$
(3.4)
The same idea can be implemented in the case of walks on a finite segment.

With a few hours of computation on a workstation, using the MACSYMA software, $`𝒩_k`$ has been evaluated for $`k`$ up to $`457`$, and $`𝒩_{N,k}`$ for $`Nk125`$. Note that $`𝒩_{N,k}=𝒩_k`$ for $`Nk`$. It is therefore possible to perform accurate investigations of the number of walks and of related quantities, such as the distribution of the maximum height (see sections 6 to 8).

It is obvious that $`𝒩_k16^k`$, the total number of walks of $`4k`$ steps. It is easy to generalize this (poor) bound to see that $`f_{2k}(\alpha )2^\alpha 16^k`$. So, the radius of convergence of $`f(t,\alpha )`$ in $`t^2`$ is at least $`1/16`$. In fact, for real positive $`\alpha `$, it is exactly $`1/16`$, as can be seen from the above equations as follows. Suppose $`\mathrm{}1`$. From $`f(t,2\mathrm{})=1+_{m1}f(t,2m)t^{2m}\left(\genfrac{}{}{0pt}{}{2\mathrm{}+2m1}{2m}\right)`$, we see that the $`t^2`$-expansion of $`\left(1t^2\mathrm{}\left(\genfrac{}{}{0pt}{}{4\mathrm{}1}{2\mathrm{}}\right)\right)f(t,2\mathrm{})1=_{m0,\mathrm{}}t^{2m}\left(\genfrac{}{}{0pt}{}{2\mathrm{}+2m1}{2m}\right)f(t,2m)`$ has only non-negative coefficients. Hence the same is true of the $`t^2`$-expansion of $`f(t,2\mathrm{})\left(1t^2\mathrm{}\left(\genfrac{}{}{0pt}{}{4\mathrm{}1}{2\mathrm{}}\right)\right)^1`$. Inserting this bound for $`\mathrm{}=m`$ on the right-hand side of eq. (3.1) shows that, for $`\alpha >0`$, the $`t^2`$-expansion of $`f(t,\alpha )1_{m0}t^{2m}\left(1t^{2m}\left(\genfrac{}{}{0pt}{}{4m1}{2m}\right)\right)^1\mathrm{\Gamma }(\alpha +2m)/(\mathrm{\Gamma }(\alpha )(2m)!)`$ has only non-negative coefficients. Note that the radius of convergence (in powers of $`t^2`$) of $`\left(1t^{2m}\left(\genfrac{}{}{0pt}{}{4m1}{2m}\right)\right)^1`$ is very close to $`1/16`$ for large $`m`$. This implies that $`f(t,\alpha )`$ is singular at $`t^2=1/16`$, giving the announced value for the radius of convergence. We have not been able to refine the above bounds in a useful way, starting from the weak-disorder expansion. In the next sections we obtain much better estimates by a different method.

## 4 Invariant measure and escape probability

In order to get a feeling for the kind of singularities that appear at $`t^2=1/16`$, we first evaluate a quantity apparently unrelated to the counting of even-visiting random walks, namely the probability for the random variable $`F(t;\epsilon _1,\epsilon _2,\mathrm{})`$ to be larger than $`2`$:

$$P_{\mathrm{esc}}(t)=\mathrm{Prob}(F>2).$$

This quantity, referred to as the escape probability, vanishes identically for $`t^21/16`$, while it is non-zero for $`t^2>1/16`$. It is quite analogous to the integrated density of states in the Anderson model and similar disordered spectra .

For the time being, consider a fixed real positive $`t`$, and define the (complementary) distribution function

$$R_N(x)=\mathrm{Prob}(F_N>x).$$

The recursion relation (2.4) implies the following relation between the distribution functions of $`F_N`$ and $`F_{N1}`$:

$$2R_N(x)=R_{N1}\left(\frac{x1}{tx}\right)R_{N1}\left(\frac{1x}{tx}\right)+R_{N1}\left(\frac{1}{t}\right)R_{N1}\left(\frac{1}{t}\right)+2\mathrm{\Theta }(x),$$

where $`\mathrm{\Theta }`$ is Heaviside’s the step function:

$$\mathrm{\Theta }(x)=\{\begin{array}{cc}0\hfill & \mathrm{for}x<0,\hfill \\ 1\hfill & \mathrm{for}x0.\hfill \end{array}$$

For large $`N`$, $`R_N(x)`$ approaches a limiting distribution function $`R(x)`$, which defines an invariant measure $`\mathrm{d}R(x)`$, and obeys the Dyson-Schmidt equation

$$2R(x)=R\left(\frac{x1}{tx}\right)R\left(\frac{1x}{tx}\right)+R\left(\frac{1}{t}\right)R\left(\frac{1}{t}\right)+2\mathrm{\Theta }(x).$$
(4.1)
The fractional linear mappings involved in this equation are the reciprocals of those involved in eq. (2.4).

### 4.1 The case $`0<t<1/4`$

In this situation, both mappings involved in eq. (2.4) \[or in eq. (4.1)\] are hyperbolic. Let $`M(t)=(1\sqrt{14t})/(2t)`$ be the smallest fixed point of the mapping $`x1/(1tx)`$, and let $`m(t)=1/(1+tM(t))`$. So, $`M(t)=1/(1tM(t))`$ and $`1/t<0<m(t)<1<M(t)<2<1/t`$. The mappings $`x1/(1tx)`$ and $`x1/(1+tx)`$ are respectively increasing and decreasing (on their intervals of continuity). Hence, if $`x`$ lies in the interval $`I(t)=]m(t),M(t)[`$, then the same is true of $`1/(1tx)`$ and $`1/(1+tx)`$. From $`F_0(t)=1`$, an inductive argument shows that $`F_N(t,\epsilon _1,\mathrm{},\epsilon _N)`$ takes its values in $`I(t)`$ for all $`N`$. Hence the support of the invariant measure is contained in $`I(t)`$. In particular, we have $`R(x)=0`$ for $`xM(t)`$ and $`R(x)=1`$ for $`xm(t)`$, so that eq. (4.1) can be further simplified for $`xI(t)`$ to

$$2R(x)=R\left(\frac{x1}{tx}\right)R\left(\frac{1x}{tx}\right)+1.$$
(4.2)
The graph of $`R(x)`$ is a (decreasing) devil’s staircase, with plateaus at the dyadic numbers $`R=m/2^n`$, with $`n=1,2,\mathrm{}`$ and $`1m`$ (odd) $`2^n1`$. The support of the invariant measure $`\mathrm{d}R(x)`$, i.e., the closure of the set of points $`x`$ around which $`R(x)`$ is not a constant, is a Cantor set of measure zero. The heights of the plateaus do not depend on $`t`$, but their size monotonically increases with $`t`$. When $`t0^+`$, $`R(x)`$ tends to the step function $`\mathrm{\Theta }(1x)`$. The other, more interesting limiting case $`t=1/4`$ is investigated in next subsection.

### 4.2 The case $`t=t_c=1/4`$

In this borderline case, the mapping $`x1/(1t_cx)`$ is parabolic: it has a degenerate fixed point $`x_c=2`$, which coincides with the upper bound of the interval $`I(1/4)=[2/3,2]`$.

The behavior of the distribution function $`R(x)`$ near this upper bound will be our first example of an exponentially small, Lifshitz-like singularity. For $`x`$ close to $`2`$, eq. (4.2) simplifies to

$$2R(x)=R\left(\frac{4(x1)}{x}\right).$$
(4.3)
The above equation is actually an identity for $`x`$ greater than $`6/7`$, the image of the lower bound $`2/3`$ by the mapping $`x1/(1+t_cx)`$. We set

$$y=\frac{2}{2x},x=2\frac{2}{y},$$
(4.4)
so that the mapping $`x4(x1)/x`$ corresponds to $`yy1`$. The general solution of eq. (4.3) therefore reads

$$R(x)=2^yA(y),$$
(4.5)
where $`A(y)`$ is a bounded periodic function of its argument, with unit period. Periodic amplitudes are quite frequent in the realm of one-dimensional disordered systems . The present situation is analogous to that of the invariant measure of the Anderson model right at a band edge . Figure 1 shows a plot of the periodic amplitude $`A(y)`$, obtained from exactly iterating the Dyson-Schmidt equation (4.2) a large enough number of times, starting from the initial data $`R_0(x)=\mathrm{\Theta }(1x)`$.

### 4.3 The case $`t>1/4`$

In this situation, the mapping $`x1/(1tx)`$ is elliptic: its two fixed points acquire an imaginary part. With the parametrization

$$t=\frac{1}{4\mathrm{cos}^2\varphi }(0<\varphi <\pi /2),$$
(4.6)
the fixed points read $`x_\pm =1+\mathrm{exp}(\pm 2\mathrm{i}\varphi )`$. As a consequence, the support of the invariant measure is unbounded. A fraction of the measure escapes above $`2`$, hence the name ‘escape probability’ for the quantity $`P_{\mathrm{esc}}=R(2)`$.

When $`t1/4\varphi ^2/4`$ is small, the fixed points $`x_\pm 2\pm 2\mathrm{i}\varphi `$ are close to $`x_c=2`$, so that the escape phenomenon takes place through a narrow channel. It is therefore legitimate to deform the result (4.5) as follows. Instead of eq. (4.4), we set

$$\mathrm{e}^{2\mathrm{i}\psi }=\frac{xx_+}{xx_{}},x=2\mathrm{cos}\varphi \frac{\mathrm{sin}(\psi \varphi )}{\mathrm{sin}\psi }(0<\psi <\pi ),$$
(4.7)
so that the elliptic mapping $`x(x1)/(tx)`$ corresponds to $`\psi \psi \varphi `$. The angle ratio $`\psi /\varphi `$ is therefore the appropriate deformation of the co-ordinate $`y`$ (up to an additive constant). There is actually no such constant, as eqs. (4.4) and (4.7) yield $`\psi =y\varphi +𝒪(\varphi ^3)`$. Eq. (4.5) therefore implies

$$R(x)2^{\psi /\varphi }A(\psi /\varphi ),$$
(4.8)
in the regime where $`\varphi `$ is small, while the angular variable $`\psi `$ is arbitrary. In particular, $`x=x_c=2`$ corresponds to $`\psi =\pi /2+\varphi `$, hence the following prediction for the escape probability:

$$P_{\mathrm{esc}}2^{\pi /(2\varphi )1}A\left(\pi /(2\varphi )\right)\frac{1}{2}\mathrm{exp}\left(\frac{\pi \mathrm{ln}2}{4\sqrt{t1/4}}\right)A\left(\frac{\pi }{4\sqrt{t1/4}}\right),$$
(4.9)
which is the result we were after. The escape probability is therefore exponentially small in $`(t1/4)^{1/2}`$, with an asymptotically periodic amplitude $`A`$, already entering eq. (4.5), and shown in Figure 1. The result (4.9) is fully analogous to the Lifshitz tail of the integrated density of states in the Anderson model and similar disordered spectra .

We have not been able to obtain by direct means estimates of the (more interesting) singularities that govern the asymptotics of the number of even-visiting random walks. However, the result (4.9) for the escape probability is quite suggestive, especially because of its formal analogy with the Lifshitz tail for the integrated density of states of electron or phonon spectra. The crucial ingredient in this result is that the transformation $`x1/(1tx)`$ turns from hyperbolic to elliptic at $`t=t_c=1/4`$. Equivalently, the singularity (4.9) is due to the occurrence in the recursive formula (2.4) of regions made of a large number, of order $`\pi /(2\varphi )`$, of consecutive positive $`\epsilon `$’s, which are one-dimensional analogues of Lifshitz spheres. The other transformation $`x1/(1+tx)`$ has no impact on the exponentially small form of the singularity. Its influence is only felt on the periodic modulation.

In the next section, we introduce a simpler model, the random-mirror model, which keeps the essential feature of the original even-visiting walk problem, i.e., the transformation $`x1/(1tx)`$, but for which quantities of interest are directly computable by elementary means. The random-mirror model is similar to the binary random harmonic chain where a finite fraction of the atoms has an infinite mass, first considered by Domb et al. . In this limiting case, the system splits into an infinite collection of independent finite molecules, so that the integrated density of states, among many other quantities, can be evaluated by simple algebra. The same simplification occurs for diffusion in the presence of perfectly absorbing sites . By analogy with these situations, we claim that the random-mirror model and the even-visiting walk problem have similar exponentially small singularities. The simple and explicitly computable periodic amplitudes of the first model are just replaced by (admittedly complicated) unknown periodic functions. In the leading terms of large-order $`t`$-expansions, this replacement is just responsible for an absolute prefactor.

## 5 A case study: the random-mirror model

If the distribution of the variables $`\epsilon _i`$ is modified, so that they take value $`0`$ or $`1`$ with respective probabilities $`p`$ or $`1p`$ (before they took value $`1`$ or $`1`$ with probability $`1/2`$), the equations obtained previously for even-visiting walks are modified in a straightforward way. We assume $`p0,1`$. The physical interpretation of the present model is that any site $`i=1,2,\mathrm{}`$ is totally reflecting with probability $`p`$, so that the particle freely moves between the origin and the first reflecting site, hence the name ‘random-mirror model’.

From a more mathematical viewpoint, this amounts to replacing the map $`x1/(1+tx)`$ by the constant map $`x1`$, without changing the second map $`x1/(1tx)`$. We shall also consider the case of a general constant map $`xX`$ (instead of $`x1`$), where $`X`$ is an indeterminate (most of the time specialized to a real number $`X<2`$).

We shall compare in detail the random-mirror model and the even-visiting walk model, so that it is useful to change notation: for the random-mirror model we use $`G`$, $`g`$, and $`S`$, instead of $`F`$, $`f`$, and $`R`$.

### 5.1 Weak-disorder expansion

Let us first fix $`X=1`$. From the recursion

$$G_N(t,\epsilon _1,\mathrm{},\epsilon _N)=\frac{1}{1t\epsilon _NG_{N1}(t,\epsilon _1,\mathrm{},\epsilon _{N_1})},$$

we get

$$g_N(t,\alpha )=\overline{G_N^\alpha (t)}=1+p\underset{n1}{}g_{N1}(t,n)t^n\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!},$$

which in turn leads to a recursion relation for the coefficients of the expansion $`g_N(t,\alpha )=_ng_{N,n}(\alpha )t^n`$, namely

$$g_{N,m}(\alpha )=p\underset{n=1}{\overset{m}{}}g_{N1,mn}(n)\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!},$$

with initial conditions $`g_{N,0}(\alpha )=1`$ for $`N0`$. For $`N\mathrm{}`$, we get

$$g(t,\alpha )=1+p\underset{n1}{}g(t,n)t^n\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!},$$
(5.1)
and

$$g_m(\alpha )=p\underset{n=1}{\overset{m}{}}g_{mn}(n)\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!}$$
(5.2)
for $`m1`$, with the initial condition $`g_0(\alpha )=1`$. We can repeat the argument given in the previous section for $`f(t,\alpha )`$ to show that for positive $`\alpha `$, the $`t`$-expansion of $`g(t,\alpha )`$ has radius of convergence $`1/4`$. The above formalism is easily extended to generic values of $`X`$. For example, eq. (5.1) becomes

$$g(t,\alpha )=1p+pX^\alpha +p\underset{n1}{}g(t,n)t^n\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!},$$

while eq. (5.2) remains unchanged, with the initial condition being modified to $`g_0(\alpha )=1p+pX^\alpha `$.

### 5.2 Escape probability

We now evaluate the escape probability by more direct means than for the case of even-visiting walks. The random variable $`G`$ is equal to $`Q_n`$ with probability $`(1p)p^n`$, where the $`Q_n`$ are defined recursively as follows:

$$Q_{n+1}=\frac{1}{1tQ_n},Q_0=X.$$
(5.3)
This definition gives a meaning to $`Q_n`$ for any integer $`n`$, positive of negative, and shows that $`Q_n`$ is a linear fractional function of $`X`$. The distribution function $`S(x)=\mathrm{Prob}(G>x)`$ reads explicitly

$$S(x)=(1p)\underset{n0}{}p^n\mathrm{\Theta }(Q_nx).$$
(5.4)
For the random-mirror model stricto sensu, we have $`X=1`$. Suppose $`t>1/4`$. With the parametrization (4.6), we infer from the recursion (5.3) the explicit expression

$$Q_n=2\left(1\mathrm{sin}\varphi \frac{\mathrm{cos}(n+1)\varphi }{\mathrm{sin}(n+2)\varphi }\right).$$
(5.5)
Then $`Q_n2`$ is equivalent to either $`\mathrm{sin}(n+2)\varphi >0`$ and $`\mathrm{cos}(n+1)\varphi 0`$, or $`\mathrm{sin}(n+2)\varphi <0`$ and $`\mathrm{cos}(n+1)\varphi 0`$. The smallest $`n`$ such that $`Q_n2`$ belongs to the first case, and reads<sup>2</sup><sup>2</sup>2The symbol $`x`$ denotes the smallest integer larger than or equal to $`x`$, while $`x`$ denotes the largest integer smaller than or equal to $`x`$. For $`x`$ not an integer, we have $`x=x+1`$, and the fractional part of $`x`$ is defined as $`Frac(x)=xx=xx+1`$. $`n_0=\pi /(2\varphi )1`$. When $`t`$ is close to $`1/4`$, i.e., $`\varphi `$ is small, many $`n`$’s following $`n_0`$, up to $`n_12n_0`$, have $`Q_n2`$. As a consequence, $`P_{\mathrm{esc}}=S(2)`$ can be estimated as $`(1p)_{n_0nn_1}p^n`$, leading to the result

$$P_{\mathrm{esc}}p^{\pi /(2\varphi )1},$$
(5.6)
up to a correction of order $`𝒪(p^{\pi /\varphi })`$.

When the initial point is not $`1`$, but a generic $`X`$, eq. (5.5) generalizes to

$$Q_n=2\left(1\mathrm{sin}\varphi \frac{\mathrm{cos}(n+1)\varphi +(1X)\mathrm{cos}(n1)\varphi }{\mathrm{sin}(n+2)\varphi +(1X)\mathrm{sin}n\varphi }\right).$$
(5.7)
The smallest $`n`$ such that $`Q_n2`$ reads $`n_0(X)=\pi /(2\varphi )\mathrm{\Delta }(\varphi ,X)`$, with

$$\mathrm{\Delta }(\varphi ,X)=\frac{1}{\varphi }\mathrm{arctan}\left(\frac{X}{2X}\mathrm{tan}\varphi \right).$$
(5.8)
Hence

$$P_{\mathrm{esc}}(X)p^{\pi /(2\varphi )\mathrm{\Delta }(\varphi ,X)},$$

again up to a correction of order $`𝒪(p^{\pi /\varphi })`$. For small $`\varphi `$, this result can be recast in a form similar to eq. (4.9), namely

$$P_{\mathrm{esc}}(X)p^{\pi /(2\varphi )}A_X\left(\pi /(2\varphi )\right),$$
(5.9)
where the periodic amplitude $`A_X(x)`$ has the explicit expression

$$A_X(x)=p^{Frac(x\mathrm{\Delta }_0(X))+1\mathrm{\Delta }_0(X)},$$

with

$$\mathrm{\Delta }_0(X)=\mathrm{\Delta }(0,X)=\frac{X}{2X}.$$

This result shows that the random-mirror problem at $`p=1/2`$ and even-visiting walks have the same leading singularities (compare eq. (5.9) and eq. (4.9)): the escape probability of the random-mirror problem ($`X=1`$), of the generalized one ($`X`$ generic), and of even-visiting walks only differ by their periodic modulation. The same property is known in the case of spectra of disordered systems: details of the distribution of random masses, random site potentials, and so on, only enter the shape of the periodic amplitudes of Lifshitz tails .

## 6 Number of even-visiting walks

The purpose of this section is to establish the asymptotic estimate (6.11) for the total number of even-visiting walks starting and ending at the origin, using eqs. (3.2) and (3.4).

The escape probability $`P_{\mathrm{esc}}`$ was easy to estimate in the even-visiting walk problem, because it is a purely singular quantity, which vanishes as the parameter $`t`$ approaches $`t_c=1/4`$. To the contrary, the singularities of the functions $`f(t,\alpha )`$ or $`g(t,\alpha )`$ at $`t=1/4`$ are less easily grasped, because the latter quantities also have a regular part.

We shall first consider the random-mirror model with $`X=1`$. In this case, we shall see that there is a simple relation between the discontinuity of $`g(t,1)`$ for $`t>1/4`$ and the escape probability. A similar correspondence is expected to hold for the discontinuity of $`f(t,1)`$, up to a periodic modulation. These discontinuities provide a direct way to estimate the large-order behavior of $`t`$-expansions. The latter will lead us to estimate the number of even-visiting random walk (eq. (6.11)).

### 6.1 Discontinuity of $`g(t,1)`$

Let us consider first the random-mirror model in the simpler situation where $`X=1`$. By analogy with eq. (5.4), we have the closed formula

$$g(t,\alpha )=(1p)\underset{n0}{}p^nQ_n^\alpha .$$
(6.1)
Their recursive definition (5.3) shows that the $`Q_n(t)`$ are rational functions of $`t`$. Their explicit expression (5.5) yields the decomposition

$$Q_n(t)=1+\frac{4}{n+2}\underset{m=1}{\overset{(n+1)/2}{}}\left(\frac{1}{14t\mathrm{cos}^2\frac{\pi m}{n+2}}1\right)\mathrm{sin}^2\frac{\pi m}{n+2}.$$
(6.2)
It is therefore clear that all the poles of the rational function $`Q_n(t)`$ lie on the half-line $`]1/4,+\mathrm{}[`$. We normalize the discontinuity of a real-analytic function $`f(t)`$ on the real axis as $`D_f(t)=(1/\pi )Imf(t\pm \mathrm{i0})`$, in such a way that the discontinuity of $`1/t`$ is exactly $`\delta (t)`$. Let $`t`$ be real and such that $`t1/4`$, and consider the integrated discontinuity

$$D(n,t)=_{1/4}^tD_{Q_n}(u)du.$$

Eq. (6.2) yields

$$D(n,t)=\frac{1}{n+2}\underset{m=1}{\overset{(n+1)/2}{}}\mathrm{\Theta }\left(t\frac{1}{4\mathrm{cos}^2\frac{\pi m}{n+2}}\right)\mathrm{tan}^2\frac{\pi m}{n+2}.$$
(6.3)
Eq. (6.1) implies that the integrated discontinuity of $`g(t,1)`$ reads

$$D_1(t)=(1p)\underset{n0}{}p^nD(n,t).$$

Eq. (6.3) leads to the following explicit result, in terms of $`\varphi `$:

$$D_1(\varphi )=\frac{1p}{p^2}\underset{n2}{}\frac{p^n}{n}\underset{m=1}{\overset{(n1)/2}{}}\mathrm{\Theta }\left(\varphi \frac{\pi m}{n}\right)\mathrm{tan}^2\frac{\pi m}{n},$$

which we reorganize as

$$D_1(\varphi )=\frac{1p}{p^2}\underset{m1}{}\underset{n2m1}{}\frac{p^n}{n}\mathrm{\Theta }\left(\varphi \frac{\pi m}{n}\right)\mathrm{tan}^2\frac{\pi m}{n}.$$

Finally, because $`0<\varphi <\pi /2`$, the condition $`\pi m/n\varphi `$ is always more stringent than $`n2m1`$, hence

$$D_1(\varphi )=\frac{1p}{p^2}\underset{m1}{}\underset{n\pi m/\varphi }{}\frac{p^n}{n}\mathrm{tan}^2\frac{\pi m}{n}.$$

For small $`\varphi `$, the contribution of $`m=1`$ is exponentially larger than the other ones, so that the asymptotic behavior of $`D_1(\varphi )`$ is

$$D_1(\varphi )\frac{\varphi ^3}{\pi }p^{\pi /\varphi 2}.$$
(6.4)
The above formula entirely comes from the pole of $`Q_n(t)`$ closest to the critical point $`t_c=1/4`$. This gives the clue for treating the case of a general constant map. In that case, $`Q_n`$, as given by eq. (5.7), has poles for $`\mathrm{tan}(n+1)\varphi /\mathrm{tan}\varphi =X/(2X)`$. Let $`\varphi _1(X,n)`$ be the smallest solution of the latter equation. One checks that, for large $`n`$,

$$\varphi _1(X,n)=\frac{\pi }{n}\frac{2\pi }{(2X)n^2}+𝒪(n^3).$$

The contributions of the closest poles for each $`n`$ to the total integrated discontinuity are

$$(1p)\underset{n}{}\frac{\mathrm{sin}^2\varphi _1}{\mathrm{cos}^3\varphi _1}p^n\frac{\mathrm{cos}(n+1)\varphi _1+(1X)\mathrm{cos}(n1)\varphi _1}{(n+2)\mathrm{cos}(n+2)\varphi _1+(1X)n\mathrm{cos}n\varphi _1}\mathrm{\Theta }(\varphi \varphi _1).$$

This expression looks formidable, but it simplifies drastically at small $`\varphi `$. From the definition of $`\varphi _1(X,n)`$, the smallest $`n`$ giving a non-vanishing contribution to the above sum is such that

$$\frac{\mathrm{tan}(n+1)\varphi }{\mathrm{tan}\varphi }\frac{X}{2X}<\frac{\mathrm{tan}n\varphi }{\mathrm{tan}\varphi },$$

so that $`n=\pi /\varphi \mathrm{\Delta }(\varphi ,X)1`$, with the definition (5.8). Hence

$$D_X(\varphi )\frac{\varphi ^3}{\pi }p^{\pi /\varphi \mathrm{\Delta }(\varphi ,X)1},$$

i.e.,

$$D_X(\varphi )\frac{\varphi ^3}{\pi }p^{\pi /\varphi }B_X(\pi /\varphi ),$$
(6.5)
with

$$B_X(x)=p^{Frac(x\mathrm{\Delta }_0(X))\mathrm{\Delta }_0(X)}=\frac{A_X(x)}{p}.$$
(6.6)
The exponentially small factor $`p^{\pi /\varphi }`$ in the result (6.5) for the discontinuity of $`g(t,1)`$ is just the square of the corresponding factor in the result (5.9) for the escape probability. This correspondence, which holds for any value of $`X`$, is expected to hold as well for the even-visiting walk problem. The simple relationship between the periodic amplitudes $`A_X`$ and $`B_X`$ is, however, a peculiarity of the random-mirror problem, that we do not expect to hold in general. We quote for further reference the constant Fourier component of the periodic function $`B_X(x)`$,

$$B_X^{(0)}=_0^1B_X(x)dx=\frac{1p}{\left|\mathrm{ln}p\right|}p^{\mathrm{\Delta }_0(X)1}=\frac{1p}{\left|\mathrm{ln}p\right|}p^{2/(2X)}.$$
(6.7)
### 6.2 Asymptotic results

We again consider first the generalized random-mirror model, with arbitrary $`X`$. The knowledge of the integrated discontinuity $`D_X(t)`$ allows to compute the asymptotics of the $`t`$-expansion of $`g(t,1)`$. By definition, the coefficient of order $`k`$ in this expansion reads

$$g_k(1)=\frac{\mathrm{d}t}{2\pi \mathrm{i}t^{k+1}}g(t,1)=_{1/4}^+\mathrm{}\frac{\mathrm{d}t}{t^{k+1}}\frac{\mathrm{d}D_X}{\mathrm{d}t}=_0^{\pi /2}(4\mathrm{cos}^2\varphi )^{k+1}\frac{\mathrm{d}D_X}{\mathrm{d}\varphi }d\varphi .$$
(6.8)
For large $`k`$, the last integral is dominated by small values of $`\varphi `$, where $`\mathrm{cos}\varphi `$ can be expanded to quadratic order and exponentiated, while the asymptotic result (6.5) holds for $`D_X(\varphi )`$. It turns out that only the constant Fourier component $`B_X^{(0)}`$ of the periodic amplitude $`B_X(x)`$ matters (see below eq. (6.10)). We thus obtain

$$g_k(1)4^{k+1}\left|\mathrm{ln}p\right|B_X^{(0)}_0^{\mathrm{}}\varphi \mathrm{e}^{k\varphi ^2\pi \left|\mathrm{ln}p\right|/\varphi }d\varphi .$$
(6.9)
This integral is then evaluated by the saddle-point approximation, with the saddle-point value of $`\varphi `$ being $`\varphi _c=(\pi \left|\mathrm{ln}p\right|/(2k))^{1/3}`$. We are thus left with the explicit asymptotic result

$$g_k(1)4\left|\mathrm{ln}p\right|^{4/3}B_X^{(0)}\left(\frac{\pi }{2k}\right)^{5/6}\mathrm{exp}\left(\frac{3}{2}\left(2\pi ^2(\mathrm{ln}p)^2k\right)^{1/3}\right)4^k,$$
(6.10)
which only depends on $`X`$ through $`B_X^{(0)}`$ (see eq. (6.7)).

It is now clear why only the constant component $`B_X^{(0)}`$ matters, just as in trapping problems . Consider the Fourier expansion $`B_X(x)=_kB_X^{(k)}\mathrm{e}^{2\pi \mathrm{i}kx}`$. The contribution of the Fourier component $`B_X^{(k)}`$ is obtained by replacing in eq. (6.10) $`\left|\mathrm{ln}p\right|^{2/3}`$ by $`(\mathrm{ln}p2\pi \mathrm{i}k)^{2/3}`$, a number with a strictly greater real part, for any $`k0`$. The oscillations of the result (6.5) therefore get damped by exponentials of $`k^{1/3}`$ in the large-$`k`$ asymptotics of quantities such as $`g_k(1)`$.

Coming back to the even-visiting walk problem, using the correspondence between the singularities at $`t=1/4`$ of $`f(t,1)`$ and of $`g(t,1)`$ with $`p=1/2`$, we arrive to the following asymptotic expression for the number $`𝒩_k=f_{2k}(1)`$ of closed even-visiting walks of $`4k`$ steps:

$$𝒩_kB^{(0)}(2\mathrm{ln}2)^{4/3}\left(\frac{\pi }{k}\right)^{5/6}\mathrm{exp}\left(\frac{3}{2}\left((2\pi \mathrm{ln}2)^2k\right)^{1/3}\right)2^{4k}.$$
(6.11)
Note that the case of even-visiting walks has the symmetry $`tt`$, so that $`f(t,1)`$ has two cuts, at $`]\mathrm{},1/4]`$ and $`[1/4,+\mathrm{}[`$, which cancel each other at odd orders in $`t`$, and add up constructively at even orders.

In the expression (6.11), $`B^{(0)}`$ stands for the constant component of the periodic amplitude $`B(x)`$ of the even-visiting walk problem, which cannot be predicted analytically. The exact enumeration procedure described in section 3 has been carried out up to $`k=457`$, i.e., $`4k=1828`$ steps. The number of these walks (exactly evaluated as an integer) reads $`𝒩_{457}3.36829575\times 10^{535}`$. Figure 2 shows a plot of the ratio $`𝒩_k`$ to its asymptotic behavior $`𝒩_k^{(\mathrm{as})}`$, defined as the result (6.11) without its prefactor $`B^{(0)}`$. The data exhibit a smooth convergence in $`k^{2/3}`$. This form of the leading finite-$`k`$ correction is indeed expected for all asymptotic estimates, as it corresponds to regular corrections to leading critical behavior, of relative order $`tt_c`$, via the expression of the saddle-point $`\varphi _c`$ given below eq. (6.9). We thus obtain the accurate estimate

$$B^{(0)}1.760.$$
(6.12)
## 7 Distribution of the maximum height

The purpose of this section is to investigate the distribution of the maximum height reached by an even-visiting walk of $`4k`$ steps, starting and ending at the origin.

The generating function for even-visiting walks with maximum height $`M`$ is the difference $`f_M(t,1)f_{M1}(t,1)`$, and its coefficient of order $`2k`$ is $`f_{M,2k}(1)f_{M1,2k}(1)`$. Hence the probability that the maximum height in an even-visiting walk of $`4k`$ steps be $`M`$ is

$$\mathrm{\Pi }_{M,k}=\frac{f_{M,2k}(1)f_{M1,2k}(1)}{f_{2k}(1)}.$$
(7.1)
The generating function of the cumulants of $`M`$ reads therefore

$$Z_k(z)=\underset{j1}{}M^j\frac{z^j}{j!}=\mathrm{ln}\left(\underset{M0}{}\mathrm{\Pi }_{M,k}\mathrm{e}^{zM}\right)=\mathrm{ln}\left(\frac{1\mathrm{e}^z}{f_{2k}(1)}\underset{M0}{}f_{M,2k}(1)\mathrm{e}^{zM}\right).$$
(7.2)
Let us again make a detour through the random-mirror model. It is straightforward to extend eq. (7.2) to the latter situation, by simply replacing $`f`$’s by $`g`$’s, and $`Z`$ by $`Y`$, even if the probabilistic interpretation of the formula thus obtained is delicate.

The identity

$$g_N(t,1)=(1p)\underset{n=0}{\overset{N}{}}p^nQ_n(t)$$

implies

$$Y_k(z)=\mathrm{ln}\frac{c_k(p\mathrm{e}^z)}{c_k(p)},$$

with

$$c_k(p)=\frac{\mathrm{d}t}{2\pi \mathrm{i}t^{k+1}}\underset{n0}{}p^nQ_n(t)=\frac{g_k(1)}{1p}.$$

The large-$`k`$ estimate (6.10) for $`g_k(1)`$ leads to

$$Y_k(z)\frac{3}{2}(2\pi ^2k)^{1/3}\left((\mathrm{ln}p)^{2/3}(\mathrm{ln}pz)^{2/3}\right).$$

Coming back to the even-visiting walk problem, using the correspondence between both problems for $`p=1/2`$, we obtain

$$Z_k(z)\frac{3}{2}(4\pi ^2k)^{1/3}\left((\mathrm{ln}2)^{2/3}(\mathrm{ln}2z)^{2/3}\right).$$

Hence all the cumulants of the distribution of the maximum height $`M`$ of $`4k`$-step walks are asymptotically proportional, namely

$$M^ja_jk^{1/3},$$
(7.3)
with explicit prefactors

$$a_j=\frac{\mathrm{\Gamma }(j2/3)}{\mathrm{\Gamma }(1/3)}(2\pi )^{2/3}(\mathrm{ln}2)^{2/3j},$$
(7.4)
i.e., $`a_1=3.847495`$, $`a_2=1.850254`$, $`a_3=3.559136`$, $`a_4=11.981080`$, and so on.

As the number of steps becomes large, the distribution of the maximum height $`M`$ therefore gets more and more peaked around its mean value $`a_1k^{1/3}`$. The bulk of this distribution is a narrow Gaussian, with standard deviation $`a_2^{1/2}k^{1/6}`$. This result is in sharp contrast with usual random walk, for which the largest positive extent $`M`$ scales as $`t^{1/2}`$, and the ratio $`\xi =M/t^{1/2}`$ admits a non-trivial limiting probability law .

The above analytical results have been checked against data obtained from the exact enumeration procedure of section 3. Figure 3 shows a plot of the first four reduced cumulants $`M^j/a_j`$ $`(j=1,\mathrm{},4)`$. The clear linear behavior in $`k^{1/3}`$, with unit slope, demonstrates a quantitative agreement between the exact data and the asymptotic results (7.3), (7.4).

## 8 Amplitude ratios

Up to now, we have only considered averages of Riccati variables with respect to the quenched disorder, such as $`f(t,1)`$ and $`g(t,1)`$, which have a direct combinatorial interpretation. However, the even-visiting walk model and the random-mirror model also display an interesting dependence in $`\alpha `$, which will be investigated in this last section.

To start with, let us consider the random-mirror model with $`X=1`$, and assume that $`\alpha =\mathrm{}`$ is a positive integer. We can use eq. (6.2) to obtain the pole structure of $`Q_n^{\mathrm{}}(t)`$, which governs the asymptotics of its $`t`$-expansion. As before, it is clear that the pole closest to the critical point ($`m=1`$) gives the leading contribution. So, we write

$$Q_n^{\mathrm{}}(t)=\underset{j=0}{\overset{\mathrm{}}{}}\frac{c_j(n,\mathrm{})}{\left(14t\mathrm{cos}^2\frac{\pi }{n+2}\right)^j}+\mathrm{},$$
(8.1)
where the dots stand for a part which is regular at $`t=1/(4\mathrm{cos}^2\frac{\pi }{n+2})`$. It turns out that the coefficients $`c_j(n,l)`$ have a simple behavior when $`n`$ is large. Indeed, let us write $`\varphi =\pi /n\pi \lambda /n^2+𝒪(1/n^3)`$, so that $`t=1/(4\mathrm{cos}^2\frac{\pi }{n})\pi ^2\lambda /(2n^3)+𝒪(1/n^4)`$, and re-express $`Q_{n2}(t)`$ in terms of $`\lambda `$. Using eq. (5.5), we can check that $`Q_{n2}=2(1+1/\lambda )+𝒪(1/n)`$, so that

$$Q_{n2}^{\mathrm{}}=2^{\mathrm{}}\left(1+\frac{1}{\lambda }\right)^{\mathrm{}}+𝒪(1/n).$$

On the other hand, a direct substitution of the above expression for $`t`$ into eq. (8.1) yields, for large $`n`$,

$$Q_{n2}^{\mathrm{}}\underset{j=0}{\overset{\mathrm{}}{}}c_j(n2,\mathrm{})\left(\frac{2\pi ^2\lambda }{n^3}\right)^j.$$

Comparing the above two formulæ, we obtain, again for large $`n`$,

$$c_j(n,\mathrm{})2^{\mathrm{}}\left(\genfrac{}{}{0pt}{}{\mathrm{}}{j}\right)\left(\frac{2\pi ^2}{n^3}\right)^j.$$

Now, the $`t`$-expansion of eq. (8.1) is straightforward. A pole of order $`j`$ gives an extra factor $`(j+k1)!/((j1)!k!)`$ compared to a pole of order $`1`$. The sum over $`n`$ is then performed by the saddle-point approximation, with the saddle-point value of $`n`$ being $`n_c(2\pi ^2k/\left|\mathrm{ln}p\right|)^{1/3}`$. We are thus left with the following explicit values of the amplitude ratios:

$$𝒢(\mathrm{})\underset{k\mathrm{}}{lim}\frac{g_k(\mathrm{})}{g_k(1)}=2^\mathrm{}1\underset{j=1}{\overset{\mathrm{}}{}}\frac{\left|\mathrm{ln}p\right|^{j1}}{(j1)!}\left(\genfrac{}{}{0pt}{}{\mathrm{}}{j}\right).$$
(8.2)
It turns out that this result is independent of $`X`$, but we shall not go into any more detail.

Quite unexpectedly, the amplitude ratios $`𝒢(\mathrm{})`$ are also determined as the solution of an eigenvalue problem. The starting point is equation (5.2), which we divide by $`g_m(1)`$. Because the radius of convergence of the $`t`$-expansion of $`g(t,\alpha )`$ is $`1/4`$, we expect that, for large $`m`$ and fixed $`n`$, $`g_{mn}(n)4^ng_m(n)`$. So, we get

$$𝒢(\alpha )=p\underset{n1}{}𝒢(n)4^n\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!}.$$

In particular, fixing $`\alpha =\mathrm{}`$, a positive integer,

$$𝒢(\mathrm{})=p\underset{n1}{}𝒢(n)4^n\frac{(\mathrm{}+n1)!}{(\mathrm{}1)!n!}.$$
(8.3)
This equation implies that the (infinite) matrix $``$ with positive entries $`_{\mathrm{},n}=(\mathrm{}+n1)!/(4^n(\mathrm{}1)!n!)`$ for $`\mathrm{},n1`$ has an eigenvector with eigenvalue $`1/p`$ and positive components $`𝒢(\mathrm{})`$, for every $`p]0,1]`$. We have therefore an explicit example of why the Perron–Frobenius theorem (stating that a finite positive matrix has exactly one eigenvector with positive components) cannot extend to infinite matrices without further hypotheses. Note that the matrix $``$ is highly asymmetric, has unbounded entries, and a divergent trace.

It is interesting to compare the linear system (8.3) with eq. (5.1) at $`t=1/4`$, which reads

$$g(1/4,\alpha )=1+p\underset{n1}{}g(1/4,n)4^n\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }(\alpha )n!}.$$

In particular, fixing $`\alpha =\mathrm{}`$, a non-negative integer,

$$g(1/4,\mathrm{})=g(1/4,0)+p\underset{n1}{}g(1/4,n)4^n\frac{(\mathrm{}+n1)!}{(\mathrm{}1)!n!}.$$

So, despite the fact that analyticity for $`|t|<1/4`$ and continuity at $`t=1/4`$ fix $`g(1/4,\alpha )`$ unambiguously, the $`𝒢(\mathrm{})`$ are a (positive) ambiguity of the solution of eq. (5.1) restricted to $`t=1/4`$. This is again incompatible with the Perron–Frobenius property.

Because of these strange properties, it is reassuring to show that the explicit formula (8.2) for the $`𝒢(\mathrm{})`$ indeed satisfies (8.3), which was derived only heuristically. The proof relies on the following identity:

$$\frac{1}{(\mathrm{}1)!}\underset{j0}{}\frac{(j+\mathrm{})!}{j!(j+1)!}s^j=\mathrm{e}^s\mathrm{}!\underset{j=1}{\overset{\mathrm{}}{}}\frac{s^{j1}}{(j1)!j!(\mathrm{}j)!},$$

which holds for $`\mathrm{}`$ a positive integer, and can be checked by an explicit expansion of its right-hand side in powers of $`s`$. We thus obtain the explicit formula

$$𝒢(\alpha )=p\frac{2^{\alpha 1}}{\mathrm{\Gamma }(\alpha )}\underset{j0}{}\frac{\mathrm{\Gamma }(j+\alpha +1)}{j!(j+1)!}\left|\mathrm{ln}p\right|^j.$$
(8.4)
When $`\alpha `$ is large and positive, at fixed $`p1`$, the sum over $`j`$ can be evaluated by the saddle-point approximation, the saddle-point value of $`j`$ being $`j_c(\left|\mathrm{ln}p\right|\alpha )^{1/2}+\left|\mathrm{ln}p\right|/2`$, yielding

$$𝒢(\alpha )\frac{1}{4}\left(\frac{p^2\alpha }{\pi ^2\left|\mathrm{ln}p\right|^3}\right)^{1/4}\mathrm{exp}\left(2(\left|\mathrm{ln}p\right|\alpha )^{1/2}\right)\mathrm{\hspace{0.17em}2}^\alpha .$$
(8.5)
In the present case of the random-mirror model, we have been lucky enough to obtain the exact amplitude ratios (8.2). Had it not been the case, eq. (8.3) would provide a convenient tool to get accurate numerical values for these ratios. In fact, this alternative approach leads to a better understanding of how the Perron–Frobenius theorem is bypassed. If the matrix $`_{\mathrm{},n}`$ is truncated to $`\mathrm{}`$, $`n\mathrm{}_{\mathrm{max}}`$, then, whatever $`\mathrm{}_{\mathrm{max}}`$ is, the spectrum does not contains any eigenvalue larger than unity. So, the right trick is to truncate $`_{\mathrm{},n}`$ to $`\mathrm{}\mathrm{}_{\mathrm{max}}`$, $`n\mathrm{}_{\mathrm{max}}+1`$, to fix the normalization $`𝒢(1)=1`$, and to solve for $`𝒢(2),\mathrm{},𝒢(\mathrm{}_{\mathrm{max}}+1)`$.

Coming back to the even-visiting walk problem, we shall follow the latter strategy. We expect, by analogy with the previous situation, that $`(\alpha )=lim_k\mathrm{}\left(f_{2k}(\alpha )/f_{2k}(1)\right)`$ are well-behaved limits. Then, starting from eq. (3.3), we can repeat the above argument, to show that

$$(\alpha )=\underset{n1}{}(2n)16^n\frac{\mathrm{\Gamma }(\alpha +2n)}{\mathrm{\Gamma }(\alpha )(2n)!},$$

and, fixing again $`\alpha =\mathrm{}`$, a non-negative integer,

$$(\mathrm{})=\underset{n1}{}(2n)16^n\frac{(\mathrm{}+2n1)!}{(\mathrm{}1)!(2n)!}.$$
(8.6)
This result implies that eq. (3.1) that determines $`f(t,\alpha )`$ is ambiguous at $`t^2=1/16`$. By truncating the system (8.6) as explained above, we can obtain the amplitude ratios $`(\mathrm{})`$ with very high accuracy, namely $`(2)=5.591806`$, $`(3)=22.285850`$, $`(4)=77.059126`$, and so on.

These predictions have been checked against exact data for the ratios $`f_{2k}(\mathrm{})/f_{2k}(1)`$, obtained by means of the enumeration approach of section 3. Figure 4 demonstrates a quantitative agreement, for $`\mathrm{}=2`$ to 4. The data smoothly converge to the predicted limiting values $`(\mathrm{})`$, with (small) $`k^{2/3}`$ corrections. This reinforces our assertion that even-visiting walks and random-mirror models have very similar asymptotic behaviors, up to multiplicative prefactors.

## 9 Discussion

In this paper we have reconsidered the problem of even-visiting random walks, and obtained several kinds of exact or asymptotic results in one dimension.

We have mapped the even-visiting walk problem onto a non-Hermitian Anderson model (2.5). The weak-disorder expansion of the latter model provides very efficient numerical tools to enumerate and characterize even-visiting walks. We have thus been able to evaluate exactly, among other quantities of interest, the total number $`𝒩_k`$ of closed even-visiting walks up to $`k=457`$, i.e., $`4k=1828`$ steps, going thus far beyond previous works. Indeed, our result is to be compared with the exact result  up to $`80`$ steps, and with the approximate numerical simulation  up to around $`1000`$ steps.

The mapping to the non-Hermitian Anderson model (2.5) also makes many concepts and techniques of one-dimensional disordered systems available . The analogy with Lifshitz tails, which was already clearly apparent from Derrida’s argument  in any dimension, has been corroborated at the level of analytical tools in the one-dimensional situation. The escape probability, investigated in section 4, is fully analogous to the integrated density of states in disordered spectra , with an exponential behavior, modulated by an oscillatory amplitude, which is asymptotically periodic in the relevant variable. Periodic amplitudes are actually ubiquitous in one-dimensional disordered systems, when the quenched disorder has a discrete distribution .

In order to investigate more complex quantities, such as the generating functions of even-visiting walks, we make an extensive use of the random-mirror model. The latter model is similar to the random harmonic chain where a finite fraction of the atoms has an infinite mass, first considered by Domb et al. : quantities of interest are directly computable by elementary means. The results thus obtained still hold true for the even-visiting walk problem, the only difference being that the simple and explicit periodic amplitudes of the first model are replaced by unknown periodic functions. The same phenomenon is well established in the case of spectra of one-dimensional disordered systems: details of the distribution of random masses, random site potentials, and so on, only enter the periodic amplitudes of Lifshitz tails .

In the present case, just as in trapping problems , the quantities of most interest just involve the constant Fourier component of the model-specific periodic amplitudes. For the total number of even-visiting $`4k`$-step walks, the amplitude (6.12) has been determined very accurately by comparing the data of the exact enumeration procedure to the analytical asymptotic estimate (6.11).

Our analytical investigations have also led to the prediction (7.3) that all the cumulants of the maximum height $`M`$ reached by an $`n`$-step even-visiting walk scale as $`n^{1/3}`$, with known prefactors. Usual random walks (for which $`M/n^{1/2}`$ has a non-trivial limit law) and even-visiting ones (for which $`M/n^{1/3}`$ becomes more and more peaked as $`n`$ gets large) are therefore very different in that respect. We have also investigated amplitude ratios, associated with higher moments of the Riccati variables, which exhibit a highly non-trivial dependence in the variable $`\alpha `$.

Finally, the asymptotic result (6.11) demonstrates that the outcome (1.2) of Derrida’s original argument is very accurate in the one-dimensional case, as it just misses a power of the number of steps (in one dimension we have $`\mathrm{\Omega }=2`$ and $`j=\pi /2`$). The analogy with Lifshitz tails also suggests that the situation is less under control in higher dimension: going beyond the leading estimate (1.2), which involves only the volume of the Lifshitz sphere, indeed remains a difficult open problem .

### Acknowledgments

It is a pleasure for us to thank Bernard Derrida for interesting discussions, and Philippe Di Francesco and Emmanuel Guitter for making us aware of this problem, for useful discussions at early stages of this work, and for a critical reading of the manuscript.


text:
# Photocount statistics of chaotic lasers

## Abstract

We derive the photocount statistics of the radiation emitted from a chaotic laser resonator in the regime of single–mode lasing. Random spatial variations of the resonator eigenfunctions lead to strong mode–to–mode fluctuations of the laser emission. The distribution of the mean photocount over an ensemble of modes changes qualitatively at the lasing transition, and displays up to three peaks above the lasing threshold.

Nonlinear interactions in active optical media can drastically affect the statistics of the photon field. The textbook example is the single–mode laser : Below the laser threshold, losses outweigh the gain by linear amplification such that nonlinear saturation is negligible. Then the photon statistics is well represented by a thermal distribution, and the probability $`P_n`$ for the occupation of the mode with $`n`$ photons decays as the power law $`P_n[\overline{n}/(1+\overline{n})]^n`$ ($`\overline{n}`$ is the mean photon number). Above threshold nonlinear interactions stabilize the field intensity. The (relative) intensity fluctuations are strongly reduced below the value found for a thermal distribution, and far above threshold $`P_n`$ approaches the Poissonian distribution $`P_n=[\overline{n}^n/n!]\mathrm{exp}(\overline{n})`$ characteristic for a coherent state. These facts have been known since the development of the quantum theory of lasing in the 1960s.

Interest in the photon statistics of amplifying media was recently renewed by experiments on artificially fabricated random media. In random media radiation is scattered in an irregular, chaotic way. Practical realizations include laser dye solutions  and semiconductors  with randomly fluctuating dielectric constant as well as chaotic resonators with irregularly shaped boundaries . Even in the absence of pumping random media show interesting interference effects such as coherent backscattering or localization of light . The theoretical investigation of the quantum optical properties of random media started  with the demonstration of deviations from blackbody radiation in the case of one dimensional scattering. In a pioneering work  Beenakker generalized these results to multidimensional chaotic scattering by establishing a general relationship between the photocount statistics of a linear random medium and its scattering matrix. Using that relationship, Beenakker and coworkers predicted an excess photon noise due to multiple scattering , computed the photon noise power spectrum , and investigated the effect of photon localization on the photocount statistics .

In contrast to the wealth of results for linear random amplifiers, little is known about the photon statistics of random media above the lasing threshold. This is an important open problem since the experimental signatures of random lasing are still controversial . In the present paper we address that problem for a chaotic laser resonator in the regime of single–mode lasing. We show that the chaotic nature of the cavity modes gives rise to fluctuations of the photocount on top of the quantum optical fluctuations known from laser theory. Chaos–induced fluctuations are found when a single–mode photodetection is performed over an ensemble of modes. The ensemble may be obtained from a single resonator upon varying suitable parameters or from different resonators with small variations in shape. The factorial moments of the photocount display strong fluctuations from one mode to another. However, we show that the probability density for each factorial moment depends only on general symmetries of the ensemble and on four parameters describing the laser dynamics; these are the coefficients of linear gain, mean escape loss, absorption loss and saturation of the amplifying medium.

A sketch of the system is shown in Fig. 1: A chaotic resonator with an irregularly shaped boundary is homogeneously filled with an active medium. A partly transmitting mirror reflects impinging radiation with mean probability $`R`$ back into the resonator. The fraction $`T=1R`$ of radiation is transmitted and injected into a waveguide supporting $`M`$ transverse modes. A frequency selective photodetector counts the transmitted photons with efficiency $`1`$. We consider the case of a quasi discrete cavity spectrum; then the mean cavity escape rate $`\overline{\mathrm{\Gamma }}=MT\mathrm{\Delta }\omega `$ is smaller than the spacing $`\mathrm{\Delta }\omega `$ of adjacent resonator modes. Accordingly, the mirror transmission coefficient $`T`$ must be small, $`MT1`$. The counting time $`t`$ will be assumed large enough that the radiation from individual cavity modes can be resolved, $`t\mathrm{\Delta }\omega 1`$.

As for the active medium, we allow for $`N`$ active atoms. The $`\mu `$th atom “sees” the field mode through a coupling constant $`g_\mu `$ which is proportional to both the atomic dipole matrix element and the value of the lasing mode at the location of the atom. We only consider the simplest situation in which the characteristic times for atomic pump and losses are short compared to the mean life time of a photon in the cavity. The effect of the atoms on the field mode can then be represented by three parameters $`𝒜`$, $``$, $`𝒞`$ that characterize the linear gain, the nonlinear saturation, and the total loss, respectively . The first two of these depend on the coupling constants, $`𝒜_{\mu =1}^Ng_\mu ^2`$, $`_{\mu =1}^Ng_\mu ^4`$. Inasmuch as the resonator mode represents wave chaos, it varies quasi-randomly on the scale of the optical wavelength $`\lambda `$ . Coupling constants $`g_\mu `$ associated with atoms separated by more than $`\lambda `$ therefore behave as independent random numbers. The coefficients $`𝒜`$ and $``$, however, become sharp in the limit $`N1`$, due to the central limit theorem. The total loss rate $`𝒞=\mathrm{\Gamma }+\kappa `$ is the sum of the photon escape rate $`\mathrm{\Gamma }`$ and the absorption rate $`\kappa `$. Here, the absorption rate may be considered fixed while the escape rate is a random quantity, as will be revealed presently.

To calculate the photocount statistics we use the input–output theory of Gardiner and Collet . Waves entering and leaving the waveguide are described by $`M`$–component annihilation and creation operators $`a^{\mathrm{in}}`$, $`a^{\mathrm{in}}`$, $`a^{\mathrm{out}}`$, $`a^{\mathrm{out}}`$, which obey the commutation relations

$`[a_p(t),a_q^{}(t^{})]=\delta _{pq}\delta (tt^{}),[a_p(t),a_q(t^{})]`$ $`=`$ $`0.`$ (1)
Here, $`p,q=1,\mathrm{},M`$ label the transverse modes and $`a=a^{\mathrm{in}}`$ or $`a=a^{\mathrm{out}}`$. The boundary condition

$`a_p^{\mathrm{in}}(t)+a_p^{\mathrm{out}}(t)=\gamma _pb(t),`$ (2)
connects incoming and outgoing radiation in each transverse mode with the annihilation operator $`b`$ of the cavity mode. While the boundary condition (2) results from assuming a linear coupling between the waveguide and the cavity field, no restriction is imposed on the intracavity dynamics. Therefore, the input-output relation (2) holds both below and above threshold. The outcoupling amplitudes $`\gamma _p`$ must be considered independent random quantities with Gaussian statistics, since they represent the local fluctuations of the resonator mode across the outcoupling mirror. In fact, their randomness constitutes the principle effect of the wave chaos within the resonator on the ensemble fluctuations of the laser output.

In the linear regime below threshold, one can eliminate the cavity operators from Eq. (2) and express the outgoing radiation in terms of the incoming radiation, the intracavity noise, and the S–matrix . One thus obtains the photocount fluctuations  through the S–matrix statistics known for chaotic scattering.

We proceed to the nonlinear regime near and above threshold. For the case of a vacuum input, the boundary condition (2) allows to express all variances of the output field entirely in terms of the variances of the cavity field and the Gaussian statistics of the $`\gamma _p`$. We characterize the photocount statistics through the factorial moments

$`\mu _r={\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}m(m1)\mathrm{}(mr+1)p(m)`$ (3)
of the probability $`p(m)`$ that $`m`$ photons are counted in a time interval $`t`$. That probability is given by

$`p(m)`$ $`=`$ $`{\displaystyle \frac{1}{m!}}:W^me^W:,`$ (4)
$`W`$ $`=`$ $`{\displaystyle _0^t}𝑑t^{}{\displaystyle \underset{p=1}{\overset{M}{}}}a_p^{\mathrm{out}}(t^{})a_p^{\mathrm{out}}(t^{}).`$ (5)
Here $`\mathrm{}`$ denote the quantum steady–state average and the colons demand normal and a certain time ordering. Combining Eqs. (3)-(5) with the input–output relation (2) for a vacuum input we get

$`\mu _r=:W^r:=\mathrm{\Gamma }^r:I^r:,`$ (6)
where $`\mathrm{\Gamma }=_{p=1}^M|\gamma _p|^2`$ is the escape rate and $`I^r`$ the $`r`$th power of the integrated cavity field intensity

$`I={\displaystyle _0^t}𝑑t^{}b^{}(t^{})b(t^{}).`$ (7)
Clearly, the first moment alias mean photocount, $`\mu _1=\mathrm{\Gamma }tb^{}b`$, is a purely static quantity as it is proportional to the stationary mean photon number inside the cavity.

The factorial moments (6) specify the output statistics of a single cavity mode with the random escape rate $`\mathrm{\Gamma }`$. The distribution of $`\mathrm{\Gamma }`$ over an ensemble of modes is the $`\chi _\nu ^2`$ distribution

$`P(\mathrm{\Gamma })=A_\nu \mathrm{\Gamma }^{\nu /21}\mathrm{exp}(\nu \mathrm{\Gamma }/2\overline{\mathrm{\Gamma }})`$ (8)
well–known from random–matrix theory . Here, $`\nu =\beta M`$ is an integer and $`A_\nu `$ a normalization constant. The value of the parameter $`\beta `$ depends on whether the system is time–reversal invariant ($`\beta =1`$) or whether time–reversal invariance is broken ($`\beta =2`$) . The special case $`M=\beta =1`$ is known as the Porter–Thomas distribution. Together with $`\mathrm{\Gamma }`$ the factorial moments $`\mu _r`$ become random quantities. Their distribution is given by

$`𝒫(\mu _r)={\displaystyle }d\mathrm{\Gamma }P(\mathrm{\Gamma })\delta (\mu _r\mathrm{\Gamma }^r:I^r:).`$ (9)
Note that the right hand side involves a twofold average, the quantum optical average (represented by the brackets $`\mathrm{}`$) and the ensemble average over the cavity modes (represented by the integral with the probability distribution $`P(\mathrm{\Gamma })`$). We emphasize that the quantum optical average $`:I^r:`$ depends on $`\mathrm{\Gamma }`$ through the total loss rate $`𝒞`$. We now discuss the result (9) in various limiting cases.

We first consider the case $`M1`$ of many transverse modes in the waveguide. The diameter of the waveguide is then much larger than the wavelength of the cavity mode. A simple saddle point argument shows that $`P(\mathrm{\Gamma })`$ for large $`M`$ approaches a Gaussian distribution with mean $`\overline{\mathrm{\Gamma }}M`$ and standard deviation $`\mathrm{\Delta }\mathrm{\Gamma }=\overline{\mathrm{\Gamma }}/\sqrt{\beta M/2}`$. The relative fluctuations are small, $`\mathrm{\Delta }\mathrm{\Gamma }/\overline{\mathrm{\Gamma }}1/\sqrt{M}0`$ for $`M\mathrm{}`$, and the same is true for the fluctuations of all factorial moments. For cavities with large outcoupling mirrors we thus recover the sharp factorial moments one is used to from non-chaotic resonators.

Second, we investigate the limit of vanishing photocount, $`\mu _r0`$. According to Eq. (6), this is the weak–coupling limit $`\mathrm{\Gamma }0`$ for which the photons in the cavity field can hardly escape into the waveguide. The laser dynamics becomes independent of the outcoupling loss as the total cavity loss $`𝒞`$ is fully dominated by the absorption loss. All cavity field moments become independent of $`\mathrm{\Gamma }`$, and Eq. (6) reduces to $`\mu _r\mathrm{\Gamma }^r`$. Substitution into Eq. (9) yields the power–law behavior

$`𝒫(\mu _r)\mu _r^{(\beta M)/(2r)1},`$ (10)
for $`\mu _r0`$. A special case of this result is $`M=\beta =1`$ for which the distribution of the first and second moment diverge as $`\mu _1^{1/2}`$ and $`\mu _2^{3/4}`$ for $`\mu _1`$, $`\mu _20`$, respectively.

The third case is the short–time regime $`tt_c`$ where $`t_c`$ is the correlation time of the intensity fluctuations of the cavity mode. During counting intervals that short the field intensity cannot vary appreciably. Therefore, the integrated intensity $`I`$ becomes the product of $`t`$ with the instantaneous photon number $`\widehat{n}=b^{}b`$ in the cavity. The factorial moments of the output take the simple form

$`\mu _r=(\mathrm{\Gamma }t)^r:\widehat{n}^r:,`$ (11)
and one easily verifies that the normally ordered moment $`:\widehat{n}^r:`$ occuring here reduces to the $`r`$th factorial moment of the stationary photon number distribution of the lasing mode. The latter has the form

$`P_n=𝒩{\displaystyle \frac{[𝒜n_s/𝒞]^{n+n_s}}{(n+n_s)!}}.`$ (12)
Note that $`P_n`$ depends on the three laser parameters $`𝒜`$, $``$, and $`𝒞`$; in particular, the nonlinearity $``$ enters through the so–called saturation photon number $`n_s=𝒜/`$. The symbol $`𝒩`$ represents a normalization constant.

The calculation of the mean photocount $`\mu _1`$ reduces to a steady–state average for arbitrary counting time. We calculated $`\mu _1`$ numerically as a function of $`\mathrm{\Gamma }`$ using the distribution (12). Substituting the result $`\mu _1(\mathrm{\Gamma })`$ into Eq. (9) and carrying out the integration over $`\mathrm{\Gamma }`$, we obtained $`𝒫(\mu _1)`$. The results for time–reversal invariant cavities connected to a single–mode waveguide are plotted in Fig. 2, for $`\beta M=1`$ and four different sets of parameters. The two distributions in Fig. 2(a) correspond to $`𝒜>\overline{𝒞}`$, i.e. to lasers above threshold in the ensemble average. By contrast, $`𝒜<\kappa `$ for the distributions of Fig. 2(b); all lasers of those ensembles are below threshold irrespective of the escape rate $`\mathrm{\Gamma }`$. We note that all distributions are strongly non–Gaussian. They are all peaked as $`\mu _1^{1/2}`$ at small photocount, in accord with the general argument presented above for the asymptotics at $`\mu _10`$. Two further features spring to the eye and demand explanation: First, above threshold but not below we encounter an additional peak at maximum photocount; second, for certain cases (dashed lines in Fig. 2) the distribution $`𝒫(\mu _1)`$ displays a shoulder for sub–maximal $`\mu _1`$.

The origin of these structures lies in the $`\mathrm{\Gamma }`$–dependence of the mean photocount $`\mu _1`$. That dependence is depicted in the insets and seen to be qualitatively different above and below threshold. While $`\mu _1`$ increases monotonically with $`\mathrm{\Gamma }`$ in the below–threshold case, it develops a maximum at an intermediate value of $`\mathrm{\Gamma }`$ when the laser is above threshold. This behavior can be understood from the simple analytic expressions

$`\mu _1/t=\{\begin{array}{cc}\mathrm{\Gamma }n_s\frac{𝒜𝒞}{𝒞}\hfill & \text{for }𝒞=\mathrm{\Gamma }+\kappa <𝒜,\hfill \\ \frac{\mathrm{\Gamma }𝒜}{𝒞𝒜}\hfill & \text{for }𝒞=\mathrm{\Gamma }+\kappa >𝒜.\hfill \end{array}`$ (13)
that follow from the photon number distribution (12) sufficiently far from threshold. According to Eq. (13), $`\mu _1`$ rises linearly with $`\mathrm{\Gamma }`$ out of the origin and approaches an asymptotic plateau for very large $`\mathrm{\Gamma }`$, as visible in Fig. 2. The rise to the plateau is monotonic when $`𝒜<\kappa `$ since the below-threshold case of (13) applies for all values of $`\mathrm{\Gamma }`$; the maximum photocount $`\mu _{\mathrm{max}}`$ is then just the plateau value. The laser of Fig. 2(a) is above threshold for small $`\mathrm{\Gamma }`$ but below for large $`\mathrm{\Gamma }`$. The maximum photocount $`\mu _{\mathrm{max}}`$ then arises for an intermediate value $`\mathrm{\Gamma }^{}`$. A simple argument can be employed to determine the border between the below– and the above–threshold ensemble. The argument follows from the observation that $`\mu _1(\mathrm{\Gamma })`$ approaches its plateau–value from above for the above–threshold ensemble while the plateau is approached from below when the ensemble is below threshold. From the lower case of Eq. (13), which becomes exact for $`\mathrm{\Gamma }\mathrm{}`$, this yields the threshold condition $`𝒜=\kappa `$. To estimate the value of $`\mathrm{\Gamma }^{}`$ we may employ the above–threshold case of Eq. (13) and find $`\mathrm{\Gamma }^{}=\sqrt{𝒜\kappa }\kappa `$.

Based on this understanding of $`\mu _1(\mathrm{\Gamma })`$ one can appreciate the above–threshold peak of $`𝒫(\mu _1)`$ at $`\mu _1=\mu _{\mathrm{max}}`$. Substituting $`\mu _1(\mathrm{\Gamma })=\mu _{\mathrm{max}}+\frac{1}{2}\mu ^{\prime \prime }(\mathrm{\Gamma }^{})(\mathrm{\Gamma }\mathrm{\Gamma }^{})^2`$ in the photocount distribution (9) we find that $`𝒫(\mu _1)`$ has a square–root singularity $`|\mu _1\mu _{\mathrm{max}}|^{1/2}`$ which is precisely the peak depicted in Fig. 2(a). Clearly, no such peak can arise in the below–threshold case of Fig. 2(b) as the photocount increases monotonically with $`\mathrm{\Gamma }`$. Inspecting Eqs. (8), (9) in the limit $`\mathrm{\Gamma }\mathrm{}`$ we rather find that $`𝒫(\mu _1)`$ vanishes when $`\mu _1`$ approaches $`\mu _{\mathrm{max}}`$.

The shoulders at sub–maximal photocount are caused by amplified spontaneous emission below the laser threshold. Formally, the shoulders arise from the asymptotic plateau of $`\mu _1`$ for large $`\mathrm{\Gamma }`$. The definition (9) immediately implies enhanced probability for photocounts $`\mu _1`$ in the vicinity of that plateau. Note that the shoulder is invisible for the distribution shown as a solid line in Fig. 2(a), since for this distribution large values of $`\mathrm{\Gamma }`$ are strongly suppressed by the small mean value $`\overline{\mathrm{\Gamma }}`$. Further note that the shoulders for both curves in Fig. 2(b) lie closer to $`\mu _{\mathrm{max}}`$ than for the dashed curve in Fig. 2(a) since $`\mu \mu _{\mathrm{max}}`$ coincides with $`\mathrm{\Gamma }\mathrm{}`$ in the regime below threshold.

In contrast to the mean photocount which can be expressed in terms of the stationary distribution of the laser, all $`\mu _r`$ with $`r2`$ involve dynamical information through correlation functions with $`r1`$ time arguments. We defer the discussion of the higher factorial moments and their distribution to a separate publication .

We finally compare our results with related fluctuation phenomena in other areas of physics and discuss possible experimental tests. In nuclear physics the Porther–Thomas distribution describes level–width fluctuations in neutron scattering . The amplitude fluctuations of Coulomb blockade oscillations in semiconductor quantum dots are also of the Porther–Thomas type . In both cases, as well as for the chaotic lasers studied in this paper, the fluctuations result from the chaotic nature of wave functions. However, in chaotic lasers new interesting features arise due to the interplay of wave chaos with the nonlinear dynamics of the laser. As a consequence of this interplay, the distribution of the mean photocount can strongly deviate from the Porter–Thomas distribution. To test the predicted mode–to–mode fluctuations experimentally, one must study the photocount statistics for an ensemble of chaotic modes. It seems feasible to generate such ensembles in tunable lasers e.g. by shape variations, or by the injection and displacement of artificial scatterers in the case of microwave cavities.

Support by the Sonderforschungsbereich “Unordnung und große Fluktuationen” der Deutschen Forschungsgemeinschaft is gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

The thermodynamics of one-dimensional solvable models is generally determined by the solution to a set of socalled thermodynamic Bethe ansatz (TBA) equations . Some lattice spin models such as the XXZ chain, XYZ chain have been treated also by the quantum transfer matrix (QTM) method , see also chapters 17 and 18 of . Correlated electron systems such as t-J model and Hubbard model, have also been treated by TBA and QTM methods .

The equations obtained by the QTM approach are quite different from those of TBA. However, the numerical results of the two methods for the free energies are the same. Mathematically the non-linear integral equations of  and  share similarities insofar as they can be interpreted as equations for dressed energies of elementary particles of magnon and spinon type, respectively.

Recently, from the stand point of TBA one of the authors (MT)  derived in the case of the XXZ chain a simple integral equation for just one unknown function. This integral equation is completely different in structure from those mentioned above. Here we aim at a derivation of this equation in the Quantum Transfer Matrix approach providing a more explicit as well as unified understanding of the structures and involved functions.

To be definite, we first consider the region $`\mathrm{\Delta }1`$,

$``$ $`=`$ $`J{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{S_i^xS_{i+1}^x+S_i^yS_{i+1}^y+\mathrm{\Delta }(S_i^zS_{i+1}^z{\displaystyle \frac{1}{4}})\right\}2h{\displaystyle \underset{i=1}{\overset{N}{}}}S_i^z.`$ (1)
The thermodynamic Bethe ansatz equations for this model at temperature $`T`$ are called Gaudin-Takahashi equations,

$`\mathrm{ln}\eta _1(x)={\displaystyle \frac{2\pi J\mathrm{sinh}\varphi }{T\varphi }}𝐬(x)+𝐬\mathrm{ln}(1+\eta _2(x)),`$
$`\mathrm{ln}\eta _j(x)=𝐬\mathrm{ln}(1+\eta _{j1}(x))(1+\eta _{j+1}(x)),j=2,3,\mathrm{},`$
$`\underset{l\mathrm{}}{lim}{\displaystyle \frac{\mathrm{ln}\eta _l}{l}}={\displaystyle \frac{2h}{T}}.`$ (2)
Here we put

$$\mathrm{\Delta }=\mathrm{cosh}\varphi ,Q\pi /\varphi ,𝐬(x)=\frac{1}{4}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{sech}\left(\frac{\pi (x2nQ)}{2}\right),𝐬f(x)_Q^Q𝐬(xy)f(y)dy.$$
(3)
The free energy per site is

$$f=\frac{2\pi J\mathrm{sinh}\varphi }{\varphi }_Q^Q𝐚_1(x)𝐬(x)dxT_Q^Q𝐬(x)\mathrm{ln}(1+\eta _1(x))dx,𝐚_1(x)\frac{\varphi \mathrm{sinh}\varphi /(2\pi )}{\mathrm{cosh}\varphi \mathrm{cos}(\varphi x)}.$$
(4)
From this equation MT  derived

$`u(x)=2\mathrm{cosh}({\displaystyle \frac{h}{T}})+{\displaystyle _C}{\displaystyle \frac{\varphi }{2}}(\mathrm{cot}{\displaystyle \frac{\varphi }{2}}[xy2i]\mathrm{exp}[{\displaystyle \frac{2\pi J\mathrm{sinh}\varphi }{T\varphi }}𝐚_1(y+i)]`$
$`+\mathrm{cot}{\displaystyle \frac{\varphi }{2}}[xy+2i]\mathrm{exp}[{\displaystyle \frac{2\pi J\mathrm{sinh}\varphi }{T\varphi }}𝐚_1(yi)]){\displaystyle \frac{1}{u(y)}}{\displaystyle \frac{\mathrm{d}y}{2\pi i}},`$ (5)
where the free energy is given by

$$f=T\mathrm{ln}u(0).$$
(6)
The contour $`C`$ is an arbitrary closed loop counterclockwise around $`0`$ where $`2nQ,n0`$ and $`\pm 2i+2nQ`$ should lie outside of this loop. Furthermore this loop should not contain zeros of $`u(y)`$. It is expected that $`u(y)`$ has no zero in the region $`|\mathrm{}y|1`$. We show that these equations can be derived in the quantum transfer matrix approach.

## 2 Quantum transfer matrix and fusion hierarchy models

The quantum transfer matrix for this model is equivalent to that of the diagonal-to-diagonal transfer-matrix of the six-vertex model which is a staggered or inhomogeneous row-to-row transfer matrix, see below. The partition function $`Z\mathrm{Tr}\mathrm{exp}(/T)`$ is given by

$$Z=\underset{\{\sigma \}}{}\underset{j=1}{\overset{N}{}}\underset{i=1}{\overset{M}{}}A(\sigma _{2i+j,j}\sigma _{2i+j+1,j};\sigma _{2i+j,j+1}\sigma _{2i+j+1,j+1}),$$

$$A(\sigma _1\sigma _2;\sigma _1^{}\sigma _2^{})=\left[\begin{array}{cccc}a& 0& 0& 0\\ 0& c& b^{}& 0\\ 0& b& c& 0\\ 0& 0& 0& a\end{array}\right].$$

$`a`$ $`=\mathrm{exp}({\displaystyle \frac{J\mathrm{\Delta }}{2MT}})\mathrm{sinh}({\displaystyle \frac{J}{2MT}}),b=\mathrm{exp}({\displaystyle \frac{h}{MT}}),`$
$`b^{}`$ $`=\mathrm{exp}({\displaystyle \frac{h}{MT}}),c=\mathrm{exp}({\displaystyle \frac{J\mathrm{\Delta }}{2MT}})\mathrm{cosh}({\displaystyle \frac{J}{2MT}}).`$ (7)
Then in the case $`N=2M\times \mathrm{integer}`$, we have

$$Z=\mathrm{Tr}𝐓^N,𝐓(\sigma _1,\sigma _2,\mathrm{},\sigma _{2M};\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{2M}^{})$$

$$A(\sigma _1\sigma _2;\sigma _{2M}^{}\sigma _1^{})A(\sigma _3\sigma _4;\sigma _2^{}\sigma _3^{})\mathrm{}A(\sigma _{2M1}\sigma _{2M};\sigma _{2M2}^{}\sigma _{2M1}^{}).$$

The eigenvalue problem of this transfer matrix is a special case of the inhomogeneous six-vertex model on the square lattice.

Consider an inhomogeneous six-vertex model with the following column dependent Boltzmann weights:

$`a_l=\rho _l𝗁(v+v_l+\eta )`$
$`b_l=\rho _l\omega ^1𝗁(v+v_l\eta )`$
$`b_l^{}=\rho _l\omega 𝗁(v+v_l\eta )`$
$`c_l=\rho _l𝗁(2\eta ),l=1,\mathrm{},L.`$ (8)
Here $`L`$ is the number of columns, $`𝗁(u)`$ is $`u`$, $`\mathrm{sin}(u)`$ or $`\mathrm{sinh}(u)`$ depending on the anisotropy parameter. The transfer matrix $`𝐓(v)`$ acts in a $`2^L`$ dimensional space,

$`𝐓=\mathrm{Tr}𝐑_1(\sigma _1,\sigma _1^{})𝐑_2(\sigma _2,\sigma _2^{})\mathrm{}𝐑_L(\sigma _L,\sigma _L^{}),`$
$`𝐑_l(++)=\left(\begin{array}{cc}a_l& 0\\ 0& b_l\end{array}\right),𝐑_l(+)=\left(\begin{array}{cc}0& 0\\ c_l& 0\end{array}\right),`$
$`𝐑_l(+)=\left(\begin{array}{cc}0& c_l\\ 0& 0\end{array}\right),𝐑_l()=\left(\begin{array}{cc}b_l^{}& 0\\ 0& a_l\end{array}\right).`$ (9)
The space is divided into subspaces characterised by the number of down spins $`k`$. Without loss of generality we can put $`kL/2`$. In this subspace we can construct Bethe-ansatz wave functions with $`k`$ parameters $`u_1,\mathrm{},u_k`$,

$`|\mathrm{\Psi }={\displaystyle f(y_1,y_2,\mathrm{},y_k)\sigma _{y_1}^{}\sigma _{y_2}^{}\mathrm{}\sigma _{y_k}^{}|0},`$
$`f(y_1,y_2,\mathrm{},y_k)={\displaystyle \underset{P}{}}A(P){\displaystyle \underset{j=1}{\overset{k}{}}}F(y_j;u_{Pj}),`$
$`F(y;u)\omega ^y{\displaystyle \underset{l=1}{\overset{y1}{}}}𝗁(u+v_l+\eta ){\displaystyle \underset{l=y+1}{\overset{L}{}}}𝗁(u+v_l\eta ),`$
$`A(P)=ϵ(P){\displaystyle \underset{j<l}{}}𝗁(u_{Pj}u_{Pl}2\eta ).`$ (10)
Imposing periodic boundary conditions the Bethe ansatz equations (BAE) take the form

$`{\displaystyle \frac{\phi (u_j+\eta )}{\phi (u_j\eta )}}=\omega ^L{\displaystyle \underset{m=1}{\overset{k}{}}}{\displaystyle \frac{𝗁(u_ju_m+2\eta )}{𝗁(u_ju_m2\eta )}},`$
$`\phi (v)={\displaystyle \underset{l=1}{\overset{L}{}}}\rho _l𝗁(v+v_l).`$ (11)
The corresponding eigenvalue of the transfer matrix is given by

$`\mathrm{T}_1(v)=\omega ^{L+k}\phi (v\eta ){\displaystyle \frac{\mathrm{Q}(v+2\eta )}{\mathrm{Q}(v)}}+\omega ^k\phi (v+\eta ){\displaystyle \frac{\mathrm{Q}(v2\eta )}{\mathrm{Q}(v)}},`$
$`\mathrm{Q}(v)={\displaystyle \underset{j=1}{\overset{k}{}}}𝗁(vu_j).`$ (12)
In order to solve the diagonal-to-diagonal transfer matrix we have to consider an inhomogeneous six-vertex model the Boltzmann weights of which are given by

$`a_l=c_l=1,b_l=b_l^{}=0\mathrm{for}\mathrm{even}l,`$
$`a_l=\mathrm{exp}\left({\displaystyle \frac{J\mathrm{\Delta }}{2MT}}\right)\mathrm{sinh}\left({\displaystyle \frac{J}{2MT}}\right),b_l=\mathrm{exp}({\displaystyle \frac{h}{MT}}),`$
$`b_l^{}=\mathrm{exp}\left({\displaystyle \frac{h}{MT}}\right),c_l=\mathrm{exp}\left({\displaystyle \frac{J\mathrm{\Delta }}{2MT}}\right)\mathrm{cosh}\left({\displaystyle \frac{J}{2MT}}\right)\mathrm{for}\mathrm{odd}l.`$ (13)
The conditions (8) are satisfied if we put

$$L=2M,\omega =\mathrm{exp}\left(\frac{h}{MT}\right),v=0,\frac{𝗁^{}(2\eta )}{𝗁^{}(0)}=\frac{\mathrm{sinh}(\frac{J\mathrm{\Delta }}{MT})}{\mathrm{sinh}(\frac{J}{MT})},$$
(14)
and

$`\rho _l=1/𝗁(2\eta ),v_l=\eta \mathrm{for}\mathrm{even}l`$
$`\rho _l={\displaystyle \frac{\sqrt{bb^{}}}{𝗁(v_l\eta )}},{\displaystyle \frac{𝗁(v_l+\eta )}{𝗁(v_l\eta )}}={\displaystyle \frac{a}{\sqrt{bb^{}}}}=\mathrm{exp}\left({\displaystyle \frac{J\mathrm{\Delta }}{2MT}}\right)\mathrm{sinh}\left({\displaystyle \frac{J}{2MT}}\right)\mathrm{for}\mathrm{odd}l.`$ (15)
Putting $`\eta +v_1=2\alpha _M`$ we have

$$\phi (v)=\left(\frac{𝗁(v+\eta )𝗁(v+2\alpha _M\eta )}{𝗁(2\eta )𝗁(2\alpha _M2\eta )}\right)^M.$$
(16)
The largest eigenvalue belongs to the $`k=M`$ sector. The Bethe-ansatz equation for $`u_j,j=1,\mathrm{},M`$ are

$$\frac{\phi (u_j+\eta )}{\phi (u_j\eta )}=e^{2h/T}\underset{m=1}{\overset{M}{}}\frac{𝗁(u_ju_m+2\eta )}{𝗁(u_ju_m2\eta )}.$$
(17)
The corresponding eigenvalue is given by

$$\mathrm{T}_1(v)=e^{h/T}\phi (v\eta )\frac{\mathrm{Q}(v+2\eta )}{\mathrm{Q}(v)}+e^{h/T}\phi (v+\eta )\frac{\mathrm{Q}(v2\eta )}{\mathrm{Q}(v)}.$$
(18)
Due to the BAE (17), the eigenvalue $`\mathrm{T}_1(x)`$ is an entire function in the complex plane. The free energy per site is given by

$$f=T\underset{M\mathrm{}}{lim}\mathrm{ln}\mathrm{T}_1(0).$$
(19)
The matrix $`\mathrm{T}_1(v)`$ can be embedded into a more general family of matrices provided by the fusion hierarchy ,

$`\mathrm{T}_j(v)`$ $``$ $`{\displaystyle \underset{l=0}{\overset{j}{}}}e^{(j2l)h/T}\phi (v(j2l)\eta ){\displaystyle \frac{\mathrm{Q}(v+(j+1)\eta )\mathrm{Q}(v(j+1)\eta )}{\mathrm{Q}(v+(2lj+1)\eta )\mathrm{Q}(v+(2lj1)\eta )}}.`$ (20)
The eigenvalues $`\mathrm{T}_j(v)`$ as functions of $`v`$ are all entire in the complex plane. It is easily seen that the following functional relations hold

$`\mathrm{T}_j(v+\eta )\mathrm{T}_j(v\eta )=\phi (v+(j+1)\eta )\phi (v(j+1)\eta )+\mathrm{T}_{j+1}(v)\mathrm{T}_{j1}(v),`$
$`\mathrm{T}_0(v)\phi (v).`$ (21)
## 3  Derivation of Gaudin-Takahashi equation

For $`\mathrm{\Delta }>1`$ we put

$`𝗁(u)=\mathrm{sin}u,\eta =i\stackrel{~}{\varphi }/2,\stackrel{~}{\varphi }=\mathrm{cosh}^1\left({\displaystyle \frac{\mathrm{sinh}(J\mathrm{\Delta }/2MT)}{\mathrm{sinh}(J/2MT)}}\right),`$
$`\alpha _M={\displaystyle \frac{i}{2}}\mathrm{tanh}^1\left(\mathrm{tanh}\stackrel{~}{\varphi }\mathrm{tanh}{\displaystyle \frac{J\mathrm{\Delta }}{2MT}}\right),`$ (22)
In the limit of $`M\mathrm{}`$ we have

$$\stackrel{~}{\varphi }=\varphi ,M\alpha _M=iJ\mathrm{sinh}\varphi /(4T).$$
(23)
We transform the parameter $`v`$ to $`xiv/\eta `$. Then equations (16) and (20) turn into

$`\mathrm{Q}(x)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{M}{}}}\mathrm{sin}{\displaystyle \frac{\stackrel{~}{\varphi }}{2}}(xx_j),\phi (x)=\left({\displaystyle \frac{\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x+i)\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x(12u_M)i)}{\mathrm{sinh}\stackrel{~}{\varphi }\mathrm{sinh}\stackrel{~}{\varphi }(1u_M)}}\right)^M,`$
$`u_M`$ $`=`$ $`\alpha _M/\eta ,x_j=iu_j/\eta .`$ (24)
$`\mathrm{T}_j(x)`$ $``$ $`{\displaystyle \underset{l=0}{\overset{j}{}}}e^{(j2l)h/T}\phi (x(j2l)i){\displaystyle \frac{\mathrm{Q}(x+(j+1)i)\mathrm{Q}(x(j+1)i)}{\mathrm{Q}(x+(2lj+1)i)\mathrm{Q}(x+(2lj1)i)}}.`$ (25)
These functions are all entire in the complex plane. Now we introduce a modified eigenvalue of $`\mathrm{T}_j(x)`$

$`\stackrel{~}{\mathrm{T}}_j(x)`$ $``$ $`\mathrm{T}_j(x)\left({\displaystyle \frac{\mathrm{sinh}(\stackrel{~}{\varphi })\mathrm{sinh}\stackrel{~}{\varphi }(1u_M)}{\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x+(j+1)i)\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x(j+12u_M)i)}}\right)^M.`$ (26)
In contrast to the entire function $`\mathrm{T}_j(x)`$, $`\stackrel{~}{\mathrm{T}}_j(x)`$ has poles of order $`M`$ at

$`x=2nQ+u_Mi\pm (1+ju_M)i`$. On the other hand, it has constant asymptotics

$$\stackrel{~}{\mathrm{T}}_j(\pm i\mathrm{})=\frac{\mathrm{sinh}(j+1)h/T}{\mathrm{sinh}h/T}.$$
(27)
From (21), we can find the following functional relation for $`\stackrel{~}{\mathrm{T}}_j(x)`$

$$\stackrel{~}{\mathrm{T}}_j(x+i)\stackrel{~}{\mathrm{T}}_j(xi)=𝐛_j(x)+\stackrel{~}{\mathrm{T}}_{j1}(x)\stackrel{~}{\mathrm{T}}_{j+1}(x),$$
(28)
where we have defined

$`𝐛_j(x)=\left({\displaystyle \frac{\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x+(j+2u_M)i)\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(xji)}{\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x+ji)\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x(j2u_M)i)}}\right)^M.`$ (29)
Note that $`\stackrel{~}{\mathrm{T}}_0(x)=1`$ and $`𝐛_j(x),\stackrel{~}{\mathrm{T}}_j(x)`$ has poles at $`x=2nQ+u_Mi\pm (ju_M)i`$ and

$`x=2nQ+u_Mi\pm (j+1u_M)i`$, respectively.

We define

$$Y_j(x)=\frac{\stackrel{~}{\mathrm{T}}_{j1}(x)\stackrel{~}{\mathrm{T}}_{j+1}(x)}{𝐛_j(x)},j=1,2,\mathrm{}.$$
(30)
For these functions the following relations stand

$`Y_1(xi)Y_1(x+i)=1+Y_2(x),`$
$`Y_j(x+i)Y_j(xi)=(1+Y_{j1}(x))(1+Y_{j+1}(x)),j=2,3,\mathrm{},`$
$`\underset{l\mathrm{}}{lim}{\displaystyle \frac{\mathrm{ln}Y_l(x)}{l}}={\displaystyle \frac{2h}{T}}.`$ (31)
As $`Y_j(x),j=2,3,\mathrm{}`$ has no pole or zero in $`1\mathrm{}x1`$, we find

$$\mathrm{ln}Y_j(x)=𝐬(\mathrm{ln}(1+Y_{j1})+\mathrm{ln}(1+Y_{j+1})),j2.$$
(32)
For $`Y_1(x)`$ one must be careful that it has poles in $`i,(12u_M)i`$. Using

$$\stackrel{~}{\mathrm{T}}_2(x+i)\stackrel{~}{\mathrm{T}}_2(xi)=𝐛_2(x)(1+Y_2(x)),$$
(33)
and $`\stackrel{~}{\mathrm{T}}_2(x)`$ has no zero or pole at $`1\mathrm{}x1`$, we have

$$\mathrm{ln}\stackrel{~}{\mathrm{T}}_2(x)=𝐬(\mathrm{ln}𝐛_2(x)+\mathrm{ln}(1+Y_2(x))).$$
(34)
Using $`Y_1(x)=\stackrel{~}{\mathrm{T}}_2(x)/𝐛_1(x)`$ we have

$$\mathrm{ln}Y_1(x)=\mathrm{ln}𝐛_1(x)+𝐬\mathrm{ln}𝐛_2(x)+𝐬\mathrm{ln}(1+Y_2(x)).$$
(35)
In the limit of $`M\mathrm{}`$, the function $`𝐛_j(x)`$ can be simplified

$`𝐛_j(x)`$ $`=`$ $`\underset{M\mathrm{}}{lim}\mathrm{exp}\left[M\mathrm{ln}{\displaystyle \frac{\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x+(j+2u_M)i)\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(xji)}{\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x+ji)\mathrm{sin}\frac{\stackrel{~}{\varphi }}{2}(x(j2u_M)i)}}\right]`$ (36)
$`=`$ $`\mathrm{exp}\left({\displaystyle \frac{2\pi J\mathrm{sinh}\varphi }{\varphi T}}𝐚_j(x)\right),𝐚_j(x){\displaystyle \frac{\varphi \mathrm{sinh}j\varphi /(2\pi )}{\mathrm{cosh}j\varphi \mathrm{cos}(\varphi x)}},`$
which has singularities at $`x=2nQ\pm ji`$. In the limit of $`M\mathrm{}`$ equations (35),(32),(31) are identical to (2). Substituting

$$\mathrm{ln}\stackrel{~}{\mathrm{T}}_1(x)=𝐬\mathrm{ln}[(1+Y_1(x))/𝐛_1(x)]$$
(37)
into (19) we have (4). Then the Gaudin-Takahashi equations are derived from the quantum transfer matrix method. (See also the treatment in  for related models).

Consider the $`M\mathrm{}`$ limit of the functions $`𝐛_j(x),\stackrel{~}{\mathrm{T}}_j(x)`$ as

$$u_j(x)\underset{M\mathrm{}}{lim}\stackrel{~}{\mathrm{T}}_j(x),$$
(38)
Then from (28) we have the relation

$$u_1(x+i)u_1(xi)=𝐛_1(x)+u_2(x).$$
(39)
Note also the asymptotics $`u_1(\pm i\mathrm{})=2\mathrm{cosh}h/T`$. We may assume the functions $`u_1(x)`$ and $`u_2(x)`$ have similar singularities at $`x=2nQ\pm 2i`$ and $`x=2nQ\pm 3i`$, respectively. If we write (39) as

$$u_1(x+i)=𝐛_1(x)/u_1(xi)+u_2(x)/u_1(xi),$$
(40)
the LHS has singularities at $`x=i,3i`$ in the fundamental region ($`|\mathrm{}x|Q`$). The first term of the RHS has singularities at $`x=i,i,3i`$ and the second term at $`x=3i,3ii`$. Then following the method in , we get an integral equation for $`u_1(x)`$,

$`u_1(x)`$ $`=`$ $`2\mathrm{cosh}h/T`$
$`+`$ $`{\displaystyle _C}{\displaystyle \frac{\varphi }{2}}\left(\mathrm{cot}{\displaystyle \frac{\varphi }{2}}[xy2i]𝐛_1(y+i)+\mathrm{cot}{\displaystyle \frac{\varphi }{2}}[xy+2i]𝐛_1(yi)\right){\displaystyle \frac{1}{u_1(y)}}{\displaystyle \frac{\mathrm{d}y}{2\pi i}}.`$
From the explicit expression of $`u_1(x)`$ (36), we see that the integral equation (LABEL:eq.integral) is identical to the one obtained in . The free energy is given by

$$f=T\mathrm{ln}u_1(0).$$
(42)
## 4 Case $`\mathrm{\Delta }<1`$

In this case we have

$`𝗁(u)=\mathrm{sinh}u,\eta =i\stackrel{~}{\theta }/2,\stackrel{~}{\theta }=\mathrm{cos}^1\left({\displaystyle \frac{\mathrm{sinh}(J\mathrm{\Delta }/2MT)}{\mathrm{sinh}(J/2MT)}}\right),`$
$`\alpha _M={\displaystyle \frac{i}{2}}\mathrm{tanh}^1\left(\mathrm{tan}\stackrel{~}{\theta }\mathrm{tanh}{\displaystyle \frac{J\mathrm{\Delta }}{2MT}}\right).`$ (43)
In the limit of $`M\mathrm{}`$ we have

$$\stackrel{~}{\theta }=\mathrm{cos}^1\mathrm{\Delta },M\alpha _M=iJ\mathrm{sin}\theta /(4T).$$
(44)
Putting $`x=iv/\eta `$ we obtain

$$\mathrm{Q}(x)=\underset{j=1}{\overset{M}{}}\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(xx_j),\phi (x)=\left(\frac{\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x+i)\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x(12u_M)i)}{\mathrm{sin}\stackrel{~}{\theta }\mathrm{sin}\stackrel{~}{\theta }(1u_M)}\right)^M.$$
(45)
Kuniba, Sakai and Suzuki  succeeded in deriving the Takahashi-Suzuki equations  for the thermodynamics of the XXZ model at $`h=0,|\mathrm{\Delta }|<1`$. The functions

$`\stackrel{~}{\mathrm{T}}_j(x)`$ $``$ $`\mathrm{T}_j(x)\left({\displaystyle \frac{\mathrm{sin}(\stackrel{~}{\theta })\mathrm{sin}\stackrel{~}{\theta }(1u_M)}{\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x+(j+1)i)\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x(j+12u_M)i)}}\right)^M.`$ (46)
are all periodic with periodicity $`2p_0i`$. We have relations for $`\stackrel{~}{\mathrm{T}}_1(x)`$ and $`\stackrel{~}{\mathrm{T}}_2(x)`$

$$\stackrel{~}{\mathrm{T}}_1(x+i)\stackrel{~}{\mathrm{T}}_1(xi)=b_1(x)+\stackrel{~}{\mathrm{T}}_2(x),$$
(47)
with

$`b_1(x)=\left({\displaystyle \frac{\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x+(1+2u_M)i)\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(xi)}{\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x+i)\mathrm{sinh}\frac{\stackrel{~}{\theta }}{2}(x(12u_M)i)}}\right)^M.`$ (48)
$`\stackrel{~}{\mathrm{T}}_1(x)`$ satisfies

$$\stackrel{~}{\mathrm{T}}_1(\pm \mathrm{})=2\mathrm{cosh}h/T.$$
(49)
By these two equations we can determine $`\stackrel{~}{\mathrm{T}}_1(x)`$ in the limit of $`M\mathrm{}`$. In this limit $`b_1(x)`$ is

$$b_1(x)=\mathrm{exp}\left(\frac{2\pi J\mathrm{sin}\theta }{\theta T}a_1(x)\right),a_1(x)\frac{\theta \mathrm{sin}\theta /(2\pi )}{\mathrm{cosh}(\theta x)\mathrm{cos}\theta }.$$
(50)
We can assume that $`\stackrel{~}{\mathrm{T}}_1(x)`$ is expanded as follows

$$\stackrel{~}{\mathrm{T}}_1(x)=2\mathrm{cosh}(\frac{h}{T})+\underset{j=1}{\overset{\mathrm{}}{}}\underset{n}{}\frac{c_j}{(x2np_0i2i)^j}+\underset{j=1}{\overset{\mathrm{}}{}}\underset{n}{}\frac{\overline{c_j}}{(x2np_0i+2i)^j}.$$
(51)
Consider the contour integral around $`x=i`$ giving the coefficients $`c_j`$

$$c_j=\frac{(xi)^{j1}b_1(x)}{\stackrel{~}{\mathrm{T}}_1(xi)}\frac{\mathrm{d}x}{2\pi i}=\frac{y^{j1}b_1(y+i)}{\stackrel{~}{\mathrm{T}}_1(y)}\frac{\mathrm{d}y}{2\pi i}.$$
(52)
The first sum of the r.h.s. of (51) is

$`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n}{}\frac{b_1(y+i)}{(x2np_0i2i)^j}\frac{y^{j1}}{\stackrel{~}{\mathrm{T}}_1(y)}\frac{\mathrm{d}y}{2\pi i}}={\displaystyle \underset{n}{}\frac{b_1(y+i)}{xy2np_0i2i}\frac{1}{\stackrel{~}{\mathrm{T}}_1(y)}\frac{\mathrm{d}y}{2\pi i}}`$
$`={\displaystyle \frac{\theta }{2}\mathrm{coth}\frac{\theta }{2}(xy2i)\mathrm{exp}[\frac{2\pi J\mathrm{sin}\theta }{T\theta }a_1(y+i)]\frac{1}{\stackrel{~}{\mathrm{T}}_1(y)}\frac{\mathrm{d}y}{2\pi i}}.`$ (53)
The second sum is calculated in a similar way. Thus we find

$`u(x)=2\mathrm{cosh}({\displaystyle \frac{h}{T}})+{\displaystyle _C}{\displaystyle \frac{\theta }{2}}(\mathrm{coth}{\displaystyle \frac{\theta }{2}}[xy2i]\mathrm{exp}[{\displaystyle \frac{2\pi J\mathrm{sin}\theta }{T\theta }}a_1(y+i)]`$
$`+\mathrm{coth}{\displaystyle \frac{\theta }{2}}[xy+2i]\mathrm{exp}[{\displaystyle \frac{2\pi J\mathrm{sin}\theta }{T\theta }}a_1(yi)]){\displaystyle \frac{1}{u(y)}}{\displaystyle \frac{\mathrm{d}y}{2\pi i}},`$ (54)
and the free energy is given by

$$f=T\mathrm{ln}u(0).$$
(55)
Apparently these equations are analytical continuations of (5) and (6) if we replace $`\varphi `$ by $`i\theta `$. Then equation (5) treats the thermodynamics in a unified way.

This research is supported in part by Grants-in-Aid for the Scientific Research (B) No. 11440103 from the Ministry of Education, Science and Culture, Japan and Japan-German cooperation program of scientific research ”Theory of low-dimensional quantum spin systems and correlated electron systems” by JSPS and DFG.