Physics and Mathematics of Disordered Systems

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

4 Responses

Subscribe to comments with RSS.

  1. […] 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 […]

  2. Great Post!
    Do you have a clear exposition of the convergence to Tracy Widom law? And about finite size effects?
    Super blog you made!


    January 28, 2018 at 11:00 pm

    • Hi Guillaume,

      regarding the convergence to Tracy-Widom, there are a number of papers studying the short-time (or arbitrary-time) distribution/correlation functions. The mathematics is quite a bit more involved than for the stationary state (you need to work with Fredholm determinants) but nevertheless exact formulae are possible in many cases. Have a look at:
      – Sasamoto, Spohn 2010: “Exact height distributions for the KPZ equation with narrow wedge initial condition”, arXiv:1002.1879
      – Gueudré, Le Doussal, Rosso, Henry, Calabrese 2012, “Short time growth of a KPZ interface with flat initial conditions”, arXiv:1207.7305
      – Le Doussal, Majumdar, Rosso, Schehr 2016: “Exact short-time height distribution in the one-dimensional Kardar-Parisi-Zhang equation and edge fermions at high temperature”, arXiv:1603.03302

      And as for finite-size effects, what exactly do you mean?


      February 5, 2018 at 7:51 pm

      • I guess I meant twice the same thing, that is how does the decreasing error in \chi scales with the size of the random matrix (t) which largest ev (that can be mapped to \chi afaiu) ? For instance, if it is relevant, is the convergence uniform on open intervals (any problem with the tails )?
        I might want to dig more into the literature though, and come back maybe later on.
        Thanks a lot (what is your name ?)!


        February 6, 2018 at 10:02 am

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: