text:
# Kinetic Equations for Diffusion in the Presence of Entropic Barriers

## 1 Introduction

The free energy landscape of a complex system, when represented as a function of an order parameter or reaction coordinate, presents an intricate aspect consisting of multiple local minima, defining metastable states, separated by barriers. The nature of the barriers depends on which thermodynamic potential varies when passing from one well to the other and their presence plays an important role in the dynamics of the system.Whereas energy barriers are more frequent in problems of solid state physics, entropy barriers are often encountered in soft condensed matter and biological systems. These barriers may also appear when coarsening the description of a complex system in order to simplify its dynamics. The elimination of some coordinates can be performed by introducing a configurational entropy accounting for the degeneracy (i.e. the number of equivalent configurations) of the state defined by the remaining coordinates. This procedure has been used for instance in the context of glasses .

The study of the dynamic properties of the system at the mesoscopic level requires the knowledge of the probability distribution function of the relevant degrees of freedom which evolves according to a kinetic equation of the Fokker-Planck type. This equation has been usually derived from kinetic theory in the diffusion limit -, from a master equation , and from projection operator techniques .

Our purpose in this paper is to propose a new and simple theoretical framework to analyze the mesoscopic dynamics of systems in the presence of energetic and entropic barriers. The procedure we follow is based on the application of mesoscopic nonequilibrium thermodynamics -, and has been previously used to analyze transport and activated processes in systems governed by a dynamics of the Fokker-Planck type. The kinetic equation follows after obtaining the probability current, occurring in the corresponding continuity equation in the configuration space, from the entropy production accounting from dissipation in that space related to the underlying diffusion process of the probability density. That equation applies irrespectively of the nature of the equilibrium state of the system. When the system is in contact with a heat bath and its volume remains unaltered, in which case the proper thermodynamic potential is the free energy, our equation describes the dynamics in the presence of energetic and entropic barriers. The validity of our method, however, goes beyond this standard situation and embraces other cases of interest characterized by other thermodynamic potentials and statistical ensembles.

The paper is distributed in the following way. In Section II, we present a general derivation of the kinetic equation in the framework of mesoscopic nonequilibrium thermodynamics. That equation governs the evolution of the probability density in the presence of barriers of any nature. In Section III we discuss an example in which the presence of an entropic barrier manifests: diffusion in a channel of varying cross section. The corresponding kinetic equation has been referred to as the Fick-Jacobs equation. Section IV is devoted to the derivation of that equation from mesoscopic nonequilibrium thermodynamics and to the formulation of a scaling law for the diffusion coefficient which is compared to previous results. In Section V we extend the results of Section IV to the case in which an external field acts on the diffusing particle. Finally, in Section VI we present our main conclusions.

## 2 Kinetic equations in the presence of potential barriers

Our purpose in this section is to present a general derivation of the kinetic equations describing the evolution in time of the probability density when the statics of the system is characterized by a thermodynamic potential landscape.

We consider that the equilibrium state of the system is described by the probability distribution function

$$P_{eq}\mathrm{exp}\left(\frac{\mathrm{\Delta }𝒲(\underset{¯}{x})}{kT}\right),$$
(1)
where $`𝒲(\underset{¯}{x})`$ is the minimum reversible work required to change the state of the system , which is related to the maximum useful work which can be extracted from it, sometimes referred to as available energy  or availability ; $`k`$ is Boltzmann’s constant, $`T`$ is the temperature of the medium and $`\underset{¯}{x}`$ an arbitrary set of coordinates which may represent the velocity of a particle, the orientation of a spin, the size of a macromolecule or whatever coordinate or order parameter whose values defines the state of the system in a phase space. Variations of the minimum work for a thermodynamic system are expressed as

$$\mathrm{\Delta }𝒲=\mathrm{\Delta }UT\mathrm{\Delta }S+p\mathrm{\Delta }V\mu \mathrm{\Delta }N+y\mathrm{\Delta }Y+\mathrm{},$$
(2)
where $`U`$ is the internal energy, $`S`$ the entropy, $`V`$ the volume, and $`N`$ the number of particles of the system, whereas $`T`$, $`p`$, and $`\mu `$ are the temperature, pressure and chemical potential of the environment. The term $`y\mathrm{\Delta }Y`$ represent other kinds of work (electric, magnetic, surface work…) performed on the system, being $`y`$ the intensive parameter and $`Y`$ its conjugate extensive variable . For the case of constant temperature, volume and number of particles, the minimum work corresponds to the Helmholtz free energy $`F`$. In general, that quantity reduces to the thermodynamics potentials by imposing to expression (2) the corresponding constraints  .

Having specified the statics of the system we will derive the kinetic equation describing the evolution of the probability density. To this purpose, we will use the framework of mesoscopic nonequilibrium thermodynamics. This theory applies the scheme of nonequilibrium thermodynamics to describe the dynamics of mesoscopic degrees of freedom. The treatment of a diffusion process in the framework of nonequilibrium thermodynamics is extended to the case in which the relevant quantity is a probability density, defined in phase space, instead of a mass density. The starting point is then the formulation of the Gibbs equation

$$\delta S=\frac{1}{T}\mu (\underset{¯}{x})\delta P(\underset{¯}{x},t)𝑑\underset{¯}{x},$$
(3)
which resembles the corresponding law proposed in nonequilibrium thermodynamics in terms of the mass density of particles. Here $`\mu (\underset{¯}{x})`$ is a generalized chemical potential whereas $`P(\underset{¯}{x},t)`$ is the probability density.

This expression is compatible with the Gibbs entropy

$$S(t)=kP(\underset{¯}{x},t)\mathrm{ln}\frac{P(\underset{¯}{x},t)}{P_{eq}(\underset{¯}{x})}d\underset{¯}{x}+S_{eq},$$
(4)
where $`S_{eq}`$ is the entropy when the system and the heat bath are at equilibrium. Effectively, taking variations in Eq. (4) and using the expression of the chemical potential per particle

$$\mu (\underset{¯}{x})=kT\mathrm{ln}\frac{P(\underset{¯}{x},t)}{P_{eq}(\underset{¯}{x})}+\mu _{eq},$$
(5)
where $`\mu _{eq}`$ is the chemical potential at equilibrium, one arrives at expression (3).

From the Gibbs equation we can obtain the entropy production related to the underlying diffusion process in $`\underset{¯}{x}`$-space. Calculating the time derivative of the entropy from Eq. (4) one has

$$\frac{S}{t}=k\mathrm{ln}\frac{P}{P_{eq}}\frac{P}{t}d\underset{¯}{x}.$$
(6)
Quite generally, we may assume that the probability density satisfies the continuity equation

$$\frac{P}{t}=\frac{}{\underset{¯}{x}}\underset{¯}{J},$$
(7)
where $`\underset{¯}{J}(\underset{¯}{x},t)`$ is a current defined in $`\underset{¯}{x}`$-space which has to be determined. Substituting now the time derivative of Eq. (7) into Eq. (6) one obtains the entropy production

$$\sigma =k\underset{¯}{J}(\underset{¯}{x},t)\frac{}{\underset{¯}{x}}\mathrm{ln}\frac{P}{P_{eq}}d\underset{¯}{x}.$$
(8)
To arrive at Eq. (8) we have also performed a partial integration and assumed that the current vanishes at the boundaries of the system in $`\underset{¯}{x}`$-space.

In the nonequilibrium thermodynamics scheme, the entropy production (8) consists of contributions of products between the current $`\underset{¯}{J}(\underset{¯}{x},t)`$ and its conjugated thermodynamic force $`k\frac{}{\underset{¯}{x}}\mathrm{ln}\frac{P}{P_{eq}}`$, for each value of the coordinate $`\underset{¯}{x}.`$ Assuming locality in $`\underset{¯}{x}`$-space, for which only currents and forces at the same value of $`\underset{¯}{x}`$ are coupled, one obtains the linear law

$$\underset{¯}{J}=kL(\underset{¯}{x})\frac{}{\underset{¯}{x}}\mathrm{ln}\frac{P}{P_{eq}},$$
(9)
in which the phenomenological coefficient $`L(\underset{¯}{x})`$ may in general depend on $`\underset{¯}{x}`$.

The resulting kinetic equation then follows by substituting Eq. (9) into the continuity equation (7)

$$\frac{P}{t}=\frac{}{\underset{¯}{x}}\left(DP\frac{}{\underset{¯}{x}}\mathrm{ln}\frac{P}{P_{eq}}\right),$$
(10)
where we have defined the diffusion coefficient

$$D(\underset{¯}{x})=\frac{kL(\underset{¯}{x})}{P}.$$
(11)
This equation, which in view of Eq. (1) can also be written as

$$\frac{P}{t}=\frac{}{\underset{¯}{x}}\left(D\frac{P}{\underset{¯}{x}}+\frac{D}{kT}\frac{\mathrm{\Delta }𝒲}{\underset{¯}{x}}P\right),$$
(12)
is the Fokker-Planck equation accounting for the evolution of the probability density in $`\underset{¯}{x}`$-space.

Under the conditions for which $`𝒲=F=UTS`$, this equation transforms into the Fokker-Planck equation for a system in the presence of energy and entropy barriers. One then obtains

$$\frac{P}{t}=\frac{}{\underset{¯}{x}}\left(D\frac{P}{\underset{¯}{x}}+\frac{D}{kT}\frac{\mathrm{\Delta }U}{\underset{¯}{x}}P\frac{D}{k}\frac{\mathrm{\Delta }S}{\underset{¯}{x}}P\right),$$
(13)
where the drift term consists of contributions due to an external potential and to variations of the entropy. When the barrier is purely energetic, only the first contribution remains, and the previous equation reduces to the well-known Fokker-Planck equation

$$\frac{P}{t}=\frac{}{\underset{¯}{x}}\left(D\frac{P}{\underset{¯}{x}}+\frac{D}{kT}\frac{\mathrm{\Delta }U}{\underset{¯}{x}}P\right).$$
(14)
If the nature of the barrier is purely entropic, Eq. (13) then corresponds to the Fick-Jacobs equation

