Isomorphs in nanoconfined liquids

Throughout the study, we use two different sets of nondimensionalized units: One set is based on the microscopic parameters of the LJ potential with length scale $\sigma_{\rm AA}$ and energy scale $\epsilon_{\rm AA}$ of the larger (A) particle, which is standard in computer simulations, and another set of nondimensionalized units using macroscopic quantities with length given in units $\rho^{-1/3}$, energy in units of $k_{\rm B}T$, and time in units of $\rho^{-1/3}\sqrt{m/k_{\rm B}T}$ ($m$ is the particle mass) as applied in isomorph scaling paper4 . We refer to macroscopic nondimensionalized units as reduced units and use a tilde above the variable name to indicate a reduced quantity; otherwise LJ units are implicitly assumed.

Nanoconfinement is modelled using a slit-pore geometry in which the 100 surface of an fcc crystal is exposed to the liquid. The two crystal planes are in registry (i.e., out of sync). The distance between the two walls is denoted by $H$, measured from the centers of the confining fcc particles (see Fig. 1).

The liquid-wall pair interactions are also described by the LJ potential in Eq. (1). For the SCLJ liquid we use the parameters $\sigma_{\rm AW}$ = 1 and $\epsilon_{\rm AW}$ = 1, in which $W$ denotes a wall particle. For the KABLJ mixture we derived the interaction parameters from Lorentz-Berthelot mixing rules with $\sigma_{\rm AW}$ = 0.97, $\epsilon_{\rm AW}=1$, $\sigma_{\rm BW}$ = 0.91, and $\epsilon_{\rm BW}$ = 0.707, where we used $\sigma_{\rm WW}$ = 0.94 and $\epsilon_{\rm WW}$ = 1 to calculate these numbers. The cutoff of the liquid-wall pair interaction is $r_{c}$ = 2.50$\sigma_{\rm\alpha W}$. The density of the fcc walls is for the majority of the simulations kept fixed at $\rho_{\rm W}$ = 1, but we study also the effect of varying this parameter. The walls consist of around $N$ = 1600 particles.

For the SCLJ liquid we use the density $\rho$ = 0.85 as a reference with $T$ = [0.70, 10] and $H$ = [2, 10]. For KABLJ we focus on $\rho=1.20$ with $T$ = [0.7, 12] and $H$ = [4, 10]. These densities are standard for studies of bulk liquids. The simulations use a Nosé-Hoover NVT thermostat on the liquid particles whereas, for simplicity, the wall particles are frozen in place, i.e., $T$ = 0 for the walls. Thermostatting the walls and the nature of the exposed surface are known to have observable effects on the structure and dynamics of confined liquids tox1981 ; schoen1987 ; schoen1988 . The phase behavior of this particular type of confinement has been studied in detail elsewhere see, e.g., Ref. dominguez2003 .

In sections IV and V, we compare results of the topological cluster classification malins2013tcc for the confined system with the bulk system at the same state points (density and temperature). For these simulations we use standard Monte-Carlo simulation with $N=4000$ particles for the models specified above.

This section studies how the above mentioned quantities are affected by changing the following parameters related to the confinement: the distance between the two walls $H$, the strength of the liquid-wall interaction $\epsilon_{\rm LW}$, and the surface roughness $\rho_{\rm W}$.

Figure 2 shows how $R$ and $\gamma$ (Eqs. (2) and (4)) vary with temperature for several slit-pore widths in the range $H$ = [2, 10]. The confinement thus ranges from almost a single layer of liquid particles to more bulk-like conditions. The densities of the liquid and wall are kept constant with $\rho_{\rm L}=0.85$ and $\rho_{\rm W}=1$, respectively. For simplicity, we here and henceforth define the liquid density as $\rho_{\rm L}\equiv N/AH$, where $A$ is the exposed surface area of the crystal. We thus do not take into account any excluded volume near the walls when calculating the confined liquid density truskdumbbellbulk . As a reference, the bulk liquid has $R\approx 0.96$ at the chosen liquid density.

The wall separation $H$ = 2 shows markedly different behavior from the other slit-pore widths having $R<0.90$ for all temperatures. This is consistent with an earlier study that also observed breakdown of R-simple liquids for wall separations close to a single particle layer ingebrigtsen2013 . However, as the slit-pore width is increased $R$ also increases and already at $H$ = 4, which still corresponds to a strongly confined liquid, we find good correlation between the virial and potential energy, i.e., $R>0.90$. Depending on $H$, we find an increase in $R$ with temperature as close encounters start to play a bigger role. We find also some irregularities in these observations for low temperatures due to crystallisation in the layer closest to the wall.

