Structural behavior and dynamics of an anomalous fluid between solvophilic and solvophobic walls: templating, molding and superdiffusion

Confinement can modify the dynamics, the thermodynamics and the structural properties of liquid water, the prototypical anomalous liquid. By considering a general anomalous liquid, suitable for globular proteins, colloids or liquid metals, we study by molecular dynamics simulations the effect of a solvophilic structured and a solvophobic unstructured wall on the phases, the crystal nucleation and the dynamics of the fluid. We find that at low temperatures the large density of the solvophilic wall induces a high-density, high-energy structure in the first layer ("tempting"effect). In turn, the first layer induces a"molding"effect on the second layer determining a structure with reduced energy and density, closer to the average density of the system. This low-density, low-energy structure propagates further through the layers by templating effect and can involve all the existing layers at the lowest temperatures investigated. Therefore, although the high-density, high-energy structure does not self-reproduce further than the first layer, the structured wall can have a long-range effect thanks to a sequence of templating, molding and templating effects through the layers. We find dynamical slowing down of the solvent near the solvophilic wall but with largely heterogeneous dynamics near the wall due to superdiffusive liquid veins within a frozen matrix of solvent. Hence, the partial freezing of the first hydration layer does not correspond necessarily to an effective reduction of the channel section in terms of transport properties.

nucleation and the dynamics of the fluid. We find that at low temperatures the large density of the solvophilic wall induces a high-density, high-energy structure in the first layer ("templating" effect).
In turn, the first layer induces a "molding" effect on the second layer determining a structure with reduced energy and density, closer to the average density of the system. This low-density, lowenergy structure propagates further through the layers by templating effect and can involve all the existing layers at the lowest temperatures investigated. Therefore, although the high-density, highenergy structure does not self-reproduce further than the first layer, the structured wall can have a long-range effect thanks to a sequence of templating, molding and templating effects through the layers. We find dynamical slowing down of the solvent near the solvophilic wall but with largely heterogeneous dynamics near the wall due to superdiffusive liquid veins within a frozen matrix of solvent. Hence, the partial freezing of the first hydration layer does not correspond necessarily to an effective reduction of the channel section in terms of transport properties.

