Elastocapillary coalescence of plates and pillars

Elastocapillary coalescence of plates and pillars


When a fluid-immersed array of lamellae or filaments that is attached to a substrate is dried, evaporation leads to the formation of menisci on the tips of the plates or pillars that bring them together. Similarly, when hair dries it clumps together due to capillary forces induced by the liquid menisci between the flexible hairs. Building on prior experimental observations, we use a combination of theory and computation to understand the nature of this instability and its evolution in both the two-dimensional and three-dimensional setting of the problem. For the case of lamellae, we explicitly derive the interaction torques based on the relevant physical parameters. A Bloch-wave analysis for our periodic mechanical system captures the critical volume of the liquid and the 2-plate-collapse eigenmode at the onset of instability. We study the evolution of clusters and their arrest using numerical simulations to explain the hierarchical cluster formation and characterize the sensitive dependence of the final structures on the initial perturbations. We then generalize our analysis to treat the problem of pillar collapse in 3D, where the fluid domain is completely connected and the interface is a surface with the uniform mean curvature. Our theory and simulations capture the salient features of both previous experimental observations and our own in terms of the key parameters that can be used to control the kinetics of the process.

liquid meniscus, surface tension, instability, cluster formation

I Introduction

While assembly of complex macromolecular structures on microscopic lengthscales Desiraju1989 () is driven by van der Waals interactions, dispersive forces and chemical interactions between constituents, on mesoscopic length scales of the order of microns to millimeters in the context of colloids and larger particles, other surface forces such as those due to capillarity play an important role Bowden1997 (); Rothemund2000 (); Vella2005 (). When particles interact with each other via capillary forces, the resulting patterns are a function of the size and shape of the constituents, and any constraints on their movement. Capillary coalescence is a natural consequence of this and occurs when free particles aggregate at an interface Grzybowski2000 (); Grzybowski2001 (), and also when extended objects such as filaments and lamellae are brought together by interfacial forces which drive aggregation Pokroy2009 (); Bico2004 (); Kim2006 (); Tanaka1993 (). Sometimes, these systems coarsen indefinitely leading to a single cluster, while at other times elastic deformations eventually arrest the process leading to many finite sized clusters. In a typical experimental setting characterizing the latter, fibers, filaments or lamellae are fully immersed in a liquid which is subsequently evaporated so that capillary forces at the liquid-gas interfaces bring the constituents together. Fibers and lamellae are long and soft objects and can easily bend. Thus a competition between actuating capillarity and resisting elasticity selects the structures formed, which can be used to construct substrates with tunable wetting and adsorption properties Bernardino2010 (); Bernardino2012 (). A variety of experimental systems that fall into this category include millimeter-scaled macroscopic brush hairs Bico2004 (), micrometer-scaled mesoscopic polymeric surface mimicking gecko foot hairs Geim2003 (), as well as nanometer-scaled carbon nanotube forests Lau2003 (); Chakrapani2004 ().

Elastocapillary interaction have been well characterized in two-body systems Kim2006 (); Duprat2011 (); Taroni2012 (). However, for the collective behavior of many elastic fibers bundling together, typically static energy minimization arguments have been employed to estimate the expected finally assembled bundle size Bico2004 (); Chandra2009 (); Zhao2006 (), although recent work Gat2013 () uses stability analysis of the unpatterned base state to predict the ordered hierarchical coalescence structure of an array of clamped parallel elastic sheets which are partially immersed in liquid. Complementing these discrete approaches, a phenomenological continuum mean-field approach has been used to model the arrested coarsening Boudaoud2007 (), while a recent continuum theory based on microscopic physics for the coupled dynamics of drying and coalescence explains the kinetics and refinement seen in experimental studies Wei2014 (). Unlike in ergodic systems where the state space is stochastically sampled due to thermal excitations, the structures in elastocapillary systems are often not selected by energetics alone. Instead, they depend critically on the dynamics of the drying process. This leads to a path dependence caused by the strong coupling of the geometry of the air-liquid interface to the local evaporation when multiple unconnected liquid domains are formed. Additional aspects, such as pinning and contact angle hysteresis, as well as the permanent adhesion of contacts formed in intermediate structures, make purely energetic arguments unable to have any predicative power. Thus, to predict and control the assembled structures in capillarity-driven self-assembly experiments, we need to (1) follow the dynamics over time as irreversible effects associated with evaporation, contact line motion and adhesion, and (2) account for non-linear effects due to large deviations from the background state and go beyond a linearized analysis.

