Effects of confinement on pattern formation in two dimensional systems with competing interactions

Abstract

Template-assisted pattern formation in monolayers of particles with competing short-range attraction and long-range repulsion interactions (SALR) is studied by Monte Carlo simulations in a simple generic model [N. G. Almarza et al., J. Chem. Phys., 2014, 140, 164708]. We focus on densities corresponding to formation of parallel stripes of particles and on monolayers laterally confined between straight parallel walls. We analyze both the morphology of the developed structures and the thermodynamic functions for broad ranges of temperature T and the separation L₂ between the walls. At low temperature stripes parallel to the boundaries appear, with some corrugation when the distance between the walls does not match the bulk periodicity of the striped structure. The stripes integrity, however, is rarely broken for any L₂. This structural order is lost at T = T_K(L₂) depending on L₂ according to a Kelvin-like equation. Above the Kelvin temperature T_K(L₂) many topological defects such as breaking or branching of the stripes appear, but a certain anisotropy in the orientation of the stripes persists. Finally, at high temperature and away from the walls, the system behaves as an isotropic fluid of elongated clusters of various lengths and with various numbers of branches. For L₂ optimal for the stripe pattern the heat capacity as a function of temperature takes the maximum at T = T_K(L₂).

---

1 Introduction

Competing interactions, in particular short-range attraction and long-range repulsion (SALR), are present in many biological and soft matter systems^1–4 as well as in magnetic films with competing ferromagnetic and dipolar forces.^5,6 The repulsion is often of electrostatic origin, and the attraction between the particles comes from the van der Waals forces, or it is induced by the solvent. For example, the solvophobic attraction is present between globular proteins in water,^2,7 or the depletion attraction is present between colloid particles when the solvent contains nonadsorbing polymer with small radius of gyration.^2,4

Despite the different origins of the SALR potential in the above examples, the structural and thermodynamic properties of the systems with competing interactions are very similar.^5,8–18 Moreover, a self-assembly into patterns that are very similar to those observed in amphiphilic systems when the range of attraction is sufficiently large and the repulsion is sufficiently strong takes place.^8–10 At sufficiently low temperature the following sequence of ordered phases was predicted for increasing volume fraction of particles by theory and simulation:^10,18–23 spherical clusters, cylindrical clusters, gyroid network of particles, layers of particles, gyroid network of voids, cylindrical voids and finally spherical voids. Heating of the system first destroys the gyroid phase. Further increase of temperature leads to disordered distribution of the clusters or layers. At low volume fractions in the disordered cluster phase spherical clusters are present. Upon increase of the volume fraction the clusters become elongated and finally percolate.^4,20 Much higher temperature is necessary for disassembly of the clusters. Structure factor, equation of state, specific heat etc. in the disordered cluster phase are significantly different than in the homogeneous fluid.^24–27

In this work we focus on two dimensional (2D) patterns. Such patterns are formed in thin magnetic films,^5,6 in thin layers of block-copolymers,^28–30 by particles adsorbed at solid surfaces or at liquid interfaces,^1,31,32 on elastic membranes or embedded in lipid bilayers.^3,33 For increasing chemical potential of the particles, hexagonally ordered clusters, next stripes and finally hexagonally ordered voids occur for various versions of the SALR potential.^11–17 The 2D self-assembly in continuous SALR models was studied by the density functional theory¹⁵ and by Monte Carlo (MC)^12–14,34 methods. Recently lattice SALR models (LSALR) were introduced and studied in the mean-field approximation (MF) and by MC simulations.^11,16,17,25 Here we investigate the effects of confinement between two parallel boundaries on the degree of order in monolayers of particles interacting through the SALR potential, using the lattice model.^11,16

The effects of confinement on surfactant solutions were investigated experimentally and theoretically already long ago.^35,36 Recently the attention focuses mainly on confined block-copolymers.^37–40 The effects of confinement depend strongly on the surfactant concentration in the amphiphilic systems, or on the volume fraction of particles in the SALR systems. The confinement effects for different densities in 2D continuous SALR models were studied in ref. 15 and 34. In this work we limit ourselves to the lamellar or stripe phase. In the surfactant systems the parallel orientation of the lamellas was stable for any wall-separation for large-period phases.^35,36 For phases with intermediate periods confined between weakly selective boundaries, parallel or perpendicular orientation of lamellas was found when the wall separation and the period of the bulk structure were or were not commensurate, respectively. In the case of thin lamellas and selective boundaries a zig-zag structure of layers was obtained in the case of the incommensurability.³⁶

Slit-like confinement in the case of the 2D SALR system has been studied for rather thick stripes.^14,15 Parallel orientation of layers was observed when the separation between the boundaries was commensurate with the bulk periodic structure. For incommensurate slit-widths, layers perpendicular to the boundaries for the whole slit or its central part were formed.^14,15 This similarity between confined amphiphilic and SALR systems could be expected based on the similarity between properties of these systems in the bulk. By analogy, because in the amphiphilic systems thick and thin lamellas respond to confinement in a significantly different way, we may expect that the effects of confinement in the SALR system self-assembling into thick and thin stripes are different. However, in the case of thin stripes the effect of confinement has not been investigated yet.

The thickness of the stripes depends on the range of the attractive and the repulsive part of the potential. To study the effects of confinement in the case of thin stripes, we shall consider the 2D LSALR model introduced and studied in the bulk in ref. 11 and 16, since in this model thin stripes are formed. In ref. 11 a rather detailed phase diagram was obtained by MC simulations. It was found that at low temperature, T, the lamellar (L) phase with both orientational and translational order of stripes is stable for the reduced density ρ ∼ 1/2. This phase melts in two steps when heated. First the translational order is lost. The phase possessing only the orientational order was named “molten lamella” (ML). In the ML the translational order is lost due to topological defects, such as fracture, branching or junction of the stripes.^29,30 The system is not isotropic, however. On the triangular lattice studied in ref. 11 the total length of stripe segments along the main lattice directions is not the same, i.e. there exists preferred orientation of the stripes. The unit vectors normal to the stripe segments are analogous to a discretized version of the director field in liquid crystals, therefore the ML phase is analogous to the nematic phase. Further heating leads to the transition from the ML to the isotropic fluid (F), when the separation between the topological defects becomes comparable with the period of the local lamellar order. In the F phase the total length of the stripe segments parallel to each principal lattice direction is the same. The transition between the ML and F phases was found to be continuous for a density interval centered at the reduced density ρ = 1/2 that is optimal for the lamellar structure.¹¹ The line of continuous transitions terminates at two tricritical points, one with the density ρ_tcp ≈ 0.4, and the other one with ρ_tcp ≈ 0.6, where the transition becomes first-order.¹¹ The ML phase is analogous to the nematic phase of stripes found in magnetic films.⁵ Recent studies of the 2D stripe-forming systems⁴¹ indicate that the transition between the anisotropic and isotropic phases can be of the second order for Coulomb-like repulsion, or of the Kosterlitz–Thouless type for dipolar interactions. Thus, the order of the transition depends on the range of the repulsion, and in general it is a very subtle issue. In this work we shall limit ourselves to ρ ∼ 1/2 in a system confined between two parallel boundaries, where no true phase transitions occur. The confined system is in contact with the bulk reservoir at fixed T and chemical potential.

