text:
# Eigenfunction approach to the persistent random walk in two dimensions

## I Introduction

Today, the concept of the random walk is an ubiquitous tool in the physical sciences, and a multitude of different mathematical and physical problems have been subsumed under this catch-all phrase Shlesinger and West (1984); Weiss (1994). Originally, the term has been coined by Pearson Pearson (1905) who was interested in the probability distribution $`w_P(R)`$ for the end-to-end distance in a chain of $`N`$ randomly chosen steps of equal length $`l`$ in a plane. This article concerns a random walk process that is closely related to the initial formulation of the problem: While Pearson’s general setup is maintained, we introduce an angular bias into the relative orientation of successive steps that is governed by a probability distribution $`p(\theta )`$. Known as the continuous persistent random walk in the plane, this modification continues to attract considerable theoretical interest Weiss (1994); Masoliver et al. (1993); Tojo and Argyrakis (1996); Larralde (1997); Wu et al. (2000). Practical applications of the planar walk cover a wide range including light scattering, where sums of random complex numbers are studied Barakat (1974), its related use in crystal structure analysis Shmueli and Weiss (1990); Weiss (1994); Shmueli and Weiss (1995), and surface diffusion processes. The model of the persistent walk examined here is of particular interest to the configurational statistics of polymers Flory (1969), where the inherent stiffness of chemical bonds restricts the relative orientation of adjacent monomers, and obviously for two-dimensional locomotion problems in biology and medicine Nossal and Weiss (1974a, b); Nossal (1983); Wu et al. (2000).

In several of these publications Tojo and Argyrakis (1996); Larralde (1997); Wu et al. (2000); Flory (1969), only moments of the probability distribution function $`w(R,\varphi )`$ were calculated, while it appeared desirable to the authors to have access to $`w(R,\varphi )`$ beyond the diffusion limit $`N\mathrm{}`$, where this function ultimately takes on a simple Gaussian shape. Unfortunately, the search for explicit solutions even for Pearson’s original problem has remained largely elusive. While some recent contributors resorted to numerical Monte Carlo simulation Tojo and Argyrakis (1996), in a pioneering work Barakat Barakat (1973) realized that expansion into a discrete Fourier series offers a superior alternative in computing the radial distribution function $`w(R)`$ of isotropic random walks. The corresponding Fourier–Bessel series solution for the Pearson random walk is equally due to Barakat Barakat (1974), and subsequently has been advantageously applied in crystallographical analysis Shmueli and Weiss (1990); Weiss (1994). In this article, we extend the Fourier analysis to the persistent random walk problem characterized by an anisotropic end-to-end distribution function $`w(R,\varphi )`$ and find excellent numerical convergence for reasonably smooth angular distribution functions $`p(\theta )`$, except for the shortest walks. This renders the method particularly suited for the study of intermediate-size chains whose behaviour conspicuously deviates from the diffusion limit. Our derivation of the series representation rests on rather elementary methods. (A more complicated general analysis based on the characteristic function $`c(\omega ,\chi )`$ has been put forward by Weiss and Shmueli Weiss and Shmueli (1987). When applied to the persistent random walk, it implicitly yields the same result.)

Let us briefly outline the organization of this article: In Section II, we develop the theory of the persistent planar walk. Starting with moments of the distribution, we proceed to derive the characteristic function $`c(\omega ,\chi )`$ in form of a matrix product, which we use to deduce an integral representation of $`w(R,\varphi )`$ reminiscent of Kluyver’s solution to Pearson’s problem Kluyver (1906). Finally, we establish the Fourier–Bessel eigenfunction expansion of this distribution. Numerical results obtained with our method are presented in Section III, where we study the convergence behaviour of the series and display the frequently distinctive shapes of the distribution function obtained for several choices of the angular weight $`p(\theta )`$. Some possible extensions of the theory are discussed in Section IV. Finally, a proof of the “projection theorem” that forms the basis of our proposed series expansion is sketched in Appendix A, while Appendix B covers the calculation of moments in the Fourier–Bessel scheme.

## II Theory of the 2D persistent walk

We start out with some basic conventions. For simplicity of notation, in this paper we deal exclusively with a planar persistent random walk built of $`N`$ steps with uniform length $`l`$, but the generalization of our theory to incorporate individually different step lengths is straightforward. The walk is then fully characterized by the set of relative angles $`\{\theta _1,\theta _2,\mathrm{},\theta _N\}`$ at the links of successive segments, where the initial angle $`\theta _1`$ is measured relative to the $`x`$–axis of our coordinate system (Figure 1). In practice, the absolute angles $`\stackrel{~}{\theta }_\nu `$ formed between the $`\nu `$.th step and the $`x`$–axis as direction of reference are equally important, and they form an alternative description of the chain. As easily inferred from Figure 1, both sets of angles are interrelated via:

$$\stackrel{~}{\theta }_\nu \stackrel{~}{\theta }_{\nu 1}=\theta _\nu \stackrel{~}{\theta }_\nu =\underset{k=1}{\overset{\nu }{}}\theta _k.$$
(1)
Since we are mainly interested in the end-to-end distribution function $`w(R,\varphi )`$, we start the random walk from the origin of our coordinate system. We denote its endpoint by $`(L_x,L_y)`$ in cartesian and $`(R,\varphi )`$ in polar coordinates, respectively.

For the persistent walk, the probability of a certain configuration $`p(\theta _1,\theta _2,\mathrm{},\theta _N)`$ depends only on the relative angles $`\theta _\nu `$ formed by adjacent steps. Assuming that these angles are not correlated, and that all links are equivalent (including the initial “joint” between the first step and the $`x`$–axis), this probability reduces to a product of independent factors: $`p(\theta _1,\mathrm{},\theta _N)=p(\theta _1)\mathrm{}p(\theta _N)`$, where the angular probability distribution $`p(\theta )=p(\theta +2\pi )`$ is a periodic function of the change of direction $`\theta `$. This suggests that a Fourier expansion of the function $`p(\theta )`$ will be convenient Nossal and Weiss (1974a, b):

$$p(\theta )=\frac{1}{2\pi }\underset{\nu }{}\mathrm{e}^{\mathrm{i}\nu \theta }p_\nu p_\nu =_0^{2\pi }d\theta \mathrm{e}^{\mathrm{i}\nu \theta }p(\theta ).$$
(2)
Proper normalization of $`p(\theta )`$ requires $`p_0=1`$; the norm of all other (generally complex) angular Fourier components is then less than unity, $`|p_\nu |<1`$ for $`\nu 0`$. (Equality may occur only for rigidly linked chains with $`p(\theta )=\delta (\theta \theta _0)`$ which we exclude here.) With these definitions, we are now set to determine quantities of interest for the walk.

### II.1 Persistence length and diffusion constant

Our first investigation concerns the average position $`(L_x_N,L_y_N)`$ of the chain end and its mean square displacement $`L_x^2+L_y^2_N`$ as a function of the number of steps $`N`$. Both quantities are accessible in a simple calculation that introduces some techniques that will be commonplace below. It is convenient to interpret vectors in the $`xy`$–plane as complex numbers; with this convention, the vector formed by the $`\nu `$.th element $`(l_x^{(\nu )},l_y^{(\nu )})`$ is characterized by its absolute angle $`\stackrel{~}{\theta }_\nu `$ with respect to the $`x`$–axis (Figure 1): $`l_x^{(\nu )}+\mathrm{i}l_y^{(\nu )}=l\mathrm{e}^{\mathrm{i}\stackrel{~}{\theta }_\nu }`$. The position of the chain end is simply the sum of all constituent step vectors, $`L_x+\mathrm{i}L_y=l_{\nu =1}^N\mathrm{e}^{\mathrm{i}\stackrel{~}{\theta }_\nu }`$, and its properly weighed average over all configurations of a $`N`$–element walk correspondingly reads:

$`L_x+\mathrm{i}L_y_N`$ $`=l{\displaystyle \underset{\nu =1}{\overset{N}{}}}{\displaystyle _0^{2\pi }}d\theta _1p(\theta _1)\mathrm{}{\displaystyle _0^{2\pi }}d\theta _Np(\theta _N)\mathrm{e}^{\mathrm{i}\stackrel{~}{\theta }_\nu }`$
$`=l{\displaystyle \underset{\nu =1}{\overset{N}{}}}p_0^{N\nu }p_1^\nu =l{\displaystyle \frac{p_1\left(1p_1^N\right)}{1p_1}}.`$ (3)
Here, we successively used the resolution of the absolute angle $`\stackrel{~}{\theta }_\nu `$ into relative angles $`\theta _k`$ (1), the definition of the angular Fourier components in (2), and finally the normalization condition $`p_0=1`$. We infer that the average displacement of the chain end depends only on the first Fourier component of $`p(\theta )`$, and if $`p_1`$ vanishes, the distribution is centered around the origin for all $`N`$. Unless the chain is absolutely rigid ($`|p_1|=1`$), the mean end-to-end vector approaches a well-defined limit as $`N\mathrm{}`$:

$$\underset{N\mathrm{}}{lim}L_x+\mathrm{i}L_y_N=l\frac{p_1}{1p_1},$$
(4)
which generalizes the notion of the persistence length of the chain put forward in Ref. Flory (1969). For comparison, we also state this result using real parameters. Introducing $`p_1=\alpha \mathrm{e}^{\mathrm{i}\beta }`$, where $`\beta `$ represents the average “slant” angle in the angular distribution $`p(\theta )`$ Tojo and Argyrakis (1996); Larralde (1997), one finds in the long-chain limit:

$`L_x_{\mathrm{}}`$ $`={\displaystyle \frac{l\alpha (\mathrm{cos}\beta \alpha )}{12\alpha \mathrm{cos}\beta +\alpha ^2}},`$ (5)
$`L_y_{\mathrm{}}`$ $`={\displaystyle \frac{l\alpha \mathrm{sin}\beta }{12\alpha \mathrm{cos}\beta +\alpha ^2}},`$ (6)
which yields for the average absolute displacement $`R=(L_x^2+L_y^2)^{1/2}`$ as $`N\mathrm{}`$:

$$R_{\mathrm{}}=\frac{l\alpha }{\sqrt{12\alpha \mathrm{cos}\beta +\alpha ^2}}.$$
(7)
This simple expression might be useful in the interpretation of experimental data Wu et al. (2000).

The determination of the mean square displacement $`R^2_N=L_x^2+L_y^2_N`$ is only slightly more involved. In the complex notation, $`R^2=(L_x+\mathrm{i}L_y)(L_x\mathrm{i}L_y)=l^2_{\nu ,\mu =1}^N\mathrm{e}^{\mathrm{i}(\stackrel{~}{\theta }_\nu \stackrel{~}{\theta }_\mu )}`$, and the same series of arguments we used in (3) yields:

$`R^2_N`$ $`=l^2{\displaystyle \underset{\nu ,\mu =1}{\overset{N}{}}}{\displaystyle _0^{2\pi }}d\theta _1p(\theta _1)\mathrm{}{\displaystyle _0^{2\pi }}d\theta _Np(\theta _N)\mathrm{e}^{\mathrm{i}(\stackrel{~}{\theta }_\nu \stackrel{~}{\theta }_\mu )}`$
$`=Nl^2+l^2{\displaystyle \underset{\delta =1}{\overset{N}{}}}(N\delta )\left(p_1^\delta +p_1^\delta \right),`$ (8)
where we exploited $`\stackrel{~}{\theta }_\nu \stackrel{~}{\theta }_\mu =\theta _{\mu +1}+\mathrm{}+\theta _\nu `$ $`(\nu >\mu )`$ and sorted the terms in the double sum according to their “distance” $`\delta =|\nu \mu |`$. The arithmetic-geometric progression in (8) is easily evaluated Gradshteyn and Ryzhik (1980):

