Isotropic finite-difference discretization of stochastic conservation laws preserving detailed balance

TL;DR

This paper introduces an isotropic finite-difference discretization method for stochastic conservation laws that preserves detailed balance, ensuring accurate and faithful simulation of thermally fluctuating systems.

Contribution

The authors develop a finite-difference scheme that maintains detailed balance and isotropy, with improved spatial accuracy and stable temporal integration for stochastic conservation laws.

Findings

Method accurately reproduces analytical results for model B dynamics.

Discretization preserves detailed balance and isotropy on the lattice.

Benchmark tests show excellent agreement with theoretical predictions.

Abstract

The dynamics of thermally fluctuating conserved order parameters are described by stochastic conservation laws. Thermal equilibrium in such systems requires the dissipative and stochastic components of the flux to be related by detailed balance. Preserving this relation in spatial and temporal discretization is necessary to obtain solutions that have fidelity to the continuum. Here, we propose a finite-difference discretization that preserves detailed balance on the lattice, has spatial error that is isotropic to leading order in lattice spacing, and can be integrated accurately in time using a delayed difference method. We benchmark the method for model B dynamics with a $\phi^{4}$ Landau free energy and obtain excellent agreement with analytical results.

Equations

$${\partial }_t\varphi =\nabla \cdot \left(M\nabla \frac{\delta \mathcal{F}}{\delta \varphi }\right)+\nabla \cdot \mathbf{\xi },$$

$$\mathcal{F}=\int d^d\mathbf{x}\left[f\left(\varphi \right)+\frac{1}{2}K(\nabla \varphi (\mathbf{x},t){)}^2\right].$$

$$⟨\mathbf{\xi }(\mathbf{x},t)\mathbf{\xi }({\mathbf{x}}^{\prime },t^{\prime })⟩=2k_{B}TM\mathbf{I}\delta (t-t^{\prime })\delta ({\mathbf{x}}^{\prime }-\mathbf{x}),$$

$$P[\varphi (\mathbf{x},t)]=Z^{-1}\exp \left(-\beta \mathcal{F}\right),$$

$$⟨\varphi (\mathbf{q})\varphi (-\mathbf{q})⟩=\frac{k_{B}T}{A+K{q}^2}.$$

$${\partial }_t\varphi (\mathbf{x},t)=M{\nabla }^2(A-K{\nabla }^2)\varphi (\mathbf{x},t)+\nabla \cdot \mathbf{\xi }(\mathbf{x},t).$$

$${\partial }_t\varphi (\mathbf{x},t)=M{\tilde{\nabla }}^2(A-K{\tilde{\nabla }}^2)\varphi (\mathbf{x},t)+\tilde{\nabla }\cdot \mathbf{\xi }.$$

$$⟨\varphi (\mathbf{q})\varphi (-\mathbf{q})⟩=\left(\frac{{\tilde{\nabla }}_{\mathbf{q}}\cdot {\tilde{\nabla }}_{\mathbf{q}}}{{\tilde{\nabla }}_{\mathbf{q}}^2}\right)\frac{k_{B}T}{A-K{\tilde{\nabla }}_{\mathbf{q}}^2},$$

$$\mathcal{R}(\mathbf{q})=\frac{{\tilde{\nabla }}_{\mathbf{q}}\cdot {\tilde{\nabla }}_{\mathbf{q}}}{{\tilde{\nabla }}_{\mathbf{q}}^2}$$

$$\tilde{\nabla }\varphi (x)=\frac{\varphi (x+\delta x)-\varphi (x-\delta x)}{2\delta x},$$

$${\tilde{\nabla }}^2\varphi (x)=\frac{\varphi (x+\delta x)-2\varphi (x)+\varphi (x-\delta x)}{(\delta x{)}^2}.$$

$$\mathcal{R}(\mathbf{q})=\frac{\cos (2q\delta x)-1}{4\cos (q\delta x)-4}.$$

$${\tilde{\nabla }}^2\varphi =\tilde{\nabla }\cdot \tilde{\nabla }\varphi =\frac{\varphi (x+2\delta x)-2\varphi (x)+\varphi (x-2\delta x)}{4(\delta x{)}^2},$$

$${\tilde{\nabla }}_{\mathbf{q}}^2=\frac{\cos (2q\delta x)-1}{2(\delta x{)}^2}=-\frac{{\sin }^2(q\delta x)}{(\delta x{)}^2}$$

