# On the solution of the multidimensional Langer's equation using the proper generalized decomposition method for modeling phase transitions

Hajer Lamari, Amine Ammar, Adrien Leygue, Francisco Chinesta

► To cite this version:

Hajer Lamari, Amine Ammar, Adrien Leygue, Francisco Chinesta. On the solution of the multidimensional Langer's equation using the proper generalized decomposition method for modeling phase transitions. Modelling and Simulation in Materials Science and Engineering, 2012, 20 (1), 10.1088/0965-0393/20/1/015007. hal-01061890

HAL Id: hal-01061890

https://hal.science/hal-01061890

Submitted on 8 Sep 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Science Arts & Métiers (SAM)

is an open access repository that collects the work of Arts et Métiers ParisTech researchers and makes it freely available over the web where possible.

This is an author-deposited version published in: http://sam.ensam.eu
Handle ID: http://hdl.handle.net/10985/8479

## To cite this version :

Hajer LAMARI, Amine AMMAR, Adrien LEYGUE, Francisco CHINESTA - On the solution of the multidimensional Langer's equation using the proper generalized decomposition method for modeling phase transitions - Modelling and simulation in materials science and engineering - Vol. 20, n°1, p.015007 - 2012
# On the solution of the multidimensional Langer's equation using the proper generalized decomposition method for modeling phase transitions

H Lamari¹, A Ammar², A Leygue¹ and F Chinesta¹

¹ EADS Corporate Foundation International Chair, GEM, UMR CNRS-Centrale Nantes,
1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France

² Arts et Métiers ParisTech, 2 Boulevard du Ronceray, BP 93525, F-49035 Angers cedex 01,
France

E-mail: Hajer.Lamari@ec-nantes.fr, Adrien.Leygue@ec-nantes.fr,
Francisco.Chinesta@ec-nantes.fr and Amine.Ammar@ensam.eu

## Abstract

The dynamics of phase transition in a binary mixture occurring during a quench is studied taking into account composition fluctuations by solving Langer's equation in a domain composed of a certain number of micro-domains. The resulting Langer's equation governing the evolution of the distribution function becomes multidimensional. Circumventing the curse of dimensionality the proper generalized decomposition is applied. The influence of the interaction parameter in the vicinity of the critical point is analyzed. First we address the case of a system composed of a single micro-domain in which phase transition occurs by a simple symmetry change. Next, we consider a system composed of two micro-domains in which phase transition occurs by phase separation, with special emphasis on the effect of the Landau free energy non-local term. Finally, some systems consisting of many micro-domains are considered.

## 1. Introduction

In the last few decades considerable attention has been paid to second-order transition by spinodal decomposition or nucleation and growth occurring in a binary mixture quench [1]; and especially to the dynamics of order parameter fluctuations that increase when approaching the critical point. The theoretical basis of these studies is the mean field theory proposed first by Landau to study first and second order phase transitions occurring by symmetry change in magnetic materials [2]. The theory describes the medium in terms of a continuous
order parameter $\eta$ which measures the symmetry degree, the limiting cases are a completely symmetric phase for which $\eta = 0$ and a completely non-symmetric one for which $\eta = 1$.

Landau proposed an expression for the free energy

$$f(T, \eta) = f_0 - \frac{1}{2}\gamma\eta^2 + \frac{1}{4}\eta^4, \qquad (1)$$

where $\gamma$ denotes the transition parameter that depends on the temperature $T$, and $f_0$ measures the free energy in the one phase region for $\eta = 0$.

For binary mixtures, the order parameter is a conservative quantity when considering phase transitions and it measures the difference in volume fraction between the two components $\eta = \phi_A - \phi_B$, the limiting cases are $\eta = -1$ for space region containing only B molecules and $\eta = +1$ for space region containing only A molecules. In the phase diagram the coexistence curve separates the one phase region which occupies the domain of high temperature in which no separation occurs and where the medium is symmetric from the two phases region in which segregation occurs and where the medium is non-symmetric. The two phases region contains the metastable region and the unstable region which have different separation kinetics and which are separated by the spinodal curve.

For binary mixtures, we can derive the bulk free energy $f(\eta)$ from the regular solution model which gives the mixing enthalpy and the mixing entropy terms. We define the interaction parameter or Flory parameter $\chi$, which controls the phase separation. The entropic effects are predominant in the one phase region $\chi \ll \chi_{cr}$ whereas enthalpic effects are dominant in the two phases region $\chi \gg \chi_{cr}$ so that segregation occurs. The bulk free energy reads

$$f(\eta) = \frac{k_B T}{v} \left[ \left(\frac{1-\eta}{2}\right) \ln\left(\frac{1-\eta}{2}\right) + \left(\frac{1+\eta}{2}\right) \ln\left(\frac{1+\eta}{2}\right) + \chi \left(\frac{1-\eta}{2}\right) \left(\frac{1+\eta}{2}\right) \right], \quad (2)$$

where $k_B$ is the Bolztmann constant, and $v$ an elementary volume in which the order parameter is uniform.

A more general expression of the free energy can be given using the general formula of the free energy expansion as sum of powers of the order parameter [3].

For a non-uniform system, the total free energy $F$ is obtained as the integral over the volume of $f$ and an additive non-local term to take into account non-local interactions between molecules in a way to obtain the Landau–Lifshitz functional:

$$F[\eta] = \int_V (f(\eta) + K(\nabla\eta)^2) dV \qquad (3)$$

with $K > 0$ to ensure stability.

The first attempt to understand phase separation kinetics using Landau theory was proposed by Cahn and Hilliard who suggested a diffusion equation for blends in which the drift governing the separation is expressed in terms of the gradient of the chemical potential, obtained by a functional derivative of the Landau free energy [4–7].

$$\frac{\partial \eta}{\partial t} = \nabla \left[ M \nabla \left( \frac{\partial f(\eta)}{\partial \eta} - 2K \nabla^2 \eta \right) \right], \qquad (4)$$

where $M$ denotes the mobility coefficient. The linear Cahn–Hilliard theory describe rigorously the initial stages of spinodal decomposition by giving an expression of the amplification factor and the effective diffusivity coefficient as functions of the wave number of the fluctuation which reproduces accurately experimental measurements. It fails however in reproducing the behavior at very early instants especially in the vicinity of the critical point where the deviation from experimental data becomes important. The weakness that exhibits this equation is that it does not take into account thermal fluctuations, the factor which initiates separation, and which increases when approaching the critical point even for later stages of separation. This follows
from the mean field approximation on which the Cahn-Hilliard model is based. In fact, the order parameter corresponds by definition to the minimum of the free energy assuming that it is unique. This assumption does not hold in the immediate critical point vicinity where the slope of free energy curves decreases and the minimum is no more unique. Thus one cannot neglect order parameter fluctuations, and mean field theory is no longer valid.

