Classical density functional theory for a two-dimensional isotropic ferrogel model with labeled particles
In this study, we formulate a density functional theory (DFT) for systems of labeled particles, considering a two-dimensional bead-spring lattice with a magnetic dipole on every bead as a model for ferrogels. On the one hand, DFT has been widely studied to investigate fluid-like states of materials, in which constituent particles are not labeled as they can exchange their positions without energy cost. On the other hand, in ferrogels consisting of magnetic particles embedded in elastic polymer matrices, the particles are labeled by their positions as their neighbors do not change over time. We resolve such an issue of particle labeling, introducing a mapping of the elastic interaction mediated by the springs onto a pairwise additive interaction (pseudo-springs) between unlabeled particles. We further investigate magnetostriction and changes in the elastic constants under altered magnetic interactions employing the pseudo-spring potential. It is revealed that there are two different response scenarios in the mechanical properties of the dipole-spring systems: while systems at low packing fractions are hardened as the magnetic moments increase in magnitude, at high packing fractions softening due to diminishing effects from the steric force, associated with increases in the volume, is observed. Validity of the theory is also verified by Monte-Carlo simulations with both real and pseudo-springs. We expect that our DFT approach may shed light on an understanding of materials with particle inclusions.
In statistical mechanics, indistinguishability of particles and consequently the correct Boltzmann counting play an essential role, see, e.g., Refs. Pathria (1996); Kardar (2007); Huang (2009). The “Gibbs paradox” is well-known in this regard and the extensivity of entropy is recovered by introduction of the factor, which corrects the number of microstates by the number of permutations of particles. Commonly, the factor is regarded as a remnant of quantum mechanics in the classical limit, in which identical particles are inherently indistinguishable. In contrast to such a point of view, however, the factor should be consistently interpreted based on an “informatic” definition of entropy Cates and Manoharan (2015). Accordingly, a modified term, “undistinguished” particles Sethna (2006), has also been proposed. Consequently, even though the classical particles such as colloidal particles are undoubtedly distinguishable, the statistical mechanics with the correction describes the macroscopic behaviors of such systems successfully Swendsen (2008); Frenkel (2014), as far as one ignores detailed differences between particles Jaynes (1957, 1992) and leaves the particles unlabeled.
At the very microscopic level, i.e., at the atomic scale, a system consisting of identical particles is invariant under permutations of the particles and the particles are unlabeled in principle. If one considers a mesoscopic length scale and employs a coarse-grained description Doi and Edwards (1986); Menzel (2019), however, it may become necessary to distinguish between particles that are physically identical. Such a scenario can emerge if the particles are permanently localized with respect to their neighbors, thereby rendering the system non-ergodic on the relevant energy scale. It is thus a major challenge for statistical mechanical theories to describe a model with labeled particles Cremer et al. (2017) or to keep track of a single localized particle Schindler et al. (2019).
Here we develop a statistical description for a new class of composite materials, which consist of magnetic particles and an elastic polymer matrix Filipcsei et al. (2007); Ilg (2013); Menzel (2019). In these materials, called ferrogels or magnetorheological elastomers, the dynamical trajectories of the magnetic particles are frequently strongly constrained by the polymeric environment Gundermann and Odenbach (2014); Landers et al. (2015). Such a magneto-mechanical coupling can even be enhanced by directly anchoring the polymers on the surface of magnetic particles Frickel et al. (2011); Messing et al. (2011); Ilg (2013); Roeder et al. (2015). Hence, the elastic properties of the materials can be tuned from outside by non-invasive applications of magnetic fields Stepanov et al. (2008); Stolbov et al. (2011); Stolbov and Raikher (2019). As a further consequence of this coupling the particles cannot exchange their positions due to the fixation by the elastic medium.
Various studies have been conducted to theoretically understand the behavior of such ferrogels with different description levels from the microscopic scale resolving the individual polymer particles to the macroscopic hydrodynamic/thermodynamic theory. For many practical purposes, one may neglect the thermal motions of the magnetic particles Menzel (2019). In particular, a mesoscopic dipole-spring model has been adopted to study the elastic and dynamical properties of ferrogels Pessot et al. (2016, 2018); Goh et al. (2018). The matrix-mediated interaction between magnetic particle inclusions has also been revealed in terms of continuous elastic backgrounds Biller et al. (2014, 2015); Cremer et al. (2015, 2016); Puljiz et al. (2016); Puljiz and Menzel (2017). Furthermore, microscopic descriptions of ferrogels via coarse-grained molecular dynamics simulations enable us to probe the role of thermal motions of the magnetic particles explicitly Weeber et al. (2012, 2015, 2019).
Here, we merge several of the aspects mentioned above. Our goal is to formulate a statistical mechanical theory for ferrogels in a dipole-spring model with thermal fluctuations taken into account. The most challenging problem that has to be addressed along the way arises from the fact that the particles in ferrogels are strictly labeled by their positions as in lattice systems, for instance, the classical Ising or XY models Goldenfeld (1992) and harmonic crystals Ashcroft and Mermin (1976). Accordingly, correcting the number of microstates by the factor does not apply to a statistical description of ferrogels and a permutation of particles will cause a change in energy (physically this results in strong distortions of the surrounding elastic matrix). While computational approaches such as Monte-Carlo (MC) simulations are still feasible Weeber et al. (2019), the formulation of a statistical mechanical theory is severely complicated by the inherent composite nature of the ferrogels, in contrast to, for instance, harmonic crystals. In practice, one would need to take into account the nonlinearity stemming from the steric and magnetic interactions.
One natural candidate for a statistical theory is classical density functional theory (DFT) Evans (1979); Lutsko (2010); Evans et al. (2016), which has been probed to be successful for variety of systems Löwen (2002), ranging from simple classical fluids Ebner et al. (1976) to systems showing a freezing transition Ramakrishnan and Yussouff (1979), from hard-spheres Roth (2010) to hard convex particles Wittmann et al. (2016) and also for two-dimensional systems Roth et al. (2012); Wittmann et al. (2017); Lin and Oettel (2018), including dipolar or electrostatic interactions Zimmermann et al. (2016); Roth and Gillespie (2016), and capturing the spinodal decomposition dynamics Archer and Evans (2004) in an adiabatic approximation for time-dependent systems. However, due to the particle labeling, a direct application of the machinery of DFT to ferrogels is not possible. To this end, we will map the elastic interaction onto an appropriate pairwise pseudo-potential Cremer et al. (2017) between unlabeled particles, which allows us to formulate a DFT. Ultimately, we aim at investigating the elastic properties of the dipole-spring systems within the DFT framework. By comparison to the MC simulations, the validity of the theory is confirmed.
The paper is organized as follows: In Sec. II, we introduce a two-dimensional model for ferrogels. Sec. III describes detailed procedures of the mapping and the subsequent formulation of a DFT. Combining MC simulations and DFT calculations, two distinctive response scenarios in elastic properties to the change of the magnetic moment are identified in Sec. IV. Lastly, a summary and a discussion are given in Sec. V.
Ii The dipole-spring model
We consider a bead-spring model Doi and Edwards (1986) in terms of a periodic two-dimensional hexagonal lattice, as illustrated in Fig. 1. There are identical magnetic particles and identical harmonic springs connecting the nearest neighbors. We denote the position and the dipole moment of the th particle of diameter by and , respectively. The total Hamiltonian of the model system is given by the sum of the kinetic part and the interaction Hamiltonian which consists of three parts in the form of
Among these three terms, the magnetic part and the steric part , can be written as
where . For isotropic interactions, the vector in the argument can be replaced by .
First, the two-body magnetic dipole-dipole interaction energy between magnetic particles reads
where is the vacuum permeability. Henceforth, we assume that the magnetic moment is constant in time and for all particles, i.e., regardless of and for the sake of simplicity. Moreover, we constrain ourselves to isotropic magnetic interactions in the two-dimensional plane by further assuming that is perpendicular to the lattice plane, so that the magnetic two-body interaction energy can be written in a simpler form as
The second term in Eq. (1) corresponds to the elastic energy of the harmonic springs and reads
where is the spring constant and is the rest length of the springs. Apparently, each spring connects a certain prescribed pair of particles, and our summation persistently runs only over nearest-neighboring pairs as indicated by the angular bracket. Consequently, all particles are labeled by the predetermined ordering on a lattice, the energetic memory of which being recalled by . We remark that Eq. (5) cannot be cast into the form of with a -independent function .
Lastly, the two-body steric repulsion energy is taken as a hard-core potential of the form
and completes the interaction Hamiltonian in Eq. (1). This term prevents the possible divergence of the magnetic interactions at . Following Refs. Bernard et al. (2009); Roth et al. (2012); Engel et al. (2013); Lin and Oettel (2018), we introduce a dimensionless packing fraction defined as the ratio of the space occupied by particles to the two-dimensional “volume” (the term “volume” is used for area throughout the paper) of the system. As we might consider the systems under constant pressure (see Sec. IV.1), however, the volume of our system is not necessarily a fixed variable. We therefore define a reference volume in which the springs are in their rest state. Accordingly, the packing fraction of a reference system with this volume is given by . Apart from , which will be used as the unit of length, depends only on the diameter of particles and therefore we employ it as a model parameter representing the steric repulsion. In contrast to that, the conventional packing fraction defined in terms of the “actual” volume of the system, may change as the volume increases/decreases in response to a decrease/increase of the pressure or an increase/decrease of the magnetic moment at constant pressure.
From now on, we measure lengths and energies in units of and , respectively. Accordingly, the magnetic moment and the spring constant are measured in units of and , respectively. Finally, let us note that our bead-spring model is identical to the classical harmonic crystal Ashcroft and Mermin (1976), if the magnetic interactions energy and the steric repulsion energy are neglected in Eq. (1).
Iii Mapping onto pseudo-springs
Now we address the issue of particle labeling and derive an approximate elastic Hamiltonian , which can be readily used for the density functional calculation. Putting aside the magnetic and steric parts of the Hamiltonian, i.e., , which are free from this issue, see Eq. (2), we only consider the elastic part in this section. Following Ref. Cremer et al. (2017), we consider a mapping of onto a pseudo-spring potential between unlabeled particles of the form
where is a two-body interaction between particles of center-to-center distance . , in contrast to Eq. (5), involves all particle pairs . Consequently, the Hamiltonian is invariant under permutations of the particles, i.e.,
where is a permutation operator constituting the symmetric group . While such a mapping is convenient from a technical point of view, it is physically appropriate only if we can ensure through the form of that each particle in effect only interacts with a prescribed set of other particles, i.e., the same number of nearest neighbors as in the real-spring system. To this end, we will cut each spring-like interaction (which we specify later) at larger distances to prevent interactions between particles too far away from each other. The crucial point is then to ensure that does neither introduce additional contacts nor miss actually existing ones.
In one spatial dimension, the particles described by a pseudo-spring model are automatically “labeled” through the steric hard-core repulsions Tonks (1936); Percus (1976), which fixes their mutual ordering. Hence, within a proper range of the spring length at sufficiently high density Cremer et al. (2017), where the particles usually interact with their two nearest neighbors, the mapping works perfectly. Note that there is no phase transition due to the strong thermal fluctuations in one dimension. In a more realistic two- or three-dimensional fluid, however, the particles can always find a path to bypass each other. Still, a similar fixation to a cell surrounded by the nearest neighbors can be achieved via ergodicity breaking associated with the freezing transition. In this sense, we can construct a mapping of the labeled particles in the original lattice model onto the crystalline phase of the unlabeled particles with all-to-all pairwise-additive interactions. We take this viewpoint as the inversion of the following interpretation from Ref. Cates and Manoharan (2015) for a system of unlabeled hard spheres. For hard-sphere crystalline systems, even though there are possible distinct crystals of labeled particles as the Hamiltonian is invariant under permutations, all microscopic configurations correspond to one unique macroscopic lattice structure. In this way, lattice models of labeled particles, for which the factor is omitted in the statistical counting, provide a good approximation for the crystalline state of materials, emerging from the freezing of fluids with undistinguished particles. Inversely, for ferrogels, given the lattice system of labeled particles as the original model, we restore the factor by mapping the unique lattice onto possible frozen states of a fluid-like system with unlabeled particles.
In Sec. III.1 and in Fig. 2, our argument justifying such a mapping in two dimensions is laid out in full detail, before establishing the specification of the elastic pseudo potential in Sec. III.2. Then we demonstrate in Sec. III.3 how to implement the mapped system within our DFT approach. Before we proceed, we make two further remarks. First, we assume that the pseudo-spring systems governed by have a crystalline phase. Second, we consider two-dimensional crystals Zahn and Maret (2000); Keim et al. (2004) in this study. As is well known, there is no true long-range order in two-dimensional systems in the absence of truly long-range interactions Mermin and Wagner (1966).
iii.1 Detailed justification of the mapping
Our argument is based on the idea that, at low temperatures, the energetic contribution overwhelms the entropic contribution involving the factor, enabling us to map the systems of labeled particles onto those of the unlabeled ones. Such an argument is strong enough to justify the mapping, for instance, at zero temperature, for which the minimum of completely determines the equilibrium properties. At finite temperatures, however, the elastic properties (or even the stability of the systems) critically depend on the details of the mapping and of the profile of the consequent interaction potential . This does not only apply at the minimum of , but even if the adjacent particles are not located exactly at their lattice sites due to thermal fluctuations. Accordingly, a quantitative matching of lattice structures in real- and pseudo-spring systems turns out to be essential. Therefore, in this section, we carefully describe how the factor and the fluctuations of particles confined to their lattice sites in the low-temperature crystalline phase should be addressed throughout the mapping.
Let us consider the -dimensional phase space of our system. The probability density to locate our system in this phase space reads with the inverse temperature . We then partition the phase space into subspaces that are related by permutations of the particles, analogously to the symmetry-related regions discussed in Ref. Goldenfeld (1992). To this end, we first specify the completely ordered set (subspace) as . In other words, every configuration in which each particle with a smaller label is located closer to the origin than each other particle with a larger label belongs to this subspace. (If there exist pairs of particles, the distances of which from the origin are the same, one can further order the particle pairs by comparing the angle, e.g., from the -axis.) Then, with the permutation operators (for where is the identity of the symmetric group) introduced in Eq. (8), one can generate subspaces exhausting the whole phase space by permuting particles, i.e., . For example, a subspace is generated by the permutation of particles in every element (configuration) of the completely ordered subspace , with a corresponding operator .
In the case of the real-spring systems, the probability densities of each subspace are not identical to each other, see Figs. 2(a) and (b), as an exchange of any particle pair is always accompanied with a change in energy, i.e., . However, for pseudo-spring systems, the probability densities corresponding to each subspace are identical to each other, as illustrated in Figs. 2(d) and (e), because the Hamiltonian is invariant under permutations. Therefore, the probability densities in the whole phase space of the real- and pseudo-spring systems are not equivalent to each other in general. More specifically, for real-spring systems, the probability density in a subspace should deviate from the one in another subspace as they involve permutations of particles, while, for pseudo-spring systems, they are still identical due to the symmetry under permutations.
Let us now turn to low-temperature systems, i.e., . Recall that there is no phase transition for harmonic crystals Jancovici (1967); Ashcroft and Mermin (1976), that is, for our real-spring system. Nonetheless, the corresponding probability density becomes highly localized, as depicted in Fig. 2(b). In the case of the pseudo-spring system, there might occur a freezing transition instead, which is the starting point of our mapping. First we constrain ourselves to each one of the permutation-related subspaces. Therefore, particle exchanges or, namely, permutations of particles via, e.g., vacancy hopping is ignored in the following analysis. Then we consider the ergodicity breaking of the systems due to freezing. Accordingly, the probability density in each permutation-related subspace is isolated as illustrated in Fig. 2(e). [In contrast to that, as shown in Fig. 2(a) and (d), the probability densities at high temperature are rather broad for both the real- and pseudo-spring systems and, therefore, it is not possible to construct such a mapping in our situation.] In both systems, we conclude that, below a certain temperature, the thermal vibration in the particle positions are much smaller than inter-particle distances in the crystalline state. Therefore, the trajectory of each particle in each subspace does not span the full configurational space but only a localized area. That is, particle remains localized within the corresponding Wigner-Seitz cell associated with its average position, see Fig. 2(c). In such a way, we provide a tiling of , where for .
Now let us be more explicit and consider only one subspace among all permutation-related subspaces at low temperature. In the case of the real-spring system, a particle is localized in the Wigner-Seitz cell around a lattice site with the same label , so that . Then, the partition function becomes
where is the mean thermal wavelength of the particles. As indicated by the arrow on the second line, the partition function has been rewritten in terms of the Wigner-Seitz cells. With the Wigner-Seitz cells corresponding to the pseudo-spring system in the crystalline state, the partition function involving the unlabeled particles can be similarly rewritten as
where the last equality follows from the fact that is invariant under permutations. In this way, one can cancel the counting factor, completely neglecting the exchange of particles, as illustrated in Fig. 2(f). We stress that the probability of the ignored configurations in both Eq. (III.1) and Eq. (III.1) is negligible and vanishes as and/or .
Comparing now both expressions to each other, we map the lattice of labeled particles to that of unlabeled particles via
where is an arbitrary constant which does not alter any physical properties of the mapping, indicating that the (many-body) pseudo-spring potential can be determined up to an additive constant. Even though the above equation is not easy to analyze as it still involves interactions between many particles, the issue of particle labeling has been resolved: the factor does not appear in the mapping. Moreover, it provides us with a route of how to construct an approximate elastic part of the Hamiltonian in Eq. (7) for practical calculations. In addition to the trivial condition that (i) the nearest-neighbor interaction of should be reproduced by , Eq. (III.1) indicates that (ii) the lattice of the real-spring system should be recovered with the pseudo-spring potential, namely, . Next, we describe our strategy to explicitly perform such a mapping by introducing a cut-off for the springs.
iii.2 Pair potential of pseudo-springs
We consider an elastic interaction energy via pseudo-springs with a cut-off. The latter introduces an additional degree of freedom to control the offset of the remaining part of the springs, leading to a two-body interaction in the form
where and take the same role as in Eq. (5) and the two mapping parameters and are fixed by the conditions (i) and (ii), respectively, as discussed below.
Regarding condition (i), the nature of the mapping onto the pair potential leads to an apparent violation of Eq. (III.1), because the number of nearest neighbors cannot be unambiguously imposed. Even when the potential has only a finite range, additional contacts with next-to-nearest-neighbor particles can be formed in one direction, simultaneously missing the contacts with nearest-neighbor particles in the other directions. Specifically, if a particle is located, e.g., at one of the corners of its Wigner-Seitz cell, the minimal possible distance to a particle in a non-adjacent cell is only one edge length, i.e., , while the maximal possible distance to another particle in a neighboring cell is as large as . This drawback cannot be overcome by any choice of another cell shape but becomes less severe as . Hence, we simply consider an isotropic pseudo-potential of range for computational convenience. This cut-off parameter in Eq. (12) should be determined to achieve the optimal connectivity with the six nearest-neighbor particles (see the inset of Fig. 3), as the best approximation for the original bead-spring Hamiltonian.
To this end, we first determine the cut-off radius for the reference system with . As a simple analytic estimate for , one may assume a uniform (fluid-like) density and replace the seven Wigner-Seitz cells containing a particle and its six nearest neighbors by a circle of the same total area. Such an assumption leads to the value of . As we only consider the crystal, in which the density is highly inhomogeneous, however, is expected to be smaller. As a more appropriate alternative, we extract the value from MC simulations of the real-spring system with , i.e., the target of our mapping. Specifically, we compute the isotropic pair correlation function defined as
By simply probing the distance at which takes its first minimum, we find , compare Fig. 3. As the effect from the discontinuity of at (see, the red solid line in Fig. 3) is minimized with , the convergence rate turns out to be faster than the minimization with . As the results are hardly affected by this small change of the parameter and the value seems to be sufficiently accurate, we use it from now on. For systems under constant pressure with volumes other than (see Sec. IV), we use
This is the only external input necessary during the formulation of our DFT in Sec. III.3.
The condition (ii) cannot be directly imposed, since the lattice structure is usually not an input but the result of a calculation based on a prescribed interaction. The fact that this condition is not automatically fulfilled even if the real- and pseudo-spring systems are at the same density is related to another aspect of a pair-potential system in general: the lattice may be imperfect, as indicated in Fig. 2(c). The real-spring system with labeled particles assigned to the lattice is completely free from defects, whereas the mismatch between the range of the pairwise interaction and the desired lattice structure suggests that there should be some vacancies or interstitials. Moreover, if there were such defects, Eq. (III.1) should also be corrected by additional factors addressing the number of possible defect configurations. Hence, we modify such that it yields a zero vacancy concentration, , as an equivalent requirement to the above condition of an equal cell structure. The only way to do so while leaving invariant is to tune the depth of the pseudo-spring potential in Eq. (12) by an offset value , which we understand as follows. On the one hand, if , the total elastic energy is lowered by forming additional contacts with new neighbors. This undesired effect results in the undesired formation of interstitials or aggregates. On the other hand, when , vacancies are generated. Closing this section on the mapping, we are now ready to formulate our DFT. The only remaining mapping parameter will be determined within the DFT framework.
iii.3 Density functional theory
Following previous studies van Teeffelen et al. (2008); Roth et al. (2012), we consider two unit cells of a hexagonal lattice in a rectangular base with periodic boundary conditions in and directions. The volume and the number of particles of the two unit cells are denoted by and , respectively. Our starting point is to construct a grand canonical free energy functional the value of which at its minimum corresponds to the equilibrium grand potential in the grand canonical ensemble at fixed temperature , chemical potential , and volume . The free parameters employed here are the magnetic moment , the spring constant , and the packing fraction of the reference system with as defined in Sec. II. In addition, we will complete the theory, by specifying the yet to be determined offset for the pseudo-spring potential . As we have also employed an additional condition for in the previous section, however, each resultant equilibrium density profile corresponds to a parameter set of .
where the Helmholtz free energy functional consists of the ideal gas term and the excess functional, which read
respectively. Here, denotes the (irrelevant) thermal wave length. Regarding the elastic and magnetic interactions, we employ the mean-field functional in the form
while we adopt fundamental measure theory Roth et al. (2012) for . Each excess functional is calculated with the aid of the Fourier convolution theorem. The detailed forms of the Fourier transforms of and used in our calculations are described in Appendix A.
Here, the chemical potential is obtained as an output of the calculation. Specifically, the minimization process consists of two distinct stages. At the first stage, we minimize for fixed using the Picard iteration algorithm
with a mixing parameter Roth (2010), where
As we have prescribed the bulk density not the chemical potential, is updated in each iteration step to keep the average density constant. Our convergence criterion for minimization at this stage is , where . Then we determine the equilibrium density profile by further minimizing with respect to . In practice, we vary by controlling while is fixed. Comparing the values of obtained for each , we determine the vacancy concentration and the equilibrium density profile with which the free energy per particle is minimized. These procedures are repeated until we have an accuracy of for . In the calculations, we set .
Now, we determine the values of to close our DFT. Specifically, we first perform the two-stage DFT minimization as discussed above to find corresponding vacancy concentrations for given values of . Then we choose the value of for which the vacancy concentration is zero, i.e., . As described in Eq. (18), controlling in DFT calculations involves the change in . Therefore, we also use a rescaled value of given by Eq. (14) when .
Performing first the minimization with and without a hard core repulsion (), we confirm that the pseudo-spring potential indeed admits a freezing transition, which is the prerequisite of the mapping. Then we perform MC simulations with , together with , to verify the mapping at finite packing fraction, using the value of obtained from such a DFT as an input. As manifested in Fig. 3, the agreement between the real- and the pseudo-spring systems is quantitatively excellent. Therefore, we conclude that the DFT approach to systems with distinguishable particles is successfully formulated with the pseudo-spring potential between indistinguishable particles and ready to use in the presence of magnetic interactions as a model for ferrogels.
Iv Elastic properties of ferrogels
We finally demonstrate the utility of the theory to investigate the mechanical properties of ferrogels, varying the magnetic moment and the density of the particles. In particular, we probe the system at constant pressure and determine the volume and the responses to elastic deformations , where is the displacement field. With the component of the corresponding linear strain tensor Landau and Lifshitz (1986)
the stress tensor and the stiffness tensor are defined as
where is the density of the Helmholtz free energy . For two-dimensional hexagonal lattices, the stiffness tensor
can be expressed in terms of only two independent elastic constants, namely the bulk modulus and the shear modulus because of the symmetry Chaikin and Lubensky (2000).
In the following, we describe in Sec. IV.1 how to calculate the volume and elastic constants in MC simulations for the real- and pseudo-spring systems and using pseudo-potentials from our DFT treatment we compare the results of those methods in Sec. IV.2. Henceforth, and the elastic constants, namely and , are measured in units of and , respectively.
iv.1.1 MC simulation
In the case of MC simulations, we employ the isobaric-isothermal () ensemble Wood (1968); Eppenga and Frenkel (1984) and vary the rectangular lengths and independently Boal et al. (1993). Accordingly, random walks in terms of and have been performed with the detailed balance condition Frenkel and Smit (2002)
where is the transition rate corresponding to the volume changes and are the scaled coordinates defined by for . We then compute the volumes of the systems, simply taking the average of and extract and from the fluctuations , , and Boal et al. (1993). The elastic constants are computed from the real-spring as well as the pseudo-spring MC simulations to verify the validity of the mapping. In particular, for the pseudo-spring MC simulation, we use the values of the offset obtained from the density functional calculations as inputs, while the average volumes extracted from the corresponding real-spring MC simulations are employed to determine the value of the cut-off , with the aid of Eq. (14).
iv.1.2 DFT calculation
Computation of thermodynamic quantities under the given pressure in DFT is formally not straight-forward, because the theory is based on the grand-canonical ensemble. In a finite system, for example, the structure depends significantly on the specific choice of the statistical ensemble González et al. (1997, 2004); de las Heras et al. (2016). Due to the equivalence of the ensembles in the thermodynamic limit, however, one can still regard the density profiles obtained from DFT also as minima in the isobaric ensemble as long as large systems are considered. Hence, the requirement to compare the DFT results to MC simulations at constant pressure does not represent a conceptual problem for our study. Specifically, we first compute pressures at various volumes from the relation and choose the volume at which the pressure reaches the prescribed value. The accuracy of the volume is and the vacancy concentration is fixed during the procedures by fine-tuning with an accuracy of . The detailed procedures are exemplified in Appendix B.
Since we have direct access to the free energy, we can directly compute the elastic constants deforming the system. We consider four types of deformations as follows:
Specifically, we first deform the system according to the given strain, simultaneously controlling the density. Then we obtain the equilibrium density profile of the deformed system adjusting the vacancy concentration equivalently to the undeformed system as described in Sec. III.3. We note that two types of deformation are enough to determine the two unknown and . Examining four types of deformation, we verify the consistency of the theory. However, the deformations involving a volume difference seem to involve inelastic changes: shifts of the vacancy concentration in the equilibrium step indeed imply changes in the number of particles in the unit-cell, which would not occur during genuinely elastic deformations. To minimize such inelastic contributions, we utilize the identity, , which does not involve any deformations, instead of directly computing the bulk modulus from the definition of Eq. (22). Specifically, the elastic constants are calculated from
Here the pressure of the deformed systems is computed from , where is the grand potential at which the density functional is minimized under the given constraint due to a deformation. We note that, in every case, the differences in the values of are less than 10%, compared to the results for that we have calculated in analogy to Eq. (27) via the changes in , instead of using the pressure as in Eq. (26).
Before we proceed, let us make a few technical remarks. First, the minimization of the functionals at requires large computational resources since a small value of the Picard iteration parameter in Eq. (19) should be used to guarantee the convergence of the minimization. To calculate the equilibrium density in a reasonable time scale, we chose to fix the iteration parameter at in general, at the cost of loosening the strict condition to impose a zero vacancy concentration. Values of smaller than that are only used for a few cases in which the minimization eventually fails. With this constraint, we could minimize the free-energy functional up to for and , and for . Secondly, as is well known, the numerical treatment of long-range interactions with periodic boundary condition would require sophisticated techniques Allen and Tildesley (2017); Arnold and Holm (2005). In our unit-cell DFT calculations, we obviously neglect the Fourier components the wave lengths of which are longer than the unit-cell size. To bypass such complications and provide a fair comparison, with the MC simulations we simply cut the magnetic interaction beyond the nearest-neighbor interaction as an approximation Goh et al. (2018); Daddi-Moussa-Ider et al. (2019).
We explore the following sets of parameters: , each for several magnetic moments . Again, we note here that corresponds to the packing fraction of the reference system with , conveying the information of the diameter of particles. The conventional packing fraction is denoted by and is not a fixed variable as we consider the isobaric ensemble, see the paragraph below Eq. (6) for details. In Fig. 4, we first present the results at low packing fractions, i.e., and , for which the conventional hard disk system is in the fluid phase Mak (2006); Roth et al. (2012); Engel et al. (2013). Therefore, the crystallization in this low packing fraction regime is due to the elastic interaction. As one can see, the DFT and the real- and pseudo-spring MC simulations agree well with each other, except that the DFT overestimates the volumes, especially for large values of . As we adopt the mean-field functionals, the repulsion effects from the energy seem to be exaggerated, leading to the increases in volumes. Overall, the volume , the bulk modulus , and the shear modulus increase as the magnetic moment increases.
Meanwhile, results with are shown in Fig. 5. In this case, we expect strong contributions from the steric forces. We first note that a significant deviation of the theory from the MC simulations in the bulk modulus is observed. Quantitatively, the bulk moduli obtained from the DFT are approximately one fourth of those obtained from the MC simulations. This is mostly due to the fact that we had to use a quite large value, i.e., , for the vacancy concentration because of the technical reasons related to the computational time and the consequent choice of the mixing parameter (see Sec. IV.1). Indeed, using the same vacancy concentration, we also obtain similar deviations in at or . Moreover, a good agreement between the real- and pseudo-spring MC simulations is still observed, confirming that such a deviation of the DFT does not indicate a failure of our mapping. Surprisingly, the behavior of is still well predicted by the DFT. In contrast to , computation of does not involve any changes in volume and consequently both the density and the vacancy concentration remain approximately constant during deformation. (In contrast to that, deformations involving volume changes are always accompanied by shifts of the vacancy concentration in equilibrium.) Therefore, we speculate that the vacancy concentration plays a much smaller role in the case of the volume-conserving deformation in DFT, for which remains basically the same. To conclude, except for the technical issue discussed, the DFT provides qualitatively correct trends even at high packing fractions.
Remarkably, we observe new response behaviors of and to an increase in for the high packing fraction, which are opposite to the low packing fraction. For , and decrease as increases, while is still an increasing function of . Emergence of such two different scenarios originates from the composite nature of the ferrogels: here, we are observing different types of crossovers from the hard-disk crystals (high packing fraction, small ) or the harmonic crystals (low packing fraction, small ) to the elastic-dipolar crystals (low packing fraction, large ).
To understand the phenomena in detail, we first note that always increases as increases because the magnetic interaction is repulsive. Then for low packing fractions, the magnetic repulsion simply causes additional increases in the bulk and shear moduli on top of the harmonic crystals. In contrast to that, for high packing fractions, the steric force dominates the mechanical properties and the system is very stiff (hard-disk crystals) with relatively large and for small values of . As the magnetic moment and consequently the volume increase, the packing fraction decreases and the contributions from the steric repulsion become insignificant compared to the dipolar repulsion at some values of the packing fraction around [see, e.g., the inset of Fig. 5(c)]. These values are slightly smaller than the fluid-crystal coexisting packing fractions of the hard-disk systems, which are in between 0.68 and 0.73 Bernard and Krauth (2011); Engel et al. (2013); Roth et al. (2012); Thorneywork et al. (2017); Lin and Oettel (2018). Once has decreased enough, as in the low packing fraction regime, the signature of the elastic-dipolar crystals should be recovered, the elastic properties of which are governed by the combination of the spring and dipolar interactions. Indeed for , the volume and the elastic constants (compare Figs. 4 and 5) are quantitatively similar among =0.3, 0.5, and 0.8. As expected, the values of and for increase in the regime of large magnetic moments, i.e., for where ( 0.68 and 0.60 for 5 and 8, respectively).
V Summary and outlook
In this study, we have formulated a DFT for a two-dimensional ferrogel model. We have replaced the labeled particles in a state of strictly permanently connected neighboring particles by the unlabeled particles in a fluid-like state to map the elastic part of the associated energy onto a pseudo-potential invariant under permutations. In particular, we have shown that the mapping provides a plausible approximation for the considered systems and their response to magnetic interactions, even though the mapping still leads to some deviations in the calculated response of the systems. These deviations have been minimized by fine-tuning the mapping parameters. Lastly, it has been demonstrated that the elastic properties of ferrogels can be successfully investigated in this framework and two scenarios have been identified for the response mechanism of the dipole-spring system, depending on the packing fraction. Notably, our DFT approach may also provide a clue for the scale-bridging Menzel (2014) between mesoscopic dipole-spring models and the macroscopic description of ferrogels Jarkova et al. (2003); Bohlius et al. (2004).
A strong feature of the presented DFT approach relying on the pseudo-spring model in two dimensions is that it works well independently of the density. First of all, the mapping is generally well justified even in the low-density regime. Moreover, the mean-field DFT, which we employed here, provides a good approximation even in the high-density regime where we found good agreement between our theoretical treatment and MC simulations, especially for the shear modulus. In contrast to that, for a one-dimensional model, the prediction of unphysical freezing within mean-field DFT restricts the parameter space to intermediate densities and weaker elasticities Cremer et al. (2017), where the pseudo-spring approximation itself is less accurate due to large possible gaps between neighboring particles. We stress that, in general, mean-field approximations become more accurate with increasing dimensionality. Besides, true long-range order generally exists in three dimensions. Therefore, we expect our DFT approach to work even better in three dimensions.
Regarding possible experiments to test our theoretically predicted scenarios, the high packing fraction of needs to be discussed. As mentioned in Ref. Huang et al. (2016), for certain systems there exists a region of increased mechanical stiffness around the magnetic particles. Effectively, this could be equivalent to an increased size of the particles. In addition, we note that in three dimensions, the fluid-crystal coexisting densities are reported as low as Hoover and Ree (1968); Roth et al. (2002). Indeed, volume concentrations of approximately 50% have already been reported in three-dimensional samples of ferrogels Liu et al. (2008), implying that the steric force should be taken into account explicitly. Therefore, decreasing elastic moduli in response to an increase of the magnitude of the magnetic moments might be observed in real three-dimensional systems. Furthermore, the pressure value employed in this study seems to be small: for instance, Pa for the system studied in Ref. Messing et al. (2011) with nm at room temperature. Therefore, for direct comparison with experiments, a broad range of pressures should be explored in future studies with more realistic settings.
There are several remaining issues. First, aligning the magnetic dipoles within the two-dimensional plane and the resulting in-plane anisotropy will give rise to additional phenomena not observed here for perpendicular dipoles. To list a few, the aspect ratio of width to height of the unit cell will vary depending on the in-plane orientation of the magnetic dipole moments, the volume may change and the magnetic particles could touch each other, leading to an abrupt change in elastic moduli even at a relatively low packing fraction Annunziata et al. (2013). As mentioned above, an extension to three dimensions would even strengthen the correspondence to real materials. Additionally, the mean-field functionals can be replaced by more sophisticated ones, especially to consider the long-range nature of the magnetic interaction van Teeffelen et al. (2008). Apart from that, the response dynamics to external magnetic fields represents another topic of interest. Dynamical DFT Marconi and Tarazona (1999); Archer and Evans (2004); Stopper et al. (2018) should be the obvious candidate to study these phenomena in the present context.
We thank Urs Zimmermann and Christian Hoell for providing numerical codes that were useful for the initiation of this study. We also thank Shang-Chun Lin, Martin Oettel, Benno Liebchen, Mehrdad Golkia, Saswati Ganguly, and Robert Evans for helpful discussions. This work was supported by funding from the Alexander von Humboldt-Stiftung (S.G.) and from the Deutsche Forschungsgemeinschaft through the SPP 1681, grant nos. ME 3571/3 (A.M.M) and LO 418/16 (H.L.).
Appendix A Fourier transform
which are numerically implemented with the aid of the convolution theorem. First, the Fourier transform of the elastic energy can be computed as follows:
For the magnetic interaction, the Fourier transformation can be performed as follows:
where is the generalized hypergeometric function. For ,
We finally obtain
The generalized hypergeometric functions were implemented using the Arb library Johansson (2017).
Appendix B DFT minimization
In this section, we describe the procedure of computing the equilibrium profile at prescribed pressure within the grand-canonical DFT, which is exemplified in Table 1. First, the volume and the offset are fixed (the first and the second column in Table 1. We then compute the free energy per particle , varying the vacancy concentration (the third column). At this stage, we control by changing while fixing the density , see Eq. (18). We choose the value of where is minimized at the prescribed value of the vacancy concentration (for the parameter set in Table 1, ). At this stage, the pressure corresponding to the fixed volume is determined simultaneously (the fourth column). Now, we change the volume, varying the density while fixing the number of particles , and repeat the above procedure to calculate corresponding pressures. Finally, we obtain the volume and corresponding density profile, comparing pressures with the prescribed pressure value, i.e., .
- Pathria (1996) R. K. Pathria, Statistical Mechanics (Butterworth-Heinemann, Boston, 1996) pp. 22–26.
- Kardar (2007) M. Kardar, Statistical Physics of Particles (Cambridge University Press, 2007).
- Huang (2009) K. Huang, Introduction to Statistical Physics (Chapman and Hall/CRC, 2009).
- Cates and Manoharan (2015) M. E. Cates and V. N. Manoharan, Soft Matter 11, 6538 (2015).
- Sethna (2006) J. Sethna, Statistical Mechanics: Entropy, Order Parameters, and Complexity, Vol. 14 (Oxford University Press, 2006).
- Swendsen (2008) R. H. Swendsen, Entropy 10, 15 (2008).
- Frenkel (2014) D. Frenkel, Mol. Phys. 112, 2325 (2014).
- Jaynes (1957) E. T. Jaynes, Phys. Rev. 106, 620 (1957).
- Jaynes (1992) E. T. Jaynes, “The Gibbs paradox,” in Maximum Entropy and Bayesian Methods (Springer, 1992) pp. 1–21.
- Doi and Edwards (1986) M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, Oxford, 1986).
- Menzel (2019) A. M. Menzel, Arch. Appl. Mech. 89, 17 (2019).
- Cremer et al. (2017) P. Cremer, M. Heinen, A. M. Menzel, and H. Löwen, J. Phys.: Condens. Matter 29, 275102 (2017).
- Schindler et al. (2019) T. Schindler, R. Wittmann, and J. M. Brader, Phys. Rev. E 99, 012605 (2019).
- Filipcsei et al. (2007) G. Filipcsei, I. Csetneki, A. Szilágyi, and M. Zrínyi, Adv. Polym. Sci. 206, 137 (2007).
- Ilg (2013) P. Ilg, Soft Matter 9, 3465 (2013).
- Gundermann and Odenbach (2014) T. Gundermann and S. Odenbach, Smart Mater. Struct. 23, 105013 (2014).
- Landers et al. (2015) J. Landers, L. Roeder, S. Salamon, A. M. Schmidt, and H. Wende, J. Phys. Chem. C 119, 20642 (2015).
- Frickel et al. (2011) N. Frickel, R. Messing, and A. M. Schmidt, J. Mater. Chem. 21, 8466 (2011).
- Messing et al. (2011) R. Messing, N. Frickel, L. Belkoura, R. Strey, H. Rahn, S. Odenbach, and A. M. Schmidt, Macromolecules 44, 2990 (2011).
- Roeder et al. (2015) L. Roeder, P. Bender, M. Kundt, A. Tschöpe, and A. M. Schmidt, Phys. Chem. Chem. Phys. 17, 1290 (2015).
- Stepanov et al. (2008) G. V. Stepanov, D. Y. Borin, Y. L. Raikher, P. V. Melenev, and N. S. Perov, J. Phys.: Condens. Matter 20, 204121 (2008).
- Stolbov et al. (2011) O. V. Stolbov, Y. L. Raikher, and M. Balasoiu, Soft Matter 7, 8484 (2011).
- Stolbov and Raikher (2019) O. V. Stolbov and Y. L. Raikher, Arch. Appl. Mech. 89, 63 (2019).
- Pessot et al. (2016) G. Pessot, H. Löwen, and A. M. Menzel, J. Chem. Phys. 145, 104904 (2016).
- Pessot et al. (2018) G. Pessot, M. Schümann, T. Gundermann, S. Odenbach, H. Löwen, and A. M. Menzel, J. Phys.: Condens. Matter 30, 125101 (2018).
- Goh et al. (2018) S. Goh, A. M. Menzel, and H. Löwen, Phys. Chem. Chem. Phys. 20, 15037 (2018).
- Biller et al. (2014) A. M. Biller, O. V. Stolbov, and Y. L. Raikher, J. Appl. Phys. 116, 114904 (2014).
- Biller et al. (2015) A. M. Biller, O. V. Stolbov, and Y. L. Raikher, Phys. Rev. E 92, 023202 (2015).
- Cremer et al. (2015) P. Cremer, H. Löwen, and A. M. Menzel, Appl. Phys. Lett. 107, 171903 (2015).
- Cremer et al. (2016) P. Cremer, H. Löwen, and A. M. Menzel, Phys. Chem. Chem. Phys. 18, 26670 (2016).
- Puljiz et al. (2016) M. Puljiz, S. Huang, G. K. Auernhammer, and A. M. Menzel, Phys. Rev. Lett. 117, 238003 (2016).
- Puljiz and Menzel (2017) M. Puljiz and A. M. Menzel, Phys. Rev. E 95, 053002 (2017).
- Weeber et al. (2012) R. Weeber, S. Kantorovich, and C. Holm, Soft Matter 8, 9923 (2012).
- Weeber et al. (2015) R. Weeber, S. Kantorovich, and C. Holm, J. Magn. Magn. Mater. 383, 262 (2015).
- Weeber et al. (2019) R. Weeber, P. Kreissl, and C. Holm, Arch. Appl. Mech. 89, 3 (2019).
- Goldenfeld (1992) N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group (Addison-Wesley, Urbana-Champaign, 1992) Chap. 2.
- Ashcroft and Mermin (1976) N. Ashcroft and N. Mermin, Solid State Physics (Holt, Rinehart and Winston, 1976).
- Evans (1979) R. Evans, Adv. Phys. 28, 143 (1979).
- Lutsko (2010) J. F. Lutsko, “Recent developments in classical density functional theory,” in Adv. Chem. Phys. (John Wiley & Sons, Ltd, 2010) Chap. 1, pp. 1–92.
- Evans et al. (2016) R. Evans, M. Oettel, R. Roth, and G. Kahl, J. Phys.: Condens. Matter 28, 240401 (2016).
- Löwen (2002) H. Löwen, J. Phys.: Condens. Matter 14, 11897 (2002).
- Ebner et al. (1976) C. Ebner, W. F. Saam, and D. Stroud, Phys. Rev. A 14, 2264 (1976).
- Ramakrishnan and Yussouff (1979) T. V. Ramakrishnan and M. Yussouff, Phys. Rev. B 19, 2775 (1979).
- Roth (2010) R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010).
- Wittmann et al. (2016) R. Wittmann, M. Marechal, and K. Mecke, J. Phys.: Condens. Matter 28, 244003 (2016).
- Roth et al. (2012) R. Roth, K. Mecke, and M. Oettel, J. Chem. Phys. 136, 081101 (2012).
- Wittmann et al. (2017) R. Wittmann, C. E. Sitta, F. Smallenburg, and H. Löwen, J. Chem. Phys. 147, 134908 (2017).
- Lin and Oettel (2018) S.-C. Lin and M. Oettel, Phys. Rev. E 98, 012608 (2018).
- Zimmermann et al. (2016) U. Zimmermann, F. Smallenburg, and H. Löwen, J. Phys.: Condens. Matter 28, 244019 (2016).
- Roth and Gillespie (2016) R. Roth and D. Gillespie, J. Phys.: Condens. Matter 28, 244006 (2016).
- Archer and Evans (2004) A. J. Archer and R. Evans, J. Chem. Phys. 121, 4246 (2004).
- Bernard et al. (2009) E. P. Bernard, W. Krauth, and D. B. Wilson, Phys. Rev. E 80, 056704 (2009).
- Engel et al. (2013) M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P. Bernard, and W. Krauth, Phys. Rev. E 87, 042134 (2013).
- Tonks (1936) L. Tonks, Phys. Rev. 50, 955 (1936).
- Percus (1976) J. K. Percus, J. Stat. Phys. 15, 505 (1976).
- Zahn and Maret (2000) K. Zahn and G. Maret, Phys. Rev. Lett. 85, 3656 (2000).
- Keim et al. (2004) P. Keim, G. Maret, U. Herz, and H. H. von Grünberg, Phys. Rev. Lett. 92, 215504 (2004).
- Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
- Jancovici (1967) B. Jancovici, Phys. Rev. Lett. 19, 20 (1967).
- van Teeffelen et al. (2008) S. van Teeffelen, H. Löwen, and C. N. Likos, J. Phys.: Condens. Matter 20, 404217 (2008).
- Oettel et al. (2010) M. Oettel, S. Görig, A. Härtel, H. Löwen, M. Radu, and T. Schilling, Phys. Rev. E 82, 051404 (2010).
- Landau and Lifshitz (1986) L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Elsevier, Oxford, 1986).
- Chaikin and Lubensky (2000) P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics (Cambridge University Press, Cambridge, 2000).
- Wood (1968) W. W. Wood, J. Chem. Phys. 48, 415 (1968).
- Eppenga and Frenkel (1984) R. Eppenga and D. Frenkel, Mol. Phys. 52, 1303 (1984).
- Boal et al. (1993) D. H. Boal, U. Seifert, and J. C. Shillcock, Phys. Rev. E 48, 4274 (1993).
- Frenkel and Smit (2002) D. Frenkel and B. Smit, Understanding Molecular Simulation, 2nd ed. (Academic Press, San Diego, 2002).
- González et al. (1997) A. González, J. A. White, F. L. Román, S. Velasco, and R. Evans, Phys. Rev. Lett. 79, 2466 (1997).
- González et al. (2004) A. González, J. A. White, F. L. Román, and S. Velasco, J. Chem. Phys. 120, 10634 (2004).
- de las Heras et al. (2016) D. de las Heras, J. M. Brader, A. Fortini, and M. Schmidt, J. Phys.: Condens. Matter 28, 244024 (2016).
- Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, 2017).
- Arnold and Holm (2005) A. Arnold and C. Holm, “Efficient methods to compute long-range interactions for soft matter systems,” in Advanced Computer Simulation Approaches for Soft Matter Sciences II, edited by C. Holm and K. Kremer (Springer, Berlin Heidelberg, 2005) pp. 59–109.
- Daddi-Moussa-Ider et al. (2019) A. Daddi-Moussa-Ider, S. Goh, B. Liebchen, C. Hoell, A. J. T. M. Mathijssen, F. Guzmán-Lastra, C. Scholz, A. M. Menzel, and H. Löwen, J. Chem. Phys. 150, 064906 (2019).
- Mak (2006) C. H. Mak, Phys. Rev. E 73, 065104(R) (2006).
- Bernard and Krauth (2011) E. P. Bernard and W. Krauth, Phys. Rev. Lett. 107, 155704 (2011).
- Thorneywork et al. (2017) A. L. Thorneywork, J. L. Abbott, D. G. A. L. Aarts, and R. P. A. Dullens, Phys. Rev. Lett. 118, 158001 (2017).
- Menzel (2014) A. M. Menzel, J. Chem. Phys. 141, 194907 (2014).
- Jarkova et al. (2003) E. Jarkova, H. Pleiner, H.-W. Müller, and H. R. Brand, Phys. Rev. E 68, 041706 (2003).
- Bohlius et al. (2004) S. Bohlius, H. R. Brand, and H. Pleiner, Phys. Rev. E 70, 061411 (2004).
- Huang et al. (2016) S. Huang, G. Pessot, P. Cremer, R. Weeber, C. Holm, J. Nowak, S. Odenbach, A. M. Menzel, and G. K. Auernhammer, Soft Matter 12, 228 (2016).
- Hoover and Ree (1968) W. G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 (1968).
- Roth et al. (2002) R. Roth, R. Evans, A. Lang, and G. Kahl, J. Phys.: Condens. Matter 14, 12063 (2002).
- Liu et al. (2008) T.-Y. Liu, S.-H. Hu, K.-H. Liu, D.-M. Liu, and S.-Y. Chen, J. Control. Release 126, 228 (2008).
- Annunziata et al. (2013) M. A. Annunziata, A. M. Menzel, and H. Löwen, J. Chem. Phys. 138, 204906 (2013).
- Marconi and Tarazona (1999) U. M. B. Marconi and P. Tarazona, J. Chem. Phys. 110, 8032 (1999).
- Stopper et al. (2018) D. Stopper, A. L. Thorneywork, R. P. A. Dullens, and R. Roth, J. Chem. Phys. 148, 104501 (2018).
- Johansson (2017) F. Johansson, IEEE Trans. Comput. 66, 1281 (2017).