From nonlinear reaction-diffusion processes to permanent microscale structures

Dynamic systems far from the thermodynamic equilibrium are a treasure trove of complex spatiotemporal patterns.^1,2 Living systems have evolved to utilize these patterns to perform tasks requiring the relay of signals over macroscopic distances as exemplified by nonlinear waves of action potentials in neuronal and cardiac tissue. In chemistry, related dissipative patterns include traveling excitation waves and Turing patterns.^1,3,4 These have typically been observed in open systems such as chemical reactors that continuously provide fresh reactants and remove products, thus establishing non-transient, far-from-equilibrium conditions. Once the externally controlled mass exchange is stopped, both chemical and biological patterns cease to exist, and the lack of driving forces causes the system to approach their unremarkable equilibrium state.

Profound opportunities for the synthesis of novel materials, various engineering applications, and the discovery of physical phenomena could result from the ability to transform the complex dissipative patterns into permanent structures that persist after the end of the chemical self-organization. To date, only a very few experimental systems are accomplishing this “fossilization” of spatiotemporal complexity. The most prominent examples in chemistry are Liesegang patterns.^5,6 These precipitation and crystallization structures are typically macroscopic bands or helical structures and form in gel systems from inorganic ions such as silver cations and chromate anions that slowly diffuse into each other. Similar precipitation reactions—i.e., reactions that form a solid, insoluble product—create the so-called chemical gardens, which consist of hollow tubes that grow out of metal salt particles placed into a sodium silicate solution.⁷ 

Another chemical system that has gained considerable interest is the reaction of aqueous solutions (optimal pH 10–12) containing alkaline earth metals (Ba$2+$⁠, Sr$2+$⁠, and Ca$2+$⁠) and silicate ions with atmospheric CO$2$ or added carbonate.^8–13 This simple reaction medium produces smoothly curved, polycrystalline aggregates known as “biomorphs.” These entirely inorganic microstructures take the shape of leaflike sheets, helices, funnels, or patterns reminiscent of corals.¹⁴ The typical size of these solid reaction products ranges from 10 $μ$m to 1 mm. In addition, they have intriguing nanoscale architectures that are based on thousands of coaligned nanorods.¹² These crystalline nanorods consist of metal carbonates (e.g., witherite in the case of barium) and have diameters of 5–50 nm and lengths of up to 500 nm. The overall alignment direction of the nanorods is parallel to the direction of the crystallization front that created the polycrystalline solid.

Some of the shapes observed in the biomorph system—especially the coral-like globules—serve as a motivation for this work. Unfortunately, the processes that give rise to biomorph growth are not well understood, partly due to the complicated chemical interaction of silicate species and metal carbonates and partly due to the complex reaction-diffusion coupling. Kaplan et al.¹⁵ modified a geometrical approach originally developed by Brower et al.¹⁶ for dendritic growth. In 2017, two of us proposed the use of nonlinear reaction-diffusion models for the description of biomorph growth.¹⁷ For the example of two-dimensional sheet growth, this approach yielded excellent agreement between the experimentally observed, leaflike shapes and the solid products formed in the wake of traveling excitation front. It also offered an explanation for the observed stalling of the crystallization front in terms of front defects resulting from local propagation failures.

Here, we continue our investigation of biomorph growth in the context of nonlinear reaction-diffusion processes with the aim to develop a simple model capable of producing three-dimensional structures reminiscent of biomorph corals and radial spikes. Unlike our earlier study,¹⁷ which used an unspecific, cubic nonlinearity in the reaction terms, we here propose and analyze a model that allows for a discussion in terms of a small number of simplified chemical reactions. We reemphasize that our basic assumptions are that the smoothly curved biomorph shapes are not directly determined by crystal growth phenomena but rather by reaction-diffusion processes and front instabilities. We suggest that this approach is justified by the clearly non-euhedral nature of the observed structures. Furthermore, other smoothly curved precipitation and mineralization structures have been successfully explained by reaction-diffusion models, specifically spiral and target patterns in the reaction of AlCl$3$ with NaOH^18,19 and on the nacreous surfaces of bivalve mollusks.^20,21

We adjust a 15 ml solution of BaCl$2$ (19.1 mM) and Na$2$SiO$3$ (8.2 mM) to a pH of 11.5 and place it in a 100 ml beaker. Then, a glass microscope slide is inserted at a tilt angle of $27°$ with respect to the gravity and the system is covered with a lid resting on top of the beaker. After 4 h, we terminate the biomorph precipitation by extracting the microscope slide and rinsing with de-ionized water. The microstructures on the slide are dried under ambient conditions and inspected using an optical microscope. In addition, samples are sputter-coated with iridium and then studied with a FEI Nova 400 field emission scanning electron microscope (SEM) operating at 7 kV.

