## Archive for the ‘**Technical Tricks**’ Category

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

*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 of the following Laplace transform:

. (1)

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

. (2)

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

.

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:

.

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

.

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

The Laplace inverse of a simple pole is just an exponential, . By linearity of the Laplace transform, we can invert each summand individually, and obtain an infinite sum representation for :

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 and . This finally gives the solution as claimed above:

. (2)

Have fun!

## Mean-field solution of the Random Field Ising Model

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 fully connected Ising spins , , subject to local random fields , is:

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

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

As usual, is the inverse temperature of the system.

In order to see the phase transition, we need to consider the free energy , averaged over all random fields :

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:

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

is the probability distribution of the random fields, assumed to be Gaussian with mean zero and variance :

Although the expression above for the free energy is written in terms of , strictly speaking I will show it only for the case where has this Gaussian form. Whether the formula for the free energy also holds for other 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

The self-consistent equation for the magnetization given above has a trivial solution for all values of the parameters . Similar to what happens in the mean-field solution for the standard Ising model, if is sufficiently large, the slope is larger than one and there is another non-trivial solution (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

One of the three parameters can be eliminated by rescaling . For example, we can eliminate and express the phase transition line in terms of the dimensionless variables , , :

Another way is to scale away , setting , .

Then one has the following expression for the phase boundary

With this expression, we see more clearly that there is a well-defined zero-temperature limit. When , we get a transition at

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

More generally, we can easily plot the phase diagram in terms of 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,

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 :

is the partition sum of copies, or *replicas*, of the system, feeling the same random fields . Indexing these copies by greek indices , we get

Since the are independent and Gaussian (with mean 0 and variance ), . Then

We have performed the disorder averaging, so what remains is the partition sum of 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 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 . Thus, this procedure is exact in the thermodynamic limit .

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

and then use

In order to have an exponential scaling homogeneously with the system size, we set :

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

where is the partition sum for all replicas on a single site:

In the thermodynamic limit , the integral in our expression for is dominated by its saddle point (due to the fact that the exponential is linear in ). At the saddle point, by symmetry we can assume all the to be equal. Setting for all , the exponent of the integrand in the above expression for becomes

.

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

where

We see that effectively, the saddle-point location corresponds to a fixed average magnetization , for the replicated spins on a single lattice site, subject to the “Hamiltonian” . Rewriting everything in terms of the magnetization , we get

is fixed by the implicit (self-consistent) equation

Now, to further simplify the result, we can compute the sums over all spin configurations by performing another Hubbard-Stratonovich transformation:

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

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

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

The magnetization is fixed by the self-consistent equation at :

By the substitution , 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!!

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

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 given an initial condition at time . 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

This is the density for the final value at time , given an initial value at time . Stationarity of the process implies that is a function of the time difference only. Define the Laplace transform of the free propagator by

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

is the density for the final value at time , given an initial value at time , *without touching the line between and *. **The upshot of this post is the following simple relationship between the Laplace transforms of the free propagator and of the propagator with absorbing boundary :**

(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 . Its stochastic differential equation is

where is white noise with .

The free propagator is a simple Gaussian

Its Laplace transform is given by:

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

The inverse Laplace transform of this expression is:

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

for some action .

Now, to add an absorbing boundary at we weigh each crossing of by a factor in the path integral, and then let . This effectively retains only those trajectories in the path integral, which never touch . We thus have:

.

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

.

Note that here I cancelled the factor , coming from the Taylor expansion of , against the possibilities of ordering the times . This is why I can assume them to be in ascending order above. Also, the functions actually set all .

Taking Laplace transforms, the time integrals (which are a -fold convolution) simplify to a simple product:

.

Now, one easily obtains the claimed result:

.

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!

## Tricks for inverting a Laplace Transform, part IV: Substitutions

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 so that

, (1)

where , and 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:

.

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

The fact that our Laplace transform is written entirely in terms of the function suggests applying a variable substitution. Instead of performing the Fourier-Mellin integral over , we will integrate over . We then have:

. (2)

Solving the definition of for , we get

and

.

Inserting this into (2), we have

.

The only singularity of the integrand is at , for . We can thus choose the contour to go parallel to the imaginary axis, from to . `Mathematica`

then knows to how to evaluate the resulting integral, giving a complicated expression in terms of Kummer’s confluent hypergeometric function . However, a much simpler expression is obtained if one introduces an auxiliary integral instead:

.

`Mathematica`

knows a simple expression for it:

.

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

.

Simplifying the resulting expressions, one obtains our **final result for **:

(3)

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

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 , and a numerical evaluation of , with given by (3). You can see that they match perfectly!

## A quick introduction to the Martin-Siggia-Rose formalism

**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:

,

where is a deterministic function and is Gaussian noise with an arbitrary correlation function :

.

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

,

where

.

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

,

where is Gaussian white noise. The resulting action for the path integral is local in time:

.

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:

,

where (else the process blows up) and is white noise:

.

In this section, I will show that the **generating functional** for this Ornstein-Uhlenbeck process is given by:

.

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

.

This path integral can be solved exactly. The key observation is that the exponent is linear in , and hence one can integrate over this field. Using , this gives (up to a normalization factor ) a -functional:

.

The factor remains since it is the only term in the action independent of .

The -functional fixes in terms of :

.

From this we get

The normalization factor is seen to be by imposing that for we must have .

Plugging this in gives the generating functional, as claimed above:

.

One simple consequence is the **exponential decay of the connected correlation function**:

.

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 -functional:

.

Now, we rewrite the -functional using the formula . Since the -functional is a product of a -function at all points in space, effectively we are introducing an auxiliary field :

.

Now, the only thing that depends on the noise is the factor . Since is Gaussian, the average of this factor can be simply expressed as

Altogether we obtain

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 .

## Tricks for inverting a Laplace Transform, part III

*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:

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

Now, using

,

and by rescaling , one finds that is solved by:

Have fun!

## An interesting integral equation

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 equation system in the limit of .

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

,

where is known and 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 , which, however, turn to to be competely unrelated to those of ).

One angle of attack to our equation is the observation that plugging in a polynomial gives a polynomial (and vice versa). The mapping of coefficients is linear, and by inverting its matrix one can explicitely write the solution for any given polynomial . 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 ):

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

In the problem which motivated that equation, we had given for and wanted to obtain for . In this case, the contour integral formula can be rewritten as

This integral is assumped to be regularized by taking the Cauchy principal value. Computing this (e.g. numerically) now only requires for , so it is no longer necessary to assume that can be extended analytically to . A similar expression can also be found for where .

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