$$R^2_N=\frac{Nl^2}{2}\frac{1+p_1}{1p_1}l^2\frac{p_1(1p_1^N)}{(1p_1)^2}+\text{c.c.}$$
(9)
(A very similar result is obtained for the mean square displacement in the three-dimensional freely rotating chain Flory (1969).) Since the complex conjugate term is formally obtained by replacing $`p_1`$ with $`p_1`$ in (9), the mean square displacement, like the absolute displacement (3), depends only on the first angular Fourier component $`p_{\pm 1}`$. For $`p_1=0`$, (9) reduces to $`R^2_N=Nl^2`$, the trivial result obtained in Pearson’s walk. In the limit $`N\mathrm{}`$, we arrive at the following diffusion constant (characteristic ratio) $`D_N=R^2_N/N`$ for the persistent walk:

$$D_{\mathrm{}}=\frac{l^2(1p_1p_1)}{(1p_1)(1p_1)}=\frac{l^2(1\alpha ^2)}{12\alpha \mathrm{cos}\beta +\alpha ^2}.$$
(10)
The asymmetry parameters $`\alpha `$ and $`\beta `$ of the angular distribution $`p(\theta )`$ thus determine the persistence length $`R_{\mathrm{}}`$ (7) and the diffusion constant $`D_{\mathrm{}}`$ of the walk and vice versa. We point out that for $`\beta 0`$, with increasing $`\alpha `$ $`D_{\mathrm{}}`$ starts to diminish from a certain stiffness of the joint and actually vanishes as $`\alpha 1`$. This property reflects the fact that the walk then tends to move around a circle, and no net transport takes place. (An alternative derivation of (10) has been given by Larralde Larralde (1997).) We conclude this section by noting that higher momenta of the distribution function $`w(R,\varphi )`$ evidently are available from similar calculations. However, in practice they are more easily evaluated from a Fourier–Bessel series representation (Appendix B).

### II.2 The characteristic function

In the previous section, we were able to determine the moments $`L_x+\mathrm{i}L_y_N`$ (3) and $`R^2_N`$ (8) of the probability distribution $`w(R,\varphi )`$ by summation of some geometric series. We now condense the ideas implicitly introduced there into a single, coherent formalism and employ it towards the corresponding characteristic function $`c(\omega ,\chi )=\mathrm{exp}(\mathrm{i}𝝎𝐋)_N`$ of the persistent walk.

So far, we evaluated averages of the form $`\mathrm{exp}(\mathrm{i}\stackrel{~}{\theta }_\nu )`$ and $`\mathrm{exp}[\mathrm{i}(\stackrel{~}{\theta }_\nu \stackrel{~}{\theta }_\mu )]`$, respectively. Both expressions may be viewed as special cases of the more general average $`\mathrm{exp}(\mathrm{i}_\nu \lambda _\nu \stackrel{~}{\theta }_\nu )`$, where $`\{\lambda _1,\mathrm{},\lambda _N\}`$ is a set of $`N`$ integers and the sum covers the range $`1\nu N`$. In view of the relation $`\stackrel{~}{\theta }_\nu =\theta _1+\mathrm{}+\theta _\nu `$ (1) between absolute and relative angles, its expression in terms of the angular Fourier coefficients $`p_\lambda `$ (2) is straightforward:

$$\mathrm{exp}\left(\mathrm{i}_\nu \lambda _\nu \stackrel{~}{\theta }_\nu \right)=p_{\lambda _N}p_{\lambda _N+\lambda _{N1}}\mathrm{}p_{\lambda _N+\mathrm{}+\lambda _2+\lambda _1}.$$
(11)
The significance of this development is as follows. Consider a functional $`Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)`$ of the random walk that depends only on the absolute angles $`\stackrel{~}{\theta }_\nu `$ of the segments. In order to be meaningful, such a function must be periodic in each of those angles, $`Q(\mathrm{},\stackrel{~}{\theta }_\mu +2\pi ,\mathrm{})=Q(\mathrm{},\stackrel{~}{\theta }_\mu ,\mathrm{})`$, which in turn implies that it may be represented by a discrete Fourier series:

$$Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=\underset{\mu _1\mathrm{}\mu _N}{}Q_{\mu _1,\mathrm{},\mu _N}\mathrm{e}^{\mathrm{i}(\mu _1\stackrel{~}{\theta }_1+\mathrm{}+\mu _N\stackrel{~}{\theta }_N)},$$
(12)
where integration yields the Fourier coefficient $`Q_{\mu _1,\mathrm{},\mu _N}`$:

$$\begin{array}{c}Q_{\mu _1,\mathrm{},\mu _N}=\frac{1}{(2\pi )^N}_0^{2\pi }d\stackrel{~}{\theta }_1\mathrm{}_0^{2\pi }d\stackrel{~}{\theta }_N\hfill \\ \hfill \times \mathrm{e}^{\mathrm{i}(\mu _1\stackrel{~}{\theta }_1+\mathrm{}+\mu _N\stackrel{~}{\theta }_N)}p(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N).\end{array}$$
(13)
(It proves convenient here to employ a Fourier expansion scheme that differs from the angular decomposition (2) in its sign convention and normalization.) By combination of (11) and (13), we find the average value for the functional $`Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)`$ in terms of Fourier coefficients:

$`Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)`$ $`={\displaystyle \underset{\mu _1\mathrm{}\mu _N}{}}Q_{\mu _1,\mathrm{},\mu _N}p_{\mu _N}\mathrm{}p_{\mu _N+\mathrm{}+\mu _1}`$
$`=\begin{array}{cc}\hfill {\displaystyle \underset{\nu _1\mathrm{}\nu _N}{}}& Q_{\nu _1\nu _2,\nu _2\nu _3,\mathrm{},\nu _{N1}\nu _N,\nu _N}\hfill \\ & \times p_{\nu _N}p_{\nu _{N1}}\mathrm{}p_{\nu 1}.\hfill \end{array}`$ (14)
In the second line of this equation, we employed an alternative indexing scheme $`\{\nu _1,\mathrm{},\nu _N\}`$ defined via:

$$\nu _k=\mu _k+\mathrm{}+\mu _N\mu _k=\nu _k\nu _{k+1},$$
(15)
$`(\mu _N=\nu _N)`$. Its advantages will become obvious below.

We are particularly interested in the special case where the functional $`Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)`$ reduces to a product of $`N`$ equivalent factors $`q(\stackrel{~}{\theta }_\nu )`$: $`Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=_\nu q(\stackrel{~}{\theta }_\nu )`$. Under these circumstances, the Fourier components $`Q_{\mu _1,\mathrm{},\mu _N}`$ will factorize as well, and we have $`Q_{\mu _1,\mathrm{},\mu _N}=q_{\mu _1}\mathrm{}q_{\mu _N}`$, where:

$$q(\stackrel{~}{\theta })=\underset{\mu }{}\mathrm{e}^{\mathrm{i}\mu \stackrel{~}{\theta }}q_\mu q_\mu =\frac{1}{2\pi }_0^{2\pi }d\stackrel{~}{\theta }\mathrm{e}^{\mathrm{i}\mu \stackrel{~}{\theta }}q(\stackrel{~}{\theta }).$$
(16)
We now exploit these properties in order to simplify the average (14). Upon inserting (16), we rearrange the summations involved:

$$\begin{array}{c}Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=\underset{\nu _N}{}q_{\nu _N}p_{\nu _N}\underset{\nu _{N1}}{}q_{\nu _{N1}\nu _N}p_{\nu _{N1}}\hfill \\ \hfill \times \mathrm{}\underset{\nu _2}{}q_{\nu _2\nu _3}p_{\nu _2}\underset{\nu _1}{}q_{\nu _1\nu _2}p_{\nu _1}.\end{array}$$
(17)
Thus, we may rewrite $`Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)`$ in the form of a trace over the $`N`$.th power of a matrix $`𝒵`$:

