Deriving phase field crystal theory from dynamical density functional theory: consequences of the approximations

TL;DR

This paper critically examines the derivation of phase field crystal (PFC) theory from dynamical density functional theory (DDFT), highlighting the impact of common approximations on the accuracy of phase diagrams and crystal stability predictions.

Contribution

The authors identify key approximations in deriving PFC from DDFT that lead to significant inaccuracies, and suggest that deriving PFC models for the logarithm of the density may improve accuracy.

Findings

Neglecting certain terms affects crystal stability predictions.

Standard polynomial approximations introduce spurious phase behavior.

A one-mode approximation for log density yields surprisingly accurate results.

Abstract

Phase field crystal (PFC) theory, extensively used for modelling the structure of solids, can be derived from dynamical density functional theory (DDFT) via a sequence of approximations. Standard derivations neglect a term of form $\nabla\cdot[n\nabla L n]$, where $n$ is the scaled density profile and $L$ is a linear operator. We show that this term makes a significant contribution to the stability of the crystal, and dropping this term from the theory forces another approximation, that of replacing the logarithmic term from the ideal gas contribution to the free energy with its truncated Taylor expansion, to yield a polynomial in $n$. However, the consequences of doing this are the presence of an additional spinodal in the phase diagram, so the liquid is predicted first to freeze and then to melt again as the density is increased; and other periodic structures are erroneously predicted to be thermodynamic equilibria. A second approximation is to replace $L$ by a gradient expansion. This leads to the possibility of solutions failing to exist above a certain value of the average density. We illustrate these conclusions with a simple model two-dimensional fluid. The consequences of the PFC approximations are that the phase diagram is both qualitatively incorrect, in that it has a stripe phase, and quantitatively incorrect (by orders of magnitude) regarding the properties of the crystal. Thus, although PFC models are successful as phenomenological models of crystallisation, we find it impossible to derive the PFC model as an accurate approximation to DDFT, without introducing spurious artefacts. However, making a simple one-mode approximation for the logarithm of the density distribution is surprisingly accurate, which gives a tantalising hint that accurate PFC-type theories may instead be derived for the field $\log(\rho(x))$, rather than for the density profile itself.

Tables

Table 1. Table 1: Various versions of DDFT and PFC, in order of appearance, along with references to the equations defining the dynamics, the free energy and the chemical potential. We also give the quadratic ( Q 𝑄 Q ), cubic ( C 𝐶 C ) and quartic ( R 𝑅 R ) coefficients. PFC- α 𝛼 \alpha  is the phenomenological model, also known as the conserved Swift–Hohenberg equation; PFC- ϵ italic-ϵ \epsilon  is equivalent to  PFC- α 𝛼 \alpha .

Name	Truncate at $𝒪\left(c^{\left(4\right)}\right)$	LDA (31) for $c^{\left(3\right)}$ and $c^{\left(4\right)}$	RY/RPA: $c^{\left(3\right)}=c^{\left(4\right)}=0$	Gradient expansion of $ℒ$ (44)	Constant mobility, expand $\log \left(1+n\right)$	Dynamics	Free energy	Chemical potential	$Q$, $C$, $R$
PFC-$\alpha$	Yes	N/A	Yes	Yes	Yes	(2)	(1)	—	$Q=0$ $C=-1$
DDFT-0						(5)	(10)	(7)	—
DDFT-1	Yes					(24), (26)	(25)	(28), (30)	—
DDFT-2	Yes	Yes				(34)	(32)	(36)	(35)
DDFT-3	Yes	N/A	Yes			(41)	(40)	(42)	$Q=\frac{1}{2}$ $C=0$ $R=0$
DDFT-4	Yes	N/A		Yes		(46)	as (32)	as (36)	(35)
DDFT-5	Yes	N/A	Yes	Yes		(48)	(47)	(49)	$Q=\frac{1}{2}$ $C=0$ $R=0$
PFC-$\beta$	Yes	Yes			Yes	(53), (54)	(52)	(55)	(56)
PFC-$\gamma$	Yes	N/A	Yes		Yes	(58)	(57)	(59)	$Q=\frac{1}{2}$ $C=-\frac{1}{3}$
PFC-$\delta$	Yes	N/A		Yes	Yes	as (54)	as (52)	as (55)	(56)
PFC-$ϵ$	Yes	N/A	Yes	Yes	Yes	(60)	as (57)	as (59)	$Q=\frac{1}{2}$ $C=-\frac{1}{3}$

Table 2. Table 2: Linear theory for the GEM-4 potential in one, two and three dimensions. Solving d ​ σ d ​ k ​ ( 1 ) = 0 d 𝜎 d 𝑘 1 0 \frac{{\rm d}\sigma}{{\rm d}{k}}(1)=0 and σ ​ ( 1 ) = 0 𝜎 1 0 \sigma(1)=0 gives R 𝑅 R and  ρ 0 ​ β ​ ϵ subscript 𝜌 0 𝛽 italic-ϵ \rho_{0}\beta\epsilon , while γ 𝛾 \gamma and σ ​ ( 0 ) 𝜎 0 \sigma(0) are computed from d 2 ​ σ d ​ k 2 ​ ( 1 ) superscript d 2 𝜎 d superscript 𝑘 2 1 \frac{{\rm d}^{2}\sigma}{{\rm d}{k^{2}}}(1) and ( 85 ) respectively. We also give the values of E B subscript 𝐸 𝐵 E_{B} for the EOF in Eq. ( 76 ).

Dimension	$R$	${\rho }_0\beta ϵ$	$\gamma$	$\sigma \left(0\right)$	$E_B$
1	4.5918	1.1629	3.2051	$-10.680$	$7.475$
2	5.0962	0.2455	4.3692	$-18.752$	$14.383$
3	5.5719	0.0455	5.6889	$-31.305$	$25.616$

Equations

$${\mathcal{F}}_{\text{α}}[n]=\int \left(\frac{1}{2}n\left((k_s^2+{\nabla }^2{)}^2-r\right)n+\frac{1}{4}{n}^4\right)\mathrm{d}\mathbf{x},$$

$$\frac{\partial n}{\partial t}={\nabla }^2\left(\frac{\delta {\mathcal{F}}_{\text{α}}}{\delta n}\right)=-{\nabla }^2\left(rn-(k_s^2+{\nabla }^2{)}^{2}n-n^3\right),$$

$$\frac{\partial n}{\partial t}=-\frac{\delta {\mathcal{F}}_{\text{α}}}{\delta n}=rn-(k_s^2+{\nabla }^2{)}^{2}n-n^3,$$

$$\frac{\partial \rho }{\partial t}=-\nabla \cdot \mathbf{j},$$

$$\frac{\partial \rho }{\partial t}=\nabla \cdot \left[\beta M(\rho )\nabla \frac{\delta \mathcal{F}}{\delta \rho }\right],$$

$$\frac{\mathrm{d}\mathcal{F}}{\mathrm{d}t}=-\int \beta M(\rho ){\left|\nabla \frac{\delta \mathcal{F}}{\delta \rho }\right|}^2\mathrm{d}\mathbf{x}\le 0,$$

$$\frac{\delta \mathcal{F}}{\delta \rho }=\text{constant}=\mu ,$$

$$\Omega [\rho ]=\mathcal{F}[\rho ]-\mu N=\mathcal{F}[\rho ]-\mu \int \rho (\mathbf{x})\mathrm{d}\mathbf{x}.$$

$$\frac{\delta \Omega }{\delta \rho }=\frac{\delta \mathcal{F}}{\delta \rho }-\mu ,$$

$$\mathcal{F}[\rho ]=k_{B}T\int \rho \left(\log ({\Lambda }^d\rho )-1\right)\mathrm{d}\mathbf{x}+{\mathcal{F}}_{\text{ex}}[\rho ]+\int \rho U_{\text{ext}}\mathrm{d}\mathbf{x},$$

$$\beta \frac{\delta \mathcal{F}}{\delta \rho }=\log \rho +\beta \frac{\delta {\mathcal{F}}_{\text{ex}}}{\delta \rho },$$

$$\beta \nabla \frac{\delta \mathcal{F}}{\delta \rho }=\frac{1}{\rho }\nabla \rho +\dots ,$$

$$\frac{\partial \rho }{\partial t}={\nabla }^2\rho +\dots ,$$

$$\begin{array}{rl}{\mathcal{F}}_{\text{ex}}[\rho ]& ={\mathcal{F}}_{\text{ex}}[{\rho }_0]-k_{B}T\int c^{(1)}({\mathbf{x}}_1)\Delta \rho ({\mathbf{x}}_1)\mathrm{d}{\mathbf{x}}_1\\ & -\frac{k_{B}T}{2!}\int c^{(2)}({\mathbf{x}}_1,{\mathbf{x}}_2)\Delta \rho ({\mathbf{x}}_1)\Delta \rho ({\mathbf{x}}_2)\mathrm{d}{\mathbf{x}}_1\mathrm{d}{\mathbf{x}}_2\\ & -\frac{k_{B}T}{3!}\int c^{(3)}({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3)\times \\ & \Delta \rho ({\mathbf{x}}_1)\Delta \rho ({\mathbf{x}}_2)\Delta \rho ({\mathbf{x}}_3)\mathrm{d}{\mathbf{x}}_1\mathrm{d}{\mathbf{x}}_2\mathrm{d}{\mathbf{x}}_3\\ & +\text{similar fourth order term}+\dots .\end{array}$$

$$c^{(n)}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n)=-\beta \frac{{\delta }^n{\mathcal{F}}_{\text{ex}}}{\delta \rho ({\mathbf{x}}_1)\dots \delta \rho ({\mathbf{x}}_n)}[{\rho }_0].$$

$$c^{(n)}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n)=c^{(n)}({\Delta \mathbf{x}}_2,{\Delta \mathbf{x}}_3\dots ,{\Delta \mathbf{x}}_n),$$

$$\rho (\mathbf{x},t)={\rho }_0(1+n(\mathbf{x},t)).$$

$$\frac{\partial n}{\partial t}={\nabla }^{2}n-{\nabla }^2{c}^{(1)}-\nabla \cdot \left[n\nabla c^{(1)}\right]+\dots$$

$$\begin{split}\frac{\partial n}{\partial t}=&\nabla^{2}n-\rho_{0}\nabla^{2}\!\!\int\!c^{(2)}({\bm{x}},{\bm{x}}_{2})n({\bm{x}}_{2}){\rm d}{\bm{x}}_{2}\\ &{}-\rho_{0}\nabla\cdot\left[n\nabla\!\!\int\!c^{(2)}({\bm{x}},{\bm{x}}_{2})n({\bm{x}}_{2}){\rm d}{\bm{x}}_{2}\right]\\ &{}-\frac{\rho_{0}^{2}}{2}\nabla^{2}\!\!\int\!c^{(3)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3})n({\bm{x}}_{2})n({\bm{x}}_{3}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}\\ &{}-\frac{\rho_{0}^{2}}{2}\nabla\cdot\left[n\nabla\!\!\int\!c^{(3)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3})n({\bm{x}}_{2})n({\bm{x}}_{3}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}\right]\\ &{}-\frac{\rho_{0}^{3}}{6}\nabla^{2}\!\!\int\!c^{(4)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3},{\bm{x}}_{4})\times{}\\ &\qquad\qquad\qquad n({\bm{x}}_{2})n({\bm{x}}_{3})n({\bm{x}}_{4}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}{\rm d}{\bm{x}}_{4}\\ &{}-\frac{\rho_{0}^{3}}{6}\nabla\cdot\bigg{[}n\nabla\!\!\int\!c^{(4)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3},{\bm{x}}_{4})\times{}\\ &\qquad\qquad\qquad n({\bm{x}}_{2})n({\bm{x}}_{3})n({\bm{x}}_{4}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}{\rm d}{\bm{x}}_{4}\bigg{]}+\dots\\ \end{split}$$

$$\frac{\partial n}{\partial t}=-{\nabla }^2\mathcal{L}n,$$

$$\mathcal{L}n(\mathbf{x})=-n(\mathbf{x})+{\rho }_0\int c^{(2)}(\mathbf{x},{\mathbf{x}}_2)n({\mathbf{x}}_2)\mathrm{d}{\mathbf{x}}_2.$$

$$\mathcal{L}{e}^{i\mathbf{k}\cdot \mathbf{x}}=\sigma (\mathbf{k})e^{i\mathbf{k}\cdot \mathbf{x}},$$

$$\begin{array}{rl}\sigma (k)& =-1+{\rho }_0\int c^{(2)}(|{\mathbf{x}}_2-\mathbf{x}|)e^{i\mathbf{k}\cdot ({\mathbf{x}}_2-\mathbf{x})}\mathrm{d}{\mathbf{x}}_2\\ & =-1+{\rho }_0{\hat{c}}^{(2)}(k),\end{array}$$

$$\begin{split}\frac{\partial n}{\partial t}=&-\nabla^{2}\left({\cal L}n+\tfrac{1}{2}n^{2}\right)-\nabla\cdot\left[n\nabla{\cal L}n\right]\\ &{}-\frac{\rho_{0}^{2}}{2}\nabla^{2}\!\!\int\!c^{(3)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3})n({\bm{x}}_{2})n({\bm{x}}_{3}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}\\ &{}-\frac{\rho_{0}^{2}}{2}\nabla\cdot\left[n\nabla\!\!\int\!c^{(3)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3})n({\bm{x}}_{2})n({\bm{x}}_{3}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}\right]\\ &{}-\frac{\rho_{0}^{3}}{6}\nabla^{2}\!\!\int\!c^{(4)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3},{\bm{x}}_{4})\times{}\\ &\qquad\qquad\qquad n({\bm{x}}_{2})n({\bm{x}}_{3})n({\bm{x}}_{4}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}{\rm d}{\bm{x}}_{4}\\ &{}-\frac{\rho_{0}^{3}}{6}\nabla\cdot\bigg{[}n\nabla\!\!\int\!c^{(4)}({\bm{x}},{\bm{x}}_{2},{\bm{x}}_{3},{\bm{x}}_{4})\times{}\\ &\qquad\qquad\qquad n({\bm{x}}_{2})n({\bm{x}}_{3})n({\bm{x}}_{4}){\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}{\rm d}{\bm{x}}_{4}\bigg{]}\\ \ \end{split}$$