Here we use a combination of theory and computation to understand the instability and the hierarchical evolution of cluster formation for both two- and three-dimensional elastocapillary systems driven by drying. Those systems represent generic situations, allow for a thorough theoretical treatment and can be validated by well controlled experiments. For the two dimensional case, we describe the drying induced collapse of an one-dimensional array of evenly spaced lamellae immersed in the evaporating liquid (Figure 1(a)(b)). For this system all relevant physical forces are considered, allowing us to derive the interaction potentials and forces, analyze the linear stability of the system, and compute the nonlinear dynamics associated with pattern coarsening. We demonstrate that different dynamical paths through the system’s state space indeed lead to different final structures which are observed in experiments.

In the three-dimensional case, our theory focuses on explaining experiments associated with the bundling of a regular square grid of fluid-immersed elastic posts anchored to a substrate Pokroy2009 () (Figure 1(c)-(e)). We determine the constant mean curvature surface for the air-liquid interface subjected to the global liquid volume constraint on the multi-connected domain. This allows us to compute the interaction forces and thence the primary unstable mode. Finally we use a numerical method to compute the aggregation dynamics.

In II we describe the experimental observations for the two-dimensional case, derive a discrete two-plate model for the deformation of the plates driven by capillary forces, and carry out a linear stability analysis of the base state. We then study the coalescence dynamics of a collection of plates in the nonlinear regime. In III we describe the experimental observations for the three-dimensional case, and introduce a dynamical model that allows us to simulate the morphology of clustered pillars, and compare these results with experiments. In IV, we conclude with a description of open problems in this rich area.

Figure 1: Dynamics of coalescence of arrays of (a-b) lamellae and (c-e) filaments driven by drying. (a) The lamellae are of thickness , depth and height , and spaced a distance apart. As the liquid evaporates, a front of dimerization coalescence initiates and propagates from right to left, after which a front of quadrimerization moves through the system in the same direction. After the liquid dries out, the quad-bundles separate and the dimers persist. (b) The lamellae are of depth , and other geometric parameters are the same as those in (a). Due to the imperfections in the system, e.g. the lamella geometry, roughness of the surface and et al., coalescence initiates simultaneously here. The 2-lamella mode still appears to be dominant at the onset of instability, and irregular bundles with a size distribution from 2-5 arise thereafter and persist even after the liquid dries out. (c) The final structures of the two-dimensional array of pillars of diameter after the liquid dries out form four-fold rhombic clusters. (d) Here the tips of the fourfold clusters form squares. (e) A larger domain shows that the four fold clusters are more asymmetric, and the tips form both rhombi and squares. The geometric parameters of the pillars and the material properties are the same in all three images; the pillars have a diameter , height and spacing . The only difference is the size of imperfections in the system.

Ii Collective dynamics of elastic lamellae

ii.1 Experimental observations

For the two-dimensional case (Figure 2), we consider a one-dimensional array of elastic micro-lamellae with height , thickness and uniform spacing , Young’s modulus and Poisson’s ratio respectively. Each lamella is assumed to be free at one end and anchored at the other on a substrate Tanaka1993 (). The lamella array is wetted by a liquid of surface tension , density and viscosity , which is confined between neighboring lamellae, defining a cell. The contact line slips from the tips as the liquid evaporates. When the system is completely immersed in the liquid, the stable configuration is a uniform array of non-interacting vertical lamellae. However, when the liquid evaporates, it is not necessarily locally stable any more: capillary forces associated with the liquid-air menisci between the free ends of the soft lamellae may cause them to deflect laterally and adhere together. In an experimental system with small imperfections, we observe a regular cascade of successive sticking events that leads to a hierarchical bundling pattern: every two neighboring lamellae incline towards each other to form a dimer first, which then collapses into quadrimers (Figure. 1(a)). The process repeats until the bending deformation induced elastic resistance eventually becomes large enough to prevent further coarsening. In the system with large imperfections, irregular bundles can and do arise but the 2-lamella-collapse mode still appears to be dominant right after the instability(Figure. 1(b)). After the liquid dries out, bundles separate if the adhesion in contact is not strong to counterbalance the elastic forces; else they persist. In our experiments, the liquid used is isopropyl alcohol (IPA), and the lamellae and the substrate are made of polydimethylsiloxane (PDMS). Throughout this entire section, we use the following experimental parameters , , , , , , , , and gravity .

ii.2 Mechanics: coupling plate bending to fluid interface shape

In our experiments, the lamellae are short and stiff, and remain almost straight as they are deflected by capillary forces, bending primarily in the neighborhood of the base. Therefore, we can approximate each lamella as a rigid plate and integrate all the bending response into an elastic hinge at the base Bernardino2012 (). This simplifies our analysis relative to the case that must account for the inhomogeneous bending and buckling (APPENDIX A) of the individual lamella Gat2013 (); Neukirch2007 (). The hinge elastic constant can be approximately derived from the bending response of a short cantilever by a transverse force at its free end, which is given by


where is the deflection at the free end, is the tilting angle between the straight plate and the horizontal direction (Figure 2(a)), is the shear modulus, and is the shear coefficient. The second term of the right hand side of Eq.(1) is due to shear deformations in the so-called Timoshenko beam theory Timoshenko1972 (), because the slender beam condition is violated in our experimental setup. We have taken and as an approximation for the rectangular cross section.

