Correlation between ordering and shear thinning in confined liquids

Despite the extensive research that has been conducted for decades on the behavior of confined liquids, detailed knowledge of this phenomenon, particularly in the mixed/boundary lubrication regime, remains limited. This can be attributed to several factors including the difficulty of direct experimental observations of the behavior of lubricant molecules under non-equilibrium conditions, the high computational cost of molecular simulations to reach steady state, and the low signal-to-noise ratio at extremely low shear rates corresponding to actual operating conditions. To this end, we studied the correlation between the structure formation and shear viscosity of octamethylcyclotetrasiloxane confined between two mica surfaces in a mixed/boundary lubrication regime. Three different surface separations corresponding to two-, three-, and five-layered structures were considered to analyze the effect of confinement. The orientational distributions with one specific peak for $n=2$ and two distributions, including a parallel orientation with the surface normal for $n>2$, were observed at rest. The confined liquids exhibited a distinct shear-thinning behavior independent of surface separations for a relatively low sliding velocity, $V_{\rm x}\lesssim 10^{-1}\,{\rm m/s}$. However, the shear viscosities at $V_{\rm x}\lesssim 10^{-1}\,{\rm m/s}$ depended on the number of layered structures. Newtonian behavior was observed with a further increase in the sliding velocity. Furthermore, we found a strong correlation between the degree of molecular orientation and the shear viscosity of the confined liquids. The magnitude of the shear viscosity of the confined liquids can primarily be determined by the degree of molecular orientation, and shear-thinning originates from the vanishing of specific orientational distributions with increasing sliding velocity.


