History-dependent Stress Relaxation of Liquids Under High-confinement: A Molecular Dynamics Study

When liquids are confined into nanometer-scale slit, the induced layering-like film structure allows the liquid to sustain non-isotropic stresses and thus being load-bearing. Such anisotropic characteristics of liquid under confinement arise naturally from the liquids’ wave number dependent compressibility that does not need solidification to take place as a prerequisite. In other words, liquids under confinement can still remain fluidity with molecules being (sub-)diffusive. However, the extensively prolonged structural relaxation time can cause hysteresis of stress relaxation of confined molecules in response to the motions of confining walls and thereby yield the quasi-static stress tensor history-dependent. In this work, by means of molecule dynamics, the discrepancy of stress tensor of a highly confined key base-oil component, i.e. 1-decene trimer, is captured after its relaxation from being compressed and decompressed. The results indicate that among the effects (e.g. confinement, molecular structure, and film density) that can potentially affect confined stress tensor, the ordering status of the confined molecules plays a predominant role.


Introduction
Liquid under confinement is a subject of central importance in diverse fields, such as tribology and lubrication, porous media and biomembrane, and electrochemistry, that has been studied extensively during the past several decades [1][2][3][4]. When liquids are confined into a nanometer-scale slit, they often exhibit layered structures [5,6], much enhanced shear viscosity [7][8][9], and oscillatory solvation forces while being continuously compressed [5,10,11]. Those characteristics, deviating from the bulk, are typically due to perturbation (i.e., the presence of confining surfaces) induced molecular ordering, which is inherently a natural consequence of liquid wavenumber dependent compressibility [12,13]. In the framework of a density-functional theory based description [12][13][14][15], one may discern that one-particle density of a confined liquid holds a similar functional form as that of the two-particle density of a bulk liquid. Such a similarity has been confirmed by recent simulation works [16,17].
Compared to the Brownian motion of bulk liquids as described by the Rouse model [18], mobility of liquid under high-confinement (with only a few molecular layers) can be largely impeded. Instead of random walk, motions of confined molecules turn to be collective (known as molecular ordering or orientational correlation) that leads to a prolonged structural relaxation time up to several orders of magnitudes [19,20]. In the case of large surface area to film thickness ratio, e.g., when using surface force apparatus (SFA) [21], the confining/shear rate typically, if not significantly, exceeds the inverse structural relaxation time with the discrepancy increasing with the number of confined monolayers decreases. Effective shear viscosity (η eff ) increases correspondingly given such an inverse dependence on the relaxation time (τ relax ), i.e., η eff ∝ τ −1 relax , beyond which liquid flow obeys Newton's law. Reversible phase transformation can take place when a liquid is close to a liquid-solid or liquid-glass transformation [4,[22][23][24], for which, however, one should be aware that confinement alone [24] cannot be served as the sole explanation since, for example, the effects of thermodynamic anomalies are omitted [25]. In fact, behavior of a liquid confined within a finite space shall be affected by a series of cross effects including molecular structure and chemical composition [8,24,26], confining/shear rate [27], the nature of confining walls [8,11,28,29] and dimensions [30], contamination [31], and so on.
Solidification is not a necessity for a confined liquid to manifest solid-like behavior (e.g., a finite yield point and critical shear stress [32]), particularly when the liquid has complex molecular structures. It has been reported [17,19,33,34] that highly confined liquids can remain fluidity in a quasi-static state with molecules being sub-diffusive. The heterogeneity of liquids under confinement, in light of the oscillatory film density, allows them to possess non-isotropic stress tensor and thereby to sustain shear and load irrelevant to the surface tension [35]. Our previous work [17], for example, shows that a triple-layer boundary film (made of 1-decene trimer) can sustain normal loads on the order of 100 MPa over sufficiently long times that requires neither crystallization nor a glass transition to occur. Note, merely 'anisotropy' cannot guarantee a sustainable load-bearing ability as flow may still occur in response to anisotropic stresses. In fact, for many liquids, such a load-bearing ability builds up rather continuously with a decreasing number of the confined monolayers [5]. Interestingly, we also found that mobility of confined molecules can rise (against an overall trend) with increasing confinement right after a monolayer molecules is expelled from the central confined region [17]. This observation suggests that, in addition to confinement, confining pressure can matter as well [7]. It too explains why a less dense film is able to endure higher loads than a denser film [34], which is far-fetched to be simply understood as a competition of disorderly statistics of polymer chain and specific surface-segment interaction.
Drainage dynamics of liquid under confinement stems from the change of confining geometry induced molecular restructuring, which can be history-dependent given the prolonged structural relaxation time [27]. As an extension of our previous work [17], the focus herein is to investigate the hysteresis of confined stress after the liquid relaxed from being compressed (from n + 1 to n monolayers) or decompressed (from n − 1 to n monolayers) by a blunt tip using molecular dynamics (MD). In the two cases studied, the initial pressure differences between the central confined region and the ambient drive the liquid to flow in the opposite directions, to which the aim is to clarify to what extent the stress anisotropy persists and whether the relaxed stresses end up comparable. Connections between the stress status and molecular morphology is established to clarify the origin of confinement-dependent stress discrepancy.

