inordinatum

Physics and Mathematics of Disordered Systems

Archive for the ‘Technical Tricks’ Category

Tricks for inverting a Laplace Transform, part V: Pole Decomposition

with 4 comments

This is a continuation of my articles on methods for inverting Laplace transforms. You can find the previous parts here: part I (guesses based on series expansions), part II (products and convolutions), part III, part IV (substitutions).

1. Result

In this post I will explain how to find the inverse P(x) of the following Laplace transform:

\displaystyle \int_0^\infty \mathrm{d}x\, e^{-sx}P(x) = \frac{\tanh \sqrt{2s}}{\sqrt{2s}}.   (1)

The solution is given in terms of the Jacobi theta function \theta_2 (which can be re-expressed in terms of Jacobi elliptic functions, though I won’t discuss that here in detail):

\displaystyle P(x) = \frac{1}{2}\theta_2 \left(0; e^{-x\pi^2/2}\right).   (2)

By taking derivatives or integrals with respect to s, which are equivalent to multiplication or division of P by x, one can obtain many other Laplace transform identities, including for example

\displaystyle \int_0^\infty \mathrm{d}x\, e^{-sx}\frac{1}{2x}\theta_2 \left(0; e^{-x\pi^2/2}\right) = \ln \cosh \sqrt{2s}.

If anyone manages to find a systematic list of those, I’d be very grateful. But for now let’s just see how one obtains (2).

2. Derivation

First, we use the classic decomposition of trigonometric functions in infinite products due to Euler:

\displaystyle   \begin{array}{rl}  \sinh z &= z \prod_{n=1}^{\infty}\left(1+\frac{z^2}{\pi^2 n^2}\right) \\  \cosh z &= \prod_{n=1}^{\infty}\left(1+\frac{z^2}{\pi^2 \left(n-\frac{1}{2}\right)^2}\right)  \end{array}  .

From the second identity, we can obtain a partial fraction decomposition of \tanh z (following this post on StackExchange):

\displaystyle   \begin{array}{rl}  \tanh z = & \partial_z \ln \cosh z = \partial_z \sum_{n=1}^{\infty} \ln \left(1+\frac{z^2}{\pi^2 \left(n-\frac{1}{2}\right)^2}\right)   \\  = & \sum_{n=1}^{\infty} \frac{\frac{2z}{\pi^2 \left(n-\frac{1}{2}\right)^2}}{1+\frac{z^2}{\pi^2 \left(n-\frac{1}{2}\right)^2}} \\  = & 2z \sum_{n=1}^{\infty} \frac{1}{\frac{\pi^2(2n-1)^2}{4}+z^2}.  \end{array}  .

Applying this to the right-hand side of (1), we obtain a sum over simple poles:

\displaystyle \int_0^\infty \mathrm{d}x\, e^{-sx}P(x) = \frac{\tanh \sqrt{2s}}{\sqrt{2s}} = 2\sum_{n=1}^{\infty} \frac{1}{\frac{\pi^2(2n-1)^2}{4}+2s}

The Laplace inverse of a simple pole is just an exponential, \int_0^\infty \mathrm{d}x\, e^{-sx}\,e^{-p x}=\frac{1}{p+s}. By linearity of the Laplace transform, we can invert each summand individually, and obtain an infinite sum representation for P(x):

\displaystyle P(x) = \sum_{n=1}^{\infty} \exp\left(-\frac{\pi^2(2n-1)^2}{8}x\right)

This sum can now be evaluated with Mathematica‘s Sum command, or by hand using the representation of theta functions in terms of the nome, for argument z=0 \Leftrightarrow w=1 and q=e^{-\pi^2 x/2}. This finally gives the solution as claimed above:

\displaystyle P(x) = \frac{1}{2}\theta_2 \left(0; e^{-x\pi^2/2}\right).   (2)

Have fun!

Advertisements

Written by inordinatum

April 15, 2013 at 10:27 pm

Mean-field solution of the Random Field Ising Model

with 7 comments

A happy new year to all the blog readers out there!

Today I will discuss a classical model of a disordered system that is extremely simple to write down: The Ising model with quenched (i.e. time-invariant) local random fields.
This model, also dubbed the Random-Field Ising Model (RFIM), is highly non-trivial. Above two dimensions, depending on the strength of the local random fields, it exhibits a phase transition between an ordered (ferromagnetic) and a disordered (paramagnetic) phase. When applying an additional homogeneous time-varying magnetic field, the local random fields prevent some regions from flipping. This leads to interesting phenomena like avalanches and hysteresis. The avalanche sizes, and the shape of the hysteresis loop, vary depending on the location in the phase diagram. Experimentally, this phenomenology – including the disorder-induced phase transition, avalanches, and hysteresis – describes very well e.g. helium condensation in disordered silica aerogels.

In this blog post, I will discuss the solution of the RFIM in the fully connected limit. This means that the ferromagnetic interaction between any two spins, no matter how far apart, has the same strength. This limit, also called the mean-field limit, is supposed to describe the model accurately for a sufficiently high spatial dimension. Just how high is “sufficient” is still a matter of debate – it may be above four dimensions, above six dimensions, or only when the dimension tends to infinity. Nevertheless, the mean-field limit already captures some of the interesting RFIM phenomenology, including the ordered-disordered phase transition. This is what I will discuss in this blog post. I will follow closely the original paper by T. Schneider and E. Pytte from 1977, but I’ve tried to extend it in many places, in order to make the calculations easier to follow.

1. The mean-field RFIM Hamiltonian and its free energy

The energy of N fully connected Ising spins S_i=\pm 1, i=1...N, subject to local random fields h_i, is:

\displaystyle H[S]=-\frac{J}{N}\sum_{i,j}S_i S_j -\sum_i S_i h_i

The prefactor 1/N in the first term (the ferromagnetic interaction energy) is chosen so that the scaling of this term with N is the same as the scaling of the second term (containing the random fields).

The corresponding partition sum Z is the thermal average over all spin configurations:

\displaystyle Z=\sum_{\{S_i^\alpha=\pm 1\}} \,e^{-\beta H[S]}