A corrective term was added by Cook on the right-hand side of the Cahn-Hilliard equation, the so-called Cahn-Hilliard-Cook equation, of stochastic nature whose stochastic noise amplitude scales with $k_B T$ and $M$ and satisfies the fluctuation dissipation theorem. This equation has shown more reliable results [8].

However, a full description of composition fluctuations based on a statistical description was achieved by Langer [9]. He adopted a discrete description of the medium as a lattice divided into micro-domains which are communicating. Thermal fluctuations and perpetual Brownian motion state that each one exchanges molecules continuously with the other. A correct description of the composition evolution should take into account such random fluctuations within a statistical mechanics framework.

Thus the order parameter being stochastic, it cannot be associated with a deterministic variable. An appealing way to describe it lies on the introduction of a probability distribution function. Langer assigned to each micro-domain $i$ a coordinate $\eta_i$ which fluctuates in a way to obtain a multivariate distribution function $\rho(\{t, \eta\}) = \rho(t, \eta_1, \dots, \eta_d)$ to describe the whole medium. He derived from the chemical master equation a Fokker-Planck type equation governing the evolution of $\rho(t, \{\eta\})$. The interest of this equation is that it is established on the basis of a statistical model of the mechanism by which composition fluctuations are produced by interactions with a thermal reservoir inducing transition from configuration $\{\eta\}$ to configuration $\{\eta'\}$. The Cahn-Hilliard equation is simply the first moment of Langer's equation whereas the Cahn-Hilliard-Cook equation derives from its second moment.

In this paper, we first introduce Langer's equation pointing out its main drawback, its high dimensionality. To circumvent the curse of dimensionality we introduce the proper generalized decomposition (PGD) that is then used for discretizing the high-dimensional Langer's equation. Different scenarios are then treated considering a one-dimensional physical space consisting of a different number of micro-domains.

### 1.1. Langer's equation

We consider a non-uniform incompressible binary mixture of small molecules A and B. Let $\phi_A$ and $\phi_B$ denote the volume fraction of each species. We define the variable $\eta$ to be

$$ \eta = \phi_A - \phi_B. \qquad (5) $$

We adopt a discrete description of the alloy, we divide the physical domain in $d$ micro-regions, $\eta_i$ denotes the value of the variable $\eta$ in the $i$th micro-domain, and we define $\{\eta\} = \{\eta_1, \eta_2, \dots, \eta_d\}$. Thermal fluctuations induce a continuous exchange of molecules, so that the variable $\{\eta\}$ does not take a unique value at equilibrium, but it is rather described in a statistical way by introducing a probability density function $\rho(t, \{\eta\})$. The time evolution of this quantity is governed by Langer's equation which contains a drift induced by the gradient of the chemical potential and a Brownian term accounting for fluctuations. We limit exchanges to immediate neighbor micro-domains. The dimensionless form of Langer's equation reads

$$ \frac{\partial \rho}{\partial t} = \sum_{i=1}^{d} \frac{\partial}{\partial \eta_i} \left[ \left( m_i \frac{\partial F}{\partial \eta_i} - \sum_{j \in C(i)} \frac{\partial F}{\partial \eta_j} \right) \cdot \rho + \frac{1}{2} \sigma^2 \left( m_i \frac{\partial \rho}{\partial \eta_i} - \sum_{j \in C(i)} \frac{\partial \rho}{\partial \eta_j} \right) \right], \qquad (6) $$
where $m_i$ represents the number of micro-domains directly connected to the $i$th micro-domain. For example in the one-dimensional case $m_i = 2$ for $i = 2, \dots, d-1$ and $m_i = 1$ for $i = 1$ and $i = d$. On the other hand, $C(i)$ is the set of index related to such neighbor micro-domains. Thus, for an internal micro-domain $i$ this set results $C(i) = \{i-1, i+1\}$.

According to the discrete description, the free energy integral taking into account that $\eta_i$ is assumed constant in the $i$th micro-domain, $f(\eta_i)$ will be also constant. The non-local term is written using a finite difference formula. Thus, in the one-dimensional case the free energy integral reads

$$F(\{\eta\}) = \mathcal{V} \left( \sum_{i=1}^{d} f(\eta_i) + \sum_{i=1}^{d-1} K \left( \frac{\eta_{i+1} - \eta_i}{h} \right)^2 \right), \quad (7)$$

where $\mathcal{V}$ represents the micro-domains volume assumed all of them of equal volume. The form given to the advective and diffusive fluxes assures the conservation balance. If we define the Langer's equation residual $\mathcal{R}(\rho(t, \{\eta\}))$ from

$$\frac{\partial \rho}{\partial t} = \sum_{i=1}^{d} \frac{\partial}{\partial \eta_i} \left[ \left( m_i \frac{\partial F}{\partial \eta_i} - \sum_{j \in C(i)} \frac{\partial F}{\partial \eta_j} \right) \cdot \rho + \frac{1}{2} \sigma^2 \left( m_i \frac{\partial \rho}{\partial \eta_i} - \sum_{j \in C(i)} \frac{\partial \rho}{\partial \eta_j} \right) \right] \quad (8)$$

Langer's weak form results:

$$\int_{\mathcal{I} \times \Omega} \rho^* \cdot \mathcal{R}(\rho(t, \{\eta\})) = 0, \quad (9)$$

$$\mathcal{R}(\rho(t, \{\eta\})) = \frac{\partial \rho}{\partial t} - \sum_{i=1}^{d} \frac{\partial}{\partial \eta_i} \left[ \left( m_i \frac{\partial F}{\partial \eta_i} - \sum_{j \in C(i)} \frac{\partial F}{\partial \eta_j} \right) \cdot \rho + \frac{1}{2} \sigma^2 \left( m_i \frac{\partial \rho}{\partial \eta_i} - \sum_{j \in C(i)} \frac{\partial \rho}{\partial \eta_j} \right) \right]. \quad (10)$$

Despite its fine description of fluctuation dynamics, Langer's equation is still not yet solved and this is due to the high dimensionality it involves. In fact, each coordinate $\eta_i$ takes values into the interval $[-1, 1]$. Thus, the probability distribution function $\rho(t, \{\eta\})$ is defined $\mathcal{I} \times \Omega$, $\mathcal{I}$ being the time interval considered. When the number of dimensions $d$ increases a mesh of the whole domain $\Omega$ is unaffordable and then, standard discretization techniques inapplicable. We discuss the curse of dimensionality illness in the next section.

## 1.2. Curse of dimensionality of Langer's equation

Problems defined in high-dimensional spaces abound. The direct solution of these problems in scientific computing has long been thought intractable in view of the so-called curse of dimensionality. Consider for example Langer's equation related to a system containing $d$ micro-domains. This equation defines a transient problem to be solved in a space of dimension $d$. A typical grid-based discretization with $M$ nodes for each coordinate would yield a total number of discrete unknowns of order $M^d$. A rather coarse discretization ($M = 10^3$) of a systems containing $d = 30$ micro-domains would thus involve $10^{90}$ unknowns and this is a gigantic number indeed.

In two recent papers [10, 11], we have proposed a technique able to circumvent, or at least alleviate, the curse of dimensionality. This method is based on the use of separated representations. It basically consists in constructing by successive enrichment an approximation of the solution in the form of a finite sum of $N$ functional products involving $d$ functions of each coordinate. In contrast to the shape functions of finite elements
approximations, these individual functions are unknown a priori. They are obtained by introducing the approximate separated representation into the weak formulation of the original problem and solving the resulting non-linear equations. If M nodes are used to discretize each coordinate, the total number of unknowns amounts to N × M × d instead of the M^d degrees of freedom of classical mesh-based discretization techniques. Thus, the complexity of the method grows linearly with the dimension d of the space wherein the problem is defined, in large contrast with the exponential growth of classical mesh-based techniques.

In [10], for example, this new technique has allowed us to compute solutions of the Fokker-Planck equation in configuration spaces of dimension 20 using the multi-bead-FENE spring model of dilute polymer solutions.

The method was later coined PGD, as in many cases the number N of terms in the separated representation needed to obtain an accurate solution is found to be close to that of the optimal decomposition obtained by applying a posteriori the proper orthogonal decomposition to the problem solution.

## 2. The PGD at a glance

Consider a problem defined in a space of dimension d for the unknown field u(x₁, ..., x₉). Here, the coordinates xᵢ denote any usual coordinate (scalar or vectorial) related to physical space, time, or conformation space, for example, but they could also include problem parameters such as boundary conditions or material parameters. We seek a solution for (x₁, ..., x₉) ∈ Ω₁ × ... × Ω₉.

The PGD yields an approximate solution in the separated form

$$u(x_1, \dots, x_d) \approx \sum_{i=1}^{N} F_i^1(x_1) \times \dots \times F_i^d(x_d). \qquad (11)$$

The PGD approximation is thus a sum of N functional products involving each a number d of functions Fᵢⱼ(xⱼ) that are unknown a priori. It is constructed by successive enrichment, whereby each functional product is determined in sequence. At a particular enrichment step n + 1, the functions Fᵢⱼ(xⱼ) are known for i ≤ n from the previous steps, and one must compute the new product involving the d unknown functions Fᵢⱼⱼ(xⱼ). This is achieved by invoking the weak form of the problem under consideration. The resulting discrete system is non-linear, which implies that iterations are needed at each enrichment step. A low-dimensional problem can thus be defined in Ωⱼ for each of the d functions Fᵢⱼⱼ(xⱼ).

If M nodes are used to discretize each coordinate, the total number of PGD unknowns is N × M × d instead of the M^d degrees of freedom involved in standard mesh-based discretizations. Moreover, all numerical experiments carried out to date with the PGD show that the number of terms N required to obtain an accurate solution is not a function of the problem dimension d, but it rather depends on the regularity of the exact solution. The PGD thus avoids the exponential complexity with respect to the problem dimension.

In many applications studied to date, N is found to be as small as a few tens, and in all cases the approximation converges towards the solution associated with the complete tensor product of the approximation bases considered in each Ωⱼ. Thus, we can be confident about the generality of the separated representation (11), but its optimality depends on the solution regularity. When an exact solution of a particular problem can be represented with enough accuracy by a reduced number of functional products, the PGD approximation is optimal. If the solution is a non-separable function for the particular coordinate system used, the PGD solver proceeds to enrich the approximation until all the elements of the functional space are
included, i.e. the $M^d$ functions involved in the full tensor product of the approximation bases in each $\Omega_j$.

Let us now consider in more detail the specific example related to Langer's equation.

### 3. PGD of Langer's equation

In what follows we consider Langer's equation, previously introduced. For the sake of clarity we consider the one-dimensional physical domain containing $d$ micro-domains and we set $\Omega = \bigcup_{i=1}^{i=d} \Omega_i$

Langer's equation reads

$$ \frac{\partial \rho}{\partial t} = \frac{\partial}{\partial \eta_1} (a_1 \cdot \rho) + \frac{\partial}{\partial \eta_2} (a_2 \cdot \rho) + \dots + \frac{\partial}{\partial \eta_d} (a_d \cdot \rho) \\ + \frac{\sigma^2}{2} \left( \frac{\partial^2 \rho}{\partial \eta_1} - 2 \frac{\partial^2 \rho}{\partial \eta_1 \eta_2} + 2 \frac{\partial^2 \rho}{\partial \eta_2} - 2 \frac{\partial^2 \rho}{\partial \eta_2 \eta_3} + \dots + \frac{\partial^2 \rho}{\partial \eta_d} - 2 \frac{\partial^2 \rho}{\partial \eta_{d-1} \eta_d} \right) \quad (12) $$

with

$$ a_i = \begin{cases} \frac{\partial F}{\partial \eta_1} - \frac{\partial F}{\partial \eta_2} & \text{if } i=1, \\ 2 \frac{\partial F}{\partial \eta_i} - \frac{\partial F}{\partial \eta_{i-1}} - \frac{\partial F}{\partial \eta_{i+1}} & \text{if } 1 < i < d, \\ \frac{\partial F}{\partial \eta_d} - \frac{\partial F}{\partial \eta_{d-1}} & \text{if } i=d. \end{cases} \quad (13) $$

Because of the high-dimensionality of Langer's equation, mesh-based discretization techniques fail in solving it. The PGD strategy introduced in the previous sections allows for circumventing the curse of dimensionality. Thus, the solution is searched in the separated form

$$ \rho(t, \eta_1, \dots, \eta_d) \approx \sum_{i=1}^{i=N} F_i^0(t) \cdot \prod_{j=1}^{j=d} F_j^i(\eta_j). \quad (14) $$

The constructor of such separated representation proceeds by successive enrichment steps until reaching the desired accuracy. Thus, if we assume at iteration $n$ the separated representation of the solution given by

$$ \rho^n(t, \eta_1, \dots, \eta_d) = \sum_{i=1}^{i=n} F_i^0(t) \cdot \prod_{j=1}^{j=d} F_j^i(\eta_j) \quad (15) $$

the solution at the next enrichment step $\rho^{n+1}(t, \eta_1, \dots, \eta_d)$ writes

$$ \rho^{n+1}(t, \eta_1, \dots, \eta_d) = \rho^n(t, \eta_1, \dots, \eta_d) + F_{n+1}^0(t) \cdot \prod_{j=1}^{j=d} F_{n+1}^j(\eta_j) \quad (16) $$

requires the calculation of the $d+1$ unknown functions: the time function $F_{n+1}^0(t)$ and the $d$ functions $F_{n+1}^j(\eta_j)$, $j=1, \dots, d$.

The resulting enrichment problem is non-linear and a linearization strategy is compulsory to solve it. The simplest procedure consists in using alternating directions fixed point algorithm that assuming known $d$ of the $d+1$ functions involved in (16), compute the unknown one. The iteration procedure starts by considering arbitrary functions $F_{n+1}^j$, $j=0, \dots, d$ that are
denoted by $F_{n+1}^{j,0}$, $j=0, \dots, d$ and then iterate until reaching convergence according to

* While $\| \sum_{j=0}^{j=d} (F_{n+1}^{j,r} - F_{n+1}^{j,r-1})^2 \| > \epsilon$, with $\epsilon$ a small enough parameter:

* Compute $F_{n+1}^{0,r+1}$ from $[F_{n+1}^{1,r}, \dots, F_{n+1}^{d,r}]$ assumed known.

* $\forall j \in [1, \dots, d-1]$ compute $F_{n+1}^{j,r+1}$ from:

$[F_{n+1}^{0,r+1}, \dots, F_{n+1}^{j-1,r+1}, F_{n+1}^{j+1,r}, \dots, F_{n+1}^{d,r}]$ assumed known.

* Compute $F_{n+1}^{d,r+1}$ from $[F_{n+1}^{0,r+1}, \dots, F_{n+1}^{d-1,r+1}]$ just computed.

* $r \leftarrow r + 1$

Thus, at each iteration of the non-linear solver we should solve $d+1$ one-dimensional linear problems.

### 3.1. Illustrating the separated representation construction

In order to illustrate in more detail the solution of each one of the problems concerned by the procedure just described, we are considering, without loss of generality, a system involving two micro-domains. In this case Langer's equation involves the time and the two configuration coordinates $\eta_1 \in \Omega_1$ and $\eta_2 \in \Omega_2$. Despite the fact that in this case the solution of the resulting two-dimensional transient problem can be obtained by invoking any discretization technique (finite elements, finite differences, finite volumes, etc) we consider the PGD solution because it is the only one able to address systems containing several micro-domains.

Thus, in a system containing two micro-domains Langer's equation reduces to

$$ \frac{\partial \rho}{\partial t} = \frac{\partial}{\partial \eta_1} (a_1 \cdot \rho) + \frac{\partial}{\partial \eta_2} (a_2 \cdot \rho) + \frac{\sigma^2}{2} \left( \frac{\partial^2 \rho}{\partial \eta_1} - 2 \frac{\partial^2 \rho}{\partial \eta_1 \eta_2} + \frac{\partial^2 \rho}{\partial \eta_2} \right) \quad (17) $$

with $(t, \eta_1, \eta_2) \in \mathcal{I} \times \Omega_1 \times \Omega_2 (\mathcal{I} \subset \mathbb{R}, \Omega_1 \subset \mathbb{R}$ and $\Omega_2 \subset \mathbb{R})$.

Now, assuming the first $n$ sums of the separated representation already computed, that is

$$ \rho^n(t, \eta_1, \eta_2) = \sum_{i=1}^{n} F_i^0(t) \cdot F_i^1(\eta_1) \cdot F_i^2(\eta_2) \quad (18) $$

the solution at the next enrichment step $\rho^{n+1}(t, \eta_1, \dots, \eta_d)$ reads

$$ \rho^{n+1}(t, \eta_1, \eta_2) = \rho^n(t, \eta_1, \eta_2) + F_{n+1}^0(t) \cdot F_{n+1}^1(\eta_1) \cdot F_{n+1}^2(\eta_2). \quad (19) $$

For alleviating the notation we introduce the following notation: $T(t) \equiv F_{n+1}^0(t)$, $R(\eta_1) \equiv F_{n+1}^1(\eta_1)$ and $S(\eta_2) \equiv F_{n+1}^2(\eta_2)$.

If we define Langer's equation residual from

$$ \mathcal{R}(\rho(t, \eta_1, \eta_2)) = \frac{\partial \rho}{\partial t} - \frac{\partial}{\partial \eta_1} (a_1 \cdot \rho) - \frac{\partial}{\partial \eta_2} (a_2 \cdot \rho) - \frac{\sigma^2}{2} \left( \frac{\partial^2 \rho}{\partial \eta_1} - 2 \frac{\partial^2 \rho}{\partial \eta_1 \eta_2} + \frac{\partial^2 \rho}{\partial \eta_2} \right) \quad (20) $$

Langer's equation results in

$$ \mathcal{R}(\rho(t, \eta_1, \eta_2)) = 0 \quad (21) $$

Because Langer's equation linearity it follows at the enrichment step at iteration $n+1$:

$$ \mathcal{R}(\rho^{n+1}) = \mathcal{R}(\rho^n) + \mathcal{R}(T \cdot R \cdot S) = 0 \quad (22) $$

that can be rewritten as

$$ \mathcal{R}(T \cdot R \cdot S) = -\mathcal{R}(\rho^n) \quad (23) $$
with

$$ \mathcal{R}(T \cdot R \cdot S) = \frac{dT}{dt} \cdot R \cdot S - \left( \frac{\partial a_1}{\partial \eta_1} + \frac{\partial a_2}{\partial \eta_2} \right) \cdot (T \cdot R \cdot S) - a_1 \cdot T(t) \cdot \frac{dR}{d\eta_1} \cdot S - a_2 \cdot T(t) \cdot R \cdot \frac{dS}{d\eta_2} - \frac{\sigma^2}{2} \cdot \left( T \cdot \frac{d^2 R}{d\eta_1^2} \cdot S - 2T \cdot \frac{dR}{d\eta_1} \cdot \frac{dS}{d\eta_2} + T \cdot R \cdot \frac{d^2 S}{d\eta_2^2} \right). \quad (24) $$

The weak form of equation (23) reads

$$ \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} \rho^* \cdot \mathcal{R}(T \cdot R \cdot S) \, dt \cdot d\eta_1 \cdot d\eta_2 = - \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} \rho^* \cdot \mathcal{R}(\rho^n) \, dt \cdot d\eta_1 \cdot d\eta_2, \quad (25) $$