$$Q(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=\underset{\nu _1}{}\left[𝒵^N\right]_{0\nu _1}=\widehat{e}_0^T𝒵^N𝐧,$$
(18)
where $`[\widehat{e}_0]_k=\delta _{0k}`$ denotes a unit vector, and $`𝐧=(\mathrm{},1,1,1,\mathrm{})^T`$ is composed of unit entries. The matrix $`𝒵`$ itself is defined via:

$$𝒵_{jk}=q_{kj}p_k,$$
(19)
and may be interpreted as the product of a circulant matrix $`𝒬_{jk}=q_{kj}`$ with a diagonal matrix $`𝒫_{jk}=p_j\delta _{jk}`$. (Occasionally, it is helpful to split off a scaling factor $`\lambda `$ from the coefficients $`q_\nu `$: $`\stackrel{~}{q}_\nu =\lambda ^\nu q_\nu `$, e. g. to accomodate functionals in a rotated frame of reference $`Q_\alpha (\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=Q(\stackrel{~}{\theta }_1\alpha ,\mathrm{},\stackrel{~}{\theta }_N\alpha )`$, where $`\stackrel{~}{q}_\nu =\mathrm{e}^{\mathrm{i}\nu \alpha }q_\nu `$ (16). A brief inspection of (17) shows that a slight generalization of (18) will cover this extension:

$$Q_\lambda (\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=\underset{\nu _1}{}\lambda ^{\nu _1}\left[𝒵^N\right]_{0\nu _1}.)$$
(20)
The expressions (18) and (20) become especially handy when only a finite number of angular Fourier components $`p_k`$ in (2) are non-vanishing. If $`p_\nu 0`$ for all $`|\nu |>m`$, only the coefficients $`q_{2m},\mathrm{},q_{2m}`$ will enter the calculation (irrespective of the nature of $`q(\stackrel{~}{\theta })`$), and the problem in (18) is effectively reduced to finding the $`N.`$th power of a $`(2m+1)\times (2m+1)`$–matrix. Eigenvalue methods then offer an efficient algorithm to calculate the average for any $`N`$. In conclusion, we remark that our matrix formulation of the planar persistent random walk problem closely resembles the treatment of polymer chains in the rotational isomeric state (RIS)–approximation Flory (1969), and bears similarity to the Ising model of a one-dimensional ferromagnetic system Kramers and Wannier (1941). The mathematical aspects of the representation of averages in Markoff chains through powers of a matrix have been summarized by Montroll Montroll (1947) who incidentally applied the method to discuss the probability distribution function for a persistent random walk in a single dimension.

As an important example for this formalism we examine the characteristic function for the end-to-end probability distribution Weiss (1994) $`c(\omega ,\chi )=\mathrm{exp}(\mathrm{i}𝝎𝐋)_N`$ (where $`𝝎=(\omega \mathrm{cos}\chi ,\omega \mathrm{sin}\chi )^T`$). The projection of the vector $`𝐋=(L_x,L_y)^T`$ onto the axis $`𝝎`$ is given by $`𝝎𝐋=\omega l_k\mathrm{cos}(\stackrel{~}{\theta }_k\chi )`$. Therefore, the function $`q(\stackrel{~}{\theta })`$ (16) reads in this case:

$$q(\stackrel{~}{\theta })=\mathrm{e}^{\mathrm{i}\omega l\mathrm{cos}(\stackrel{~}{\theta }\chi )}=\underset{\nu }{}\mathrm{i}^\nu \mathrm{e}^{\mathrm{i}\nu \chi }J_\nu (\omega l)\mathrm{e}^{\mathrm{i}\nu \stackrel{~}{\theta }},$$
(21)
since $`q(\stackrel{~}{\theta })`$ coincides with the generating function for the Bessel functions of integer index Watson (1944). We take advantage of the scaling symmetry (20) to identify $`\lambda =\mathrm{ie}^{\mathrm{i}\chi }`$ and $`q_\nu =J_\nu (\omega l)`$, and consequently obtain for the characteristic function:

$$c(\omega ,\chi )=\underset{\nu }{}\mathrm{i}^\nu \mathrm{e}^{\mathrm{i}\nu \chi }\left[𝒵(\omega l)^N\right]_{0\nu },$$
(22)
where $`𝒵(\omega l)_{jk}=J_{kj}(\omega l)p_k`$ (19). (In form of a recursive relation, this result is also stated in Ref. Larralde (1997).) Equation (22) may be interpreted as the angular Fourier series for the characteristic function $`c(\omega ,\chi )`$. Remarkably, if $`p_m`$ is the non-vanishing Fourier component of largest index $`|m|`$ in the series for $`p(\theta )`$ (2), then also the angular Fourier expansion of $`c(\omega ,\chi )`$ will be a finite series limited to the range $`|\nu ||m|`$. In this sense, the angular decomposition of the end-to-end distribution directly reflects the properties of the angular bias $`p(\theta )`$ between the individual segments that make up the walk.

Finally, we note that our treatment covers Pearson’s problem Pearson (1905) as the isotropic case $`p(\theta )=1/2\pi `$, or $`p_k=\delta _{0k}`$ (2). Then, the characteristic function (22) is also isotropic and reduces to a simple power of a Bessel function, $`c_P(\omega )=J_0(\omega l)^N`$ Weiss (1994).

### II.3 The probability distribution function

It has been pointed out by various authors Weiss (1994); Flory (1969) that the probability distribution function $`w(R,\varphi )=\delta [𝐑𝐋(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)]_N`$ is connected to the characteristic function $`c(\omega ,\chi )=\mathrm{exp}(\mathrm{i}𝝎𝐋)_N`$ via a Fourier transform in $`𝝎`$–space:

$$w(R,\varphi )=\frac{1}{(2\pi )^2}\mathrm{d}^2𝝎\mathrm{e}^{\mathrm{i}𝝎𝐑}c(\omega ,\chi ).$$
(23)
Since $`𝝎𝐑=\omega R\mathrm{cos}(\chi \varphi )`$, we may expand the exponential in (23) into the series of Bessel functions (21). Inserting the matrix-type expression for $`c(\omega ,\chi )`$ (22) renders the integration over the polar angle $`\chi `$ in (23) trivial, and one obtains an angular Fourier series for the end-to-end distribution function:

$$w(R,\varphi )=\frac{1}{2\pi }\underset{\nu }{}\mathrm{e}^{\mathrm{i}\nu \varphi }_0^{\mathrm{}}\omega d\omega J_\nu (\omega R)\left[𝒵(\omega l)^N\right]_{0\nu }.$$
(24)
We may interpret this result as the generalization of Kluyver’s formula Kluyver (1906) towards the persistent random walk, to which it reduces in the isotropic case $`p_k=\delta _{0k}`$:

$$w_P(R)=\frac{1}{2\pi }_0^{\mathrm{}}\omega d\omega J_0(\omega R)J_0(\omega l)^N.$$
(25)
As for $`c(\omega ,\chi )`$, we again infer that the range of non-vanishing angular Fourier components in $`w(R,\varphi )`$ equals the corresponding range for the local weight function $`p(\theta )`$ concerning the individual angles between steps.

While a rather elegant expression, (24) is nevertheless not well suited towards the accurate determination of $`w(R,\varphi )`$, since this pursuit involves the numerical evaluation of an oscillating integrand (a sum over products of $`N+1`$ individual Bessel functions, see (17)). Indeed, since $`w(R,\varphi )=0`$ for $`R>Nl`$ (a fact that we will take advantage below), the distribution is not an analytic function of $`𝐑`$, and the asymptotics of the integrand in (24) suggests that $`w(R,\varphi )`$ shows singular behavior whenever the lengths $`\pm R`$, $`\pm l`$ involved sum up to zero Flory (1969). (For identical step lengths $`l`$, this occurs at regular intervals $`R=(N2k)l`$, where $`k=0,1,2,\mathrm{},N`$.) From the practical point of view, the eigenfunction expansion of $`w(R,\varphi )`$ that we are going to develop now is much superior. On the other hand, the analytical expressions for the integral (25) known for $`N3`$ provide a useful testing ground for any numerical attempt to evaluate the probability distribution function (Section III).

### II.4 The eigenfunction method

Since the end-to-end distance $`R`$ of a walk of $`N`$ steps, regardless of its details, cannot exceed the contour length $`Nl`$ of the chain, the distribution function $`w(R,\varphi )`$ vanishes outside the disc $`RNl`$. This renders a technique attractive that is familiar from the theory of field equations under boundary conditions, viz., the expansion into eigenfunctions of a self-adjoint operator. For a function confined within a circle, the modes of a drum, i. e. solutions of the equation $`(\mathrm{\Delta }+\lambda _{mk}^2)f_{mk}(R,\varphi )=0`$, provide a natural choice. Properly normalized, these functions read Sommerfeld (1949):

$$f_{mk}(R,\varphi )=\frac{\mathrm{e}^{\mathrm{i}m\varphi }J_m(z_{mk}R/Nl)}{Nl\sqrt{\pi }J_{m+1}(z_{mk})},$$
(26)
where $`J_m(z)`$ is the Bessel function of order $`m`$, $`m`$, $`z_{mk}`$ denotes the $`k`$.th zero of $`J_m(z)`$ ($`k1`$), and $`\lambda _{mk}=z_{mk}/Nl`$ is the assigned eigenvalue. The functions (26) form an orthonormal set that is furthermore complete, i. e., with $`L_x=L\mathrm{cos}\psi `$ we find for all $`R,L<Nl`$:

$$\delta (𝐑𝐋)=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{\mathrm{}}{}}f_{mk}(L,\psi )^{}f_{mk}(R,\varphi ).$$
(27)
Now, the probability distribution function $`w(R,\varphi )=\delta [𝐑𝐋(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)]_N`$ presents just the average of (27), and we thus formally obtain Weiss and Shmueli (1987):

$$w(R,\varphi )=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{\mathrm{}}{}}\frac{F_{mk}\mathrm{e}^{\mathrm{i}m\varphi }J_m\left(z_{mk}R/Nl\right)}{(Nl)^2\pi J_{m+1}(z_{mk})^2}$$
(28)
as the Fourier-Bessel expansion of $`w(R,\varphi )`$, where the weight coefficients $`F_{mk}`$ are given by the averages:

$$F_{mk}=\mathrm{e}^{\mathrm{i}m\psi (\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)}J_m\left(\frac{z_{mk}}{Nl}L(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)\right)_N.$$
(29)
The actual calculation of the coefficients $`F_{mk}`$ follows the general line of argument presented in Section II.2, where only a slight modification is required. The weights $`F_{mk}`$ present averages of functionals of the absolute angles of the segments $`\stackrel{~}{\theta }_\nu `$ in the sense of (12), and we thus seek to establish the angular Fourier series of $`Q_{mk}(\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N)=\mathrm{e}^{\mathrm{i}m\psi }J_m\left(z_{mk}L/Nl\right)`$. Although this may appear as a formidable task at first, an easily proven “projection theorem” (51) for Bessel functions (Appendix A) will provide the desired result:

$$Q_{\mu _1,\mathrm{},\mu _N}^{(mk)}=\delta _{m,\mu _1+\mathrm{}+\mu _N}J_{\mu _1}(z_{mk}/N)\mathrm{}J_{\mu _N}(z_{mk}/N).$$
(30)
Inserting this result into (14) and switching to the indexing scheme $`\{\nu _1,\mathrm{},\nu _N\}`$ (15), we find that $`\nu _1=m`$. Rearranging the remaining summations in the spirit of (17) yields at once:

$$\begin{array}{c}F_{mk}=\underset{\nu _N}{}J_{\nu _N}\left(\frac{z_{mk}}{N}\right)p_{\nu _N}\underset{\nu _{N1}}{}J_{\nu _{N1}\nu _N}\left(\frac{z_{mk}}{N}\right)p_{\nu _{N1}}\hfill \\ \hfill \times \mathrm{}\underset{\nu _2}{}J_{\nu _2\nu _3}\left(\frac{z_{mk}}{N}\right)p_{\nu _2}J_{m\nu _2}\left(\frac{z_{mk}}{N}\right)p_m.\end{array}$$
(31)
By comparison with (18), we infer that this multiple sum can be recast again into the $`N`$.th power of a matrix $`𝒵(z_{mk}/N)`$: $`F_{mk}=\left[𝒵(z_{mk}/N)^N\right]_{0m}`$, whose matrix elements are given by:

$$\left[𝒵(z_{mk}/N)\right]_{\mu \nu }=J_{\nu \mu }\left(\frac{z_{mk}}{N}\right)p_\nu .$$
(32)
In conjunction with (28), this finally yields the Fourier–Bessel expansion of the end-to-end distribution function $`w(R,\varphi )`$ of the persistent random walk:

$$\begin{array}{c}w(R,\varphi )=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{\mathrm{}}{}}\frac{\left[𝒵(z_{mk}/N)^N\right]_{0m}}{(Nl)^2\pi J_{m+1}(z_{mk})^2}\hfill \\ \hfill \times \mathrm{e}^{\mathrm{i}m\varphi }J_m\left(\frac{z_{mk}R}{Nl}\right).\end{array}$$
(33)
Several remarks are in order here. First, we point out that the expansion coefficient $`F_{mk}`$ (31) vanishes if the corresponding internal angular Fourier component $`p_m`$ does. Thus, the angular character of the end-to-end distribution function $`w(R,\varphi )`$ and the angular bias $`p(\theta )`$ at the joints of individual segments (2) are directly related. Like in Section II.2, a terminating Fourier expansion of $`p(\theta )`$ implies a finite-size matrix $`𝒵(z_{mk}/N)`$ in the determination of the coefficients in (33). In particular, for Pearson’s isotropic walk the radial distribution function $`w_P(R)`$ (25) simplifies to:

$$w_P(R)=\underset{k=1}{\overset{\mathrm{}}{}}\frac{J_0(z_{0k}/N)^N}{(Nl)^2\pi J_1(z_{0k})^2}J_0\left(\frac{z_{0k}R}{Nl}\right),$$
(34)
a result first obtained by Barakat Barakat (1974).

Second, series similar in structure to (33) account for arbitrary derivatives of the probability distribution function $`w(R,\varphi )`$. We invoke the following differentiation relation Watson (1944) for the eigenfunctions $`f_{mk}(R,\varphi )`$ (26):

$$\begin{array}{c}\left(\frac{}{x}\pm \mathrm{i}\frac{}{y}\right)^n\left[\mathrm{e}^{\mathrm{i}m\varphi }J_m\left(\frac{z_{mk}R}{Nl}\right)\right]=\hfill \\ \hfill \left(\pm \frac{z_{mk}}{Nl}\right)^n\left[\mathrm{e}^{\mathrm{i}(mn)\varphi }J_{mn}\left(\frac{z_{mk}R}{Nl}\right)\right],\end{array}$$
(35)
(which immediately yields the eigenfunction property $`\mathrm{\Delta }f_{mk}(R,\varphi )=\lambda _{mk}^2f_{mk}(R,\varphi )`$ of these solutions to the wave equation) and apply it to (33):

