a limitation of a large-deviation approximation in the 1D Ising model

Note: I used GPT-5.6 Sol to help consolidate my notes and emails into an initial draft of this blog. The underlying work and analysis are my own, and I revised the resulting writing extensively. Some of my original, more organized, AI-free notes are available here.

TLDR: Large-deviation theory can be used to turn a rare-fluctuation problem into an equation of state (EOS) problem. In an interacting 1D chain, a low-order Padé approximation to that equation of state reproduces every limit used to construct it, but becomes nonmonotone at stronger coupling. The resulting rate function is nonconvex, which is thermodynamically forbidden. A Maxwell construction can restore convexity, but in this case it would (probably?) repair an artifact of the approximation rather than highlighting 'interesting physics.'

This note grew out of conversations over a couple of months that I had with Professor Richard Stratt, the Newport Rogers Professor of Chemistry at Brown University. The subject of these conversations started with simplified models of polymer folding, and eventually progressed to using large-deviation theory to infer rare orientational fluctuations. Understanding the motivation for using large-deviation theory is straightforward: important physical events, particularly in glassy liquids, are typically controlled by coordinated rearrangements of molecules. These happen very rarely, and, as such, are hard to observe in simulation or with low-order statistics. Large-deviation theory provides natural tools to study these events using pretty accessible information, eliminating the need for prohibitively long brute-force simulations. The interpolation strategy discussed here is from Mainas and Stratt, J. Chem. Phys. 162, 024501 (2025). I worked through this 1D test case almost entirely, though I consulted Professor Stratt frequently. He also derived these results independently, albeit slightly differently.

The paper’s central idea is neat, and in my opinion, very hard to dislike: combine information that is easy to obtain near equilibrium with analytic information in extreme-field limits, then interpolate between them to reconstruct the full fluctuation distribution.

The method works well in the liquid-crystal systems studied in the paper. The point of this note is to stress-test the generality of the approach by probing simple cases; in the simplest interacting model Professor Stratt and I discussed, the same kind of asymptotic matching can generate a globally impossible EOS even though every imposed local constraint is correct.

the basic statmech picture

A statmech system has an enormous number of microscopic configurations. We usually do not care about which exact configuration the system occupies, but rather about a coarse variable such as the magnetization, density, or, in this note, the extension of a chain. Let’s denote the intensive chain extension by $x$.

Some values of $x$ can be produced by many microscopic configurations, while others require the system to organize itself in a very specific way. So the probability distribution $P_N(x)$ contains thermodynamic information. Near its maximum are the ordinary fluctuations seen all the time in a simulation. Far into its tails are rare, highly organized fluctuations. To study these organized fluctuations, it is useful to rewrite the probability as an effective free energy:

$$F_N(x)=-k_{\mathrm{B}}T\log P_N(x)+C$$

For a system made of $N$ identical “pieces”, the free energy cost of holding the whole system at an atypical value usually grows approximately in proportion to $N$. Large-deviation theory expresses this as

$$P_N(x)\asymp e^{-N I(x)}$$

The rough intuition is that an atypical global value has to be maintained across the whole system: piece 1 must do something weird, piece 2 must also do something weird, etc. For weakly correlated blocks, those probabilities roughly multiply, so the total probability becomes exponentially small in $N$, while the corresponding free energy becomes proportional to $N$.¹1.The symbol $\asymp$ denotes asymptotic equivalence on the exponential scale: $P_N(x)=\exp[-NI(x)+o(N)]$. Equivalently, for $I(x)>0$, $\log P_N(x)\sim-NI(x)$ as $N\to\infty$. The argument in the main text is only intuition.

The function $I(x)$ is called the rate function. It is essentially the free energy cost per particle, measured in units of $k_{\mathrm{B}}T$. Its minimum is the equilibrium value of $x$. Its curvature near that minimum controls ordinary Gaussian fluctuations, while its shape farther away controls rare events.²2.If $x_0$ minimizes $I(x)$, then $I(x)\approx I(x_0)+\frac12I''(x_0)(x-x_0)^2$ nearby. Substituting this into $P_N(x)\asymp e^{-NI(x)}$ gives a Gaussian with variance approximately $[NI''(x_0)]^{-1}$, which is the usual $N^{-1/2}$ central-limit scaling of fluctuations about the mean. For the 1D chain examined below, $x=0$ means no net extension and $x=\pm1$ means every link points in the same direction. A large value of $I(x)$ means that extension is exponentially unlikely in the chain length.

the Gärtner–Ellis theorem

The difficult object is $I(x)$, because measuring it directly requires sampling the rare tails of $P_N(x)$. The Gärtner–Ellis theorem gives a complementary description in terms of a field that biases the system toward the value of $x$ we want to study.