Self-assembly into parallel stripes of particles might be a first step in the formation of ordered patterns of practical importance. From the point of view of some applications, it would be advantageous to have thin parallel self-assembled stripes without topological defects. In some other cases interconnected stripes forming porous structure would be preferable. The number of topological defects in the ML phase is significantly larger than in the L phase. Motivated by the phenomenon of capillary condensation in simple fluids, capillary lamellarization in surfactant solutions and capillary smectization in liquid crystals,^42–44 we shall study if the ML structure can be transformed into the L structure by decreasing the separation between the boundaries at fixed temperature and chemical potential. We shall also investigate how the structural transformations are reflected in properties of thermodynamic functions such as the heat capacity and the solvation pressure.

In the next section the model and the simulation method are described, and the thermodynamic and structural functions are introduced. In Section 3 we present our results. The ground state (T = 0) is discussed in Section 3.1. Results for the heat capacity and the solvation pressure are described in Section 3.2. In Section 3.3 we discuss the structure. Typical configurations, and results for the density profiles across the slit and for the correlation function along the slit are presented for several values of T and for wall–wall separations L₂ commensurate and incommensurate with the period of the bulk structure. Finally, the orientational order parameter as a function of T is obtained for a few values of L₂. Section 4 contains our conclusions.

2 Model and simulation procedures

We consider a surface in equilibrium with a bulk reservoir and assume that the particles can occupy the sites of a triangular lattice with a lattice constant, σ, comparable with the diameter of the adsorbed particles. The lattice sites are denoted by x = x₁e₁ + x₂e₂ ≡ (x₁,x₂), where e₁, e₂ and e₃ = e₂ − e₁ are the unit lattice vectors on the triangular lattice, i.e. |e₁| = |e₂| = |e₁ − e₂| = 1 (in σ-units), and x_i are integer. We consider the triangular lattice, because it allows for close packing of the particles.

We introduce the slit-like confinement by assuming that the particles are located between two lines parallel to e₁, and separated by L₂ layers of sites along the direction e₂. Thus, 1 ≤ x₂ ≤ L₂. In the direction e₁ we assume periodic boundary conditions (PBC). In simulations we assume 1 ≤ x₁ ≤ L₁, where L₁ ≫ L₂, and x₁ = L₁ + 1 is identified with x₁ = 1. The system under consideration is shown in Fig. 1.

Fig. 1 Sketch of the two-dimensional triangular lattice with the slit boundary conditions. L₁ is the length of the system in direction e₁ (in this direction periodic boundary conditions are assumed), whereas L₂ is the number of layers between the confining walls. The walls interact only with the particles which occupy the first neighboring row of the wall, and the dimensionless interaction energy is denoted by h* = h/J₁, where J₁ is the nearest-neighbor attraction and h is the interaction energy between the particle in the boundary layer and the confining wall.

The interaction of the particles with the binding sites on a solid substrate, or with the lipids in the membrane plays analogous role as the chemical potential of the particles, and we introduce μ that is the sum of the binding energy and the chemical potential. We assume that the nearest-neighbors attract each other (SA), the interaction changes sign for the next-nearest neighbors, becomes repulsive for the third neighbors (LR), and vanishes for larger separations. The nearest-neighbor attraction is the standard assumption in the lattice-gas models. In the case of charged particles in electrolytes the assumed range of repulsion should be of order of the Debye screening length, 2.5σ ∼ λ_D. Since in various solvents with weak ionic strength λ_D ∼ 1–100 nm, the model is suitable for charged molecules, nanoparticles or globular proteins. The interaction between the occupied sites x and x + Δx of the form described above is given by

where −J₁ and J₂ represent the attraction well and the repulsion barrier respectively.

The thermodynamic Hamiltonian has the form

where denotes the summation over all lattice sites, the microscopic density at the site x is [small rho, Greek, circumflex](x) = 1(0) when the site x is (is not) occupied, and h denotes the interaction energy between the particle in the boundary layer and the confining wall.

The microscopic state of the whole system is specified by indicating which sites are occupied, and is denoted by {[small rho, Greek, circumflex](x)}. The probability of a particular microscopic state has the form

p[{[small rho, Greek, circumflex](x)}] = Ξ⁻¹ exp(−βH[{[small rho, Greek, circumflex](x)}]), (3)

where β = 1/(k_BT) and k_B is the Boltzmann constant. The grand potential is expressed in terms of the grand statistical sumin the standard way

Ω = −k_BT ln Ξ. (5)

We choose J₁ as the energy unit, and use the notation X* = X/J₁ for any quantity X with the dimension of energy, except from the repulsion to attraction ratio, denoted by J = J₂/J₁. Reduced temperature is defined as: T* = k_BT/J₁. Along the paper we have used the same simulation techniques as used for this model in the bulk,¹¹ that is: Metropolis Monte Carlo, thermodynamic integration and parallel tempering technique. More details about the methodology are given in the Appendix.

2.1 The ground state

In the previous work¹⁶ we determined the ground state (GS) of the model and the MF phase diagram. Later we computed¹¹ the phase diagram by means of MC simulations for strong repulsion, J = 3. In this work we limit ourselves to J = 3, and to the chemical potential μ* = 6, for which the L phase is stable.^11,16 In the most of the cases we shall assume h* = −1, which is a moderate wall–fluid interaction.

The grand potential for T = 0 reduces to the minimum of the grand canonical Hamiltonian per site H̄₀ ≡ H[{[small rho, Greek, circumflex](x)}]/(L₁L₂). The ground states for μ* = 6 and different wall–wall separations L₂ and values of h* were computed as follows: for fixed values of μ*, L₂, and h*, we consider different values of the lateral length L₁, and determine by means of parallel tempering MC simulations their corresponding minima, H̄₀(L₁). The GS configuration and the GS energy, at the given conditions (μ*,L₂,h*), are taken as those corresponding to the value L₁ which minimizes H̄₀.

2.2 Thermodynamic properties

The heat capacity. We start the investigation of the thermodynamic properties of the model at T* > 0 by computing the heat capacity defined as:where the angular brackets represent the average of the corresponding quantity for given conditions in the grand canonical ensemble. H is the grand canonical Hamiltonian given by (2).

The solvation pressure. The impact of the system geometry and the particle–wall interactions on thermodynamic properties is described by the excess grand potential

ω_ex ≡ (Ω − Ω_bulk)/L₁ = γ₁ + γ₂ + Ψ(L₂), (7)