As usual, \beta = T^{-1} is the inverse temperature of the system.
In order to see the phase transition, we need to consider the free energy F, averaged over all random fields h:

\displaystyle \left\langle F\right\rangle_h = -T\left\langle \ln\, Z \right\rangle_h

Through a somewhat lengthy calculation (which I explain in detail in the last section 3. of this post), one obtains the following expression for the disorder-averaged free energy:

\displaystyle \left\langle F\right\rangle_h = N\left[J m^2 - T \int \mathrm{d} h\, p(h) \ln 2 \cosh \beta \left(2 J m + h\right) \right]

where the typical magnetization m is fixed through the implicit (self-consistent) equation

\displaystyle m = m_{SC}(m):=\int \mathrm{d} h\, p(h) \tanh \beta \left(2 J m + h\right)

p(h) is the probability distribution of the random fields, assumed to be Gaussian with mean zero and variance \sigma^2:

\displaystyle p(h) = \frac{1}{\sqrt{2\pi \sigma^2}}\exp\left(-\frac{h^2}{2\sigma^2}\right)

Although the expression above for the free energy is written in terms of p(h), strictly speaking I will show it only for the case where p(h) has this Gaussian form. Whether the formula for the free energy also holds for other p(h) is a good question which I’d like to understand better – if any expert happens to be reading this and can clarify, it would be great!

Now, let’s look at what information on the phase diagram of the model we can deduce from these expressions.

2. Phase diagram of the mean-field RFIM

Self-consistent equation for the magnetization of the RFIM

Self-consistent equation for the magnetization of the RFIM. The dashed line indicates the consistency condition m_{SC}(m)=m. The red line is an example of m_{SC}(m) in the ferromagnetic phase (here J=1, T=0.5, \sigma=1) It has a non-trivial value m>0 where it satisfies the self-consistency condition, i.e. intersects the dashed line. On the other hand, the yellow line is an example of m_{SC}(m) in the paramagnetic phase (here J=1, T=0.5, \sigma=2). There, the only point satisfying the self-consistency condition is the trivial solution m=0.


The self-consistent equation for the magnetization m given above has a trivial solution 0=m_{SC}(0) for all values of the parameters \sigma, J, T. Similar to what happens in the mean-field solution for the standard Ising model, if \beta is sufficiently large, the slope m_{SC}'(0) is larger than one and there is another non-trivial solution m>0 (see figure on the right). This indicates the existence of ferromagnetic order. Thus, the phase boundary between the paramagnetic and the ferromagnetic phase is given precisely by the condition

\displaystyle   \begin{array}{rl}  m_{SC}'(m=0) =& 1 \\  \Leftrightarrow 2\beta J \int \mathrm{d} h\, \frac{p(h)}{\left(\cosh \beta h\right)^2} = &1   \end{array}

One of the three parameters \beta, J, \sigma can be eliminated by rescaling h. For example, we can eliminate \beta and express the phase transition line in terms of the dimensionless variables \tilde{J}:=\beta J, \tilde{\sigma}:= \beta \sigma, \tilde{h}:=\beta h:

\displaystyle 2 \tilde{J}(\tilde{\sigma}) = \left[\int{\frac{\mathrm{d}\tilde{h}}{\sqrt{2\pi\tilde{\sigma}^2}}}\frac{e^{-\frac{1}{2}\tilde{h}^2/\tilde{\sigma}^2}}{\left(\cosh \tilde{h}\right)^2}\right]^{-1}

Another way is to scale away \sigma, setting \hat{J}:= J/\sigma, \hat{\beta}:=\sigma \beta.

Then one has the following expression for the phase boundary

\displaystyle 2\hat{J}(\hat{\beta}) = \left[\int{\frac{\mathrm{d}\tilde{h}}{\sqrt{2\pi}}}\frac{e^{-\frac{1}{2}\tilde{h}^2/\hat{\beta}^2}}{\left(\cosh \tilde{h}\right)^2}\right]^{-1}

With this expression, we see more clearly that there is a well-defined zero-temperature limit. When \hat{\beta}\to\infty, we get a transition at

\displaystyle 2\hat{J}(\infty) = \left[\int{\frac{\mathrm{d}\tilde{h}}{\sqrt{2\pi}}}\frac{1}{\left(\cosh \tilde{h}\right)^2}\right]^{-1} = \sqrt{\pi/2}

This is quite interesting, since it indicates a phase transition driven purely by the fluctuations of the quenched disorder, without any thermal noise!

Phase diagram of the mean-field RFIM

Phase diagram of the mean-field RFIM

More generally, we can easily plot the phase diagram in terms of \hat{J}, \hat{\beta} using Mathematica; see the figure on the left.

Of course, the expression for the free energy given above can also be used to deduce more elaborate observables, like the susceptibility, etc. For this, I refer to the original paper by T. Schneider and E. Pytte. For example, one can find that near the phase transition the magnetization with a power-law,

\displaystyle m(T) \sim (T-T_c)^{1/2}

This is the typical mean-field power law, which is also found when studying the phase transition in the pure Ising model on the mean-field level.

Now that we’re convinced that these formulas for the free energy and the magnetization (boxed in orange above) are useful for understanding the physics of the mean-field RFIM, I will explain in detail how they were derived.

3. Calculating the free energy using the replica trick

Since averages over the logarithm of a random variable are not easy to compute, we use the replica trick to re-express the free energy in terms of moments of Z:

\displaystyle \left\langle F\right\rangle_h = -T \frac{\partial}{\partial n}\Big|_{n=0}\left\langle Z^n \right\rangle_h

Z^n is the partition sum of n copies, or replicas, of the system, feeling the same random fields h. Indexing these copies by greek indices \alpha=1...n, we get

\displaystyle \left\langle Z^n \right\rangle_h = \sum_{\{S_i^\alpha=\pm 1\}} e^{\frac{\beta J}{N} \sum_\alpha \sum_{i, j}S_i^\alpha S_j^\alpha}\left\langle e^{\beta \sum_i \left(\sum_\alpha S_i^\alpha\right) h_i}\right\rangle_h

