Symmetry Breaking in Chemical Systems: Engineering Complexity Through Self‐Organization and Marangoni Flows
Sangram Gore, Binaya R. Paudyal, Duarte Rocha, Mohamed Ali, Nader Masmoudi, Albert J. Bae, Christian Diddens, Detlef Lohse, Oliver Steinbock, Azam Gholami
TL;DR
This study shows how chemical waves can form flower-like patterns around obstacles due to fluid flows, enabling control over chemical reactions for microfluidic devices.
Contribution
The study demonstrates how solutal Marangoni flows can be used to engineer specific wave patterns in chemical reactions.
Findings
Marangoni flows destabilize circular chemical waves into flower-like patterns around hydrophilic obstacles.
Petal number increases linearly with obstacle diameter above a threshold.
Solutal Marangoni forces dominate thermal ones in influencing wave patterns.
Abstract
Far from equilibrium, chemical and biological systems can form complex patterns and waves through reaction‐diffusion coupling. Fluid motion often interferes with these self‐organized concentration patterns. This study examines the influence of Marangoni‐driven flows inside a thin layer of fluid ascending the outer surfaces of hydrophilic obstacles on the spatio‐temporal dynamics of chemical waves in the modified Belousov–Zhabotinsky reaction. These observations reveal that circular waves originate nearly simultaneously at the obstacles and propagate outward. In a covered setup, where evaporation is minimal, the wavefronts maintain their circular shape. However, in an uncovered setup with significant evaporation and resulting Marangoni flows, the interplay between surface tension‐driven Marangoni flows and gravity destabilizes the wavefronts, creating distinctive flower‐like patterns…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14- —Tamkeen under the NYU Abu Dhabi Research Institute (SITE research center)
- —Netherlands Organisation for Scientific Research (NWO)
- —Ministry of Economic Affairs (EZ); Canon Production Printing Netherlands B.V.; University of Twente; and Eindhoven University of Technology (project TKI HTSM ‐ CANON ‐ P1 ‐ PRINTHEAD & DROPLET FORMATI
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsNonlinear Dynamics and Pattern Formation · Fluid Dynamics and Thin Films · Micro and Nano Robotics
Introduction
1
Symmetry breaking is a fundamental concept in physics, biology and chemistry where small perturbations can lead to the emergence of complex structures and patterns from initially uniform states.^[^ 1, 2 ^]^ In fluid dynamics, this phenomenon is exemplified by the transition from simple, laminar flows to intricate turbulent structures, driven by instabilities that disrupt the symmetry of the system. In biology, Alan Turing's reaction‐diffusion theory proposes that the initial symmetry in embryos can be broken by the interaction of two diffusible molecules, whose dynamic interplay leads to the formation of complex patterns.^[^ 3 ^]^ Similarly, in chemical systems like the Belousov–Zhabotinsky (BZ) reaction, symmetry breaking can give rise to a rich spectrum of spatio‐temporal structures, such as rotating spirals and target waves as well as Turing patterns.^[^ 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 ^]^ These patterns, resulting from the interplay between nonlinear reaction kinetics and diffusion, offer profound insights into the mechanisms of self‐organization in both natural and synthetic systems.^[^ 1, 24, 25 ^]^
The interaction between chemical reactions and convective fluid motion, particularly through mechanisms such as Marangoni‐driven flows, profoundly influences the formation and breakdown of complex spatio‐temporal patterns.^[^ 26, 27, 28, 29, 30 ^]^ In pioneering studies on the BZ reaction,^[^ 31, 32, 33, 34, 35 ^]^ it was observed that in confined environments, chemical instabilities manifest themselves as propagating waves. In contrast, in systems with a free surface, the presence of Marangoni flows lead instead to the formation of cellular structures.^[^ 12, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46 ^]^ These flows are initiated by gradients in surface tension,^[^ 45, 47 ^]^ which might arise from either temperature variations or differences in solute concentrations across the fluid's surface, causing the fluid to flow from regions of lower to higher surface tension. For instance, in the ferroin‐catalyzed BZ reaction, the Fe^2 +^ complex has a higher surface activity than the Fe^3 +^ complex which is attributed to the differing hydrophilic properties of the catalyst's oxidized (more hydrophilic) and reduced forms. This difference leads to a decrease in surface tension of approximately 2 mN/m when the system transitions from the reduced to the oxidized state.^[^ 47 ^]^ In an uncovered setup, where evaporation is significant, the interplay of chemical waves, buoyancy‐driven flows due to density gradients, and thermal and solutal Marangoni flows due to surface tension gradients collectively influences the resulting chemo‐hydrodynamic patterns.^[^ 30, 32, 35, 41, 42, 48, 49, 50 ^]^
Evaporation‐induced Marangoni flows can be intensified by introducing hydrophilic vertical or inclined obstacles within the reactive medium. The thin film rising along the outer surfaces of these obstacles evaporates more quickly than the bulk fluid, creating thermal and/or solutal surface tension gradients by evaporative cooling and/or different volatilities of the multicomponent liquid. The delicate balance between upward‐driving surface tension forces and gravity‐induced downward motion leads to undulations at the fluid interface, potentially developing into finger‐like patterns if the fluid becomes sufficiently unstable.^[^ 51, 52, 53, 54, 55, 56, 57 ^]^ The obstacle's geometry also affects the dynamics at the droplet‐gas interface of thin films.^[^ 57, 58 ^]^ A cylindrical obstacle introduces unique conditions, such as azimuthal curvature and periodic boundary conditions, absent in flat obstacles.^[^ 57 ^]^ This periodicity ensures that an integer number of fingers can form along a cylinder, a condition that is not required for a flat obstacle, and the smallest number of fingers that can form along the contact line of a cylinder is one.^[^ 57 ^]^
Building on these groundworks, our study introduces a novel approach to explore the intricate interaction between evaporation‐driven Marangoni flows and chemical waves in the modified BZ reaction around hydrophilic obstacles. In a covered setup—where evaporation and therefore the evaporation‐induced surface tension differences are minimal—these obstacles act as focal points, emitting synchronous circular waves that retain their shape over multiple cycles. By contrast, in an uncovered setup with pronounced evaporation, thermal and solutal Marangoni flows trigger a hydrodynamic instability, yielding characteristic flower‐like wave patterns. The number of “petals” scales with obstacle diameter, and a critical diameter threshold is required for the instability to emerge. We also show that in our system, solutal Marangoni forces are more relevant than thermal ones. To elucidate the driving mechanism, we perform numerical simulations of a simplified binary fluid (e.g., water–ethanol) subject to solutal Marangoni instabilities, due to selective evaporation of one component (e.g., ethanol as compared to water). These simulations qualitatively reproduce the observed symmetry‐breaking patterns and lay the groundwork for future quantitative modeling.
Results
2
Experimental Results
2.1
Our quasi 2D setup consists of a periodic array of millimeter‐sized PDMS‐based cylindrical obstacles (see Figure 1A) that influence the spatio‐temporal dynamics of chemical waves in the modified BZ reaction (see Experimental Section). This version of the BZ reaction, known as CHD‐BZ, employs 1,4‐cyclohexanedione as a bubble‐free organic substrate and ferroin as a high‐absorbance redox indicator.^[^ 59, 60, 61, 62, 63, 64, 65, 66 ^]^ The addition of the ferroin redox indicator to the reaction mixture renders the oscillations visible, with periodic color changes between blue and red. The PDMS has undergone plasma treatment to make it hydrophilic, which allows the fluid to climb up around these obstacles, as illustrated in Figure 1B. We observe circular oxidation waves in blue color (high value of ferriin concentration Fe^3 +^) that initiate almost synchronously at the pillars and propagate inside the reduced state in red color (high concentration of ferroin Fe^2 +^) at a speed of approximately 7 mm/min and recurring every 45 sec (Figure 1C–L). When these waves collide, they mutually annihilate each other. The time evolution of these waves is strongly influenced by whether the setup is covered or uncovered (as illustrated in Figure 2B), which impacts the evaporation rate of the liquid (Figure 2C). Below, we discuss the results for both setups.
Wave dynamics in a covered versus an uncovered setup. A) A standard experimental setup featuring a periodic arrangement of PDMS‐based cylindrical obstacles designed to fit within a 10 cm Petri dish. B) A side view depicting the fluid profile around an exemplary hydrophilic pillar. C,D) Top views of concentric chemical waves in a covered setup, captured at two distinct times, show the waves initiating almost simultaneously at the pillars and spreading outward. The waves retain their circular shapes over time. E,F) Top views of chemical waves captured at two distinct times in an uncovered setup: initially, concentric waves form around the obstacles (E), and later, the wavefronts fracture to create flower‐like patterns (F). G) A space‐time plot along the orange line indicated in (C) demonstrates that the pillars serve as the wave centers. Black bars show the position of the obstacles. H) A space‐time plot traced along the cyan‐colored circle shown in (C) illustrates that the circular waves retain their shape through multiple wave cycles. I) The space‐time plot along the yellow line shown in (E). J) The space‐time diagram along the green circle indicated in (E) illustrates that after two cycles of circular waves, the wavefronts break. In this experiment, where pillars are 1 mm in diameter, 3 mm tall, and spaced 5 mm apart, most pillars develop three petals. K,L) Images displaying the entire Petri dish from two separate experiments: one covered and the other uncovered. The color‐coded wave period around each pillar is also shown in seconds. The timestamps reflect the elapsed time since the start of recording, which began shortly after the chemical solution was introduced into the setup. These images are typical of at least 50 experiments.
Number of petals depends on the pillar diameter and can vary over time. A) A PDMS mold with alternating pillar diameters of 1 mm and 1.5 mm is shown. B) A sketch illustrates the covered versus uncovered experimental setup. C) In the uncovered setup, the evaporation rate is approximately 2.5× 10−3 mg/(min mm2), which is about seven times higher than in the covered setup. D,E) Top views of chemical waves in an uncovered setup captured before and after the onset of instability in an experiment where pillar diameters alternate between rows, with 1.5 mm pillars in the top row and 1 mm pillars in the second row. Pillars with a diameter of 1.5 mm generally develop four petals, whereas those with a diameter of 1 mm typically produce three petals. F–H) Space‐time plots along the highlighted circles in panel (D) illustrate the progression of the instability. I) An example of a pillar that initially developed three petals, but over time, two of the petals merged.
Covered Setup: Circular Waves Initiate Synchronously at the Pillars and Maintain Their Circular Shape
2.1.1
In a covered setup, chemical waves begin nearly simultaneously on the pillars, with a wave period of about 45 sec, and maintain their circular shape across multiple wave cycles (Figure 1C,D; Video S1, Supporting Information). A space‐time plot along an arbitrary line crossing the obstacles illustrates that they serve as wave centers (Figure 1G). Additionally, to visualize the wave dynamics over several cycles, we compiled the light intensity data along a circle surrounding an arbitrary pillar, confirming that the wavefronts consistently retain their circular form (Figure 1H). Figure 1K and Video S2 (Supporting Information) present a separate experiment where the entire Petri dish is imaged, demonstrating similar wave dynamics. Fourier analysis of the waves emanating from the obstacles indicates that the majority have a periodicity of 45 sec (see Figure S1A– D).
Uncovered Setup: Circular Waves Begin Synchronously at the Pillars and Eventually Evolve into Flower‐Like Patterns
2.1.2
In an uncovered setup, where evaporation is more significant (Figure 2C), the obstacles continue to act as centers that emit synchronous concentric waves (Figure 1E,I; Video S3, Supporting Information). However, after approximately two cycles of these circular waves, an instability in the wavefront becomes apparent: the circular wavefronts fragment, leading to the formation of striking flower‐like patterns (see Figure 1F, J). Figure 1L and Video S4 (Supporting Information) present a separate experiment where the entire Petri dish is imaged, demonstrating similar wavefront instabilities. Consistent with the covered setup, Fourier analysis shows that the waves originating from the obstacles typically have a periodicity of about 45 sec (see Figure S1E– H).
Our experimental findings indicate that the number of petals is influenced by the diameter of the pillars. Typically, pillars with a diameter of 1 mm form three petals (Figure 1F), while those with a diameter of 1.5 mm often produce four petals (Figures 2 and 3), and pillars measuring 3.8 mm in diameter generate six to seven petals (Figure 4; Video S5, Supporting Information). Figure 2 presents an illustrative experiment where the pillar diameters alternate between 1 mm and 1.5 mm across neighboring rows, with the top row featuring pillars of 1.5 mm in diameter. In this setup, the 1 mm pillars typically form patterns with three petals, whereas the 1.5 mm pillars display flower patterns consisting of four petals (Figure 2E–I; Video S6, Supporting Information).
Wavefront instability at the boundary in an uncovered setup. A,B) Pillars with a diameter of 1.85 mm typically form four petals. The color‐coded phase map, generated using Hilbert transform, indicates that the pillars simultaneously act as wave centers. C) A space‐time plot along the cyan‐highlighted circle in panel (A) shows the development of the instability, with a noticeable evolution in the number of petals over time. The wave period is 50 sec. On the right, these snapshots at three different times (126 sec, 180 sec and 233 sec) are shown. D) A space‐time plot along the orange‐highlighted circle demonstrates the wavefront instability that also manifests at the PDMS wall on the periphery. This experiment is conducted in a small Petri dish.
Effect of the pillar diameter on the number of petals. A–D) As the pillar diameter decreases, so does the number of petals formed. For instance, pillars with a diameter of 3.8 mm typically develop around 7 petals, while those with a diameter of 2.8 mm produce about 6 petals. This number drops to 4 petals for pillars with a diameter of 1.8 mm, and further reduces to 3 petals for those with a diameter of 1 mm. E) Pillars with a diameter of 0.5 mm exhibit no wavefront instability, which verifies that there is a minimum pillar perimeter required for such instability to develop. F) The number N of petals increases linearly with the diameter d of the pillar. Below a critical diameter, wavefronts do not break and remain circular. The inset shows the average lateral extension πd/N of the petals.
Interestingly, for a small subset of the pillars, we observed that the number of petals can vary over time, as illustrated in Figure 2I. In a related observation, Figure 3 depicts an experiment where, after two cycles of circular wavefronts, the wavefronts around the pillars break. Initially, four petals form, but as time progresses, a fifth petal develops (Figure 3C). The synchronization between the waves initiated at different pillars is demonstrated through a phase map calculated using the Hilbert transform (see Experimental Section), as shown in Figure 3A,B. In this experiment, where we have captured images of the entire small Petri dish, the waves initiated at the PDMS boundary also exhibit wavefront instability, similar to those observed at the pillars (see Figure 3D; Video S7, Supporting Information).
To unravel the mechanisms (i) governing the formation of synchronous circular waves around the pillars and (ii) contributing to the instability of the wavefront in an uncovered setup, we began by examining the hypothesis that the permeability of PDMS to bromine gas (and potentially other chemicals) facilitates wave formation on the pillars, an idea supported by literature.^[^ 67, 68, 69 ^]^ To test this possibility, we conducted experiments using PDMS pillars coated with a thin layer of gold (approximately 200 nm thick, see Figure S2A, Supporting Information). Contrary to expectations, the pillars continued to function as centers for synchronous wave activity even with the gold coating, which blocks chemical absorption (Figure S2B and Video S8, Supporting Information). Additionally, we conducted experiments using obstacles made of acrylic (PMMA) with plasma treatment and observed that these pillars also functioned as centers for wave activity (Figure 5A–C). These results confirm that the absorption of chemicals into PDMS is not the primary mechanism driving the formation of concentric waves around the obstacles.
Experiments with PMMA‐made obstacles. A) Side view of the cylindrical obstacle array made from PMMA, with walls forming a square boundary. The fluid climbs up the plasma‐treated hydrophilic obstacles as well as the wall boundary. B) Images are subtracted every 3 seconds to better visualize the wavefronts initiating in the thin film region, wetting the pillars and the walls, and traveling downward. The obstacles act as wave centers, as shown in the space‐time plot in the red box, which stacks up light intensity along the red line. C) Initially, the wavefronts are planar along the wall and circular around the pillars, but they rapidly deform, developing into finger‐like patterns. This wavefront instability on the wall is clearly visible in the space‐time plot in the green box, created by stacking up the light intensity along the green line on the wall boundary. D) A top view of the waves around the obstacles, illustrating the breaking of the wavefront and the emergence of flower‐like patterns. These images are representative of at least 10 experiments.
Next, we explored the impact of the thin liquid layer ascending the PDMS‐made pillars, which are hydrophilic due to plasma treatment. Our investigation involved several experimental approaches. Initially, we conducted experiments with pillars that had not undergone plasma treatment and observed that most of the pillars did not act as centers for wave generation (see Figure 6A,B; Video S9, Supporting Information). Subsequently, we decreased the height of the pillars from 3 to 1 mm, ensuring that these shorter pillars were submerged to the same fluid level. We noted that the submerged pillars also did not act as wave centers (Video S10, Supporting Information).
Hydrophilic versus hydrophobic obstacles: An overview of a partitioned, uncovered Petri dish taken at two distinct times: before (A) and after (B) the emergence of instability. The upper portion of the dish has been plasma‐treated to render PDMS hydrophilic, whereas the lower section has been left untreated, maintaining its hydrophobic properties (see Video S9, Supporting Information). A zoomed‐in side view of the fluid near the pillars is also presented, emphasizing the rising fluid around the hydrophilic pillar. C–E) Side‐view imaging of the fluid climbing the hydrophilic pillars confirmed that the contact line remains pinned. As the liquid evaporates, the contact angle gradually decreases (see Video S12, Supporting Information). The cyan curves represent the air‐liquid interface at t=0; See the Supporting Information for meniscus analysis and curvature calculations.
Finally, to assess whether the formation of circular waves around the pillars was influenced by the proximity of adjacent pillars, we conducted experiments with a single isolated pillar. As shown in Figure 7 and Video S11 (Supporting Information), this solitary pillar also acted as a wave center and exhibited wavefront instability in an uncovered setup. Notably, under uncovered conditions, only waves originating from the pillar and the PDMS wall at the periphery exhibited wavefront instability, while off‐center waves retained their circular shape without fragmentation. These findings collectively highlight the crucial role of the climbing liquid layer on the pillars in both the generation of synchronous circular waves and the onset of wavefront instability.
A single pillar also acts as the wave center. A–C) An isolated pillar serves as the wave center and exhibits wavefront instability in an uncovered Petri dish. Additionally, we observe wavefront instability at the PDMS walls around the periphery. D) shows a flower‐like pattern with six or seven petals at different time points that formed around the solitary central pillar.
Experimental Evidence of a Fingering Instability
2.1.3
As previously discussed, our experiments reveal that the climbing film on the external surfaces of cylindrical pillars is essential for the formation of synchronous circular waves and the onset of wavefront instability. We observe that wavefront instability predominantly occurs in uncovered experimental setups, where evaporation is more significant compared to covered setups. Notably, in our control experiments, where the CHD‐BZ solution is covered with a thin layer of hexane to suppress evaporation, the obstacles continue to function as wave centers, but no wavefront instability is observed (Figure S2D,E, Supporting Information).
Moreover, side‐view imaging of the fluid climbing around the pillars confirms that the contact line remains pinned (Figure 6C–E; Video S12, Supporting Information), indicating the presence of upward capillary flows that counteract fluid loss caused by evaporation. Within this ascending liquid layer, the thinner downstream precursor film is expected to evaporate faster than the thicker upstream film, creating slight temperature and/or concentration gradients, due to evaporative cooling and/or selective evaporation of one or more components of the multicomponent liquid. These gradients generate surface tension differences, driving thermal and/or solutal Marangoni flows that ultimately disrupt the wavefront. Our hypothesis that waves originate from the thin ascending fluid layer is supported by our observations in Figure 5, using a PMMA mold with square wall boundaries. The side view of the waves confirms that they initiate at the thin film in the contact region, which wets the square wall boundaries, and then travel downward (Video 13, Supporting Information). Importantly, in an uncovered setup, the downward‐traveling wave along the wall boundary (as well as the waves initiated around the cylindrical obstacles) initially maintain a planar (circular) form but eventually break into fingers at later times (Figure 5B–D). This fingering instability is similarly visible in Figure 8E and Video S14 (Supporting Information) with PMMA‐based triangular obstacles within square wall boundaries. Lastly, we note that the initiation of waves at the ascending thin film is consistent with a separate experiment in which the Petri dish was slightly tilted.^[^ 38 ^]^ In this setup, we observed that synchronization waves initially appeared in regions with lower fluid thickness and subsequently propagated toward areas with higher fluid thickness (see Figure S5A,B and Video S15, Supporting Information).
Experiments with rounded rectangular obstacles and triangular obstacles. A) In an experiment examining the impact of evaporation, PDMS‐made obstacles with a rectangular body and rounded ends were used. B) The sequence of wave pattern of the highlighted pillar in panel (A) indicates that instability mainly arises at the curved ends. C) Oscillations of light intensity averaged over the red box shown in panel (B). D) A space‐time plot along the circular path marked in cyan in panel (B) illustrates the progression of wavefront instability. E) A side view of an experiment featuring triangular PMMA obstacles within a square boundary, where waves predominantly initiate at the sharp edges. Space‐time plots along the green and orange lines are also shown. The intensity averaged over the red box shows oscillations of 60 sec. Notably, the chemical fingers are clearly visible along the wall.
Literature reveals that Marangoni‐driven thin films on planar or curved substrates exhibit a complex array of dynamics,^[^ 54, 57, 70, 71 ^]^ dependent on the flat upstream film thickness trailing the advancing front. When the upstream film thickness is set by the competition effect of meniscus curvature and the surface tension gradient, the film is sufficiently thin for the advancing front to form a stationary compressive or Lax shock, which are well known to be unstable to contact line perturbations.^[^ 54, 55, 57, 70 ^]^ It is important to note that surface tension gradients, which drive Marangoni flows, can result from temperature gradients, concentration gradients, or a combination of both. Concentration gradients often develop due to the evaporation of more volatile chemical components within the ascending thin film, similar to the mechanism observed in the tears of wine phenomenon.^[^ 72, 73 ^]^ The interplay between Marangoni flows and gravity can disrupt the contact lines, evolving into a finger pattern. Specifically, the linear stability analysis in ref. [57] correlates linearly the number of fingers with the cylindrical pillar's perimeter, aligning with our findings (see Figure 4). This analysis also suggests a minimum pillar diameter threshold below which fingering instability does not occur. To validate this, we reduced the pillar diameter to 0.5 mm in our uncovered setup experiments, observing that the wavefronts remained nearly circular (see Figure 4E and Figure S3, Supporting Information). This resembles the situation of Rayleigh–Bénard convection where, for given driving strength, the linear instability cannot develop for containers with a too small aspect ratio (width/height).^[^ 74 ^]^
Moreover, we tested the impact of evaporation on fingering instability using obstacles with a rectangular geometry and circular ends (see Figure 8A–D; Video S16, Supporting Information). Given the increased surface area at the curved boundaries, we anticipated a higher evaporation rate here, which likely triggers the fingering instability predominantly at these locations, a hypothesis supported by our observations. A similar pattern was noted in our experiments with obstacles of triangular cross‐section, where waves primarily started at the sharp boundaries (Figure 8E, Video S14, Supporting Information). Additionally, we conducted experiments with both inverted tapered and tapered pillars, confirming that as long as a thin wetting layer surrounds the obstacles, they act as synchronous wave centers (Figure 9; Videos S17–S20, Supporting Information). In the experiments with tapered pillars, we observed that wavefront instability first appeared in non‐tapered pillars (our standard straight pillars with a 0° taper angle) and, after a time delay, in the tapered pillars as the angle gradually increased from 0° to approximately 26° (Figure 10). We attribute this effect to the higher evaporation rate in standard pillars compared to tapered pillars, due to their greater available surface area at the fluid‐pillar contact line. Interestingly, we found that the average number of petals remains unchanged across different taper angles. However, as the perimeter of the pillars at the fluid‐pillar contact line decreases with larger taper angles (see yellow circles in Figure 10B), the wavelength also decreases as the taper angle increases, ensuring that the same number of petals is accommodated.
Experiments with both inverted tapered and tapered pillars. A) Conical pillars with base diameters ranging from 0.5 to 2 mm, with the top row displaying standard cylindrical pillars of 3 mm diameter. B) Side view of centered waves observed before and after the onset of instability. The space‐time plot within the green box represents the light intensity along the green line over time, revealing that the waves originate from the climbing thin film and propagate downward. C,D) A close‐up side view of two neighboring cones with a base diameter of 1 mm. The cyan curve illustrates the air–liquid interface at t=0. E) Top view of the centered waves prior to the instability. F) Tapered pillars with a 3 mm base diameter that gradually decrease to smaller diameters toward the top. G–J) The side and top views of the centered waves are shown, along with the contact angle of the CHD‐BZ solution around the hydrophilic conical structures at t=0. Over time, as evaporation progresses, the contact angle decreases (t=300 sec). The cyan curves in panel (H) represent the liquid‐air interface at t=0; see Videos S17–S20 (Supporting Information).
A) Side‐view of tapered pillars used for data processing of waves centered around them in panels (B–D), with the last pillar at the bottom representing the standard non‐tapered pillar. B–D) The top view highlights the circles along which the space‐time plot in panel (C) is generated. The colored squares indicate the areas where light intensity is averaged and plotted over time in panel (D). Both the space‐time plot and light intensity oscillations show that after the first cycle of nearly synchronous circular waves, the second round of waves with a broken wavefront first appears in the non‐tapered pillar at the bottom, followed by a time delay in the tapered pillars. The two tapered pillars at the top, both tapered at approximately 26°, exhibit the longest delay. The orange lines serve as visual guides (see Video S18, Supporting Information).
Furthermore, we observed that during multiple cycles of covering and uncovering the setup, the wavefronts, which exhibit flower‐like patterns when uncovered, revert to their circular form once the setup is covered, and vice versa (Videos S21–S22, Supporting Information). To further investigate the role of boundaries, we conducted a control experiment using an uncovered, plasma‐treated glass Petri dish without any obstacles (Figure S4; Video S23, Supporting Information). Interestingly, we observed wavefront instability at the glass boundary, closely resembling the phenomena observed at the PDMS wall in Figure 3 and the PMMA wall in Figures 5 and 8E. A similar control experiment was performed in an uncovered Petri dish with a PDMS substrate (no pillars) and a PDMS wall boundary, where the wavefront instability was only observed at the PDMS wall (see Video S24, Supporting Information). Overall, these experiments confirm that the thin layer of climbing fluid at the boundary is critical for initiating the waves, and that evaporation significantly contributes to the wavefront instability.
Our final experiment that clearly proves the important role of evaporation in spontaneous breaking of azimuthal‐symmetry is an evaporating CHD‐BZ droplet. The preferential evaporation of more volatile chemical components at the rim can contribute to the segregation of other chemicals near the rim and eventually break the azimuthal symmetry. Figure 11A–D and Video S25 (Supporting Information) show spot formation near the boundary of a droplet. A similar effect has been reported in water/1,2‐hexanediol binary drops in Refs. [75, 76, 77, 78], where 1,2‐hexanediol micro‐droplets nucleate at the rim of the droplet. We note that in our evaporating droplet experiments, where micron‐sized beads were added as tracer particles to track the flow, we observed the coffee‐stain effect (see Video S26, Supporting Information; Figure 11E‐H.). This effect is characterized by the formation of a ring‐like structure of deposited beads after the droplet evaporates. The phenomenon arises from a combination of the non‐uniform evaporation rate along the droplet interface and the presence of a pinned contact line, which generates a continuous capillary flow from the center of the droplet toward the contact line to compensate for liquid loss during evaporation.^[^ 79, 80 ^]^ Furthermore, we identified a secondary instability at the droplet boundary (Figure 11H), which requires further investigation in future studies.
Symmetry breaking in an evaporating CHD‐BZ droplet and coffee‐stain effect. A–D) A droplet of volume 2.5 µl is deposited on a hydrophobic glass substrate. After about 500 sec, the azimuthal symmetry breaks and spots form at the rim evolving over time (T=22°C, RH=40%, θ0=12°). The images are captured in the reflection mode of a confocal microscope; see Video S25 (Supporting Information). E–H) In a separate experiment imaging the entire droplet, micron‐sized beads were introduced to track the flow. During the first 5 min of evaporation, the droplet boundary remains pinned to the substrate, leading to the accumulation of tracer particles at the rim due to capillary flows, clearly exhibiting the coffee‐stain effect, as depicted in panels (F–H). As evaporation continues and the reaction progresses, a secondary instability develops at the rim of the shrinking droplet, as shown in panel (H) (see Video S26, Supporting Information).
Estimation of Dimensionless Numbers and Timescale Comparison
2.2
In our experiments, the chemical fingers forming on the pillars or on the walls of the PDMS/PMMA setups exhibit a wavelength of approximately λ ≈ 2 mm, see Figure 4F, inset. We hypothesize that the fluid fingers have a similar wavelength. Assuming the climbing fluid depth is in the range h = 20 − 100 µm,^[^ 73, 81 ^]^ we calculate the Marangoni stress τ from the balance with gravity, τ = 2ρgh/3,^[^ 70 ^]^ which represents the balance between the competing Marangoni stress and gravitational forces. This calculation yields a surface stress range of τ = 0.13 − 0.67 N/m ^2^. Evaluating the corresponding solutal Marangoni number for τ = 0.67 N/m ^2^ and h = 20 µm yields:
where μ is the viscosity of the CHD‐BZ solution and D is molecular diffusivity of chemical components. The corresponding thermal Marangoni number, defined as:
is two orders of magnitude smaller due to the ratio of molecular diffusivity to thermal diffusivity, given by the inverse Lewis number
Indeed, in general, for evaporating multicomponent liquids with different volatilities of their components, solutal Marangoni effects are more relevant than thermal ones.^[^ 82, 83 ^]^ We now proceed by considering a dimensional analysis similar to ref. [73] that examines the interaction between a deformable free surface and the steady‐state concentration profile, which is governed by diffusion, evaporation, and convective transport along the pillars. Evaporation reduces the concentration of the more volatile chemical component with lower surface tension near the surface; however, under the assumption of a quasi‐steady state, fresh fluid from the reservoir continuously replenishes the depleted region, establishing a stable vertical concentration gradient. We note that given the complexity of the CHD‐BZ reaction, where multiple chemical components are continuously produced and consumed, the validity of this assumption requires further extensive investigation in future studies. If the free surface experiences an upward perturbation, the concentration profile must adjust before the interface returns to its original shape. If the concentration realigns more rapidly than the interface relaxes, the perturbed region develops a lower concentration than its surroundings, generating a surface tension gradient that draws fluid into the bump and amplifies the perturbation. Similarly, depressions in the free surface create local surface tension minima, pulling fluid away and reinforcing the deformation. This feedback loop sustains convection and enhances surface fluctuations, making the film susceptible to instability driven by Marangoni forces. The system's stability depends on the relative timescales of diffusion and interface relaxation. If the diffusive timescale is significantly shorter than the relaxation time, the system is unstable, i.e., T diff/T relax ≪ 1. The relaxation time T relax is estimated by balancing pressure gradients with viscous stresses,
where V _ f _ ∼ 0.5 mm/s is the estimated speed of fluid movement up the pillars. The restoring pressure gradient is determined by either surface tension or gravity:
For h = 20 µm and surface tension γ = 65 mN/m,^[^ 84 ^]^ the corresponding timescale comparison yields the dimensionless ratio for surface tension restoring pressure as:
while for gravitational restoring pressure, it is
Since both values are significantly less than one, we conclude that the system is susceptible to Marangoni‐driven convection, where surface tension gradients enhance surface deformations and sustain instability. We note that accurately estimating the characteristic timescales and length scales in our study requires simultaneous experimental measurements of film thickness and concentration gradients driving Marangoni flows. Such simultaneous measurements are experimentally challenging, and efforts to address this are currently ongoing in our laboratory.
Numerical Simulations
2.3
Comprehensive numerical modeling of the CHD–BZ reaction—including its full chemical kinetics—is beyond the scope of this study and will be pursued in future work. To capture the Marangoni‐driven hydrodynamic instability qualitatively, we instead employ a simplified two‐component fluid mixture exhibiting Marangoni instabilities – e.g., a binary ethanol–water mixture, in which ethanol is both more volatile and has lower surface tension than water. We then reproduce the experimental geometry by placing this mixture around a pillar (Figure 12) and considering diffusion‐limited evaporation of both components. The contact line remains pinned to the pillar walls, and a symmetry boundary condition is imposed on the right. Our evaporation model follows the approach of Section 2.1 in ref. [83], suitably adapted for this configuration (see Experimental Section 4).
A) Geometry of the problem setup based on the experiments at time t = 599 sec (see Figure 6E). The liquid layer (magenta) is pinned to the pillar walls, forming an angle θ = 52° with the pillar walls and it is at a height of 2.5 mm at the walls. The liquid evaporates through diffusion of vapor in the gas phase, as schematically represented by the green arrows. Note the longer arrows close to the contact line, representing the larger evaporation flux close to the contact line singularity. B) Schematics of the Marangoni instability phenomena in a box with periodic side boundary and an evaporating top surface, where only one of the two components evaporates. When the evaporating component has a lower surface tension, the interface increases its concentration on higher surface tension fluid. May a perturbation locally reduce the surface tension somewhere in the interface, it will onset a Marangoni flow, which in turn draws fluid away from the perturbed region (yellow arrows in second panel). Due to continuity, the displaced fluid must be replenished, which brings fluid from the bulk toward the perturbed region (green arrow in second panel). This draws lower surface tension fluid to the perturbed region, further decreasing the surface tension and creating a positive feedback loop that destabilizes the system – leading to the so‐called Marangoni instability.
First, we compute the axisymmetric quasi‐stationary base state as a function of the Marangoni number, Ma. Next, we carry out a linear stability analysis with respect to azimuthal perturbations of wavenumber m. Our results reveal that the radial width of the liquid layer strongly influences stability: wider layers support longer‐wavelength modes. As illustrated schematically in Figure 12B (see also Figure 1 of ref. [83]), a local drop in surface tension drives Marangoni flow away from the perturbed region. The resulting pressure gradient pulls fluid – enriched in the low‐surface tension component – from the bulk back into that region, creating a positive feedback loop that can trigger the instability. The resulting positive feedback can destabilize the system, leading to a Marangoni instability. A wider liquid layer reduces flow restrictions, thereby allowing longer‐wavelength perturbations to develop.
To specifically assess the role of the Marangoni number in this confined geometry and avoid varying geometry effects, the liquid width is fixed at 3.5 mm. We then vary Ma and analyze the stability of the base compositional state with respect to perturbations of different azimuthal wavenumbers. Figure 13 shows the real eigenvalues of the system as functions of both m and Ma. The results reveal that the system remains stable for low Ma values (Ma⪅100) but becomes unstable as Ma increases. At low Ma, the dominant instabilities arise from perturbations with lower m, while higher Ma shifts the instability toward larger m.
Dispersion relation of the system for different Ma values. The real part of the nondimensionalized eigenvalues is shown as a function of the azimuthal wavenumber m for Ma values of 100, 2500, 5000, 7500, and 10000. The Re(λ) = 0‐line indicates the stability threshold, where the real part of the eigenvalue becomes positive, indicating instability.
To visualize a fully 3D compositional field, we can plot
where ϵ is a small perturbation amplitude, y ^0^ is the axisymmetric base state, and y ^ m ^ is the eigenfunction for azimuthal mode m, see Figure 14 showing the full flow (panel A) and the mode m = 5 (panels B–D). This representation assumes a single dominant eigenmode—which may not strictly hold here—but serves to qualitatively illustrate the emerging instability pattern.
*A) Nondimensional axisymmetric base solution at Ma = 2500 for liquid's velocity magnitude (left), liquid's composition (right) and gas vapor's concentration. In the liquid‐gas interface, the arrows shown represent the evaporation rate. B) Eigenfunction with azimuthal mode m = 5 for the liquid's (r, z)‐velocity magnitude (left) and liquid's composition (right). C) Top view of the 3D expansion y(r, z, ψ) = y 0(r, z) + ϵy
m (r, z)e
imψ + c.c. for the composition. D) Isometric view of the same expansion.*
Conclusion and Outlook
3
In this study, we introduce a periodic array of hydrophilic cylindrical obstacles into the CHD‐BZ reaction, which rapidly become synchronized centers of wave activity. The thin layer of chemical solution ascending the obstacle surfaces is essential for enabling these structures to function as wave sources, thus precisely determining their spatial positions (Figure 1C). Although an isolated obstacle can also act as a wave center (Figure 7), an array arrangement effectively prevents the formation of random wave centers elsewhere in the solution. Under covered conditions, where evaporation is minimized, wavefronts maintain stable circular shapes through multiple cycles. In contrast, uncovered conditions—with significant evaporation—induce Marangoni‐driven instabilities, causing wavefront breakage and the emergence of striking flower‐like patterns (Figure 1F). The number of petals in these patterns scales linearly with obstacle diameter (Figure 4), and there is a minimum diameter required for this instability to manifest (Figure S3, Supporting Information). Additionally, control experiments where the weight of PDMS molds was monitored after each five‐minute wave recording show negligible weight changes, conclusively ruling out any influence of PDMS swelling on our experimental results (Figure S7, Supporting Information).
We performed numerical simulations of a simplified binary fluid mixture exhibiting Marangoni instabilities (for example, an ethanol–water system) to investigate the evaporation‐driven hydrodynamic instability. Convective flows in this system arise from selective evaporation of one component (ethanol) and the resulting surface‐tension gradients at the air–liquid interface. In these qualitative simulations—where we assume that the contact line of the fluid mixture is pinned around a single pillar—we show that, in principle, such a simplified ethanol–water mixture can also exhibit a “compositional flower” pattern surrounding the pillar. We stress that these simulations are intended only for qualitative insight and do not quantitatively reproduce the experiments, which involve many additional complexities. Moreover, they do not capture thin‐film dynamics or fingering instabilities that could emerge once the contact line moves; those phenomena will be addressed in future work.
The behavior of a film climbing the exterior of a vertical or an inclined wall, influenced by the opposing forces of a thermally/chemically driven surface tension gradient and gravity, is extensively explored in the literature through numerical analysis using a thin‐film model (see review articles Refs. [51, 52, 53] and references therein). In the case of a cylindrical obstacle, these models depend on three parameters: the cylinder radius, the upstream film thickness, and the downstream precursor film thickness. As the cylinder radius approaches infinity, the model converges to the scenario of a Marangoni‐driven film ascending a vertical plate.^[^ 54, 55, 56, 73 ^]^ ref. [57] investigates the stability of the contact line by introducing an infinitesimal 2D perturbation to the traveling wave solution, characterized by sinusoidal variations in the azimuthal direction. Through linear stability analysis, the study determines that the traveling wave (the so called Lax shock) is unstable to long‐wavelength perturbations at the contact line. Given the periodic nature of the contact line in the azimuthal direction, only an integer number of fingers can form along it. The analysis identifies the most unstable mode and a neutrally stable cutoff mode, beyond which all modes are stable. These modes increase linearly with the cylinder radius, indicating that the number of fingers formed along the contact line correlates with the cylinder's circumference. There exists a critical cylinder radius below which the contact line remains stable, with no fingers forming. Our experimental results provide direct validation of these theoretical predictions. In future work, we will extend the thin‐film framework of Smolka et al.^[^ 57 ^]^ to rounded‐rectangular obstacles—introducing asymmetric evaporation rates and surface‐tension gradients—and couple those thin‐film equations with the reaction–diffusion equations,^[^ 6 ^]^ to capture the fully coupled fluid–chemical dynamics in the thin film around obstacles of various geometries.
In our experiments, by tilting the Petri dish, it is possible to selectively control which array of obstacles initially becomes the focal point for wave generation (Figure S5A,B). As reported in ref. [32] and consistent with our experiments, in a tilted Petri dish, waves originate in the thinner areas of the fluid, where convective flows are negligible, and progress toward the thicker regions, where convective flows can develop due to evaporation and the resulting Marangoni forces. The authors of that paper assumed that these convective flows generate a heterogeneous medium where the refractory tail of a preceding wave annihilates part of the following wavefront, resulting in segmented waves that curl inward to form spirals. Over time, a chaotic wave pattern develops. In another study to understand the role of the fluid thickness in an uncovered Petri dish, Rossi et al.^[^ 41 ^]^ have demonstrated that for fluid layers thinner than 1.0 mm, the BZ system behaves predominantly as a pure reaction‐diffusion system, and convective effects are negligible. In layers with an intermediate thickness of 1.0–3.0 mm (the fluid thickness in our experiments was about 1.5 mm), Marangoni instabilities dominate, allowing rippled waves to propagate without breaking, and segmented patterns are not observed. However, for layers thicker than 3.0 mm, buoyancy‐driven flows become significant, and wave segmentation processes occur. This resembles the competition between buoyancy‐driven flow and solutal Marangoni in evaporating binary droplets.^[^ 80 ^]^ Our experiment using an uncovered, plasma‐treated glass Petri dish with a fluid thickness of 2 mm also revealed rippled wavefronts and Marangoni‐driven mosaic patterns (Figure S7A and Video S23, Supporting Information).^[^ 43, 44, 45, 46 ^]^ We notice that at the boundary of the Petri dish, we observed a Marangoni‐driven fingering instability, similar to the phenomena we see near the hydrophilic obstacles.
Contact angle measurements on both tapered and non‐tapered pillars confirm that the contact line remains pinned (Figure 6C–E and 9C,H), indicating the presence of compensatory upward capillary flows. Additionally, experiments with tapered and inverse‐tapered pillars demonstrate that as long as a wetting fluid layer surrounds the obstacle, it continues to function as a wave center. In tapered pillars, we observed that wavefront instability initially appears in straight pillars and, after a time delay, in pillars tapered at approximately 26°. This delay is attributed to the higher evaporation rate in straight pillars due to their greater exposed surface area at the contact line. Interestingly, in a separate experiment investigating the effect of surface roughness, we 3D‐printed conical pillars with a wavy surface (Figure S8, Supporting Information). The results were largely consistent with our standard experiments without surface roughness, except for increased variability in petal size.
Finally, in our experiments with an evaporating CHD‐BZ droplet, we observed the formation of spots at the rim. To elucidate the underlying mechanism, identify the more volatile components, and determine which chemical elements are segregated at the boundary, additional comprehensive experiments are necessary. Tracer beads added to the CHD‐BZ solution in our evaporating droplet experiment reveal a pinned contact line and the coffee‐stain effect, further suggesting the presence of compensatory capillary flows. Interestingly, after approximately one hour, we observed a complex secondary instability at the droplet boundary, which warrants further investigation. Ongoing studies in our laboratory are focused on exploring this phenomenon in more detail.
Our study on the spatio‐temporal dynamics of chemical waves interacting with a periodic array of obstacles highlights several potential applications. By demonstrating how these obstacles can synchronize and guide wave propagation, we introduce a strategy for engineering specific wave patterns that could be especially beneficial in advanced microfluidic devices. We further propose that these guided patterns can be harnessed not only for exploring fundamental reaction–diffusion dynamics, but also for practical applications such as chemical computing, where controlled wave interactions serve as the basis for logic operations and information processing.^[^ 85 ^]^ Moreover, our findings on manipulating Marangoni‐driven flows around hydrophilic obstacles to shape wave patterns can be extended to colloidal droplet evaporation. By employing similar obstacle‐based or surface‐modification techniques, one could direct evaporative flows and steer colloidal self‐assembly, opening new possibilities for creating uniform, tunable deposition patterns in printing, photonics, and related applications.^[^ 86 ^]^ Finally, by adjusting parameters like obstacle geometry or evaporation conditions, one could similarly guide chemical wave propagation within gels, thereby influencing their mechanical oscillations. This capability paves the way for designing advanced biomimetic actuators or “smart” materials where desired wave‐driven shape changes (or gel motions) can be precisely programmed or directed.^[^ 87 ^]^ Overall, our results provide a foundation for further manipulation of chemical wave behaviors, offering exciting opportunities for technology and manufacturing.^[^ 21, 88 ^]^
Experimental Section
4
Chemical Preparation
The following chemicals and materials were used in this study. Sodium Bromate (NaBrO_3_), Sodium Bromide (NaBr), Conc. Sulfuric acid (Conc. H_2_SO_4_), 1,4‐Cyclohexanedione, Ferrous sulfate, 1,10‐phenanthroline was provided by sigma and used without any further purification. Ferroin was synthesized by mixing Ferrous sulfate hexahydrate and 1,10‐phenanthroline in 1:3 molar ratio. For CHD‐BZ experiments different concentration of the stock solutions were prepared in dH_2_O. Solution A: (494.58 mM solution of sodium bromate) Solution B: (891.9 mM solution of 1,4‐Cyclohexanedione) Solution C: (971.9 mM solution of NaBr). A mixture was prepared in a beaker by combining 9 mL of solution A with 1.5 mL of solution B, which was then homogenized using a magnetic stirrer. Following this, 1.5 mL of the Ferroin indicator was added to the mix. Under a chemical fume hood, 0.8 mL of solution C was cautiously introduced to the reaction. The stirring continued for an estimated 1 h and 45 min. Subsequent to this phase, an additional 1.5 mL of Ferroin was incorporated into the mixture. The stirring persisted at a controlled temperature of 23°C for a further 30 minutes, after which it was halted. Attention was then given to the mixture for any signs of color transformation, signaling the start of wave propagation, typically observed 3.5–4 h post the initial mixing process. Following the onset of color change, the CHD‐BZ reaction mixture was transferred to a plasma‐treated PDMS mold featuring an array of macro‐pillars, as specified in Figure 1A. These pillars, typically 1 mm in diameter and 3 mm high, are arrayed in a square grid with 5 mm spacing. In experiments with variable pillar dimensions and spacings, the volume of the CHD‐BZ solution was adjusted to ensure a liquid height of about 1.5 mm, thus avoiding submergence of the pillars. Upon pouring, the solution was swirled gently until the mixture achieved a brown color. The Petri dish was then positioned undisturbed beneath a camera to record the wave propagation, under the illumination of white LED light. For each experiment, we record the waves for 5–10 min, observing until the waves from the Petri dish boundary interact with those centered around the obstacles. Afterward, we stop the recording, swirl the solution in the Petri dish to reset the system, and start recording again. Typically, the reagents are fully depleted after about 2 h. In some experiments, a smaller, modified 50 mm Petri dish was used, and the volume of the CHD‐BZ solution was correspondingly reduced to 1.5 mL to maintain the same liquid height. It's important to note that in our standard setups, unless specified otherwise, the pillars do not submerge in the chemical solution. The hydrophilic nature of plasma‐treated PDMS causes a thin layer of liquid to ascend the pillars (Figure 1B), altering the liquid surface near the obstacles as depicted in Figure S2C (Supporting Information). For some control experiments, the surface of the PDMS was coated with a 20 nm titanium layer and then with a pure gold layer. This was done using a Physical Vapor Deposition (PVD) technique with the Lesker Proline PVD 75 apparatus, rotating the substrate four times to achieve uniform coverage on the PDMS pillars' surfaces.
Image Processing
The wave period was roughly 45 sec, and images of the chemical waves were captured at a rate of 1 frame per second across multiple wave cycles. To minimize spatial noise, we average the captured images over several wave cycles, typically 6–8, and then subtract each individual image from the time‐averaged image. After this, a Gaussian blur filter was applied and a customized MATLAB code was used to perform the Hilbert transform, which helps to calculate the phase map. This code is provided in the Supporting Information.
Numerical Simulations
Our numerical simulations qualitatively capture the hydrodynamic instabilities observed in CHD‐BZ experiments. Our aim is to demonstrate that solutal Marangoni–driven instabilities, induced by evaporation, are sufficient to explain the symmetry‐breaking patterns observed experimentally. We do not attempt a quantitative fit to the data, as that would require a substantially more complex model.
We model the system using the experimental geometry: a binary mixture exhibiting Marangoni instabilities—one example being an ethanol–water system—is confined around a cylindrical pillar of radius R = 1 mm, height of 5 mm, and center‐to‐center spacing of 7 mm (see Figure 12). At the evaluation time t = 599 sec—recorded from our experiments—the liquid layer meets the pillar walls at a contact angle θ = 52° and has a height of 2.5 mm there. The evaporation is assumed to be limited by vapor diffusion, a reasonable assumption for liquids sufficiently below the boiling point. The contact line remains pinned to the pillar walls, and a symmetry boundary is imposed at the midpoint between adjacent pillars. Evaporation is modeled using the quasi‐stationary framework of Rocha et al.,^[^ 83 ^]^ suitably adapted to this geometry.
We track the mass fraction y _ A _ of the more volatile component A (e.g., ethanol), with y _ B _ = 1 − y _ A _ denoting the mass fraction of the less volatile component B (e.g., water). For slowly evaporating liquids, compositional variations remain small, so we decompose
where y _ A, 0_(t 0) is the spatially averaged mass fraction at the evaluation time t 0 and y(x, t) captures local deviations.
Although y _ A, 0_(t) evolves over time, the flow can be treated as quasi‐steady at each drying stage, provided the evaporation rate remains low. All material properties are held constant except for surface tension, which is a function of the mass fraction variations. We assume a linear dependence
The system is nondimensionalized as follows:
where D 0 and μ_0_ are the diffusion coefficient and viscosity at the mean composition,
is the equilibrium vapor concentration (with cαpure the pure‐component saturation, γ_α_ the activity coefficient, and x α the mole fraction), and cα∞ is the far‐field vapor concentration.
Assuming small local variations in y, we treat cαeq as constant.^[^ 80 ^]^ Under this scaling, the dimensionless vapor field c∼ in the gas phase satisfies the Laplace equation
subject to the boundary conditions c∼=0 at the liquid–vapor interface and c∼=1 in the far field.^[^ 80, 89, 90 ^]^ The solution c∼ then maps directly back to the physical vapor concentration c α. The dimensionless evaporation rate j∼ is given by the diffusive mass flux
with n being the outward normal vector at the liquid–vapour interface. Evaporation drives variations in the liquid composition through the flux boundary condition:
where Ev is the evaporation number, which was originally defined by [80] as
where DAvap and DBvap are the vapour diffusion coefficients of components A and B, respectively, and ρ_0_ is the liquid density evaluated at the spatially averaged composition.
In the liquid phase, we consider the incompressible Stokes equations (neglecting inertia and gravity):
In order to reduce the number of parameters in the problem, we define ξ = y/Ev, such that Equation (15) becomes:
The evolution of ξ in the liquid phase is then given by the advection‐diffusion equation:
Marangoni stresses are considered at the liquid‐vapour interface via the boundary condition:
with t being the tangential vector at the liquid–vapour interface and Ma being the solutal Marangoni number, defined as
It was considered that the shape of the liquid layer does not change due to surface tension effects (i.e., capillary number Ca → 0). At the liquid–gas interface, it was considered that evaporation is slow, i.e., u _ I _ ≪ u, where u _ I _ is the velocity at the interface and u is the velocity in the bulk liquid. This allows to assume that the interface is stationary, i.e., the boundary condition u∼·n=0 was imposed at the liquid–vapour interface, which was enforced via a locally adjusted normal traction. The volume of liquid in the system is fixed to a value such that the contact angle is θ = 52° at the pillar walls and the liquid is pinned to the pillar walls. At the substrate, we impose a no‐slip boundary condition u∼=0 and a no‐penetration boundary condition ∇ξ·n=0. The spatial average of ξ was constraint to zero, in accordance with the decomposition (9).
The system of equations is solved using the finite element method with the open‐source library pyoomph,^[^ 91 ^]^ which is based on oomph‐lib ^[^ 92 ^]^ and GiNaC.^[^ 93 ^]^ All domains are discretised using triangular elements. The operators in the system of equations are given in a cylindrical coordinate system. Linear basis functions are used for ξ and p∼, whilst quadratic basis functions are used for the velocity field u∼. We follow the methodology outlined in ref. [83] to compute the steady base state of the system and its stability against perturbations characterised by different azimuthal wavenumbers m. This allows to assess whether the system's symmetry is broken and forms patterns with a specific azimuthal wavenumber, like those observed in the experiments.
The stability of the simplified system depends on the Marangoni number Ma, on the distance between pillars D, on the angle θ at the pillar walls, and on the liquid height at the pillar walls. In the following, for simplicity, the liquid height was fixed at 2.5 mm, D = 7 mm, θ = 52°, and vary Ma.
Conflict of Interest
The authors declare no conflict of interest.
Author Contributions
S.G. and B.P. conducted the experiments, D.R., C.D. and D.L. designed and carried out the numerical simulations. A.G. conceived the idea, designed the experiments, analyzed the experimental data, and wrote the first draft of the manuscript. All authors discussed the results and revised the manuscript.
Supporting information
Supporting Information
Supplemental Video 1
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1R. B. Hoyle , Pattern Formation: An Introduction to Methods , Cambridge University Press, 2006.
- 2H. Meinhardt , Rep. Prog. Phys. 1992, 55, 797.
- 3A. M. Turing , Bull. Math. Biol. 1990, 52, 153.2185858 10.1007/BF 02459572 · doi ↗ · pubmed ↗
- 4A. Zaikin , A. Zhabotinsky , Nature 1970, 225, 535.10.1038/225535 b 016056595 · doi ↗ · pubmed ↗
- 5I. R. Epstein , J. A. Pojman , An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos , Oxford University Press, 1998.
- 6R. J. Field , R. M. Noyes , J. Chem. Phys. 1974, 60, 1877.
- 7J. J. Tyson , The Belousov‐Zhabotinskii Reaction , vol. 10, Springer Science & Business Media, 2013.
- 8O. Steinbock , V. Zykov , S. C. Müller , Nature 1993, 366, 322.