where the dependence of $T$, $R$ and $S$ on $t$, $\eta_1$ and $\eta_2$, respectively, is omitted for the sake of clarity.

As in our former works we consider the test function $\rho^*$ given by

$$ \rho^* = T^* \cdot R \cdot S + T \cdot R^* \cdot S + T \cdot R \cdot S^*. \quad (26) $$

In what follows we are illustrating the sequential construction of functions $T(t)$, $R(\eta_1)$ and $S(\eta_2)$:

* Calculation of $T(t)$ assuming known $R(\eta_1)$ and $S(\eta_2)$.

In this case the test function reduces to

$$ \rho^* = T^* \cdot R \cdot S \quad (27) $$

that introduced into the weak form gives:

$$ \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} T^* \cdot R \cdot S \cdot \mathcal{R}(T \cdot R \cdot S) \, dt \cdot d\eta_1 \cdot d\eta_2 \\ = - \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} T^* \cdot R \cdot S \cdot \mathcal{R}(\rho^n) \, dt \cdot d\eta_1 \cdot d\eta_2 \quad (28) $$

As all the functions depending on $\eta_1$ and $\eta_2$ are assumed known at the present step, we could integrate the weak form in the domains $\Omega_1$ and $\Omega_2$ to obtain a problem depending on the time coordinate:

$$ \int_{\mathcal{I}} T^* \cdot \left( \alpha \frac{dT}{dt} + \beta T \right) dt = \int_{\mathcal{I}} T^* \cdot f(t) dt, \quad (29) $$

where coefficients $\alpha$ and $\beta$ come from the integrations of the different functions involving the $\eta_1$ and $\eta_2$ coordinates and $f(t)$ comes from all the integrals involved in the right-hand side of equation (28).

This weak form could be integrated using an appropriate discretization technique taking into account its hyperbolic character (e.g. discontinuous Galerkin-DG-, streamline upwind-SU-, etc). Another possibility taking into account that equation (29) works for any test function $T^*$ is to come back to the strong formulation

$$ \alpha \frac{dT}{dt} + \beta T = f(t) \quad (30) $$

and then integrating it using an appropriate technique accounting for its first order (e.g. backward Euler, Runge-Kutta, etc)

* Calculation of $R(\eta_1)$ assuming known $T(t)$ and $S(\eta_2)$.

Now, using the just computed $T(t)$ and the former $S(\eta_2)$ we will compute function $R(\eta_1)$.

In the present case the test function reduces to

$$ \rho^* = T \cdot R^* \cdot S \quad (31) $$
that introduced into the weak form gives:

$$ \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} T \cdot R^* \cdot S \cdot \mathcal{R}(T \cdot R \cdot S) \, dt \cdot d\eta_1 \cdot d\eta_2 \\ = - \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} T \cdot R^* \cdot S \cdot \mathcal{R}(\rho^n) \, dt \cdot d\eta_1 \cdot d\eta_2 \qquad (32) $$