$${\partial }_t\varphi =D{\partial }_x^2\varphi$$

$${\varphi }_i^{n+1}={\varphi }_i^n+\frac{\alpha }{(m+1{)}^2}\left[{\varphi }_{i+m+1}^n+{\varphi }_{i-m-1}^n-2{\varphi }_i^n\right],$$

$$\lambda =1-\frac{4\alpha }{(m+1{)}^2}{\sin }^2\left[\frac{(m+1)q\delta x}{2}\right]$$

$$0\le \alpha \le \frac{(m+1{)}^2}{2}.$$

$$\frac{\partial \varphi }{\partial t}+\frac{\Delta t}{2}\frac{{\partial }^2\varphi }{\partial t^2}=D\left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{(m+1{)}^2\delta x^2}{12}\frac{{\partial }^4\varphi }{\partial x^4}\right).$$

$$\frac{{\partial }^2\varphi }{\partial t^2}=D^2\frac{{\partial }^4\varphi }{\partial x^4}+O({ϵ}^2).$$

$$\frac{\partial \varphi }{\partial t}=D\frac{{\partial }^2\varphi }{\partial x^2}+\frac{D\delta x^2}{2}\frac{{\partial }^4\varphi }{\partial x^4}{\mathcal{I}}^{CD2}(m,\alpha ),$$

$${\mathcal{I}}^{CD2}(m,\alpha )=\frac{(m+1{)}^2}{6}-\alpha .$$

$${\varphi }_i^{n+1}={\varphi }_i^n+\frac{\alpha }{(m+1{)}^2}\left[{\varphi }_{i+m+1}^{n-m}+{\varphi }_{i-m-1}^{n-m}-2{\varphi }_i^{n-m}\right],$$

$${\lambda }^{m+1}-{\lambda }^m+\frac{4\alpha }{(m+1{)}^2}{\sin }^2\left(\frac{(m+1)q\delta x}{2}\right)=0,$$

$$\lambda =\frac{1}{2}\left(1\pm \sqrt{1-4\alpha {\sin }^2\left(q\delta x\right)}\right).$$

$$\alpha \le 1,$$

$$\frac{\partial \varphi }{\partial t}+\frac{\delta t}{2}\frac{{\partial }^2\varphi }{\partial t^2}=D\left(\frac{{\partial }^2\varphi }{\partial x^2}-m\delta t\frac{{\partial }^3\varphi }{\partial t\partial x^2}+\frac{(m+1{)}^2\delta x^2}{12}\frac{{\partial }^4\varphi }{\partial x^4}\right)+O({ϵ}^3)$$

$$\frac{\partial \varphi }{\partial t}=D\frac{{\partial }^2\varphi }{\partial x^2}+\frac{D\delta x^2}{2}\frac{{\partial }^4\varphi }{\partial x^4}{\mathcal{I}}^{\mathrm{Delayed}}(m,\alpha )+O({ϵ}^3),$$

$${\mathcal{I}}^{\mathrm{Delayed}}(m,\alpha )=\frac{(m+1{)}^2}{6}-\alpha (2m+1).$$

$${\partial }_t\varphi =D{\partial }_x^2\varphi +\frac{D\delta x^2}{2}\left[\frac{2}{3}-3\alpha \right]{\partial }_x^4\varphi .$$

Full text

1 Introduction

Thermal fluctuations are an essential part of complex phenomena as diverse as Brownian motion in colloidal suspensions[1, 2], concentration fluctuations in semi-dilute polymer solutions [3, 4, 5], capillary waves at fluctuating interfaces[6, 7, 8], critical dynamics in binary mixtures [9, 10, 11, 12] and pattern formation [13, 14]. In these systems, stochastic partial differential equations in the form of conservation laws are used to describe the time evolution of conserved densities. Examples include the fluctuating Cahn-Hilliard-Cook equation [15, 16, 17], the fluctuating Navier-Stokes equations of Landau and Lifshitz[18, 19], the fluctuating lubrication equation [20, 21] and models of electrohydrodynamic instabilities in electrospinning experiments [22].

The flux, in the conservation laws, which are coarse-grained expressions of microscopically reversible Hamiltonian dynamics, must satisfy an important constraint: the stochastic part of the flux cannot be chosen independently but must be related through a detailed balance condition to the irreversible part of the flux [23]. This relation is necessary to ensure that the coarse-grained dynamics yields a stationary state with a Gibbs distribution, or in other words, that the dynamics expressed by the conservation law is consistent with micro-reversibility.