$$\begin{array}{c}\left(_x+\mathrm{i}_y\right)^r\left(_x\mathrm{i}_y\right)^nw(R,\varphi )=(1)^n\hfill \\ \hfill \times \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{\mathrm{}}{}}\left(\frac{z_{mk}}{Nl}\right)^{n+r}\frac{\left[𝒵(z_{mk}/N)^N\right]_{0m}}{(Nl)^2\pi J_{m+1}(z_{mk})^2}\\ \hfill \times \mathrm{e}^{\mathrm{i}(m+nr)\varphi }J_{m+nr}\left(\frac{z_{mk}R}{Nl}\right).\end{array}$$
(36)
All derivatives of the form $`_x^\mu _y^\nu w(R,\varphi )`$ are then available as properly chosen linear combinations of the Fourier–Bessel series (36), where $`n+r=\mu +\nu `$. The representation (36) is particularly useful in the random walk model of polymers Flory (1969), where the gradient $`\mathbf{}\mathrm{log}w(R,\varphi )`$ is related to the average force exerted by the polymer chain.

Finally, also all physically sensible moments $`x^\mu y^\nu _N`$ ($`\mu ,\nu _0`$) of the probability distribution function in the persistent walk can be represented as fairly simple sums involving the zeroes $`z_{mk}`$ of Bessel functions (58), (59). Since the evaluation of these moments involves an integral over Bessel functions that is not contained in standard tables Gradshteyn and Ryzhik (1980), the details of the calculation have been assembled in Appendix B.

## III Example calculations

In the following, we probe different practical aspects of the Fourier–Bessel expansion (33) derived in the previous section. The numerical convergence of the method is the first topic to be examined in detail, and we present analytical and numerical results obtained on Pearson’s original walk problem below (Section III.1). Some of the distinctive properties of persistent walks are highlighted in Section III.2 where we will study how the underlying angular bias $`p(\theta )`$ influences the shape of the probability distribution function. Our examples include a “chiral” walk problem Larralde (1997) that is biased towards left or right turns, and a model illustrating the transition from diffusive to elastic behavior as the stiffness of the chain is increased.

### III.1 The Pearson random walk

Obviously, particularly simple results are obtained when the Fourier–Bessel method of Section II.4 is applied towards the isotropic Pearson random walk with $`p(\theta )=1/2\pi `$ or $`p_\nu =\delta _{0\nu }`$ (2). The angular sum in (33) is then absent, and the matrix power in the Fourier coefficient reduces to a simple power, which allows convenient access to error estimates for the truncated sums necessarily encountered in a numerical evaluation of $`w_P(R)`$ via the expansion (34).

Unfortunately, the probability distribution function $`w_P(R)`$ for a random walk in two dimensions, unlike its three-dimensional counterpart Flory (1969); Treloar (1946), generally defies evaluation in closed form. However, analytical expressions can be found for the shortest chains with $`N=2`$ and $`N=3`$. A formula for the two-segment chain distribution is available for arbitrary persistent walks (see Section III.2), and reads here $`w_P(R)=[\pi ^2R\sqrt{4l^2R^2}]^1`$ (45). For $`N=3`$, the integral representation (25) can be evaluated in closed form Watson (1944). Denoting the complete elliptic integral of the first kind Gradshteyn and Ryzhik (1980) by $`K(k)`$, we find:

$$\begin{array}{c}w_P(R)=\frac{2}{\pi ^3(l+R)\sqrt{(3lR)(l+R)}}\hfill \\ \hfill \times K\left[\frac{4l}{l+R}\sqrt{\frac{lR}{(3lR)(l+R)}}\right],\end{array}$$
(37)
in the interval $`0R<l`$, and

$$w_P(R)=\frac{1}{2\pi ^3l\sqrt{lR}}K\left[\frac{l+R}{4l}\sqrt{\frac{(3lR)(l+R)}{lR}}\right],$$
(38)
for $`l<R3l`$. Note that $`w_P(R)`$ possesses a logarithmic singularity at $`R=l`$ and a steplike discontinuity at the maximum extension $`R=3l`$, where it approaches the value $`w_P(3l)=1/(4\sqrt{3}\pi ^2l^2)`$. (The division of $`w_P(R)`$ into distinct intervals of length $`2l`$ was already mentioned at the end of Section II.3.) Its exact shape is displayed in Figure 2.

These formulas provide a useful testing ground for the Fourier–Bessel expansion (34). Clearly, this series cannot be convergent in a strict sense since $`w_P(R)`$ diverges as $`Rl`$. For an approximate analysis, we estimate the asymptotic value of the sum terms $`S_\mu (N,R)`$ with large index $`\mu `$. Since $`z_{0\mu }\pi (\mu \frac{1}{4})`$ and $`J_0(\lambda z_{0\mu })\sqrt{2/\pi \lambda z_{0\mu }}\mathrm{cos}(\lambda z_{0\mu }\pi /4)`$ holds as $`\mu \mathrm{}`$ Watson (1944), we find for $`RNl/\mu `$:

$$\begin{array}{c}S_\mu (N,R)\frac{1}{\pi Nl^2}\left(\frac{2N}{\pi ^2\mu }\right)^{\frac{N1}{2}}\sqrt{\frac{l}{R}}\hfill \\ \hfill \times \mathrm{cos}\left[\frac{\pi (4\mu N1)}{4N}\right]^N\mathrm{cos}[\frac{\pi R(4\mu 1)}{4Nl}\frac{\pi }{4}].\end{array}$$
(39)
Thus, the series terms overall decay like $`S_\mu (N,R)\mu ^{(N1)/2}`$. Equation (39) allows us to examine the absolute error $`ϵ(k,R)`$ introduced by truncating the series (34) to terms with index $`\mu k`$. This quantity is the sum of the neglected terms, $`ϵ(k,R)=_{\mu >k}S_\mu (N,R)`$, and its asymptotics depend on the behavior of the oscillating terms in the second line of (39). If $`NR/l`$ is not en even integer, the average of the oscillating term vanishes, the series is of the alternating type, and for an estimate of $`ϵ(k,R)`$ we may multiply $`S_k(N,R)`$ with the largest period present in the product of cosine functions. These periods turn out to be proportional to the inverse distance of $`R`$ to the interval boundaries $`(N2\nu )l`$ (see Section II.3). Indeed, for $`N=3`$ the maximum error $`ϵ(k,R)`$ incurred in the truncation of the series is approximated well by the empiric formula:

$$ϵ(k,R)\frac{0.0076}{kl\sqrt{lR}}\left(\frac{1}{2}+\frac{3l}{|lR|}+\frac{2l}{3lR}\right).$$
(40)
In the case that $`NlR`$ is an even multiple of $`l`$, the Fourier series contains a strictly positive (non-oscillating) part that dominates the remainder. Its summation yields the estimate $`ϵ(k,R)k^{(N3)/2}`$ for the interval boundaries, except for $`N=3`$, where it reproduces the logarithmic singularity at $`R=l`$ (37). Examples for the numerical convergence of the Bessel–Fourier series for $`N=3`$ are shown in Figure 2 and Figure 3, where the relative error $`ϵ_{\mathrm{rel}}(k,R)=ϵ(k,R)/w_P(R)`$ is displayed together with the estimate (40) for various values of $`k`$.

Clearly, for fixed $`N`$ the absolute convergence of the expansion (34) is slowest at the origin $`R=0`$. We may then use arguments similar to (39) to establish an absolute error bound $`ϵ(k)`$ for the Fourier–Bessel series, which then yields $`k`$ as a function of the desired accuracy $`ϵ`$:

$$k(ϵ)\frac{4N}{\pi ^2}\left[\frac{16\mathrm{\Gamma }\left(\frac{N+2}{4}\right)^2}{\pi ^4(N4)\mathrm{\Gamma }\left(\frac{N+2}{2}\right)ϵ}\right]^{\frac{2}{N4}},$$
(41)
Note that for short chains, the absolute series convergence at $`R=0`$ ($`kϵ^{2/(N4)}`$) is much worse than its typical convergence ($`kϵ^{2/(N1)}`$). In particular, (41) predicts the divergence of $`w_P(R)`$ for $`N=4`$ as $`R`$ approaches the origin. Invoking the general relation $`w_P(N+1,R=0)=w_P(N,R=l)`$ that follows at once from (25), the analytic expressions (37), (38) immediately confirm this assertion. With increasing number of segments $`N`$, the number of sum terms $`k(ϵ)`$ drops rapidly, and finally becomes of the order of the chain length itself, in which case the estimate (41) becomes too conservative. In the long-chain regime $`\pi k(ϵ)N`$ the argument of the Bessel function in (34) tends towards zero, and the approximation $`J_0(z_{0k}/N)^N(1z_{0k}^2/4N^2)^N\mathrm{exp}(z_{0k}^2/4N)`$ instantly yields the asymptotic behavior:

$$k(ϵ)\frac{2}{\pi }\sqrt{N\mathrm{log}(\pi Nϵ)}.$$
(42)
Therefore, for long chains, the numerical complexity of the Fourier–Bessel series approach grows merely with $`(N\mathrm{log}N)^{1/2}`$. Figure 4 depicts numerical results for $`k(ϵ)`$ for various values of the absolute accuracy $`ϵ`$, together with the estimates (41) and (42). We conclude that these expressions closely reproduce $`k(ϵ)`$ within their respective regimes of validity.

Following this excursion on the numerical aspects of the Fourier-Bessel approach, we next study the resulting probability distributions $`w_P(R)`$ (34) for the Pearson random walk. As a random process, $`w_P(R)`$ must be normally distributed in the limit $`N\mathrm{}`$, and thus is completely determined by the diffusion constant $`D_P=l^2`$ (10). After proper normalization, $`w_P(R)`$ asymptotically behaves like:

$$w_P(R)\frac{1}{\pi Nl^2}\mathrm{exp}\left(\frac{R^2}{Nl^2}\right),$$
(43)
assuming that $`RNl`$. For short walks of less than $`N=10`$ segments, deviations from the Gaussian limit become clearly apparent, and the separation of $`w_P(R)`$ into distinct functions for every interval of length $`2l`$ reveals itself. The nonanalytic behavior of $`w_P(R)`$ at the interval boundaries $`R=(N2\nu )l`$ is manifested by the presence of divergences ($`N=3,4`$) and kinks ($`N=4,5,6`$) in the distribution functions, as shown in Figure 5.

Finally, we examine the force function $`F_N(R)`$ assigned to the random walk problem. Persistent random walks have been proposed as a model for the configuration of polymer chains Flory (1969). In this statistical interpretation, the bias function $`p(\theta )`$ (2) is related to the angular potential between adjacent segments $`U(\theta )`$ via (properly normalized) Boltzmann factors $`p(\theta )\mathrm{exp}[U(\theta )/k_BT]`$ which renders the Pearson walk equivalent to the uncorrelated freely joined chain (FJC) model of the polymer. For a fixed end position $`𝐑`$ of the polymer chain, the molecule exerts an entropic, temperature-dependent force whose average is given by the logarithmic derivative of the distribution function Flory (1969): $`𝐅(𝐑)=k_BT\mathbf{}\mathrm{log}w(𝐑)`$, and its absolute value reduces to $`F_N(R)=k_BTw_P(R)^{}/w_P(R)`$ for the isotropic walk. Like the distribution function itself, its derivatives are conveniently calculated using the Fourier–Bessel formalism (36), and the results of this effort are displayed in Figure 6 for the two-dimensional FJC model.

