Bubbles determine the amount of alcohol in Mezcal

Mezcal is a traditional Mexican spirit, obtained from the distillation of fermented agave juices. Its preparation has been conducted for centuries in an artisanal manner. The method used to determine the correct alcohol content is of particular interest: a stream of the liquor is poured into a small vessel to induce surface bubbles. These bubbles, known as pearls by the Mezcal artisans, remain stable for tenths of seconds only if the alcohol content is close to 50%. For higher or lower alcohol content, the bubbles burst rapidly. The long bubble lifetime is the result of surfactant-induced surface tension changes. However, the precise mechanism and its relation to alcohol content remain unexplained. In this investigation, the extended lifetime of pearls was studied both experimentally and numerically. It was found that changes in surface tension, density, viscosity (resulting from mixing ethanol and water), and the presence of surfactants are all relevant to extend the bubble lifetime. The dimensionless bubble lifetime was found to reach its maximum value when the Bond number was close to unity, corresponding to 2 mm Mezcal bubbles. These findings show that the traditional empirical method does work. Beyond this, the understanding of the process provides physical insight to many other natural and industrial problems for which the stability of surface bubbles is of importance, such as bio-foams, froth floatation, and volcanic flows.

Physical mechanisms that produce surface tension gradients. Three physical mechanisms may induce gradients of surface tension: surfactant transport, heterogenous evaporation, and local changes in temperature. All of these mechanisms affect the time the film takes to drain and, therefore, bubble lifetime. Let us first consider the effect of alcohol concentration. The addition of 10% of mass of ethanol to water, for example, decreases the surface tension of the mixture more than 30% as seen in the study by Vazquez et al. 26 , where variations in temperature were shown to have a small effect -an increase of 30 C resulted in a surface tension reduction of less than 6%.
In ethanol-water mixtures, the alcohol evaporates at a constant rate regardless of the initial concentration 27 . Therefore, it would be expected for the evaporation rate to be uniform along the film, avoiding gradients of surface tension. However, the evaporation rate does vary with the thickness of the film as in Mezcal pearls. Such non-uniform evaporation naturally results in a decrease of bubble lifetime. Hence, we do not discard the effect of evaporation. However, given the long lifetimes observed experimentally, we consider this effect as secondary. We can thus argue that surface tension gradients result primarily from the changes in the surfactant concentration. To fully resolve this issue, a numerical analysis would have to account for both evaporation and thermal effects. Such a scheme is beyond the scope of the present work. Instead, based on the arguments above, we account only for the surfactant-induced Marangoni effect. The details are explained in the following sections. www.nature.com/scientificreports www.nature.com/scientificreports/ The aim of our investigation is to demonstrate that the traditional technique to assess the ethanol content in Mezcal actually works and to explain the relationship between the lifetime of pearls and alcohol content. Although modern methods of determining alcohol content are accessible and reliable, even for small remote rural communities, the fundamental understanding of the physical mechanisms behind bubble and foam stability is relevant due to their importance to many other fields. In this study, we conducted experiments to determine the lifetime of surface bubbles of many different types of Mezcal and water-ethanol mixtures. To gain insight into the film drainage process, we also conducted numerical simulations and analyzed all results in terms of dimensionless quantities.
Dimensional analysis. As proposed by Barenblatt 28 , dimensional analysis can be used to understand a physical phenomenon more deeply. There are several ways in which the relevant dimensionless numbers can be identified. We consider two alternatives below.
Dimensionless groups by simple inspection. Drainage results from a balance between viscous forces and either gravitational or surface tension forces. The ratio that determines the relative importance between gravitational and surface tension effects is the Bond number, 2 where ρ is the liquid density, g is the gravitational acceleration, D is the bubble diameter, and σ is the surface tension. All of these quantities are attainable via laboratory experiments. Here, D is determined by analyzing images of bubbles from a top view (Methods Section). It should be noted that this metric gives an apparent bubble size rather than an equivalent bubble diameter, which measures the diameter of a spherical bubble with equivalent volume.
Considering surface tension effects, the bubble lifetime can be scaled by incorporating viscous effects leading to: life life where μ is the liquid viscosity. One can propose the following functional relationship between ⁎ T life and Bo: life Now, as explained above, the ever-presence of surfactants plays an important role in the film drainage process. If one defines Δσ as the change in surface tension resulting from surfactant gradients along the interface, and the relative variation as one can therefore argue that Here, we aim to understand the effect of each group by considering that In other words, we shall consider that the effect of surfactants can be de-correlated from the drainage. This fact is not obvious and represents an important assumption of our study. However, as discussed below, the assumption appears to be reasonable. Furthermore, the analysis indicates that the surfactant-induced (Marangoni) stresses should be considered to determine an appropriate timescale.
Dimensionless groups from the analysis of the flow type. The paper by Champougny et al. 23 discusses the unsolved problem of bubble drainage in the presence of surfactants. They propose an ad-hoc model based on an extrapolation length λ that allows describing the transition from stress-free interfaces (λ = ∞) to no-slip interfaces (λ = 0). The flow between stress-free interfaces is a plug flow at leading-order, dominated by extensional viscous stresses. Debregeas et al. 29 have shown that the film thickness of such a 'bare' bubble decays exponentially in time and that the film always ruptures at the apex. Contrarily, the flow between no-slip (rigid) interfaces is a Poiseuille-like flow dominated by viscous shear stresses (as described, for instance, in 30 ). Now, surfactant-induced Marangoni stresses have been shown to rigidify, at least partially, the interfaces, such that the flow is essentially a shear flow (e.g., Champougny et al. 31 for film formation or Atasi et al. 32 for confined bubbles in microchannels, both in the presence of surfactants). Even traces of surfactants, such as impurities, have been found to induce dominant Marangoni stresses 33 . Additionally, local thinning mechanisms of the film have been observed in soap bubbles 18 , reminiscent to marginal regeneration 34 , and eventually leading to the rupture of the film away from the apex. Based on these arguments, we conjecture that viscous shear stresses dominate the film drainage. This hypothesis is corroborated below via numerical simulations. www.nature.com/scientificreports www.nature.com/scientificreports/ In absence of a simple model available that accounts for the presence of surfactants and with the aim to find an appropriate timescale, let us consider a drainage flow dominated by viscous shear, balanced either by gravity or capillarity. Starting with a gravity-driven drainage, the momentum balance at leading order in the lubrication approximation writes yy where the cross-stream coordinate y in the film can be scaled with the critical thickness for rupture h rup . The choice of this length scale is justified by the fact that the drainage dynamics is not significantly influenced by the initial condition since the main contribution to the lifetime occurs at the later stage, i.e. when the film is the thinnest, which is also when the viscous dissipation dominates. Note, however, that h rup is not a quantity that can be obtained directly from the experiments. By considering the speed at which the film retracts after it ruptures, as explained in the Methods Section below, an estimate of this thickness can be inferred. For simplicity, we consider that the rupture always occurs at the same thickness value. Considering that u ~ D/t g , the timescale for drainage driven by gravity can be obtained as: Similarly, the momentum balance for a capillary-driven drainage is given by where the streamwise coordinate x can be scaled with D. If we assume that the pressure P results from a Laplace effect, we can write P ~ σ/D. By estimating that in this case u ~ D/t c , the timescale for drainage driven by capillarity is given by Comparing these two timescales leads to the Bond number, namely = Bo t t / c g , exactly as defined in Eq. (1). It is interesting to note that the Bond number plays two roles in this problem, a static and a dynamic one. On the one hand, it allows evaluating the static shape of the interfaces, namely a spherical bubble underneath an almost undeformed liquid surface for  Bo 1, and a deformed bubble under a deformed liquid surface for  Bo 1. On the other hand, as explained above, it provides a way to evaluate the dominant driving force for drainage, namely the capillary force for  Bo 1 and the gravity force for  Bo 1. The striking feature of the present problem is that ∼ Bo 1 indicates a transition which, as shown below, coincides with the maximum dimensionless bubble lifetime.
A different dimensionless lifetime can, therefore, be defined in terms of the capillary timescale, t c : where ⁎ T life is defined in Eq. (2). We also introduce ε =  h D / 1 rup , a measure of the slenderness of the film. Equation (11) indicates that the viscous forces act perpendicularly to the capillary/gravity forces, as reminiscent to shear-dominated flow. Conversely, forces in the case of an extensional flow would act along in the flow direction affecting the drainage; for such a case ε would be close to unity. Consequently, ⁎ T life would be the appropriate group for an extensional drainage. In the present problem, ε − 10 2 4 , therefore ⁎⁎ T life is more appropriate in the present context.
As will be shown later, the data shows a clear transition in trend at around Bo ~ 1. Using the ⁎⁎ T life scale, we can derive the trends in both limits of large and small Bond numbers. In the limit of large Bond number the drainage is governed by gravity, therefore, using Eq. (8), one can write On the other hand, when the Bond number is small the drainage is dominated by capillarity, which, considering Eq. (10), implies that life ⁎⁎  These two trends will be contrasted against the experimental results below.
Review of relevant models to predict the bubble lifetime. As stated in the Introduction, many studies have addressed the issue of bubble lifetime. Here, we summarize some relevant studies. Note that we do not aim to comprehensively summarize the state-of-the-art of the subject. To draw comparisons with our experiments and among the different models, we recast the predictions in terms of the dimensionless groups discussed above and recognize that for a given liquid, varying the Bond number amounts to varying the bubble size.
www.nature.com/scientificreports www.nature.com/scientificreports/ Large bubbles. Bo ≫ 1 In the limit of large Bond numbers, the bubbles are large and the film near the apex is essentially uniform and the drainage is driven by gravity.
For the case of rigid interfaces, an analytical solution exists for the dimensionless lifetime as follows 35 Notably, this prediction underestimates most of the experimental times shown later in the Results section, which indicates that the Marangoni stress in the experiments is large enough to retard the drainage longer than what rigid interfaces would do. This is only possible if the mean flow is upwards, i.e. towards the apex, entrained by Marangoni stresses. Such gradient-tension-induced stresses would need to built up prior to the drainage phase 36 . The simulations shown below show this type of behavior.
In the case of stress-free interfaces, Kocarkova et al. 37 calculated the evolution of the film thickness considering an extensional flow in the film, assumed to be uniform and axisymmetric. They observed an exponential thinning, from which the dimensionless lifetime can be shown to be: Champougny et al. 23 extended the analysis of Kocarcova et al. 37 considering certain levels of surface rigidification of the interfaces to account for the presence of surfactants. They also found an exponentially thinning of the film with time, which leads to the same functional relation of Eq. (15) but the proportionality constant was larger in the case of surfactants. Interestingly, they also found that the puncture of the bubble film changed from the apex to the foot of the bubble as the amount of surfactants increased.
Finally, Lhuissier and Villermaux 18 considered the case of bubbles in water. They argued that the scaling of a viscous drainage was not appropriate for water due to the phenomenon known as marginal regeneration 38 . In short, the balance of capillary pressure from the film curvature, the meniscus at the foot of the bubble and surface tension gradients lead to a non-uniform film thickness which causes the appearance of a localized pinching. Moreover, Lhuissier and Villermaux 18 recognized that the film breakup could also be influenced by the Bénard-Marangoni convection flows within the film. Considering the fluctuations from marginal regeneration convection cells and probabilistic arguments of the puncture breakup mechanism, they found which also indicates a decreasing lifetime with an increase of Bo, as qualitatively obtained in Eq. (12).
Small bubbles. Bo ≪ 1 In the limit of small Bond numbers, no simple model exists for the bubble lifetime because in this limit the drainage is driven by capillary forces which also induce surface deformations. The model of Howell 39 , considers this case but assumes an extensional flow, which leads to where k 1 is a constant that depends on the initial and rupture thicknesses. This expression indicates that the rupture time increases with Bo: the larger the bubble, the longer it will take to burst. This calculation assumes that the film drains uniformly and axisymmetrically without surfactants. The small Bo limit indicates that the bubble is nearly spherical and is mostly immersed below the liquid free surface. Note that this trend is different from that obtained considering a shear flow (Eq. 13).

