text:
# First passage time exponent for higher-order random walks: Using Levy flights

## I Introduction

We investigate the general random walk $`x(t)`$ obeying the equation of motion

$$\frac{d^nx(t)}{dt^n}=\eta (t),$$
(1)
where $`\eta `$ is white noise with zero mean and unit variance and $`x(0)=x_o`$.

Let us begin with $`n=1`$. The *first passage time* is the time it takes for the walk to reach zero. When the walk has no bias, as above, there is no definite time to expect such an event and the distribution is a power law. To find the first passage time distribution $`\rho (t)dt`$, start an ensemble of random walkers at $`x_o>0`$ and at time $`t=0`$. Whenever a random walk reaches zero, it is removed from the ensemble. Let $`P_t(x)`$ be the number density of walks at time $`t`$ and positive $`x`$. $`P_t(x)`$ is a solution of the diffusion equation with absorbing boundary condition, $`P_t(0)=0`$. More precisely,

$`P_t(x)={\displaystyle \frac{1}{\sqrt{2\pi t}}}(e^{\frac{(xx_o)^2}{2t}}e^{\frac{(x+x_o)^2}{2t}})`$ (2)
$`2x_o{\displaystyle \frac{d}{dx}}({\displaystyle \frac{1}{\sqrt{2\pi t}}}e^{\frac{x^2}{2t}})`$ (3)
at long times.

Let $`g(t)`$ be the integral of this probability distribution over positive values of $`x`$. This is the probability that a random walk at time $`t`$ has not crossed zero. The first passage time distribution is then given by

$$\rho (t)dt=\frac{dg}{dt}dt=\frac{x_o}{\sqrt{2\pi t}}\frac{dt}{t}.$$
(4)
Although we will not show it, this result is universal for all symmetric walks. Any random walk which is equally likely to move forward as backward by a given amount has a first passage time distribution with the same asymptotic form. This is the Sparre-Anderson theorem . Note that $`x_o^2`$ sets the time scale.

The previous discussion is one of the few first passage time problems where an exact solution may be easily found by solving a Fokker-Planck equation with absorbing boundaries. Extensions of this method to more complicated random walks, such as the second-order random walk described by

$$\frac{d^2x}{dt^2}=\eta (t),$$
(5)
exist . However, we investigate the first passage time distribution for this walk, and for any $`n`$, in a different way.

At long times, the first passage time distribution for these processes is a power law,

$$\rho _n(t)dt\frac{1}{t^{\beta _n}}\frac{dt}{t}.$$
(6)
As computed above, $`\beta _1=\frac{1}{2}`$. $`\beta _2`$ is known to be $`\frac{1}{4}`$ . The others are most likely not exact fractions and are nontrivial critical exponents of certain statistical models. We first present numerical results for $`n=2,3,4,5`$. Then, after presenting a heuristic derivation of $`\beta _2`$, we make a quantitative estimate for the shift in first passage exponent for $`n3`$. Our analytical results will draw heavily from the theory of Levy flights . We also address separately the $`n=\mathrm{}`$ limit .

Finally, within zero temperature mean field theory, we demonstrate that the first passage time exponent in the $`n=2`$ case is the avalanche size exponent for a dynamically elastic extended object, like a crack front or interface, being pulled through a random medium by a uniform applied force at a special point in the parameter space . Generically, for these nonequilibrium systems there is a transition from an overall stationary phase to an overall moving phase as the applied force is increased toward a critical value. Increasing the applied force from the static side triggers *local* motion of the interface for some finite amount of time as long as the applied force is held fixed from the time the toppling starts until the time the toppling stops. The distribution of amount the interface has moved during an “avalanche” event, i.e. the avalanche size, gives us information about the continuity (or discontinuity) of the depinning transition and, therefore, has been a focus of study over the years .

## II Levy Flights

A *Levy variable* $`q_i`$ is a random variable with a power law distribution for large $`q`$,

$$P(q)dq\frac{1}{q^\beta }\frac{dq}{q},$$
(7)
with $`0<\beta <2`$. The variance $`<q^2>`$ is infinite for these distributions and $`\beta `$ is called the Levy exponent of $`q`$.

A *Levy flight* is a random walk with each step a Levy variable. It is the sum of many independently drawn Levy variables. Because of the infinite variance, a Levy flight is an irregular walk dominated by the few largest jumps .

Let $`L_N(Q)`$ be the distribution of $`Q=_{i=1}^Nq_i`$. The Fourier transform $`\stackrel{~}{L}_N(k)`$ is the $`N`$th power of the Fourier transform of $`P(q)`$. If $`\stackrel{~}{P}(k)`$ were twice differentiable at zero, then $`\stackrel{~}{L}_N(k)(1bk^2+𝒪(k^3))^Ne^{bNk^2}`$ for $`b`$ real and positive. This is the central limit theorem, and it applies when the second moment of a distribution is finite. For a Levy variable, the second moment of the distribution is infinite, and $`\stackrel{~}{P}(k)`$ has a cusp at zero for $`0<\beta <2`$. The form of the cusp may be determined as follows: for small $`k`$, the Fourier transform is approximately the integral of the distribution over the first wavelength $`\lambda =\frac{2\pi }{k}`$, or

$$_1^\lambda \frac{C}{q^\beta }\frac{dq}{q}=1C^{}\frac{1}{\lambda ^\beta }=1C^{\prime \prime }k^\beta .$$
(8)
When there is a cusp, the limiting distribution is not a Gaussian, but has the following form

$$\stackrel{~}{L}_N(k)=e^{bN|k|^\beta }k>0$$
(9)
$`\stackrel{~}{L}_N(k)=e^{b^{}N|k|^\beta }k<0.`$

By rescaling $`k`$, $`b`$ can be made into a pure phase. For the case where $`P(q)`$ is symmetric about zero, the Fourier transform is real and $`b=1`$. It is clear that this distribution has the correct form near $`k=0`$, and increasing $`N`$ is equivalent to a rescaling of $`k`$. The limiting distribution gets wider without changing shape — so it is a fixed point of convolution.

We will need one special nonsymmetric distribution, the distribution of a flight composed of only positive steps. In this case, $`L_N(Q)`$ is zero for $`Q`$ negative, therefore $`\stackrel{~}{L}_N(k)`$ must be analytic in the lower half-plane. In other words, $`\stackrel{~}{L}(k)`$ for $`k<0`$ is the analytic continuation of $`\stackrel{~}{L}(k)`$ for $`k>0`$. As $`k`$ is rotated to $`k`$ through negative imaginary values, $`k^\beta `$ acquires a phase $`e^{\pi i\beta }`$. We conclude that

$$b^{}=be^{i\pi \beta }$$
(10)
$`b=e^{i\pi \beta /2}.`$

For the special case $`\beta =\frac{1}{2}`$, $`b=\frac{1+i}{\sqrt{2}}`$.

And finally, if we were to compute the first passage time distribution for the process,

$$\frac{dx}{dt}=\xi (t),$$
(11)
where $`\xi (t)`$ is Levy noise (a Levy variable symmetrically distributed about zero), then we would find a first passage time exponent of $`\frac{1}{2}`$. Even though the Levy flight is irregular, the position does not cross zero any faster than it does in the nearest-neighbor random walk. This is the Sparre-Anderson theorem once more, which we refer the reader to Ref.  for the details.

## III Numerical Results

To obtain the first passage time exponent $`\beta _n`$ numerically is not as easy as it might appear. Direct numerical integration of the equation of motion becomes more cumbersome with increasing $`n`$. To efficiently simulate the equation, we have calculated the free-space propagation kernels for $`n=2,3,4,5`$ (Appendix A). In this context, the propagation kernel is the probability distribution of $`x(t_o+t_s)`$, given initial values $`x(t_o)`$, $`x^{}(t_o)`$, … , $`x^{(n1)}(t_o)`$ and final values $`x^{}(t_o+t_s)`$, … $`x^{(n1)}(t_o+t_s)`$. We first generate the highest-order time derivative, $`x^{(n1)}(t_o+t_s)`$, and then new values for each lower-order time derivative until $`x(t_o+t_s)`$ is updated.

The next time step $`t_s`$ is chosen so that the variance for the next value of $`x`$ will be smaller. The ratio of the new variance to the old is $`\chi `$. $`\chi `$ should be small so that the walk does not become negative then positive within one time step. The probability of this occurring is exponentially small in the inverse of $`\chi `$. With this algorithm, if $`x(t_o)`$ is large, the time steps are large as well. The number of time steps required for the simulation only grows logarithmically as a function of the first passage time. Near the end of the simulation, when $`x`$ is small, the time steps become small as well, and accuracy is not sacrificed.

The first passage time distributions for $`n=2,3,4,5`$ are shown in Fig. 1 on a log-log scale. Each plot contains approximately $`10^7`$ samples in bins of doubling size. Table I contains the linear regression values of the exponents. The results are independent of $`\chi `$ over the range $`[0.025,0.005]`$, indicating that $`\chi `$ is small enough.

## IV The First Passage Time Exponent for $`n=2`$

For the $`n=2`$ case, $`x(t)`$ is the integral of a random walk. In other words, the variable that executes the random walk is $`x^{}(t)`$, the velocity. To exploit this fact, we introduce a new time counter $`i`$ and divide the time axis into intervals $`\mathrm{\Delta }t_i`$ between the points where the velocity crosses zero. Then the interval sizes $`\mathrm{\Delta }t_i`$ are first passage times for an ordinary random walk.

There is one complication. Referring back to Eq. (3), the time scale until a zero-crossing is the initial value squared. So right after a zero-crossing, the initial value is zero. Since the probability distribution is singular at zero, the random walk reaches zero again instantly, and then again, infinitely many times. This is a well known property of random walks — they jiggle about every value before moving on. Since this affects only the distribution of the smallest times, we cut off the distribution of times near zero by imagining the system on a lattice. Now we have a finite and discrete set of time intervals $`\mathrm{\Delta }t_i`$ with each interval distributed with the Levy exponent $`\beta _1=\frac{1}{2}`$. These time intervals are independently distributed because the velocity undergoes a simple random walk.

To compute the first passage time for the position, we must compute the total time $`t=_{i=1}^N\mathrm{\Delta }t_i`$ until the position becomes negative. $`N`$ is the first zero-crossing of the velocity which happens at a negative value of the position so we are actually slightly overestimating $`t`$ by summing $`N`$ intervals. We can control for this by summing $`N1`$ intervals instead, which would be a slight underestimate. This estimate of $`t`$, be it over or under, should not affect the first passage time exponent since we do not expect the last position step (corresponding to the partial time step) to be arbitrarily large.

To continue we need to know two things — 1) the distribution of $`N`$, $`f(N)`$, and 2) given $`N`$, the probability distribution of $`t`$.

The first problem is actually no problem at all. Since the velocity is just as likely to go up at the start of a time interval as it is to go down, the position is as likely to be positive as negative. Therefore, the position is a symmetric walk on the zero-crossings (actually, a symmetric Levy flight with independent steps distributed with Levy exponent 1/3 using scaling arguments). By the Sparre-Anderson theorem , the distribution of the first-passage-$`N`$ is asymptotically the same as for an unbiased, simple random walk, and

$$f(N)\frac{1}{N^{\frac{3}{2}}}$$
(12)
for large $`N`$.

And now for the second problem. The distribution of the sum of $`N`$ independent $`\mathrm{\Delta }t`$’s, $`L_N(t)dt`$, is a Levy distribution. The fact that the steps are distributed with Levy exponent $`\beta _1=\frac{1}{2}`$ and that all of the $`\mathrm{\Delta }t`$’s are positive fixes the Fourier transform of the Levy distribution. The inverse Fourier transform can be computed exactly in this case. For large t,

$$L_N(t)dt\frac{N}{t^{1/2}}e^{\frac{N^2}{2t}}\frac{dt}{t}.$$
(13)
This distribution is small for $`N^2>t`$, and it may be approximated by $`\frac{N}{t^{\frac{1}{2}}}\frac{dt}{t}`$ for $`N^2<t`$. Using scaling, analogous results may be derived for $`\beta \frac{1}{2}`$, where a closed form does not exist.