$$\frac{P}{t}=\frac{}{\underset{¯}{x}}\left(D\frac{P}{\underset{¯}{x}}\frac{D}{k}\frac{\mathrm{\Delta }S}{\underset{¯}{x}}P\right),$$
(15)
first derived in the context of diffusion of Brownian particles in a channel of non-constant cross section .

The general form of Eq. (12), in which the equilibrium distribution function does not need to be specified and is given in general by Eq. (1), makes that result applicable to a great diversity of situations. For example, for a Brownian particle, for which the minimum work is simply its kinetic energy, that equation correspond to the usual Fokker-Planck equation. When $`\underset{¯}{x}`$ is an order parameter or a reaction coordinate, Eq. (12) provides the kinetic equation in the presence of barriers. The method used in this section then offers a common formalism able to analyze the dynamics of a system in the presence of energy and entropy barriers. In this paper, we will mainly focus our analysis to the case in which the potential is strictly entropic.

## 3 Diffusion through a channel of varying cross section: the Fick-Jacobs equation

The influence that the presence of an entropic barrier exerts on the dynamics of a system can be illustrated by deriving the kinetic equation governing Brownian diffusion through a channel of varying cross-section. This equation was first proposed by Jacobs  and subsequently rederived by Zwanzig  on the basis of more fundamental arguments.

### 3.1 Jacobs’ derivation

Jacobs in his book *Diffusion Processes*  provided an heuristic derivation of the equation governing diffusion in a symmetric tube whose cross section $`A(x)`$ varies along the $`x`$ axis, defined by the center line of the tube. The argument runs as follows.

Consider an elementary volume of thickness $`dx`$ perpendicular to the axis of the tube. The total amount of particles in this slice, at $`x`$ and time $`t`$ is $`C(x,t)dx`$, which is the integral of the concentration over the volume $`A(x)dx`$. The rates of entrance and exit of the diffusing substance into this volume are given by the Fick’s law $`J=D_0\frac{C/A}{x}`$, where $`D_0`$ is the diffusion coefficient, and $`C/A`$ is the local volume concentration. In this case, both rates are different not only because the concentration gradient $`\frac{C}{x}`$ changes with the distance, but also due to the variation of the cross-section of the channel. Explicitly, the rate of entrance into this elementary volume is $`D_0A\frac{C/A}{x}`$ , whereas the rate of exit is given by $`D_0\left(A\frac{C/A}{x}+\frac{}{x}\left(A\frac{C/A}{x}\right)dx+\mathrm{}\right)`$. The difference between both rates provides the rate of change of the substance in the elementary volume, which can also be expressed as $`\frac{C}{t}dx`$. Neglecting quadratic terms in $`dx`$, one easily arrives at the equation governing the diffusion in the channel

$$\frac{C}{t}=D_0\frac{}{x}\left(A\frac{}{x}\left(\frac{C}{A}\right)\right)=D_0\frac{}{x}\left(\frac{C}{x}\frac{1}{A(x)}\frac{dA(x)}{dx}C\right)$$
(16)
which is referred to as the Fick-Jacobs equation.

### 3.2 Zwanzig’s derivation

Zwanzig reported in Ref.  a more general and rigorous derivation of the Fick-Jacobs equation. For the sake of simplicity, let us review the 2D case. The method consists of performing a reduction in the number of coordinates from the 2D Smoluchowski equation to a 1D description.

The starting point is the two-dimensional Smoluchowski equation for diffusion through a general potential $`U(x,y)`$

$`{\displaystyle \frac{c(x,y,t)}{t}}=D_0{\displaystyle \frac{}{x}}e^{\beta U(x,y)}{\displaystyle \frac{}{x}}e^{\beta U(x,y)}c(x,y,t)`$ (17)
$`+D_0{\displaystyle \frac{}{y}}e^{\beta U(x,y)}{\displaystyle \frac{}{y}}e^{\beta U(x,y)}c(x,y,t)`$ (18)
where $`\beta =1/kT`$, and $`c(x,y,t)`$ is the concentration. To perform the reduction to 1D, this equation is integrated over the variable $`y`$, leading to

$$\frac{C(x,t)}{t}=D_0\frac{}{x}e^{\beta U(x,y)}\frac{}{x}e^{\beta U(x,y)}c(x,y,t)𝑑y$$
(19)
where the reduced concentration is defined as

$$C(x,t)=c(x,y,t)𝑑y.$$
(20)
The key point of the derivation is the assumption of equilibration in the transverse direction. Under this assumption, one can define an averaged $`x`$-dependent free energy $`F(x)`$ through the expression

$$e^{\beta F(x)}=e^{\beta U(x,y)}𝑑y$$
(21)
from which one can define a normalized conditional probability distribution

$$\rho (x;y)=\frac{e^{\beta U(x,y)}}{e^{\beta F(x)}}.$$
(22)
Then, under the local equilibrium approximation, the concentration $`c(x,y,t)`$ factorizes as follows

$$c(x,y,t)C(x,t)\rho (x;y)$$
(23)
Taking these considerations into account, one finally obtains

$$\frac{C(x,t)}{t}\frac{}{x}D_0e^{\beta F(x)}\frac{}{x}e^{\beta F(x)}C(x,t)$$
(24)
which constitutes the generalization of the Fick-Jacobs equation (16) to the case of a two dimensional potential $`U(x,y)`$.

In the previous analysis, we have not taken into account the fact that the normal flux must vanish at the boundaries. The role played by this zero normal flux conditions, can be replaced by the confining potential $`U(x,y)`$. In fact, if $`U(x,y)`$ is a box-like potential, i.e. zero inside the tube and infinite outside, it is clear that $`e^{\beta F(x)}=2w(x)`$, being $`w(x)`$ the half width of the tube. In this case, the barrier is purely entropic, and one recovers the usual Fick-Jacobs equation.

The extension to the three-dimensional case is quite straightforward by taking into account that in that case the integration in the transverse coordinates involves variables $`y`$ and $`z`$, and the width of the 2D tube must be replaced by the transverse section $`\pi w(x)^2`$.

Zwanzig’s analysis clearly manifests that the accuracy of the Fick-Jacobs equation is conditioned to the existence of local equilibrium in the transverse direction. The author analyzed the effect of the deviations from the local equilibrium related to the variations of the density $`\delta c(x,y,t)c(x,y,t)C(x,t)\rho (x;y)`$. He derived an equation for the evolution of this deviations, which suggested that the accuracy of the Fick-Jacobs equations is restricted to situations verifying $`\left|w^{}(x)\right|<1`$, that is, when the section of the tube varies smoothly. In addition, he showed that the range of validity of the one-dimensional Fick-Jacobs description could be extended by introducing a position-dependent effective diffusion coefficient

$$D_Z(x)=D_0\frac{1}{1+\gamma w^{}(x)^2}$$
(25)
where $`D_0`$ is the molecular diffusion coefficient and the parameter $`\gamma `$ depends on the dimensionality. The explicit expression for this coefficient is obtained through an expansion in powers of $`w^{}(x)`$ . The result is

$$D(x)=D_0(1\gamma w^{}(x)^2+\mathrm{})$$
(26)
where

$$\gamma =\frac{\nu ^2e^{\beta V(\nu )}𝑑\nu }{e^{\beta V(\nu )}𝑑\nu },$$
(27)
being $`\nu `$ the transverse coordinate scaled by the function $`w(x)`$ (i.e. $`\nu =\frac{y}{w(x)}`$ in 2D, and $`\nu =\frac{r}{w(z)}`$ for a 3D tube of radius $`r`$ with cylindrical symmetry), and $`V(\nu )=U(x,y)`$. For a purely confining potential, the results obtained by Zwanzig are $`\gamma =1/3`$ for the 2D case, and $`\gamma =1/2`$ for a 3D tube with cylindrical symmetry. From the expansion (26) he infers the expression (25) based merely on the fact that this reconstruction of the series improves the agreement with the exact results.

## 4 The Fick-Jacobs equation from mesoscopic nonequilibrium thermodynamics

For the case of an enclosure of varying cross-section, the concept of entropic barrier is remarkably simple. At equilibrium, the density of diffusing material $`\rho _0`$ is constant. If we contract the 3D description retaining the coordinate $`x`$, the resulting 1D equilibrium distribution is

$$\rho _{eq}(x)=\rho _0𝑑y𝑑z=\rho _0A(x).$$
(28)
In this case, the diffusing particle constitute an isolated system. Therefore its corresponding thermodynamic potential obtained form Eq. (2) is simply $`\mathrm{\Delta }𝒲=T\mathrm{\Delta }S`$, where the entropic barrier is, in accordance with Eq. (1),

$$\mathrm{\Delta }S(x)=k\mathrm{ln}A(x).$$
(29)
The previous equation clearly manifests that the entropic barrier originates from the variation of the space available for the diffusing particles. Notice that Eq. (29) corresponds to the usual microcanonical definition of the entropy in terms of the number of states which in this case is simply proportional to the area of the tube.

Once identified the equilibrium distribution (or the entropic barrier), the kinetic equation describing diffusion in the presence of the entropic barrier follows from our general scheme developed in Section 2. According to Eq. (15), the kinetic equation is

$$\frac{P}{t}=\frac{}{x}\left(D(x)\frac{P}{x}\frac{D(x)}{A(x)}\frac{A(x)}{x}P\right)$$
(30)
which has the same structure as the Fick-Jacobs equation, but with a spatial dependent diffusion coefficient.

At this point, it is worth to analyze the role played by this coefficient. In Ref.  it was shown that the validity of the 1D description could be extended if we use an effective spatial-dependent diffusion coefficient. However, the agreement one can achieve with the expression (25) proposed by Zwanzig is not always satisfactory. In particular, for the 3D example discussed in Ref. , results for the current of particles obtained by using the expression of the diffusion coefficient (25) present severe discrepancies with the corresponding expression coming from the exact result of the 3D Smoluchowski equation, when the section of the tube changes abruptly.

The appearance of an entropic barrier originates from the reduction of the space to a single coordinate. But this reduction may also have implications on the form of the diffusion coefficient. The molecular diffusion coefficient in the real space, $`D_0`$, gives information about the dispersion of the displacement. In 2D, we can estimate

