Thermodynamic and kinetic anisotropies in octane thin films

Confinement alters the thermodynamic and kinetic properties of matter. Such changes can partly arise from quantum effects at the nanoscale.^1–4 They can also emerge in purely classical systems, simply due to the presence of confinement-induced interfacial regions. On a fundamental level, confinement breaks the translational isotropy of a bulk material, making all its physical properties functions of position.⁵ The extent of such anisotropy is, however, variable, and can depend on both the thermodynamic conditions (temperature, density, composition) as well as the nature of the interface(s), such as boundary conditions (free interface vs. substrate), interactions (attractive, repulsive, inert), and dimensionality (pores, channels, slabs). Regardless, the behavior of confined matter can deviate significantly from the bulk in many different ways. Examples include the emergence of new phases,^6–8 changes in thermodynamic properties such as melting and glass transition temperatures,^9–13 and in kinetic properties such as nucleation rates,^14–16 viscosities,^17–20 diffusivities^21–23 and elastic constants.^18,24,25

The translational anisotropy that is induced as a result of confinement can be utilized for modulating the structural and functional properties of materials. A widely known example is heterogeneous nucleation²⁶ in which a substrate enhances a particular first-order phase transition in its vicinity by decreasing the associated free energy barriers. On a microscopic level, this is mediated by the formation of a subsurface region that is more conducive to nucleation than the homogeneous bulk material. The nucleating potency of a substrate is thus related to the microstructure of the material that is in its proximity, i.e., the anisotropy that it induces therein.²⁷ The relationship between the anisotropy of the subsurface region and the facilitated nucleation can be rigorously explained for phase transitions such as hydrophobic evaporation,^14,28 but is far more difficult to discern for more complex phase transitions, such as crystallization.¹⁶ 

A more recent example of the interesting behavior arising in anisotropic systems is provided by ultrastable glasses obtained by depositing the vapor of a glass-forming liquid onto a cold substrate. By tuning the substrate temperature, it is possible to form glasses that are far more stable than the ordinary glasses obtained by rapidly quenching the liquid.²⁹ In addition to their superior thermal and mechanical stability, such ultrastable glasses can be structurally anisotropic as evident from refractive index measurements and X-ray scattering.^30–32 Possible structural differences between ordinary and ultrastable glasses are manifest in the different scaling of their heat capacities in the limit of T → 0 K.³³ It has even been suggested that a first-order liquid-liquid transition might exist between the ordinary and ultrastable amorphous states.³⁴ Vapor deposition has been successfully used for making organic,^29,35–41 polymeric⁴² and metallic⁴³ ultrastable glasses. There is, however, a significant gap in understanding why these vapor-deposited glasses are ultrastable. It has been argued that the enhanced mobility of molecules at the vapor-liquid interface enables them to explore the potential energy landscape more efficiently, thereby giving rise to structures that reside far deeper in the potential energy landscape.⁴⁴ Such enhanced mobility has been observed in experimental^20,45–47 and computational^5,48–50 studies of thin films. However, no causal relationship has been unambiguously established between the existence of such accelerated regions and the formation of ultrastable glasses.

A systematic way of identifying a possible link between enhanced surface mobility and increased thermodynamic and kinetic stability is to probe thermodynamic and kinetic anisotropies in films of different compositions, and to assess their sensitivity to changes in thermodynamic conditions. It is of particular interest to study the role of substrates in inducing such anisotropies, and in modifying the properties of the resulting glasses. This is a question that has not been extensively studied, especially with regard to the vapor-deposited ultrastable glasses. Two particular parameters that are relevant in this quest are the temperature and the interaction parameter, with the latter quantifying the strength of (attractive) interactions between the substrate and the liquid. In our earlier publication,⁵ we performed a detailed analysis of such anisotropies in thin films of a model atomic glass-forming liquid,⁵¹ and we discovered two distinct dynamical regimes: accelerated dynamics near loosely attractive substrates, and decelerated dynamics near sticky substrates. We also established correlations between oscillations in density and in normal stress, and between oscillations in relaxation time and in lateral stress.

In this work, we revisit the main findings of Ref. 5, but now in the context of molecular thin films of n-octane. We choose n-octane as a prototypical chain molecule. Alkane thin films are particularly interesting as they can undergo a process known as surface freezing^52,53 in which a frozen layer emerges at the vapor-liquid interface at temperatures exceeding the equilibrium melting temperature. Also, alkanes constitute one of the most important components of crude oil, and their presence can result in interesting phase separation phenomena that are typically governed by confinement.^54–57 It is therefore worthwhile to inspect the microstructure of interfacial regions of alkane films with an eye towards understanding the mechanism(s) and predicting the kinetics of such phase transition processes.