Since the h_i are independent and Gaussian (with mean 0 and variance \sigma^2), \left\langle e^{\lambda h_i}\right\rangle_h = \int \mathrm{d}h\,p(h)\,e^{\lambda h} = e^{\sigma^2\lambda^2/2}. Then

\displaystyle \left\langle Z^n \right\rangle_h = \sum_{\{S_i^\alpha=\pm 1\}} \exp\left[ \frac{\beta J}{N} \sum_\alpha \sum_{i, j} S_i^\alpha S_j^\alpha + \frac{\beta^2\sigma^2}{2}\sum_i\left(\sum_\alpha S_i^\alpha\right)^2\right]

We have performed the disorder averaging, so what remains is the partition sum of n copies/replicas of a non-disordered fully-connected Ising model. They are interacting via the second term in the exponential, which came from the disorder averaging. We will now compute the partition sum of this replicated Hamiltonian. The plan for this is the following:

  • Decouple the individual sites i using a Hubbard-Stratonovich transformation.
  • Apply the saddle-point method to compute the integral. Due to the decoupling, the saddle-point is controlled by the number of sites N. Thus, this procedure is exact in the thermodynamic limit N\to\infty.

To apply the Hubbard-Stratonovich transformation, we first rewrite the ferromagnetic interaction as

\displaystyle \sum_\alpha \sum_{i, j} S_i^\alpha S_j^\alpha = \sum_\alpha \left(\sum_i S_i^\alpha \right)^2

and then use

\displaystyle \exp\left[\frac{\beta J}{N} \sum_\alpha \left(\sum_i S_i^\alpha \right)^2\right] = \frac{1}{(2\pi)^{n/2}}\int\prod_\alpha\mathrm{d}x_\alpha\,\exp\left[\sum_\alpha \left(x_\alpha \sqrt{2\frac{\beta J}{N}}\sum_i S_i^\alpha - \frac{1}{2}x_\alpha^2\right)\right]

In order to have an exponential scaling homogeneously with the system size, we set x = \sqrt{N} \tilde{x}:

\displaystyle \left\langle Z^n \right\rangle_h = \left(\frac{N}{2\pi}\right)^{n/2}\sum_{\{S_i^\alpha=\pm 1\}} \int\prod_\alpha\mathrm{d}\tilde{x}_\alpha\,\exp\left[-\frac{N}{2}\sum_\alpha\tilde{x}_\alpha^2 + \sqrt{2\beta J}\sum_i \sum_\alpha \tilde{x}_\alpha S_i^\alpha  + \frac{\beta^2\sigma^2}{2}\sum_i\left(\sum_\alpha S_i^\alpha\right)^2\right]

The exponential is now a linear sum over all sites i, without any interaction between them. Thus, the partition sum is the N-th power of the partition sum of a single site:

\displaystyle \left\langle Z^n \right\rangle_h = \left(\frac{N}{2\pi}\right)^{n/2} \int\prod_\alpha\mathrm{d}\tilde{x}_\alpha\,\exp\, N\left[-\frac{1}{2}\sum_\alpha\tilde{x}_\alpha^2 + \ln Z_1(\tilde{x})\right],

where Z_1(\tilde{x}) is the partition sum for all replicas on a single site:

\displaystyle Z_1(\tilde{x}) = \sum_{\{S^\alpha=\pm 1\}} \exp\left[\sqrt{2\beta J}\sum_\alpha \tilde{x}_\alpha S^\alpha  + \frac{\beta^2\sigma^2}{2}\left(\sum_\alpha S^\alpha\right)^2\right]

In the thermodynamic limit N\to\infty, the x integral in our expression for \left\langle Z^n \right\rangle_h is dominated by its saddle point (due to the fact that the exponential is linear in N). At the saddle point, by symmetry we can assume all the x_\alpha to be equal. Setting \tilde{x}_\alpha := x for all \alpha, the exponent of the integrand in the above expression for \displaystyle \left\langle Z^n \right\rangle_h becomes

\displaystyle N\left[-\frac{n}{2}x^2 + \ln Z_1(x)\right], \quad \text{where} \quad    Z_1(x) = \sum_{\{S^\alpha=\pm 1\}} \exp\left[\sqrt{2\beta J}x \sum_\alpha S^\alpha  + \frac{\beta^2\sigma^2}{2}\left(\sum_\alpha S^\alpha\right)^2\right].

Hence, the saddle-point value x=x_m is fixed by the implicit equation

\displaystyle \partial_x \big|_{x=x_m}\left[-\frac{n}{2}x^2 + \ln Z_1(x)\right]=0 \Leftrightarrow n\, x_m = \frac{\partial_{x}|_{x=x_m}Z_1(x)}{Z_1(x_m)} = \sqrt{2\beta J}\frac{\sum_{\{S^\alpha=\pm 1\}} \left(\sum_\alpha S^\alpha\right)e^{A[S,x_m]}}{\sum_{\{S^\alpha=\pm 1\}} e^{A[S,x_m]}},

where

\displaystyle A[S,x] = \sqrt{2\beta J}\sum_\alpha \tilde{x}_\alpha S^\alpha  + \frac{\beta^2\sigma^2}{2}\left(\sum_\alpha S^\alpha\right)^2.

We see that effectively, the saddle-point location x_m corresponds to a fixed average magnetization m := x_m / \sqrt{2\beta J}, for the replicated spins S^\alpha on a single lattice site, subject to the “Hamiltonian” A[S,x]. Rewriting everything in terms of the magnetization m, we get

\displaystyle \left\langle Z^n \right\rangle_h \approx \exp\, N\left[-n\beta J m^2 + \ln Z_1(m)\right],

\displaystyle   \begin{array}{rl}  Z_1(m) = & \sum_{\{S^\alpha= \pm 1\}}e^{A[S,m]} \\  A[S,m] = & 2\beta J m\sum_\alpha S^\alpha  + \frac{\beta^2\sigma^2}{2}\left(\sum_\alpha S^\alpha\right)^2 \end{array}

m is fixed by the implicit (self-consistent) equation