$$\begin{split}\beta{{\cal F}_{1}}[n]&=\int\!\big{(}[1+n({\bm{x}}_{1})]\log[1+n({\bm{x}}_{1})]-n({\bm{x}}_{1})\big{)}{\rm d}{\bm{x}}_{1}\\ &\quad{}-\frac{1}{2}\!\int\!\left(n^{2}({\bm{x}}_{1})+n({\bm{x}}_{1}){\cal L}n({\bm{x}}_{1})\right){\rm d}{\bm{x}}_{1}\\ &\quad{}-\frac{\rho_{0}^{2}}{6}\!\!\int\!c^{(3)}({\bm{x}}_{1},{\bm{x}}_{2},{\bm{x}}_{3})\times{}\\ &\qquad\qquad\quad n({\bm{x}}_{1})n({\bm{x}}_{2})n({\bm{x}}_{3}){\rm d}{\bm{x}}_{1}{\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}\\ &\quad{}-\frac{\rho_{0}^{3}}{24}\!\int\!c^{(4)}({\bm{x}}_{1},{\bm{x}}_{2},{\bm{x}}_{3},{\bm{x}}_{4})\times{}\\ &\qquad\qquad\quad n({\bm{x}}_{1})n({\bm{x}}_{2})n({\bm{x}}_{3})n({\bm{x}}_{4}){\rm d}{\bm{x}}_{1}{\rm d}{\bm{x}}_{2}{\rm d}{\bm{x}}_{3}{\rm d}{\bm{x}}_{4},\end{split}$$

$$\frac{\partial n}{\partial t}=\nabla \cdot \left[\beta (1+n)\nabla \frac{\delta {\mathcal{F}}_1}{\delta n}\right].$$

$$\frac{\delta \mathcal{F}}{\delta \rho }=\frac{\delta {\mathcal{F}}_1}{\delta n}.$$

$$\frac{\delta {\mathcal{F}}_1}{\delta n}=\Delta \mu ,$$

$${\mu }_0=k_{B}T\log {\Lambda }^d{\rho }_0-k_{B}T{c}^{(1)}[{\rho }_0],$$

$$\begin{array}{rl}\beta \frac{\delta {\mathcal{F}}_1}{\delta n}& =\log \left(1+n(\mathbf{x})\right)-n(\mathbf{x})-\mathcal{L}n(\mathbf{x})\\ & -\frac{{\rho }_0^2}{2}\int c^{(3)}(\mathbf{x},{\mathbf{x}}_2,{\mathbf{x}}_3)n({\mathbf{x}}_2)n({\mathbf{x}}_3)\mathrm{d}{\mathbf{x}}_2\mathrm{d}{\mathbf{x}}_3\\ & -\frac{{\rho }_0^3}{6}\int c^{(4)}(\mathbf{x},{\mathbf{x}}_2,{\mathbf{x}}_3,{\mathbf{x}}_4)\times \\ & n({\mathbf{x}}_2)n({\mathbf{x}}_3)n({\mathbf{x}}_4)\mathrm{d}{\mathbf{x}}_2\mathrm{d}{\mathbf{x}}_3\mathrm{d}{\mathbf{x}}_4.\end{array}$$

Full text

I Introduction

The phase field crystal (PFC) theory for matter is widely used and has been successfully applied to describe a broad range of phenomena, including, for example, grain boundary dynamics Elder et al. (2002); Elder and Grant (2004), crystal nucleation Backofen and Voigt (2010); Toth et al. (2010), crystal growth Archer et al. (2012), glass formation Berry et al. (2008), crack propagation Elder and Grant (2004) and many other properties of condensed matter. For more background and examples of situations to which the PFC theory has been applied, see the excellent review Emmerich et al. (2012). The PFC theory was originally proposed, in the spirit of ‘regular’ phase field theory (PFT), as a diffuse-interface theory for the time evolution of an order parameter field Elder et al. (2002). The equations of PFT are obtained via symmetry, thermodynamic and other arguments and the result is a theory that is widely used in materials science and other disciplines to model the structure of materials. For more background on PFT see for example Ref. Boettinger et al. (2002) and references therein.

The central and original idea in extending PFT to arrive at PFC theory is to incorporate aspects of the microscopic structure of the material into the model Elder et al. (2002). The result is a theory that operates on atomic length scales and diffusive time scales Emmerich et al. (2012). By this we mean that PFC theory is a theory for a field that exhibits numerous maxima, each of which is identified as the average location of the atoms (or more generally ‘particles’) in the system. This idea is powerful because, by building into the theory more information about the underlying material structure, it enables the inclusion of much more of the physics coming from particle correlations to be incorporated. Over the years several variants of PFC theory have been developed that are able to describe a range different crystalline (and even quasicrystalline) structures Jaatinen et al. (2009); Pisutha-Arnond et al. (2013); Wu et al. (2010); Barkan et al. (2011); Achim et al. (2014); Subramanian et al. (2016); Jiang et al. (2017); Savitz et al. (2018).

Thus, the original PFC Elder et al. (2002) may be viewed as the simplest partial differential equation model one can conceive of for a conserved order parameter exhibiting peaks arranged with crystalline ordering. It is obtained from a (scaled) free energy {\cal F_{\text{\alpha}}} that is a functional of the dimensionless order parameter $n$:

[TABLE]

where $n({\bm{x}},t)$ is a field that depends on position in space ${\bm{x}}$ and on time $t$, and $k_{s}$ is an inverse length scale that determines the lattice spacing of the crystal. The parameter $r$ defines how near the system is to freezing. The time evolution of the conserved field $n$ is given by the dynamics

[TABLE]

where \frac{{\delta{\mkern-1.0mu}{\cal F_{\text{\alpha}}}}}{{\delta{\mkern-1.0mu}n}} is the functional derivative of {\cal F_{\text{\alpha}}} with respect to $n({\bm{x}})$.

Given the ingredients in the model, it is therefore no surprise that PFC theory is closely related to the Swift–Hohenberg equation Swift and Hohenberg (1977):

[TABLE]

which is one of the archetypal equations in pattern formation theory. As one can see above, both the Swift–Hohenberg equation and PFC theory can be expressed as a different type of dynamics based on the same free energy functional Emmerich et al. (2012). The Swift–Hohenberg equation (3) is based on an underlying dynamics that seeks to minimise the free energy over time, whilst the PFC dynamics (2), which also decreases the free energy over time, in addition enforces a conservation of the average value of the order parameter in the system. Thus, the PFC equation (2) is sometimes referred to as the conserved Swift–Hohenberg equation Thiele et al. (2013); Sagui and Desai (1994); Knobloch (2015); Matthews and Cox (2000); Emmerich et al. (2012).

In the years after PFC theory was originally proposed it was realised that it could be derived from classical dynamical density functional theory (DDFT) Elder et al. (2007); van Teeffelen et al. (2009); Huang et al. (2010); Emmerich et al. (2012); Archer et al. (2012), via a sequence of several different approximations. Below, we say much more on what these approximations are. DDFT is a theory for the time evolution of the ensemble average one-body (number) density profile $\rho({\bm{x}},t)$, for a non-equilibrium system of interacting classical particles. DDFT is based on equilibrium density functional theory (DFT) Evans (1979, 1992); Hansen and McDonald (1992) and for an equilibrium system, DDFT is equivalent to DFT. DDFT was originally developed as a theory for Brownian particles with over-damped stochastic equations of motion Marconi and Tarazona (1999, 2000); Archer and Rauscher (2004); Archer and Evans (2004), but it has also been extended to describe the dynamics of under-damped systems and atomic or molecular systems where the particle dynamics is governed by Newton’s equations of motion Archer (2006, 2009); Goddard et al. (2012, 2013); Durán-Olivencia et al. (2017); Schmidt (2018). This body of work shows that when such systems are not too far from equilibrium, then the dynamics predicted by the original DDFT is still often correct in the long-time limit where the particle dynamics is dominated by diffusive processes. This is because DDFT corresponds to a dynamics given by the continuity equation

[TABLE]

where the current $\mathbf{j}\propto-\nabla\mu({\bm{x}},t)$, with $\mu({\bm{x}},t)$ a local (non-equilibrium) chemical potential Marconi and Tarazona (1999, 2000); Archer and Rauscher (2004); Archer and Evans (2004). Eq. (4) is of course expected since the total number of particles in the system $N=\int\rho({\bm{x}},t){\rm d}{\bm{x}}$ is a conserved quantity.

Refs. Elder et al. (2007); van Teeffelen et al. (2009); Huang et al. (2010); Emmerich et al. (2012); Archer et al. (2012) give various different derivations of the PFC model, starting from DFT and/or DDFT. Here, starting from DDFT, we systematically show how all the various different theories are related and we identify and highlight the significance of each of the approximations that are made in the derivation of PFC theory. We show that there is a particular term of the form $\nabla\cdot[n\nabla{\cal L}n]$, where ${\cal L}$ is a linear operator, that is almost universally neglected because it is ‘of higher order’ Huang et al. (2010), but this term is actually important for stabilizing crystalline structures: its contribution is of the same order as some of the terms that are retained. As we explain in detail, neglecting this term essentially forces one to make the Taylor expansion of the ideal gas logarithmic term in the free energy in order to recover something physically reasonable. We show that neglecting this term, as is done in PFC theory, and the subsequent replacement of the logarithm by its Taylor series, leads to the spurious appearance in the phase diagram of an extra spinodal and also alters the relative stabilities of the crystal state compared to a stripe phase and also other phases, leading in two dimensions (2D) to the stripe phase becoming the global free energy minimum state for certain parameter values. Essentially, all this behaviour originates because the function $\log(1+n)$ has one root, but when it is replaced by a truncated Taylor expansion, the resulting polynomial generally has two roots. Our arguments also directly apply in three dimensions to explain why lamellar phases occur as equilibrium phases in PFC theory. Recall that most PFC theories predict that as one moves in the phase diagram away from the region where there is coexistence between the liquid and the crystal, moving deeper into the crystalline portion of the phase diagram, such stripe/lamellar phases appear as equilibrium structures and are global minima of the free energy Emmerich et al. (2012). Of course, particles with isotropic pair interactions generally never ‘freeze’ to form striped phases, unless they have an unusual and special form for the pair potential between the particles Imperio and Reatto (2004, 2006); Archer and Wilding (2007). DDFT, taken together with a reliable approximation for the Helmholtz free energy functional of course does not predict such stripe phases for crystallisation from simple liquids.

The linear operator ${\cal L}$ has the form of a non-local convolution involving the pair direct correlation function plus another simpler term (see Eq. (20) below). Another approximation that is often made in deriving PFC theories is to approximate ${\cal L}$ by a gradient expansion involving derivatives. We show below that if one makes this approximation whilst simultaneously retaining the logarithmic term from the ideal gas free energy, this results in a theory that still predicts reasonably accurately the freezing transition, but as one increases the average density, moving deeper into the region of the phase diagram where the crystal phase occurs, there is a point where $\rho({\bm{x}})\to 0$ at the points in space ${\bm{x}}$ between the density peaks, where the density is a minimum. On increasing the average density beyond this point in the phase diagram, there is no solution to the theory. We analyse in detail this singular behaviour. As $\rho({\bm{x}})\to 0$ we have $\log\rho({\bm{x}})\to-\infty$, of course. In the equation for the equilibrium density profile this divergence is initially balanced by the term involving the fourth derivative, $\partial^{4}\rho/\partial x^{4}\to+\infty$. However, when the average density in the system is increased beyond the value at which this divergence occurs, we find there is no solution.

We illustrate these conclusions by finding the predicted structures and phase diagram for the 2D version of the GEM-4 (Generalised Exponential Model of index 4) Mladek et al. (2006); Prestipino and Saija (2014), chosen because DDFT based on a simple approximation (the so-called random phase approximation (RPA) Likos (2001)) for the Helmholtz free energy functional can be very accurate for predicting the equilibrium structures formed in this model and also the thermodynamics Mladek et al. (2007); Prestipino and Saija (2014); Archer et al. (2014). At higher temperatures, the 2D GEM-4 system exhibits just a single fluid phase and at higher densities a single crystal phase. At lower temperatures, where the RPA DDFT is no longer accurate, there is a hexatic phase and multiple crystalline phases as the density is increased Prestipino and Saija (2014). Here we do not consider this regime, restricting ourselves to the regime where there is just one fluid and one crystal phase, which are predicted accurately by the RPA DDFT. This enables us to investigate the effect of making subsequent approximations to the DDFT, including those made to derive PFC theory. We find that the PFC type theories spuriously predict three additional phases that are in reality not present in the phase diagram (i.e., are not thermodynamically stable). These are (i) a stripe phase, (ii) what we refer to as ‘down hexagons’ (in contrast to the true crystal structure, which we refer to as ‘up hexagons’) and then at even higher densities a melting to form (iii) another uniform liquid phase. We show how the approximations made in deriving the PFC result in these structures being predicted.

The final contribution of this paper is to show that there is a very simple and accurate ansatz one can make for the form of the equilibrium crystal density profile in DDFT (and so also for DFT, of course). The ansatz is $\rho({\bm{x}})=\rho_{0}e^{\phi({\bm{x}})}$, where $\rho_{0}$ is a constant and the field $\phi({\bm{x}})$ is approximated by a sinusoid of the form $\phi({\bm{x}})\approx\phi_{0}+\phi_{1}e^{i{\bm{k}}\cdot{\bm{x}}}+$complex conjugate (in one dimension), plus other similar terms (in higher dimensions), where $\phi_{0}$ and $\phi_{1}$ are constants. The results presented here are for the GEM-4 model and show why this approximation is unexpectedly accurate: the approximation is able to replicate almost exactly the numerical solution to the DDFT problem, from small to arbitrarily large amplitude density variations. We expect this ansatz also to be reliable for other systems. This form of one-mode theory gives a hint for future directions to develop accurate PFC-type theories, since using a one-mode approximation in PFC is often fairly accurate.

This paper is structured as follows: In Sec. II we present our systematic step-by-step derivation of PFC, starting from DDFT. After each approximation, we carefully state the model, i.e., we give the corresponding free energy functional and also the expression for the chemical potential, which is a quantity that is a constant at all points in space for equilibrium states. In order to keep track of the different orders in which the approximations can be made, we give each model a name, starting with PFC-$\alpha$ for the original PFC model in Eq. (2) above, and with DDFT-0 for the original formulation of DDFT below. The different DDFT approximations result in five different versions, DDFT-1 to DDFT-5. Similarly, we explain the various different approximations that can be made to each of these, leading to a corresponding PFC theory, which we name PFC-$\alpha$ to PFC-$\epsilon$. Note that the criterion we use here for distinguishing between whether we refer to a theory as a DDFT or a PFC is based on whether the free energy which is minimised by the dynamical equations (i.e., the Lyapunov functional) has the logarithmic ideal gas term or not: if it does not have the logarithm, we refer to it as a PFC. Table 1 below is there to help the reader navigate the various models and the approximations made in each one. Sec. II concludes with a summarising discussion. In Sec. III we present results for the GEM-4 system comparing predictions for the density profiles and thermodynamics of equilibria for two of the different DDFT theories and also two of the PFC theories. In this section we also present the phase diagrams for the GEM-4 system predicted by these different DDFT and PFC theories. By comparing all of these we are able to assess the accuracy of the different theories and the validity of the various approximations. In Sec. IV we discuss the implications of the main two approximations and analyse the singular behaviour displayed by some models. In Sec. V we introduce the ansatz $\rho({\bm{x}})=\rho_{0}e^{\phi({\bm{x}})}$ and derive the new one-mode approximation for DDFT. We draw our conclusions in Sec. VI. The paper includes two appendices in which we describe the numerical (continuation) methods we use to calculate the density profiles.