In the absence of added carbonate ions, biomorph growth depends directly on the influx of CO$2$ from the gas phase. In the alkaline medium, CO$2$ reacts to carbonate ions and ultimately forms the barium carbonate structures. Along the tilted microscope slides, the carbonate concentration decreases with increasing solution depths, thus, providing different growth conditions. The biomorphs closest to the meniscus are star-shaped structures with three to five rodlike arms that extend radially away from a common starting point. The number of arms increases with increasing depths, and in addition, the structures change from nearly planar, substrate-bound units to three-dimensional, hemispherical objects. A typical example for the latter biomorphs is shown in Fig. 1(a), which features about 40 rodlike extensions.

Deeper into the solution, biomorphs develop sheetlike subunits and rods turn into tubes, urns, and cones [Figs. 1(b) and 1(c)]. For the investigated reaction time of 4 h and the employed conditions, we find the limit of biomorph formation at a depth of about 6.5 mm. Near this boundary, the structures [Fig. 1(d)] are reminiscent of corals and best described as folded, labyrinthine sheets within a hemispherical envelope. The spacing or wavelength of the sheets within these structures also varies with depth.

For SrCO$3$-based biomorphs grown in silica gels, Terada et al.²² further distinguished between coral-like and “petal-like” structures. The latter develop from biomorphs similar to the one in Fig. 1(d) at later reaction times and might simply be examples of extreme pattern wavelengths that ultimately terminate the growth. Noorduin et al.²³ produced coral-like structures in microfluidic reactors and noted an interesting thickening of the side walls. This feature is also discernible in Fig. 1(d) and reminiscent of curling sheet edges.^8,9,17

We note that all of these hemispherical biomorphs form in the presence of a vertical gradient in carbonate concentration; however, due to the small size of the structures, the actual concentration differences across the individual biomorphs are essentially negligible. Accordingly, no noticeable differences are observed when the top (closer to air) and bottom (farther from air) sides of a given biomorph are compared. Nonetheless, this particular class of experiments clearly unfolds an important chemical parameter (i.e., the carbonate concentration) in space and thus provides a convenient cross section through the system’s morphological parameter space.

Our model is based on three reaction steps that describe the conversion of a reactant species $A$⁠, through an intermediate $M$⁠, to a solid product $P$

The model species A is closely related to BaCO$3$(aq) that prior to biomorph growth has formed a supersaturated solution, although the complex involvement of silicate species complicates this simple chemical interpretation.¹¹ We currently consider the intermediate species $M$ to be nanoparticles that aggregate in front of or attach to the sheet edge, but other assignments within the still unresolved chemical processes are possible. After the intermediate $M$ has formed in an autocatalytic step, Eq. (1), the formation of the product occurs via two reaction channels: a simple, linear deposition equation (2) and a second-order process equation (3) that we interpret as surface-induced deposition, i.e., the attachment of nanoparticles onto or a molecular growth mechanism on the solid structure. Equations (1) and (2) also define the reaction scheme of the well-studied Gray–Scott model.^24,25 However, in our model, the simple precipitation step [Eq. (2)] is controlled by the solubility (or aggregation) threshold $m∗$ according to the Heaviside function $Θ(m−m∗)$⁠, where $m$ denotes the local concentration of species $M$⁠. Likewise, $a$ and $p$ describe the concentrations of $A$ and $P$⁠, respectively, and we obtain

Here $t$⁠, $DA$⁠, and $DM$ are time and the diffusion coefficients of $A$ and $M$⁠, respectively. The solid product $P$ is immobile but does not affect the diffusion of the other species. The rate constants $ka$⁠, $kp$⁠, and $ks$ as well as the threshold value $m∗$ are positive constants.

Note that the model considers neither an influx nor an outflux of matter, in analogy to biomorph growth, which can occur in closed systems. As such, it is concerned with the conversion of $A$ to $P$⁠. The backbone of our nonlinear reaction-diffusion model for biomorph formation is the Gray–Scott kinetics supplemented with a threshold-dependent conversion of the intermediate species $m$ into the precipitate $p$ as well as an enhanced precipitation step at places were solidification has already occurred [Eqs. (5) and (6)].

To reduce the number of parameters in Eqs. (4)–(6), an appropriate scaling is introduced by referencing all concentrations to the initial concentration $a0$ and by defining dimensionless parameters

Using this non-dimensionalized system, we perform a set of three-dimensional simulations that produce structures with diverse patterns depending on the chosen parameter values.

We emphasize the role of lateral wave instabilities for biomorph morphology at the micrometer scale. These instabilities, well-known for the Gray–Scott model, cause periodic rupture and segmentation of traveling M-waves. In detail, we will study this aspect in a separate paper. Here, we only note that the kinetic parameter $k$ and a sufficiently small ratio $δ$ of the diffusion coefficients of $A$ and $M$ set the characteristic length scale for the instability of traveling M-waves and hence the basic microscopic length scale of the biomorph to be calculated from the Kuramoto–Sivashinsky equation.²⁶ 