\displaystyle n \, m = \frac{1}{Z_1(m)}\sum_{\{S^\alpha= \pm 1\}}\left(\sum_\alpha S^\alpha\right)e^{A[S,m]}

Now, to further simplify the result, we can compute the sums over all spin configurations \{S^\alpha\} by performing another Hubbard-Stratonovich transformation:

e^{A[S,m]} = \int \frac{\mathrm{d} s}{\sqrt{2\pi}}\exp\left[-\frac{1}{2}s^2 + (2\beta J m + \beta \sigma s)\sum_\alpha S^\alpha\right]

Then, the single-site partition sum can be computed as following:

\displaystyle   \begin{array}{rl}  Z_1(m) = &\sum_{\{S^\alpha= \pm 1\}}e^{A[S,m]} \\  = &\int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, e^{-\frac{1}{2}s^2} \prod_{\alpha}\sum_{S^\alpha= \pm 1}e^{(2\beta J m + \beta \sigma s)S^\alpha} \\  = & \int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, e^{-\frac{1}{2}s^2}\left[2\cosh(2\beta J m + \beta \sigma s)\right]^n   \\  = &\int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, \exp\left[-\frac{1}{2}s^2 + n \ln 2\cosh(2\beta J m + \beta \sigma s)\right]   \end{array}

Similarly, the self-consistent equation for the mean magnetization simplifies to:

\displaystyle   \begin{array}{rl}  n \, m = &\frac{1}{Z_1(m)}\sum_{\{S^\alpha= \pm 1\}}\left(\sum_\alpha S^\alpha\right)e^{A[S,m]} \\  =& \frac{1}{Z_1(m)}\int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, e^{-\frac{1}{2}s^2} \frac{1}{2\beta J }\frac{\partial}{\partial m}\prod_{\alpha}\sum_{S^\alpha= \pm 1}e^{(2\beta Jm + \beta \sigma s)S^\alpha}    \\  =& \frac{1}{Z_1(m)}\int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, e^{-\frac{1}{2}s^2} \frac{1}{2\beta J}\frac{\partial}{\partial m} \left[2\cosh(2\beta Jm + \beta \sigma s)\right]^n   \\  =& \frac{1}{Z_1(m)}\int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, \exp\left[-\frac{1}{2}s^2 + n \ln 2\cosh(2\beta Jm + \beta \sigma s)\right]\,n\, \tanh (2\beta Jm + \beta \sigma s) \\  \Rightarrow m = & \frac{1}{Z_1(m)}\int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, \exp\left[-\frac{1}{2}s^2 + n \ln 2\cosh(2\beta Jm + \beta \sigma s)\right]\tanh (2\beta Jm + \beta \sigma s)  \end{array}

In order to obtain the disorder-averaged free energy, as discussed above we need to take the derivative with respect to n at n=0:

\displaystyle   \begin{array}{rl}  \left\langle F\right\rangle_h = &-T \frac{\partial}{\partial n}\Big|_{n=0}\left\langle Z^n \right\rangle_h \\  = & N\left[Jm^2 -\frac{\partial}{\partial n}\Big|_{n=0} Z_1(m) \right] \\  = & N\left[Jm^2 - \int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, e^{-\frac{1}{2}s^2} \ln 2\cosh(2\beta Jm + \beta \sigma s)\right]   \end{array}

The magnetization m is fixed by the self-consistent equation at n=0:

\displaystyle   m = \int \frac{\mathrm{d} s}{\sqrt{2\pi}}\, e^{-\frac{1}{2}s^2}\tanh (2\beta Jm + \beta \sigma s)

By the substitution h := \sigma s, these formulae reduce to those shown in section 1, and then used to analyze the phase diagram in section 2.

As always, I’m happy about any questions or comments!!

Written by inordinatum

January 20, 2013 at 7:16 pm

How does an absorbing boundary modify the propagator of a stochastic process?

with one comment

A recurring theme throughout my blog are stochastic processes, widely used to model quantities under the influence of both deterministic evolution laws and random fluctuations. A complete solution for a stochastic process is its propagator: The probability density for the final state at some time t_f given an initial condition at time t_i. The propagator depends on the boundary conditions imposed on the stochastic process. For example, an absorbing boundary is useful for describing extinction in a stochastic model of population dynamics: When the population of a given species reaches zero, it cannot be re-introduced.
Typically, however, the “free” propagator for a stochastic process — i.e. the propagator without any boundaries — is much easier to determine than the propagator with any given boundary condition.

However, common sense tells us that the free propagator for a stochastic process contains already all there is to know about it. So I have wondered for quite a while if and how one can recover e.g. the propagator with an absorbing boundary, knowing only the free propagator. I found the solution only quite recently in this preprint on arXiv, and think that it is interesting enough to explain it here briefly.

The result

Consider a one-dimensional, stationary, Markovian stochastic process. Let us denote its free propagator by

\displaystyle G(x_i,x_f|t_f-t_i).

This is the density for the final value x_f at time t_f, given an initial value x_i at time t_i. Stationarity of the process implies that G_0 is a function of the time difference \tau := t_f-t_i only. Define the Laplace transform of the free propagator by

\displaystyle \hat{G}(x_i,x_f|\lambda) := \int_{0}^{\infty} \mathrm{d}\tau\, e^{\lambda \tau}G(x_i,x_f|\tau).

Now let us add an absorbing boundary at x=a. Then, the propagator

\displaystyle G_a(x_i,x_f|t_f-t_i)

is the density for the final value x_f at time t_f, given an initial value x_i at time t_i, without touching the line x=a between t_i and t_f. The upshot of this post is the following simple relationship between the Laplace transforms of the free propagator G and of the propagator with absorbing boundary G_a :

\displaystyle \hat{G}_a(x_i,x_f|\lambda) = \hat{G}(x_i,x_f|\lambda)-\frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\hat{G}(a,a|\lambda)}.      (1)

In case the denominator on the right-hand side is zero, one may have to perform some limiting procedure in order to make sense of this formula. We can now check how nicely it works for simple cases, like Brownian motion with drift, where the result is known by other means.

Example: Brownian motion with drift