Now, the time intervals we are adding are not really independent. This is because we are restricting attention to Levy flights that end on the $`N`$th step and this an atypical sample of all Levy flights with $`N`$ steps. However, the properties we actually use from the previous distribution — namely that it is small for $`N^2>t`$ and the appropriate power law for large $`t`$ — are the scaling laws for the sum of almost any collection of $`N`$ Levy variables, where $`N`$ is large. So, although we do not prove it, we assume that these properties hold for the correct $`L_N(t)dt`$, or the distribution of the time elapsed for those Levy flights that end on the $`N`$th step. This is why we refer to our derivation as heuristic.

All the ingredients are now in place for computing $`\beta _2`$. There is a probability $`f(N)`$ for any value of $`N`$, and given $`N`$, we know the probability distribution of $`t`$. So to find the total distribution of $`t`$, we sum up the conditional distributions $`L_N(t)`$ weighted by the probability of $`N`$, or

$$\rho _2(t)dt=\underset{N=0}{\overset{\mathrm{}}{}}f(N)L_N(t)dt.$$
(14)
Approximating the sum with an integral and using the approximate scaling form for $`L_N(t)`$,

$$\rho _2(t)dt_0^{\mathrm{}}𝑑Nf(N)L_N(t)𝑑t_0^\sqrt{t}\frac{1}{N^{\frac{3}{2}}}\frac{N}{t^{\frac{3}{2}}}𝑑N.$$
(15)
We obtain the following asymptotic form,

$$\rho _2(t)dt=\frac{1}{t^{\frac{1}{4}}}\frac{dt}{t},$$
(16)
so that $`\beta _2=\frac{1}{4}`$. This result is unchanged when $`N`$ is shifted by one unit, so that the last step is of no consequence for the exponent as there are typically many zero-crossings, making our derivation self-consistent.

Y. G. Sinai has rigorously computed $`\beta _2`$ ; we have used some similar methods.

## V The first passage time exponent for $`n3`$

The next random walk we consider is the random surge process, governed by

$$\frac{d^3x}{dt^3}=\frac{da}{dt}=\eta (t).$$
(17)
The acceleration is now an ordinary random walk. After a time t, the acceleration, velocity and position scale as $`a_{\mathrm{typ}}t^{\frac{1}{2}}`$, $`v_{\mathrm{typ}}t^{\frac{3}{2}}`$, and $`x_{\mathrm{typ}}t^{5/2}`$, respectively. We proceed as we did for the $`n=2`$ case. Once again, we divide the time axis by the zero-crossings of the velocity. The time intervals $`\mathrm{\Delta }t_i`$ are now distributed with the Levy exponent $`\beta _2=\frac{1}{4}`$. The sum of $`N`$ independent $`\mathrm{\Delta }t`$’s is approximately zero for $`t<N^4`$, and for larger $`t`$, is approximated by

$$L_N(t)dt\frac{N}{t^{\frac{1}{4}}}\frac{dt}{t}.$$
(18)
Proceeding casually, we may think that, as before, the distribution $`f(N)`$ has Levy exponent of $`\frac{1}{2}`$ since that result is universal for all symmetric walks. We perform an integral analogous to Eq. (14) and obtain a first passage time Levy exponent $`\beta _3=\frac{1}{8}`$. However, our numerical simulations yield a Levy exponent of $`0.220\pm 0.001`$. The exponent is closer to $`\frac{1}{4}`$ than $`\frac{1}{8}`$. Clearly, there is something wrong.

The method fails because there are now *correlations* between the $`\mathrm{\Delta }x_i`$ steps, where $`\mathrm{\Delta }x_i`$ is the change in $`x`$ during a $`\mathrm{\Delta }t_i`$ step. In the previous $`n=2`$ case, the velocity undergoes a simple random walk and so its sign is random after a zero-crossing. In this $`n=3`$ case, the velocity is the integral of a random walk. It is a differentiable function and so almost certainly changes sign when it hits zero. Consequently, the $`\mathrm{\Delta }x_i`$ steps alternate in sign. In addition, as the acceleration drifts about, the $`\mathrm{\Delta }x_i`$ steps increase in size. The acceleration at the beginning of each next time interval is larger on average than the previous one, making the next $`\mathrm{\Delta }t_i`$ and, therefore, the next $`\mathrm{\Delta }x_i`$ larger than the last.

A large number of alternating, increasing steps will reach zero quickly. From the numerical simulation, with each oscillation there is a definite probability of reaching zero which is almost constant. In other words, the distribution $`f(N)`$ is not a power law as in the $`n=2`$ case, but an exponential,

$$f(N)e^{\lambda _3(2N+1)}\mathrm{with}\lambda _3=0.79\pm \mathrm{\hspace{0.17em}0.02}.$$
(19)
See Table II and the corresponding Fig. 2. The number of velocity zero-crossings must be odd because the position can become negative only after the velocity has become negative. With this exponential distribution, the number of time intervals one must add up is actually quite small. However, adding a finite number of independently drawn Levy variables distributed with Levy exponent $`\frac{1}{4}`$ produces a flight with Levy exponent $`\frac{1}{4}`$, so it is quite surprising that adding $`N`$ of them together, where $`N`$ has a finite mean, leads to anything other than an exponent of $`\frac{1}{4}`$.

And yet, adding a small number of *correlated* Levy variables does shift the exponent. To quantitatively describe the correlation in magnitude, we need to determine how the initial acceleration sets the time scale for each $`\mathrm{\Delta }t_i`$ interval. We saw from Eq. (3) for the $`n=1`$ random walk that the time scale is the initial value squared. Therefore, the time scale between two acceleration zero-crossings is the initial acceleration squared. Because the time between velocity zero-crossings is comprised of many acceleration zero-crossings, the time scale between velocity zero-crossings is fixed in the same way.

This argument is somewhat general, so it is nice to verify that it is correct. The coefficient in the formula for the first passage time distribution was computed in Ref.  for the case of the $`n=2`$ random walk by solving the Fokker-Planck equation with the appropriate boundary conditions.

$$\rho _2(t;v_o,x_o\mathrm{\hspace{0.17em}0})dt\frac{\sqrt{v_o}}{t^{1/4}}\frac{dt}{t}=(\frac{v_o}{\sqrt{t}})^{\frac{1}{2}}\frac{dt}{t}$$
(20)
for large t. Note that $`t_{\mathrm{typ}}v_o^2`$. Since the acceleration in the $`n=3`$ case is exactly the same as the velocity in the $`n=2`$ case, we can translate this result into $`n=3`$ language by replacing $`v`$ with $`a`$. Therefore in the $`n=3`$ case, the square of the initial acceleration determines the time scale until the next zero-crossing.

So the square of the initial acceleration at each time step determines the scale of the next time step. Since the acceleration is an ordinary random walk, its square has a size proportional to the total elapsed time. Therefore, $`\mathrm{\Delta }t=T\times q`$ where $`q`$ is a unit Levy variable and $`T`$ is the total elapsed time. To describe this correlated process, we set the units of time so that $`\mathrm{\Delta }t_1`$ is a unit size Levy variable $`q_1`$. The next time step, $`\mathrm{\Delta }t_2`$, is no longer unit sized, but has a magnitude determined by the square of the acceleration, or, equivalently, the total elapsed time. So to find $`\mathrm{\Delta }t_2`$ we take $`\mathrm{\Delta }t_1`$ and multiply it by separate, independent, unit Levy variable $`q_2`$. To find $`\mathrm{\Delta }t_3`$, we take the total elapsed time, $`\mathrm{\Delta }t_1+\mathrm{\Delta }t_2`$, which is the expected square acceleration, and multiply by $`q_3`$. In equations,

$$\mathrm{\Delta }t_1=q_1$$
(21)
$`\mathrm{\Delta }t_2=(\mathrm{\Delta }t_1)q_2`$

$`\mathrm{\Delta }t_3=(\mathrm{\Delta }t_1+\mathrm{\Delta }t_2)q_3`$

$`\mathrm{\Delta }t_4=(\mathrm{\Delta }t_1+\mathrm{\Delta }t_2+\mathrm{\Delta }t_3)q_4.`$

The total time $`t`$ that has elapsed after $`N`$ steps is then

$$t=q_1\underset{i=2}{\overset{N}{}}(1+q_i).$$
(22)
The asymptotic distribution is extracted when $`q_i>>1`$. This means we need to determine the distribution of products of Levy variables. Taking the natural logarithm, the product becomes a sum. Define $`z_2=ln(\mathrm{\Delta }t_2)`$, $`f_1=ln(q_1)`$, and $`f_2=ln(q_2)`$. We use the following approximate distribution for both Levy variables,

$$P(q)dq=\frac{\alpha c_0^\alpha }{q^\alpha }\frac{dq}{q}$$
(23)
$`P(q)dq=0q<c_0.`$

with $`c_0=1`$. Keep in mind, in order to compute $`\beta _3`$, the Levy exponent $`\alpha `$ of the $`q_i`$ variables equals $`\frac{1}{4}`$. In terms of the transformed variables, the distribution of $`f_1`$ and $`f_2`$ is given by

$$P(f)=\alpha e^{\alpha f}.$$
(24)
If $`q_1`$ and $`q_2`$ had different Levy exponents, $`\alpha _1>\alpha _2`$, the distribution of the sum would be

$$_0^{z_2}𝑑f_1P(f_1)P(z_2f_1)=\frac{\alpha _1\alpha _2}{(\alpha _2\alpha _1)}(e^{\alpha _1z_2}e^{\alpha _2z_2}).$$
(25)
For large $`z_2`$, $`P(z_2)\alpha _2e^{\alpha _2z_2}`$, i.e. the smaller exponent dominates the tail. However, when $`\alpha _1=\alpha _2=\alpha `$, we have

$$P(z_2)=\alpha ^2z_2e^{\alpha z_2}.$$
(26)
In terms of $`q_1`$ and $`q_2`$, the probability distribution of $`q_1q_2`$ is

$$P(\mathrm{\Delta }t_2)d(\mathrm{\Delta }t_2)=\alpha ^2\frac{ln(\mathrm{\Delta }t_2)}{\mathrm{\Delta }t_2^\alpha }\frac{d(\mathrm{\Delta }t_2)}{\mathrm{\Delta }t_2}.$$
(27)
In general, the probability distribution for $`i`$ independent Levy variables multiplied together is

$$P(z_i)=\frac{\alpha ^i}{(i1)!}z_i^{i1}e^{\alpha z_i}.$$
(28)
The $`(i1)!`$ is required for normalization. Transforming back to the original $`\mathrm{\Delta }t_i`$ variables,

$$P(\mathrm{\Delta }t_i)d(\mathrm{\Delta }t_i)=\frac{\alpha ^i}{(i1)!}\frac{(ln\mathrm{\Delta }t_i)^{i1}}{\mathrm{\Delta }t_i^\alpha }\frac{d(\mathrm{\Delta }t_i)}{\mathrm{\Delta }t_i}.$$
(29)
We substitute $`i=N`$ and $`\alpha =\frac{1}{4}`$ and this distribution becomes $`L_N(t)`$— the probability distribution of the $`N`$th passage time of the velocity. Unlike the $`N`$th passage time of the acceleration, $`L_N(t)`$ is not a Levy stable distribution because of the correlations. More specifically, $`L_N(t)`$ is the distribution of a product of $`N`$ Levy variables $`q_i`$. The expression for the $`N`$th passage time we derived earlier is $`(q_i+1)`$, which differs in subleading terms from the previous expression.

Note that $`L_N(t)`$ coincides with the $`N`$th passage time distribution for the position in the $`n=2`$ case. Both the velocity for $`n=3`$ and the position for $`n=2`$ are integrals of an $`n=1`$ random walk. Observe that the distribution $`L_2(t)`$ is not the same as the distribution of the sum of two variables distributed as $`L_1(t)`$— there is an additional logarithm. We determined $`L_2(t)`$ numerically and it agrees with this form (Fig. 3). For higher $`N`$, there are $`N1`$ factors of $`ln(t)`$.