where Ω_bulk = −k_BT ln Ξ_bulk is the grand potential of the bulk system with the same volume, i.e. with PBC in both directions, γ_i is the wall–fluid surface tensions at the i-th wall, and Ψ(L₂) is the interaction between the confining walls per unit length induced by the confined fluid (solvation pressure). In order to compute Ψ(L₂) we have applied the MC thermodynamic integration procedure for fixed chemical potential, which was previously used for the bulk system.¹¹ The grand potential per lattice site was computed using the following relationwhere H is the grand canonical Hamiltonian, in the case of the slit given by (2). In the case of the bulk, the interactions with the walls are neglected in (2), and PBC are assumed in the direction e₂. The term −ln 2 results from the assumption that at infinite temperature our model acts like an ideal lattice gas. To obtain Ψ(L₂) out of ω_ex = (Ω − Ω_bulk)/L₁, we note that these two quantities at fixed temperature differ only by a constant (see (7)). Since , we have estimated the constant by the value of ω_ex for which the curve reaches a plateau.

2.3 Structural description

The density profile. In order to quantify the description of the structure across the slit, we introduced the density profile averaged over the longitudinal direction e₁,

The longitudinal correlation function. The structure in the longitudinal direction, e₁, is studied by defining a correlation function along each line parallel to the boundary, and average it over all the lines located at 1 ≤ x₂ ≤ L₂. This global correlation function along the confining walls allows us to compute the average period and decay length in direction e₁. We define the global density correlation function along the direction e₁ aswhere [small rho, Greek, vector](x⁰₁) ≡ ([small rho, Greek, vector](x⁰₁,1),…,[small rho, Greek, circumflex](x⁰₁,L₂)), is a vector containing the densities of the sites belonging to a compact line of sites in the direction e₂ (see Fig. 2), · is the standard scalar product, and 〈⋯〉 indicates an average over different system configurations. From g₀(x₁) we can define a normalized correlation function:

g(x₁) ≡ g₀(x₁)/g₀(0). (11)

Fig. 2 Sketch of the definition of the density vectors, [small rho, Greek, vector], used in the definition of the global correlation function g₀(x₁). Here the vectors are oriented in the e₂ direction. The red rings connect the sites which belong to one density vector. In the scheme two density vectors are depicted: [small rho, Greek, vector](x⁰₁) = (1,1,0,0,0,0,1,1) and [small rho, Greek, vector](x¹₁) = (1,1,0,0,1,1,0,0), with x₁ = x¹₁ − x⁰₁.

The orientational order parameter. We have made use of an orientational order parameter O_L similar to that of the bulk system.¹¹ In short, we count the total lengths (S₁, S₂, S₃) of the segments of the stripes oriented in each of the three main directions of the lattice, and the order parameter is defined as:In the isotropic phase we expect S₁ ≃ S₂ ≃ S₃, and therefore O_L should vanish in the thermodynamic limit, whereas in the bulk L phase there are lamellar segments only in two directions,¹¹ and therefore O_L = 1. In order to focus on the behavior of the region not directly influenced by the walls, only sites with coordinates 2 < x₂ < L₂ − 1 were considered in the calculation.

3 Results

3.1 The ground state

Long-range order of the system for given μ and L₂ is found for T = 0. Before studying finite temperatures with thermally induced defects, we present how the bulk GS is modified in the slit confinement. Let us first briefly remind the GS in the bulk. For strong repulsion the energy in the bulk assumes a minimum when neither intra-cluster nor inter-cluster repulsion is present, and as many nearest-neighbors as possible are occupied. The energy competes with −μN, where N is the number of occupied sites. In this model the Hamiltonian takes the same value for straight stripes, and for zig-zag stripes that in one of the main lattice directions have the thickness 2.^11,16 The segments of the zig-zag stripes can be parallel to two of the three main lattice directions. Thus, in the bulk the GS is strongly degenerated in the stability region of the L phase.

If the confining surfaces strongly attract particles, then energetically favorable configurations have double layers (stripes) of particles adsorbed onto the walls. For fixed J the configuration of the particles between the adsorbed stripes will depend upon the distance between the confining walls, L₂, and on the wall–fluid interaction h*.

If L₂ = 4n + 2, where n is a positive integer, then the distance is commensurate with the period of the bulk structure and the straight lamella structure has the lowest energy. Note that in the case of the straight lamella structure no degeneracy of the GS is present. For L₂ incommensurate with the period of the bulk structure geometrical defects appear in the stripes adsorbed at the walls (see Fig. 3 for h* = −1). The GS depends on both L₂ and h*, and there are many GS configurations in the (L₂,h*) plane. To fix attention we assume h* = −1, except from a few cases illustrating the effect of h*. We call the defects resulting from the incommensurability between the lamellar period and L₂ at T = 0 “geometrical”, because the continuity of the stripe is not broken. The stripe next to the wall becomes locally thicker or thinner, or undulates to bypass a void at the wall. The first case can be seen for example in the first column in Fig. 3. The second and the third cases are shown in the two top panels in the third and the central column in Fig. 3, respectively. In order to avoid excess repulsion, the nearest stripe (second from the wall) has to bypass the defect. The turns of the nearest stripe entail turns of the other stripes inside the slit in such a way that it is energetically favorable to build a defect at the opposite wall. Such construction ensures that the distance between every two defects is the same, hence three consecutive defects in the e₁ direction form an equilateral triangle (Fig. 4).

Fig. 3 The ground state (T = 0) configurations for L₂ = 4n − 1 (left column) L₂ = 4n (center column) and L₂ = 4n + 1 (right column), where n = 2, 3, 4, 5. The confinement attracts the particles with h* = − 1. L₁ is the period of the structure in direction parallel to the confining walls. For the sake of figure clarity the walls are not shown.

Fig. 4 The defects are in the vertices of the equilateral triangle.

We found that for slit widths L₂ > 13 the distance L₂ determines the shape and the position of the defects, hence having L₂ one can also tell the size of the unit cell in the e₁ direction, L₁, which is given by the following formula

We cannot tell, however, if for large wall–wall separations (e.g. L₂ > 30) L₁ is still determined uniquely by L₂.

In contrast to the bulk, the segments of the lamella can be parallel not only to the two main lattice directions, but to all the three directions. Moreover, the presence of the confinement removes the degeneracy and the residual entropy of the GS.

In order to illustrate the effect of the wall–fluid interaction, we present in Fig. 5 the GS for L₂ = 16 and the neutral, h* = 0, and strongly repulsive, h* = 3, boundaries. Note that the period in direction e₁ is L₁ = 12 for h* = −1 (see Fig. 3), and L₁ = 16 for h* = 0 (see Fig. 5), but the general structure of the system is quite similar in the case of the attractive and neutral walls. In the case of strongly repulsive walls the particles are expelled from the bottom and top rows of the slit, and in the GS configuration four parallel straight stripes that do not interact with the walls occur (right panel of Fig. 5). Hence having the wall separation L₂ = 4n one can remove the defects by changing the particle–wall interaction from attractive to strongly repulsive.

