format: text

text:
# Strong field limit of black hole gravitational lensing

## I Introduction

The simplicity of the theory of gravitational lensing derives from some basic assumptions which are satisfied by most physical situations . Through the weak field and the thin lens approximations, the whole theory was built and succeeded in explaining the phenomenology risen up to now.

Yet, gravitational lensing must not be conceived as a weak field phenomenon, since high bending and looping of light rays in strong fields is one of the most well-known and amazing predictions of general relativity. The importance of gravitational lensing in strong fields is highlighted by the possibility of testing the full general relativity in a regime where the differences with non-standard theories would be manifest, helping the discrimination among the various theories of gravitation . However, the complexity of the full mathematical treatment and the difficulties of experimental or observational evidences are obstacles for these studies.

Without the linear approximations, we have trascendent equations which are hard to be handled even numerically.

Several studies about light rays close to the Schwarzschild horizon has been lead: for example, Viergutz  made a semi–analytical investigation about geodesics in Kerr geometry; in refs.  the appearance of a black hole in front of a uniform background was studied. Recently, Virbhadra & Ellis  faced the simplest strong field problem, represented by deflection in Schwarzschild space–time, by numerical techniques. The existence of an infinite set of relativistic images has been enlightened and the results have been applied to the black hole at the centre of the Galaxy. Later on, by an alternative formulation of the problem, Frittelli, Kling & Newman  attained an exact lens equation, giving integral expressions for its solutions, and compared their results to those by Virbhadra & Ellis.

All these studies are affected by the high complexity in the investigation of the geodesics in the neighbourhood of the horizon and must resort to numerical methods to return valuable results. This issue prevents any kind of general and systematic investigation of light bending in this region, since fair analytical formulae for the interesting quantities are absolutely missing.

Starting from the black hole lens equation of ref. , we perform a set of expansions exploiting the source–lens–observer geometry and the properties of highly deflected light rays. In this way, we manage to solve the lens equation and find analytical expressions for the infinite set of images formed by the black hole. This approach leads to extremely simple formulae which allow an immediate comprehension of the problem and a straightforward application to the physically interesting situations. This strong field approach can be surely considered as the direct counterpart of the weak field limit for its striking simplicity and reliability. On the other hand, this approach allows to investigate where the strong field limit effects become relevant.

This paper is structured as follows: In Sect. 2, the black hole lens equation is presented and the first basic approximations are explained. Sect. 3 contains the calculation of the deflection angle. In Sect. 4 the deflection angle is then plugged in the black hole lens equation to derive the relativistic images, their amplification and the critical curves. Sect. 5 contains a discussion of these results. Conclusions are drawn in Sect. 6.

## II The black hole lens equation

The geometrical configuration of gravitational lensing is shown in Fig. 1. The light emitted by the source $`S`$ is deviated by the black hole $`L`$ and reaches the observer $`O`$. Here, as a black hole, we mean any compact object having a radius comparable to its Schwarzschild radius, so that even very compact object which has not undergone a full gravitational collapse would work in the same way; $`\beta `$ is the angular position of the source with respect to the optical axis $`OL`$ and $`\theta `$ is the angular position of the image seen by the observer. It is important to stress that the closest approach distance $`x_0`$ does not coincide with the impact parameter $`b`$, unless in the limit of vanishing deflection angle $`\alpha `$.

By inspection of Fig. 1, it is possible to write a relation among the source position, the image position and the deflection angle $`\alpha `$.

$$\mathrm{tan}\beta =\mathrm{tan}\theta \frac{D_{LS}}{D_{OS}}\left[\mathrm{tan}\theta +\mathrm{tan}\left(\alpha \theta \right)\right].$$
(1)
This is what is called the full lens equation. Given a source position $`\beta `$, the values of $`\theta `$, solving this equation, give the position of the images observed by $`O`$.

In the weak field limit, several standard approximations are performed. The tangents are expanded to the first order in the angles since they are, at most, of the order of arcsec. The weak field assumption reduces the deflection angle to $`{\displaystyle \frac{4GM}{c^2x_0}}`$. Then the lens equation can be solved exactly and two images are found: one on the same side of the source and one on the opposite. Their separations from the optical axis are of the order of the Einstein angle

$$\theta _E=\sqrt{\frac{4GM}{c^2}\frac{D_{LS}}{D_{OS}D_{OL}}}.$$
(2)
It is easy to see that, in most relevant cases, these images are formed by light rays passing very far from the event horizon, justifying the weak field approximation .

Now we turn to the study of gravitational lensing in strong field. From now on, all lengths will be expressed in units of the Schwarzschild radius $`{\displaystyle \frac{2GM}{c^2}}`$. The deflection angle $`\alpha `$ contains the physical information about the deflector and must be calculated through the integration of the geodesic of the light ray . Its integral expression, as a function of $`x_0`$, is

$$\alpha \left(x_0\right)=\underset{x_0}{\overset{\mathrm{}}{}}\frac{2}{x\sqrt{\left(\frac{x}{x_0}\right)^2\left(1\frac{1}{x_0}\right)\left(1\frac{1}{x}\right)}}𝑑x\pi .$$
(3)
The next section is completely devoted to the achievement of a manageable expression for this quantity as a function of $`\theta `$.

When the light ray trajectory gets closer to the event horizon, the deflection increases. At some impact parameter, $`\alpha `$ becomes higher than $`2\pi `$, resulting in a complete loop of the light ray around the black hole. Decreasing further the impact parameter, the light ray winds several times before emerging. Finally, for $`b=(3\sqrt{3})/2`$, corresponding to $`x_0=3/2`$, the deflection angle diverges and the light ray is captured by the black hole. For each loop we add to the light ray geodesic, there is one particular value of the impact parameter such that the observer is reached by the light coming from the source. So there will be an infinite sequence of images on each side of the lens.

We shall put our attention on situations where the source is almost perfectly aligned with the lens. In fact, this is the case where the relativistic images are most prominent. In this case, we are allowed to expand $`\mathrm{tan}\beta `$ and $`\mathrm{tan}\theta `$ to the first order. Some more words are needed for the term $`\mathrm{tan}\left(\alpha \theta \right)`$. Even if $`\theta `$ is small, $`\alpha `$ is not small in the situations of our interest. However, if a ray of light emitted by the source $`S`$ is going to reach the observer after turning around the black hole, $`\alpha `$ must be very close to a multiple of $`2\pi `$. Writing $`\alpha =2n\pi +\mathrm{\Delta }\alpha _n`$, with $`n`$ integer, we can perform the expansion $`\mathrm{tan}\left(\alpha \theta \right)\mathrm{\Delta }\alpha _n\theta `$. The lens equation becomes

$$\beta =\theta \frac{D_{LS}}{D_{OS}}\mathrm{\Delta }\alpha _n.$$
(4)
As defined by (3), $`\alpha `$ is a positive real number, corresponding to a clockwise winding in Fig. 1. We have inserted it in (4) without giving a sign. Taking a positive $`\beta `$, this equation describes only images on the same side of the source ($`\theta >0`$). To obtain the images on the opposite side, we can solve the same equation with the source placed in $`\beta `$. Taking the opposite of these solutions, we obtain the full set of the secondary images.

## III The deflection angle

The deflection angle can be evaluated exactly, but its expression does not allow the resolution of (4). Anyway, it is possible to make some simplifying but very general approximations which reduce the deflection angle to an expression easier to handle. By this strategy, we shall calculate this fundamental quantity. All approximations are essentially based on the proximity of the closest approach distance $`x_0`$ to its minimum value which is $`3/2`$.

The integral in Eq.(3) gives the following result

$$\alpha =\pi 4F(\varphi _0,\lambda )G\left(x_0\right)$$
(5)
where

$$G\left(x_0\right)=\sqrt{\frac{x_0\left(3+3x_0\sqrt{3+2x_0+x_0^2}\right)}{\left(32x_0\right)\left(1x_0+\sqrt{3+2x_0+x_0^2}\right)}}$$
(6)
and

$$F(\varphi _0,\lambda )=\underset{0}{\overset{\varphi _0}{}}\left(1\lambda \mathrm{sin}^2\varphi \right)^{1/2}𝑑\varphi $$
(7)
is an elliptic integral of first kind. The parameters $`\varphi _0`$ and $`\lambda `$ themselves are functions of $`x_0`$:

$$\varphi _0=\mathrm{arcsin}\sqrt{\frac{3+x_0\sqrt{3+2x_0+x_0^2}}{2\left(3+2x_0\right)}},$$
(8)
and

$$\lambda =\frac{3x_0\sqrt{3+2x_0+x_0^2}}{3x_0+\sqrt{3+2x_0+x_0^2}}.$$
(9)
Fig. 2 shows a plot of the deflection angle. This expression is, of course, too complicated to use in Eq. (4). However, we can see that $`\alpha `$ diverges when $`x_03/2`$. On the other hand, we are interested just into small closest approaches, since they correspond to the high deflection angles producing relativistic images. If we let $`x_0=3/2+ϵ`$, we can search for the leading order term in the divergence when $`ϵ0`$.

After some non–trivial expansions, we find that the leading order of the deflection angle is logarithmic in $`ϵ`$, that is

$$\alpha 2\mathrm{log}\frac{\left(2+\sqrt{3}\right)ϵ}{18}\pi .$$
(10)
In Fig. 3, we plot the ratio of the exact deflection angle and its leading order expansion. Of course, the two functions coincide in the limit $`x_01.5`$. Using our expansion up to $`x_0=1.55`$, the error we commit is about $`1\%`$. For our purposes, this is largely sufficient as we shall see later.

