## Abstract

When a fluid-immersed array of supported plates or pillars is dried, evaporation leads to the formation of menisci on the tips of the plates or pillars that bring them together to form complex patterns. 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- and three-dimensional setting of the problem. For the case of plates, we explicitly derive the interaction torques based on the relevant physical parameters associated with pillar deformation, contact-line pinning/depinning and fluid volume changes. A Bloch-wave analysis for our periodic mechanical system captures the window of volumes where the two-plate eigenvalue characterizes the onset of the coalescence instability. We then study the evolution of these binary clusters and their eventual elastic arrest using numerical simulations that account for evaporative dynamics coupled to capillary coalescence. This explains both the formation of hierarchical clusters and the sensitive dependence of the final structures on initial perturbations, as seen in our experiments. We then generalize our analysis to treat the problem of pillar collapse in three dimensions, where the fluid domain is completely connected and the interface is a minimal surface with the uniform mean curvature. Our theory and simulations capture the salient features of experimental observations in a range of different situations and may thus be useful in controlling the ensuing patterns.

## 1. Introduction

While the assembly of complex macromolecular structures on microscopic length scales [1] is driven by van der Waals interactions, dispersive forces and chemical interactions, on mesoscopic length scales of the order of micrometres to millimetres in the context of colloids and larger particles, other surface forces such as those due to capillarity play an important role [2–4]. When particles interact with each other via capillary forces, the resulting patterns are a function of the size and shape of the constituents, and the constraints on their movement. Capillary coalescence is a natural consequence of this and occurs when free particles aggregate at an interface [5,6], and also when extended objects such as pillars and plates are brought together by interfacial forces which drive aggregation [7–10]. 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, soft, slender pillars or plates are fully immersed in a liquid which is subsequently evaporated so that capillary forces at the liquid–gas interfaces bring the constituents together. The competition between capillarity and elasticity selects the structures formed, which can be used to construct substrates with tunable wetting and adsorption properties [11,12]. A variety of experimental systems that fall into this category include millimetre-scaled macroscopic brush hairs [8], micrometre-scaled mesoscopic polymeric surface mimicking gecko foot hairs [13], as well as nanometre-scaled carbon nanotube forests [14,15].

Elastocapillary interaction has been well characterized in two-body systems [9,16,17]. Building on this, aspects of collective bundling behaviour of elastic pillars and plates were studied using one-dimensional static energy minimization arguments to estimate the assembled bundle size [8,18,19]. Recent work has delved into the dynamics of elastocapillary coalescence of an array of elastic elements separated by thin fluid films, where the fluid flow is described by a low-dimensional lubrication-type model. In particular, Gat & Gharib [20] have shown that pairwise clustering is the fastest growing mode by using a linear stability analysis of the unpatterned base state; iterating this pairwise collapse results in an ordered hierarchical coalescence structure. Singh *et al.* [21] have also captured the same most unstable pairwise mode in terms of a spring-block model, but their numerical results suggest that the dynamics arising from different perturbations can be considerably more complex than pairwise collapse. Complementing these discrete approaches, a phenomenological continuum mean-field approach has been used to model the arrested coarsening [22], 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 [23].

However, 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 predictive power.