Methodology
The simulation model is shown in Fig. 1, where the liquid is sandwiched between a blunt ridge and a flat substrate. Here we continue to focus on a branched hydrocarbon, namely 1-decene trimer (C 30 H 62 , a snapshot of a single molecule is shown as the inset of Fig. 1), as in our previous work [17] since it is an important base-oil component and less likely to order compared to linear alkanes. Both the tip and the substrate are made of iron with (100) surface and are atomically flat. Using a flat punch instead of a curved tip simplifies the problem and improves statistics, but it can also lead to a slightly increase of the free energy barrier when a monolayer is expelled. Molecules can be squeezed into (or sucked from) a reservoir with a nominally constant chemical potential [36] from the sideways along the x-direction. The in-plane dimensions of the substrate and the bottom of the ridge are 20.0×5.1 nm 2 and 9.9×5.1 nm 2 , respectively. There is a total of 450 liquid molecules in the model system, giving rise to an adsorbing liquid film of an approximately six-monolayer thickness in the absence of a counter surface. The substrate and the bottom of the ridge (in mauve) are strongly adhesive to the liquid molecules whereas other solid walls (in gray) are repulsive with a cutoff of 6 √ 2σ (σ is the distance at which the pairwise potential energy is zero). Atoms from the solid walls remain fixed relative to their centerof-mass velocities throughout the simulations provided that deformation of the solids is negligible. The inter-and intra-molecular interactions of the liquid are described by a unitedatom variant-the transferable potentials for phase equilibria (TraPPE) [37,38] force-field, in which monomers, i.e., CH 3 , CH 2 , and CH, are represented by single interaction sites (pseudo-atoms). This force-field can faithfully reproduce experimentally measured thermophysical properties [17,39] and significantly lower the computational cost compared to all-atom variants [40,41]. The intra-molecular (bonded) interactions contain harmonic energy expressions for bond stretching, bond angle bending, and torsional motions. As suggested [42], bond stretching is allowed in order to facilitate parallel computing in MD, where a reasonable spring constant is adopted from OPLS force-field. While, the intermolecular (non-bonded) interactions are governed by 12-6 Lennard-Jones (LJ) potential with a cutoff of 1.0 nm. Long-range Coulomb interactions are not taken into account because the pseudo-atoms are electronically neutral. The Lorentz-Berthelot (L-B) mixing rules [43,44] are applied for unlike atom pairs, including both the liquid-liquid and liquidsolid interactions [45]. Detailed force-field parameters and functional forms can be found in Refs [37,38,45]. Note, the purpose here is prone to construct a boundary condition to mimic macroscopic observations rather than to accurately simulate iron-hydrocarbon interactions at the interface, otherwise the iron surfaces would have been passivated with oxygen [46].
There are two sets of simulations carried out: (1) compress/decompress the liquid at a constant rate, and (2) relax the liquid at fixed separations. Regarding (1), the displacing rate of the indenter is 0.2 m/s (along z-axis), which is the lowest rate used in our previous work [17]. In the case of decompression, unloading of the indenter takes place either right after it reaches the lowest position (at which the separation between the tip and the substrate is ∼0.8 nm) or after halting the indenter still for 100 ns. Regarding (2), the liquid is relaxed at a fixed separation d (d=1.68, 1.48, and 1.38 nm) for a total duration of 105 ns after being compressed or decompressed. At each separation, a film with a total number of n=3 monolayers is retained. The relaxation duration is chosen based on the fact that all the confining stresses reach ostensibly stable (quasi-static) states at the end of the simulation. Periodic boundary conditions are applied in the xy-plane; in the z-direction, atom movements are confined by the solid walls. The system is maintained at 300 K using a Langevin thermostat with a damping time constant of 100 fs. The thermostat is only applied in parallel to the y-component of the velocities of liquid molecules so as not to interfere the flow dynamics and molecular diffusion. The simulation time-step is 2 fs. Results reported in this work are each averaged over three independent simulations in terms of different seeds. All the simulations are carried out using the open-source code LAMMPS [47].