Discretization in space and time is a necessary step in seeking numerical solutions to stochastic conservation laws. The discretization of the spatial part of the conservation law commonly requires discrete analogs of the vector differential operators - the gradient, divergence, and curl - and of the Laplacian. It is well-known that the discrete differential operators do not always inherit all properties of the continuum operators. In particular, special care is needed to preserve the continuum “div-grad-curl” identities like the vanishing of the discrete curl of a discrete gradient or of the discrete divergence of a discrete curl. The discretization of the stochastic flux is more involved compared to Langevin equations as it is the divergence of a Gaussian random field. For detailed balance in a stochastic conservation law, it is necessary to ensure that the discrete divergence of a discrete gradient is identical to the discrete Laplacian (see below). Common discrete representations of these operators, constructed in the setting of deterministic conservation laws, typically do not satisfy this last property. Their use in stochastic conservation laws results in a violation of the detailed balance condition [24, 25, 26].

Even when the above “mimetic” property of the discrete spatial derivative operators is ensured, detailed balance can be broken by the temporal discretization: the probability distribution in the stationary state may acquire a dependence on the time step (and, hence, on the kinetic coefficients), when detailed balance would explicitly rule out such a dependence. To the best of our knowledge, much work is needed on the temporal discretization of the Cahn-Hilliard-Cook equation to yield a stationary distribution that is independent of the temporal time step [27] so as to minimize the violation of detailed balance at the discrete level, if not to eliminate it altogether.

Here, we propose a finite-difference discretization of a stochastic conservation law following a semi-discretization strategy and illustrate our general results with the specific example of model B [19]. The stochastic partial differential equation is first discretized in space to yield a set of coupled stochastic ordinary differential equations. In this, we use discretizations of the gradient and divergence that ensure isotropic truncation errors to leading order in lattice spacing [28, 29]. We then define the Laplacian to be a composition of these discrete operators, so that the identity $\nabla^{2}=\boldsymbol{\nabla}\cdot\boldsymbol{\nabla}$ is satisfied by construction. The resulting Laplacian operator is negative semi-definite with a trivial null space: the only eigenvector with a zero eigenvalue is the constant. The support of this Laplacian, though, is larger than those of comparable accuracy commonly used in finite-difference methods.

The temporal discretization of the semi-discrete system is completed using a novel delayed difference scheme. A distinct advantage of this integration scheme, deriving from the trivial null space of the Laplacian, is that it does not produce spurious checker-board modes in two- and three-dimensional spaces. The isotropy and trivial null space of the “mimetic” discretization of the Laplacian, together with the delayed time integrator, yields a numerical method that is stable, accurate and efficient. Our numerical results are for two-point correlation functions and order parameter distributions are in excellent agreement with the well-known analytical results for model B.

To contextualize the contributions of this manuscript, we briefly survey previous work on the topic. The necessity of a consistent discretization of the vector differential operators for satisfying in model H was first pointed out by one of the authors in an unpublished report. Similar ideas were expressed in the subsequent work of Delgado et. al.[30] and Garcia et. al.[31] on the finite-volume discretization of the compressible isothermal fluctuating hydrodynamics at nanoscale and the same for stochastic conservation law obtained from a large-volume expansion of the chemical master equation for reacting and diffusing species respectively. Thampi et. al.[32] presented both spectral and finite-volume discretizations of the order parameter equation in model H, using fluctuating discrete kinetic theory to describe momentum conservation. A systematic study of finite-volume discretization schemes preserving detailed balance for a variety of stochastic conservation laws whose evolution is generated by Poisson and dissipation brackets has recently been initiated. Several alternatives to explicit temporal integration has been explored by Torre et. al.[33]. The purpose of this (possibly incomplete) survey is to emphasize that our work here focuses on the combination of finite-difference spatial discretizations, which are both simple and popular, with delayed temporal integrators which mitigate some of the drawbacks of using finite-difference spatial discretizations.