Transition-sized bubbles. Bo ~ 1
From the brief review shown above, an important conclusion can be reached: the trends are opposite for small and large Bond numbers. This fact indicates that there is a transition at a certain value of Bo. From our arguments, confirmed below by the experiments, the critical Bo-value is around unity, where capillary and gravity forces are of comparable magnitude to drive the drainage. And this transition corresponds to a maximum dimensionless lifetime, as suggested by the cross-over between the models given above.
experimental Results pearl formation. Figure 2(a) shows an image of the formation process of pearls (see also Video S2, Supplementary Information section), considering the controlled experiment described in the Methods Section. The fluid jet fragments into droplets, which continuously splash against the surface of the liquid. The continuous splashing induces the formation of air cavities that, in turn, form bubbles that eventually rise to the surface. However, only bubbles of a certain size (of approximately 2 mm in diameter) remain on the surface for longer times.
The process of air entrainment from a plunging jet has been studied in detail [40][41][42] , due to its relevance for many flow phenomena such as aeration of the ocean and other water reservoirs. It is believed to be well understood. Hence, we did not pursue a more complete study of the formation process. Instead, we focus on the residence time of bubbles on the surface. Figure 2(b) shows a time sequence of a typical Mezcal bubble at the moment of bursting (see also Video S3, Supplementary Information section), produced with the system described in the Methods Section. It was observed that bursting initiates with a small hole rupturing the film, which opens quickly until it www.nature.com/scientificreports www.nature.com/scientificreports/ retracts completely. The puncture does not typically appear at the apex of the bubble; the thinning of the film is not uniform arguably due to the partial regeneration mechanism 18 . Several experiments were conducted under the same nominal conditions (same fluid and same bubble diameter) to calculate the typical value of the rupture thickness, according to the scheme described in the Methods  www.nature.com/scientificreports www.nature.com/scientificreports/ Section. It was found that the thickness was h rup = 23.5 ± 4.8 μm for the image shown in Fig. 2(b). This value is relatively large, compared to what has been measured for the case of water and seawater 18 . This value was used as a reference for the numerical simulations and is also discussed in the scaling of the bubble lifetime.