Solvation force
Solvation forces (F z ), measured as the normal forces on the indenter, are shown as functions of wall-wall separations (d) in Fig. 2, where results from the two cases appear contrasting. First, forces from prompt decompression (red squares) are predominantly attractive as opposed to the solely repulsive forces during compression (black circles). Given the opposite flow directions and viscous nature, this is mainly because the inertia forces from molecules resist tip motion while the molecules being either squeezed or sucked. Second, undulation, as an indicative of transition of monolayers from (n±1) to n, appears only when the molecules are expelled (during compression). Such collective behavior of molecules in response to the external forces, often being described as layer-by-layer transition [16,48], is, however, not detected in the case of decompression even after the system pre-equilibrated for 100 ns before unloading (blue dashed curve). This simply suggests that the current displacing rate (0.2 m/s) is much too fast for molecules to flow and then to restructure into analogous arrangements within the confining geometry. Also irrelevant to the lowest tip position, discrepancies between the curves in Fig. 2 will remain exist as long as the deforming rate surpasses the inverse of structural relaxation time of molecules confined. In that sense, the transient solvation forces will be dependent on the moving history of the tip.

Stress relaxation
To study stress anisotropy and its sustainability, evolution of the stress tensor (σ) of liquid at the central confined region is monitored at fixed separations (d 1−3 ). Using the virial estimator [49], the stress tensor components are calculated right after the indenter halted from either a compression or a decompression run with the number of monolayers evolved from n±1 to n (n=3). The virial estimator includes kinetic energy contribution and virial contribution that consists both the intra-and inter-molecular interactions. Although the virial theorem is known to be rigorous for homogeneous systems in thermal equilibrium [50], it predicts, according to our tests, fairly close mean stresses of a multi-layer film as those from the method of planes (MOP) [51]. Also, the normal stresses predicted using virial estimator agree reasonably well (the errors are less than 5%) with the mechanical definition of the normal pressure: P z = F z /A, where F z denotes the total normal force on the indenter and A the contact area. Stress anisotropy, detected on both simple [52] and complex [17] molecules, can be attributed to the heterogeneity of liquids under confinement. Compared to being compressed (as demonstrated in our previous work [17]), the anisotropic characteristics of liquid after decompression appear less pronounced (not show) due to the lack of molecular ordering (discuss next). Figure 3 shows a comparison of normal stress (σ zz ) evolution after compression (solid symbols) and decompression (open symbols) runs, where stresses from each pair of cases first rapidly approach each other from the opposite ends, then gradually evolve into nominally stable states after approximately 50 ns. Notably, the discrepancy between the quasi-static pair stresses (at ∼100 ns) rises from trivial to somewhat significant as the confining gap decreases. In specific, according to the fitting curves, the stress differences read 6.1, 28.6, and 77.1 MPa for the gaps of d 1 , d 2 , and d 3 , respectively. Those discrepancies, as one can imagine, would slowly decay with time until finally disappear once the corresponding structural relaxation times are exceeded, provided that fluidity is remained within the confined zone [17]. However, the corresponding time-span can easily beyond the typical time-scale of MD [8], which is assuredly out of the scope of the current study. Nevertheless, before reaching such a threshold, the ability of a confined liquid to sustain non-isotropic stresses persists and ought to be dependent on both the confinement and the moving history of the tip.

Radius of gyration
The aforementioned quasi-static history dependence of stress tensor and earlier claimed stress anisotropy [17] are inherently determined by the detailed confined molecular arrangement, which is not readily accessible by the experimental techniques. To elaborate such correlations, radius of gyration (R g ) of central confined molecules is calculated based α . In these calculations, only molecules located within the mid-monolayer are taken into consideration for reasons to exclude boundary effects from the confining walls. In the expression, α denotes Cartesian coordinates, δr is the relative position between monomer j and the center-of-mass of molecule i it belongs, N mol and N mon are the total numbers of embraced molecules and monomers in each single molecule, respectively, angle brackets denote an ensemble average. The principle moments of R 2 g at d 1−3 during the last 25 ns relaxation are shown in Fig. 4.
Within the confining space, liquid molecules are highly deformed and being parallel to the confining walls, as opposed to those in a bulk phase of which the three principle components of R 2 g are comparable (R 2 αα ≈ 10.8 Å 2 for 1-decene trimer) [17]. As shown in   xx and R 2 yy ) are at least two-fold as that in the normal direction (R 2 zz ), indicating a sheet-like structure that complements the observed oscillatory density profile. In the case of compression, the extension of molecules along the flow direction (x-axis) appears more sensitive to the confinement change given a rapidly increased pressure gap between the initial and stable states (see Fig. 3). High pressure induced strong tendency to flow can ease molecular alignment along the streaming direction and thereby yield stretched molecular appearances. In contrast, the bottleneck effect, particularly at high confinement, would hamper the travelling of disordered molecules from an open space into a narrow channel (decompression case), resulting in a relatively low film density in the central void (discuss next). Under low confining pressure, the inert response of molecular conformation to the confinement change suggests that a much longer time is required for molecules to replenish the central void and become orderly arranged thereafter. This explains why there is no sign of collective behavior (e.g., undulations on the solvation force curve) in the course of continuous decompression, as shown in Fig. 2.

