Can one predict a drop contact angle?

The study of wetting phenomena is of great interest due to the multifaceted technological applications of hydrophobic and hydrophilic surfaces. The theoretical approaches proposed by Wenzel and later by Cassie and Baxter to describe the behaviour of a droplet of water on a rough solid were extensively used and improved to characterize the apparent contact angle of a droplet. However, the equilibrium hypothesis implied in these models means that they are not always predictive of experimental contact angles due to strong metastabilities typically occurring on heterogeneous surfaces. A predictive scheme for contact angle is thus urgently needed both to characterise a surface by contact angle measurements and to design superhydrophobic and oleophobic surfaces with the desired properties, e.g., contact angle hysteresis. In this work a combination of Monte Carlo simulation and the string method is employed to calculate the free energy profile of a liquid droplet deposited on a pillared surface. For the analyzed surfaces, we show that there is only one minimum of the free energy that corresponds to the superhydrophobic wetting state while the wet state can present multiple minima. Furthermore, when the surface roughness decreases the amount of local minima observed in the free energy profile increases.


Introduction
Contact angle is perhaps the most immediate way to characterise the wettability of a given surface, due to the ease of depositing a drop and estimate the angle formed by it. Indeed, the Young equation allows one to relate the surface tensions with the contact angle on an ideally smooth surface [1]. However, the crucial question is whether this observable is a reliable one in the case of actual surfaces, in which heterogeneities are inevitably present.
Surface heterogeneities can dramatically change the contact angle as compared to the case of an ideal, homogeneous surface. In particular, topographical heterogeneities tend to magnify the intrinsic wettability of the surface, making them "superhydrophobic" or "superhydrophilic" [2]. Since the pioneering work of Wenzel [3] and Cassie and Baxter [4], the concept of apparent contact angle has been used, which is meant to predict the angle measured at the actual surface accounting for the effect of both wettability and heterogeneity. Since then, a continued effort has been spent to provide accurate models for predicting the apparent contact angles based on the chemical and topographical characteristics of the surface. The main ambition of this line of research is finding a rapid and inexpensive way to characterise surface wettabilities via contact angle measurements.
The standard approach to develop wetting models for heterogeneous surfaces consists in assuming a wetting state of the surface and predict the contact angle based on the homogenised surface energy of such composite interface. The task becomes more involved when a surface allows multiple visibly different states, e.g., for superhydrophobic surfaces in which the roughness can be fully wet (Wenzel, W) or dry (Cassie-Baxter, CB).
Importantly, major difficulties arise in relating the apparent contact angle and surface characteristics, most importantly due to static contact angle hysteresis (CAH), i.e., a scatter of the contact angle values which depends on the experimental procedure or history. Already in 1964, Johnson and Dettre [5] showed that an idealised sinusoidal roughness can be an important source of CAH; they further showed that the W state is connected to larger CAH than the CB state. More recent work further showed that CAH emerges from micro and nanoscale phenomena which are beyond the scope of Cassie and Wenzel equations [6]. A heated debate originated from the question whether "Cassie and Wenzel were wrong" [7,8,9]. The main merit of such discussion was to highlight that the experimental contact angle may be determined by features localised at the threephase contact line where the drop meets the surface. In fact, the classical CB and W models assume an equilibrium hypothesis in which the droplet is allowed to relax to the absolute free energy minimum. In other words, due to the rough free energy landscape featuring multiple minima, the measured contact angle may significantly depart from the value predicted based on averaged surface characteristics and the commonly used wetting models [10,11,12].
Several attempts have been made to relate the characteristics of single defects to CAH, starting from the pioneering work of Joanny and de Gennes [13,14,15,16]. However, densely packed heterogeneities [17], together with the finite size of the drop and the multiple wetting states possible on superhydrophobic surfaces [18,19,20,21,22,23] considerably complicates the matter. Indeed to date, notwithstanding the theoretical and practical importance, the analysis of the free energy landscape connected to wetting of a realistic surface by a droplet is still lacking. Thus making a connection between measured contact angles and available wetting models is highly needed in the community. At the same time, the difficulties connected to the exploration of rough free energy landscapes [24] represent a major obstacle to develop the computational tools required for engineering textured surfaces with the desired wetting properties. In this work, we fill such a gap, providing a generic framework for reconstructing the free energy connected to wetting of a pillared superhydrophobic surface based on the combination of the cellular Potts model [25] and the (zero temperature) string method for rare events [26,27,28]. This approach allows us to assess the variation of CAH with the pillar spacing, to connect CAH with the existence of multiple free energy minima, and to verify the energetic origin of the superhydrophobicity of the CB state. In addition, the presented approach provides useful guidelines for contact angle measurements on actual surfaces and for designing surfaces with tailored wetting properties, notably superhydrophobic ones.
The paper is organised as follows. In Section 2 we present the numerical setup used: a Monte Carlo scheme combined with the string method. Section 3 contains the results of the free energy calculations and the physical origin of its saddle points. In Section 4, we discuss the origin of contact angle hysteresis and compare the CB and W models with the simulations. The paper is then concluded in Section 5.