Fig. 5 Ground state configuration for L₂ = 16 for different interactions with the walls. Left panel: h* = 0, right panel: h* = 3.

We conclude that in the case of thin stripes the incommensurability between the slit width and the lamellar period leads to zig-zag or undulating stripes, and to a periodic order in the longitudinal direction. Very similar behaviour was found previously for thin lamellas in water–surfactant mixtures.³⁶ The detailed shape of the undulating stripes and the geometrical defects at the near-surface layers both depend on the wall–wall separation, on the wall–fluid interactions and on the interparticle potential. For strong interparticle repulsion the continuity of the stripes is not broken at T = 0, i.e. the topological properties of the stripes do not change, when the wall–wall separation increases a little. With further increase of L₂ a new stripe is introduced into the system.

3.2 Thermodynamics

The heat capacity. As in the GS, we take h* = −1 and μ* = 6, but we limit ourselves to the pore widths L₂ = 6, 10, 14, 18,…,38 (i.e. those sizes whose GS present straight stripes). In Fig. 6 we show c_μ as a function of the temperature. For the narrowest slits only one maximum of c_μ(T) appears, whereas for L₂ ≥ 30 there are two maxima. One of them is located at T* ≃ 0.7. The other one appears at lower temperature that decreases significantly with increasing L₂.

Fig. 6 Reduced heat capacity, c_μ*(T) ≡ c_μ(T)/k_B for μ* = 6 and h* = −1, and for different pore widths. c_μ(T) exhibits two maxima for L₂ > 30, and an inflection point for L₂ = 30.

We do not expect true thermodynamic phase transitions in slit systems at T* > 0, given the fact that the system is virtually macroscopic in only one direction of the space. However, in the 1D LSALR model pseudo-phase transitions were observed,²⁵ and we can expect such pseudo-phase transitions, or crossovers between structures with different degree of order, in the 2D slit.

In the bulk, the continuous F–ML transition was found¹¹ for μ* = 6 at T* ≃ 0.75 that is close to the location of the high-T maximum of c_μ. On the other hand, MC simulations¹¹ for L₁ = L₂ = 120 show that in the bulk the ML–L transition occurs at T_bulk* ≃ 0.25 that is much lower than the temperature at which c_μ takes the low-T maximum for L₂ ≤ 38. Note, however that in simple fluids the phenomenon of capillary condensation leads to the shift of the temperature at which the condensation occurs in thin capillaries.⁴² The thinner the capillary, the higher the temperature of the capillary condensation. The phenomenon is described by the Kelvin equation T_K*(L₂) = T_bulk* + a/L₂, where T_bulk* is the transition temperature in the bulk, and a is associated with the difference in densities and entropy per particle in the coexisting phases.^42,45 Similar behavior, i.e. capillary lamellarization or capillary smectization was observed in surfactant solutions or in liquid crystals.^43,44

In order to verify if the low-T* maximum of the specific heat is associated with the capillary lamellarization, we plot in Fig. 7 the temperature T_max*(L₂) corresponding to the maximum of c_μ. In addition, we plot the best fit of the data to the equation

T_K*(L₂) = T_b* + a/L₂. (13)

We obtain a fair agreement between T_max*(L₂) and (13), with T_b* similar to T_bulk*. Thus, we conclude that the solid line shown in Fig. 7 represents the pseudo-transition between the structure with and without the periodic order (capillary melting of the ordered L structure upon heating). Since there are no thermodynamic phase transitions in our slit, we expect a crossover between the L and the ML or F structures that occurs for some temperature range around T_K*(L₂). The shaded region in Fig. 7 is a very rough estimation of the crossover between the L and ML or F structures based on the width of the c_μ peak. Recall that in the bulk there are very few topological defects in the L phase, and many topological defects in the ML and F phases. Based on the above thermodynamic considerations, we expect small and large number of the topological defects on the low-T and the high-T sides of the crossover region respectively.

Fig. 7 Structure diagram based on the features of the heat capacity and particle distribution considerations. Symbols represent temperature corresponding to: (i) the low-T global maximum, c^max_μ (bullets), (ii) the high-T maximum or the inflection point of the heat capacity (small symbols connected by the dashed line), (iii) the boundary between the low-T and high-T regions where the structure does and does not dependent on the commensurability between the slit size and the period of the stipe pattern (big triangles). The reduced chemical potential and wall–fluid interaction are μ* = 6 and h* = −1, and the pore width is L₂ = 4n + 2 with integer n (see Fig. 6). The solid line is the best fit to the Kelvin-like equation T_K*(L₂) (see (13)), where the fitting parameters are T_b* = 0.28 and a = 7.422. Note that T_b* is close to T_bulk* ≃ 0.25 associated with the transition between the L and ML phases in the unrestricted monolayer. The shaded region around T_K*(L₂) corresponds to the width of the peaks of c_μ(T*), estimated based on temperatures for which c_μ(T*) ≈ (c^max_μ + c^p_μ)/2, where c^p_μ is the high-temperature value of the heat capacity (the plateau in Fig. 6). L denotes the periodically ordered structure. ML and F are the anisotropic and isotropic structures without the periodic order respectively.

The Kelvin-like equation was fitted to the maxima of the heat capacity only for L₂ commensurate with the periodic structure, because in ref. 43 it was observed that the capillary lamellarization in surfactant mixtures is delayed if the slit width is incommensurate with the periodic structure. By analogy, deviations from (13) for incommensurate L₂ are expected in the SALR system. We verified that the geometrical defects present for L₂ incommensurate with the period of the striped structure influence the heat capacity too. For this reason we did not consider the values of L₂ for which the simple eqn (13) cannot be valid.

The solvation pressure. A typical shape of Ψ(L₂) is shown in Fig. 8 for T* = 0.5. The shape of Ψ(L₂) is very similar to that obtained for the 1D LSALR model⁴⁶ by exact calculations, and to the Ψ measured experimentally in surfactant mixtures.³⁵ It can be observed that the effective potential between the walls exhibits a significant decay length, with an oscillatory behavior superimposed on an attractive background. In addition, slit widths that fulfill L₂ = 4n + 2, with n = 1, 2,…,5 correspond to local minima of Ψ(L₂) for L₂ ≤ 22. We can anticipate that this effect is due to the good fitting of the straight lamellar stripe structures to those widths. Note, however, that for L₂ > 22 the minima are located at L₂ = 31 and L₂ = 35, which do not correspond to the wall separation optimal for the straight lamella structure. For T* > 0, however, the lamellar period may increase due to thermally induced defects. In fact in the 1D lattice model the period of the oscillatory decay of Ψ(L₂) obtained by exact calculations is a real number.⁴⁶