This paper is organized as follows. Section II A provides details on the computational setup as well as the force-fields employed. Technical details of the simulations performed are given in Section II B. The procedures utilized for computing spatial profiles of thermodynamic and kinetic properties are given in Section II C. Section III A discusses the qualitative behavior of the octane films, particularly, with regard to surface freezing. Anisotropies in thermodynamic and kinetic properties are discussed in Sections III B and III C, respectively. Finally, Section IV is reserved for concluding remarks.

A schematic representation of the n-octane films considered in this work is depicted in Fig. 1(b). Octane molecules are represented with the NERD force-field⁵⁸ in which each octane molecule is comprised of eight united-atom interaction sites: two for CH₃ (green in Fig. 1(a)) and six for CH₂ (light blue in Fig. 1(a)). Our choice of a united-atom potential such as NERD is due to its simplicity, as we only consider n-alkanes as prototypical chain molecules. Indeed, the increased quantitative accuracy obtained from using more realistic, but computationally costly, all-atom force-fields is unlikely to impact the main findings of this work, and can only shift the n values or temperatures for which the reported phenomena are observed. In the NERD potential, the non-bonded interactions between these sites are modeled through the Lennard-Jones potential,⁵⁹ 

The bonded interactions include bond-stretching, bond-bending and torsional terms:

All model parameters are given in Table I. Lennard-Jones interactions are shifted and truncated at r_c = 1.38 nm.

In addition to the n-octane film, there are two substrates in the system: the attractive substrate (C atoms, steel blue in Fig. 1(b)) and the repulsive substrate (D atoms, light pink in Fig. 1(b)). The constituent atoms of both substrates are arranged into a face-centered cubic (fcc) lattice. They interact with themselves and the united-atom LJ sites via the Lennard-Jones potential with ϵ_AC = ϵ_AD = ϵ_Sϵ_AA, ϵ_BC = ϵ_BD = ϵ_Sϵ_BB, and σ_AC = σ_BC = σ_AD = σ_BD = 0.4 nm. Here, the interaction parameter, ϵ_S, is used for tuning the strength of attractive (repulsive) interactions between the substrates and octane molecules. All substrate-octane interactions are truncated and shifted at r_c = 1.2 nm and r_c = 0.44 nm for the attractive and repulsive substrates, respectively. The inclusion of a second repulsive substrate is to assure that the octane molecules evaporating from the film would never redeposit onto the opposite side of the substrate as a result of periodic boundary conditions. Simulations are carried out for three distinct values of ϵ_S = 0.5, 1.0, and 3.0, with the temperature range 200 K ≤ T ≤ 290 K considered at each ϵ_S value. As will be shown in Section III, this temperature range covers the states that are both above and below freezing.

All Molecular Dynamics (MD) simulations are performed using LAMMPS⁶⁰ in the isochoric (NVT) ensemble. We integrate Newton’s equations of motion using the velocity Verlet algorithm⁶¹ with a time step of 2 fs, and we control temperature using the Nosé-Hoover thermostat^62,63 with a time constant of τ = 0.2 ps. All simulations are carried out in cuboidal boxes that are periodic in all dimensions. The simulation box is always longer along the z direction in order to assure the lack of correlation between the liquid film and its periodic images.

In all simulations, the starting configuration is comprised of an n-octane film in which all molecules are arranged into a simple cubic lattice. This configuration is initially heated at T = 300 K for 200 ps in order to melt the crystal. The melted film is then gradually quenched to T_f at a cooling rate of 2.5 × 10¹² K/s. The arising configuration is equilibrated at T_f for 4 ns. This equilibration time is far larger than the structural relaxation time in the bulk for the temperatures considered in this work. Throughout the entire process, a separate thermostat is applied to the atoms in the substrate, always maintained at a temperature T_f. After this initial equilibration stage, the production runs are carried out for 80 ns.