For short chains, the singularities of $`w_P(R)`$ at the interval boundaries $`R=(N2\nu )l`$ confer a conspicious irregular shape to the corresponding force curves $`F_N(R)`$. As $`N`$ increases, $`F_N(R)`$ becomes smoother, and as a function of the scaled end-to-end distance $`R/Nl`$, approaches a well-defined limit $`F_{\mathrm{}}(R)`$. Its analytic form may be obtained by considering the FJC in the Gibbs ensemble, where the force $`𝐅`$ is kept constant, and $`𝐑`$ is allowed to vary Flory (1969). Its partition function is given by $`Z(𝐅)=\mathrm{exp}(𝐋𝐅/k_BT)`$, and the average position of the chain end $`𝐑`$ under the applied force $`𝐅`$ follows at once from $`𝐑(𝐅)_N=k_BT\mathbf{}\mathrm{log}Z(𝐅)`$. Using the technique presented in Section II.2, the function $`Z(𝐅)`$ can be cast as the trace over a matrix power (18) for an arbitrary persistent walk. Indeed, $`Z(𝐅)`$ may be interpreted as the characteristic function $`c(𝝎)`$ for imaginary argument $`𝝎=\mathrm{i}𝐅/k_BT`$ (22), and thus reduces for the isotropic walk to $`Z_P(F)=I_0(lF/k_BT)^N`$, where $`I_0(z)`$ denotes the modified Bessel function Watson (1944). Hence, for the FJC model we find $`R(F)_N=NlI_1(u)/I_0(u)`$, where $`u=lF/k_BT`$. In the limit $`N\mathrm{}`$, we may identify the average position $`R(F)_N`$ with the actual position $`R`$, and likewise $`F_N(R)`$ with $`F`$. This allows to revert the relation to express $`F_{\mathrm{}}(R)`$ as a function of $`R/Nl`$:

$$F_{\mathrm{}}(R)=\frac{k_BT}{l}𝒰(R/Nl),$$
(44)
where $`𝒰(z)`$ is the inverse function to $`I_1(z)/I_0(z)`$. It is defined for $`|z|<1`$ and behaves like $`𝒰(z)=2z+z^3+\frac{5}{6}z^5+\mathrm{}`$ for small $`|z|`$, while its Laurent series near $`z=1`$ is given by $`𝒰(1z)=\frac{1}{2z}+\frac{1}{4}+\frac{3}{8}z+\mathrm{}`$

### III.2 Properties of persistent walks

Analytical methods towards the calculation of the probability distribution function $`w(R,\varphi )`$ for the more general persistent walk problem are usually limited to the combination of only two segments ($`N=2`$). In this case, for a fixed endpoint $`𝐋=(R\mathrm{cos}\varphi ,R\mathrm{sin}\varphi )`$ only two possible locations for the intermediate joint exist, and simple geometrical arguments lead to the following result:

$$\begin{array}{c}w(R,\varphi )=\frac{2}{R\sqrt{4l^2R^2}}\hfill \\ \hfill \times \left[p(2\gamma )p(\varphi \gamma )+p(2\gamma )p(\varphi +\gamma )\right],\end{array}$$
(45)
where we used $`\gamma =\mathrm{arccos}(R/2l)`$ to denote the angle between the end-to-end vector and the individual segments. (For the uncorrelated Pearson walk, the term in brackets reduces to $`1/(2\pi ^2)`$.) Likewise, the limit of very extended walks ($`N\mathrm{}`$) is amenable to an analytical treatment. Here, the end-to-end probability function $`w(R,\varphi )`$ must be again normally distributed, and is completely determined by the persistence lengths $`L_x_{\mathrm{}}`$ (5), $`L_y_{\mathrm{}}`$ (6), and the diffusion constant $`D_{\mathrm{}}`$ (10):

$$w(𝐑)\frac{1}{\pi ND_{\mathrm{}}}\mathrm{exp}\left[\frac{(𝐑𝐋_{\mathrm{}})^2}{ND_{\mathrm{}}}\right].$$
(46)
We point out that the distribution $`w(R,\varphi )`$ for sufficiently large $`N`$ depends merely on the single angular Fourier component $`p_1`$. Little can be said a priori about the behavior of $`w(R,\varphi )`$ for intermediate chain lengths $`N`$. Fortunately, the Fourier–Bessel expansion (33) provides an efficient means to determine the probability distribution function numerically.

In the following, we select normalized angular bias functions $`p(\theta )`$ of the general form:

$$p(\theta )=\frac{2^{2\xi 1}\mathrm{\Gamma }(\xi +1)^2}{\pi \mathrm{\Gamma }(2\xi +1)}\mathrm{cos}^{2\xi }\left(\frac{\theta \beta }{2}\right),$$
(47)
where $`\xi 0`$ and $`\beta `$ are real parameters. The corresponding angular Fourier components $`p_\nu `$ (2) are given by Gradshteyn and Ryzhik (1980):

$$p_\nu (\xi ,\beta )=\frac{\mathrm{\Gamma }(\xi +1)^2\mathrm{e}^{\mathrm{i}\nu \beta }}{\mathrm{\Gamma }(\xi +\nu +1)\mathrm{\Gamma }(\xi \nu +1)}.$$
(48)
In practice, the coefficients $`p_\nu (\xi ,\beta )`$ are best calculated by the simple recursion formula $`p_{\nu +1}=\mathrm{e}^{\mathrm{i}\beta }(\xi \nu )p_\nu /(\xi +\nu +1)`$ ($`\nu 0`$), where normalization requires $`p_0=1`$. We find in particular $`p_1=\mathrm{e}^{\mathrm{i}\beta }\xi /(\xi +1)`$ and thus:

$`𝐋_{\mathrm{}}`$ $`={\displaystyle \frac{l\xi }{1+4\xi (\xi +1)\mathrm{sin}^2\frac{\beta }{2}}}\left(\begin{array}{c}12(\xi +1)\mathrm{sin}^2\frac{\beta }{2}\\ (\xi +1)\mathrm{sin}\beta \end{array}\right),`$ (49)
$`D_{\mathrm{}}`$ $`={\displaystyle \frac{l^2(2\xi +1)}{1+4\xi (\xi +1)\mathrm{sin}^2\frac{\beta }{2}}},`$ (50)
in (46). Note that the angle $`\beta `$ in (47) controls the chiral bias of the walk Larralde (1997); a symmetrical distribution will ensue if $`\beta =0`$ is chosen. (Then, (49) simplifies to $`L_x_{\mathrm{}}=l\xi `$, $`L_y_{\mathrm{}}=0`$, and $`D_{\mathrm{}}=l^2(2\xi +1)`$.) Finally, we remark that the series in (47) terminates whenever $`\xi `$ is an integer. (Otherwise, $`p(\theta )`$ is not analytic at its zero $`\theta _0=\pi +\beta `$.)

We start out with a study of the dependence of $`w(R,\varphi )`$ on the number of steps $`N`$ for the fixed angular weight $`p(\theta )=\frac{4}{3\pi }\mathrm{cos}^4[\frac{1}{2}(\theta \frac{\pi }{6})]`$, which is of the form (47) with $`\xi =2`$ and $`\beta =30^{}`$. The angular bias evidently favors a “spiral” configuration of the chain in which adjacent segments are joined at an angle of $`\beta `$, which translates to an absolute angle $`\varphi =(N+1)\beta /2`$ of the end-to-end vector with respect to the $`x`$–axis. However, configurations that end up near the origin constitute the vast majority in the phase space $`\{\theta _1,\mathrm{},\theta _N\}`$ of the walk. For sufficiently large $`N`$, the latter contribution always dominates and ultimately yields the Gaussian distribution (46), while the shape of $`w(R,\varphi )`$ is governed by the intrinsic stiffness of the chain for short walks.

In this regime, large deviations from the ideal intersegment angle $`\beta `$ are suppressed by the angular bias function $`p(\theta )`$, so walks that significantly contribute to $`w(R,\varphi )`$ show small variations of their bond angles $`\theta _k`$ around $`\beta `$. While these variations may significantly alter the total turning angle $`\varphi `$ of the persistent walk, they usually have little effect on the end-to-end distance $`R`$, at least when the preferred configuration is mostly stretched (i. e., $`N\beta `$ is small). The situation then resembles the elastic response of a thin rod (or, in the case of chiral bias $`\beta 0`$, a circular spring), with one end fixed at the origin: While easily yielding to bending forces, the rod (spring) is resilient to compression, and resists changes of its length. Under applied force, the condition of constant curvature thus requires that the rod bends into a circular arc of fixed length $`\mathrm{\Lambda }`$, and a simple geometrical argument shows that the free end of the rod moves on the curve $`R=\mathrm{\Lambda }\mathrm{sin}\varphi /\varphi `$. The shape of this curve is displayed in Figure 7.

A numerical calculation of the probability distribution function $`w(R,\varphi )`$ for the aforementioned angular bias $`p(\theta )`$ confirms these deliberations. In Figure 8, we present $`w(R,\varphi )`$ for walks of length $`N=4`$, $`N=7`$, and $`N=10`$. For the short walk, $`w(R,\varphi )`$ is concentrated along a narrow groove whose shape is reminiscent of the geometrical curve presented in Figure 7, while for the extended chain, the distribution is largely of Gaussian shape, albeit shifted from the origin, as predicted in (49). In the transitory case $`N=7`$, the end-to-end distribution function conspicuously combines the typical features of both limiting cases, as both the central hump and a protruding ridge are clearly visible. (Note that the isolated spike in the center of the left image is not a numerical artifact; the actual distribution for $`N=4`$, like the corresponding uncorrelated Pearson walk (Figure 5) contains a logarithmic singularity at the origin. In the calculation on display, we chose to truncate the Fourier–Bessel series (33) after $`k=10^4`$ terms.)

In a second simulation, we investigate the probability distribution function $`w(R,\varphi )`$ in a persistent walk of $`N=20`$ segments as a function of the stiffness of the joints, i. e., the parameter $`\xi `$ in the angular bias function $`p(\theta )`$ (47). In this example, we chose a symmetric distribution ($`\beta =0`$) and calculated the average entropic force $`F(𝐑)=k_BT|\mathbf{}\mathrm{log}w(𝐑)|`$ (see Section III.1) for the three stiffness parameters $`\xi =2`$, $`\xi =4`$, and $`\xi =6`$. (Under these circumstances, convergence of the Fourier–Bessel series (33) and (36) was achieved after summing $`k=250`$ terms.) Lines of constant force are plotted in Figure 9. For the softest chain ($`\xi =2`$), the pattern resembles a set of equidistant, concentric rings that has been shifted from the origin. This observation implies that the corresponding distribution $`w(R,\varphi )`$ has approached the long-chain limit (46) characterized by an isotropic harmonic force field $`F(𝐑)`$ directed towards a minimum $`𝐑_{\mathrm{min}}`$ of the free energy $`U(𝐑)=k_BT\mathrm{log}w(𝐑)`$. Its position is given by the persistence length vector $`𝐑_{\mathrm{min}}=𝐋_{\mathrm{}}=(l\xi ,0)^T`$ (49), and the elastic constants $`ϵ_{jk}=_j_kU(𝐑_{\mathrm{min}})`$ there read $`ϵ_{xx}=ϵ_{yy}=2k_BT/ND_{\mathrm{}}=2k_BT/[(2\xi +1)Nl^2]`$, and $`ϵ_{xy}=0`$. The numerically determined values $`x_{\mathrm{min}}=2.3224l`$, $`ϵ_{xx}=1.52210^2k_BT/l^2`$, and $`ϵ_{yy}=1.43810^2k_BT/l^2`$ still deviate considerably from their counterparts in the diffusive limit, $`x_{\mathrm{min}}=2l`$ and $`ϵ_{xx}=ϵ_{yy}=0.02k_BT/l^2`$, however.