Let us consider an example, the Brownian motion with drift \mu. Its stochastic differential equation is

\displaystyle  \partial_t X(t) = \mu + \xi(t),

where \xi is white noise with \overline{\xi(t)\xi(t')} =2\sigma\delta(t-t').

The free propagator is a simple Gaussian

\displaystyle  G(x_f,x_i|\tau) = \frac{1}{\sqrt{4\pi\sigma\tau}}\exp\left[-\frac{(x_f-x_i-\mu \tau)^2}{4\sigma \tau}\right].

Its Laplace transform is given by:

\displaystyle  \begin{array}{rl}  \hat{G}(x_f,x_i|\lambda) =& \int_{0}^{\infty} \mathrm{d}\tau\, e^{\lambda \tau}G(x_i,x_f|\tau) \\  =&\frac{1}{\sqrt{\mu^2-4\lambda\sigma}}\exp\left[\frac{\mu}{2\sigma}(x_f-x_i) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-x_i|\right].  \end{array}

To compute G_a, let us assume x_i > a, x_f > a. Note this is without loss of generality, since if x_i and x_f are on opposite sides of the absorbing boundary the propagator G_a is zero. Applying the result (1), we have

\displaystyle   \begin{array}{rl}  \hat{G}_a(x_i,x_f|\lambda) =& \hat{G}(x_i,x_f|\lambda)-\frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\hat{G}(a,a|\lambda)}  \\  \quad=& \frac{1}{\sqrt{\mu^2-4\lambda\sigma}}\exp\left[\frac{\mu}{2\sigma}(x_f-x_i) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-x_i|\right] - \\   & \frac{1}{\sqrt{\mu^2-4\lambda\sigma}}\exp\left[\frac{\mu}{2\sigma}(a-x_i) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|a-x_i| \right. \\  &\quad \left.+\frac{\mu}{2\sigma}(x_f-a) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-a|\right]  \\  \quad=& \frac{1}{\sqrt{\mu^2-4\lambda\sigma}}e^{\frac{\mu}{2\sigma}(x_f-x_i)}\left[e^{-\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-x_i|} - e^{-\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}(x_i+x_f-2a)}\right].    \end{array}

The inverse Laplace transform of this expression is:

\displaystyle   \begin{array}{rl}  G_a(x_i,x_f|\tau) =& G(x_i,x_f|\tau)-e^{\frac{\mu}{\sigma}(a-x_i)}G(2a-x_i,x_f|\tau) \\  =& \frac{1}{\sqrt{4\pi\sigma\tau}}e^{-\frac{\mu^2 \tau}{4\sigma} + \mu\frac{x_f-x_i}{2\sigma}}\left[e^{-\frac{(x_f-x_i)^2}{4\sigma\tau}}-e^{-\frac{(2a-x_f+x_i)^2}{4\sigma\tau}}\right].  \end{array}

This is, of course, exactly the same result (the so-called Bachelier-Levy formula) which I had already discussed earlier in a guest blog post at Ricardo’s blog, using slightly different methods (a variant of the Cameron-Martin-Girsanov theorem). Note that some signs which were wrong in the original post are corrected here.

The proof

Now let us go back to the case of a general stochastic process, and discuss how formula (1) can be proven. Let us assume that the stochastic process in question has a path-integral representation, i.e. that

\displaystyle G(x_i,x_f|\tau) = \int_{x(t_i)=x_i}^{x(t_f)=x_f} \mathcal{D}[x(t)]\,e^{-S[x(t)]}

for some action S.

Now, to add an absorbing boundary at x=a we weigh each crossing of x=a by a factor e^{-\gamma} in the path integral, and then let \gamma \to \infty. This effectively retains only those trajectories in the path integral, which never touch x=a. We thus have:

\displaystyle G_a(x_i,x_f|\tau) = \lim_{\gamma \to \infty}\int_{x(t_i)=x_i}^{x(t_f)=x_f} \mathcal{D}[x(t)]\,e^{-S[x(t)]-\gamma\delta(x(t)-a)}.

Now, to evaluate the path integral, we expand it in a series in \gamma (or, in other words, in the number k of crossings of x=a). We get:

\displaystyle   \begin{array}{rl}  \lefteqn{G_a^{(\gamma)}(x_i,x_f|\tau) = G(x_i,x_f|\tau) + }&\\  & + \sum_{k=1}^{\infty} (-\gamma)^k \int_{t_i}^{t_f}\mathrm{d}t_1\, G(x_i,x_1|t_1-t_i) \int_{t_1}^{t_f}\mathrm{d}t_2\, G(x_1,x_2|t_2-t_1) \cdots \\  &  \cdots \int_{t_{k-1}}^{t_f}\mathrm{d}t_k\,G(x_{k-1},x_k|t_k-t_{k-1})\,G(x_k,x_f|t_f-t_k)  \end{array}  .

Note that here I cancelled the factor \frac{1}{k!}, coming from the Taylor expansion of e^{-\gamma \delta(x(t)-a)}, against the k! possibilities of ordering the times t_1\cdots t_k. This is why I can assume them to be in ascending order above. Also, the \delta functions actually set all x_1=x_2=\cdots=x_k=a.
Taking Laplace transforms, the time integrals (which are a k-fold convolution) simplify to a simple product:

\displaystyle  \begin{array}{rl}  \hat{G}_a^{(\gamma)}(x_i,x_f|\lambda) =& \hat{G}(x_i,x_f|\lambda) + \\  & + \sum_{k=1}^{\infty} (-\gamma)^k \hat{G}(x_i,a|\lambda)\hat{G}(a,a|\lambda)^{k-1} \hat{G}(a,x_f|\lambda) \\  =& \hat{G}(x_i,x_f|\lambda) - \frac{\gamma \hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{1+\gamma \hat{G}(a,a|\lambda)}\\  =& \hat{G}(x_i,x_f|\lambda) - \frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\gamma^{-1}+ \hat{G}(a,a|\lambda)}  \end{array}.

Now, one easily obtains the claimed result:

\displaystyle  \begin{array}{rl}  \hat{G}_a(x_i,x_f|\lambda) =& \lim_{\gamma \to \infty}\hat{G}_a^{(\gamma)}(x_i,x_f|\lambda) \\  =& \hat{G}(x_i,x_f|\lambda) - \frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\hat{G}(a,a|\lambda)}  \end{array}  .

This computation can, of course, be generalized to more complicated situations like several boundaries, higher dimensions, etc. For details see again the original preprint by C. Grosche. Have fun and let me know if you have questions!

Written by inordinatum

November 24, 2012 at 11:50 pm

Tricks for inverting a Laplace Transform, part IV: Substitutions

with 9 comments

After a few less technical posts recently, I now continue the series of articles on tricks for the inversion of Laplace transforms. You can find the other parts here: part I (guesses based on series expansions), part II (products and convolutions), part III, and part V (pole decomposition).

Today’s trick will be variable substitutions. Let’s say we need to find P(x) so that

\displaystyle \int_0^\infty \mathrm{d}x\, e^{\lambda x}P(x) = \hat{P}(r(\lambda))=e^{a r(\lambda)}\left[1-r(\lambda)\right]^b,   (1)

where r(\lambda) = \frac{1}{2}\left(c-\sqrt{c^2-\lambda}\right), and a,b,c are constants. This example is not as absurd as it may seem; it actually came up recently in my research on a variant of the ABBM model.

There is a general formula permitting one to invert the Laplace transform above, in terms of an integral in the complex plane:

\displaystyle P(x) = \int_\gamma \frac{\mathrm{d} \lambda}{2\pi i} e^{-\lambda x} \hat{P}(r(\lambda)).

This so-called Fourier-Mellin integral runs along a contour \gamma from -i \infty to +i\infty. \gamma can be chosen arbitrarily, as long as the integral converges and all singularities (poles, branch cuts, etc.) of \hat{P} lie on the right of it. Note that my convention for the Laplace variable \lambda is the convention used for generating functions in probability theory, which has just the opposite sign of \lambda of the “analysis” convention for Laplace transforms.

The fact that our Laplace transform \hat{P} is written entirely in terms of the function r(\lambda) suggests applying a variable substitution. Instead of performing the Fourier-Mellin integral over \lambda, we will integrate over r(\lambda). We then have:

\displaystyle P(x) = \int \frac{\mathrm{d}r}{2\pi i} \frac{\mathrm{d}\lambda(r)}{\mathrm{d}r} e^{-\lambda(r) x + a r}\left(1-r\right)^b.  (2)

Solving the definition of r(\lambda) for \lambda, we get

\displaystyle \lambda(r) = 4 (c r -r^2)

and

\displaystyle \frac{\mathrm{d}\lambda}{\mathrm{d}r} = 4 c - 8r.

Inserting this into (2), we have

\displaystyle P(x) = 4\int_{\gamma'} \frac{\mathrm{d}r}{2\pi i} (c-2r)e^{4 r^2 x - (4cx-a)r}\left(1-r\right)^b.

The only singularity of the integrand is at r=1, for b \in \mathbb{R} \backslash \mathbb{N}. We can thus choose the contour \gamma' to go parallel to the imaginary axis, from 1-i \infty to 1+i\infty. Mathematica then knows to how to evaluate the resulting integral, giving a complicated expression in terms of Kummer’s confluent hypergeometric function \,_1 F_1(a,b,z) = M(a,b,z). However, a much simpler expression is obtained if one introduces an auxiliary integral instead:

\displaystyle g(q,c,d):= \int_{-\infty}^{\infty} \mathrm{d}h\, e^{-q h^2}\left(c-i h\right)^d.

Mathematica knows a simple expression for it:

\displaystyle g(q,c,d) = \sqrt{\pi} q^{-(d+1)/2}U(-d/2;1/2;c^2 q).

U is Tricomi’s confluent hypergeometric function, which is equivalent to Kummer’s confluent hypergeometric function but gives more compact expressions in our case.
Using this auxiliary integral, (2) can be expressed as

\displaystyle P(x) = \left[\frac{4}{\pi}g\left(4x,1+\frac{a-4cx}{8x},b+1\right)+\frac{c-4}{\pi}g\left(4x,1+\frac{a-4cx}{8x},b\right)\right]\exp\left[-\frac{(a-4cx)^2}{16x}\right].

Simplifying the resulting expressions, one obtains our final result for P(x):

\displaystyle   \begin{array}{rcl} P(x) & = & \frac{4}{\sqrt{\pi}}(4x)^{-1-b/2}\left[U\left(-\frac{1+b}{2};\frac{1}{2};\frac{(a+4x(2-c))^2}{16x}\right) + \right. \\  & & \left.+(c-2)\sqrt{x}U\left(-\frac{1+b}{2};\frac{1}{2};\frac{(a+4x(2-c))^2}{16x}\right)\right]\exp\left[-\frac{(a-4cx)^2}{16x}\right].  \end{array} (3)

Equivalently, the confluent hypergeometric functions can be replaced by Hermite polynomials:

\displaystyle    P(x) =  \frac{8}{\sqrt{\pi}}(8x)^{-1-b/2}\left[H_{1+b}\left(\frac{a+4x(2-c)}{4\sqrt{x}}\right) + 2(c-2)\sqrt{x}H_b\left(\frac{a+4x(2-c)}{4\sqrt{x}}\right)\right]\exp\left[-\frac{(a-4cx)^2}{16x}\right].

For complicated Laplace transforms such as these, I find it advisable to check the final result numerically. In the figure below you see a log-linear plot of the original expression (1) for \hat{P}(\lambda), and a numerical evaluation of \int_0^\infty \mathrm{d}x\, e^{\lambda x}P(x), with P(x) given by (3). You can see that they match perfectly!

Numerical verification of the Laplace transform. Yellow line: Original expression \hat{P}(\lambda) as given in (1). Red crosses: Numerical evaluation of \int_0^\infty \mathrm{d}x\,e^{\lambda x}P(x) for the final expression P(x) given in (3).

Written by inordinatum

November 4, 2012 at 2:14 pm

A quick introduction to the Martin-Siggia-Rose formalism

with 20 comments

The Martin-Siggia-Rose (MSR) formalism is a method to write stochastic differential equations (SDEs) as a field theory formulated using path integrals. These are familiar from many branches of physics, and are simple to treat perturbatively (or, in some cases, even to solve exactly).

We start from an SDE of the form:

\displaystyle \partial_t x(t) = F(x(t),t) + \xi(x(t),t),

where F is a deterministic function and \xi(x(t),t) is Gaussian noise with an arbitrary correlation function G:

\displaystyle \left\langle \xi(x,t) \xi(x',t') \right\rangle = G(x,t,x',t').

MSR recognized that observables, averaged over solutions of this SDE, can be written as a path integral with a certain action:

\left\langle \mathcal{O}[x(t)] \right\rangle = \int \mathcal{D}[x,\tilde{x}]\,\mathcal{O}[x(t)]\,e^{-S[x,\tilde{x}]},
where
\displaystyle S[x,\tilde{x}] = \int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t),t)\right] + \frac{1}{2}\int_{t,t'} G(x(t),t,x(t'),t')\tilde{x}(t)\tilde{x}(t').