Figure 2: A unit cell confined by two adjacent lamellae/plates. (a) Each lamella is modeled as a rigid plate elastically hinged at the base. The plate is of height and thickness . is the angle of plate with respect to the horizontal. The plate spacing is and the distance between the two free tips is . is the volume per unit depth of the liquid confined in the cell. The red curve represents the air-liquid interface, an arc of a circle of radius , with being the half angle subtended by the meniscus arc. is the critical contact angle at which the meniscus slides down from the tip. is the wetting length on the left side of the plate and is the wetting length on the right side of the plate. (b) 1-8 are the 8 possible cases of menisci, showing that the menisci can be pinned on both tips, or slip down from one or both tips. 7 and 8 show that when the two free tips are so close that no circular arcs exist to connect them, in which case the arc is replaced by a line. These situations are solely used to prevent simulation failure in rare cases.

To explicitly derive the torques due to capillarity, we consider a unit cell consisting of two plates, with liquid confined in between, and air outside (Figure 2(a)). The pressure field inside the liquid is nonuniform due to the effects of gravity and flow. However, in our system, gravity can be neglected, because the Bond number . Comparing the viscous moment due to flow (APPENDIX C) with that caused by the surface tension , we find that for , which indicates that pressure effects due to fluid flow can be neglected. Therefore, the air-liquid interface is always a segment of a circle because of the uniform pressure in each cell. For each plate, the moment results from both the line tension at the contact line, and the pressure gradient across the plate caused by the curvature difference of two neighboring menisci. To calculate the moment we need to find the wetting length, the contact angle, and the pressure difference across the air-liquid interface due to its local curvature. We assume the contact angle can take any value equal to or larger than a critical value , consistent with the fact that the meniscus can be pinned at the tip of the lamella while being concave-up (Figure 2(b1)), flat (Figure 2(b2)), or concave-down (Figure 2(b3)). We assume that once the contact angle reaches the critical value , it remains constant as the meniscus starts to slide down from the tip (Figure 2(b4)-(b6)). Consequently, the meniscus profile and the resulting moment can be calculated for any given tilting angles , and the liquid volume per unit depth . Scaling lengths by , volumes by and moments by leads to the following results for the different cases that we summarize here for our subsequent linear stability analysis (for details, see APPENDIX B).

  1. The meniscus is pinned on both tips (Figure 2(b1)). The half angle subtended by the meniscus arc is determined by solving


    for given , and , where and tip separation is


    must satisfy , where is the critical angle at which the meniscus starts to slide down from at least one lamella. when the meniscus concaves down, when the meniscus is flat, and when the meniscus concaves up. The moments on the and plates are given respectively by

  2. The meniscus has slipped down from both tips (Figure 2(b3)). The contact angle is fixed at . When , the meniscus radius is independent of , and . is determined by solving


    The wetting length on the right side of the plate and that on the left side of the plate are given respectively by


    The moments on the and plate are given respectively by


To get the total moment on a plate, we must add the contributions from adjacent cells. For example, we can obtain the full expression of simply by replacing by in Eq. (5) and add up to Eq. (4) for the case when the meniscus is pinned at both tips, which readily yields .

The dynamics of the plate neglecting inertia (APPENDIX C) follows the overdamped first order equation of motion


where is the damping coefficient, is defined in Eq. (1) and the dimensionless moment is due to capillarity. (We use Eqs. (B3) and (B4), Eqs. (B11) and (B12), Eqs. (B17) and (B18), Eqs. (B23)-(B26), and Eqs. (B29) and (B30) for different situations to obtain the full expressions of as explained above. Please see APPENDIX B for more details.) To estimate , we need to account for both the internal viscosity of the solid and the external viscosity of the fluid, and find that the former dominates, which gives (APPENDIX C). Together with Eq. (11) and the dynamics of drying that will be discussed later, we can now understand the dynamics of the lamella array completely.

ii.3 Onset of bundling: linear stability of the uniform base state

For the base state with all lamella being vertical and a uniform meniscus associated with constant liquid volume in each cell is decreased, there is a potential for instability as the curvature of the menisci increase.

To understand this, we consider a periodic domain of plates and the same volume of liquid confined in each cell, when all plates being vertical () is an equilibrium state. To determine the stability of this state, we study the perturbations in the moment as a function of variations in the angles linearized around the current state,


where and are the stiffness matrices due to the elasticity of the plate and the geometrical change of menisci respectively. is a diagonal matrix with all elements being the dimensionless hinge constant , and is a tridiagonal matrix with two additional elements of on the upper-right and lower-left corners reflecting the periodic boundary condition. At the base state with , is expressed as