The force distributions obtained for the stiffer walks with $`\xi =4`$ and $`\xi =6`$ bear little resemblance with the diffusive circular pattern observed for $`\xi =2`$. For $`\xi =4`$ (center panel in Figure 9), the force field is distinctively flat over a wide range of end-to-end vectors $`𝐑`$, and the free energy $`U(𝐑)`$ actually features two local minima $`x_<=3.9808l`$, $`x_>=15.160l`$ on the $`x`$–axis. The latter represents the maximum of the probability density with $`w(x_>,0)1.9110^3l^2`$, but their difference in free energy is only $`\mathrm{\Delta }U=0.2336k_BT`$. The minima are separated by a shallow maximum near $`x_{\mathrm{max}}=9.15l`$ ($`\mathrm{\Delta }U=0.2846k_BT`$). It is instructive to compare the elastic constants $`ϵ_{jk}`$ at both minima. For $`x_<`$, we find an almost isotropic response ($`ϵ_{xx}=7.42110^3k_BT/l^2`$, $`ϵ_{yy}=7.00910^3k_BT/l^2`$), while the free energy surface near $`x_>`$ is strongly prolate ($`ϵ_{xx}=9.52910^2k_BT/l^2`$, $`ϵ_{yy}=2.55710^3k_BT/l^2`$). Clearly, this anisotropy indicates the onset of elastic behavior.

For $`\xi =6`$ (right panel), elastic response becomes the dominant feature of the distribution which is predominantly concentrated along a narrow arc, while a rather flat plateau in its center still persists. The maximum of $`w(R,\varphi )`$ has shifted further out to $`x_>=16.790l`$. (An extremely shallow local minimum at $`x_<=2.4924l`$ ($`\mathrm{\Delta }U=1.222k_BT`$) with a barrier of only $`0.009k_BT`$ is contained within the plateau structure.) The anisotropy of the free energy surface near $`x_>`$ has become even more marked, where the bending mode ($`ϵ_{yy}=3.23210^3k_BT/l^2`$) is almost 75 times softer than the response towards compression ($`ϵ_{xx}=2.40510^1k_BT/l^2`$). Empirically, the shape of the trough structure in $`U(𝐑)`$ is accurately reproduced by the elastic spring model (with $`\mathrm{\Lambda }=Nl`$) put forward in Figure 7, if we allow for a displacement of its origin so that the curve $`R(\varphi )=Nl\mathrm{sin}\varphi /\varphi `$ passes through $`x_>`$ (dashed line in Figure 9). The elastic compressive response is limited to small deviations from the equilibrium point $`x_>`$: For applied loads greater than a critical force $`F_{\mathrm{crit}}=0.195k_BT/l`$ (at $`x_{\mathrm{crit}}=13.90l`$), the stretched chain will buckle and collapse into the more globular shapes making up the plateau. This again resembles the behavior of an elastic rod under stress Feynman et al. (1964).

## IV Discussion

In this article, we examined persistent random walks in two dimensions, where the relative orientation of successive segments is correlated through an angular bias function $`p(\theta )`$. We found that the eigenfunction expansion into a Fourier–Bessel series provides a convenient and efficient means to calculate the end-to-end probability distribution function $`w(𝐑)`$ encountered in a planar random walk except for extremely short chains. Related series expansions also cover all derivatives and moments of $`w(𝐑)`$. We used Pearson’s original random walk problem to assess the numerical aspects of the Fourier–Bessel series and found that the number of sum terms required to achieve a desired accuracy initially sharply decreases with the chain length $`N`$ and, for very extended walks, grows like $`(N\mathrm{log}N)^{1/2}`$. This property renders the method attractive in studying the properties of chains of intermediate length whose distribution function $`w(𝐑)`$ may strongly deviate from the diffusive behavior expected in the limit $`N\mathrm{}`$. Even for uncorrelated planar walks, the singularities contained in $`w_P(R)`$ cause conspicuous features in the probability distribution and force curves that only gradually fade away as $`N`$ increases. Naturally, correlated walks generate a much richer set of phenomena that we only have begun to explore in this work. Here, we briefly studied “chiral” walks with asymmetric angular bias $`p(\theta )`$ that have a tendency to loop, and inquired into the dependence of $`w(𝐑)`$ on the stiffness of the intersegment joints. Clear evidence of a transition from diffusive character to a form of elastic response familiar from the elastomechanics of a stiff rod, was obtained.

The latter observation raises questions that deserve further examination. It appears worthwhile to study the dependence of the crossover from diffusive to elastic behavior as a function of $`N`$ and the bias $`p(\theta )`$. Characteristic parameters of the distribution, like the location of the minima $`𝐑_{\mathrm{min}}`$ of the free energy, the elastic moduli $`ϵ_{jk}`$ and the critical force $`F_{\mathrm{crit}}`$, may obey scaling laws. The existence of several minima of $`U(𝐑)`$ suggests that the average extension of a stiff chain under an applied external force shows hysteresis. The formalism developed in this article directly applies to these problems.

Clearly, for many physical applications, e. g., in polymer theory Flory (1969), realistic models will require the random walk to take place in three spatial dimensions. While the analysis of the uncorrelated random walk problem encounters no major obstacles Barakat (1973); Treloar (1946), the extension of the eigenfunction expansion technique to persistent walks in three-dimensional space presents a challenging problem, since the angular bias $`p(\theta ,\varphi )`$ then generally will describe the relative orientation of a segment and its two predecessors, unless $`p(\theta ,\varphi )`$ is independent of the dihedral angle $`\varphi `$. Alternative methods, like the transfer operator formalism Livadaru et al. (2002), may prove to be better suited for the numerical determination of the probability density function $`w(𝐑)`$ in this case.

###### Acknowledgements.

The author appreciates useful discussions with H.–J. Kreuzer. This work was financially supported by the Alexander von Humboldt foundation and the Killam trust.
* * *

## Appendix A The projection theorem

This appendix provides a brief mathematical derivation of the projection theorem for Bessel functions which forms the basis of the Fourier–Bessel expansions derived in Section II.4. If $`(R,\varphi )`$ denotes the location of the chain end in polar coordinates for a walk consisting of $`N`$ steps of individual lengths $`\{l_1,\mathrm{},l_N\}`$ characterized by absolute angles $`\{\stackrel{~}{\theta }_1,\mathrm{},\stackrel{~}{\theta }_N\}`$ (see Figure 1), this proposition states that:

$$\begin{array}{c}\mathrm{e}^{\mathrm{i}\mu \varphi }J_\mu \left(z\frac{R}{\mathrm{\Lambda }}\right)=\underset{\lambda _1+\mathrm{}+\lambda _N=\mu }{}\mathrm{e}^{\mathrm{i}(\lambda _1\stackrel{~}{\theta }_1+\mathrm{}+\lambda _N\stackrel{~}{\theta }_N)}\hfill \\ \hfill \times J_{\lambda _1}\left(zl_1/\mathrm{\Lambda }\right)\mathrm{}J_{\lambda _N}\left(zl_N/\mathrm{\Lambda }\right),\end{array}$$
(51)
where $`\mathrm{\Lambda }=l_1+\mathrm{}+l_N`$ represents the contour length of the walk, $`\mu `$ is an integer, and $`z`$ denotes an arbitrary complex number. (For the case of equal step lengths $`l_\nu =l`$ primarily considered in this article, the arguments of the Bessel functions simplify to $`zR/Nl`$ on the l.h.s. and $`z/N`$ throughout the r.h.s. of this equation, respectively.)

We now sketch the proof of (51). From Figure 1, it is clear that the projection $`\mathrm{\Lambda }_x`$ of the end-to-end vector on the $`x`$–axis is simply the sum of the correponding projections of the steps that make up the walk: $`R\mathrm{cos}\varphi =_\nu l_\nu \mathrm{cos}\stackrel{~}{\theta }_\nu `$ ($`\nu =1,\mathrm{},N`$). It is convenient to perform this projection also in an alternative reference frame which is rotated by the angle $`\alpha `$, i. e. with respect to an arbitrary direction. Since this operation does not affect the lengths involved and shifts the angles $`\varphi `$ and $`\stackrel{~}{\theta }_\nu `$ merely by the constant value $`\alpha `$, we find the relation $`R\mathrm{cos}(\varphi \alpha )=_\nu l_\nu \mathrm{cos}(\stackrel{~}{\theta }_\nu \alpha )`$, or in exponential form:

$$\mathrm{e}^{\mathrm{i}uR\mathrm{cos}(\varphi \alpha )}=\underset{\nu =1}{\overset{N}{}}\mathrm{e}^{\mathrm{i}ul_\nu \mathrm{cos}(\stackrel{~}{\theta }_\nu \alpha )},$$
(52)
for all $`u`$ and $`\alpha `$. We employ the generating function expansion for the Bessel functions $`J_\nu (u)`$, $`\mathrm{e}^{\mathrm{i}u\mathrm{cos}\alpha }=_k\mathrm{i}^k\mathrm{e}^{\mathrm{i}k\alpha }J_k(u)`$ (21) on both sides of (52) with $`u=z/\mathrm{\Lambda }`$ to obtain the equality:

$$\begin{array}{c}\underset{\lambda }{}\left(\mathrm{i}\mathrm{e}^{\mathrm{i}\alpha }\right)^\lambda \mathrm{e}^{\mathrm{i}\lambda \varphi }J_\lambda \left(z\frac{R}{\mathrm{\Lambda }}\right)=\underset{\lambda _1,\mathrm{},\lambda _N}{}\left(\mathrm{i}\mathrm{e}^{\mathrm{i}\alpha }\right)^{\lambda _1+\mathrm{}+\lambda _N}\hfill \\ \hfill \times \mathrm{e}^{\mathrm{i}(\lambda _1\stackrel{~}{\theta }_1+\mathrm{}+\lambda _N\stackrel{~}{\theta }_N)}J_{\lambda _1}(zl_1/\mathrm{\Lambda })\mathrm{}J_{\lambda _N}(zl_N/\mathrm{\Lambda }).\end{array}$$
(53)
The Fourier decomposition of (53) with respect to the angle $`\alpha `$, i. e. an integration $`_0^{2\pi }d\alpha \mathrm{e}^{\mathrm{i}\mu \alpha }\mathrm{}`$, then immediately yields the projection theorem (51).

## Appendix B Evaluation of moments

The Fourier–Bessel series (33) for $`w(R,\varphi )`$ permits fast and accurate numerical evaluation of all moments of the probability distribution in the persistent planar walk. For this two-dimensional problem, the moments are averages of the general form $`L_x^jL_y^k_N`$. Since it is easily verified that they may be expressed as a finite combination of the complex averages $`R^{m+2n}\mathrm{e}^{\mathrm{i}m\varphi }_N`$, where $`j+k=m+2n`$ and $`m,n_0`$, we prefer to use the radial definition, as it leads to shorter expressions. Examples are encountered in Section II.1, where we calculated the average position $`L_x+\mathrm{i}L_y_N=R\mathrm{e}^{\mathrm{i}\varphi }_N`$ (11), corresponding to $`m=1`$ and $`n=0`$, and the mean square displacement $`L_x^2+L_y^2_N=R^2_N`$ (17), where $`m=0`$ and $`n=1`$. (We also note in advance that for $`m=n=0`$, conservation of probability demands $`R^0_N=1`$.)