$$D_0\frac{(\mathrm{\Delta }𝐫)^2}{t}=\frac{\mathrm{\Delta }x^2\left(1+\left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)^2\right)}{t}$$
(31)
which manifests that the diffusion coefficient involves the displacement in both coordinates. Once we have contracted the $`y`$ coordinate, the resulting effective diffusion coefficient only depends on the dispersion in the remaining $`x`$ coordinate, that is qualitatively $`D_{eff}\frac{(\mathrm{\Delta }x)^2}{t}`$. Therefore, from Eq. (31) we can infer the behaviour $`D_{eff}\frac{D_0}{\left(1+\left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)^2\right)}`$. Previous equation provides a hint about the fact that the reduction of coordinates may involve not only the appearance of an entropic barrier, but also a scaling of the diffusion coefficient. Following this heuristic reasoning, we will propose the scaling law

$$D_s(x)=D_0\frac{1}{(1+y^{}(x)^2)^\alpha }$$
(32)
for the diffusion coefficient appearing in Eq. (30), where $`y^{}(x)=\frac{dy}{dx}`$, and $`y(x)=w(x)`$ defines the shape of the enclosure. The objective is then to test if the use of this expression is able to extend the range of validity of the 1D Fick-Jacob-type description.

The values of the scaling exponents can be estimated from the calculations performed in Ref. . We can expand our expression for the diffusion coefficient $`D_s(x)`$ in terms of $`w^{}(x)^2`$ and compare the result with the one obtained by Zwanzig (Eq. 26). From that comparison, one obtains that a reasonable choice of the scaling exponents is $`\alpha =1/3`$ for the 2D case and $`\alpha =1/2`$ for the case of a 3D tube.

In the remaining of this section, we will test this result for the two cases discussed in Ref.  which admit an exact solution. We will see that the agreement with the exact result, when we use our effective diffusion coefficient given by Eq. (32), is significantly improved.

### 4.1 Effective Diffusion Coefficient in a 2D periodic symmetric channel

The first case under scrutiny is the diffusion in a 2D periodic channel defined by

$$x=u+a\mathrm{cosh}v\mathrm{sin}u;y=v+a\mathrm{sinh}v\mathrm{cos}u$$
(33)
where diffusion occurs in the region comprised between $`0<u<2\pi `$ and $`V<v<V`$. The parameter $`\lambda a\mathrm{cosh}V`$ must be smaller than $`1`$, to avoid double valuation of the walls. Fig. 1 illustrates the shape of the tube for $`V=0.5`$ and $`a=0.5`$.

At very long times, diffusion in periodic channels can be described through an effective diffusion coefficient $`D^{}`$, accounting for the effects of passing through many constrictions in the channel. This effective coefficient can be evaluated as

$$\frac{1}{D^{}}=\frac{1}{D(x)w(x)}w(x)$$
(34)
where $``$ denotes average over one period of the potential.

For this particular geometry, an exact expression for $`D^{}`$ was derived in Ref.

$$\left(\frac{D^{}}{D_0}\right)_{exact}=\left[1+\frac{\lambda ^2}{2}\frac{\mathrm{tanh}V}{V}\right]^1$$
(35)
and was compared to the results obtained with the Fick-Jacobs equation using a constant diffusion coefficient $`D_0`$, and with the position-dependent diffusion coefficient $`D_Z(x)`$.

In Fig. 2 we have represented the results for $`D^{}`$ obtained using the different approaches mentioned above, plus the ones corresponding to our expression $`D_s(x)`$. One can easily realize that the results obtained with the expression of $`D_s(x)`$ we propose improve the ones when using $`D_Z(x)`$ and constitute an excellent approximation of the exact result.

### 4.2 Steady-state flux through a 3D hyperboloidal cone

The improvement of the results one can achieve with $`D_s(x)`$ is more remarkable for the 3D example discussed in Ref. . In this case, the shape of the tube corresponds to a 3D hyperboloidal cone as the one depicted in Fig. 3. It can be conveniently described by using oblate spheroidal coordinates $`(\xi ,\eta ,\phi ),`$ related to the cylindrical coordinates $`(r,z,\phi )`$ by

$$r^2=a^2(\xi ^2+1)(1\eta ^2),z=a\xi \eta $$
(36)
In this coordinate system, $`\xi =0`$ corresponds to the small hole and $`\xi =\mathrm{}`$ to the far end of the tube; $`\eta =1`$ is the z-axis, $`\eta =0`$ is the (x,y)-plane, whereas $`\eta =\eta _0`$ corresponds to the surface of the tube. Diffusion then takes place in the region

$$0<\xi <\mathrm{},\eta _0<\eta <1,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<\phi <2\pi $$
(37)
The steady-state 3D diffusion equation can be solved by using the boundary conditions: $`C=0`$ for $`\xi =0`$, $`C=C_0`$ for $`\xi =\mathrm{}`$, and the condition that the normal flux vanishes at the walls $`\eta =\eta _0`$. In this situation, the exact steady-state flux through the exit hole $`\xi =0`$ is

$$J_{exact}=4Da(1\eta _0)C_0$$
(38)
Alternatively, one can use the Fick-Jacobs equation derived in the previous section to calculate that flux, yielding

$$J=C_0\left[_0^{\mathrm{}}𝑑z\frac{1}{D(z)A(z)}\right]^1$$
(39)
The results are depicted in Fig. 4, corresponding to the “exact”, the unmodified Fick-Jacobs equation, in which $`D(z)=D_0`$ is a constant, and to the choices $`D_Z(z)`$ and $`D_s(z)`$.

From the examples we have discussed, it becomes clear that our proposal of a scaling law for the spatial-dependent diffusion coefficient leads to a better agreement with the exact results.

## 5 Diffusion in a channel of varying cross-section under the influence of a field

In the previous section, we have discussed the effect of entropic barriers in the dynamics. The presence of that barriers is sometimes accompanied by energy barriers. The theory we have developed admits a generalization to the case in which both an entropic and an energetic barrier coexist. This will be the purpose of this section.

An illustrative example of that situation is the problem of diffusion in a 2D channel of non-constant cross-section in the presence of gravity. The situation is depicted schematically in Fig. 5.

If gravity is involved, the 2D equilibrium distribution is given by the usual Boltzmann’s factor

$$P_{eq}(x,y)e^{\frac{mgy}{kT}}$$
(40)
Considering the normalization condition

$$1=P_{eq}(x,y)𝑑x𝑑y=P_{eq}(x)𝑑x$$
(41)
and using the expression (40), one can identify the reduced 1D equilibrium distribution as

$$P_{eq}(x)_{y_1(x)}^{y_2(x)}𝑑y\mathrm{exp}\left(\frac{mgy}{kT}\right)=\frac{kT}{mg}\mathrm{exp}\left(\frac{mgy_1(x)}{kT}\right)\left(1\mathrm{exp}\left(\frac{2mgw(x)}{kT}\right)\right)$$
(42)
where $`y_1(x)`$ and $`y_2(x)`$ define the boundaries of the channel and $`2w(x)=y_2(x)y_1(x)`$ is its width. The dynamic equation then follows after introducing $`P_{eq}(x)`$ in Eq. (10), yielding

$$\frac{P}{t}=\frac{}{x}\left(D(x)\frac{P}{x}+\frac{D(x)mg}{kT}y_1^{}(x)\mathrm{coth}\left(\frac{mgw(x)}{kT}\right)P\right)$$
(43)
for the case of a symmetric channel for which $`y_2(x)=y_1(x)`$.

From the previous expressions one can distinguish different limiting cases:

1. In the limit $`ϵ\frac{mgw(x)}{kT}1`$, we can expand the reduced equilibrium distribution (42) in terms of $`ϵ`$

$$P_{eq}(x)\frac{kT}{mg}(1+O(ϵ))\left(11+\frac{2mgw(x)}{kT}+O(ϵ^2)\right)2w(x),$$
(44)
recovering the case of purely entropic barrier discussed in the previous section. Eq. (43) reduces then to the Fick-Jacobs equation (30).
2. For the case $`ϵ1`$, gravity dominates, $`\mathrm{coth}ϵ1`$, and we recover the dynamics for diffusion along a 1D purely energetic barrier.

$$\frac{P}{t}=\frac{}{x}\left(D(x)\frac{P}{x}+\frac{D(x)mg}{kT}y_1^{}(x)P\right)$$
(45)
3. If the boundaries are flat, that is $`y_1(x)=y_2(x)=cte`$, the equilibrium distribution $`P_{eq}(x)`$ is then constant, and the dynamics can be properly described by the 1D Smoluchowski equation.

$$\frac{P}{t}=\frac{}{x}\left(D(x)\frac{P}{x}\right)$$
(46)

The evolution of this confined system in the presence of gravity dictated by Eq. (43) presents some peculiarities. On one hand, since the energy barrier depends on the coordinate we have eliminated, there exists a coupling between entropy and energy barriers and not a mere superposition of the drifts related to each of them, which is reflected by the presence of the term $`\mathrm{coth}\left(\frac{mgw(x)}{kT}\right)`$ in the drift of Eq. (43). On the other hand, it is important to highlight that Eq. (43) does not satisfy the detailed balance condition: the mobility $`b\frac{D(x)}{kT}\mathrm{coth}\left(\frac{mgw(x)}{kT}\right)`$ and the diffusion $`D(x)`$ are not related through the usual Einstein relation $`D(x)=kTb`$. That is the signature that the fluctuation-dissipation theorem, which holds at equilibrium, may lose its validity when we perform a reduction of the variables describing the state of the system.

Concerning the effective diffusion coefficient, a scaling law as the one proposed through Eq. (32), which is valid in absence of external forces, will be in general not correct when an external field is present. The value of the scaling exponent $`\alpha `$ varies in the presence of external potentials, since the scaling exponents changes due to interactions. The range of validity of the Fick-Jacobs description and the modifications in the scaling of the effective diffusion coefficient under the presence of external fields requires a more elaborated treatment.

## 6 Conclusions

In this paper we have presented a theory to describe the kinetics of a system whose equilibrium state is characterized by a given landscape of an unspecified thermodynamic potential. The theory is based on mesoscopic non-equilibrium thermodynamics, which uses the scheme of nonequilibrium thermodynamics  at the mesoscopic level of description in which the pertinent fields are probability densities.