where and are the stiffness of the effective spring connecting two neighboring plates due to capillarity, where is defined in Eq. (4) or Eq. (9), and is defined in Eq. (5) or Eq. (10). For the case when the meniscus is pinned at both tips,


where is determined by the volume


which follows from Eq. (2) and Eq. (3) by substituting . Here we have omitted the subscript for the translationally invariant base state. Note that and in this case are not necessarily equal, because while the wetting length is constant, the contact angles change by different amount on the two neighboring plates as they are deflected except when the meniscus is flat. For the case when the meniscus is no longer at both tips,


with is the radius of the meniscus, and the wetting length is determined by the volume


which follows from Eqs. (6)-(8) by substituting . and in this case are always identical, because the contact angle keeps constant, and the change of wetting length is the same on both plates when they are deflected. From Eqs. (14)-(17), we see that for any given geometric parameter and critical contact angle , (). The inset of Figure 3(a) shows that switches sign as decreases, corresponding to the case when the system becomes unstable. We note that the discontinuity in corresponds to the transition from a static meniscus to a dynamic one that slides down from the tip.

Figure 3: Stability analysis of a vertical array of plates. (a) The curve of following Eqs. (14) and (16) as a function of given by Eqs. (15) and (17) intersects with the line . The interval between the two intersections indicates the region where the system is unstable as according to Eq. (18), . For typical experiments, any leads to coalesce given infinitesimal perturbations to the system. Inset figure shows the effective dimensionless spring constants (); switches signs from positive to negative as decreases, while the discontinuity corresponds to the meniscus sliding down from the tip. (b) shows the positive bifurcation branch using both simulations of Eq. (11) and the asymptotic result in Eq. (20). When , , and when , , which is a clear indication of supercritical bifurcation. The inset illustrates the dimer mode. We simulate 2 plates by solving Eq. (11) with periodic boundary conditions, and display them for visualization. The black dots correspond to stable dimers, while the red dots correspond to unstable dimers that will further coalescence.

To study the instability of the system, we investigate the eigenmodes of the stiffness matrix which resembles a discrete Laplacian, with eigenvectors


and the corresponding eigenvalues


must be positive definite to ensure stability, which is equivalent to requiring the smallest eigenvalue positive. When , the meniscus is pinned at both tips, and the smallest eigenvalue is , which follows from Eqs. (14) and (19). In our experiments, , so that the array of vertical plates is stable. When , the smallest eigenvalue is , so that stability is controlled by the competition between elasticity and capillarity. As (), the condition sets the range of in which the system is unstable (Figure 3(a)). The primary eigenmode corresponds to , and the eigenvector thus is from Eq. (18), corresponding to dimerization of the lamellae. This calculation shows that to avoid the lamella collapse, we need to keep or reduce when and increase . Practical approaches to implement this include the use of liquid with low surface tension and contact angle close to , at which is at its maximum positive value, the use of stiff solids, and/or a proper choice of geometric parameters, e.g. a large aspect ratio .

To explain the nature of the instability transition, we make use of the fact that the fastest growing mode is the dimer mode and assume . In the vicinity of the critical volume at which the instability happens, the dynamics of a lamella is governed by the equation


where is the dimensionless damping coefficient, and is an algebraically lengthy coefficient of the cubic term. Eq. (20) is derived by expanding Eq. (11) in powers of ; even orders of do not appear in Eq. (20) due to the reflection symmetry inherent in the system. When the volume in a cell reaches the critical value (equivalently ), , while when , a Taylor series expansion of Eq. (14) in the neighborhood of gives us , while a similar expansion of Eq. (15) in the neighborhood of yields , and hence . From Eq. (20), we see that the stable equilibrium state has two branches of solutions for and only one solution otherwise, suggesting that the bifurcation is supercritical. Figure 3(b) shows both the positive branch of the asymptotic solution , and the results of simulation for lamellae obtained by solving Eq. (11) with periodic boundary conditions. The results agree well, and confirm that the instability is supercritical and leads to lamellar dimerization, as shown in the inset.

Figure 4: Elastocapillary coalescence of plates. (a) A periodic domain of 32 plates is simulated by solving Eq. (11) for different prescribed volumes . The tilting angle of the first plate is perturbed by from as the initial condition. The left column shows a sequence of transient states and the right column shows the steady states. (b) The standard deviation of deflection angle as a function of time. is time and is the time scale for mechanical relaxation. (c) Complementing (b), these plots show the evolution of hierarchical formation of bundles for different prescribed volumes. The time resolution of the plots is and plates are considered as a bundle if the tilting angle gradient is positive and larger than the perturbation amplitude.

ii.4 Nonlinear dynamics: drying, coarsening and refining

Controlled liquid volume