I. INTRODUCTION
The study of the properties of liquids confined at the nanometer scale is a topic of high interest for its technological, experimental and theoretical implications . Furthermore, confinement plays an important role in hydrated biological systems and organic solvents [22][23][24][25]. Structural, thermodynamical and dynamical properties of a liquid can change near an interface (solid, liquid, etc.) [26]. When the surface-to-volume ratio is large, at least along one direction as for the slit pore geometry, the effect of the confining surfaces has to be taken into account. Experiments and simulations on nanoconfined fluids show that molecules arrange in layers parallel to the surface. The effect becomes stronger for decreasing temperature or increasing density, until the fluid eventually solidifies. The nature of the solid, as an amorphous or a crystal, can depend on the interparticle potential and the confinement conditions [27][28][29][30]. In different cases the role of the interfaces can results in complex behaviors, e.g. a persisting fluid mono-layer around a spherical impurity while the rest of the system is in a polycrystal or glassy state [31], or persisting amorphous water mono-layer near an hidrophilic disordered surface [32]. In a recent experiment [33], Kaya and coworkers found that a thin-film water on a BaF 2 (111) surface remains in a high density liquid form for temperatures ranging from ambient (300K) to supercooled (259K). The result is unexpected because, based on thermodynamics arguments [34,35], Kaya and coworker would expect that the templating effect of BaF 2 (111) on the structure of the water film should promote the tetrahedral structure of the low-density liquid. Moreover, the presence of an interface can promote the heterogeneous nucleation of the crystal. However, recent experiments and simulations based on crystallographic analysis of liquids on crystals [36], or on colloidal self-assembly [37], show that the traditional theories need to be revised in these cases. Confined fluids, under suitable conditions of density and temperature, can spontaneously develop patterns, e.g. stripes, and different mesophases or more complex structures [38,39]. The confinement affects also the dynamics of the liquid. It has been found that near an interface there is a reduction of the local diffusivity of the liquid [40]. Computer simulations can help in interpreting the experimental results for nanoconfined fluids that are difficult to understand [33,[41][42][43][44]. Different numerical approaches based on first principles simulations can give detailed informations, but are limited by their high computational cost. Classical molecular dynamics of empirical fluid models employ parameters tested for the bulk case that not necessarily hold in confinement. The difficulty to adopt these bulk fluid models to the case of confinement leads to a variety of simulation results concerning the aggregation state of the fluid near the surfaces that in principle are model-dependent.
It is, therefore, useful to develop coarse-grained models that allow for analytic calculations [45][46][47][48][49][50][51][52][53][54] and more efficient simulations [55], like isotropic pairwise core-softened potentials, and that could allow us to better understand common features of fluids under confinement and the basic mechanisms of complex phenomena emerging in these systems, like pattern formation, e.g. stripes and different mesophases. This has been confirmed in a recent work [56] in which the authors claim that the origin of quasi-crystals could be understood in the context of a coarse-grained model by the competing effect of the hard and soft core radius of interacting particles.
Here we focus on the study of nanoconfined anomalous fluids, relevant for biological and technological applications [57,58], by means of molecular dynamics simulations of a system of identical particles interacting through the continuous shouldered well potential (CSW), an isotropic pairwise core-softened potential with a repulsive shoulder and an attractive well [59,60]. The CSW fluid is confined in a slit pore obtained by a solvophilic structured wall, and a solvophobic wall with no structure, as described in Sec.II B. The CSW model is suitable for studying globular proteins in solution [61], specific colloids [62][63][64], and liquid metals [65], and displays water-like anomalies [66]. In particular the CSW reproduces density, diffusion and structure anomalies following the water hierarchy [66] and displays a liquid-gas and a liquid-liquid (LL) phase transition, both ending in critical points [60].
Here we focus on structural and dynamical properties of the CSW model in confinement and show the role that the characteristic length scales of the inter-particle potential have in the structuring of the fluid.

A. The model
We consider a system of N identical particles interacting by means of the CSW potential confined between two parallel walls. The CSW potential is defined as [59] where a is the diameter of the particles, R A and R R are the distance of the attractive minimum and the repulsive radius, respectively, U A and U R are the energies of the attractive well and the repulsive shoulder, respectively, δ 2 A is the variance of the Gaussian centered in R A and ∆ is the parameter which controls the slope between the shoulder and the well at R R .
The parameters employed are the same as in Refs. [59,60,67]: In order to reduce the computational cost, we impose a cutoff for the potential at a distance r c /a = 3. In the present simulations we use ∆ = 15 that allows to better evidence the anomalies in density, diffusion and structure [59,60].

B. Theoretical and simulation details
In our simulations we consider a NV T ensemble system composed of N = 1024 particles at fixed temperature T and volume V . The temperature of the thermal bath T is kept constant by rescaling the velocity of the particles at each time step by a factor (T /T ) 1/2 , where T is the instantaneous kinetic temperature (Allen thermostat) [68]. Pressure, temperature, density and diffusion constant are all expressed in internal units: P * ≡ P a 3 /U A , T * ≡ k B T /U A , ρ * ≡ ρa 3 , and D * ≡ D(m/a 2 U A ) 1/2 , respectively. The equation of motion are integrated by means of the velocity Verlet method [68], using the time-step dt * = 0.0032 defined in units of (a 2 m/U A ) 1/2 (that corresponds to ∼ 1.7 · 10 −12 s for water-like molecules and to ∼ 2.1 · 10 −12 s for argon-like atoms [67]). We performed the same check as in Ref. [67] in order to verify that the value used for dt is small enough to satisfy the energy conservation of the system.
The confining parallel walls are placed along the z axis at a separation distance L z . The solvophilic wall is composed of a triangular lattice of CSW particles quenched with the position of the centers placed at z phil = 0. The lattice constant is d = a. The solvophobic wall has no structure and is obtained by imposing the repulsive potential U phob (z) ≡ (σ/z) 9 where z ≡ |z particle − z phob |, and z phob = L z [69]. In the following we consider the parameter σ = 1. We adopt periodic boundary conditions in the x and y directions. In order to compute the effective density ρ ef f , we need to compute the effective volume V ef f accessible to particles V ef f = L ef f z A, where A ≡ L x L y is the section of the simulation box and L ef f z is the effective distance between the plates. By considering the quenched particles forming the philic wall and the repulsive strength of the potential of the phobic wall, we obtain To explore different densities for the confined system in the NV T ensemble, we change L x and L y . We keep L z and N constant to exclude finite-size effects when we compare results for different densities. For each density, we equilibrate the system by annealing from T * = 4. For each temperature, the system is equilibrated during 10 6 time steps. We observe equilibrium after 10 4 time steps for the range of ρ and T considered here. All our results are averaged over 10 independent samples.
In order to check the stability of the system, we verify that the energy and the pressure are equilibrated. To compute the pressure we follow the approach of Refs. [69,70]. Due to the inhomogeneous nature of the system, the pressure is a tensor P(r) ≡ P K (r)+P U (r)+P W (r) where P K (r) ≡ k B T ρ(r)1 is the kinetic contribution as in an ideal gas,1 is the unit tensor, P U (r) is the potential contribution due to the interparticle interaction, and P W (r) is the contribution of the walls. At equilibrium the system is mechanically stable if ∇ · P = 0.
Considering the planar symmetry of the system, the normal (or orthogonal) and tangential (or parallel) component of the pressure can be written respectively as: P ⊥ (z) = P zz = const.
To verify the stability of the system, we computed the normal component of the pressure as a function of z The kinetic part is P K zz (z) = k B T ρ(z). For the potential part we used the Todd, Evans and Daivis formulation of the pressure [71] where F z ij is the z-component of the interaction force between particles i and j. The products of Heaviside step functions, Θ, select couple of particles that lie in different semispace respect to the plane parallel to the walls with coordinate z. The walls contribution to the pressure can be computed from Eq. 3, with z phil < z i < z phob for i = 1, ..., N, as [69] where F z i,phil ≡ j∈phil F z ij and F z phob,i ≡ −dU phob (z)/dz are the interaction forces of the particle i with the philic wall and the phobic wall, respectively. The different components compensate each other in order to keep the normal pressure constant (Fig.1). In particular, the philic wall contributes with a positive term P W phil from z = 0 to z = R R /a = 1.6 and with a negative term from z = R R /a = 1.6 to z = r c /a = 3, due to the interaction force between the particles of the fluid and those of the philic wall. The phobic wall contributes with a positive term P W phob from z ≃ L z − (1/T ) 1/9 to z = L z due to the repulsive potential. The kinetic pressure P K (z) is a signature of the density profile, while the potential pressure P U (z) can contribute with a positive or negative term depending which part of the interparticle interaction, repulsive or attractive, respectively, dominate. After having verified that the system is stable, i.e., P ⊥ = P zz = const., we compute the normal pressure where z * 1 and z * 2 are z-coordinates sufficiently near the philic wall and the phobic wall, respectively, for which the kinetic P K zz and the potential P U zz part of the pressure are zero. In particular, P W phil is the contribution of the philic wall, P W phob is the contribution of the phobic wall, P K is the kinetic contribution as in an ideal gas, P U is the potential contribution due to the interparticle interaction.

III. RESULTS AND DISCUSSION
A. Density profile and aggregation state As pointed out in the introduction (Sec.I), the confinement can modify the structure of a fluid resulting in an inhomogeneous density profile. In a slit pore geometry, near the confining walls, particles form layers parallel to the walls as the temperature is decreased or the density is increased, as shown by our calculations for the density profile ρ(z) (Fig. 2).
To establish the aggregation state for each layer, we compute, layer by layer, the lateral radial distribution function g (r | ); the 2d Voronoi tessellation of each layer; the mean square displacement (MSD) and the survival probability (SP) of molecules in each layer. All these quantities together, as we discusse in the following, allow us to identify the presence and coesistence of the solid, heterogeneous fluid and homogeneous fluid.
Our analysis shows that the system organizes forming layers parallel to the solvophilic wall at any temperature T and average density ρ. At high T and low ρ we find only one well defined layer and no layering near the solvophobic wall, while the whole system is in the fluid state. By decreasing T and increasing ρ, the number of layers increases up to eleven. Furthermore, at higher value of ρ, the layers appear also near the solvophobic wall.
However, for high enough T and ρ we observe that away from the walls the system is in the fluid state. Nonetheless, the walls affect the density profile, changing system density and aggregation state, over an extension that in our case can be up to 11 layers that depends on temperature and average density of the system.

B. Radial distribution function and spatial configuration analysis
The in-layer radial distribution function g n (r ) for the n-th layer (with n = 1, 2, ..., 11) is computed as where ρ n and z n are the density of particles and the z-coordinate of the layer n, respectively, r || = (x 2 + y 2 ) 1/2 is the transverse distance between two particles in the same layer (and in internal units is r * = r /a). The Heaviside step functions, Θ, select couple of particles that lie in the layer n of width δz. The g n (r ) is proportional to the probability of finding a molecule in the layer n at a distance r from a randomly chosen molecule of the same layer n. The definition of the layer n in which lies a particle i is once for all established according to the value of the z-coordinate of the particle i as: for the first layer 0 < z n=1 i < (3/2)δ z , with δ z = R R /a, and (j − 1/2)δ z < z n=j>1 i < (j + 1/2)δ z for the others. This is a natural choice because particles at low temperatures or high densities tend to stratify in layers whose interdistance is approximately equal to δ z .
At low density (ρ * = 0.11, Fig.3), we observe that the system is in a fluid state for high temperature (T * = 1.4) in any layer. The g n (r ) of the layer n = 1 shows a first peak around the shoulder radius (r * ≃ 1.6) and a second peak around the attractive well radius (r * ≃ 2), while for the other layers only the second peak is present. We interpret this difference between the first and the other layers as the consequence of a "templating" effect of the solvophilic wall that at high T * is observed only on the first layer.
At T * = 0.7 the system shows the same behavior observed for higher temperatures, except that the layer near the solvophilic wall develops patterns. This behavior is reminiscent of what has been observed in monolayers with an interparticle potential composed by a hard core and a soft repulsive shoulder [38].
At T * = 0.6 the layer near the solvophilic wall is still showing patterns, while the layer n = 2 is forming crystal patches. This is evident from the analysis of the g n=2 that goes to zero for r * ≃ 2.75 and 4.75 at T * = 0.6, consistent with an incipient triangular crystal with lattice step given by the interaction potential attractive distance R * A = 2. The other layers are in a fluid state.
At T * = 0.5 the first layer is forming a hexagonal crystal. Although the hexagonal crystal in n = 1 does not overlap exactly with the triangular wall structure, the comparison of the g n=1 and g n=0 of the wall shows a strong correlation between the two structures, suggesting a "templating" effect. This effect due to the attraction to the solvophilic wall is so strong at T * = 0.5 that is forcing particles to be at their repulsive distance. The high-energy cost of the resulting honeycomb lattice forming in the layer n = 1 is compensated by the large number of attractive interactions between the particles at n = 1 and those of the wall (n = 0) from one hand, and between the particles themeselves at n = 1 from the other hand. This free energy minimization process is analyzed in Sec.III F in the discreet potential approximation to understand the stripe phase formation. The triangular structure that was incipient for n = 2 at high T * , for T * = 0.5 is well defined for n = 2 and n = 3, with defects in the layer n = 3. This triangular structure is the dual lattice of the n = 1 hexagonal layer and its formation is the consequence of a "molding" effect of the layer n = 1 onto the layer n = 2. Note that while the wall (n = 0) layer has a templating effect on the n = 1 layer, the n = 1 layer has a molding effect on the n = 2 layer. The difference between the two cases is due to the smaller density of the n = 1 layer with respect to that of the wall. The smaller density does not allow to compensate the high energy cost of the propagation of the hexagonal crystal to the layer n = 2. On the other hand, the triangular crystal of the n = 2 layer is energetically favorable, because the particles are all at the attractive distance, and at this temperature can propagate to the n = 3 layer again with a "templating" effect. The layer n = 4 is made of a few triangular crystallites immersed in the fluid, while the other layers are in a fluid state.
At T * = 0.3 both the templating and the molding effect are stronger. In particular the template of the n = 2 layer propagates over all the six layers that are formed at this density and temperature. Fig.4), for T * = 1.4 and T * = 0.7 we observe the same qualitative behavior as for the low density case. For T * = 0.6 the layer n = 1 has less tendency to form patterns respect to the low density case, and the layer n = 2 to order in a crystal structure. Therefore, the confined system is more fluid at this density than at lower density. We understand this result as a consequence of the larger hydration at higer density.
At T * = 0.5 the first layer has partially crystallized in the hexagonal and partially in the triangular structure following the template of the wall. Therefore, the templating effect is now stronger then the corresponding case at lower density. The hexagonal crystal shows now a preferred direction of symmetry. This direction propagates to the layer n = 2, where we observe stripes along the preferred direction. The stripes propagate up to n = 4 layer, while the other layers are in a fluid state. The peak of g n (r ) at r * ≃ 2.1 (that corresponds to the average second nearest neighbor distance) is a signature of the stripe phase formation.
At T * = 0.3 the preferred direction in the deformation of the hexagonal crystal for n = 1 is more evident and we observe a clear stripe phase for the layers from n = 2 to n = 4, with a peak of g n (r ) at r * ≃ 2.1 more pronounced than the case at T * = 0.5. The other layers form a triangular crystal at the attractive distance.
At high density (ρ * = 0.30, Fig.5), for T * = 1.4 we observe the same qualitative behavior as for the lower density cases. At T * = 0.7 we found that the only difference with the lower density case is that the layer n = 1 is forming crystallites following the template of the wall.
At T * = 0.6 the layer n = 1 has a different and incipient crystal structure (Kagome lattice) with defects that is better defined at lower T . This is evident from the analysis of the g n=1 that goes to zero for r * ≃ 1.6 and 2.75 at T * = 0.6. The layer n = 2 shows patterns very close to the stripe configuration. The corresponding g n=2 goes to zero for r * ≃ 1.5 and shows a peak for r * ≃ 2.1. These characteristics of the g n are typical of a stripe phase. The layer n = 3 and n = 4 still show patterns close to the stripe phase, but in a less pronounced way.
From the layer n = 5 to the n = 10 the pattern is vanishing. The layer n = 11 is showing an incipient triangular crystal with lattice step given by the interaction potential attractive distance R * A = 2. This is evident from the analysis of the g n=11 that approaches zero for r * ≃ 2.75 and 4.75 at T * = 0.6.
At T * = 0.5 the layer n = 1 is forming a kagome crystal with defects. The layers from n = 2 to n = 10 show a stripe phase and the layer n = 11 is forming a triangular crystal with defects.
At T * = 0.3 the Kagome crystal of layer n = 1 has no defects. The layers from n = 2 to n = 10 are in a stripe phase and the layer n = 11 is forming a well defined triangular crystal.
As discussed above, the layer n = 1 close to the solvophilic wall (n = 0) is subjects to the templating effect for all densities at low temperatures. In Fig.6 we compare the g n=1 (r ) and the snapshots of the first layer for T * = 0.3 at several densities. In order to compare layers, that correspond to different densities, between them, we considered a portion of each layer of the same size (L x xL y of the system at ρ * = 0.30). For densities ρ * = 0.11, 0.13, 0.16 the first layer is forming a distorted hexagonal lattice characterized by a g n=1 (r ) with a first peak at r * /a ∼ 1.15 and vanishing for r * ≃ 1.6 and 2.75. By increasing ρ we observe a progressive shift to higher values of r of all the peaks of g n=1 (r ), but the first that, instead, is becoming more pronounced as a consequence of a better local order. For densities between ρ * = 0.20 and ρ * = 0.25 the first layer shows a polycrystal phase with coexistence of triangular and square lattices. This corresponds to an intermediate stage toward the well defined Kagome lattice that is formed for ρ * ≥ 0.27, as showed by the splitting of the second peak of g n=1 (r ) into two close peaks at r * ≃ 1.9 and 2.2.

