Orientational order induced mode switching at coupled interfaces of a nematic–isotropic free bilayer

Abstract

Instabilities at the coupled deformable interfaces of a nematic–isotropic free bilayer can engender rich morphological patterns governed by an intricate interplay of nematic elasticity, surface and interfacial energetics, gravitational forces, and intermolecular forces. This work presents a theoretical investigation of the instability in a bilayer composed of a liquid crystal film resting on a water layer having a free deformable nematic–air surface and a confined nematic–isotropic interface. Utilizing a long-wave hydrodynamic model, the formulation couples the Ericksen–Leslie equations with the Navier–Stokes equations to quantify growth rates of the instability across NLC thickness and director orientations. The formulation examines the competing roles of gravitational, van der Waals, and elastic forces in governing the deformations at the bilayer interfaces. For micro-thick films, density contrast between the layers initiates classical Rayleigh–Taylor instabilities (RTIs), while in nano-thin regimes, disjoining pressure drives mode dynamics. In both scenarios, we identify a pair of primary instability modes—the confined interfacial mode (CIM) and the free surface mode (FSM)—whose relative dominance is dictated by film thickness, anchoring configuration, and the spreading coefficients at the surface or interface. Systematic variation of director orientations and three canonical surface anchoring arrangements reveals critical transitions and instability mode responses, especially under asymmetrical spreading conditions. The results uncover mode switching controlled by NLC orientation and interfacial energetics, offering fundamental insights into the spatiotemporal structuring of layered soft materials. The framework provides an elementary direction for designing reconfigurable morphologies in NLC-based bilayer systems relevant to wetting, templating, and self-patterning applications.

---

1. Introduction

The interfacial dynamics of unstable thin liquid films have inspired research enthusiasts for ages because of their prominence in a wide range of applications, which include paints,^1,2 coatings,^2,3 flotation,⁴ artificial skins,⁵ mucous linings,⁶ packaging,⁷ emulsification,⁸ and phase transition,⁹ among many others. It is now well explored that the physics behind the interfacial instability of a thick film is fundamentally different from that of an ultra-thin (<100 nm) film.¹⁰ For example, while an ultra-thin liquid film can self-organize an array of correlated holes when the van der Waals interactions overpower the stabilizing surface tension forces,^10–12 a relatively thicker film can dewet a surface following the defect-induced heterogeneous nucleation pathway.^12,13 In such a scenario, in the late stages of the dewetting process, the competition between the stabilizing and destabilizing components of the surface tension forces triggers Rayleigh–Plateau instability^14,15 to facilitate the formation of the droplets. Although the breaking of thin films is an unwanted phenomenon for the paints or coatings, the self-organized dewetted patterns with a large-area ordering¹⁶ find important applications in modern display units,¹⁷ microreactors,¹⁸ optofluidic devices,¹⁹ photodetectors,²⁰ photovoltaic cells,²¹ and soft lithography.²²

Interestingly, when another immiscible layer confines a liquid film, the interface may also develop interesting morphologies under varied conditions based on the contrast of body forces – Rayleigh–Taylor instability (RTI),^23,24 viscous forces – Saffman–Taylor instability (STI),²⁵ surface tension forces – Marangoni (MI),^26,27 or inertial forces – Kelvin–Helmholtz (KHI)²⁸ instability. Importantly, such systems can largely be classified into two types: (i) a confined bilayer – a pair of immiscible fluids with a single confined interface²⁹ and (ii) a free bilayer – a pair of films with confined and free interfaces.^30,31 In particular, the dynamics of the free bilayers are far more complex owing to the presence of the twin liquid–liquid and liquid–air interfaces, which help develop fascinating flow and droplet morphologies under varied conditions.³² For example, when the films are thicker and a denser fluid is dispensed on a rarer fluid, the destabilizing gravitational field overpowers the stabilizing surface tension to manifest an RTI.^24,33,34 On the other hand, when both the films are ultra-thin, a different type of RTI may be expected wherein the destabilizing van der Waals forces overpower the stabilizing surface tension forces to engender micro- or nanoscale patterns.^12,30 Furthermore, the kinetic factors, such as the contrast of the viscous, elastic, or viscoelastic forces across the interface, play key roles in deciding the morphologies of the patterns decorated.^31,35 Previous studies suggested that using external fields in such systems can significantly influence the length and time scales of such instabilities.^36–39

The spontaneous rupture and pattern formation in thin polymer films over thick liquids have been widely studied experimentally and theoretically over the past few decades. The seminal experiments conducted by Lambooy et al.⁴⁰ and Sferrazza et al.⁴¹ demonstrated dewetting and instability in thin polymer bilayers, revealing how van der Waals and capillary forces govern interfacial modulation and rupture kinetics. More recently, Das et al.⁴² and Bhandaru and Mukherjee⁴³ have reported ordering transitions and regular pattern formation in polymer bilayers during dewetting, showing that the interplay of viscosity contrast, interfacial tension, and confinement leads to tunable morphologies. These studies collectively underline that when two immiscible polymeric layers are coupled, instabilities arise through competing capillary and intermolecular interactions, resulting in coupled-mode deformations. Of late, investigations have highlighted ordering dynamics in polymer bilayers and composite films.^29,44 These works emphasize how hydrodynamic coupling and substrate interactions can dictate the wavelength and symmetry of dewetted patterns, offering parallels to the coupled-interface phenomena explored in the present work.

Importantly, the free surface instabilities of the isotropic viscous, elastic, or viscoelastic films are vastly different from those of the anisotropic liquid crystal (LC) films. The inherent orientational order of LCs quantified by a ‘director’ field^45,46 adds an additional handle to tune the length and time scales of such instabilities. Previous studies showed that the primary focus in this regard has been on studying the phase transition of ultrathin LC films^47–49 or the contact line instabilities of LC droplets.^50,51 The literature indicates that the elastic rigidity of the LCs originating from their director field may help control the deformation dynamics of the interfaces in a free bilayer.^45,52 In particular, if a nematic liquid crystal (NLC) film is placed on a liquid layer to form a free bilayer, the orientation of the LC molecules at the boundaries has the capacity to control the dynamics of both the NLC–air and NLC–liquid interfaces to develop interesting flow morphologies.^53–55

On the theoretical front, a number of pioneering articles^56,57 published previously have explored theoretically the free surface dynamics of the NLC droplets while considering an additional stress term originating from Frank's bulk elasticity. In this regard, different approaches^58,59 unveil the influence of the weak to strong elastic effects of the NLC droplets or films alongside considering the effects of anchoring conditions of NLC directors at the boundaries, e.g., homeotropic, planar, and angular.^45,46 In this regard, the electrohydrodynamic convection of nematic films has been studied using the Ericksen–Leslie governing equation,^60–62 considering the Frank–Oseen elastic energy.^63,64 The deformation of the NI interface has also been studied using uniform and non-uniform AC electric fields⁶⁵ and magnetic fields⁶⁶ employing long-wave linear stability analysis. Recently, general linear stability analysis (GLSA)⁶⁷ involving the Ericksen–Leslie governing equations coupled with the anisotropic Maxwell stresses for a free NLC film showed that the air-to-NLC filling ratio between the electrodes, dielectric anisotropy of the film, Frank elasticity, director orientations, Leslie coefficients, and boundary conditions play pivotal roles in selecting instability modes.^60,62–64

In view of this background, the instabilities of a free bilayer have been explored, which is composed of an NLC upper layer and a water layer resting on a rigid surface, as schematically shown in Fig. 1. The system is composed of a free NI surface and confined nematic–water (NW) interface, wherein the NLC molecules are capable of anchoring at both these places. For the NLC layer, the Ericksen–Leslie governing equations, coupled with Frank bulk elasticity, have been considered to capture the dynamics, while Navier–Stokes equations have been used to model the water layer. The set of equations with the appropriate boundary conditions are linearized (LWLSA) following the lubrication approximations to estimate these instabilities' time and length scales. A mathematical formulation has been developed in such a manner that the bilayer composed of micro-thick (μm thickness) and nano-thin (≤100 nm) films can be analyzed simultaneously in the asymptotic limits.

Fig. 1 Schematic representations of (a) an NLC–water bilayer system bound by the air and (b) director orientation configurations. In the schematic (a) h₁(x, t) and h₂(x, t) are the local thicknesses with h₂₀ and h₁₀ as the mean thicknesses of NLC and water films, respectively. ρ₂ and ρ₁ are the densities of NLC and water films, and s stands for the substrate location, respectively. The θ is the angle that the director (n) makes with the positive z-axis direction as displayed on the right side zoomed circle in the schematic (b) and has components n_x = sin θ and n_z = cos θ. The schematic (b) shows the different director orientation configuration combinations used in the present theoretical analysis, wherein θ₁ and θ₂ are the director angles at the NW and NI interfaces, respectively. The designated abbreviations in b(i)–(v) stand for homeotropic–planar (H–P, θ₁ = 0, θ₂ = π/2), homeotropic–homeotropic (H–H, θ₁ = 0, θ₂ = 0), planar–homeotropic (P–H, θ₁ = π/2, θ₂ = 0), planar–planar (P–P, θ₁ = π/2, θ₂ = π/2), and planar–angular (P–A, θ₁ = π/2, θ₂ = π/4), respectively. Here, g is the gravity acting in the negative z-direction.