Particularly, we have established the kinetic equation for a system in the presence of entropy barriers. The barriers may be inherent to the intimate structure of the system or may emerge as a consequence of the elimination procedure of some coordinates when one tries to simplify its dynamic description.

An illustrative example treated in the literature is the diffusion of a particle in an region of irregular geometry. The governing equation for the probability density is known as the Fick-Jacobs equation and has been derived heuristically by Jacob  and directly from the proper coordinate reduction procedure by Zwanzig . In our derivation, the entropic barrier directly comes out from Boltzmann’s entropy, which follows from the proper accounting of the number of accessible states of the system. In the kinetic equation we have obtained, the dependency of the diffusion coefficient on the coordinate follows from the general dependency of the Onsager transport coefficients on the state variable, according to the rules of nonequilibrium thermodynamics. This dependency becomes crucial in order to the reduced description accurately resembles the exact solution. We have proposed a scaling law for the diffusion coefficient reaching a very good agreement with the exact solution even in the case when the original Fick-Jacobs equation does not provide a good approximation. In this way, we have put the validity of a Fick-Jacob-like description in a more general context.

The theory we have proposed can also be applied to cases in which an energy barrier is also present. We have obtained the corresponding Fokker-Planck equation describing its dynamics. Interesting characteristics of this equation are the coupling of the entropy and energy barriers and the apparent violation of the Fluctuation Dissipation Theorem, resulting from the elimination of variables.

The theory we have presented is applicable to a wide variety of systems of different nature. Apart from the diffusion problem in the presence of an entropic barrier we have discussed in this paper, we can quote protein folding , glassy systems , transport of ions  and macromolecules - through membranes or channels, motion of polymers subjected to rigid constraints , protein binding kinetics , drug release , nucleation  or polymer crystallization , to mention just a few examples in which the presence of entropic barriers becomes relevant.

There still remain open questions whose answers goes beyond the scope of this paper. The main pitfalls arise from the coordinate reduction procedure. The question is how to proceed with the coarsening in more complex situations. For instance, it will be desirable to analyze the effect of the asymmetries either intrinsic to the nature of the landscape or due to the presence of external forces which impose a preferred direction, which occurs in many situations of practical interest. Another interesting point will be to analyze the behaviour when the slope diverges, giving rise to regions in which the particle may become trapped thus breaking the ergodicity of the system. The solutions to these problems probably demand a more elaborated analysis of the scaling of the diffusion coefficient which takes into account these effects.

Just proceeding systematically along the lines we have established by means of our theory one could arrive at a more complete theoretical framework providing the dynamic description of systems exhibiting entropic barriers.

## 7 Acknowledgments

We want to acknowledge P. Mazur, J.M. Vilar, P. Reimann and G. Schmidg for their comments and suggestions. D. Reguera wishes to thank Generalitat de Catalunya for financial support. This work has been partially supported by DGICYT of the Spanish Government under grant PB98-1258.


format: text

text:
# On non-Adiabatic Holonomoic Quantum Computer

1National Centre of High Energy Physics,Bogdanovich Str.153, Minsk 220040, Belarus
2Institute of Physics National Academy of Sciences Skoryna av.68, Minsk 220072, Belarus

Abstract
Non-adiabatic non-Abelian geometric phase of spin-$`3/2`$ system in the rotating magnetic field is considered. Explicit expression for the corresponding effective non-Abelian gauge potential is obtained. This formula can be used for construction of quantum gates in quantum computations.
PACS: 03.65.V; 03.67
Keywords: quantum computation; quantum gate; non-adiabatic non-Abelian
geometric phase

A great deal of investigations in quantum computers (see  for a bibliographic review on this subject) refreshed some interest on Berry phase effect  in quantum mechanics. The idea of using unitary evolution operators produced by the non-Abelian Berry phase  as quantum calculations is proposed in  and realized in  in a concrete optical model of holonomic quantum computer, in  in the solid state model, in the model of network of Josefson junctions , in  for the ion-trap model and in  for that of QED cavity. Calculation aspects of this model are considered in . For other references where also Abelian Berry phase is considered in the context of quantum computer see e.g.  \- . On the other hand non-adiabatical Berry phase can exist and be measured if transitions in a given statistical ensemble do not lead to loose of coherency . Therefore it is also possible to use the corresponding unitary operators to perform quantum calculations. This fact has been noticed in ,  for the NMR model of QC and in  for the ion-trap model. In this paper we show a realization of quantum gates for a concrete 4-level quantum system driven by external magnetic field. Let us consider a spin-$`3/2`$ system with quadrupole interaction. Physically it can be thought of as a single nucleus with the spin above. A coherent ensemble of such nuclei manifest geometric phase when placed in rotating magnetic field. This phase is non-Abelian due to degenerate energy levels with respect to the sign of the spin projection. Depending on the experiment setup the phase can be both adiabatic as in Rb experiment by Tycko  and non-adiabatic as in Xe experiment by Appelt et al . This non- Abelian phase results in mixing of $`\pm 1/2`$ states in one subspace and $`\pm 3/2`$ in another one and thus can be regarded as a 2-qubit gate. The gate is generated by a non-Abelian effective gauge potential $`𝐀`$ that is the subject of computation in this paper.
Let us consider a parametric quantum system described by a Hamiltonian $`H(𝑹)`$, where $`𝑹(t)`$ is a set of evolving parameters. We suppose that evolution of the Hamiltonian is determined by a unitary rotation

$$H(t)=U^{}(t)H_0(t)U(t),$$
(1)
where $`H_0(t)`$ is diagonal in some basis in the Hilbert space. Looking for solutions of the Schrödinger equation

$$i\frac{\psi }{t}=H(t)\psi (t),$$
(2)
we take rotating frame by assigning $`\stackrel{~}{\psi }(t)=U(t)\psi (t)`$ and get in such a way

$$i\frac{\stackrel{~}{\psi }}{t}=(H_0(t)iU^{}(t)\dot{U}(t))\stackrel{~}{\psi }(t).$$
(3)
Of course, this transformation generally does not help to solve equation (2) due to the fact that the algebraic structure of the coupling term $`iU^{}\dot{U}`$ can appear to be rather complicated and the last generally does not commute with $`H_0`$. However if a receipt is known how to evaluate the last term in (3), further solution of this equation is straightforward:

$$\psi (t)=e^{i\varphi _n(t)}𝒯\mathrm{exp}\left(i_0^tU^{}(\tau )\dot{U}(\tau )𝑑\tau \right)\psi (0),$$
(4)
where $`\varphi _n(t)=_0^tE_n(\tau )𝑑\tau `$ is so called dynamical phase, $`𝒯`$ denotes time-ordering and $`E_n`$ are elements of $`H_0`$ that is by definition diagonal. Of course if there is no way to find $`U^{}(\tau )\dot{U}(\tau )`$, expression (4) is useless, but there is at least one sufficiently general case in which this term can be independently computed. It is the case of a parametric quantum system with adiabatically evolving parameters. The corresponding class of solutions was discovered by Berry  for non-degenerate $`E_n`$ and was generalized by Wilczek and Zee  for the degenerate case. In this case

$$(U^{}(t)\dot{U}(t))_{ab}=(𝑨_n)_{ab}(𝑹)=i<\phi _b^{(n)}|_𝑹\phi _a^{(n)}>,$$
(5)
where $`𝑹=(R_1,\mathrm{},R_N)`$ is a set of externally driven parameters and $`\phi _a^{(n)}a=1,\mathrm{},\kappa _n`$ ($`\kappa _n`$ is the degeneration rate for the energy level $`E_n`$) are solutions of the eigenvalue problem

$$\widehat{H}_0(𝒓,𝑹(t))\phi _a^{(n)}=E_n(𝑹(t))\phi _a^{(n)}(𝒓,𝑹(t)).$$

Here $`𝒓`$ denotes the set of dynamic variables affected by $`\widehat{H}_0`$. To our knowledge other classes of quantum systems that admit such a transparent evaluation of the term $`U^{}\dot{U}`$. In this paper we consider a particular case of 4-level system investigated experimentally in . For this system a constructive derivation of the above term is possible. Evolution of the magnetic field governing the the system can be generally non-adiabatic so the solution is of more general type (4) but not of (5). It is important to realize that although the solution part corresponding to the second factor in the right hand side of formula (4) can also be called ”geometric phase”, it can not be computed by means of formula (5) because of the system’s non-adiabaticity.

We assume the condition of the $`{}_{}{}^{131}Xe`$ NMR experiment to be held so one does not need to trouble about the coherency in the system. The last is described by the following Hamiltonian in the frame where the magnetic field is parallel to the z-axis ($`\mathrm{}=1`$)

$$H_0=\omega _0(J_3^21/3j(j+1))$$
(6)
For a spin-$`3/2`$ system we choose the third projection of the angular momentum in the form

$$J_3=\left(\begin{array}{cccc}3/2& 0& 0& 0\\ 0& 3/2& 0& 0\\ 0& 0& 1/2& 0\\ 0& 0& 0& 1/2\end{array}\right)=\left(\begin{array}{cc}3/2\sigma _3& 0\\ 0& 1/2\sigma _3\end{array}\right),$$
(7)
Then two other projection operators are

$$J_1=\left(\begin{array}{cccc}0& 0& \sqrt{3}/2& 0\\ 0& 0& 0& \sqrt{3}/2\\ \sqrt{3}/2& 0& 0& 1\\ 0& \sqrt{3}/2& 1& 0\end{array}\right)=\left(\begin{array}{cc}0& \frac{\sqrt{3}}{2}\hfill \\ \frac{\sqrt{3}}{2}\hfill & \sigma _1\end{array}\right),$$
(8)
$$J_2=\left(\begin{array}{cccc}0& 0& 0& \sqrt{3}/2\\ 0& 0& \sqrt{3}/2& 0\\ 0& \sqrt{3}/2& 0& i\\ \sqrt{3}/2& 0& i& 0\end{array}\right)=\left(\begin{array}{cc}0& i\frac{\sqrt{3}}{2}\sigma _3\\ i\frac{\sqrt{3}}{2}\sigma _3& \sigma _2\end{array}\right).$$
(9)
In the laboratory frame the Hamiltonian takes the form

