inordinatum

Physics and Mathematics of Disordered Systems

Archive for the ‘Maths’ Category

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!

Advertisements

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

Random matrix theory and the Coulomb gas

with 2 comments

Today I have the pleasure of presenting you a guest post by Ricardo, a good physicist friend of mine in Paris, who is working on random matrix theory. Enjoy!

After writing a nice piece of hardcore physics to my science blog (in Portuguese, I am sorry), Alex asked me to come by and repay the favor. I am happy to write a few lines about the basis of my research in random matrices, and one of the first nice surprises we have while learning the subject.

In this post, I intent to present you some few neat tricks I learned while tinkering with Random Matrix Theory (RMT). It is a pretty vast subject, whose ramifications extend to nuclear physics, information theory, particle physics and, surely, mathematics as a whole. One of the main questions on this subject is: given a matrix M whose entries are taken randomly from a known distribution, what would the distribution for its eigenvalues be?

The applications of such a question are countless. We might be interested in the energy levels of a random hamiltonian, the information transmitted through a random noisy channel, the transmission probability of an electron through an open quantum system; every eigenvalue problem, when the effects of the environment are taken into account, must answer this most fundamental question in RMT. To do so, first we restrict our analysis to symmetric (resp. hermitian) matrices. We have two reasons for it: firstly, most applications in physics involve this kind of matrix (hamiltonians, S-matrices), secondly, they are much easier to analyze, since they can always be diagonalized by a orthogonal (resp. unitary) matrix.

To enable a complete analysis, I must introduce an extra assumption: our matrix must have a probability density function (p.d.f) invariant with respect to orthogonal (resp. unitary) conjugation. This is not such a narrow class of matrices as one might think, we can construct this kind of matrix with nearly any infinitely divisible probability. Roughly, we want to be able to write:

\displaystyle P(M) = P(UDU^{-1}) = P(D),

where D is a diagonal matrix.

However, we should not get too excited with this property, as to write P(D) = P(\{\lambda_1,\ldots ,\lambda_n\}) right away. We are always talking about a probability density function, the computation of the actual probability requires an integral over a certain space, and the passage from the matrix to its eigenvalues requires a variable change, hence we must calculate the Jacobian of such a transformation.

Let us take the simplest and perhaps most useful example: the Gaussian real (resp. complex) symmetric (resp. hermitian) matrix whose entries are i.i.d. normal variables with zero mean and variance 1. We can verify that this matrix is invariant with respect to orthogonal conjugation (actually, we can prove that the Gaussian matrix is the only matrix invariant with respect to orthogonal conjugation and whose entries are independent). Hence, we expect the eigenvalue probability to be:

\displaystyle P(\{\lambda_1,\ldots ,\lambda_n\})=C_N \prod_{j=1}^N e^{-\frac{\beta}{2}\lambda^2}.J(\{\lambda_1,\ldots ,\lambda_n\}).

The value for \beta will be discussed soon, C_N is a normalization constant. I will not derive the Jacobian term and I will refer to [1] for the complete demonstration, not only of this expression, but also of the uniqueness theorem mentioned above. Even though the entries of our matrix are independent, the p.d.f. of its eigenvalues its highly coupled due to the Jacobian:

\displaystyle J(\{\lambda_1,\ldots ,\lambda_n\})=\prod_{j\neq i}|\lambda_i-\lambda_j|^\beta

\beta is one of three possible values: if our initial matrix is real, \beta=1. If it is a complex matrix, \beta=2. If M is a quaternion matrix, \beta=4. This is called Dyson’s threefold way.

Since we have the eigenvalues p.d.f., it should be easy to calculate all quantities involved in a eigenvalue problem. We could be interested in the probability of having an eigenvalue larger than a threshold [2], or the chance of not finding any eigenvalue at all inside an interval. The calculation necessary would be simply the N-dimensional integral of the p.d.f. multiplied by a function, but, oh dear, what a probability function!

The Jacobian will make most variables changes futile, all attempts of uncoupling the equations are faced with the product (\lambda_j-\lambda_i). One standard approach would be to interpret the Jacobian as a Vandermonde determinant (up to a sign), rearrange rows to compose polynomials and adjust them as to create orthogonal functions, whose integral would preserve only the “diagonal” terms of our calculations. I will not explain this technique in detail, even though it is quite beautiful. I will travel another path.

Let us write the whole probability as an exponential. We already have the Gaussian term as an exponential, let us raise the Jacobian to the exponent in the most standard way, we write:

\displaystyle P(\{\lambda_j\})=C_N e^{-\frac{\beta}{2}\sum_{j=1}^N\lambda^2}\prod_{j\neq i}|\lambda_i-\lambda_j|^\beta=C_Ne^{-\beta\mathcal{H}},