For the micro-thick films, a denser NLC layer on a rarer water film is expected to stimulate RTI due to the density difference, whereas for the nano-thin films, the long-range van der Waals interactions facilitate the instabilities suppressing the classical RTI due to the density difference.^68–70 In both situations, the surface and interfacial tensions at the interfaces, the NLC film's nematic elasticity, and the NLC molecules' anchoring at the surface/interfaces may impart a stabilizing influence to the systems. The results reported in this work uncover the possible length and time scales of the interfacial instabilities recently encountered in an NLC–water system.^53,71 Furthermore, this work also reports the possibilities in the case of similar systems at the nanoscale.

2. Mathematical problem formulation

Fig. 1 schematically shows a bilayer comprised of water (lower, 1) and NLC (upper, 2) films resting on a solid substrate (s), depicting relative deformation at the respective interfaces. The schematic pictorially shows the typical rod-like NLC molecules oriented inside the NLC film while resting on the water layer. The NLC film with long-range orientational order directors is exposed to the bounding medium – an inviscid gas such as air. The long-range orientational order of NLC molecules is described by a director vector n = (n_x, n_z). The formulation utilizes a two-dimensional (2D) Cartesian coordinate system with the origin at the solid substrate interface. Here, x and z represent the coordinates parallel and perpendicular to the substrate, respectively, while the variable t denotes time. Boldface variables signify vectors and tensors. The present formulation assumes isothermal and incompressible films for both layers.

A. NLC description

The director field deformation contribution to the free energy, which, in the context of the NLC, is referred to as the Frank–Oseen elastic energy (E_f) as represented in eqn (1), wherein, K₁, K₂, and K₃ represent the elastic splay, twist, and bend constants. The saddle-splay constant (K₂ + K₄) is later eliminated due to no contribution to governing equations for the case of strong anchoring, as elaborated in the SI. These modifications turn into⁴⁶

E_f = 0.5K₁(∇·n)² + 0.5K₂(n·∇ × n)² + 0.5K₃(n × ∇ × n)² + 0.5(K₂ + K₄)∇·((n·∇)n − (∇·n)n). (1)

The implication of one constant approximation^46,58 simplifies the E_f as follows,

E_f = 0.5K_f((∇·n)² + |∇ × n|²). (2)

where K_f= K₁ = K₂ = K₃ and K₄ = 0 and K_f is the NLC bulk elastic constant.

In order to describe the NLC dynamics, the equation of balance of couple (angular momentum balance) reads as follows, ensuring the director field constraint⁵⁸

n·n = 1, (3)

as,The Lagrange multiplier λ also confirms the director constraint in the couple balance.^58,72 A unit vector describing the director, n, in a 2D can be expressed as n (sin θ, cos θ), where θ represents the director angle with positive z-axis (Fig. 1(b)) that also satisfies eqn (3). The other variables in eqn (4) are – the body force due to any external field (G), an intrinsic director and body force (g). In our formulation, there is no external field present that leads to G = 0, and g reads as,

g = −λ₁N − λ₂(e·n). (5)

The rotation vector N and the strain tensor e in eqn (5) can be expressed as, N = ṅ − ω·n, where the spin tensor, ω = 0.5(∇u − ∇u^T) and e = 0.5(∇u + ∇u^T).

The NLC film's constitutive relation is described using a viscous stress tensor τ^V as,^56,73,74

τ^V = α₁e: nnnn + α₂nN + α₃nN + α₄e + α₅nn·e + α₆e·nn, (6)

where the component α_i (i = 1–6) represents the Leslie viscosity coefficient that later can be related to the rotational viscosity (λ₁) and irrotational torque coefficient (λ₂)^58,74 as λ₁ = α₃ − α₂, λ₂ = α₆ − α₅ and to each other using the Onsager relation as α₂ + α₃ = α₆ − α₅ respectively. However, they do not necessarily have to be all positive, and 0.5α₄ > 0 when other coefficients are zero resembles a typical Newtonian-isotropic system of common viscosity. Furthermore, the elastic (Ericksen) stress tensor, τ^E, is described as,⁵⁸

B. Momentum balance equations

The conservation of mass and momentum equations for the lower layer (1) in the vectorial form can be expressed as,

∇·u₁ = 0, (8)

where u(u, w), p₁ read as (P₁ + ϕ₁), ϕ₁, ρ₁, ∇, I and τ read as (∇u + ∇u^T) are the velocity vector, pressure, body force term, water density, gradient operator, identity tensor, and hydrodynamic stress tensor, respectively. Here, D/Dt = ∂/∂t + u·∇ is the total derivative, and t is the time. Similar common notations are used in the following equations unless mentioned.

Similarly, for the upper NLC layer, the mass and momentum equations read as,

∇·u₂ = 0, (10)

where p₂(P₂ + ϕ₂), ϕ₂, τ^V, and τ^E are pressure, body force potential, viscous, and elastic stress tensors, respectively. The detailed term-by-term simplifications of eqn (11) are mentioned in the SI and elsewhere.^56,58

C. Boundary conditions

In the bilayer, the lower layer is placed on the rigid and impervious substrate – s and the NLC film rests on top, exposed to the air medium, respectively. The no-slip and impermeability are assumed at the substrate at the location z = 0, which reads as,

u₁ = 0. (12)

The system consists of an interface – NW at z = h₁(x, t) and a surface NI at z = h₂(x, t), respectively. At these locations, the continuity of velocities, normal and tangential stress balances, and kinematic conditions are implemented as boundary conditions. Subscripts h₁ and h₂ in eqn (16) and (19) denote locations at which u₁ and u₂ are evaluated. These boundary conditions read as, at z = h₁(x, t),

u₁ = u₂ and w₁ = w₂, (13)

n₁·T₂·n₁ − n₁·T₁·n₁ = γ₁₂κ₁, (14)

n₁·T₂·t₁ − n₁·T₁·t₁ = 0, (15)

∂h₁/∂t + (u₁)_h1∂h₁/∂x = (w₁)_h1, (16)

at z = h₂(x, t),

n₂·T₀·n₂ − n₂·T₂·n₂ = γ₂κ₂, (17)

n₁·T₀·t₁ − n₁·T₂·t₁ = 0, (18)

∂h₂/∂t + (u₂)_h2∂h₂/∂x = (w₂)_h2. (19)

Subscripts 0, 1, and 2 denote air, water, and NLC layers. The components h_i, γ₁₂, γ₂, and κ_i = ∇·n_i (i = 1, 2) represent the respective film thickness, interfacial tension, surface tension, and the respective deforming interface curvature. The notations ∇_s (= I − nn), and represent the surface gradient operator, outward surface normal, and tangent vectors. At the NI interface, we assume that there exists no tangential surface tension gradient, which leads to no jump at NI, that is n₂·T₂·n₂ = −γ₂κ₂ and n₁·T₂·t₁ = 0, since T₀ = 0 (eqn (16) and (17)) for the bounding medium (air), making it appropriate for assumed strong anchoring conditions. Essentially, by maintaining the general form of normal and tangential stresses at the interface, the derived equation can examine different body and interface forces arising from specific influences such as gravity and surface tension and including disjoining pressures within the body force potential term. Furthermore, the strong anchoring boundary conditions are executed to express the director field (n) in order to achieve the planar (P) and homeotropic (H) director positions at NW and NI, respectively. In this study, the LC–water interface is modeled with boundary conditions representing two experimentally relevant anchoring types, as mentioned above, which are commonly observed in NLCs like 5CB. The polarity of water aligns the LC's polar cyano (–CN) groups toward the aqueous phase, promoting homeotropic or tilted orientation, while the hydrophobic alkyl chains orient away.⁷¹ These molecular interactions form the physical basis for the imposed boundary conditions. The enforcement of these boundary conditions at NW (z = h₁(x, t)) and NI (z = h₂(x, t)) leads to,

n·z = cos θ₁, (20)

n·t_i = cos θ₂. (21)

3. Long-wave hydrodynamics

A. Thin film equation

In this section, we show the equations that describe the temporal evolution of the deforming NW interface and NI-free surface. This is achieved in the framework of the lubrication or long-wave approximation (ε = 2πh_i0/λ_w ≪ 1), considering the general properties of the layers involved. In the framework of lubrication approximation, we determined that the x-direction and time scale as O(1/ε), while the z-direction scales as O(1). This considers the wavelength for the x-directional length scale and the mean film thickness (h_i0) for the z-directional length scale. Additionally, the velocities (u, w), total pressure, director orientation (θ),⁶⁷ and surface or interfacial tension (γ) scale as follows – O(1), O(ε), O(1/ε), O(1), and O(1/ε)³. Utilizing this, we simplify all governing equations and boundary conditions. We conduct an order analysis on each equation, retaining only the leading-order terms and disregarding any higher-order terms. The evolution equation is subsequently derived from these simplified equations. The following equations outline all the reduced governing equations obtained from eqn (8)–(11) and boundary conditions from eqn (12)–(19).