Thus, to predict and control the assembled structures in capillarity-driven self-assembly experiments, we need to (i) include the irreversible effects associated with evaporation, contact line motion and adhesion, that have not been treated so far, (ii) account for geometrically nonlinear effects due to large deviations from the background state that go beyond a linearized analysis, and finally, (iii) generalize our study to the three-dimensional problem of pillar coalescence, that also remains unexplained.

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 plates 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, analyse 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 [7] (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 §2, 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 §3, 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 §4, we conclude with a description of open problems in this rich area.

## 2. Collective dynamics of elastic plates

### (a) Experimental observations

For the two-dimensional case of plate coalescence (figure 2), we consider a one-dimensional array of elastic micro-plates with height *L*, thickness *h* and uniform spacing *D*, Young's modulus *E* and Poisson's ratio *ν*, respectively, wetted by a liquid of surface tension *σ*, density *ρ* and viscosity *μ*. Each plate is assumed to be free at one end and anchored at the other end to a substrate [10]. Between neighbouring plates, defining a cell, we assume that there is a liquid reservoir. When the system is completely immersed in the liquid, the stable configuration is a uniform array of non-interacting vertical plates. 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 plates 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 neighbouring plates tilt towards each other to form a dimer, which then collapses into quadrimers (figure 1*a*). The process repeats until the elastic resistance owing to bending eventually becomes large enough to prevent further coarsening. In a system with large imperfections, irregular bundles can and do arise but the two-plate-collapse mode still appears to be dominant at the onset of instability (figure 1*b*). After the liquid dries out, bundles separate if the adhesion between pillars or plates is not strong to counterbalance the elastic forces; else they persist. In our experiments, the liquid used is isopropyl alcohol (IPA), and the plates and the substrate are made of polydimethylsiloxane (PDMS). Throughout this entire section, we use the following experimental parameters *σ*=0.022 N m^{−1}, *ρ*=786 kg m^{−3}, *μ*=1.96×10^{−3} Pa *s*, *E*=1.5 MPa, *ν*=0.5, *h*=10 μm, *L*=40 μm, *D*=10 μm and gravity *g*=9.8 m s^{−2}.

### (b) Mechanics: coupling plate bending to fluid interface shape

In our experiments, the plates are short and stiff, and remain almost straight as they are deflected by capillary forces, bending primarily in the neighbourhood of the base. Therefore, we approximate each plate as a rigid element with a localized bending response in an elastic hinge at the base [12]. This simplifies our analysis relative to the case of inhomogeneous bending and buckling (electronic supplementary material, appendix A) of individual plates [20,24]. The hinge elastic constant can be approximately derived from the bending response of a short cantilever by a transverse force *F* at its free end, with
*δ*=4*FL*^{3}(1−*ν*^{2})/*Eh*^{3}+*α*_{s}*FL*/*Gh* is the deflection at the free end, *θ* is the tilting angle between the straight plate and the horizontal direction (figure 2*a*), *G* is the shear modulus and *α*_{s} is the shear coefficient. The second term of the right-hand side of equation (2.1) is due to shear deformations in the so-called Timoshenko beam theory [25] for beams with moderate thickness. We have taken *ν*=0.5 and *α*_{s}=10(1+*ν*)/(12+11*ν*) 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 inside the liquid can be non-uniform owing to the effects of gravity and flow. However, in our system, gravity can be neglected, because the Bond number *Bo*=Δ*ρgD*^{2}/*σ*∼10^{−5} is small. Comparing the viscous moment owing to flow *M*_{σ}∼*σL*^{2}/*D*, we find that *M*_{μ}/*M*_{σ}∼10^{−4}≪1 for *a*.), which indicates that pressure due to fluid flow can be neglected. Then, the air–liquid interface is always an arc of a circle because of the uniform pressure in each cell. For each plate, the bending moment results from both the line tension at the contact line, and the pressure jump across the plate caused by the curvature difference between two neighbouring menisci. To calculate the moment we need to find the wetting length, the contact angle, and the pressure difference across the air–liquid interface owing 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 plate while being concave-up (figure 2*b*1), flat (figure 2*b*2) or concave-down (figure 2*b*3). 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*b*4–*b*6). The meniscus profile and the resulting moment per unit depth *M*_{n} can be calculated for any given tilting angles *θ*_{n−1}, *θ*_{n}, *θ*_{n+1} and the liquid volumes per unit depth *V* _{n−1} and *V* _{n}. All the lengths, volumes *V* _{n} and moments *M*_{n} henceforth in this entire section are dimensionless and have been scaled by *L*, *L*^{2} and *σL*, respectively.

Then, the dynamics of the *n*th plate, neglecting inertia (electronic supplementary material, appendix C), follows the overdamped first order equation of motion
*C* is the damping coefficient, *k* is defined in equation (2.1) and the dimensionless moment *M*_{n} is due to capillarity (electronic supplementary material, appendix B). To estimate *C*, 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 *C*≈*τ*_{m}*k* (electronic supplementary material, appendix C), where *τ*_{m} is the time scale for the plate to relax mechanically due to the internal material viscosity of PDMS.

### (c) Onset of bundling: linear stability of the uniform base state

The base state consists of 2*N* vertical plates with *θ*_{n}=*π*/2, *n*=1,2,…,2*N* and a uniform meniscus associated with constant liquid volume *V* in each cell that separates the plates. As this volume is decreased, there is a potential for instability as the curvature of the menisci increase. To determine the stability of this state, we study the perturbations in the bending moment as a function of variations in the angles *θ*_{n} linearized around the current state
*N*×2*N* stiffness matrices owing to the elasticity of the plate and the geometrical change of menisci, respectively. *k*_{b}=(3*σL*^{2}/*Eh*^{3}+6*σ*/7*Gh*)^{−1}, and *θ*_{n}=*π*/2, is given by
*k*_{1}=∂*M*_{n}/∂*θ*_{n+1}=∂*M*_{n+1}/∂*θ*_{n} and *k*_{2}=−∂*M*_{n}/∂*θ*_{n}=−∂*M*_{n+1}/∂*θ*_{n+1} are the effective spring constants associated with the two neighbouring plates owing to capillarity and *M*_{n} are the effective capillary bending moments (see electronic supplementary material, appendix B for a detailed derivation).

For the case when the meniscus is pinned pinned at both tips,
*β* is determined by the volume (electronic supplementary material, appendix B)
*n* for the translationally invariant base state. Note that *k*_{1} and *k*_{2} in this case are not necessarily equal, because while the wetting length is constant, the contact angles change by different amounts on the two tilted neighbouring plates.

For the case when the meniscus is no longer pinned at both tips,
*L*_{w} the wetting length determined by the volume (see electronic supplementary material, appendix B)
*k*_{1} and *k*_{2} are identical, because the contact angle is constant, and the change of wetting length is the same on both plates when they are deflected. From equations (2.5)–(2.8), we see that for any given geometric parameter *d* and critical contact angle *α*, *k*_{i}=*k*_{i}(*V*) (*i*=1,2). The inset of figure 3*a* shows that *k*_{i} switches sign as *V* decreases, corresponding to the case when the system becomes unstable. We note that the discontinuity in *k*_{i} 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 *k*_{1}≥0, the meniscus is pinned at both tips, and the smallest eigenvalue is *k*_{1}<0, the smallest eigenvalue is λ(*π*)=2*k*_{1}+2*k*_{2}+*k*_{b}, so that stability is controlled by the competition between elasticity and capillarity. As *k*_{i}=*k*_{i}(*V*) (*i*=1,2), the condition λ(*π*)<0 sets the range of *V* in which the system is unstable (figure 3*a*). The primary eigenmode corresponds to *f*=*π*, and the eigenvector thus is [1,−1,1,−1,…,1,−1] from equation (2.9), corresponding to dimerization of the plates.

This result is consistent with the predictions of spring-block models for the same process [20,21] that show that a pair of blocks should be unstable whenever the dimensionless spring stiffness is below a critical value, which in our case corresponds to the critical volume of the liquid. However, by correctly accounting for the geometry of the deforming pillars and the changes in fluid volume induced by their deflection as well a minimal model for contact line pinning and depinning, we also see there is both a lower and an upper limit on the volume of liquid for which instability can be observed, resulting in a finite window of instability. This is because at relatively large volumes, the meniscus is pinned at the plate or pillar tips and can have either a concave or convex shape with a relatively small curvature so that the capillary stiffness is positive and the capillary pressure jump across the plates is able to stabilize any perturbations. This calculation also shows that to avoid the plate collapse, it suffices to have *k*_{1}>0 or increase *k*_{b} when *k*_{1}<0. Practical approaches to implement this include the use of liquid with low surface tension and contact angle close to 90°, so that *k*_{1} 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 *h*/*L*.

To explain the nature of the instability, we note that the fastest growing mode is the dimer mode and assume *θ*_{n−1}−*π*/2=*π*/2−*θ*_{n}=*θ*_{n+1}−*π*/2=*ϕ*. In the vicinity of the critical volume *V* * at which the instability happens, the dynamics of a plate is determined by expanding equation (2.2) in powers of *ϕ* and yields the equation
*c* is the dimensionless damping coefficient, and *g*(*V* *)>0 is an algebraic coefficient of the cubic term (electronic supplementary material, appendix C). In (2.11), even orders of *ϕ* do not appear due to the reflection symmetry *ϕ*→−*ϕ* inherent in the system. When the volume in a cell reaches the critical value *V* =*V* * (equivalently *β*=*β**), 2*k*_{1}+2*k*_{2}+*k*_{b}=0, while when 0<*V* *−*V* ≪1, a Taylor series expansion of equation (2.5) in the neighbourhood of *β** gives us 2*k*_{1}+2*k*_{2}+*k*_{b}∼*β**−*β*, while a similar expansion of equation (2.6) in the neighbourhood of *β** yields *V* *−*V* ∼*β*−*β**, and hence 2*k*_{1}+2*k*_{2}+*k*_{b}∼*V* −*V* *. From equation (2.11), we see that the stable equilibrium state has two branches of solutions *V* *>*V* and only one solution *ϕ*=0 otherwise, suggesting that the bifurcation is supercritical. Figure 3*b* shows both the positive branch of the asymptotic solution *ϕ*_{1}, and the results of simulation for two plates obtained by solving equation (2.2) with periodic boundary conditions. The results agree well, and confirm that the supercritical instability leads to plate dimerization, as shown in the inset.

### (d) Nonlinear dynamics: drying, coarsening and refining

#### (i) Controlled liquid volume

For an initially translationally invariant system although the two-plate-collapse mode is the fastest growing mode. However, the cluster of dimers is not necessarily the final stable state (dots in figure 3*b*). We note 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 §2c breaks down. Therefore, we numerically integrate equation (2.2) directly to follow the hierarchical dynamics by which an array of vertical individual plates first forms dimers, then quadrimers, and so on until the very large bundles are arrested due to elastic effects. Figure 4*a* shows snapshots of dynamical coarsening for different values of the control parameter *V* , in a periodic domain of 32 plates, triggered by tilting one plate by 0.1% from 90° as the initial condition. As expected, for volumes outside the range *V*/*d*∈(0.5073,0.9393), e.g. for the cases of *V*/*d*=0.95 and 0.4, 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 *V* . As examples, we see that for *V*/*d*=0.92, quadrimers sweep through the domain right after the dimers form and the system reaches the stable equilibrium, for *V*/*d*=0.80 the dimers persist for a while before they eventually collapse to quadrimers, and for *V*/*d*=0.70 the dimers stay stable, as shown in figure 4*c*. The propagating waves of coalescence that we see as these patterns evolve are similar to those in earlier models [21] wherein a localized perturbation generates a propagating front that leaves behind clusters of a finite size. In the context of our model, when *V* <*V* *, the springs have a negative stiffness. Then, a localized perturbation results in a single plate element escaping from its equilibrated vertical configuration. Rotation of one plate causes the effective stiffness of its adjacent neighbours to change, with one increasing and another decreasing due to the changing meniscus levels and contact angles even as the liquid volume of a cell is conserved. The local change in stiffness then breaks the symmetry between two adjacent plates and either pushes them apart or brings them together to form pairs. Repeating this pair-wise clumping at the front where the deflected and undeflected elements meet results in a propagating wave of disturbance from the localized perturbation in both directions. The clumped pairs behave as effectively new elements. Iterating this process until the elastic resistance prevents further coalescence selects a well-defined bundle size.

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

#### (ii) Coupling capillary coalescence to drying dynamics

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 plates. 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 [8,18,19] or renormalization analysis [20].

A minimal evaporation model that is sufficient to capture the qualitative features of hierarchical bundle formation is given by
*c*_{e} is a constant, *τ*_{m} and *τ*_{e} are time scales for mechanical relaxation and evaporation, respectively, *r*_{n} is the radius of the meniscus and 2*β*_{n} is the angle subtended by the meniscus arc. Equations (2.2) and (2.12) coupled together determine the dynamics of plate coalescence driven by evaporating liquid, with three regimes. When *τ*_{m}≪*τ*_{e}, *V* _{n} decreases quasi-statically so that equation (2.11) relaxes to a static state for a prescribed value of *V* _{n}. When *τ*_{m}∼*τ*_{e}, the evolution of the cell volumes *V* _{n} and the plate configuration *θ*_{n} are coupled. When *τ*_{m}≫*τ*_{e}, the evaporation is so fast that the plate array does not coarsen. In our typical experiments, *τ*_{m}≈16 ms and *τ*_{e} is of the order of seconds, and we choose *c*_{e} in equation (2.12) so that the liquid dries out in approximately 5 s. Figure 5 shows the simulation results for a periodic domain of 100 plates with two 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 *V* _{99} is smaller by 2% 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. If adhesion between contacting plates is neglected, capillary forces are not sufficiently large to hold the plates together, and the bundles separate symmetrically when the liquid volume falls below 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*, we show the results associated with an initial uniform random perturbation with maximum relative amplitude of 5% that is applied to all the tilting angles and liquid volume in all cells. Irregular bundles of size ranging from two to five plates form transiently, and separate as the liquid evaporates. Figure 5*d* presents the number count of bundles of different size for figure 5*c*; although the dimer is still the dominant mode in the early stages, bundle formation can be irregular. This is because 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 et al. [7] and Gat & Gharib [20]. All these features also agree qualitatively well with our own experimental observations shown in figure 1*a*,*b*.

## 3. Collective dynamics of a two-dimensional array of pillars

### (a) Experimental observations

We now generalize our study of the one-dimensional dynamics of plate 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 [7,26]. The dynamics of pillars is different from that of plates in three qualitative 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 neighbours. Second, the three-dimensional geometry allows the pillars to bend in two principal directions, and also twist after the pillars are in contact to further reduce the surface energy. For pillars with a circular cross-section, the twist must be a constant. If we neglect friction between adhering pillars, the twist must identically vanish, and here we will assume that this is the case. Third, experimentally we see that the segmented, ‘wormlike’ geometry of specially treated pillars increases pinning of the receding contact line by reentrant curvature [7]; here will neglect the reentrant curvature and consider the menisci to be pinned at the pillar tips for simplicity.

As in the case of plates, 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 [26]. 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 *L*=4.5 μm, the pillar radius *R*=150 nm, the pillar spacing *D*=2 μm, the Young's modulus *E*=0.2 GPa, the surface tension of the liquid *σ*=0.022 N m^{−1} and the density of the liquid *ρ*=786 kg m^{−3}. Unlike in the one-dimensional array of plates, where the dimer is the primary unstable mode, for pillars the quadrimer is the primary unstable mode. As the liquid evaporates, this mode leads to a hierarchical growth of larger assemblies which eventually get arrested elastically.

### (b) Mechanics: coupling pillar deformation to fluid interface shape

As in the plate 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 ** X**(

*x*,

*y*) (figure 6

*a*), and is given by

*c*is the drag coefficient,

*F*_{b}is the elastic bending resistance force at the tip,

*F*_{σ}is the capillary driving force owing to surface tension

*σ*and

*V*is the liquid volume in the system.

To estimate the drag coefficient *c* and the time scale in equation (3.1), we consider both the internal contribution from viscoelasticity of the solid material and the external contribution from the viscous fluid and find that viscoelastic effects dominate. Indeed, AFM characterization of the epoxy micropillars [27] show that the force versus displacement curve has a clear hysteresis loop, which reveals that the epoxy micropillars are viscoelastic. For a pillar of radius *R*=1.5 μm, length *L*=9 μm and bending spring constant *k*=3*πER*^{4}/4*L*^{3}≈3.3 N m^{−1}, where *E*=0.2 GPa is the Young's modulus, the measured hysteresis accounts for 32% of the total work done at max deflection *x*=1.5 μm and rate of deflection *v*=6 μm *s*^{−1}. Thus, the effective damping coefficient *c*_{1} from internal viscoelasticity can be calculated from 2*c*_{1}*vx*≈0.32×0.5×*kx*^{2}, based on the Kelvin–Voigt model. To calculate the damping coefficient *c*_{2} from the external viscous fluid, we consider a circular cylinder of radius *R* moving with velocity *U* normal to its axis at small Reynolds number *Re*=2*RρU*/*μ*≈7×10^{−6}, in which case, the drag of magnitude *μ*≈10^{−3} Pa s is the viscosity of the fluid. As the centre of the pillar is moving at a velocity of almost *v*/2 when the cantilever tip is deformed at a velocity of *v* in the AFM test, we take *c*_{1}/*c*_{2}∼10^{7} for the experimental parameters listed above, so that external damping from viscous fluid can be neglected and *c*≈*c*_{1}. Therefore, the time scale for the pillar to relax mechanically is *τ*_{m}=*c*/*k*∼10^{−2} *s* in equation (3.1), and *τ*_{m} is a material property and does not depend on the geometric dimensions of the pillar.

To compute the bending resistance force *F*_{b} in the horizontal direction, we use the theory of the elastica for the inextensional, unshearable deformation of thin pillar. Letting *ϑ* be the angle of the tangent to the pillar centreline with the vertical, with *s*∈[0,*L*] the arc length coordinate, and |*F*_{b}| the force amplitude at the tip, equilibrium implies that
** X**, |

*F*_{b}| is uniquely determined, and

*F*_{b}=|

*F*_{b}|

**/|**

*X***|. To obtain**

*X*

*F*_{σ}, we need to determine the shape of the air–liquid interface. Since the Bond number

*Bo*∼10

^{−6}is small, gravity can be neglected. Moreover, the time scale for the fluid to equilibrate in the porous brush

*t*

_{f}is much smaller than that for the pillars to relax mechanically

*t*

_{m}, which is much smaller than that for the evaporation

*t*

_{e}, i.e.

*t*

_{f}∼10

^{−3}

*s*≪

*t*

_{m}∼10

^{−2}

*s*≪

*t*

_{e}∼10

^{0}

*s*(electronic supplementary material, appendix D). Assuming the ambient pressure to be zero, the pressure in the liquid domain can be regarded as uniform and the air–liquid interface

*z*=

*S*(

*x*,

*y*) will thus be a surface of uniform mean curvature that satisfies the equation

*H*is the mean curvature of the interfacial surface. Volume conservation in the whole domain yields

*p*. Here

*A*is the projected domain of the air–liquid interface to the horizontal plane (meshed area in figure 6

*a*), and

*h*

_{i}is the elevation of

*i*th pillar tip. Since the menisci are always pinned on the pillar tips, we need to solve equation (3.3) 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*)