II Derivation of the Phase Field Crystal model from DDFT

In this section we progress from the original formulation of DDFT (which we call DDFT-0) through a series of approximations (DDFT-1, …, DDFT-5), as listed in Table 1. Our main starting point is DDFT-2. From this point, there are three main approximations that can be made (or not made): (i) the Ramakrishan–Yussouff (RY) or the random phase approximation (RPA) for the free energy; (ii) the gradient expansion of the convolution term; and (iii) the Taylor expansion of the logarithmic term. Making (or not making) the first two of these approximations results in DDFT-3, DDFT-4 and DDFT-5. Then, making the third approximation from DDFT-2 results in PFC-$\beta$, from DDFT-3 results in PFC-$\gamma$, and so on up to PFC-$\epsilon$. The PFC-$\epsilon$ model can be rescaled to recover the original version of PFC, PFC-$\alpha$, see Eqs. (1) and (2). The various models are summarised in Table 1. Amongst the models we present below, DDFT-5 is equivalent to the model derived by Huang et al. Huang et al. (2010) and advocated by van Teeffelen et al. van Teeffelen et al. (2009) (named PFC1 in that paper), and DDFT-3 and PFC-$\epsilon$ are equivalent to the models named DDFT and PFC2 by van Teeffelen et al. van Teeffelen et al. (2009).

II.1 Dynamic Density Functional Theory: DDFT-0

The starting point for all of our derivations is the key DDFT equation Marconi and Tarazona (1999, 2000); Archer and Rauscher (2004); Archer and Evans (2004):

[TABLE]

where $\beta=(k_{B}T)^{-1}$ (with $k_{B}$ being Boltzmann’s constant and $T$ being temperature), $M(\rho)$ is the positive $\rho$-dependent mobility. The Helmholtz free energy ${\cal F}[\rho]$ depends on the density profile $\rho({\bm{x}},t)$ integrated over space; hence ${\cal F}[\rho]$ depends on time but not on position Archer and Rauscher (2004). The expression ${{\delta{\mkern-1.0mu}{\cal F}}}/{{\delta{\mkern-2.0mu}\rho}}$ is the functional derivative of ${\cal F}$ with respect to $\rho({\bm{x}},t)$, which therefore depends on both time and on position. DDFT usually takes $M(\rho)=D\rho$, i.e., the mobility is proportional to density Marconi and Tarazona (1999, 2000); Archer and Rauscher (2004); Archer and Evans (2004), where $D$ is the diffusion coefficient. We henceforth scale time so that $D=1$. With boundary conditions that do not allow material to enter or leave the system, $N=\int\!\rho({\bm{x}}){\rm d}{\bm{x}}$ (or equivalently, the mean density) is a constant of the motion and is the total number of particles in the system.

With suitable boundary conditions, one can readily show that the Helmholtz free energy decreases monotonically with time:

[TABLE]

so (assuming that ${\cal F}[\rho]$ is bounded below) the system typically evolves to a (local) minimum of ${\cal F}$, which is a dynamically stable equilibrium of (5). Here, ‘dynamically stable’ means that small perturbations away from the equilibrium decay, and ‘equilibrium’ means that ${\partial\rho}/{\partial t}=0$ and ${{\rm d}{\cal F}}/{{\rm d}{t}}=0$. Owing to the dynamics being governed by a continuity equation (4), such perturbations cannot change the mean density. Local minima of ${\cal F}$ that are not the global minimum are thermodynamically metastable. The system can also have dynamically unstable equilibria, for which ${\cal F}$ is a saddle or maximum. From (6), we see that all equilibria of (5) satisfy $\nabla({\delta{\mkern-1.0mu}{\cal F}}/{\delta{\mkern-2.0mu}\rho})=0$, so

[TABLE]

where $\mu$ is the chemical potential of the equilibrium. This is of course the Euler–Lagrange equation for the problem of finding stationary points of the functional ${\cal F}[\rho]$, subject to the constraint of fixed mean density. Note however that when evolving (5) forward in time from an arbitrary initial condition, $\mu$ is not necessarily known a priori.

The theory can also be cast in terms of the grand potential (also called the Landau free energy) functional Evans (1979, 1992); Hansen and McDonald (1992):

[TABLE]

From this it follows that the functional derivative of $\Omega$ is

[TABLE]

and that this is zero at equilibrium: equilibria are extreme values of $\Omega$. Like the Helmholtz free energy, the grand potential decreases monotonically with time, since Eq. (6) is also true if one replaces ${\cal F}$ by $\Omega$. Therefore, for two phases to coexist, they must have the same specific grand potential (i.e., the same pressure) and the same chemical potential. Thus, the global minimum of $\Omega$ for a given $\mu$ and $T$ is the thermodynamic equilibrium state of the system Evans (1979, 1992); Hansen and McDonald (1992).

Following the usual approach in DFT, we separate the Helmholtz free energy into three parts: the ‘ideal gas’ contribution, which is proportional to the temperature but takes no account of particle interactions, an excess (${\cal F_{\text{ex}}}$) over the ideal gas contribution arising from the particle interactions, and the contribution due to an external potential $U_{\text{ext}}({\bm{x}})$, as follows Evans (1979, 1992); Hansen and McDonald (1992):

[TABLE]

where the integral is taken over the volume $V$ in three dimensions ($d=3$) (or the area in 2D, $d=2$) and where $\Lambda$ is the thermal de Broglie wavelength. Since for our purposes here the value of $\Lambda$ is irrelevant (changing $\Lambda$ will shift the values of ${\cal F}$ and $\mu$ by constants), we henceforth set $\Lambda=1$. We also consider bulk systems and so we assume that $U_{\text{ext}}=0$. With the separation in Eq. (10), we have

[TABLE]

which gives

[TABLE]

where on the right hand side we only explicitly write the contribution from the ideal gas part of the free energy. Inserting this into Eq. (5) with $M=D\rho$ we obtain

[TABLE]

in which the coefficient in front of the term $\nabla^{2}\rho$ is $D$, but our choice of time scaling has $D=1$. Note that this term is linear in $\rho$, in spite of it originating from a nonlinear logarithmic contribution to the free energy.

We refer to the model up to this point as DDFT-0.

II.2 Expansion of ${\cal F_{\text{ex}}}$: DDFT-1

To proceed, we must have an expression for the excess Helmholtz free energy functional ${\cal F_{\text{ex}}}[\rho]$. We use a functional Taylor expansion, which is also that used in all derivations of PFC theory. This gives the free energy functional of the system of interest in terms of properties of a reference system, which are assumed to be known. The reference system that is chosen is a uniform liquid, with constant density $\rho_{0}$. The density profile of the system of interest may be varying in space and with an average density that may be different from $\rho_{0}$. The functional Taylor series expansion of the excess free energy can be written in terms of the density difference ${\Delta{\mkern-1.0mu}\rho}({\bm{x}},t)=\rho({\bm{x}},t)-\rho_{0}$ as follows Evans (1992); Hansen and McDonald (1992):

[TABLE]

The expressions $c^{(n)}$ in the equation above are proportional to the first and higher functional derivatives of ${\cal F_{\text{ex}}}$ with respect to density, all evaluated at $\rho=\rho_{0}$:

[TABLE]

These functions $c^{(n)}$ are known as direct correlation functions Evans (1979, 1992); Hansen and McDonald (1992), and are related to $n$-point density correlation functions. In the two-point case, $c^{(2)}$ is the pair direct correlation function and is related to the pair correlation function (i.e., the radial distribution function) through the Ornstein–Zernike equation Evans (1979, 1992); Hansen and McDonald (1992). These direct correlation functions depend on our choice of $\rho_{0}$ and depend directly on temperature through the linear factor of $\beta$ in the definition (15) and also indirectly via the fact that the correlations in a liquid change with temperature. Note also that $c^{(1)}[\rho_{0}]$ is a constant when $\rho_{0}$ is a constant.

For a homogeneous liquid with distant (or periodic) boundaries, these functions depend on displacements but not on absolute position, so (through a slight abuse of notation) we also write

[TABLE]

where ${\Delta{\mkern-1.0mu}{\bm{x}}}_{j}={\bm{x}}_{j}-{\bm{x}}_{1}$ Hansen and McDonald (1992). We also take the liquid to be isotropic.

We are considering density perturbations away from the liquid state, so it is convenient to write

[TABLE]

We do not assume that $n$ is small, but it is often the case that the average of $n({\bm{x}},t)$ over the whole system is small. Note also that $\rho({\bm{x}},t)=\rho_{0}$ is a stationary solution of (5). Substituting Eq. (14) into Eq. (5) and writing only the terms up to $c^{(1)}$, we get:

[TABLE]

That the uniform liquid state is an equilibrium of (5) implies that $n=0$ is a solution of equation (18): all terms not written down involve ${\Delta{\mkern-1.0mu}\rho}$ and so are zero for the uniform liquid with density $\rho_{0}$. Recall that $c^{(1)}[\rho_{0}]$ is a constant, which means terms involving gradients of this can be dropped. Whilst this constant term does not influence the structure (density profile) both in and out of equilibrium, it does affect the thermodynamics (i.e., free energy value) and so also mechanical properties Wang et al. (2018). With this, we can write the equation for the time evolution of $n({\bm{x}},t)$ (up to ${\cal O}(c^{(4)})$) as:

[TABLE]

where we have suppressed writing the time dependence of $n$ throughout and the ${\bm{x}}$ dependence of $n$ when it is not inside an integral. We have written this equation so that the first line is linear in $n$, the next two lines are quadratic in $n$, the fourth and fifth lines are cubic in $n$, and the last line is quartic in $n$.

Since the first line is linear in $n$, and both terms involve a Laplacian, we can write the linearised version of (19) in terms of the negative Laplacian of a linear operator ${\cal L}$:

[TABLE]

where

[TABLE]

The non-local operator ${\cal L}$ is most conveniently considered in terms of its Fourier transform, or equivalently, in terms of how it acts on modes of the form $\exp(i{\bm{k}}\cdot{\bm{x}})$. If

[TABLE]

then $\sigma({\bm{k}})$ is the eigenvalue of ${\cal L}$ with eigenfunction $\exp(i{\bm{k}}\cdot{\bm{x}})$. With this, the linear equation (20) can readily be solved in terms of linear combinations of functions like $\exp\left({k^{2}\sigma({\bm{k}})t+i{\bm{k}}\cdot{\bm{x}}}\right)$, where $k^{2}\sigma({\bm{k}})$ is the growth rate for a mode with wavevector ${\bm{k}}$, and $k=|{\bm{k}}|$. If $\sigma({\bm{k}})$ is negative for all ${\bm{k}}$, all small amplitude density modulations decay to zero, and the liquid state is dynamically stable.

Recall that for a bulk liquid, $c^{(2)}({\bm{x}},{\bm{x}}_{2})=c^{(2)}({\Delta{\mkern-1.0mu}{\bm{x}}}_{2})$, with ${\Delta{\mkern-1.0mu}{\bm{x}}}_{2}={\bm{x}}_{2}-{\bm{x}}$, and for spherically symmetric (isotropic) particles, $c^{(2)}({\bm{x}},{\bm{x}}_{2})=c^{(2)}(|{\Delta{\mkern-1.0mu}{\bm{x}}}_{2}|)$. Therefore, in this case $\sigma({\bm{k}})=\sigma(k)$, i.e., $\sigma$ depends only on the wavenumber $k=|{\bm{k}}|$. The eigenvalue $\sigma(k)$ can be expressed as:

[TABLE]

where ${\hat{c}}^{(2)}$ is the Fourier transform of $c^{(2)}$. Recall from (15) that $c^{(2)}$ is proportional to $\beta$, so if ${\hat{c}}^{(2)}$ has any positive Fourier components, decreasing the temperature (increasing $\beta$) can be expected to lead to a range of wavenumbers with positive growth rates, and the liquid being dynamically unstable to modes with wavenumbers centered around the maximum of $k^{2}\sigma(k)$, see Fig. 1. We have scaled lengths so that the maximum growth rate occurs at wavenumber $k=1$. This argument, of course, assumes that the product $\beta^{-1}c^{(2)}$ is independent of temperature. This is not true in general, but for some systems it is a good approximation (at least over a limited range of temperatures) – see Ref. Somerville et al. (2018) for a recent discussion on this for a particular colloidal system. Recall too that for an equilibrium liquid the static structure factor $S(k)=[1-\rho_{0}{\hat{c}}^{(2)}(k)]^{-1}$. $S(k)$ is proportional to the Fourier transform of the radial distribution function Hansen and McDonald (1992). So, for the stable uniform liquid, we have $\sigma(k)=-1/S(k)$.

We refer to the state point at which the uniform liquid becomes linearly unstable to density modulations with wavenumber $k\neq 0$ as the spinodal point, in keeping with the terminology of Trudu et al. (2006). The more common usage of the term ‘spinodal’ relates to the onset of the zero-wavenumber phase separation instability of liquid–liquid or gas–liquid phase separation Hansen and McDonald (1992); Archer and Rauscher (2004). At the spinodal point, the density and temperature are such that the liquid is dynamically marginally stable, that is, the maximum of $k^{2}\sigma(k)$ is zero. Therefore, at higher temperatures, small amplitude density modulations decay, and at lower temperatures, small amplitude density modulations grow. For a given fixed value of $\rho_{0}$, the spinodal temperature is $T_{s}$, with a corresponding $\beta_{s}=(k_{B}T_{s})^{-1}$. Similarly, either increasing the density $\rho_{0}$ of the liquid or increasing the chemical potential $\mu$ can also lead to crossing the spinodal.

With (21), we can eliminate $c^{(2)}$ in favour of ${\cal L}$ in (19), and obtain (truncating at ${{\cal O}}(c^{(4)})$):

[TABLE]

where we have used the result $\nabla\cdot\left[n\nabla n\right]=\frac{1}{2}\nabla^{2}n^{2}$. For an ideal gas, with ${\cal L}n=-n$ and $c^{(2)}=c^{(3)}=c^{(4)}=0$, the first line of the equation above reduces to the diffusion equation, $\frac{\partial n}{\partial t}=\nabla^{2}n$, similar to (13).