where \mathcal{H}=\frac{1}{2}\sum_j \lambda_j^2 - \sum_{i\neq j}\log|\lambda_i-\lambda_j|. If we write C_N=\frac{1}{Z}, the analogy is complete, we can treat this probability density function as a system analyzed by the Boltzmann Canonical Ensemble! Hence, we import all our techniques and tools from statistical mechanics to solve this ugly, nasty piece of probability density function.

It is interesting to remark that this Hamiltonian has a natural physical interpretation. If you recall your electrodynamics course, you might find, deep in your memory, the derivation of the Coulomb potential for two point charges in 2D: if we say their position is \lambda_i and \lambda_j, we have V \propto \log|\lambda_i-\lambda_j|. Our Hamiltonian is the description of a 2D gas under a quadratic potential, or, as we call it, Dyson’s Coulomb gas. We can actually infer the equilibrium position of these electrons, if we consider the quadratic potential and the mutual logarithmic repulsion:

It is no surprise that, for the large N limit, the p.d.f. of the electrons, and hence that of the eigenvalues, is the Wigner Semicircle Law. To deduce it, the next natural steps are (I am not going to detail them any further): replace \lambda_j by its density: \rho(x)=\frac{1}{N}\sum_j \delta(x-\lambda_j), write the Hamiltonian in terms of the density, replace the integral over the eigenvalues by a functional integral over the possible densities \rho, use a saddle-point approximation to argue that only the “best” density will matter for our calculations, derive functionally the Hamiltonian and solve a Tricomi equation to obtain the optimal value for \rho, the equilibrium position of the electrons. Surely, we suppose that N is large for the saddle-point approximation, results for small N require a different approach.

To wrap it up, we can avoid treating a highly coupled integral by sending the nasty coupled part to the exponent, interpret the system as a Canonical Ensemble, and treat it in the thermodynamical limit to find the density of eigenvalues. This beautiful treatment was popularized by Dyson and revisited many times in the last decade. RMT is a very active subject, for physicists and mathematicians,  and I hope that, with these few lines, you have been able to see a little but of what I see on it, and to appreciate, for a change, a place where physics has come to help mathematics on a hard day of work.

[1] M. L. Mehta, Random Matrices, Third Edition, Elsevier, 2004.

[2] S. N. Majumdar, C. Nadal, A. Scardicchio, P. Vivo, How many eigenvalues of a Gaussian random matrix are positive?, Phys. Rev. E 83, 041105 (2011)

Written by Ricardo Marino

October 19, 2012 at 11:42 am

The Alessandro-Beatrice-Bertotti-Montorsi model

with 2 comments

When a magnet is submitted to a slowly varying external magnetic field, its magnetization changes not smoothly, but in discrete jumps. These avalanches can be made audible using an induction coil. The resulting crackling signal is called Barkhausen noise. By analysing various features of this signal one can deduce information on material properties, for example residual stress or defect sizes, which is important for applications such as non-destructive testing. In this post, I will discuss a simple model describing the physics of Barkhausen noise. I will explain some of its predictions, including the stationary signal distribution, and sizes and durations of the avalanches.

Stochastic differential equation of the ABBM model

Example of signal generated by ABBM model

Example of signal generated by ABBM model

As you probably know, a ferromagnetic material with zero net magnetization consists of many magnetic domains. Inside one domain, the spins are oriented in parallel (thus each domain has a non-vanishing magnetization), however the magnetizations of different domains are randomly aligned and cancel out on average.
We will be interested in so-called soft magnets. In these materials domain walls can move quite freely, until they encounter a defect. This means they have a wide hysteresis loop. The dominant mechanism for magnetization is the motion of domain walls (and not the changing of magnetization inside one domain, as for hard magnets).

Alessandro, Beatrice, Bertotti and Montorsi (for details, see ref. [1]) model the change in magnetization under an external field through the motion of a single domain wall transverse to the magnetic field. They propose the following stochastic differential equation for the domain wall position u(t):

\displaystyle \Gamma \partial_t u(t) = H(t) -k u(t) + F(u(t))

Here, H(t) is the external applied magnetic field, and \Gamma is the domain wall relaxation time. The term -ku(t) is the so-called demagnetizing field, which keeps the domain wall from moving indefinitely. F(u) is a random pinning force, which depends on the position of the particle u(t). In the ABBM model, it is assumed to be a realization of Brownian motion. For more details on the motivation of the individual contributions, see e.g. the review in ref. [2].

Note that the random pinning force is quenched, i.e. depends on the particle position u(t) and not directly on the time t. A time-dependent random force would be a model for thermal noise (instead of localized defects).

Simulating the stochastic differential equation above yields a trajectory (see the figure on the right) which is very similar to the results of Barkhausen noise measurements. Due to the specific properties of Brownian motion, the ABBM model is easy to treat analytically. I will now discuss several observables which can be computed analytically: The stationary distribution of domain wall velocities, and the distributions of avalanche sizes and durations.

Stationary domain wall velocity distribution