C. Mean square displacement analysis
In order to characterize space-dependent diffusion properties of our system, we compute the mean square displacement (MSD) associated to each layer of the slit. We observe that, except for low temperatures, a particle can visit different layers in which the aggregation state can change from homogeneous to heterogeneous liquid and vice versa. For this reason we calculate the MSD only for those particles that remain in a layer over the entire time interval under consideration and we average over all possible time interval. Therefore, the MSD associated to each layer n is defined as where τ ≃ t−t 0 is the time spent in the layer n by a particle that entered in the layer at time t 0 . In according to the standard definition of the MSD, lim τ →∞ (∆r n (τ )) 2 = 4D τ α where D is the lateral, or parallel, diffusion coefficient and α the diffusion exponent. The value α = 0 means that the system is arrested, as in a solid state where particles can only vibrate around theirs equilibrium positions; for 0 < α < 1 the system is subdiffusive corresponding in general to particles diffusing in complex structures (with non trivial microscopic disorder); α = 1 is the standard diffusive behavior as in a normal fluid state. For α > 1 the system is superdiffusive. On the other hand, for early times free diffusion we expect the ballistic regime with α = 2. Our analysis (Fig.7) shows that the in-layer MSD has always a ballistic regime for t * ≤ 10. The corresponding mean displacement is approximately half particle diameter at high T and low ρ and weakly decreases for increasing ρ and decreasing T corresponding to the expected decrease of the mean free path of the particles.
For T * ≥ 1.4 all the layers reach the diffusive (α = 1) behavior for long times. By decreasing the temperature the behavior of the layers becomes more heterogeneous. In particular, we observe that the layer n = 1 near to the solvophilic wall slows down in a sensible way with respect to the layers at T * ≤ 0.7 and becomes arrested for T * ≤ 0. instead, at T * = 0.5 and ρ * = 0.30 we observe that for all the layers but the one near the solvophobic wall (n = 11), after the plateau, the dynamics enters in a superdiffusive regime with 1 < α < 2. This effect is related to the presence of defects and of a nonuniform stress field, as discussed in Sec.III F. For T * = 0.3 we observe that all the layers are arrested at ρ * = 0.11. At this low density the slit is only partially filled (Fig.7). At T * = 0.3 and ρ * = 0.18 and ρ * = 0.22 also the layer n = 11 near the solvophobic wall is present and it is characterized by a larger MSD with respect to the other layers and by a diffusive regime at long times. The other layers have an arrested dynamics. At T * = 0.3 and ρ * = 0.30 all the layers from n = 1 to n = 10 have a superdiffusive regime at long times, while the layer n = 11 reaches the diffusive log-time regime. However, its MSD is smaller than that of the other layer for very long times (t * ≥ 100).