At this point, we have made no approximations beyond expanding the free energy in Eq. (14) and truncating at ${{\cal O}}(c^{(4)})$. We refer to the model at this stage, truncated in this way, as DDFT-1. In the new variables, and incorporating $c^{(2)}$ into ${\cal L}$, the Helmholtz free energy ${\cal F}$ can be expressed (up to fourth order) in terms of a scaled free energy ${{\cal F}_{1}}={\cal F}/\rho_{0}$, where

[TABLE]

and where we have also dropped terms that do not contribute to (24). In these variables, the DDFT that leads to the dynamics (24) is

[TABLE]

Note that, because of the $\log(1+n)$ term in (25), $n$ is constrained so that $1+n$ is always non-negative. Also, because of Eq. (17), we have

[TABLE]

Moreover, states that satisfy

[TABLE]

where ${\Delta{\mkern-0.5mu}\mu}=\mu-\mu_{0}$ and where [see (7), (10) and (14)]

[TABLE]

are equilibrium solutions of (24), or equivalently, extrema of ${{\cal F}_{1}}$. Henceforth, we redefine $\mu$ to be $\beta{\Delta{\mkern-0.5mu}\mu}/\rho_{0}$, which is a shifted and rescaled chemical potential. For the free energy in (25), we have

[TABLE]

At equilibrium, this expression (the rescaled chemical potential $\mu$) does not vary in space. The reference liquid $n=0$ has ${{\cal F}_{1}}=0$ and $\mu=0$. The zero value for ${{\cal F}_{1}}$ arises (in part) from dropping ${\cal F_{\text{ex}}}[\rho_{0}]$ from (14), while the zero value for the rescaled chemical potential is a consequence of (29), which is equivalent to dropping $c^{(1)}$ from (14) and setting $\Lambda=1$ in (10).

II.3 Simplification of $c^{(3)}$ and $c^{(4)}$: DDFT-2

As the next step, Huang et al. Huang et al. (2010) kept only the zero-wavenumber components of $c^{(3)}$ and $c^{(4)}$, or equivalently, they took

[TABLE]

where $c^{(3)}_{0}$ and $c^{(4)}_{0}$ are constants (our sign convention is opposite to that of Huang et al. (2010)). This is equivalent to making a local density approximation (LDA) Evans (1992) for these terms in the free energy. We could in principle include terms involving $c^{(5)}$ and higher as well, treated in the same way: these would contribute a more general function of $n$ in the free energy, treated with the LDA. However, since we are investigating the effect of approximations that have not yet been discussed, we keep as simple a free energy as possible at this point, consistent with truncating at ${\cal O}(c^{(4)})$. With this, the free energy in (25) becomes

[TABLE]

and the four terms involving $c^{(3)}$ and $c^{(4)}$ in (24) become

[TABLE]

Using $\nabla\cdot\left[n\nabla n^{2}\right]=\frac{2}{3}\nabla^{2}n^{3}$ and $\nabla\cdot\left[n\nabla n^{3}\right]=\frac{3}{4}\nabla^{2}n^{4}$, Huang et al. Huang et al. (2010) combined (33) and (24) to get

[TABLE]

where

[TABLE]

We also have a chemical potential

[TABLE]

which does not vary in space at equilibrium. Up to this point, we refer to the model as DDFT-2.

Here, we retain the $n^{4}$ term (as did Huang et al. Huang et al. (2010)), because otherwise the dynamics in (34) would not be consistent with the free energy (32) and the DDFT dynamics (26) (with ${\cal F}_{2}$ instead of ${\cal F}_{1}$).

The next three models involve making (or not making) two approximations: (i) assuming the Ramakrishan–Yussouff or random phase approximation, which leads to a quadratic excess Helmholtz free energy functional, and (ii) making a gradient expansion of the linear operator ${\cal L}$.

II.4 Quadratic excess free energy: DDFT-3

Often, the free energy functional in (14) is truncated at ${\cal O}({\Delta{\mkern-1.0mu}\rho}^{2})$. This is known as the Ramakrishan–Yussouff (RY) approximation Ramakrishnan and Yussouff (1979); van Teeffelen et al. (2009); Emmerich et al. (2012), which effectively sets $c^{(3)}=c^{(4)}=0$. A mathematically equivalent approximation arises in the treatment of soft purely repulsive particles modelling soft matter, namely the RPA Likos (2001). Here, two soft isotropic particles at ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$ separated by a distance $x_{12}=|{\bm{x}}_{1}-{\bm{x}}_{2}|$ interact through a potential energy $u(x_{12})$, which depends only on the magnitude of the distance and is finite for all values of $x_{12}$. The excess free energy [c.f. Eq. (14)] is then

[TABLE]

This amounts to setting $c^{(3)}=c^{(4)}=0$ and

[TABLE]

in DDFT-2. The eigenvalues $\sigma(k)$ can thus be related to the Fourier transform of $u$ through (23) Likos (2001); Archer et al. (2012):

[TABLE]

Setting $c^{(3)}=c^{(4)}=0$ implies from (35) that $Q=\frac{1}{2}$, $C=0$ and $R=0$, and results in a free energy

[TABLE]

With this choice of free energy, the dynamics in (34) becomes:

[TABLE]

along with an analogous version of (36), for the chemical potential:

[TABLE]

We refer to this model as DDFT-3; it is equivalent to DDFT-1 with the RY approximation, and to DDFT-0 with ${\cal F_{\text{ex}}}$ given by the RPA approximation.

Before moving on to make further approximations, it is worth noting a useful property that DDFT-3 and the subsequent theories derived from it possess. If the pair potential $u(x_{12})$ in Eq. (38) can be written as $u(x_{12})=\epsilon\psi(x_{12})$, where $\epsilon$ is a parameter that controls the overall strength of the potential, then from Eqs. (21), (38) and (42) we obtain:

[TABLE]

The consequence of this is that for a given $\psi$, the behaviour of the model depends only on the combination of parameters $\rho_{0}\beta\epsilon$ and the value of $\mu$. If one changes the value of the reference density $\rho_{0}$ to some other value, then this is entirely equivalent to solving the system with the original reference density $\rho_{0}$ at a different value of $\beta\epsilon$. We should emphasize that this is only true if $\psi$ does not change with density, which in general is not true, but is approximately the case for some systems.

II.5 Gradient expansion of the linear term: DDFT-4

Returning to DDFT-2, Huang et al. Huang et al. (2010) (following Elder and Grant (2004)) expanded ${\cal L}$ in powers of the gradient operator $\nabla$, replacing ${\cal L}$ by the simplest linear operator that allows a positive growth rate for modes with a wavenumber $k_{s}$. Scaling lengths so that $k_{s}=1$ results in:

[TABLE]

so $\sigma(k)=r-\gamma(1-k^{2})^{2}$ from (22). This approximation is equivalent (within scaling) to a local gradient expansion of (21), expanding the Fourier transform of $c^{(2)}$ about its maximum:

[TABLE]

where the function $\rho_{0}{\hat{c}}^{(2)}(k)$ and its second derivative evaluated at $k=1$ are $1+r$ and $-8\gamma$, respectively. Here, $r$ is a parameter, notionally increasing with $\beta$ (and with $\rho_{0}$) and equal to zero at the spinodal point, when $\beta=\beta_{s}$. This parameter controls the growth rate of waves with wavenumber $1$: effectively, $r$ is the height of the maximum at $k=1$ in the growth rate curve in Fig. 1. The second parameter $\gamma$ can be used to fit the curvature of ${\hat{c}}^{(2)}(k)$ at $k=1$.

With this gradient expansion, the dynamics is

[TABLE]

We refer to this model as DDFT-4: ${{\cal L}_{\text{grad}}}$ is now a (local) partial differential operator and (46) is a partial differential equation. The free energy and chemical potential can be found from (32) and (36), setting ${\cal L}={{\cal L}_{\text{grad}}}$. The lower bound $n\geq-1$ is still respected. This model is equivalent to that written down by Huang et al. (2010).

Higher powers (or other functions) of the Laplacian can be retained in ${{\cal L}_{\text{grad}}}$, to improve the accuracy of the match between the eigenvalues of ${\cal L}$ and ${{\cal L}_{\text{grad}}}$, as done for example by Jaatinen et al. (2009); Pisutha-Arnond et al. (2013), or to introduce additional unstable length scales, as done for example by Wu et al. (2010); Barkan et al. (2011); Achim et al. (2014); Subramanian et al. (2016) and others. See also Eq. (76) below and the associated discussion.

II.6 RY and gradient expansion: DDFT-5

Finally, we can make both the RY/RPA approximation ($c^{(3)}_{0}=c^{(4)}=0$) and replace the linear operator ${\cal L}$ by ${{\cal L}_{\text{grad}}}$ to get the model advocated in Ref. van Teeffelen et al. (2009). The free energy and evolution equation are

[TABLE]

and

[TABLE]

along with an analogous version of Eq. (42) for the chemical potential:

[TABLE]

This model is named PFC1 in van Teeffelen et al. (2009), but here we call it DDFT-5 for consistency.

II.7 PFC models

The final simplification that can be made (or not made) is to discard the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ (or $\nabla\cdot\left[n\nabla{{\cal L}_{\text{grad}}}{n}\right]$) term from the dynamical equations for the four DDFT models DDFT-2, …, DDFT-5, resulting in four PFC models PFC-$\beta$, …, PFC-$\epsilon$. Huang et al. Huang et al. (2010) justify making this simplification on the grounds that this term is not truly quadratic in $n$: the presence of ${\cal L}{n}$ in the expression means that it is effectively of higher order. However, we show below that this term does in fact make an important contribution to the free energy: at least as important as the $c^{(3)}$ term.

In addition, dropping this term implies significant changes to the DDFT dynamics, the mobility and the nonlinear terms in the free energy. In fact, the $(1+n)$ factor in the mobility in (26), the logarithm in the ideal gas free energy in (10) and the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term in (24) are inextricably linked. This can be seen in the progression from (10) to (13): the functional derivative of the ideal gas term in (10) (the first term on the right hand side) leads to the $\log\rho$ term in (11), the gradient of this leads to $\rho^{-1}\nabla\rho$ in (12), and the mobility being $M=D\rho$ cancels the $\rho^{-1}$, leading to a diffusion equation in (13). If the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term is dropped from (34), the equation for $n$ becomes of the form $\frac{\partial n}{\partial t}=\nabla^{2}\frac{{\delta{\mkern-1.0mu}{\cal G}}}{{\delta{\mkern-1.0mu}n}}$ for some functional ${\cal G}[n]$. We can see the implications of this by returning to (5) and taking the steps needed to get to this modified version of (34). Clearly the mobility in (5) has been taken to be constant. If we now think of the ideal gas part of the free energy in (10) and (11), but with a constant mobility in the dynamical equation, we end up with the ideal gas term contribution to the equation for $\rho$ being the form

[TABLE]

instead of the diffusion equation (13). This unlikely equation can be avoided, and the diffusion equation recovered at leading order, by expanding the logarithm in (32) in a Taylor series. Thus, dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term is equivalent to taking constant mobility and expanding the logarithm.

It is because of these substantial changes that we opt to use the term ‘DDFT’ for all models based on free energies that have the logarithmic ideal gas term, the non-constant mobility and the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term retained. In contrast, we use the term ‘PFC’ for models based on expanding the logarithm, having a constant mobility and the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term dropped. One consequence of expanding the logarithm up to ${\cal O}(n^{4})$, as is done in most PFC derivations Emmerich et al. (2012), is that the ideal gas part of the free energy contributes cubic and quartic (as well as quadratic) terms to the free energy, so going from DDFT-2 to PFC-$\beta$ turns out not to be just a matter of dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term.

So, a consistent free energy–dynamics derivation Elder and Grant (2004); van Teeffelen et al. (2009) involves going back to DDFT-2 and replacing the logarithm in (32) by:

[TABLE]

resulting in a free energy

[TABLE]

where we have suppressed writing the ${\bm{x}}_{1}$ dependency of $n({\bm{x}}_{1})$. Taking the mobility $M(\rho)$ in (5) to be a constant ($M=D\rho_{0}$) implies (after scaling)

[TABLE]

similar to (2). This leads to the PFC dynamical equation:

[TABLE]

and to a chemical potential

[TABLE]

where $Q$ is as in (35) but $C$ is different:

[TABLE]

We refer to this model as PFC-$\beta$, and recall that the factor of $\beta$ in front of {\cal F_{\text{\beta}}} is the inverse temperature.

The end result here is that PFC-$\beta$ (54) is not the same as DDFT-2 (34) with the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term removed: the cubic coefficient $C$ is different and the quartic contribution $Rn^{4}$ in (34) is absent. For the cubic coefficient, the contribution proportional to $c^{(3)}_{0}$ in (35) comes from the non-constant mobility, while the $-\frac{1}{3}$ term in (56) comes from expanding the logarithm. The contribution to $C$ proportional to $c^{(4)}_{0}$ is the same. Moreover, the $\frac{1}{2}$ in $Q$ in (35) and (56), while having the same numerical value, arises for two different reasons: non-constant mobility versus expanding the logarithm. An additional difference between the DDFT and PFC models is that in the PFC models, the constraint that $n\geq-1$ (i.e., $\rho\geq 0$) is not enforced.

As in the DDFT derivations, we can now make (or not make) the RY/RPA approximation and the gradient expansion. We consider first the RY/RPA approximation, setting $c^{(3)}=c^{(4)}=0$ in PFC-$\beta$. The free energy is

[TABLE]

the dynamics is

[TABLE]

and the chemical potential is

[TABLE]

We refer to this model as PFC-$\gamma$, and it is effectively the same as PFC-$\beta$ but with $Q=\frac{1}{2}$ and $C=-\frac{1}{3}$.

Finally, the gradient expansion can be made, replacing ${\cal L}$ by ${{\cal L}_{\text{grad}}}$ in all expressions in this subsection, resulting in PFC-$\delta$ (without RY/RPA) and PFC-$\epsilon$ (with RY/RPA).

We refer to these models collectively as the PFC models, and have chosen the names PFC-$\alpha$ etc. to distinguish these from the PFC1 and PFC2 models of Ref. van Teeffelen et al. (2009). The quadratic term in the dynamics ($Qn^{2}$) can be removed (provided $C\neq 0$) by adding a constant to $n({\bm{x}})$, but we choose not to do this as it implies a change to what was meant by $\rho_{0}$ in the reference liquid. In addition, a negative $C$ can be scaled to $-1$. With these changes, PFC-$\epsilon$ is equivalent to the original PFC-$\alpha$ model (2) of Elder et al. (2002); Elder and Grant (2004):

[TABLE]

where we have written out ${{\cal L}_{\text{grad}}}$ explicitly, and $Q=\frac{1}{2}$ and $C=-\frac{1}{3}$ (or $Q=0$ and $C=-1$ after scaling and adding a constant to $n$ – returning to the conserved Swift–Hohenberg equation).