The determination of the radial moments involves an integral over Bessel functions that apparently is not listed in the major tables Gradshteyn and Ryzhik (1980). Hence, we present a short sketch of its derivation. We start out with a recursion formula stated by Watson Watson (1944). Introducing the notation:

$$I_{mk}^{(n)}=_0^{z_{mk}}duu^{m+2n+1}J_m(u),$$
(54)
where $`m,n0`$ are integers, and $`z_{mk}`$ denotes the $`k`$.th zero of the Bessel function $`J_m(u)`$, this recurrence reads:

$$I_{mk}^{(n)}=z_{mk}^{m+2n+1}J_{m+1}(z_{mk})4n(m+n)I_{mk}^{(n2)}.$$
(55)
Since the recursion terminates for $`n=0`$, we may add up the terms and rearrange them as the product of $`J_{m+1}(z_{mk})`$ with a polynomial in the zero $`z_{mk}`$:

$$I_{mk}^{(n)}=J_{m+1}(z_{mk})\underset{\nu =0}{\overset{n}{}}\frac{(4)^{n\nu }n!(m+n)!}{\nu !(m+\nu )!}z_{mk}^{m+2\nu +1}.$$
(56)
From the series expansion (33), we now calculate the radial moments $`R^{m+2n}\mathrm{e}^{\mathrm{i}m\varphi }_N`$. Since the angular integration is trivial, using the abbreviation (54) in the first step we obtain:

$`\mathrm{}_N`$ $`={\displaystyle _0^{Nl}}dRR^{m+2n+1}{\displaystyle _0^{2\pi }}d\varphi \mathrm{e}^{\mathrm{i}m\varphi }w(R,\varphi )`$ (57)
$`=2(Nl)^{m+2n}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left[𝒵(z_{mk}/N)^N\right]_{0m}I_{mk}^{(n)}}{z_{mk}^{m+2n+2}J_{m+1}(z_{mk})^2}},`$
where the elements of the matrix $`𝒵`$ are given by $`𝒵_{jl}=J_{lj}(z_{mk}/N)p_l`$ (see Section II.4). Using (56), we rearrange this sum into a series over powers of $`z_{mk}`$. It is convenient to define a set of auxiliary sums $`B_{mr}^{(N)}`$ involving the zeroes of the Bessel functions $`J_m(u)`$:

$$B_{mr}^{(N)}=\underset{k=1}{\overset{\mathrm{}}{}}\frac{\left[𝒵(z_{mk}/N)^N\right]_{0m}}{z_{mk}^{2r+1}J_{m+1}(z_{mk})}.$$
(58)
The final result for the radial moments then can be cast in the following form:

$$\begin{array}{c}R^{m+2n}\mathrm{e}^{\mathrm{i}m\varphi }_N=2^{2n+1}(Nl)^{m+2n}\hfill \\ \hfill \times \underset{\nu =0}{\overset{n}{}}\frac{(1)^{n\nu }n!(m+n)!}{2^{2\nu }\nu !(m+\nu )!}B_{m,n\nu }^{(N)}.\end{array}$$
(59)
For not too small $`N`$, the series $`B_{mr}^{(N)}`$ (58) converge rapidly, thus permitting efficient numerical evaluation of the moments using (59). (For practical purposes, we note that the average $`R^{m+2n}\mathrm{e}^{\mathrm{i}m\varphi }_N`$ (59) is a function of the angular Fourier coefficients $`p_0,p_{\pm 1},\mathrm{},p_{\pm (m+n)}`$ only, as the method of direct evaluation exposed in Section II.1 shows. In the calculation of $`B_{mr}^{(N)}`$ (58), it is thus convenient to set $`p_\mu =0`$ for $`|\mu |>m+r`$.)

Finally, we present several examples to illustrate our result for the moments of $`w(R,\varphi )`$. For $`n=0`$, the series (59) reduces to a single term: $`(L_x+\mathrm{i}L_y)^m_N=2(Nl)^mB_{m0}^{(N)}`$. In the case $`m=0`$, conservation of probability then requires $`B_{00}^{(N)}=1/2`$, irrespective of $`N`$ and $`p(\theta )`$. Similarly, comparison with (3) yields for $`m=1`$:

$$B_{10}^{(N)}=\frac{p_1(1p_1^N)}{2N(1p_1)}.$$
(60)
We may use these relations to keep track of the accuracy in the summation of the Fourier-Bessel series (33). Similarly, from (59) we obtain all purely radial moments $`R^{2n}_N`$ of even order. They read for $`n3`$:

$`R^2_N=`$ $`(Nl)^2\left(18B_{01}^{(N)}\right),`$ (61)
$`R^4_N=`$ $`(Nl)^4\left(132B_{01}^{(N)}+128B_{02}^{(N)}\right),`$ (62)
$`R^6_N=`$ $`(Nl)^6\left(172B_{01}^{(N)}+1152B_{02}^{(N)}4608B_{03}^{(N)}\right),`$ (63)
where we already inserted $`B_{00}^{(N)}=1/2`$. Clearly, an expression for $`B_{01}^{(N)}`$ in terms of the angular bias Fourier components $`p_{\pm 1}`$ akin to (60) is available from our earlier result for the mean square displacement (10). From there, we also infer that $`R^2_N`$ in the limit $`N\mathrm{}`$ only grows linearly with $`N`$ (18). Hence, asymptotically the Bessel sum coefficient $`B_{01}^N`$ in (61) approaches the value $`B_{01}^{\mathrm{}}=\frac{1}{8}`$. Similar considerations for the higher moments in conjunction with (58) then show that:

$$B_{0r}^{(\mathrm{})}=\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{z_{0k}^{2r+1}J_1(z_{0k})}=\frac{a_r}{2^{2r+1}(r!)^2},$$
(64)
where the coefficients $`a_r`$ form a sequence of integers: $`a_0=1`$, $`a_1=1`$, $`a_2=3`$, $`a_3=19`$, etc.


format: text

text:
# INTEGRAL and XMM-Newton observations of the weak GRB 030227 1footnote 11footnote 1Based on observations with INTEGRAL, an ESA project with instruments and science data centre funded by ESA member states (especially the PI countries: Denmark, France, Germany, Italy, Switzerland, Spain), Czech Republic and Poland, and with the participation of Russia and the USA, and XMM-Newton an ESA science mission with instruments and contributions directly funded by ESA member states and the USA.

## 1 Introduction

Our understanding of Gamma Ray Bursts (GRBs), one of the mysteries of high-energy astrophysics for more than 25 years, advanced dramatically after the discovery of their X–ray, optical and radio afterglows (see, e.g., van Paradijs, Kouveliotou & Wijers 2000). Currently, the multi-wavelengths study of GRBs is providing a wealth of results relevant for several branches of astrophysics. Besides the study of the prompt emission, a rapid and accurate localization is crucial to pursue these objectives and fully exploit the information on the GRB afterglows and host environment. Such a capability was first achieved by the BeppoSAX satellite, which located more than 45 GRBs from 1997 to 2002 (see, e.g., Amati et al. 2002), and is currently provided by HETE II, a satellite specifically devoted to this task (Ricker et al. 2002).

The INTEGRAL satellite (INTErnational Gamma Ray Astrophysical Laboratory, Winkler et al. 1999) was successfully launched on October 17, 2002. The INTEGRAL payload consists of two $`\gamma `$-ray instruments operating in the $``$15 keV - 10 MeV range: IBIS (Ubertini et al. 1999) and SPI (Vedrenne et al. 1999). Both are coded mask telescopes222With the coded mask technique it is possible to obtain images at energies where photons cannot be easily reflected by interposing a partially absorbing mask between the source and a position sensitive detector (see, e.g., Dean 1983). optimized for angular and spectral resolution respectively.

Although not specifically devoted to GRB studies, thanks to the large field of view and good imaging capabilities of its $`\gamma `$-ray instruments, INTEGRAL is able to localize GRBs at an expected rate of 1-2 per month (Mereghetti, Cremonesi & Borkowski 2001). The data are immediately transmitted to the ground and it is possible to derive the coordinates of detected GRBs with very small delays. Automatic software, the INTEGRAL Burst Alert System (IBAS, Mereghetti et al. 2001), has been developed at the INTEGRAL Science Data Center (ISDC, Courvoisier et al. 1999) to exploit these capabilities. As soon as INTEGRAL telemetry reaches the ISDC, the IBAS software screens the data in real time looking for the presence of potential GRBs, performs a rapid imaging analysis of the candidates, and eventually distributes the positions of the GRBs via the Internet.

In line with pre-launch expectations, five GRBs have been detected to date in the field of view of the INTEGRAL instruments: GRB021125 (Bazzano & Paizis 2002), GRB021219 (Mereghetti, Götz & Borkowski 2002), GRB030131 (Borkowski et al. 2003), GRB030227 (Götz et al. 2003a), and GRB030320 (Mereghetti, Götz & Borkowski 2003). Their positions have been determined with errors in the 2 to 4 range. Here we report on the observations of GRB030227, the first GRB for which a prompt localization with INTEGRAL has led to the successful identification of X–ray (Loiseau et al. 2003) and optical afterglows (Castro-Tirado et al. 2003a, Soderberg et al. 2003).

## 2 INTEGRAL observation

IBIS detected GRB 030227 with its low energy detector ISGRI (Lebrun et al. 2001), an array of 128$`\times `$128 CdTe crystals sensitive in the energy range from $``$15 keV to $``$300 keV. ISGRI has an effective area of the order of 1000 cm2 and provides an angular resolution of $``$12 over a 29$`{}_{}{}^{}\times `$29 field of view. Bright sources can be located with good accuracy (for example, the 90% confidence level error radius for a source with signal to noise ratio of $``$10 is as small as 1).

The SPI instrument observes with a coarser angular resolution ($``$2) the same region of sky covered by IBIS, providing a better energy resolution (FWHM=2.5 keV at 1 MeV). The detector consists of 19 Ge crystals cooled to 85 K, surrounded by a thick anticoincidence shield of BGO scintillation crystals.

GRB 030227 was detected by the IBAS program based on the analysis of the IBIS/ISGRI data, consisting of positional, timing and energy information of each detected photon. The IBAS alert message with the preliminary position of the burst (R.A. = 4h 57m 24s, Dec.= +20 28 24′′, J2000) was delivered to the members of the IBAS Team at 08:42:38 UT of February 27, 2003, only 35 s after the start of the burst (most of this delay was due to buffering of the telemetry on board the satellite and to data transmission between the ground station and the ISDC). Unfortunately, the Internet message with these coordinates could not be distributed in real time333The detection of GRB 030227 occurred during a calibration observation of the Crab Nebula. Since the instrument configuration during these observations caused some spurious IBAS triggers, the automatic delivery of the alerts to external clients had been temporarily disabled. . Nevertheless, this information was distributed within less than one hour, after the GRB had been confirmed by an interactive quick look analysis of the data (Götz et al. 2003a). Further analysis resulted in an improved localization at R.A. = 4h 57m 32.2s, Dec. = +20 29 54′′ (Götz et al. 2003b), with an error of 3 dominated by systematic uncertainties. This position differs by only 50′′ from that of the optical transient (Castro-Tirado et al. 2003).