Fig. 8 The effective wall–wall interaction potential as a function of the wall–wall separation, Ψ(L₂), at temperature T* = 0.5 for L₁ = 120, μ* = 6 and attractive walls with h* = −1.

Since in a lattice model Ψ(L₂) is meaningful only for integer L₂, we did not try to fit Ψ(L₂) to a smooth curve. In the 1D case the exact Ψ(L₂) exhibits exponentially damped oscillations for large L₂, and similar asymptotic behavior, i.e. exponentially damped oscillations for large L₂, can be expected in this case. However, in 1D the asymptotic form of Ψ(L₂) deviates significantly from the exact solution for relatively large L₂, and by analogy we expect that for L₂ < 30 no simple analytical formula can describe Ψ(L₂) in Fig. 8.

At large separations L₂ the interaction between the confining walls follows from the mismatch between L₂ and the period of the periodic structure, therefore it should vanish in the absence of the periodic order, i.e. beyond the crossover between the L and ML or F structure. In particular, for T* = 0.5 from Fig. 8 we can estimate this range as ≃40, which is close to L₂ ≃ 38 at the boundary of the shaded region in Fig. 7 for this temperature. We conclude that the range of Ψ(L₂) increases significantly with decreasing temperature, and its very rough estimate is given by the high-T boundary of the shaded region in Fig. 7.

3.3 The structure

We are interested in the structure in the slit for temperatures below the first and above the second peak of the heat capacity, and between the two peaks. We shall verify if below the first and above the second peak the structure resembles the L and F phases, and if between the peaks the ML-like structure is present, as suggested by the thermodynamic considerations.

Representative configurations. It is instructive to compare the snapshots for different temperatures for the commensurate and incommensurate slit widths. In Fig. 9 representative configurations for slit widths L₂ = 30, 31, 32, 33, 34, with L₁ = 120, and at T* = 0.25 are shown. For L₂ = 30 and L₂ = 34 the system presents respectively eight and nine straight stripes, which are parallel to the walls. The periodicity of the system in direction perpendicular to the walls is four, and the top and bottom stripes are located close to the walls. For 30 < L₂ < 34 the structure of the system can be described as composed by eight stripes, but the stripes that are not close to the walls are no longer straight, showing some cooperative corrugation. By inspection of the configurations presented in Fig. 9 it seems that at low T* the system self-assembles in an almost periodic structure in the e₁ direction, similar to the GS. Note that the chosen temperature corresponds to very small values of the heat capacity (see Fig. 6). The effect of the temperature is shown for L₂ = 30 and L₂ = 32 in Fig. 10. The bottom panels correspond to the temperature T* = 0.403 that is lower than the low-T* maximum of the heat capacity, T_max*(30) ≈ 0.53. In the case of L₂ = 30 almost no geometrical defects are present, whereas for L₂ = 32 more geometrical defects than in the GS appear. They lead to small distortions of the ordered zig-zag structure. For temperatures T* < T_K*(L₂) (see (13)) the system seems to exhibit large correlation at relatively large distances in the e₁ direction. The topological defects such as the fracture or branching of the stripe remain almost absent. The shown configurations agree with the thermodynamic considerations suggesting stability of the L-type structure for T* < T_K*(L₂). The second and third panels from the bottom correspond to temperatures between the two maxima of c_μ. For T* = 0.556 the topological defects appear in both slits, but the GS-like undulations and pieces of straight stripes are still visible for L₂ = 32 and L₂ = 30 respectively. The anisotropic structure with topological defects characterizes the ML phase. For T* ≈ 0.625 that is close to the minimum of the heat capacity one cannot distinguish the structure for L₂ = 30 and L₂ = 32 by visual inspection. This indicates that the topological defects that appear independently of L₂ dominate. The two top panels correspond to temperatures T* = 0.806, and T* = 1.00 that are above the high-T* maximum of c_μ. Note that for T* = 1 the structure looks isotropic, whereas for T* = 0.806 some degree of anisotropy is visible. Note, however that the very broad peak of the heat capacity may indicate large crossover region between the anisotropic ML and the isotropic F structures.

Fig. 9 Representative configurations of the systems with attractive walls (h* = −1), at low temperature T* = 0.25, and different pore widths. From top to bottom: L₂ = 30, L₂ = 31, L₂ = 32, L₂ = 33, and L₂ = 34. Only the system sites are plotted.

Fig. 10 Representative configurations of the system with attractive walls (h* = −1) and L₁ = 120 for L₂ = 30 (left column) and L₂ = 32 (right column) and several temperatures. From top to bottom T* = 1.00, T* = 0.806, T* = 0.625, T* = 0.556, and T* = 0.403. The bottom panels correspond to temperature below the global maximum of the heat capacity, i.e. below T* = T_K*(30) (see (13)). T* = 0.556 is slightly above T_K*(30). Temperature at the central panels, T* = 0.625, is close to the minimum of the heat capacity. The two top panels correspond to temperature above the high-T maximum or the inflection point of c_μ*, i.e. above the dashed line in Fig. 7.

The above inspection of the representative configurations shows consistency of the thermodynamic and structural properties. Before looking at the structural functions, ρ(x₂) and g(x₁), let us consider a single line parallel to the boundary. To fix attention let us choose the 9th line from the bottom. At low T* it is completely filled with particles for L₂ = 30, and there are alternating occupied and empty segments for L₂ = 32, because of the stripe corrugation. The density averaged over the longitudinal direction for x₂ = 9 is smaller for L₂ = 32 than for L₂ = 30. Further away from the boundary this difference in the average density increases. Moreover, for L₂ = 32 the alternating occupied and empty segments along the direction e₁ lead to an oscillating correlation function. The period of the oscillations at very low T* is determined by L₂, as in the GS. For T* ≃ T_max*(30) ≈ 0.53 the topological defects lead to additional alternating occupied and empty segments in each line parallel to the boundaries for L₂ = 32, and to appearance of such segments for L₂ = 30. The average length of the occupied and empty segments decreases with increasing temperature, because the number of the topological defects increases. This leads to smaller density averaged over the longitudinal direction in lines that were occupied at low T*, and to a smaller period of oscillations of the correlation function g(x₁). The above qualitative observations concerning representative configurations are quantified in the functions representing the average structure in the transverse and in the longitudinal direction, and in the orientational order parameter discussed below.

