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

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 L2 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 L2. This structural order is lost at T = TK(L2) depending on L2 according to a Kelvin-like equation. Above the Kelvin temperature TK(L2) 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 L2 optimal for the stripe pattern the heat capacity as a function of temperature takes the maximum at T = TK(L2).


Introduction
Competing interactions, in particular short-range attraction and long-range repulsion (SALR), are present in many biological and soft matter systems [1][2][3][4] as well as in magnetic films with competing ferromagnetic and dipolar forces. 5,6The 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][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,20Much higher temperature is necessary for disassembly of the clusters.][26][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][29][30] by particles adsorbed at solid surfaces or at liquid interfaces, 1,31,32 on elastic membranes or embedded in lipid bilayers. 3,33][13][14][15][16][17] The 2D self-assembly in continuous SALR models was studied by the density functional theory 15 and by Monte Carlo (MC) [12][13][14]34 methods. Rectly lattice SALR models (LSALR) were introduced and studied in the mean-field approximation (MF) and by MC This journal is © The Royal Society of Chemistry 2016 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 blockcopolymers.[37][38][39][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 wor 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.36 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 r B 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,30The 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 r = 1/2 that is optimal for the lamellar structure. 11The line of continuous transitions terminates at two tricritical points, one with the density r tcp E 0.4, and the other one with r tcp E 0.6, where the transition becomes first-order. 11The ML phase is analogous to the nematic phase of stripes found in magnetic films. 5Recent studies of the 2D stripe-forming systems 41 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 r B 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][43][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 2 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 2 .Section 4 contains our conclusions.

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, s, comparable with the diameter of the adsorbed particles.The lattice sites are denoted by x = x 1 e 1 + x 2 e 2 (x 1 ,x 2 ), where e 1 , e 2 and e 3 = e 2 À e 1 are the unit lattice vectors on the triangular lattice, i.e. |e 1 | = |e 2 | = |e 1 À e 2 | = 1 (in s-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 1 , and separated by L 2 layers of sites along the direction e 2 .Thus, 1 r x 2 r L 2 .In the direction e 1 we assume periodic boundary conditions (PBC).In simulations we assume 1 r x 1 r L 1 , where L 1 c L 2 , and x 1 = L 1 + 1 is identified with x 1 = 1.The system under consideration is shown in Fig. 1.
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 m 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 nearestneighbor 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.5s B l D .Since in various solvents with weak ionic strength l D B 1-100 nm, the model is suitable for charged molecules, nanoparticles or globular proteins.The interaction between the occupied sites x and x + Dx of the form described above is given by VðDxÞ ¼ ÀJ 1 for jDxj ¼ 1; ðnearest neighborsÞ þJ 2 for jDxj ¼ 2; ðthird neighborsÞ 0 otherwise; where ÀJ 1 and J 2 represent the attraction well and the repulsion barrier respectively.The thermodynamic Hamiltonian has the form where P x denotes the summation over all lattice sites, the microscopic density at the site x is r(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 {r(x)}.The probability of a particular microscopic state has the form where b = 1/(k B T) and k B is the Boltzmann constant.The grand potential is expressed in terms of the grand statistical sum in the standard way We choose J 1 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 2 /J 1 .Reduced temperature is defined as: T* = k B T/J 1 .Along the paper we have used the same simulation techniques as used for this model in the bulk, 11 that is: Metropolis Monte Carlo, thermodynamic integration and parallel tempering technique.More details about the methodology are given in the Appendix.

The ground state
In the previous work 16 we determined the ground state (GS) of the model and the MF phase diagram.Later we computed 11 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 m* = 6, for which the L phase is stable. 11,16In the most of the cases we shall assume h* = À1, which is a moderate wallfluid interaction.
The grand potential for T = 0 reduces to the minimum of the grand canonical Hamiltonian per site % H 0 H[{r(x)}]/(L 1 L 2 ).The ground states for m* = 6 and different wall-wall separations L and values of h* were computed as follows: for fixed values of m*, L 2 , and h*, we consider different values of the lateral length L 1 , and determine by means of parallel tempering MC simulations their corresponding minima, % H 0 (L 1 ).The GS configuration and the GS energy, at the given conditions (m*,L 2 ,h*), are taken as those corresponding to the value L 1 which minimizes % H 0 .
Fig. 1 Sketch of the two-dimensional triangular lattice with the slit boundary conditions.L 1 is the length of the system in direction e 1 (in this direction periodic boundary conditions are assumed), whereas L 2 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 1 , where J 1 is the nearest-neighbor attraction and h is the interaction energy between the particle in the boundary layer and the confining wall.
This journal is © The Royal Society of Chemistry 2016

Thermodynamic properties
The heat capacity.We start the investigation of the thermodynamic properties of the model at T* 4 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 where O bulk = Àk B T ln X bulk is the grand potential of the bulk system with the same volume, i.e. with PBC in both directions, g i is the wall-fluid surface tensions at the i-th wall, and C(L 2 ) is the interaction between the confining walls per unit length induced by the confined fluid (solvation pressure).
In order to compute C(L 2 ) we have applied the MC thermodynamic integration procedure for fixed chemical potential, which was previously used for the bulk system. 11The grand potential per lattice site was computed using the following relation bOðb; mÞ where 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 2 .The term Àln 2 results from the assumption that at infinite temperature our model acts like an ideal lattice gas.To obtain C(L 2 ) out of o ex = (O À O bulk )/L 1 , we note that these two quantities at fixed temperature differ only by a constant (see (7)).Since lim we have estimated the constant by the value of o ex for which the curve reaches a plateau.

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 1 , The longitudinal correlation function.The structure in the longitudinal direction, e 1 , is studied by defining a correlation function along each line parallel to the boundary, and average it over all the lines located at 1 r x 2 r L 2 .This global correlation function along the confining walls allows us to compute the average period and decay length in direction e 1 .
We define the global density correlation function along the direction e 1 as where r(x 0 1 ) (r(x 0 1 ,1),. ..,r(x 0 1 ,L 2 )), is a vector containing the densities of the sites belonging to a compact line of sites in the direction e 2 (see Fig. 2), Á is the standard scalar product, and hÁ Á Ái indicates an average over different system configurations.From g 0 (x 1 ) we can define a normalized correlation function: g(x 1 ) g 0 (x 1 )/g 0 (0). ( The orientational order parameter.We have made use of an orientational order parameter O L similar to that of the bulk system. 11In short, we count the total lengths (S 1 , S 2 , S 3 ) 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 1 C S 2 C S 3 , and therefore O L should vanish in the thermodynamic limit, whereas in the bulk L phase there are lamellar segments only in two directions, 11 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 o x 2 o L 2 À 1 were considered in the calculation.

The ground state
Long-range order of the system for given m and L 2 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 ÀmN, 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 2 , and on the wall-fluid interaction h*.
If L 2 = 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 2 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 2 and h*, and there are many GS configurations in the (L 2 ,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 2 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 1 direction form an equilateral triangle (Fig. 4).
We found that for slit widths L 2 4 13 the distance L determines the shape and the position of the defects, hence having L 2 one can also tell the size of the unit cell in the e direction, L 1 , which is given by the following formula We cannot tell, however, if for large wall-wall separations (e.g.L 2 4 30) L 1 is still determined uniquely by L 2 .
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 2 = 16 and the neutral, h* = 0, and strongly repulsive, h* = 3, boundaries.Note that the period in direction e 1 is L 1 = 12 for h* = À1 (see Fig. 3), and L 1 = 16 for h* = (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 2 = 4n one can remove the defects by changing the particle-wall interaction from attractive to strongly repulsive.
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. 36The 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 2 a new stripe is introduced into the system.

Thermodynamics
The heat capacity.As in the GS, we take h* = À1 and m* = 6, but we limit ourselves to the pore widths L 2 = 6, 10, 14, 18,. ..,38 (i.e.those sizes whose GS present straight stripes).In Fig. 6 we show c m as a function of the temperature.For the narrowest slits only one maximum of c m (T) appears, whereas for L 2 Z 30 there are two maxima.One of them is located at T* C 0.7.The other one appears at lower temperature that decreases significantly with increasing L 2 .
We do not expect true thermodynamic phase transitions in slit systems at T* 4 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, 25 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 11 for m* = 6 at T* C 0.75 that is close to the location of the high-T maximum of c m .On the other hand, MC simulations 11 for L 1 = L 2 = 120 show that in the bulk the ML-L transition occurs at T bulk * C 0.25 that is much lower than the temperature at which c m takes the low-T maximum for L 2 r 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. 42The thinner the capillary, the higher the temperature of the capillary condensation.The phenomenon is described by the Kelvin equation T K *(L 2 ) = T bulk * + a/L 2 , 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,45Similar behavior, i.e. capillary lamellarization or capillary smectization was observed in surfactant solutions or in liquid crystals. 43,44n 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 2 ) corresponding to the maximum of c m .In addition, we plot the best fit of the data to the equation We obtain a fair agreement between T max *(L 2 ) 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 2 ).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 m 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.
The Kelvin-like equation was fitted to the maxima of the heat capacity only for L 2 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 2 are expected in the SALR system.We verified that the geometrical defects present for L 2 incommensurate with the period of the striped structure influence the heat capacity too.For this reason we did not consider the values of L 2 for which the simple eqn ( 13) cannot be valid.
The solvation pressure.A typical shape of C(L 2 ) is shown in Fig. 8 for T* = 0.5.The shape of C(L 2 ) is very similar to that obtained for the 1D LSALR model 46 by exact calculations, and to the C measured experimentally in surfactant mixtures. 35It 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 2 = 4n + 2, with n = 1, 2,. ..,5 correspond to local minima of C(L 2 ) for L 2 r 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 2 4 22 the minima are located at L 2 = 31 and L 2 = 35, which do not correspond to the wall separation optimal for the straight lamella structure.For T* 4 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 C(L 2 ) obtained by exact calculations is a real number. 46ince in a lattice model C(L 2 ) is meaningful only for integer L 2 , we did not try to fit C(L 2 ) to a smooth curve.In the 1D case the exact C(L 2 ) exhibits exponentially damped oscillations for large L 2 , and similar asymptotic behavior, i.e. exponentially damped oscillations for large L 2 , can be expected in this case.However, in 1D the asymptotic form of C(L 2 ) deviates significantly from the exact solution for relatively large L 2 , and by analogy we expect that for L 2 o 30 no simple analytical formula can describe C(L 2 ) in Fig. 8.
At large separations L 2 the interaction between the confining walls follows from the mismatch between L 2 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 C40, which is close to L 2 C 38 at the boundary of the shaded region in Fig. 7 for this temperature.We conclude that the range of C(L 2 ) increases significantly with decreasing temperature, and its very rough estimate is given by the high-T boundary of the shaded region in Fig. 7.

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 2 = 30, 31, 32, 33, 34, with L 1 = 120, and at T* = 0.25 are shown.For L 2 = 30 and L 2 = 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 o L 2 o 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 1 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 2 = 30 and L 2 = 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) E 0.53.In the case of L 2 = 30 almost no geometrical defects are present, whereas for L 2 = 32 more geometrical defects than in the GS appear.They lead to small distortions of the ordered zig-zag structure.For temperatures T* o T K *(L 2 ) (see ( 13)) the system seems to exhibit large correlation at relatively large distances in the e 1 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* o T K *(L 2 ).The second and third panels from the bottom correspond to temperatures between the two maxima of c m .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 2 = 32 and L 2 = 30 respectively.The anisotropic structure with topological defects characterizes the ML phase.For T* E 0.625 that is close to the minimum of the heat capacity one cannot distinguish the structure for L 2 = 30 and L 2 = 32 by visual inspection.This indicates that the topological defects that appear independently of L 2 dominate.The two top panels correspond to temperatures T* = 0.806, and T* = 1.00 that are above the high-T* maximum of c m .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.
The above inspection of the representative configurations shows consistency of the thermodynamic and structural properties.Before looking at the structural functions, r(x 2 ) and g(x 1 ), 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  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 m *, i.e. above the dashed line in Fig. 7.
completely filled with particles for L 2 = 30, and there are alternating occupied and empty segments for L 2 = 32, because of the stripe corrugation.The density averaged over the longitudinal direction for x 2 = 9 is smaller for L 2 = 32 than for L 2 = 30.Further away from the boundary this difference in the average density increases.Moreover, for L 2 = 32 the alternating occupied and empty segments along the direction e 1 lead to an oscillating correlation function.The period of the oscillations at very low T* is determined by L 2 , as in the GS.For T* C T max *(30) E 0.53 the topological defects lead to additional alternating occupied and empty segments in each line parallel to the boundaries for L 2 = 32, and to appearance of such segments for L 2 = 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 1 ).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 2 = 32), and for the slit commensurate with the lamellar period and containing 9 straight stripes (L 2 = 34).For T* o 0.5 the maxima of r(x 2 ) for x 2 approaching the center of the slit decrease for L 2 = 32 but not for L 2 = 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 2 = 32 the profile very weakly depends on T* for T* r 0.5.Note that these temperatures correspond to T* o T max *(L 2 = 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 2 = 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 2 = 34) E 0.495, and for T* = 0.5 4 T max *(L 2 = 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.
Interestingly, the density profiles averaged over the longitudinal direction resemble strongly the density profiles of the onedimensional cluster phase. 25he correlations in the longitudinal direction.The correlation function g(x 1 ) is shown in the upper panel of Fig. 12 for L 2 = 32 at several temperatures.It is interesting that g(x 1 ) 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.
We have carried out a Fourier-like analysis of the correlation function g(x 1 ).The Fourier integral applied to the lattice was taken as: Fig. 11 The density profiles in the direction perpendicular to the walls, averaged over the longitudinal direction (see ( 9)), for several temperatures.This journal is © The Royal Society of Chemistry 2016 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 1 = 2p/q max , which characterizes the typical wavelength of the oscillations in g(x 1 ).In the caption of Fig. 12 we collect the corresponding values for slits with L 2 = 32 and different temperatures.The results for x p 1 are consistent with the position of the second maximum of g(x 1 ).From the trends of the correlation functions g(x 1 ) shown in Fig. 12 one could have anticipated an ordered structure at the ground state, with the periodicity in the direction e 1 , as it was already presented in Fig. 3.
In Fig. 13 we compare the correlation function g(x 1 ) for L 2 = 30 and L 2 = 32 for temperatures that correspond to the L, ML and F pseudo-phases (see Fig. 7).As expected, in the L phase g(x 1 ) quickly vanishes for L 2 = 30, and for L 2 = 32 exhibits damped oscillations with a very large decay length (see Fig. 12 for low T*).When the temperature approaches T K *(L 2 ), damped oscillations are present in both cases, but the decay length and the period are both larger for L 2 = 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* C 0.6 the shape of g(x 1 ) becomes nearly the same for L 2 = 30 and L 2 = 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 2 , as long as the same number of stripes is present between the walls.
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 2 = 30 and L 2 = 32 together with the O L obtained for the bulk. 11For T* o 0.5 and L 2 = 30 we find O L C 1, indicating the perfect lamellar order.This behavior is very similar to the bulk L phase.For L 2 = 32 the zig-zag structure of the lamellar stripes leads to O L E 0.8 for T* o 0.5.
For T* 4 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* o 0.725, and then decreases rapidly to zero with the inflection point at T* E 0.75.In contrast, in the slit it decreases significantly for increasing T* when 0.5 o T* o 0.9, and for T* 4 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 E 1 (see (12)).In a slit, however, the geometrical defects in the nearboundary stripes induce a similar number of segments in directions e 2 , and e 3 ; this number is smaller that the number or lamellar segments in direction e 1 , 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* Z 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 2 and L 2 + 2 coincide above certain temperature, in agreement with the identical shapes of the global correlation function g(x 1 ) (see Fig. 13).
Fig. 13 The correlation function g(x 1 ) (see (10) and ( 11)) for L 2 = 30 (solid lines) and L 2 = 32 (dashed lines).The length of the system in the longitudinal direction in the MC simulations was L 1 = 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 2 r 30, and the beginning of the L-ML crossover for L 2 = 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.This temperature is marked by big triangles in Fig. 7.This is another confirmation that the mismatch between L 2 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 m , 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* C 1.4 that is much higher than in the system with PBC, i.e. at T* C 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* C 0.7 that is somewhat lower than the temperature of the ML-F transition in the bulk, estimated as T* C 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.

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 2 ).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 2 ), 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 2 .
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 o T K (L 2 ) 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 o T K (L 2 ) should help.If isotropic or anisotropic porous medium is desired, then T 4 T K (L 2 ) should be chosen.We have verified that the low-T maxima of the heat capacity measured for L 2 commensurate with the period of the stripe pattern lie on the T = T K (L 2 ) line, and that L 2 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.

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 ri = 0,1: with probabilities proportional to exp[Àbh i (r 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.

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 After completion of each simulation cycle for all the n b cases, we perform attempts of configuration interchange between (n b À 1)/2 pairs of neighbor states, alternating between pairs (2k, 2k + 1) and (2k + 1, 2k + 2), with k = 0, 1,. ..,(n b À 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 2 (slit width), and L 1 (length in the direction where periodic boundary conditions are applied).
The overall GS for a given slit width L 2 , is attained by considering the lateral length L 1,min which minimizes H[{r(x)}]/ (L 1 L 2 ).The result is checked by testing that the same minimum energy and configuration is obtained when considering L 1 = 2L 1,min .
The required simulation lengths to achieve reliable estimates of the energy minima for each system size L 1 Â L 2 , ranges from B5 Â 10 4 to B10 6 cycles.In general, more cycles are required as the slit width increases.

Thermodynamic properties
The thermodynamic properties were computed using the parallel tempering technique for systems with different slit widths: The simulations runs comprised typically 10 7 cycles of equilibration plus 10 7 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 2 sites, using n b = 501, db* = 0.01, with about 3.2 Â 10 6 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.

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 1 = 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 2 = 32 and T* r 0.50, the simulation runs comprised B2 Â 10 7 cycles of equilibration, and the same length for the sampling part.Shorter simulations (by one order of magnitude) could be used for T* 4 0.50, and for L 2 = 30 at all T*.Data for further computation of the correlation functions and the density profiles were collected from configurations every 10 3 cycles.

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 o x 2 o L 2 À 1 are considered in the calculation.
The total length of so defined lemallar segments is S = S 1 + S 2 + S 3 , 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 E 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.

Fig. 3
Fig.3The ground state (T = 0) configurations for L 2 = 4n À 1 (left column) L 2 = 4n (center column) and L 2 = 4n + 1 (right column), where n = 2, 3, 4, 5.The confinement attracts the particles with h* = À 1. L 1 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
Fig.4The defects are in the vertices of the equilateral triangle.

Fig. 6
Fig. 6 Reduced heat capacity, c m *(T) c m (T)/k B for m* = 6 and h* = À1, and for different pore widths.c m (T) exhibits two maxima for L 2 4 30, and an inflection point for L 2 = 30.

Fig. 7
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 m (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 m* = 6 and h* = À1, and the pore width is L 2 = 4n + 2 with integer n (see Fig. 6).The solid line is the best fit to the Kelvin-like equation T K *(L 2 ) (see (13)), where the fitting parameters are T b * = 0.28 and a = 7.422.Note that T b * is close to T bulk * C 0.25 associated with the transition between the L and ML phases in the unrestricted monolayer.The shaded region around T K *(L 2 ) corresponds to the width of the peaks of c m (T*), estimated based on temperatures for which c m (T*) E (c max m + c p m )/2, where c pm 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.

Fig. 10
Fig. 10 Representative configurations of the system with attractive walls (h* = À1) and L 1 = 120 for L 2 = 30 (left column) and L 2 = 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 m *, i.e. above the dashed line in Fig.7.
Fig.11The density profiles in the direction perpendicular to the walls, averaged over the longitudinal direction (see (9)), for several temperatures.Left panel: L 2 = 32, right panel: L 2 = 34.The color code is the same in both panels.The walls are attractive with h* = À1.

Fig. 14
Fig. 14 The orientational OP O L computed for L 1 = 180 and for the commensurate (L 2 = 30, L 2 = 34) and incommensurate (L 2 = 32, L 2 = 36) slit widths as a function of temperature.The thick solid line is the O L computed for the bulk system. 11 our procedure we considered typically n b = 151 different values of the reverse temperature, b, as: b k * = k Â db*, with k = 0, 1, 2,. ..,n b À 1, with db* = 0.10; whereas the chemical potential was fixed at m* = 6.
L 2 = 2, 3, 4,. ..,40.We took for all the cases L 1 = 120.For systems with L 2 = 4n + 2 (with n being an integer number) we considered n b = 401 values of the inverse temperature b, whereas for the other slit widths we took n b = 501.In both cases the particular values of b are given b k * = kdb*, with db* = 0.01 and k = 0, 1,. ..,n b À 1.This journal is © The Royal Society of Chemistry 2016