Bubble lifetime.
Experiments were initially conducted using Mezcal with the 'correct' ethanol volume fraction, identified as M1 (see Table 1). The mean lifetime in this case was 28.1 ± 12.5 s, depicted in Fig. 3(a) by the filled red circle. To vary the amount of ethanol in this liquid, either pure water or ethanol was added. The measured lifetimes for this ethanol-adjusted Mezcal are shown Fig. 3(a), denoted as liquid MA1 from Table 1. The figure shows the measured bubble lifetime as a function of ethanol volume fraction. As the amount of ethanol increases or decreases, the bubble lifetime decreases. A clear sharp maximum was observed for the Mezcal sample with 55% ethanol. It should be noted that even though the amount of other components was diluted (by adding water or ethanol), increasing or decreasing the alcohol content led to a measurable variation of the lifetime. Therefore, based on these observations, we conclude that the traditional technique does work: the bubble lifetime shows a maximum value for a certain intermediate value of ethanol in Mezcal. Although the experiments were conducted in a controlled environment (standard laboratory conditions), a large variability of the measured lifetimes was found. Large scatter is often reported in experimental reports of the lifetime of bubbles, especially with water and other surface-contaminated liquids 13,18 . To our knowledge single bubble lifetime measurements for water-alcohol mixtures have not been reported to date, but similarly large scatter can be expected due to evaporation effects 19 . In the present case, the large uncertainty bars in our measurements may also result from the use of unfiltered artisanal Mezcal. According to the norm 7 , particulate matter is expected to be present. The presence of particles, and maybe fibers, would also explain the relatively large bursting thickness of the film.  Table 1. www.nature.com/scientificreports www.nature.com/scientificreports/ To evaluate the effect of surfactants and other components, a second set of experiments was conducted. In this set, different amounts of a 55% water-ethanol mixture were added to the sample M1, also shown in Fig. 3(a) and denoted as liquid MA2 in Table 1. In this manner, the amount of ethanol remains approximately constant but the quantity of surfactants is diluted with respect to the original sample. The result is clear: the pearl lifetime is significantly reduced compared to the original sample, even though the volumetric concentration of ethanol is nearly the same. We can, therefore, argue that the amount of surfactants is an important factor to determine the lifetime of pearls.
To begin elucidating the underlying physical mechanism of the process, experiments were performed considering mixtures of water and ethanol, varying the volume fraction of ethanol ranging from pure water to pure ethanol, as shown in Fig. 4. It was observed that even when no surfactants are present, the maximum pearl lifetime occurs at approximately 55% of ethanol fraction. Note, however, that the lifetimes for these mixtures are one order of magnitude smaller than those observed for Mezcal. Assuming that the dimensionless bubble lifetime (Eq. 2) is only a function of the Bond number (Eq. 1), we can argue that the lifetime of pearls is directly proportional to the fluid viscosity. While density and surface tension of water-ethanol mixtures have a monotonic behavior with ethanol fraction (from the pure water to the pure ethanol values), viscosity does not: its value increases from that of water, reaching a maximum at around 55% of ethanol, to then decrease to reach the value of pure ethanol 22 . The dotted line in Fig. 4 shows how viscosity varies (plotted on the right scale of the figure) with the volume of ethanol. Therefore, we can argue that both fluid properties and surfactant concentration affect the bubble lifetime. Consequently, it may be inferred that the maximum lifetime observed at the ethanol content at 55% is related to the fact that the liquid viscosity has a maximum value at that concentration.
More experiments were conducted with other Mezcal types, from different regions and agave species. For two cases, the alcohol content was also modified in the same manner as described before. The results are shown in Fig. 3(b), where the bubble lifetime for all liquids is shown as a function of alcohol content. The result is essentially the same, the bubble lifetime shows a clear maximum at around 55% of ethanol volume fraction.
numerical Results. From the results shown above, it is clear that the bubble lifetime is both dependent on the fluid properties and the amount of surfactants. However, to elucidate the details of the phenomena, we conducted a numerical study. Details of the numerical scheme can be found in the Methods Section.
First, the residence time of bubbles was computed for water-ethanol mixtures, with the aim to validate the numerical method and to gain some insight into the drainage process. Figure 4 shows the numerically obtained bubble lifetime, considering two bubble sizes which are close to the experimental values. The numerical lifetime was calculated by following the procedure described in the Methods Section. The simulations were obtained for a small amount of surfactant, as it is suspected to be the case for laboratory-clean aqueous mixtures. In very close agreement with the experiments, the bubble lifetime increases as the ethanol volume fraction increases, reaching a maximum for concentrations close to 55% of ethanol, to then decrease for higher concentrations. Given the variability of bubble size in the experiments, the quantitative agreement between experiments and numerical results is remarkable.
From the experiments, it is clear that the amount and type of surfactants affect the bubble lifetime. We conducted numerical simulations with increasing amounts of surfactants, as shown in Fig. 5(a). The bubble lifetime increases as the ethanol volume fraction increases, reaching a maximum for concentrations close to 55% of ethanol, to then decrease for higher concentrations, in close agreement with the experiments. Most importantly, the maximum lifetime increases with the amount of surfactants; the lifetime increases sharply with the concentration www.nature.com/scientificreports www.nature.com/scientificreports/ but appears to not increase significantly beyond a certain value of surfactants. These results are also in good qualitative agreement with the experiments. We have conducted many more simulations and analyses to further understand this effect. These results can be found in Atasi et al. 43 .
An interesting and challenging aspect of the Mezcal case is that, in addition to the amount of surfactants, the specific type of it is not known. Furthermore, there may be changes of the surfactant type depending on the species of Agave and production region for each Mezcal. To gain some general understanding in this respect, we also conducted numerical simulations varying the surface elasticity, E. This quantity measures the interface resistance to deformation 44 . Fig. 5(b) shows the bubble lifetime as a function of the surface elasticity parameter, considering a fixed value of ethanol content. Clearly, the lifetime increases significantly with the elasticity of the surface. It is important to note that this effect has not been discussed in-depth in the literature. Some initial tests and discussions can be found in Atasi et al. 43 .
The numerical simulations can also provide significant insight into the complex drainage process of bubbles in Mezcal. In addition to the effects on the lifetime, details of the flow within the film can also be analyzed. Figure 6 shows the fluid velocity (tangent to the interface) as a function of position within the film, at a certain angle from the apex (θ = 0.35 radians) at different times, with and without surfactants. An important difference between the two cases is that speed is significantly larger and uniform (across the thickness) for the case of a bubble without surfactants, in comparison to the case with surfactants. The reason for this reduction stems, first, from the immobilization of the surfaces due to surfactants. Moreover, when surfactants are present, for the angle shown, the fluid at the inner edge of the film (s = 0) can move in the opposite direction to that of the gravitational drainage, even for early times (t = 37 ms). For later times, t = 102 ms the fluid moves at a very small downward speed all across the gap. Note that the simulations consider the effect of surfactant transport but do not contemplate evaporation www.nature.com/scientificreports www.nature.com/scientificreports/ or temperature gradients. Despite this simplification, the lifetimes obtained numerically are in good agreement with those found experimentally. Therefore, we can argue that surfactant transport dominantly affects the drainage process.
Additionally, thanks to the unprecedented access to other properties of the flow from the numerical simulations, it is also possible to track the deformability of the interface. We have identified that the film becomes non-uniform along the bubble during the drainage. In particular, a neck region appears for a sufficiently large Bond number, associated to a local thinning of the film as shown in Fig. 7(b,c), which are used to validate the numerical results. Also, the position of the neck changes position from the top of the bubble toward the liquid bath, as the Bond number increases. The appearance of the neck region has been documented for the case of films 38 , but not as extensively for the case of surface bubbles 18 . This thinning zone leads to the rupture of the film which, as shown in Fig. 2(b), generally does not occur at the bubble apex. The results indicate that the position where the rupture starts is also a function of the Bond number.

