Elastocapillary coalescence of plates and pillars
Abstract
When a fluidimmersed 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 twodimensional and threedimensional setting of the problem. For the case of lamellae, we explicitly derive the interaction torques based on the relevant physical parameters. A Blochwave analysis for our periodic mechanical system captures the critical volume of the liquid and the 2platecollapse 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.
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 liquidgas 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 millimeterscaled macroscopic brush hairs Bico2004 (), micrometerscaled mesoscopic polymeric surface mimicking gecko foot hairs Geim2003 (), as well as nanometerscaled carbon nanotube forests Lau2003 (); Chakrapani2004 ().
Elastocapillary interaction have been well characterized in twobody 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 meanfield 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 airliquid 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 capillaritydriven selfassembly 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 nonlinear 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 threedimensional 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 onedimensional 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 threedimensional case, our theory focuses on explaining experiments associated with the bundling of a regular square grid of fluidimmersed elastic posts anchored to a substrate Pokroy2009 () (Figure 1(c)(e)). We determine the constant mean curvature surface for the airliquid interface subjected to the global liquid volume constraint on the multiconnected 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 twodimensional case, derive a discrete twoplate 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 threedimensional 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.
Ii Collective dynamics of elastic lamellae
ii.1 Experimental observations
For the twodimensional case (Figure 2), we consider a onedimensional array of elastic microlamellae 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 noninteracting vertical lamellae. However, when the liquid evaporates, it is not necessarily locally stable any more: capillary forces associated with the liquidair 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 2lamellacollapse 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
(1) 
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 socalled 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.
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 airliquid 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 airliquid 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 concaveup (Figure 2(b1)), flat (Figure 2(b2)), or concavedown (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).

The meniscus is pinned on both tips (Figure 2(b1)). The half angle subtended by the meniscus arc is determined by solving
(2) for given , and , where and tip separation is
(3) 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
(4) (5) 
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
(6) The wetting length on the right side of the plate and that on the left side of the plate are given respectively by
(7) (8) The moments on the and plate are given respectively by
(9) (10)
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
(11) 
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,
(12) 
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 upperright and lowerleft corners reflecting the periodic boundary condition. At the base state with , is expressed as
(13) 
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,
(14) 
where is determined by the volume
(15) 
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,
(16) 
with is the radius of the meniscus, and the wetting length is determined by the volume
(17) 
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.
To study the instability of the system, we investigate the eigenmodes of the stiffness matrix which resembles a discrete Laplacian, with eigenvectors
(18) 
and the corresponding eigenvalues
(19) 
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
(20) 
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.
ii.4 Nonlinear dynamics: drying, coarsening and refining
Controlled liquid volume
For an initially translationally invariant system, the twoplatecollapse 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 airliquid 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 ().
A minimal evaporation model that is sufficient to capture the qualitative features of hierarchical bundle formation is given by
(21) 
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 quasistatically 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 steplike 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 twodimensional 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 twodimensional array of epoxy nanopillars 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 multiconnected 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 threedimensional geometry allows the pillars to bend in two principal directions and also twist. For pillars with a circular crosssection, 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 noninteracting 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 onedimensional 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
(22) 
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
(23) 
Geometry implies that , so given , is uniquely determined, and .
To obtain , we need to determine the shape of the airliquid 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 airliquid interface is thus a surface of uniform mean curvature and satisfies the equation
(24) 
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
(25) 
and serves to determine . Here is the projected domain of the airliquid 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 airliquid 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 shortrange 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.
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 symmetryrelated 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 customcoded 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 multiplyconnected 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 airliquid 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 threedimensional case has a much longer range than in the twodimensional 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 onedimensional 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 elastocapillary 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 selforganization 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 multiplyconnected 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 3dimensional coalescence remains an open question.
Acknowledgment
We thank Sam Ocko for many discussions that helped to sharpen and clarify our arguments. We thank the HarvardMRSEC 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
(26) 
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
(27) 
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 gasliquid 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
(28) 
for given , and , where and half of the tip distance is
(29) 
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
(30)  
(31) 
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
(32) 
The wetting length on the right side of the plate and that on the left side of the are given respectively by
(33)  
(34) 
When , the meniscus radius is independent of , and . is determined by solving
(35) 
The wetting length on the right side of the plate and that on the left side of the plate are given respectively by
(36)  
(37) 
In either case, the moment on the plate and that on the plate are given respectively by
(38)  
(39) 
Case 3: when the meniscus slides down from the plate and is pinned on the one. In this case and .
(40) 
(41) 
From the condition that contact angle on the plate is , we can get the following relations,
(42) 
and
(43) 
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
(44)  
(45) 
Case 4: when the meniscus slides down from the plate and is pinned on the one. In this case and .
(46) 
(47) 
From the condition that contact angle on the plate is , we can get the following relations,
(48) 
and
(49) 
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
(50)  
(51) 
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
(52)  
(53) 
where is determined by solving from the given volume
(54) 
and is
(55) 
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
(56)  
(57) 
where is determined by solving from the given volume
(58) 
and is
(59) 
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.