Physics and Mathematics of Disordered Systems

Archive for the ‘Directed Polymers’ Category

Directed Polymers on hierarchical lattices – Tails of the free energy distribution

leave a comment »

A directed polymer in a random medium is a popular model for studying how an elastic interface in a random energy landscape becomes rough on large scales. On a cartesian lattice in 1+1 dimensions, the roughness exponent is known exactly (see e.g. my earlier introduction to directed polymers) and, in some cases, even the distribution of the free energy can be computed (see my earlier post on simulating the Tracy-Widom distribution). In higher dimensions, analytical results on the scaling exponents of directed polymers are scarce. However, numerical simulations indicate that non-trivial scaling exponents persist above 1+1 dimensions (see ref. [1] below and many more).

Instead of 1+d dimensional cartesian space, one can consider a directed polymer on any other lattice. In particular, hierarchical lattices (as introduced e.g. in ref. [2] below) turn out to be interesting. They exhibit some of the phenomena known from directed polymers on cartesian lattices (e.g. a phase transition between a rough and a high-temperature phase for sufficiently high dimensions), and can be approached analytically. In the following, I’ll show (following ref. [2]) the exact recursion relation for the free energy distribution of a directed polymer. I’ll then deduce some exponent relations for the stretched exponential tails of the free energy.


Figure 1: Construction of a hierarchical lattice with branching ratio b

Free energy of a directed polymer on a hierarchical lattice

Consider an hierarchical lattice with branching ratio b, constructed recursively following Figure 1. At each step of the recursion, every bond of the lattice is replaced by b independent branches consisting of two independent bonds each. The branching ratio b plays a role similar to the dimension of cartesian space – as ref. [2] shows, it controls the scaling exponents and the existence of phase transitions.

Let us denote by f^{(L)} the free energy of a directed polymer spanned between the bottom and the top points of the hierarchical lattice of size L. At zero temperature, the polymer will choose the least-energy path between the bottom and the top point. Since the lattice of size 2L consists of b independent branches, f^{(2L)} is the maximum of the b branches, each consisting of two independent sub-systems of size f^{(L)}. Hence,

\displaystyle f^{(2L)} = \text{max}_{j=1...b} \left(f^{(L)}_{j,1}+f^{(L)}_{j,2}\right).

Let us now denote by F the cumulative probability distribution of f, i.e. F(x) = P(f \leq x). Then the recursion above translates into