For an initially translationally invariant system, the two-plate-collapse mode is the fastest growing mode, yet the cluster of dimers is not necessarily the final stable state (red dots in Figure 3(b)). We notice that for a range of sufficient liquid volumes in each cell, both the deflection of the plates and the geometric change of the liquid menisci are large so that the linear approximation in Section II(II.3) breaks down. Therefore, we numerically integrate Eq. (11) directly to follow the hierarchical dynamics by which an array of vertical individual plates first forms dimers, then quadrimers, until eventually forming large bundles that are limited by elastic effects. Figure 4(a) shows snapshots of dynamical coarsening for different control parameters , in a periodic domain of 32 plates, triggered by tilting one plate by from as the initial condition. As expected, for volumes outside the range , e.g. for the cases of and , the array is stable to perturbations and remains uniformly vertical. In the unstable parameter range, the primary mode corresponds to two plates collapsing into dimers (Figure 3(a)(c)). These dimers may further collapse into quadrimers or stay as the final stable state with different amplitudes of deformation angles as shown in Figure 4(a)(b). For the same initial perturbation, the dynamical path of successive bundle aggregation depends on the control parameter . As examples, we see that for , quadrimers sweep through the domain right after the dimers form and the system reaches the stable equilibrium, for the dimers persist for a while before they eventually collapse to quadrimers, and for the dimers stay stable, as shown in Figure 4(c).

Having seen how the system coarsens when the volume is controlled, we now turn to the more realistic case when the volume itself evolves dynamically.

Coupling to dynamics of drying

For an evaporation dominated situation, the rate at which the liquid volume in each cell is reduced depends on the local surface area of the air-liquid interface and thus on deflection angles of the adjacent lamellae. Consequently, the drying dynamics is coupled to the evolution of the geometric configuration, and new instabilities associated with inhomogeneous cell volumes are expected, which cannot be captured by either energy minimization Bico2004 (); Chandra2009 (); Zhao2006 () or renormalization analysis Gat2013 ().

Figure 5: Drying induced elastocapillary coalescence of plates assuming evaporation rates to depend on the surface area in each cell. A periodic domain of 100 plates is simulated by solving the coupled Eqs. (11) and (21), of which only 32 are displayed. (a) Snapshots of the array of plates and menisci at different times. The initial conditions are except that is smaller by and . (b) Number count of bundles of different size as a function of time for (a). (c) Same as (a) but with uniform random initial perturbations with the maximum amplitude of on both and . (d) Number count of bundles of different size as a function of time for (c). Note the 3- and 5-plate bundles.

A minimal evaporation model that is sufficient to capture the qualitative features of hierarchical bundle formation is given by


where is a constant, and are time scales for mechanical relaxation and evaporation respectively, is the radius of the meniscus and is the angle subtended by the meniscus arc. Eq. (11) and Eq. (21) coupled together determine the dynamics of lamella coalescence driven by evaporating liquid, with three regimes. When , decreases quasi-statically so that Eq. (20) relaxes to a static state for a prescribed value of . When , the evolution of the cell volumes and the lamellar configuration are coupled. When , the evaporation is so fast that the lamella array does not coarsen. In our typical experiments, and is of the order of seconds, so we choose correspondingly in Eq. (21) so that the liquid dries out in about seconds. Figure 5 shows the simulation results for a periodic domain of 100 plates with 2 different initial perturbations.

In Figure 5(a), all plates are perfectly vertical, and the liquid volume in each cell is constant for all cells except that is smaller by to mimic boundary effects in an experimental system. The evolution of the system indicates that a front of dimer coalescence propagates from the imperfection site and sweeps through the entire domain, followed by a successive front of quadrimer coalescence. The largest transient bundles have 4 lamellae. If adhesion between contacting lamella is neglected, capillary forces are not sufficiently large to hold the plates together, and the bundles separate symmetrically when the liquid volume falls bellow a second threshold, and a perfectly vertical configuration is restored. Figure 5(b) shows the number of bundles of different sizes as a function of time corresponding to the configuration shown in Figure 5(a), and highlights the fact that bundle formation/separation is perfectly hierarchical and regular.

In Figure 5(c), uniform random perturbations with maximum relative amplitude of are applied to all the tilting angles and liquid volume in all cells. Irregular bundles of size ranging from 2 to 5 form transiently, and separate as the liquid evaporates. Figure 5(d) presents the number count of bundles of different size for Figure 5(c), which shows that although the dimer is still the dominant mode in the early stages, bundle formation can be irregular. This is because the the jump in the location of the contact line when the contact angle reaches a critical value leads to a sudden decrease in the effective stiffness as shown in Figure 3(a), and the uniform random initial perturbation generates multiple sites from which the front of dimer coalescence starts propagating. Consequently the sites are not necessarily separated by an even number of plates. Therefore, trimers and pentamers also arise in addition to quadrimers. Bundles of larger size do not appear because of the large elastic energy associated with their formation.