The density-scaling exponent $\gamma$ displays the opposite trend in Fig. 2(b) with a monotonic decrease with temperature for all $H$. $\gamma$ is noted to increase with slit-pore width, but no theory currently exists for how $\gamma$ should depend on $H$, as is the case for bulk liquids paper2 ; thermoscl . We find a possible slit-pore-width dependent plateau for $\gamma$ between 4 and 5, which signifies that $\gamma$ is also dependent on $H$ in confined liquids. More investigations are nevertheless needed to determine if this $H$-dependent plateau truly exists and if it has a physical significance.

We now consider the effect of varying the attraction between the wall and the confined liquid particles. Figure 3 shows $R$ and $\gamma$ as a function of $\epsilon_{\rm LW}$ at $T$ = 2 and $\rho$ = 0.85, again for several $H$. As the simulations use $\epsilon_{\rm AA}=1$, for $\epsilon_{\rm LW}>1$ we have an “attractive” wall and when $\epsilon_{\rm LW}<1$ it becomes a “repulsive”’ wall with respect to the interactions between the liquid particles.

Figure 3(a) shows that $R$ depends significantly on $\epsilon_{\rm LW}$. Depending on $H$, the correlation coefficient $R$ may both decrease and increase when the wall becomes more attractive or repulsive. The interplay between the liquid-liquid and liquid-wall interactions is thus highly nontrivial. The density-scaling exponent in Fig. 3(b) displays a behavior that mimics that of the correlation coefficient, but with the maximum displaced to lower values of $\epsilon_{\rm LW}$.

The highest value of $R$, and therefore the maximum, is expected to appear when the wall particles are most similar to the liquid particles, which means here around $\epsilon_{LW}=1$. In spite of this, we observe that the maximum occurs somewhat to the right around $\epsilon_{LW}=1.5$ for $R$ and the opposite for $\gamma$. We currently have no explanation for why this is the case.

Another means to probe the coupling between the confined liquid and the walls is to change the density of the fcc walls $\rho_{\rm W}$. In Fig. 4 we show the correlation coefficient and density-scaling exponent as a function of the ”surface roughness” of the crystal (i.e., as a function of $\rho_{\rm W}$). The liquid density is kept constant at $\rho_{\rm L}$ = 0.85 and $H$ = 6.

We observe for low $\rho_{\rm W}$, i.e., when the surface roughness is high, that $R$ decreases when the density of the wall particles is reduced, though it remains above 0.90. This effect can be attributed to particles penetrating into the walls (not shown). For high $\rho_{\rm W}$ the correlation cofficient remains virtually constant. Almost no effect on $\gamma$ is noted in Fig. 4(b) for both low and high $\rho_{\rm W}$.

We now turn our attention to isomorphs in the nanoconfined system. To this end, we consider the behaviour along an isomorph and contrast it with an isochore at a density $\rho_{\rm L}=0.85$ as the isomorph concept is approximate. Isomorphs in the SCLJ liquid were generated by the direct-isomorph-check method (see Sec. III). We consider two different slit-pore widths, one with $H\approx$ 6 and one with $H\approx$ 10; recall that $H$ is adjusted with the liquid density along isomorphs. These two distances span a strong and medium confined liquid.

Density profiles perpendicular to the walls along an isomorph with $H\approx$ 10 are shown in Fig. 5(a). For comparison, results for an isochore with the same temperature variation are given in Fig. 5(b). From this point forward a yellow figure background denotes data obtained along an isomorph and a pink figure background denotes data obtained along an isochore.

We find that the density profiles to a good approximation are invariant along the isomorph with a 24% density increase, whereas this is seen not to be the case along the isochore in Fig. 5(b). For the isochore, significant changes in all peak heights with temperature are observed, in particular for the layer closest to the wall indicating pre-crystallisation even at $H$ = 10. This pre-crystallization is, however, nicely preserved on the isomorph.

Figure 6 displays the in-plane density profiles in the layer closest to the walls, i.e., for the layer around $|\widetilde{z}|\approx$ 4, for the first and the last state points of Fig. 5. The density profile is shown for one unit cell of the crystalline walls. The isomorph shows excellent scaling with the liquid particles situated in between the wall particles and exhibits very little change in the density profile whereas the isochore displays a density field that is increasingly smeared out as $T$ is increased, confirming the pre-crystallisation. Similar behavior is observed the remaining layers (not shown).