A useful special case is an automonous (time-independent) SDE with multiplicative noise:

\displaystyle \partial_t x(t) = F(x) + H(x)\xi(t),

where \overline{\xi(t)\xi(t')} = 2\sigma \delta(t-t') is Gaussian white noise. The resulting action for the path integral is local in time:

\displaystyle S[x,\tilde{x}] = \int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t))\right] + \frac{\sigma}{2} H(x(t))^2\tilde{x}(t)^2.

I will now show a few concrete examples, and then explain how this mapping is derived.

Example: Ornstein-Uhlenbeck process

The Ornstein-Uhlenbeck process is one of the few stochastic processes treatable analytically. Its SDE is:
\displaystyle \partial_t x(t) = -\alpha x(t) + \xi(t),
where \alpha > 0 (else the process blows up) and \xi(t) is white noise:
\left\langle \xi(t)\xi(t') \right\rangle = 2\sigma \delta(t-t').

In this section, I will show that the generating functional for this Ornstein-Uhlenbeck process is given by:
\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \exp \left[ \frac{\sigma}{2\alpha}\int_{-\infty}^\infty \mathrm{d} s_1 \int_{-\infty}^\infty \mathrm{d} s_2 \,e^{-\alpha|s_1-s_2|}\, \lambda(s_1) \lambda(s_2) \right] .

To show this, we apply the MSR formula to express the generating functional as a path integral:

\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \int \mathcal{D}[x,\tilde{x}] \exp\left\{-\int_t i\tilde{x}(t)\left[\partial_t x(t) +\alpha x(t)\right] - \sigma\int_{t} \tilde{x}(t)^2 + \int_t \lambda(t) x(t)\right\}.

This path integral can be solved exactly. The key observation is that the exponent is linear in x(t), and hence one can integrate over this field. Using \int_x e^{i x q} \propto \delta(q), this gives (up to a normalization factor \mathcal{N}) a \delta-functional:

\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \mathcal{N}\int \mathcal{D}[\tilde{x}] e^{ \sigma\int_{t} \tilde{x}(t)^2 } \delta\left[i\left(\partial_t - \alpha \right) \tilde{x}(t)+\lambda(t)\right].
The e^{ \sigma\int_{t} \tilde{x}(t)^2 } factor remains since it is the only term in the action independent of x(t).

The \delta-functional fixes \tilde{x}(t) in terms of \lambda:
\displaystyle \tilde{x}(t) = i \int_{t}^\infty \mathrm{d} s \,e^{\alpha(t-s)}\, \lambda(s) .

From this we get
\displaystyle \int_t \tilde{x}(t)^2 = - \int_t \int_{t}^\infty \mathrm{d} s_1 \,e^{\alpha(t-s_1)}\, \lambda(s_1) \int_{t}^\infty \mathrm{d} s_2 \,e^{\alpha(t-s_2)}\, \lambda(s_2) \\  =  - \frac{1}{2\alpha}\int_{-\infty}^\infty \mathrm{d} s_1 \int_{-\infty}^\infty \mathrm{d} s_2 \,e^{-\alpha|s_1-s_2|}\, \lambda(s_1) \lambda(s_2)

The normalization factor \mathcal{N} is seen to be 1 by imposing that for \lambda=0 we must have \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = 1.
Plugging this in gives the generating functional, as claimed above:
\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \exp\left[ \frac{\sigma}{2\alpha}\int_{-\infty}^\infty \mathrm{d} s_1 \int_{-\infty}^\infty \mathrm{d} s_2 \,e^{-\alpha|s_1-s_2|}\, \lambda(s_1) \lambda(s_2) \right] .

One simple consequence is the exponential decay of the connected correlation function:
\displaystyle \left\langle x(t)x(t')\right\rangle^c = \frac{\delta^2}{\delta \lambda(t) \delta \lambda(t')}\big|_{\lambda=0} \ln \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \frac{\sigma}{2\alpha}e^{-\alpha|s_1-s_2|}.

But of course, the full generating functional tells us much more – in some sense, it is the exact solution of the SDE!

Note added 01/12/2012: Of course, the generating functional can be written down immediately once one knows the two-point correlation function and the fact that the Ornstein-Uhlenbeck process is Gaussian. However, the latter fact is not immediately obvious (and was not used in our derivation), in fact one can see our MSR calculation as one way to prove Gaussianity for the OU process.

Derivation of the Martin-Siggia-Rose formula