∂u₁/∂x + ∂w₁/∂z = 0, (22)

∂u₂/∂x + ∂w₂/∂z = 0, (23)

∂p₁/∂x = μ₁∂²u₁/∂z², (24)

∂p₁/∂z = 0, (25)

∂/∂x(p₂ + 0.5K_f(∂θ/∂z)²) = ∂/∂z(c₂∂u₂/∂z), (26)

∂/∂z(p₂ + 0.5K_f(∂θ/∂z)²) = 0, (27)

−K_f(∂²θ/∂z²) = 0, (28)

P₁ + ϕ₁ − P₂ − ϕ₂ − K_f(∂θ/∂z)² = −γ₁₂(∂²h₁/∂x²); at (z = h₁(x, t)), (29)

c₂∂u₂/∂z = μ₁∂u₁/∂z; at (z = h₁(x, t)), (30)

P₂ + ϕ₂ − p₀ + K_f(∂θ/∂z)² = −γ₂(∂²h₂/∂x²); at (z = h₂(x, t)), (31)

∂u₂/∂z = 0; at (z = h₂(x, t)), (32)

∂h₁/∂t + (u₁)_h1∂h₁/∂x = (w₁)h₁, (33)

∂h₂/∂t + (u₂)_h2∂h₂/∂x = (w₂)_h2. (34)

The term c₂ in eqn (26) and (30) is one of the coefficients (the remaining are given in the SI) that belongs to the hydrodynamic stress, which survives after scaling analysis and is expressed as, c₂ = 0.5(2α₁ sin² θ cos² θ + α₄ + (α₂ + α₅)sin² θ + (α₆ − α₃)cos² θ). Solving the aforementioned governing equations and boundary conditions and substituting the obtained terms in the kinematic condition result in the evolution equations for the respective layers. For the NLC layer, the Leibnitz formula is utilized to transform eqn (19) into a mathematically more convenient form, . The resulting evolution equations, detailed below, describe the stability, dynamics, and morphology of both NI and NW.where, ; ; ; ; and . Here, we call P_2M = P₂ + 0.5K_f(∂θ/∂z)² as a modified pressure, and the NLC bulk viscosity is defined by μ₂ = 0.5α₄, where, α₄ is the fourth Leslie coefficient.

B. Micro-thick and nano-thin cases

(i) Gravity and surface tension. In the present study, the bilayer comprises an NLC film on the top of a water layer, which introduces the density difference, as NLC is physically a higher-density fluid than water. In order to examine the system dynamics and interfacial stability under this condition, we incorporate the inertial effects in the body force terms appearing in the momentum balances and normal stresses for the individual layers as follows.

ϕ₁ = ρ₁gz and ϕ₂ = ρ₂gz, (37)

where g is the gravitational constant, which is in the negative z direction in the present case. The linear stability analysis is performed after solving momentum balance and normal stress simplification for the interface locations utilizing the normal linear modes (h_i = h_i0 + h̃_ie^ωt+ikx), where, h_i0, ω, k and h̃_i are the mean or base thickness, linear growth rate of disturbance, wavenumber, and infinitesimal initial surface perturbation amplitudes. The present longwave formulation utilizes strong anchoring conditions, where surface and interfacial tensions are assumed constant. Therefore, it is essential to understand that constant interfacial tensions correspond to the fixed base-state director angles θ₁ and θ₂. Moreover, if an explicit γ(θ) was included, the surface tension would vary locally with the perturbation amplitude. The resulting dispersion relation (see Bandyopadhyay et al. (2005) for a detailed derivation of the isotropic bilayer)³⁰ from eqn (33) and (34) is expressed as follows:

h̃₁ω + h̃₂A + h̃₁B = 0, (38)

h̃₂ω + h̃₂C + h̃₁D = 0, (39)

with

(ii) Disjoining pressure. The long-range intermolecular forces can significantly influence the nano-thin film stability by affecting the intermolecular interactions at the interfaces.^12,30 These forces are particularly critical in determining the equilibrium configurations and potential dewetting phenomena in thin liquid films. The van der Waals effects are incorporated into the stability analysis (detailed steps in the SI) through an additional pressure term that represents the disjoining pressure following a few previous studies.^47,75–78 The effective Hamaker constants (A_l) are related to the binary Hamaker constants (A_lm) appearing in the ΔG^LW_E as follows,

A_l = A_mm − A_lm = −12πd₀²S_lm; A_m = A_ll + A_sm − A_sl − A_lm = −12πd₀²S_slm;

and

A_n = A_sm − A_lm = −12πd₀²S_slm + 12πd₀²S_sl,

where l, m, and n represent subscripts 1, 2, and 3. The macroscopic equilibrium spreading coefficients utilised in primed coefficients of eqn (46) below, denoted by the respective subscripts where n is a gas medium. The letter ‘a’ is used to represent the bounding medium (air) in Section 4(C) while writing the respective spreading coefficients in place of n just for simplification. Therefore, S_lm = γ_l − γ_m − γ_lm, and d₀(∼0.157 nm) is the van der Waals cut-off distance between the two surfaces, determined by incorporating the short-range Born repulsion.^30,79 The excess free energy (ΔG^LW_E) for the present bilayer configuration can be expressed as,^47,76,78The term h₃, representing the NLC layer thickness, is the same as h₂−h₁ and alternatively used in the rest of the manuscript. Thus, the disjoining pressure (Π = −∂/∂h(ΔG^LW_E)) is negative of the van der Waals conjoining pressure (Φ) and for a particular interface that modifies the pressure distribution within the film and is given by,

Here, is a constant and obtained from the angular momentum balance equation. The total pressures (P₁ = p₀ − γ₂∂²h₂/∂x² − γ₁₂∂²h₁/∂x² − Π₁) and (P₂ = p₀ − K_f(∂θ/∂z)² − γ₂∂²h₂/∂x² − Π₂) at the NW and NI interfaces are double differentiated (details in the SI) with respect to x and then linearized as follows,Eqn (35) and (36) can be modified by substituting eqn (42) and (43) and performing linear stability analysis utilizing normal modes, which leads to the following forms,

h̃₁ω + h̃₂A′ + h̃₁B′ = 0, (44)

h̃₂ω + h̃₂C′ + h̃₁D′ = 0, (45)

where,



The obtained quadratic equations (eqn (40) and (46)) provide the wavenumber spectrum for disturbances where the real part of ω(Re{ω} > 0) is positive, indicating the condition for instability. We focus on the largest positive root of this equation for a given wavenumber, as it represents the maximum growth rate of the disturbance. In the unstable range of wavenumbers, the most rapidly growing or dominant wavenumber (k_m) is determined by the condition ∂ω/∂k = 0. This wavenumber corresponds to the dominant wavelength (λ_m) of the instability, where λ_m = 2π/k_m. The associated growth rate, ω_m, at k_m is the highest initial growth rate for small-amplitude perturbations, which defines the characteristic time scale of the instability.

4. Results and discussion

A. Validation

In this section, the formulations for both micro-thick and nano-thin cases are validated in the asymptotic limits with existing analytical solutions, as illustrated in Fig. 2. Following a previous⁸⁰ formulation, which focuses on the bilayer RTI in isotropic films, the result is reproduced for the dispersion relation as displayed in plot (a) in the creeping flow limit. In order to validate these existing results, the NLC terms are eliminated from the present formulation, and as per their criterion, equal viscosities (μ₁ = μ₂) are maintained with upper film density (in their case – ρ₁) ρ₁ = 1100 kg m⁻³ (NLC – ρ₂ in the present case) and interfacial tension γ₁₂ = 0.0036 N m⁻¹ (γ₁₂ = 0.0075 N m⁻¹ in the present case) at the lower layer mean thickness h₀ = 0.1 μm (h₁₀ in the present case). Note here that the present formulation is in good agreement with their work. However, the slight deviation in the dispersion curve at higher wavenumbers might be attributed to the fact that the existing work has solved the equations at the interface in the semi-infinite upper film thickness.