The implication of dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term in the dynamics (41) for DDFT-3 is now apparent: without this term, Eq. (41) reduces to (58) but with the cubic term removed. The absence of the cubic term here implies a free energy as in (57) that is not bounded below, i.e., a free energy that is non-physical, and so the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term can have the effect of stabilizing patterns. In addition, dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term is consistent only with a theory with a constant mobility.

II.8 Summary

To summarise, we have carefully laid out the various approximations made in the progression from the DDFT-0 starting point (5) to the final PFC (2,58) written down in Elder et al. (2002); Elder and Grant (2004). We have largely followed earlier derivations van Teeffelen et al. (2009); Huang et al. (2010); Emmerich et al. (2012), seeking to clarify the approximations that are made. Along the way, we have identified four intermediate versions of DDFT, listed for clarity in Table 1. The change in name from DDFT to PFC could be made at any point in this progression, but we prefer to make the name change at the point where the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term is dropped (along with all the other changes that are implied by this), since removing this term marks a considerable alteration to the free energy expression and to the dynamics.

The PFC model (54) is appealing in its simplicity, and it gives insight into a variety of crystallisation phenomena, but the derivations of the model from DDFT presented here, as well as the derivation from Ref. Huang et al. (2010), are both problematic. Just dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term, as done by Huang et al. Huang et al. (2010), means that the dynamics is not equivalent to a DDFT with mobility proportional to $\rho$. On the other hand, the alternative is to expand the logarithm up to ${\cal O}(n^{4})$ in (51) in order to provide a nonlinear stabilizing term ($\frac{1}{12}n^{4}$) in the free energy (57). However, in the original formulation, the logarithm comes from the ideal gas term in (10), and leads to a linear diffusive term in the dynamics. The stabilizing nonlinear terms in (34) are provided by $c^{(3)}_{0}$, $c^{(4)}_{0}$ (in DDFT-2) and by the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term – these are all absent in PFC-$\gamma$.

Indeed, all these models only make physical sense if their free energies are bounded below. The free energies for DDFT-0 and DDFT-1 are too general to make any comment, but that for DDFT-2 (32) etc. can be discussed. The $(1+n)\log(1+n)-n$ term is bounded below by zero, and the $-\tfrac{1}{2}\int{n}{\cal L}{n}{\rm d}{\bm{x}}$ term is bounded below because the eigenvalues of ${\cal L}$ are bounded above:

[TABLE]

where $\sigma_{\text{max}}$ is the maximum over $k$ of $\sigma(k)$ (we have in mind a $\sigma(k)$ as in Fig. 1). In any case, this term, along with the other quadratic and cubic terms, is dominated by the quartic in $n$, which is bounded below provided $c^{(4)}_{0}<0$. If $c^{(4)}_{0}=0$, then $c^{(3)}_{0}<0$ will do, recalling that $n\geq-1$. For DDFT-3, with the RY/RPA approximation $c^{(3)}_{0}=c^{(4)}=0$, the boundedness of the free energy (40) depends on the $n^{2}+n{\cal L}{n}$ combination. From (21), the relevant term is

[TABLE]

In general, this is not bounded below, but it is in certain circumstances. For example, it is if $\sigma_{\text{max}}<-1$, and it is if $c^{(2)}({\bm{x}}_{1},{\bm{x}}_{2})\leq 0$ (or $u(|{\bm{x}}_{1}-{\bm{x}}_{2}|)\geq 0$ for RPA) for all ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$, which is the case in the numerical examples below. The PFC models are not constrained to have $n\geq-1$, but {\cal F_{\text{\beta}}} (52) is bounded by the $n^{4}$ term as long as its coefficient is positive; {\cal F_{\text{\gamma}}} (57) is always bounded below, because the expansion of the logarithm in (51) was truncated after an even powered term.

Throughout we have made the simplest choices in the approximations, but other authors have made many other choices. For example, the original PFC paper Elder et al. (2002), as well as later papers Elder et al. (2007); Huang et al. (2010); Robbins et al. (2012); Alster et al. (2017); Elder et al. (2018), included a two-component (binary) version of the PFC model. Recently, Wang et al. Wang et al. (2018) took a much closer look at $c^{(3)}$ and $c^{(4)}$, expressing these in terms of isotropic tensors and so allowing these functions to introduce bond angle dependence into the free energy. Some choices of $c^{(3)}$ and $c^{(4)}$ lead to nonlinear terms that include gradients, which can affect the selection of the final stable crystal Wu et al. (2010). The gradient expansion approach has been generalised in two ways: (i) higher order terms or rational functions were considered by Jaatinen et al. (2009); Pisutha-Arnond et al. (2013) in order to improve the fit between the functional form and the Fourier transform of $c^{(2)}$, and (ii) PFC models with two unstable length scales have been put forward by several authors Pisutha-Arnond et al. (2013); Wu et al. (2010); Barkan et al. (2011); Achim et al. (2014); Subramanian et al. (2016); Jiang et al. (2017), since these allow more complex crystals (face-centered cubic, icosahedral quasicrystals, …) to be stabilized. We discuss the model of Jaatinen et al. (2009) in more detail below. Alternative approaches involving weighted densities are also possible Jaatinen and Ala-Nissila (2010).

III Comparison of DDFT and PFC

We are interested in the effects of the approximations made in going from DDFT to PFC. A full assessment of the validity of the RY/RPA approximation for ${\cal F_{\text{ex}}}$, which in itself constitutes a major simplification, is beyond the scope of the present study. The general conclusion on the validity of the RY/RPA approximation is that it depends on the nature of the interactions between the particles; there are examples in the literature where this approximation is reliable and others where it works badly – see for example the discussion in Refs Evans (2010); Löwen (2010) and references therein.

Here, we consider one particular system where the RPA is accurate and then we focus on the effects of approximating ${\cal L}$ by ${{\cal L}_{\text{grad}}}$ and of making the suite of other approximations inherent in going from DDFT to PFC: expanding the logarithm, assuming constant mobility and dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term. To this end, we start with DDFT-3 and solve (42), rewritten here as:

[TABLE]

The system that we consider is particles interacting via the GEM-4 Mladek et al. (2006); Prestipino and Saija (2014) potential: this is a model for soft-matter particles and in particular for dendrimers and other polymers in suspension, treating the polymers via an effective pair potential between their centers of mass. This potential is soft, i.e., finite for all values of $x_{12}$ Mladek et al. (2006); Prestipino and Saija (2014); Likos (2001, 2006); Lenz et al. (2012), and is

[TABLE]

where the parameter $\epsilon$ controls the strength of the potential and $R$ controls its spatial range. We consider here the system in 2D Prestipino and Saija (2014); Archer et al. (2014). As long as the temperature and density are high enough that the particle cores regularly overlap (the regime in which the system freezes), the RPA approximation (37) is known to be rather accurate for the GEM-4 system and gives a good account of the phase diagram and the structure of the liquid and solid phases Mladek et al. (2007); Prestipino and Saija (2014); Archer et al. (2014).

From Eqs. (21), (38) and (64) we obtain the linear operator ${\cal L}$:

[TABLE]

Recall from (22) that ${\cal L}$ has eigenvalues $\sigma(k)$ with eigenfunctions $e^{i{\bm{k}}\cdot{\bm{x}}}$. We can choose the combined parameter $\rho_{0}\beta\epsilon$ and soft-particle radius $R$ so that the maximum in $\sigma(k)$ occurs at $k=1$ when the system is at the linear stability threshold, i.e., this is a maximum with $\sigma(1)=0$, similar to Fig. 1. In 2D, to satisfy this condition we must have $\rho_{0}\beta\epsilon=0.2455$ and $R=5.0962$ – see Appendix A for details.

With this choice of parameters, the eigenvalue $\sigma(k)$ is shown as a solid line in Fig. 2. The figure also shows (dashed line) the eigenvalue for the gradient expansion of ${\cal L}$ around $k=1$:

[TABLE]

where $\gamma=4.37$ is chosen to match the second derivative $\frac{{\rm d}^{2}\sigma}{{\rm d}{k}^{2}}$ at $k=1$, as done for example in Refs. Elder and Grant (2004); Wu and Karma (2007); Jaatinen et al. (2009). The dotted line in Fig. 2 is the eigenvalue for (76), the eighth-order fitting model proposed in Jaatinen et al. (2009) and discussed in more detail below.

In what follows we compare solutions of (63) for DDFT-3 with solutions of the analogous equations for DDFT-5, PFC-$\gamma$ and PFC-$\epsilon$:

[TABLE]

See Appendix B for details of the pseudo-arclength continuation numerical method we use for solving these equations. Also, in the supplementary material we include a Matlab code for solving DDFT-3. Note that throughout what follows, we refer to the quantity $1+n({\bm{x}})$ as the ‘density’.

Since the DDFT-5, PFC-$\gamma$ and PFC-$\epsilon$ represent different forms of Taylor expansion around the reference state with density $\rho_{0}$, there are a variety of ways comparison between solutions can be made. Here, we opt to fix ${\cal L}$ and ${{\cal L}_{\text{grad}}}$ as in (65) and (66) with the specified values of $\rho_{0}\beta\epsilon$, $R$ and $\gamma$. This implies that at $\mu=0$ the reference state with $n=0$ is at the spinodal point and is marginally unstable to modes with wavenumber $k=1$. We then vary $\mu$ starting from $\mu=0$ and follow the liquid, stripe and hexagonal solutions of (63) and (67)–(69) in appropriately sized two-dimensional domains. For a given value of $\mu$ the different solutions have different values for the mean density $1+\bar{n}=1+\frac{1}{A}\!\!\int\!n({\bm{x}}){\rm d}{\bm{x}}$, where $A$ is the area of the domain. For each state we calculate the specific grand potential:

[TABLE]

where ${\cal F}$ is ${{\cal F}_{3}}$, ${{\cal F}_{5}}$, {\cal F_{\text{\gamma}}} or {\cal F_{\text{\epsilon}}}, as appropriate. We also minimise $\Omega/A$ with respect to the domain size $A$ by applying the approach described in Appendix B. For a given value of the chemical potential $\mu$ and the combined parameter $\rho_{0}\beta\epsilon$, the thermodynamic equilibrium state is that with the minimum value of $\Omega/A$. Note that equilibria with the same $\mu$ do not necessarily have the same value of $\bar{n}$, which is important when considering which equilibria might result from initial conditions via the dynamics.

The solution corresponding to the uniform density liquid state with $n({\bm{x}})=n_{\text{liq}}$ can readily be found. In this case we have ${\cal L}{n_{\text{liq}}}=\sigma(0)n_{\text{liq}}$, and so we must solve the following algebraic equations for $n_{\text{liq}}$:

[TABLE]

recalling that the value of $\sigma(0)$ depends on whether or not the gradient expansion is carried out (see Fig. 2). Finding $n_{\text{liq}}$ for a given value of $\mu$ is done easily using Newton’s method, and the resulting $n_{\text{liq}}$ and specific $\Omega_{\text{liq}}$ are shown in Fig. 3. In all cases, we see that $n_{\text{liq}}$ is an increasing function of $\mu$, while $\Omega_{\text{liq}}/A$ is a decreasing function of $\mu$. The figure shows that the specific grand potential for the liquid state predicted by all four models are similar close to $\mu=0$, but the predicted liquid state densities are rather different away from $\mu=0$. This difference originates from the different values of $\sigma(0)$ ($-18.75$ for DDFT-3 and PFC-$\gamma$, in contrast to $-4.37$ for DDFT-5 and PFC-$\epsilon$). We see from Fig. 3 that the density of the liquid is erroneously predicted to increase too rapidly as $\mu$ is increased by the gradient expansion theories (DDFT-5 and PFC-$\epsilon$). This is because these get the value of the isothermal compressibility $\chi_{T}$ to be too large Jaatinen et al. (2009). This compressibility is related to $\sigma(0)$ via $\chi_{T}=-\beta/[\sigma(0)\rho_{0}(1+n_{\text{liq}})]$ Hansen and McDonald (1992): see Eq. (23) and following discussion. Expanding the logarithm makes relatively little difference over this range of densities.

Since crystallisation occurs at higher densities, we expect a transition from the liquid to the crystal to occur as $\mu$ increases. At the spinodal the uniform liquid becomes linearly unstable and the patterned state solution branches bifurcate from the liquid at this point. To find these states, we seek a solution of the form

[TABLE]

where near the bifurcation point ${\delta{\mkern-1.0mu}n}\ll 1$, and ${\delta{\mkern-1.0mu}n}$ is of the form $e^{i{\bm{k}}\cdot{\bm{x}}}$, so that ${\cal L}{\delta{\mkern-1.0mu}n}=\sigma(k){\delta{\mkern-1.0mu}n}$. Expanding Eqs. (63) and (67–69) in powers of ${\delta{\mkern-1.0mu}n}$ we find that the ${\cal O}(1)$ equations to solve are just those for finding the liquid state density, Eqs. (71)–(72). The ${\cal O}({\delta{\mkern-1.0mu}n})$ equations are

[TABLE]

The spinodal point for DDFT-3,5 or for PFC-$\gamma$,$\epsilon$ is where there are solutions of the equation with ${\delta{\mkern-1.0mu}n}\neq 0$.

Since we are looking for a change in stability, we take the extreme value of $\sigma(k)$, i.e., $\sigma(k)=0$ (see Fig. 2). Then, Eq. (74) is solved (with ${\delta{\mkern-1.0mu}n}\neq 0$) only for $n_{\text{liq}}=0$, which leads to $\mu=0$ from (71). In contrast, Eq. (75) with $\sigma(k)=0$ has two solutions, $n_{\text{liq}}=0$ and $n_{\text{liq}}=1$, leading to $\mu=0$ and $\mu=-\tfrac{1}{6}-\sigma(0)$ from (72). The implication of this is that the PFC has two spinodal points: the liquid loses stability at $n_{\text{liq}}=0$ as $\mu$ increases through [math], but it regains stability at $n_{\text{liq}}=1$, which gives $\mu=18.58$ for PFC-$\gamma$ and $\mu=4.20$ for PFC-$\epsilon$. This prediction that the liquid regains stability for higher $\mu$ is a consequence of expanding the logarithm, or equivalently of Taylor expanding the $1/(1+n_{\text{liq}})$ term in (74) and is confirmed by direct computation of the crystal solutions below. Of course, this prediction is erroneous, since the simulation results for the GEM-4 system Mladek et al. (2006); Prestipino and Saija (2014) show no sign of a second spinodal point or the associated stable second liquid in the equilibrium system phase diagram.