I. INTRODUCTION
Energy-saving techniques such as friction reduction have become increasingly important as countermeasures against global warming. 1 Particularly, an improvement in the friction reduction 2 for mechanical systems, including automobile systems, is essential because the internal combustion engine is still used as the primary power source for most vehicles, and various tribological contacts such as electric motors, wheel bearings, and steering systems are still present in electric vehicles. 3 In the analysis of lubrication problems, the behavior of the lubricant molecules in a wide range of fluid lubrication to boundary lubrication must be considered based on the practical operating conditions of mechanical systems (e.g., engine oil experiences an engine starting/idling state (boundary lubrication) and a high-speed operating state (fluid lubrication)). The Stribeck curve 4 is a fundamental concept of tribology. This concept describes the friction coefficient as a function of the Sommerfeld number (viscosity × sliding velocity/load). In the fluid lubrication region, the frictional force is determined by the bulk properties of the lubricant. Therefore, the use of a lower bulk viscosity can effectively reduce the friction force. A transition from fluid lubrication to a mixed/boundary lubrication region is observed with the decrease in the sliding velocity and/or increase in the load, and lubricant films are formed under extreme confinement at the nanometer scale. Therefore, obtaining a better understanding of the bulk properties along with the structural and shear behavior of the lubricant under nanoscale confinement is a) Electronic mail: kobayashi@mech.keio.ac.jp crucial.
Confined liquids have attracted considerable attention in the past few decades owing to their potential applications in electronics 5 , engineering 6 , and biomedical 7 fields. The molecular arrangement and orientation in the static state, and the dynamic properties can be measured on a length scale from µm to nm owing to the progress in experimental techniques such as atomic force microscopy (AFM) 8 , surface force apparatus (SFA) 9 , and resonance shear measurement (RFM) 10 . Confined liquids exhibit distinct phases 11,12 and flow properties 13,14 not exhibited in the bulk solution even for low-molecular-weight solutions such as water owing to the effect of spatial constraints and solid-liquid interfaces. Confined water is an interesting system because it is ubiquitous in nature; examples include biochannels 15,16 and water films confined between clay mineral surfaces 17,18 . Pure water and lubricating oil with trace quantities of water are also potential candidates for future lubricants due to their eco-friendly nature. 14,19,20 A recent experiment 20 conducted on hexadecane with trace quantities of water confined between mica walls demonstrated that the nucleation of water nanodroplets at the surface is induced by an increase in the temperature of the system, and the addition of amphiphilic molecules such as oleic acids to the lubricant suppresses this nucleation. Additionally, research conducted on the viscosity behavior of four different phenyl ether lubricants under nanoslit confinement using RSM indicated a significant increase in the viscosity with the decrease in the surface separation ( 2 nm), and the magnitude of the viscosity correlation between the bulk and the confined systems was reversed. 21 Despite the extensive research conducted on the development of measurement techniques for the properties of confined lubricants, obtaining a better understanding of the correlation be-tween the structure and flow properties of lubricants under nanoconfinement remains a challenge owing to the difficulty of obtaining direct experimental observations of the behavior of lubricant molecules, particularly under non-equilibrium conditions. Molecular dynamics (MD) simulation is a powerful tool which can be used to obtain microscopic insights into the properties of confined liquids behind the experimental results and to guide future research applications. Several studies have been conducted using both equilibrium and non-equilibrium MD simulations on confined liquids such as Lennard-Jones (LJ) particles 22-24 , water 11,14,25,26 , and hydrocarbons [27][28][29] . A recent simulation conducted by Ta et al. 27 of a hexadecane film confined between iron oxide surfaces demonstrated that an increase in the shear rate can result in a decrease in the degree of ordering of the hexadecane film and an increase in the velocity slip at the solid-liquid interfaces. A combination of tribometer experiments and MD simulations 30 demonstrated the effect of molecular structure on the friction behavior. The study reported that flexible and broadly linear hydrocarbon molecules exhibit low friction coefficients which increase with the strain rate and pressure, whereas inflexible molecules that contain cycloaliphatic groups have high friction coefficients that are almost independent of the strain rate and surface pressure. However, the task of computing non-equilibrium MD simulations in the shear rate range corresponding to those used in practical operating conditions is extremely difficult, which hinders the possibility of obtaining a detailed understanding of the underlying physics despite the extensive simulation research conducted on this topic. Therefore, most of the previous simulation efforts have been focused on friction and flow behavior in the elastohydrodynamic lubrication regime.
This study aims to provide new molecular-level insights into the correlation between the structural and lubrication properties in the mixed/boundary lubrication regime. Octamethylcyclotetrasiloxane (OMCTS) and cyclohexane are among the most widely used lubricant molecules in practical lubrication systems; these model lubricants are more frequently used in tribometer experiments because of their globular and non-polar nature. Therefore, these molecules have also been used as model lubricants in molecular simulation studies conducted on confinement. Hitherto, force-field parameters for OMCTS [31][32][33] and cyclohexane 34,35 have been developed to reproduce the layered structure and period of oscillatory force observed in experiments. In addition to the development of potential models, analyses of the transport properties 36 and the possibility of freezing transition 37 under nanoscale confinement have been conducted via equilibrium MD simulations. A confinementinduced diffusion slowdown was observed for OMCTS confined between mica surfaces, and the mechanism of this phenomenon was explained using the concept of activated diffusion theory 36 . As an exception, a recent simulation of cyclohexane confined between mica surfaces by Xu et al. 38 presented the first molecular evidence for the structural changes and dynamics of stick-slip motion in boundary lubrication.
In this study, we performed MD simulations to demonstrate the correlation between structure formation and shear viscosity of a lubricant film under nanoscale confinement in a mixed/boundary lubrication regime. We focused on the OMCTS-mica system, given the progress made in lubricant models to accurately reproduce a realistic view of the microscopic structure of the confined system. By using this system, our work determines the physical mechanism of the behavior of confined liquids in the mixed/boundary lubrication regime from a molecular viewpoint.