As all the functions depending on $t$ and $\eta_2$ are assumed known at the present step, we could integrate the weak form in the domains $\mathcal{I}$ and $\Omega_2$ to obtain a problem depending on the $\eta_1$ coordinate:

$$ \int_{\Omega_1} R^* \cdot \left( \gamma_1 \frac{d^2 R}{d\eta_1^2} + \mu_1 \frac{dR}{d\eta_1} + \xi_1 R \right) d\eta_1 = \int_{\Omega_1} R^* \cdot g(\eta_1) d\eta_1, \qquad (33) $$

where coefficients $\gamma_1$, $\mu_1$ and $\xi_1$ come from the integrations of the different functions involving the $t$ and $\eta_2$ coordinates and $g(\eta_1)$ comes from all the integrals involved in the right-hand side of equation (32).

This weak form could be integrated using an appropriate discretization technique. Another possibility taking into account that equation (33) works for any test function $R^*$ is to come back to the strong formulation

$$ \gamma_1 \frac{d^2 R}{d\eta_1^2} + \mu_1 \frac{dR}{d\eta_1} + \xi_1 R = g(\eta_1) \qquad (34) $$

and then integrating it by applying any collocation technique accounting for its second order derivatives (e.g. finite differences, pseudo-spectral collocation techniques, etc).

* Calculation of $S(\eta_2)$ assuming known $T(t)$ and $R(\eta_1)$.

Finally from the just updated $T(t)$ and $R(\eta_1)$ we will compute function $S(\eta_2)$. In the present case the test function reduces to

$$ \rho^* = T \cdot R \cdot S^* \qquad (35) $$

that introduced into the weak form gives:

$$ \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} T \cdot R \cdot S^* \cdot \mathcal{R}(T \cdot R \cdot S) \, dt \cdot d\eta_1 \cdot d\eta_2 \\ = \int_{\mathcal{I} \times \Omega_1 \times \Omega_2} T \cdot R \cdot S^* \cdot \mathcal{R}(\rho^n) \, dt \cdot d\eta_1 \cdot d\eta_2 \qquad (36) $$

As all the functions depending on $t$ and $\eta_1$ are assumed known at the present step, we could integrate the weak form in the domains $\mathcal{I}$ and $\Omega_1$ to obtain a problem depending on the $\eta_2$ coordinate:

$$ \int_{\Omega_2} S^* \cdot \left( \gamma_2 \frac{d^2 S}{d\eta_2^2} + \mu_2 \frac{dS}{d\eta_2} + \xi_2 S \right) d\eta_2 = \int_{\Omega_2} S^* \cdot h(\eta_2) d\eta_2, \qquad (37) $$