The remainder of the paper is organized as follows. In section 2 we present model B of the Halperin-Hohenberg classification of dynamical critical phenomena [19] and list those analytical results used later in our benchmarking. Section 3 describes the different topological properties of the discrete operators preserving the FDR at lattice level in connection with the non-interacting order parameter dynamics of a system subjected to a single phase bulk free energy potential. Section IV depicts the asynchronous time discretization method which is crucial to achieve the required stability and accuracy for the FDR preserving discrete operators. In section V, we investigate the cases of the interacting order parameter dynamics for a system which is subjected to both single and two phase equilibrium potentials. In all these cases we find an excellent agreement between the present method and the analytical and pseudo-spectral results. In section VI, we conclude with the multidimensional generalization of the present work and draw a comparison of our method to the well-known cell-dynamical method of Oono and Puri[34, 35], showing that isotropic differences and delay-difference integrators provide an independent formulation of going beyond the Oono-Puri method.

2 Model B

Model B is a stochastic partial differential equation for a conserved scalar order parameter field $\phi(\mathbf{x},t)$ whose dynamics is driven by a competition between deterministic thermodynamic forces and stochastic forces of thermal origin [19, 15, 16]. The equation of motion is

[TABLE]

where $M$ is the order parameter mobility and ${\cal F}$, the Landau free energy, is a functional of the order parameter

[TABLE]

The local part of the free energy density is here taken to be $f(\phi)=\frac{1}{2}A\phi^{2}+\frac{1}{4}B\phi^{4}$ , where $A$ can be either positive or negative but $B$ is always positive. The positive coefficient $K$ in the non-local part is related to the energy cost for gradients in the order parameter. The stochastic flux $\boldsymbol{\xi}(\mathbf{x},t)$ is a zero-mean Gaussian random field whose correlation is local in both space and time,

[TABLE]

where $k_{B}$ is the Boltzmann constant, $T$ is the temperature and $\mathbf{I}$ is the identity matrix in the space of Cartesian indices. This fluctuation-dissipation relation for the stochastic flux ensures that the stationary probability distribution is

[TABLE]

where $Z$ is the partition function.

Model B has a Gaussian fixed point, corresponding to the parameters $A,K>0$ and $B=0$. It also has a non-trivial Wilson-Fisher fixed point when $A=0$ and $B,K>0$. In addition, it allows for two-phase coexistence between the phases $\phi=\pm 1$ when $A<0$ and $B,K>0$. The “domain wall” between these two phases is described by the $\phi^{4}$ soliton. Correlation functions in all three cases can be calculated in closed form.

Here, our principal interest is in the Gaussian fixed point. The free energy is quadratic and the order parameter distribution, consequently, is Gaussian. The two-point correlation determines all remaining correlation functions. In Fourier space, it is

[TABLE]

A first check on the accuracy of the discrete numerical method is provided by a comparison with the two-point correlation function. Away from the Gaussian fixed point, a more stringent check is provided by a comparison with the distribution of the order parameter. Below, we use both these checks to validate our numerical methods.

3 Spatial discretization and detailed balance

In this section we discretize Model B in space to show how naive discretizations break detailed balance and how detailed balance can be restored by a redefinition of the discrete Laplacian. This analysis most illustrative without the additional complication of order parameter non-linearity and therefore, we shall restrict ourselves to the Gaussian phase, though the results obtained will be generally applicable. In the Gaussian phase the equation of motion is linear and, for a constant mobility, takes the form

[TABLE]

Let us denote the discrete gradient, divergence and Laplacian by $\tilde{\boldsymbol{\nabla}}$, $\tilde{\boldsymbol{\nabla\cdot}}$ and $\tilde{\nabla}^{2}$ respectively. It follows that the equation of motion of the discretely sampled field is

[TABLE]

It is then a straightforward exercise to show, using the fluctuation-dissipation relation for the random flux, that the two-point correlation function of the Fourier modes of the order parameter is given by

[TABLE]

where the subscripts indicate the Fourier transforms of the respective operators. Comparing with the two-point correlation function of the continuum theory, Eq.4, it is evident that the discrete two-point correlation function will contain the factor

[TABLE]

which, generally, will differ from unity. To ensure this “equilibrium ratio” to be unity for all wave numbers requires that the Fourier transform of the discrete gradient, divergence and Laplacian operators be related exactly as in the continuum.

To illustrate this analysis with a simple example, consider the standard central-difference stencils in one dimension, for which

[TABLE]

[TABLE]

Using Fourier transform of Eq.9 and Eq.10 and using Eq.8 it can be shown that,

[TABLE]