$$H=\omega _0((𝑱𝒏)^21/3j(j+1))=e^{i\phi J_3}e^{i\theta J_2}H_0e^{i\phi J_2}e^{i\theta J_3}.$$
(10)
Rotation around the z-axis means that $`\phi =\omega _1t`$ and one should perform the unitary transformation

$$|\psi >=U_1|\stackrel{~}{\psi }>,U_1=e^{i\omega _1tJ_3}.$$
(11)
In the rotating frame we get

$$H_1=e^{i\theta J_2}(\omega _0J_3^2\omega _1\stackrel{~}{J}_3)e^{i\theta J_2}\frac{5\omega _0}{4},$$
(12)
where

$$\stackrel{~}{J}_3=e^{i\theta J_2}J_3e^{i\theta J_2}$$

Expression (12) is equivalent to

$$H_1=\left(\begin{array}{cc}\omega _0\frac{3}{2}\omega _1\mathrm{cos}\theta \sigma _3& \omega _1\frac{\sqrt{3}}{2}\\ \omega _1\frac{\sqrt{3}}{2}& \omega _0\frac{1}{2}\omega _1\mathrm{cos}\theta \sigma _3+\omega _1\mathrm{sin}\theta \sigma _1\end{array}\right)$$
(13)
It is convenient to diagonalize this matrix in two steps. First we get rid of $`\sigma _1`$ in the last matrix element by performing of the block-diagonal transformation

$$U_2=diag(1,e^{i\alpha \sigma _3}),$$
(14)
where $`\mathrm{tan}\alpha =2\mathrm{tan}\theta `$. Thereafter the Hamiltonian $`H_1`$ reads

$$H_1=\left(\begin{array}{cc}\omega _0\frac{3}{2}\omega _1\mathrm{cos}\theta \sigma _3& \omega _1\frac{\sqrt{3}}{2}\\ \omega _1\frac{\sqrt{3}}{2}& \omega _0\frac{1}{2}\omega _1\frac{\mathrm{cos}\theta }{\mathrm{cos}\alpha }\sigma _3\end{array}\right).$$
(15)
At the second step we apply the transformation

$$\left(\begin{array}{cc}\beta _1& \beta _2\\ \beta _2^{}& \beta _1^{}\end{array}\right),$$
(16)
where $`\beta _1`$, $`\beta _2`$ are diagonal $`2\times 2`$ matrices that must obey the unitarity condition

$$|\beta _1|^2+|\beta _2|^2=1.$$
(17)
Supposing $`\beta _{1,2}`$ to be real and performing transformation (16) we come to the diagonalization condition in the form

$$\xi (\beta _1^2\beta _2^2)+(\lambda _1\lambda _2)\beta _1\beta _2=0,$$
(18)
where $`\lambda _1,\lambda _2`$ are $`2\times 2`$ diagonal matrices and $`\xi `$ is a parameter

$`\lambda _1`$ $`=`$ $`\omega _03/2\omega _1\mathrm{cos}\theta \sigma _3,`$ (19)
$`\lambda _2`$ $`=`$ $`\omega _01/2\omega _1{\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{cos}\alpha }}\sigma _3`$ (20)
$`\xi `$ $`=`$ $`\omega _1\sqrt{3}/2\mathrm{sin}\theta `$ (21)
Assuming $`\beta _2=\mu \beta _1`$ where $`\mu `$ is a diagonal $`2\times 2`$ matrix as well we come to the following expressions for the matrix elements of $`\mu `$

$$\mu _i=k_i+\sqrt{1+k_i^2},$$
(22)
where

$$k_i=\frac{\mathrm{\Delta }\lambda _i}{2\xi },\mathrm{\Delta }\lambda _i=\lambda _{1i}\lambda _{2i}.$$

Finally we get for the matrix elements of $`\beta _{1,2}`$

$`\beta _{1i}^2`$ $`=`$ $`1/2(1+k_i^2)^{1/2}\left(k_i+\sqrt{1+k_i^2}\right)^1`$ (23)
$`\beta _{2i}^2`$ $`=`$ $`1/2\left(1+{\displaystyle \frac{k_i}{\sqrt{1+k_i^2}}}\right)`$ (24)
Now one can evaluate the connection 1-form. It is convenient to represent it as follows:

$$𝑨=iU^+dU=Ad\varphi =\left(\begin{array}{cc}A_{3/2}& A^{tr}\\ \stackrel{~}{A}^{tr}& A_{1/2}\end{array}\right)d\varphi ,$$
(25)
where all matrix elements of $`A`$ denote $`2\times 2`$ matrix-valued blocks, $`U=U_1U_2U_3`$ and $`U_i`$ are determined by (11), (14), (16) correspondingly. Here tilde denotes a transposed matrix. After some algebra we get for the matrix elements of (25)

$`A^{tr}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\beta _1\beta _2(3\mathrm{cos}\alpha )\sigma _3+{\displaystyle \frac{1}{2}}\mathrm{sin}\alpha \beta _2\sigma _1\beta _1,`$ (26)
$`A_{3/2}`$ $`=`$ $`(a_{3/2}+b_{3/2}\sigma _3+c_{3/2}\sigma _1)d\varphi ,`$ (27)
$`a_{3/2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(3\beta _{11}^23\beta _{12}^2+\beta _{21}^2\mathrm{cos}\alpha \beta _{22}^2\mathrm{cos}\alpha \right),`$ (28)
$`b_{3/2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(3\beta _{11}^2+3\beta _{12}^2+\beta _{21}^2\mathrm{cos}\alpha +\beta _{22}^2\mathrm{cos}\alpha \right),`$ (29)
$`c_{3/2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}\alpha \beta _{21}\beta _{22},`$ (30)
$`A_{1/2}`$ $`=`$ $`(a_{1/2}+b_{1/2}\sigma _3+c_{1/2}\sigma _1)d\varphi ,`$ (31)
$`a_{1/2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(3\beta _{21}^23\beta _{22}^2+\beta _{11}^2\mathrm{cos}\alpha \beta _{12}^2\mathrm{cos}\alpha \right),`$ (32)
$`b_{1/2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(3\beta _{21}^2+3\beta _{22}^2+\beta _{11}^2\mathrm{cos}\alpha +\beta _{12}^2\mathrm{cos}\alpha \right),`$ (33)
$`c_{1/2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}\alpha \beta _{11}\beta _{12},`$ (34)
where $`d\varphi =\omega _1dt`$. Note that as $`A`$ does not depend on time, the final solution does not require $`𝒯`$-ordering. It should be also emphasized here that in the non-adiabatic case we discuss the term $`A_{3/2}`$ contains non-diagonal terms that is not the case when the adiabaticity condition is held . Now the solution of the problem takes a particular form of (4):

$$\psi (t)=e^{i\varphi _n(t)}e^{i\omega _1tA}\psi (0),$$
(35)
Formula (35) solves the problem of the evolution control for the system under consideration. The resulting quantum gate is entirely determined by $`A`$ and the evolution law of the magnetic field, i.e. by a contour in the parameter space. Of course it is always possible to choose the parameters so that $`A`$ turns out to generate a 2-qubit transformation that produces a superposition of basis states. For this reason the gate can be thought of as a universal one . Of course, a suitable speed of the parameters evolution can not be reached by rotation of the sample as it took place in the experiment by authors of . Nevertheless it is clear that this manner of control is not principle and one could imagine a situation where the parameters evolution is provided by the controlling magnetic field by adding a non-stationary transverse component.


text:
# The s–wave implosion and explosion of free particles

## Abstract

A free quantum particle in two dimensions with vanishing angular momentum (s–wave) in the form of a ring–shaped wave packet feels an attractive force towards the center of the ring, leading first to an implosion followed by an explosion.

PACS numbers: 03.65.-w,67.40.Vs,03.40.Gc

The spreading of a wave packet is an important topic in quantum mechanics . It arises due to different momenta contained in the packet. As pointed out by Max Born  even a classical ensemble of particles spreads. The only quantum effect lies in the connection between the position and the momentum distributions, that is, the spreads in position and momentum are not independent. In the present paper we predict a wave effect of the free particle that manifests itself in the dynamics: A ring–shaped wave packet in two dimensions of vanishing angular momentum feels an attractive force towards the center of the ring. This effect is unique to two dimensions and manifests itself most strongly for vanishing angular momentum.

In order to illustrate this surprising phenomenon we first consider the free motion of wave packets in $`N`$ dimensions. We then concentrate on two and three dimensions, where we propagate a circular shell and a spherical shell, respectively. All wave packets have initially vanishing angular momentum and vanishing average radial momentum. We show that the circular shell first implodes and then explodes. In contrast, the spherical shell does not implode but only explodes. We demonstrate this wave effect unique to two–dimensional space using (i) analytical expressions for the average radial position and momentum, (ii) the Wigner function approach  towards quantum mechanics, (iii) analytical expressions  for the radial Green’s function in two dimensions and (iv) the numerical propagation of wave packets.

Our analytical approach considers the propagation of the initial wave function

$$\mathrm{\Phi }_0^{(N)}(x_1,\mathrm{},x_N)=\frac{𝒩}{\sqrt{S^{(N)}}}r^\gamma \mathrm{exp}\left[\frac{1}{2}\left(\frac{r}{\delta r}\right)^2\right]$$
(1)
in the absence of any potential. Here $`𝒩\left[2/\mathrm{\Gamma }\left(N/2+\gamma \right)\right]^{1/2}\delta r^{N/2\gamma }`$ denotes the radial normalization constant and $`r\left(x_1^2+x_2^2+\mathrm{}x_N^2\right)^{1/2}`$ is the radial variable in $`N`$ space dimensions. Moreover, $`S^{(N)}2\pi ^{N/2}/\mathrm{\Gamma }\left(N/2\right)`$ is the total solid angle. The quantity $`\delta r`$ defines a characteristic length scale of the packet and is a measure of the radial width .

