A density functional approach to ferrogels

TL;DR

This paper introduces a density functional approach to model ferrogels, capturing the interplay of magnetic and elastic interactions, and predicts their elastic properties through a novel pseudopotential mapping and simulations.

Contribution

It presents a new classical density functional framework for ferrogels, using a pseudopotential to simplify interactions and enabling predictions of bulk elastic properties.

Findings

The pseudopotential mapping is validated by Monte-Carlo simulations.

The approach predicts the bulk elastic modulus under various conditions.

It extends to inhomogeneous situations using external pseudopotentials.

Abstract

Ferrogels consist of magnetic colloidal particles embedded in an elastic polymer matrix. As a consequence, their structural and rheological properties are governed by a competition between magnetic particle-particle interactions and mechanical matrix elasticity. Typically, the particles are permanently fixed within the matrix, which makes them distinguishable by their positions. Over time, particle neighbors do not change due to the fixation by the matrix. Here we present a classical density functional approach for such ferrogels. We map the elastic matrix-induced interactions between neighboring colloidal particles distinguishable by their positions onto effective pairwise interactions between indistinguishable particles similar to a "pairwise pseudopotential". Using Monte-Carlo computer simulations, we demonstrate for one-dimensional dipole-spring models of ferrogels that this mapping is justified. We then use the pseudopotential as an input into classical density functional theory of inhomogeneous fluids and predict the bulk elastic modulus of the ferrogel under various conditions. In addition, we propose the use of an "external pseudopotential" when one switches from the viewpoint of a one-dimensional dipole-spring object to a one-dimensional chain embedded in an infinitely extended bulk matrix. Our mapping approach paves the way to describe various inhomogeneous situations of ferrogels using classical density functional concepts of inhomogeneous fluids.

Equations

$$U=U_{\text{h}}+U_{\text{m}}+U_{\text{e}},$$

$$U_{\text{h}}=\sum _{i=0}^{N+1}\sum _{j>i}{u}_{\text{h}}(x_{ij});u_{\text{h}}(x)=\{\begin{array}{ll}\infty & \text{for}x<d,\\ 0& \text{for}x\ge d,\end{array}$$

$$U_{\text{m}}=\sum _{i=0}^{N+1}\sum _{j>i}{u}_{\text{m}}(x_{ij});u_{\text{m}}(x)=-2\frac{{\mu }_0}{4\pi }\frac{m^2}{x^3},$$

$$U_{\text{e}}=\frac{k}{2}\sum _{i=0}^N(x_{i,i+1}-\ell {)}^2$$

$$\begin{array}{c}{\tilde{U}}_{\text{e}}=\sum _{i=0}^{N+1}\sum _{j>i}{\tilde{u}}_{\text{e}}(x_{ij});\\ {\tilde{u}}_{\text{e}}(x)=\{\begin{array}{ll}\frac{k}{2}\left[(x-\ell {)}^2-(2d-\ell {)}^2\right]& \text{for}x<2d,\\ 0& \text{for}x\ge 2d.\end{array}\end{array}$$

$$\varphi =\frac{Nd}{L-d},$$

$$\Omega[\rho]=\mathcal{F}[\rho]+\int_{0}^{L}\rho(x)\big{(}U_{\textrm{ext}}(x)-\mu\big{)}\,dx$$

$$\mathcal{F}_{\textrm{id}}[\rho]=k_{B}T\int_{0}^{L}\rho(x)\bigg{(}\ln\big{(}\Lambda\rho(x)\big{)}-1\bigg{)}\,dx\,$$

$$\begin{gathered}\begin{aligned} \mathcal{F}_{\textrm{ex}}^{\textrm{P}}[\rho]=-k_{B}T\int_{0}^{L}&\frac{\rho(x+d/2)+\rho(x-d/2)}{2}\\ &\quad\times\,\ln\big{(}1-\eta(x)\big{)}\,dx\,,&\end{aligned}\\ \textrm{with }\eta(x)=\int_{x-d/2}^{x+d/2}\rho(x^{\prime})\,dx^{\prime}\,.\end{gathered}$$

$${\mathcal{F}}_{\text{ex}}^{\text{MF}}[\rho ]$$

$$\times g(|x-x^{\prime }|)\rho (x)\rho (x^{\prime })d{x}^{\prime }dx,$$

$$\mathcal{F}[\rho ]={\mathcal{F}}_{\text{id}}[\rho ]+{\mathcal{F}}_{\text{ex}}^{\text{P}}[\rho ]+{\mathcal{F}}_{\text{ex}}^{\text{MF}}[\rho ].$$

$$U_{\text{ext}}(x)=u(L-x)+u(x),$$

$$\frac{\delta \Omega [\rho ]}{\delta \rho (x)}=\frac{\delta \mathcal{F}[\rho ]}{\delta \rho (x)}+U_{\text{ext}}(x)-\mu \stackrel{!}{=}0,$$

$$p$$

$$=\frac{k_{B}T}{\Delta L}\ln ⟨{\left(\frac{L+\Delta L}{L}\right)}^N\exp (-\Delta U/k_{B}T)⟩.$$

$${\mathbf{F}}_{\text{mp}}(\mathbf{U})=-\frac{12\pi (1-\nu )Gd}{5-6\nu }\Delta \mathbf{R},$$

$$U_{\text{mp}}(\Delta \mathbf{R})=\frac{1}{2}{k}_{\text{mp}}(\Delta \mathbf{R}{)}^2$$

$$\left(\begin{array}[]{c}\mathbf{F}_{1}\\[2.84544pt] \mathbf{F}_{2}\end{array}\right)=\left(\begin{array}[]{cc}-k_{\textrm{mp}}&\frac{k_{\textrm{mp}}^{2}}{4\pi G}\frac{1}{r}\\[2.84544pt] \frac{k_{\textrm{mp}}^{2}}{4\pi G}\frac{1}{r}&-k_{\textrm{mp}}\end{array}\right)\cdot\left(\begin{array}[]{c}\Delta\mathbf{R}_{1}\\[2.84544pt] \Delta\mathbf{R}_{2}\end{array}\right)\,.$$

$$u_{\text{pp}}(\pm \Delta \mathbf{r})=\frac{1}{2}{k}_{\text{pp}}(r)(\Delta \mathbf{r}{)}^2,$$

$$k_{\text{pp}}(r):=\frac{k_{\text{mp}}^2}{8\pi G}\frac{1}{r}=\frac{3}{2}\frac{1-\nu }{5-6\nu }{k}_{\text{mp}}\frac{d}{r}.$$

$$U_{\text{mp}}=\sum _{i=1}^N\frac{1}{2}{k}_{\text{mp}}(x_i-i\ell {)}^2.$$

$$U_{\textrm{pp}}=\sum_{i=1}^{N}\sum_{j>i}\frac{1}{2}\frac{k_{\textrm{pp}}}{|j-i|}\big{(}x_{ij}-|j-i|\ell\big{)}^{2}\,,$$

$$U_{\text{ext}}(x)=\underset{i=1,\dots ,N}{\min }\left\{\frac{1}{2}{k}_{\text{mp}}(x-i\ell {)}^2\right\}$$

$$\begin{array}{c}\begin{array}{rl}{\tilde{u}}_{\text{e}}(x)=& \frac{k_{\text{pp}}}{2}\left[(x-\ell {)}^2-\frac{{\ell }^2}{4}\right]\text{mathds}{1}_{[0,\frac{3}{2}\ell ]}(x)\\ +\sum _{i=2}^{\infty }& \frac{k_{\text{pp}}}{2i}\left[(x-i\ell {)}^2-\frac{{\ell }^2}{4}\right]\text{mathds}{1}_{[(i-\frac{1}{2})\ell ,(i+\frac{1}{2})\ell ]}(x),\end{array}\\ \text{with }\text{mathds}{1}_{[a,b]}(x)=\{\begin{array}{ll}1& \text{for}x\in [a,b],\\ 0& \text{else}.\end{array}\end{array}$$

$$\mathcal{L}[\alpha ]$$

$$\displaystyle\qquad-\lambda\left(\int_{0}^{L}\exp\big{(}\alpha(x)\big{)}\,dx-\langle N\rangle\right)\,$$

$$\frac{d}{dt}\frac{\delta \mathcal{L}[\alpha ]}{\delta \dot{\alpha }(x)}-\frac{\delta \mathcal{L}[\alpha ]}{\delta \alpha (x)}\stackrel{!}{=}0$$

$$\ddot{\alpha }(x)=-\rho (x)\left(\frac{\delta \Omega [\rho ]}{\delta \rho (x)}+\lambda \right).$$

$$\begin{gathered}\frac{d^{2}}{dt^{2}}\int_{0}^{L}\exp\big{(}\alpha(x)\big{)}\,dx=0\quad\Leftrightarrow\\ \int_{0}^{L}\rho(x)\left(\dot{\alpha}(x)^{2}+\ddot{\alpha}(x)\right)\,dx=0\quad\Leftrightarrow\\ \lambda=\int_{0}^{L}\rho(x)\left(\dot{\alpha}(x)^{2}-\rho(x)\frac{\delta\Omega[\rho]}{\delta\rho(x)}\right)dx\,\bigg{/}\int_{0}^{L}\rho(x)^{2}\,dx\,.\end{gathered}$$

$$\gamma(x)=\bar{\rho}\int\limits_{-\infty}^{\infty}dx^{\prime}\big{(}\gamma(x^{\prime})+c(x^{\prime})\big{)}c(x-x^{\prime})$$

$$c(x)=\frac{e^{-\beta u_{\text{r}}(x)}\left[f(x)-1+e^{f(x)\left(\gamma (x)-\beta u_{\text{a}}(x)\right)}\right]}{f(x)}-\gamma (x)-1,$$

Full text

I Introduction

Classical density functional theory for inhomogeneous fluids is nowadays used for many-body systems governed by a pair potential (such as hard or soft spheres) and has found widespread applications for phase separation, freezing and interfacial phenomena, for reviews see Refs. Evans et al., 2016; Roth, 2010; Tarazona et al., 2008; Wu and Li, 2007. In a one-component system, though classical, these particles are indistinguishable in principle according to standard statistical mechanics Hansen and McDonald (2006), which implies that the interaction between two particles is the same for any pair of particles provided they are at the same separation. This standard assumption breaks down for particles embedded in an elastic polymeric gel, if the particles are anchored to the surrounding gel matrix and/or cannot diffuse or propagate through it. In this case, the particles can be labeled according to their position in the matrix with their interaction energy persistently depending on the labeling. Thus, they are distinguishable. As a basic example, this situation is encountered for a simple bead-spring model, where the springs represent the elasticity and connectivity provided by the matrix and the beads represent the particles.

Particle distinguishability leads to a different combinatorial prefactor in the classical partition function and therefore affects the entropy Gibbs (1957). However, at high density, a fluid of indistinguishable particles typically undergoes a freezing transition into a crystal. At low temperature, this crystal can be modeled by a harmonic solid Chaikin and Lubensky (2000); Español (1996); Jancovici (1967), where the neighboring particles are connected by springs. In fact, this effective model of distinguishable particles provides a good approximation as the free energy of the system is dominated by the particle interactions that overwhelm the combinatorial contribution.

In the present paper, we exploit this idea to introduce a density functional approach for ferrogels and related systems. Such ferrogels consist of magnetic colloidal particles that are embedded in a polymeric gel matrix Filipcsei et al. (2007); Ilg (2013); Menzel (2015); Odenbach (2016); Lopez-Lopez et al. (2016). Examples of similar materials are given by magnetic elastomers, magnetorheological elastomers, or magnetosensitive elastomers, where often these terms are used interchangeably. Remarkably, the structural properties and rheological behavior of these materials are governed by a competition between the magnetic particle–particle interactions and the mechanical elasticity of the embedding polymeric gel matrix. As a consequence, it is possible to tune their properties during application by modifying the magnetic interactions via external magnetic fields. Therefore, these magnetorheological systems have many prospective and promising applications, such as tunable dampers Sun et al. (2008) or vibration absorbers Deng et al. (2006).

The theoretical description of these materials is challenging. While the specific properties arise on the mesoscopic colloidal particle scale, for practical applications one is interested in the overall macroscopic response. To connect these scales in simulations, large numbers of individual particles need to be covered. For this purpose, recent work has focused on simplified minimal models. Starting on the microscale, at most a few individual polymer chains are resolved by coarse-grained bead-spring models Weeber et al. (2012); Ryzhkov et al. (2015); Weeber et al. (2015). In still more reduced mesoscopic dipole-spring models, the elasticity of the matrix is directly represented by effective spring-like interactions between the particles, combined with long-ranged magnetic dipolar interactions between them Annunziata et al. (2013); Pessot et al. (2014); Tarama et al. (2014); Ivaneyko et al. (2015); Pessot et al. (2016). More explicit approaches treat the matrix directly by continuum elasticity theory, yet at the price of reduced accessible overall particle numbers Han et al. (2013); Cremer et al. (2015, 2016); Metsch et al. (2016). A kind of compromise between the two concepts can be found in Refs. Biller et al., 2014 and Biller et al., 2015. Previous analytical approaches to link the different scales often relied on substantially simplifying idealizations concerning the positional particle configurations Ivaneyko et al. (2011); Menzel (2014); Zubarev et al. (2016). Therefore, it is desirable to develop statistical means that allow for a more profound connection between the different scales in the future. As a step in this direction, we now suggest to employ the framework of classical density functional theory for a characterization of these complex materials.

Here we mainly follow the dipole-spring concept of distinguishable particles often used for the description of ferrogels. In order to keep the models simple, we study effective one-dimensional set-ups. Such a situation is realized, for instance, for elongated magnetic particle chains embedded into an elastic matrix Huang et al. (2016), but also for magnetic filaments Sánchez et al. (2013) made, e.g., of magnetic colloidal particles connected by DNA polymer strands Dreyfus et al. (2005). We map this system with its particle-distinguishing connectivity onto another one with an effective connectivity and indistinguishable particles Denton (2007). Based on the considerations above, one expects a good agreement between real and effective connectivity at least for strong particle–particle interactions. We use Monte-Carlo computer simulations of both situations and confirm that the results agree at high packing fractions and/or strong particle interactions. This opens the way to employ statistical-mechanical theories like classical density functional theory to also describe systems of particles that are, in principle, distinguishable. For the one-dimensional model, we use the exact Percus free-energy functional for hard rods Percus (1976) combined with a mean-field theory for the elastic and dipolar interactions and minimize the resulting grand canonical free energy functional with respect to the equilibrium one-body density field.

We study two different models. In the first one, the elastic matrix is represented by harmonic springs between nearest-neighboring particles. Including thermal fluctuations, such a simple one-dimensional bead-spring model cannot show a phase transition Van Hove (1952); Cuesta and Sánchez (2004). Thermal fluctuations have a strong impact in one spatial dimension and fuel the Landau-Peierls instability Landau (2008); Peierls (1934); Chaikin and Lubensky (2000), which impedes periodic ordering. This fact is captured by our Monte-Carlo simulations, which take all contributions by thermal fluctuations into account exactly. Our mean-field density functional theory, however, introduces an artificial crystallization at low temperatures. Still, at higher temperatures we can obtain qualitative agreement between density functional theory and Monte-Carlo simulations for the pressure and compressibility of the system. These provide key material properties for the practical use of ferrogels.

Later in this paper, we turn to an extended model, including an additional external elastic pinning potential for the colloidal magnetic particles. Such pinning potentials arise when the particles are embedded in a three-dimensional elastic bulk matrix Huang et al. (2016). Similarly, the displacement of one embedded particle results in a matrix-mediated force on all other particles, so that we have additional long-ranged elastic particle–particle interactions. In combination with the pinning potential, this suppresses the Landau-Peierls instability even though our model is effectively one-dimensional. Consequently, our density functional theory immediately shows much better agreement with Monte-Carlo simulations for both, the model with real connectivity and the version mapped towards indistinguishable particles. During the synthesis of ferrogels such permanent straight one-dimensional magnetic particle chains embedded within three-dimensional ferrogel blocks are readily generated by applying external magnetic fields during the manufacturing process Collin et al. (2003); Varga et al. (2003); Günther et al. (2012); Borbáth et al. (2012); Gundermann et al. (2013); Huang et al. (2016); Gundermann et al. (2017).

We remark at this stage that the problem of mapping from distinguishable to indistinguishable particles also occurs in density functional descriptions of polymeric bead models Chandler et al. (1986a, b). Typically, in tangential bead models for hard spheres Gu et al. (2003); Słyk et al. (2016), one neglects the linking constraints of the chain and maps the excess free energy of the system onto an unconstrained hard-sphere fluid.

Our analysis paves the way to a future application of density functional theory of freezing also to two- and three-dimensional ferrogel models. There, we anticipate thermal fluctuation effects to be of less influence, leading to a better agreement with simulations. It will further be useful to characterize other particulate systems embedded in a permanent elastic matrix such as electrorheological elastomers An and Shaw (2003); Allahyarov et al. (2015) or possibly even drug carriers and compartments within biological tissue Tietze et al. (2013).

This paper is organized as follows. In Sec. II, we describe our first one-dimensional dipole-spring model and offer a method to map the real connectivity to an effective one. Next, we describe the various methods used to study our systems in Sec. III. These methods are mainly mean-field density functional theory and canonical Monte-Carlo simulations. The supplemental material sup contains an additional treatment using the Zerah-Hansen liquid-integral equation. Then we complete the discussion of the first dipole-spring model in Sec. IV, where we present results from our density functional theory and compare them to Monte-Carlo simulations. Subsequently, we proceed to our extended model in Sec. V. Following a motivation of this extended model, we again compare results from density functional theory and Monte-Carlo simulations, showing their improved agreement. Finally, in Sec. VI, we revisit our overall approach and discuss prospective uses and extensions beyond the one-dimensional models discussed here.

II Dipole-spring model

We consider the following one-dimensional dipole-spring model, which is sketched in Fig. 1. There are two outer particles at a fixed distance $L$, forming the system boundary, and $N$ mobile particles in between. All particles have a hard core of diameter $d$, which limits the closest approach of two particle centers to this distance. Additionally, all particles carry magnetic dipole moments of magnitude $m$ that all point in the same direction aligned with the system axis. Finally, each particle is connected to its nearest neighbors by a harmonic spring of spring constant $k$ and equilibrium length $\ell$.

As we will discuss below, the connectivity introduced by the harmonic springs renders the particles distinguishable. We label the particles with indices $i=0,\dots,N+1$ according to their position $x_{i}$ in ascending order. The indices $i=0$ and $i=N+1$ are used for the left and right boundary particles, respectively. The total potential energy of the system consists of three contributions

[TABLE]

i.e., the hard core repulsion $U_{\textrm{h}}$, the magnetic dipolar interaction $U_{\textrm{m}}$, and the elastic interaction $U_{\textrm{e}}$. We can write the former two as sums over the interactions between all particle pairs $i,j$ with $j>i$

[TABLE]

[TABLE]

where $\mu_{0}$ is the vacuum permeability and $x_{ij}:=|x_{j}-x_{i}|$ is the distance between a pair of particles. Since the pair interactions $u_{\textrm{h}}(x)$ and $u_{\textrm{m}}(x)$ do not depend on any particular labeling of the particles, the total interactions $U_{\textrm{h}}$ and $U_{\textrm{m}}$ are invariant under a relabeling of all particles. In contrast to that, the elastic interactions between nearest-neighbors

[TABLE]

persistently depend on the labeling and therefore render the particles distinguishable.

To facilitate a description of these systems with the tools of statistical mechanics, we map it onto a system of indistinguishable particles. This can be achieved by replacing the elastic interaction (4) with an approximative potential $\tilde{U}_{\textrm{e}}$ that can be decomposed into pairwise interactions $\tilde{u}_{\textrm{e}}(x)$. Ideally, such an approximative potential should still affect only nearest neighbors and provide the same result as Eq. (4) under realistic circumstances. We make the choice

[TABLE]

The “pseudo-spring” pair potential $\tilde{u}_{\textrm{e}}(x)$ is illustrated in Fig. 2 and consists of a harmonic well of spring constant $k$ centered around a distance $\ell<2d$. The harmonic well is cut and shifted to zero potential strength at a distance $x\geq 2d$. For two particles at a distance $x<2d$, this potential acts as a common harmonic spring. Beyond this distance, the spring “breaks” leading to zero interaction. The combination with the hard-core repulsion $u_{\textrm{h}}(x)$ in Eq. (2) limits the possible harmonic interaction to pairs of nearest-neighbors. Only nearest neighbors can be at a distance $x<2d$. Next-nearest neighbors are always at a greater distance and, thus, excluded from the interaction.

In the following, we refer to this potential as “pseudo-spring” interaction as opposed to the “real-spring” permanent connectivity between nearest neighbors. Good agreement between the real-spring system and its mapped version using pseudo-springs can be expected in situations where the pseudo-springs do not break. First, this is the case at high packing fraction, when the confinement enforces small distances between nearest neighbors. The packing fraction in our finite system is defined as

[TABLE]

because $L-d$ is the system length enclosed between the two hard boundary particles. At $\phi>(L-2d)/(L-d)$, the distance between nearest neighbors is smaller than $2d$ everywhere, such that breaking of pseudo-springs becomes impossible. Another limit is reached at high potential strength (large value of $k$) and moderate packing fraction of $\phi\gtrsim d/\ell$. Under these conditions, the harmonic well, as illustrated in Fig. 2, is deep compared to the thermal energy $k_{B}T$ and the system is sufficiently filled such that all particles are effectively trapped in the harmonic wells created by their nearest neighbors.

For fixed values of $L$, $N$, and $d$, the physical input parameters determining all interactions are the magnetic moment $m$, the spring constant $k$, and the spring equilibrium length $\ell$. From now on, we measure all energies in units of $k_{B}T$ and all lengths in units of the particle diameter $d$. This implies to measure the spring constant in units of $k_{0}=k_{B}T/d^{2}$ and the magnetic moment in units of $m_{0}=\sqrt{\frac{4\pi}{\mu_{0}}k_{B}Td^{3}}$, while the pressure and the compression modulus are given in units of $p_{0}=K_{0}=k_{B}T/d$.

III Methods

We use three different methods to study our dipole-spring model. The first and most notable one is our density functional theory (DFT) description, for which we use the pseudo-spring approximation to make particles indistinguishable. Second, we perform canonical Monte-Carlo (MC) simulations for real springs as well as for pseudo-springs as a benchmark to test our DFT results. Finally, we have also solved the Zerah-Hansen liquid-integral equation to show that our pseudo-spring approximation is meaningful beyond the scope of DFT, see the supplemental material sup for results and a description of the method.

III.1 Density functional theory

The central statement of classical DFT is that for a fixed temperature $T$ and interparticle pair potential $u(x)$, the Helmholtz free energy $\mathcal{F}[\rho]$ is a unique functional of the one-body density distribution $\rho(x)$. Likewise, there is a unique grand canonical free energy functional $\Omega[\rho]$ describing the system when it is exposed to an external potential $U_{\textrm{ext}}(x)$ and a particle reservoir at chemical potential $\mu$. This grand canonical free energy functional has the form Tarazona et al. (2008)

[TABLE]

and is minimized by the equilibrium one-body density profile $\rho_{\textrm{eq}}(x)$. The minimum $\Omega[\rho_{\textrm{eq}}]$ corresponds to the thermodynamic grand canonical free energy in equilibrium.

Unfortunately, the exact free energy functional $\mathcal{F}[\rho]$ is usually unknown, so that one has to resort to approximations. These approximations usually start by splitting the free energy functional $\mathcal{F}[\rho]=\mathcal{F}_{\textrm{id}}[\rho]+\mathcal{F}_{\textrm{ex}}[\rho]$ into the exact free energy for the ideal gas

[TABLE]

with $\Lambda$ the thermal de Broglie wavelength, plus an excess contribution $\mathcal{F}_{\textrm{ex}}[\rho]$. For some special problems in one spatial dimension, the exact excess contribution can be derived Tutschka and Kahl (2000). One such example is the Percus excess functional Percus (1976) for the one-dimensional hard-rod fluid,

[TABLE]

It takes one-dimensional hard repulsions exactly into account and, thus, provides a good starting point for the construction of a functional describing our dipole-spring model. Here, we combine it with an approximate mean-field excess functional accounting for the soft pair interactions consisting of our pseudo-spring pair potential $\tilde{u}_{\textrm{e}}(x)$ and the magnetic dipolar pair interaction $u_{\textrm{m}}(x)$,

[TABLE]

where the distribution function $g(x)$ satisfies the no-overlap condition $g(x)=0$ for $x<d$. The mean-field approximation assumes that the pair potentials are soft enough to regard the particle positions as basically uncorrelated Tarazona et al. (2008). Here we make the simplifying assumption that $g(x)=1$ for all distances $x>d$. In total, our free energy functional is given by

[TABLE]

The boundary of our finite systems consists of the leftmost and rightmost particles, which are fixed but otherwise identical to the enclosed particles, see Fig. 1. Their influence on the enclosed density profile enters via an external potential

[TABLE]

where $u(x)=u_{\textrm{h}}(x)+u_{\textrm{m}}(x)+\tilde{u}_{\textrm{e}}(x)$. This completes our grand canonical free energy functional $\Omega[\rho]$.

Functional derivation of Eq. (7) leads to the Euler-Lagrange equation

[TABLE]

which can be used to determine the equilibrium density profile minimizing $\Omega[\rho]$. In practice, however, we numerically calculate our equilibrium density profile $\rho(x)$ by performing a dynamical relaxation of $\Omega[\rho]$ Löwen et al. (1993). This scheme fixes the average particle number $\langle N\rangle$ instead of the chemical potential $\mu$ and is described in detail in the appendix.

After the relaxation, we have access to the grand canonical free energy $\Omega$. This enables us to calculate a pressure p=-\frac{\partial\Omega}{\partial L}\big{|}_{\langle N\rangle,T} and a compression modulus K=-L\frac{\partial p}{\partial L}\big{|}_{\langle N\rangle,T} by varying the system length $L$ at fixed average particle number $\langle N\rangle$ and probing the corresponding change in $\Omega$.

III.2 Monte-Carlo simulation

We perform canonical Monte-Carlo (MC) simulations at fixed particle number $N$, system length $L$, and temperature $T$ Frenkel and Smit (2002). After equilibrating the systems, we sample the pressure $p$, the compression modulus $K$, and the equilibrium density profile $\rho(x)$. To sample the pressure, we affinely deform the system by a factor $(L+\Delta L)/L$ and probe the corresponding change in Helmholtz free energy $F(N,L,T)$. $\Delta L$ is a small change in system length. It can be shown that the pressure is related to the acceptance ratio of such volume moves by Harismiadis et al. (1996)

[TABLE]

$\Delta U$ is the change in system energy associated with the volume move and $\langle\cdot\rangle$ denotes the ensemble average. In order to capture the pressure contributions of the hard repulsions in our systems, the volume moves must be compressive ($\Delta L<0$). Given the pressure, the compression modulus can be calculated using K=-L\frac{\partial p}{\partial L}\big{|}_{N,T}.

IV Results

In the following, we present results for our one-dimensional dipole-spring model. First we concentrate on a non-magnetic system to test the feasibility of the mapping onto indistinguishable particles. Then we add the magnetic interaction and discuss how this affects the density profile and the pressure in our systems. Finally, we turn to the thermodynamic compression modulus, which is a key quantity to characterize ferrogel systems as it can be controlled by changing the magnetic properties.

IV.1 Non-magnetic system

First of all, we confirm within our MC-simulations that the pseudo-spring pair potential is an appropriate replacement for real springs between nearest neighbors. Figure 3 compares the equations of state $p(\phi)$ for both situations in a system of length $L=51d$. The spring parameters $k=40k_{0}$ and $\ell=1.5d$ are the same as in Fig. 2.

Using these parameters, we can confirm that at packing fractions $\phi\gtrsim d/\ell=\frac{2}{3}$ the mapping to indistinguishable particles using pseudo-springs works well.

Let us now compare MC and DFT results using the same parameters. Figure 4 shows three density profiles $\rho(x)$ at a packing fraction $\phi=\frac{2}{3}$, one from the real-spring MC, one from the pseudo-spring MC, and one from DFT. While the two MC density profiles expectedly agree with each other and display a liquid-like behavior near rigid boundaries, the DFT density profile is qualitatively different and resembles a crystal.

This crystalline appearance displayed by the DFT is unphysical. Our system is one-dimensional, all particle interactions are short-ranged and there are no external fields. For such systems, the existence of a phase transition can be ruled out Cuesta and Sánchez (2004); Dyson (1969); Ruelle (1968); Van Hove (1952). This fact is accurately captured by our MC-simulations that display a liquid phase even for this high value of $k=40k_{0}$, as they explicitly include all effects of thermal fluctuations. In one spatial dimension, thermal fluctuations have a particularly strong effect. They can escalate into long-ranged fluctuations scaling in amplitude with the system size, capable of destroying periodic ordering. This is the well-known Landau-Peierls instability Landau (2008); Peierls (1934); Chaikin and Lubensky (2000). Within the DFT, some thermal fluctuations are introduced by the ideal gas term [see Eq. (8)], which pushes the system towards disorder. However, the mean-field term [see Eq. (10)] excludes other contributions by fluctuations. We conclude that this term is responsible for the unphysical crystallization. Our conjecture is supported by setting $k=0$ and $m=0$, i.e. setting the mean-field term to zero. Then, we recover the hard-rod fluid also for the DFT and find perfect agreement with MC-simulations, see Fig. 5.

The Landau-Peierls instability is well-known to be most prominent in one spatial dimension. In future studies in two and three dimensions, we therefore expect a significantly weaker effect of the thermal fluctuations, which should lead to a better agreement between simulations and mean-field DFT. For now, we achieve qualitative agreement between DFT and MC by raising the temperature (which means decreasing $k$) until the DFT system enters the liquid state.

Figure 6 shows a comparison between density profiles as well as equations of state for the same systems as in Fig. 4, but with a ten times lower spring constant $k=4k_{0}$.

The depth of the harmonic well in the pseudo-spring potential is now of the order of $k_{B}T$ so that breaking of pseudo-springs is a common event. This renders the density profile in DFT more liquid-like, which improves the agreement with the pseudo-spring MC substantially. At the same time though, the pseudo-spring mapping becomes a bad approximation for the real connectivity. Only at high packing fractions, where the confinement prevents pseudo-spring breaking, we can reach agreement between the pseudo-spring and real-spring MC again.

IV.2 Influence of magnetic interactions

We now activate the magnetic dipolar interactions and discuss the resulting changes for our systems. Figure 7 demonstrates that increasing $m$ increases the amplitudes of all peaks in the DFT, whereas in the pseudo-spring MC only the first peak is affected.

Again, our mean-field DFT seems to overestimate the tendency to form a patterned structure because of its incomplete representation of thermal fluctuations. The reason is that, effectively, the particles do not fluctuate as much around their average positions and do not come as close to each other, where the pseudo-spring interaction and the dipolar interaction increase (the latter with inverse cubic distance). As a consequence, the DFT underestimates the averaged strength of the pair interactions in the system. This becomes apparent in the equation of state, where the MC predicts a much stronger downwards shift when increasing the magnetic moments (see the insets of Fig. 7).

As a liquid-state approach alternative to DFT, we have also solved the Zerah-Hansen liquid-integral equation. Corresponding results in comparison to MC-simulations can be found in the supplemental material sup .

IV.3 Thermodynamic compression moduli

Finally, we evaluate the elastic moduli of the DFT and pseudo-spring MC systems for various magnitudes $m$ of the magnetic moment. We present them as a function of packing fraction $\phi$ in Fig. 8. The DFT predicts only a very slight downward shift of the compression modulus when increasing the magnetic moment. In contrast to that, the shift is significantly more pronounced in the pseudo-spring MC. Additionally, the overall value of the compression modulus at high packing fractions is lower in the DFT.

These observations are in line with our earlier results. The mean-field DFT overestimates the tendency to form patterned structures. It therefore underestimates both, contributions by magnetic and elastic pair potentials. If the fluctuations of the particle positions were more pronounced, there would be more emphasis on configurations with strong elastic and magnetic interactions and their influence on the compression modulus would be stronger.

V Embedding into the elastic matrix

So far, we have considered a simple one-dimensional dipole-spring model. There, the elastic matrix is solely represented by springs between nearest-neighbor magnetic particles. Now we turn to an extended model, explicitly describing a single linear chain of magnetic particles that is embedded in a three-dimensional elastic matrix.

V.1 Dipole-spring model for a linear embedded chain

We begin by constructing an effective pinning potential $U_{\textrm{mp}}$ for the embedded particles within the three-dimensional matrix as well as an effective pair interaction $u_{\textrm{pp}}$ between two embedded particles mediated by the matrix. Subsequently, we translate these potentials into a network of springs describing the overall elastic interactions.

If a single spherical particle of diameter $d$ embedded in an infinitely extended homogeneous elastic matrix is displaced by a vector $\Delta\mathbf{R}$, it distorts the elastic environment. Then the restoring force $\mathbf{F}_{\textrm{mp}}$ that the matrix exerts onto the particle is given by Phan-Thien (1993); Phan-Thien and Kim (1994); Puljiz et al. (2016); Puljiz and Menzel (2016)

[TABLE]

where $\nu$ is the Poisson ratio, that equals $\nu=1/2$ for incompressible matrices, and $G$ is the shear modulus. The force can be connected via $\mathbf{F}_{\textrm{mp}}=-\nabla U_{\textrm{mp}}$ to a harmonic potential

[TABLE]

with the spring constant k_{\textrm{mp}}:=\big{(}12\pi(1-\nu)Gd\big{)}/\big{(}5-6\nu\big{)}.

Now we consider two embedded particles, labeled as “1” and “2”, respectively. Upon displacing these particles by vectors $\Delta\mathbf{R}_{1}$ and $\Delta\mathbf{R}_{2}$, they experience the forces $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$. In our one-dimensional set-up, we only consider forces and displacements along the particle center-to-center vector $\mathbf{r}$. To first order in the particle distance, i.e.  to order $1/r$ with $r=|\mathbf{r}|$ we then obtain Puljiz et al. (2016); Puljiz and Menzel (2016); Phan-Thien and Kim (1994); Phan-Thien (1993):

[TABLE]

Here, entries on the diagonal represent the restoring pinning forces (15). The off-diagonal contributions result from the matrix distortions that are caused by the displacement of one particle but affect the other embedded particle. To construct an effective pair potential, we here only consider symmetric situations where $\Delta\mathbf{R}_{1}=-\Delta\mathbf{R}_{2}$. Then, the change in distance between the two particles is $\Delta\mathbf{r}=\Delta\mathbf{R}_{1}-\Delta\mathbf{R}_{2}=2\Delta\mathbf{R}_{1}$. Per particle, we can thus rewrite the effective matrix-mediated inter-particle interaction as a function of $\Delta\mathbf{r}$ in the form of an effective potential

[TABLE]

with a distance-dependent spring constant

[TABLE]

Using the two potentials $U_{\textrm{mp}}$ and $u_{\textrm{pp}}$ as an input, we now motivate an extended dipole-spring model for a linear, initially homogeneous chain of $N$ particles embedded into an elastic matrix, see Fig. 9. We label the particles from left to right by $i=1,\dots,N$, according to their equilibrium positions $x_{i}^{0}:=i\ell$ within the chain. The total pinning potential based on Eq. (16) then becomes

[TABLE]

To approximate the matrix-mediated particle–particle interactions between two particles $i,j$ we replace the $1/r$ dependence of the spring constant (19) by $1/|j-i|\ell$. Thus, we have for the total interaction between all pairs of particles

[TABLE]

where $k_{\textrm{pp}}:=\frac{3}{2}\frac{1-\nu}{5-6\nu}\frac{d}{\ell}k_{\textrm{mp}}$. Essentially, this means that each particle $i$ is connected to all other particles $j$ with harmonic springs of spring constant $k_{\textrm{pp}}/|j-i|$ and spring equilibrium length $|j-i|\ell$, see Fig. 9. From now on, we assume incompressibility of the elastic matrix and set $\nu=1/2$.

Of course, in this extended dipole-spring model the particles are again distinguishable by their positions. As before, it needs to be mapped to the use in our DFT. For this purpose, we replace the harmonic springs in the model by “pseudo-springs”, following the ideas outlined in Sec. II. To include the pinning potential $U_{\textrm{mp}}$, we use an external potential consisting of a series of $N$ harmonic wells

[TABLE]

as depicted in Fig. 10.

To represent the network of springs in Eq. (21) between one particle and all other particles by “pseudo-spring” interactions, we use

[TABLE]

This pair potential is illustrated in Fig. 11 and consists of a series of harmonic wells. The spring constants of the wells decay with the neighbor number $i$ from the origin just like the individual springs connecting a particle to all other particles in the dipole-spring model, see Fig. 9. Furthermore, the boundaries of the wells are shifted to zero so that we have a vanishing pair potential at infinite distance. Since the depth of the wells roughly decays as $1/x$, the interaction is long-ranged.

The width of the wells in both, the external and the pair potential, is given by $\ell$. This width should be larger than $d$ and here we choose $\ell=2d$. Thus, both potentials in principle allow more than one particles to occupy a single well. However, the external potential pinning the particles to their equilibrium positions in the matrix is relatively strong. The particles should, therefore, remain centered in their respective wells on average.

V.2 Results

We now discuss our results for our extended dipole-spring model for a magnetic particle chain embedded in a three-dimensional elastic matrix. To this end, we consider a chain of $N=40$ particles, with an equilibrium interparticle distance $\ell=2d$ and spring constants $k_{\textrm{pp}}=4k_{0}$, $k_{\textrm{mp}}=\frac{8\ell}{3d}k_{\textrm{pp}}=\frac{16}{3}k_{\textrm{pp}}$.

To evaluate the contribution of the magnetic chain to the pressure and compression modulus as well as to evaluate the DFT numerically, we address one part of the elastic matrix of length $L=100d$ that contains the magnetic chain. This choice of $L$ is arbitrary, the only requirement for $L$ is to be larger than the total equilibrium length $N\ell=80d$ of the chain by a reasonable amount.

Figure 12 shows density profiles obtained from DFT, pseudo-spring MC, and real-spring MC when setting the magnetic moment to $m=0$. In contrast to our former dipole-spring model, periodic structures appear here in the density profiles resulting from all three methods, even though this model is still effectively one-dimensional. The reason is, first, that the elastic particle–particle interaction decays only slowly with the distance and is effectively long-ranged. This applies to both, the real-spring and the pseudo-spring version. Second, we have a pinning potential suppressing large amplitude fluctuations of the particles around their pinning positions. Together, both contributions counteract the Landau-Peierls instability and can facilitate periodic structures also in one spatial dimension Cuesta and Sánchez (2004); Wong et al. (1990); Dyson (1969); Ruelle (1968); Van Hove (1952). In this way, the role of thermal fluctuations is substantially reduced, and our mean-field DFT performs much better when compared to the MC-simulations.

Let us now address the pressure and compression modulus. As before, they can be determined by probing the energetic change of the system upon deformation. However, we keep in mind that we have one system (the chain of particles) embedded into another system (the surrounding matrix, where here we present our results for one part of length $L$ of this infinitely extended matrix). When we perform a small affine deformation $\Delta L$ of the part of the matrix containing the chain, we alter the properties of the embedded system. In particular, the equilibrium distance $\ell$ between the embedded particles changes by a factor $(L+\Delta L)/L$. In our approach, this affects the pinning positions $x_{i}^{0}=i\ell$ of the particles as well as the spring constant $k_{\textrm{pp}}=\frac{3d}{8\ell}k_{\textrm{mp}}$, which is accounted for in the energetic change upon deformation. Furthermore, what we can calculate from this energetic change are only the contributions $\Delta p$ and $\Delta K$ of the embedded chain to the overall pressure and compression modulus of the composite. To obtain the overall pressure or compression modulus of the whole composite, the energetic change associated with the macroscopic deformation of the three-dimensional matrix would need to be included as well, which is beyond our particle-based approach.

Figure 13 shows the contribution $\Delta p$ to the pressure as well as the contribution $\Delta K$ to the compression modulus as a function of the magnetic dipole moment. We can observe a linear decrease with $m^{2}$ in both quantities with good agreement between the DFT and MC calculations. The linear behavior is expected, because the magnetic interaction energy scales with $m^{2}$ and the particle chain remains relatively homogeneous while increasing the magnetic moment.

As a final result, we present the contribution of the embedded chain to the stress-strain behavior of the composite material. For this purpose, we compress the surrounding matrix by $\Delta L$ and measure the pressure contribution of the embedded chain as a function of this compression. The results are shown in Fig. 14 for values of the squared magnetic moment in the range $m^{2}=0.0\,{m_{0}}^{2}-1.0\,{m_{0}}^{2}$. At vanishing magnetic moment, the pressure contribution slightly increases when compressing the system. This is probably due to entropic effects that favor an elongated chain and, thus, work against a compression combined with a slight increase in the spring constant $k_{\textrm{pp}}$ in our description. Increasing $m^{2}$, however, leads to a stronger magnetic attraction between the particles. This renders an overall compression more favorable. Since decreasing the particle distance also enhances the magnetic attraction, we have a negative pressure contribution that increases in magnitude when compressing the system further.

Again, we can observe good agreement between DFT and MC. The best agreement is observed at small compressions and low magnetic moments. Remarkably, the pseudo-spring and the real-spring MC agree exceptionally well at all considered values of $\Delta L$ and $m^{2}$. This demonstrates, that our approach to map the spring network to effective interactions between indistinguishable particles is a promising approach also beyond the scope of our mean-field DFT.

VI Conclusions

In summary, we have proposed a density functional theory to address ferrogel model systems, here evaluated in one spatial dimension. These systems are in principle non-liquid, because the particles are arrested by the elastic matrix surrounding them. To enable the investigation with statistical-mechanical theories, we map the elastic interactions onto effective pairwise interactions and, thus, make the particles indistinguishable.

The one-dimensional nature of the ferrogel model systems investigated here poses a challenge, because thermal fluctuations have a special impact in one dimension. Fluctuations can become long-ranged and destroy periodic structural order. These fluctuations, driving the Landau-Peierls instability, are not resolved within our mean-field density functional theory. Therefore, within our first dipole-spring model we observe deviations from Monte-Carlo simulations where these fluctuations are included.

In a second, more advanced approach, we explicitly model a linear particle chain embedded into a three-dimensional matrix. Within this model, the Landau-Peierls instability is counteracted by a stronger long-ranged coupling between the particles and a pinning potential that localizes the particles within the elastic matrix. Since the role of the fluctuations is therefore reduced, our density functional theory now provides results that are in good agreement with Monte-Carlo simulations. Numerous experimental realizations of such a systems exist Collin et al. (2003); Varga et al. (2003); Günther et al. (2012); Borbáth et al. (2012); Gundermann et al. (2013); Huang et al. (2016); Gundermann et al. (2017), see particularly the set-up in Ref. Huang et al., 2016.

For the future, it would be promising to extend the concept proposed here to higher spatial dimensions, that is to two-dimensional sheets of ferrogels or full three-dimensional samples. In those dimensions, the Landau-Peierls instability will be less relevant. We expect that especially for regular crystal-like particle arrangements, where the one-body density is regularly peaked, density functional theory is reliable and provides a useful framework to study the properties of these promising materials. It will be challenging to extend the present analysis to include the dynamics of the colloidal particles by using the concept of dynamical density functional theory Archer and Evans (2004); Marconi and Tarazona (1999); Español and Löwen (2009); Löwen (2017). For particles of different sizes, or different dipole moments, the same ideas can be used to map the system onto binary and multicomponent systems. Moreover, orientational degrees of freedom, such as rotating dipole moments or anisotropic particle shapes, can be tackled by density functional theory as well, both in statics Groh and Dietrich (1996); Klapp and Patey (2000) and dynamics Rex et al. (2007); Menzel et al. (2016), and can therefore be treated within the same framework proposed here.

Our results show that density functional theory can be used to describe non-liquid systems like ferrogels, still leading to reasonable results. More generally, we have established that mapping bead-spring models to effective potentials is a feasible approach to make them accessible to statistical-mechanical theories.

These theories often make use of correlation functions as an input Denton and Ashcroft (1991), which here are related to the particle distribution in the ferrogel. Experimental extraction of the particle distribution and the corresponding correlation functions is still challenging Gundermann et al. (2017). However, this route could be explored in the future once the available experimental techniques for particle detection in ferrogel materials are more advanced and particularly can address larger system sizes. These correlation functions could then help to construct effective pair potentials representing the real connectivity in the gel Hansen-Goos et al. (2006), providing a formal route for the mapping onto systems of indistinguishable particles.

Appendix A Numerical relaxation scheme to obtain the equilibrium density profile from our DFT

Here, we describe in detail our numerical relaxation scheme to obtain the equilibrium density profile $\rho(x)$ minimizing $\Omega[\rho]$ within our DFT. Instead of directly solving Eq. (13), we perform a dynamical relaxation of the Lagrange functional

[TABLE]

with respect to the logarithmic density profile \alpha(x)=\ln\big{(}\rho(x)\big{)} Ohnesorge et al. (1994). Minimizing with respect to the logarithmic density profile ensures that \rho(x)=\exp\big{(}\alpha(x)\big{)} remains positive during the relaxation. The artificial kinetic term $\int_{0}^{L}\frac{1}{2}\dot{\alpha}(x)^{2}\,dx$ drives $\alpha(x)$ and thus $\rho(x)$ towards the minimum in the grand canonical free energy $\Omega[\rho]$. The Lagrange multiplier $\lambda$ with the corresponding constraint \langle N\rangle=\int_{0}^{L}\exp\big{(}\alpha(x)\big{)}\,dx=\int_{0}^{L}\rho(x)\,dx allows us to set the average particle number $\langle N\rangle$ instead of the chemical potential $\mu$. This is more convenient for evaluating the pressure and the compression modulus defined as p=-\frac{\partial\Omega}{\partial L}\big{|}_{\langle N\rangle,T} and K=-L\frac{\partial p}{\partial L}\big{|}_{\langle N\rangle,T}, respectively.

Solving the Euler-Lagrange equation

[TABLE]

then leads to the equation of motion for $\alpha(x)$

[TABLE]

The Lagrange multiplier $\lambda$ is determined by

[TABLE]

In order to perform the numerical relaxation, we discretize the system into $n$ equally spaced sampling points such that one particle diameter $d$ is represented by $100$ points. The density profile, the potentials, and the radial distribution function are all defined on this numerical grid. All integrals appearing in the calculation of the “acceleration” $\ddot{\alpha}(x)$ can then be solved numerically, making use of Fast-Fourier-Transforms in the case of convolution integrals.

We iterate the equation of motion for $\alpha(x)$ forward in time using the standard Velocity-Verlet scheme, obtaining the “velocity” $\dot{\alpha}(x)$ and an update for the density profile $\rho(x)$ in each time step $\Delta t$. To ensure that the constraint $\langle N\rangle=\int_{0}^{L}\rho(x)\,dx$ remains fulfilled, we renormalize $\rho(x)$ after each time step. The time step is variable and increases by a factor $1.1$ up to a maximum $\Delta t_{\textrm{max}}=0.01$ when the grand canonical energy has decreased for $5$ consecutive time steps. The decrease in energy is monitored by the “power” $P=\int_{0}^{L}\dot{\alpha}(x)\ddot{\alpha}(x)\,dx$, which is supposed to be positive. If $P\leq 0$ occurs, we set $\dot{\alpha}(x)=0$ and halve the time step. We consider the density profile $\rho(x)$ sufficiently close to equilibrium when our measure for the error $\varepsilon:=\sqrt{\int_{0}^{L}\ddot{\alpha}(x)^{2}\,dx}$ becomes smaller than $10^{-6}$. At that stage, the left hand side of Eq. (26) as well as $\lambda$ are close to zero, so that we have $\frac{\delta\Omega[\rho]}{\delta\rho}(x)\approx 0$ as required by Eq. (13).

Acknowledgements.

We thank Mate Puljiz for a stimulating discussion concerning the matrix-mediated interactions. PC, AMM, and HL thank the Deutsche Forschungsgemeinschaft for support of this work through the priority program SPP 1681.