This expression tends to unity as $q$ tends to zero but is less than unity for all other values of $q$ in the first Brillouin zone $|q|\leq\pi$. Therefore, the use of such a discretization will, even for a (hypothetically) perfect temporal integrator, introduce spurious wave number dependence in the two-point correlation. This artifact of the standard discretization has been noted earlier [32].

However, if the Laplacian is defined as

[TABLE]

repeating the above exercise shows that $\mathcal{R}(\mathbf{q})$ is unity for all wave numbers and the semi-discretization ensures that the correlation function approximates that of the continuum and it is free of spurious wave number dependent contributions. The Fourier transform of this Laplacian is

[TABLE]

and is obviously negative semi-definite. The only null eigenvector in the first Brillouin zone is a constant.

4 Explicit time integrators

The use of the Laplacian in Eq.12 in conventional time discretization schemes will lead to reduced overall accuracy, which can be seen by considering a simple but illustrative example of the diffusion equation

[TABLE]

for the scalar field $\phi$. Using central differences and forward Euler for spatial and temporal discretization respectively, the resulting difference equation is

[TABLE]

where $\phi_{i}^{n}\equiv\phi(i\Delta\mathbf{x},n\Delta t)$ and $\alpha=D\,\Delta t/\delta x^{2}$ is the Courant–Friedrichs–Lewy (CFL) number. The Laplacians of Eq.10 and Eq.12 correspond to $m=0$ and $m=1$ respectively. The numerical stability of the difference scheme can be analyzed using the von Neumann method. The amplification factor for Eq.15 is easily obtained to be

[TABLE]

and thus stability, $|\lambda|<1$, requires

[TABLE]

The continuum limit of Eq.15 under the “diffusive scaling” $\delta x\sim O(\epsilon)$ and $\Delta t\sim O(\epsilon^{2})$, correct to $O(\epsilon^{4})$, is

[TABLE]

This equation can be simplified using using the Cauchy-Kowalewsky backward error analysis, where all higher order time and mixed derivatives are estimated by space derivatives obtained using the differential equation itself. For example, by taking a derivative in time of the evolution Eq.18, we can estimate

[TABLE]

Thus, the effective differential equation with an error of $O(\epsilon^{4})$ is

[TABLE]

where the transport coefficient associated with the biharmonic operator is

[TABLE]

Based on this effective differential equation, we can analyze the accuracy for different values of $m$ and thus the effect of a wider stencil on accuracy. The trade-off in using the “mimetic” stencil is now obvious: the wider stencil has a lower accuracy as $|\mathcal{I}^{{\rm CD2}}(1,\alpha)|>|\mathcal{I}^{{\rm CD2}}(0,\alpha)|$ in the parameter range $0\leq\alpha\leq 1/2$ where the schemes are stable.

5 Delayed time integrators

The above shortcoming of the combination of the spatial “mimetic” Laplacian and the temporal explicit integrator can be remedied by using a recently introduced delayed-in-time integration scheme [36]. This scheme is motivated by the observation that computing derivatives on a wider stencil, while using spatial data from earlier times, can dramatically improve both stability and accuracy [36].

The delayed integrator applied to Eq.14 gives

[TABLE]

which should be compared with Eq.15 for $m=0$. The amplification factor for this scheme obeys

[TABLE]

which can be solved for $m=1$ to obtain

[TABLE]

Hence, to satisfy $|\lambda|\leq 1$, we have the stability condition

[TABLE]

which implies a gain in stability compared to $m=0$. However, the gain is less than naive use of wider stencil (Eq.17). Here, we remind that wider stencil leads to better stability but much lower accuracy. However, the delayed scheme removes this problem of lower accuracy associated with wider stencil. This can be seen from the effective differential equation corresponding to the delayed scheme. In order to obtain such effective equation, similar to previous section, we use Cauchy-Kowalewski procedure. We write Eq.21 in differential form using Taylor series as

[TABLE]

and then replacing the time and the mixed derivatives using evolution equation (Eq.24), to obtain the effective differential equation at the leading order as

[TABLE]

where,

[TABLE]

Thus, for $m=1$, Eq.25 implies

[TABLE]

Thus, for $\alpha>0.25$ the delayed scheme has higher accuracy $\mathcal{I}^{{\rm Delayed}}(1,\alpha)$ than the standard CD2 schemes with $\mathcal{I}^{{\rm CD2}}(0,\alpha)$ as well as naive scheme with wider stencil $\mathcal{I}^{{\rm CD2}}(1,\alpha)$. This is illustrated in Fig.1, where the pre-factors $\mathcal{I}^{{\rm CD2}}(m,\alpha)$ and $\mathcal{I}^{{\rm Delayed}}(m,\alpha)$, are plotted with respect to the CFL($\alpha$)