To obtain a stationary state of domain wall motion, one ramps up the external field H linearly:
\displaystyle H(t) = c t.
Then, the instantaneous domain wall velocity \dot{u} has a stationary distribution given by
\displaystyle P(\dot{u}) = \frac{e^{-\dot{u}}\dot{u}^{-1+c}}{\Gamma(c)}.
Here I use dimensionless units k = \Gamma = 1. One way to derive this result is by solving the Fokker-Planck equation associated to the SDE above (as was done by ABBM in ref. [1]).
This distribution is interesting since it exhibits two different kinds of behaviour: For c < 1, P(\dot{u}=0)=\infty, meaning that the domain wall is pinned at near zero velocity most of the time. On the other hand, for c > 1, P(\dot{u}=0)=0 and the motion is smooth.
In the following I will focus on the case c<1, where we have intermittent avalanches.

Avalanche statistics

One way to obtain information on avalanches in the ABBM model is mapping the SDE above to a path integral (the Martin-Siggia-Rose formalism). This is done e.g. in references [3] and [4] below, where the resulting path integral is solved for any external field H(t). Probably the simplest way to define an avalanche is to apply a step-like field, H(t) = w \theta(t). The instantaneous increase at t=0 from H=0 to H=w triggers precisely one avalanche. Its size S and duration T are distributed according to (again in dimensionless units):
\displaystyle P(S) = \frac{w}{\sqrt{4\pi}S^{3/2}}e^{-\frac{(w-S)^2}{4S}},
\displaystyle P(T) = \frac{w}{(2\sinh T/2)^2}e^{-\frac{w}{e^T-1}}.
For small w and small avalanches, these distributions display a large power-law regime where one has
P(S) \sim S^{-3/2}, \quad P(T) \sim T^{-2}.
These power laws indicate that avalanches in the ABBM model are scale-free: There are both extremely small and extremely large ones (between the microscopic cutoff scale given by w and the macroscopic cutoff scale given by k=1).

Universality (or the lack of it)

The exponents of the power-law regimes in P(S) and P(T) above are universal for mean-field models of elastic interfaces. They do not depend on material properties or on details of the dynamics, but only on the fact that one has sufficiently long range interactions between different parts of the interface. These exponents are well-verified experimentally for magnets falling into the mean-field universality class.

Being universal, they also apply to other elastic interfaces with long-range interactions: Some even argue that the P(S) \sim S^{-3/2} behaviour is related to the Gutenberg-Richter distribution of earthquake moments. The avalanche size S would correspond to the earthquake moment, related to its magnitude M (the number reported in newspapers) via M \propto \log_{10} S. The exponent 3/2 for P(S) would give a Gutenberg-Richter b-value of b\approx 0.75, which is not too far off from the observed one.

On the other hand, I find it a little overly simplistic to try and find universal aspects of completely disparate physical systems. We know after all, that motion of magnetic domain walls and of earthquakes is not the same thing — so maybe the more interesting physics is in their differences, rather than their similarities.
A more detailed analysis of avalanches thus requires going beyond just power-law exponents. Several more refined observables — like mean shapes of avalanches — have been proposed to that end. It has been shown (see ref. [5]) that they are sensitive to the details of the dynamics. In my view, the interesting question (which is still not completely answered) is: What features should one look at, in order to determine if a signal is Barkhausen noise or something else? What can one learn from it about the microscopic disorder in one particular sample by listening to the Barkhausen noise it emits?

Outlook and References

If there is interest, in the future I may extend this blog post to a Wikipedia article, since I believe the model is simple but frequently used. It is still a field of active research, thus the list of references is certainly incomplete. Let me know if you want to add anything!

[1] B. Alessandro, C. Beatrice, G. Bertotti, and A. Montorsi, “Domain-wall dynamics and Barkhausen effect in metallic ferromagnetic materials. I. Theory,” J. Appl. Phys., vol. 68, no. 6, p. 2901, 1990.

[2] F. Colaiori, “Exactly solvable model of avalanches dynamics for Barkhausen crackling noise,” Adv. Phys., vol. 57, no. 4, p. 287, Jul. 2008.

[3] P. Le Doussal and K. J. Wiese, “Distribution of velocities in an avalanche,” Europhys. Lett., vol. 97, no. 4, p. 46004, Apr. 2012.

[4] A. Dobrinevski, P. Le Doussal, and K. J. Wiese, “Nonstationary dynamics of the Alessandro-Beatrice-Bertotti-Montorsi model,” Phys. Rev. E, vol. 85, no. 3, p. 18, Mar. 2012.

[5] S. Zapperi, C. Castellano, F. Colaiori, and G. Durin, “Signature of effective mass in crackling-noise asymmetry,” Nature Physics, vol. 1, no. 1, p. 46, Oct. 2005.

Written by inordinatum

October 11, 2012 at 8:37 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