The density profiles. In Fig. 11 we present the density profiles for a few temperatures for a slit incommensurate with the lamellar period and containing 8 undulating stripes (L₂ = 32), and for the slit commensurate with the lamellar period and containing 9 straight stripes (L₂ = 34). For T* < 0.5 the maxima of ρ(x₂) for x₂ approaching the center of the slit decrease for L₂ = 32 but not for L₂ = 34, since the stripe undulations are present only in the former case (see Fig. 10 and the discussion in the previous paragraph). Interestingly, for L₂ = 32 the profile very weakly depends on T* for T* ≤ 0.5. Note that these temperatures correspond to T* < T_max*(L₂ = 30), where the L-like structure is expected based on thermodynamic considerations and on the snapshots in Fig. 10. On the other hand, for L₂ = 34 the profiles for T* = 0.25 and T* = 0.4 are very similar, but for T* = 0.5 the degree of order decreases significantly. However, T_max*(L₂ = 34) ≈ 0.495, and for T* = 0.5 > T_max*(L₂ = 34) the topological defects start to destroy the order. Finally, T* = 1 is above the temperature of the second maximum of the heat capacity, and we expect the isotropic F structure based on thermodynamic considerations and based on the snapshots in Fig. 10. Indeed, for both slits the oscillations in the slit center are very weak. Moreover, the density profiles near the boundaries are nearly the same, because the wall-induced ordering of the stripes is not influenced by the second wall when the slit interior has no periodic order.

Fig. 11 The density profiles in the direction perpendicular to the walls, averaged over the longitudinal direction (see (9)), for several temperatures. Left panel: L₂ = 32, right panel: L₂ = 34. The color code is the same in both panels. The walls are attractive with h* = −1.

Interestingly, the density profiles averaged over the longitudinal direction resemble strongly the density profiles of the one-dimensional cluster phase.²⁵

The correlations in the longitudinal direction. The correlation function g(x₁) is shown in the upper panel of Fig. 12 for L₂ = 32 at several temperatures. It is interesting that g(x₁) has a very regular shape, despite rather complex structure of individual configurations (Fig. 10). In agreement with the analysis of the representative configurations, the damping of the oscillations with the distance decreases, and the apparent wavelength increases substantially on cooling the system.

Fig. 12 Upper panel: The correlation function g(x₁) for L₂ = 32. Lower panel: Fourier transforms S(q) for L₂ = 32. The pseudo-periodicity parameter x^p₁(T) ≡ 2π/q_max takes the values: x^p₁ = 23.9 (for T* = 0.25); x^p₁ = 20.3 (for T* = 0.40); x^p₁ = 15.6 (for T* = 0.50); and x^p₁ = 6.6 (for T* = 1.00).

We have carried out a Fourier-like analysis of the correlation function g(x₁). The Fourier integral applied to the lattice was taken as:

Examples of the results can be found in the lower panel of Fig. 12. Defining q_max as the value of q that maximizes S(q) we can compute a length x^p₁ = 2π/q_max, which characterizes the typical wavelength of the oscillations in g(x₁). In the caption of Fig. 12 we collect the corresponding values for slits with L₂ = 32 and different temperatures. The results for x^p₁ are consistent with the position of the second maximum of g(x₁). From the trends of the correlation functions g(x₁) shown in Fig. 12 one could have anticipated an ordered structure at the ground state, with the periodicity in the direction e₁, as it was already presented in Fig. 3.

In Fig. 13 we compare the correlation function g(x₁) for L₂ = 30 and L₂ = 32 for temperatures that correspond to the L, ML and F pseudo-phases (see Fig. 7). As expected, in the L phase g(x₁) quickly vanishes for L₂ = 30, and for L₂ = 32 exhibits damped oscillations with a very large decay length (see Fig. 12 for low T*). When the temperature approaches T_K*(L₂), damped oscillations are present in both cases, but the decay length and the period are both larger for L₂ = 32, indicating that the corrugations induced by the size incommensurability are still important for T* = 0.55 that is slightly above T_K*(L = 32). Near the temperature T* ≃ 0.6 the shape of g(x₁) becomes nearly the same for L₂ = 30 and L₂ = 32. This close similarity remains present for higher temperatures, indicating that no trace of the L structure is left. The global correlation function in the ML and F pseudo-phases is not influenced by the changes of L₂, as long as the same number of stripes is present between the walls.

Fig. 13 The correlation function g(x₁) (see (10) and (11)) for L₂ = 30 (solid lines) and L₂ = 32 (dashed lines). The length of the system in the longitudinal direction in the MC simulations was L₁ = 2520. Panels from (a) to (d) correspond to T* = 0.5, 0.55, 0.625 and 1 respectively. The thermodynamic considerations and Fig. 7, 10 and 11 indicate stability of the following phases: (i) for T* = 0.5 the L pseudo-phase when L₂ ≤ 30, and the beginning of the L–ML crossover for L₂ = 32, (ii) for T* = 0.55 the crossover between the L and ML pseudo-phases, (iii) for T* = 0.625 the ML pseudo-phase, and (iv) the F pseudo-phase for T* = 1.

The orientational order parameter. The O_L computed for the slits of widths that are commensurate and incommensurate with the lamellar period respectively is shown in Fig. 14 for L₂ = 30 and L₂ = 32 together with the O_L obtained for the bulk.¹¹ For T* < 0.5 and L₂ = 30 we find O_L ≃ 1, indicating the perfect lamellar order. This behavior is very similar to the bulk L phase. For L₂ = 32 the zig-zag structure of the lamellar stripes leads to O_L ≈ 0.8 for T* < 0.5.

Fig. 14 The orientational OP O_L computed for L₁ = 180 and for the commensurate (L₂ = 30, L₂ = 34) and incommensurate (L₂ = 32, L₂ = 36) slit widths as a function of temperature. The thick solid line is the O_L computed for the bulk system.¹¹

For T* > 0.5 the OP shown in Fig. 14 differs significantly from the OP in the bulk. In the bulk O_L(T*) very weakly decreases with temperature for T* < 0.725, and then decreases rapidly to zero with the inflection point at T* ≈ 0.75. In contrast, in the slit it decreases significantly for increasing T* when 0.5 < T* < 0.9, and for T* > 0.9 a slow decrease takes place. Comparison of Fig. 14 and 7 suggests that the O_L(T*) starts to decrease when the crossover region between the L and ML structures is entered by heating the system. Thus, in contrast to the bulk, O_L(T*) decreases with increasing temperature in the whole region that we identify with the ML pseudo-phase.

The confinement-reduced orientational order in the ML phase is counter-intuitive, and it is due to the definition of the order parameter. Our simulations in the bulk typically show stripes that consist of short segments parallel to two main lattice directions, and the number of lamellar segments in one of the main directions of the lattice is very small, which makes O_L ≈ 1 (see (12)). In a slit, however, the geometrical defects in the near-boundary stripes induce a similar number of segments in directions e₂, and e₃; this number is smaller that the number or lamellar segments in direction e₁, but by no means negligible, particularly for slit widths incommensurate with the lamellar periodicity. By visual inspection of Fig. 10 we can see that for T* ≥ 0.556 the stripes are parallel to the boundaries only close to the walls, whereas away from the boundaries the stripes consist of short segments parallel to different e_i. It remains to be verified if in the off-lattice system the orientational order in the ML phase is reduced by confinement too.