Mean-squared displacement
Diffusion of central confined molecules at d 1−3 is characterized via measuring the in-plane mean-squared displacement (MSD) of monomers in the mid-confined-layer: MSD is the position of the j-th monomer from the i-th molecule at time t (t ≥ t 0 = 80 ns). A comparison between the results from compression and decompression runs is shown in Fig. 5. For all the cases studied, the confined liquid remains fluidity with molecules being sub-diffusive given that the smallest mean displacement (∼4 Å) of monomers over a 25 ns time window exceeds the typical diffusion rate of polymers on surfaces [53] and atoms within crystalline structures, even those located along high-diffusivity paths (dislocations, grain boundaries, and free surfaces) [54]. Clearly, at a given separation, molecules relaxed from decompression are more diffusive compared to those after compression, particularly evident in the case of high-confinement (see the blue curves in Fig. 5). It is the different film densities and thereby molecular morphologies render the mobility of confined molecules confining history dependent. It ought to clarify that higher molecular mobility in the central confined region does not correspond to an easier incoming flow at the entrance of the confined channel (i.e., the edges of the tip), but the opposite instead. Moreover, as pointed out in our previous work [17], diffusivity of confined molecules may not necessarily hold a monotonic function of confinement given the fact that it is the confining pressure that plays an essential role. This claim is hereby endorsed by the results shown in Fig. 5, where molecules at a higher confinement (blue dashed curve) can be more diffusive than those at a lower confinement (red curves for example).  Figure 6 shows confined liquid density profiles (ρ) along the normal direction (z-axis) with a bin size of 0.5 Å. At each separation, three molecular layers, evolved from t=0 to 100 ns, can be identified (by counting the oscillation periods). Evidently, wider wall-wall separation (d 1 ) inhibits molecular transfer to a lesser degree given the analogous density profiles ending up after 100 ns relaxation (Fig. 6d). While, at high confinement (d 3 ), the notable discrepancy between the density peaks of the mid-layer (Fig. 6f) shall be the origin of the aforementioned dissimilarities in molecular morphology and stress tensor at rest. That being said, the lower the density peak, the less pronounced molecular ordering and thereby the lower confining stresses. Note, the correlation between confining stresses and film density is not ubiquitous. As an example, the average film density in the case of d = 1.68 nm at 100 ns (∼0.97 g/cm 3 ) is 24% higher than that of the bulk liquid density (0.78 g/cm 3 ) [17] at ambient conditions, whereas the normal stress appears attractive (the negative stress reads -27.9 MPa). Essentially, what determines the stress tensor of a liquid confined within a limited space is not merely whether the molecules are densely packed (mean film density) but rather to what extent they are orderly arranged (oscillatory peak height).

Conclusions
Heterogeneity of liquids under confinement enables them to sustain non-isotropic stresses that irrelevant to the surface tension. In this work, stress relaxation of a key baseoil component, 1-decene trimer, under nanoscale confinement was studied by means of molecular dynamics. In terms of a triple-layer film, a stress discrepancy was identified while the liquid relaxed from tip compression and decompression in nominally quasi-static states. Such a discrepancy, increasing with confinement, stems from the unalike ordering status of molecules in the central confined region that dominated by the confining pressure. Although highly relevant, factors such as confinement and film density are not intrinsic. The comparably lower normal stress in the case of decompression can be attributed to the impeding flow of disordered molecules from an open space into a narrow channel, thus a much longer time is required for molecules to evolve into an orderly arranged film structure. Such prolonged structural relaxation time can induce hysteresis of molecular restructuring in response to the external forces and thereby yields the stress tensor history dependent. Findings from this work, along with our previous study [17], may complement the understanding of drainage dynamic and mechanical response of liquids under smallscale confinement.