where coefficients $\gamma_2$, $\mu_2$ and $\xi_2$ come from the integrations of the different functions involving the $t$ and $\eta_1$ coordinates and $h(\eta_2)$ comes from all the integrals involved in the right-hand side of equation (36).

This weak form could be integrated using an appropriate discretization technique. Another possibility taking into account that equation (37) works for any test function $S^*$ is to come back to the strong formulation

$$ \gamma_2 \frac{d^2 S}{d\eta_2^2} + \mu_2 \frac{dS}{d\eta_2} + \xi_2 S = h(\eta_2) \qquad (38) $$

and then integrating it by applying any collocation technique accounting for its second order derivatives (e.g. finite differences, pseudo-spectral collocation techniques, etc).
### 3.2. Parametric Langer's equation in a single domain

Consider Langer's equation with a single domain. In that case the single conformation coordinate $η$ represents an order parameter taking values in the interval $η \in Ω = [0, 1]$. The limit values are associated with a perfect disorder $η = 0$ or a perfect order $η = 1$.

Using Landau's local free energy the resulting balance Langer's equation reads

$$ \frac{\partial \rho}{\partial t} = \frac{\partial}{\partial \eta} ((-\gamma\eta + \eta^3) \cdot \rho) + \frac{\sigma^2}{2} \frac{\partial^2 \rho}{\partial \eta^2}. \quad (39) $$

Classically, this equation should be solved for different values of the parameter $γ$. Thus, a solution must be computed for each value of the parameter $γ$ by solving equation (39). In order to avoid the solution of many problems, one for each value of the model parameter, one could introduce the parameter $γ$ as an extra coordinate in the model, looking for a general solution $ρ(t, η, γ)$ of equation (39). Thus, the resulting multidimensional model must be solved only once, and then the resulting solution can be particularized for each considered value of the parameter $γ$.

The residual of equation (39) $\mathcal{R}(\rho)$ reads

$$ \mathcal{R}(\rho) = \frac{\partial \rho}{\partial t} - (-\gamma + 3\eta^2) \cdot \rho - (-\gamma\eta + \eta^3) \cdot \frac{\partial \rho}{\partial \eta} - \frac{\sigma^2}{2} \frac{\partial^2 \rho}{\partial \eta^2} \quad (40) $$

with $(t, η, γ) = \mathcal{I} \times \Omega \times \Gamma$, where $\mathcal{I} \subset \mathbb{R}$, $\Omega \subset \mathbb{R}$ and $\Gamma \subset \mathbb{R}$.

The weak form related to equation (39) reads

$$ \int_{\mathcal{I} \times \Omega \times \Gamma} \rho^* \cdot \mathcal{R}(\rho) dt \cdot d\eta \cdot d\gamma = 0. \quad (41) $$

The separated representation of the solution at iteration $n$ writes

$$ \rho^n(t, \eta, \gamma) = \sum_{i=1}^{i=n} F_i^0(t) \cdot F_i^1(\eta) \cdot F_i^2(\gamma) \quad (42) $$

from which the following enrichment steps read

$$ \rho^{n+1}(t, \eta, \gamma) = \rho^n(t, \eta, \gamma) + T(t) \cdot R(\eta) \cdot \Upsilon(\gamma). \quad (43) $$

Again, we consider the test function $ρ^*$ given by

$$ \rho^* = T^* \cdot R \cdot \Upsilon + T \cdot R^* \cdot \Upsilon + T \cdot R \cdot \Upsilon^*. \quad (44) $$

Because the linearity of equation (39) the weak form (41) can be written as

$$ \int_{\mathcal{I} \times \Omega \times \Gamma} \rho^* \cdot \mathcal{R}(T \cdot R \cdot \Upsilon) dt \cdot d\eta \cdot d\gamma = - \int_{\mathcal{I} \times \Omega \times \Gamma} \rho^* \cdot \mathcal{R}(\rho^n(t, \eta, \gamma)) dt \cdot d\eta \cdot d\gamma \quad (45) $$

Now, we can follow the procedure described in the previous section for computing sequentially the functions $T(t)$, $R(η)$ and $Υ(γ)$ until reaching the fixed point of the alternating directions iteration:

* Computing $T(t)$ from $R(η)$ and $Υ(γ)$.

In this step the test function reduces to

$$ \rho^* = T^* \cdot R \cdot \Upsilon \quad (46) $$

that introduced in the weak form (45) gives

$$ \int_{\mathcal{I} \times \Omega \times \Gamma} T^* \cdot R \cdot \Upsilon \cdot \mathcal{R}(T \cdot R \cdot \Upsilon) dt \cdot d\eta \cdot d\gamma \\ = - \int_{\mathcal{I} \times \Omega \times \Gamma} T^* \cdot R \cdot \Upsilon \cdot \mathcal{R}(\rho^n(t, \eta, \gamma)) dt \cdot d\eta \cdot d\gamma \quad (47) $$
Because all the functions depending on the coordinates $η$ and $γ$ are assumed known at this step, we can integrate the weak form (47) in the domain $\Omega \times \Gamma$ to obtain

$$ \int_{\mathcal{I}} T^* \cdot \left( \alpha \frac{dT}{dt} + \beta T \right) dt = \int_{\mathcal{I}} T^* f(t) dt \quad (48) $$

or the equivalent strong form:

$$ \alpha \frac{dT}{dt} + \beta T = f(t) \quad (49) $$

* Computing $R(η)$ from $T(t)$ and $Υ(γ)$.
In this step the test function reduces to

$$ \rho^* = T \cdot R^* \cdot \Upsilon \quad (50) $$

that introduced in the weak form (45) gives

$$ \begin{aligned} & \int_{\mathcal{I} \times \Omega \times \Gamma} T \cdot R^* \cdot \Upsilon \cdot \mathcal{R}(T \cdot R \cdot \Upsilon) \, dt \cdot d\eta \cdot d\gamma \\ & = - \int_{\mathcal{I} \times \Omega \times \Gamma} T \cdot R^* \cdot \Upsilon \cdot \mathcal{R}(\rho''(t, \eta, \gamma)) \, dt \cdot d\eta \cdot d\gamma \end{aligned} \quad (51) $$

Because all the functions depending on the coordinates $t$ and $γ$ are assumed known at this step, we can integrate the weak form (51) in the domain $\mathcal{I} \times \Gamma$ to obtain

