https://orcid.org/0000-0002-3069-7967
Figures
Abstract
Over the last 40 years, applied mathematicians and physicists have proposed a number of mathematical models that produce structures exhibiting a fractal dimension. This work has coincided with the discovery that objects with fractal dimension are relatively common in the natural and human-produced worlds. One particularly successful model of fractal growth is the diffusion limited aggregation (DLA) model, a model as notable for its simplicity as for its complex and varied behavior. It has been modified and used to simulate fractal growth processes in numerous experimental and empirical contexts. In this work, we present an alternative fractal growth model that is based on a growing mass that bonds to particles in a surrounding medium and then exerts a force on them in an iterative process of growth and contraction. The resulting structure is a spreading triangular network rather than an aggregate of spheres, and the model is conceptually straightforward. To the best of our knowledge, this model is unique and differs in its dynamics and behavior from the DLA model and related particle aggregation models. We explore the behavior of the model, demonstrate the range of model output, and show that model output can have a variable fractal dimension between 1.5 and 1.83 that depends on model parameters. We also apply the model to simulating the development of polymer thin films prepared using spin-coating which also exhibit variable fractal dimensions. We demonstrate how the model can be adjusted to different dewetting conditions as well as how it can be used to simulate the modification of the polymer morphology under solvent annealing.
Citation: Mulder K, Lee SM, Chen W (2024) A triangular model of fractal growth with application to adsorptive spin-coating of polymers. PLoS ONE 19(2): e0298916. https://doi.org/10.1371/journal.pone.0298916
Editor: Malgorzata J. Krawczyk, AGH University of Science and Technology Faculty of Physics and Applied Computer Science: Akademia Gorniczo-Hutnicza im Stanislawa Staszica w Krakowie Wydzial Fizyki i Informatyki Stosowanej, POLAND
Received: August 30, 2023; Accepted: January 31, 2024; Published: February 23, 2024
Copyright: © 2024 Mulder et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and its Supporting Information files.
Funding: SL was funded by the National Science Foundation, grant number DMR-1807186 (nsf.gov). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
In nature as well as in human-determined realms, a surprising number of structures exhibit self-similar and scale-invariant features. Historically, such structures sparked the interest of mathematicians as they defied our conventional notion of dimension. Dubbed “fractals” by Benoit Mandelbrot in his seminal book The Fractal Geometry of Nature [1], rather than playing by the normal rules of Euclidean geometry in which idealized objects (e.g. lines, planes, and spheres) have integer dimension, these structures exhibit a fractional dimension. (In this work, fractal dimension will generally refer to Minkowski-Bouligan (or box-counting) dimension.) Since the emergence of the concept, fractal structures have been discovered in nearly every field of science [1,2]. Many physical materials exhibit a fractal structure [3–5], and bacteria and other organisms follow fractal growth patterns [6,7] as do lightning [8] and river networks [9]. Fractal patterns have been discovered on Mars [10] and in patterns of urban migration and development [11]. Dendritic growth is a particularly common form of fractal growth [2,12,13].
Due to the ubiquitous nature of fractal structures, including in the material sciences, mathematical models that produce fractal growth are of particular interest with many models focused on the aggregation of particles [2,14]. Early models include the Eden model of random growth on the boundary of a mass [15] and the Ballistic Aggregation model in which “sticky” particles are stochastically shot or dropped onto a surface [16]. Other, slightly more complicated models include the Dielectric Breakdown (DB) model [8] and the Born model of fracture growth [17]. (See [14] for an historical overview.) But arguably the most successful model of fractal growth has been the Diffusion Limited Aggregation (DLA) model of Witten and Sander [14,18].
The DLA model is notable for its stark simplicity, motivated by the physical dynamics of “sticky” particles diffusing in a medium [19]. The model begins with a centrally located particle followed by a sequence of introduced particles on random walks (Brownian motion), either on a lattice or in Euclidean space. Each particle keeps moving until it encounters a particle already incorporated into the aggregation. At this point, it stops moving and becomes part of the aggregate. Despite its simplicity, the DLA model has been successful in replicating numerous observed phenomena including Martian araneiforms [10], breast lesions [20], and thin film “islands” of one material grown on another [4].
DLA structures tend to exhibit dendritic growth and have a fractal dimension. The exhibited dimension when grown in two dimensional space is 1.71 and seems to be a constant even though the model exhibits a high degree of variability over a range of parameter values or through relaxation of some assumptions such as lowering the probability of sticking below unity [19,21]. Other, more extreme, modifications of the model can result in structures with different fractal dimensions. Modifying the nature of the random walk [22] or the probability structure of where particles enter the space [20] can produce significant variability in structure and dimension. A further modification is to allow the clusters themselves to diffuse, a so-called cluster aggregation model, which results in a significant reduction in dimension [23].
Other modifications to the DLA model have been problem specific. Wang et al. [9] modified the movement of the particles to reflect geographical constraints in order to allow the DLA model to predict river network formation. Xia [24] added specific orientation criteria to encourage production of needle-like structures. And the dielectric breakdown model can be replicated by adding an electric field to the DLA model thereby increasing the probability of attaching to an extending tip [8]. In all three cases, the general fractal structure of the DLA model is preserved while the output is modified to meet specific criteria. More importantly, the overall conceptual simplicity of the DLA model is largely preserved even if the computational complexity is increased.
Other models of fractal growth tend to be significantly more complex as they seek to incorporate continuous processes and capture a more nuanced picture of reality. Examples include models of frost formation that seek to include fluid flow and phase field dynamics (e.g. [25]) and models of bacterial growth that incorporate nutrient gradients [6]. Polimeno, Kim, and Blanchette [26] added several new factors to the DLA model with the goal of making it more consistent with observed empirical phenomena. Such models are important and often include DLA-like mechanisms, but their greater realism comes at a cost of higher complexity and consequent reduction in tractability. As noted by many authors [27,28], we often learn the most from early, simple models, and gains in understanding and management can decline rapidly with increased complexity.
In this work, we seek to make an addition to the library of fractal growth models that is of intermediate complexity. It is more complex than the DLA model, but its definition and dynamics are conceptually and mechanistically straightforward. Here we introduce a 2D model of fractal growth that we believe is novel and different from cellular or particle-based models of growth and aggregation. Instead of modeling particles, we model the bonds formed by a growing mass as it incorporates particles from the surrounding medium. Geometrically, a particle attaches to the mass by linking to the two endpoints of an existing link forming a triangle. Additionally, as the mass is forming it is also contracting. Particles are initially stationary. When they bond to the mass, they become mobile and are pulled toward any particles they are attached to. This contraction process continues until the mass solidifies and further movement ceases.
Below we introduce the model, explore its properties including a variable fractal dimension, and give an example of its application by showing how it replicates the growth of dewetted polymer structures.
Model definition
Model structure
Our model consists of two primary components—points in the plane (representing particles) and triangles consisting of three linked particles which represent bonds connecting these particles. The model is created iteratively by attaching a previously unattached particle to both ends of an existing link thereby creating a new triangular bond. Links that are candidates for forming new triangular bonds are considered “active”. When a new triangle forms, the initiating link is turned off while the other two links that form the triangle are now active.
Model initiation.
Because our model frequently exhibits radial growth, and one of our applications is to adsorptive spin-coating of polymer solution, we begin with a random, uniform distribution of particles on a disk. The number of particles is set by a parameter, and the radius of the disk is set so as to achieve a specified particle density. Throughout this paper, unless otherwise mentioned, the model has 40,000 particles with an average density of 100 particles per square unit. This yields a radius of 11.3 units. One particle is placed at the origin and a starting triangle comprised of three links is formed by connecting that particle to its two closest neighbors. The three initial links are all set to be active.
The model runs iteratively, alternating between a growth step in which active links have the potential to bond to new particles and a contraction step in which attached particles adjust their position seeking to minimize the distance to their linked neighbors. The two steps continue to alternate until either all particles have been attached or no active links remain.
Growth step.
Each time step, every active link has a specified probability (bond-rate) of seeking to bond with a previously unattached particle. This parameter determines the rate at which new bonds are formed under ideal bonding conditions. The number of links seeking to bond is determined by a random draw from a binomial distribution and the specified number of links is randomly chosen from all active links. Links seek to bond in a random order. When a link seeks to form a bond, the model searches for the unattached particle with the shortest average distance to the endpoints of the link. If the average distance to this particle is below a given threshold (bond-range), then both endpoints of the link form links with the new particle forming a new triangular bond. Each new link is then activated. Regardless of whether a bond is formed, the initiating link is deactivated and cannot be activated again. If a link is deactivated without bonding, we refer to it as a dead link. Otherwise, it becomes an internal link.
Contraction step.
After the growth step is complete and all selected links have been deactivated, the model experiences a contraction step as if the links between particles were ideal springs. When a particle is connected to the network, it is mobile for a specified number of time steps (mobility) after which it remains fixed in location. During the contraction phase, all mobile particles are pulled closer to their link neighbors (particles they are linked to) by experiencing a displacement given by where N is the set of all particles linked to the target particle and
is the vector extending from the target particle to particle i. We consider this to be a discrete approximation to a continuous contraction process.
When an attached particle has experienced contraction for a number of time steps equal to the mobility, it then becomes fixed in its location though it will continue to exert a force on any of its link neighbors that are still mobile. A stylized representation of the initial growth and contraction steps is shown in Fig 1. When no new triangular bonds are possible because no unattached particles remain within the bonding range, the model then continues to experience contraction for a number of time steps equal to the mobility thereby allowing recently attached particles to experience full contraction. Model parameters are reported in Table 1.
Each connection step, active links are randomly selected (colored to blue) and connected to the nearest particle if the particle is within bonding range. The “dead link” in the last window attempted to bond but had no unattached particles within the bonding range. In the contraction step, all attached particles that are mobile experience a force in the direction of the particles they are attached to which leads to the network contracting. The amount of contraction has been exaggerated for demonstration purposes.
Example models
Figs 2 and 3 show a series of intermediate stages and the final output for two different runs of the model. (All model images reported in this paper will show the links connecting particles since it is the bonds that are being modeled. The particles themselves are hidden.) For demonstration purposes, these models are smaller (5,000 particles) than the other models shown in the paper (40,000 particles) allowing the entire growth process to be seen. The masses generally grow radially, although some branches in Fig 2 extend back in toward the center to reach unattached particles. As the model grows, branches closer to the center complete their contraction while branches at the outer edge display the denser structure that precedes full contraction. The maximum bonding distance is greater in Fig 3 which results in fewer dead links and more direct, radial growth.
Parameter values are bond-rate = 0.05, mobility = 50 and bond-range = 0.2. Output is shown at times 150, 300, 450, 600 and 963 (the final time). The final image (lower right) is the output at time 963 zoomed in to demonstrate the triangular nature of the structure.
Parameter values are bond-rate = 0.02, mobility = 50 and bond-range = 0.5. Output is shown at times 150, 300, 450, 600 and 766 (the final time). The final image (lower right) is the output at time 766 zoomed in to demonstrate the triangular nature of the structure.
Model behavior
Determinants of behavior.
Model behavior is primarily dictated by two factors. The first is the average growth rate in active links as a function of the number of attached particles. This is directly related to the probability that an active link forms a new bond when it is selected. In order for a link to bond, an unattached particle must be within the maximum bonding distance (bond-range). This is in part determined by the ratio of bond-range to the average distance between particles at the start of the model. In this analysis, the latter is kept constant by keeping the particle density constant, thus bond-range is the sole parameter governing this ratio. However, each newly attached particle produces two active links, an exponential growth process that outpaces the growth in the number of particles within a given range as the structure grows. Thus, as closer unattached particles become linked, the required bonding distance steadily increases and will reach a point at which some links fail to bond and become dead links.
Contraction of the model also causes the distance between active links and unattached particles to increase as mobile particles and their links are pulled toward the center. This will increase the probability of a dead link forming. In particular, how many contractions the endpoints of a link experience before that link seeks to bond determines how far toward the center and away from unattached particles a link moves. The number of contractions a link experiences before attempting to bond is determined by how long it waits to be selected for bonding (inversely proportion to bond-rate) and the mobility (which determines the total number of contractions a particle experiences). If mobility is low, then active links do not move very far regardless of the bond-rate. But if mobility is high and the bond-rate is low, then links may experience significant contraction before seeking to bond potentially moving beyond the bond-range.
The second determinant of model behavior is the amount of contraction that takes place after bonding which is determined by mobility. Models with high mobility experience significant contraction which leads to shorter side branches and lower density in the final structure. Models with low mobility have longer side branches and higher density.
Together these two determinants lead to a heuristic for predicting model growth. First, the degree of curliness and spatial variability caused by dead links can be adjusted by increasing or decreasing the active link growth rate. This can generally be effected by adjusting the bond-range although changing the bond-rate can also affect this if the mobility is sufficiently high. Once the proper growth rate has been determined, the relative density and average length of branches in the structure can be adjusted by varying the mobility.
Descriptive model statistics.
We calculate five statistics that describe model growth and final structure. First, we directly estimate the active link growth rate. If there are no limits to bonding, then every active link bonds and each newly bonded particle adds two active links. However, generally one or more links fail for each link that succeeds in attaching to a new particle. We calculate the ratio of the change in active links divided by the change in attached particles between 5000 and 10,000 timesteps. Generally, we see an average increase in active links per attached particle of between 0.025 and 0.2. However, if too many links are unable to bond, the growth rate can be negative and the structure will fail to expand into all areas and some particles may fail to attach. Consequently, as a second statistic, we calculate the percentage of initial particles that become attached, referred to as percent attached.
Third, since unattached particles lie outside the area of current growth, when the growth rate in active links is high, we expect the growth to have a strong, radial pattern. Conversely, when more deadlinks appear, we expect growth to show greater variability in direction and thus be less radial. One way to measure radial growth is to estimate the growth rate in the number of links intersecting a concentric circle of the disk as a function of the radius of the circle. This will generally be a linear function of the radius, but the slope of this line will be lower for networks that are curvier with less direct radial motion. Plotting the number of intersecting branches as a function of radius for a variety of models supports this claim. We refer to this statistic as the branch-slope.
Fourth, we calculate the variability in network structure by calculating the standard deviation in the number of links each particle is connected to, referred to as the standard deviation of the network degree. Finally, to assess the overall level of contraction, we report the average branch length. Higher mobility should lead to shorter branch lengths. However, this should also be influenced by bond-range with greater bonding ranges leading to longer links.
Example models.
This leads to two primary axes describing model behavior. The first axis is bond-range, which dictates how direct the radial growth is. Low radial growth results in curvier branches and a higher proportion of dead links. The second axis is mobility, which controls the relative density and how long the side branches are. Fig 4 shows nine examples of model output arranged according to these two axes. As the maximum bond range increases, growth becomes straighter and more radial. As the mobility increases, the side branches become shorter and the overall structure becomes less dense.
The vertical axis controls the mobility which dictates the degree of contraction the structure experiences as it grows. The horizontal axis controls the maximum bonding range which dictates the growth rate in active links and subsequent radial growth.
Fig 5 shows heat maps of five model statistics and fractal dimension for five values of mobility and seven values of bond-range. Bond-rate was set to 0.1. The maximum bonding range increases from left to right and the mobility increases from top to bottom in each map. Parameter values were chosen to give provide a thorough coverage of parameter phase space based on exploratory analysis.
Bond-rate = 0.1. From left to right in each map, bond-range values are 0.165, 0.175, 0.185, 0.2, 0.25, 0.3 and 0.5. From top to bottom in each map, mobility values are 25, 50, 100, 300, and 500. Heat maps are color coded according to the range of each variable with the lowest values shaded green and the highest values shaded red. The values are the averages for 50 model runs for each combination of parameters. Data are available in S1 File.
All five model statistics are significantly impacted by bond-range. Only branch-slope and active branch length are affected by mobility. The growth rate, the proportion of particles attached, and the variability in the network degree of each particle all increase with bond-range and are unaffected by mobility. The branch-slope increases with bond-range as growth becomes more radial and diminishes with mobility as the side branches become shorter. The same is true for average branch length. Similar figures were calculated for bonding rates of 0.02 and 0.5 (not shown). The same patterns with regard to bond-range and mobility were observed. The growth rate, branch-slope, and average branch length all increased with the bonding rate. The proportion of particles attached increased when bond-rate went from 0.02 to 0.1 but did not change when it was increased further to 0.5. The variability in network degree increased with bond-rate for lower levels of bond-range and decreased at higher levels of bond-range.
When the growth rate in active links is too low, the model has the potential to die out before utilizing all (or even most) of the particles, and even when the space is fully colonized, the path taken from the center to a branch tip can be very circuitous. Fig 6 shows two such outcomes, one with a lower mobility (higher density) and one with a higher mobility (lower density). In both situations, the average growth rate is barely positive leading to a high proportion of dead links. As a result, growth is often more similar to a random-walk than it is to radial growth. Indeed, in the model with higher mobility, initial growth was restricted to just two branches and as a result the network is comprised of two distinct clusters. The model with lower mobility is composed of only three clusters although this is more difficult to discern visually.
The average growth rate is ~ 0.015 which is near the lower end for viable growth or, on average, a slight increase in the number of active lengths for each new particle attached. Higher mobility (300) in the model on the left leads to a lower-density network while lower mobility (25) in the model on the right leads to a higher-density network. Parameter values for both runs are bond-rate = 0.1 and bond-range = 0.16.
Fractal nature of model growth and variability in the fractal dimension
Image analysis using the box-counting method for estimating fractal dimension [29] shows model output to consistently have a fractional dimension. Because model output consists of triangular bonds, dimension could not be estimated using the coordinates of the particles (the vertices of the bonds). Instead, model output was converted to a binary image and dimension was estimated using the Fraclac plugin [30] in the image processing software ImageJ [31]. To ensure dimension was comparable across models, image scale and resolution were kept constant.
Analysis across models shows the mean and variance of the fractal dimension both vary with parameter settings, a feature somewhat unusual for a relatively simple growth model. Observed fractal dimensions range from 1.50 to 1.84 with most structures having a dimension greater than 1.70. Average fractal dimensions for seven different bonding ranges and five different mobilities are shown in Fig 5. Fractal dimension increases with bond-range and decreases with mobility. At lower levels of bond-range, dimension increases with bond-rate.
Lower dimensions occur when the average growth rate is close to zero, which can be achieved generally by decreasing bond-range (Fig 5). When the maximum bonding range is greater than 0.2, fractal dimension is generally greater than 1.75. Lower fractal dimensions are generally characterized by a high proportion of unattached particles. Fig 7 shows three representative model outputs with different fractal dimensions.
From left to right, the fractal dimensions using the box-counting method are 1.65, 1.75, and 1.79. Parameter values from left to right are as follows: {bond-rate = 0.64, mobility = 456 and bond-range = 0.15}, {bond-rate = 0.01, mobility = 31 and bond-range = 0.25}, and {bond-rate = 0.02, mobility = 19 and bond-range = 0.3}.
It is of some interest that the same set of model parameters can lead to significant variation in dimension. Fig 8 shows the range of fractal dimensions achieved over multiple runs when all parameters are kept fixed. For the runs shown, only bond-rate was varied, and multiple runs were conducted for each setting of the bonding rate. The average fractal dimension increases with the bonding rate, but significant variability is seen for low and moderate bonding rates where a high proportion of particles remained unattached. This is likely due to the higher variability of the negative binomial distribution with lower bond-rate. Links that wait longer to attempt bonding experience more contractions and thus a higher likelihood of failing to bond. At higher bonding rates, the variability in dimension decreases significantly.
Other parameter settings are as follows: bond-range = 0.18, mobility = 100. The sample size varies as some runs were discarded due to a high proportion of unattached particles. Sample size was 7–13 for bonding rates less than 0.025 and 20 for the other groups. Data are available in S2 File.
Simulation of polymer development
Fractal polymer growth
Model output for our triangular growth model produces patterns very similar to those observed in the morphologies of polymer thin films prepared by adsorptive spin-coating. In particular, Qi et al. [32], by adjusting the amount of poly(vinyl alcohol) (PVOH99%H: 89−98 kDa and 99% hydrolyzed) solution deposited on a polydimethylsiloxane substrate (PDMS49k: 49 kDa), produced polymer patterns very similar to our model output both in terms of structure and fractal dimension (see Fig 9). Poly(vinyl alcohol) is synthesized by hydrolyzing poly(vinyl acetate) to convert a large portion of acetate groups to alcohol groups, -OH. On average ~99% of the repeat units in PVOH99%H consists of -OH groups, which are capable of hydrogen bonding interactions and crystallization. The remaining ~1% of the repeat units contains acetate groups, which can engage in hydrophobic interactions. The formation of the fractal features in the PVOH99%H thin films is driven by crystallization during the drying process [33]. PDMS49k is a hydrophobic, flexible support that facilitates the dewetting of the PVOH99%H thin films [32–34]. In the study by Qi et al., the PVOH99%H polymer thickness (amount of polymer per unit area) was the same under various experimental conditions, similar to the constant particle density in our model.
5−50 μL PVOH99%H solutions spin-coated on PDMS49k substrates (1.4 cm x 1.4 cm) at 6000 rpm. Each image (12.5x) was converted to binary after background subtraction so that fractal dimension (D) and lacunarity (L) values could be determined. The scale bars are 2 mm in length. Reprinted from Qi et al. [32] under a CC BY license, with permission from the American Chemical Society, original copyright 2019.
In Fig 9, there is a clear relationship between drop size and fractal dimension as well as with density. The primary cause of this relationship is the fact that smaller drops have a much longer drying time due to a greater proportion of the mass experiencing weaker centrifugal force. This gives particles more time to experience contracting forces prior to solidification. In the model, increased drying time is simulated by increased mobility, and we do indeed see a similar relationship between mobility and dimension/density (see Figs 4 and 5).
Moreover, Qi et al. observed an inverse relationship between the fractal dimension and the lacunarity of the images. Lacunarity (λ) was introduced by Mandelbrot [1] as a complementary measure to dimension for describing fractals. It captures the texture of a fractal by measuring the spatial variability in image density. To better compare our model output to the spin-coated polymers, we estimated lacunarity for a range of parameters using the Fraclac plugin [30] in ImageJ [31]. Our model also produces an inverse relationship between dimension and lacunarity that closely matches values reported by Qi et al. (Fig 10). Generally, we see a negative linear correlation between lacunarity and dimension. Model parameters varied across their full range with the exception of bond-range. When the bonding range was 0.2 or greater, output tended to cluster around a dimension of 1.77 to 1.79 and lacunarity between 0.28 and 0.3, values beyond the range of those observed in the polymer data in Fig 9. For the parameters chosen, while the sample size from Qi et al. is small, all six samples lie within the general range of values seen in Fig 10.
Model parameters were randomly drawn for each run. Bond-rate varied between 0.01 and 1.0, mobility between 15 and 500, and bond-range between 0.16 and 0.20. Two of the empirical samples had identical values. Data are available in S3 File.
Additional aspects of polymer growth
There are additional aspects of polymer pattern development that the model is able to capture. Fig 11 demonstrates two distinct annular regions in polymer pattern formation that result from a difference in spin rates and consequent drying times. The inner region experiences a longer drying time due to reduced centrifugal force leading to a circular boundary between a looser, lower dimensional inner core and an outer band displaying a tighter, higher dimensional pattern that results from a shorter drying time. Faster spin rates result in the region with shorter drying times extending farther in toward the center. Fig 12 shows a similar pattern produced by our model. The change is the result of a difference in mobility between particles that are above and below a fixed distance from the center. In particular, the radius of the figure is 11.3 units, and particles closer than 6.5 units to the center have a higher mobility than those more than 6.5 units from the center. Thus, the center models a longer drying time meaning more ability for the polymer to contract. The similarity between our model and the polymer system suggests that the mobility parameter simulates the experimental drying time.
A 200 μL PVOH99%H solution (drop diameter: 9.2 ± 0.1 mm) spin-coated on PDMS49k substrates (1.4 cm x 1.4 cm) at 900−6000 rpm. Each image was stitched together from sectional optical images (50x) using Autopano Giga. Reprinted from Qi et al. [32] under a CC BY license, with permission from the American Chemical Society, original copyright 2019.
The radius of the disk is 11.3 units. Particles within 6.5 units of the center have mobility of 200. Those further away have mobility of 20. This simulates a change in drying time during polymer pattern formation. Longer drying times cause more extensive polymer aggregation and contraction resulting in a lower fractal dimension. Remaining model parameters are: bond-range = 0.25, and bond-rate = 0.1.
The other phenomenon our model replicates with regard to polymer morphology is the rearrangement of polymer thin films when they are exposed to water and then allowed to reform, a process known as solvent annealing. The first two panels in Fig 13 show a polymer before and after solvent annealing. The third panel is an overlay of the two images for comparison. The applied water appears to have liberated a portion of the particles, modifying the original structure. When the water is removed, the liberated polymer reforms. The new growth is less curvy than the original structure suggesting an increase in bond-range as shown in Fig 4. The soluble portion of the polymer is less crystalline and differs from the insoluble portion. This difference most likely gives rise to different bond-range values between the original and the new structures.
Optical images (500x) of dewetted patterns obtained from adsorptive spin-coating of 400 μL PVOH99%H solution on PDMS49k substrates (1.4 cm x 1.4 cm) at 900 rpm before and after solvent annealing. Solvent annealing involved soaking samples in Milli-Q water for 1 h followed by drying under nitrogen gas. Overlaid optical image depicts polymer rearrangement as the result of solvent annealing—before image in magenta and after image in green. Fractal dimensions were 1.58 (+/- 0.07) before solvent annealing and 1.64 (+/- 0.07) after solvent annealing.
We model annealing in three stages. First, the normal model was run with the parameters set to replicate a polymer with the dimension and structure seen in the first panel of Fig 13 (before annealing). Second, we removed all of the bonds as well as a fixed proportion of randomly selected particles from the structure and allowed the original structure to reform. The removed particles were randomly placed on the original disk. This simulates the liberating of some of the polymer into solution while some of the original structure is retained. Finally, all external links were set to active, and the structure was allowed to finish growing. Two growth stages were necessary after annealing as it is hypothesized that many of the original bonds are maintained during annealing and many of the particles have limited mobility. Fig 14 shows the results of such a simulation. A similar deformation of the original polymer combined with higher-dimensional new growth is apparent. Fractal dimensions for the polymer before and after annealing were a bit higher than the average values observed experimentally, but still capture an increase in fractal dimension after annealing (see Fig 13).
Parameters for the left model were bond-range = 0.2, bond-rate = 0.02. and mobility = 100. 30% of the particles were liberated, and after annealing bond-range was increased to 0.5. The disks were cropped to replicate the images in Fig 13. Fractal dimensions were 1.71 and 1.75, slightly higher than those in Fig 13. Coloration of the before and after components in the overlaid image is the same as in Fig 13.
Discussion
Similarities and differences with other models of fractal growth
There are many mathematical growth models that exhibit a fractal dimension and produce structures that resemble empirical phenomena. Some are quite simple, such as the Eden model [15] and the DLA model [19], others are of medium complexity such as the Born model of fracture growth [17], the dielectrical breakdown (DB) model [8], and the model of bacterial growth of Ben-Jacob et al. [6], while others are numerically if not conceptually quite complex such as the mathematical model for Hele-Shaw cells [35,36] or phase models of ice development (e.g. 25). We feel that our model is comparable to the DLA model in presenting a conceptually straightforward model of growth. Computationally, the contraction step adds in some complexity but still remains significantly simpler than the numerical approximations required for the Born model or the DB model.
The variability in output of our model is comparable to other models. The bacterial growth model of Ben-Jacob et al. [6] also has two parameters that dictate model output, producing a range of behaviors not dissimilar from ours (see Fig 3 in [6]). And variations to the DLA model in terms of the probability distribution of incoming particles and random walk step size also produces a range of behaviors (see Fig 5 in [20]). Variability in these models is comparable to ours in that much of the variation concerns the number of branches and how dendritic the growth is in appearance. However, neither model is able to produce the kind of random walk behavior seen in our model output with low dimension.
One significant difference between our model and most fractal growth models is that mathematically, our structure is a triangular network with vertices of zero mass. In other words, our model produces bonds between particles rather than modeling the aggregation of particles themselves. This is not entirely unique. The Born model output consists of the bonds that break as the fracture spreads. L-system models produce geometric structures by converting a sequence of strings into a series of commands for drawing straight line segments, in effect producing a network [37,38]. In terms of polymer pattern development, the target of simulation for our model is the same as if the DLA model were being utilized, but the structure is being modeled as an expanding mass of bonds rather than an aggregating mass of particles.
Our triangular bond model simulates multi-body interactions in forming a realistic fractal network. A new particle in the vicinity of a growing front experiences multiple, attractive, intermolecular forces and can get pulled in to bond to the structure. Bond-range is adjustable to scale with the nature of intermolecular forces. Formation of solid structures, such as those prepared by adsorptive spin-coating of polymer solution, often involves solvent evaporation. Mobility simulates variation in drying time and extent of structure contraction. The key features in our model–multi-body bonding, bond-range, and mobility–are drawn upon physical parameters in real applications.
Finally, we would note that while the model we present here is inherently a 2D, planar model, it is readily extendable to higher dimensions. In three dimensions, bonds are formed by connecting a new particle to each vertex of an external triangular face on the growing mass. Contraction occurs precisely as with the 2D model, and initial runs we have conducted also show fractal growth. Extensions to higher dimensions would be analogous.
Models with variable dimension
The variability of the fractal dimension of our growth model is an important feature, especially since the dimension can be tuned with the model parameters. While the original DLA model appears to have a fractal dimension that is consistently 1.71 [19,21], variations of the model do exhibit other dimensions, and some exhibit dimensions that vary with the parameters. The generalized DLA model with different dimensional random walks achieves dimensions between 1.71 and 1.95 [22], and allowing the aggregates to also diffuse (cluster aggregation) reduces the dimension to approximately 1.5 [4,23]. DB models, which can be considered to be an extension of the DLA model, exhibit a variable and tunable fractal dimension between 1.6 and 2 dependent upon the strength of the interaction with the electrical field [8].
Variable dimension often comes at the cost of increased model complexity, both conceptually and computationally. Ehrl, Soos, and Lattuada [39] were able to produce clusters in 3D of any desired dimension from 2.2 to 3 by using an aggregation algorithm specifically focused on producing the desired dimension. The algorithm has value for modeling observed aggregations, however it is significantly more complex than aggregation models such as DLA and the DB model, and it lacks a direct connection to a known aggregation process [40]. Similarly, the bacterial growth model of Ben-Jacob et al. [6] successfully replicates bacterial growth patterns and has variable dimension less than 1.7. But again, model dynamics are quite complex. Perhaps the simplest model with tunable dimension is the DB model, but this model requires an iterative solution to the Laplace equation after each particle docks.
By comparison, our model is conceptually simple and also has a tunable dimension although with significant variability between model runs at low and medium dimensions. We are working on algorithms that would enable runs with a greater number of particles to determine the degree to which the stochasticity in dimension is inherent in the model. Having the ability to adjust dimension increases our ability to fit the model to specific phenomena as well as to test hypotheses related to those phenomena. This is demonstrated above by our ability to adjust dimension to match polymers with different drying times and consequently different fractal dimensions. Further, by associating such changes with specific model parameters, we can craft hypotheses about what elements of the system cause the change in dimension.
Comparable applications to polymers and other phenomena
Fractal growth models like ours have been used to model a wide variety of phenomena [2]. Zhang and Liu [12] modeled viscous fingering in liquids with varying fractal dimension using a modified DLA model. Other modifications of the DLA model have been used to model chemical vapor deposition [41], urban migration and growth patterns [11], general dendritic growth [13], and the evolution of river networks [9]. L-systems have been used to model many biological systems [37] including the growth of mycorrhizal fungi [38].
Polymer pattern development in particular has been modeled by several different fractal growth models. Gao et al. [3] demonstrated a similarity between fullerene doped polymer thin films and the cluster diffusion model [23]. Similarly, Amir, Ali, and Mohamed [42] modeled fractals in polymer electrolyte films using the DLA model. Morozova et al. [43] used a clustering model for charged nanoparticles that has good agreement with observed colloidal clusters and networks. On a finer scale, Monte Carlo based methods have been used to model individual polymer chains as lattice animals, a set of connected vertices on a hypercubic lattice. [44,45].
In a similar vein, we have demonstrated the ability of our model to create structures very similar to those produced by the spin coating of polymer thin films and subsequent dewetting. We are not only able to match the dimension and general shape of these polymers, but also, by making reasonable changes to model parameters, we are able to account for differences in polymer structure caused by variations in drying time or the application of a solvent and the resultant annealing. This suggests a deeper connection between our model and the dynamics of polymer crystallization and rearrangement, something we hope to explore in future work. Ideally, our model will be able to make predictions for polymer features under different conditions as well as lead to new hypotheses regarding changes to polymer dynamics under different experimental treatments.
More importantly, given the fundamental differences in approach between our model and other fractal growth models, it is our hope that this model will allow for the relatively simple and efficient simulation of experimental systems that are not amenable to being modeled with simple, aggregation-based models. We would add that our model can easily be extended to three dimensions where the activated “links” are triangles on the surface of the growing structure which connect to unattached particles by forming bonds between the particle and all three vertices of the triangle forming a tetrahedron. The structure of these 3D fractals is the subject of ongoing work.
Conclusion
The power of the diffusion limited aggregation model of Witten and Sander [18] lies in its tractability and simplicity. It is a model that can be effectively described to a wide audience and easily implemented in many different platforms. And yet, its behavior is varied and complex and able to simulate a wide variety of empirical phenomena. In this work, we have introduced a model of comparable (albeit greater) complexity that takes a fundamentally different approach to fractal growth, what might be dubbed a “reach-and-grab” model as opposed to a “sit-and-wait” model. Our model also exhibits a range of complex behaviors that appear similar to a variety of real-world phenomena. Model output demonstrates variability in fractal dimension that is correlated with model parameters. This leads us to hope that for specific situations like the development of polymer fractals using spin-coating, the model can be tuned to simulate specific features of the system and ideally lead to predictive models for yet-to-be-produced fractal structures.
Acknowledgments
We wish to thank two reviewers whose comments and critiques improved the work presented here.
References
- 1. Mandelbrot BB. The fractal geometry of nature. Vol. 1. WH freeman New York; 1982.
- 2.
Vicsek T. Fractal Growth Phenomena [Internet]. 2nd ed. WORLD SCIENTIFIC; 1992 [cited 2023 Jul 20]. Available from: https://www.worldscientific.com/worldscibooks/10.1142/1407.
- 3. Gao HJ, Xue ZQ, Wu QD, Pang SJ. 2D fractal pattern in fullerene doped polymer. Solid State Commun. 1996 Feb 1;97(7):579–82.
- 4. Hwang RQ, Schröder J, Günther C, Behm RJ. Fractal growth of two-dimensional islands: Au on Ru(0001). Phys Rev Lett. 1991 Dec 2;67(23):3279–82. pmid:10044692
- 5. Wang WJ, Xu X, Shao Y, Liao JW, Jian HX, Xue B, et al. Fractal Growth of Giant Amphiphiles in Langmuir-Blodgett Films. Chin J Polym Sci. 2022 Jun 1;40(6):556–66.
- 6. Ben-Jacob E, Schochet O, Tenenbaum A, Cohen I, Czirók A, Vicsek T. Generic modelling of cooperative growth patterns in bacterial colonies. Nature. 1994 Mar;368(6466):46–9. pmid:8107881
- 7. Wang J, Li X, Kong R, Wu J, Wang X. Fractal morphology facilitates Bacillus subtilis biofilm growth. Environ Sci Pollut Res. 2022 Aug 1;29(37):56168–77. pmid:35325386
- 8. Niemeyer L, Pietronero L, Wiesmann HJ. Fractal Dimension of Dielectric Breakdown. Phys Rev Lett. 1984 Mar 19;52(12):1033–6.
- 9. Wang S, Ji H, Li P, Li H, Zhan Y. Growth diffusion-limited aggregation for basin fractal river network evolution model. AIP Adv. 2020 Jul 22;10(7):075317.
- 10. Portyankina G, Hansen CJ, Aye KM. How martian araneiforms get their shapes: morphological analysis and diffusion-limited aggregation model for polar surface erosion. Icarus. 2020 May 15;342:113217.
- 11. Murcio R, Rodríguez-Romo S. Colored diffusion-limited aggregation for urban migration. Phys Stat Mech Its Appl. 2009 Jul 1;388(13):2689–98.
- 12. hua Zhang J, hua Liu Z. Study of the relationship between fractal dimension and viscosity ratio for viscous fingering with a modified DLA model. J Pet Sci Eng. 1998 Sep 1;21(1):123–8.
- 13. Ruzhitskaya DD, Ryzhikova YuV, Ryzhikov SB. Algorithms for Analyzing the Characteristics of Dendritic Structures. Bull Russ Acad Sci Phys. 2018 Nov 1;82(11):1375–8.
- 14. Meakin P. A Historical Introduction to Computer Models for Fractal Aggregates. J Sol-Gel Sci Technol. 1999 Aug 1;15(2):97–117.
- 15. Eden M. A two-dimensional growth process. Dyn Fractal Surf. 1961;4:223–39.
- 16. Meakin P, Ramanlal P, Sander LM, Ball RC. Ballistic deposition on surfaces. Phys Rev A. 1986 Dec 1;34(6):5091–103. pmid:9897896
- 17. Yan H, Li G, Sander LM. Fracture Growth in 2d Elastic Networks with Born Model. Europhys Lett. 1989 Sep;10(1):7.
- 18. Witten TA, Sander LM. Diffusion-Limited Aggregation, a Kinetic Critical Phenomenon. Phys Rev Lett. 1981 Nov 9;47(19):1400–3.
- 19. Witten TA, Sander LM. Diffusion-limited aggregation. Phys Rev B. 1983 May 1;27(9):5686–97.
- 20. Rashidnasab A, Elangovan P, Yip M, Diaz O, Dance DR, Young KC, et al. Simulation and assessment of realistic breast lesions using fractal growth models. Phys Med Biol. 2013 Jul;58(16):5613. pmid:23892735
- 21. Meakin P. Diffusion-controlled cluster formation in 2—6-dimensional space. Phys Rev A. 1983 Mar 1;27(3):1495–507.
- 22. Hentschel H. Fractal dimension of generalized diffusion-limited aggregates. Phys Rev Lett. 1984;52(3):212.
- 23. Meakin P. The structure of two-dimensional Witten-Sander aggregates. J Phys Math Gen. 1985 Aug;18(11):L661.
- 24. Xia Z. The growth simulation of pine-needle like structure with diffusion-limited aggregation and oriented attachment. RSC Adv. 2022;12(35):22946–50. pmid:36106001
- 25. Kim T, Henson M, Lin MC. A hybrid algorithm for modeling ice formation. In: Proceedings of the 2004 ACM SIGGRAPH/Eurographics symposium on Computer animation—SCA ‘04 [Internet]. Grenoble, France: ACM Press; 2004 [cited 2023 Jul 20]. p. 305. Available from: http://portal.acm.org/citation.cfm?doid=1028523.1028564.
- 26. Polimeno M, Kim C, Blanchette F. Toward a Realistic Model of Diffusion-Limited Aggregation: Rotation, Size-Dependent Diffusivities, and Settling. ACS Omega. 2022 Nov 15;7(45):40826–35. pmid:36406481
- 27. Lee DB Jr. Requiem for large-scale models. J Am Inst Plann. 1973;39(3):163–78.
- 28. Zellner ML, Milz D, Lyons L, Hoch CJ, Radinsky J. Finding the Balance Between Simplicity and Realism in Participatory Modeling for Environmental Planning. Environ Model Softw. 2022 Nov 1;157:105481.
- 29. Smith TG, Lange GD, Marks WB. Fractal methods and results in cellular morphology—dimensions, lacunarity and multifractals. J Neurosci Methods. 1996 Nov;69(2):123–36. pmid:8946315
- 30. Karperien A. FracLac for ImageJ [Internet]. 1999. Available from: http://rsb.info.nih.gov/ij/plugins/fraclac/FLHelp/Introduction.htm.
- 31. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012 Jul;9(7):671–5. pmid:22930834
- 32. Qi Y, Nguyen H, Lim KSE, Wang W, Chen W. Adsorptive Spin Coating To Study Thin-Film Stability in Both Wetting and Nonwetting Regimes. Langmuir. 2019 May 28;35(21):6922–8. pmid:31082251
- 33. Karki A, Nguyen L, Sharma B, Yan Y, Chen W. Unusual Morphologies of Poly(vinyl alcohol) Thin Films Adsorbed on Poly(dimethylsiloxane) Substrates. Langmuir. 2016 Apr 5;32(13):3191–8. pmid:27002807
- 34. Jiang Y, Minett M, Hazen E, Wang W, Alvarez C, Griffin J, et al. New Insights into Spin Coating of Polymer Thin Films in Both Wetting and Nonwetting Regimes. Langmuir. 2022 Oct 18;38(41):12702–10. pmid:36201003
- 35. Saffman PG. Viscous fingering in Hele-Shaw cells. J Fluid Mech. 1986 Dec;173:73–94.
- 36. Ben-Jacob E, Garik P. The formation of patterns in non-equilibrium growth. Nature. 1990 Feb;343(6258):523–30.
- 37. Lindenmayer A. Mathematical models for cellular interactions in development II. Simple and branching filaments with two-sided inputs. J Theor Biol. 1968 Mar 1;18(3):300–15.
- 38. Schnepf A, Leitner D, Schweiger PF, Scholl P, Jansa J. L-System model for the growth of arbuscular mycorrhizal fungi, both within and outside of their host roots. J R Soc Interface. 2016 Apr;13(117):20160129. pmid:27097653
- 39. Ehrl L, Soos M, Lattuada M. Generation and Geometrical Analysis of Dense Clusters with Variable Fractal Dimension. J Phys Chem B. 2009 Aug 6;113(31):10587–99. pmid:19594146
- 40. Thouy R, Jullien R. A cluster-cluster aggregation model with tunable fractal dimension. J Phys Math Gen. 1994 May;27(9):2953.
- 41. Li J, Chen M, Zhang C, Dong H, Lin W, Zhuang P, et al. Fractal-Theory-Based Control of the Shape and Quality of CVD-Grown 2D Materials. Adv Mater. 2019;31(35):1902431. pmid:31265203
- 42. Amir S, Ali SAH, Mohamed NS. Implementation of a diffusion-limited aggregation model in the simulation of fractals in PVDF-HFP/PEMA–NH4CF3SO3–Cr2O3 nanocomposite polymer electrolyte films. Phys Scr. 2011 Sep;84(4):045802.
- 43. Morozova SM, López-Flores L, Gevorkian A, Zhang H, Adibnia V, Shi W, et al. Colloidal Clusters and Networks Formed by Oppositely Charged Nanoparticles with Varying Stiffnesses. ACS Nano [Internet]. 2023 Jul 17 [cited 2023 Jul 21]; Available from: pmid:37459253
- 44. Hsu HP, Grassberger P. A Review of Monte Carlo Simulations of Polymers with PERM. J Stat Phys. 2011 Aug 1;144(3):597–637.
- 45. Binder K, Paul W. Recent Developments in Monte Carlo Simulations of Lattice Models for Polymer Systems. Macromolecules. 2008 Jul 1;41(13):4537–50.