*F*_{σ}, 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

*c*represents the tip circle and

**is its unit outward normal. Note that for very small values of**

*n**V*, 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 a short-range repulsion force when the discs representing pillar tips come close enough (10% of the pillar diameter), but do not consider the elastic deformation of the cross section owing to contact.

### (c) 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 finite element code. For a given pillar tip displacement *X*_{i}, equation (3.2) is solved with the geometric condition *F*_{bi} and tip elevation of each pillar *h*_{i}. Given *X*_{i} and *h*_{i}, equation (3.3) is solved using the finite element method to determine the interface *S*(*x*,*y*) on the multiply-connected domain with the liquid volume constraint equation (3.4), which determines *p* and thence *F*_{σi}. Then equation (3.1) is integrated in time explicitly to update the pillar tip positions, and the domain is re-meshed accordingly at every time step.

Figure 6*d*–*g* show the simulated evolution of an array of 14 by 14 pillars collapsing into fourfold bundles for a prescribed liquid volume of *V*/*V* _{flat}=0.90, where *V* _{flat} 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 centre 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 plate collapse. However, here the primary eigenmode of instability has a fourfold symmetry (figure 6*e*) independent of initial perturbations as observed in experiments (figures 6*g* and 7). An intuitive way to understand this is to recognize that the fourfold bundles have two principal directions, along each of 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 the effective spring constants in the two principal directions are decoupled, which is true in the linearized regime, the dimer is the primary mode in each direction as in the one-dimensional plate 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 favourable 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*).

## 4. Conclusion

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 plate collapse, we explicitly derived the conditions for the primary dimerization instability in terms of the state variables—tilting angles *θ*_{n} and liquid volumes *V* _{n}, and the relevant geometrical and physical parameters. Our analysis correctly accounts for the geometry of the system and consequently predicts a finite window of liquid volumes over which the vertical configuration is unstable, consistent with observations. Complementing our linear analysis, we also present full numerical simulations that show how 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 depends 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 plates, 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 plate coalescence [23], but the question for three-dimensional coalescence remains an open question.

## 5. Funding statement

We thank the Harvard-MRSEC DMR-0820484, the MacArthur Foundation (LM), and the NRF of Korea (Grant No. 2013034978, H.-Y.K.) for support.

## 6. Author contributions

Z.W. contributed to the formulation, did the numerical simulations and drafted the manuscript; T.S. contributed to the formulation and drafted the manuscript; J.K., H.Y.K. and J.A. contributed to the experiments; L.M. conceived of the study, designed the study, contributed to the formulation of the problem, and drafted the manuscript. All authors gave final approval for publication.

## 7. Conflict of interests

The authors have no competing interests.

## Acknowledgements

We thank Sam Ocko for many discussions that helped to sharpen and clarify our arguments.

- Received August 4, 2014.
- Accepted December 22, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.