Dynamical properties such as the mean-squared displacement (MSD) or the intermediate scattering function are also invariant along an isomorph and we now examine how well this behaviour holds in nanoconfinement. To do so, we consider the reduced MSD parallel and normal to the walls as a function of reduced time in Fig. 7 along the same isomorph and isochore as before. The MSD is averaged over the entire slit-pore. We find excellent invariance along the isomorph and visible variation for the isochore. For the bulk liquid, at these high temperatures, one would see a similar scaling in comparison to the isochore; the differences becoming more pronouced with the degree of supercooling.

Two-point spatial correlation functions, such as the radial distribution function, have been shown to a good approximation to be invariant along isomorphs in simulations paper4 ; ingebrigtsen2013 ; ingebrigtsen2015 ; exp1 . Higher-order structural correlations have been investigated to a lesser extent, and could be less invariant than the two-body correlation functions as the isomorph theory is approximate. In particular, geometric motifs, so-called locally favoured structures (LFS) such as the bicapped square antiprism coslovich2007 ; malins2013fara have been seen to vary by a factor of two along isomorphs in supercooled KABLJ malins2013iso .

As a final probe of isomorphs for the confined SCLJ liquid we investigate in Figs. 8 and 9 invariance of higher-order structures. Now the liquids here are not supercooled much and thus rather than a locally favoured structure, we instead consider minimum energy clusters of $5\leq m\leq 13$ particles for these systems taffs2010 ; malins2013fara . These minimum energy clusters typically include the locally favoured structure, although the latter typically exhibits a specific symmetry royall2015a while the minimum energy clusters for each system exhibit a range of symmetries. In addition, we consider populations of the hcp and fcc crystalline structures.

The minimum energy structures are identified by the topological cluster classification malins2013tcc . The topological cluster classification (TCC) algorithm carries out a Voronoi decomposition and seeks structures topologically identical to geometric motifs of particular interest. The eight minimum energy structures of the SCLJ system are depicted in the figures (top left in each panel) and which are minimum energy structures of the SCLJ system ben2018 ; malins2013tcc ; taffs2010 .

We find that the distribution of clusters is, to an excellent approximation, invariant along the isomorph but not along the isochore. For all structures on the isochore the variation is around a factor of two with hardly any visible variation along the isomorph. Very minor deviations are, however, noted for the 9B and 11C structures along the isomorph. We emphasize that the probing of these structures does not imply relevance of these structures for the dynamics of the liquid, but is merely used for testing invariance of higher-order structures along an isomorph.

While the structures considered exhibit very little change along the isomorph, the changes between the behaviour of the different clusters is notable in itself. We can identify three regimes: small amorphous clusters, larger amorphous clusters and crystalline structures. As shown in Figs. 8(a) and (c), the smaller amorphous clusters, the triangular bipyramid 5A and octahedron 6A largely follow the density profiles illustrated in Fig. 5. Larger amorphous clusters, beginning with the pentagonal bipyramid 7A have a degree of fivefold symmetry, and their population is suppressed close to the wall (see Figs. 8(e) and (g) and Figs. 9(a) and (c)). Interestingly, this is not observed in the case of a free interface where the cluster population is rather slaved to the density profile godonoga2011 . A rather different behaviour is found for the crystalline structures, where the layer by the wall has a high population of particles in a crystalline environment but the population in the middle of the slit is very small (Figs. 9(e) and (g)).

We compare these results with bulk populations of the same clusters for the same state points (temperature and density) as shown by the colored data points in Figs. 8 and 9. The data points are plotted for the isochore data, but may be taken to be representative of the isomorph data on the left hand side of the figure. Even in the centre of the slit, the results show a strong enhancement of cluster population with respect to the bulk in all cases except for the hcp and fcc crystals, whose population in the centre of the slit is negligible. This is remarkable, given that for a free liquid–vapour interface, cluster populations reach their bulk value with around a diameter from the interface godonoga2011 . Further work is called for to understand this unexpected increase in higher-order structure in confinement.

In the Supplementary Material (SM) we present results for an isomorph with $H\approx$ 6. We find also here excellent scaling of density profiles and mean-square displacements along the isomorph. To conclude on the SCLJ liquid, we find that isomorphs survive into heavily confined systems with fcc crystalline walls and have excellent scaling even for higher-order structures.