The dynamics of coalescence in these two cases shows how the number of plates per cluster varies in a step-like manner, very similar to the experimental data reported by Pokroy Pokroy2009 () and Gat Gat2013 (). All these features also agree qualitatively well with our own experimental observations shown in Figure 1(a)(b).

Iii Collective dynamics of a two-dimensional array of pillars

iii.1 Experimental observations

We now generalize our study of the one dimensional dynamics of plate or lamella aggregation driven by capillarity to the coalescence of a two-dimensional array of epoxy nano-pillars immersed in an evaporating wetting liquid as reported in detail in previous work Pokroy2009 (); Kang2010 (). The dynamics of pillars is different from that of plates in 3 major ways. First, fluid can flow freely around the multi-connected domains associated with pillars, so that the interaction between them occurs over much longer ranges rather than being limited to just nearest neighbors. Secondly, the three-dimensional geometry allows the pillars to bend in two principal directions and also twist. For pillars with a circular cross-section, the twist must be a constant. If we neglect friction between adhering filaments, the twist must identically vanish, and here we will assume that this is the case. Thirdly, experimentally we see that a segmented, “wormlike” geometry of the specially treated pillars increases pinning of the receding contact line by reentrant curvature Pokroy2009 (); here will neglect this effect for simplicity.

As in the case of lamellae, the uniform array of non-interacting straight pillars loses stability as the liquid evaporates. The dynamics of the ensuing structures is a result of the competition between elasticity and capillarity, and the morphology of the final assembly is determined by intrapillar elasticity and interpillar adhesion Kang2010 (). Figure 1(c)-(e) show the scanning electron microscopy (SEM) images of the assembly into uniform periodic fourfold clusters of nanopillars, in which the pillar height , the pillar radius , the pillar spacing , the Young’s modulus , the surface tension of the liquid , and the density of the liquid . Unlike in the one-dimensional array of lamellae, where the dimer is the primary unstable mode, for pillars the quadrimer is the primary unstable mode. As the liquid evaporates, this mode gives way to hierarchically grow into larger assemblies which eventually get arrested by the increase in the elastic resistance.

iii.2 Mechanics: coupling filament deformation to fluid interface shape

As in the lamellar case, inertial effects can be neglected here, so that the dynamics of each pillar tip can be characterized by its displacement vector relative to its base (Figure 6(a)), and is given by


where is the drag coefficient, is the elastic bending resistance force at the tip, is the capillary driving force due to surface tension , and is the liquid volume in the system. The dominant contribution to is from the internal damping of the viscoelastic solid similar to the lamellar case, and (APPENDIX D), where is the time scale for the fiber to relax mechanically.

To compute the bending resistance force in the horizontal direction, we use the theory of the elastica for the inextensional, unshearable deformation of thin filament. Letting be the angle of the pillar centerline tangent with the vertical direction, with is the arc length coordinate, and the force amplitude at the tip, equilibrium implies that


Geometry implies that , so given , is uniquely determined, and .

To obtain , we need to determine the shape of the air-liquid interface. Since the Bond number , gravity can be neglected. Moreover, the time scale for the fluid to equilibrate in the porous brush is much smaller than that for the pillars to relax mechanically , which is much smaller than that for the evaporation , i.e. (APPENDIX E). Therefore, the pressure throughout the liquid domain can be regarded as uniform and the air-liquid interface is thus a surface of uniform mean curvature and satisfies the equation


where is the mean curvature of the interfacial surface. Without loss of generality, we set the ambient pressure to zero, and let be the pressure inside the liquid. Volume conservation in the whole domain yields


and serves to determine . Here is the projected domain of the air-liquid interface to the horizontal plane (meshed area in Figure 6(a)), and is the elevation of pillar tip. Since the menisci are always pinned on the pillar tips, we need to solve Eq. (24) on a multiply connected domain, in which pillar tips are regarded as solid circles with the identical radius (Figure 6(a)) and the height of the surface is fixed at the elevation of pillar tips (Figure 6(b)) calculated from Eq. (23). As the pillars are effectively immersed in liquid, the integration of pressure over the lateral surface of the cylinder does not contribute to , and the only active contribution is the line tension at the contact line. For a given air-liquid interface (Figure 6(c)), the angle of the meniscus tangent with the horizontal direction on the circular boundary of the tip is known, which we denote as . To calculate the capillary driving force on each pillar, we integrate the interfacial force over the contact line contour at the tip and determine the component in the horizontal direction that contributes to the deflection, with , where the subscript represents the tip circle and is its unit outward normal. Note that for very small values of , the assumption of the surface being pinned at the pillar tips breaks down; however clusters form well before this assumption is violated, so that this is not a cause for concern.

To prevent penetration upon collision between pillars, we treat each pillar as a rod with finite radius with an artificial short-range repulsion force when the discs representing pillar tips come close enough ( of the pillar diameter), but do not consider the elastic deformation of the cross section due to contact.

