# Vacancy diffusion in colloidal crystals as determined by dynamical density functional theory and the phase field crystal model

###### Abstract

A two-dimensional crystal of repulsive dipolar particles is studied in the vicinity of its melting transition by using Brownian dynamics computer simulation, dynamical density functional theory and phase-field crystal modelling. A vacancy is created by taking out a particle from an equilibrated crystal and the relaxation dynamics of the vacancy is followed by monitoring the time-dependent one-particle density. We find that the vacancy is quickly filled up by diffusive hopping of neighbouring particles towards the vacancy center. We examine the temperature dependence of the diffusion constant and find that it decreases with decreasing temperature in the simulations. This trend is reproduced by the dynamical density functional theory. Conversely, the phase field crystal calculations predict the opposite trend. Therefore, the phase-field model needs a temperature-dependent expression for the mobility to predict trends correctly.

###### pacs:

82.70.Dd, 64.70.D-, 81.10.-h, 81.10.Aj## I Introduction

Most of the mechanical properties of crystals depend crucially on the presence of crystalline defects. For material processing it is therefore of principal importance to understand and control the defect concentration and dynamics. The nature and dynamics of defects is much easier to classify for crystalline sheets in two spatial dimensions. In this case, it is known for a long time, that the formation and unbinding of topological defects provides an efficient way of melting according to the two-stage scenario proposed by Kosterlitz-Thouless-Nelson-Halperin-Young (KTNHY) Strandburg:88 (). Defects can also be accumulated near edges of crystalline sheets and do occur for two-dimensional crystals on more complicated manifolds Nelson:10 (); Irvine:10 ().

Defects are highly dynamic: Whereas the structure of a crystal is static over long time scales, defects undergo diffusion in the crystalline background. The diffusive dynamics of individual point defects were observed directly in two-dimensional colloidal suspensions of charged micro-spheres by video microscopy Pertsinidis:01 (); Pertsinidis:05 () and they were confirmed and further analyzed by computer simulations Libal:07 (). The dynamics of defects were also explored by real-space methods in a two-dimensional crystal of weakly damped dust particles in a plasma by Nosenko and coworkers Nosenko:11 (); Nosenko:09 (); Nosenko:07 (), see also Book ().

Describing the defect dynamics by a microscopic theory is still a formidable challenge, in particular close to melting. Such a theory should contain the solid and fluid phase and give a reliable picture on the defect concentration and its dynamics. Recent progress in this respect has been made by classical density functional approaches of freezing Evans:79 (); Singh:91 (); Loewen:94 (); EmmerichLowen2012AdvPhy (). For Brownian particles, the density functional approach can be generalized towards dynamics Marconi:99 (); Archer:04 (); EspanolL2009 () and the dynamics of solidification has been examined in two dimensions SvenPRL (); SvenRainerPFC (). Similar in line, a more coarse-grained phase-field-crystal model has been proposed to describe crystal growth ElderKHG2002 (); ElderG2004 (); ElderPBSG2007 (); Emmerich2009 () and the defect structure and dynamics, for various applications see e.g. YuLV2009 (); JaatinenAEAN2010 (); AchimRKEGNY2009 (); RamosAEYN2009 (); AchimKGYN2006 (); JaatinenAEAN2009 (); TegzeGTPJANP2009 (); HuangE2008 (); MellenthinKP2008 (); McKennaGV2009 (); WuV2009 (). However a systematic exploration of defect dynamics by such a density functional theory has not yet been performed nor has the reliability of the phase-field-crystal model been systematically checked as far as the trends of defect dynamics are concerned.

Here, we study the dynamics of vacancies in a two-dimensional colloidal crystal by using Brownian dynamics computer simulations, dynamical density functional theory and the phase-field crystal approach and thereby test the ability of the theoretical approaches to qualitatively reproduce the observations made in the simulations. The model system we use here is a two-dimensional suspension of dipolar colloids. This system has been realized experimentally as superparamagnetic particles at an air-water interface Maret:09 (). When exposing superparamagnetic particles to an external magnetic field perpendicular to the plane, their induced magnetic dipole moment leads to an effective repulsive interaction whose amplitude can be tuned by the magnetic field strength. At sufficiently high field strength, the system crystallizes into a two-dimensional triangular (i.e. hexagonal) crystal. This system has been studied extensively by computer simulations and by the aforementioned dynamical density functional theory SvenPRL () and phase field crystal models SvenRainerPFC ().

Out of a perfect triangular crystal, some particles are removed and the relaxation of the resulting defect and its mobility are extracted by monitoring the one-particle density as a function of time. We confirm by simulation that the defect mobility is increasing with increasing temperature, as was already observed for charged particles by Líbal et al Libal:07 (). This behavior is qualitatively and semi-quantitatively reproduced in dynamical density functional theory based on a static Ramakrishnan-Yussouff-like density functional Rogers:84 (); Teeffelen:06 (). The phase-field-crystal model, on the other hand, fails to predict the trend of the temperature-dependence of the mobilities. This is mainly attributed to the constant kinetic prefactor involved in the phase-field-crystal approach. For predicting the trend, a temperature-dependent corrective mobility is needed for the phase-field-crystal model.

This paper is organized as follows: in section II, we briefly propose the different approaches used in this paper. Results are presented in section III and we conclude in section IV.

## Ii Theoretical models

Dynamical density functional theory and the more coarse-grained phase-field-crystal model describe the overdamped Brownian dynamics in terms of a continuity equation for the deterministic, time-dependent, and ensemble averaged one-particle density . Note that in many applications of the phase-field-crystal model is interpreted as a fluctuating density field that changes for different realizations of the dynamical evolution even under the same initial conditions. Here, we will not take this approach but regard as a purely deterministic quantity. For a more thorough discussion of the physical interpretation of in the phase-field-crystal model see SvenRainerPFC (); Rauscher ().

The temporal evolution of the density field according to dynamical density functional theory Marconi:99 (); Marconi:00 () is given by

(1) |

with the single-particle diffusion constant and the thermal energy. The Helmholtz free energy functional is provided by classical density functional theory Evans:79 (). In the crystal, the driving current obviously assumes the hexagonal symmetry of the underlying crystal. In the interstitial regions, this current is small since the density itself is small.

Note that Eq. (1) can be derived from first principles Marconi:99 (); Marconi:00 (), i.e., from the microscopic Langevin equations of motion or from the Smoluchowski equation for the time evolution of their respective probability distribution (for a review see SvenRainerPFC ()). Here, we employ the same approximation to the density functional theory of Ramakrishnan and Yussouff Ramakrishnan:79 () already introduced in Refs. SvenPRL (); SvenRainerPFC (). Eq. (1) then reads

(2) |

where is the time-dependent external potential. and are the excess free energy and the direct correlation function of the reference fluid of density , respectively.

In this work we consider a two-dimensional system of magnetic dipoles that are oriented perpendicular to the 2D plane. The pair potential of two dipoles in the plane is given by

(3) |

where is the interaction strength. The thermodynamics and structure depend only on one dimensionless coupling parameter , where is the average one-particle density and is the thermal energy.

The two-particle direct correlation function of the fluid hansen-mcdonald:86 () has been obtained for a large range of coupling constants from liquid-state integral equation theory as previously described SvenJPCM (); Teeffelen:06 (); SvenRainerPFC ().

In order to measure the diffusion of defects, Eq. (1) is numerically solved on a rectangular periodic box of a fine grid with grid points per nearest-neighbor distance . A finite difference method with variable time step is applied. The convolution integrals are solved using the method of fast Fourier transform.

For the more coarse-grained phase-field-crystal model we employ the two different versions termed PFC1 and PFC2 in Ref. SvenRainerPFC (). The PFC1 model constitutes an approximation of dynamical density functional theory, as introduced above. The last term in Eq. (2) is replaced by its gradient expansion. The dynamical equation then reads

(4) |

The parameters , , and are the fit parameters of a parabola to the first maximum of the Fourier transform of the two-particle correlation function . The coefficient is artificially introduced to match the melting points of PFC and DDFT.

The second, more coarse-grained model termed PFC2 in Ref. SvenRainerPFC (), which is frequently used in the phase-field-crystal literature, can be obtained from dynamical density functional theory by assuming a constant mobility, in front of the functional derivative in Eq. (1) and a gradient expansion. The model equation then reads

(5) |

with the dimensionless density modulation.

The equation of motion is solved using the finite difference method with a semi-implicit time integration Chen1999 ().

## Iii Setup and results

In the following subsections we qualitatively compare the temperature dependence of defect diffusion as obtained by computational Brownian dynamics simulations and as predicted by the dynamical density functional theory and the phase-field-crystal model 2 (PFC2).

### iii.1 Brownian dynamics computer simulation

As a reference for the theoretical models we use Brownian dynamics computer simulations Ermak:75 () to quantify the diffusion of vacancies for different coupling strengths (well above the melting point at Haghgooie:05 ()). Following Líbal et al. Libal:07 (), we equilibrate a perfectly hexagonal crystal of particles in a rectangular, almost square box employing periodic boundary conditions while keeping one particle tagged at the origin (see Fig. 1).

Subsequently, the particle at the origin is removed and the diffusion of the defect is followed over time. In this setup the vacancy concentration is typically relatively small (of the order ). The position of the defect is determined by a Voronoi construction: The vacancy appears as one of several configurations of multiple particles with more or less than six neighbors (see Fig. 2).

The center of mass of these particles is considered as the position of the vacancy.

As was already observed for Yukawa particles in 2D by Líbal et al. Libal:07 (), we also find that the defect undergoes diffusion and that the diffusion constant increases with increasing temperature, corresponding to a decreasing coupling strength (Fig. 3). The diffusion constant of the vacancy ranges between , for , and , for .

### iii.2 Theory

In dynamical density functional theory and in the phase-field-crystal models, crystals appear as strongly modulated density fields that have the symmetry of the corresponding crystal SvenPRL (); SvenRainerPFC (). These density fields are mechanically and thermodynamically stable at low temperature or high coupling strength Teeffelen:06 (). In an equilibrium density field, the integrated density field over one Wigner-Seitz cell is equal to the probability to find a particle at the corresponding lattice site. In the approximation to the density functional theory by Ramakrishnan and Yussouff for magnetic dipoles in 2D described above this number is very close to . A number smaller than can be interpreted as a finite probability to find a vacancy, a number larger than one as a finite probability to find an interstitial, respectively.

Whereas we have addressed the short-time relaxation dynamics of
crystals in a previous paper SvenPRL (), we are here concerned
with the long-time dynamics of vacancies. For ease of computation and
to assess larger time scales we thus start with a slightly different
setup, as compared to the setup in the computer simulations: Instead
of introducing a vacancy with probability at the center of a large
two-dimensional crystal, which constitutes a problem of cylindrical
symmetry, we study the relaxation dynamics of the
quasi-one-dimensional setup sketched in Fig. 1: Our
reference state
is an equilibrium crystalline density field in a periodic rectangular
box of dimensions , where
is the equilibrium lattice constant. At time we reduce the
integrated density of all density peaks lying on an infinite row of
neighboring crystal sites along the axis by a small amount of only percent, thus rendering the problem quasi-one-dimensional (see Fig. 1). Specifically, low-amplitude Gaussians with the same center positions and similar width as the equilibrated density peaks were subtracted from the density field. Thus, the integrated density of each altered peak is smaller by a few percent than its equilibrium value. The temporal behavior of this setup is expected to be qualitatively similar to that of a cylindrical setup at late times, i.e., once the vacancy has diffused a large distance from its initial position ^{1}^{1}1Ideally one could have started with an initial density profile which describes an ideal single vacancy, but this causes strong numerical problems in the ideal-gas entropy term.. Strictly speaking the long-time diffusion constant should be calculated from the self-part of a van-Hove function in the presence of interactions (e.g. by using the test particle approach to the DDFT HopkinsJCP2010 ()). For a low vacancy density, however, the collective and self-dynamics are expected to be similar such that we have avoided the significant additional effort to implement the test-particle approach.

The outcome of the dynamical density functional theory is summarized in Fig. 4 and based on the -averaged density field

(6) |

Fig. 4a) displays the equilibrium averaged density field
, i.e., before the introduction of vacancies,
for two different coupling strength, and ,
that are close and far from the freezing transition at
, respectively ^{2}^{2}2We note that while there are large differences between the freezing point in DFT and the melting point in BD, investigating the fluid-solid transition is beyond the scope of this work. Here we simply quenched the system deep enough into the solid state where the details of the equilibrium melting process do not matter much.. The higher coupling strength
corresponds to higher and more pronounced density peaks. Removing a
fraction of particles from the row of peaks at leads to a
restoring density current from neighboring particle rows towards the
origin. This is represented by the difference between the
-averaged density fields of the perturbed and the unperturbed
systems at their original -positions (Fig. 4b):

(7) |

For small initial perturbations and at long times the envelope of approaches a Gaussian function, which broadens over time. The variance is plotted in Fig. 4c. As expected, shows a linear dependence of at long times. The long-time slope translates into a diffusion constant given by (Fig. 4d). In agreement with the Brownian dynamics simulations the diffusion constant is higher for low coupling constant of () than for the high coupling constant of (().

The same setup is studied in the PFC1 and PFC2 models. The temporal evolution of the width of the Gaussian envelope function describing the relaxing density field is shown for the PFC2 model in Figure 5. After a fast transient relaxation process its dependence is linear in time. Remarkably, the corresponding slope which is presented in the inset of Figure 5 is increasing for increasing conversely to what has been found before in simulation and DDFT.

The diffusion constant increases linearly with from corresponding to to and for equal to and . The PFC1 model gives the same incorrect trend as the PFC2 models. Again, after a transient process the system reaches a state where the relaxation towards equilibrium is getting diffusive but the slope increases with increasing . The physical reason for the incorrect trend in both variants of PFC models is that the PFC has a smoothened density profile and therefore allows for a quick diffusive current of particles from one lattice site to another. In DDFT, on the other hand, the full density profile is kept and the density decays to very small values in the interstitial region between the density peaks. Thus, particle currents between lattice sites are strongly reduced. For increasing coupling , the interstitial density drops even further down and therefore reduces vacancy diffusion more. This effect is not contained in both PFC models. Here rather the only remaining trend is set by the increasing interactions which leads to larger inter-particle forces and therefore accelerate the dynamics of vacancy diffusion. Therefore, to account for the correct defect dynamics, the temperature dependence of the mobility prefactor in the PFC models needs a proper fitting to reference data. Though this reduces the predictive power of the PFC model, it may still be useful for a fast explorative numerical study for various dynamical effects in solids provided this fitting is before a priori.

## Iv Conclusions

In conclusion, we have used dynamical density functional theory, phase-field-crystal theory and particle-resolved Brownian dynamics computer simulations to calculate the diffusion of defects in a two-dimensional crystal of repulsive dipoles. The typical diffusion coefficient of defect motion is expected to decrease with increasing system temperature as confirmed by the simulation data. This trend is reproduced in the dynamical density functional theory but not in the phase field crystal calculations. These findings show that the PFC model requires a fitting of the kinetic mobilities as a function of the thermodynamic parameters if a realistic description of trends is required to predict material properties. Given the fact that the efficiency of the PFC model is in general achieved by an optimal fitting procedure, also for structural predictions Oettel_Nestler_et_al_PRE_2012 (), such a phenomenological input of the mobility is an acceptable fact. However, it shows that clearly the dynamical density functional theory is more appropriate to predict the microscopic time evolution as a first-principle theory for Brownian systems.

Future work should address the dynamics of dipolar mixtures of colloidal particles with different dipole moments reviewAdvancesin physicalChemistry () where an equilibrium density functional for a binary mixture Kruppa_JCP_2012 () is needed. These mixtures show more complex possibilities of mixed crystals as a function of the asymmetry in their dipole moments AssoudEPL2007 (). In this case one will expect different diffusion coefficients for different defects topologies. Moreover, one should do a similar calculation for hard discs, for which a very accurate functional based on fundamental measure theory R._RothJPCMreview2010 () was proposed recently Oettel_Roth_JCP2012 (). A similar comparison can be performed in three dimensions, e.g. for hard spheres, where the phase field crystal model has been tested against density functional theory recently Oettel_Nestler_et_al_PRE_2012 (). Finally it would be interesting to consider particles with orientational degrees of freedom which form liquid crystals with interesting defect structures. In fact, phase-field-crystal models for liquid crystals have been developed LowenJPCM2010 (); WittkowskiPRE2010 () and applied AchimWittkowski () to freezing recently which opens the way to describe defect dynamics in liquid crystalline states.

###### Acknowledgements.

We thank Alexandros Chremos for many inspiring discussions. This work has been supported by the DFG through the DFG priority program SPP 1296.## References

- (1) K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988).
- (2) A. M. Turner, V. Vitelli, and D. R. Nelson, Rev. Mod. Phys. 82, 1301 (2010).
- (3) W. T. M. Irvine, V. Vitelli, and P. M. Chaikin, Nature 468, 947 (2010).
- (4) A. Pertsinidis and X.S. Ling, Nature 413, 147 (2001).
- (5) A. Pertsinidis and X. S. Ling, New J. Phys. 7, 33 (2005).
- (6) A. Libál, C. Reichhardt, and C. J. Olson Reichhardt, Phys. Rev. E 75, 011403 (2007).
- (7) V. Nosenko, G. Morfill, and P. Rosakis, Phys. Rev. Lett. 106, 155002 (2011).
- (8) V. Nosenko, S. Zhdanov, A. Ivlev, C. Knapek, and G. Morfill, Phys. Rev. Lett. 103, 015001 (2009).
- (9) V. Nosenko, S. Zhdanov, and G. Morfill, Phys. Rev. Lett. 99, 025002 (2007).
- (10) A. V. Ivlev, H. Löwen, G. E. Morfill, C. P. Royall, Complex plasmas and colloidal dispersions: particle-resolved studies of classical liquids and solids, World Scientific, Singapore (2012).
- (11) R. Evans, Adv. Phys. 28, 143 (1979).
- (12) Y. Singh, Phys. Rep. 207, 351 (1991).
- (13) H. Löwen, Phys. Rep. 237, 249 (1994).
- (14) H. Emmerich, H. Löwen, R. Wittkowski, T. Gruhn, G. I. Tot́h, G. Tegze, L. Gránásy, Advances in Physics, in press (2012).
- (15) U. M. B. Marconi and P. Tarazona, J. Chem. Phys. 110, 8032 (1999).
- (16) A. J. Archer and R. Evans, J. Chem. Phys. 121, 4246 (2004).
- (17) P. Español and H. Löwen, J. Chem. Phys. 131, 244101 (2009).
- (18) S. van Teeffelen, C. N. Likos, and H. Löwen, Phys. Rev. Lett. 100, 108302 (2008).
- (19) S. van Teeffelen, H. Löwen, R. Backofen, and A. Voigt, Phys. Rev. E 79, 051404 (2009).
- (20) K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, Phys. Rev. Lett. 88, 245701 (2002).
- (21) K. R. Elder and M. Grant, Phys. Rev. E 70, 051605 (2004).
- (22) K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, Phys. Rev. B75, 064107 (2007).
- (23) H. Emmerich, J. Phys.: Condens. Matter 21, 4103 (2009).
- (24) Y. Yu, B. Liu, and A. Voigt, Phys. Rev. B 79, 235317 (2009).
- (25) A. Jaatinen, C. V. Achim, K. R. Elder, T. Ala-Nissila, Technische Mechanik 30, 169-176 (2010).
- (26) C. V. Achim, J. A. P. Ramos, M. Karttunen, K. R. Elder,E. Granato, E. , Ala-Nissila, T. and Ying, S. C., Phys. Rev. E 79, 011606 (2009).
- (27) J. A. P. Ramos, C. V. Achim, K. R. Elder, E. Granato, S. C. Ying and T. Ala-Nissila, Phys. Rev. E 78, 031109 (2008).
- (28) C. V. Achim, Mikko Karttunen, K. R. Elder, E. Granato, S. C. Ying and T. Ala-Nissila, Phys. Rev. E 74, 021104 (2006); C. V. Achim, Mikko Karttunen, K. R. Elder, E. Granato, S. C. Ying and T. Ala-Nissila, Journal of Physics: Conference Series 100, 072001 (2008).
- (29) A. Jaatinen, C. V. Achim, K. R. Elder, and T. Ala-Nissila, Phys. Rev. E 80, 031602 (2009).
- (30) G. Tegze, L. Gránásy, G. Tóth, F. Podmaniczky, A. Jaatinen, T. Ala-Nissila, and T. Pusztai, Phys. Rev. Lett. 103, 035702 (2009).
- (31) Z. Huang and K. R. Elder, Phys. Rev. Lett. 101, 158701 (2008).
- (32) J. Mellenthin, A. Karma, and M. Plapp, Phys. Rev. B 78, 184110 (2008).
- (33) I. M. McKenna, M. P. Gururajan, and P. W. Voorhees, J. Mat. Sci. 44, 2206 (2009).
- (34) K. Wu and P. W. Voorhees, Phys. Rev. B 80, 125408 (2009).
- (35) F. Ebert, P. Dillmann, G. Maret, and P. Keim, Rev. Sci. Instr. 80, 083902 (2009).
- (36) F. J. Rogers and D. A. Young, Phys. Rev. A 30, 999 (1984).
- (37) S. van Teeffelen, C. N. Likos, N. Hoffmann, and H. Löwen, Europhys. Lett. 75, 583 (2006).
- (38) A. Archer and M. Rauscher, J. Phys. A.: Math. Gen. 37, 9325 (2004).
- (39) U. M. B. Marconi and P. Tarazona, J. Phys.: Condens. Matter 12, A413 (2000).
- (40) T. V. Ramakrishnan and M. Yussouff, Phys. Rev. B 19, 2775 (1979).
- (41) G. Tegze, G. Bansel, G. I. Tóth, T. Pusztai, Z. Fan, and L. Gránásy, J. Comput. Phys. 228, 1612 (2009); B. P. Vollmayr-Lee and Andrew D. Rutenberg, Phys. Rev. E 68, 066703 (2003); J. Zhu, L.-Q. Chen, J. Shen, and V. Tikare, ibid. 60, 3564 (1999).
- (42) J.-P. Hansen and I. R. McDonald, Theory of simple liquids, Academic Press, London, 3rd edition (2006).
- (43) S. van Teeffelen, H. Löwen, and C. N. Likos, J. Phys.: Condens. Matter 20, 404217 (2008).
- (44) D. L. Ermak, J. Chem. Phys. 62, 4189 (1975); 64 4197 (1975).
- (45) R. Haghgooie and P. S. Doyle, Phys. Rev. E 72, 011405 (2005).
- (46) P Hopkins, A. Fortini, A. J. Archer, and M. Schmidt, J. Chem. Phys. 133, 224505 (2010).
- (47) H. Löwen, E. C. Oguz, L. Assoud, R. Messina, Advances in Chemical Physics 148, 225-249, (2012).
- (48) T. Kruppa, T. Neuhaus, R. Messina, H. Löwen, J. Chem. Phys. 136, 134106 (2012).
- (49) L. Assoud, R. Messina, H. Löwen, Europhys. Letters 80, 48001 (2007).
- (50) R. Roth, K. Mecke, and M. Oettel, J. Chem. Phys. 136, 081101 (2012).
- (51) R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010).
- (52) M. Oettel, S. Dorosz, M. Berghoff, B. Nestler, and T. Schilling, Phys. Rev. E 86, 021404 (2012).
- (53) H. Löwen, J. Phys.: Condensed Matter 22, 364105 (2010).
- (54) R. Wittkowski, H. Löwen, H. R. Brand, Phys. Rev. E 82, 031708 (2010).
- (55) C. V. Achim, R. Wittkowski, H. Löwen, Phys. Rev. E 83, 061712 (2011).