Before estimating $`\beta _3`$, we must make one last argument in the form of an approximation. Up until this point, $`L_N(t)dt`$ assumes that each step is free to be larger than the previous one. This is because we have not yet taken into account the conditional dependence of the correlation on the position being positive for the walk to persist. In fact, the generic situation when the steps are ever-increasing and alternating in sign is for there to be only one step. This means that walks that survive more than one step are special. For these walks, the second step is an exceptional pick from the distribution. The *easiest* way for the walk to survive is for the second step to be approximately the same size as the first step. In other words, the size of the first step cuts off the distribution of the second step. The result is that the second step shares the distribution of the first. For the third step, the position is positive, and now its time interval is distributed the same as the second step would have been. Generalizing our argument, any two *complete* time intervals can simply be lumped together as one and we conclude that the distribution of the total after $`N`$ steps is the same as the distribution of the step sizes after $`\frac{N}{2}`$ steps. In other words, $`L_N(t)dt`$ is not the correct distribution for $`N`$ steps, but $`L_{N/2}(t)dt`$ is. We present evidence for the even-numbered complete steps being typically small in Fig. 4. For $`n=3`$, the second step is, on average, about five times smaller than the first time interval. For $`n=4`$, the second time interval is typically twice as small as the first.

So we now have all the ingredients to estimate the Levy exponent $`\beta _3`$. We arrive at

$$\rho _3(t)dte^{\lambda _3}\frac{1}{t^{5/4}}\underset{N=0}{\overset{\mathrm{}}{}}\frac{1}{4^N}\frac{1}{N!}(ln(t))^Ne^{2\lambda _3N}dt$$
(30)
$$=e^{\lambda _3}\frac{1}{t^{5/4}}e^{\frac{e^{2\lambda _3}}{4}ln(t)}dt,$$
(31)
which gives $`\beta _3^{\mathrm{est}}=\frac{1}{4}\frac{e^{2\lambda _3}}{4}`$. Substituting the numerical value $`\lambda _3=0.79`$, we find the Levy index $`0.199`$. This result is a shift in the right direction from our original estimate, but it is still too large.

Figure 4 suggests that our estimation scheme may work better for $`n>3`$ and so we extend the method to $`n=4,5,\mathrm{}`$. As before, we define the two Levy flights of the time and position intervals between velocity zero-crossings. For $`n>3`$, the size correlation factor for the time interval Levy flight turns out to be the same factor of the elapsed time as in the $`n=3`$ case with each step retaining only the memory of the highest-order derivative. To be more precise, the initial acceleration at the start of each time interval for the $`n=3`$ case is replaced with the highest-order random variable. We ignore any other memory effects. Then, for the $`n`$th-order process, the step-size distribution for the $`\mathrm{\Delta }t_i`$ intervals are governed by the passage time exponent for the $`(n1)`$th-order process. Referring back to Eq. (28), $`i=N`$ and $`\alpha =\beta _{n1}`$. In addition, the correlation in the sign of the $`\mathrm{\Delta }x_i`$ still persists for the $`n`$th-order random walk so the distribution $`f(N)`$ remains exponential with a decay constant that we determine numerically. Therefore, the $`n\mathrm{th}`$ Levy exponent will be a small perturbation about the $`(n1)`$th first passage time exponent with a shift of $`\beta _{n1}e^{2\lambda _n}`$. For $`n=4`$ and $`n=5`$, we observe that our estimation scheme yields more accurate results. See Table II.

As in the $`n=2`$ case, for $`n3`$ we consider the last partial time step to be a full one. Given that $`L_{N/2}(t)dt`$ is the more appropriate distribution, the bulk of the first passage time for a given $`N`$ is taken up during the last two steps (up to subleading corrections). So the these last partial time steps should be distributed just as the first passage time. We have numerically verified this to be the case.

Majumdar et. al.  have constructed a completely different approximation scheme to compute the first passage time exponent for these higher-order random walks. While we refer the reader to Ref.  for the details, their scheme is rather accurate for $`\beta _2`$ and $`\beta _3`$.

## VI The $`n=\mathrm{}`$ Limit

As $`n`$ increases, we also expect the decay constant $`\lambda _n`$ to grow. In fact, in the limit that $`n`$ becomes large, the first passage time exponent reaches a limiting value. We consider the solution to Eq. (1) with initial conditions $`x(0)=x^{}(0)=\mathrm{}=x^n(0)=0`$. It is

$$x(t)=_0^t𝑑t_n_0^{t_n}𝑑t_{n1}\mathrm{}_0^{t_2}𝑑t_1\eta (t_1),$$
(32)
which is a linear function of $`\eta `$ and therefore is the inner product of $`\eta `$ with a kernel. In this case,

$$x(t)=_0^t𝑑t^{}\frac{(tt^{})^{(n1)}}{(n1)!}\eta (t^{}).$$
(33)
It is easy to see that this expression satisfies Eq. (1).

When solving for the time when $`x(t)`$ first becomes zero, it is permitted to rescale $`x(t)`$ by any finite number, even if that number is a function of time. This does not change the set of points where $`x(t)`$ is zero. Rescaling by $`t^n`$ so that $`y(t)=\frac{x(t)}{t^n}`$ and dropping constant factors,

$$y(t)=\frac{1}{t}_0^t𝑑t^{}(1\frac{t^{}}{t})^{n1}\eta (t^{}).$$
(34)
In the large $`n`$ limit, $`(1\frac{t^{}}{t})^n`$ is indistinguishable from an exponential in the region where it has the most weight ($`t^{}<nt`$). The result is that

$$y(t)=\frac{1}{t}_0^t𝑑t^{}e^{\frac{nt^{}}{t}}\eta (t^{}).$$
(35)
Majumdar it et al.  have noted that the $`n=\mathrm{}`$ random walk has the same first passage time exponent as a very different system — the solution of the two-dimensional diffusion equation with a random initial condition,

$$\frac{d\rho }{dt}=^2\rho $$
(36)
$`\rho (x,y,0)=\eta (x,y).`$

In polar coordinates, the value of $`\rho (0,t)`$ is given by

$$\rho (0,t)=\frac{1}{2\pi t}d^2\stackrel{}{r}e^{\frac{r^2}{2t}}\frac{1}{\sqrt{r}}\eta (r,\theta ),$$
(37)
where we have used $`\delta (x)\delta (y)=r^1\delta (r)\delta (\theta )`$. The angular integration leaves the integral in terms of a new random variable, $`\eta ^{}(r)=\frac{1}{\sqrt{2\pi }}𝑑\theta \eta (r,\theta )`$. We now have

$$\rho (0,t)=\frac{1}{\sqrt{2\pi }t}𝑑r\sqrt{r}e^{\frac{r^2}{2t}}\eta ^{}(r)$$
(38)
Transforming variables to $`u=r^2`$, we arrive at

$$\frac{1}{2}\frac{1}{t\sqrt{\pi }}𝑑ue^{\frac{u}{2t}}\eta (u),$$
(39)
where we used the fact that $`\delta (r)=\delta (u)2\sqrt{u}`$ so that $`\eta ^{}(r)=\eta (u)\sqrt{2}u^{\frac{1}{4}}`$.

Comparing this kernel with the kernel for the $`n=\mathrm{}`$ random walk, we see that they are indeed the same. As a result, the distribution of times for the diffusive field to change sign at one point is the same as the distribution of the first passage time for the $`n=\mathrm{}`$ random walk for long times, where the initial conditions are irrelevant.

## VII Physical Systems

Consider a planar crack front moving through rough, brittle material. In the quasistatic regime, the planar crack front is a line with long-range static elasticity. There is a uniform stress applied to the material, which drives the crack front forward. But there are also local, random variations in the fracture toughness which the crack front needs to overcome to move forward. These randomly varying forces are pinning forces. The competition between the global pulling and local pinning forces as mediated by the elasticity of the crack front determines its dynamics. For small values of the applied force, the pinning forces eventually dominate and the crack front remains stuck. For some finite value of the applied force, the pinning forces are not strong enough to keep the crack front in place and it moves forward. The boundary between these two phases is the depinning transition.

Approaching the depinning transition from the static side, a small increase in the applied force causes a tiny portion of the crack front to move forward. This portion causes a few more regions to move forward as well, only to be eventually stopped by more strongly pinned regions. The only role of the applied force here is to induce the initial motion of the tiny portion. It remains fixed until the local motion of the crack front stops. Near the depinning transition, there exists a sequence of these discrete, localized *avalanche events* with distribution of sizes $`\tau `$. The avalanche size is defined as the total amount of the crack front that has moved forward during these events. A power law distribution of avalanche sizes as the force approaches the critical force indicates a continuous, second-order-like depinning transition.

Near the depinning transition, the motion of the crack front is jerky. We therefore model its dynamics in terms of discrete space and time. Within an infinite-range model, where each crack front segment is coupled equally to every other, the spatial degrees of freedom along the crack front are averaged out. We define $`y_t`$ to be the total number of segments of the crack front that hop forward at discrete time $`t`$. Given that we are only concerned with the few segments that are on the verge of hopping forward, $`y_t`$ is Poissonian distributed. In this infinite-range limit, and in the dissipative regime, the mean of $`y_t`$ is determined by the number of segments that have hopped at the previous time step only. In addition, the mean of $`y_t`$ fluctuates statistically until it reaches zero, whereupon the motion stops. In the limit of large avalanches where things are changing slowly with time, the equation of motion for the avalanches at the depinning transition is

$$\frac{dy}{dt}=y^{\frac{1}{2}}\eta (t)$$
(40)
in the Ito representation (the time derivatives are forward differences). When $`y(t)`$ reaches zero, the avalanche is over. We refer the reader to Ref.  for a more detailed derivation.

So far, we have neglected any effects like elastic waves that are indeed present along the crack front. If a piece of the crack front moves forward, that motion creates an extra transient force on its neighbors as the elastic wave propagates along the crack front. The extra transient stress is called a stress overshoot. This effect is different than pure inertia, which would be an extra transient stress on the segment itself after it jumps. However, we demonstrate that the two effects are similar in this infinite-range model for a particular value of the stress overshoot . If we were to take into account elastic waves along the crack front, the mean of $`y_t`$ depends on several previous time steps. In other words, there are higher-order corrections to Eq. (39) that are usually irrelevant in the long time limit. The most significant terms are the first and second time derivatives, and the first time derivative is the only relevant term. If we fine-tune the coefficient of $`\frac{dy}{dt}`$ to zero, we arrive at a tricritical point. Here, the second-order derivative is the most relevant term and we obtain the equation of motion

$$\frac{d^2y}{dt^2}=y^{\frac{1}{2}}\eta (t).$$
(41)
The variable we are interested in is

$$\tau =_0^{t_f}𝑑ty(t),$$
(42)
the avalanche size. We emphasize that $`t_f`$ is determined by the boundary condition $`y(t_f)=0`$.

In order to find the distribution of $`\tau `$, we use the size of the avalanche as a new path-dependent time coordinate. We transform time on each particular path in the path-integral differently, in a way that depends on the history of $`y(t)`$. This transformation alters the density of $`y(t)`$’s at any given time, mixing together $`y(t)`$’s from different times so that they appear simultaneous. Because of this property, it is difficult to perform the transformation in the Fokker-Planck equation; but the transformation is easy in the path integral.

A reparameterization of time, even a path dependent one, does not affect the answers to questions that do not involve the time explicitly. The probability for a stochastic walker to arrive at point B from point A is independent of the global time, but the probability for a stochastic walker to be at point B at time $`t`$ is not. Similarly, the avalanche size of a given walk does not involve the time lapse of the walk as the upper cutoff of the integral is determined by the boundary condition on $`y(t)`$ and not the time explicitly.

Choosing the new time coordinate $`\tau `$ to tick at the new rate,