The bulk of the methodology that is used for computing the spatial profiles of thermodynamic and kinetic properties is discussed in our earlier publication.⁵ For thermodynamic quantities that are time-invariant, spatial profiles can be rigorously determined via simple time averaging of the corresponding property in thin cuboidal slices of the simulation box. Here, profiles of potential energy, atomic and molecular density, lateral and normal stress and lateral radial distribution function are computed in slabs that are 0.025 nm thick, and the binning of all molecular properties, including molecular density, and radial distribution function, is performed based on the centers of mass of the molecules. We also compute inherent structures using the FIRE algorithm,⁶⁴ and the contribution to the average inherent structure potential energy of a slice is from those atoms that were originally located in that slice prior to energy minimization. We also compute the orientational distribution function (ODF), f(θ, z), of n-octane molecules defined as follows. First, 𝒢_i the gyration tensor of each n-octane molecule is computed from

Here, m_i,j and r_i,j are the mass and the position of the jth united atom site of the ith molecule and r_i,CM is the center of mass of the ith molecule. Subsequently, v_i, the longest principal axis of the ith molecule is determined, which is the unit eigenvector corresponding to the largest eigenvalue of 𝒢_i. Due to the inversion symmetry of 𝒢_i, both ±v_i will be equally valid choices. In order to avoid any ambiguity, we choose the v_i for which v_i ⋅ n_S ≥ 0. Here, n_S is a unit vector perpendicular to the substrate, and pointing towards the liquid. The orientational distribution function is then defined as

Note that $∫ 0 π / 2 f ( θ , z ) sinθdθ=1$⁠. f(θ, z) is utilized for computing further orientational order parameters that are introduced and discussed in Section III A.

The time averaging procedure that is used for computing profiles of thermodynamic properties cannot be utilized for kinetic properties such as relaxation times as the latter are obtained from ensemble averages of autocorrelation functions. This difficulty stems from the ambiguity of defining autocorrelation functions in open systems. In our earlier publication,⁵ we present a thorough discussion of different heuristics that can be used for defining autocorrelation functions in confined systems. Here, we take the same convention as used in that work, and define atomic and molecular translational and rotational relaxation times as follows. For translation relaxation, we compute the z-dependent self-intermediate scattering function

Here, X = a, m, with a and m corresponding to atomic and molecular self-intermediate scattering functions, respectively. $Δ r i ‖ ( t )$ is the lateral displacement of entity i over time t. 𝔇(z₁, z₂, z) = δ(z₁ − z) δ(z₂ − z) is an indicator that assures that entity i is present at z both in the beginning and at the end of the time window, and $N X ( z , t ) =〈 ∑ i = 1 N X D ( z i ( t ) , z i ( 0 ) , z ) 〉$ is the average number of entities that contribute to the sum for a particular slice in Eq. (5). Similar to thermodynamic properties, contributions to molecular autocorrelation functions are based on the centers of mass of n-octane molecules. For octane films, we use a q = 16 nm⁻¹, which corresponds to the first peak of S(q), the structure factor, in the bulk n-octane liquid computed in an NpT simulation at T = 300 K and p = 0 bar. We observe no significant change in the locus of the maximum of S(q) in bulk octane with temperature. Relaxation times are determined from F_S,X(q, z, t = τ) = 0.2. The atomic and molecular translational relaxation times are denoted by τ_tr,a and τ_tr,m, respectively.

In order to quantify molecular rotational relaxation time profiles, we compute the following z-dependent orientational auto-correlation function:

with P₂(cosθ) = (3/2)cos²θ − (1/2) the second Legendre polynomial. Note that P₂(⋅) is the natural choice considering the inversion symmetry of the gyration tensor used for determining v_i. Analogously, the relaxation of orientations of individual bonds is characterized using h_b(z, t), the bond autocorrelation function, defined as follows:

Here, b_i is the unit vector in the direction of bond i that connects two united-atom sites on an n-octane molecule, and z_i is the z coordinate of the center of the bond. Since a bond is essentially a directed entity with no inversion symmetry, we use P₁(cosθ) = cosθ, the first Legendre polynomial, to quantify its rotational relaxation, as using P₂(⋅) will lead to loss of relevant orientational information. Similar to translation relaxation times, molecular and bond rotational relaxation times are determined from h_m(z, t = τ_rot,m) = 0.2 and h_b(z, t = τ_rot,b) = 0.2, respectively.

Since it is more difficult to obtain suitable statistics for the z-dependent autocorrelation functions introduced above, we use slices that are 0.1 nm thick, four times thicker than the slices used for computing z-dependent time averages of thermodynamic properties.