D. Survival probability function analysis
As the time proceeds, the average in Eq. 6 for the MSD is performed on a decreasing number of particles because some of them can leave the layer. A priori this reduction of the statistics is not homogeneous, that means that in general there can be a correlation between particles that leave the layer and theirs properties, as theirs velocity components. Therefore, for T * ≥ 0.5, low enough ρ, and for the most diffusive layers there is a time, τ max , after which the in-layer MSD (Fig.7) is not well defined. To estimate τ max as function of T and ρ for different layers we analyse the population relaxation of particles in each layer. In particular, we compute the survival probability (SP) function, S i (τ ), which is the probability that a given particle stay in the layer i for a time interval τ . The SP can be calculated as where N i (t) is the number of particle in the layer i at time t and N i (t, t + τ ) is the number of particle that do not leave the layer i during the time interval [t, t + τ ]. The SP give an indication of the time interval τ max over which the MSD is well defined. We observe that S(τ ) has an exponential decay in our simulations (Fig.8). We, therefore, define τ max as the characteristic decay time S(τ ) ∼ e −τ /τmax . This choice is consistent with the observation that the MSD in Fig.7 is well defined when S(τ ) ≥ 1/e. We observe that for T * ≥ 0.6 and all the densities and for T * = 0.5 and ρ * ≤ 0.22, the SP decay is slower for the layer n = 1 near the solvophilic wall and becomes faster for the layers away from the two walls. When the layer n = 11 near the solvophobic wall is present, we observe that it has a decay in SP slower than those layers that are farther away from the wall. At T * = 0.3 for all densities, and at T * = 0.5 for ρ * = 0.30, there is no decay in SP, consistent with the crystallization of the layers. The non monotonic behavior of τ max is reported in Fig.9 as a function of layers for different densities and temperatures. By comparing Fig.9 and Fig.2 we observe that τ max increases when the layers are more structured in the z direction. Hence, both walls facilitate the stratification of the fluid, although the structureless phobic wall does it in a less strong way with respect to the structured solvophilic wall.