Numerical Model
The main aim of our simulation campaign is connecting the apparent contact angles attainable on a given surface with its physical and chemical characteristics. We focus on a hydrophobic surface decorated with a three-dimensional array of pillars with the geometry shown in Figure 1-a. A three-dimensional spherical droplet with volume V 0 = 4/3πR 3 0 is placed over the surface. We use a mesoscale wetting model, the cellular Potts model [25,29,30], which has been successfully used in studying the wetting of structured surfaces [31,32,33,34,35]. Its low computational cost as compared to, e.g., molecular dynamics, allows computations of relatively large 3D drops via Monte Carlo minimisation.
We combine the MC minimisation with the string method [26,28] in order to calculate the most probable path leading from the superhydrophobic CB to the fully wet W state or vice versa. The string method, combined with atomistic simulations, diffuse interface methods, or sharp interface models, has proven instrumental to clarify the origin of metastability in wetting of structured surfaces [23,36,37,38]. We will employ this approach to overcome the free energy barriers, to establish the mechanism of breakdown/recovery of the CB state, and to identify intermediate minima.

The cellular Potts model
The Cellular Potts model (CPM) assumes a three state system in a simple cubic lattice, with each site representing a different state: solid, liquid, and gas. We employ this technique to simulate a three-dimensional droplet in contact with air and a solid surface. The Hamiltonian of the system is given by: where s i ∈ {0, 1, 2} is the state of a given site in the lattice, gas (G), liquid (L), or solid (S). The first term on the right hand side describes the surface interaction between the three states, with E s i ,s j = J ij · a L , where J ij is the site-site coupling between i and j type sites, a L = L 2 is the unit area and L is the lattice size. The summation ranges over the n neigh = 26 nearest and next-nearest lattice sites of the 3D simple cubic lattice. Notice that the energy contribution will be non-zero only when s i = s j . Gravitational contributions to the energy are neglected due to small size of the droplets hereby investigated [39,34,33]. The second term in Equation (1) corresponds to a harmonic bias used to drive the system across metastabilities allowing to more effectively explore phase space [27] using general descriptors, the collective variables z, which are function of the lattice sites. Biasing is used because, in an unrestrained simulation (κ = 0), drops relax and remain trapped in the closest minimum state. Simulation modeling a droplet deposited on a textured surface without any bias were previously studied in Refs. [33,40], where the authors indeed showed that the final contact angle depends on the initial wetting configuration.
The collective variables z that appear in Equation (1) should be chosen carefully [41] in order to correctly discriminate among relevant configurations of the system. More in detail z T represents the target value for the collective variable, while κ is stiffness of the harmonic constraint. In this work we shall identify z l with liquid-phase occupation of N interpillar grooves. The usage of such discrete density indicators is customary when dealing with wettability problems [42].
The energy scale is set once J ij is defined. In this work spin-spin couplings are chosen to render the same ratios between the interfacial surface tensions as in Ref. [43], in which a droplet of water is studied on a surface composed by Polydimethylsiloxane (PMDS) with the following surface tensions σ ij between the gas (G), liquid (L), and solid (S): σ GL = 70 × 10 −3 N/m, σ SG = 25 × 10 −3 N/m, and σ SL = 53 × 10 −3 N/m; the last value was calculated by Young's equation, σ GL cos θ Y = σ SG − σ SL with θ Y = 114 • . The data will be presented in dimensionless form using the lattice size L as unit length, a L = L 2 as unit area, and J GL · a L as unit energy, referred as e.u. In such reduced units, J GL = 1, J SG = 0.36, and J SL = 0.76. For the forcing parameter we use κ 0.004, which drives the system to the target configuration while allowing some fluctuation. More details about the choice of κ is reported in the Supporting Information (SI). In order to convert the dimensionless units to physical ones, one can for instance assume L = 1 µm and the energy for the gas-liquid interface of water, which yields σ GL a L = 7 × 10 −14 J.
To evolve the system we use the Metropolis-Hastings algorithm, which consists in changing the state of two random sites at the gas-liquid interface with an acceptance rate equal to min{1, exp[−β∆H]}, where β = 1/T is the inverse of the effective temperature of the CPM, which introduces some noise, allowing for a more effective exploration of the phase space. We set T = 4.8 throughout the article, which allows the system to fluctuate with an acceptance rate of approximately 9%. The attempted MC moves consist in swaps between liquid and gas sites, which guarantees that the volume of the droplet is constant throughout the simulation.

The string method
In order to attempt a simulation of the wetting of rough surfaces it is crucial to tackle the challenge associated with the presence of metastabilities. Previous methods to sample rough free energy landscapes required projecting to a low dimensional representation of the free energy as a function of a handful of parameters and reconstructing the full landscape as a function of such collective parameters. It is easy to understand how this task requires a computational effort which is exponential in the number of collective variables thus limiting the applicability of this class of methods typically to two/three variables. Such low dimensional representation of the energetic landscape often results in a poor description of the phenomena [41]. The (zero temperature) string method in collective variables was first introduced by Vanden-Eijnden and collaborators [27] and is a path method that requires only the computation of the local gradient of the free energy landscape at certain points along the path; its computational cost thus only grows linearly with the number of collective variables.
The string method allows for a fast and convenient identification of the minimum free energy path (MEP) connecting two metastable regions of the rough landscape, along with providing a measure of the free energy along the path, without requiring to fully sample a high dimensional variable space. This is achieved by iteratively refining a discretized guess path (i.e., the string). In this work the collective variables (CV) represent the liquid occupation of each cavity, as represented by the blue region in Figure 1-b. The available volume in the cavity is given by (d 2 − w 2 )h. The algorithm proceeds as follows: 1. We run N r Monte Carlo (MC) replicas, each biased to explore the vicinity of a z T l,i point in the collective variable space. Each replica corresponds a string point in the collective variable space.
2. The N r replicas run for 10 5 MC moves allowing to sample the mean biasing forces and the local free energy gradient at each string point.
3. The string is updated by a convenient gradient descent of the string points and reparametrized [28] so that one has a new set of z T l,i . The process is iterated from point 1.
Refinement of the initial string guess is repeated until convergence is reached, as evaluated by a suitable threshold in the changes to the path/path energy. In all runs convergence was obtained within 20 iterations. At convergence, numerical integration of the free energy gradients allows for the reconstruction of the free energy profile along the string: with ∆F j the free energy at point z T l,j , the index l running over the N collective variables, i up to the current replica j, with j ≤ N r − 1, and ∆z T l,i = z T l,i+1 − z T l,i . The results shown in this work are averages over the distinct realizations of the simulation and . . . represents the average over MC steps.

Results
We consider three substrates with same pillar height h = 10L and width w = 5L and different distances between pillars: substrate referred to as S 1 presents a low interpillar distance (a = 5L), S 2 an intermediate value (a = 8L), while S 3 has the largest interpillar distance (a = 11L). The volume of the spherical droplet is defined by imposing an initial radius of R 0 = 50L. We first show the rough free energy landscape connected with wetting of the three surfaces, characterizing the droplet configuration at the minima and maxima of the free energy. We then introduce a possible physical explanation of the roughness in the free energy and end this section discussing the connection between the free energy landscape and contact angle hysteresis. Figure 2-a reports, for substrates S 1 , S 2 , and S 3 , the free energy of the drop as a function of the total filling of the cavities below the drop, which is defined as the sum of the target filling of each cavity divided by the total volume of the droplet: f = l z T l /V 0 . The minima in the free energy correspond to stable (global) and metastable (local) states of the system. The free energy extrema are numerically identified from the data; minima are indicated in Figure 2-a as black circles, while maxima by red triangles. Since some minima are very shallow and may be affected by numerical accuracy we performed numerical MC minimisation using different initial conditions close to the putative minima, to confirm that the system can indeed reside in such states; empty symbols denote shallower minima which will not be further discussed. Figure 2-b shows the fraction of water wetting the bottom area. The vertical dotted line in the figure roughly separates two regimes: on the left, the liquid does not touch the bottom of the substrate, indicating that the corresponding configurations of the droplet are associated to the superhydrophobic CB state, while, on the right of the line, the liquid reaches to the bottom of the substrate and the droplet configurations are associated with the wet W state(s). Importantly, the free energy landscape appears rough for all considered surfaces, with markedly different trends for the three interpillar distances: monotonically growing for S 1 , almost at the CB-W coexistence, with additional high free energy minima, for S 2 , and with a significant W basin with multiple local minima for S 3 .

Fraction of bottom
area wetted  Each point on the free energy profiles corresponds to a droplet configuration, whose sequence thus defines a wetting path, see Figure 2-c; in particular, these paths represent the most probable way in which the transition from CB to one of the various W states (or vice versa) occurs. Other paths may exist connecting minima, especially in such complex landscape, see e.g. Ref. [44]. Figure 2-c presents a lateral view of the 3D droplet in correspondence of the minima and maxima of substrate S 3 . Minimum I corresponds to the CB state, with air trapped between the pillars -this numbering is the same for the three substrates. As the droplet starts infiltrating the substrate, it touches the bottom for the first time at point II, where free energy increases to a maximum. A similar behaviour was also reported in molecular dynamics simulations [23,45,46,47]. As the filling level increases, the free energy of the substrate presents several local minima which correspond to the progressive infiltration of liquid within the pillars: lateral views of the droplet in Figure 2-c clearly show full wetting of 3 lines of cavities at minimum III, 4 lines at V, 5 lines at VII. For filling levels above 7 lines of cavities, the free energy does not present any other minima for this droplet size.
The maxima separating the mentioned minima are found to be associated to an incomplete wetting of some of the cavities at the drop perimeter; we will further analyse their origin in Section 3.2. We note that the wetting of the substrate in 3D is a more complex problem than what can be inferred from a lateral view, which is however convenient to picture the main features of the process; we thus make available in the Supporting Information videos in 2D and 3D of the droplet wetting the 3 substrates.
The inset of Figure 2-a defines the barrier H L , which is the difference in ∆F between a minimum and the consecutive maximum on its left and H R , being the difference in ∆F between a minimum and the first maximum on its right. Using these definitions, we found that the barrier H R and H L are typically of the same order of magnitude for substrate S 3 , while for S 1 H L H R . We will come back to this point later. These observations raise several questions, which are addressed in the next section: why some substrates present global minima at the CB state and others in the W configuration? What is the physical origin and significance of the local minima and the intervening maxima?

Physical origin of the minima and maxima of the free energy
The wetting state corresponding to the global free energy minimum is different for the considered substrates and can be rationalized using a continuous global energy model [33]. This model takes into account the energy of creating interfaces between the liquid, the gas, and the solid of a spherical droplet with fixed volume V 0 placed on a surface. The droplet is allowed to display two wetting states: it either stays in the superhydrophobic CB state or it homogeneously wets the surface in the W state. We refer to these configurations as ideal Cassie-Baxter, CB I and ideal Wenzel, W I , states. The difference in energy of a system with and without the droplet on the surface can be written as [33]: where S s = 2πR s2 [1 − cos(θ s C )] is the surface of the spherical cap in contact with air, B s = R s sin(θ s C ) is the base radius, R s the radius of the droplet, and θ s C its contact angle in the state s. φ S = w 2 /d 2 is the fraction of solid surface area wet by the liquid (or pillar density) and r = 1 + 4wh/d 2 is the surface roughness ratio.
The main difference between the model defined by Equation (4) and the Wenzel one [3] is that the former one considers a droplet of finite size, which means that the interface of the cap and the lower interface in contact with the substrate compete in the minimisation of the total free energy. We will return to this point in Section 4. Concerning the CB state, it is instead found that results are equivalent when one explicitly considers the cap as in Equation (3) or simply homogenises the surface free energy as in the original CB theory [4]. are not in perfect agreement for the W state, due to the local pinning at pillars, which is not captured in Equation (4), see below for details. On the other hand, the model predicts correctly that W is the global free energy minimum and the free energy values are in reasonable accord, see the Supporting Information. Besides the fair agreement at the global minimum, the main difference between the two approaches is that the simulated free energy profile presents several local minima which are not accounted for in any homogenised model. This difference is particularly important for the W state, which helps to explain the difference in the global W minimum. Indeed, an important simplification introduced in the model (4) is that the W state is achieved by homogeneous wetting: below the droplet base, there is a perfect cylinder filling the cavities; individual wetting of pillars is thus disregarded; the cost to create the interface between the liquid and air is also neglected in the model. However, simulations show that the infiltration of the substrate is not homogeneous, that pinning at individual pillars may occur, and that the interface between the liquid and the gas below the droplet is quite rough, as exemplified in Figure 3-c for one minimum and one maximum: the configurations are particularly far from being a cylinder when ∆F is at a maximum.
The nontrivial shape of the droplet in contact with the surface drove us to investigate the contribution to the free energy of each interface. We propose a putative free energy defined as: where ∆P can be interpreted as the Laplace Pressure, ∆P = −ασ GL with α being the mean curvature of the drop and A ij are the interfaces between two different phases ij (solid, liquid, and gas). We then use the Young equation σ SG = σ SL + σ GL cos θ Y to rearrange the terms and note that the total area of the substrate A tot = A LS + A VS is constant. The contribution of liquid-gas surface area A GL is split in two parts, one corresponding to the surface of the spherical cap A C GL , and other to the part below the droplet, A B GL . We can then write the difference in free energy with respect to the reference one Ω ref = σ SG A tot : Interestingly, the surface tension σ GL factorises and only geometrical quantities appear in the parenthesis on the right hand side, together with the contact angle θ Y . Figure 4-a compares, for surface S 3 , the free energy ∆F obtained in the MC simulations via Equation (2) with the putative free energy ∆Ω always computed from MC simulations but by measuring the geometrical quantities in Equation (5). In particular, the surface areas A ij and the mean curvature α = 2/R are measured from droplet configurations along the converged string. Since the free energy ∆F is computed by integration, it is known up to a constant and thus can be freely displaced along the y-axis. The only free parameter left to compare ∆Ω to ∆F is σ GL , which should be of the order of the coupling parameter J ij . Figure 4-a shows ∆F and ∆Ω using σ GL = 1, which is the same value used for J ij . The two quantities show similar trends and are in semi-quantitative agreement for all three substrates analysed in this work (see the Supporting Information for the other two substrates). Interestingly, the ruggedness observed in ∆F is also present in ∆Ω.
In Figure 4-b each geometrical term in Equation (5) is shown separately and compared with ∆F. Note the different scales and units for areas (right axis) and ∆F (left axis). Most of the interfacial areas show a smooth dependence on f , the exception being A B GL that presents some maximum and minima which is thus a good candidate to explain the roughness of the free energy connected to wetting. In Figure 4-c we show A B GL and ∆F to assess the importance of the former in the free energy. Vertical dashed lines correspond to maxima of ∆F. We observe variations in the curve of A B GL that are indicated by vertical red lines. The typical size of these variations is ∆A B GL ≈ 1000L 2 . In terms of dimensionless energy, this corresponds to ∆A B GL σ GL ≈ 1000. This energy variation is comparable to the free energy barrier for S 3 , ∆A B GL σ GL ≈ H R , which suggests that the term relative to the interface between liquid and gas below the droplet plays an important role in generating local minima in the free energy. For S 1 , instead, abrupt variations of A B GL are not sufficient to generate local minima, due to the steep slope of ∆F vs f , see the Supporting Information.
In the previous section, we mentioned that local minima of ∆F correspond to pinning of the drop at the pillars edges, which gives rise to several possible minimal configurations characterised by different numbers of pillars. To connect this picture with the observation that variations of A B GL correlate with variations of ∆F, we suppose that the wet domain below the droplet can be approximated by a cylinder of height h and increasing radius B, such that A B GL = 2πhB. This approximation does not take into account the roughness of A B GL shown in Figure 3-c, but is reasonable in the case of minima. From the volume differences between neighbouring minima in Figure 4 we can thus compute the jump in droplet radius ∆B in the cylindrical approximation, ∆B ≈ ∆A B GL /(2πh) ≈ 16L. The estimated value of ∆B corresponds to the typical size of a cavity d = w + a, which is plausible and supports the idea that local minima correspond to jumps of the droplet front across discrete numbers of pillars.
To summarize, we propose that the global minimum of the free energy corresponds to configurations that minimise the total interfacial energy for a droplet of fixed volume. Local minima, instead, occur in correspondence of abrupt variations of the liquid-vapor interface below the droplet connected with the overcoming of individual pillars.

Minima of the free energy and contact angle hysteresis
We now evaluate the free energy barrier sizes that separate the local minima and the consequences on the metastability of the substrates. In Figure 2-d, the barriers H R and H L are defined. From Figure 2-a, one measures for S 1 typically H R ∈ (3000, 9000) and H L ∈ (0, 100) showing that, even when the system is initialised in a W-like state, it rapidly evolves towards the CB minimum. For substrate S 3 both barriers vary typically in the range H R ≈ H L ∈ (370, 3700), which in physical units are between 9 × 10 −12 J and 9 × 10 −11 J, assuming L = 1 µm. This is much higher than the thermal energy k B T ≈ 4.1 × 10 −21 J at ambient temperature: thermal fluctuations are not sufficient to drive the system from one minimum to the other. Only other larger sources of energy can move the drop away from local minima, e.g., mechanical vibrations.
When the system is prepared with some generic initial condition, it will fall in the closest minimum and remain trapped there; this was verified by running unrestrained MC simulations from several points along the curves in Figure 2-a. Only for the S 1 the system always returned to the superhydrophobic CB state, even though two shallow minima were identified numerically; this can either be due to the numerical accuracy of the free-energy profiles or to the size of the barriers, which is so low that the fluctuation imposed by the effective temperature of the Monte Carlo simulations are enough to bypass them. In other words, the free energy profiles in Figure 2 help understanding the origin of contact angle hysteresis in contact angle measurements, which can be rationalised in terms of a rough free energy landscape with multiple minima, separated by large free energy barriers. Furthermore, the importance of the initial conditions becomes apparent, which correspond to the preparation of the sessile drop experiment.
Finally, the profiles account for the superhydrophobic properties of the CB state, which are connected to the existence of a single minimum, i.e., with low hysteresis. On the other hand, the presence of multiple W minima explains why CAH is so pronounced in this wetting state and its stickiness [6,49,48]. Figure 5 shows the apparent contact angle θ C of the droplet measured in MC simulations together with ∆F for substrate S 3 . The fact that there is a basin of Wenzel states and that each minimum has a different contact angle allows us to define the contact angle hysteresis θ H related to the wet W state as the difference in θ C of the configuration associated to the first and the last local minimum of ∆F, see Figure 5. Table 1 summarizes the values θ H for the three substrates together with their roughness ratio r. We measured an increase of θ H when r decreases, which is in line to what was previously observed [50].

Discussion: modeling rough wetting
We have identified by free energy simulations that some pillared substrates present several local minima separated by high barriers, while others present only one minimum. Incidentally, for tall/tightly packed pillars only the superhydrophobic CB state is possible, while for more sparse pillars multiple local wet minima arise. The latter surfaces thus display a behavior which is strongly dependent on the initial conditions: if a droplet is deposited on a substrate at a random configuration, it would accommodate in a wetting state correspondent to the closest minimum which can be either superhydrophobic or sticky. The present findings also imply that theories that predict a single W state cannot be complete [2]. We measured all apparent contact angles θ C of the droplet in configurations correspondent to physical minima identified in our simulations for the three substrates and compared them with theoretical models, see Figure 6. Squares corresponds to the superhydrophobic CB minimum and circles to wet W states. Lines are solutions of the classical Cassie-Baxter and Wenzel models, whose apparent contact angles are given by: For the case where the global minimum is CB, which happens for substrate S 1 , Figure 6 shows that the prediction of the Cassie-Baxter model is almost quantitative. The CB model also has a reasonable agreement with simulated θ C in cases where it is only a local minimum, which happens for substrates S 2 and S 3 . The agreement with CB model deteriorates as the pillar distance increases, which can be explained by the increasing curvature of the menisci suspended among pillars, not accounted for in the classical CB model.
In principle, no direct comparison can be made between the W model and the MC results, mainly because the model predicts only one state, while simulations demonstrate the existence of multiple wet states. However, a fair comparison can be made considering only the global minimum, which is compatible with the minimisation procedure used in all homogenised models. It is seen that the W model prediction is far from the measured contact angle and, for S 2 , it even predicts a contact angle higher than the actual CB one. When the finite size of the droplet is taken into account, as done in Equation (3) and (4), the solution for the CB does not change but the W curve shifts to smaller values. This simple correction captures the overall trend with interpillar distance and improves the agreement with the global minima, although it is not quantitative. However, we remark that the finite size is not enough to account for the existence of multiple minima, which is due to pinning at individual pillars and is crucial to account for contact angle hysteresis. Figure 6 confirms that the extent of CAH has the following trend S 1 < S 2 < S 3 , which suggests that the more favorable the W state, the higher θ H . This trend is related to the number of local minima and to the facility of wetting the bottom of the surface for short pillars, but what limits the total number of minima still remains an open question. Continuos mode l, Eq.3 Continuos mode l, Eq.4 CB mode l, Eq.6 W mode l, Eq.7 Homogenised CB Homogenised W CB model W model Figure 6: Apparent contact angle computed at physical minima in Figure 2

Conclusions
In this work, we have used a combination of the cellular Potts model [25,33] and the string method [26,27,28] to compute numerically the free energy ∆F of a droplet with a fixed volume placed on surfaces decorated with pillars. The computed free energy landscape for wetting the surface is found to be rough, with one minimum corresponding to the superhydrophobic Cassie-Baxter state; depending on the pillar spacing, multiple local minima can exist, associated with different wet Wenzel states. Free energy barriers typically larger than the thermal energy are found, whose origin can be traced back to pinning of the liquid-gas interface below the droplet at individual pillars. Wenzel minima are characterized by an increasing number of filled cavities, corresponding to different apparent contact angles. This scenario accounts for the strong contact angle hysteresis of the W state(s), in which the final wetting state depends on the initial condition.
We have compared the apparent contact angles obtained in our in silico experiments with the predictions of simple models, including the classical Cassie-Baxter [4] and Wenzel [3] ones. Results showed that the Cassie-Baxter model has a good prediction capacity, which can be improved by considering the curvature of menisci overhanging surface textures. On the other hand, the prediction of the apparent contact angle cannot be made by simple models in the W state effectively hindering the deduction of surface characteristics from θ C measurements. Actually, the problem of predicting θ C in the W state seems ill-posed and the relevant challenge is to assess its contact angle hysteresis, which can be achieved only by detailed models.
An interesting and long-standing question is the possibility of predicting CAH [13,14,16], which was possible for individual surfaces within our approach. It is found that the wider the pillars are spaced, the larger CAH. Our analysis suggests that larger drops would have liquid enough to fill more cavities and thus the free-energy profile would show even more local minima. What limits the total number of filled cavities and how these local minima connect to CAH is however not clear and deserves to be investigated.
Finally, our findings suggest that, for a complete characterization of the wetting properties of a substrate, it is necessary to repeat sessile drop experiments several times for different initial conditions. Moreover, the fact that free energy barrier between minima are typically much larger than thermal fluctuations suggests that mechanical vibrations of the substrate are necessary to drive the droplet across different minima. The amount of mechanical energy necessary to reach different minima need to be evaluated by computing the barrier sizes, which is possible with the methodology proposed in this work.

Supporting Information
The supporting information describes the algorithm used to find the thermodynamic wetting state for droplet placed on the pillared surface. Its also presents the criteria to adjust parameter κ and description of the videos available. Then there is a discussion of the comparison between the free energy obtained from the simulation with the Cassie-Baxter and Wenzel models. Finally the areas of the components for the proposed analytical free energy are presented.
under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 803213).