The next step is to convert the dependence on $`x_0`$ into a dependence on the impact parameter $`b`$ and then on $`\theta `$. The study of light ray dynamics in Schwarzschild metric gives the relation

$$x_0^2=\left(1\frac{1}{x_0}\right)b^2.$$
(11)
For values of $`x_0`$ close to $`3/2`$, this relation can be solved perturbatively (or, equivalently, one can expand the exact root of the third degree polynomial equation in $`x_0`$ to the first non-trivial order), giving the relation

$$ϵ=\sqrt{\frac{b\frac{3\sqrt{3}}{2}}{\sqrt{3}}}.$$
(12)
Remembering that $`b\theta D_{OL}`$, we can substitute this relation into Eq. (10) to retrieve $`\alpha `$ as a function of $`\theta `$, that is

$$\alpha \mathrm{log}\left(\theta D_{OL}\frac{3\sqrt{3}}{2}\right)+A,$$
(13)
with

$$A=\mathrm{log}\frac{\left(5+3\sqrt{3}\right)}{1944}\pi =2.109.$$
(14)
Eq. (13) is the highly simplified expression for the deflection angle we were looking for. It allows reliable calculations, as we shall see in the next section, where we insert it in the lens equation.

## IV Images, magnification, critical curves

According to the considerations done at the end of Sect. 2, what enters the lens equation is not the full deflection angle but its deviation from a multiple of $`2\pi `$. First of all, we have to find the values of $`\theta `$ (denoted by $`\theta _n^0`$) such that

$$\alpha \left(\theta _n^0\right)=2n\pi .$$
(15)
Inverting this equation, with $`\alpha `$ given by (13), we have

$$\theta _n^0=\frac{3\sqrt{3}+2e^{A2n\pi }}{2D_{OL}}.$$
(16)
As expected, if we let the number of loops around the black hole tend to infinity, we have $`\theta _{\mathrm{}}^0={\displaystyle \frac{3\sqrt{3}}{2D_{OL}}}`$. $`\theta _n^0`$ are the starting points for all the successive calculations. The offsets $`\mathrm{\Delta }\alpha _n`$ can be found by expanding $`\alpha `$ to the first order in $`\mathrm{\Delta }\theta _n=\theta \theta _n^0`$:

$$\mathrm{\Delta }\alpha _n=\frac{D_{OL}}{e^{A2n\pi }}\mathrm{\Delta }\theta _n.$$
(17)
We have all the ingredients to solve the lens equation (4), which becomes

$$\beta =\left(\theta _n^0+\mathrm{\Delta }\theta _n\right)+\left(\frac{D_{OL}}{e^{A2n\pi }}\frac{D_{LS}}{D_{OS}}\right)\mathrm{\Delta }\theta _n.$$
(18)
We can observe that $`D_{OL}`$, $`D_{LS}`$, $`D_{OS}`$ are all much greater than unity (remember that all distances are measured in Schwarzschild radii) and $`e^{A2\pi }=3.25\times 10^3`$. This means that the last term in Eq.(18) prevails on the $`\mathrm{\Delta }\theta _n`$ at the second place in the rhs. Neglecting this term, we finally get the position of the $`n^{th}`$ image as

$$\theta _n\theta _n^0+\frac{e^{A2n\pi }\left(\beta \theta _n^0\right)D_{OS}}{D_{LS}D_{OL}}.$$
(19)
We see that, when $`\beta `$ equals $`\theta _n^0`$, there is no correction to the position of the $`n^{th}`$ image, that remains in $`\theta _n^0`$ simply. In this particular case, the image position coincides with the source position. It is worthwhile to note that the second term in (19) is much smaller than the first one. For practical purposes, $`\theta _n^0`$ are already a good approximation for the position of relativistic images.

The relation (17) can help us to estimate the error in the determination of the position of the images. We do this in the least favorable case that is for the first image, which is the most external and thus the farthest from the divergence we started from. The relative error for $`\theta _1^0`$ is

$$\frac{\mathrm{\Delta }\theta _1}{\theta _1^0}=\left(\frac{e^{A2\pi }}{D_{OL}}\frac{\alpha }{\theta _1^0}\right)\frac{\mathrm{\Delta }\alpha }{\alpha }$$
(20)
For the first image we have $`\alpha 2\pi `$. In the previous section, we have estimated the relative error on $`\alpha `$ to be about one percent for $`x_01.55`$. By using $`b=\theta D_{OL}`$ and Eq. (12), we can find the closest approach distance for the first image to be $`1.545`$. In conclusion, this estimate for $`\mathrm{\Delta }\alpha /\alpha `$ can be reasonably used.

The error on $`\theta _1^0`$ is

$$\frac{\mathrm{\Delta }\theta _1}{\theta _1^0}\frac{4\pi e^{A2\pi }}{3\sqrt{3}+2e^{A2\pi }}\frac{\mathrm{\Delta }\alpha }{\alpha }8\times 10^5.$$
(21)
For the other images, the errors can be found to be orders of magnitude smaller. This estimate supports our results which thus prove to be extremely reliable.

The critical curves are defined as the points where the Jacobian determinant of the lens equation vanishes, that is

$$\frac{\beta }{\theta }\frac{\beta }{\theta }=0.$$
(22)
The term $`{\displaystyle \frac{\beta }{\theta }}`$ is always positive (there are no radial critical curves). Then the (tangential) critical curves are obtained for $`\beta =0`$. As we have already solved the lens equation for each $`\beta `$, it is sufficient to put $`\beta =0`$ in Eq.(19) to get the angles of these relativistic curves, that is

$$\theta _n^{cr}\theta _n^0\left(1\frac{e^{A2n\pi }D_{OS}}{D_{LS}D_{OL}}\right).$$
(23)
A source perfectly aligned to the black hole produces an infinite series of concentric rings with these radii.

The magnification of the images (19) is nothing else but the inverse of the modulus of the Jacobian determinant already used for the critical curves. In the region of our interest, i.e. for small source angles, it is easy to calculate all the needed quantities from the lens equation. We do this by approximating the images by the angles $`\theta _n^0`$:

$$\frac{\beta }{\theta }|_{\theta _n^0}=1+\frac{D_{OL}}{e^{A2n\pi }}\frac{D_{LS}}{D_{OS}}.$$
(24)
The second term is much higher than $`1`$ for the same reasons expressed before. Then the magnification of the $`n^{th}`$ image is

$$\mu _n=\frac{1}{\left|detJ|_{\theta _n^0}\right|}=\frac{\theta _n^0}{\beta \frac{\beta }{\theta }|_{\theta _n^0}}=e^{A2n\pi }\frac{\left(3\sqrt{3}+2e^{A2n\pi }\right)D_{OS}}{2\beta D_{OL}^2D_{LS}}.$$
(25)
This expression gives a magnification decreasing very quickly and then the luminosity of the first image dominates all the others. The amplification diverges for $`\beta 0`$, confirming the fact that the possible detection of relativistic images is maximal for a perfect alignment of the source with the lens. In our approximation, the amplification of the images on the opposite side of the source is just the same of those on the same side.

Another interesting quantity is the total magnification of the relativistic images. As these would be seldom resolved as single images, it is likely to see them as one image with a total flux equal to the sum of the partial contributions coming from each image. Then, we just have to sum up the series

$$\mu _{tot}=2\underset{n=1}{\overset{\mathrm{}}{}}\mu _n,$$
(26)
which is a geometrical one, and then

$$\mu _{tot}=\frac{e^A\left(3\sqrt{3}\left(1+e^{2\pi }\right)+2e^A\right)D_{OS}}{\left(e^{4\pi }1\right)\beta D_{OL}^2D_{LS}}.$$
(27)
To take into account the finite extension of the source, one must integrate over its luminosity profile. As a simple example, integrating the $`\frac{1}{\beta }`$ dependence in Eqs. (25) and (27) over a uniform disk with angular radius $`\beta _S`$, we get

$`{\displaystyle \frac{1}{\pi \beta _S^2}}{\displaystyle \underset{D(\beta ,\beta _S)}{}}{\displaystyle \frac{1}{\beta ^{}}}d^2\beta ^{}`$ $`=`$ $`{\displaystyle \frac{2Sign\left[\beta _S\beta \right]}{\pi \beta _S^2}}[(\beta _S\beta )E({\displaystyle \frac{\pi }{2}},{\displaystyle \frac{4\beta _S\beta }{\left(\beta _S\beta \right)^2}})+`$ (28)
$`+`$ $`(\beta _S+\beta )F({\displaystyle \frac{\pi }{2}},{\displaystyle \frac{4\beta _S\beta }{\left(\beta _S\beta \right)^2}})],`$ (29)
where we have indicated the disk with radius $`\beta _S`$ centered on $`\beta `$ by $`D(\beta ,\beta _S)`$ and

$$E(\varphi _0,\lambda )=\underset{0}{\overset{\varphi _0}{}}\left(1\lambda \mathrm{sin}^2\varphi \right)^{1/2}𝑑\varphi ,$$
(30)
is the elliptic integral of the second kind.

Then, the magnification of an extended uniform source of angular radius $`\beta _S`$ can be obtained by substituting the $`1/\beta `$ dependence in Eq. (25) and Eq. (27) with the right member of Eq. (28).

## V Discussion

The formulae just derived provide a complete characterization of the two infinite sets of relativistic images surrounding a black hole or, in general, any compact object acting as a lens whose size is comparable with its Schwarzschild radius. They can be employed in all kinds of phenomenological calculations to test their detectability.