In Fig. 4 we display examples of the three different types of periodic solutions that can be found for DDFT-3. These are (i) the crystal solution, which we refer to as ‘up hexagons’, which exhibits a triangular array of isolated density maxima surrounded by hexagonal regions where the density is close to zero. There are also (ii) ‘down hexagons’ which are the opposite, with isolated density minima and hexagonal density maxima. Finally, there is (iii) the stripe state. Depending on the state point these solutions are not necessarily linearly stable. Our naming convention to distinguish the two different hexagonal solutions originates in the convection literature Bodenschatz et al. (2000). These solutions were initiated at $\mu\approx 0$ and then continued numerically (see Appendix B) up to $\mu=10$. For DDFT-3 it is possible to go a bit higher in $\mu$, but with increasing $\mu$ (i.e., increasing average density) the peaks in the density profile get narrower and higher and so more and more grid points are required to resolve these peaks correctly. However, as we show below for some of the other models and for different reasons, it is not possible to continue the solutions this far in $\mu$. The domains on which the profiles are calculated have periodic boundary conditions, with 4 wavelengths in each direction (for stripes), or $8\times\tfrac{8}{\sqrt{3}}$ wavelengths (for hexagons). The wavelength is initially equal to $2\pi$ for $\mu=0$ and is then adjusted by up to about 2% in order to minimise the specific grand potential as $\mu$ is varied; i.e., we minimise $\Omega/A$ with respect to variations in the size of the crystal unit cell or, for the stripe phase, we minimise with respect to variations in the spacing between the stripes – see Appendix B for details.

In Fig. 5 we display a series of plots showing the maximum, minimum and average values of the density profiles $1+n$ for the stripe and hexagonal structures as a function of $\mu$. We also plot the specific grand potential $\Omega/A$ for the different structures. Recall that for a given $\mu$ the thermodynamic equilibrium phase corresponds to the global minimum of $\Omega/A$. The results for DDFT-3 are shown in Fig. 5(a–c). The (a) stripes originate in a supercritical pitchfork bifurcation at $\mu=0$, and (b) hexagons originate in a transcritical bifurcation at the same value of $\mu$. The density of the up hexagons ranges from about $2\times 10^{-5}$ up to about 50, for $\mu=10$. All of these branches can be continued to larger values of $\mu$.

DDFT-5, in Fig. 5(d–f), initially behaves in the same way, but all three branches have their minimum density heading to zero before $\mu$ gets to 10: this happens at $\mu\approx 3.37$ for (d) stripes, and for $\mu\approx 0.28$ and $\mu\approx 2.73$ for (e) the up and down branches of hexagons, respectively. The numerical method cannot continue the branches beyond these points. We argue in Sec. IV that this is not an artefact of the numerical method, rather it is a genuine feature of solutions of Eq. (67) that the density $1+n$ can go to zero. In this limit, $\log(1+n)\rightarrow-\infty$, but this is balanced by a lack of smoothness in $n({\bm{x}})$: the fourth derivative in ${{\cal L}_{\text{grad}}}n$ can go to $+\infty$ and so balance the singularity in $\log(1+n)$. Therefore, $\mu\approx 0.28$ is the limit of validity of the DDFT-5 model.

The two PFC examples are similar to each other, and it is easier to discuss PFC-$\epsilon$, in Fig. 5(j–l), first. Here, (j) stripes and (k) hexagons bifurcate from the liquid at $\mu=0$, but they rejoin the liquid at $\mu=4.20$ as explained in the discussion following Eqs. (74–75). The maximum and minimum densities for the up and down hexagon cross between the two bifurcations. The behaviour of PFC-$\gamma$, in Fig. 5(g–i), is similar, though the second bifurcation is at $\mu=18.58$, off the scale of the figure.

Figure 5(c,f,i,l) shows that the curves of the specific grand potential $\Omega/A$ as functions of $\mu$ are very close, so in Fig. 6, we plot instead $(\Omega-\Omega_{\text{liq}})/A$ versus $\mu$, where $\Omega_{\text{liq}}$ is the specific grand potential for the liquid at the same value of $\mu$. In (a) the DDFT-3 case, the up hexagons clearly have the lowest grand potential for $\mu\geq-2.8$ with the uniform liquid being the global minimum for $\mu<-2.8$, and at no point do stripes come anywhere near, as one should expect from the particle simulation results Prestipino and Saija (2014). For (b) DDFT-5, the two hexagon branches stop before the stripe branch when their minimum densities go to zero (the limit of validity), but otherwise the relative values for the hexagon and stripe grand potentials is qualitatively similar to DDFT-3. For (c,d) the PFC examples, once again it is easiest to discuss PFC-$\epsilon$ first. In Fig. 6(d), the hexagon and stripe branches bifurcate from the liquid at $\mu=0$ and rejoin the liquid at $\mu=4.20$, with stripes having the lowest grand potential for intermediate values of $\mu$, and up or down hexagons being the lowest grand potential state for smaller or larger values of $\mu$. The behaviour of (c) PFC-$\gamma$ is similar, but stretched to larger values of $\mu$ (off scale).

The insets in the four panels of Fig. 6 display magnifications that show that the behaviour near the spinodal point at $\mu=0$ is qualitatively similar in all four cases: the up hexagons start with $\Omega>\Omega_{\text{liq}}$ for negative $\mu$, but the branch changes direction, forming a cusp, close to which is the thermodynamic coexistence point (Maxwell point), where $\Omega=\Omega_{\text{liq}}$. The down hexagons start with $\Omega<\Omega_{\text{liq}}$ for positive $\mu$, and the stripes, also with $\Omega<\Omega_{\text{liq}}$ for positive $\mu$, have a value of the grand potential intermediate between the up and down hexagons. We note that the range of $\mu$ over which this behaviour occurs is about a factor of ten smaller in the DDFT-5 and PFC-$\epsilon$ cases as compared to DDFT-3, also with a roughly ten-fold drop in the overall range of values of $(\Omega-\Omega_{\text{liq}})/A$.

The observation that the bulk phase behaviour of the system depends only on $\mu$ and the value of $\rho_{0}\beta\epsilon$ if the pair potential can be written as $u(x_{12})=\epsilon\psi(x_{12})$ – see the discussion around Eq. (43) – is true for the GEM-4 system. As a consequence, having calculated the coexisting densities for a particular value of $\rho_{0}\beta\epsilon$, the linear stability threshold, these results can be scaled to give the phase diagram in the full average density $1+\bar{n}$ versus dimensionless temperature $k_{B}T/\epsilon$ plane, which is one of the usual ways the GEM-4 phase diagram is displayed Mladek et al. (2006, 2007); Prestipino and Saija (2014); Archer et al. (2014); Archer and Malijevskỳ (2016). The phase diagrams obtained from doing this are in Fig. 7. In Fig. 7(a) we display the phase diagram obtained from DDFT-3, which is identical (to within the resolution of the calculations) to that previously calculated in Refs. Archer et al. (2014); Archer and Malijevskỳ (2016). For example, when $\beta\epsilon=1$, the average densities $\rho_{0}(1+\bar{n})R^{2}$ of the coexisting liquid and the crystal are $5.41$ and $5.68$, respectively. As a result of the scaling behaviour, the coexisting densities (binodals) are two straight lines going from the origin and passing through these two points.

In Fig. 7(b) we display the phase diagram obtained for the DDFT-5. The binodals are a little closer to the linear stability threshold line than for DDFT-3, but other than that, it looks similar overall. Note however that the up hexagon branch cannot be continued beyond $\mu\approx 0.28$ (where the minimum density goes to zero): this line is indicated as the ‘limit of validity’. Beyond this line, in the bottom right region of the phase diagram, there is no up hexagon solution to the equations, for the reasons discussed in Sec. IV.

In Fig. 7(c) we display the phase diagram for PFC-$\gamma$ and in (d) for PFC-$\epsilon$. The binodals almost overlie each other, so the predicted difference between the average densities of the liquid and the crystal at coexistence are much smaller than that predicted by DDFT-3 and DDFT-5. Furthermore, on moving to higher average densities or to lower temperatures $k_{B}T/\epsilon$ one encounters the stripe phase, followed by the down hexagon phase and then finally the uniform liquid becomes stable again. The prediction of the occurrence of these later phases is of course wrong, signifying a breakdown in the accuracy of the PFC theory at even the qualitative level.

Before finishing this section, we note that it is possible to extend the gradient expansion in (66) by including higher powers of the Laplacian. For example, Ref. Jaatinen et al. (2009) proposed an eighth-order fitting (EOF), which in our notation is

[TABLE]

where $\gamma$ fits the curvature of the dispersion relation as before, and $E_{B}$ allows the eigenvalue $\sigma(0)$ of ${\cal L}$ to be matched as well, i.e., allows the model to match correctly the isothermal compressibility $\chi_{T}$. An example of the dispersion relation for this operator is shown as a dotted line in Fig. 2. This EOF version of the theory, with ${{\cal L}_{\text{grad-8}}}$ in (76), improves over the standard version, with (66), since $\sigma(0)$ for ${\cal L}$ and ${{\cal L}_{\text{grad-8}}}$ are the same. Therefore, the liquid properties of DDFT-5 match those of DDFT-3, and the liquid properties of PFC-$\gamma$ match those of PFC-$\epsilon$, once ${{\cal L}_{\text{grad}}}$ is replaced by ${{\cal L}_{\text{grad-8}}}$.

However, the drawbacks of the gradient expansion are still present. With ${{\cal L}_{\text{grad-8}}}$, the values of $\mu$ at which the DDFT-5 stripe and hexagon densities go to zero are larger, but this undesirable feature is only deferred, not eliminated. The reason is that the singularity in the logarithm is now balanced against an eighth-order derivative. Note too that introducing even higher derivatives does not cure this problem, it just pushes the singularity to higher order. In addition, with ${{\cal L}_{\text{grad-8}}}$ the second liquid spinodal in the PFC-$\epsilon$ model is still present, it is just pushed to higher values of $\mu$ (similar to the value for PFC-$\gamma$) and since this second spinodal is present, there is still a range of values of $\mu$ for which stripes have the lowest specific $\Omega$.

IV Effect of the approximations

The qualitative change in going from a DDFT to a PFC model (dropping the $\nabla\cdot\left[n\nabla{\cal L}{n}\right]$ term, assuming constant mobility and expanding the logarithm) is apparent in the phase diagrams shown in Fig. 7, comparing (a,b) to (c,d). Here we discuss in detail additional effects of the approximations.

IV.1 Expanding the logarithm

The effect of expanding the logarithm as in Eq. (51) is very significant. In Sec. III we demonstrated that this expansion leads to the liquid having a second spinodal point, illustrated in the phase diagrams in Fig. 7(c,d), with the crystal re-melting as the density is increased. The reason for this is that Eq. (75) can be solved with $n=0$ and $n=1$, while Eq. (74) is only solved by $n=0$, with $\sigma(k)=0$ in both cases.

Intimately connected to the existence of this second spinodal is the transition from stable up hexagons at the $\mu=0$ spinodal to stable down hexagons at the higher density spinodal. These connections to the spinodals mean that the free energy of the up hexagons increases again compared to the liquid state free energy as the chemical potential is increased, in order to reconnect to the liquid state at the upper spinodal. An intermediate region of stable stripes is not inevitable, but is evident in both PFC examples in Fig. 6(c,d).

Taking more terms in the expansion in Eq. (51) does not help. The highest power should be even (otherwise the free energy is not bounded below), and the improved versions of Eq. (75), which involves the second derivative of Eq. (51) with respect to $n$, also have $n=0$ and $n=1$ as (the only) real roots, regardless of how many terms are kept in the expansion of the logarithm. The exception is if only terms up to $n^{2}$ are kept in (51); in this case, $c^{(3)}$ and $c^{(4)}$ (if they are non-zero) provide the stabilizing nonlinearities and may also lead to a spurious spinodal.

IV.2 Gradient expansion of ${\cal L}$

As discussed in Sec. III above, the results in Fig. 5(d,e) suggest that for certain values of $\mu$, the DDFT-5 equation (67) has solutions for the density $1+n$ which go to zero at certain places. In contrast, the density in the PFC models appears to stay away from zero (although there is no reason for it to do so and there would be no singularity if it did), and in DDFT-3, the density minimum gets smaller and smaller as $\mu$ increases, but remains positive, without the sharp cutoff seen in DDFT-5. In this section we argue that the density reaching zero is not an artefact of numerical difficulties, rather it is a feature of the DDFT-5 equation (67). Here, we focus on singularities in the solution, not on stability. Our discussion in this section is mainly framed in terms of the stripe solution, which is unstable, but the stable up-hexagon branch has similar issues, as is illustrated in Fig. 5(e).

Figure 8(a) shows that even close to the end of the branch of DDFT-5 stripes, the density profile $1+n$ remains smooth, but since its minimum at $x=x_{\rm min}$ is very close to zero ($1+n(x_{\rm min})\approx 5\times 10^{-7}$), therefore the logarithm $\log(1+n)$ is sharply spiked towards large negative values at $x_{\rm min}$. Writing out the terms in (67) for a density profile only varying in the $x$-direction, we have:

[TABLE]

which suggests that the only way to balance a large negative contribution from $\log(1+n)$ is to have a large positive $\gamma{n_{xxxx}}$. Figure 8(b) shows that these two terms (solid lines with markers at the grid points) do indeed balance each other. The figure also shows that the other terms in (77) are well behaved and that the equation is satisfied at each grid point. Therefore this singularity is not a numerical artefact, but rather a genuine feature of DDFT-5 stripes: the minimum density goes to zero at a certain finite value of $\mu$. In Fig. 8(c,d) we show for comparison results from DDFT-3: the stripes have sharply peaked density maxima and density minima just as small as in DDFT-5, but all terms in Eq. (63) (Fig. 8d) are well behaved.

As $\mu$ is further increased, the minimum of the density in DDFT-5 gets closer to zero, so the logarithm of the density goes further towards $-\infty$ and correspondingly $\gamma{n_{xxxx}}$ goes towards $+\infty$. Fig. 5(d) shows that the density minimum gets to zero at a finite value of $\mu$. We have not been able to develop a consistent asymptotic approximation for this limit in the DDFT-5 equation. However, to illustrate that apparently smooth solutions with logarithmic singularities in their fourth derivatives can easily be found, consider for example taking $\gamma=1$ in (77) and taking a density profile that has a quadratic minimum at $x=x_{\rm min}=0$:

[TABLE]

where $A$, $B$ and $C$ are constants. For small $x$, the largest of these three terms is $Ax^{2}$, so $\log(1+n)\approx\log(Ax^{2})$, which goes to $-\infty$ as $x\rightarrow 0$. The other terms are $n_{xx}\approx 2A+{\cal O}(x^{2},x^{2}\log(x^{2}))$ and $n_{xxxx}\approx 24C\log(x^{2})+24B+100C$. Adding these three together requires $1+24C=0$ to cancel the logarithmic singularity at $x=0$, and the remaining terms are constants or go to zero as $x\rightarrow 0$.

As can be seen from Fig. 8, having an adequate resolution for our numerical calculations was a challenge but for different reasons for the different models. In DDFT-3, the density maxima can be sharply peaked while the logarithm of the density is smooth, so inadequate resolution in the density field prompts an increase in the number of grid points (we implement automatic regridding, as discussed in Appendix B). In contrast, in DDFT-5, the density field can be smooth but with minima very close to zero, so its logarithm has very sharp negative peaks. In this case, inadequate resolution in the logarithm of the density prompts regridding. The difference is that in DDFT-5, the equation involves derivatives so any problem is magnified, while in DDFT-3, the equation involves convolutions that smooth out any problems.

These arguments indicate that a singular solution to the DDFT-5 equation (67) of the type seen in Fig. 8 is possible, with the density going to zero, and that this is not a problem of inadequate numerical resolution, but rather a consequence of replacing the convolution in DDFT-3 with derivatives. A full asymptotic theory should result in a prediction for the value of $\mu$ at which the branch terminates. The stripe solutions of the DDFT-3 equation (63) can also have small density but without any singularity in the solution.

V One mode approximation for DDFT

The data displayed in Fig. 8(a,c) lead to an interesting observation: in (a) DDFT-5, the logarithm of the density is sharply negatively peaked, while the density is smooth (at least up to its second derivative), slowly varying and resembles a cosine. In contrast, in (c) DDFT-3, the density is sharply peaked, while the logarithm of the density is slowly varying and resembles a cosine. One of the attractions of PFC theory is that it has slowly varying solutions that are well represented by a few Fourier modes Elder et al. (2002); Elder and Grant (2004); Emmerich et al. (2012); Wu et al. (2010), and this carries over to some extent to DDFT-5. Such Fourier representations of the density profiles in DDFT-3 are unsatisfactory, apart from for the unstable solutions very close to the spinodal point, since any solution of reasonable amplitude is sharply peaked.

However, the data in Fig. 8 suggest that representing the logarithm of the density with a few Fourier modes should work well in DDFT-3. In this section, we elaborate how such a theory can be developed.

The key is to write Eq. (63) in terms of $\phi({\bm{x}})\equiv\log(1+n({\bm{x}}))$, and approximate $\phi$ by a few Fourier modes. As a first step, we write out ${\cal L}$ using Eq. (65) and re-write (63) as:

[TABLE]

where $\psi(|{\bm{x}}-{\bm{x}}_{2}|)=e^{-|{\bm{x}}-{\bm{x}}_{2}|^{4}/R^{4}}$. The convolution term (including the $\rho_{0}\beta\epsilon$ prefactor) in this equation is $-n({\bm{x}})-{\cal L}n({\bm{x}})$. We know the eigenvalues of ${\cal L}$: ${\cal L}\exp(ikx)=\sigma(k)\exp(ikx)$, which means that the convolution term acting on a Fourier mode $\exp(ikx)$ has eigenvalue $-(1+\sigma(k))$. We also know that for high wavenumbers, the convolution averages to zero, and indeed $\sigma(k)\rightarrow-1$, as can be seen in Fig. 2.

We focus first on stripes, which have Fourier components only at integer wavenumbers, and notice that $1+\sigma(2)$ is already very small (less than $0.01$), because $\hat{\psi}(2)$ is small. This implies that the Fourier components of the convolution term at $k=2$ and higher will be much smaller than the Fourier components at $k=0$ and $k=1$ – regardless of the spectrum of $n$ itself. The other two terms in Eq. (79) are $\mu$, which is constant ($k=0$ only), and $\log(1+n)$, so $\log(1+n)$ can only have significant Fourier components at $k=0$ and $k=1$: there is nothing to balance modes with $|k|\geq 2$. This explains why $\log(1+n)$ in the lower left panel of Fig. 8 is smooth, and it is also why approximating the logarithm of the density by a few Fourier modes must work regardless of the amplitude of the modulations in the density, or how sharply they are spiked, or of the value of $\mu$.

For the stripe phase we write

[TABLE]

where $\phi_{0}$ and $\phi_{1}$ are constants (real and complex, respectively) that we need to find. This can easily be generalised for hexagons and other periodic phases by adding more modes in (80). The $k=0$ and $k=1$ components of $\exp\left(\phi_{1}e^{ix}+\text{c.c.}\right)$, where c.c. denotes the complex conjugate, can be expressed in terms of integrals, defining two functions $f_{0}(\phi_{1})$ and $f_{1}(\phi_{1})$ given by

[TABLE]

i.e., $f_{0}$ is modified Bessel function of the first kind of order zero and $f_{1}$ is a Fourier transform generalisation of $f_{0}$. Using these functions, $n(x)=e^{\phi(x)}-1$ can be written in terms of its Fourier components as

[TABLE]

The modes with $|k|\geq 2$ in (82) are large in amplitude, but they are reduced in significance in Eq. (79) by the convolution, as explained above. The action of the convolution on modes with $|k|<2$ can be written in terms of $\sigma(k)$. Retaining only these terms, we are left with the $k=0$ and $k=1$ components of (79):

[TABLE]

Notice that the only information remaining from the GEM-4 potential is the values of $\sigma(0)$ and $\sigma(1)$, i.e., the values of $\hat{u}(0)$ and $\hat{u}(1)$. Recall also that if the reference density $\rho_{0}$ is chosen to be the value at the spinodal, then we have $\sigma(1)=0$. These equations can also be written in terms of the pair potential as

[TABLE]

The two equations in (83) can easily be solved for $\phi_{0}$ and $\mu$ in terms of $\phi_{1}$, from which the density can be reconstructed. The agreement between this and the full solutions of DDFT-3 is astonishing. Figure. 9(a) shows a 1D example at $\mu=28.9$, with the full solution as a black line and the approximate solution as a dashed line. Even though the density varies by two orders of magnitude, the two are almost indistinguishable. There is similar excellent agreement with the branch of stripe solutions (Fig. 9b) in 2D DDFT-3 (recall that the GEM-4 potential has different values of $\sigma(0)$ in 1D and 2D).

For 2D hexagons, the approach is similar, with the two $e^{\pm{ix}}$ terms in Eq. (81) replaced by six similar terms, with wavenumbers ${\bm{k}}$ that are uniformly spaced around a circle of radius 1 in ${\bm{k}}$-space (c.f. Eq. (102) in Appendix B). The agreement in this case (Fig. 9c) is also very good. If one adds a further six modes with $|{\bm{k}}|=\sqrt{3}$, then the agreement is as good as that for stripes. These approximate solutions can easily be continued up to $\mu=100$ without difficulties, where we observe very sharply peaked density maxima and extremely small but non-zero values for the density minimum.

VI Discussion and conclusions

In this paper, starting from DDFT, we have presented a step-by-step derivation of PFC theory, at each stage explaining the consequences of the approximations. The approximations can be listed under three main groupings: (i) making a truncated functional Taylor expansion approximation for the excess Helmholtz free energy, and then making the RPA/RY approximation. This leads to DDFT-3 in our classification. (ii) Neglecting the $\nabla\cdot[n\nabla{\cal L}n]$ term, which effectively also forces making a Taylor expansion of the logarithmic ideal gas term and assuming constant mobility. (iii) Replacing the nonlocal convolution in ${\cal L}$ with a local gradient expansion. The consequence of (ii) is to introduce a second spinodal into the phase diagram and to significantly alter the relative stabilities of the different periodic states, in particular making striped states to become an equilibrium phase for some state points, which is contrary to the physics. The consequence of making (iii) without first making (ii), is to generate a theory (DDFT-5) that has a no-solution region in the phase diagram, such as that displayed for the GEM-4 model in Fig. 7(b). All these consequences have been illustrated for the GEM-4 system, chosen because DDFT-3 is fairly accurate for this model for temperatures $k_{B}T/\epsilon>0.1$, allowing us to see the influence of the subsequent approximations.

Throughout, there is good quantitative agreement between DDFT-3, PFC-$\gamma$ and the EOF versions of DDFT-5 and PFC-$\epsilon$ (data not shown) only for unstable small amplitude solutions close to the spinodal point. The region of quantitative agreement agreement between DDFT-3 and PFC-$\gamma$ is circled in red in Fig. 5(b,h). Beyond this region, the agreement between the four theories is at best qualitative.

Given all these problematic consequences for PFC theory, especially the issues related to Taylor expanding the logarithmic ideal gas term, it raises the question of why then is PFC theory so successful? In our view, there are several reasons for this. The first reason is that PFC theory is qualitatively correct near to the spinodal. Therefore, it can satisfactorily describe the coexistence between the liquid and the crystal phase, which is often an important aspect in applications of the theory. Second, despite the approximations, PFC theory still incorporates some very important physics: (i) the free energy satisfies the correct symmetries, (ii) the dynamical equation gives a time evolution that decreases the free energy monotonically over time and (iii) the current is proportional to the gradient of the chemical potential. These are all important features for describing many phenomena. Also, many of the features that PFC theory is used to describe are generic, and the model parameters can be scaled to fit (for example) iron Jaatinen et al. (2009) and graphene Fan et al. (2017), but could equally well be scaled to match other materials, with similar good agreement. This universality (having the correct symmetries etc.) underlines the importance of PFC theory as a powerful model of generic features of crystalization, but it means that PFC theory will in general not be able to predict any unusual (non-generic) behaviour.

Our results in Sec. V for the GEM-4 system showing that one can derive a very accurate one-mode amplitude equation approximation for the field $\phi({\bm{x}})=\log(\rho({\bm{x}})/\rho_{0})$, rather than the density itself, gives a tantalising hint as to how PFC-type theories may more properly be derived and what the order parameter field in PFC theory really represents: Should we consider the PFC order parameter to be a scaled logarithm of the density distribution or some other similar function of the density, rather than being proportional to the density profile itself? On the basis of the work presented here, the answer to this question is ‘probably yes’, but clearly more work is required to fully address this.

Returning to theories for the density profile, in our view it is preferable to retain the logarithmic ideal gas term in the approximation used for the Helmholtz free energy functional, since this is required to have physical (i.e., positive) density profiles, and also because with this the DDFT dynamics leads to the (correct physics) linear diffusion equation term – see Eq. (13). The difficult consequence of retaining this in the free energy is that one then has to deal with terms in the dynamical equation of the form $\nabla\cdot[n\nabla{\cal L}n]$, which makes solving numerically far more difficult. However, this term also contributes to making the crystalline phase more stable than the stripe phase, which also makes it important. The crucial contribution of this term can especially be seen in DDFT-3, since in this version of the theory it is clearly required for stabilizing the crystal structure. In general, we would advocate using one of DDFT-1, DDFT-2 or DDFT-3 (depending on how ${\cal F_{\text{ex}}}$ is treated) over all existing PFC theories, for studying the properties of real materials.

It is worth noting that in all the approximations made here, we only consider those that retain the form of a dynamics that decreases the free energy monotonically over time – see Eq. (6). This is a feature of both DDFT and PFC theory. In our view this structure is important and should not be broken by any approximations made, i.e., any that might be made in future in attempting to avoid any of the above mentioned issues.

It is also worth noting that, while we have not discussed the consequences of the LDA in going from DDFT-1 to DDFT-2 (see Eq. (31)), the polynomial terms in the chemical potential (36) can potentially lead to the same problem of having a second spinodal, even while retaining the logarithm term. This applies to DDFT-4 as well.

In this paper we have largely focussed on making our arguments in two dimensions, in order to keep the presentation as simple and comprehensible as possible. However, we should emphasis that all of our arguments apply for three dimensional (3D) systems. For example, at the higher temperatures we have focussed on here, the 3D GEM-4 model exhibits at equilibrium a single fluid phase and two crystalline phases: the body-centered cubic phase at lower densities and the face-centered cubic structure at higher densities. These are all accurately predicted by the 3D version of DDFT-3 Mladek et al. (2006). There are no columnar or lamellar phases, which are the 3D equivalents of the stripe phase. On the other hand, the 3D versions of the PFC theories presented here all predict a lamellar phase at some state points Thiele et al. (2013), the 3D generalisations of the down hexagons and a second spinodal with the uniform liquid becoming the equilibrium phase at higher densities.

It would be interesting to explore how the dynamics of defects, the elasticity and the plasticity of crystals differs between DDFT and PFC, and our results are also relevant to binary systems. In the derivation in Ref. Huang et al. (2010) of a PFC theory for binary mixtures, the generalisation of the $\nabla\cdot[n\nabla{\cal L}n]$ term is retained until the last moment in the derivation, but then dropped for the same reasons that is is neglected for one-component systems. Given the importance of this term for one-component systems, it is surely also important for stabilizing crystal structures in binary systems. Note also that when determining mechanical properties such as elastic constants, the terms in the free energy that are linear in $n$ can be important Wang et al. (2018). These have been neglected here throughout since such terms do not contribute to determining density profiles.

The singularity observed for DDFT-5 as the chemical potential $\mu$ (or equivalently, the average density $1+\bar{n}$) is increased was found by continuing equilibrium solutions determined at lower values of $\mu$. One aspect that needs further investigation relates to determining the influence of this when DDFT-5 is solved for state points where the final equilibrium crystal (or stripe) solution for the density profile does not exist. For example, a situation we have in mind is that studied by van Teeffelen et al. van Teeffelen et al. (2009) consisting of a solidification front propagating into the unstable liquid. These authors compared results for this situation between (in our terminology) DDFT-3, DDFT-5 and PFC-$\epsilon$. Their results are for two-dimensional dipolar colloidal particles. From the DDFT-5 results (PFC1 in their terminology) displayed in Figs. 4 and 5 of Ref. van Teeffelen et al. (2009), it can be seen that they did not consider values of the coupling parameter large enough to encounter any of the singularities; the density profiles stay well away from zero. It would be interesting to quench deeper into the crystal phase to study the evolution of the density distribution towards the singular state. However, the numerics to resolve this accurately would surely be difficult.

One aspect of PFC theory that the derivation from DDFT highlights is that in general one is not free to independently vary the parameters $r$ and $\bar{n}$ in Eqs. (1) and (2). For example, for the GEM-4 model there are certain values of $r$ that are not generated by the mapping from DDFT. Of course, by changing the pair potential, different combinations of the PFC model parameters can become accessible. We should also recall that although we have illustrated many of our conclusions by considering the soft-core GEM-4 model, PFC theory can be derived for systems of particles with hard potentials since it is the pair direct correlation function $c^{(2)}({\bm{x}}_{1},{\bm{x}}_{2})$ that enters the theory; this quantity is finite for all values of ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$.

As a final point, we mention that our results will also be of interest to the pure mathematics community. DDFT-3 is also referred to as the McKean–Vlasov equation and in this context there are a number of recent interesting rigorous results Carrillo et al. (2018); Gomes et al. (2019). Our results for DDFT-5, showing that for a finite value of $\mu$ there is a singularity with the density profile going to zero, may well be of interest to those who study the mathematics of solutions to partial differential equations with compact support – see for example Ref. Bernis et al. (1992). For values of $\mu$ beyond the singular point where $1+n({\bm{x}})\to 0$ it may be that the solutions become complex. If one were interested to find these solutions, we believe it might require treating $\mu$ as a complex variable. Of course, all of this is out of the realm where the model represents a theory for matter.

Acknowledgements.

This work was supported in part by a L’Oréal UK and Ireland Fellowship for Women in Science (PS), by the EPSRC under grants EP/P015689/1 (AJA, DR) and EP/P015611/1 (AMR), and by the Leverhulme Trust (RF-2018-449/9, AMR). We are grateful for conversations with Tapio Ala-Nissilä, Daniele Avitabile, Ken Elder, Zhi-Feng Huang, Kai Jiang, Edgar Knobloch, Ron Lifshitz, Chris Malcotte, Daniel Read, Uwe Thiele, Steve Tobias, Gyula Tóth, Laurette Tuckerman and Joanna Tumelty. We are grateful also for constructive comments from two anonymous referees.

Appendix A Linear theory for GEM-4

In this appendix, we discuss how we compute the linear theory for the GEM-4 potential in (64). To be specific, in a two-dimensional periodic domain, the eigenvalue $\sigma(k)$ is defined by ${\cal L}e^{i{\bm{k}}\cdot{\bm{x}}}=\sigma(k)e^{i{\bm{k}}\cdot{\bm{x}}}$ and (65), which can be written as

[TABLE]

where the integral is taken over the periodic domain (the GEM-4 potential is replaced by its periodic extension). We set ${\bm{k}}=(k,0)$, we integrate from $[-N\pi/k,N\pi/k]$ in each dimension, and we choose the integer $N$ large enough that the GEM-4 exponential is effectively zero at the boundaries; $N=4$ suffices. We then scale ${\bm{x}}$ by a factor of $k$, replacing ${\bm{x}}_{2}-{\bm{x}}$ by ${\bm{x}}/k$, so that the integral becomes:

[TABLE]

where $x$ is the first component of ${\bm{x}}$. With this scaling, the limits of the integral do not depend on $k$.

We choose $R$ and $\rho_{0}\beta\epsilon$ so that $\sigma(1)=0$ and $\frac{{\rm d}\sigma}{{\rm d}{k}}(1)=0$. The derivative of $\sigma$ with respect to $k$ is

[TABLE]

Evaluating this at $k=1$ and removing the constant factor outside the integral gives a function $F(R)$:

[TABLE]

We solve the equation $F(R)=0$ using Newton’s method to give $R$. We then calculate $\rho_{0}\beta\epsilon$ by requiring that $\sigma(1)=0$ in (86). Values for $R$ and $\rho_{0}\beta\epsilon$ in one, two and three dimensions are given in table 2. We compute the GEM-4 dispersion relation in Fig. 2 in a similar way, and use the second derivative of (86) with respect to $k$, at $k=1$, to find the $\gamma$ parameter (also given in table 2) in ${{\cal L}_{\text{grad}}}$. The table also gives $\sigma(0)$, since this is useful for computing properties of the liquid, as well as $E_{B}$, the coefficient in the eighth-order model of Jaatinen et al. (2009), given in Eq. (76).

Appendix B Numerical method: continuation

We use numerical continuation to solve the four equations (63) and (67)–(69) for $n({\bm{x}})$ as the parameters vary. Our approach is based on Doedel et al. (1991) for the pseudo-arclength continuation method, and we use the approach advocated by Kelley (2003) to solve the large linear systems at each Newton step. The main parameters are the chemical potential $\mu$, the parameters in the linear operators ${\cal L}$ and ${{\cal L}_{\text{grad}}}$, and the domain size. In this discussion, we focus on $\mu$ as the parameter that is varied.

B.1 Pseudo-arclength continuation

The main idea behind pseudo-arclength continuation is to suppose that we are looking to calculate a branch of solutions $n({\bm{x}})$ depending on the parameter $\mu$. The branch may have folds (as in Fig. 5), and the method should be able to go around these. The method defines a parameter (the arclength $s$) that increases or decreases monotonically along the branch (including its folds) such that both $n({\bm{x}})$ and $\mu$ can be regarded as single-valued functions of $s$. Then the equation to be solved is

[TABLE]

where ${\cal G}$ represents the equation we are solving for $n$. Instead of thinking of $n$ as being a function of position ${\bm{x}}$, we represent $n$ as a series of values $n_{i}$ on $N$ grid points ${\bm{x}}_{i}$ ($i=1,\dots,N$), so $n$ is now a vector in $\mathbb{R}^{N}$, and ${\cal G}$ is a function from $\mathbb{R}^{N+1}\rightarrow\mathbb{R}^{N}$. Equation (89) represents $N$ equations for $N+1$ unknowns, and so it is supplemented by an orthogonality condition, that the next point on a branch should lie in a plane orthogonal to a line connecting the two previous points. It is this that allows the branch following technique to go around folds. If we have two points on the branch $(n(s),\mu(s))$ at $s_{0}$ and $s_{1}$, then we take the derivatives of $n$ and $\mu$ with respect to the arclength to be approximately

[TABLE]

with the scaling factor $S$ chosen so as to satisfy

[TABLE]

The $\frac{1}{N}$ prefactor means that the parameterization of the branch by the arclength is essentially independent of the number of grid points.

The method then proceeds in a predictor–corrector fashion. The predictor step, with a target stepsize $\Delta{s}$ provides $(n_{2},\mu_{2})$:

[TABLE]

Then, $(n_{2},\mu_{2})$ is used as an initial iterate for a Newton solver for equation (89), supplemented by the condition that the Newton iterates lie in a plane orthogonal to the line given in (92), parameterised by $\Delta{s}$. This means that we are solving ${\cal H}(n,\mu)=0$, where ${\cal H}$ is $\mathbb{R}^{N+1}\rightarrow\mathbb{R}^{N+1}$, with the first $N$ equations in ${\cal H}$ being the same as ${\cal G}$, and the last equation being

[TABLE]

To be precise, we take

[TABLE]

where $P$ is a linear preconditioner for ${\cal G}$ (see below). The $\sqrt{N}$ scaling means that the norm $\lVert{\cal H}(n,\mu)\rVert$ (the square root of the sum of the squares of its components) is independent of the number of grid points $N$, and it also means that the equations in ${\cal G}$ and the orthogonality condition (93) are given a similar weighting by the Newton method.

Solving ${\cal H}(n,\mu)=0$ results in a new point on the branch of solutions, $(n(s_{2}),\mu(s_{2}))$, where $s_{2}$ is given by

[TABLE]

with $\frac{{\rm d}{n}}{{\rm d}{s}}$ and $\frac{{\rm d}\mu}{{\rm d}{s}}$ given by (90). This last equation comes from replacing $n_{2}$ by $n(s_{2})$ and $\mu_{2}$ by $\mu(s_{2})$ in (92) and finding a $\Delta{s}=s_{2}-s_{1}$ from $1/N$ times the first equation plus the second equation. This is not quite the actual change to the arclength that was achieved in the step, and the approximation is the reason that the method is called the pseudo-arclength method.

B.2 Newton’s method

For Newton’s method, we define $X=(n,\mu)$ and solve ${\cal H}(X)=0$. We start with $X_{0}$ given by the predictor step above in (92), and follow Kelley (2003), at each step solving the linear equation

[TABLE]

where $\frac{\partial{\cal H}}{\partial{X}}$ is the $(N+1)\times(N+1)$ matrix of derivatives of ${\cal H}$ with respect to $X$, and then improving our estimate of the root by using ${\delta{\mkern-1.0mu}X}$.

The Newton method proceeds until convergence, defined by

[TABLE]

where $N_{\text{abs}}$ and $N_{\text{rel}}$ are the Newton absolute and relative convergence tolerances respectively, typically $10^{-10}$ and $10^{-8}$. We also monitored the maximum of $|{\cal H}(n({\bm{x}}),\mu)|$ across the domain, and this was typically no larger than ten times $N_{\text{abs}}$, so the equations are well satisfied at each point in space as well as in norm.

The linear equations in (96) are solved to find ${\delta{\mkern-1.0mu}X}_{n}$ using Matlab’s biconjugate gradient stabilized (l) (bicgstabl) method. This allows the matrix–vector multiplications to be evaluated using a function (rather than by explicitly computing a large matrix). The method is iterative, and proceeds until

[TABLE]

where the relative tolerance $L_{\text{rel}}$ of the linear solver is chosen so as to balance the number of Newton steps against the number of bicgstabl iterations. Based on Kelley (2003), we choose

[TABLE]

subject to the constraint that $L_{\text{rel}}$ should be no larger than $0.1$. The effect of this is that in the initial first or second of the Newton iterations, when $\lVert{\cal H}(X_{n})\rVert$ is at its largest, the linear solver is not asked to work too hard to solve (96), since any reasonably good approximate solution is likely to improve the estimate of the root, and an absolutely perfect solution isn’t going to do much better. In the middle stages of the Newton iterations, when $\lVert{\cal H}(X_{n})\rVert$ is about $10^{-6}$, the tolerance $L_{\text{rel}}$ is about $2\times 10^{-4}$, which is not good enough for quadratic convergence of Newton’s method, but is good enough to provide two or three orders of magnitude improvement to the quality of the solution at a considerably lower cost. In the final stages of the Newton iterations, $L_{\text{rel}}=0.1$, good enough for polishing the solution to the tolerance $N_{\text{abs}}$ while not attempting to solve the linear problem down to round-off error.

Also based on Kelley (2003), we implement the Armijo rule, which ensures that each Newton step results in an improvement to the solution. The idea is that the solution of the linear equation (96) should give the correct direction for improving the solution of ${\cal H}(X)=0$, but taking a full step may not actually result in an improvement, so instead we set

[TABLE]

where $j=0,1,2,\dots$ is chosen to be the smallest such that

[TABLE]

In most cases, the first ($j=0$) Armijo step satisfies (101), equivalent to the normal Newton method, but when the density is close to zero in the DDFT calculations, and small changes in density lead to large changes in its logarithm, the Armijo rule is helpful.

We do not use a preconditioner in the GEM-4 calculations (so $P$ in (94) is the identity), but in the gradient expansion calculations, a preconditioner is helpful. In Fourier space, ${{\cal L}_{\text{grad}}}$ can easily be inverted, so the preconditioner is ${{\cal L}_{\text{grad}}^{-1}}$ when the absolute value of the eigenvalue $\sigma(k)$ is greater than $1$, otherwise the preconditioner is the identity. This has the effect of reducing the number of iterations needed to solve (96) by a factor of 10 or even 100.

A sample Matlab code to solve Eq. (63) for DDFT-3 by Newton’s method (without the continuation aspect) is given in the supplementary material.

B.3 Additional considerations

We start the computation of each branch close to $n=0$ and $\mu=0$ using an approximate solution derived from weakly nonlinear theory. For example, for hexagons in DDFT-3, we take

[TABLE]

where the initial value of $\mu$ is small, $\sigma(0)$ comes from Table 2, ${\bm{k}}_{1}=(1,0)$, ${\bm{k}}_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})$, ${\bm{k}}_{3}=(-\frac{1}{2},-\frac{\sqrt{3}}{2})$, and c.c. stands for complex conjugate.

The equations are posed on periodic domains and we use $N_{x}\times{N_{y}}$ grid points, depending on the number of wavelengths in the domain and the nature of the solution. The GEM-4 hexagonal calculations require $8\times\frac{8}{\sqrt{3}}$ wavelength domains, with resolutions starting at $80\times 48$ Fourier modes close to onset. At larger amplitude, $512\times 320$ Fourier modes (or even more) are needed, especially if the density maxima are sharply peaked (in DDFT-3) or if the density minima are very close to zero (in DDFT-5). In order to accommodate the changing needs for resolution along a branch, we monitor whether the solutions are well resolved and implemented automatic regridding, so as to maintain enough grid points to resolve the solution well, regardless of what features emerge as $\mu$ is varied. Typically we require that the amplitudes of the highest-wavenumber Fourier modes be no higher than $10^{-10}$ times the largest Fourier amplitude.

We also implement automatic pseudo-arclength stepsize control: $\Delta{s}$ is increased by a factor of 1.1 (up to a maximum of 0.1) if the Newton method converges quickly (in fewer than 5 iterations), or is decreased by a factor of 2 (down to a minimum of $10^{-6}$) if it converges slowly (more than 8 iterations) or fails.

Finally, we adjust the domain size continuously along each branch so as to minimise the specific grand potential $\Omega/A$. This is done by adding an extra parameter (the domain stretch factor $K$), so that the real domain size is $KL_{x}\times{K}L_{y}$ instead of an unstretched $L_{x}\times{L_{y}}$. The number of grid points is not altered. Then the real GEM-4 potential is proportional to $\psi(|{\bm{x}}|)=\exp\left(-K^{4}|{\bm{x}}|^{4}/R^{4}\right)$, where ${\bm{x}}$ is the unstretched coordinate on the unaltered grid. Quantities like the mean value of $n$ are the sum of the values of $n$ at each grid point divided by the number of grid points, and so are not affected by the stretch factor. The only parts of $\Omega/A$ that are affected are the convolutions (for ${\cal L}$ and the GEM-4 potential) and the spatial derivatives (for ${{\cal L}_{\text{grad}}}$). In the case of the GEM-4 potential, when considering the specific grand potential $\Omega_{3}/A$ arising from (40) for example, and ${\cal L}$ given by (65), the integrals in $\Omega_{3}/A$ are proportional to $K^{2}$, with additional $K$ dependence coming from the GEM-4 potential itself. Therefore,

[TABLE]

where $\otimes$ represents the convolution integral evaluated on the unstretched grid, and the ${\bm{x}}$ integral is also on the unstretched grid. In the case of ${{\cal L}_{\text{grad}}}$, Laplacians on the real grid are a factor of $K^{-2}$ times Laplacians on the unstretched grid, so ${\rm d}(\Omega/A)/{\rm d}{K}$ is evaluated accordingly. In both cases, an extra equation (${\rm d}(\Omega/A)/{\rm d}{K}=0$) is added to ${\cal H}$ in (94), and this is solved alongside all the other equations. The Jacobian also needs to be evaluated. In practice this made little difference in the cases considered here, and the domain stretch factor largely stayed between $0.98$ and $1.02$.