$$\frac{d\tau }{dt}=y,$$
(43)
and using the property of Gaussian noise, $`\eta (\tau (t))=\frac{\eta (t)}{\sqrt{\frac{d\tau }{dt}}}`$, the (Ito) equation pair

$$\frac{dv}{dt}=y^{\frac{1}{2}}\eta (t)$$
(44)
$`{\displaystyle \frac{dy}{dt}}=v`$

become

$$\frac{dv}{d\tau }=\eta (t)$$
(45)
$`y{\displaystyle \frac{dy}{d\tau }}=v.`$

Replacing $`y`$ with $`x=\frac{y^2}{2}`$ completes the transformation. The zeros of $`y`$ are also the zeros of $`x`$ so the value of $`\tau `$ when $`x`$ is zero is then the avalanche size. We therefore find the distribution of avalanche sizes to be the distribution of first passage times $`\tau `$ for the $`n=2`$ random walk. It has a Levy exponent $`\frac{1}{4}`$.

We may perform the same transformation on the Fokker-Planck equation directly, although the motivation then becomes obscure. Consider the following probability-conserving equation,

$$\frac{P}{t}=y\frac{P}{\tau }+\frac{1}{2}y\frac{^2P}{v^2}v\frac{P}{y},$$
(46)
for a new quantity $`P_t(x,v,\tau )`$. If we integrate $`P_t(x,v,\tau )`$ over all $`\tau `$, we recover the original Fokker-Planck equation. On the other hand, if integrate $`P_t(x,v,\tau )`$ with respect to $`t`$ from zero to infinity, we obtain

$$\frac{K}{\tau }=\frac{1}{2}\frac{^2K}{v^2}\frac{v}{y}\frac{K}{y},$$
(47)
where $`K(v,y,\tau )=_0^{\mathrm{}}𝑑tP_t(v,y,\tau )`$. Defining $`x=\frac{y^2}{2}`$, as before, gives

$$\frac{K}{\tau }=\frac{1}{2}\frac{^2K}{v^2}v\frac{K}{x},$$
(48)
which is the Fokker-Planck equation for the $`n=2`$ random walk. This demonstrates the equivalence of the two problems.

It should now be clear that the avalanche size exponent for Eq. (39) is $`\frac{1}{2}`$ since the same time reparameterization may be performed in that case also. This relates the avalanche size exponent to the first passage time of the $`n=1`$ random walk.

## VIII Conclusions

We have presented a relation between the shift in the first passage time exponent and the decay rate of the probability of $`N`$ velocity zero-crossings for the $`n`$th random walk. The method can be exact for $`n=2`$, but not for $`n3`$ because of correlations. More work is needed to estimate the decay constant of $`f(N)`$ since this would determine the convergence rate of the exponents to the $`n=\mathrm{}`$ value. In addition, we do not have any bounds on the error of the first passage time exponent as the approximation is uncontrolled. We also do not know which types of correlations for Levy flights lead to shifts in exponents and which do not. A classification of correlations in terms of these two categories might be useful.

We have opted to slice up time in terms of velocity zero-crossings, instead of using the global time in the Fokker-Planck equation. This approach allows us to analyze the higher-order random walks in terms of a one-dimensional Levy process since the phase variables are the position and the time only. The remaining variables are subsumed in the $`(n1)`$th-order random walk that the velocity undergoes between the zero-crossings.

With Eq. (42) we have also reparameterized the global time, but in a different way. We use a path-dependent time transformation to analyze the avalanche statistics of a nonlinear second-order random walk in terms of the first passage time exponent of the linear second-order random walk. This allows us to give a novel physical interpretation of the first passage time exponent for the $`n=2`$ case.

Since the avalanche size does not depend on the global time, it has been known to us that one can find the stationary solution of the Fokker-Planck equation corresponding to the nonlinear random walk in Eq. (40) and then arrive at the avalanche size exponent after invoking a simple scaling argument . Since this method gives a new derivation of the avalanche size exponent, it indirectly gives a new derivation of the first passage time exponent for the $`n=2`$ case. This insight may be useful to construct a derivation of $`\beta _n`$ for $`n3`$.

Finally, these higher-order random walks could be relevant for more general forms of dynamic stress transfer along a crack front, or other extended elastic objects. For instance, one can approach a fine-tuned critical point where the coefficients of the first and second time derivative terms are zero. Then the third-order time derivative is the most relevant term, and we arrive at a higher-order critical point. While the time reparameterization given by Eq. (42) does not allow us to equate the avalanche size exponent with the first passage time exponent in this case, we remain optimistic that either exponent may be relevant for various physical systems.

## Acknowledgements

The authors would like to thank P.A. Argyres and Lev Kaplan for useful discussion. J.M.S. would also especially like to thank Daniel S. Fisher for suggesting the problem, for many helpful discussions and for useful comments on the manuscript. J.M.S. has been supported by the National Science Foundation via DMR-9630064, DMR-9976621, and DMR-9809363. R.M. has been supported by PHY95-13717.

## The Kernel Of Higher-Order Random Walks

We begin with the partition function for all higher-order random walks,

$$Z=D[\eta (t)]e^{\frac{1}{2}_0^T(\eta (t))^2}.$$
(49)
We change variables from $`\eta (t)`$ to $`x(t)`$, where $`\frac{d^nx}{dt^n}=\eta (t)`$, and use forward difference time derivatives. We change variables with the insertion $`_{x(0)=x_o}^{x(T)=x_f}D[x(t)]\delta (\frac{d^nx(t)}{dt^n}\eta (t))`$ to arrive at

$$_{x(0)=x_o}^{x(T)=x_f}D[x(t)]e^{\frac{1}{2}_0^T(\frac{d^nx(t)}{dt^n})^2}.$$
(50)
We then minimize the action to find the classical equation of motion with vanishing variations on the boundary for all derivatives up to the $`(n1)`$th-derivative. In the $`n=2`$ case we find the equation of motion,

$$\frac{d^4x(t)}{dt^4}=0.$$
(51)
We then impose the following constraints on the $`x(t)`$ trajectory, $`x(0)=x_o`$, $`x(T)=x_f`$, $`x^{}(0)=v_o`$, and $`x^{}(T)=v_f`$, yielding

$`x(t)=x_o+v_ot+{\displaystyle \frac{(3x_f3x_ov_fT2v_oT)}{T^2}}t^2`$ (52)
$`+{\displaystyle \frac{(v_f+v_o\frac{2}{T}(x_fx_o))}{T^2}}t^3.`$ (53)
The path integral is quadratic, so the partition function is proportional to the exponential of the classical action. Substituting the classical solution into the action and integrating gives the propagation kernel,

$`P_t(v_f,x_f;v_o,x_o)={\displaystyle \frac{\sqrt{3}}{\pi t^2}}e^{\frac{(v_fv_o)^2}{2t}}`$ (54)
$`e^{\frac{6((x_fx_o)\frac{1}{2}(v_f+v_o)t)^2}{t^3}}.`$ (55)
We normalize the propagation kernel to have integral unity. Notice that the kernel is the product of a Gaussian in $`v`$ and a gaussian in $`x`$. Given that $`x(t)`$ is the time integral of $`v(t)`$ and that the sum of Gaussian variables is Gaussian, if the velocity is Gaussian, so is the position. Not only does the position spread faster than the velocity, but it also drifts with time with a coefficient that depends on the velocity. This method works for higher $`n`$ giving us the following propagation kernels for $`n=3,4,5`$:

$`P_t(a_f,v_f,x_f;a_o,v_o,x_o){\displaystyle \frac{1}{t^{\frac{9}{2}}}}e^{\frac{1}{2t}(a_fa_o)^2}e^{\frac{6}{t^3}((v_fv_o)\frac{1}{2}(a_f+a_o)t)^2}`$ (56)
$`e^{\frac{360}{t^5}((x_fx_o)\frac{1}{2}(v_o+v_f)t+\frac{1}{12}(a_fa_o)t^2)^2}`$ (57)
for $`n=3`$,

$`P_t(s_f,a_f,v_f,x_f;s_o,a_o,v_o,x_o){\displaystyle \frac{1}{t^8}}e^{\frac{(s_fs_o)^2}{2t}}`$ (58)
$`e^{\frac{6}{t^3}((a_fa_o)\frac{1}{2}(s_f+s_o)t)^2}`$ (59)
$`e^{\frac{360}{t^5}((v_fv_o)\frac{1}{2}(a_f+a_o)t+\frac{1}{12}(s_fs_o)t^2)^2}`$ (60)
$`e^{\frac{50400}{t^7}((x_fx_o)\frac{1}{2}(v_f+v_o)t+\frac{1}{10}(a_fa_o)t^2\frac{1}{120}(s_f+s_o)t^3)^2}`$ (61)
for $`n=4`$ with $`\frac{da(t)}{dt}=s(t)`$, and

$`P_t(u_f,s_f,a_f,v_f,x_f;u_o,s_o,a_o,v_o,x_o){\displaystyle \frac{1}{t^{\frac{25}{2}}}}e^{\frac{(u_fu_o)^2}{2t}}e^{\frac{6}{t^3}((s_fs_o)\frac{1}{2}(u_f+u_o)t)^2}`$ (62)
$`e^{\frac{360}{t^5}((a_fa_o)\frac{1}{2}(s_f+s_o)t+\frac{1}{12}(u_fu_o)t^2)^2}e^{\frac{50400}{t^7}((v_fv_o)\frac{1}{2}(a_f+a_o)t+\frac{1}{10}(s_fs_o)t^2\frac{1}{120}(u_f+u_o)t^3)^2}`$ (63)
$`e^{\frac{12700800}{t^9}((x_fx_o)\frac{1}{2}(v_f+v_o)t+\frac{3}{28}(a_fa_o)t^2\frac{1}{84}(s_f+s_o)t^3+\frac{1}{1680}(u_fu_o)t^4)^2}`$ (64)
for $`n=5`$ with $`\frac{ds(t)}{dt}=u(t)`$.


format: text

text:
# Subvarieties of abelian varieties

## 1. Jacobians

Let $`C`$ be a smooth, complete, irreducible curve of genus $`g`$ over the complex numbers. The jacobian $`JC=Pic^0C`$ of $`C`$ is the connected component of its Picard group parametrizing degree $`0`$ invertible sheaves. For any nonnegative integer $`e`$, the choice of an invertible sheaf $``$ of degree $`e`$ on $`C`$ gives a morphism

$$\begin{array}{ccc}\hfill \varphi _{}:C^{(e)}& & JC\hfill \\ \hfill D_e& & 𝒪_C(D_e)^1\hfill \end{array}$$

where $`C^{(e)}`$ is the $`e`$-th symmetric power of $`C`$. For $`eg`$, such a morphism is surjective. When $`e=g1`$, the image of $`C^{(g1)}`$ by $`\varphi _{}`$ is a theta divisor on $`JC`$ which we will denote by $`\mathrm{\Theta }_C`$ (always well-determined up to translation). For any $`e`$ between $`1`$ and $`g`$ the image of $`\varphi _{}`$ has class $`\frac{[\mathrm{\Theta }]^{ge}}{(ge)!}`$. If $`e=1`$, the map $`\varphi _{}`$ is an embedding and its image is called an Abel curve. By a theorem of Matsusaka , the dimension $`1`$ minimal class is represented by an algebraic curve $`C`$ in $`(A,\mathrm{\Theta })`$ if and only if $`(A,\mathrm{\Theta })`$ is the polarized jacobian $`(JC,\mathrm{\Theta }_C)`$ of $`C`$. The higher-dimensional analogue of this theorem has the following counterexample: by a result of Clemens and Griffiths, the intermediate jacobian of a smooth cubic threefold in $`^4`$ is not the jacobian of a curve but it contains a surface (the image of the Fano variety of lines in the cubic threefold) whose cohomology class is the dimension $`2`$ minimal class . Debarre has shown that for any $`e`$ strictly between $`1`$ and $`g1`$ jacobians form an irreducible component of the family of ppav in which the dimension $`e`$ minimal class is represented by an algebraic subvariety .

