Complex Nanowrinkling in Chiral Liquid Crystal Surfaces: From Shaping Mechanisms to Geometric Statistics.

Surface wrinkling is closely linked to a significant number of surface functionalities such as wetting, structural colour, tribology, frictions, biological growth and more. Given its ubiquity in nature's surfaces and that most material formation processes are driven by self-assembly and self-organization and many are formed by fibrous composites or analogues of liquid crystals, in this work, we extend our previous theory and modeling work on in silico biomimicking nanowrinkling using chiral liquid crystal surface physics by including higher-order anisotropic surface tension nonlinearities. The modeling is based on a compact liquid crystal shape equation containing anisotropic capillary pressures, whose solution predicts a superposition of uniaxial, equibiaxial and biaxial egg carton surfaces with amplitudes dictated by material anchoring energy parameters and by the symmetry of the liquid crystal orientation field. The numerical solutions are validated by analytical solutions. The blending and interaction of egg carton surfaces create surface reliefs whose amplitudes depend on the highest nonlinearity and whose morphology depends on the anchoring coefficient ratio. Targeting specific wrinkling patterns is realized by selecting trajectories on an appropriate parametric space. Finally, given its importance in surface functionalities and applications, the geometric statistics of the patterns up to the fourth order are characterized and connected to the parametric anchoring energy space. We show how to minimize and/or maximize skewness and kurtosis by specific changes in the surface energy anisotropy. Taken together, this paper presents a theory and simulation platform for the design of nano-wrinkled surfaces with targeted surface roughness metrics generated by internal capillary pressures, of interest in the development of biomimetic multifunctional surfaces.

1. Introduction

Biological cholesteric liquid crystals are organic materials [1] found in human compact bone (collagen) [2,3], beetles’ shells (chitin fibre) [4,5,6], plywood (cellulose) [7,8], virus RNA [9,10] and others [11]. Biological cholesteric liquid crystals have been found as a promising material in optics [12], tribology [13,14] and wetting applications [15,16]. The manifestation of molecular or colloidal orientation, chirality, self-assembly in nature’s liquid crystals [17,18,19,20,21,22] has served as a biomimetic model, and reciprocally, liquid crystal material physics has been used to shed light on biological structures [18,23,24,25,26,27,28,29]. Liquid crystals are orientationally ordered materials that flow like viscous liquids but display anisotropy like crystals [23] and usually form from high geometric aspect ratio molecular, supramolecular, and colloidal units, such as fibre, fibrils, and rods. The energetic interaction between the surface geometry and the anisotropy liquid crystal direction generates surface forces [20,30,31,32] and novel modes of intrinsic surface pattern formation explored in this paper.

The alignment of rod-like molecules determines the symmetry of liquid crystal phases. The two most common phases are the nematic liquid crystal (NLC) and cholesteric liquid crystal (CLC). In the NLC phase (Figure 1a), all the molecules tend to align toward a certain symmetry breaking direction, known as the director field or anisotropy axis. However, in the chiral CLC phase, the molecules have a layered structure. The preferred orientation within each layer in Figure 1b has a twisting angle with respect to its adjacent layer. The helix pitch corresponds to the length required to perform a full 360-degree rotation. Therefore, the director field in NCL can be regarded as a constant, while the director field in CLC is a periodic function with respect to its position.

Figure 1

Nematic and cholesteric liquid crystal phases show different symmetry. (a) The nematic liquid crystal phase with an anisotropy axis known as the director. (b) The chiral cholesteric liquid crystal phase with a spatial periodicity .

The cholesteric phase is a preferred hierarchical structure adopted in many biological materials since it integrates mechanical, optical, tribological, and wetting functionalities. The blue colour observed on Morpho butterfly wings is the consequence of the presence of a chiral structure in the chitin fibres [33]. The beetles’ shells, on the other hand, generate a unique structural colour that changes with different polarizer orientation [34]. The chitin fibres are stacked in a chiral structure, therefore showing selectivity to light [5]. Cholesteric chitin fibres can enhance the mechanical properties of some animals’ organs as well. For example, the dactyl clubs of Odontodactylus scyllarus mantis shrimp can accelerate the dactyl heel of 65 to 104 / and can reach a peak speed of 23 / from a stationary position [35,36]. A double-spiral defect pattern is observed in the chitin fibres since it creates a larger surface area that can easily dissipate the total energy when the fracture propagates. The multifunctional cholesteric liquid crystal is therefore a promising biomimetic material that can be applied in many fields.

The self-assembly and self-organization of cholesteric liquid crystal with free surfaces or interfaces usually generates complex multi-directional and multi-wavelength periodic wrinkling patterns [37,38,39,40,41]. The origin of these wrinkles is the presence of anisotropic surface/interfacial tension due to the anisotropic axis, its spatial gradients, and its chirality. These multiscale wrinkling surface reliefs with distinct morphology generate material functionalities in diffraction gratings [42,43,44,45] and elastomers [46]. However, a lack of a fundamental understanding of wrinkling mechanisms and geometry and the challenges in experimental characterization creates barriers to the biomimetic development of liquid crystalline wrinkled surfaces with controllable and targeted geometries. Thus, the opportunity of using liquid crystal material science and self-assembly coupled with computational geometry and high-performance computing offers a new tool in this functional surfaces’ area.

In our previous work [47] on 2D surface wrinkling, using the lowest possible order anisotropic surface energy, we captured the emergence of single and doubly periodic egg carton surfaces of the type:

Furthermore, we showed [47] that in the linear small-amplitude regime, the superposition of these egg carton modes predicted by this simple model leads to surface reliefs rich in geometric diversities, such as minimal surface patches, umbilic points, zero Gaussian curvature patches and more. Encouraged by this previous work, here we extend this previous work on egg carton surfaces considering higher-order anisotropic surface energy models. Two important novel features of surface pattern formation in liquid crystal surfaces, exploited in this paper, are:

The dimensionless anisotropic-to-isotropic interfacial tension anchoring coefficient ratio (r) transforms orientation features into geometric surface features [48]. If r is relatively small, material orientation scales of the order microns generate surface nanowrinkling, otherwise submicron or micron features will arise;

The spatial orientation of the director field usually depends on current or prior self-assembly [49] or self-organization conditions [50] and hence this key input to the surface shape equation can generate surface wrinkling complexity not obtained by external stress loading of elastic surfaces.

Complex wrinkling patterns can be formed by various mechanisms and predicted by various models. For example, the nonlinear membrane elastic model [51], the dynamic approach by using the Landau–Ginzburg dynamics equations [52], the anisotropic spontaneous curvature model [53], and surface wrinkling due to temperature changes [46] or plasma treatment [54]. However, all the approaches mentioned above require an external effect that is independent of the material itself, such as the presence of an external load. In the present paper, we present an intrinsic self-assembly approach where surface wrinkling is only developed by the anisotropic interfacial tension (anchoring effect) and no other surface force fields or external loads are necessary. In other words, the coupling of surface geometry and surface orientation is the geometric morphing agent. In particular, this research focuses on the lyotropic cholesteric liquid crystals inspired by biological systems where the temperature is the ambient temperature. The concentration of the liquid. The typical concentration for this type of material is of the order of 50 mg ml −1 [55,56,57], and the environment (such as humidity) could affect the surface profile [58]. The other type of liquid crystal is the thermotropic liquid crystal, whose properties are affected by temperature, including, for example, the wetting properties [59], anchoring and optics [60], as well as coatings and surface roughness [61].

Surface roughness plays a critical role in various biological and biomimetic systems. For example, a biomimetic approach inspired by shark skin improves aerodynamic performance [62]. Multiscale roughness-induced hydrophobicity can be found in lotus and butterfly wings [63], as well as osteogenic biomimetic surfaces [64]. The colours of butterfly wings are due to their multiscale surface patterns [65], and the surface roughness is later shown to determine the imaging quality in an AR/VR display system [66]. The usual characterization of surface roughness geometry is based on the distribution function of the surface relief height and its moments. The statistical properties of the surface geometry are described by its roughness (second-order moment), the skewness (third-order moment) and the kurtosis (fourth-order moment). The average roughness and root-mean-square roughness contains the mean value and its overall deviation from the mean, while the higher order skewness and kurtosis moments capture surface symmetry breaking features that can generate surface functionalities such as static friction [67] and biological surface topography [68]. For example, a surface with positive skewness has a much higher shear strength than a negative skewness given the same average roughness (first-order moment) and root-mean-square roughness (second-order moment) [69]. The surface kurtosis quantifies features such as isolated peaks and valleys. In this research (Section 3.3), we propose an integrated analytical approach to obtain wrinkled surface profiles with rich skewness and kurtosis features.

Based on the previous discussion, Figure 2 summarizes the key issues addressed in this paper: (a) the surface wrinkling pattern formation process is due to an instability of a flat surface above the orientation director field which is unstable and tends to relax to a wrinkled surface. In this paper, we address the mechanism using an efficient shape equation, a complex director orientation and a high order interfacial energy formulation by keeping the quartic term in the surface free energy in Equation (2). (b) The three primary liquid crystal orientation gradients are known as splay, twist and bend; the surface wrinkling process is generated by bend and splay modes whose presence, absence, and relative mixing impact the wrinkling geometry. In this paper, we address how the different modes impact final shapes explicitly, which is crucial in optimizing and controlling surface reliefs. (c) The wrinkled surface can be decomposed into various fundamental wrinkling modes (egg cartons, corrugated, etc.), whose contribution to the final wrinkling surface is captured by a revealing shape -matrix (Equation (15)); the figure shows a matrix with each element corresponding to wrinkling patterns. In this paper, we found optimal ways to synthesize information for the wrinkling patterns using a shape classification based on the superposition of primary modes. (d) Finally (not shown in Figure 2), an in-depth characterization of shape statistic and scaling is provided, aiming at future potential applications in tribology, friction, structural colour and more.

Figure 2

Schematic of this research. (a) A cholesteric liquid crystal tends to have a wrinkled surface rather than a flat surface driven by the director field (short lines) and the gradient of the director field. (b) The gradient of the director field shows three primary elastic modes: splay, bend and twist. (c) The wrinkling surface is the summation of various simple wrinkling modes, whose contributions are integrated into the -matrix of fundamental wrinkling modes (shown in Equation (15)). (d) The higher-order moment of the wrinkling surface characterizes the surface roughness.

The scope of this paper is restricted to equilibrium patterns, and all time-dependent phenomena are not considered. The liquid crystalline phase is described by a director orientation field, and issues arising from molecular ordering are not considered. The director field is a known input to the overall model, and hence, there is no effect of geometry on material structure. The chosen defect-free director field is selected as it manifests a rich potential for wrinkling. The nonlinear (quartic) anisotropic interfacial energy is a well- known [70,71], widely used function of , where is the surface unit normal. Tangential Marangoni forces of any origin are not considered. The essential nature of the model is the liquid crystal shape equation, and the sole focus is on intrinsic capillary pressures; hence, bulk elastic corrections and bulk stress jumps are not considered. Future work can include all these effects as required by experimental data, biological observations and material innovation. An effort is made to highlight mechanisms, novelty and significance, while mathematical details are found in the appendices. Without ambiguity, in this paper, we refer to the anisotropic interfacial tension coefficient as the anchoring coefficient and anisotropic interfacial energy as the anchoring energy.

The organization of this paper is as follows. In Section 2.1, we present a summary of the surface energy expression and the liquid crystal director elastic deformation. The governing equation of surface wrinkling and the mechanism are explained in Section 2.2, with a linear approximation proposed in Section 2.3. The results are presented in Section 3, where we compute the morphology of a wrinkling surface (Section 3.2) and the skewness-kurtosis profiles (Section 3.3). In Section 4, we characterize surface roughness parameters in the surface anchoring (anisotropic interfacial energy) parametric space. Section 5 summarizes key results and concludes with future work. Appendix A presents the fundamental equations that are used to describe the intrinsic geometry, and Appendix B provides details on the wrinkling scaling laws (Equation (19)).

2. Theory and Model Formulation

2.1. Surface Energy and Elastic Deformation

The anisotropic surface free energy of liquid crystal is described by the generalized Rapini–Papoular equation [20,70,71] where is the isotropic surface tension, and are the temperature- and/or chemical composition-dependent anchoring coefficients. The anisotropy is presented in the product of the material director field and geometric surface unit normal . In this paper, we focus on the quartic model such that and . The dimensionless quartic model used in this paper reduces to . The signs of the material parameters are restricted as follows. The thermodynamic stability criterion within the region requires . This criterion implies that not every combination of and has physical meaning. The purple region in Figure 3 represents the possible choices in the anchoring parametric space by assuming that for every . If , Equation (2) degenerates to the quadratic model, which is discussed in [58,72,73,74,75], and it is not included or discussed here.

Figure 3

(a) Stability diagram of the Rapini–Papoular Equation (2) presented in the parametric anchoring energy space . (b) A zoom-in plot of the pink square (, ) in (a). to are critical curves separating different stable zones, where the local surface energy density shows various landscapes summarized and characterized in Table 1.

The stable regions in Figure 3 are separated into 8 different sub-regions based on the location of the local maximum and minimum of the surface energy density . The critical straight lines in Figure 3 are , and . The critical parabolic curve is . Any region that is outside the purple part in Figure 1 does not ensure that holds for every independent director and unit normal . Therefore, those regions are thermodynamically unstable. For nanowrinkling where both and are on the order of ≪, surface tension is always positive, and the stable condition is always satisfied.

Once an anchoring coefficient pair is chosen, the surface energy density profile is uniquely determined and is a function of , where and is equivalent to the system in the framework of the non-polar liquid crystal model. The local minimum or maximum of is only determined by the anchoring coefficient ratio . The position (in terms of ) of the local minimum and maximum of all 8 stable regions separated by 4 critical curves to in Figure 3, A, B1, B2, C, D, E1, E2 and F, are summarized in Table 1.

Table 1

A summary of local energy density extrema in the sub-regions presented in Figure 3.

Regions in Figure 3	Local Minimum	Local Maximum
A	n⊥k	n‖k
B1	(n·k)^2=|r|	n‖k
B2	(n·k)^2=|r|	n⊥k
C	n‖k	n⊥k
D	n‖k	n⊥k
E1	n⊥k	(n·k)^2=|r|
E2	n‖k	(n·k)^2=|r|
F	n⊥k	n‖k

Table 1 demonstrates that the local minimum or maximum of the surface energy density only occurs either at , , or , corresponding to planar (tangential director to the surface) anchoring, homeotropic (perpendicular director to the surface) anchoring and oblique (tilted director to the surface) anchoring, respectively. The local energy density under these three conditions are , 1 and , respectively. The significance of Table 1 is that it shows the parametric bubble where a spatially variable director field samples the energy extrema and the types of extrema. For example, in region B1, one can see that if in one region adjacent to , the energy landscape will drive a geometric change through surface tilting, as shown below.

For an isotropic surface tension, the only capillary pressure is the Laplace pressure generated by the mean curvature or divergence of . For a liquid crystal, we have other sources of capillary pressures [20] generated by the director gradients, captured by the splay , twist and bend deformation given by [76,77]

In partial summary, the parametric anchoring coefficient space (Figure 3a), together with the surface energy (Equation (2)) and the director field and its deformation modes (Equation (3)), are the primary information model needed to analyze nanowrinkling.

2.2. Cahn–Hoffman Capillary Vector Method

The surface shape equation is the interfacial normal stress balance equation. This balance equation is given by the sum of all the capillary forces and the bulk normal stress jump (the difference between normal stresses in the two bulk phases) [20]; see [20] for details.

The Cahn–Hoffman capillary vector is defined by [78,79,80] where is the surface projection tensor. The physical interpretation of Equation (4) is demonstrated with Figure 4; for simplicity, we only consider a quadratic anchoring model. In Figure 4a,b, the tangential component of the Cahn–Hoffman capillary vector gives the direction where the rotation of the surface unit normal around results in the fastest rate of decreasing anchoring energy.

Figure 4

Schematic of the main vectors in the nematic Cahn-Hoffman vector thermodynamics for a (a) planar easy axis () and (b) homeotropic easy axis (). The principal surface frame is selected by the director orientation [80]. Reprinted (figure) with permission from Ae-Gyeong Cheong and Alejandro D. Rey, Physical Review E, 66, 021704 (2002). Copyright (2002) by the American Physical Society.

In the absence of external loads and zero bulk stress jumps, the governing shape implies a surface divergence free , which reveals three fundamental capillary pressures by applying the chain rule in Equation (4): (1) the change in the local surface area (variation of along normal direction) causes dilation pressure; (2) the change in the local surface orientation (variation of along tangent direction) causes rotation pressure, unique to anisotropic media; (3) the director pressure induced by surface gradients in the average molecular orientation [80].

Replacing from Equation (2) to Equation (4) yields the three dimensionless fundamental pressures (scaled by ) in Equation (5) where is the mean curvature. The cancellation of the final net pressure (Equation (5)) indicates that the wrinkling phenomenon arises due to the balance of the three capillary pressures. Equations (6) to (8) show the additional capillary pressures from geometry-orientation couplings such as: , , .

2.3. Linear Approximation

It turns out that given the complexity of the shape in Equation (5), the linear regime is an indispensable tool to characterize, analyze, and shed light on the underlying principles of pattern formation. In the linear region ( and ), the rotation pressure is relatively small compared to the dilation pressure and director pressure. The following identities hold:

Equation (9) and Equation (10) demonstrate that in the director pressure (Equation (8)), the deformation is only presented by the bending elastic mode and splay elastic mode (see Equation (3)). A further approximation yields the Laplacian for the surface relief: where is expanded in the orthonormal basis , and the dimensionless variables are

This generic result shows the direct connections between the surface geometry and director splay/bend deformations. It can be seen that increasing the surface energy nonlinearity “m” with a fixed director field increases the number of shaping modes. In addition, the symmetry properties of the terms in Equation (11) shows that the normal director orientation () plays a unique role in generating curvature.

The shape Equation (11) provides the groundwork for a new way to classify and extract surface wrinkling information which would not be possible with numerical simulations alone. In the linear region, the surface profile can be analytically solved using Equation (11). The resulting is a linear combination of fundamental wrinkling modes where a director field with a doubly helical structure and an m-th order surface energy (Equation (2)) generates the following geometric matrix and basis shape vectors :

The entries of the generic -matrix are -dependent () constants whose values reflect the symmetry and periodicities of the director field and the anchoring model adopted (Equation (2)), as expected from our discussions above on geometry-orientation couplings. For a quartic anchoring model (), the -matrix is , indicating 24 fundamental modes, excluding the trivial flat mode. If a director of the following form is given, we can solve for the -matrix, resulting in a shape -matrix for the given director field

The -matrix contains 24 fundamental wrinkling modes, and the -matrix reduces the entries to 10 nonzero modes due to the symmetry present in the director .

Next, we present a specific application based on Equations (11) and (13). The wrinkling mode is denoted as for clarity. As mentioned above, we use a quartic model (Equation (2)) with two-dimensional parametric space (Figure 3). For the egg carton-type surfaces of interest in this work, we introduce the following compact nomenclature (briefly mentioned in the introduction):

uniaxial egg carton modes: or ;

equibiaxial egg carton modes: and ;

biaxial egg carton modes: and .

Hence, (1) denotes a 1D surface corrugation, (2) denotes the ideal equi-doubly periodic egg carton, and (3) denotes an egg carton with different x- and y-wavelengths, as follows. Assuming (for please consult [47]), the analytic solution () to Equation (11) can be rearranged such that

Here, the shape matrix is the specific realization of the generic -matrix (Equation (13)). The relationship between and is . It is important to note that we have expressed the coefficient matrix in terms of an r-dependent coefficient matrix multiplied with the coefficient . Equation (15) reveals that r only controls 3 modes: [1,2] (biaxial), [2,0] (uniaxial) and [2,2] (equibiaxial). The [1,2] and [2,0] modes vanish at the same time when , and the [2,2] mode does not contribute when . We can conclude that directly controls the surface profile magnitude, while r determines the surface morphology through the coefficient matrix . The nullification of fundamental modes at specific negative r-values indicates the impact of the plus or minus effect and the importance of the II and IV quadrants in the parametric space (Figure 3) to generate shape transitions and extrema.

Figure 5 shows the representation of the wrinkling modes by the generic -matrix (a) and (b) the corresponding specific -matrix. The -matrix is all of the theoretically possible 24 wrinkling modes, and the -matrix, as mentioned above, is the specific coefficient matrix in this paper (determined by the director field in Equation (13) and a quartic anchoring model ). In our adopted model (-matrix), Figure 5b shows that there are only 10 fundamental wrinkling modes: 2 uniaxial modes (cylindrical wrinkling), 2 equibiaxial modes (symmetric wrinkling) and 6 biaxial modes (asymmetric wrinkling). The rest of the wrinkling modes in the -matrix (indicated by black zones in (b)) do not contribute to the surface profile. Importantly, the anchoring coefficient ratio modulates the amplitude of the three modes (see Equation (15)).

Figure 5

Visualization of the general shape - and -matrix in Equations (13) and (15). The position of the wrinkling mode in Figure 5 exactly corresponds to the entry in the - and -matrix. The red zone, blue zone, and yellow zone demonstrate uniaxial wrinkling, equibiaxial wrinkling and biaxial wrinkling, respectively. (a) The -matrix presents all possible wrinkling modes in a quartic model; (b) the -matrix indicates that only 10 wrinkling modes are presented in this model due to the symmetry of the adopted director field.

In partial summary, symmetry analysis combined with basic liquid crystal surface physics can be systematically represented by an shape matrix containing the fundamental wrinkling modes for an m-th anchoring model. The symmetry of the adopted director field then reduces the number of fundamental modes in the -matrix, and the reduced actual shape -matrix eliminates forbidden modes. For the adopted -vector, rich and complex wrinkling modes remain active. The systematic approach formulated here can efficiently be applied to any director field and any anchoring model, as required by future experimental observations and/or computational discovery.

3. Results and Discussion

3.1. Parametric Space and Numerical Methods

The governing Equation (5) is a second-order nonlinear partial differential equation of with respect to two variables , . We replace 0 with in the governing equation, and the numerical solution converges when . We use the EPS method [81,82,83] along the time dimension ( until ) and the finite difference method on the two spatial directions with periodic boundary conditions (matrix index satisfies and , where ). In this research, , .

h shows a translational invariance in Equation (5) such that if , and does not change from Equation (6) to (8), which meets our expectation since lifting or lowering the surface should not affect the morphology. To obtain a consistent result in numerical simulation, the final surface profile is , which ensures . The numerical results were compared with analytical results for the linear regime (see Section 2.3), and in all cases, we find an essentially perfect agreement.

3.2. Surface Profile

Equation (15) implies that only controls the absolute magnitude of the wrinkling profile, while the anchoring coefficient ratio r determines the morphology. Once is fixed, r can be chosen arbitrarily. Here, we focus mainly on (see Figure 3). The numerical solution, analytic linear approximation (via solving Equation (11)) and their cross-sectional plots for , , , and are presented in Figure 6. These specific anchoring coefficient ratio values are chosen because each of them corresponding to a special case in the surface skewness and kurtosis analysis discussed later in Section 3.3. In the parametric space (Figure 3), the specific r-values correspond to: , zone F (); , and are in zone E1 (); and is in zone D (). The average difference ratio Ar is defined as where is the solution matrix with entries, and the division represents element-wise division.

Figure 6

The numerical solution (first row, (a1)–(e1)), analytic solution (second row, (a2)–(e2)) and cross-section plots along a diagonal line ((a3)–(e3)) as well as two boundaries ((a4)–(e4)) and ((a5)–(e5)) with (column (a), min skewness), (column (b), min kurtosis), (column (c), zero skewness), (column (d), max skewness) and (column (e), max kurtosis) by fixing . The percentage values shown in the second row are the average difference ratio Ar defined in Equation (16). The full and dashed-line curves for the cross-sectional plots represent numerical solution and analytic solutions, respectively.

Figure 6 presents the dimensionless surface wrinkling profile as a function of two spatial variables . The first row demonstrates the numerical solution, and the second row shows the analytic solution. The cross-section plots along the diagonal curve and the two boundaries and are shown in the third row of Figure 6. The full and dashed curves represent the numerical cross-section plot and analytic cross-section plot, respectively. The average difference ratio Ar in Figure 6 varies from 0.0114% to 0.0613%, validating the fact that the linear model fits the numerical solution accurately. The dimensionless surface profile has an order of magnitude of , implying a nanowrinkling profile. The results shown in Figure 6 demonstrates the complex wrinkling patterns that result from the superposition of fundamental egg carton shape modes given in Equations (13) and (15). Comparing Figure 6(a1,e1), we find that . This relationship verifies that the -matrix does not contain any symmetry with respect to , unlike the relationship found in the 2D surface wrinkling phenomenon [84].

In partial summary, the model and simulation are able to capture non-ideal complex wrinkling patterns whose symmetries and periodicities can be modulated by the anchoring coefficients. In practical terms, these coefficients are functions of the chemistry, concentration, and temperature, offering a number of routes to fix the value of r and hence target different egg cartons in a more predictable way.

3.3. Skewness and Kurtosis

The statistics of the surface geometry given by moments of the surface relief are required to establish potential applications that depend on roughness. Here the superposition of various egg carton will be shown to generate a rich response to the anchoring coefficient ratio “r”. In what follows, we only consider the morphology factor r in the scaled dimensionless surface profile to eliminate the effect of on , since it only linearly increases the magnitude of .

The surface roughness parameters root mean square (), skewness (Ssk) and kurtosis (Sku) are defined by (with a mean value of equal to 0) where is the metric. If , the following scaling laws hold (see Appendix B)

If , the relationship between the roughness parameter and r is given by Figure 7, and the root mean square increases when increases. If , the skewness profile presents two local minima and two local maxima due to the competition of the , and -modes (Equation (15)), shown in Figure 7(b,b1). The kurtosis profile only presents one local maximum and two local minima for the same reason. The skewness and kurtosis profiles present rich complexity for compared with . The opposite signs between the quadratic term and the quartic term result in a more complicated surface energy profile (the second and fourth quadrants in Figure 3), further creating various behaviours in the surface roughness parameters. Figure 7d presents the skewness-kurtosis plot, which is a prevalent tool in the application of lubrication [85,86], flows [87,88] and boundary layer analysis [89]. Figure 7d demonstrates a global binary plot with a single loop within , showing the complex coupling phenomenon of the two roughness parameters due to the lack of a dominant term for and in Equation (6) to Equation (8).

Figure 7

Surface roughness parameters of are dependent on the anchoring coefficient ratio r. (a) Root mean square (), (b) skewness (Ssk) and (c) kurtosis (Sku) are functions of r. (d) Sku as a function of Ssk parametrized by r. The red and blue curves demonstrate and , respectively. A square marker represents a local minimum, and a filled circle represents a local maximum. A star marker indicates a global extremum. (b1,c1,d1) are the zoom-in details of the region with squares in (b–d), respectively.

The maximum and minimum values of the skewness and kurtosis plots presented in Figure 8 within are analyzed along with the computation of the kernel density (KD) of the surface profile at the corresponding r-value, which is defined by . Here, PDF is the probability distribution function, n denotes the number of bins for histogram, and in this research for a smooth figure. The distribution KD of surface relief is illustrated by the profile in Figure 8. Figure 8a shows a right-leaning curve with a long-left tail (), which validates the fact that the surface profile under has a negative skewness. Similarly, Figure 8b has a maximum, positive skewness, where the KD is a left-leaning curve. The skewness measures the distribution asymmetry, while kurtosis captures the sharpness of the distribution profile. Figure 8c suggests a more even KD profile than Figure 8d, which verifies that the surface under has a smaller kurtosis than under . Figure 8 concludes again that the surface roughness features can be regulated by controlling the anchoring ratio r with a minimum kurtosis of (Figure 7). These results on roughness parameters can be further applied to tribology [67], topography [68] and analyzing shear mechanics [69].

Figure 8

The kernel density (KD) when (a) (minimum skewness), (b) (maximum skewness), (c) (minimum kurtosis) and (d) (maximum kurtosis).

4. Applications: Pathways to Targeted Surface Roughness Metrics

In practical applications, one would like to know the specific correlations between the most important surfaces roughness metrics and the anchoring coefficient ratio r. Figure 9 summarizes these key results in the parametric anchoring space .

Figure 9

Summary of the parametric sensitivity of the key surface roughness metrics: root mean square, skewness, and kurtosis (Equations (17) and (18)). A straight line with slope passing through the origin in the parametric space corresponds to a morphology (evaluated by -matrix), and its y-axis coordinate (purple line ) determines the amplitude of the wrinkling profile. The location of these r-points in parametric space (Figure 3) are: , zone F; , zone E2; to are in zone E1; and to are in zone D. The density of the special surface roughness metrics in the negative r-value interval indicates that zone E1 in Figure 3 is a highly sensitive area.

Figure 9 summarizes the parametric sensitivity of the surface roughness metrics (root mean square, skewness, kurtosis) as a function of the anchoring coefficient ratio r, shown as straight lines (top) in the plane. The r-ratio interval is from to (from the blue positive slope line to the red negative slope line). The constant r-lines intersect the horinzontal purple line ) at an ordered set of decreasing r-values, corresponding to the local extrema of the surface roughness metrics within the interval. As r decreases, the root mean square keeps declining until it reaches a minimum; then, it climbs up with increasing r. The skewness and kurtosis show a completely different profile. Within the range , the skewness demonstrates two maximum and two minimum values, while the kurtosis presents one maximum and two minimum values. This phenomenon verifies that r plays a complex role in 2D surface wrinkling. Importantly, the local extrema in the surface roughness metrics occur at different values of r, and the presence of several extrema in the negative r-value interval shows that sector E1 of Figure 3 is a highly sensitive parametric area that can be exploited in future applications.

Figure 9 provides a road map to design egg carton surfaces with specific surface roughness metrics. For example, a self-assembled biocompatible cholesteric liquid crystal scaffold with the lowest root mean square can be produced while maintaining the anchoring coefficient ratio . If we want to increase the surface shear strength, we need a surface with high positive skewness [69]. If r is negative, we can choose , and if r is positive, the higher r, the better. The surface skewness and kurtosis are also shown to determine the contact angle and further influence the wettability of materials [90]. As mentioned above, the desired anchoring coefficients can be tuned by controlling the chemical composition and temperature [91].

5. Conclusions

Generating multifunctional wrinkled surfaces usually requires external stress loads, erosion and removal processes, or high-energy beams. In this paper, we present a platform to create nanowrinkling using an intrinsic surface pattern formation mechanism that exists in anisotropic chiral soft matter materials such as cholesteric liquid crystals. Here, shape is determined by internal capillary pressures generated by couplings and interactions between the geometry and the material structure. Specifically, the anisotropic interfacial energy depends on the coupling between the surface unit normal and the average molecular orientation of liquid crystals. Including chirality, spatial periodicity, and basic deformation modes generates additional capillary pressures not found in isotropic materials. The higher the coupling order between the geometry and the orientation, the more complex the surface energy landscape is and the greater the potential for complex surface wrinkling, such as egg carton surfaces, as demonstrated in this paper.

Specifically, in this paper, we studied the complex nanowrinkling patterns of cholesteric liquid crystals in a quartic anchoring model parameterized by two anchoring coefficients . For a given director field, the surface has 10 fundamental wrinkling modes. In a linear analysis of the quartic model (such that is never nonzero), the wrinkling morphology (the contribution of each fundamental wrinkling mode) is determined by the anchoring coefficient ratio , while the wrinkling magnitude is determined by the higher-order anchoring only. The surface wrinkling information is captured by a shape matrix, which reflects the specific role of the director-filed symmetry as well as the anchoring model. This efficient approach can easily be generalized to higher-order wrinkling by increasing the the dimensionality of the shape matrix as well as other director fields.

For functional surfaces that depend on specific surface roughness metrics, we can condense this information by performing a parametric study in terms of the anchoring coefficient ratio r and the parametric anchoring plane. Specifically, the higher-order moment information of surface geometry such as root mean square, skewness and kurtosis are evaluated, presenting complex nonlinear effects when the anchoring coefficient ratio is between and . The linear analysis validates the analytic wrinkling profile as well as the scaling laws.

Taken together, this paper presents a theory and simulation platform for the design of nanowrinkled surfaces with targeted surface roughness metrics, generated by internal capillary pressures, of interest in the development of biomimetic multifunctional surfaces.