$$ \int_{\Omega} R^* \cdot \left( \gamma \frac{d^2 R}{d\eta^2} + \mu \frac{dR}{d\eta} + \xi R \right) d\eta = \int_{\Omega} R^* g(\eta) d\eta \quad (52) $$

or the equivalent strong form

$$ \gamma \frac{d^2 R}{d\eta^2} + \mu \frac{dR}{d\eta} + \xi R = g(\eta) \quad (53) $$

* Computing $Υ(γ)$ from $T(t)$ and $R(η)$.
In this step the test function reduces to

$$ \rho^* = T \cdot R \cdot \Upsilon^* \quad (54) $$

that introduced in the weak form (45) gives

$$ \begin{aligned} & \int_{\mathcal{I} \times \Omega \times \Gamma} T \cdot R \cdot \Upsilon^* \cdot \mathcal{R}(T \cdot R \cdot \Upsilon) \, dt \cdot d\eta \cdot d\gamma \\ & = - \int_{\mathcal{I} \times \Omega \times \Gamma} T \cdot R \cdot \Upsilon^* \cdot \mathcal{R}(\rho''(t, \eta, \gamma)) \, dt \cdot d\eta \cdot d\gamma. \end{aligned} \quad (55) $$

Because all the functions depending on the coordinates $t$ and $η$ are assumed known at this step, we can integrate the weak form (55) in the domain $\mathcal{I} \times \Omega$ to obtain

$$ \int_{\Gamma} \Upsilon^* \cdot (\theta + \varphi\gamma) \cdot \Upsilon(\gamma) \, d\gamma = \int_{\Gamma} \Upsilon^* h(\gamma) \, d\gamma \quad (56) $$

or the equivalent strong form:

$$ (\theta + \varphi\gamma) \cdot \Upsilon(\gamma) = h(\gamma). \quad (57) $$

We can note that as the original partial differential equation (39) does not contain derivatives involving the parameter $γ$, the equation allowing to update function $Υ(γ)$ is an algebraic equation.
<figure><img src="image_5.png" /><figcaption>Figure 1. Modes of the separated representation: (a) $F_i^1(\eta \in [0, 1])$, $i = 1, \dots, 10$; (b) $F_i^2(\gamma \in [-0.8, 0.8])$, $i = 1, \dots, 10$; (c) $F_i^0(t \in [0, 10])$, $i = 1, \dots, 10$.</figcaption></figure>

## 4. Numerical results

For all the algorithmic details concerning the practical application of the PGD (residual minimization for non-symmetric differential operators, introduction of boundary conditions, etc) the interested reader can refer to [12, 13] and the references therein. In what follows we are considering three different scenarios: (i) a system composed of a unique micro-domain; (ii) a system consisting of two micro-domains and finally a system composed of many micro-domains defining a high-dimensional model suffering the so-called curse of dimensionality.

### 4.1. Phase transition in a system consisting of a unique micro-domain

First, we consider a single micro-domain in which the order parameter is homogeneous. We assume that the phase transition occurs by a simple symmetry change. The associated parametric Langer's equation and its solution procedure using the PGD was described in detail in section 3.2.

In what follows we consider $(t, \eta, \gamma) \in [0, 10] \times [0, 1] \times [-0.8, 0.8]$ and $\sigma = 0.5$. The distribution function $\rho(t, \eta, \gamma)$ was searched in the separated form

$$ \rho(t, \eta, \gamma) \approx \sum_{i=1}^{i=N} F_i^0(t) \cdot F_i^1(\eta) \cdot F_i^2(\gamma), \quad (58) $$
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Time evolution of the distribution function $\rho(\eta; t_j, \gamma_{cr})$—each curve corresponds to a time instant, (b) Steady-state distribution function $\rho(\eta; t \to \infty, \gamma_j)$ for different values of the interaction parameter $\gamma_j \in [-0.8 \cdots 0.8]$ moving from one phase to two phases.</figcaption></figure>

where the one-dimensional functions $F_i^j$ ($i = 1, \dots, N$; $j = 0, 1, 2$) where approximated from a one-dimensional finite element approximation defined on a uniform mesh consisting of 101 nodes. The nodal values are denoted by $F_i^0(t_j)$, $F_i^1(\eta_j)$ and $F_i^2(\gamma_j)$, with $j = 1, \dots, 101$.

We consider as initial condition:

$$ \rho(t=0, \eta, \gamma) = G^0(t) \cdot G^1(\eta) \cdot G^2(\gamma) \qquad (59) $$

with

$$ G^0(t_j) = \begin{cases} 1 & \text{if } j = 1, \\ 0 & \text{if } j > 1 \end{cases} \qquad (60) $$

the values $G^1(\eta_j)$ taken from a Gaussian distribution with a null mean value and a unit variance, and $G^2(\gamma_j) = 1 \forall j$.

Figures 1(a)–(c) show the 10 first modes $F_i^0(t)$, $F_i^1(\eta)$ and $F_i^2(\gamma)$ in the resulting separated representation (58). The fact that time functions reach a constant value for long times indicates that the steady state is attained.

Figure 2 depicts the steady-state distribution function curves for different values of the parameter $\gamma$ around the critical point. We find that over the critical point, the mean value of the order parameter vanishes and that the fluctuation grows when approaching the critical temperature $T_{cr}$. The bifurcation behavior occurs once under $T_{cr}$ and the order parameter increases to tend to the limiting value +1.

Nevertheless, the results found above do not agree with critical theory predictions. In fact, we expected to have a divergence of the fluctuation for $\gamma = \gamma_{cr}$ and this has not been found. It was expected that $\rho(t \to \infty, \eta, \gamma_{cr})$ have a constant value. From a mathematical point of view, the equation as it has been written, does not predict any divergent behavior because the stationary distribution function is approximately $e^{-\frac{F}{k_B T}}$, so even if the free energy slope decreases around $\eta = 0$, it remains non-null far away from it.

## 4.2. Phase transition by spinodal decomposition in a system composed of two micro-domains

In this section, we consider a system of two micro-domains characterized by a distribution function $\rho(t, \eta_1, \eta_2)$. The associated Langer's equation and its solution procedure was described in section 3.1.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Transient marginal distribution functions $\rho_1(\eta_1, t)$ and $\rho_2(\eta_2, t)$: (a), (b) $\rho_1(\eta_1, t)$ and $\rho_2(\eta_2, t)$ for $\chi < \chi_{cr}$ and $\sigma = 0.1$. (c), (d) $\rho_1(\eta_1, t)$ and $\rho_2(\eta_2, t)$ for $\chi = \chi_{cr}$ and $\sigma = 0.5$. (e), (f) $\rho_1(\eta_1, t)$ and $\rho_2(\eta_2, t)$ for $\chi > \chi_{cr}$ and $\sigma = 0.1$.</figcaption></figure>