Thus, long time integration can be efficiently performed by the delayed scheme, as it enhances the stability by a factor $2.0$ as compared to that of the usual CD2 scheme, without reduction in accuracy. This is evident from the schematics (Fig.2), which shows that data taken from past requires wider stencil($2\delta x$) and thus, widening of stencil compensates for error due to the time delay.

6 Delayed time integrators in multi-dimensional space

In last section, we have shown that the delayed scheme leads to better stability and accuracy for discrtized diffusion equation in one dimension. However, for the multi-dimensional extension of the delayed scheme,

[TABLE]

the increase in accuracy and stability requires further restrictions on the form of discrete the operators. This can be seen by repeating the analysis of previous section on multi-dimensional scheme. In this case, similar to the previous section, we write Eq.21 in differential form using Taylor series as

[TABLE]

If the discrete Laplacian preserves an isotropic structure at least at the leading order, with $a$ as a stencil dependent constant, i.e,

[TABLE]

the effective differential equation for the discrete analog of diffusion equation can be written as

[TABLE]

Similar to the 1-D case, this equation for the delayed scheme also has better stability and accuracy. However, conventional discrete operators such as central difference operators do not satisfy Eq.29[28, 34, 35]. This can be seen by Taylor series expansion of central difference Laplacian operator in 3-D

[TABLE]

Thus, we use recently proposed lattice differential operators[28, 29], where discrete operators are constructed from lattice kinetic models. In this approach, the basic discrete vector operators viz. Gradient($\tilde{\boldsymbol{\nabla}}$), Divergence($\tilde{\boldsymbol{\nabla}}\cdot$) and Curl($\tilde{\boldsymbol{\nabla}}\wedge$) for a given vector field $\mathbf{\Phi}$ on this lattice are formulated as

[TABLE]