As mentioned in Section I, sufficiently long n-alkanes can undergo a process known as surface freezing.^52,53 To be more precise, alkane films behave differently at temperatures above and below $T ˜ m$⁠, their effective melting temperature. For films that are sufficiently thick, $T ˜ m ≈ T m$ the equilibrium bulk melting temperature. At $T≤ T ˜ m$⁠, the entire film freezes into a lamellar crystal. At $T ˜ m ≤T≤ T s$⁠, however, only a frozen monolayer emerges at the vapor-liquid interface. Here, T_s stands for the surface freezing temperature. Experimentally, surface freezing occurs for 15 ≤ n ≤ 50.⁶⁵ Yet, we observe surface freezing in n-octane films simulated using the NERD force-field, even though n = 8 for n-octane. For the films considered in this work, complete freezing occurs for T ≤ 230 K irrespective of the ϵ_S value (Fig. 2(a)). The existence of long-range lateral order in these “cold” films is clearly visible in the lateral radial distribution functions depicted in Fig. 2(a). For 235 K ≤ T ≤ 250 K, however, only a frozen monolayer emerges at the vapor-liquid and solid-liquid interfaces, with the center of the film remaining amorphous (Figs. 2(b) and 2(c)). Here, we define freezing in a purely structural sense, i.e., based on the presence of long-range lateral translational order. As will be further elucidated in Section III C, such frozen regions also correspond to dynamically decelerated regions in which translational and rotational degrees of freedom are decoupled.

The crossover temperature of $T ˜ m =232.5±2.5$ K is reasonably close to 216 K, the experimental melting temperature of n-octane.⁶⁶ (No computational estimate of T_m is available for the NERD force-field.) Our findings are also qualitatively consistent with earlier computational studies of n-octane,^67,68 and other closely related n-alkanes⁶⁹ using different force-fields, all reporting surface freezing.

In order to examine the microstructure of frozen monolayers, we compute the ODFs defined in Eq. (4). Fig. 3 depicts the ODFs for films at T = 240 K. The frozen monolayers emerging at vapor-liquid interfaces are characterized by a single peak in f(θ, z), corresponding to a perpendicular arrangement of octane molecules. Such an arrangement is easily visible in the films depicted in Figs. 2(a)-2(c). In the vicinity of the substrate, however, the microstructure of the frozen monolayer depends on ϵ_S. For loose substrates, i.e., ϵ_S = 0.5 (Fig. 3(a)) and ϵ_S = 1 (Fig. 3(b)), the microstructure closely resembles that of the free interface, as evident in the single peak of f(θ, z) at θ = 0 and z ≈ 1 nm. In contrast, a sticky substrate induces a different type of ordering, with the arising ODF having two distinct strong peaks at θ = 0: one at z = 0.65 nm, and one at z = 0.9 nm. The loci of these peaks are identical to the peaks in the molecular density profile depicted in the rightmost panel of Fig. 6. The existence of two peaks in f(θ, z) is a consequence of the corrugated surface of the 001 facet of fcc substrate. In other words, the molecules at the frozen monolayer can be present both in the valleys and peaks of the rough 001 surface. This changes their z value even though they all have the same orientation with θ = 0.

At higher temperatures, no frozen monolayer emerges at vapor-liquid interfaces (Fig. 4(b)). However, n-octane molecules tend to have a mild preponderance to align along the z axis, as evident in the ODF depicted in Fig. 4(a). Such a propensity can be qualitatively described as “pre-freezing,” a phenomenon previously observed in earlier computational studies of alkane films.⁶⁸ The notion of pre-freezing here is different from and should not be confused with the pre-freezing that is occasionally used for describing the interfacial freezing that precedes bulk freezing at T > T_m.⁷⁰ We can quantify the extent for such pre-freezing with the following two order parameters. The first one is $ξ= max z ≥ 6 nm f ( θ , z )$⁠, the maximum of the free-surface peak of ODF. Note that ξ = 1 for a fully isotropic film. Any deviation from unity will henceforth correspond to broken rotational symmetry. As depicted in Fig. 4(c), ξ is always significantly larger than unity, even for the films at T = 290 K. Also the peak always occurs at θ^∗ = argmax_z≥6nmf(θ, z) = 0  irrespective of T. As expected, pre-freezing becomes stronger at lower temperatures.

The second order parameter that is adopted from the liquid crystal literature is called the nematic order parameter (OP) and is essentially the second moment of the ODF,⁷¹ 