We first take the exponent $`\gamma `$ to be equal to 2. This special form of the initial wave function – polynomial in coordinates multiplied by a Gaussian – allows us to calculate analytically the time dependence of the average radial position $`r^{(N)}`$ and the average radial momentum $`p^{(N)}`$ in $`N`$ dimensions. We then turn to the general case of arbitrary $`\gamma >1`$ but concentrate on two dimensions.

In the derivation of our results we follow three independent approaches that provide identical results: (i) We consider the time evolution in Wigner phase space and calculate the averages with the help of the Wigner function; (ii) we use the method of the generating function  to propagate a modulated Gaussian wave packet of the free particle in $`N`$ dimensions, and evaluate the averages; (iii) we propagate the free wave packet with the help of the radial Green’s function of the free particle in two dimensions. Here we do not present the detailed calculations but only summarize the results .

We start by presenting first the general results for the propagation in $`N`$ dimensions and then focus on the special case of two and three dimensions. It is convenient to introduce the dimensionless time variable $`\tau \mathrm{}t/(\delta r^2M)`$, where $`M`$ denotes the mass of the particle. We find  the expression

$$r^{(N)}(\tau )=\frac{1+\tau ^2/a^{(N)}}{\sqrt{1+\tau ^2}}r_0^{(N)}$$
(2)
for the average radial position. Here we have introduced the abbreviation

$$a^{(N)}1+\frac{4N}{N^2+3}$$
(3)
and the initial radial position

$$r_0^{(N)}\frac{\mathrm{\Gamma }\left(\frac{N+5}{2}\right)}{\mathrm{\Gamma }\left(\frac{N+4}{2}\right)}\delta r,$$
(4)
in $`N`$ dimensions, where $`\mathrm{\Gamma }`$ denotes the Euler gamma function.

Equation (2) shows that for large times the average radial position increases linearly with time, corresponding to an explosion of the wave packet. This holds true for any number of dimensions.

For short times we can expand the square root in Eq. (2) which yields

$$r^{(N)}(\tau )\left[1+\left(\frac{1}{a^{(N)}}\frac{1}{2}\right)\tau ^2\right]r_0^{(N)}.$$
(5)
An interesting situation occurs when $`a^{(N)}>2`$. In this case, the average radius initially decreases, corresponding to an implosion of the wave packet. From the definition Eq. (3) of $`a^{(N)}`$ we arrive at the implosion condition

$$(N1)(N3)<0.$$
(6)
Hence an implosion occurs only for $`N=2`$ corresponding to $`a^{(2)}=15/7=2+1/7`$.

The phenomenon of the implosion can also be viewed by considering the time evolution of the average radial momentum

$$p^{(N)}(\tau )=M\frac{d}{dt}r^{(N)}=\frac{d}{d\tau }r^{(N)}\frac{\mathrm{}}{\delta r^2}.$$
(7)
determined by the time derivative of the average position $`r^{(N)}`$, which yields

$$p^{(N)}(\tau )=\frac{\tau }{(1+\tau ^2)^{3/2}}\left(a^{(N)}2\tau ^2\right)p_{\mathrm{}}^{(N)}$$
(8)
where

$$p_{\mathrm{}}^{(N)}\underset{\tau \mathrm{}}{lim}p^{(N)}(\tau )=\frac{1}{a^{(N)}}\frac{\mathrm{}r_0^{(N)}}{\delta r^2}.$$
(9)
The radial momentum is negative, corresponding to an implosion, provided $`\tau <\left[a^{(N)}2\right]^{1/2}`$. Since in two dimensions $`a^{(2)}=15/7`$ this condition reads $`\tau <\tau _{\mathrm{min}}1/\sqrt{7}`$. At that moment the momentum vanishes and the implosion turns into an explosion. The average radial position assumes a minimum $`r_{\mathrm{min}}r^{(2)}\left(\tau _{\mathrm{min}}\right)`$.

In order to highlight the special property of two dimensions we now compare the dynamics of the wave packet Eq. (1) for $`\gamma =2`$ in two and three dimensions. In Fig. 1 we show the time evolution of the average radial position and the average radial momentum in two (solid curves) and three (dashed curves) dimensions based on Eqs. (2) and (8) with $`a^{(2)}=15/7`$ and $`a^{(3)}=2`$.

For short times the average radius of a ring–shaped wave packet decreases as a function of time. This initial implosion of the circular wave is followed by a final explosion. We emphasize that this phenomenon is not due to an external classical force but occurs for a free particle moving in two dimensions . Obviously the spreading of the ring–shaped wave packet in the radial direction is asymmetric; initially the probability flux towards the center is larger than towards the outside.

In the three–dimensional case, where we propagate a spherical shell, the average position and the average momentum do not display the soft implosion effect of two dimensions. Indeed, since $`a^{(3)}=2`$, the average momentum is always positive and the average position cannot decrease but must increase.

The contraction effect is maximal at $`\tau _{\mathrm{min}}=1/\sqrt{7}`$ and according to Eq. (2) the average radial position $`r^{(2)}`$ assumes the minimal value

$$r_{\mathrm{min}}\sqrt{\frac{224}{225}}r_0^{(2)}0.9978r_0^{(2)}.$$
(10)
Hence, it is a rather soft implosion.

We note that this result is independent of $`\delta r`$, since there is only one length scale in the problem. When we choose a different form for the packet we can make the effect larger. For example, by replacing $`r^\gamma `$ in Eq. (1) by $`\mathrm{sin}(r^2/\delta r^2)`$, we obtain at $`\tau _{\mathrm{min}}1.11`$

$$r_{\mathrm{min}}0.9964r_0^{(2)}.$$
(11)
In order to study the implosion effect we concentrate on two dimensions and now consider an arbitrary value of the exponent $`\gamma >1`$ in the the wave packet Eq. (1). With the help of the corresponding radial propagator  we obtain the radial wave function which yields the expression

$$r^{(2)}(\tau )=2^\gamma \frac{\mathrm{\Gamma }^2\left(\frac{\gamma +1}{2}\right)}{\mathrm{\Gamma }\left(\gamma +\frac{1}{2}\right)}\frac{\tau ^{\gamma 1}}{\left(1+\tau ^2\right)^{\frac{\gamma +1}{2}}}I(\tau ;\gamma )r_0^{(2)}$$
(12)
for the average radial position. Here the integral