Figure 6: Elastocapillary coalescence of a square array of pillars. (a) Demonstration of triangular meshes on a domain that contains an array of 2 by 2 vertical pillars, where the white solid circles represents pillar tips. The mesh density used in the actual simulation is 4 times denser. (b) The three-dimensional air-liquid interface is obtained by solving Eq. (24) on the domain shown in (a) for a given liquid volume in Eq. (25). (c) Given a slight perturbation to the vertical state in (b), 4 pillars coalesce to form a bundle. (d)-(g) A domain of 14 by 14 pillars evolves to the steady state of fourfold clusters for a given liquid volume . The dynamics follows the coupled evolution equations (22)-(25). is the dimensionless time scaled by (see text). The gray scale shows the air-liquid interface height, scaled by the pillar spacing. White solid circles represent pillar tips, and the black open circles represent pillar bases. The dashed lines connecting bases and tips correspond to a projection of pillars viewed from the top.
Figure 7: More examples of final structures of the fourfold clusters for different prescribed volumes and domain sizes with random initial perturbations and symmetric boundary conditions. (a) 14 by 14 pillars with a given liquid volume . (b) 10 by 16 pillars with a given liquid volume . Simulations were done by solving equations (22)-(25).

iii.3 Nonlinear dynamics: simulations and comparison with experiments

To complete the formulation of the problem, we need some boundary conditions. We consider a square array of pillars inside a domain with straight vertical walls and impose symmetry-related conditions on the walls for the fluid interface and the pillar deformations. On the pillars, the contact lines are assumed to be pinned on all pillar tips. For initial conditions, the pillar bases (open circles in Figure 6(g)) are assumed to form a perfect periodic square lattice, while the pillar tips (solid white circles in Figure 6(g)) are perturbed from the vertical configuration so that the layer of pillars closest to the boundaries inclines inwards to trigger the inward motion from the boundaries, and other pillar tips are uniformly and randomly perturbed. We numerically solve the coupled evolution equations (3.1)-(3.4) using a custom-coded finite element code. For a given pillar tip displacement , Eq. (23) is solved with the geometric condition to determine the reaction force and tip elevation of each pillar . Given and , Eq. (24) is solved using the finite element method to determine the interface on the multiply-connected domain with the liquid volume constraint Eq. (25), which determines and thence . Then Eq. (22) is integrated in time explicitly to update the pillar tip positions, and the domain is remeshed accordingly at every time step.

Figure 6(d)-(g) show the simulated evolution of an array of by pillars collapsing into fourfold bundles for a prescribed liquid volume of , where is the control volume inside the system when all the pillars are vertical and the air-liquid interface is flat, using the parameters from experiments in Figure 1(c)-(e). Due to the initial boundary perturbations, coalescence is initiated from the boundary and propagates towards the center of the domain. We observe that there is a critical liquid volume above which uniform vertical pillars are stable, and below which pillar coalescence occurs, similar to the case of lamellar collapse. However, here the primary eigenmode of instability has a fourfold symmetry (Figure 6(e)) independent of initial perturbations as observed in experiments (Figure 6(g) and Figure 7). An intuitive way to understand this is to recognize that the fourfold bundles have two principal directions, along which one sees dimers. Although the interaction potential between pillars in the three-dimensional case has a much longer range than in the two-dimensional case because of connectivity, it is still monotonic and decreases with distance. Provided that within the linear analysis the effective spring constants in the two principal directions are decoupled, the dimer is the primary mode in each direction as in the one-dimensional lamella array. Beyond the linear regime, our simulations capture the coarsening that is arrested and eventually leads to a maximum bundle size. We also note that the tips in a bundle can form either rhombi or squares depending on the initial perturbations, although rhombi are more likely as they are stable against shear deformations; indeed as we neglect friction between the tips, the rhombus is a more energetically favorable configuration than the square. However, the energy difference between these two states is very small, so that contact and friction in real experiments leads to both square tips (Figure 1(d)) and rhombi tips (Figure 1(c)).

Iv Conclusions

Our study has focused on understanding the onset and evolution of elasto-capillary coalescence of plates and pillars driven by evaporation. For the case of lamellar collapse, we explicitly derived the conditions for the primary dimerization instability in terms of the state variables - tilting angles and liquid volumes , and the relevant geometrical and physical parameters. Complementing our analysis, full numerical simulations show that the final coalescent states sensitively depend on initial perturbations because of the discontinuous motion of the contact line when the contact angle reaches a critical value. This implies that the self-organization of clusters cannot be predicted by energy minimization arguments alone, but depend on the dynamics of the drying process - this is especially true when the coupling of geometry to local evaporation rates is taken into account. Our model accounts for this, and is in qualitative agreement with experimental observations of the intermittence of coalescent transitions. For the case of pillar collapse, our model correctly accounts for the multiply-connected nature of the fluid interface, and the large elastic deflections of the pillars. The analysis based on this model captures the primary fourfold eigenmode associated with the onset of collapse, consistent with experimental observations. Numerical simulations of the full dynamics allow us to follow the evolution of the clusters whose eventual size is determined by the competition between capillarity and elasticity. For both cases, our numerical results agree well with many of the salient experimental observations. In particular, we can explain the eigenmodes at the onset of instability, and the time scales on which clusters form, while providing explanations for both regular and irregular hierarchical bundling till the final state.