Figure 10 displays $R$ and $\gamma$ as a function of temperature and several slit-pore widths $H$ for the KABLJ mixture at $\rho_{\rm L}$ = 1.2 and $\rho_{\rm W}$ = 1 with $\epsilon_{\rm AW}$ = 1, and $\epsilon_{\rm BW}$ = 0.707. The bulk correlation coefficient at this density is $R\approx 0.96$. As for the SCLJ liquid we find an increase in $R$ and a decrease in $\gamma$ with temperature. Similary $R$ increases with increasing $H$ but even for $H=4$ is the confined liquid R-simple.

Next, we consider the effect of changing the liquid-wall interaction strength $\epsilon_{\rm LW}$ in Fig. 11. A maximum is again noted for both $R$ and $\gamma$ and is displaced away from the value of $\epsilon_{\rm AW}=1$; The location of the maximum seems to be more dependent on $H$ than for the SCLJ. For the KABLJ mixture the effect of surface roughness $\rho_{\rm W}$ was not investigated.

For the KABLJ mixture we also investigate two isomorphs with $H\approx$ 6 and $H\approx 10$ to facilitate comparison with the SCLJ liquid in the previous section. Figures 12(a), (c), and (e) show reduced normal density profiles for both particles in the mixture ($A$ and $B$), as well as the total density along an isomorph with 17$\%$ density increase and $H\approx 6$. Figures 12(b), (d), and (f) show the corresponding quantities along an isochore. Good invariance is noted along the isomorph but not along the isochore where again significant changes in the peak heights are noted for all density profiles, in particular for the $B$-particle density profiles (Fig. 12(f)).

Figure 13 shows reduced in-plane total density profiles for the layer closest to the wall, i.e. $|\widetilde{z}|\approx$ 2.7, for the first and last state points of the same isomorph and isochore. For the KABLJ mixture the in-plane density profile does not seem to show a strong deterioration along the isochore as found for the SCLJ liquid, which could be anticipated from the height variation of the first peak in Fig. 12(b). Nevertheless the invariance is still visually worse than the isomorph.

For the reduced normal and parallel A-particle MSDs in Fig. 14 almost perfect scaling is observed along the isomorph while approximately a decade deviation in diffusion coefficient is observed for the isochore. These deviations in MSD are similar to what is seen for supercooled bulk liquids paper4 ; ingebrigtsen2015 . We find for the B-particles a very similar scaling behaviour (not shown).

Finally, we consider invariance of selected minimum energy clusters for the KABLJ mixture malins2013tcc ; malins2013fara in Figs. 15 and 16, reaching a similar conclusion as for the SCLJ liquid with excellent invariance along the isomorph but not along the isochore showing around a factor of two variation in almost all structures. Although, the bicapped square antiprism (11A) structure has been shown to correlate reasonably to the dynamics of the KABLJ system peter2015 ; malins2013fara ; speck2012 ; turci2017 we find very few bicapped square antiprisms at the higher temperatures considered here (the onset temperature for glassy dynamics, at which the bicapped square antiprism become popular is around $T_{\mathrm{on}}\approx 1.0$) and it is therefore not included in the figures.

In particular, it is quite remarkable that the minimum energy clusters of the KABLJ system malins2013fara ; malins2013tcc show such a good invariance, even for higher-order correlations, given that previous work showed a significant discrepancy in precisely the same system malins2013iso , although at a lower temperature $T$ = 1.0. The results presented here is at a higher tempererature ($T>$ 2.0), the magnitude of the discrepancy, and its rather weak temperature dependence leads one to speculate whether the confinement may somehow influence the agreement.

In comparing with the bulk values (data points in Figs. 15 and 16 right hand side), we see that as above in the case of the SCLJ system (Figs. 8 and 9) that in many cases the cluster population even at the centre of the slit is markedly higher than the bulk liquid at the same state point. However this is not universally the case here, as the $m=7$ polytetrahedron (7K) in fact seems to sit right on the confined data for some state points, while at higher temperature the bulk population seems rather higher than the confined system of interest here. The reasons for the change in behaviour of this structure and, as noted above, why the minimum energy clusters typically have a reduced population with respect to the bulk is an interesting topic for future work.

We provide in SM figures for an isomorph with $H\approx$ 10. Similar conclusions are reached showing again good invariance along the isomorph.

There are no conflicts of interest to declare

B.M.G.D.C., J.C.D., and T.S.I. are supported by the VILLUM Foundation’s Matter Grant (No. 16515). B.M.G.D.C. acknowledges funding from EPSRC with grant code EP/L016648/1. C.P.R. acknowledges the Royal Society and European Research Council (ERC consolidator grant NANOPRS, project number 617266) for financial support.