A. Simulation model
We use the MD method to analyze the structural formation and flow properties of OMCTS confined between two mica walls in equilibrium and under shear. In the employed OMCTS model 31 , both the crystal lattice constants and liquid transport properties such as bulk density and diffusion coefficient in the range of 300-440 K are reproduced, despite only considering the methyl-methyl interactions. The atomic sites interact with each other through the standard LJ potential, which is given as: where r ij denotes the distance between atomic sites, σ ij is the length parameter, ε ij is the energy parameter, and r cut = 10.0Å 32 is the cutoff radius. All the LJ parameters used in this study are listed in Table I.

B. Simulation procedure
We performed MD simulations using a liquid-vapor (L-V) molecular ensemble proposed by Leng 39 to prepare the confined system. The system setup for the L-V MD simulations of the OMCTS-mica system was established based on previous studies 32 , as shown in Fig. S1 in the ESI † . The mica surface consists of 614 O and 288 K + atoms. The cleaved mica surface contains half of the K + ions present in the crystal. The 1100 OMCTS molecules were sandwiched between two mica walls with a surface area of 81.8Å×81.2Å, and the two mica walls were connected to smooth repulsive walls expressed by harmonic potentials. The simulation box was 400Å×81.2Å×110Å, and the periodic boundary condition was applied in all three spatial conditions. The complete details of the procedure have been provided in Ref 32 . This atomistic model for the OMCTS-mica system can reproduce the layered structure and period of oscillatory force observed in the experiments.
We calculated the normal surface pressure (force/area)-distance curve at 1Å intervals by continuously moving the upper mica with dz/dt = 0.01Å/ps to obtain confined systems with a targeted number of layers (n). The initial configurations were prepared by randomly inserting OMCTS molecules between two mica walls with H = 63Å using the Packmol software 40 . The temperature was maintained at 300 K using a Nosé-Hoover thermostat 41 with a time constant of 0.2 ps. The atoms of the upper and lower mica were frozen, and each OMCTS molecule was considered as a rigid body 42 with a time step of 2.0 fs. For each surface separation, H, an equilibrium run was carried out for at least 1 ns before performing a production run to calculate the normal force F z . We selected H from the oscillation peak of the target layer region after deriving the normal surface pressure-distance curves from the L-V MD simulations; however, the value of H depends on the degree of surface pressure (see Fig. S2 in the Supplementary Materials). For example, the surface separation for a five-layered system is determined according to the following procedure. First, we selected H from the oscillation peak of the five-layer region that corresponds to the normal surface pressure, which can be obtained from the surface force experimental measurements. A five-layered system was obtained at a normal surface pressure, P zz , of 0.1 MPa in the experiments 43 . Here, we identified two values of H at P zz = 0.1 MPa within the oscillation peak of the five-layer region. Previous MD simulations of the LJ particles confined within a nanoslit have demonstrated that the layered structures at larger values of H within the oscillation peak are energetically more stable, resulting in the increase of the freezing/melting points. 22 Therefore, a larger H value is selected in this study to reflect practical experimental conditions. Subsequently, the constraint of the z-position of the upper mica is released, surface pressure is applied, and the time evolution of the z-position of the upper mica is measured until a plateau is reached. Lastly, the OMCTS molecules outside the mica walls were removed from the system, and the equilibrium simulations were repeated. Confined systems with a targeted number of layers (n = 2, 3, and 5) were obtained using this procedure, and the equilibrium properties of the system are explained in Section III A.
After performing the equilibrium simulations, shear flow was imposed by applying a constant sliding ve- locity, V x , on the upper mica. The fluid temperature was controlled at 300 K by applying the velocity scaling method to the thermal velocity, which was subtracted from the streaming velocity. The applied V x ranged from 10 −3 m/s to 10 2 m/s, leading to shear rates,γ, in the range from 1.32 × 10 5 s −1 γ 3.57 × 10 10 s −1 for the film thickness simulated. The given low sliding velocities (V x < 10 −2 m/s) are of particular interest because this range enables an overlap with those used in the tribometer experiments. Conversely, there are significant fluctuations in friction force and velocity gradient at such low V x owing to the thermal noise. Furthermore, the calculation time required for the simulation to reach a steady state is also very long. Therefore, a moving average was applied to smooth out the data to determine the steady state and the shear viscosity, η, which is given by: where τ denotes the shear stress, which is measured by the friction force on the upper mica per unit surface area. We used the graphics processing unit (GPU) versions 44 of the large-scale atomic/molecular massively parallel simulator (LAMMPS) 45 MD package are used in this study.