with the set of connecting vectors on stencil as $\mathbf{c}_{i}$ with $i=1,\cdots N$, where, $N$ being the total number of neighbours, and the corresponding wights $w_{i}$ (normalized to one ($\sum_{i}^{N}w_{i}=1$) are chosen so as to satisfy

[TABLE]

with $\Delta_{\alpha\beta\kappa\eta}$ being the fourth order isotropic Kronecker-delta and the connecting vectors are chosen such that $\sum_{i}^{N}w_{i}c_{i\alpha}=0$. These conditions on weights ensure the isotropy of discrete operators up to the leading orders.

Let us consider a 2-D analog of the formulation depicted in[28, 29], where the computational grid is constructed through a sequence of square lattices. For a 2-D grid with connecting vectors, shown in Fig.3, one could explicitly write the discrete gradient operators as

[TABLE]

[TABLE]

Various stencils on which isotropic operators can be written are documented in [28, 29]. An implementation of these isotropic operators in discretization of PDE and SPDE are also used in [32, 37].

6.1 Factorizable sign-definite isotropic Laplacians

With this definition of gradient and Eq.12 as definition of Laplacian, we create an FDT preserving isotropic discrete Laplacian as, $\tilde{\Delta}^{{\rm iso}}=\tilde{\nabla}^{{\rm iso}}\cdot\tilde{\nabla}^{{\rm iso}}$, which allows us to write an FDT preserving discrete space-time representation of model B(Eq.5) as

[TABLE]

Now, following Eq.12 one could formulate the discrete isotropic Laplacian by using $\tilde{\nabla}\cdot\tilde{\nabla}$, which can be expressed explicitly on a grid, shown in Fig.4 as

[TABLE]

which in Fourier domain is a positive quantity

[TABLE]

7 Model B: Harmonic fluctuations

For the sake of simplicity and without loss of generality, we first consider the following simplified case of model B, describing the dynamics of a non-interacting order parameter in single phase equilibrium,

[TABLE]

with $D=A\,M$ being the diffusion coefficient.

We contrast the present approach with the traditional discretization schemes, by performing a long time integration for this setup where the steady state probability distribution can be compared with the Gaussian distribution expected from the continuous model B dynamics. The superiority of the present work is apparent in the Fig.5, where the the probability distribution of the order parameter $\phi(\mathbf{x},t)$ is plotted. It should be noted that, with a considerably larger($\alpha=0.22$) time step than that of the conventional CD2, current scheme shows much better agreement with the Gaussian. We also compare the spectra of the normalized energy distribution(equilibrium ratio) in Fig.6 for the CD2 and present scheme with the analytical one, from which the break down of FDT at discrete level for CD2 is quite apparent.

To highlight the quantitative improvements due to the present work, we present the polar plots of the normalized energy distribution(equilibrium ratio) at different wave numbers in Fig.7. We remind that for model B with non-interacting single well Landau-Ginzburg Hamiltonian of the order parameter, energy in any wave number is $0.5k_{B}T$ irrespective of the wave number, thus the normalized energy at any wave number should be one. As expected, result from current scheme is quite close to the analytical result, while central discretization shows violation of FDT even with the use of a time step which is four times smaller than that of the present scheme.

8 Model B: Anharmonic fluctuations

In this section, we extend our treatment to the case of a inhomogeneous system, where different regions are coupled via free-energy gradients of the order parameter entering in the Landau-Ginzburg Hamiltonian ${\cal F}$. Use of the aforementioned isotropic operators and delayed discretization, the discrete model B assumes the following form, with $f(\phi_{i}^{n})$ denoting the value of the respective discrete free energy density at lattice cite $i$ and time step $n$.

[TABLE]

To establish the consistency of the present work over the usual central difference type operators, we compare the spectra of the structure factors obtained from three different discrete formulation of model B, viz. isotropic, central difference and Fourier pseudo spectral in Fig.8. Here free energy densities for both single phase $f(\phi)=\frac{1}{2}A\phi^{2}$ and two phase equilibrium $f(\phi)=\frac{1}{2}A\phi^{2}+\frac{1}{4}B\phi^{4}$ are considered. The pseudo spectral, ensuring exact space derivatives, preserve the FDT at discrete level but the computation is much more expensive than the other two counterparts viz. CD2 and isotropic. On the other hand it is evident here that the failure of preserving FDT at the lattice level leads to an energy loss at the higher wave numbers for the case of CD2 as compared to that of the pseudo spectral result. The isotropic formulation does not show any such energy loss at the higher wave numbers, instead its energy spectra is very close to that of the pseudo spectral. To illustrate the quantitative aspect of the present formulation over the traditional schemes we also present the polar spectra in Fig.9 and Fig.10 for the structure factors of these three cases, which clearly bring out the anisotropy and breakdown of FDT at the discrete level for the traditional CD2 schemes as opposed to the isotropic discretization.

It should be stressed here that in all these cases the isotropic scheme operates at a time step four times larger than that of the CD2 scheme.

To characterize the diffusive behaviour of the central difference scheme over the present formulation, we investigate the dynamics of the order parameter $\phi$ through its instantaneous distribution. In Fig.11 three consecutive instantaneous states of the order parameter evolution are compared for three different cases, namely CD2, isotropic and pseudo-spectral.

It is evident from this plot that, for CD2, the system is in a little more quenched state than that of the isotropic or pseudo-spectral methods. This is indicative of the fact that the CD2 has more diffusive behaviour at higher wave numbers, which is also observed in Fig.8.

9 Outlook

To conclude, we have presented a discrete framework where the essence of phase separation dynamics in terms of fluctuation-dissipation relation is preserved. Thus, similar to the cell dynamical method, a fully self-consistent framework at discrete level is obtained. The present approach also allows the discrete framework to inherit transport properties and free energies known from PDE based formulation.

It should also be pointed out that, the present formulation can easily be extended to three dimensions(3D), as all of these aforementioned isotropic operators can also be derived in 3D. For a detailed formulation of these operators in higher dimensions one might refer to [32, 29]. Also in order to contrast the present formulation with cell dynamics, we write Eq.38 in a form analogous to Eq.(2.7,2.8) of [35] as:

[TABLE]

where,

[TABLE]

[TABLE]

Here, $N_{i},N_{j},N_{k},N_{l}$ are the nearest and next nearest neighbours and so on. It can be observed that the discrete form of the model B in Eq.42 and the discrete operators Eq.43 and Eq.44 preserve the same structure of a Cell Dynamical System(CDS)[35]. Though the form is similar, two key differences, wider stencil and use of past data are must be noted. These two differences allow us to formulate a cell dynamical system where FDT is preserved even at the discrete level and the connection with the PDE is also apparent.

Reference