**1. Introduction**

Competing attractive and repulsive interactions can be found in a wide variety of systems, ranging from block copolymers, proteins, or colloids, just to mention a few examples [1]. Even though the physical origin of the interactions are different in these systems, theory predicts that they all exhibit similar phase diagrams in which periodic microphases (cluster-crystal, hexagonal, bicontinuous gyroid and lamellar phases) are stable at low temperatures [2–6]. One might think that colloidal systems, in which the attractive and repulsive interactions can be tuned by modifying the colloidal solution, could be a good playground to experimentally study the formation of periodic microphases. Still, these periodic microphases have not ye<sup>t</sup> been experimentally observed in colloidal systems [7], which has been attributed to the particle size polidispersity [8], to the slow kinetics of the fluid [9–11], or to the inability of a simple effective potential to capture the behaviour of the colloidal solution [7]. Additionally, recent studies have suggested that apart from the strength of the interactions, the attractive and repulsive ranges play an important role and their variation can induce different phase behaviours [1,12]. In this regard, the ranges of interaction can be easily tuned in block copolymers by varying the length of the chains composed of one or another monomer [13]. Block copolymers self-assembly has many potential applications in nanotechnology and industry, such as separation and ion conduction in batteries, templating for nanomaterial synthesis and sensing [14,15].

**Citation:** Serna, H.; Gó ´zd ´z, W.T.; Noya, E.G. Structural and Dynamical Behaviour of Colloids with Competing Interactions Confined in Slit Pores. *Int. J. Mol. Sci.* **2021**, *22*, 11050. https://doi.org/10.3390/ ijms222011050

Academic Editor: Małgorzata Borówko

Received: 16 September 2021 Accepted: 10 October 2021 Published: 13 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland; hserna@ichf.edu.pl (H.S.); wtg@ichf.edu.pl (W.T.G.)

Obtaining the proper ranges in experimental colloidal systems has been challenging, but some recent approaches in which the colloidal particles are functionalised with hydrophobic molecules have shown promising results [16]. In a previous work [17], we showed how the Lennard-Jones plus Yukawa (LJY) potential, with the proper ranges and strengths of the interactions, can form ordered microphases in bulk. Here we want to stress the importance of simulations in predicting new physical phenomena. In particular, simulation is useful to guide the design of experiments that can lead finally to new discoveries. Regarding the applications, confined colloidal particles in channels of different geometries have been used to build wave-guide devices useful in sensing [18].

There are several ways in which the ordering of the periodic microphases can be induced, for example, by applying shear [19] or by confining the fluid in pores with the appropriate geometry [20]. In this work, we will explore this second route. It is known that confinement of simple and complex fluids can change the phase behaviour by shifting coexistence lines to lower or to higher temperatures than in bulk, depending on the shape and size of the pores and on the nature of the interactions of the fluid with the pore walls [21]. It can also induce significant changes on the dynamic behaviour, in some cases finding a nonmonotonous variation of the diffusion coefficient with the pore size [22–24]. In the particular case of systems with competing interactions, previous theory and simulation studies showed that confinement can promote or inhibit the formation of periodic microphases depending on whether the pore size is commensurate or not with the periodicity of the bulk microphase. For example, we showed in a previous work that confinement in channels with triangular and hexagonal cross-sections favour the formation of the hexagonal phase, as well as by introducing wedges in pores with cylindrical crosssections (which otherwise promote the formation of helical structures) [25,26]. We also found that new phases that are not stable in bulk can be stabilised when confined by the appropriate pore geometry. In this way, cluster-crystals with different symmetries were obtained by confinement in bicontinuous porous materials [27]. Surprisingly, the study of confinement in simple geometries, such as a slit pore, has not been sufficiently explored. Indeed we are only aware of a few studies in which fluids with competing interactions confined between parallel plates were studied, but those were restricted to two- and one-dimensional cases [20,28,29]. The literature on the dynamic behaviour of fluids with competing interactions in confinement is also scarce [30].

In this work we undertake a simulation study to investigate the effects of confinement on the structural and dynamic behaviour of fluids with competing interactions in narrow slit pores as a function of the pore width. The study is performed under conditions at which the cluster-crystal, the hexagonal and the lamellar phases are stable in bulk.

Our goal is to determine whether the formation of periodic microphases can be thermodynamically and/or kinetically favoured by confinement.

#### **2. The Model and the Simulation Method**

The colloidal particles interact with each other via an effective short-range attraction long-range repulsion (SALR) model potential resulting from the addition of a Lennard-Jones potential plus a Yukawa repulsive term:

$$u\_{SALR}(r\_{ij}) = 4\epsilon \left[ \left(\frac{\sigma}{r\_{ij}}\right)^{2a} - \left(\frac{\sigma}{r\_{ij}}\right)^a \right] + \frac{A}{(r\_{ij}/\zeta)} \exp\left(-r\_{ij}/\zeta\right) \tag{1}$$

The parameters of the model were assigned the same values as in our previous work in which the bulk phase diagram was investigated [17]. In particular, we chose = 1.6, *σ* = 1.0, *α* = 6, *A* = 0.65, and *ζ* = 2.0. For computational efficiency, the potential is truncated and shifted at *rc* = 4.0 *σ*. In what follows, all the magnitudes are reduced taking *σ* and as units of length and energy, respectively.

The confinement is implemented along the *z* direction by placing two parallel walls at *zw* = ± *W*/2, so that the separation between them is *W*. Periodic boundary conditions are imposed along the *x* and *y* directions. The walls are structureless and repulsive. Particles interact with the walls via a Lennard-Jones model truncated and shifted at the energy minimum:

$$\mathcal{V}\_{z\_{\rm av}}(z\_{\rm inv}) = \begin{cases} 4\epsilon\_{\rm uv} \left[ \left( \frac{\sigma\_{\rm uv}}{z\_{\rm inv}} \right)^{12} - \left( \frac{\sigma\_{\rm uv}}{z\_{\rm inv}} \right)^{6} \right] + \epsilon\_{\rm uv} & : z\_{\rm inv} < 2^{1/6}\sigma\_{\rm uv} \\ 0 & : z\_{\rm inv} \ge 2^{1/6}\sigma\_{\rm uv} \end{cases} \tag{2}$$

where *w* = 1.0, *σw* = 1.0, and *ziw* is the distance from particle *i* to the pore wall *zw*. The Lennard-Jones plus Yukawa interaction potential used to model the interactions between the particles and the truncated Lennard-Jones potential that accounts for the interactions between the particles and the walls are plotted in Figure 1.

**Figure 1.** The Lennard-Jones plus Yukawa potential used to model the interactions between colloidal particles and the Lennard-Jones potential truncated and shifted at the energy minimum used to model the interactions between the slit walls and the particles

Thus, the total energy of the system is given by:

$$\mathcal{U}\_{\text{tot}} = \sum\_{i=1}^{N-1} \sum\_{j>i}^{N} \mu\_{SALR}(r\_{ij}) + \sum\_{i=1}^{N} (\mathcal{V}\_{W/2}(z\_{i:w}) + \mathcal{V}\_{-W/2}(z\_{i:w})).\tag{3}$$

The slit width is given by the centre–centre separation between the confining walls, *W*. Taking into account that the energy of a particle becomes very repulsive for distances to the wall shorter than *σw*, the available width for the particle volume is actually 2*σw* smaller than the pore width. The edges of the simulation box were set to *Lx* = *Ly* = 40*σ*, and *Lz* = *W*, with *W*∗ = *W*/*σ* = 5.0, 7.0, 9.0 and 11.0. Thus, the number density is calculated as *ρ*∗ = (*<sup>N</sup>σ*<sup>3</sup>)/(*LxLyW*). Given the large dimensions of the simulation box along the *x* and *y* directions, we expect that finite-size effects will be small.

The phase behaviour of the confined SALR fluid was explored by performing a series of Monte Carlo (MC) simulations in the grand canonical ensemble at *T*∗ = *kBT*/ = 0.30 for several values of the chemical potential within −1.2 ≤ *μ*∗ = *μ*/ ≤ 0.5. For each wall separation, we chose three states at which the cluster-crystal, the hexagonal and the lamellar phases exhibit the most ordered structure (as compared to those obtained at other chemical potentials). The numbers of particles confined in the slit pore at each considered state are given in Table 1.


**Table 1.** Average number of particles confined in the slit pores at which the fluid organises into ordered structures at *T*∗ = 0.3 at densities at which the bulk fluid assembles into a cluster-crystal, a cylindrical and a lamellar phase. Note that the chemical potential for the more dense lamellar phase might not be reliable due to the low acceptance probability of the insertion/deletion MC moves. In any case, we only used these simulations to generate the initial configurations for the NVT MD runs.

∗denotes reduced variables.

Starting from these configurations, the confined fluid was then heated and cooled using Molecular Dynamics (MD) simulations in the canonical ensemble (*NVT*). The MD simulations were performed with the LAMMPS simulation package [31], in which the truncated and shifted SALR model described above was implemented in an external subroutine coded by us. The time step was set to *dt* = 0.005√*m<sup>σ</sup>*2/. Temperature was controlled with the Nose-Hoover thermostat with a relaxation time of 100*dt*. Simulations were evolved for 10<sup>6</sup> MD steps for equilibration, followed by another 10<sup>6</sup> MD steps for taking averages.