Discussion
Let us now recast the bubble lifetimes in dimensionless terms. Considering the scaling from simple inspection (also arising from the extensional flow analysis), the dimensionless lifetime can be calculated according to Eq. (2). Figure 8(a) shows the bubble lifetime, ⁎ T life , as a function of the Bond number, Bo, defined in Eq. (1). Despite the large dispersion of the data, several key features can be readily identified. First, the dimensionless lifetime of the unaltered Mezcal is higher than any of the other cases, when either water or ethanol was added, or even for other Mezcal samples with other ethanol content. The dimensionless lifetime appears to increase for small values of the Bond number, to reach a maximum at around ≈ Bo 1, and then decrease as the value of Bo continues to increase. The Bond number allows us to evaluate the driving force for drainage, namely the capillary force for  Bo 1 and the gravity force for  Bo 1. The striking feature of the results in Fig. 8(a) is that Bo 1 indicates a transition, which coincides with the maximum dimensionless bubble lifetime. The Bond number compares the two main driving forces for drainage: gravitational drainage for large Bo and capillary-induced drainage for small Bo. Interestingly, the maximum dimensionless lifetime is reached when both effects are of the same order of magnitude.
To further support the existence of a critical Bo, we compared our results with predictions from the literature, shown in Fig. 8(a). For Bo < 1 the prediction by Howell, Eq. (17), shown by the solid line, shows an increasing  www.nature.com/scientificreports www.nature.com/scientificreports/ lifetime with Bond number. On the other hand, for Bo > 1, Eqs. (15) and (16), dashed and dotted line respectively, predict a decreasing trend of the bubble lifetime with Bo. These contrasting trend predictions corroborate the existence of a value of Bo at which ⁎ T life has a maximum value. The numerical results, in dimensionless form, are also shown in Fig. 8, along with the experiments. The numerical results show an increasing trend for small Bo, but at around ≈ Bo 1, a clear change of behavior is observed. Note that all simulations, conducted for two different bubble diameters and different ethanol concentrations, collapse into a single band when presented in dimensionless terms. The lifetime of the numerical results is shorter than that in the experiments resulting from the assumptions considered in the simulations. Nevertheless, the qualitative agreement is noteworthy. The cross-over for which gravity and capillary effects are of the same order of magnitude shows a non-monotonous behavior for which a maximum dimensionless lifetime has been observed, both numerically and experimentally. Strikingly, ≈ Bo 1 corresponds to a bubble size of about 2 mm, i.e. the size for which the lifetime is the longest in the experiments.
One important feature of the dimensionless bubble lifetime, ⁎ T life , is its large value, in the order of 10 4 -10 5 time-scale units. This is an indication that the characteristic time μ σ D/ is orders of magnitude smaller than the bubble lifetime. Hence, its magnitude cannot be used to characterize the phenomena properly. Therefore, we use ⁎⁎ T life instead, as defined by Eq. (11). Furthermore, if we consider a compensated dimensionless lifetime, given by the product ⁎⁎ T Bo life 3/2 , the data dispersion can be reduced by eliminating the dependence on the bubble diameter. www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 8(b) shows the compensated dimensionless lifetime as a function of the Bond number, Bo. The time scaling considers the rupture thickness of the film. As discussed above, the factors that determine the rupture thickness are not completely understood even for pure liquids 18 . In the present case, the presence of particles and surfactants, pose additional difficulties 45 . Hence, since we cannot reasonably assess how the rupture thickness changes for different liquids, we assume that h rup = 23.5 ± 4.8 μm is the same for all cases. This value corresponds to the measured film thickness for the case of unaltered Mezcal. The data presented in this form shows that the dimensional lifetime increases up to Bo of order unity, which is in good agreement with the scaling argument of Eq. (13). For bubbles for > Bo 1, the data clearly shows that the compensated dimensionless lifetime is relatively constant. In contrast, the prediction from Eq. (12, shows a slightly increasing trend. Note also that the magnitude of the dimensionless times is of O(10) for all experimental points, which indicates that the characteristic time t c (Eq. 10) is more representative of the flow phenomena in the draining film. Again, this argument supports the idea that the flow in the film is shear-dominated.
The numerical results are also shown in Fig. 8(b). As in the previous case, the lifetimes obtained from the simulations are shorter for all cases. However, and most importantly, the numerically obtained lifetimes collapse into a single band of data indicating that the scaling is correct. Furthermore, the trend found with the numerical results matches closely with the scaling arguments: the compensated dimensionless time increases for Bo < 1 and is roughly constant for Bo > 1. One of the advantages of the numerically obtained times is that the data scattering is greatly reduced, which allows us to observe the trend of results more clearly. From these results, it is further demonstrated that there is a critical Bo number, around unity, that characterizes a long bubble lifetime.
In summary, we conducted an experimental and numerical investigation to elucidate the phenomena behind the 'pearls of Mezcal' . By observing the extended lifetime of superficial bubbles, artisanal Mezcal makers have been able to identify the correct amount of ethanol in this distilled spirit without the need for any other measuring device. Intrigued by this empirical technique, we conducted experiments to test the technique under controlled conditions. Single bubbles were generated into a small container, and the time from reaching the surface until bursting spontaneously was measured. It was found that, indeed, bubbles in samples with ethanol volumetric content close to 55% had a characteristically long lifetime. We have found that, in dimensionless terms, the lifetime of a surface bubble increases with Bond number, reaching a maximum value at ≈ Bo 1, and then decreases for large values of Bo. We found that both an increase of the liquid viscosity and the presence of surfactants are needed to observe a long lifetime of bubbles. These factors explain why pearls in Mezcal have a particularly long life at a certain concentration of alcohol and a certain size. Clearly, the artisanal technique is the result of observation, empiricism and tradition. The insight gained in investigating this technique may help elucidate the fundamental mechanisms of bubble formation and stability in many other contexts. For example, the lifetime of surface bubbles could be used as a diagnostic tool to infer the presence of surfactants in a liquid: if the lifetime is larger than that expected for a pure/clean liquid, then the liquid is most likely contaminated. For instance, the bubble nests of some tropical fish have long lifetimes thanks to the presence of bio-surfactants 15 ; short-lived nest could indicate environmental contamination or health problems of the fish. Similarly, the contamination levels of water bodies could be correlated to long-lived foam and bubbles, as discussed by Zenit and Rodriguez-Rodriguez 46 . Furthermore, the lifetime of surface bubbles in applications such as froth floatation or anti-bubble formation could be adjusted by choosing the 'correct' surfactant to tune particle or gas capture rate 16,47 . Finally, recent reports indicate that the rupture process of surface bubbles plays a significant role in the spread of infectious diseases 48 ; hence, it is important to fully understand the nature of this process. The production of aerosols resulting from bursting bubbles could be reduced if the effects of surfactants are fully understood. The results presented in this investigation may shed some light into such processes of current importance.

Materials and methods informed consent.
The person who appears in Fig. 1(a) and in the S1 Video in the Supplementary Information section has given informed consent to be shown in the paper. He understands that the images may lead to identification.
Test fluids and physical properties. The physical properties of all the test liquids are reported in Table 1.
The origin of the Mezcal samples is also given. According to the Mexican norm, they all come from the Oaxaca state region. Seven Mezcal types were considered; three cases (M1, M6 and M7) were purposely altered to change their water or alcohol content. Also, pure water/ethanol mixtures were tested.
The surface tension was measured with a tensiometer (DynaTester, SITA). The density was measured with a 25 ml pycnometer. The alcohol content was inferred from the density measurement, considering room temperature, T room = 23 °C. The viscosity was not measured; considering the alcohol content, its value was assumed to correspond to that of a water-ethanol mixture and was obtained from tables 22 . Visualization of bubble formation process. The traditional technique to evaluate the alcohol content is shown in S1 Video (Supporting Information) and Fig. 1(a). The Maestro Mezcalero issues a stream of fluid from a long reed, with a round opening of about 2 mm in diameter. The jet impinges onto a small gourd cup, of about 10 cm in diameter. As the cup fills up, a pool of fluid is formed of several centimeters in depth. To reproduce the process in a controlled manner, the traditional reed was replaced by a glass pipette of 25 ml, with approximately the same exit diameter as the reed (2 mm). The sample of Mezcal (or other test liquids) of about 20 ml was placed in the pipette and was held fixed in a vertical position with a laboratory holding bracket. Below the exit of the pipette, a 10 × 10 cm 2 square transparent container was placed, filled with the same liquid to a depth of 2 cm. The process was filmed with a high-speed camera (FASTCAM-APX, Photron) using diffuse back-lighting at 1500 fps.
www.nature.com/scientificreports www.nature.com/scientificreports/ Measurement of lifetime. To accurately measure the lifetime of single bubbles, a simple experimental arrangement was used. A short cylindrical glass container of 1.6 cm in diameter and 1.2 cm in height was filled with the test liquid; the rim of the container was roughened to make it slightly hydrophobic. The amount of liquid was slightly larger than the volume of the container, such that a convex meniscus was formed. The arrangement was placed inside a closed container to minimize environmental disturbances during the experiment. At the side of the container, a needle was inserted through the wall via a plastic stopper plug. Bubbles were formed by slowly pushing air through the needle with a syringe. Since the free surface was slightly convex, when a bubble reached the surface it moved to the center of the container where it could be filmed. The bubble traveled approximately 1 cm through the liquid, lasting approximately 0.03 s from its formation until it reached the surface; the bubble moved slightly on the surface quickly reaching a stationary position. The bubble lifetime was measured from the moment when the bubble no longer moved and until it burst.
Experiments were performed with needles of different gauges ranging from 159 μm to 210 μm of internal diameter, to produce bubbles with slightly different diameters. The setup was illuminated from the bottom with a LED light. The bubble diameter was measured from the image obtained from the top; this measurement overestimates the bubble diameter of an equivalent spherical bubble by approximately 13% for the highest Bond numbers 49 . The process of rupture was measured with two synchronized high-speed cameras (Phantom SpeedSense), one aligned from the top and the other from the side. Different recording rates were considered: since the rupture time could take tens of seconds, an ordinary 30 fps was used to determine the bubble lifetime. To measure the rupture speed of the film, the recording rate was as high as 5,000 fps. All experiments were performed under standard laboratory conditions. The container was thoroughly rinsed with distilled water prior to each experiment.
Thickness of the film during bursting. The thickness of the film can be inferred from the speed at which the leading piercing rim moves as the bubble bursts. Considering the Taylor-Culick velocity 50,51 : where σ and ρ are the surface tension and liquid density, respectively. h rup is the thickness of the film. Experiments were conducted using fluid M1 and bubbles with the same diameter (D = 1.9 mm). numerical simulations. Direct numerical simulations were conducted by solving the Navier-Stokes equations coupled with the Level-Set method. We refer the reader to Atasi et al. 32 and Abadie et al. 52 for a detailed description of the method and its validation. Briefly, Navier-Stokes equations are solved for two Newtonian and incompressible fluids using the finite volume method (second-order accurate in time and space). Continuity is ensured through a projection method, and the capillary contribution is considered through the classical Continuum Surface Force method. The interface position is tracked using the Level-Set method where the transport of the signed distance to the interface is controlled through the re-distancing techniques. The novel aspect of the numerical scheme presented here is its ability to account for the surfactant concentration both in the liquid and on the interface as given by 53 : a I d where k a and k d are the adsorption and desorption kinetic constants, respectively, and C I is the surfactant concentration in the liquid in contact with the interface. It is assumed that the surface tension depends on the surfactant concentration on the interface according to an equation of state derived from the Langmuir adsorption isotherm: 0 where Γ σ = ∞ E RT / 0 is the surface elasticity parameter, R is the ideal gas constant, T is the absolute temperature, σ 0 is the surface tension of the clean interface and Γ ∞ is the maximum packing concentration of surfactants on the interface. The surface elasticity, which relates the surfactant concentration on an interface to its surface tension, has a great influence on the drainage of the liquid film atop of a bubble. In particular, for given flow www.nature.com/scientificreports www.nature.com/scientificreports/ conditions, the elasticity parameter dictates the strength of the Marangoni stress induced by the surfactant concentration gradient on the interface and the consequent immobilization of this interface. Typical values of the elasticity parameter for soluble surfactants in aqueous solutions are E = 0.05-0.25 55 .
The numerical solution of these equations is extensively described in Atasi et al. 32 . The proper implementation of each term, in particular, the surfactant transport on the gas-liquid interface and the computation of the resulting Marangoni stress, has been verified by adapted validation cases and a free rising bubble situation was compared with results from the literature.

Mesh.
Simulations were performed with a non-uniform axisymmetric orthogonal mesh characterized by 400 and 200 cells in the vertical and radial direction, respectively. The mesh size was refined in the vicinity of the interface to be able to properly capture both the film drainage and its rupture. A mesh convergence test was conducted by simulating the same case with increasing number of grid cells, as shown in Fig. 9(a). Clearly, the draining behavior is captured correctly for by the mesh size chosen (400 by 200) and higher.
We have observed that the numerical rupture of the film occurs when the film thickness is about where Δ is the grid size. From the experiments (see above), the film thickness at rupture is h rup ≈ 24 μm. Due to computing limitations, the minimum grid size we could apply in the numerical simulations was Δ = 10 μm, i.e.  . For the case of film drainage between two rigid interfaces, and in the limit  Bo 1, the thinning at the apex follows ∝ − h t 1/2 35 . Consequently, a decrease by a factor two for the critical film rupture should imply an increase by a factor four of the lifetime. Though this correction factor is certainly a good approximation, it was not applied since the drainage dynamics with surfactants can significantly differ from the one between rigid interfaces, especially at ≈ Bo 1 and for  Bo 1. To reconcile the calculated numerical bubble lifetimes with those measured experimentally, an extrapolation scheme was adopted, shown schematically in Fig. 9(b). Since the evolution of the film thickness showed a power-law dependence, the normalized film thickness was fitted to Validation. The numerical results were extensively validated in Atasi et al. 32 . One additional validation was conducted for this study, considering the shape of bubbles floating on the surface as the film drains. The shape that the bubble adopts while resting on the free surface changes with the value of the Bond number as shown in Fig. 7. The predictions are in close agreement with recent experiments 49 . More importantly, the simulations do capture the retarded drainage due to surface tension gradients induced by concentration gradients of surfactants along the interfaces, as discussed in the paper.