Some general considerations can be done on the features expected by a candidate lensing system. First, we need a very compact massive object, possibly a black hole, in order to gain access to the region of strong field immediately outside of the event horizon. Moreover, the matter surrounding this compact object should be transparent to the wavelength of the radiation emitted by the lensed source, otherwise the photons would be absorbed before the light rays complete their loops and would not reach the observer.

We recall that the amplification of each of the weak field images can be expressed (with distances in Schwarzschild units) as:

$$\mu _{wfi}=\frac{1}{\beta }\sqrt{\frac{2D_{LS}}{D_{OL}D_{OS}}}$$
(31)
for the small $`\beta _s`$ we are interested in. We can observe that the dependence on $`\beta `$ is the same of the weak field images. Then the relative importance of the relativistic images to the weak field ones is constant for high alignments of the source to the lens.

The ratio of this quantity and the total magnification of the relativistic images (27) is

$$\frac{\mu _{wfi}}{\mu _{tot}}=\frac{\sqrt{2}}{B}\left(\frac{D_{OL}D_{LS}}{D_{OS}}\right)^{3/2}$$
(32)
where $`B`$ is a numerical coefficient:

$$B=\frac{e^A\left(3\sqrt{3}\left(1+e^{2\pi }\right)+2e^A\right)}{\left(e^{4\pi }1\right)}=0.017$$
(33)
Then we must expect relativistic images to be always very faint with respect to the weak field images, since this ratio goes as an astronomical distance elevated to $`3/2`$. This fact rules out microlensing as a method for detection of relativistic images.

The separation between the two sets of relativistic images is of the order of the Schwarzschild diameter of the compact object. Depending on the specific situation, the angle corresponding to this length is generally very small. A very massive black hole would surely help to separate these images each other. The mass of the black hole at the center of our Galaxy  is believed to be about $`2.8\times 10^6M_{}`$. Its Schwarzschild angle would be about $`15`$ microsecs, which could become accessible by VLBI experiments . In this case, at least the two sets of relativistic images of some rear source would be distinguishable. With their very low magnification, depending on the alignment degree, the relativistic images of any rear source would cast a very hard but not impossible challenge for observational astronomy.

Possible better candidates could be black holes at the centers of other galaxies lensing some compact source (e.g. Quasars) on the background, or black holes at the centers of globular clusters.

## VI Conclusions

The expansion we have presented in this paper can be considered as the strong field gravitational lensing limit, in opposition to the usual weak field limit. Its starting point is the critical null geodesic followed by a massless particle captured by the black hole. We have expanded the deflection angle in this limit and used its leading order in the full lens equation. We have solved this equation and found two infinite sets of relativistic images. We have given a complete analytical description of these images with simple formulae for their position, magnification and for the critical curves.

The striking importance of these images lies on the fact that they could provide a profound test of general relativity in its full regime. In this situation, it would be surely possible to distinguish among relativistic theories of gravitation (e.g. Brans-Dicke, induced gravity, etc). For this reason, it is imperative to perform similar expansions within alternative pictures and evaluate the differences with the standard results we have found.

The scheme we have used for this calculation encourages its application to more complicated problems of general relativity such as charged rotating black holes or other relevant astrophysical objects where strong field general relativity is involved.

Acknowledgments

Work supported by fund 60% D.P.R. 382/80.


text:
# A Price Dynamics in Bandwidth Markets for Point-to-point Connections11footnote 1 ToN paper 01-1009. Submitted Jan. 17, 2001 to IEEE Trans. on Networking, revised Feb. 8, 2001

## 1 Introduction

To be able to provide guaranteed quality of service, QoS, in a network, a user needs to be able to reserve capacity, or ’bandwidth’, in congested routers. The reservation scheme should be efficient in the sense that one reservation should not unnecessarily block other reservations, and it should not require extensive negotiation. One scheme that fits these requirements is to trade router capacity in spot markets. The assumption is that someone reserves capacity for a connection by buying the capacity in the routers along the cheapest path between the source and destination node. When the capacity is no longer needed, it is sold to someone else. Increased demand increases price, so alternative paths, if they exist, may become competitive, and the users will tend to move their bandwidth usage away from congested routers.

We ultimately want to be able to price contingent claims on resources, formalized as options or futures. The working hypothesis is that adding a suitable set of such claims will improve the efficiency of the resource allocation, as prices will better reflect all available information. Indeed, the purpose of contingent claims on future value, such as options, futures, or more generally, derivative securities, is to construct a market for trading the expected future value. For the problem studied here, this means trading in expected future demand and supply. One way to price derivatives is to use statistical modeling of the price dynamics, since under suitable assumptions, which we will cover in a forthcoming separate contribution, derivative prices are functions of current market prices and the statistical model\[Black73\]\[Avellaneda99\]. This paper addresses the necessary preliminary issue of how to estimate parameters in two stochastic models from observed market prices, and also presents measures of the efficiency of the market resource allocation scheme. These questions are also of independent interest, as the efficiency is to be used to compare different market system with each other, and with other allocation schemes.

Using artificial markets for resource allocation in distributed systems dates back to the mid 80ies, ranging over markets for storage capacity \[Kurose89\], CPU time \[Ferguson88\] \[Waldspurger92\], and network capacity \[Kurose85\] \[Sairamesh95\]. The emphasis has been on evaluating the efficiency of the resource allocation, rather than understanding the resulting price dynamics. More recent work stresses the *agent* aspect, i.e. that the trading parties are locally optimizing entities \[Faratin00\]. Combinatorial markets, i.e. trading of bundles of distinct resources is yet a relatively new area. Somewhat related is the area of combinatorial auctions \[Rassenti82\] \[Rothkopf98\] \[Sandholm99\]. The use of derivatives for network admission control has been used in \[Lazar98\].

Previous bandwidth market models usually only include a primary market, in which end-users can buy and sell capacity only from the router owner. In the presented model end-users trade amongst each other, i.e. the capacity is traded on secondary markets.

This paper is organized as follows: in section 2 we present the model, of which a central ingredient is a detailed mechanism of price formation, with very low overhead in user-user communication. In section 3 we model mathematically the resulting price processes. The main tools used are stochastic differential equations (Itô-processes), their associated Fokker-Planck equations, and the stationary distribution of those. In section 4 we sum up an discuss our results.

## 2 Method

### 2.1 Market model

We simulate a network consisting of $`N`$ inter-connected routers and $`M`$ network users concurrently making reservations of router capacity for point-to-point connections. There is one spot market per router, in which users trade router capacity.

We use Farmer’s non-equilibrium market dynamics \[Farmer00\] as a prescription of how prices change due to trading (see also sec. 4). Farmer’s model is based on the assumption of a market maker that guarantees liquidity at all times, and that buying causes prices to increase, while selling causes prices to decrease. Farmer assumes that the price dynamics is such that it is impossible to move the market by performing a sequence of trades, where the net traded volume is zero, and that the relative price changes are independent of the current value of the price. From these assumptions, Farmer derives a formula that the price per unit in a transaction of $`\omega `$ units is $`\stackrel{~}{S}(S,\omega )=Se^{\omega /\lambda }`$, where $`S`$ is the unit price in the previous trade, and $`\lambda `$ the market depth or liquidity, i.e. the rate at which the price is changed by trading. For a derivation, see Appendix A. This model is useful since we do not have to simulate details of the order-book in each market. Instead, we can calculate the price change caused by trading directly from the last price and the size of the trade.

### 2.2 Simulation setup

At each time interval, $`m`$ new demands are generated. A demand is a 7-tuple,

$$d<id,uid,src,dst,cap,dur,max>$$

The user identities are chosen independently, so one user may receive zero or one or several new demands. If the current time is $`t,`$ demand $`id`$ specifies that the user $`uid`$ demands $`cap`$ units in each node on a path from $`src`$ to $`dst`$, starting at time $`t`$ and ending at time $`t+dur`$ if it costs less than $`max`$ to obtain the resources. If there are several paths, a choice will be made, see below. During the simulation a user reserves capacity in a router by buying that capacity, and sells excess capacity that is no longer in need.

User $`i`$ owns $`r_{i,j}`$ units of $`j`$. Initially, none of the simulated users are assigned any resources or money, i.e. $`r_{i,j}=0,cash_i=0`$ for all $`i`$ and $`j`$, nor do they have any resource demand, $`\omega _{i,j}=0`$. When a user manages to satisfy a demand, its capital is increased by the amount $`max`$, and when resources are bought and sold, the capital is decreased and increased, respectively, by the cost of the resource. The simulation is run in $`L`$ time steps from time $`0`$ to $`T`$ with time increments $`\mathrm{\Delta }t`$. At each time step where the current time is $`t`$:

* Generate $`m`$ new demands. A demand is specified by a unique demand number $`id`$, a user $`uid`$ randomly drawn from the set of $`M`$ users, a source node $`src`$ and a destination node $`dst`$, both randomly drawn from the set of $`N`$ nodes, and the required capacity per node $`cap`$, which is randomly drawn as the ceiling of $`e^{K\xi }`$ where $`\xi `$ is a uniform random variable between 0 and 1. If some demand is more important than others, a user is willing to pay more for that resource. In this simulation, $`max`$ is a linear function of the required capacity, i.e. the maximum total cost $`max=C_{unit}cap`$, where $`C_{unit}`$ is a simulation parameter. The duration (number of time steps) $`dur`$ is randomly drawn integer between $`1`$ and $`D`$.
* Calculate $`\omega _{i,j}`$, the net change of resource $`j`$ of user $`i=uid`$ in the following way:

+ For each new demand $`d`$, user $`i`$ looks at the last known transaction prices $`S_j(t)`$ and decides to buy $`cap`$ resources along *the least cost path.* Since the user will not know the actual cost of a buying and later selling resources along a path, its decision on whether to buy resources or not, is determined by a parameter $`C_{max}`$. The user decides to buy the resources if and only if the estimated total cost to buy the resources is less than $`C_{max}max`$. The resulting demand $`cap`$ is then added to the $`\omega _{i,j}`$ for all the resources $`j`$ on the least cost path, and the amount $`max`$ is added to $`cash_i`$.
    + For each satisfied demand $`d`$ terminating at $`t`$, decide to sell the resources that were allocated (if any) to the demand, i.e. subtract the resulting supply of liberated capacity on the least cost path if the demand was satisfied and router capacity was bought.
* All the demands on capacity in single routers (i.e. $`\omega _{i,j}`$) are effectuated one by one in random order, and for each trade prices are updated to $`S_je^{\omega _{i,j}/\lambda _j}`$. The $`cash_i`$ is decreased by $`\omega _{i,j}S_j`$ and the owned amount of resource is increased by $`\omega _{i,j}`$. According to the price formation formula, the trades are made at prices which depend on the actual order. The prices payed by the users are thus not the same as those used for determining the least cost path. The next known price is the one that pertains after all the trades, and is independent of the order.
* When all trading is done for this time-step, log the last transaction prices, $`\widehat{S}(t)`$, the number of satisfied demands, and repeat.

### 2.3 Simulation Parameters

The simulation was run for $`L=1000`$ time-steps ($`\mathrm{\Delta }t=0.01)`$ using the network in (fig. 1). There were $`10`$ routers and 10 users, so $`N=10`$ and $`M=10`$. The liquidity in all Farmer markets was chosen to be $`\lambda _i=10`$. The maximal cost per route users were willing to accept was determined by $`C_{unit}=100`$ and $`C_{max}=1`$ for all users. Initial prices were set to $`\widehat{S}_i(0)=10`$ in all markets. Every time step, $`m=10`$ new demands were generated. The duration $`dur`$ was uniformly distributed between $`1`$ and $`10`$, i.e. $`D=10`$. The call duration sets the time scale in the simulation, as we will see below. The required amount of capacity was determined by $`K=2`$ (see above).

## 3 Results

### 3.1 Statistical modeling of the price process

#### Inhomogeneous drift, additive noise, constant coefficients

The price as shown in (fig. 2) does not appear to be drifting freely. Instead it appears to return towards the same area. Assuming the price process is an Ornstein-Uhlenbeck process, the dynamics for the price of router $`i`$ would be

$$dS_i(t)=\alpha _i\left(\mu _iS_i(t)\right)dt+\sigma _idW_i(t)$$
(1)
where $`W_i(t)`$ is a Wiener process and the correlation between two processes $`i`$ and $`j`$ is $`Corr[dW_j(t),dW_j(t)]=\rho _{i,j}`$. Recall that a Wiener process has independent normal distributed increments with mean $`0`$ and variance $`ts`$, i.e. $`W_i(t)W(s_i)N[0,\sqrt{ts}]`$ where $`t>s`$. For improved readability, we omit the indices in $`S_i(t)`$, etc., when they are irrelevant for the understanding.

In (eq. 1) the drift term (the $`dt`$ term) detracts when $`S`$ is bigger than $`\mu `$, and adds when $`S`$ is less than $`\mu `$. The amount of the increase is determined by $`\alpha `$. The diffusion (the $`dW`$ term) is independent of $`S`$.

When the simulation has run for sufficiently long time, the price $`S_i(t)`$ becomes independent of the starting state of the system, and reaches a stationary distribution $`P(t)`$ of $`S(t)`$. Using the Fokker-Planck equation (see Appendix A) gives

$$P(s)=C_0e^{\frac{1}{2}\left(\frac{s\mu }{\sigma /\sqrt{2\alpha }}\right)^2}$$
(2)
where $`C_0=(\pi \sigma ^2/\alpha )^{\frac{1}{2}}`$, which is the density function for a normal distribution $`N[\mu ,\frac{\sigma }{\sqrt{2\alpha }}]`$. We note first that the normal distribution is non-zero in all of $`(\mathrm{},\mathrm{})`$, meaning that $`S(t)`$ can take on negative values, something that is not possible in a market with Farmer’s dynamics. Second, the normal distribution has many good properties, such as that a weighted sum of normal distributed variables is also normal distributed. If prices are far from zero, this dynamics may therefore be a convenient approximation to the true distribution. Third, the stationary distribution only depends on the ratio $`\frac{\sigma ^2}{\alpha }`$ and can therefore not distinguish between separate variations in $`\sigma ^2`$ and $`\alpha `$, which has to be done by other means.

The observations $`\widehat{S}(i)`$, $`i[1,\mathrm{},L]`$ of the process $`S`$ in its stationary state are regularly spaced with distance $`\mathrm{\Delta }t`$. Since $`E[S]=\mu ,`$ we estimate $`\mu `$ with

$$\widehat{\mu }=\widehat{E}[S]=\frac{1}{L}\underset{i=1}{\overset{L}{}}\widehat{S}(i)$$

We estimate $`\sigma `$ after noting that $`E[(dS)^2]=\sigma ^2dt+O(dtdW)`$. Therefore, for small $`\mathrm{\Delta }t`$,

$`\widehat{\sigma }^2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }t}}\widehat{E}[(dS)^2]`$
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }t}}{\displaystyle \frac{1}{L1}}{\displaystyle \underset{i=1}{\overset{L1}{}}}\left(\widehat{S}(i+1)\widehat{S}(i)\right)^2`$
Since $`Var[S]=\frac{\sigma ^2}{2\alpha }`$, the unbiased estimate is

$`{\displaystyle \frac{\widehat{\sigma }^2}{2\widehat{\alpha }}}`$ $`=`$ $`{\displaystyle \frac{L}{L1}}\left(\widehat{E}[S^2]\widehat{E}[S]^2\right)`$
so the estimate of $`\alpha `$ is

$$\widehat{\alpha }=\frac{1}{2\mathrm{\Delta }t}\frac{_{i=1}^{L1}\left(\widehat{S}(i+1)\widehat{S}(i)\right)^2}{\left(_{i=1}^L\widehat{S}(i)^2\left(_{i=1}^L\widehat{S}(i)\right)^2\right)}$$

To verify that the observed estimate is indeed correct, we can plot the histogram of the observed data together with the estimated density function $`P(s)`$, see (fig 3). This estimation coincides with the maximum likelihood estimate of the variables. A least square fit to the values in the histogram gives a better-looking curve. However, this fit only gives $`\mu `$ and the ratio $`\frac{\sigma ^2}{\alpha }`$, as noted above.

Having estimates of $`\alpha `$, $`\mu `$ and $`\sigma `$ we are able to estimate the correlation $`\rho _{i,j}`$ between the random sources for two price processes $`S_i`$ and $`S_j`$ by solving $`dW_i`$ and $`dW_j`$ in (eq. 1) and using the estimates

$$d\widehat{W}_i(k)=\frac{\widehat{S}_i(k+1)\widehat{S}_i(k)\widehat{\alpha }_i\left(\widehat{\mu }_i\widehat{S}_i(k)\right)\mathrm{\Delta }t}{\widehat{\sigma }_i}$$

so that the estimate becomes

$`\widehat{\rho }_{i,j}`$ $`=`$ $`{\displaystyle \frac{\widehat{\text{C}}\text{ov}[dW_i,dW_j]}{^{\widehat{\text{V}}\text{ar}[dW_i]\widehat{\text{V}}\text{ar}[dW_j]}}}`$
$`=`$ $`{\displaystyle \frac{\frac{1}{L1}_{k=1}^{L1}d\widehat{W}_i(k)d\widehat{W}_j(k)}{(\mathrm{\Delta }t)^2}}`$
where we have used $`E[dW]=0`$ and $`Var[dW]=dt`$.

#### Inhomogeneous drift, multiplicative noise, constant coefficients

The dynamics in (sec. 3.1) can generate negative prices, something which is should not be possible in a well-functioning market. By asserting a mean-reverting dynamics with multiplicative noise we get a process that is strictly positive. Assume the dynamics for the price of router $`i`$ to be

$$dS_i(t)=\alpha _i\left(\mu _iS_i(t)\right)dt+\sigma _iS_i(t)dW_i(t)$$
(3)
where $`W_i(t)`$ is a Wiener process and the correlation between two processes $`i`$ and $`j`$ is $`Corr[dW_j(t),dW_j(t)]=\rho _{i,j}`$. As before, we we omit the indices in $`S_i(t)`$, etc., for readability.

The stationary distribution $`P(s)`$ of $`S`$ is

$$P(s)=\frac{(\gamma \mu )^\gamma \mu }{\mathrm{\Gamma }(\gamma )}e^{\frac{\gamma \mu }{s}}\left(\frac{1}{s}\right)^{\gamma +2}$$
(4)
where $`\gamma \frac{2\alpha }{\sigma ^2}`$, and $`\mathrm{\Gamma }(z)`$ is the gamma function. In (eq. 4) $`s`$ takes only positive values, which is consistent, since the dynamics in (eq. 3) does not move an $`s`$ from positive to negative.

The first moment of this distribution, $`E[S]`$ is $`\mu `$ (see Appendix A). As for the Ornstein-Uhlenbeck process above, we estimate $`\mu `$ with

$$\widehat{\mu }=\widehat{E}[S]=\frac{1}{L}\underset{i=1}{\overset{L}{}}\widehat{S}(i)$$

We cannot easily estimate $`\sigma `$ from $`E[(dS)^2]`$ since $`S`$ is in the drift term, but note that the process $`X(t)=\mathrm{log}(S(t))`$ in (eq. 8) has additive noise. Therefore, we get for small $`\mathrm{\Delta }t`$,

$`\widehat{\sigma }^2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }t}}\widehat{E}[(dX)^2]`$
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }t}}{\displaystyle \frac{1}{L1}}{\displaystyle \underset{i=1}{\overset{L1}{}}}\left(\mathrm{log}{\displaystyle \frac{\widehat{S}(i+1)}{\widehat{S}(i)}}\right)^2`$
To estimate $`\alpha `$, we try two approaches. First we use the conditional expectation $`E[S(t+\tau )|S(t)]`$ for $`S`$ in the stationary state of the system. Taking the partial derivative w.r.t $`\tau `$ gives a first order ODE with the solution