Fig. 2 Asymptotic validation plots for the present formulation via dispersion relations comparing the linear growth rate (ω, s⁻¹) as a function of wave number (k, m⁻¹). All the solid color hollow circle symbols denote the previous studies, and the solid (cyan) line represents the NLC–water bilayer system of the present work. Plot (a) shows the validation with Stergios and Brian's (1989)⁸⁰ formulation to compare the bilayer RTI in the limit of isotropic films, wherein ρ₁ (upperlayer) > ρ₂ and ρ₁ is the density of the NLC layer in the present case. Plot (bii) shows the validation result with a single NLC film on the solid substrate studied by Ben and Cumming (2001)⁵⁶ and validated in the limit of h₁₀ → 0 and using P–H director orientation, as shown in plot (bi). Plot (c) shows the isotropic bilayer validation with Bandyopadhyay et al. (2005)³⁰ in the limit of K_f → 0, θ₁ → 0 and θ₂ → 0 considering the nano-thin film. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Plot (bii), on the other hand, is for the case where a previous work⁵⁶ has dealt with the single layer of NLC on the solid substrate and modeled the system using lubrication approximation. The study has considered two cases – gravity in perpendicular and parallel directions to the film, wherein they found that instability may arise for the former case, but instability always prevails for the latter case. However, in the present work, we find different modes of instability for the case, pointing towards the different NLC behavior when in contact with the liquid substrate instead of the solid counterpart. For this case, we validate our results within the limit of h₁₀ → 0 (water film thickness) for the P–H director orientation, as shown in plot (bi).

In order to validate the nano-thin film bilayer configuration and the characteristic instability emerging from the van der Waals interactions between films, we verify the present formulation with the literature.³⁰ The result in plot (c) shows an agreement in the isotropic limit with case 3 for the surface tension combination γ₂ > γ₁ > γ_s. The bilayer instability of NLC film atop a water film is examined through a theoretical framework based on the thickness ranges of the films, as mentioned in the problem formulation section. The forces that drive instabilities are highly sensitive to the thickness of the films, so to achieve a thorough understanding of the instability dynamics, the present study is divided into two sections: micro-thick and nano-thin bilayers. Analyzing these two regimes allows us to capture each case's distinct instability mechanisms accurately.

B. Micro-thick films

For the micro-thick bilayers, the gravitational forces arising from the density difference between the NLC and water layers play a dominant role in destabilizing the system. For example, a higher-density 5-cyano-bisphenyl (5CB) film on a water layer may manifest an RTI. The length and time scales of the instability have been decided by the competition between gravitational forces with the surface and interfacial tensions at the interfaces, along with the elastic rigidity of the NLC film. Thus, the focus of discussion in this section has been on parameters such as the density difference, film thickness, surface and interfacial tensions, and Frank's elasticity of the NLC layer. As mentioned before, the free bilayer system considered here has a pair of coupled deformable interfaces – NI and NW, respectively. Thus, there is a possibility that both interfaces may assume a similar length or time scale when they are fully coupled or may assume different length and time scales if they are not. Previous studies indicate that while the coupled interfaces may deform following a single length and time scale, the decoupled ones may lead to bimodality, wherein the interfaces may deform following different lengths and time scales.^12,30 In the latter case, the NI surface is expected to assume a longer wavelength owing to a larger surface tension, whereas the NW interface may assume a shorter one owing to a weaker stabilizing interfacial tension.

In this direction, Fig. 3 shows some interesting results wherein such free bilayers are considered with two different director orientations at the boundaries of the NLC layer, homeotropic–homeotropic (H–H – θ₁ = 0, θ₂ = 0) and planar–homeotropic (P–H – θ₁ = 0, θ₂ = π/2), as shown schematically in the first column of the images. The P–P and H–P configurations at the boundaries are not explored because of their remote possibility. Furthermore, angular anchoring is also not considered because such director orientation shows results similar to those of homeotropic anchoring.

Fig. 3 The dispersion curves and neutral stability for the NLC–water bilayer configuration under varying layer thickness and director orientations. Schemes (i) and (ii) show NLC director (solid green) orientations under P–H and H–H conditions, respectively. Plots (a) and (c) show the growth rate (ω) versus wavenumber (k) for varying water layer thickness (h₁₀, μm) for the respective director configurations at constant h₃₀ = 10 μm, indicating increased instability at higher thicknesses. Plots (b) and (d) depict the corresponding (to plots (a) and (c)) contour of the critical wavenumber (k_c, m⁻¹) on the plane k_c–h₁₀ showing the influence of thickness on stability. Plots (e) and (g) show the ω for varying NLC layer thickness (h₃₀, μm) in both configurations at constant h₁₀ = 10 μm, showing a stabilizing effect with increasing thickness. Plots (f) and (h) depict contour plots of k_c for varying h₃₀ on the plane k_c–h₃₀, highlighting the stabilizing role of the NLC layer in both configurations. Plots (c), (e) and (g) show bimodality, indicating the presence of FSM and CIM instability modes in the bilayer system shown on the respective curves. Arrows (solid color with dash-dotted lines) on the particular curves indicate increasing and/or decreasing trends. Besides contour plots, color bars indicate the respective contour curve values. The neutral stability curves (ω = 0) are represented by the dash-dotted lines (bold black) on the respective contour plots demarcating the region of stability (ω < 0) from the unstable (ω > 0) region. Schemes (j)–(l) illustrate the in-phase (coupled) and out-of-phase (decoupled) deformation modes of a two-layer thin film system supported on a solid substrate, highlighting the mechanistic origin of unimodal and bimodal responses. The other dash-dotted color lines (green) present the respective ω values. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Table 1 Typical parameter values to perform calculations for bilayer films^45,46,53,67

Parameters	Units	Values
α 1	Pa s	0.0065
α 2	Pa s	−0.007
α 3	Pa s	−0.001
α 4	Pa s	0.085
α 5	Pa s	0.05
γ 1	N m^−1	0.072
γ 2	N m^−1	0.033
γ s	N m^−1	0.0265
ρ 1	kg m^−3	1000
ρ 2	kg m^−3	1020
K f	pN	5
μ 1	Pa s	0.001
μ 2 = α4/2	Pa s	0.0425
h 10	μm (thick films) or nm (thin films)	5–100
h 20	μm (thick films) or nm (thin films)	5–100
g	m s^−2	9.81

---

Fig. 3(a), (c), (e) and (g) shows the variations in the growth rate (ω) and wavenumber (k) for the P–H and H–H systems. The plot (a) shows that, for the P–H system, both the NI surface and NW interface deform by assuming a single mode with a distinct single maximum, ω_m, when the films are of similar thicknesses. The plots also show that, as the lower layer thickness, h₁₀, increases under a thin upper layer, h₃₀ = 30 μm, a reduction in the frictional influence at the bottom water layer rapidly increases ω while the length scale shifts more towards the shorter wavelength regime (larger wavenumber). On the other hand, plot (c) shows that when an H–H configuration is considered, the interfaces may assume a bimodal behavior with a pair of distinct maxima in place. The emergence of bimodality can be elucidated by examining schemes (j)–(l), which reveal that a unimodal response arises when both layers deform cooperatively, i.e., in-phase.^81–84 In contrast, when the layers deform out-of-phase,^81–84 exhibiting relative displacement or antagonistic curvature, the system transitions to a bimodal configuration. In such a scenario, the NI surface may assume an instability mode with a longer wavelength (smaller wavenumber), and the confined NLC–water interface may assume a shorter wavelength (larger wavenumber) mode. However, the system may eventually assume the length and time scale corresponding to the dominant mode among the two. In order to distinguish the pair of modes in a bimodal behavior, we term them as the free surface (FSM) and confined interface (CIM) modes in this work. The plot (c) also suggests that when the h₁₀ is smaller, the dominant FSM mode of instability stays at the NI surface, showing a larger wavelength (smaller wavenumber). However, with the increase in the h₁₀, a mode switching is observed wherein the CIM becomes the dominant one at the NW interface.

The plots (b) and (d) show the contours of ω > 0 in the k vs. h₁₀ plane wherein we clearly distinguish the zones of CIM and FSM dominance. Unlike in plots (a) and (c), where the gravitational forces drive the instability with the increasing h₁₀, in plot (e), the instability is more influenced by the NLC layer thickness. For example, the plot suggests that when both films are thinner and under the P–H configuration, a bimodal behavior may be observed, with the CIM being the dominant mode. However, as h₃₀ increases, the increasing effect of Frank's elasticity at the confined interface may cause the reduction of ω_m. Plot (g) shows a situation wherein an H–H configuration may show the emergence of a dominant FSM mode with increased NLC layer thickness. In this situation, we have employed the term associated with the Ericksen stress as the destabilizing term (K_f(∂θ/∂z)²) in eqn (29) and (31) on the left-hand side, which eventually helps fuel the FSM mode to become the dominant one over the CIM.^{46,47,58,67,85}Fig. 3(b), (d), (f) and (h) shows the situations discussed in their respective dispersion relation plots through contours for ω on the h₁₀–k_c and h₃₀–k_c planes when h₃₀ and h₁₀ are kept constant in the alternative case. Plots (b) and (d) reveal that the CIM mode can be fueled up merely by switching the director orientation from P–H to H–H and increasing h₁₀. In contrast, both the modes appear with increasing NLC-layer thickness, as shown in the plot (f) for the values ∼h₃₀ > 20 μm. Furthermore, plot (h) shows that changing the boundary conditions from P–H to H–H results in a dominant FSM mode with increasing h₃₀.