We consider first only the local part of the free energy. Figure 3 shows the time evolution of the marginal distribution functions $\rho_1(\eta_1, t)$ and $\rho_2(\eta_2, t)$ for three values of the interaction parameter: $\chi < \chi_{cr}$, $\chi = \chi_{cr}$ and $\chi > \chi_{cr}$. The two functions are identical because of the symmetric role of $\eta_1$ and $\eta_2$. We note also that we reproduce the bifurcation but we do not find any divergent behavior. When plotting $\rho(\eta_1, \eta_2, t)$ we note the complementarity of $\eta_1$ and $\eta_2$.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Time evolution of distribution function $\rho(\eta_1, \eta_2, t)$ during a quench from one phase region to the instable region for $t = 0, 2, 4, 6, 8$ and $16 (\sigma = 0.1)$.</figcaption></figure>

Figure 4 depicts the distribution function at different times for $\chi \gg \chi_{cr}$. Mass conservation has been checked. The distribution evolves from the symmetric distribution localized around the $(\eta_1, \eta_2) = (0, 0)$ towards the strong segregated distribution localized around $(\eta_1, \eta_2) = (-1, 1)$ and $(\eta_1, \eta_2) = (1, -1)$. The resulting distribution verifies the mass conservation.

We are now introducing the non-local term $K(\nabla\eta)^2$ in the free energy expression. We consider different values: $K = 0, 0.025, 0.05, 0.075, 0.1$. We notice in figure 5 that the degree of segregation reduces as $K$ increases. The segregation is canceled for high values of $K$. We find here a classical result which estimates that phase separation and interface creation occur
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Steady-state distribution function $\rho(\eta_1, \eta_2)$ in the unstable region ($\chi = 2.5$, $\sigma = 0.1$) for $K = 0, 0.025, 0.05, 0.075, 0.1$. Starting from $K = 0$, for which spinodal decomposition occurs, phase separation goes down as $K$ increases until being canceled for $K = 0.1$.</figcaption></figure>

only when the decrease in the free energy resulting from separation compensate the increase
due to surface energy creation.

4.3. Phase transition by spinodal decomposition in a system consisting of many micro-domains

We consider a system composed of five micro-domains. We solve the transient Langer's equation using the PGD. We consider the following values of the model parameters: $\chi = 2.5$
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Time evolution of the marginal distribution function in a system composed of five micro-domains system ($\chi = 2.5$ and ($\sigma = 0.1$): (a), (b) $\rho_1(\eta_1, t)$ and $\rho_5(\eta_5, t)$. (c), (d) $\rho_2(\eta_2, t)$ and $\rho_4(\eta_4, t)$. (e) $\rho_3(\eta_3, t)$.</figcaption></figure>

and $\sigma = 0.1$. We compute $\rho(t, \eta_1, \eta_2, \eta_3, \eta_4, \eta_5)$. We consider 60 time instants and 101 values uniformly distributed along each domain $\Omega_i = [-1, 1]$.

Figure 6 shows the marginal distribution functions relative to each coordinate. The simulation was also done for 7, 8, 10 and even 20 micro-domains partitions.

**5. Conclusion**

In this work, we have studied phase separation by spinodal decomposition under a quench from one phase region, and we have focused especially on composition fluctuations dynamics
increase when passing through the critical region. This was done by solving Langer's equation and calculating the transient probability distribution function $\rho(t, \{\eta\})$. Because the resulting equation is highly multidimensional, mesh-based discretization techniques fail to solve that equation. PGD allows circumventing the curse of dimensionality by constructing a separated representation of the distribution function. It allows also the introduction of some parameters as extra coordinates. Thus, from the solution of the resulting multidimensional parametric model, we can particularize the solution for each choice of the model parameters. In this paper, we introduced the interaction parameter $\chi$ as an extra coordinate proving the potentialities of the PGD for addressing multidimensional parametric models.

Finally, we analyzed different scenarios ranging from systems composed of a unique micro-domain to the ones involving many micro-domains. The obtained results were in agreement with the expectations, proving the ability of PGD to cope with the modeling of spinodal decomposition within a statistical mechanics framework.

## References

[1] Gunton J D, Sam Miguel M and Sahni P S 1983 *Phase Transitions and Critical Phenomena* vol 8 ed C Domb and J L Lebowitz (New York: Academic) pp 267–99

[2] Landau L and Lifshitz E 1984 *Physique Theorique* vol 5 (Paris: Mir) p 3

[3] De Gennes P G 1979 *Scaling Concepts in Polymer Physics* (Ithaca, NY: Cornell University Press)

[4] Cahn J W and Hilliard J E 1958 Free energy of a nonuniform system: I. Interfacial free energy *J. Chem. Phys.* **28** 258–67

[5] Cahn J W 1959 Free energy of a nonuniform system: II. Thermodynamic basis *J. Chem. Phys.* **30** 1121–4

[6] Cahn J W and Hilliard J E 1959 Free energy of a nonuniform system: III. Nucleation in a two-component incompressible fluid *J. Chem. Phys.* **31** 688–99

[7] Cahn J W 1965 Phase separation by spinodal decomposition in isotropic systems *J. Chem. Phys.* **42** 93–9

[8] Cook H E 1970 Brownian motion in spinodal decomposition *Acta Metall.* **18** 297–306

[9] Langer J S 1971 Theory of spinodal decomposition in alloys *Ann. Phys.* **65** 53–85

[10] Ammar A, Mokdad B, Chinesta F and Keunings R 2006 A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids *J. Non-Newtonian Fluid Mech.* **139** 153–76

[11] Ammar A, Mokdad B, Chinesta F and Keunings R 2007 A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part II: Transient simulation using space-time separated representation *J. Non-Newtonian Fluid Mech.* **144** 98–121

[12] Chinesta F, Ammar A and Cueto E 2010 Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models *Arch. Comput. Methods Eng.* **17** 327–50

[13] Gonzalez D, Ammar A, Chinesta F and Cueto E 2010 Recent advances in the use of separated representations *Int. J. Numer. Methods Eng.* **81** 637–59