All three-dimensional simulations were performed using the forward Euler method for time integration and a discretized Laplacian with a 19-point stencil on a regular cubic lattice (256 $×$ 256 $×$ 256 lattice points) with an edge length of 2.56 space units. The space and time steps are set to $Δx$ = $Δt$ = 0.01. The diffusion ratio is kept constant at $δ$ = 0.02 with the specific diffusion rates chosen as $DA=10−4$ and $DM=2×10−6$⁠. The initial condition is the homogeneous state $(a~,m~,p~)$ = (1,0,0) locally perturbed to (0.7, 0.3, 0.0) in a small rectangular volume measuring $0.05×0.05×0.03$ cubic space units and located in the bottom center of the lattice. In order to simulate a closed reaction system, homogeneous Neumann boundary conditions are imposed.

The three-dimensional structures we find in simulations with our model show good qualitative agreement with the experimentally observed silica-carbonate precipitation structures with coral- and rodlike morphology. Simulation and experiment share the hemispheric envelope and similar morphological details as, for example, funnel-like substructures. The rate constant $k$ and the ratio of diffusion coefficients $δ$ define the critical wave number of the lateral wave instability, which is of the order of the segmentation length and the minimum distance between neighboring $m~$ pulses.

Within the model, effects of the quadratic surface deposition and the threshold-dependent linear process—described by parameters $s$ and $mthres$⁠, respectively—on the morphology of the pattern can be studied separately. Figure 2 compares four patterns obtained after equal simulation time for identical $k$ and $δ$ values taken from parameter intervals, where lateral wave instabilities are supported in two space dimensions.

Numerical simulations show that stronger surface-assisted deposition, i.e., enhanced precipitation at locations where precipitate has already formed, leads to a finer segmentation of the $m~$-front and thus results in more delicate $p~$-structures [Figs. 2(b) and 2(d)].

The higher the solubility threshold $mthres$⁠, the more frequent morphological substructures differentiate in the pattern that are also observed in the real biomorph formation. Figures 2(c) and 2(d), for example, display a funnel-shaped substructure arising in response to an increase in the threshold. Comparison between Figs. 2(b) and 2(d) or 2(a) and 2(c) reveals a more rapid formation of funnels for larger $mthres$ values. We emphasize that funnel-like structures do not form in the Gray–Scott model, Fig. 2(a), after equal simulation times.

Figure 3 shows numerical simulations for different parameter regimes that display distinct morphological features. In panel (a), rodlike arms extend from the area of local perturbation and grow out radially. In panel (b), in contrast, large bulges extend from the perturbation and the surface subsequently splits up into growing sheets and cones. Both cases maintain the hemispherical envelope that is also seen in Fig. 2.

Biomorphs grown in simple chemical experiments exhibit various morphological features like rods, cones, funnels, and sheets. Qualitative considerations suggest a possible connection to certain instabilities in nonlinear reaction-diffusion systems. For example, biomorph rods could result from a localized reaction zone traveling outwards while leaving behind a solid reaction product. Other features are reminiscent of Turing and labyrinthine patterns. With few exceptions,²⁷ the experimental analysis of such instabilities has focused on systems that are confined to two space dimensions and conditions that maintain the pattern via a global inflow of reactants and outflow of products. In contrast, biomorph growth is three-dimensional and does not depend on this exchange of matter as the structures can form in thermodynamically closed systems²⁸ (although in our experiments CO$2$ does enter the solution).

Here, we have proposed a simple three-step reaction model that, inspired by the above considerations, allows numerical simulations of nonlinear processes that form an immobile product by the consumption of a finite amount of reactant. Accordingly, the complex precipitation patterns arise in a closed system in response to a supercritical perturbation of a metastable state. Therefore, these patterns can be viewed as frozen transients of the intermediate $M$ during relaxation to a stable equilibrium state.

In addition, our study reproduces many of the experimentally observed features in three-dimensional numerical simulations, e.g., the funnels in Fig. 2(c). Other similarities include the structures with rodlike morphologies [compare Figs. 1(a) and 3(a)] and coral-like morphologies [compare Figs. 1(d) and 3(b)]. Future work should further investigate such similarities and the role of lateral wave instabilities. We believe that in this context, comparative studies to the well-established Gray–Scott model will provide a valuable starting point.

This paper is dedicated to Ken Showalter on the occasion of his 70th birthday. The material is based upon work supported by the National Science Foundation (NSF) under Grant No. 1609495 to O.S. and a Graduate Research Fellowship No. 1449440 to P.K. The Condensed Matter and Material Physics (CMMP) User Facility at Florida State University provided access to the SEM instruments. A.-K.M. and H.E. acknowledge funding by the German Science Foundation in the framework of Sfb 910 and GRK 1558.