For a fully isotropic fluid, S_zz = 0, while the values of +1 and −0.5 correspond to perfect alignment along and perpendicular to the z axis, respectively. Fig. 5 depicts S_zz(z) profiles for different films. Note that the nematic OP is very close to unity for surface-frozen monolayers at low temperatures. Surface pre-freezing of high-temperature films is manifest in positive peaks of the nematic OP corresponding to weak alignment along the z direction. Similar to ξ that decreases with T, the heights of these peaks diminish as T increases.

Similar to the free interface, a mildly pre-frozen monolayer emerges in the vicinity of loose substrates (ϵ_S = 0.5, 1.0), with characteristic mild peaks in the nematic OP profiles (Fig. 5). For sticky substrates, however, a frozen monolayer emerges at all temperatures, with its ODF resembling the one depicted in Fig. 3(c). The nematic order parameter also demonstrates two strong peaks close to unity, consistent with perfect alignment along the z direction spatially modulated by the corrugated surface of the substrate.

Profiles of atomic and molecular densities are depicted in Fig. 6. Both profiles are peaked at the frozen monolayers that emerge at the vapor-liquid and solid-liquid interfaces. In the case of free interfaces and loosely attractive substrates, a single large peak is observed. For sticky substrates, however, the peak splits into two as a result of spatial modulation induced by the corrugated substrate. The loci of these peaks correspond to the maxima of ODF (Fig. 4) and nematic OP (Fig. 5). Such large peaks are always followed by distinct valleys that correspond to the crystal/liquid interface. As evident in Fig. 5, the nematic OP is negative at these valleys as the octane molecules have a (weak) propensity to align parallel to the substrate.

A very interesting feature of octane films is the stratification of octane molecules at free interfaces even in the absence of surface freezing, i.e., at T > 250 K. This behavior is contrary to what is observed in atomic thin films in which density drops monotonically across a free interface.^5,48 A single mild peak emerges in the atomic and molecular density profile, with its amplitude closely correlating with the pre-freezing order parameters defined in Section III A. A linear correlation between ξ (Fig. 4(c)) and ρ_max,free, the amplitude of the molecular density peak, is depicted in Fig. 7. Henceforth, the nonmonotonic behavior of density is a result of pre-freezing that creates local high-density regions in the interfacial region. Density oscillations at a free interface have been previously observed in experimental and computational studies of not only alkanes^68,72 but also alkali metals^73,74 and ionic liquids.⁷⁵ The emergence of density oscillations has been attributed to a wide range of factors, including large separation between the critical and triple-point temperatures⁷⁶ and building block anisotropy.⁷⁷ Unlike the atomic thin films, the oscillatory density profiles in octane films cannot be fitted to the commonly used hyperbolic tangent functional form,⁷⁸ henceforth making the determination of the width of the free interface based on density nontrivial.

Close to the substrates, both atomic and molecular density profiles are oscillatory. This is consistent with the traditional picture of confinement in which a wall induces structure in the surrounding liquid.^79–81 The stratification of the octane liquid is, however, an interplay between the structuring induced by the substrate and the natural propensity of the interfacial region to pre-freeze. Density oscillations in the vicinity of a loosely attractive substrate (ϵ_S = 0.5) are very similar to the profiles at the free interface, with the substrate inducing very little structure in the liquid. This is very similar to theoretical predictions obtained from density functional theory for associating fluids near a hard wall.⁸¹ As ϵ_S increases, however, stratification becomes more pronounced. For instance, multiple peaks in density emerge at ϵ_S = 1.0, even though the solid-liquid interface has the same orientational fingerprints as the vapor-liquid interface.

Profiles of 〈U〉, the average potential energy, are depicted in Fig. 8. 〈U〉 exhibits a maximum at the center of a frozen monolayer. This apparent energetic penalty can be attributed to the lamellar microstructure of a monolayer in which the united-atom CH₂ sites cluster in the middle. Since |ϵ_CH2,CH2| < |ϵ_CH2,CH3| < |ϵ_CH3,CH3|, such CH₂-rich regions will correspond to higher potential energies. This energetic penalty is, however, offset by the CH₃-rich regions at two ends of a monolayer, corresponding to minima in 〈U〉. It is also noteworthy that z^∗, the locus of the maximum of 〈U〉, matches the corresponding z^∗ for density, f(θ, z) and S_zz(z), further confirming that all these oscillations are manifestations of the same phenomenon, i.e., interfacial freezing. Similar to density, 〈U〉 profiles are oscillatory close to the substrate, even in the absence of interfacial freezing, with the amplitudes and number of oscillations increasing with 1/T and ϵ_S. However, we do not observe any nonmonotonicity in the pre-frozen free interfaces, suggesting that the propensity of a surface to pre-freeze has a purely entropic origin.

