Buckling in Armored Droplets

The issue of the buckling mechanism in droplets stabilized by solid particles (armored droplets) is tackled at a mesoscopic level using dissipative particle dynamics simulations. We consider spherical water droplet in a decane solvent coated with nanoparticle monolayers of two different types: Janus and homogeneous. The chosen particles yield comparable initial three-phase contact angles, chosen to maximize the adsorption energy at the interface. We study the interplay between the evolution of droplet shape, layering of the particles, and their distribution at the interface when the volume of the droplets is reduced. We show that Janus particles affect strongly the shape of the droplet with the formation of a crater-like depression. This evolution is actively controlled by a close-packed particle monolayer at the curved interface. On the contrary, homogeneous particles follow passively the volume reduction of the droplet, whose shape does not deviate too much from spherical, even when a nanoparticle monolayer/bilayer transition is detected at the interface. We discuss how these buckled armored droplets might be of relevance in various applications including potential drug delivery systems and biomimetic design of functional surfaces.

Pickering emulsions [1], i.e. particle-stabilized emulsions, have been studied intensively in recent years owning to their wide range of applications including biofuel processing [2] and food preservation [3,4]. They have also been developed as precursors to magnetic particles for imaging [5] and drug delivery systems [6]. Even with their widespread use, they remain, however, underutilized. In Pickering emulsions, particles and/or nanoparticles (NPs) with suitable surface chemistries adsorb at the droplet surfaces, with an adsorption energy of up to thousands of times the thermal energy. The characteristics of Pickering emulsions pose a number of intriguing fundamental physical questions including a thorough understanding of the perennial lack of detail about how particles arrange at the liquid/liquid interface. Other not completely answered questions include particle effects on interfacial tension [7], layering [8], buckling [9][10][11] and particle release [8,12].
In some important processes that involve emulsions, it can be required to reduce the volume of the dispersed droplets [9,[13][14][15]. The interface may undergo large deformations that produce compressive stresses, causing localized mechanical instabilities. The proliferation of these localized instabilities may then result in a variety of collapse mechanisms [8,10,11]. Despite the vast interest in particle-laden interfaces, the key factors that determine the collapse of curved particle-laden interfaces are still subject of debate. Indeed, although linear elasticity describes successfully the morphology of buckled particle-laden droplets, it is still unclear whether the onset of buckling can be explained in terms of classic elastic buckling criteria [16,17], capillary pressure-driven phase transition [9], or interfacial compression phase transition [18]. Numerous experiments have been conducted to link the rheological response of particle-laden interfaces to the stability of emulsions and foams. However, their results could be dependent on the method chosen for preparing the interfacial layer. Due to their inherent limited resolution, direct access to local observables, such as the particles three-phase contact angle distribution, remains out of reach [19]. This crucial information can be accessed by numerical simulations sometimes with approximations. All-atom molecular dynamics (MD) simulations have become a widely employed computational technique. However, all-atom MD simulations are computationally expensive. Moreover, most phenomena of interest here take place on time scales that are orders of magnitude longer than those accessible via all-atom MD. Mesoscopic simulations, in which the structural unit is a coarse-grained representation of a large number of molecules, allow us to overcome these limitations. It is now well established that coarse-grained approaches offer the possibility of answering fundamental questions responsible for the collective behaviour of particles anchored at an interface [20].
We employ here Dissipative Particle Dynamics (DPD) [21] as a mesoscopic simulation method. We study the shape and buckling transitions of model water droplets coated with spherical nanoparticles and immersed in an organic solvent. The procedure and the parametrisation details are fully described in prior work [22][23][24] and in the Supporting Information (SI). The particles are of two different types: Janus and homogeneous. They are chosen so that the initial three-phase contact angles (≈ 90 • ) result in maximum adsorption energy. The volume of the droplets is controllably reduced, pumping randomly a constant proportion of water molecules out of the droplet (more details in the SI). At every stage we remove 10 percent of the water from the droplet. Throughout this letter, E i refers to the i th removal of water, with E 0 corresponding to the initial configuration and E 20 to the final configuration. We seek FIG. 1. Sequence of simulation snapshots representing buckling processes of water in oil droplets armored with 160 spherical Janus (top) and homogeneous (bottom) nanoparticles after successive removals of water. The number of water beads removed increases from left to right with Ei refering to the i th removal. Cyan and purple spheres represent polar and apolar beads, respectively. Pink spheres represent water beads. The oil molecules surrounding the system are not shown for clarity.
to determine whether the NPs at the droplet interface buckle, causing the droplets to deviate from the spherical shape. We show that Janus particles affect strongly the shape of the droplet via the formation of a crater-like depression. This evolution is actively controlled by a closepacked particle monolayer at the curved interface. On the other hand, homogeneous particles follow passively the volume reduction of the droplet. The shape of the droplet remains approximately spherical with a nanoparticle monolayer/bilayer transition, with some NPs desorbing in water. We discriminate the two mechanisms with the evolution of their respective nanoparticle threephase contact angle distributions. While for Janus particles the distribution remains unimodal, albeit skewed when the droplet is significantly shrinked, for homogeneous particles, the evolution of the contact angle distribution becomes bimodal with some particles becoming more/less immersed in the aqueous phase. We consider a system initially made by a spherical water droplet immersed in oil, and stabilized by a sufficiently dense layer of NPs [24]. The initial shape of the droplet is spherical. The only difference between the two systems is the NP chemistry, i.e. the distribution and proportion of polar and apolar beads around the spherical particles and their efficiency in interacting with the two fluids at the interface. Janus and homogeneous NPs are designed to present comparable three-phase contact angles, θ c = (91.6 ± 2.0) • and θ c = (88.7 ± 3.5) • , respectively (cf. SI for details). We consider throughout this study the same NP density on the droplets. We calculate the radius of gyration, R GYR , and the asphericity, A s , for the droplet covered by either Janus or homogeneous NPs (cf. SI for details). For the initial configurations, we ob-tain R GYR = 13.837 ± 0.003 and R GYR = 13.860 ± 0.003, and A s = 0.156±0.05 and A s = 0.153±0.05, respectively, expressed in R C units (cf. SI for details).
In Fig.1 we show representative snapshots obtained during the simulations for systems containing Janus NPs (top panels) and homogeneous NPs (bottom panels). Visual inspection of the simulation snapshots highlights some fundamental differences between the two buckling processes. We start with spherical initial droplets (E 0 ). When the water droplet is coated with Janus particles (top), the system starts developing dimples as moderate amount of water is removed (E 2 ). The morphology then becomes more crumpled with increasing numbers of dimples (E 5 ). For stronger removal, the droplet geometry evolves to a large and smooth curved shape, yielding a crater-like depression to minimize the interfacial energy of the system (E 8 and E 20 ). During this evolution, Janus NPs remain strongly adsorbed at the interface, forming a close-packed monolayer between the two fluids.
The buckling process is fundamentally different when the water droplet is stabilized with homogeneous NPs (bottom). When the volume of the droplet is reduced, the shape of the system evolves smoothly and does not present any sharp transition to morphologies showing dimples and cups, nor crater-like depressions. Instead, the NPs reorganize progressively into a bilayer, presumably to minimize the system energy. Unlike Janus NPs, homogeneous NPs either protrude exceedingly towards the decane solvent, or recede into the water droplet with some particles even desorbing into the water phase (from E 2 to E 20 ). For reference, we recall that the change in energy accompanying desorption of a spherical particle from the oil-water interface to either bulk phase is ap- proximated by ∆E = πr 2 γ ow (1 ± cos θ) 2 , in which r is the particle radius, γ ow is the bare oil-water interfacial tension, and the plus (minus) sign refers to desorption into oil (water) [19]. Even if this expression assumes the oil-water interface remains planar up to the contact line with the particle, it can give a rough approximation of the energy at play. Considering the system parameters given in the SI, we obtain ∆E ≈ 85 k B T in our systems when one NP desorbs.
These two different behaviours are quantitatively investigated in Fig. 2, where we show the temporal evolution of the radius of gyration (left panel), R GYR , and the asphericity (right panel), A s , of the two droplets as a function of the dimensionless parameter ∆N W ≡ N W /N 0 W , with N W the number of water beads remaining in the droplet, and N 0 W the initial number of water beads in the droplet. When N W > 0.6, the radius of gyration of two systems follow the same evolution, regardless the chemistry of the NPs (Janus or homogeneous). For one droplet coated with Janus NPs, R GY R then departs from its linear trend when N W < 0.6. This departure corresponds to the evolution from E 5 to E 8 in Fig. 1, i.e. the transition from a droplet interface made of dimples and cups to the formation of the crater-like depression. During this transition, the size of one dimple increases when the system relaxes after evaporation. This local evolution yields a larger depression, which causes the progressive coalescence of the small dimples. This transition is consistent with the surface model numerical analysis from Ref. [17], which study the shape evolution of a spherical elastic surface when the volume it encloses is decreased. This model, which has long been considered as valid to describe the deformation of thin shells [16,25], showed that a thin shell with a single dimple has lower energy than a shell containing multiple dimples. This occurs because elastic energy mainly concentrates in dimple edges as bending energy. Dimples coalescence lowers the total elastic energy. Below ∆N W ≈ 0.6, R GY R increases as the droplet can be described as half-sphered. Let us note that this evolution is coherent with the temporal evolution of the radial distribution function of the NPs, g(r), with r the distance between the centers of the NPs, given in the SI. In contrast, the droplet coated with homogeneous NPs shrinks isotropically when the volume reduces even below ∆N W ≈ 0.6. This evolution yields continuous decrease of R GY R and a relatively low A s in Fig. 2. Eventually, the NP concentration becomes too high and some NPs move into the droplet. When N W < 0.25, the number of water beads that remain in the droplet is not sufficient to define unambiguously the droplet volume. This limitation impacts the system shape and the evolution of R GY R and A s for Janus and homogeneous NPs.
We also quantify the NP layer properties to see if the NPs actively influence or passively follow the evolution of the droplet geometry. In Fig. 3, we compare the three phase contact angle distribution of Janus (left panel) and homogeneous (right panel) NPs at the initial stage E 0 , where the shape of the droplet is spherical, and the final stage, E 20 . The initial distributions, fitted with continuous lines, can be described with Gaussian distributions for both NPs. The values of the respective means, µ J and µ H , and variances, σ J and σ H , differ due to the NPs chemistry. We obtain µ J = 91.6 • and µ H = 88.6 • , and σ J = 2.0 • and σ J = 3.4 • for Janus and homogeneous NPs, respectively. When the droplet coated with Janus NPs shrinks, the contact angle distribution evolves to a skewed one, but it remains unimodal, with a single peak centered at the same value as the one measured for the initial configuration. The emergence of the skewness of the distribution is linked to the decrease of the NP-NP distance when the droplet volume is reduced. It is due to the major role played by steric effects. As discussed earlier, to minimize its interfacial energy, the system must deform its shape, eventually forming a crater-like depression. We conclude this transition is achieved through the active role played by the Janus NPs. In the final structure, some NPs are forced to deviate from their original contact angle, increasing the skewness of the distribution on both sides of the peak.
The evolution of the system is different when homogeneous NPs are present. As the droplet volume is reduced, the contact angle distribution firstly evolves as a monolayer interface with a single peak (cf. SI). As the droplet shrinks further, and the distance between the NPs decreases, the distribution becomes bimodal, with two distinct peaks emerging on both sides of the initial equilibrium contact angle. This feature is characteristic of a particle bilayer. Indeed, homogeneous NPs are more weakly attached to the interface than Janus NPs. In the case of buckling mechanism studied here, the homogenous NPs mainly follow the volume reduction, sharing the interfacial area, either receding into the water droplet or protruding towards the organic solvent. Unlike Janus NPs, homogeneous NPs do not drive the evolution of the droplet shape, which does not differ too much from the spherical geometry. The behaviour just described is characteristic of the passive role played by the homogeneous NPs, which mainly follow the volume reduction, only modulating the droplet shape due to the steric constraints.
The curved shape obtained when the droplet is coated with Janus NPs can also be characterized by the wettability associated with the local arrangement of the NPs at the interface. The particles with a contact angle θ c < 85 • , i.e. the blue ones in Fig. 3 (left panel), can be found in the crater-like depression. These particles have receded into the water droplet due to the concave local geometry of the interface. The particles with θ c > 100 • , i.e. the red ones in Fig. 3 (left panel), can be found at the transition between the concave and convex areas of the interface, where they are likely to protrude towards the solvent. The shape deformation of the droplet is achieved through the active role played by the Janus NPs. Their specific chemistry causes them to create an interface with excess wettability (θ c < 85 • ) in a pocket delimited by the crater-like depression, and surrounded by a cup with low wettability (θ c > 100 • ).
Our results are consistent with experiments reporting buckling and crumpling of nanoparticle-coated droplets [9][10][11]. In particular, we observe a close analogy to the experimental work of Datta et al. [11], who studied water-in-oil droplets of varying sizes. In these experiments the dispersed phase is slightly soluble in the continuous phase. The volume reduction was controlled with the addition of a fixed amount of unsatured continous phase. As shown in Fig. 4, Datta et al. observed droplet shapes including dimples, cups, and folded configurations, in agreement with our simulations (cf. experimental details in the caption in Fig. 4). Unlike our mesoscopic analysis, Datta et al. [11] do not have access to the particle three-phase contact angle distribution. This information provides a deeper understanding of the organisation of the NPs at the interface, and allows us to decipher the active or passive role of the NPs.
As explained in the SI, the layering properties of the particles depend strongly on the numerical protocol. For example, decreasing the relaxation time between successive water removals can induce NP release from the interface, which is in agreement with experiments [12]. The results presented here seem to be due to the chemistry of the nanoparticles simulated (i.e. Janus vs. homogeneous). It is however possible that homogeneous NPs with large adsorption energy become active and yield buckled armored droplets similar to those observed when Janus NPs are simulated here. The new physical insights discussed in this letter could be useful for a variety of applications. For example, controlling the positions of the solid particles with respect to the interface could help in heterogeneous catalysis [26]. In biomimetic design, where the identification and evaluation of surface binding-pockets is crucial, the ability of controlling pockets such as those created by the craterlike depression in the presence of Janus NPs, could play a central role in designing structures with a defined geometry [27]. The analogy between Fig. 1 and the shape of protein active site might play an important role for ligand docking [28,29]. Finally, buckled armored droplets might also be of relevance as potential drug delivery systems [6]. Over the last decade, nanoscale droplets have been used for instant real-time ultrasound imaging of specific organs [5]. Superparamagnetic solid NPs provide a means of manipulating the droplets using an external magnetic field [5]. One of the main limitations in such applications is droplet coalescence, which can happen before droplets reach the target. The specific shapes obtained with buckled armored droplets might prevent coalescence. Indeed, the NP arrangements on the droplets show increased packing, which reduces significantly the NPs mobility. The particle layers would then provide enough mechanical resistance to guarantee the droplet stability.

Buckling in Armored Droplets
Supporting Information

MD SIMULATION METHOD
The Dissipative Particle Dynamics (DPD) simulation method [21] was implemented within the simulation package LAMMPS [30]. The procedure and the parametrisation details are fully described in prior work [22,23]. The system simulated here is composed of water, oil (decane), and nanoparticles (NPs). One "water bead" (w) represents 5 water molecules and a reduced density of one DPD bead is set to ρ = 3. One decane molecule is modeled as two "oil beads" (o) connected by one harmonic spring of length 0.72 R c and spring constant 350 k B T /R c [31], where R c is the DPD cutoff distance. The initial size of the simulation box is L x × L y × L z ≡ 72 × 72 × 78 R 3 c , where L i is the box length along the i th direction. Periodic boundary conditions are applied in all three directions. The NPs are modelled as hollow rigid spheres and contain polar (p) and nonpolar (ap) DPD beads on their surface. One DPD bead was placed at the NP center for convenience, as described elsewhere [22,23]. Hollow models have been used in the literature to simulate NPs, and hollow NPs can also be synthesized experimentally [32]. We considered spherical NP of the same volume, 4/3πa 3 0 , where a 0 is the radius of the sphere. We imposed a 0 = 2R c ≈ 1.5 nm. All types of beads in our simulations have reduced mass of 1. We maintain the surface bead density on the NPs sufficiently high to prevent other DPD beads (either decane or water) from penetrating the NPs (which would be unphysical), as it has already been explained elsewhere [23]. To differenciate every NPs, we report the nonpolar fraction of the NP surface beads and the NP type. For example, 75HP (55JP ) indicates that 75% (55%) of the beads on the NP surface are nonpolar, and that we consider an homogeneous (Janus) NP.
The interaction parameters shown in Table I are used here. These parameters were adjusted to reproduce selected atomistic simulation results, as explained in prior work [22]. By tuning the interaction parameters between polar or nonpolar NP beads and the water and decane beads present in our system, it is possible to quantify the effect of surface chemistry on the structure and dynamics of NPs at water-oil interfaces. Specifically, the interaction parameters between NP polar and nonpolar beads were adjusted to ensure that NPs are able to assemble and disassemble without yielding permanent dimers at the water/oil interface [22]. All simulations were carried out in the NVE ensemble [22]. The scaled temperature was 1, equivalent to 298.73 K. The DPD time scale can be gauged by matching the selfdiffusion of water. As demonstrated by Groot and Rabone [31], the time constant of the simulation can be calculated as τ = NmDsimR 2 c Dwater , where τ is the DPD time constant, D sim is the simulated water self-diffusion coefficient, and D water is the experimental water self-diffusion coefficient. When a w−w = 131.5 k B T /R c (cf. Tab. I), we obtained D sim = 0.0063 R 2 c /τ . For D water = 2.43 × 10 −3 cm 2 /s [33], we finally obtain τ = 7.6 ps.

SYSTEM CHARACTERISATION: DROPLETS AND NANOPARTICLES
In our simulations, the initial size of the droplet was fixed. At the beginning of each simulation, the solvent (oil) beads were distributed within the simulation box forming a cubic lattice. One water droplet of radius ≈ 20 R c was generated by replacing the oil beads with water beads within the volume of the spherical surface. A number of spherical NPs were placed randomly at the water-decane interface with their polar (nonpolar) part in the water (oil) phase to reach the desired water-decane interfacial area per NP. Following previous work [24,34], the NPs considered in this study are spherical and of two different types: Janus and homogeneous. The emulsion systems are stabilized by a sufficiently dense layer of NPs [24]. We considered water droplets coated with 160 spherical Janus and homogeneous nanoparticles of type 55JP and 75HP , respectively. Considering the NP surface coverage, φ, defined in Ref. [22], we obtain φ ≈ 0.9. Considering the results obtained in Ref. [22] for a flat interface, this yields an interfacial tension γ ow ≈ 6.8 k B T /R 2 c for both Janus and homogeneous NPs. The initial configuration obtained was simulated for 10 6 timesteps in order to relax the density of the system and the contact angle of the nanoparticles on the droplet. The system pressure and the three-phase contact angles did not change notably after 5000 simulation steps. We let then run the system for an additional 2 × 10 6 timesteps to generate two new intial configurations after 2 × 10 6 and 3 × 10 6 timesteps, respectively. We then repeated the buckling simulation with these different initial configurations to test the reproducibility of the simulation.
The surface area of the droplets is slowly diminished, pumping randomly a constant proportion, i.e. 10 percent, of water molecules out of the droplet and letting the system pressure and the three-phase contact angles equilibrate at constant density. By slowly, we mean we do not create any hollow volume in the droplet that would strongly drive the system out-of-equilibrium. Doing so, the three-phase contact angle distribution of the NPs evolves sufficiently smoothly when the droplet buckles and becomes nonspherical, thereby preventing particles to be artifactually realeased. This numerical protocol can be similarly compared with an emulsion system where the dispersed phase is slightly soluble in the continuous phase [11]. By adding a fixed amount of unsatured continuous phase, the volume of the droplets can then be controllably reduced.
Three-phase contact angle. To estimate the three phase contact angle on the droplets we calculate the fraction of the spherical NP surface area that is wetted by water [35], where A w is the area of the NP surface that is wetted by water and R is the radius of the NP. The ratio A w /4πR 2 is obtained by dividing the number of NP surface beads (ap or p), which are wetted by water, by the total number of beads on the NP surface (192 for spherical NP). One surface bead is wet by water if a water bead is the solvent bead nearest to it. One standard deviation from the average is used to estimate the statistical uncertainty.
Radius of gyration. The description of the geometrical properties of complex systems by generalized parameters such as the radius of gyration or principal components of the gyration tensor has a long history in macromolecular chemistry and biophysics [24,36,37]. Indeed, such descriptors allow an evaluation of the overall shape of a system and reveal its symmetry. Considering, e.g., the following definition for the gyration tensor, where the summation is performed over N atoms and the coordinates x, y, and z are related to the geometrical center of the atoms, one can define a reference frame where T GY R can be diagonalized: In this format we obey the convention of indexing the eigenvalues according to their magnitude. We thus define the radius of gyration R 2 GY R ≡ S 2 1 + S 2 2 + S 2 3 , and the asphericity A s ≡ S 1 − 1 2 (S 2 + S 3 ), which measures the deviation from the spherical symmetry. To determine the properties of a droplet, we calculate R GY R and A s using the centers of the water beads.
In this letter, Janus and homogeneous NPs present similar three-phase contact angle θ c = (91.6 ± 2.0) • and θ c = (88.7±3.5) • , respectively. The radius of gyration, R GYR , and the asphericity, A s , for the Janus and homogeneous initial configurations are R GYR = 13.837 ± 0.003 and R GYR = 13.860 ± 0.003, and A s = 0.156 ± 0.05 and A s = 0.153 ± 0.05, respectively, expressed in R C units.

EVOLUTION OF THE NANOPARTICLE RADIAL DISTRIBUTION FUNCTION
The transition from dimples and cups to crater-like depression observed when Janus nanoparticles cover the droplet can be reflected in the temporal evolution of the radial distribution function of the NPs, g(r), with r the distance between the centers of the NPs. We first consider the initial spherical configuration, and extract the list of nearest neighbours of each NP within a shell of radius r < 12R c . This treshold value defines the first and second neighbouring shells [24,38]. As the system is densely packed at the interface, we follow the temporal evolution of g(r) considering this subset of particles, which corresponds to the nearest neighbours shell. When the volume of the droplet reduces, we first see in Fig. S1 (left panel) the emergence of a new peak around r ≈ 4.75R C < r DPD , where r DPD = 5R C represents the distance between the centers of NPs above which the NPs do not interact with each other through the DPD non-bonded force. As the droplet volume further reduces, the heigh of the new peak increases. This evolution continues until ∆N W ≈ 0.6 (cf. Fig. 2 in the main text). Below that value, the heigh of the first peak decreases, together with the increase of the heigh of the second peak. This corresponds to the transition from dimples and cups to crater-like depression. The evolution continues until ∆N W ≈ 0.25 where the number of water beads which remain in the droplet is not sufficient to define unambiguously the droplet volume.
As the g(r) evolution between the homogeneous particles is concerned, we see in Fig. S1 (right panel), a constant increase of the heigh of the first peak, centered at r ≈ 4.75R C . This first peak located below r DPD is already present in the initial spherical configuration, and is linked to the difference between the NP chemistry. Indeed, Janus NP are made of two hemispherical faces, one hydrophobic in contact with the hydrocarbon beads, and one hydrophilic in contact with the water beads. The NP hydrophobic faces which protrude in the solvent interact strongly and repeal each other to minimize the interaction potential energy. Unlike Janus NPs, both hemispherical faces of homogeneous NPs are made with hydrophobic and hydrophilic beads. This allows the NPs to interact smoothly with each other and to fluctuate more. Hence they can come closer than r DPD , keeping spherical the droplet interface.

EVOLUTION OF THE THREE-PHASE CONTACT ANGLE DISTRIBUTION
In Fig. S2, we show the evolution of the three phase contact angle distribution of Janus (left panel) and homogeneous (right panel) NPs from the initial stage E 0 (cf. Fig. 1 in the main text), where the shape of the droplet is spherical, to the final stage E 20 . The initial distributions, fitted with continuous lines, can be described with Gaussian distributions for both the Janus and homogenous emulsion systems. The values of the respective means, µ J and µ H , and variances, σ J and σ H , differ according to the chemistry of the stabilizers. We obtain µ J = 91.6 • and µ H = 88.6 • , and σ J = 2.0 • and σ J = 3.4 • for Janus and homogeneous NPs, respectively.
When the droplet is coated with Janus NPs, the system evolves to a skewed distribution as the droplet shrinks, but remains unimodal with a single peak at the same value as the one measured for the spherical initial configuration. The emergence of the skewness of the distribution is linked to the decrease of the NP-NP distance when the droplet volume is reduced due to the major role played by the steric effect. The evolution is different when homogeneous NPs cover the droplet. As the volume is reduced, the contact angle distribution firstly evolves as a monolayer interface with a single peak (from E 0 to E 5 , Fig. 1 in the main text). As the distance between the NPs decreases further, the distribution becomes bimodal as two distinct peaks emerge on both sides of the original equilibrium contact angle. This fundamental difference is characteristic of a particle bilayer.
FIG. S3. Snapshots of the final configurations E20 of the buckling processes of water droplets armored with 160 spherical Janus (left) and homogeneous (right). The proportion of water beads removed is around 10 percent in each evaporation. The equilibration time after each evaporation is reduced to 2 × 10 3 timesteps (instead of 10 5 as in the main text).

HOW NUMERICAL ALGORITHMS AFFECT THE DROPLET EVOLUTION
In our mesoscopic analysis, we controllably reduced the volume of the droplet, removing a constant proportion, i.e. 10 percent, of water beads from the droplet. As a consequence, the three-phase contact angle distribution of the NPs evolves smoothly when the droplet buckles, thereby preventing particles to be artifactually released. However, the particles behaviour strongly depends on the numerical protocol.
In Fig. S3, we show qualitatively, for comparison, the evolution of the droplet volume when the equilibration time after each water bead removal is reduced. The proportion of water beads removed remains in 10 percent in each stage. When the droplet is coated with homogeneous NPs (right panel), we observe that the shape of the droplet remains spherical, with some NPs desorbed in the organic solvent. This evolution is representative of the passive role played by the homogeneous NPs and is in agreement with experiments [12]. When the droplet is coated with Janus NPs (left panel), we observe a significant curved-shape deformation of the droplet, along with the abscence of NP release. This evolution is representative of the active role played by the Janus NPs. The morphology of the droplet becomes noticeably crumpled, with large dimples, and no transition to crater-like depression is observed. This results are consistent with the surface model numerical analysis from Ref. [17], where more than one dimple may nucleate if the evaporation is rapid, leading to metastable multi-indented shapes. Experimentally, the term rapid may correspond to kinetic barriers, which prevent thermally activated coalescence between adjacent dimples. Our result highlights the central role played by the relaxation time of the system after each evaporation in the evolution of the interface geometry.