A. Equilibrium behavior
We investigated the equilibrium properties of the confined systems with layered structures (n = 2, 3, and 5). Figure 1 presents the density profiles of the centers-ofmass of the OMCTS along the direction normal to the mica surface for different surface separations. We can clearly observe the layered structures can be clearly observed in the confined system for all values of H. Zero density is observed between the layers at H = 2.2 (n = 2) and 2.8 (n = 3). However, a non-zero density is observed at H = 3.8 (n = 5), indicating a loose layered structure. Previously conducted simulation studies 32,36,37 have reported similar behavior. Table II lists the resulting equilibrium properties for each H value, including the surface pressure P zz , number densities ρ z , and diffusion coefficient D. In a previously conducted molecular simulation study of the OMCTS-mica system 36 , the reason for selecting a specific H within a single layer region was ambiguous. We used the value of P zz obtained from the surface force measurements 43 to reproduce the realistic state of the confined system. These settings can be used to capture phenomena in practical systems more accurately, and the obtained results can be connected with the behavior observed in experiments. The diffusion coefficient, D, was calculated from the slopes of the time-averaged mean-square displacement curves. A previous study 36 reported a diffusion slowdown induced by confinement, where D was approximately two orders of magnitude lesser than the bulk value 31 . A considerable decrease was observed for H = 2.2 nm, which was attributed to the contact layers. For a confined system, the existence of contact layers causes a significant decrease in the diffusivity. 36 It is observed that D increased despite the decrease of the surface separation from H = 3.8 nm to H = 2.8 nm. Previous simulations 36 have also reported that the diffusion curve exhibits a non-monotonous decrease when considered as a function of surface separation. We found that this result is primarily attributed to the number density of the confined systems (see Table II).
The orientational distributions, P (cos θ), of the confined OMCTS, where θ is the angle between the Si and Si vector and normal to the mica surface, are computed to obtain the detailed structural information. A perpendicular orientation with the surface normal corresponds to cos θ = 0 and a parallel orientation with the surface normal corresponds to cos θ = 1.0. Figure 2(a) presents P (cos θ), of confined OMCTS for various surface separations. For H = 2.2 nm, a distinct peak is observed at cos θ ≈ 0.7. This specific orientation was also observed in the previous simulation, but the distribution shifted to cos θ = 0 when the surface separation was decreased further. 37 For larger separations with contact layers and middle layers (H = 2.8 nm (n = 3) and H = 3.8 nm (n = 5)), we observed a different behavior, where two distinct peaks appeared at cos θ = 0 and cos θ ≈ 0.7. The peak at cos θ ≈ 0.7 gradually weakens when compared to the case of H = 2.2 nm. A broader distribution was observed with the increase of H, indicating that the OMCTS molecules did not have any preferred orientation corresponding to the surface normal. A layer-by-layer analysis of P (cos θ) is performed for H = 3.8 nm (n = 5) to obtain a more detailed understanding of the difference between n = 2 and others, as shown in Fig. 2(b). It can be observed from this plot, that the first contact layers (L1) had sharper peaks when compared to the second layer (L2) and middle layer (L3). The peaks clearly became weak for L2 and L3, and a more uniform distribution is observed in the central region (from L2 to L3). The formation mechanism of the layered structures in contact layers that have two distributions is attributed to the molecular pair of contact layers and inner layers induced by the confinement effect, as explained by a previous study conducted by Matsubara et al. 32 . These characteristic structures induced by nanoconfinement signifi-cantly affect both the equilibrium properties and the flow properties, which will be discussed in the following section.