Note that the O_L(T*) lines computed for L₂ and L₂ + 2 coincide above certain temperature, in agreement with the identical shapes of the global correlation function g(x₁) (see Fig. 13). This temperature is marked by big triangles in Fig. 7. This is another confirmation that the mismatch between L₂ and the local lamellar period is irrelevant when there is no periodic order.

The ML–F phase transition in the bulk was identified with the inflection point of the O_L(T*). The shape of O_L in the slit is much different. There is no clearly distinguished temperature marking the transition between the anisotropic and isotropic structures when the OP changes gradually in a large temperature interval. We identify the pseudo phase transition between these phases in the slit with the second maximum of c_μ, but in fact the crossover region between the two structures is very broad. The high-T plateau of the O_L is reached at T* ≃ 1.4 that is much higher than in the system with PBC, i.e. at T* ≃ 0.8. This orientational ordering in the confined F-like structure is associated with the layering near the walls. As shown in Fig. 11, three maxima of the density profile in the transverse direction are present near each wall at T* = 1, and this periodic near-surface structure contributes to O_L.

The pseudo-transition ML–F occurs at the temperature T* ≃ 0.7 that is somewhat lower than the temperature of the ML–F transition in the bulk, estimated as T* ≃ 0.75. Note that this transition in the bulk is continuous. An apparent continuous transition in a slit corresponds to the temperature at which the bulk correlation length and the slit width are comparable and it happens before the bulk transition occurs. Thus, lower temperature of the pseudo-transition in the slit agrees with the expected behavior.

4 Conclusions

Despite the fact that in a two-dimensional slit there are no true thermodynamic phase transitions, we have found significantly different structures for different parts of the temperature – slit width structure-diagram, (T*,L₂). The structures with and without the periodic order lie on the low- and the high-temperature side of the Kelvin-like line T = T_K(L₂), respectively. Thus, the melting of the periodic lamellar structure occurs at the temperature that is much higher than in the unrestricted monolayer, and increases for decreasing L₂.

Importantly, in the periodic structure the stripes integrity is rarely broken. When the distance between the boundaries and the period of the structure are commensurate or incommensurate, the stripes are straight or show collective corrugation, respectively. In contrast, in the absence of the periodic order there are many topological defects. In this case changes of the distance between the boundaries from commensurate to incommensurate have no effect on the self-assembled structure.

Our results show that for T < T_K(L₂) the template in the form of parallel lines has a significant ordering effect on the structure occurring spontaneously on the length scale smaller than the wall separation. Thus, if ordered patters with few defects are needed, parallel boundaries such that T < T_K(L₂) should help. If isotropic or anisotropic porous medium is desired, then T > T_K(L₂) should be chosen. We have verified that the low-T maxima of the heat capacity measured for L₂ commensurate with the period of the stripe pattern lie on the T = T_K(L₂) line, and that L₂ at this line is comparable with the range of the solvation pressure. We hope that our findings can help in designing and interpreting experimental studies of systems self-assembling into stripe patterns in 2D.

5 Appendix

5.1 Grand canonical Monte Carlo simulation

We have extensively used the grand canonical Monte Carlo simulation through the paper. The elementary step for these simulations in our lattice system was:

(i) Choosing randomly one of the sites, i of the system

(ii) Choose one of the two possible states for that site [small rho, Greek, circumflex]_i = 0,1:

with probabilities proportional to exp[−βh_i([small rho, Greek, circumflex]_i)], where h_i includes the contributions of the site i to the Hamiltonian given in eqn (2).

We define a cycle as a number of elementary steps equal to the number of sites of the system.

5.2 Computation of the ground state configurations

In order to compute the ground state of the system, we have made use of the parallel tempering technique. In our procedure we considered typically n_β = 151 different values of the reverse temperature, β, as:

β_k* = k × δβ*, with k = 0, 1, 2,…,n_β − 1, with δβ* = 0.10;

whereas the chemical potential was fixed at μ* = 6.

After completion of each simulation cycle for all the n_β cases, we perform attempts of configuration interchange between (n_β − 1)/2 pairs of neighbor states, alternating between pairs (2k, 2k + 1) and (2k + 1, 2k + 2), with k = 0, 1,…,(n_β − 3)/2.

After a given number of simulation cycles (that depends on the system lengths), the energies of the lowest temperatures reach a clear plateau, with just occasional elementary excitations. The configurations in this region correspond to minimum energy for a given system size, given in terms of the lengths L₂ (slit width), and L₁ (length in the direction where periodic boundary conditions are applied).

The overall GS for a given slit width L₂, is attained by considering the lateral length L_1,min which minimizes H[{[small rho, Greek, circumflex](x)}]/(L₁L₂). The result is checked by testing that the same minimum energy and configuration is obtained when considering L₁ = 2L_1,min.

The required simulation lengths to achieve reliable estimates of the energy minima for each system size L₁ × L₂, ranges from ∼5 × 10⁴ to ∼10⁶ cycles. In general, more cycles are required as the slit width increases.

5.3 Thermodynamic properties

The thermodynamic properties were computed using the parallel tempering technique for systems with different slit widths: L₂ = 2, 3, 4,…,40. We took for all the cases L₁ = 120. For systems with L₂ = 4n + 2 (with n being an integer number) we considered n_β = 401 values of the inverse temperature β, whereas for the other slit widths we took n_β = 501. In both cases the particular values of β are given β_k* = kδβ*, with δβ* = 0.01 and k = 0, 1,…,n_β − 1.

The simulations runs comprised typically 10⁷ cycles of equilibration plus 10⁷ cycles of sampling. These relative long simulations were required to get a reliable sampling of the properties at low temperature.

For the bulk system, whose properties were required to compute the solvation pressure, we used also the parallel tempering scheme on a system of 120² sites, using n_β = 501, δβ* = 0.01, with about 3.2 × 10⁶ cycles for equilibration, and equal simulation length for the sampling.

Replica exchanges were attempted for both, confined and bulk systems, every five simulation cycles following the same scheme described previously.

5.4 Correlation functions

The density profiles and the correlation functions presented in Section 3.3 were computed using standard GCMC simulations. In order to minimize the possible effects of the periodic boundary conditions, we considered large system size in the longitudinal direction, namely L₁ = 2520.

These simulations required some care to ensure that the systems were equilibrated before computing the correlation functions. The number of cycles to reach a good equilibration of the systems depends on the temperature and the slit width.

For L₂ = 32 and T* ≤ 0.50, the simulation runs comprised ∼2 × 10⁷ cycles of equilibration, and the same length for the sampling part. Shorter simulations (by one order of magnitude) could be used for T* > 0.50, and for L₂ = 30 at all T*. Data for further computation of the correlation functions and the density profiles were collected from configurations every 10³ cycles.

5.5 The orientational order parameter