Fig. 4 shows that the influence of K_f in P–H and H–H configurations arises from the fundamental differences in the director orientation at the interfaces. The asymmetry in director orientation results in an imbalance in the elastic response at the two interfaces. When perturbations occur, the NLC's ability to resist these disturbances is strongly influenced by K_f,⁶⁷ mainly because the elastic forces at the planar interface are weaker.

Fig. 4 The dispersion curves and neutral stability for the NLC–water bilayer configuration under varying NLC bulk elastic constant (K_f) and density differences (Δρ). Plots (a), (c), (e) and (g) depict ω as a function of k for various K_f and Δρ at P–H and H–H director configurations, respectively, at constant film thicknesses h₁₀ = 20 μm, h₃₀ = 10 μm for K_f and h₁₀ = 10 μm, h₃₀ = 60 μm for Δρ, respectively. Plots (b) and (d) depict corresponding (to plots (a) and (c)) neutral stability curves on the plane (k_c–K_f) for P–H and H–H director configurations, respectively. Similarly, plots (f) and (g) depict neutral stability curves (corresponding to plots (e) and (g)) on the plane k_c–Δρ for H–P and H–H director configurations, respectively. Arrows (dash-dotted solid color) on the respective curves indicate increasing and/or decreasing trends. Besides contour plots, colour bars indicate individual contour curve values. The neutral stability curves (ω = 0) are represented by the dash-dotted lines (bold black) on the respective contour plots, demarcating the region of stability (ω < 0) from the unstable (ω > 0) region. The other dash-dotted colour lines (green) present the respective ω values. Insets inside plots (e) and (g) are meant to show distinct curves. An inset inside plot (f) shows the zoomed-in section of the plot representing the values marked inside the rectangular box (pink color). Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

As K_f increases, the overall director field elastic energy becomes larger, especially at the planar interface, where the director is more flexible. This additional elastic energy provides a restoring force that suppresses perturbations,⁶⁷ particularly those with short wavelengths (high k), which are more sensitive to local distortions in the director field. The K_f effectively reduces the ω (plot (a)), especially at higher wavenumbers, by providing more resistance to director distortions. However, the H–H alignment at both interfaces is resistant to deformations; increasing K_f does not significantly alter the stability of the system. The elastic forces are uniformly applied across both interfaces, so even a high K_f value (50 pN) does not introduce any new restoring forces or additional stability, as shown in plot (c). Furthermore, plots (b) and (d) collectively show that irrespective of the K_f magnitude, the system becomes more stable in the H–H director configuration. Nevertheless, an unstable system attempts to nullify the destabilizing effect of increasing K_f indicating its stabilizing role in conjunction with the surface tension. The stabilizing contributions of K_f (the fourth term on the left side) and surface tension (right-hand side term) are also evident in eqn (31). Moreover, we observe that K_f can potentially reduce the CIM arising in addition to the prevailing FSM at higher NLC (h₃₀ ≥ 40 μm) and lower water film (h₁₀ = 10 μm) thicknesses⁶⁷ in contrast to the results shown in plot (a) for h₁₀ = 20 and h₃₀ = 10 μm.

The stability analysis of the P–H and H–H configurations under varying density differences (Δρ = ρ₂–ρ₁) reveals interesting insights into the underlying physics governing these systems. In both the P–H and H–H configurations (plots (e) and (g)), we again observe a bimodal behavior with a pair of distinct maxima at the smaller (CIM) and larger (FSM) wavelengths. Interestingly, when the magnitude of Δρ is smaller, the longer wavelength FSM modes show a dominant behavior. However, with the increase in Δρ, mode switching is observed as the dominant mode of instability shifts from the higher to lower wavelength regime when the CIM becomes the dominant one. Furthermore, the neutral stability contour plots (f) and (h) corresponding to plots (e) and (g) reveal that the mode switching can simply be achieved by tuning the Δρ in the case of P–H as well as H–H configurations. Concisely, Fig. 4 summarizes that the elasticity of the NLC film and the density contrast between the NLC and liquid layer can not only play a key role in altering the length scales from larger to smaller wavelengths but also may cause switching of the dominant mode from the free surface to the confined interface.

Fig. 5 and 6 summarize the predictions of the dominant growth rate – ω_m and wavelength – λ_m for the various parametric combinations discussed in Fig. 3 and 4. Plots (a)–(d) depict the variations of ω_m and λ_m, respectively, with h₃₀ for various h₁₀ (10–30 μm for P–H and 10–50 μm for H–H configurations). Similarly, plots (e)–(h) depict the same results with h₁₀ for various h₃₀ (10–20 μm for both P–H and H–H configurations). As discussed previously, the dominant instability mode is indeed a smaller wavelength CIM, wherein λ_m reduces with the increasing h₁₀ (b) and (d) for both P–H and H–H configurations, whereas the same decreases in the case of increasing h₃₀.

Fig. 5 Variation of the maximum growth rate (ω_m) and the dominant wavelength (λ_m) for parametric variations at different director orientation configurations. Schematics (i) and (ii) show director orientation configurations P–H and H–H, respectively. Plots (a)–(d) depict ω_m and λ_m for h₁₀ variations at different constants h₃₀ for P–H and H–H orientations. Similarly, plots (e)–(h) depict ω_m and λ_m for h₃₀ variations at different constants h₁₀ for P–H and H–H orientations. The inset in (ii-g) shows the prominent mode switching starting from h₁₀ = 16 μm and onward. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Fig. 6 Plots (j)–(n) depict ω_m and λ_m for K_f variations at different h₃₀ while maintaining constant h₁₀ = 20 μm for P–H and H–H orientations. Plots (p) and (q) depict ω_m and λ_m for Δρ variations at different h₃₀ while maintaining constant h₁₀ = 10 μm. For Δρ variations, both P–H and H–H orientations show nearly similar trends without significant differences. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

The increased ω_m may be attributed to the lower frictional influence of the increase in the film thicknesses. On the other hand, the plots also suggest that the presence of the thinner water layers favors longer wavelengths, while thicker layers produce shorter wavelengths for the CIM. Overall, the plots highlight that there can be a pair of instability modes, FSM and CIM, and the transition of the same depends on both NLC and water layer thicknesses.

Fig. 6 displays the impact of K_f and Δρ on ω_m and λ_m for the P–H and H–H configurations. The plots suggest that ω_m decreases and λ_m increases with the increase in the K_f. The rapid decline in the growth rate with increasing K_f suggests that as the NLC elasticity increases, it resists deformation more effectively, thereby slowing the growth of any perturbations. Thinner NLC films are less affected, implying that thinner films are less sensitive to elastic effects, while thicker films stabilize more as elasticity increases. Interestingly, plot (m) shows a slight crossover at h₃₀ between 10 and 20 μm. This crossover is attributed to the FSM dominance, depicting slightly higher growth rates. Beyond the crossover thickness, the system switches to the CIM, where a thinner water layer stabilizes the interface. Consequently, the growth rate trend reverses, making this configuration a stable one. Furthermore, the plots suggest that the bimodality and mode switching can be fueled by the variation in the Δρ. The ω_m increases monotonically with the increase in Δρ for both P–H and H–H configurations with the possibility of transition from the FSM to CIM and switching of the dominant mode from the free surface to the confined interface. In summary, Fig. 3–6 show the key possibilities of different modes of instabilities for the RTI of a free bilayer composed of a nematic and a water film.

C. Nano-thin films

Apart from instigating RTI arising due to the density contrast across a soft interface, the instabilities in the NLC–water free bilayers can also be triggered by employing other body forces. For example, the nano-thin films (<100 nm) can also show the instabilities originating from the van der Waals forces across the films³⁰ because of the insignificance of the gravitational effects in such systems. In the present study, we consider three different varieties of NLC–water free bilayers based on the different spreading coefficients that are most common in the experimental systems – Case (I) S_12s > 0 and S_21a < 0 when the NLC layer has the highest surface tension (γ₂ > γ₁ > γ_s); Case (II) S_12s < 0 and S_21a > 0 when the water layer has the highest surface tension (γ₁ > γ₂ > γ_s); and Case (III) S_12s < 0 and S_21a < 0 when the substrate has the highest surface tension (γ_s > γ₂ > γ₁). Here S_lmn (A_lmn) is the spreading coefficient (effective Hamaker constant) of the layer l sandwiched between the bounding layers m and n, wherein S_lmn < 0 and A_lmn > 0 (S_lmn > 0 and A_lmn < 0) signify an unstable (stable) ‘lth’ layer. The Hamaker constant values used in nano-thin film analysis are calculated from formulas enlisted in Section 3(B-ii).

Case I (S_12s > 0 and S_21n < 0; γ₂ > γ₁ > γ_s): the spreading coefficients chosen in this case indicate that the free bilayer has a stable liquid layer while the NLC layer can be unstable when the van der Waals forces overcome the stabilizing influences arising from the surface/interfacial tensions at the surface/interfaces and the NLC elastic forces.^47,76,77Fig. 7 shows some interesting results regarding the P–H configuration. The plots (a) and (b) suggest that when the K_f is vanishingly small, we observe a bimodal behavior for this case, as has been seen previously in the literature for a purely viscous free bilayer.³⁰ Physically, the bimodality indicates the presence of competing instability modes at both the NI surface and NW interface, in which the shorter-wavelength CIM is the dominant one at the thickness combination chosen for the bilayer. However, with the increase in the K_f, as the deformability of the NLC layer reduces, the dominant mode's length scale progressively switches towards the FSM while suppressing the CIM. The phenomenon resembles the mode-switching behavior observed for the micro-thick films. Fig. 7(c) and (d) shows that such mode-switching behaviors, in this case, can also be achieved by changing either h₁₀ or h₃₀, which is signified by the discontinuities in the plots. A thinner upper layer on a thicker lower layer may ensure a smaller wavelength and faster CIM, while a thicker upper layer may assume a larger wavelength FSM.

Fig. 7 Dispersion relation results, the maximum growth rate (ω_m), and dominant wavelength (λ_m) for the variation of the NLC bulk elastic constant (K_f) considering Case I (γ₂ > γ₁ > γ_s) at the P–H director orientation configuration. In plot(a), the variation of the linear growth rate (ω, s⁻¹) as a function of wave number (k, μm⁻¹) at different K_f values for fixed h₁₀ = h₃₀ = 10 nm upper NLC and solvent film thicknesses. Plot (b) displays the expanded view of the short wavenumber in the dispersion curve for K_f variation. Plots (c) and (d) are for the variation of the ω_m and λ_m with the K_f at h₃₀ = 5–25 nm NLC on h₁₀ = 10 nm solvent film. Similarly, plots (e) and (f) display the variation of ω_m and λ_m with K_f for a h₁₀ = 5–50 nm solvent layer, under a 10 nm thick NLC film. The different symbols in plots (e) and (f) merely depict the distinction between the overlapping trends. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Fig. 8 highlights the role of the film thickness for a given K_f. Fig. 8(a)–(c) shows that for a thin and fixed h₁₀, at the smaller values of h₃₀ the smaller wavelength CIM is the dominant mode of instability, but with increasing h₃₀, mode switching from the CIM to the FSM takes place with λ_m shifting towards the longer wavelength regime. On the other hand, Fig. 8(d)–(f) shows another interesting scenario of mixed-mode instability, wherein a thin h₃₀ on a thicker h₁₀ shows a smaller λ_m, which progressively increases with the increase in h₁₀. Fig. 9 shows the same systems as in Fig. 7 and 8, however, with P–A and P–P configurations at the interfaces, as depicted schematically. A comparison between the results of Fig. 7 and 8 help conclude that P–P, P–H, and P–A configurations essentially show very similar behaviors in terms of mode-switching if other parameters are kept very similar. The figures also suggest that the conditions for mode switching and the associated length (λ_m) and time scales (ω_m) can be altered simply by tuning the director orientation at the boundaries.

Fig. 8 Dispersion relation results, the maximum growth rate (ω_m), and dominant wavelength (λ_m) for the variation of the respective layer thicknesses considering Case I (γ₂ > γ₁ > γ_s) at the P–H configuration. Plot (a) shows the variation of the linear growth rate (ω, ms⁻¹) as a function of wave number (k, μm⁻¹) at a fixed solvent layer thickness h₁₀ = 10 nm with h₃₀ = 5–15 nm NLC layer thickness. Plots (b) and (c) show the variations of ω_m and λ_m with h₃₀ for a h₁₀ = 10–30 nm thick solvent layer. Plot (d) shows the variation of ω with k for a fixed NLC layer thickness of h₃₀ = 10 nm with h₁₀ = 5–20 nm solvent layer thicknesses. Plots (e) and (f) show the variations of ω_m and λ_m with the h₃₀ = 9–17 nm thick NLC layer. The different symbols in the plot (c) merely depict the distinction between the overlapping trends. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Fig. 9 Dispersion relation results, the maximum growth rate (ω_m), and dominant wavelength (λ_m) for the variation of layer thicknesses considering Case I (γ₂ > γ₁ > γ_s.) at P–A (plots (a)–(f)) and P–P (plots (j)–(f)) director configurations, respectively. Plot (a) shows the variation of the ω, ms⁻¹ as a function of k, μm⁻¹ at a fixed lower layer thickness h₁₀ = 10 nm with varying NLC layer thickness in the range h₃₀ = 5–15 nm. Plots (b) and (c) show the variations of ω_m and λ_m with a h₃₀ = 10–30 nm thick solvent layer. Plot (d) shows the variation of ω with k for a h₃₀ = 10 nm fixed upper NLC layer and h₁₀ = 5–20 nm lower layers. Plots (e) and (f) show the variations of ω_m and λ_m with the h₁₀ = 8–22 nm thick upper NLC layer. Plots (g) and (h) show the variations of ω_m and λ_m with h₃₀ = 10–30 nm thick solvent layers. Plots (i) and (j) show the variations of ω_m and λ_m with the h₁₀ = 8–22 nm thick NLC layers. The different symbols in the plot (c) merely depict the distinction between the overlapping trends. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Interestingly, Fig. 8 and 9 highlight the critical influence of orientational anisotropy in governing the transition from a mixed-mode response to distinct individual mode-switching, observed upon varying the orientational order while keeping all other system parameters fixed.

Case II (S_12s < 0 and S_21a > 0; γ₁ > γ₂ > γ_s): the characteristics of a bilayer system under the spreading conditions γ₁ > γ₂ > γ_s, wherein the NLC layer is placed between the water layer and the air, illustrated in Fig. 10. The resulting negative spreading coefficient (S_12s < 0) in this case indicates the unstable configuration. When the NLC layer comes into contact with the substrate, it also results in an unstable case – S_2sa < 0, that is the reason bimodality appears. However, in the presence of a water layer, the NLC layer shows stable behavior. On the other hand, in the absence of a water layer, NLC may destabilize; as a result, the upper interface may show instability.

Fig. 10 Dispersion relation results, the maximum growth rate (ω_m), and dominant wavelength (λ_m) for the variation of layer thicknesses considering Case II (γ₁ > γ₂ > γ_s) at P–H and P–A director configurations, respectively. Furthermore, plots (a) and (d) show variations of ω, ms⁻¹ as a function of k, μm⁻¹ for h₁₀ and h₃₀ variations, respectively. Plots (b), (c), (e) and (f) illustrate the variation of ω_m and λ_m for the P–H (i) configuration. Similarly, plots (g)–(k) display the variation of ω_m and λ_m for the P–A (ii) configuration. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

Fig. 10(a)–(k) summarizes the prominent features of Case II when the layer thicknesses are varied at P–H and P–A director orientations. The increasing h₁₀ initially shows a unimodal behavior at small water layer thickness, which is a longer wavelength (shorter wavenumber) FSM. As the h₁₀ increases, the mode switches towards a shorter wavelength (larger wavenumber), that is, the CIM slowly becomes dominant (plot (a)). On the other hand, for thinner NLC layers, the FSM initially shows dominance, driven by van der Waals near the NI surface. However, as the h₃₀ increases, mode-switching is observed, where the CIM begins to dominate due to increased destabilization at the NW interface (plot (d)). Furthermore, ω_m and λ_m indicate that the thicker NLC layers correspond to the instability of the water layer and exclusively depend on the water layer thicknesses, as shown in plots (b), (c), (g) and (h) for both P–H and P–A, respectively. In the case of small h₃₀, ω_m (plots (b) and (g)) and λ_m (plots (c) and (h)) show the mode switching behavior as the unstable water layer largely controls the instability. This instability grows stronger by decreasing h₁₀, as shown in plots (e), (f), (j) and (k). This is attributed to the fact that the increasing h₁₀ results in weaker intermolecular disjoining pressure, and NLC upper layer instability weakens because the NLC layer on a thick water layer configuration is stable.

Case III (S_12s < 0 and S_21a < 0; γ_s > γ₂ > γ₁): this case corresponds to the surface tension ordering γ_s > γ₂ > γ₁, wherein both the layers are unstable, closely resembling Case-I, as shown in Fig. 8 and 9. In this scenario, a thin water lower layer shows stable behavior in the absence of the upper NLC layer. However, it is unstable in the presence of a thicker NLC layer as the spreading coefficient becomes negative, S_12s < 0. Moreover, if the water layer is absent, the NLC layer is stable but becomes unstable on a thicker lower layer. Given this background, it can be expected that the overall dominating instability mechanism is essentially controlled by the relatively thin NLC upper layer unstable behavior. In contrast, the lower layer is confined between two bounding media, and the water layer instability indeed controls the system for a comparatively thick NLC layer.

In view of this, Fig. 11 summarizes the linear stability analysis for Case III bilayer configuration and surface tension arrangements. The dominant mode of instability is certainly CIM (shorter wave), and the wavelength decreases with either of the layer thickness (h₁₀ or h₃₀) elevation for both P–P and P–A director configurations. However, the present surface tension and layer thickness arrangement depict the presence of bimodality through the occurrence of FSM, which increases with either of the layer thickness, as shown in plots (a) and (d). Moreover, plots (b)–(f) and (g)–(k) depict that ω_m and λ_m gradually become nearly independent of the water layer thickness, h₁₀, in the presence of relatively smaller h₃₀, indicating that the upper NLC thickness of a bulk of water layers largely controls the instability. On the other hand, plots also show that ω_m and λ_m gradually become independent of the h₃₀ when the water layer thickness, h₁₀, is small enough. As a result, the instability in this case is entirely controlled by the lower water layer, which is confined between the top and bottom bulk media. Interestingly, the figure also indicates that the mode switching and associated ω_m and λ_m are the cumulative outcomes of the net influences of film thickness tuning and director arrangements. However, for physically realistic parameter regimes, the long-wavelength (short-wavenumber) maxima remain subdued in magnitude, as can be seen in plots (a) and (d). Consequently, a clear manifestation of mode-switching from FSM to CIM is not observed in ω_m and λ_m plots (b)–(f) and (g)–(k) for the P–P and P–A director configurations. This suppression of distinct mode transitions may be attributed to director anchoring, which reduces the effective elasticity role within the layer.

Fig. 11 Dispersion relation results, the maximum growth rate (ω_m), and dominant wavelength (λ_m) for the variation of layer thicknesses considering Case III (γ_s > γ₂ > γ₁) for k_f = 0.05 pN considering P–P and P–A director configurations, respectively. Plots (a) and (d) show variations of the ω, ms⁻¹ as a function of k, μm⁻¹ for h₁₀ and h₃₀ variations, respectively. Furthermore, plots (b), (c), (e) and (f) illustrate the variation of the ω_m and λ_m for the P–H (i) configuration. Similarly, plots (g)–(k) display the variation of ω_m and λ_m for the P–A (ii) configuration. Unless otherwise stated, the other essential parameters used for the plots are displayed in Table 1.

As the dispersion relation intrinsically depends on the density contrast, viscosity, and interfacial tension between the two layers, variations in these parameters can significantly influence the system's stability characteristics if water is replaced with another fluid. For example, (i) a more viscous subphase (e.g., glycerol) would suppress short-wavelength perturbations, favoring long-wave-dominated instability; (ii) a less viscous and less dense subphase (e.g., a light oil) would diminish gravitational effects and enhance capillary and elastic contributions, thereby promoting finite-wave modes; and (iii) modifications in interfacial tension would alter the stabilizing term in the dispersion relation, consequently shifting the growth rate and the wavelength of the most unstable mode. However, the qualitative nature of the instability, particularly the coexistence and possible transition between the CIM and the FSM, is expected to remain robust. This robustness arises from the fundamental competition between viscous–gravitational destabilization and elastic–capillary stabilization that governs the bilayer dynamics. Hence, while quantitative aspects such as growth rates and characteristic wavelengths may vary with the choice of subphase, the essential physical mechanisms delineated in this study would largely persist. Furthermore, the anisotropy in the present system manifests through elastic and viscous stresses (eqn (11)), rather than through solid-like bending rigidity as happens in transversely isotropic elastic films. Thus, no direct quantitative comparison can be drawn between the two systems, although the presence of anisotropy conceptually connects them.

5. Conclusions

We investigate two asymptotic regimes – micro-thick and nano-thin bilayer systems revealing two distinct forms of RTIs – one driven by density differences and the other by intermolecular force contrasts. Mode switching from the free surface – the FSM to the confined interface – the CIM governed by the interplay of elasticity, surface tension combinations, and film thickness is observed. The analysis is performed using linear stability analysis and coupled Ericksen–Leslie and Navier–Stokes equations across different director configurations. In micro-thick films, mode selection is influenced by layer elasticity, density contrast, and thickness, revealing bimodality emergence for out-of-phase deformations and unimodal behavior for in-phase relative deformation. For nano-thin bilayers, three cases are classified based on spreading coefficient magnitudes. These cases, combined with stability analysis, assist in understanding the instability criteria across physical settings and asymptotic film behaviors when regimes where one layer dominates due to vanishing or excessive thickness. Notably, Case I shows the existence of the bimodality even at vanishingly small K_f values, with mode switching emerging as K_f increases due to reduced NLC layer deformability in the P–H configuration. Case II highlights the distinctive mode switching at P–H and P–A, indicating the predominance of layer thickness that can be harnessed to control the instability. Case III underlines the subtle role of spreading asymmetry and director orientation at boundaries to modulate the instability pathways. Comparative analysis of cases I–III demonstrates that the layer thickness, director orientations, and spreading condition serve as effective handles to modulate interfacial coupling and mode switching in bilayer systems. The identified linear modes have far-reaching consequences on the morphological evolution of thin films, representing the main uniqueness of the present investigation. The findings offer a predictive and tunable framework for precise control of the length and time scales of the instability-driven morphologies in layered NLC systems, with broader relevance to soft matter design, microfluidic patterning, and anisotropic functional material engineering involving thin films. However, fully nonlinear simulations of the present system could reveal complex interfacial morphologies arising from the nonlinear coupling of the two modes. A comprehensive nonlinear nemato-hydrodynamic analysis thus represents a promising direction for future investigation.

Author contributions

VBV and DB conceived the project. VBV, NY, and AG have performed derivations. VBV and NY performed all the simulations and analyzed and validated the results. VBV, NY, DB, and SD wrote the manuscript. DB acquired funding, supervised the progress and revised the final draft. TKM supervised the progress.

Conflicts of interest

The authors do not have any conflicts to declare.

Data availability

All the data utilized in this work are incorporated in the main manuscript. The detailed derivations concerning the manuscript are given in the supplementary information (SI). Supplementary information: the detailed derivation steps of longwave linear stability analysis and dispersion relations – Appendix A, B, and C. See DOI: https://doi.org/10.1039/d5sm00634a.

Acknowledgements

We thank the Government of India for providing financial support through MeitY Grant No. 5(1)/2022-NANO and ICMR Grant No. 5/3/8/20/2019-ITR. VBV is incredibly thankful to Dr Siddharth Thakur for engaging in discussions about the experimental relevance of the theoretical NLC bilayer systems.

References

1. H. M. van der Kooij and J. Sprakel, Soft Matter, 2015, 11, 6353–6359.
2. J. Marthelot, E. F. Strong, P. M. Reis and P.-T. Brun, Nat. Commun., 2018, 9, 4477.
3. S. J. Weinstein and K. J. Ruschak, Annu. Rev. Fluid Mech., 2004, 36, 29–53.
4. J. Ralston, S. S. Dukhin and N. A. Mishchuk, Adv. Colloid Interface Sci., 2002, 95, 145–236.
5. C. Jiang, X. Wang, R. Gunawidjaja, Y.-H. Lin, M. K. Gupta, D. L. Kaplan, R. R. Naik and V. V. Tsukruk, Adv. Funct. Mater., 2007, 17, 2229–2237.
6. A. Sharma and E. Ruckenstein, J. Colloid Interface Sci., 1986, 113, 456–479.
7. A. Sorrentino, in Nanocoatings and Ultra-Thin Films, ed. A. S. H. Makhlouf and I. Tiginyanu, Woodhead Publishing, 2011, pp. 203–234.
8. P. Tchoukov, F. Yang, Z. Xu, T. Dabros, J. Czarnecki and J. Sjöblom, Langmuir, 2014, 30, 3024–3033.
9. S. R. Nagel, Rev. Mod. Phys., 2017, 89, 025002.
10. A. Sharma and R. Khanna, Phys. Rev. Lett., 1998, 81, 3463.
11. A. Sharma, Langmuir, 1993, 9, 861–869.
12. D. Bandyopadhyay and A. Sharma, J. Chem. Phys., 2006, 125, 054711.
13. R. Konnur, K. Kargupta and A. Sharma, Phys. Rev. Lett., 2000, 84, 931.
14. C. Zhao, J. E. Sprittles and D. A. Lockerby, J. Fluid Mech., 2019, 861, R3.
15. C. Zhao, Y. Zhang and T. Si, J. Fluid Mech., 2023, 954, A46.
16. R. Mukherjee and A. Sharma, Soft Matter, 2015, 11, 8717–8740.
17. R. Martinez-Cuenca, G. Saavedra, M. Martinez-Corral and B. Javidi, J. Disp. Technol., 2005, 1, 321–327.
18. S. Maity, T. Bhuyan, R. Bhattacharya and D. Bandyopadhyay, ACS Appl. Mater. Interfaces, 2021, 13, 19430–19442.
19. X. Fan and I. M. White, Nat. Photonics, 2011, 5.
20. X. Meng, Y. Lu, B. Yang, G. Yi and J. Jia, ACS Appl. Mater. Interfaces, 2010, 2, 3467–3472.
21. S. Thakur, S. Rarotra, M. Bhattacharjee, S. Mitra, G. Natu, T. K. Mandal, A. K. Dasmahapatra and D. Bandyopadhyay, Phys. Rev. Appl., 2018, 10, 064012.
22. Y. Xia and G. M. Whitesides, Annu. Rev. Mater. Sci., 1998, 28, 550–575.
23. N. Sabet, H. Hassanzadeh, A. De Wit and J. Abedi, Phys. Rev. Lett., 2021, 126, 094501.
24. V. B. Vanarse, S. Thakur, A. Ghosh, P. R. Parmar and D. Bandyopadhyay, Phys. Fluids, 2024, 36, 022115.
25. P. Tabeling, G. Zocchi and A. Libchaber, J. Fluid Mech., 1987, 177, 67–82.
26. O. K. Matar and S. M. Troian, Chaos, 1999, 9, 141–153.
27. S. S. Usgaonkar, C. J. Ellison and S. Kumar, Langmuir, 2021, 37, 6660–6672.
28. D. I. Pullin, J. Fluid Mech., 1982, 119, 507–532.
29. N. Bhandaru, A. Das and R. Mukherjee, Nanoscale, 2015, 8, 1073–1087.
30. D. Bandyopadhyay, R. Gulabani and A. Sharma, Ind. Eng. Chem. Res., 2005, 44, 1259–1272.
31. D. Bandyopadhyay, A. Sharma and V. Shankar, J. Chem. Phys., 2008, 128, 154909.
32. J. Léopoldès and P. Damman, Nat. Mater., 2006, 5, 957–961.
33. A. R. Piriz, J. J. L. Cela, O. D. Cortázar, N. A. Tahir and D. H. H. Hoffmann, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 056313.
34. Y. Chao, L. Zhu and H. Yuan, Phys. Rev. Fluids, 2021, 6, 064001.
35. A. Pototsky, M. Bestehorn, D. Merkt and U. Thiele, J. Chem. Phys., 2005, 122, 224711.
36. D. Bandyopadhyay, A. Sharma, U. Thiele and D. S. Reddy, Langmuir, 2009, 25, 9108–9118.
37. D. Bandyopadhyay, A. Sharma and V. Shankar, EPL, 2010, 89, 36002.
38. D. Bandyopadhyay, P. D. S. Reddy, A. Sharma, S. Joo and S. Qian, Theor. Comput. Fluid Dyn., 2012, 26, 23–28.
39. D. Bandyopadhyay, P. D. S. Reddy and A. Sharma, Phys. Fluids, 2012, 24, 074106.
40. P. Lambooy, Phys. Rev. Lett., 1996, 76, 1110–1113.
41. M. Sferrazza, Phys. Rev. Lett., 1998, 81, 5173–5176.
42. A. Das, A. B. Dey, G. Manna, M. K. Sanyal and R. Mukherjee, Macromolecules, 2022, 55, 1657–1668.
43. N. Bhandaru and R. Mukherjee, Macromolecules, 2021, 54, 4517–4530.
44. S. Roy, D. Biswas, N. Salunke, A. Das, P. Vutukuri, R. Singh and R. Mukherjee, Macromolecules, 2013, 46, 935–948.
45. S. Chandrasekhar, Liquid Crystals, Cambridge University Press, 1992.
46. P. G. D. Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, Oxford, 1993.
47. F. Vandenbrouck, M. P. Valignat and A. M. Cazabat, Phys. Rev. Lett., 1999, 82, 2693.
48. B. Ravi, R. Mukherjee and D. Bandyopadhyay, Soft Matter, 2015, 11, 139–146.
49. S. Žumer and J. W. Doane, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 34, 3373–3386.
50. C. Poulard and A. M. Cazabat, Langmuir, 2005, 21, 6270–6276.
51. A. D. Rey and E. E. Herrera-Valencia, Soft Matter, 2014, 10, 1611–1620.
52. P. G. D. Gennes, Rev. Mod. Phys., 1985, 57, 827.
53. S. Thakur, A. K. Dasmahapatra and D. Bandyopadhyay, ACS Appl. Mater. Interfaces, 2021, 13, 60697–60712.
54. V. B. Vanarse, S. Thakur and D. Bandyopadhyay, Langmuir, 2025, 41, 5933–5946.
55. R. Bolleddu, S. Chakraborty, M. Bhattacharjee, N. Bhandaru, S. Thakur, P. S. Gooh-Pattader, R. Mukherjee and D. Bandyopadhyay, Ind. Eng. Chem. Res., 2020, 59, 1902–1913.
56. M. Ben Amar and L. J. Cummings, Phys. Fluids, 2001, 13, 1160–1166.
57. L. J. Cummings, Eur. J. Appl. Math., 2004, 15, 651–677.
58. T. S. Lin, L. J. Cummings, A. J. Archer, L. Kondic and U. Thiele, Phys. Fluids, 2013, 25, 082102.
59. T.-S. Lin, L. Kondic, U. Thiele and L. J. Cummings, J. Fluid Mech., 2013, 729, 214–230.
60. J. L. Ericksen, Arch. Ration. Mech. Anal., 1962, 9, 371–378.
61. J. L. Ericksen, Trans. Soc. Rheol., 1967, 11, 5–14.
62. F. M. Leslie, Continuum Mech. Thermodyn., 1992 DOI:10.1007/BF01130288.
63. C. W. Oseen, Trans. Faraday Soc., 1933, 29, 883–899.
64. F. C. Frank, Discuss. Faraday Soc., 1958 10.1039/DF9582500019.
65. H. Yokoyama, S. Kobayashi and H. Kamei, Mol. Cryst. Liq. Cryst., 1985, 129, 109–126.
66. V. A. Raghunathan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1995 DOI:10.1103/PhysRevE.51.896.
67. K. Mondal, A. Ghosh, J. Chaudhuri and D. Bandyopadhyay, J. Fluid Mech., 2018, 834, 464–509.
68. D. H. Sharp, Phys. D, 1984 DOI:10.1016/0167-2789(84)90510-4.
69. L. Rayleigh, Philos. Mag., 1882, 14, 184–186.
70. G. Taylor, Proc. R. Soc. London, Ser. A, 1950, 201, 192–196.
71. P. Dhara and R. Mukherjee, London, Edinburgh Dublin Philos. Mag. J. Sci., 2020, 124, 1293–1300.
72. M. J. Stephen and J. P. Straley, Rev. Mod. Phys., 1974, 46, 617–704.
73. S. J. Tavener, T. Mullin, G. I. Blake and K. A. Cliffe, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 63, 011708.
74. L. J. Cummings, T.-S. Lin and L. Kondic, Phys. Fluids, 2011, 23, 043102.
75. D. van Effenterre, M. P. Valignat and D. Roux, Europhys. Lett., 2003, 62, 526–532.
76. M.-P. Valignat, D. Duca, R. Barberi, A.-M. Cazabat and R. Bartolino, Mol. Cryst. Liq. Cryst. Sci. Technol., Sect. A, 1997, 301, 151–156.
77. M. P. Valignat, S. Villette, J. Li, R. Barberi, R. Bartolino, E. Dubois-Violette and A. M. Cazabat, Phys. Rev. Lett., 1996, 77, 1994–1997.
78. E. Perez, J. E. Proust and L. Ter-Minassian-Saraga, Colloid Polym. Sci., 1978, 256, 784–792.
79. C. J. Van Oss, M. K. Chaudhury and R. J. Good, Chem. Rev., 1988, 88, 927–941.
80. S. G. Yiantsios and B. G. Higgins, Phys. Fluids, 1989, 1, 1484–1501.
81. C. Maldarelli, R. K. Jain, I. B. Ivanov and E. Ruckenstein, J. Colloid Interface Sci., 1980, 78, 118–143.
82. K. D. Danov, V. N. Paunov, N. Alleborn, H. Raszillier and F. Durst, Chem. Eng. Sci., 1998, 53, 2809–2822.
83. K. D. Danov, V. N. Paunov, S. D. Stoyanov, N. Alleborn, H. Raszillier and F. Durst, Chem. Eng. Sci., 1998, 53, 2823–2837.
84. V. N. Paunov, K. D. Danov, N. Alleborn, H. Raszillier and F. Durst, Chem. Eng. Sci., 1998, 53, 2839–2857.
85. S. Herminghaus, K. Jacobs, K. Mecke, J. Bischof, A. Fery, M. Ibn-Elhaj and S. Schlagowski, Science, 1998, 282, 916–919.

---