Let $M=Nx$ be the corresponding extensive quantity, and define a dimensionless field $f$. Weighting each configuration by $e^{fM}$ favors configurations with large $M$ when $f>0$ and small $M$ when $f<0$.³3.This reweighting is equivalent to adding the term $-k_{\mathrm{B}}TfM$ to the Hamiltonian, since $e^{-\beta H}e^{fM}=e^{-\beta(H-k_{\mathrm{B}}TfM)}$. At finite $N$, we define the scaled cumulant-generating function $\lambda_N(f)=N^{-1}\log\left\langle e^{fM}\right\rangle_0$, where $\langle\cdot\rangle_0$ denotes an average in the unbiased ensemble. Its thermodynamic limit form is $\lambda(f)=\lim_{N\to\infty}\lambda_N(f)$.

Physically, $\lambda(f)$ is minus $\beta$ times the change in free energy per particle caused by the field. Differentiating it gives the average value selected by that field, i.e. $x(f)=\lambda'(f)$. This identity is standard in statmech and easy to see before taking the large-$N$ limit. Differentiating the logarithm pulls down a factor of $M$, so

$$\lambda_N'(f) = \frac1N \frac{\left\langle M e^{fM}\right\rangle_0} {\left\langle e^{fM}\right\rangle_0} = \langle x\rangle_f$$

The subscript $f$ means an average in the biased ensemble. A second derivative gives $\lambda_N''(f)=\operatorname{Var}_f(M)/N=N\operatorname{Var}_f(x)\ge0$, so $x(f)$ must be nondecreasing. In particular, increasing the field cannot decrease the equilibrium response. Convexity thus follows directly from the fact that variance is nonnegative.⁴4.I am skipping much of the algebraic manipulation here—it is worthwhile to go through it to convince yourself that these identities do hold. The key point is that the numerator and denominator together are just an expectation value in the field-biased ensemble.

Now for the actual theorem! In its simplest form, the Gärtner–Ellis theorem says that if the large-$N$ limit $\lambda(f)$ is sufficiently “regular”—most importantly, differentiable over the relevant range—then the rate function is its Legendre–Fenchel transform:

$$I(x)=\sup_f\left[fx-\lambda(f)\right]$$

At a smooth optimum, differentiating the expression inside the brackets with respect to $f$ gives $x=\lambda'(f)$, and if this relationship can be inverted, the maximizing field is $f(x)$ and

$$f(x)=I'(x), \qquad I(x)-I(0)=\int_0^x f(u)\,du$$

There is a relatively simple physical way to understand this. Without a field, realizing $x$ costs $I(x)$. Applying $f$ rewards that extension by $fx$, so the biased system minimizes $I(x)-fx$. The selected value therefore satisfies $I'(x)=f$.

The theorem itself is not the approximation studied here. For the exact 1D model, it gives the correct convex rate function. The approximation enters later, when we do not know the full function $f(x)$, so we try to reconstruct it from limited information.

the EOS shortcut

Now, near $f=0$, the response of $x$ is roughly linear, i.e. $x(f)\approx x(0)+\chi f$, where $\chi=x'(0)$ is the zero field susceptibility. For a symmetric system (as ours will be and many simplified folding models are), $x(0)=0$. In practice, the weak field slope can be estimated from simulation; in the exactly solvable model below, we will just derive it by expanding $x(f)$ around the origin.

At the opposite extreme, a very large field controls the leading approach to saturation, making that limit analytically pretty simple. The interactions can still enter via subleading corrections. In bounded problems, reaching perfect order generally requires a divergent field, but the leading asymptotic behavior is typically known. This is great news for us!

With this in mind, the strategy used by Mainas and Stratt can be summarized as follows:

1. Obtain the weak field susceptibility from normal, unbiased simulations.
2. Derive the strong field saturation behavior analytically.
3. Construct a minimal interpolation between the found limits.
4. Integrate that EOS to estimate the rate function.

If it works, information from the center of a distribution predicts its tails without directly sampling them or running a dense set of field-biased simulations. This would be awesome—rare events happen... rarely. The question in this post is whether matching the two endpoint regimes is generally sufficient to provide a valid function between them.

the exact 1D test case

Consider a chain of $N$ links, each pointing left or right:

$$s_i\in\{-1,+1\}, \qquad x=\frac1N\sum_{i=1}^Ns_i, \qquad -1\le x\le1$$