E. Liquid veins
In Sec.III C, analysing the MSD layer by layer, we have seen (Fig.7) that for high densities and low temperatures, after a plateau, the dynamics can enter in a superdiffusive regime with 1 < α < 2. In this section we show how this behavior is due to the formation of liquid "veins" in such layers. nanometer-sized ice crystal [72]. In polycrystalline systems, the liquid is found along intergranular junctions, as grain boundaries (see [72] and references therein for the case of water). Residual stress in these polycrystal structures can be localized along integranular junctions, and can results in an effective force that acts on fluid particles present in these junctions. The origin of the residual stress in our system is due to the fact that when the fluid solidifies as the temperature is decreased, the minimization process of the free energy take place locally, instead of globally. In glass forming liquids, this effect is caused by a fast cooling, while in our system it is due to the layering of the fluid caused by the confinement.
In Fig.10 we show the spatial configuration of the first three layers of the system close to the philic wall, for three different runs (i.e. for three different realization of initial conditions) at T * = 0.3 and ρ * = 0.30. We observe that particles in the first layer (n = 1) are characterized by the same MSD, while particles in other layers (n = 2, 3) can have different MSD. In particular in some configurations, we observe veins with mobility higher than the rest of the system (Fig.12). We analyzed the trajectories of the particles of these specific realizations of the system (Fig.12). For these cases we find that these particles with a MSD higher than the majority belong to the same stripe and diffuse along the stripe itself. We observe that the majority of particles in the layer n = 2 (and in a less evident way for the layer n = 3 (Fig.12b), remain spatially localized during the entire simulation, except those belonging to two stripes moving in the same direction as along stripe veins (Fig.12a). We observe a similar situation for the layer n = 3, but here all the particles are more mobile and the particles in the veins move in opposite directions. Further analysis, that goes beyond the goals of the present work, is necessary to understand the effect of the vicinity of the solvophilic wall and if the veins are related to point-like defects as seems to be suggested by  FIG. 11. Single-particle MSD for layers and runs that are in Fig.10. For Run #1, red and green colors are used for those particles belonging to veins performing a dynamics different from the rest of the particles in the layer.
FIG. 12. In-layer trajectories for particles in the second (a) and third (b) layer of Run #1 in Fig.10. In red and green we show the trajectories of representative particles belonging to stripes veins. The MSD of these particles is represented with the same color code in Fig.11. We apply periodic boundary conditions for y * = 0, L y where L y ≃ 12.1 (dashed lines) and for x * = 0, L x where L x = 15.

F. Voronoi tessellation and structural analysis of the solid state
Our MSD and SP analysis show that at low T and high ρ there are layers taht behave as a solid. However, by looking only at the MSD and SP is not possible to establish if a solid layer is in an amorphous, crystal or polycrystal state [73]. In order to better understand the structure of solid layers, we computed the standard 2d Voronoi tessellation (useful to identify defects present in the crystal structures, as vacancies, Frenkel-like, dislocations and grain boundaries), and a modified version of it (suitable to identify distorted crystal structures). With this analysis we can also disentagle the role that the three relevant length scales (the diameter of the particles a, the repulsive radius R R , and the attractive minimum R A ), giving rise to two competing length scale R R /a and R A /a, play in the determination of layer's structure. Indeed, the interdistance between two adjacent layers is ∼ R R , while when stripes form in a specific layer for intermediate densities, just to consider a specific case, particles within a stripe are compressed at a distance ∼ a, while the distance between stripes depends on ∼ R R and ∼ R A , as discussed in the last part of this section.
In the standard Voronoi tessellation we construct polygons centered around particles forming a lattice whose edges are crossed in their middle point by the edges of the Voronoi cells. This procedure garantees that each Voronoi cell represent the proper volume of each particle. To better visualize the result, we represent Voronoi cells having a different number of edges with different colors. To reduce the noise in our analysis we adopt also a modified version of the Voronoi tessellation in which we associate a color to a polygon in according to the number of edges of the polygon that have a length 10% larger than that of the average edge lenght calculated over the specific polygon itself. This procedure allows us to better visualize polycrystal structures despite the presence of small lattice deformations.
We compute the Voronoi tessellation for low density (ρ * = 0.11) and high density

3/2
A /(k B am 1/2 ). At low density (ρ * = 0.11), the first layer at T * = 0.3 (Fig.13) is in a frustrated solid state that by annealing toward T * = 0.0005 (Fig.14) becomes a frustrated polycrystal. The very low-T polycrystal has point and line defects as grain boundaries dividing a deformed honeycomb lattice (the deformed green triangles) from a stripe phase (the stretched hexagonal cyan polygons). For both considered T * the other layers are organized in a triangular lattice (where each particle is surrounded by a hexagonal cyan polygon). We only observe defects, such as dislocations (Run # 2 in Fig.13a and Fig.16a) that are not eliminated by annealing (Run # 2 in Fig.13b).
At high density (ρ * = 0.30), the first layer at T * = 0.3 (Fig.14a) is in a polycrystal state with defects. At T * = 0.0005 (Fig.14b) we observe two principal crystal grains: a triangular lattice (cyan polygons) and a Kagome lattice with defects (blue rhombouses). At T * = 0.3, the layers n = 2, 3 present a zigzagging stripe structure with orientation and angles that can change from run to run. At T * = 0.0005 the stripe structure of these layers becomes more regular. These observations emphasize that the increase of density induces an increase of disorder in the solid layers, propagating from the layer n = 1 to the other layers and up to the layer n = 11. The formation of crystal defects during the annealing, and the fact that they are different for different initial conditions, indicate that the system cannot reach easely the global minimum of the free energy landscape, corresponding to the crystal configuration, but is trapped in local minima due to the slowing down of the dynamics and the templating effect of the solvophilic wall. In particular, the mismatch of the wall structure with the bulk crystal structure induces a frustrating effect that is more evident near the wall (in layers n = 1 and n = 2) for increasing density. To understand the formation of stripes, we follow the same approach as in Refs [74,75] to show that, if the principal contribution to the minimization of the free energy comes from the energetic therm, in the discreet potential approximation, under suitable conditions of density and temperature, particles organize in straight or zigzagging stripes. In the discreet potential approximation the energetic terms can be reduced to the soft core (U R ) and the attractive well (U A ). As a consequence the energetic contribution coming from the interaction with particles of the adjacent layers is approximately constant when the layer density is fixed. Therefore, the energetic cost of stripes formation is determined only by the contribution of the in-layer particle interactions. In particular, if a layer has a triangular structure, as the stable configuration of layers n =2, ...,11 at low density (Fig.13a,b), then for sufficiently high density the layer will prefer to form stripes.
FIG. 15. 2d schematic representation of particles composed by an hard core (in dark blue) of size a, a soft corona (in cyan) of size R R and an attractive external corona (in green) that extends uo to the potential cutoff r C . The equilateral and isosceles triangles represent the unitary cell of the triangular and straight-stripes lattice, respectively. An example of (maximal) zigzagging-stripes of the same density of straight-stripes is also shown.
Consider our fluid made of particles with a hard core surrounded by a soft corona and an external attractive corona (Fig.15). If d is the lattice constant of the triangular structure at the soft-corona distance, then the density of the layer is where N l is the number of particles present in this layer. If we allow the triangular lattice to deform in order to minimize the energy of the layer, the new unit cell will be composed by the isosceles triangle in which one side is equal to x and the other two are equals to y with x < d < y. The fact that the density doesn't change implies that y = (3d 4 )/(4x 2 ) + x 2 /4 = 1/(ρ 2 x 2 ) + x 2 /4 = y(ρ l , x).
The energy per particle of the layer is For N l ≫ 1 it becomes U l ≃ U x + 2U y , where U x and U y are the energy associated to the interaction between the particle along the x and y side of the triangle respectively.
In the discreet potential approximation it is U x = U R for a < x < R R , U x = −U A = −U R /2 for R R < x < r C and U x = 0 for x > r C (note that we obtain the same result if instead of r C we consider any value between R R and r C . Indeed, the present approach has been applyed to show the stability of stripes cnofiguration for a pure repulsive potential model [74,75]). The same holds for U y substituting x with y.
For sufficiently high density, i.e. for d < R R or ρ l > 2/( √ 3R 2 R ), the energies per particle associated to a layer formed by equilateral or isosceles triangles are U equi l = 3U R or U iso l = 0, respectively. Therefore, under these conditions, the layer will prefer to form stripes.
Furthermore, from geometric consideration it is possible to conclude that the zigzaggingstripe lattice can be obtained as a deformation of the straight-stripe lattice without changing the density and keeping the energies per particle U iso l ≃ 0 [74,75]. In general, many different zigzagging stripe lattices are possible all with comparable energy per particle (Fig.15).
Considering the stripes that form in the layers n = 2, ..., 9 for system density ρ * = 0.30 and temperature T * = 0.3, the resulting average in-layer density is ρ l ≃ 0.53, and x ≃ 1.05.
Hence, from the previous equation for y = y(ρ l , x), we find y ≃ 1.87. In view of all the approximation made, we consider this value consistent with y ≃ 2 of the distance between the closest particles belonging to two adjacent stripes in the same layer.

IV. CONCLUSIONS
By considering many layers of a confined anomalous fluid [21,[76][77][78] we show that the effect of the structured solvophilic wall can extend up to the entire slit pore. In particular, we study structural and dynamical properties of a monocomponent anomalous liquid under confinement. The fluid has two characteristic distances and can be considered as a coarsegrained model for globular proteins [61], colloidal systems [62][63][64] or, to some extent, liquid metals with water-like anomalies [65]. We perform molecular dynamics simulations of the fluid in a slit pore with a solvophilic wall and a solvophobic wall. The solvophilic wall has structure while the solvophobic one has no structure. We observe that the molecules organize in an inhomogeneous way, forming layers that are parallel to the surfaces, with higher density near the solvophilic surface with respect to the center of the slit pore. For sufficiently high densities, for which the fluid occupy entirely the pore, we observe an increase of density also close to the solvophobic surface, but in a less prominent way. These results are consistent with experimental and theoretical works for nanoconfined fluids.
At low temperature we observe coexistence between the homogeneous liquid, heteroge-neous liquid and solid phase of the fluid. The influence of the structured solvophilic surface on the solid layers can extend as far as the sixth hydration layer at low T and high ρ. In particular, we find a strong correlation between the structure of the solvophilic surface and that of the first layer suggesting a "templating" effect. Indeed, the large density of solvophilic surface particles allows the formation of a first layer at high density. The high energy cost of this first layer is compensated by the large number of attractive interactions between the particles of the first layer and those of the surface. Further energy gain comes from an extra energy term due to the first in-layer particle interactions.
Moving further from the wall, we find that the first layer has a "molding" effect on the second layer. This is because the density of the first layer is smaller than that of the surface and is not high-enough to propagate its template. Nevertheless, the low-density second layer is in condition to template the third layer replicating its structure and inducing a long-range effect that can, eventually, involve the whole system.
From the calculation of the mean square displacement and the Voronoi tessellation we conclude that at low temperature the first layer close to the solvophilic surface is a polycrystal with two competing phases that generate low-energy states with high degeneracy and very slow dynamics. At low densities the two competing phases are stripes and honeycomb lattice, while at high densities are triangular and kagome lattice. We understand this result as a consequence of the high density (triangular) structure of the solvophilic wall with a lattice step that corresponds to the hard repulsive distance of the solvent, and the strong wallsolvent attractive interaction. These properties of the wall generate in the first solvent layer local regions with density and energy that are higher than the average of the layer. As a consequence, other regions within the layer have density and energy below the average, giving rise to a competing crystal structure. Our results remind us of some recent experiments and simulations for a thin-film of water on BaF 2 (111) surface for which the authors found a very high density first interfacial layer for all temperatures, while they would expect, from thermodynamic arguments, a lower density liquid at supercooled conditions [33]. In other recent experiments and simulations [36,37], the authors pointed out that it is necessary to revise the theory of heterogeneous nucleation when the crystalization induces a non-zero entropy at zero temperature and the system initially is far from equilibrium.
Apart from the layers close to the two surfaces, we observe that the structure of each layer mainly depends on its density. In the case of the stripe phase, using simple geometrical and energetic considerations, we find that straight and zigzagging stripes are the stable configurations for intermediate densities.
Furtheremore, analysing the mean square displacement layer by layer, we observe layers at high densities and low temperatures with a caging-like behavior characterized by a ballistic dynamics followed by an arrested state (plateau) and a superdiffusive regime with a diffusion exponent 1 < α < 2. Our analysis shows that this behavior is due to the formation of liquid veins within the stripe phase. In particular we observe that each vein can behave differently from the others diffusing in one of the two possible directions along the stripes and having a different diffusion exponent 1 < α < 2. We rationalize the different possible values of α as a consequence of the presence of residual stress that could introduce an effective force acting on the fluid. Under suitable conditions, e.g. a constant effective force along the stripe, the particles in the vein could perform a biased one-dimensional random walk characterized by an exponent of the MSD that approaches the ballistic value (α = 2). The behavior of these veins can be analyzed in a more quantitative way by computing, for example, the temporal autocorrelation function, the intermediate scattering function [79,80], the relative displacement of nearest neighbors or the particles displacements following the Lindemann criterion [81]. We will present this analysis in future works.
Our results show that the dynamical slowing down of the anomalous solvent near the solvophilic wall does not imply by necessity the complete freezing of the first hydration layers, because at low T and high ρ we observe largely heterogeneous dynamics in three layers with the formation of liquid veins within a frozen matrix of solvent. Therefore, under these considerations the partial freezing of the first hydration layer does not correspond necessarily to an effective reduction of the channel section in terms of transport properties, at variance with the conclusions of Ref. [82].