\displaystyle F^{(2L)}(x) = \left[\int \mathrm{d}y \, F'^{(L)}(y) F^{(L)}(x-y)\right]^b.

To see this, note that \int \mathrm{d}y \, F'(y) F(x-y) is the cumulative probability distribution of f_1 + f_2 and the maximum of b independent branches is equivalent to taking the b-th power of the cumulative distribution function.

We expect that for large L, the distribution F should converge to a fixed point of the form

\displaystyle F^{(L)}(x)=F\left(\frac{x-L c}{L^\theta}\right),

where c is the typical free energy per unit length and L^\theta is the typical scale for free energy fluctuations. \theta is the scaling exponent for the anomalous free energy fluctuations which we are interested in.

Inserting this scaling form into the recursion relation for F^{(2L)}, some simple transformations yield the following nonlinear integral equation for the fixed point F:

\displaystyle F(px)=\left[\int \mathrm{d}y\, F'(y) F(x-y)\right]^b.


In contrast to the previous recursion relation, the F on both sides are now the same function (namely the fixed point distribution), and p = 2^{-\theta} is the scaling factor indicating how the free energy fluctuations increase from one iteration to the next.

Eq. (1) looks astonishingly simple, but finding its fixed points and the corresponding values of p is equivalent to solving exactly the directed polymer problem on an hierarchical lattice. In ref. [2], an expansion around the Gaussian case b=1 is performed. This turns out to be solvable, and the corresponding expansion for p and F can be obtained analytically. Other than that, not much is known on the solutions of eq. (1) – I would be very much interested in hearing of any approaches to solving or approximating it.

Stretched exponential tails

Eq. (1) puts some interesting constraints on the behaviour of the free energy distribution F(x) for very small and very large x. In analogy to the known solution for the directed polymer in 1+d cartesian dimensions, let us assume that in both of these limits the probability distribution is a stretched exponential, i.e.

\displaystyle F(x) \sim \exp \left[-c_1\,(-x)^\alpha\right] for x \to -\infty

\displaystyle 1-F(x) \sim \exp \left[-c_2\,x^\beta\right] for x \to \infty

Inserting these asymptotics into eq. (1), the right-hand side integral can be approximated by Laplace’s method (“saddle-point expansion”) and is dominated in each case by y \approx x/2. In particular, for large x also y is large, so the approximation is consistent. We obtain:

\displaystyle \exp \left[-c_1\,(-px)^\alpha\right] \sim \exp \left[-2b\,c_1\,(-x/2)^\alpha\right] for x \to -\infty

\displaystyle \exp \left[-c_2\,(px)^\beta\right] \sim \exp \left[-2c_2\,(x/2)^\beta\right] for x \to \infty

Equating the exponents, we see (as first noted in ref. [3]) that the exponents of the tails can be expressed in terms of p and b:

\displaystyle \alpha = \frac{\log 2b}{\log 2p} (tail for x \to -\infty)

\displaystyle \beta= \frac{\log 2}{\log 2p} = \frac{1}{1-\theta} (tail for x \to \infty)

As a simple check, note that these equations are fulfilled by the Gaussian limit b=1, \alpha=\beta=2, \theta=1/2.

It might be possible to extend the saddle-point expansion to higher orders, in order to obtain more precise results on the pre-exponential factors (let me know if you manage to!). Nevertheless, it does not appear to put any specific constraints on p (or, equivalently, \theta).

So, a challenge for the inclined reader remains – how is the scaling factor p in eq. (1) fixed in terms of the branching ratio b? While it is possible that determining p is just as hard as finding the actual fixed-point distribution F, I’ve still not given up hope that a simpler approach or approximation might exist…


[1] Marinari, E., Pagnani, A., Parisi, G., & Rácz, Z. (2002). Width distributions and the upper critical dimension of Kardar-Parisi-Zhang interfaces. Physical Review E, 65(2), 026136.

[2] Derrida, B., & Griffiths, R. B. (1989). Directed polymers on disordered hierarchical lattices. EPL (Europhysics Letters), 8(2), 111.

[3] Monthus, C., & Garel, T. (2008). Disorder-dominated phases of random systems: relations between the tail exponents and scaling exponents. J. Stat Mech (2008) P01008.


Written by inordinatum

October 30, 2017 at 8:31 pm

Simulating Directed Polymers and the Tracy-Widom Distribution

with 4 comments

A few posts ago I briefly talked about the problem of a directed polymer in a random medium. In that post, I discussed the scaling behaviour of the polymer fluctuations due to the disordered environment, i.e. its roughness exponent, in 2 dimensions.

As with many other statistical physics problems, for gaining an intuitive understanding it is also useful to perform some numerical simulations. In this post I will show how the directed polymer partition sum can be simulated using Wolfram Mathematica. I will also compare the results to the known exact result for the distribution of the free energy, the Tracy-Widom distribution.

The transfer-matrix algorithm

The partition sum Z of the directed polymer is defined as the sum over all directed paths \Gamma of the Boltzmann weights e^{-\beta E_\Gamma} of each path. The energy E_\Gamma of the path \Gamma is the sum over the random energies E_i of the sites i which are touched by \Gamma. \beta = \frac{1}{T} is the inverse temperature.

In any space dimension, this partition sum can be evaluated rapidly by using a transfer-matrix algorithm. Let us order the lattice along the longitudinal direction of the directed polymer t. Then, the partition sum at site x, t is obtained from the partition sums at this site’s neighbours y\wedge x in the preceding slice t-1:

\displaystyle Z(x,t) = e^{-\beta E_{x,t}}\sum_{y \wedge x} Z(y,t-1).

For the free energy F(x,t) := -\frac{1}{\beta}\ln Z(x,t), we obtain

\displaystyle F(x,t) = E_{x,t} - \frac{1}{\beta} \log \sum_{y \wedge x} e^{-\beta F(y,t-1)}.

Let us focus on the case of very low temperatures, \beta \to \infty. Then the logarithm is dominated by the neighbour y with the minimal (i.e. most negative) free energy F(y,t-1). Hence we obtain

\displaystyle F(x,t) = E_{x,t} + \min_{y \wedge x} F(y,t-1).

This also makes sense intuitively: At zero temperature, the partition sum of the directed polymer is dominated by the minimal-energy configuration. The “cost” of getting to site x is the cheapest “cost” of getting to one of its neighbours, augmented by the (fixed) cost E(x,t) of going to site x.

The Mathematica simulation code

The following short Mathematica routine returns the free energy at the final time t for a directed polymer

fe = 
  Compile[{{wd, _Integer}, {len, _Integer}}, 
   Module[{pot1, pot2, fval}, 
    pot1 = RandomVariate[NormalDistribution[], {len, wd}];
    pot2 = RandomVariate[NormalDistribution[], {len, wd}];
    fval = ConstantArray[0., wd];
    Do[fval = (fval - (fval - RotateLeft[fval]) UnitStep[
           fval - RotateLeft[fval]]) + pot1[[k]]; 
     fval = (fval - (fval - RotateRight[fval]) UnitStep[
           fval - RotateRight[fval]]) + pot2[[k]];, {k, 1, len}]; 

k is the longitudinal coordinate of the polymer, and runs from 1...len. The vector fval has wd components and contains the free energy values at the current value of k, for all transversal coordinates x=1...wd.
For each k I perform two steps forward, one where site i neighbours on i+1 and i, and one where site i neighbours on i and i-1. This ensures that in the end one obtains a standard square lattice of size 2len\timeswd, where the directed polymer can choose between two symmetric alternatives at each time step. pot1 and pot2 are len\timeswd matrices containing the random energies for each site, for the even and odd values of k, respectively.

Since Mathematica‘s built-in Min function does not operate componentwise on vectors, I emulate it by using the componentwise UnitStep function:

a-(a-b)UnitStep[a-b] = Min[a,b].

At t=0, I use a free initial condition: F(x,t=0)=0 \Leftrightarrow Z(x,t=0)=1, i.e. the polymer end at t=0 can be anywhere. Modifying the line

fval = ConstantArray[0., wd];,

one can change this e.g. to a fixed initial condition where the polymer end at t=0 is fixed.
The returned value is the vector F(x,t) where x=1...wd and t=2len is the final time.

For example, to generate 10^4 samples of the free energy for a directed polymer on a 256\times 256 lattice, with the final end fixed, and the initial end free, you would call

Table[fe[256,128][[1]], {10000}]

Note that the [[1]] picks out the free energy for a fixed final site. This simulation takes about a minute on my laptop.

The Tracy-Widom distribution of the directed polymer free energy

Theory predicts the following form for the free energy F of the directed polymer model in 2 dimensions:

\displaystyle F(t) = ct - A t^{1/3}\chi + ... ,

where ... indicates terms subdominant in t. c is a constant and c t is the nonuniversal, deterministic part of the free energy. The fluctuations of the free energy due to randomness of the disorder, described by the second term, are universal: The random variable \chi is distributed according to a Tracy-Widom distribution, and only its amplitude A depends on discretisation details.

The precise distribution depends on the symmetries imposed by the initial condition.

  • For a polymer free at t=0, P(\chi \leq s) = F_2(s), where F_2(s) is the Tracy-Widom distribution for the Gaussian Unitary Ensemble (GUE).
  • For a polymer fixed at t=0, P(\chi \leq s) = F_1(s), where F_1(s) is the Tracy-Widom distribution for the Gaussian Orthogonal Ensemble (GOE).

The functions F_1 and F_2 are solutions of Painlevé differential equations. Tables of their numerical values can be found with very high precision e.g. on the web page of M. Prähofer.

Tracy-Widom distribution and directed polymer free energy

Tracy-Widom distribution and directed polymer free energy.

Blue line: Tracy-Widom distribution F_1 of the GOE ensemble.
Red line: Tracy-Widom distribution F_2 of the GUE ensemble.
Yellow line: Gaussian distribution.
Histogram: 10^5 simulations of a 256\times 256 directed polymer with free (blue) and fixed (red) initial condition, respectively.

All distributions are normalized so that mean=0 and variance=1.

In the figure on the right you can see how this theoretical prediction compares to numerical simulations using the code above. To get rid of the nonuniversal constants c and A, I consider the “normalized” free energy

f_n := \left(F-\left\langle F \right\rangle\right)/\sqrt{\left\langle F^2 \right\rangle -\left\langle F \right\rangle^2}.

We see that its distribution observed in directed polymer simulations has a skewed shape, which agrees very well with the Tracy-Widom prediction (clearly much better than with the symmetric Gaussian distribution). In the tails, one can even “almost” distinguish the two universality classes corresponding to different initial conditions. For the difference to be really significant, one would need to run some longer simulations though.

Further optimizations

Naturally one would like to make the simulation even faster, in order to have better statistics and reach higher system sizes. One may think that implementing the code above in a low-level programming language (e.g. C++) would yield a large performance gain. However, it turns out that a naive C++ implementation (where the operations on the free-energy vector are implemented by a simple loop over the transversal coordinate x) is actually about three times slower than the Mathematica code above. This is since such array-wise operations can be optimized by using modern processor SIMD extensions, which are utilized in Mathematica‘s routines.

Of course, one can also take advantage of this in C++. One can either hand-optimize the loop over transversal coordinates using the processor-dependent SIMD routines, or use an array math library. For example, I implemented the transfer-matrix algorithm above in C++ using the vector routines from Intel’s Math Kernel Library, and was able to achieve a speed-up by a factor of five compared to the Mathematica code. Still, considering the ease of programming and debugging, I find that the simple Mathematica implementation fares amazingly well!

If you’re interested in more details, references, higher-dimensional versions, the C++ code, or anything related, go ahead and drop me a comment below — I’ll extend this post if there is interest 😉

Written by inordinatum

January 22, 2014 at 11:27 pm

A very brief introduction to directed polymers

with 10 comments

A classical problem in the field of disordered systems is an elastic string (or elastic polymer) stretched between two points in a random environment. One is interested in knowing how the polymer geometry (roughness), ground state energy, and excitations change due to the disordered environment. A fundamental question is if the polymer becomes rough for arbitrarily weak disorder, or if there is a phase transition between a weakly disordered (smooth) and a strongly disordered (rough) phase.

In the following I will try explain some basic results on the simple case of a two-dimensional directed polymer (i.e. a “stretched” polymer which does not contain loops or overhangs). I will explain why for large polymer lengths L, moments of the partition sum Z grow as \ln \overline{Z^n} \propto n^3 L. I will show why this indicates a rough polymer with roughness exponent \zeta=\frac{2}{3}. This analysis holds for arbitrarily weak disorder, meaning that in two dimensions disorder is always relevant, and the polymer is always in the rough phase.

Although none of these observations are new, the existing results are often buried in technical literature (there is not even a Wikipedia page describing the directed polymer problem). I hope this self-contained summary will be useful for people interested in a quick introduction to the subject. Some original references where more details can be found are cited at the end of the post.

The directed polymer partition sum, and its moments

We can consider a polymer configuration as a trajectory for a particle moving between two points. We call the polymer directed if the particle always moves from the initial towards the final point, i.e. if the trajectory does not contain loops or overhangs.

Formally, this means that the polymer configuration in 1+d dimensions can be parametrized by giving d transversal coordinates (x_1...x_d)=\vec{x} as a function of the longitudinal coordinate t (distance along the path). The partition sum is then given by a path integral

Z = \int \mathcal{D}[\vec{x}]e^{-\int_0^L \mathrm{d}t \, \left[(\dot{\vec{x}})^2 + \eta(\vec{x},t)\right]}

(\dot{\vec{x}})^2 is the elastic energy term, which is minimized for straight paths. \eta(\vec{x},t) is the energy it costs the polymer to pass through the site (\vec{x},t) (and which is minimized for rough paths picking out the minimal energy sites). We model a disordered environment by taking \eta to be uncorrelated Gaussian white noise:

\displaystyle \overline{\eta(\vec{x}_1,t_1)\eta(\vec{x}_2,t_2)} = 2c\delta^{(d)}(\vec{x}_1-\vec{x}_2)\delta(t_1-t_2)

Applying the Feynman-Kac formula, the partition sum Z(\vec{x},t) satisfies the following stochastic differential equation:

\displaystyle \partial_t Z(\vec{x},t) = \nabla^2 Z(\vec{x},t) +\eta(\vec{x},t) Z(\vec{x},t)

Z is a random quantity, because \eta is random. Since computing the distribution of Z is difficult, we will consider its moments, defined as

\displaystyle \Psi_n(\vec{x}_1...\vec{x}_n,t) := \overline{Z(\vec{x}_1,t)\cdots Z(\vec{x}_n,t)}

The time evolution of \Psi_n is determined by applying the Itô formula to the equation for Z:

\displaystyle \partial_t\Psi_n(\vec{x}_1...\vec{x}_n,t) = \left[\sum_{j=1}^n \nabla_{\vec{x}_j}^2 + 2c\sum_{j<k}\delta(\vec{x}_j-\vec{x}_k)\right] \Psi_n(\vec{x}_1...\vec{x}_n,t) =: H_n \Psi_n(\vec{x}_1...\vec{x}_n,t)

This is, up to an overall sign, the Schrödinger equation for n bosons interacting through pairwise pointlike attraction of strength c. In one dimension, this is called the Lieb-Liniger model. It is easy to determine the ground state of this system exactly, which will allow us to obtain the leading behaviour of Z^n for large t.

Solution in one spatial dimension: Bethe ansatz ground state

In d=1, the ground state of the Hamiltonian H_n defined above has energy E_n = \frac{c^2}{12}n(n^2-1) and is given by

\displaystyle \Psi_n(\vec{x}_1...\vec{x}_n) \propto \exp\left(-\frac{c}{2}\sum_{j<k} |x_j-x_k|\right)

This is a special case of the general Bethe ansatz ubiquitious in modern theoretical physics. Its application to the directed polymer problem, as discussed here, was recognized by Kardar (see references below).

To check that this ansatz for \Psi gives, indeed, an eigenstate, let us take two derivatives with respect to x_j:

\displaystyle \begin{array}{lcl} \partial_{x_j}^2 \Psi_n & = & \partial_{x_j} \left[-\frac{c}{2}\sum_{j \neq k} \text{sgn}(x_j-x_k) \Psi_n\right] \\ & = & \left[-c\sum_{j \neq k} \delta(x_j-x_k)\right] \Psi_n + \left[-\frac{c}{2}\sum_{k \neq j} \text{sgn}(x_j-x_k)\right]^2 \Psi_n \end{array}

Taking the sum over all j we get

\displaystyle \sum_{j=1}^n \partial_{x_j}^2 \Psi_n = -2c\sum_{j < k} \delta(x_j-x_k) \Psi_n + \frac{c^2}{4}\sum_{j=1}^n \left[\sum_{j \neq k} \text{sgn}(x_j-x_k)\right]^2 \Psi_n

To compute the second term, note that it is symmetric under permutations of the x_j. Let us assume (without loss of generality) that x_1 < x_2 < ... < x_n. Then

\displaystyle \begin{array}{lcl} \sum_{j=1}^n \left[\sum_{j \neq k} \text{sgn}(x_j-x_k)\right]^2 & = & \sum_{j=1}^n \left[\sum_{k>j} (-1)+\sum_{k<j} (+1)\right]^2 \\  & = & \sum_{j=1}^n \left[-(n-j) + j-1\right]^2 \\  & = & \sum_{j=1}^n \left(2j-n-1\right)^2 = \frac{n(n^2-1)}{3} \end{array}

Thus, we see that as claimed,

\displaystyle H_n \Psi_n = E_n \Psi_n

with the energy E_n = \frac{c^2}{12}n(n^2-1). To see that this is, indeed, the ground state, is beyond the scope of this post (essentially it would require one to take a finite system volume V, construct all excited states using the Bethe ansatz, and check that they all have higher energy. Since one knows the number of eigenstates in a finite volume, this would exclude any lower-energy state).

Believing that this is, indeed, the ground state our argument shows that for large polymer length L, the moments of the partition sum Z grow as:

\displaystyle \overline{Z^n (L)} \propto e^{\frac{c^2}{12}n(n^2-1) L}

Implications for the string geometry

The asymptotics of \overline{Z^n} derived above is consistent with a free energy that scales with polymer length L as \ln Z = F =: f L^\frac{1}{3}. Furthermore, the distribution P(f) of the rescaled free energy should have a tail decaying as

\displaystyle P(f) \propto e^{-f^{\frac{3}{2}}} \quad \text{for } f \gg 1

To see this, let us compute \overline{Z^n} assuming this form for P(f):

\displaystyle \overline{Z^n} = \int \mathrm{d}f e^{n f L^{\frac{1}{3}}-f^{\frac{3}{2}}}

For large n, we expect f to be large; thus the integral can be approximated by the maximum of the exponent. It occurs at

f \propto \left(n L^{\frac{1}{3}}\right)^2

and thus we get

\displaystyle \ln \overline{Z^n} \propto L n^3,

consistent with the large-n behaviour of the explicit formula for \overline{Z^n} derived above.
Thus, the asymptotics of \overline{Z^n} allowed us to determine the tail of the free energy distribution P(f) and the scaling F = f L^\frac{1}{3}. The scaling of the typical displacement x with polymer length L follows:

\displaystyle L^\frac{1}{3} \sim F \sim \int_0^L \mathrm{d}t (\dot{x})^2 \sim L \left(\frac{x}{L}\right)^2 \Rightarrow x \sim L^\frac{2}{3}

Thus, the polymer is rough in the sense that the typical fluctuations don’t remain bounded, but increase with polymer length as

\displaystyle \overline{\left[x(L)-x(0)\right]^2}\sim L^{2\zeta}

with the so-called roughness exponent \zeta= \frac{2}{3}.

This whole analysis remains valid at an arbitrarily weak coupling strength c. Thus, the polymer is rough for arbitrarily weak disorder (or, equivalently, arbitrarily high temperatures). In order to observe a true phase transition between weak-disorder and strong-disorder behaviour one needs to go to higher dimensions.

Higher dimensions and phase transitions

As we have now seen, in two dimensions (one longitudinal + one transversal) the directed polymer is always in a rough phase. It turns out (but it would take too long to explain here) that the same holds in three (1+2) dimensions. However, starting from one longitudinal + three transversal dimensions, weak disorder is insufficient to make the polymer rough. In other words, there is a high-temperature phase where the moments of Z scale as \ln \overline{Z^n} \propto n^2, indicating Gaussian fluctuations, i.e. purely thermal roughness.


The original argument for obtaining the moments of Z from the Bethe ansatz solution of the Lieb-Liniger model is due to Kardar (see e.g. Nucl. Phys. B 290, 1989 and the references therein). The saddle-point argument for the tail of the free energy distribution is due to Zhang (see e.g. PRB 42, 1990). If you feel some other reference should be included please let me know!

Written by inordinatum

July 29, 2012 at 9:42 pm