$$E[S(t+\tau )|S(t)]=e^{\alpha \tau }(S(t)\mu )+\mu $$

showing that the expected value approaches $`\mu `$ exponentially with $`\tau `$ and $`\alpha `$ determines the speed of the return. Rearrange to keep $`e^{\alpha \tau }`$ (which is independent of $`S(t)`$) on one side, take the logarithm and take the expected value of both sides. Let $`\tau =k\mathrm{\Delta }t`$. We now have an estimate for $`\alpha `$,

$$exp(\widehat{\alpha }k\mathrm{\Delta }t)=\frac{1}{Lk}\underset{i=1}{\overset{Lk}{}}\frac{\widehat{S}(i+k)\widehat{\mu }}{\widehat{S}(i)\widehat{\mu }}$$

As we can see from the left hand side of the equation above, plotting the right hand side as a function of $`k`$ should result in a straight line if $`\alpha `$ is constant. However, as we can see in (fig. 4), the line is straight up to some $`k`$ that depends on the simulation parameters and error. The estimation is very sensitive to errors in $`\widehat{\mu }`$, especially in the denominator if $`\widehat{S}`$ is near $`\widehat{\mu }`$, which results in a flattened out jagged line, since $`e^{\widehat{\alpha }k\mathrm{\Delta }t}`$ is less than the simulation error for larger $`k`$. It eventually flattens out which means that $`\alpha `$ appears to decrease for larger $`\tau `$.

Another way to estimate $`\alpha `$ is to estimate $`\mu `$ and $`\sigma `$ as above, and then fit the observed distribution to the distribution (eq. 4) using the least square method. See (fig. 5) for a plot of the model fit.

Having estimates of $`\alpha `$, $`\mu `$ and $`\sigma `$ we are able to estimate the correlation $`\rho _{i,j}`$ between the random sources for two price processes in the same way as above, by solving $`dW_i`$ and $`dW_j`$ in (eq. 3) and using the estimates

$$\mathrm{\Delta }\widehat{W}_i(k)=\frac{\widehat{S}_i(k+1)\widehat{S}_i(k)\widehat{\alpha }_i\left(\widehat{\mu }_i\widehat{S}_i(k)\right)\mathrm{\Delta }t}{\widehat{\sigma }_i\widehat{S}_i(k)}$$

so that the estimate becomes

$$\widehat{\rho }_{i,j}=\frac{\frac{1}{L1}_{k=1}^{L1}d\widehat{W}_i(k)d\widehat{W}_j(k)}{(\mathrm{\Delta }t)^2}$$

where again we have used $`E[dW]=0`$ and $`Var[dW]=dt`$.

Plotting the histograms of $`\mathrm{\Delta }\widehat{W}`$ shows us if the model is good. If that is the case, the $`\mathrm{\Delta }\widehat{W}`$ should behave like samples of a Wiener process, i.e. be normal distributed. As can be seen in (fig. 6), the model (eq. 3) fits well if the market liquidity is high. However, if market liquidity is low then the assumed price model does not provide a good fit.

#### Price Correlations