The IBIS light curve of GRB 030227 is shown in Fig. 1. The burst, which started approximately at 08:42:03 UT, had a typical FRED (Fast Rise and Exponential Decay) profile, with a rise to the peak in 2 s and a decay well described by an exponential with e-folding decay time of 11$`\pm `$1 s. The peak flux, integrated over 1 s, is 1.1 photons cm-2 s-1 in the 20-200 keV energy range. The fluence, in the same range and assuming the average spectrum discussed below, is 7.5$`\times `$10-7 erg cm-2.

GRB 030227 was detected with a signal to noise ratio of $``$20 at off-axis angles Z=5.3, Y=6.9. Since a fully calibrated response matrix valid at this off-axis angles is not yet available, we derived the spectrum of GRB 030227 by comparing its count rate in different energy bins to the corresponding values obtained from the Crab Nebula observed at a similar position in the field of view. With this procedure we extracted the average burst spectrum for the time interval 08:42:04 to 08:42:26 UT, which is well described by a single power law with photon index 1.85$`\pm `$0.2 over the 20-200 keV energy range.

GRB 030227 was also detected by SPI with a significance of 7.7 $`\sigma `$ in the 20-200 keV range and localized (R.A.= 74.766, Dec.= +20.531) only 0.36 off the IBIS position, owing to the less-accurate location precision of SPI. A spectrum was extracted for a time interval of 18 s starting at 08:42:04 UT. The background was estimated from a time interval of 35 min around the GRB, but excluding the time span of the event itself. A power law with photon index of 1.9$`\pm `$0.3 and 20-200 keV flux of 4.7$`\times `$10-8 erg cm-2 s-1 gives a good fit to these data. These values are consistent with those obtained with IBIS, confirming that the method used in the IBIS spectral analysis does not introduce important systematic effects.

To study the spectral evolution as a function of time, we defined a hardness ratio HR = (H-S)/(H+S) based on the background subtracted count rates in the ranges H=40-100 keV and S=20-40 keV. From the HR values of both SPI and IBIS we found evidence for a hard to soft evolution. A time resolved spectral analysis gave the power law indices shown in Fig. 2, which confirm the indication of a possible spectral softening during the decaying part of the event.

The overall veto count rate of SPI’s anticoincidence shield (ACS) showed no evident count rate increase which could be associated with the GRB event. This is consistent with the fluxes quoted above and the low effective area of the ACS for directions corresponding to the SPI field of view (von Kienlin et al. 2001).

## 3 XMM-Newton observation

XMM-Newton observed the position of GRB 030227 for $``$13 hours, starting on February 27 at 16:58 UT. Due to the presence of a high background induced by low energy protons of solar origin, the observation had to be interrupted twice, resulting in net exposure times of 33 and 36 ks respectively in the PN and MOS cameras of the EPIC instrument (Strüder et al. 2001, Turner et al. 2001). All the cameras operated in Full Frame mode and with the thin optical blocking filter. The data were processed using SAS version 5.4.1.

A bright and variable source was clearly detected at R.A. = 04h57m 33.1s, Dec. = +20 29 05′′ (error radius of 4′′), well within the IBIS error circle of GRB 030227 . The flux variability immediately suggested that this source was the GRB afterglow (Loiseau et al. 2003), as it was later confirmed by the discovery of a fading optical transient within its small error region. Fig. 3 shows the background subtracted PN light curve in the 0.2-10 keV energy range. The figure also shows the first flux measurement (at $``$17:00 UT) obtained with the MOS camera, which started the observation earlier than the PN. The X-ray flux decay is well described by a power law function, F $``$ t, with $`\delta `$=0.97$`\pm `$0.07.

The source spectra were extracted from a circle of radius 40′′ and rebinned to have at least 20 counts per channel and to oversample by a factor three the instrumental energy resolution. For the PN the background was extracted from a nearby rectangular region (1.5$`\times `$2) in the same chip as the source and for the MOS cameras from an annulus (inner and outer radii of 1 and 2) centered on the source. All the errors on the spectral parameters given below are at the 90% confidence level.

After checking that no significant variations in the best fit parameters, except for the flux value, occurred during the observation, we performed a spectral fit of the whole PN data set. A model consisting of an absorbed power law resulted in an unacceptable fit ($`\chi ^2`$/dof=272/209) with photon index $`\mathrm{\Gamma }`$ = 2.04$`\pm `$0.05 and NH = (3.6$`{}_{0.2}{}^{+0.3})\times 10^{21}`$ cm-2, larger than the galactic absorption in this direction (NH=1.75$`\times `$1021 cm-2, Hartmann & Burton 1997). An acceptable fit ($`\chi ^2`$/dof=235/208, Fig. 4) could be obtained by fixing NH to the galactic value and adding a redshifted neutral absorption, NHz. This resulted in $`\mathrm{\Gamma }`$ = 1.94$`\pm `$0.05, NHz =(6.8$`{}_{3.8}{}^{+1.8})\times 10^{22}`$ cm-2, $`z`$=3.9$`\pm `$0.3, and 0.2-10 keV observed flux Fx=8.5$`\times 10^{13}`$ erg cm-2 s-1. As shown in Fig. 5, the redshift is not strongly constrained, its value being correlated with that of NHz. The best fit values are only slightly dependent on the assumed galactic NH value: by varying it by $`\pm `$30% acceptable fits were always obtained with $`\mathrm{\Gamma }`$2, while the 90% c.l. range of $`z`$ and NHz varied by less than 35%.

The fit residuals of Fig. 4 suggest the presence of possible lines in the spectrum. We therefore tried to improve the fit by adding to the above model gaussian lines at different energies and fixed widths, smaller than the instrumental resolution. The only possibly significant line was found at observed energy of 1.67$`{}_{0.03}{}^{}{}_{}{}^{+0.01}`$ keV. According to an F-test this line is significant at the 3.2 $`\sigma `$ level. However, such a probability value should be used with caution (Protassov et al. 2002), and we give it just to allow a comparison with other possible detections of lines reported for previous GRB’s.

We also tried thermal models, both for the whole observation and for shorter time intervals, replacing the power law component with a redshifted optically thin plasma model (MEKAL in XSPEC, with $`z`$ linked to that of NHz). Keeping the elemental abundances fixed to the solar values resulted in unacceptable fits. Letting the abundances be free parameters, we obtained slightly better results ($`\chi ^2`$/dof=260/207), with values of NHz and $`z`$ similar to the power law case, and a temperature of $``$25 keV. However, the corresponding $`\chi ^2`$/dof values were always greater than in the power law model.

Entirely consistent results were obtained from the MOS spectra. As mentioned above, a short initial exposure was carried out only with the MOS, from 16:58 to 17:16 UT. These data do not provide enough statistics for a detailed spectral analysis, but they are consistent with the spectral parameters obtained for the rest of the observation. The corresponding flux is 1.34$`\times 10^{12}`$ erg cm-2 s-1.

## 4 Discussion

The INTEGRAL data indicate that GRB 030227 belongs to the class of long GRBs. Its peak flux and fluence, converted to the 50-300 keV energy range used in the BATSE catalog, are respectively $``$0.6 photons cm-2 s-1 and $``$7$`\times `$10-7 erg cm-2 (these values slightly depend on the extrapolation; for example, a break in the spectrum at 120 keV to a slope of 2.5 would reduce them by 15%). About three quarter of the GRBs in the 4th BATSE Catalog (Paciesas et al. 1999) have a higher peak flux, indicating that GRB 030227 is quite faint. The possible evidence for a hard to soft spectral evolution in GRB 030227 suggests that such a trend, already observed in samples of brighter GRBs (Ford et al. 1995, Preece et al. 1998, Frontera et al. 2000), may also apply to relatively faint bursts.

GRB spectra usually show a curvature requiring more complex models than a single power law. A broken power law, with a smooth transition between the low energy ($`\alpha `$) and high energy ($`\beta `$) slope provides in most cases an adequate empirical description of the data (Band et al. 1993). By fitting such a function to our data no improvement of the fit was obtained. A spectral break, if any, could occur outside the energy range over which we detected GRB 030227 , most likely above $``$100-200 keV as seen in most GRBs. Our photon index $``$1.9 lies in the lower tail of the distribution of values of $`\alpha `$ (Preece et al. 2000), indicating that GRB 030227 is relatively soft. An extrapolation of its spectrum to lower energies, leads to a ratio of fluences in the 6-30 to 30-400 keV ranges of 0.45 (or more if there is a high-energy break). This would qualify GRB 030227 as an X–ray rich GRB (see, e.g., Barraud et al. 2003). In this context it is interesting to note that its optical afterglow is particularly faint, R$``$23.3 at t-to = 12 hr (Gorosabel et al. 2003), comparable to the dimmest GRB afterglows found so far.

The fading behavior of the source detected with EPIC, and its positional coincidence with the optical transient, qualify it as the X–ray afterglow of GRB 030227 . Among the afterglows detected so far by XMM-Newton , this one is the brightest and the one observed with the shortest delay after the prompt event (t-to=8 hr). The EPIC spectrum reported here has the highest statistical quality among the afterglows detected to date with XMM-Newton . A power law with a redshifted absorber provides a statistically acceptable fit to the spectrum. We found tentative evidence for an emission line at 1.67 keV. If this line is due to Fe, as it has been suggested for similar features observed in other afterglows (see, e.g., Piro 2002), the implied redshift of $``$2.7-3 (depending on the Fe ionization state) would be consistent with the value derived from the absorption. Other interpretations in terms of lighter elements are of course possible, but somewhat arbitrary in the lack of an independent measure of $`z`$ and/or other statistically significant lines.

Low energy emission lines found in the initial part (t-t$`{}_{o}{}^{}`$11-13 hr) of the GRB 011211 afterglow (Reeves et al. 2002, but see also Borozdin and Trudolyubov 2003) provide evidence for thermal emission from a plasma enriched in metals. Evidence for thermal emission in two other afterglows observed with XMM-Newton , GRB 010220 (t-t$`{}_{o}{}^{}`$15 hr) and GRB 001025A (t-t$`{}_{o}{}^{}`$45 hr), was reported by Watson et al. (2002a). On the other hand, the XMM-Newton spectra of the afterglow of GRB 020322 starting at t-t$`{}_{o}{}^{}`$15 hr (Watson et al. 2002b) were adequately fit with a power law. These results suggest that the presence of thermal components is not an ubiquitous property of all GRBs and/or it might be related only to some short duration phases of the afterglows.

An upper limit to the redshift of $`z`$3.5 has been estimated from the absence of the onset of the Lyman forest blanketing in the optical data (Castro-Tirado et al. 2003b). This allows us to constrain the $`z`$ values derived by the X–ray spectral fitting (see Fig. 5), which is in any case $`>`$1, supporting a non local origin of GRB 030227 . The corresponding value of the local absorption NHz of the order of a few times 1022 cm-2 supports the scenarios involving the occurrence of GRBs in regions of star formation (e.g. Galama & Wijers 2001).

Finally we note that the rapid IBAS localization of GRB 030227 leading to the successful detection of its X–ray and optical afterglows, as well as the rate of $``$1 GRB per month in the IBIS field of view observed so far, demonstrate that INTEGRAL will efficiently complement other satellites specifically devoted to GRB studies. We also expect that particularly interesting results will be obtained for a few bursts falling in the central part of the field of view, which is covered also by the INTEGRAL monitors operating in the optical (V band) and in the X-ray range (4.5 - 35 keV).

We thank the XMM-Newton staff for the prompt execution of the TOO observation. The IBAS development has been funded by the Italian Space Agency. This research has been partially supported by the Spanish Program AYA2002-0802 (including FEDER funds) and by Polish grant 2P03C00619p02 from KBN.