Inverse patchy colloids with small patches: fluid structure and dynamical slowing down

Inverse Patchy Colloids (IPCs) differ from conventional patchy particles because their patches repel (rather than attract) each other and attract (rather than repel) the part of the colloidal surface that is free of patches. These particular features occur, .e.g., in heterogeneously charged colloidal systems. Here we consider overall neutral IPCs carrying two, relatively small, polar patches. Previous studies of the same model under planar confinement have evidenced the formation of branched, disordered aggregates composed of ring-like structures. We investigate here the bulk behavior of the system via molecular dynamics simulations, focusing on both the structure and the dynamics of the fluid phase in a wide region of the phase diagram. Additionally, the simulation results for the static observables are compared to the Associative Percus Yevick solution of an integral equation approach based on the multi-density Ornstein-Zernike theory. A good agreement between theoretical and numerical quantities is observed even in the region of the phase diagram where the slowing down of the dynamics occurs.


Introduction
In the last ten years, patchy colloidal particles [1,2] have emerged as one of the most promising solutions to the problem of designing new materials via a bottom-up strategy. Due to the specificity of their bonding patterns, patchy particles have indeed proven to be perfect candidates as nano-and micro-scale building blocks for materials with desired symmetries and physical properties [3]. The interest in these systems is triggered by the wide range of potential applications, e.g. in photonics, biomimetic materials synthesis or electronics, and relies on the impressive advances in the production of colloidal particles with tunable chemical/physical surface designs [1]. Even though experimental studies on assembly processes of patchy particles are so far rather rare [4][5][6][7], there is no doubt about the great potentialities of these colloids as self-assembling units of completely new macroscopic structures. In fact, numerical and theoretical studies have opened the way to a surprising multitude of both ordered and disordered phases [2] that significantly extends the many possibilities already offered by isotropically interacting colloidal particles.
A further direction of investigation in the field of patchy particles has been recently put forward by the introduction of patchy systems with charged surface regions; these systems are also referred to as Inverse Patchy Colloids (IPCs) [8]. While conventional patchy particles are built by adding attractive patches on the top of otherwise repulsive spherical particles, inverse patchy particles have patches that repel each other, while they attract the rest of the colloidal surface that is free of patches. IPC model systems are thus suitable to study self-assembling units with a non-homogeneous surface charge distribution. Within this broad class of systems, we focus here on heterogeneously charged colloids with two charged polar patches and an oppositely charged equatorial belt. The analytical description of the underlying microscopic system has been provided in Ref. [8]; it features a complex expression in terms of Bessel functions and Legendre polynomials for the pair potential between two of such IPCs. The related coarsegrained description, also presented in Ref. [8], allows to deal with a much simpler model that maintains the key physical features of the microscopic systems. The model was originally developed to describe the interaction between complex units emerging from the adsorption of charged polyelectrolyte stars onto the surface of oppositely charged colloids, but it works equally well to describe heterogeneously charged particles synthesized in the lab [9]. The only prerequisite of the applicability of our description -analytical as well as coarse-grained -is that the center of charge of a given surface region does not lie outside the colloidal particle.
The equilibrium behavior of IPCs is dominated by a non-trivial interplay between attractive and repulsive directional interactions. A very recent example of exotic structures emerging in IPC systems in the bulk has been reported in Ref. [12]. There, slightly overall charged IPCs with two, relatively extended and long-ranged, polar patches were selected and their phase diagram was evaluated by means of Monte Carlo simulations and free energy calculations [10,11]. In addition to the fluid phase and to the fcc (ordered and plastic) solids, the resulting phase diagram features a wide region in temperature and density where a crystal formed by parallel monolayers is stable; the stability of this laminar phase was also investigated on changing the interaction parameters [13]. The effect of the interaction range and of the patch extension was also studied in the bulk fluid phase: the structural and thermodynamic properties of a selection of overall neutral IPC systems have been studied by means of Monte Carlo simulations and theoretical approaches. In particular, the multi-density Ornstein-Zernike integral equation, supplemented with the Associative Percus Yevick (APY) approximated closure [15][16][17][18], has been extended to describe IPCs with an arbitrary number of patches. From the comparison between the simulation data and other theoretical approaches, such as the Reference Hypernetted Chain description [19,20] and the Barker-Henderson thermodynamic perturbation theory [21,22], it emerged that the APY approach is a better candidate to describe the IPC class of systems, performing a higher accuracy in the results and a wider convergence ability when lowering the system temperature.
In the present paper, we focus on overall neutral IPCs with relatively small patches and a short interaction range. Under planar confinement conditions and at intermediate densities, the selected system did not form ordered structures, but rather disordered, branched aggregates composed of particle rings [23]. In contrast, IPCs with the same interaction range and particle charge, but bigger patch size tend to form planar colloidal monolayers with a broad variety of translational and orientational ordering [23,24].
In this contribution, we consider the relatively small patch case and we focus on its bulk behavior in the fluid phase, while the bulk counterpart of the other models investigated under planar confinement [23,24] will be addressed in future publications. We use Molecular Dynamics (MD) simulations to study a selected IPC system, focusing on the structural and dynamical properties of the fluid phase in a wide range of temperatures and densities, including the region where the dynamical slowing down occurs. Additionally, we compare the numerical results for the static properties at different state points with the APY theoretical description.
Our goal is twofold: on one side we want to check how far the APY theory can be extended to the high-density and/or low-temperature regions; on the other side we want to study the behavior of the system in the region where the dynamical slowing down occurs.
The paper is organized as follows. In section 2.1 we introduce the IPC model. In section 2.2 we describe the investigation methods, with particular attention to the way the interaction potential was implemented in the MD simulation code. In section 3 we describe our results, namely the static properties are discussed in section 3.1, while the dynamic ones are reported in section 3.2. Our final considerations are presented in section 4.