B. Flow behavior
This section presents a detailed description of the structure formation and viscosity behavior to elucidate the underlying physical phenomena of flow behavior for low sliding velocities, reflecting practical operating conditions. Figure 3 illustrates the shear viscosity, η, of OM-CTS confined between two mica surfaces as a function of the sliding velocity, V x , on the upper mica. For a relatively low sliding velocity, V x 10 −1 m/s, the confined liquids exhibited a distinct shear-thinning behavior for all the systems. When a higher sliding velocity, V x 10 −1 m/s, was applied, a Newtonian-like behavior was observed with η ≈ constant which corresponds to a slightly smaller zero-shear viscosity of the bulk liquid, η bulk = 2.13 mPa · s. Additionally, it was observed that the shear viscosities at V x 10 −1 m/s were largely dependent on the number of layered structures, n. η increases with the increase in H, particularly for H = 2.2 (n = 2).
The orientational distributions of the confined OM-CTS under shear, P (cos θ), are analyzed to better understand the mechanism behind this viscosity change, as shown in Fig. 4. The shear-induced rearrangements, such as enhancement of the orientation degree and breakup of the oriented structures, are observed for all H values. It can be observed that a relatively weak shear (V x 10 −1 m/s) observed at H = 2.2 nm (n = 2) leads to an increase in the distribution at cos θ ≈ 0.7 (see Fig. 4(a)). This is interesting because of a possible scenario for the increase in shear viscosity at low sliding velocities, which reflects the actual operating conditions. Additionally, the zero-shear viscosity data is also plotted for H = 2.2 nm using the Green-Kubo relation in Fig. 4, as follows: where P αβ denotes the off-diagonal components of the pressure tensor. While applying the Green-Kubo relation to confined systems is a controversial debate [46][47][48] because of the inability to adapt to heterogeneous systems, we computed the zero-shear viscosity in parallel (defined by axial pressure components P xy ) and perpendicular (defined by radial components P xz and P yz ) directions and then averaged them. The viscosity computed using the Green-Kubo formula in this study was almost an order of magnitude smaller than the bulk viscosity, which is similar to that reported in a previous study 48 .
The calculated results cannot be guaranteed in a confined system since both the liquid-liquid interactions and the solid-liquid boundary condition affect the thermal fluctuations of the shear stress in a confined liquid. However, the layered structures under low V x have a strong degree of orientation at cos θ ≈ 0.7 when compared to that in an equilibrium state, as shown in Fig. 4(a); the shearinduced rearrangement can cause an enhanced shear viscosity. It will be interesting to reveal the correlation between the shear viscosity and structure formation, such as the degree of order under extremely low sliding velocities (V x < 10 −3 m/s), in future work. Additionally, the development of theoretical prediction methods, as well as the issues of both large thermal noise and high computational cost to reach a steady state, will be the objective of future studies. We continue our discussion with the confined systems for n > 2. An increase in the distribution at cos θ ≈ 0.7 was observed at H = 2.8 nm (n = 3), even for a relatively high sliding velocity, V x = 10 0 m/s. In addition to the peak at cos θ ≈ 0.7, an increase in the distribution is observed at cos θ = 0; however, the peak intensity at cos θ = 0 was weaker than that at cos θ ≈ 0.7 (see Fig. 4(b)). The same behavior was observed at H = 3.8 nm (n = 5) as in the case of n = 3, but η was smaller than that in the case of n = 3 for the shear-thinning regime (V x < 5.0 × 10 −1 m/s). The three-layered system had a higher η value for the shearthinning regime (V x < 5.0 × 10 −1 m/s) than that in the case of n = 3 although the self-diffusion coefficient, D, of the three-layered system was higher than that of the five-layered system (see Table II). This indicates that the formation of oriented structures in a confined liquid is one of the factors which significantly affects the shear viscosity.
Subsequently, the origin of the shear-thinning of the confined liquids is discussed. A clear correlation can be obtained between the gradual breakup of the oriented structures and shear-thinning behavior by comparing Fig. 3 and Fig. 4. This can be more clearly observed from the result of H = 2.2 nm. The intensity of the peak at cos θ ≈ 0.7 gradually decreased along with shear-thinning as the sliding velocity increased to V x = 10 −1 m/s. The confined liquid exhibits a Newtonian-like behavior even though the peak at cos θ ≈ 0.7 remains at V x = 10 0 m/s. We anticipate the existence of a threshold peak intensity below which there is no significant effect on the viscosity behavior of confined liquids. For larger separations (n > 2), the slope of η becomes less steep when compared to the case of n = 2, as shown in Fig. 4. The effect of the specific orientation distribution of the contact layers on the shear viscosity reduces with the increase in H because of the existence of middle layers which have random orientations. From the peak intensity at cos θ ≈ 0.7 for all the analyzed H ( Fig. 4(a-c)) values, the threshold for the molecular orientation which causes the shear-thinning of the confined liquids was observed to be approximately 2 %. When V x was further increased (V x 10 0 m/s), we observed a broader orientational distribution for all the H values, indicating that the OMCTS molecules do not have any preferred orientation. These confined liquids continued to exhibit a Newtonian response to shear while breaking up the oriented structures. The correlation between the shear viscosity and structure formation corresponds to all the analyzed H values, and the obtained results indicate that the shear viscosity of simple nonpolar lubricants under nanoslit confinement is primarily determined by the degree of molecular orientation.