Although potential energy profiles provide valuable information about the relative energetic stability of different regions of a confined material, they convey very little about the ability of a material to explore its confinement-induced potential energy landscape. A systematic way of characterizing the latter is to obtain an ensemble of inherent structures, and to compute the amount of decrease in potential energy for different regions of the material.⁴⁸ Here, z binning is performed based on the initial (pre-minimization) positions of the individual atoms. Spatial profiles of this quantity, which we denoted as dive profiles in our earlier publication,⁵ are depicted in Fig. 9. For all the films, whether they go through surface freezing or pre-freezing, no noticeable overall difference is detected between the landscape depth in the bulk and at the free interface. The ability of a free interface to modify the potential energy landscape accessible to a liquid is, therefore, far more limited in the case of n-octane. This is unlike the behavior observed in the atomic thin films considered in Refs. 5 and 48 in which a vapor-liquid interface increases the depth of the potential energy landscape accessible to the liquid at its vicinity. Such a heightened access to the minima of the potential energy landscape appears to be key in the ability of a material to form ultrastable glasses upon vapor deposition.⁸² By this token, the ability of n-alkanes to form ultrastable glasses will be limited considering the lack of elevated access to the potential energy landscape of the liquid at the free interface.

Of all the thermodynamic quantities considered in this work, none oscillates more strongly than normal and lateral stress (Fig. 10). Such oscillations are particularly strong in frozen monolayers at both interfaces, and extend well beyond the frozen region. Indeed, sharp interfaces are maintained as a result of surface tension that stems directly from anisotropies of the stress tensor.⁸³ In the case of frozen monolayers, in particular, two such sharp interfaces are present: one between the substrate (or vapor) and the frozen monolayer, and the other between the monolayer and the liquid. Consequently, oscillations in stress occur across a relatively wider region of the film.

It can be argued that stress anisotropy is among the most systematic ways of defining what constitutes an interfacial region. Indeed, a region with anisotropic stress responds differently to mechanical stress, and will thus have mechanical properties distinct from the bulk. In this work, we define the bulk region as the largest connected region in which the difference between normal and lateral stress does not exceed 20 bars. Choosing this cutoff, which is always less than 10 per cent of the largest difference between lateral and normal stress at T = 290 K, is necessary considering the difficulty of converging virial-based estimates of stress in molecular simulations. Applying this criterion yields an interfacial width that increases with temperature, from ≈2.30 nm at 260 K to ≈2.45 nm at 290 K. The interfacial width, as determined from the anisotropy in stress tensor, is very difficult to measure in experiments. Nevertheless, the temperature scaling of the dynamical length scales in thin films of polystyrene glasses reveals a similar dependence of temperature, with the width of the highly mobile region increasing with temperature.⁸⁴ It can indeed be argued that a correlation must exist between the thermodynamic length scale (obtained from stress tensor) and the dynamical length scales (defined as per mobility).

In our earlier work,⁵ we established a strong correlation between the substrate-induced profiles of normal density and normal stress in amorphous regions of atomic thin films. By inspecting the density and normal stress profiles, we explore the existence of such a correlation in octane films. Here, the quantity of relevance is the atomic—and not the molecular—density, as the positions of individual molecules are not sufficient for determining the stress tensor. We observe a correlation between atomic density and normal stress in amorphous regions of octane films. For loosely attractive substrates (Fig. 11(a)), the correlation is weaker to the extent that the first maximum in normal stress only corresponds to a shoulder in atomic density. In other words, the solid-liquid interfacial region is very similar to a free interface, as the substrate exerts minimal effect in ordering and stratifying the liquid in its vicinity. In the case of stickier substrates (Fig. 11(b)), however, the correlation is much more pronounced, and the peaks and valleys of atomic density and normal stress follow one another more closely. It is, of course, necessary to emphasize the absence of such correlations in parts of the film that have long-range translational order, as crystals can respond to stress in nontrivial ways.

Figs. 12 and 13 depict profiles of translational and rotational relaxation times. In surface-frozen monolayers, a sharp increase in relaxation times is observed, with both the atomic and molecular relaxation times increasing by at least one order of magnitude with respect to their bulk values. We observe an even larger increase in molecular and bond orientational relaxation times. We are unable to quantify the extent of such an increase as the orientational autocorrelation functions given in Eqs. (6) and (7) never relax to zero in the time scale of our simulations. This decoupling of translational and orientational degrees of freedom is a direct consequence of layering at the surface, and is common in phases with long-range orientational order such as liquid crystals.⁸⁵ In such systems, the relaxation of translational degrees of freedom proceeds through a slipping mechanism that does not require orientational relaxation.