As the above suggests, not every ppav is a jacobian. In fact, the moduli space of ppav of dimension $`g`$ has dimension $`\frac{g(g+1)}{2}`$ whereas the moduli space of curves of genus $`g2`$ has dimension $`3g3`$. As soon as $`g4`$, we have $`\frac{g(g+1)}{2}>3g3`$ and so not all ppav of dimension $`g`$ are jacobians. The question then becomes how else can one parametrize ppav? The first step of a generalization of the notion of jacobian is the construction of a Prym variety which we describe below.

## 2. Pryms

Suppose that $`C`$ is a smooth, complete and irreducible curve of genus $`g+1`$ with an étale double cover $`\pi :\stackrel{~}{C}C`$. Then $`\stackrel{~}{C}`$ has genus $`\stackrel{~}{g}:=2g+1`$. Let $`\sigma :\stackrel{~}{C}\stackrel{~}{C}`$ be the involution of the cover $`\pi `$. The involution $`\sigma `$ acts on the jacobian $`J\stackrel{~}{C}`$ and the Prym variety $`P`$ of $`\pi `$ is an abelian variety of dimension $`g`$ defined by

$$P:=im(\sigma 1)J\stackrel{~}{C}.$$

The principal polarization of $`J\stackrel{~}{C}`$ induces twice a principal polarization on $`P`$ which we will denote $`\mathrm{\Xi }`$.

A priori, there are two ways of obtaining subvarieties of a Prym variety: by projection and intersection. Since $`P`$ is a quotient of $`J\stackrel{~}{C}`$ via $`\sigma 1`$, we can take the images of subvarieties of $`J\stackrel{~}{C}`$ in $`P`$. For $`e=1`$, we can take the images in $`P`$ of Abel curves in $`J\stackrel{~}{C}`$. We obtain in this way embeddings of $`\stackrel{~}{C}`$ in $`P`$ whose images are called Prym-embedded curves. The class of a Prym-embedded curve is $`2\frac{[\mathrm{\Xi }]^{g1}}{(g1)!}`$. Welters classfied all curves of class twice the minimal class in a ppav . The list is short but it is not limited to Prym-embedded curves. Therefore the analogue of Matsusaka’s theorem is false for curves representing twice the minimal class. For $`e>1`$, using Pontrjagin product, it is easy to see that the projection in $`P`$ of the image of $`\varphi _{}`$ has class $`2^e\frac{[\mathrm{\Xi }]^{ge}}{(ge)!}`$.

Secondly, since the Prym variety is also a subvariety of $`J\stackrel{~}{C}`$, we can intersect the images of the symmetric powers of $`\stackrel{~}{C}`$ with it. Since $`\mathrm{\Theta }_{\stackrel{~}{C}}`$ induces twice $`\mathrm{\Xi }`$ on $`P`$, for $``$ of degree $`e`$, the intersection of the image of $`\varphi _{}`$ with $`P`$ has class $`2^{\stackrel{~}{g}e}\frac{[\mathrm{\Xi }]^{\stackrel{~}{g}e}}{(\stackrel{~}{g}e)!}`$ , provided that the intersection is proper. There is at least one case in which one can do better: there are a finite number of theta divisors in $`J\stackrel{~}{C}`$ whose intersection with $`P`$ is $`2\mathrm{\Xi }`$ as a divisor (see for instance ). In particular, we have a nice way of parametrizing the theta divisor of the ppav $`(P,\mathrm{\Xi })`$.

It can easily happen that the intersection $`\varphi _{}(\stackrel{~}{C}^{(e)})P`$ is not proper. In such a case the cohomolgy class of the resulting subvariety needs to be determined by other means. Such subvarieties appear in the work of Recillas , Donagi , Clemens-Griffiths  and Beauville . Following Beauville, we shall call them special subvarieties. A different way of defining a special subvariety of a Prym variety which allows one to compute its cohomology class is as follows , ,  and . Let $`g_e^r`$ be a complete linear system of dimension $`r`$ and degree $`e`$ on $`C`$. Let $`L`$ be the corresponding invertible sheaf on $`C`$ and let $``$ be an invertible sheaf on $`\stackrel{~}{C}`$ whose Norm is $`Nm()=L`$ (i.e., if $`𝒪_{\stackrel{~}{C}}(D)`$, then $`L𝒪_C(\pi _{}D)`$). Consider $`g_e^r`$ as a subvariety of $`C^{(e)}`$, isomorphic to $`^r`$. Assuming that $`g_e^r`$ contains reduced divisors, the inverse image of $`g_e^r`$ in $`\stackrel{~}{C}^{(e)}`$ is reduced. It splits as the union of two connected components whose images in $`J\stackrel{~}{C}`$ by $`\varphi _{}`$ are contained in $`P`$ and the translate $`P^{}`$ of $`P`$ such that $`PP^{}`$ is the kernel of the Norm $`Nm:J\stackrel{~}{C}JC`$. Therefore, after translating one of these subvarieties, we obtain two subvarieties of $`P`$. They are isomorphic if $`e`$ is odd but not if $`e`$ even. They both have cohomology class $`2^{e2r1}\frac{[\mathrm{\Xi }]^{gr}}{(gr)!}`$ at least when $`1e2g+1`$ and $`e>2r`$ .

To see that these special subvarieties are indeed the nonproper intersections that we mentioned above, one needs to note that the fibers of the map $`J\stackrel{~}{C}JC`$ are translates of $`PP^{}`$. The special subvarieties are intersections of these fibers with images of maps $`\varphi _{}`$. Note that when the $`g_e^r`$ is nonspecial, i.e., $`h^1(g_e^r)=0`$, the special subvariety is in fact a proper intersection and Beauville’s cohomology class is equal (as it should be) to the cohomology class of the proper intersection above. As we shall see below, looking at special subvarieties as such nonproper intersections allows us to define them for arbitrary ppav.

For $`g5`$, all ppav are Prym varieties (in the generalized sense of Beauville ). For $`g5`$ Prym varieties of dimension $`g`$ depend on the same number of moduli as curves of genus $`g+1`$, meaning $`3g`$ moduli. Therefore, for $`g6`$, a general ppav is not a Prym variety. So we need to find a different way to parametrize a ppav.

## 3. Prym-Tjurin varieties

Again, one would want to use a construction involving curves. Higher degree coverings $`\stackrel{~}{C}C`$ do not yield general ppav because the dimension of the Prym variety (defined as the connected component of $`𝒪_{\stackrel{~}{C}}`$ of the kernel of the Norm (or pushforward on divisors) $`J\stackrel{~}{C}JC`$) is too high and therefore the families of ppav that one would obtain are too small, their dimensions being the dimensions of the moduli spaces for the bottom curves $`C`$. Looking back at a Prym variety, we note that it was defined as a special type of abelian subvariety of a jacobian. An abelian subvariety $`A`$ of a jacobian $`(JX,\mathrm{\Theta }_X)`$ such that $`\mathrm{\Theta }_X`$ induces $`m`$-times a principal polarization $`\mathrm{\Theta }`$ on $`A`$ is called a Prym-Tjurin variety. Welters has proved that all ppav are Prym-Tjurin varieties .

To say that $`\mathrm{\Theta }_X`$ induces $`m`$-times $`\mathrm{\Theta }`$ on $`A`$ is equivalent to saying that the class of the image of an Abel curve in $`A`$ is $`m`$-times the minimal class for curves . Here we are taking the image of an Abel curve by the composition

$$JX\stackrel{}{}\widehat{JX}\widehat{A}\stackrel{}{}A$$

where $`\widehat{}`$ denotes the dual abelian variety, the first and the last map are induced by the polarizations $`\mathrm{\Theta }_X`$ and $`\mathrm{\Theta }`$ respectively, and the middle map is the transpose of the embedding of $`A`$ in $`JX`$.

Therefore, finding a structure of Prym-Tjurin variety on a ppav is equivalent to finding a reduced and irreducible curve $`\overline{X}`$ in $`A`$ representing $`m`$-times the minimal class and such that $`A`$ embeds in the jacobian of the normalization $`X`$ of $`\overline{X}`$. Given such a structure, we can find subvarieties of $`A`$ as in the case of Prym varieties: by projection and intersection. Since the class of the image of an Abel curve is $`m`$-times the minimal class, Pontrjagin product shows that the projection in $`P`$ of the image of $`\varphi _{}`$ has class $`m^e\frac{[\mathrm{\Theta }]^{ge}}{(ge)!}`$. Since $`\mathrm{\Theta }_X`$ induces $`m`$-times $`\mathrm{\Theta }`$ on $`A`$, for $``$ of degree $`e`$, the intersection of the image of $`\varphi _{}`$ with $`P`$ has class $`m^{g_Xe}\frac{[\mathrm{\Theta }]^{g_Xe}}{(g_Xe)!}`$, provided that the intersection is proper ($`g_X`$ is the genus of $`X`$). However, unlike Prym varieties, it is not clear whether one can find translates of $`\mathrm{\Theta }_X`$ whose intersection with $`A`$ is $`m`$-times a theta divisor. Kanev has shown that this is possible under a restrictive hypothesis which we explain below.

Any abelian subvariety of $`JX`$ is the image of an endomorphism of $`JX`$ (which is not unique). The datum of an endomorphism of $`JX`$ is equivalent to the datum of a correspondence, i.e., a divisor in $`X\times X`$, up to addition and subtraction of the fibers of the two projections. This is best seen as follows. Start with a divisor $`DX\times X`$. To $`D`$ one can associate an endomorphism of $`JX`$ in the following way

$$\begin{array}{cccc}\psi _D:& \hfill JX& & JX\hfill \\ & \hfill 𝒪_X(E)& & 𝒪_X(p_{2}^{}{}_{}{}^{}((p_1^{}E)D))\hfill \end{array}$$

where $`p_1`$ and $`p_2`$ are the two projections $`X\times XX`$. If $`D`$ is linearly equivalent to a sum of fibers of $`p_1`$ and $`p_2`$, then $`\psi _D`$ is the zero endomorphism. If we exchange the roles of $`p_1`$ and $`p_2`$ in the above definition then $`\psi _D`$ is replaced by its image under the Rosati involution. The correspondence $`D`$ is said to be symmetric if there are (not necessarily effective) divisors $`a`$ and $`b`$ on $`X`$ such that $`DD^t`$ is linearly equivalent to $`p_1^{}(a)+p_2^{}(b)`$, where $`D^t`$ is the transpose of $`D`$, i.e., the image of $`D`$ under the involution exchanging the two factors of $`X\times X`$. So $`D`$ is symmetric if and only if $`\psi _D`$ is fixed by the Rosati involution. We shall assume that this is the case. This is not restrictive since any abelian subvariety of $`JX`$ is always the image of an endomorphism which is fixed by the Rosati involution (see e.g. ).

Kanev  has shown that if the endomorphism can be represented by a symmetric fixed-point-free correspondence $`D`$ (i.e., the support of $`D`$ does not intersect the diagonal of $`X\times X`$), then one can find theta divisors $`\mathrm{\Theta }_X`$ such that $`\mathrm{\Theta }_X|_A=m\mathrm{\Theta }`$ as divisors. Furthermore, fixing an invertible sheaf $`_0Pic^{g_X1}X`$, an invertible sheaf $`PJX`$ is on $`\mathrm{\Theta }`$ if and only if $`h^0(_0)m`$ and $`\mathrm{\Theta }`$ if and only if $`h^0(_0)=0`$. This gives a nice parametrization of $`\mathrm{\Theta }`$ and even allows one to analyze the singularities of $`\mathrm{\Theta }`$. It is not known however, whether every ppav is a Prym-Tjurin variety for a (symmetric) fixed-point-free correspondence. In addition, two correspondences could induce the same endomorphism of $`JX`$ while one is fixed-point-free and the other is not. In general it is difficult to determine whether a given endomorphism can be induced by a fixed-point-free correspondence.