IV. CONCLUSIONS
We studied the correlation between the shear viscosity, η, and the structure formation of OMCTS confined between two mica walls using the MD method. Three different surface separations, H, reproducing two-(n = 2), three-(n = 3), and five-layered (n = 5) structures, were prepared by performing liquid-vapor molecular dynamics simulations. A confinement-induced diffusion slowdown is observed at rest, particularly for n = 2. For larger separations (n > 2), the self-diffusion coefficient, D, of n = 5 was higher than that of n = 3, and we found that this was because of the difference in the number density of the confined systems. In addition to the transport properties, the orientational distributions of the confined OMCTS molecules were also analyzed. We observed distinct distribution with one specific peak for n = 2 and two distributions, including a parallel orientation with the surface normal for n > 2. These distributions for n > 2 became more uniform when moving towards the central region in the confined system.
We then applied sliding velocities to the upper mica wall and generated shear flow covering the actual operating conditions in the range of 10 −3 m/s ≤ V x ≤ 10 2 m/s. The confined liquids demonstrated a distinct shear-thinning behavior for all the analyzed H values for a relatively low sliding velocity, V x 10 −1 m/s. Furthermore, it is observed that η at V x 10 −1 m/s depends on the number of layered structures, contrary to the results of the diffusion coefficient in equilibrium. Newtonian behavior was observed when a higher sliding velocity, V x 10 −1 m/s, was applied. A strong correlation was identified between the degree of molecular orientation and shear viscosity of confined liquids. By applying a low V x , we observed the shear-induced enhancement of orientation degree compared to the equilibrium state. The magnitude of η in the shear-thinning regime was found to be dependent on the orientation intensity. Furthermore, it was observed that shear-thinning originated from the vanishing of specific orientational distributions, corresponding to the gradual breakup of the oriented structures. The findings of this study provide new insights into the correlation between the structural and lubrication properties in the mixed/boundary lubrication regime and support breaking away from the trial-anderror empirical approaches for developing and improving lubricants. ogy (MEXT) as Research and Development of Next-Generation Fields. The computation was partially carried out using the computer resources offered under the category of general projects by the Research Institute for Information Technology, Kyushu University. The authors are grateful to Prof. K. Kurihara (Tohoku University) and Prof. M. Mizukami (Tohoku University) for helpful discussions and useful comments.
Supplementary Materials for "Correlation between ordering and shear thinning in confined liquids"