Now that we’ve seen what the MSR path integral formulation is useful for, let us see how it can be derived. The average over an observable can be written as a path integral with the SDE enforced through a \delta-functional:
\displaystyle \left\langle \mathcal{O}[x(t)] \right\rangle = \left\langle\int \mathcal{D}[x]\,\mathcal{O}[x(t)]\,\delta\left[\partial_t x(t) - F(x(t),t) - \xi(x(t),t)\right]\right\rangle.
Now, we rewrite the \delta-functional using the formula \delta(q) \propto \int_{\tilde{x}} e^{i \tilde{x} q}. Since the \delta-functional is a product of a \delta-function at all points in space, effectively we are introducing an auxiliary field \tilde{x}(t):
\displaystyle \left\langle \mathcal{O}[x(t)] \right\rangle = \left\langle\int \mathcal{D}[x,\tilde{x}]\,\mathcal{O}[x(t)]\,\exp\left\{-\int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t),t) - \xi(x(t),t)\right]\right\}\right\rangle.
Now, the only thing that depends on the noise \xi(x(t),t) is the factor e^{\int_t i\tilde{x}(t)\xi(x(t),t)}. Since \xi is Gaussian, the average of this factor can be simply expressed as
\displaystyle \left\langle e^{\int_t i\tilde{x}(t)\xi(x(t),t)} \right\rangle = \exp\left[-\frac{1}{2}\int_{t,t'}\tilde{x}(t)\tilde{x}(t')\left\langle\xi(x(t),t)\xi(x(t'),t')\right\rangle\right] =  \exp\left[-\int_{t,t'}\tilde{x}(t)\tilde{x}(t')G(x(t),t,x(t'),t') \right]
Altogether we obtain
\displaystyle \left\langle \mathcal{O}[x(t)] \right\rangle = \int \mathcal{D}[x,\tilde{x}]\,\mathcal{O}[x(t)]\,\exp\left\{-\int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t),t)\right]-\int_{t,t'}\tilde{x}(t)\tilde{x}(t')G(x(t),t,x(t'),t') \right\}
which is just the MSR formula claimed above.

Amazingly, I did not yet find a clear derivation of this path integral representation on Wikipedia or anywhere else on the web… So I hope this will be useful for someone starting to work with the formalism!

If there is interest, I may write another blog post in the future with more examples and explanations e.g. on the physical role of the auxiliary field \tilde{x}.

Written by inordinatum

September 27, 2012 at 11:03 pm

Tricks for inverting a Laplace Transform, part III

with 6 comments

This is a continuation of the series of articles on Laplace transforms… You can also have a look at part I (guesses based on series expansions), part II (products and convolutions), part IV (substitutions), and part V (pole decomposition).

This time we will cover the following inverse Laplace problems:

\displaystyle \int_0^\infty e^{\lambda x}P(x)\,dx = e^{\frac{1}{a+b \lambda}}

The basic result that can be used for such problems is the Laplace transform of the Bessel function:

\int_0^\infty e^{\lambda x}\left[\frac{1}{\sqrt{x}}I_1(2\sqrt{x})+\delta(x)\right]dx = e^{\frac{1}{\lambda}}

Now, using

\displaystyle LT\{f(x)e^{cx}\}(s)=\int_0^\infty \mathrm{d}x\,e^{s x}f(x)e^{c x} = LT\{f(x)\}(s+c),

and by rescaling \lambda, one finds that \int_0^\infty e^{\lambda x}P(x)\,dx = e^{\frac{1}{a+b \lambda}} is solved by:

\displaystyle P(x)= \left[\frac{1}{\sqrt{bx}}I_1\left(2\sqrt{\frac{x}{b}}\right)+\delta(x)\right]e^{-\frac{a}{b}x}

Have fun!

Written by inordinatum

October 15, 2011 at 11:17 am

An interesting integral equation

leave a comment »

Recently, I discussed an interesting problem with a friend of mine. Basically, the objective was to find an approximation for the solution of a linear N\times N equation system in the limit of N\to\infty.

It turned out that in this limit the solution can be obtained by solving the following integral equation

\displaystyle h(x)=\int_{-1}^1 \frac{f(x)-f(y)}{x-y}\,\mathrm{d}y,

where h(x) is known and f(x) is to be found. Now, such integral equations are notoriously difficult to solve (in particular, much more difficult than differential equations). An excellent resource (which contains most known exact solutions) is the EqWorld page on integral equations (and the “Handbook of Integral Equations” by EqWorld’s authors, which you can find by googling). Unfortunately, it does not help in this case (the closest thing one finds there are the eigenfunctions for the integral operator \int_{-1}^1 \frac{f(x)-f(y)}{|x-y|}\,\mathrm{d}y, which, however, turn to to be competely unrelated to those of \int_{-1}^1 \frac{f(x)-f(y)}{x-y}\,\mathrm{d}y).

One angle of attack to our equation is the observation that plugging in a polynomial f(x) gives a polynomial h(x) (and vice versa). The mapping of coefficients is linear, and by inverting its matrix one can explicitely write the solution f(x) for any given polynomial h(x). Looking at the resulting coefficients polynomials of small degree one can guess the following expression for the solution (which is then easily verified to hold for general holomorphic h(x)):

\displaystyle f(x) = \frac{1}{2\pi i} \oint_\gamma \frac{h(z)}{2 \text{arctanh}\frac{1}{z}}\frac{\mathrm{d}z}{x-z}

Here, the contour \gamma has to be chosen in such a way that all singularities of the integrand are inside \gamma. For example, one can choose a circle with radius \max 1,|2x|. For non-holomorphic h(x), this formula still works but the choice of integration contour is more complicated.

In the problem which motivated that equation, we had h(x) given for x\in [-1;1] and wanted to obtain f(x) for x\in [-1;1]. In this case, the contour integral formula can be rewritten as

\displaystyle f(x)=\frac{h(x)\log \frac{1+x}{1-x}}{\left(\log \frac{1+x}{1-x}\right)^2 + \pi^2} + \int_{-1}^1 \frac{h(y)}{\left(\log \frac{1+x}{1-x}\right)^2 + \pi^2}\frac{\mathrm{d}y}{x-y}

This integral is assumped to be regularized by taking the Cauchy principal value. Computing this (e.g. numerically) now only requires h(x) for x\in [-1;1], so it is no longer necessary to assume that h(x) can be extended analytically to \mathbb{C}. A similar expression can also be found for f(x) where |x|>1.

If somebody finds another use or application for this pretty result, please drop me a comment below 😉

Written by inordinatum

October 1, 2011 at 8:03 pm