As we noted above, we can generalize the notion of special subvarieties to Prym-Tjurin varieties by defining them to be non-proper intersections of $`A`$ with images of symmetric powers of $`X`$. It would be interesting to compute the cohomology classes of these special subvarieties and see whether the analogue of Beauville’s formula holds, meaning, the cohomology class of a special subvariety of dimension $`r`$ is $`m^{e2r1}\frac{[\mathrm{\Theta }]^{gr}}{(gr)!}`$.

Welters showed that every principally polarized abelian variety is a Prym-Tjurin variety . Birkenhacke and Lange showed that every principally polarized abelian variety is a Prym-Tjurin variety for an integer $`m3^g(g1)!`$ (see  page 374 Corollary 2.4)111Their proof uses $`3`$-theta divisors. Using the fact that a general $`2`$-theta divisor is smooth, the exact same proof would give $`m2^g(g1)!`$. For abelian varieties with a smooth theta divisor, the same proof would give $`m(g1)!`$. One needs the Lefschetz hyperplane theorem which also works for mildly singular theta divisors, see e.g.  Chapter $`2`$..

## 4. Deforming curves

The question is to find the smallest integer $`m`$ for which $`m\frac{[\mathrm{\Theta }]^{g1}}{(g1)!}`$ can be represented by an algebraic curve. This naturally defines a stratification of the moduli space $`𝒜_g`$ of ppav. Using results of Kanev, Debarre shows that if $`(A,\mathrm{\Theta })`$ is the Prym-Tjurin variety for a symmetric fixed-point-free correspondence, then either $`Sing(\mathrm{\Theta })`$ is empty or its dimension is at least $`g2m2`$. Since the theta divisor of a general ppav is smooth, this suggests that, for a general ppav $`A`$, the smallest integer $`m`$ for which there is a curve of class $`m\frac{[\mathrm{\Theta }]^{g1}}{(g1)!}`$ in $`A`$ which in addition gives $`A`$ a structure of Prym-Tjurin variety should be at least $`\frac{g1}{2}`$. It is unlikely however that this bound is effective. Debarre has proved in  that the smallest integer $`m`$ for which $`m\frac{[\mathrm{\Theta }]^{g1}}{(g1)!}`$ is the class of an algebraic curve is at least $`\sqrt{\frac{g}{8}}\frac{1}{4}`$ if $`(A,\mathrm{\Theta })`$ is general.

The difficulty is to produce curves in ppav in nontrivial ways. One approach that we have considered is to deform curves in jacobians of curves out of the jacobian locus. More precisely, let $`C`$ be a curve of genus $`g`$ with a $`g_d^1`$ (a pencil of degree $`d`$). Define

$$X_e(g_d^1):=\{D_e:DC^{(de)}\text{ such that }D_e+Dg_d^1\}C^{(e)}$$

(for the precise scheme-theoretical definition see  when $`e=2`$ and  for $`e>2`$). If $`de+1`$, the restriction of a given morphism $`\varphi _{}`$ to $`X_e(g_d^1)`$ is nonconstant and so we can map $`X_e(g_d^1)`$ to $`JC`$. The cohomology class of the image of $`X_e(g_d^1)`$ in $`JC`$ is $`m`$-times the minimal class with $`m=\left(\genfrac{}{}{0pt}{}{d2}{e1}\right)`$. Given a one-parameter infinitesimal deformation of the jacobian of $`C`$ out of the jacobian locus $`𝒥_g`$ we ask when the curve $`X_e(g_d^1)`$ deforms with it. Infinitesimal deformations of $`JC`$ are parametrized by $`H^1(T_{JC})`$ where $`T_{JC}`$ is the tangent sheaf of $`JC`$. The principal polarization $`\mathrm{\Theta }_C`$ provides an isomorphism between $`H^1(T_{JC})`$ and the second tensor power $`H^1(𝒪_C)^2`$. Under this isomorphism the globally unobstructed deformations of the pair $`(JC,\mathrm{\Theta }_C)`$ are identified with the symmetric square $`S^2H^1(𝒪_C)`$. Therefore any quadric in the canonical space of $`C`$ defines a linear form on the space of these infinitesimal deformations. When we say that an infinitesimal deformation $`\eta S^2H^1(𝒪_C)`$ is in the annihilator of a quadric, we mean that it is in the kernel of the corresponding linear form. We prove the following in

Theorem Suppose $`C`$ nonhyperelliptic and $`d4`$. If the curve $`X_2(g_d^1)`$ deforms out of $`𝒥_g`$ then

1. either $`d=4`$
2. or $`d=5`$, $`h^0(g_5^1)=3`$ and $`C`$ has genus $`5`$ or genus $`4`$ and only one $`g_3^1`$.

In the case $`g=5`$ if $`X_2(g_5^1)`$ deforms in a direction $`\eta S^2H^1(𝒞_C)`$ out of $`𝒥_5`$, then $`\eta `$ is in the annihilator of the quadric $`_{Dg_5^1}D`$.

Here $`D`$ denotes the span of the divisor $`D`$ in the canonical space of $`C`$. For $`d=3`$ the image of $`X_2(g_3^1)`$ in $`JC`$ is an Abel curve and so by the result of , the curve $`X_2(g_3^1)`$ cannot deform out of $`𝒥_g`$. For $`d=4`$, it follows from the theory of Prym varieties that $`X_2(g_4^1)`$ deforms out of $`𝒥_g`$ (into the locus of Prym varieties): in fact $`X_2(g_4^1)`$ is a Prym-embedded curve . For $`d=5`$, $`h^0(g_d^1)=3`$ and $`g=4`$ (with only one $`g_3^1`$) or $`g=5`$ we believe that $`X_2(g_5^1)`$ deforms out of $`𝒥_g`$ but we do not have a proof of this. An interesting question is what are these deformations of $`(JC,\mathrm{\Theta }_C)`$ into which $`X_2(g_5^1)`$ deforms. Can one describe them in a concrete geometric way.

For $`e>2`$, the analogous result would be the following. The curve $`X_e(g_d^1)`$ deforms out of $`𝒥_g`$ if and only if

* either $`e=h^0(g_d^1)`$ and $`d=2e`$
* or $`e=h^0(g_d^1)1`$ and $`d=2e+1`$.

We expect this to be true most of the time. There could, however, be special pairs $`(C,g_d^1)`$ for which the curve $`X_e(g_d^1)`$ deforms out of the jacobian locus but $`g_d^1`$ does not verify the above conditions. For instance, so far my calculations  seem to indicate that if there is a divisor $`DX`$ with $`h^0(D)2`$, then $`X`$ might deform in directions $`\eta `$ whose images in the projectivization $`(S^2H^1(𝒪_C))`$ are in the span of the image of $`_{D^{}|D|}D^{}`$. Finally, note that a standard Brill-Noether calculation shows that for general curves of genus $`7`$, the smallest $`d`$ for which they can have $`g_d^1`$’s satisfying $`d=2h^0(g_d^1)`$ is $`d=2g4`$. In such a case the class of $`X_e(g_d^1)`$ is $`m`$-times the minimal class with $`m=\left(\genfrac{}{}{0pt}{}{2g6}{g3}\right)`$ which is then what we would find for a general ppav. We address the case $`e>2`$ in .

## 5. The genus

The cohomology class is one discrete invariant that one can associate to a curve in a ppav. Another discrete invariant is the genus of the curve. We refer the reader to the nice paper by Bardelli, Ciliberto and Verra  for a discussion of this.


text:
# The 𝑈(1)-invariant field theories with normal field operators

## 1.  INTRODUCTION

In quantum field theory there are today many general theorems which give limitations to possible field theories, e.g. the CPT theorem, the spin and statistics theorem .

In this paper we investigate the implications of the $`U(1)`$-invariance of a field theory. Our assumptions will be very limited, namely we consider a normal field operator $`\mathrm{\Phi }`$ (an operator which commutes with its adjoint) such that the vacuum state is $`U(1)`$-invariant on the algebra generated by $`\mathrm{\Phi }`$ and $`\mathrm{\Phi }^{}`$. Particularly we make no assumptions concerning the structure of the spacetime.

We show that from such an invariance follows the existence of particles and antiparticles associated to the same physical state. The field operator turns out to be a sum of a particle creation and an antiparticle annihilation operators. We study in details the case when creation and annihilation operators factorize and show that such operators are closely related to $`q`$-deformed commutation relations.

Our method bases on an analysis of orthogonal polynomials in two variables .

## 2.  $`U(1)`$-INVARIANT FIELD THEORIES

We consider a field theory with a Hilbert space $``$ and a field operator $`\mathrm{\Phi }:`$. We assume that there is a state $`\varphi :𝒜𝐂`$ (where $`𝒜`$ is an algebra generated by operators $`\mathrm{\Phi }`$ and $`\mathrm{\Phi }^{}`$) which is $`U(1)`$-invariant, i.e. if in any expression in $`𝒜`$ we replace operators $`\mathrm{\Phi }`$ by $`e^{is}\mathrm{\Phi }`$ and $`\mathrm{\Phi }^{}`$ by $`e^{is}\mathrm{\Phi }^{}`$, the value of the state should not change for any $`s𝐑`$.