$`I(\tau ;\gamma )`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑\xi \xi ^2\mathrm{exp}\left({\displaystyle \frac{1}{1+\tau ^2}}\xi ^2\right)`$ (14)
$`\left|{}_{1}{}^{}F_{1}^{}[{\displaystyle \frac{\gamma 1}{2}};1;{\displaystyle \frac{i}{2\tau }}{\displaystyle \frac{1}{1+i\tau }}\left({\displaystyle \frac{r}{\delta r}}\right)^2]\right|^2`$
containing the confluent hypergeometric function $`{}_{1}{}^{}F_{1}^{}`$ can be expressed in terms of the generalized hypergeometric function $`{}_{3}{}^{}F_{2}^{}`$. However, we found it more convenient to evaluate this integral numerically. For each value of $`\gamma `$ we obtain an implosion but the ratio $`r_{\mathrm{min}}/r_0^{(2)}`$ is always of the same order of magnitude.

We now turn to the explanation of the implosion effect using the Wigner function . The free particle is one of the few quantum mechanical systems in which the classical phase space dynamics is identical to the quantum mechanical one. Indeed, the Wigner function satisfies the classical Liouville equation. Hence, we find the Wigner function at time $`t`$ from the Wigner function at time $`t=0`$ by replacing the position $`\stackrel{}{r}`$ by $`\stackrel{}{r}\stackrel{}{p}t/M`$. We emphasize that now $`\stackrel{}{r}`$ and $`\stackrel{}{p}`$ denote $`N`$–dimensional vectors. The corresponding probability distributions in position and momentum follow by integration over the conjugate variable.

Since for the free particle classical and quantum phase space dynamics are identical the quantum effect of implosion rather than explosion is stored in the initial Wigner function. To be more precise, it is in the correlations between position and momentum. A classical phase space distribution enjoys a probability interpretation and therefore always has to be positive. In contrast, the Wigner function can assume negative values. Indeed, both the Wigner functions of the circular shell and of the spherical shell contain domains in phase space where they become negative. Since the Wigner function is normalized the total volume of the Wigner function is unity. We can calculate numerically the volumes $`V_{}^{(2)}`$ and $`V_{}^{(3)}`$ corresponding to the negative parts. We find in the two dimensional case $`V_{}^{(2)}0.27`$, whereas in three dimension $`V_{}^{(3)}0.23`$. Hence, there is a slightly larger contribution of negative parts in two dimensions than in three dimensions.

In the phase space description of quantum mechanics quantum effects are stored in the negative parts of the initial Wigner function. Since a classical ensemble of particles cannot have a negative phase space distribution it cannot display the implosion effect. Needless to say a classical ensemble with a ring–shaped initial position distribution can show an implosion provided there exists a non–vanishing negative initial radial momentum. However, in our quantum mechanical wave packet the average radial momentum vanishes.

For the choice Eq. (1) of the initial wave packet the radial probability density

$`W^{(N)}(r)dr`$ $`=`$ $`\mathrm{\Phi }_0^{(N)}(r)^2r^{N1}drS^{(N)}`$ (15)
$`=`$ $`𝒩^2r^{N+2\gamma 1}e^{(r/\delta r)^2}dr`$ (16)
depends on the number of dimensions. Hence, one might argue that the implosion is a result of the slightly different initial radial wave functions. In order to exclude this argument we have performed a numerical integration of the Schrödinger equation

$$i\mathrm{}\frac{}{t}u_0^{(N)}(r,t)=\left[\frac{\mathrm{}^2}{2M}\frac{^2}{r^2}+V^{(N)}(r)\right]u_0^{(N)}(r,t)$$
(17)
for the radial wave function $`u_0^{(N)}`$ of the free particle in two and three space dimensions, corresponding to zero angular momentum. Here we always start with the same initial radial wave function. In this way we guarantee identical radial probabilities densities in two and three dimensions.

The radial potential

$$V^{(N)}(r)\frac{\mathrm{}^2}{2M}\frac{(N1)(N3)}{4r^2}$$
(18)
in $`N`$ dimensions involves the product $`(N1)(N3)`$ familiar from the implosion condition Eq. (6).

In three dimensions the potential $`V^{(3)}`$ vanishes. However, for $`N=2`$ the quantum anti–centrifugal potential

$$V_Q(r)V^{(2)}(r)\frac{\mathrm{}^2}{2M}\frac{1}{4r^2}$$
(19)
does not vanish and is negative, that is, attractive . This potential is responsible for the implosion.

A large class of wave packets with vanishing angular momentum displays this effect of first imploding and then exploding . Indeed, in the insert of Fig. 1 we show the short time evolution of the average radial position and the average radial momentum in two and three dimensions for an initial radial Gaussian wave packet $`u(r,t=0)\stackrel{~}{𝒩}\mathrm{exp}\left[\left(r\rho \right)^2/\delta r^2\right]`$ displaced by an amount $`\rho `$ from the origin. As in the analytical results we find an implosion of a circular but not of a spherical shell.

So far we have concentrated on a single particle but the effect of implosion manifests itself also in the relative motion of two free particles. We can envision the total wave function of the two particles to be of the form Eq.(1) corresponding to an entangled state. When we assume that each particle is allowed to move in a $`D`$–dimensional space the total wave function lives in a $`N=2D`$ dimensional configuration space. The physical situation of a single particle imploding in two dimensions corresponds to two particles constrained to one dimension. When they move along a line their average relative separation first decreases and then increases. It looks like the particles first attract and then repel each other. However, when they move in two or three dimensions, their relative distance always increases.

We conclude by briefly outlining an experiment to demonstrate the implosion effect. At first sight, one might think of cold atoms around a wire familiar from atom chips . However, it might be easier to make use of the analogy between the Schrödinger equation and the paraxial wave equation of classical electrodynamics. We shine monochromatic light onto an opaque screen with an annular aperture. In order to avoid diffraction from sharp edges, we introduce apodization, that is, a position dependence transmission coefficient. At the edges of the ring almost no light is transmitted, whereas at the center of the ring almost all light passes through. This arrangement creates an intensity distribution after the screen confined mainly to a ring. The subsequent free propagation in space is analogous to the free propagation in time of the Schrödinger wave function. The time dependent probability distribution then corresponds to the transverse intensity distributions at different locations downstream from the screen .

We thank I. Sh. Averbukh, M. V. Berry, D. Bouwmeester, H. Carmichael, D. Greenberger, D. Kobe, R. F. O’Connell, W. C. Schieve, F. Straub, S. Varro and K. Wódkiewicz for many fruitful discussions. I. B.–B., J. P. D. and M. F. gratefully acknowledge the support of the Humboldt foundation and the great hospitality enjoyed at the Abteilung für Quantenphysik. The work of W. P. S. is partially supported by DFG.


text:
# Numerical analysis of the master equation

## Abstract

Applied to the master equation, the usual numerical integration methods, such as Runge-Kutta, become inefficient when the rates associated with various transitions differ by several orders of magnitude. We introduce an integration scheme that remains stable with much larger time increments than can be used in standard methods. When only the stationary distribution is required, a direct iteration method is even more rapid; this method may be extended to construct the quasi-stationary distribution of a process with an absorbing state. Applications to birth-and-death processes reveal gains in efficiency of two or more orders of magnitude.

KEYWORDS: Markov process; master equation; numerical integration

The master equation is the basic tool for describing Markovian stochastic processes on a discrete state space, in continuous time . Despite its central role in stochastic analysis in the physical and biological sciences, only a handful of exactly soluble examples are known. Thus the need to integrate the master equation numerically arises frequently, and with it the issue of computational efficiency.

Consider a Markov process with transition rates $`W_{n,m}`$ from state $`m`$ to state $`n`$. In the master equation (ME),

$$\dot{p}_n=\underset{n^{}}{}W_{n,n^{}}p_n^{}p_n\underset{n^{}}{}W_{n^{},n}(n=0,\mathrm{},N),$$
(1)
the factor $`_n^{}W_{n^{},n}w_n`$ multiplying $`p_n`$ can become large in a typical birth-and-death process, where $`W_{n^{},n}n`$ or some higher power of $`n`$. This poses a problem for numerical integration via the usual discretization schemes, such as the Runge-Kutta method (RKM). Since instabilities appear when $`|w_n\mathrm{\Delta }t|1`$, we must use a small time increment, $`\mathrm{\Delta }t1/(\mathrm{max}_nw_n)`$, and convergence is slow. With some (but not all) of the transition rates large, the ME is, in effect, a stiff system of differential equations, requiring special numerical treatment . In this note I introduce numerical integration schemes for the ME that are efficient when the transition rates vary over a wide range, and an iteration method that eliminates the need for step-by-step integration, when only the stationary (or quasi-stationary) distribution is required. Simple but detailed examples are used to illustrate the methods; further applications will be reported elsewhere . (We consider stationary Markov processes, i.e., time-independent rates, although the method is not limited to this class of problem.)

We begin by writing the ME in the form

$$\dot{p}_n=w_np_n+r_n,$$
(2)
where

$$r_n(t)=\underset{n^{}}{}W_{n,n^{}}p_n^{}(t).$$
(3)
Note that $`p_n`$ does not appear in the sum for $`r_n`$, since $`W_{n,n}0`$. Both $`w_n`$ and $`r_n`$ are nonnegative. (In fact $`w_n`$ is zero only if state $`n`$ is absorbing.) Integrating Eq. (2) we have

$$p_n(t)=e^{w_nt}p_n(0)+_0^t𝑑t^{}e^{w_n(tt^{})}r_n(t^{}).$$
(4)
This is only a formal solution, since we need the $`p_n(t)`$ to evaluate $`r_n`$, but it is a useful starting point for approximations. If we adopt a time increment $`\mathrm{\Delta }t`$ such that $`r_n(t)`$ is approximately constant over this interval, then we have

$$p_n(\mathrm{\Delta }t)e^{w_n\mathrm{\Delta }t}p_n(0)+\left[1e^{w_n\mathrm{\Delta }t}\right]\frac{r_n(0)}{w_n}.$$
(5)
This relation can then be iterated: we use $`p_n(\mathrm{\Delta }t)`$ to evaluate $`r_n`$ at the start of the second interval, and thereby find $`p_n(2\mathrm{\Delta }t)`$, and so on. This simple formal integration (FI) scheme, analogous to Euler’s method for direct numerical integration, already represents a significant advantage over the usual approaches when some of the $`w_n`$ are large. The reason is that the exponential factor is already included in the solution, whereas in the usual discretizion methods it has to be constructed term by term, in powers of $`w_n\mathrm{\Delta }t`$.

Suppose that the ME of interest possesses a unique, stable stationary solution $`\overline{p}_n`$. For a stationary solution, $`\overline{p}_n=\overline{r}_n/w_n`$, where $`\overline{r}_n`$ is given by Eq. (3) with $`p_n=\overline{p}_n`$. Setting $`p_n=\overline{p}_n`$ on the r.h.s. of Eq. (5) we immediately see that $`\overline{p}_n`$ is a stationary solution of our iteration method. To investigate the linear stability of the stationary solution, let $`p_n=\overline{p}_n+\delta _n`$. According to the ME,

$$\dot{\delta }_n=\underset{n^{}}{}W_{n,n^{}}\delta _n^{}\delta _n\underset{n^{}}{}W_{n^{},n},$$
(6)
or, in matrix notation,

$$\frac{d}{dt}\delta =S\delta $$
(7)
where $`S_{n,m}=W_{n,m}`$ for $`nm`$ and $`S_{n,n}=w_n`$. Since, by hypothesis, the stationary solution $`\delta =0`$ is stable, $`S`$ has one zero eigenvalue and all others negative.

If we insert $`p_n=\overline{p}_n+\delta _n`$ in the r.h.s. of Eq. (5) (so that $`r_n=\overline{r}_n+_n^{}W_{n,n^{}}\delta _n^{}`$), we find

$$\delta (t+\mathrm{\Delta }t)=T\delta (t)$$
(8)
where $`T_{n,m}=(1e^{w_n\mathrm{\Delta }t})W_{n,m}/w_n`$ for $`nm`$, and $`T_{n,n}=e^{w_n\mathrm{\Delta }t}`$. Expanding to first order in $`\mathrm{\Delta }t`$, we see that Eq. (7) applies in this case as well, so the stationary solution is stable if it is so for the original ME.

We refer to the integration method embodied in relation Eq. (5) as a “first-order FI scheme” since the error in the solution $`p_n(t)`$ is proportional (for fixed $`t`$) to $`\mathrm{\Delta }t`$. \[In each step, the error incurred in treating $`r_n`$ as constant is $`𝒪(\mathrm{\Delta }t)^2`$ (since in fact $`r_n(\mathrm{\Delta }t)=r_n(0)+r_n^{}(0)\mathrm{\Delta }t+\mathrm{}`$), and we require $`N=t/\mathrm{\Delta }t`$ steps.\] An obvious way to improve accuracy, analogous to going from Euler’s method to the midpoint (or second-order Runge-Kutta) method, is to replace the assumption of a constant $`r_n`$ on the interval $`[t,t+\mathrm{\Delta }t]`$ with a linear approximation. To do this we first use Eq. (5) to find $`p_n(t+\mathrm{\Delta }t)`$, and then estimate the $`r_n(t+\mathrm{\Delta }t)`$ using these values. Then we form the linear approximation

$$r_n(s)r_n(t)+(st)r_n^{},(tst+\mathrm{\Delta }t),$$
(9)
where

$$r_n^{}=\frac{r_n(t+\mathrm{\Delta }t)r_n(t)}{\mathrm{\Delta }t}$$
(10)
Inserting this in the formal solution, Eq. (4), we find

$$p_n(t+\mathrm{\Delta }t)e^{w_n\mathrm{\Delta }t}p_n(t)+\left[1e^{w_n\mathrm{\Delta }t}\right]\frac{r_n(t)}{w_n}+\left[\mathrm{\Delta }t\frac{1e^{w_n\mathrm{\Delta }t}}{w_n}\right]\frac{r_n^{}}{w_n}$$
(11)
(Note that we do not use the expression for $`dr_n/dt`$ that follows from the master equation in place of $`r_n^{}`$, since it does not, in general, represent the behavior of $`r_n`$ over a substantial time interval.) The total error associated with the second-order FI procedure is $`(\mathrm{\Delta }t)^2`$. Further improvement is clearly possible, by introducing higher order polynomial approximations to the $`r_n`$, but the second-order scheme is quite adequate for the applications considered here. The error estimates apply for a fixed, finite time; the error for $`t\mathrm{}`$ is zero, if the system possesses a unique stationary distribution.

A tremendous simplification and speedup is possible if one only requires the stationary probability distribution. The relation $`\overline{p}_n=\overline{r}_n/w_n`$ suggests an iterative procedure of the form

$$p_n^{}=ap_n+(1a)\frac{r_n}{w_n}$$
(12)
where $`0<a<1`$ is a parameter. We expect $`p_n^{}`$ to converge to $`A\overline{p}_n`$, where the factor $`A`$ depends on $`a`$ and on the initial distribution.

Some insight into the choice of $`a`$ is afforded by the simple example of a two-state system with transition rates $`W_{1,0}=\lambda `$ and $`W_{0,1}=\mu `$. The evolution of the probability distribution is given (in matrix notation) by $`p^{}=Mp`$, where

$$M=\left(\begin{array}{cc}a& (1a)\frac{\mu }{\lambda }\\ & \\ (1a)\frac{\lambda }{\mu }& a\end{array}\right)$$
(13)
Matrix $`M`$ has eigenvalues 1 and $`2a1`$. Thus the iterative procedure converges to a multiple of the stationary distribution, $`(\mu ,\lambda )`$ for $`0<a<1`$, and is instantaneous for $`a=1/2`$.

A similar analysis of the three-state processes ($`n=0,1,2`$), with rates

a) $`W_{n+1,n}=\gamma `$, $`W_{n1,n}=1`$,

and

b) $`W_{n+1,n}=\gamma `$, $`W_{2,0}=\gamma ^2`$, and $`W_{n,m}=1`$ for $`n<m`$,

yields $`a=1/3`$ as the optimal choice. (By optimal choice we mean the value that minimizes $`\mathrm{max}_i|\omega _i|`$, the $`\omega _i`$ being the eigenvalues of $`M`$, excluding, of course $`\omega _01`$, associated with the stationary distribution.) These results suggest that for processes with a large number of states, even smaller values of $`a`$ may be advantageous. The numerical examples discussed below support this notion.

A simple modification of the iterative scheme, Eq. (12), generates the quasi-stationary distribution (if such exists) of a Markov processs with an absorbing state. Consider a process on the states $`n=0,1,2,\mathrm{}N`$ with $`n=0`$ absorbing, i.e., $`W_{n,0}0`$ for all $`n`$, while $`W_{0,n}`$ for at least one $`n>0`$. (It seems reasonable to exclude manifestly transient states, in other words, we assume that for each $`n`$, $`W_{n,m}>0`$ for at least one $`m`$.) In this case the stationary state is $`\overline{p}_n=\delta _{n,0}`$, but it is possible that the probability distribution, conditioned on survival, attains a stationary form, that is, for long times $`p_n(t)C(t)q_n`$ ($`n>0`$), where the $`q_n`$ are time-independent. Such quasi-stationary distributions arise in birth-and-death processes with saturation, for example the Malthus-Verhulst process or the contact process .

The defining feature of the quasi-stationary distribution is that relative probabilities $`p_n(t)/p_m(t)`$ ($`n`$, $`m>0`$) are constant, or equivalently, $`\dot{p}_n/p_n=\kappa `$, constant and independent of $`n`$ for $`n1`$. Suppose the probability distribution has attained a quasi-stationary form at some time $`t`$, and that at this instant the distribution is normalized so: $`_{n1}p_n(t)=1`$. Using Eq. (2) we have $`r_n=(\kappa +w_n)p_n`$, and summing on $`n1`$ yields

$$\kappa =\underset{n1}{}\dot{p}_n=\underset{n1}{}W_{0,n}p_nr_0.$$
(14)
($`r_0`$ is the decay rate of the survival probability in the quasi-stationary regime.) Thus we find that in the quasi-stationary state,

$$p_n=\frac{r_n}{w_nr_0},(n1).$$
(15)
This relation suggests that we iterate as follows:

$$p_n^{}=ap_n+(1a)\frac{r_n}{w_nr_0}.$$
(16)
The new distribution $`p^{}`$ should be normalized after each iteration, since Eq. (15) assumes this property. \[If $`r_0=0`$ the process of course possesses a true stationary distribution and Eq. (16) reduces to Eq. (12).\] We have verified that this scheme converges to the quasi-stationary distribution much more rapidly than via integration of the ME.

As a first example we consider the nonlinear one-step process (on states $`n=1,2,3,\mathrm{}`$), with nonzero transition rates:

$$W_{n1,n}=n(n1),W_{n+1,n}=\lambda $$
(17)
This represents the coagulation process $`A+AA`$ with addition of species $`A`$ at rate $`\lambda `$, in a well-stirred system. In Fig. 1 we compare the mean population size $`n_t`$ as furnished by integration of the ME via a fourth-order RKM (solid line) with the second-order FI scheme, Eq. (11), (open circles). \[Here $`\lambda =100`$, $`p_n(0)=\delta _{n,1}`$, and we have set $`p_n(t)0`$ for $`n>100`$.\] The results are nearly identical, but the RKM requires and integration step of $`\mathrm{\Delta }t2\times 10^4`$ for stability, whereas the FI integration uses $`\mathrm{\Delta }t=0.002`$. The latter yields $`n_t`$ with a relative error of less than $`0.8\%`$, the largest deviation occurring at relatively short times. (As noted above, the stationary values are identical in all cases.) To maintain fidelity to the solution at short, as well as long times, one may use the FI with an adjustable time step ($`+`$ in Fig. 1), such that $`r^{}\mathrm{\Delta }t`$ is small, where $`r^{}=\mathrm{max}_n\{r_n^{}\}`$. For example, using the second-order FI with a time step of $`\mathrm{\Delta }t=0.0005+.02/(1+r^{})`$, the relative error is reduced to $`0.06\%`$, while the average time step over the region of interest is 0.0024.

Fig. 2 shows the stationary probability distribution for $`\lambda =400`$. The solid line represents the RKM result, while the points come from iteration of Eq. (12). The distributions are identical to within one part in $`10^5`$. In this case, the distribution converges (such that the error in $`n_t`$ is $`<10^{10}`$), after about 900 iterations, when we use $`a=1/2`$. As $`a`$ is reduced, the number of iterations required falls steadily; in fact the procedure converges even for $`a=0`$, after only 460 steps.

Our second example is a multi-step generalization of the Malthus-Verhulst process, with transition rates:

$$W_{m,n}=\{\begin{array}{cc}n[1+\mu (n1)]e^{\gamma (mn+1)},\hfill & m<n\\ \lambda ne^{\gamma (nm+1)},\hfill & m>n\\ 0,\hfill & m=n\end{array}$$
(18)
for $`n`$, $`m=0,1,2,\mathrm{}`$. (Thus transition rates between states $`n`$ and $`m`$ fall off $`\mathrm{exp}[\gamma |mn|]`$. We choose this example to show that the methods proposed here are not limited to one-step processes, which admit a relatively simple analysis .) In this case $`n=0`$ is absorbing. Figure 3 compares the quasi-stationary distribution for $`\lambda =3`$, $`\mu =0.1`$, and $`\gamma =1`$, obtained via the RKM (with $`\mathrm{\Delta }t=10^4`$; larger values produce instability), with iteration of Eq. (16). As is clear from the plot, the results are identical in every detail. The RKM requires about 1270 seconds of cpu on a DEC-alpha workstation to converge to the quasi-stationary distribution; the FI method requires about 50 sec. Iteration of Eq. (16), (with $`a=0.5`$), by contrast, converges in 9 seconds (about 1600 steps). Reducing parameter $`a`$ to 0.1 yields convergence in just 900 steps (5 sec. cpu), and in 815 steps for $`a=0`$.

The steady decrease in computation time as we reduce $`a`$ leads one to ask whether it is possible to use $`a`$ negative. Iteration of Eq. (16) with $`a<0`$ does not work, as negative probabilities are generated, but we can exclude these by writing

$$p_n^{}=\mathrm{max}[0,ap_n+(1a)\frac{r_n}{w_nr_0}],$$
(19)
This functions even for negative $`a`$, converging, in the present case, after about 630 steps when $`a=0.3`$. For $`a0.4`$, however, the scheme does not converge. On the other hand, in the first example \[i.e., rates given by Eq. (17)\], using $`a<0`$ does not offer any advantage over $`a=0`$.

We conclude from these and other examples with large numbers of states , that $`a0`$ generally yields rapid convergence, but that some amount of experimenting (for example, with negative $`a`$ values), may reduce computation time further, in particular cases. But even without such optimization, the iterative schemes of Eqs. (12) and (16) yield economies in computation time of two or more orders of magnitude, compared with conventional integration methods. As the number of variables and/or the disparity in the magnitudes of the transition rates increases, direct iteration and the FI scheme become ever more valuable.

Acknowledgment

I thank Dani ben-Avraham for helpful comments. This work was supported by CNPq, Brazil.

FIGURE CAPTIONS

FIG. 1. Mean population $`n`$ versus time in the coagulation process, Eq. (17), with $`\lambda =100`$, using Runge-Kutta integration (solid line), the second-order FI scheme with a fixed time step (circles), and with a variable time step, as described in the text ($`+`$).

FIG. 2. Stationary probability distribution in the coagulation process with $`\lambda =400`$. Solid line: RKS; $`+`$: direct iteration of Eq. (12).

FIG. 3. Quasi-stationary probability distribution in the generalized Malthus-Verhulst process, Eq. (18), with $`\lambda =3`$, $`\mu =0.1`$, and $`\gamma =1`$. Solid line: RKS; $`+`$: direct iteration of Eq. (16).