The scaled extension $x$ is the total chain displacement per link. Neighboring links interact through the dimensionless coupling $K=\beta J$, and a dimensionless pulling field $f$ couples to the total extension. Positive $K$ makes neighboring links prefer the same direction, so the chain develops correlated domains. Positive $f$ favors right-pointing links. At zero field the two directions remain symmetric and the mean extension is zero, but increasing $K$ makes the chain much easier to polarize; once a small field tips one link, its neighbors prefer to align.

$$\beta H = -K\sum_{i=1}^Ns_is_{i+1} -f\sum_{i=1}^Ns_i$$

This is the nearest-neighbor 1D Ising model in a field.⁵5.Calling this a polymer model is slightly silly, but it has the variable we care about: a bounded extension produced by locally interacting links. More importantly, it is exactly solvable. As many readers likely know, this model is analytically solvable, and a common solution is via the transfer matrix method. A transfer matrix packages the Boltzmann weights of each neighboring pair into a $2\times2$ matrix. Multiplying this matrix along the chain sums over all link configurations, and in the thermodynamic limit its largest eigenvalue determines the free energy. That eigenvalue is

$$\Lambda_+(K,f) = e^K\left[ \cosh f+\sqrt{\sinh^2f+e^{-4K}} \right]$$

so differentiating the corresponding free energy with respect to the field gives the mean extension per link:

$$x(f) = \frac{\partial}{\partial f}\log\Lambda_+(K,f) = \frac{\sinh f} {\sqrt{e^{-4K}+\sinh^2f}}$$

Inverting this relation gives the exact EOS

$$f_{\mathrm{exact}}(x) = \operatorname{arcsinh}\!\left( \frac{e^{-2K}x}{\sqrt{1-x^2}} \right)$$

which we’ll use to test our approximate forms. Its derivative is

$$g_{\mathrm{exact}}(x) \equiv f'_{\mathrm{exact}}(x) = \frac{e^{-2K}} {(1-x^2)\sqrt{1-(1-e^{-4K})x^2}}$$

Every factor is positive for $|x|<1$. Therefore, $f_{\mathrm{exact}}(x)$ is strictly increasing and $I_{\mathrm{exact}}(x)$ is strictly convex at every finite $K$. Another way to see this is by using the fact that $\lambda''(f)$ is a variance, since $g(x)=df/dx=1/(dx/df)=1/\lambda''(f)$, so $g(x)$ must be nonnegative. This condition is important, and we’ll revisit it as such!

constructing the Padé approximation

Now, we use the large-deviation approach of Mainas and Stratt. The field diverges logarithmically as $x\to\pm1$, whereas its derivative has a simple pole. It is therefore convenient to interpolate $g(x)$ instead of $f(x)$.⁶6.In the liquid-crystal problem studied by Mainas and Stratt, the field itself has the pole used in their Padé form. Here $f(x)$ diverges logarithmically, so differentiating turns the endpoint divergence into a simple pole. This is therefore a closely related adaptation of their interpolation strategy, rather than literally the same ansatz.

The symmetry $(x,f)\mapsto(-x,-f)$ makes $f$ odd and $g$ even. The minimal ansatz with the endpoint poles and enough coefficients to match one weak field condition and two strong field conditions is⁷7.A Padé approximant is a rational-function approximation. Here the denominator enforces poles at $x=\pm1$, the numerator is even so that $g(x)$ is even, and the three coefficients are fixed by one weak field and two strong field conditions.⁸8.It would also be reasonable to use two conditions from the weak field expansion and only the leading strong field pole. Mainas and Stratt discuss the analogous alternative in the supplementary material and explain why it is less practical.

$$g_{\mathrm{P}}(x;K) = \frac{A(K)\left[1+B(K)x^2+C(K)x^4\right]} {1-x^2}$$

where the P subscript indicates the Padé approximant. I will also usually suppress the explicit $K$ dependence below.

weak field

Expanding the exact response around $f=0$ gives $x(f)=\chi f+O(f^3)$ with $\chi=e^{2K}$. Thus $f(x)=\chi^{-1}x+O(x^3)$, so $g(0)=\chi^{-1}=e^{-2K}$, which fixes $A=e^{-2K}$.

strong field

The expansions here are a little tricky and took me a while to get right.⁹9.I am skipping a lot of algebra here. To match the strong field behavior, define $\epsilon=1-x$ and $u=e^{-2f}$. Both $\epsilon$ and $u$ approach zero as $f\to+\infty$. Expanding the exact force-extension curve in small $u$ gives

$$\epsilon = 2e^{-4K}u + \left(4e^{-4K}-6e^{-8K}\right)u^2 + O(u^3)$$

Inverting this series to write $u$ in terms of $\epsilon$ gives