Model
Our two-patch IPCs are modeled as hard spheres of radius σ c carrying three interaction sites: one corresponding to the colloid center and two corresponding to the patch centers; the latter are placed at a distance e from the first, in opposite directions (see Figure 1). The bare colloid interaction sphere has radius σ c + δ/2, where δ is the resulting interaction range, while both the patch interaction spheres have radius σ p ; 2σ c is assumed to be the unit of length. By construction the following relations hold [8] σ c + δ/2 = e + σ p and cos γ = where γ is half of the patch opening angle. By virtue of these constraints, the model is characterized only by two geometrical parameters: the interaction range δ and the patch extension γ.
Beyond the hard core repulsion, the pair potential between two IPCs at distance r is given by [ where i (j) specifies one of the three interaction sites of the first (second) IPC, w ij is the overlap volume of the corresponding interaction spheres, and u ij is the energy strength of the ij interaction. We note that, while the u ij are constants, the w ij -as well as the potential U -depend on both the inter-particle distance and the relative orientation of the two IPCs; we drop the arguments for brevity. For two identical patches, the set of the u ij reduces to three independent energy constants: u cc , u pp and u cp . The coarse-grained model parameters δ, γ andū = (u cc , u pp , u cp ) for a specific microscopic system can be fixed by taking advantage of the Debye-Hückel analytical description of the underlying model [8].
In the present contribution, we consider IPCs under relatively high screening conditions, namely we fix κσ c = 5, where κ −1 is the Debye screening length, and κδ = 2; the corresponding interaction range is δ = 0.4σ c . The polar patches are chosen to be relatively small, namely γ ≈ 30 • , corresponding to the choice e = 0.64σ c and σ p = 0.56σ c . The arrayū results from considering overall neutral colloids of diameter 60nm, dispersed in water at room temperature. The energy minimum, obtained when two particles at contact are in a T-shape configuration, sets the unit of energy m . These features lead toū = (0.8631, 241.8, −26.62) in units of m . In the following we use reduced units for the temperature, T * = k B T / m , the energy per particle, u * = U/(N m ), the chemical potential µ * = µ/ m and the density, ρ * = ρ(2σ c ) 3 .

Methods
The bulk fluid phase of our IPC system was investigated with MD simulations in the NVE ensemble using the velocity Verlet integration scheme [26] and the RATTLE algorithm for geometric constraints [27]. Since the conventional RATTLE procedure is singular for linear particles [28], we describe the patch-center-patch complex reducing the number of degrees of freedom according to the scheme introduced in Ref. [28]. We proceed as follows: the two patches are treated as two particles separated by a fixed distance and move according to effective forces that include both the inertia of the central particle and the forces acting on it. As a consequence of these effective forces, the central particle automatically satisfies the constraint of being in the middle point of the line joining the two patches and its trajectory is automatically provided by the knowledge of the patch trajectories. This approach reduces the number of equations of motion from 18N to 12N for a system of N IPCs.
In order to integrate the equations of motion, we approximated the hard core repulsion via a continuous harshly repulsive soft sphere interaction with k = 15 and A = 500. We checked the consistency between the continuous and the hard core [8] versions of the model by comparing energies, pair distribution functions and static structure factors at several state points with those obtained via Monte Carlo simulations of the original model [23]. No significative differences were observed. Since in the microcanonical ensemble the temperature fluctuates, our temperature values are calculated as time averages over the kinetic temperature.
The algorithm requires to specify a mass for each interaction center. We chose to assign the same mass, m 0 , to the central particle as well as to each of the two patches. This choice fixes the total mass of an IPC to 3m 0 and its inertia tensor is thus I = diag(2m 0 e 2 , 2m 0 e 2 , 0). Defining m 0 as the mass unit also fixes the time unit, i.e.
The timestep of the simulation was chosen to be 2.5×10 −3 in units of t 0 ; this choice ensured that the total energy was conserved within 10 −2 % of its mean value.
Most results shown in this contribution were obtained with a system of 2048 particles; on few occasions ensembles with 4000 were considered. To check for possible size effects, simulations for ρ = 0.40 state points were carried out for both ensemble sizes. All observables agreed within numerical accuracy: for instance the internal energy showed differences only in the third digit. The systems were equilibrated for over 10 5 t 0 , and in the smallest-density case for even over 10 6 t 0 .
The gas-liquid critical point of our system was evaluated via grand-canonical Monte Carlo (GCMC) simulations and using the histogram reweighting technique [29]. We consider a cubic box with side L = 11 and looked for values of temperature (T ) and chemical potential (µ) at which the system exhibited large fluctuations in both the number of particles (N ) and the total internal energy of the system (U ). At the state points where such fluctuations were found, we started ten independent MC runs -each one extending over more than 10 6 MC steps -and evaluated the probability distribution of the order parameter M ∼ N + sU , s being the mixing parameter [29]. The precise location of the critical point was then obtained by using the histogram reweighting technique [29] and searching for those (T, µ, s) values at which the probability distribution of M matched the one of the 3D Ising model [29].
The theoretical description of the static properties of the model was carried out using the multi-density integral-equation theory developed in Ref. [14]. The pair distribution functions, the structure factors and the internal energy were calculated using numerical solution of the corresponding multi-density Ornstein-Zernike equation supplemented by the APY closure relations. To avoid unnecessary repetition we refer the reader to the original publication [14], where the theory is described in details. In contrast to Ref. [14] where the reference model has always a HS contribution, here we consider both a HS and a soft sphere core; the latter is the same used in MD simulations. In the following, theoretical results for the soft sphere case are reported only if we observed significative differences with respect to the HS case.

Results
In an effort to obtain a better overview over the system behavior we have first located the critical point of the gas-liquid transition via GCMC simulations, using the histogram reweighting technique [29]; this point was found to be located at (T * c , µ * c , ρ * c ) = (0.122, −0.477, 0.266). In matching the distribution of the order parameter M to the binodal curve of the Ising model, the best fit was achieved for s ≈ 0.
MD simulations were carried out on a grid of state points, defined by six densities (ρ * = 0.1, 0.25, 0.4, 0.55, 0.7 and 0.8) and six temperature values (T * = 0.10, 0.11, 0.12, 0.13, 0.15 and 0.30); ρ * = 0.8 was the highest density value for which we could equilibrate the system in the liquid phase. Temperatures were chosen to span from super-critical to sub-critical values in such a way that the isotherms better sample the region around the critical point.

Static structure and thermodynamics
To describe the static structure of the system we have calculated the pair distribution function, g(r), defined by [30] and the static structure factor [30] where r is the distance between two IPC centers, r i with i = 1, · · · , N represents the position of particle i and k is the wavevector. S(k) and g(r) were calculated directly during the simulation, i.e., "on the fly".
While in the simulations g(r) and S(k) are obtained according to the aforementioned definitions, within the APY theory g(r) is obtained from the Fourier transform of S(k).
In Figure 2 data for g(r) at the selected state points are reported. Each panel refers to a different thermodynamic state: the left column corresponds to low density cases, i.e. ρ * = 0.10, the central column displays data at an intermediate density, i.e. ρ * = 0.40, and the right column shows results for the highest density state points, i.e. ρ * = 0.80; temperature decreases from the top to the bottom, dropping from T * = 0.30 to T * = 0.10. Since we report here only equilibrium data, two low-temperature panels are missing in Figure 2 because they correspond to state points located within the binodal.
Throughout, a good, sometimes even excellent, agreement between theoretical and simulation results can be observed; this holds in particular for distances covering the first neighbor shell, i.e. for r ≤ 2σ c + δ. Since the APY approach -as all integral equation approaches -allows to calculate the thermodynamic properties by integrating g(r) over the interaction range of the potential, a reliable evaluation of the pair distribution function within the first neighbor shell is indispensable for an accurate calculation of these properties. In addition, a reliable prediction of the contact value of the pair distribution function, namely g(2σ c ), is responsible for a reliable estimate of the pressure in systems with hard-core interactions. At low densities and high temperatures (upper-left panel) the numerical g(r) reproduces the characteristic features of a hard sphere system, characterized by a pronounced peak at contact, i.e., for r = 2σ c , and a fast, structureless decay to 1 on increasing r. Upon lowering the temperature and/or increasing the density, a second nearest neighbors peak gradually emerges at r = 4σ c . The overimposed theoretical data show that APY predictions are in good agreement with the simulation results not only at contact, but also around the second peak of g(r). Only for the high-density, low-temperature state point (bottom-right panel) small discrepancies can be observed: the theoretical data are not able to reproduce the complex shape of the peak at r 4σ c . In subsection 3.2 we will provide evidence that at this state point the system dynamically slows down.
In Figure 3 we show the corresponding S(k), using the same arrangement of panels as in Figure 2. Similar to the pair distribution function the agreement between theory and simulations is extremely good, even though data are slightly noisy due to substantial differences in the statistics. While at high temperatures and/or high densities the theoretical and numerical data fall essentially on the top of each other, at low temperatures and low/intermediate densities slight discrepancies can be observed at k = 0. The value of the static structure factor at k = 0 is of particular relevance since S(0)/ρk B T is the isothermal compressibility. An increase of S(0) indicates the onset of a liquid-gas phase separation. In our case at T * = 0.15, i. e. at a temperature slightly above the gas-liquid transition temperature, we observe the onset of an increase in S(0) in the simulation data; the theoretical approach is able to reproduce this increase of the compressibility (with slight discrepancies). In contrast, far from the gas-liquid phase separation region, theoretical predictions are very accurate, even at very high density.
Finally, in Figure 4 we display the internal energy per particle versus the density of the system along four isotherms; data along the lowest temperature isotherm are incomplete in a density range where the system undergoes phase separation. At the state points where APY converges, the agreement between the theoretical and numerical data is good for super-critical temperatures, while it is less convincing for sub-critical temperatures. We also report APY data for the soft sphere core. A small difference can be observed between the two reference models: the soft sphere provides a slightly more accurate agreement only at high temperatures, while at low temperatures the two theoretical sets of data reproduce simulation results with the same accuracy.

Dynamics
We study the dynamics of the system exclusively via simulations by calculating the autocorrelation functions of two particular observables as well as the diffusion coefficient at several state points.
To be more specific, we have calculated the autocorrelation function of the particle velocities v(t), defined as as well as the autocorrelation function of the particle orientation C n (t) = n(t) · n(0) , with n(t) being a unit vector pointing from one patch to the other. Both functions are single particle correlation functions; thus the angular brackets denote a time average along the simulation as well as an average over all particles. Both correlation functions are shown in Figure 5 for three selected densities (considering a low, an intermediate and a high ρ-value) covering a range of temperatures from 0.3 down to 0.1.
C v (t) features the expected behavior: (i) a monotonous decay at low/intermediate densities and high temperatures, and (ii) a characteristic minimum, signaling a rebound mechanism, as the density increases and/or the temperature decreases.
In contrast, C n (t) shows a rather complex, strongly state-dependent behavior. At low temperatures and for all densities this correlation function is slowly decaying, indicating a strong persistence in the orientational motion of the particles. As the temperature increases the decay of the correlation function is still monotonous, indicating that the particles prefer to maintain their original orientation; C n (t) drops more rapidly with time as the density decreases. Only at rather high temperatures and low/intermediate densities, the correlation function shows a pronounced minimum; this feature signals that the particles start to freely rotate in space.
Another feature of C n (t) deserves particular attention, even though a thorough discussion of this phenomenon requires a much more comprehensive and systematic investigation that we postpone to a future contribution [25]. In particular for low temperature states the orientational correlation function shows characteristic shortrange oscillations in time, whose frequency increases with density. This feature stems from fast oscillations of the particle orientation around a preferred spatial axis and is undoubtedly a finger-print of the strong coupling between the particles which occurs preferentially at low temperatures (where particles are tightly bonded) and/or high densities (where the average distance between particles is rather small). Also the fact that the frequency of the spatial oscillations decreases with increasing density can be understood via a stronger spatial coupling of the particles as ρ grows.
Finally, we determined the diffusion coefficient, D, at all the state points investigated. This quantity is calculated via the mean square displacement of the particles, using the well-known relation at sufficiently large t values, when the system is in the diffusive regime.
Results for D are shown in the left panel of Figure 6 as a function of ρ along several isotherms. The diffusion coefficient shows a monotonous decrease of D on decreasing both T and/or ρ. From these data, we were able to trace isodiffusivity lines, selecting D values in a range from 10 −3 to 7.5 × 10 −5 in units of (2σ c ) 2 /t 0 . In the right panel of Figure 6 we report the resulting lines in the temperature vs. density plane. In the figure, the state points are characterized by different symbols to specify if they correspond to an equilibrium fluid phase or if they are characterized by a large value of S(0), indicating a phase-separated system.
It is worth to observe that many of the simulated state points are found to be in the part of the phase diagram where dynamic slowing down occurs, characterized by a diffusion constant less than 10 −4 . We also observe that the simulation data suggest a rather broad gas-liquid phase separation region, extending, e.g., from ρ * = 0.10 to ρ * = 0.40 at T * = 0.10. We note that MD simulations are consistent with the estimate of the critical point provided by GCMC simulations.

Conclusions
In the present contribution we considered overall neutral inverse patchy colloids (IPCs) with a charged equatorial belt and two oppositely charged polar patches defined by a relatively small opening angle. While the same system has been investigated under planar confinement in Ref. [23], we study here its bulk behavior in the fluid phase by means of numerical and theoretical approaches. On one side, we performed extensive Molecular Dynamics (MD) simulations over a wide range of temperatures and densities, including extremely low temperatures, where a gas-liquid phase separation takes place, as well as high densities, where a dynamical slowing down of the system occurs. On the other side, we applied an integral equation approach developed to describe heterogeneously charged systems with two or more patches [14]. We found that this theory is able to predict the static and the thermodynamic properties of the systems with high accuracy within the investigated temperature-density range; this holds even at state points characterized by the slowing down of the dynamics.
Additionally, we numerically characterized the dynamics of the system by calculating the velocity and orientation autocorrelation function as well as the diffusion constant. Finally, we traced a phase diagram where the gas-liquid phase separation region, the equilibrium fluid phase and the dynamically arrested state points are highlighted. We note that the observed bulk behavior, characterized by a wide region where the fluid phase is stable, is consistent with what observed in quasi two-dimensions. Indeed, under planar confinement, the system was shown to form an extended network of bonded particles with no well-defined spatial order. In contrast, IPCs with relatively bigger patches, also studied under confinement in Ref. [23], were shown to assemble into (spatially and orientationally) ordered planar arrays. Thus, we expect to observe a completely different bulk phase diagram for the latter system [25]. Figure 1. Schematic representation of the IPC model with two patches. The dark gray sphere is the hard sphere (HS) colloid with radius σ c . The yellow dots inside the particle at distance e from the particle center are the interaction sites, while the yellow spheres, partially located inside the HS core, represent the interaction zone of the patch with radius σ p ; only the cap of the yellow sphere that emerges from the HS core can interact with other particles, thus defining the patch angular semiamplitude γ. The external black circle of radius σ c + δ/2 represents the colloid interaction sphere. Geometric constraints impose that the interaction range is given by δ/2 = e + σ p − σ c .  Figure 3. Static structure factor S(k) of the system at different state points (as labeled). Solid green lines correspond to APY results, crosses refer to MD data. In the missing panels the system is phase separating.  . Potential energy per particle as a function of the density of the system along several isotherms (as labeled). Crosses are MD results, solid lines refer to APY data for the system with the HS core, and dashed lines correspond to APY data for the system with the soft sphere core.   Left: MD results for the diffusion coefficient along several isotherms (as labeled). Right: phase diagram of the system. The points highlighted in the temperature-density plane correspond to the numerically investigated state points: a red denotes a phase separating system, a blue corresponds to a homogeneous fluid, the black is the critical point. Along the gray curves the value of the diffusion coefficient D is constant and decreases from 10 −3 to 7.5 · 10 −5 in units of (2σ c ) 2 /t 0 .