The orientational order parameter is computed following the same strategies introduced in ref. 11. Firstly we establish which sets of triangles of occupied nearest-neighbor sites fulfill the definition of lamellar triangles as defined in ref. 11, and determine their corresponding orientation on the lamellar stripe. Then, we count how many lamellar triangles are oriented in each of the three main lattice directions.

In order to focus on the ordering in the region not directly influenced by the presence of the walls, only sites with coordinates 2 < x₂ < L₂ − 1 are considered in the calculation.

The total length of so defined lemallar segments is S = S₁ + S₂ + S₃, where S_i is the number of lamellar triangles with orientation e_i (see Fig. 1). The order parameter is defined as O_L = 1 − 3 min (S_i/S) as it was done in the bulk (ref. 11). If the three directions are almost equiprobable we will have O_L ≈ 0, whereas if there is a preferential orientation for the lamellar stripes, O_L will get a large value (relatively close to one) on approaching the ground state that will depend on the corrugation of the inner lamellar stripes for the corresponding slit width.

Acknowledgements

JP acknowledges the financial support by the National Science Center under Contract Decision No. DEC-2013/09/N/ST3/02551. JP received funding for the preparation of the doctoral dissertation from the National Science Center in the funding of PhD scholarships on the basis of the Decision No. DEC-2014/12/T/ST3/00647. AC acknowledges the financial support by the National Science Center grant 2015/19/B/ST3/03122. Financial support from the Dirección General de Investigación Cientfica y Técnica under Grant FIS2013-47350-C5-4-R is also acknowledged.

References

1. R. P. Sear and W. M. Gelbart, J. Chem. Phys., 1999, 110, 4582.
2. A. Stradner, H. Sedgwick, F. Cardinaux, W. Poon, S. Egelhaaf and P. Schurtenberger, Nature, 2004, 432, 492.
3. S. Rozovsky, Y. Kaizuka and J. T. Groves, J. Am. Chem. Soc., 2005, 127, 36.
4. A. I. Campbell, V. J. Anderson, J. S. van Duijneveldt and P. Bartlett, Phys. Rev. Lett., 2005, 94, 208301.
5. D. G. Barci and D. A. Stariolo, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 075437.
6. D. G. Barci, A. Mendoza-Coto and D. A. Stariolo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062140.
7. A. Shukla, E. Mylonas, E. D. Cola, S. Finet, P. Timmins, T. Narayanan and D. I. Svergun, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 5075.
8. M. Seul and D. Andelman, Science, 1995, 267, 476.
9. A. Ciach, J. Pȩkalski and W. T. Góźdź, Soft Matter, 2013, 9, 6301.
10. A. Ciach and W. T. Góźdź, Condens. Matter Phys., 2010, 13, 23603.
11. N. G. Almarza, J. Pȩkalski and A. Ciach, J. Chem. Phys., 2014, 140, 164708.
12. A. Imperio and L. Reatto, J. Phys.: Condens. Matter, 2004, 18, S2319.
13. A. Imperio and L. Reatto, J. Chem. Phys., 2006, 124, 164712.
14. A. Imperio and L. Reatto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 040402.
15. A. J. Archer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 031402.
16. J. Pȩkalski, A. Ciach. and N. G. Almarza, J. Chem. Phys., 2014, 140, 114701.
17. B. Chacko, C. Chalmers and A. J. Archer, J. Chem. Phys., 2015, 143, 244904.
18. Y. Zhuang, K. Zhang and P. Charbonneau, Phys. Rev. Lett., 2016, 116, 098301.
19. A. Ciach, W. T. Góźdź and A. Perera, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 021203.
20. Y. Zhuang and P. Charbonneau, J. Phys. Chem. B, 2016, 120, 6178.
21. H. Shin, G. M. Grason. and C. D. Santangelo, Soft Matter, 2009, 5, 3629.
22. B. A. Lindquist, R. B. Jadrich and T. M. Truskett, Soft Matter, 2016, 12, 2663.
23. M. Edelmann and R. Roth, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2016, 93, 062146.
24. A. Ciach and W. T. Gozdz, J. Phys.: Condens. Matter, 2016, 28, 244004.
25. J. Pȩkalski, A. Ciach and N. G. Almarza, J. Chem. Phys., 2013, 138, 144903.
26. A. Ciach and W. T. Góźdź, J. Phys.: Condens. Matter, 2016, 28, 244004.
27. D. F. Schwanzer and G. Kahl, J. Phys.: Condens. Matter, 2010, 22, 415103.
28. H. Hasegawa and T. Hashimoto, Polymer, 1992, 33, 475.
29. A. Horvat, G. J. A. Sevink, A. V. Zvelidovsky, A. Krekhov and L. Tsarkova, ACS Nano, 2008, 2, 1143.
30. M. R. Hammond, E. Cochran, G. H. Fredrickson and E. J. Kramer, Macromolecules, 2005, 38, 6575.
31. G. M. Whitsides and M. Boncheva, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 4769.
32. M. Rey, M. A. Fernández-Rodríguez, M. Steinacher, L. Scheidegger, K. Geisel, W. Richtering, T. M. Squires and L. Isa, Soft Matter, 2016, 12, 3545.
33. A. Saric and A. Cacciuto, Soft Matter, 2013, 9, 6677.
34. A. Imperio, L. Peatto and S. Zapperi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 021402.
35. D. Antelmi, P. Kékicheff and P. Richetti, J. Phys. II, 1995, 5, 103–112.
36. M. Tasinkevych and A. Ciach, J. Chem. Phys., 2001, 115, 8705.
37. B. Yu, P. Sun, T. Chen, Q. Jin, D. Ding, B. Li and A.-C. Shi, Phys. Rev. Lett., 2006, 96, 138306.
38. B. Yu, Q. Jin, D. Ding, B. Li and A.-C. Shi, Macromolecules, 2008, 41, 4042–4054.
39. P. Chi, Z. Wang, B. Li and A.-C. Shi, Langmuir, 2011, 27, 11683–11689.
40. S. Perkin and J. Klein, Soft Matter, 2013, 9, 10438.
41. L. Nicolao, A. Mendoza-Coto and D. A. Stariolo, J. Phys.: Conf. Ser., 2016, 686, 012005.
42. R. Evans, J. Phys.: Condens. Matter, 1990, 2, 8989.
43. V. Babin, A. Ciach and M. Tasinkevych, J. Chem. Phys., 2001, 114, 9585.
44. D. de las Heras, E. Velasco and L. Mederos, Phys. Rev. Lett., 2005, 94, 017801.
45. A. Ciach, Prog. Colloid Polym. Sci., 2004, 129, 40.
46. J. Pȩkalski, A. Ciach and N. G. Almarza, J. Chem. Phys., 2015, 142, 014903.

---

Footnotes

† Noé G. Almarza passed away unexpectedly during preparation of this manuscript for publication. This paper is dedicated to his memory.

‡ Deceased.

---