In Section III B, we discuss non-monotonicities in several thermodynamic quantities, including ξ, S_zz, and density, across free interfaces of the films that do not undergo surface freezing. We do not detect any such non-monotonicity in relaxation time profiles. Yet, we are unable to rule out its possibility for the following reasons. First of all, the uncertainties associated with computing dynamical quantities such as relaxation times are generally larger than the uncertainties in estimating thermodynamic quantities such as density. After all, the former are computed from auto-correlation functions while the latter are obtained from simple time averaging. Second, the convention used for defining z-dependent autocorrelation functions will inevitably lead to some mixing between neighboring slices as the molecules that contribute to the autocorrelation function of a particular z slice can leave and re-enter that slice. Such a procedure can mask the existence of subtle spatial differences in relaxation time, a situation that is completely plausible in the pre-frozen interfaces considered in this work.

The dynamical features of a film in the vicinity of a substrate depend on ϵ_S. Similar to the free interface, the surface-frozen monolayers close to loose substrates demonstrate a deceleration of dynamics and a decoupling of translational and rotational degrees of freedom. This slow-down happens in the vicinity of a sticky substrate as well. However, it is not clear whether translational and rotational degrees of freedom decouple as the corresponding relaxation times cannot be computed in the time scale of our simulations. For the films that do not undergo surface freezing, however, two distinct dynamical regimes are observed. For sticky walls (ϵ_S = 3.0), the dynamics is decelerated close to the substrate, even in amorphous regions of the film. Indeed, relaxation times can be as much as two orders of magnitude larger than the corresponding bulk values. In the vicinity of loose substrates (ϵ_S = 0.5), however, the dynamics is accelerated, as in atomic thin films. For ϵ_S = 1.0, the relaxation time profiles are more oscillatory with the overall dynamics becoming only slightly faster with respect to the bulk. This is consistent with the behavior observed in atomic thin films that also demonstrate two distinct dynamical regimes in the vicinity of sticky and loose substrates.

Fig. 14 depicts the temperature dependence of bulk relaxation times for ϵ_S = 0.5. Those are obtained from averaging the relaxation time profiles of Figs. 12 and 13 over z[nm] ∈ [2.5, 7.5]. An identical behavior is observed for other ϵ_S values, with the results not shown for conciseness. All relaxation times demonstrate an Arrhenius-type dependence on temperature, with no sign of fragility. This is not unusual considering the earlier viscosity measurements for n-octane,^86,87 which predict a glass transition temperature of T_g = 85 K, markedly lower than the temperatures considered here. Note that atomic translational relaxation times are always lower than their molecular counterparts, due to the relative ease of relaxing atomic degrees of freedom. It is, however, not possible to systematically compare the bond and molecular orientational relaxation times, since they correspond to relaxation features of different Legendre polynomials. Therefore, the fact that the bond relaxation times are higher does not have any physical meaning. Similarly, it is not meaningful to compare translational and orientational relaxation times, as they have also been computed from decays of different autocorrelation functions.

In our earlier work,⁵ we established a close correlation between lateral relaxation time and lateral stress profiles in amorphous regions of atomic thin films. This is consistent with what happens in simple liquids in which diffusivity increases upon decreasing pressure.⁸⁸ We find a similar correlation between the profiles of atomic translational relaxation time and lateral stress in octane films. Fig. 15 shows such a correlation for two representative films, with the peaks and valleys of τ_tr,a and lateral stress closely following one another in the amorphous region (determined from RDF). It is interesting to note that such a correlation can also be rationalized by the fact that the compressibility of a simple fluid increases upon decreasing pressure, and this enhances structural relaxation by facilitating local density fluctuations.

In this work, we study the effect of substrate on thermodynamic and kinetic anisotropies in n-octane thin films. We observe complete freezing at temperatures below 232.5 ± 2.5 K. For 235 K ≤ T ≤ 250 K, however, a frozen monolayer emerges at vapor-liquid interfaces, as well as in the vicinity of loosely attractive substrates. Such frozen monolayers correspond to large peaks in several thermodynamic and kinetic properties, such as density, potential energy, nematic order parameter, and translational and rotational relaxation times. Also, translational and rotational degrees of freedom are decoupled in the monolayer due to long-range orientational ordering of the n-octane molecules.