Usually the state is simply a vacuum expectation: $`\varphi (S)=\mathrm{\Omega }|S\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ denotes the vacuum. In this case the $`U(1)`$-invariance of a state means that the vacuum is $`U(1)`$-invariant as well, so the $`U(1)`$-symmetry of the theory remains unbroken.

Assume that $`\mathrm{\Phi }`$ is a normal operator, i.e. $`\mathrm{\Phi }\mathrm{\Phi }^{}=\mathrm{\Phi }^{}\mathrm{\Phi }`$. This assumption is usually fulfilled in the theories of bosonic particles .

By standard methods of functional analysis we can assign to $`\mathrm{\Phi }`$ a measure $`\mu `$ with support on the complex plane $`𝐂`$ such that

$$\varphi (\mathrm{\Phi }^k\mathrm{\Phi }^l)=_𝐂z^k\overline{z}^l𝑑\mu (z).$$

The $`U(1)`$-invariance of the state $`\varphi `$ implies that the measure $`\mu `$ is rotation invariant.

Let us consider polynomials on the complex plane of two variables: $`z`$ and $`\overline{z}`$. The following theorem proven by Zygmunt  holds:

###### Theorem 1.

Let polynomials $`P_{k,l}(z,\overline{z})`$ form a system of orthonormal polynomials given by Gram–Schmidt procedure applied to the sequence

$$1,z,\overline{z},z^2,z\overline{z},\overline{z}^2,\mathrm{},z^n,z^{n1}\overline{z},\mathrm{},\overline{z}^n,$$

with respect to the rotation invariant measure $`\mu `$ on $`𝐂`$ (the polynomial $`P_{k,l}`$ has the leading term proportional to $`z^k\overline{z}^l`$).

For all nonnegative integers $`k,l`$ they fulfill the following recurrence relations:

$$zP_{k,l}=\alpha _{k,l}P_{k+1,l}+\alpha _{l1,k}P_{k,l1},$$

$$\overline{z}P_{k,l}=\alpha _{l,k}P_{k,l+1}+\alpha _{k1,l}P_{k1,l}.$$

The coefficients $`\alpha _{k,l}`$ are positive for $`k,l0`$ and $`\alpha _{1,l}=0`$. Furthermore for all nonnegative integers $`k,l`$ they fulfill the following relations:

(1)
$$\alpha _{k,l}\alpha _{l,k+1}=\alpha _{l,k}\alpha _{k,l+1},$$

(2)
$$\alpha _{k,l}^2+\alpha _{l1,k}^2=\alpha _{l,k}^2+\alpha _{k1,l}^2.$$

Using the notation of this theorem let us denote

$$\psi _{k,l}=P_{k,l}(\mathrm{\Phi },\mathrm{\Phi }^{})\mathrm{\Omega }.$$

We have

(3)
$$\psi _{k,l}|\psi _{m,n}=_𝐂\overline{P_{k,l}}P_{m,n}𝑑\mu (z)=\delta _{k,m}\delta _{l,n},$$

(4)
$$\mathrm{\Phi }\psi _{k,l}=\alpha _{k,l}\psi _{k+1,l}+\alpha _{l1,k}\psi _{k,l1},$$

(5)
$$\mathrm{\Phi }^{}\psi _{k,l}=\alpha _{l,k}\psi _{k,l+1}+\alpha _{k1,l}\psi _{k1,l}.$$

We see that $`\psi _{k,l}`$ is a family of orthonormal vectors in $``$ indexed by a pair nonnegative numbers. We can think that the index $`k`$ represents the number of particles and the index $`l`$ the number of antiparticles which are in the physical state corresponding to the operator $`\mathrm{\Phi }`$.

Equations (4) and (5) tell us that the operator $`\mathrm{\Phi }=K^{}+\mathrm{\Lambda }`$ is a sum of two operators: $`K^{}`$—a creator of a particle and $`\mathrm{\Lambda }`$—an annihilator of an antiparticle, while $`\mathrm{\Phi }^{}=K+\mathrm{\Lambda }^{}`$ is a sum of a annihilator of a particle and a creator of an antiparticle.

## 3.  FACTORIZABLE CREATION AND ANNIHILATION OPERATORS

It natural to restrict our considerations to the case when each of the operators $`K,K^{},\mathrm{\Lambda },\mathrm{\Lambda }^{}`$ is a product of two operators, each depending either on the number of particles or the number of antiparticles.

This means that we are solving the equation $`\alpha _{k,l}=f_kg_l`$. From the equation (1) it follows that

$$g_lg_{k+1}=g_kg_{l+1},$$

so $`(g_i)`$ form a geometric sequence. We may take it in the form $`g_i=q^i`$ since the constant factor can be moved to the $`(f_i)`$ sequence. From the assumption $`\alpha _{k,l}>0`$ follows that $`q>0`$.

Now the equation (2) takes the form

$$[(f_k)^2(f_{k1})^2]q^{2k}=[(f_l)^2(f_{l1})^2]q^{2l},$$

therefore there exists a constant $`c`$ such that

$$[(f_k)^2(f_{k1})^2]q^{2k}=c^2$$

from which follows that

$$\alpha _{k,l}=c\sqrt{1+q^2+\mathrm{}+q^{2k}}q^l=\{\begin{array}{cc}c\sqrt{\frac{1q^{2k+2}}{1q^2}}q^l& \text{for }q1\hfill \\ & \\ c\sqrt{k+1}& \text{for }q=1\hfill \end{array}$$

We see that the case of $`q=1`$ corresponds to $`\mathrm{\Phi }`$ being equal (up to a constant) to a sum of a classical bosonic particle creator and antiparticle annihilator.

In the general case we have $`K=cA_kq^{N_l}`$ and $`\mathrm{\Lambda }=cA_lq^{N_k}`$ where $`N_k`$, $`N_l`$ denote the number of particles or number antiparticles operators respectively and $`A_k`$, $`A_l`$ denote the particle and antiparticls annihilation operators, $`A_k^{}`$, $`A_l^{}`$ denote the particle and antiparticle creation operators.

Operators concerning particles and antiparticles commute:

$$A_kA_l=A_lA_k,$$

$$A_k^{}A_l^{}=A_l^{}A_k^{},$$

$$A_k^{}A_l=A_lA_k^{}.$$

Creation and annihilation operators fulfill the $`q^2`$-deformed commutation relation :

$$A_kA_k^{}q^2A_k^{}A_k=1$$

$$A_lA_l^{}q^2A_l^{}A_l=1$$


text:
# 1 INTRODUCTION

## 1 INTRODUCTION

More and more supporting observations (Sahai & Trauger 1998; Kwok, Su, & Hrivnak 1998; Hrivnak, Kwok, & Su 1999; Kwok, Hrivnak, & Su 2000; Huggins et al. 2001) and theoretical considerations (Soker 1990; Soker 2001) are accumulated in support of the view that significant shaping of the circumstellar material takes place just before, after, and during the transition from the asymptotic giant branch (AGB) to the planetary nebula (PN) phase. Both the wind and radiation properties are significantly changed during these stages. As the star is about to leave the AGB the mass loss rate increases substantially, up to $`10^510^4M_{}\mathrm{yr}^1`$. This wind was termed superwind by Renzini (1981); because of confusion with other winds termed superwinds, hereafter we’ll refer to this wind as FIW, for final intensive wind. After the star leaves the AGB the mass loss rate decreases down to $`10^8M_{}\mathrm{yr}^1`$, and its velocity increases from $`10\mathrm{km}\mathrm{s}^1`$ to few$`\times 10^3\mathrm{km}\mathrm{s}^1`$ at the PN phase. Simultaneously, the effective temperature increases, and the post-AGB star starts to ionize the nebula around it, when by definition the PN phase starts. The changes in the mass loss rate and velocity of the wind are accompanied by a change in the wind geometry, e.g., jets and bipolar structures are formed (Soker 1990; Sahai & Trauger 1998; Kwok et al. 2000). The change in the wind geometry may result from an intrinsic processes in the AGB and post-AGB mass-losing star, or from processes in an accreting companion (e.g., Soker 2001). The increase in the wind velocity leads to collision of winds (Kwok, Purton, & Fitzgerald 1978), which can lead to instabilities during the post-AGB (proto-PN) stage (Dwarkadas & Balick 1998). As the central star starts ionizing the nebula, an ionization front propagates outward, and plays a significant role in shaping the nebula, both in the radial direction (e.g., Mellema & Frank 1995; Chevalier 1997; Schönberner & Steffen 2000), and in the transverse directions (e.g., Mellema 1995; Soker 2000b).

In the present paper we examine yet another process which may occur during the transition from the AGB to the PN stage, i.e., after the FIW ceases and before ionization starts. This is a backflow of a fraction of the dense wind toward the central star(s). This backflowing material may be accreted by the central star and/or its companion. The idea of accreting backflowing material during the post-AGB phase was raised before to explain and account for: a possible mechanism for the formation of jets from an accretion disk (Bujarrabal, Alcolea, & Neri 1998); a slower evolution along the post-AGB track (Zijlstra et a;. 2000); and post-AGB stars depleted of refractory elements which compose the dust particles (e.g., Van Winckel et al. 1998) by accretion of a dust-depleted circumstellar gas (Waters, Trams & Waelkens 1992), most likely in binary systems (Van Winckel 1999).

The goal of the present paper is to explore the conditions required for a backflow to occur such that it plays a non-negligible role in the post-AGB evolution. The conditions are derived in $`\mathrm{\S }2`$, while the implications for the processes mentioned above, as well as other processes, together with a short summary, are in $`\mathrm{\S }3`$.

## 2 CONDITIONS FOR A BACKFLOW

In this section we derive the conditions for a backflow to occur during the post-AGB phase. We assume that a very slow flow exists along some directions, e.g., in the equatorial plane, and consider the conditions for some of this material to flow inward and be accreted by the central system, before ionization starts. After being ionized, any dense cool gas will expand and will be pushed outward by radiation and ram pressure (see below). We do not consider the deceleration of the slowly outward moving mass element, but simply assume that if there is a slowly outward moving mass element, it will reach zero radial velocity at some radius.

We therefore consider a mass-element $`M_i`$ with a constant density $`\rho _i`$ within a solid angle $`\mathrm{\Omega }`$ and a radial extention $`\mathrm{\Delta }r`$ at a distance $`r\mathrm{\Delta }r`$ from the central star, such that

$$M_i=\rho _i\mathrm{\Omega }r^2\mathrm{\Delta }r.$$
(1)
The mass element can also be a shell where $`\mathrm{\Omega }=4\pi `$. If the sound crossing time $`\mathrm{\Delta }r/c_s`$, where $`c_s`$ is the sound speed, is shorter than any other time scale in the process, we can take the shell to move more or less coherently. Three relevant forces are acting on the mass element in the radial direction. The gravitational force of the central star(s)

$$f_g=\frac{GMM_i}{r^2},$$
(2)
where $`M`$ is the total mass of the central system, a binary system or a single star. The fast wind blown by the central star and its radiation push outward. The force due to the fast wind’s ram pressure, assuming it is much faster than the slow wind velocity, is given by

$$f_w=\rho _w\left(r\right)v_w^2\mathrm{\Omega }r^2=\dot{M}_wv_w\beta ,$$
(3)
where $`\rho _w\left(r\right)=\dot{M}_w/\left(4\pi r^2v_w\right)`$, $`v_w`$ and $`\dot{M}_w`$ are the density, velocity and mass loss rate (defined positively) of the fast wind, and $`\beta \mathrm{\Omega }/4\pi `$. The radiation imparts a force of

$$f_r=\frac{L_{}}{c}\beta \left(1e^\chi \right),$$
(4)
where $`L_{}`$ is the luminosity of the central system, $`c`$ the speed of light, and

$$\chi =\rho _i\kappa \mathrm{\Delta }r=7\left(\frac{M_i}{0.01M_{}}\right)\left(\frac{\kappa }{10\mathrm{cm}^2\mathrm{g}^1}\right)\left(\frac{r}{100\mathrm{AU}}\right)^2\beta ^1$$
(5)
is the optical depth of the mass element, and we used equation (1) for the density. We scale the opacity $`\kappa \left(r\right)`$ with a typical value for AGB stars (Jura 1986; Winters et al. 2000), and assume that the fast wind inner to the mass element absorbs a negligible fraction of the radiation. For convenience we define three dimensionless variables. The ratio of maximum radiation pressure to the fast wind ram pressure

$$q\frac{L_{}/c}{\dot{M}_wv_w}=\left(\frac{L_{}}{5000L_{}}\right)\left(\frac{\dot{M}_w}{10^6M_{}\mathrm{yr}^1}\right)^1\left(\frac{v_w}{100\mathrm{km}\mathrm{s}^1}\right)^1,$$
(6)
which does not depend on $`r`$. We scale the fast wind mass loss and velocity as appropriate for a post-AGB star before it starts ionizing the nebula. The ratio of the force due to radiation and wind to that of gravity depends on $`r`$, both through the dependence of optical depth $`\chi `$ and gravity on $`r`$. We define $`\eta `$ to be this ratio at a scaling radius $`r_0`$

$$\eta =\frac{\dot{M}_wv_w}{GMM_i}k\left(r_0\right)r_0^2,$$
(7)
where

$$k\left(r_0\right)\beta \left[1+q\left(1e^\chi \right)\right].$$
(8)
The equation of motion for the mass element can be written as

$$\frac{d^2r}{dt^2}=\frac{GM}{r_0^2}\left(\eta \frac{r_0^2}{r^2}\right).$$
(9)
We take $`r_0`$ to be the radius at which the radial velocity of the mass element is zero. When there is only gravity, the free-fall time from $`r=r_0`$ to the center $`r=0`$, with $`v\left(r_0\right)=0`$, is

$$t_{ff}=\frac{\pi }{2^{3/2}}\frac{r_0^{3/2}}{\left(GM\right)^{1/2}}=1400\left(\frac{r_0}{400\mathrm{AU}}\right)^{3/2}\left(\frac{M}{1M_{}}\right)^{1/2}\mathrm{yr}.$$
(10)
We define the dimensionless variables

$$\tau =t/t_{ff}\mathrm{and}x=r/r_0,$$
(11)
and write the equation of motion (9) in the form

$$\frac{d^2x}{d\tau ^2}=\frac{\pi ^2}{8}\left(\eta x^2\right).$$
(12)
Assuming that $`\eta `$ is constant and does not depend on $`x`$ allows us to integrate once the last equation to give the velocity

$$v\frac{dx}{d\tau }=\frac{\pi }{2}\left(\eta x+x^11\eta \right)^{1/2},$$
(13)
where we substituted the initial condition $`v\left(1\right)=0`$. This can be integrated analytically for constant values of $`\eta =1`$ and $`\eta =0`$. The case $`\eta =0`$ gives the free fall solution. The time left for the object to fall from $`x`$ to $`x=0`$ is given by

$$\tau _f\left(x\right)=\frac{2}{\pi }\left(\mathrm{sin}^1x^{1/2}\left[x\left(1x\right)\right]^{1/2}\right).$$
(14)
As expected from our scaling $`\tau \left(1\right)=1`$, i.e., the free fall time from $`r=r_0`$ to $`r=0`$ is $`t_{ff}`$. For $`\eta =1`$ and with $`v\left(1\right)=0`$, the fall time from $`x=1`$ is infinite. This is because the wind and radiation outward-acceleration equals the gravity inward-acceleration. However, the time left to fall from $`x`$ to $`x=0`$, with $`v\left(1\right)=0`$, is finite, and is given by

$$\tau _f\left(x\right)=\frac{2}{\pi }\left(\mathrm{ln}\left[\left(1+x^{1/2}\right)/\left(1x^{1/2}\right)^1\right]2x^{1/2}\right).$$
(15)
Comparing the last two equations, we find as expected, that radiation pressure and wind’s ram pressure slow down the inflow. For example, for $`\eta =0`$, i.e., no wind and radiation pressure, the time left to fall from $`x=0.9`$ to $`x=0`$ is $`\tau =0.604`$, while for $`\eta =1`$ it is $`\tau =1.11`$. In both cases the initial condition is $`v\left(1\right)=0`$. These cases are less interesting, since for $`\eta =1`$ the flow actually stagnates at $`x=1`$. More interesting is the case of $`v\left(1\right)=0`$ and $`0<\eta <1`$, since for $`\eta >1`$ it will be push away from the center. We numerically integrated equation (13) for these conditions and for constant values of $`\eta `$. We find the backflow time to be $`\tau _f\left(1\right)=1.28`$, $`1.52`$, and $`2.08`$, for $`\eta =0.5`$, $`0.7`$, and $`0.9`$, respectively. For $`\eta =0.8`$, for example, the fall back times from $`x_i=1`$, $`0.9`$ and $`0.8`$ are $`\tau _f=1.72`$, $`1.23`$ and $`0.925`$, respectively, where $`v\left(x_i\right)=0`$ in these cases. For $`\eta =0`$ the backflow times for the same initial conditions are $`1`$, $`0.85`$, and $`0.72`$, respectively. The conclusion from the numerical values cited above is that the typical backflow time from $`rr_0`$ is the free fall time $`t_{ff}`$ at $`r_0`$, but because of the radiation and wind pressures the region from which this is the backflow time is much larger than the $`\eta =0`$ cases, extending from $`r_0`$ down to $`0.70.8r_0`$ for $`\eta 0.8`$. So the question is what is the value of $`r_0`$ for which $`\eta =1`$. Below and close to this radius the gas falls back in a time $`t_{ff}`$, while it is accelerated away for larger radii. From equation (7) we find

$$r_0=\left(\frac{\eta GMM_i}{k\dot{M}_wv_w}\right)^{1/2}=430\mathrm{AU}\left[\left(\frac{\mathrm{M}}{1\mathrm{M}_{}}\right)\left(\frac{\mathrm{M}_\mathrm{i}}{0.01\mathrm{M}_{}}\right)\left(\frac{\dot{\mathrm{M}}_\mathrm{w}}{10^6\mathrm{M}_{}}\right)^1\left(\frac{v_w}{100\mathrm{km}\mathrm{s}^1}\right)^1\left(\frac{\mathrm{k}}{0.1}\right)^1\eta \right]^{1/2}$$
(16)
For the backflowing mass to influence the evolution significantly, we required the fall back time to be $`1000\mathrm{yr}`$. For the typical parameters used in equations (5) and (6) we find from equation (8) that $`k\left(r\right)2\beta `$; using the time scale given by equation (10) in equation (16) gives the desired condition

$$\frac{M_i}{\beta }0.1M_{}$$
(17)
where as before, $`\beta \mathrm{\Omega }/4\pi `$, and $`\mathrm{\Omega }`$ is the solid angle covered by the dense backflowing material. The last condition is limited by a maximum density, since a large value of the backflowing mass $`M_i`$ and small value for $`\beta `$ means a very high density. We now estimate a reasonable value for the density. We assume a very slow equatorial flow, with a speed of $`v_s1\mathrm{km}\mathrm{s}^1`$ and with a mass loss rate per unit solid angle of $`\dot{m}_s=\dot{M}_s/4\pi `$. The density of the mass elements formed by this wind is

$$\rho _{iw}=\frac{\dot{m}_s}{r^2v_s}=1.4\times 10^{15}\left(\frac{\dot{M}_s}{10^3M_{}\mathrm{yr}^1}\right)\left(\frac{r}{400\mathrm{AU}}\right)^2\left(\frac{v_s}{1\mathrm{km}\mathrm{s}^1}\right)^1\mathrm{g}\mathrm{cm}^3.$$
(18)
The minimum density is that for which the fast wind compresses the dense cool wind such that the ram pressure $`\rho _wv_w^2`$ equals the thermal pressure of the cool gas. For a molecular gas, we find this density to be

$$\rho _{ip}=\rho _w\frac{v_w^2}{c_i^2}=10^{16}\left(\frac{\dot{M}_w}{10^6M_{}}\right)\left(\frac{v_w}{100\mathrm{km}\mathrm{s}^1}\right)\left(\frac{r}{400\mathrm{AU}}\right)^2\left(\frac{T_i}{300\mathrm{K}}\right)^1\mathrm{g}\mathrm{cm}^3.$$
(19)
where $`c_i`$ and $`T_i`$ are the isothermal sound speed and temperature, respectively, of the cool gas. We find that a density of $`\rho _i10^{15}\mathrm{g}\mathrm{cm}^3`$ is reasonable. Using equation (1) and the definition of $`\beta `$ gives

$$\frac{M_i}{\beta }=0.14\left(\frac{\rho _i}{10^{15}\mathrm{g}\mathrm{cm}^3}\right)\left(\frac{r}{400\mathrm{AU}}\right)^3\left(\frac{\mathrm{\Delta }r}{0.1r}\right)M_{}.$$
(20)
Condition (17) is met for the density given by equation (18), but this requires a very high mass loss rate per unit solid angle. If $`\beta =0.1`$ this requires a total mass loss rate of $`\dot{M}=10^4M_{}\mathrm{yr}^1`$, but concentrated in particular directions, probably in the equatorial plane. All these considerations strongly suggest an equatorial dense and slow flow, such as expected in a close binary system (Mastrodemos & Morris 1999; Soker 2000a). A very fast rotation can also form such a wind (Bjorkman & Cassinelli 1993), but then a binary companion is needed to substantially spin-up the envelope.

## 3 IMPLICATIONS AND SUMMARY

Despite the assumptions and simplifications in deriving condition (17), we feel that the results obtained in the previous section and the implications discussed in this section are quite general. We find that for a backflow to occur on a time scale of $`t_{\mathrm{acc}}10^3\mathrm{yr}`$ after the termination of the AGB, so that it has a non-negligible role in the post-AGB evolution, the following conditions should be met: (1) The total backflowing mass should be larger than the combined mass lost in the wind and that burned in the core. For a post-AGB mass loss rate of $`10^6M_{}\mathrm{yr}^1`$ the nuclear burning is negligible, and the total required mass is $`M_{\mathrm{acc}}10^3(t_{\mathrm{acc}}/1,000\mathrm{yr})\mathrm{M}_{}`$. (2) The backflowing mass should have a very low, $`1\mathrm{km}\mathrm{s}^1`$, terminal velocity, so that eventually it will be decelerated to zero velocity, and flow back. (3) Condition (17) on the ratio of the mass of the mass-element and the solid angle it covers $`\beta =\mathrm{\Omega }/4\pi `$ should be met. (4) For reasonable densities (eq. 18), and the required mass, we find (eq. 20) that the mass range is $`M_i0.110^3M_{}`$, and the appropriate solid angle covered by the backflowing mass is $`1<\beta 10^2`$. We can take the typical values to be $`M_i0.01M_{}`$ and $`\beta =0.1`$.

These mass loss properties required that ($`i`$) the flow be concentrated along particular directions, and ($`ii`$) have an inefficient acceleration by the stellar radiation. In a previous paper (Soker 2000a) two mechanisms which lead to such a flow were discussed. In the first mechanism proposed by Soker (2000a) magnetic cool spots (as in the Sun) are formed on the surface of slowly rotating AGB stars. The lower temperature above the spots enhances dust formation, which during the final intensive wind (FIW) may lead to an optically thick wind. This means an inefficient radiative acceleration, hence a slow flow above the spot. If the spot is small, material from the surrounding flows into the shaded region and accelerates the slow flow. However, if the spot is large, material from the surroundings will not accelerate the flow much, and it will stay slow. In the present paper we find that the condition for the slowly moving material to flow back is that the spot has $`\beta 0.01`$, which for a circular spot means a radius of $`R_{\mathrm{spot}}0.2R_{}`$, where $`R_{}`$ is the stellar radius. The required mass in the slow flow is given by equation (17).

In the second mechanism a high density in the equatorial plane is formed by a binary interaction, where the secondary star is close to, but outside the AGB envelope. This mechanism is supported by observations, e.g., slowly moving equatorial gas is found around several binary post-AGB stars (Van Winckel 1999, and references therein). In the process proposed by Soker (2000a) the strong interaction between the two stars forms a dense equatorial outflow which is optically thick, leading to an inefficient radiative acceleration and a very slow equatorial flow. For a very massive and significant backflow, with a mass of $`M_i0.01M_{}`$ and $`\beta 0.1`$, to occur, a binary mechanism is required. For the cool spots model to form such a massive slow flow, several large spots are required. This means a strong magnetic activity, which probably requires the AGB star to be spun-up by a stellar companion. We therefore argue that in both mechanisms a binary companion is required to cause a massive flow, such that it may last for $`t_{\mathrm{acc}}10^3\mathrm{yrs}`$, possibly by as long as $`10^4\mathrm{yrs}`$ in extreme cases.

Such a backflow may have the following effects. If it has a large specific angular momentum, as expected in strongly interacting binary systems, the backflowing material may form an accretion disk around one or two of the two stars. The disk(s) may blow jets or collimated fast winds (CFWs), which will play a significant role in shaping the circumstellar material. Such a possibility was briefly mentioned by Bujarrabal et al. (1998) for the proto-PN M1-92. The accreted mass may slow down the post-AGB evolution, as suggested by Zijlstra et al. (2000) for some OH/IR stars. Zijlstra et al. (2000) termed these stars retarded stars, and argue for a delay as long as $`10^4\mathrm{yrs}`$ by accretion from a near-stationary reservoir. They bring supporting observations, but don’t consider the formation of such a reservoir of mass. The results of the present paper put the idea of Zijlstra et al. (2000) on a more solid ground. Another effect attributed to backflow accretion is the formation of post-AGB stars depleted of refractory elements which compose the dust particles (e.g., Van Winckel et al. 1998; Waters, Trams & Waelkens 1992; Van Winckel 1999). The separation of dust from the gas was not considered in the present paper.

Finally, we speculate on another plausible effect of the dense backflowing gas. The central star’s wind is shocked when it hits the dense material. If some dense backflowing blobs survive long into the PN phase, when the central star’s wind velocity is $`10^3\mathrm{km}\mathrm{s}^1`$, then there will be a hard X-ray emission from the post-shock fast wind material. The dense blob will be close to the central star, making the hard X-ray emitting region hard to resolve. It is not clear if this compact hard X-ray emitting region can explain the recent Chandra observations of a “point source” in the centers of the Helix (NGC 7293) and Cat’s Eye (NGC 6543) PNe (Guerrero et al. 2000).

ACKNOWLEDGMENTS: I thank Benzion Kon for his help with the mathematics. This research was supported by a grant from the US-Israel Binational Science Foundation.