The resource prices in a network of resources depend on the prices of other resources, since they are traded in groups. Using the parameter estimates of the price dynamics above, we get the correlation matrix in (tab. 1) for a simulation with high liquidity markets ($`\lambda _i=100)`$.

We find that prices generally are positively correlated. Comparing with the network graph(fig. 1), one can see that prices of neighboring nodes often are strongly positively correlated with an average correlation coefficient of around 0.4, compared to 0.03 for nodes that are not connected. The most significant exception is nodes 3 and 7. Looking at the network graph it is clear that no *least cost path* can contain both 3 and 7, except for the path from 3 to 7.

### 3.2 Efficiency of the market based resource allocation

The efficiency of a market based resource allocation scheme depends on how well prices reflect the available information about resource demands. Two sources for bad performance is that price quotes are outdated, or that they do not reflect knowledge about the price behavior, e.g. periodic price fluctuations, or future prices. To measure efficiency, we can measure the resource utilization to compare which markets are able to capture most information.

#### Communication costs

The number of messages sent in the negotiation phase of this kind of system is negligible, since all users operate on old price quotes and place bids at market, i.e. accept the price whatever it may be. They do not update their price quotes for every bid. Therefore, less than $`M`$ messages per trade come in to a user (the potential quote update from each market). One message per trade (the bid) go out from each user to the markets that contain resources that the user have chosen. No messages need to be communicated between the end-users.

#### Successful connection ratio

In the Farmer market model, the less liquid the markets are, the less valid are the price quotes. To the left in (fig. 7) the ratio of successful connections is plotted as a function of the market liquidity $`\lambda _i`$ for simulations with the parameters described in (sec. 2.3). The graphs corresponds to different values of the decision parameter $`C_{max}`$. Low liquidity causes large price fluctuations, making the prices higher than the limit $`C_{max}max`$ (see sec. 2.2) which inhibits many connections. Note that a connection is considered successful even if the net cost (after releasing the resources) is higher than $`max`$.

#### Net profit

To the right in (fig. 7) we plot the average profit as a function of the liquidity for a number of values of $`C_{max}`$. Large values of $`C_{max}`$ causes the users to buy resources when they are expensive. If the liquidity is high, the users sell resources at approximately the same price, but with low liquidity, prices will move significantly downwards when the resources are sold, causing a net loss to the user.

The particular way the $`max`$ cost for a connection is determined of course very much determines which of the different kinds of traffic that is promoted in the network. Different schemes could for instance promote short or long paths, high or low capacity connections, etc.

#### Average load

The average load, or reserved capacity in a router is with the Farmer dynamics a direct function of the market price, $`\omega (t)=\lambda \mathrm{log}\frac{S(t)}{S(0)}`$. If two simulations that allocate the same connections differ in load, the higher load depends on inefficient routing that does not choose efficient path. If the simulations differ in allocated connections, the comparison is less straight forward, and depends on the choice of metric above. We intend to return to a longer discussion on efficiency, in particular for load balancing, in a future contribution.

## 4 Discussion

In our simulation, Farmer’s market dynamics can be interpreted either as prescriptions of how a rational market maker should modify the market prices, or as a model of the aggregate behavior of market price changes during some period of time.

With the former interpretation, the Farmer dynamics can be used to implement a market maker program that brokers trading between end-users. The dynamics is derived from the assumption that it is impossible to change prices by ’trade in circles’. This assumption is a necessary condition for any market maker strategy, since otherwise anyone can exploit the market and gain an unlimited amount of money from the market maker.

The latter interpretation can be used when we want to simulate a part of a market with many concurrent trades. Since trading causes prices to change means that it is impossible to have updated price information at the time they place their bid. It is only possible to have completely updated price info if the trading is synchronous, something that severely reduces the number of bids a market can handle per time unit if communication delays are taken into account.

The market liquidity (or market depth) parameter $`\lambda _i`$ determines the speed at which prices change. When the dynamics interpreted as prescriptions for a market maker, it is up to the market maker to adjust $`\lambda _i`$ in order to reduce the risk of running out of resources. If we model an existing bandwidth market, $`\lambda _i`$ is an observable parameter which must be determined from the distribution of price jumps in relation to traded volume. The more resources in relation to the order size, the larger $`\lambda _i`$. A bandwidth market that is trading the capacity of a high capacity router can, everything else being equal, be expected to have a higher $`\lambda _i`$ than a similar market for a congested or inefficient router.

In the mean reverting processes investigated here, $`\alpha `$ determines the speed with which the process returns towards its statistical average $`\mu `$. The speed of return depends on the characteristic time length of the system, determined by the simulation parameter $`D`$, the call length. With a large $`D`$, the prices will be correlated with previous prices over a longer time period, which results in a small $`\alpha `$. With a small $`D`$, the price effect of previous calls will soon be forgotten, resulting in a high $`\alpha `$.

The traditional statistical models used for telephony have been found to inadequately model data network traffic. Data communication has been found to show a very bursty or fractal behavior as one communication event often generates a burst of more communication to other parts of the network. Data communication is generally short lived and often with strong latency bounds, due to the increasing use of computer networks for interactive communication. For short lived communication with low transfered amount, traditional switched best effort networks will probably continue to provide a very efficient solution for a long time. However, with increasing demand of streaming real-time data such as high quality video-telephony, video-on-demand, etc., it is necessary to be able to reserve capacity. The alternatives are large buffers, which has bad latency performance, or migrating data (intelligent replication), which can reduce load for one-to-many communication. Neither is suited for point-to-point communication service guarantees.

In the ’bandwidth market’ presented in this paper, end users buy the required resources themselves. The demand for updated price quotes results generates additional network traffic. In an extended model, we could allow risk neutral middle men to sell options on the resources. Since updated quotes will then only be required by the (few) middle men doing the actual trading, the additional traffic would be greatly reduced. We intend to return to these topics in future work.

## Acknowledgement

This work was supported by the Swedish National Board for Industrial and Technical Development (NUTEK), project P11253, COORD, in the Promodis program, and by the Swedish Research Institute for Information Technology (SITI), project eMarkets, in the Internet3 program. We want to thank Sverker Janson and Bengt Ahlgren for helpful discussions and comments.

## Appendix A

### Stationary Distribution

$`S(t)`$ is an Itô process with the dynamics $`dS(t)=a(t,S(t))dt+b(t,S(t))dW(t)`$, and $`P(s,T;s_0,t)`$ is the contingent probability distribution of $`S`$ at time $`T`$ given $`S(t)=s_0`$. $`P`$ obeys the Fokker-Planck equation a.k.a the forward Kolmogorov equation,

$$_TP_s\left(a(T,s)P\right)+_{ss}\left(\frac{b(T,s)^2}{2}P\right)=0$$
(5)
An important class of stochastic processes are those that are stationary. For those, if $`Tt`$ is sufficiently large, the contingent distribution $`P`$ no longer depends on $`s_0`$, and can be substituted with the stationary distribution $`P(s)`$, which does not depend on $`s_0`$, $`t`$ or $`T`$. Two ways in which a stochastic process can fail to be stationary is if the coefficients, e.g. $`a`$ and $`b`$, are explicitly time-dependent, or if there is no distribution $`P(s)`$. The latter happens for instance in ordinary random diffusion, in which the probability gradually spreads out without reaching a limit.

Assuming that a stationary distribution has been reached, the time derivative drops out of (eq. 5), and the equation can be integrated once to be

$$_sP(s)\frac{2a(s)2b(s)_sb(s)}{b(s)^2}P(s)=0$$
(6)
which is an ODE. One further integration of (eq. 6) gives

$$P(s)=Ce^{_0^s\frac{2a(u)}{b(u)^2}𝑑u}\frac{1}{b(s)^2}$$
(7)
where $`C`$ is a normalization constant determined by $`P(s)𝑑s=1`$.

### The Additive Noise Process

#### Stationary density function, mean and variance

Assume the price dynamics

$$dS(t)=\alpha (\mu S(t))dt+\sigma dW(t)$$

Then the conditional stationary probability distribution of $`S`$ is, using (eq. 7),

$`P(s)`$ $`=`$ $`C_0e^{\frac{2\alpha }{\sigma ^2}_0^s(\mu u)𝑑u}`$
$`=`$ $`C_0e^{\frac{1}{2}\left(\frac{s\mu }{\sigma /\sqrt{2\alpha }}\right)^2}`$
where $`C_0=(\pi \sigma ^2/\alpha )^{\frac{1}{2}}`$. $`P(s)`$ can be identified as the density function for a normal distribution with mean $`\mu `$ and variance $`\sigma ^2/2\alpha `$.

### The Multiplicative Noise Process

#### Stationary density function

Assume the price dynamics

$$dS(t)=\alpha (\mu S(t))dt+\sigma SdW(t)$$

Using (eq. 7)

$`P(s)`$ $`=`$ $`C_1e^{\frac{2\alpha }{\sigma ^2}\left(_{s_0}^s\frac{\mu }{u^2}𝑑u\stackrel{=\mathrm{log}s\mathrm{log}s_0}{\stackrel{}{_{s_0}^s\frac{1}{u}𝑑u}}\right)}\left({\displaystyle \frac{1}{\sigma s}}\right)^2`$
$`=`$ $`{\displaystyle \frac{1}{C_2}}e^{\frac{2\alpha }{\sigma ^2}\left(\frac{\mu }{s}\right)}\left({\displaystyle \frac{1}{s}}\right)^{\frac{2\alpha }{\sigma ^2}+2}`$
An alternative derivation is to let $`X(t)=\mathrm{log}S(t)`$. Itô’s lemma gives that

$$dX=\left(\alpha \mu e^X(\alpha +\frac{\sigma ^2}{2})\right)dt+\sigma dW$$
(8)
Then the stationary probability distribution $`Q(x)`$ of $`X`$ is , using (eq. 7),

$`Q(x)`$ $`=`$ $`C_0e^{_0^x\left(\frac{2\left(\alpha \mu e^X(\alpha +\frac{\sigma ^2}{2})\right)}{\sigma ^2}\right)𝑑u}{\displaystyle \frac{1}{\sigma ^2}}`$
$`=`$ $`C_1e^{\gamma \mu e^x+(\gamma +1)x}`$
where $`C_1`$ is a normalization constant, and $`\gamma \frac{2\alpha }{\sigma ^2}`$ for readability. The stationary density function $`P(s)`$ for $`S`$ is found by recalling $`X=\mathrm{log}(S)`$ and

$`P(s)ds`$ $``$ $`Q(\mathrm{log}(s))d(\mathrm{log}(s))`$
$`=`$ $`{\displaystyle \frac{1}{C_2}}e^{\gamma \mu \frac{1}{s}}\left({\displaystyle \frac{1}{s}}\right)^{\gamma +1}{\displaystyle \frac{1}{s}}ds`$
The constant $`C_2`$ is determined by normalization and identifying the integral as a gamma function, $`\mathrm{\Gamma }(z+1)=_0^{\mathrm{}}t^ze^t𝑑t=z\mathrm{\Gamma }(z)`$. After the substitution $`t=\gamma \mu \frac{1}{s}`$ we get

$`C_2`$ $`=`$ $`{\displaystyle _{\mathrm{}}^0}e^t\left({\displaystyle \frac{1}{\gamma \mu }}t\right)^{\gamma +2}(\gamma \mu {\displaystyle \frac{1}{t^2}})𝑑t`$
$`=`$ $`\left({\displaystyle \frac{1}{\gamma \mu }}\right)^{\gamma +1}{\displaystyle _0^{\mathrm{}}}e^tt^\gamma 𝑑t`$
$`=`$ $`\left({\displaystyle \frac{1}{\gamma \mu }}\right)^\gamma {\displaystyle \frac{1}{\mu }}\mathrm{\Gamma }(\gamma )`$
#### Stationary mean

Use the same substitution, $`t=\gamma \mu \frac{1}{s}`$ to obtain

$`E[S]`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}sP(s)𝑑s`$
$`=`$ $`{\displaystyle _{\mathrm{}}^0}({\displaystyle \frac{1}{\gamma \mu }}t){\displaystyle \frac{1}{C_2}}e^t\left({\displaystyle \frac{1}{\gamma \mu }}t\right)^{\gamma +2}(\gamma \mu {\displaystyle \frac{1}{t^2}})𝑑t`$
$`=`$ $`{\displaystyle \frac{1}{C_2}}\left({\displaystyle \frac{1}{\gamma \mu }}\right)^\gamma \underset{=\mathrm{\Gamma }(\gamma )}{\underset{}{{\displaystyle _0^{\mathrm{}}}e^tt^{\gamma 1}𝑑t}}`$
$`=`$ $`\mu `$
#### Conditional expected value

Let $`f(\tau |S(t))`$ be the conditional expectation of $`S(t+\tau )`$ given $`S(t)`$.

$`f(\tau |S(t))`$ $``$ $`E[S(t+\tau )|S(t)]`$
$`=`$ $`S(t)+E[{\displaystyle _t^{t+\tau }}𝑑S(t^{})|S(t)]`$
$`=`$ $`S(t)+\alpha \mu \tau \alpha {\displaystyle _t^{t+\tau }}E[S(t^{})|S(t)]𝑑t^{}`$
Take the partial derivative with respect to $`\tau `$

$`_\tau f(\tau |S(t))`$ $`=`$ $`\alpha \mu \alpha E[S(\tau )|S(t)]`$
$`=`$ $`\alpha \mu \alpha f(\tau |(S(t))`$
Multiplying with $`e^{\alpha \tau }`$ and collecting the terms gives

$$_\tau (e^{\alpha \tau }f(\tau |(S(t)))=e^{\alpha \tau }\alpha \mu $$

Integrate both sides over $`\tau `$ from $`0`$ to $`\tau `$ and multiplication with $`e^{\alpha \tau }`$ gives

$$f(\tau |S(t))=e^{\alpha \tau }\left(S(t)\mu \right)+\mu $$

where we have used that $`f(0|S(t))=S(t)`$. Therefore

$$e^{\alpha \tau }=\frac{f(\tau |S(t))\mu }{S(t)\mu }$$

The left side does not depend on $`S(t)`$. Taking the expected value of both sides

$`e^{\alpha \tau }`$ $`=`$ $`E[e^{\alpha \tau }]`$
$`=`$ $`E\left[{\displaystyle \frac{E[S(t+\tau )|S(t)]\mu }{S(t)\mu }}\right]`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{E[S(t+\tau )|S(t)=x]\mu }{x\mu }}P(x)𝑑x`$
where $`P`$ is the stationary distribution of $`S`$, and we have substituted a time average with an ensemble average, which is true under suitable assumptions.

#### Estimation of parameter $`\alpha `$

We can estimate the integral from the observed data using the estimates

$`\widehat{P}(x)`$ $``$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{t}{}}\delta (x\widehat{S}(t))`$
$`\widehat{E}[S(t+\tau )|S(t)=x]`$ $``$ $`{\displaystyle \frac{_{t|\widehat{S}(t)=x}\widehat{S}(t+\tau )}{|\{x|\widehat{S}(t)=x\}|}}`$
where $`\delta (x)`$ is the Dirac delta function and $`|A|`$ is the cardinality of the set $`A`$. The estimation is

$`e^{\widehat{\alpha }\tau }`$ $``$ $`{\displaystyle _0^{\mathrm{}}}\widehat{P}(x){\displaystyle \frac{\widehat{E}[S(t+\tau )|S(t)=x]\widehat{\mu }}{x\widehat{\mu }}}𝑑x`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{t}{}}\delta (x\widehat{S}(t))\right)\widehat{E}[{\displaystyle \frac{S(t+\tau )\widehat{\mu }}{x\widehat{\mu }}}|S(t)=x]𝑑x`$
$`=`$ $`{\displaystyle \underset{x|\widehat{S}(t)=x}{}}{\displaystyle \frac{1}{N}}|\{x|\widehat{S}(t)=x\}|{\displaystyle \frac{_{t|\widehat{S}(t)=x}\left(\frac{S(t+\tau )\widehat{\mu }}{x\widehat{\mu }}\right)}{|\{x|\widehat{S}(t)=x\}|}}`$
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{t}{}}{\displaystyle \frac{\widehat{S}(t+\tau )\widehat{\mu }}{\widehat{S}(t)\widehat{\mu }}}`$
so we can estimate $`\alpha `$ with

$$\widehat{\alpha }(\tau )\frac{1}{\tau }\mathrm{log}\left(\frac{1}{N}\underset{t}{}\frac{\widehat{S}(t+\tau )\widehat{\mu }}{\widehat{S}(t)\widehat{\mu }}\right)$$

which should be a constant function.

### Farmer’s market dynamics

Let $`S`$ be the current price, $`\omega `$ the net demand and $`\stackrel{~}{S}`$ the price at which the demand can be met. We seek a functional relationship of the type $`\stackrel{~}{S}\stackrel{~}{S}(S,\omega )`$, where $`\stackrel{~}{S}`$ depends continuously on $`S`$ and $`\omega `$. Assume the price to be positive and bounded, $`\stackrel{~}{S}`$ to be an increasing function of $`\omega `$, and that prices are only changed through trading, $`\stackrel{~}{S}(S,0)=S`$. Assume furthermore that one cannot make money by trading in circles

$$\stackrel{~}{S}(\stackrel{~}{S}(\stackrel{~}{S}(S,\omega _1),\omega _2),(\omega _1+\omega _2))=S$$
(9)
and the relative price change is independent of the absolute price,

$$\frac{\stackrel{~}{S}(S,\omega )}{S}=\varphi (\omega )$$
(10)
Because of (9) we see that $`\stackrel{~}{S}(\stackrel{~}{S}(S,\omega ),\omega )=S`$, so the inverse is $`\stackrel{~}{S}^1(S,\omega )=\stackrel{~}{S}(S,\omega )`$, where $`\stackrel{~}{S}^1`$ is the inverse function of $`\stackrel{~}{S}`$. Applying $`\stackrel{~}{S}^1(,(\omega _1+\omega _2))`$ on (9) to see that

$$\stackrel{~}{S}(\stackrel{~}{S}(S,\omega _1),\omega _2)=\stackrel{~}{S}(S,\omega _1+\omega _2)$$
(11)
(10) together with (11) implies that $`\varphi (\omega _1)\varphi (\omega _2)=\varphi (\omega _1+\omega _2)`$ which implies $`\varphi (x)=e^{x/\lambda }`$ for some constant $`\lambda `$. Therefore

$$\stackrel{~}{S}(S,\omega )=Se^{\omega /\lambda }$$
(12)


text:
# The universality class of fluctuating pulled fronts

## Abstract

It has recently been proposed that fluctuating “pulled” fronts propagating into an unstable state should not be in the standard KPZ universality class for rough interface growth. We introduce an effective field equation for this class of problems, and show on the basis of it that noisy pulled fronts in d+1 bulk dimensions should be in the universality class of the (d+1)+1D KPZ equation rather than of the d+1D KPZ equation. Our scenario ties together a number of heretofore unexplained observations in the literature, and is supported by previous numerical results.

Consider spatio-temporal systems in which the important dynamics is governed by the propagation of fronts or interfacial zones separating two domains whose bulk dynamics is relatively trivial or uninteresting. In the presence of fluctuations, the theory of the stochastic behavior of such fronts or interfaces is well-developed . In particular it is known that many such fluctuating $`d`$-dimensional interfaces in d+1 bulk dimensions are described by the KPZ equation  for their height $`h`$,

$$\frac{h}{t}=\nu ^2h+\lambda (h)^2+\eta ,$$
(1)
with $`\eta `$ a random gaussian noise with correlations

$`\eta (𝐫_{},t)=0,`$ (2)
$`\eta (𝐫_{},t)\eta (𝐫_{}^{},t^{})=2ϵ\delta ^d(𝐫_{}𝐫_{}^{})\delta (tt^{}).`$ (3)
We will follow common practice to refer to this equation as the d+1D(imensional) KPZ equation, where the $`d`$ refers to the dimension of the interface and the +1 to the time dimension; $`𝐫_{}`$ denotes the coordinates perpendicular to the direction of propagation of the interface.

The fact that the scaling behavior of so many stochastic interfaces fall in the d+1D KPZ universality class, is due to the fact that (1) this equation contains all the terms in a gradient expansion which are relevant in a RG sense; and (2) that the long wavelength deterministic dynamics of many interfaces is local in space and time, i.e., of the form $`v_n=v_n(h,^2h,\mathrm{})`$, expressing that the normal velocity $`v_n`$ becomes essentially a function of the instantaneous slope (angle) and curvature of the interface only. Upon expanding in the gradients, adding noise, and retaining only relevant terms, one then arrives at (1).

The starting point of such an argument, the fact that one can integrate out the internal structure of the interface and on long length and time scales think of it as a mathematically sharp boundary with effective dynamics expressed by a boundary condition $`v_n`$=$`v_n(h,^2h,\mathrm{})`$ which is local in space and time, is appealing and usually correct. Intuitively, one associates it with the interfacial zone being sufficiently sharp on a spatial scale. Nevertheless, there have been scattered observations in the literature which indicate that there is more to it: (a) Some continuum reaction-diffusion equations have propagating planar interfaces of finite width which are stable, but which become weakly unstable for discrete particle model equivalents , contrary to what the above coarse-graining picture would suggest. (b) The empirical relation observed for the distribution of DLA fingers in a channel and the interface shape of a viscous finger could not be understood from the standard continuum model until the innocuously looking reaction term was regularized ; on hindsight, this was because the standard mean-field DLA equations do not give the appropriate “local” boundary conditions of the type $`v_n=\mu _np`$. (c) In a simple stochastic particle model with fluctuating fronts, non-KPZ scaling was observed  contrary to what one would naively have expected.

It turns out that these observations all have one common denominator , in that they are related to the existence of two classes of fronts, “pushed” and “pulled” fronts. Pushed fronts are the usual ones: their dynamics is determined by the behavior in the interfacial zone, a region of finite thickness, and their response to the bulk fields is local in space and time . Pulled fronts, on the other hand, propagate into a linearly unstable state. Although they do not differ noticeably from pushed fronts in their appearance, their dynamics is driven by the growth and spreading of perturbations about the unstable state in the semi-infinite region ahead of the front ; hence they are particularly sensitive to slight changes in the dynamics there . These important differences led two of us  to propose recently that fluctuating variants of $`d`$-dimensional pulled fronts in d+1 bulk dimensions would indeed not be in the d+1D KPZ universality class, even though pushed fronts do effectively give local boundary conditions on long length and time scales, and hence do give rise to d+1D KPZ scaling in the absence of coupling to a diffusion or laplace field in the bulk . Simulations of a simple stochastic lattice model were consistent with these arguments, and with the earlier observations of .

In this paper, we will argue that fluctuating pulled fronts are indeed in a different universality class from the usual pushed ones which show the standard KPZ behavior. Indeed, we will show that the semi-infinite region ahead of the front can not be integrated out, and effectively enhances the dimension by 1: we introduce a field equation for fluctuating pulled fronts and argue that $`d`$-dimensional fronts in d+1 bulk dimension are in the universality class of the (d+1)+1D KPZ rather than the d+1D KPZ equation. This surprising scenario, which also builds on the insight of  for the stochastic behavior of pulled fronts in one bulk dimension, is fully consistent with our earlier 2D simulations  and also with the heretofore unexplained results of  in higher dimensions. In addition, as we shall discuss, our scenario leads to a number of interesting new questions and challenges.

A stochastic equation for pulled fronts should obey two requirements: in the usual stochastic lattice models with fronts, no particles are spontaneously generated when there are none already. Secondly, the average front speed and the local fluctuations ahead of the front remain always finite in such lattice models . The field equation should be consistent with these basic facts. So when we consider a stochastic field equation for $`\varphi (x,𝐫_{},t`$ in d+1 dimensions $`(x,𝐫_{})`$ of the type

$$\frac{\varphi }{t}=D^2\varphi +f(\varphi )+g(\varphi )\eta $$
(4)
these requirements put constraints on the function $`f`$ and the noise term $`g(\varphi )\eta `$. The stochastic noise $`\eta (x,𝐫_{},t)`$ has delta correlations as in (3), and is interpreted in the Stratonovich sense, but our arguments will not rely on the distinction between Ito and Stratonovich calculus. Stochastic field equations of this type have e.g. been used already since long for studying the scaling behavior of homogeneous bulk phases like directed percolation ; investigations of noisy fronts in such equations are more recent — see e.g.  for an analysis of stochastic pulled fronts and a discussion of the applicability to various systems. Here we focus on the proper form for an effective stochastic field equation for pulled fronts. For $`f(\varphi )`$, which determines the dynamics of deterministic fronts in the absence of noise, we choose the standard form for pulled front propagation $`f(\varphi )=\varphi \varphi ^3`$, which gives saturation of the field $`\varphi `$ behind the front where $`\varphi 1`$. How should the noise term $`g(\varphi )\eta `$ look like ? The requirement that if there are no particles ($`\varphi `$=$`0`$) none are created spontaneously, implies that there should be no additive noise term, and hence that $`g(\varphi `$=$`0)=0`$. For $`\varphi `$ nonzero but small, it is natural to assume a power law behavior $`g(\varphi )\varphi ^\alpha `$; in the studies of the homogeneous bulk properties of directed percolation the choice $`\alpha =1/2`$ was been made , motivated by the idea that typical bulk fluctuations are of the order of the square root of the particle density. For pulled fronts, however, the dynamically important region is ahead of the front, where $`\varphi 0`$. Our second requirement that the relative fluctuations $`g(\varphi )\eta /\varphi `$ remain finite here shows that the natural choice is $`\alpha =1`$, i.e., $`g(\varphi )\varphi `$ for $`\varphi 1`$ . The linearity of $`g`$ for small $`\varphi `$ is sufficient for our subsequent analysis. In our numerical studies, we have actually taken $`g(\varphi )=\varphi (1\varphi ^2)`$, a form taken to suppress fluctuations behind the front. This makes it numerically easier to focus on the fluctuations of the front position itself, without affecting the essential results. Specifically, we thus propose as the generic stochastic field equation for pulled fronts

$$\frac{\varphi }{t}=D^2\varphi +(1+\eta )\varphi (1\varphi ^2).$$
(5)
Let us now turn to the analysis of the stochastic behavior of fronts which propagate along the $`x`$-direction into the linearly unstable state $`\varphi `$=$`0`$. The crucial feature of pulled fronts is that even though the full dynamics of the fronts is nonlinear, it is essentially determined in the “leading edge”, the region ahead of the front where $`\varphi `$ remains small enough that the nonlinear saturation term $`\varphi ^3`$ which limits the growth, plays no role: the linear spreading and growth of perturbations about the state $`\varphi `$=$`0`$ almost literally “pull the front along”. An important recent development has been the realization that this simple intuitive picture can be turned into a systematic scheme to calculate even the convergence of the front speed to its asymptotic value $`v^{}`$. Remarkably, this relaxation is governed by universal power laws which can be calculated exactly even for general equations . The fact that the stochastic fluctuation effects that we want to investigate are dominant relative to the deterministic velocity relaxation terms, suggests to calculate these along similar lines. For the deterministic case ($`\eta `$=$`0`$), fronts in (5) propagate with an asymptotic speed $`v^{}=2\sqrt{D}`$. In a frame $`\xi =xv^{}t`$ moving with this speed, the asymptotic front solution has an exponential fall-off $`e^{\lambda ^{}\xi }`$ with $`\lambda ^{}=1/\sqrt{D}`$ for large positive $`\xi `$. The asymptotic relaxation analysis of deterministic fronts is based on the so-called leading edge transformation $`\varphi =e^{\lambda ^{}\xi }\psi `$, which transforms (5) into

$$\frac{\psi }{t}=D^2\psi +\eta \psi (1+\eta )\psi ^3e^{2\lambda ^{}\xi }.$$
(6)
Here $`\psi `$ has been written in the frame moving with velocity $`v^{}`$ in the $`x`$-direction: $`\psi =\psi (\xi ,𝐫_{},t)`$.

In the analysis of deterministic fronts , the nonlinear term on the right hand side (which is exponentially small for $`\xi 1`$) essentially plays the role of a boundary condition for the semi-infinite leading edge region where $`\varphi `$ is small — it allows the nonlinear region to match properly to the leading edge which “pulls” the front. As explained above, this holds a fortiori for fluctuating pulled fronts: their stochastic fluctuations are essentially determined by the region where the linearized equation can be used. Now, as is well known, upon making a Cole-Hopf transformation $`\psi =e^h`$, the linearized equation transforms to

$$\frac{h}{t}=D^2h+D(h)^2+\eta ,$$
(7)
which is nothing but the (d+1)+1D KPZ equation (1) for the d+1 dimensional field $`h(\xi ,𝐫_{},t)`$!

As illustrated for two bulk dimensions in Fig. 1, the 1D fluctuations in the front position in the propagation direction are defined by tracing a line where $`\varphi =const`$, e.g., $`\varphi =1/2`$. Since $`\varphi =e^{\lambda ^{}\xi +h}`$, the front fluctuations in the $`\xi `$ direction are given by $`\xi (𝐫_{},t)=h(\xi ,𝐫_{},t)/\lambda ^{}+\xi _0`$ $``$ $`h(\xi _0,𝐫_{},t)/\lambda ^{}+\xi _0`$, where the constant $`\xi _0`$ is determined by the level curve of $`\varphi `$ which we trace to determine the front position. Thus indeed the position fluctuations of a $`d`$-dimensional pulled front in d+1 bulk dimensions map onto the height fluctuations along a line of a KPZ surface in d+1 dimensions — see Fig. 1 The growth and roughness exponents are therefore those of the (d+1)+1D KPZ universality class!

The above scenario unifies a number of different results. It can immediately be compared with the simulation results of the stochastic lattice model of . In that paper a 2D lattice model was introduced in which by changing a simple birth and death rule of particles 1D fronts could be tuned from pushed to pulled. The scaling exponents of the pushed model were found to be the standard 1+1D KPZ ones, as it should, while those of the pulled variants were close to those of the 2+1D KPZ universality class. More importantly, without any adjustable parameters, the distribution functions for the long-time saturated width of the fronts in this model for finite transverse width $`L_{}`$  are completely in accord with our scenario . Moreover, although fronts in 1D do not have transverse fluctuations, the wandering of the position of pulled fronts in one dimenension is also consistent with 1+1D KPZ scaling . Finally, the observations of Riordan et al.  that in three (and higher) bulk dimensions their fronts did not appear to show a power law growth of the front width finds a natural explanation. According to  their fronts are pulled and so they should be governed by the 3+1D KPZ equation. The free $`\lambda =0`$ fixed point in this equation is stable and has no divergent interface width. Apparently above two dimensions the model of  is in the weak-coupling limit.

On hindsight, our arguments also justify the regularization of  of the mean-field equations for DLA in a channel: the full problem involves pushed fronts but the mean-field equations have pulled front solutions. The regularization effectively cures this by making the fronts into pushed ones.

The validity of the crucial step of our derivation, the assertion that the nonlinearities in (5) or (6) can be neglected because the leading edge where $`\varphi `$$``$$`1`$ is the essential region, can be tested independently. In Fig. 2 we show simulation data of the wandering of the lines where $`\varphi =0.5`$ in Eq. (5) in 2D, both with and without the nonlinearity. Following , where the linearized version of (6) was already employed to study the 2+1D KPZ exponents numerically, we have taken parameters so as to make the dimensionless coupling $`\stackrel{~}{\lambda }=2\lambda ^2ϵ/\nu ^325`$. This value appears to be close to the fixed point value and so slow transients are minimized . We find indeed that the two datasets with and without the nonlinear term in (5) show the same growth exponent, with a value close to the one $`\beta 0.24`$ of the 2+1D KPZ equation. This gives confidence in the validity of our assertion that the nonlinear terms in the front equation are not important for the scaling behavior of pulled fronts.

The main steps of our line of argument are elegantly direct and build on various previously established ideas; at the same time our scenario also raises a number of new questions and challenges for further research:
(i) There is no systematic theory for the transition from the pushed to the pulled regime in stochastic lattice models, so it is difficult to determine a priori which models lead to the standard pushed fronts and which ones to pulled ones. E.g., fronts in the directed percolation problem are pushed and obey KPZ scaling in one special case , but it is not known whether this is generally so.
(ii) Finite size scaling of the KPZ equation is normally done for interfaces of size $`L_{}`$ in all directions. Our scenario, on the other hand, leads one to consider anisotropic scaling, since there is effectively a time-dependent cutoff in the $`\xi `$-direction . The crossover scaling is completely unexplored, but is most likely quite tricky: for fixed $`L_{}`$ the results of  for fronts in one dimension suggest that one should see subdiffusive wandering of the average front position, $`(\overline{h})^2\sqrt{t}`$ (rather than $`t`$) because the cutoff in the $`\xi `$-direction grows as $`L_\xi \sqrt{t}`$, but our simulations seem to suggest that the crossover to this regime happens at such extremely long times that it can not convincingly be seen in practice. Moreover, the crossover is likely to depend significantly on the initial conditions .
(iii) According to the results of , pulled fronts are very sensitive to finite particle effects, so that the convergence to a continuum limit is extremely slow. This crossover is poorly understood, especially in higher dimensions where it is not even known how important the effects are in practice.

In conclusion, we have put forward an effective field equation for pulled fronts and argued on the basis of it that pulled fronts in d+1 bulk dimensions are in the (d+1)+1D KPZ universality class rather than the d+1D KPZ universality class. The scenario ties together various results in the literature and brings up various new issues for future research.

WvS would like to thank Uwe Tauber and David Mukamel for stimulating discussions. Financial support from the Dutch science foundation FOM, the Spanish project BXX2000-0638-C02-02 and the TMR network ERBFMRX-CT96-0085 is gratefully acknowleged.