$$u = \frac{e^{4K}}2\epsilon - \left( \frac{e^{8K}}2-\frac{3e^{4K}}4 \right)\epsilon^2 + O(\epsilon^3)$$

Since $f=-\frac12\log u$ and $g=df/dx$, this implies

$$g_{\mathrm{exact}}(x) = \frac1{2(1-x)} + \left( \frac34-\frac{e^{4K}}2 \right) + O(1-x)$$

Expanding the Padé ansatz about the same endpoint gives

$$g_{\mathrm{P}}(x) = \frac{A(1+B+C)}{2(1-x)} + \frac{A(1-3B-7C)}4 + O(1-x)$$

Matching the pole and the constant term gives the constraints

$$A(1+B+C)=1, \qquad \frac{A(1-3B-7C)}4 = \frac34-\frac{e^{4K}}2$$

Together with $A=e^{-2K}$, we have 3 equations with 3 unknowns. These conditions give

$$B = -2+\frac52e^{2K}-\frac12e^{6K}, \qquad C = 1-\frac32e^{2K}+\frac12e^{6K}$$

Integrating $g_{\mathrm{P}}(x)$ from the origin, with $f_{\mathrm{P}}(0)=0$, gives

$$f_{\mathrm{P}}(x) = -\frac{Ax}{3}\left[3B+C(3+x^2)\right] + A(1+B+C)\operatorname{arctanh}x$$

The construction is exact at $K=0$, has the exact weak field slope for every $K$, and has the exact first two strong field terms. At small positive coupling it is also numerically very accurate. However, in between, the results are a bit disappointing.

the failure

At $K=0$, the model is noninteracting and the approximation is exact:

The approximation remains extremely accurate at weak positive coupling:

This is encouraging! The approximation has the correct slope at the origin, the correct divergence near saturation, and seems to interpolate smoothly between the two.

By $K=0.5$, the two curves are still close, but differences begin to appear in the interior:

At $K=0.75$, the Padé approximation develops a large backward-bending region:

This cannot be the exact equilibrium EOS.¹⁰10.Thermodynamic stability fixes the sign of the response between conjugate variables. Here the field couples as $-fM$, so increasing $f$ must increase the equilibrium extension: $dx/df\ge0$. For the familiar pressure-volume pair, the sign convention is reversed because $P=-\partial F/\partial V$. Stability then requires $(\partial V/\partial P)_T\le0$: increasing pressure should compress the system. A branch where increasing pressure increases volume is thermodynamically unstable in an ordinary equilibrium system, though related behavior can occur in constrained, metastable, or otherwise unusual systems. The exact $f(x)$ remains monotone, while the approximation predicts a region where increasing the extension requires a smaller field. Equivalently, $g_{\mathrm{P}}(x)=f_{\mathrm{P}}'(x)<0$ over part of the interval, so the inferred rate function satisfies $I_{\mathrm{P}}''(x)<0$ there. In other words, the approximate rate function is nonconvex.

The breakdown becomes more dramatic at stronger coupling:

The interesting part is that the approximation is still correct in every limit used to build it. It has the exact weak field slope and the exact strong field pole and constant term. The failure is entirely in the intermediate region, which the asymptotic constraints do not control.

This is the main result. Matching the local physics at both ends does not guarantee a globally valid EOS between them.

a brief note on the Maxwell construction

A standard response to a nonconvex approximate free energy is to replace it with its convex envelope. In the EOS picture, this is the Maxwell construction: the backward-bending branch is replaced by a constant-field plateau chosen using an equal-area condition.

Equivalently,¹¹11.I don't think it's obvious that this is actually equivalent, it's pretty cool. taking the Legendre–Fenchel supremum globally rather than following every local solution of $f=I'(x)$ discards the backward-bending branch and returns this convex envelope. That would make the approximate rate function convex again. However, I do not think it resolves the underlying issue here. The exact finite-$K$ 1D Ising model has no phase transition, so the Maxwell construction would (I think) just be repairing the Padé approximation rather than recovering the true model behavior.

I did not probe much further into whether there is some interesting physical interpretation of the nonconvex branch. My strong suspicion is that there is not. The exact result remains well behaved, while the pathology appears only after imposing a low-order rational interpolation.

reference

E. Mainas and R. M. Stratt, “Exceptionally large fluctuations in orientational order: The lessons of large-deviation theory for liquid crystalline systems”, J. Chem. Phys. 162, 024501 (2025)

Thanks to Professor Richard Stratt for helpful discussions and guidance with this work. Thanks to Corin Wagen, Moe Zhang, and Raphael Stone for reading drafts of this. Last updated 07/21/2026.