However, our analysis still leaves out some effects and thus cannot explain some observations. Neglecting adhesion and friction between pillars implies that we cannot explain the twisting of pillars that leads to the formation of chiral clusters often seen. This is a natural next step in the analysis. Furthermore, we have limited ourselves to a discrete theory in both cases, but in the thermodynamic limit of a large number of pillars or lamellae, one might ask what the nature of a continuum theory might be. A recent continuum theory that addresses the explicit connection of the essential geometric and physical parameters to determine the maximal size and dynamics of the assembly has been carried out for lamellar coalescence Wei2014 (), but the question for 3-dimensional coalescence remains an open question.


We thank Sam Ocko for many discussions that helped to sharpen and clarify our arguments. We thank the Harvard-MRSEC DMR -0820484, the MacArthur Foundation (LM), and the NRF of Korea (Grant No. 2013034978, H.-Y.K.) for support.

Appendix A Critical buckling length of a plate

The critical buckling height of a thin plate under compression due to surface tension is


where is the typical elastocapillary length scale. For the case when the plate is clamped at one end and axially compressed at the other, the exact expression Timoshenko1961 () for is


where we have substituted in the experimentally observed parameter values. As is much larger than the lamella height , buckling of the lamellae when they pierce the gas-liquid interface does not happen in our system.

Appendix B Moment calculation for the 2D case

The 8 possible meniscus configurations can be classified into 6 cases, where we explicitly express the moments in terms of the state variables - tilting angles and liquid volumes . We have scaled all lengths by and moment by , so all results below are dimensionless.

Case 1: The meniscus is pinned on both tips (Figure 2(b1)). The half angle subtended by the meniscus arc is determined by solving


for given , and , where and half of the tip distance is


must satisfy , where is the critical angle at which the meniscus starts to slide down from at least one lamella. when the meniscus concaves down, when the meniscus is flat, and when the meniscus concaves up. The moments on the and plates are given respectively by


Case 2: The meniscus is down from both tips (Figure 2(b3)). The contact angle is fixed at . When , the radius of the meniscus is determined by solving


The wetting length on the right side of the plate and that on the left side of the are given respectively by


When , the meniscus radius is independent of , and . is determined by solving


The wetting length on the right side of the plate and that on the left side of the plate are given respectively by


In either case, the moment on the plate and that on the plate are given respectively by


Case 3: when the meniscus slides down from the plate and is pinned on the one. In this case and .


From the condition that contact angle on the plate is , we can get the following relations,




and are determined by solving either (41) and (42) or (41) and (43).

Moment on the plate and that on the plate are given respectively by


Case 4: when the meniscus slides down from the plate and is pinned on the one. In this case and .


From the condition that contact angle on the plate is , we can get the following relations,




Similarly, and are determined by solving either (47) and (48) or (47) and (49).

Moment on the plate and that on the plate are given respectively by


For certain given , and , although one end of the meniscus is depinned from the tip, there is no arc satisfying the enforced contact angle condition. In the following two cases, the menisci are approximated by straight lines. They are used to prevent the numerics from blowing up when two plates almost contact, yet very unlikely to happen.

Case 5: when the meniscus slides down from the tip of the plate and keeps flat. The moment on the plate and that on the plate are given respectively by


where is determined by solving from the given volume


and is


Case 6: when the meniscus slides down from the tip of the plate and keeps flat. The moment on the plate and that on the plate are given respectively by


where is determined by solving from the given volume


and is


Appendix C Damping coefficient for the 2D case of plates

In order to calculate the damping coefficient, we need to account for both the viscosity of the fluid and the viscoelasticity of the solid. First we consider the contribution from the fluid. The Reynolds number or less, so that inertia of the fluid is negligible. We use lubrication theory to calculate the moment acting on the plate caused by flow although the ratio of plate height to the spacing does not strictly satisfy that . Figure 8 illustrates 3 rigid plates hinged at the base. The upper and lower one are perpendicular to the substrate, and the middle one is rotating clockwise at the angular velocity . The liquid is confined between two plates and the ambient pressure is set to be 0.

In the bottom chamber, the momentum conservation of the fluid in the and directions are


where is the pressure and is the fluid viscosity. The boundary conditions are


From (60) and (61), we can get