The structure of the fluid was identified mainly by visual inspection of local density plots. These plots were built by dividing the simulation box in small cubic cells of approximate edge length *σ*, measuring the particle density in each of these cells and averaging over 10,000 independent configurations, so that we can evaluate the local density function:

$$
\rho\_{xyz}(x,y,z) = \frac{\langle \mathcal{N}(x,y,z) \rangle}{\Delta V},\tag{4}
$$

where *<sup>N</sup>*(*<sup>x</sup>*, *y*, *z*) is the average number of particles in a cubic cell of edge *σ* and centred at the point (*<sup>x</sup>*, *y*, *<sup>z</sup>*), and Δ*V* is the volume of each small cubic cell, in our case Δ*V* = *σ*3. Isosurfaces of these density maps were visualised using OpenDX software. Clusters were also identified by performing a cluster size analysis [32], adopting the convention that two particles are nearest neighbours if the distance between them is lower than *rcut* = 1.6*σ* for the cluster-crystal and hexagonal phases and lower than *rcut* = 1.4*σ* for the lamellar phase, i.e., roughly the distance to the first minimum in the radial distribution function of each periodic microphase [17]. This information was used to calculate the cluster size distribution (CSD).

The spatial distribution of the particles along the direction perpendicular to the pore walls was investigated by measuring the density profiles, calculated by dividing the pore volume in small slabs of width Δ*z* = 0.1*σ* and averaging the number density in each of these slabs:

$$\rho\_{\overline{z}}(z) = \frac{\langle \mathcal{N}(z + \Lambda z) \rangle}{L\_x L\_y \Delta z}. \tag{5}$$

Here *<sup>N</sup>*(*z* + <sup>Δ</sup>*z*) is the ensemble average of the number of particles in the slab between *z* − Δ*z*/2 and *z* + Δ*z*/2, and *Lx* and *Ly* are the two periodic edges of the simulation box.

Following our preliminary study of the bulk system [17], we also investigated the internal ordering of the clusters at the particle scale as a function of temperature. For that purpose, for each periodic microphase, we chose an order parameter that is able to discriminate particles within local ordered environments from those within local disordered environments. As in the bulk system, the spherical and cylindrical clusters that form the cluster-crystal and hexagonal phases at low temperatures have local icosahedral

symmetries. In this case, a common neighbour analysis (CNA) [33] allows us to distinguish particles with local icosahedral environments from liquid environments. The CNA analysis was made with the OVITO visualisation tool [34], using a fixed cutoff radius of 1.6 *σ*, that corresponds to the distance to the first minimum in the pair distribution function in the bulk cluster-crystal and hexagonal phases [17]. On the contrary, in the lamellar phase, particles are arranged in stacks of hexagonally-packed layers. Therefore, it is more convenient to use the Lechner and Dellago order parameter, that is able to effectively distinguish particles in ordered local environments (i.e., solid-like, including particles in the frozen lamellae [17]) from those in disordered local environments. For the lamellar phase, first neighbours were defined using a slightly shorter cutoff distance than for the cluster-crystal and hexagonal phases of 1.4 *σ*, corresponding to the first minimum in the pair distribution function of the bulk lamellar phase [17].

Finally, we also measured the mean squared displacement (MSD) that provides information on the single particle dynamics:

$$
\left< \Delta r(t)^2 \right> = \left< \frac{1}{N} \sum\_{i=1}^{N} (r\_i(t) - r\_i(0))^2 \right>,\tag{6}
$$

where *<sup>r</sup>i*(*t*) and *<sup>r</sup>i*(0) are the positions of particle *i* at times *t* and zero, respectively. The diffusion coefficient, *D*, is estimated from Einstein's relation:

$$D = \frac{1}{2d} \lim\_{t \to \infty} \frac{\partial \langle \Delta r(t)^2 \rangle}{\partial t} \tag{7}$$

where *d* is the dimensionality of the system. For the confined systems, we calculated the diffusion coefficient in the direction parallel to the walls ( *<sup>D</sup>*), because the particle displacement in the perpendicular direction is limited by the narrow width of the pores. In this case, the MSD is calculated using only the *x* and *y* coordinates, and *d* is set to 2. For the bulk system, movements in the three directions of space are considered and *d* = 3. To calculate the diffusion coefficient, we divide the MSD data into 10 independent blocks. Then, we fit the MSD to a straight line and calculate the diffusion coefficient in each block following Equation (7). We discard the first steps in which the systems usually exhibit ballistic behaviour. We only calculate the diffusion coefficient when the MSD scales linearly with *t*, i.e., in the diffusive regime. The diffusion coefficients are averaged over the independent blocks and the errors are estimated as the standard deviation of the sample of blocks.