At higher temperatures, interfacial freezing only occurs in the vicinity of sticky substrates. At vapor-liquid interfaces, and close to loose substrates, only a weak propensity is observed for the molecules to “pre-freeze” and align perpendicular to the interface. Such a propensity manifests itself in mild peaks in atomic and molecular density and f(θ, z) and S_zz orientational order parameters.

The amorphous regions of the films are stratified in the vicinity of substrates, with the extent of stratification increasing upon decreasing temperature, or increasing the interaction parameter ϵ_S. Density oscillations are minimal close to loose substrates, and the interface almost resembles a vapor-liquid interface. Similar to density, oscillations are observed in other thermodynamic quantities such as potential energy, dive (as a result of energy minimization), stress, and orientational order parameters. In the vicinity of a substrate, we confirm the existence of correlation between atomic density and normal stress in amorphous regions of the film, consistent with our observations in atomic thin films. We also observe two distinct dynamical regimes at the solid-liquid interface. In the vicinity of loose substrates, dynamics is accelerated, while a sticky substrate decelerates dynamics in its vicinity. Finally, we are able to establish a correlation between lateral atomic translational relaxation times and lateral stress in amorphous regions of the films.

Short-chain n-alkanes are known to be poor glass formers due to their strong propensity to crystallize.⁸⁹ Our findings suggest that the vapor deposition process that yields ultrastable glasses for many materials is unlikely to be very successful in the case of n-alkanes. This is not only due to the propensity of alkane films to undergo surface freezing, but also because of the fact that the free interface does not create a more accessible potential energy landscape even when the alkane film is amorphous. It has, indeed, been recently demonstrated that there are materials that cannot form an ultrastable glass upon vapor deposition,⁹⁰ unlike earlier suggestions that this process might be universal.⁹¹ This has led to speculations that a correlation might exist between the fragility of a liquid and its ability to form an ultrastable glass.⁹² This correlation is not clear-cut as some relatively strong liquids still form ultrastable glasses upon vapor deposition.⁹³ There might, however, be other dimensions to this problem, namely, the presence of a mobile interfacial region, as well the accessibility of deeper minima of the potential energy landscape at the interface. If this picture is accurate, n-octane must not yield ultrastable glasses upon physical vapor deposition. It is therefore a worthwhile experimental endeavor to investigate vapor-deposited glasses of n-alkanes, and other molecules with aliphatic chains to determine whether this picture is accurate.

It is necessary to emphasize that, according to this picture, elevated orientational ordering at a free interface will not necessarily make a material a poor ultrastable glass former. Indeed, experimental^30–32 and computational⁵⁰ studies of ultrastable glasses have confirmed the existence of such ordering at vapor-liquid interfaces of good ultrastable glass formers.

It is necessary to emphasize that the microstructure of the liquid in the vicinity of the solid-liquid interface might depend on the structure of the substrate. For instance, the orientational order that would emerge close to a sticky structureless substrate will be different from what is observed here, as reported in Ref. 94. Since the liquid that is close to a low-ϵ_S substrate is structurally very similar to the liquid at the free interface, we expect that such differences will only be important at large values of ϵ_S. Even then, the observed deceleration of the dynamics (with respect to bulk) is not expected to disappear if a structureless sticky substrate is utilized. It is, however, interesting to study the combined effect of ϵ_S and substrate corrugation on the microstructure of the frozen monolayers, as well as the thermodynamic and kinetic anisotropies in the film.

There are other interesting questions to be asked about alkane films besides their ability or lack thereof to form ultrastable glasses. One such question is the role of substrate in inducing—or suppressing—surface freezing. This can be systematically addressed by performing rigorous free energy and/or rate calculations. Finally, the ultra-thin films of long-chain alkanes can be studied to investigate the possible interplay between solid-liquid and vapor-liquid interfaces. All these topics can be the subject of further studies.

We acknowledge the support from the National Science Foundation, NSF Grant No. CHE-123343 (to P.G.D.) and the Princeton Center for Complex Materials, a MRSEC supported by NSF Grant No. DMR-1420541. These calculations were performed on the Terascale Infrastructure for Groundbreaking Research in Engineering and Science (TIGRESS) at Princeton University. We gratefully acknowledge A. Panagiotopoulos, M. Ediger, D. Sussmann, M. Tylinski, Y. E. Altabet, and R. Singh for useful discussions.