Predator-prey dynamics stabilised by nonlinearity explain oscillations in dust-forming plasmas

Dust-forming plasmas are ionised gases that generate particles from a precursor. In nature, dust-forming plasmas are found in flames, the interstellar medium and comet tails. In the laboratory, they are valuable in generating nanoparticles for medicine and electronics. Dust-forming plasmas exhibit a bizarre, even puzzling behaviour in which they oscillate with timescales of seconds to minutes. Here we show how the problem of understanding these oscillations may be cast as a predator-prey problem, with electrons as prey and particles as predators. The addition of a nonlinear loss term to the classic Lotka-Volterra equations used for describing the predator-prey problem in ecology not only stabilises the oscillations in the solutions for the populations of electrons and particles in the plasma but also explains the behaviour in more detail. The model explains the relative phase difference of the two populations, the way in which the frequency of the oscillations varies with the concentration of the precursor gas, and the oscillations of the light emission, determined by the populations of both species. Our results demonstrate the value of adopting an approach to a complex physical science problem that has been found successful in ecology, where complexity is always present.

no model has successfully explained this aspect of the behaviour. It is notable that whilst the ion current and dust particle mean charge, and therefore size, oscillate, the electron temperature may rise continuously, not correlating with cyclic behaviour 25 . Therefore, it is appropriate in the interests of simplicity to neglect electron temperature when modelling cyclic behaviour, and concentrate on the oscillatory aspects of two interacting populations.
The predator-prey problem is addressed in the classic works of Lotka 27 , Volterra 28 , and May 29 . Behaviour in which the predator and prey populations oscillate indefinitely is the solution of the original problem as set out by Lotka and Volterra [29][30][31][32] . A foundational study linked oscillations in the populations of the snowshoe hare and its predator the lynx as measured over 200 years in Canada, by noting periodically recurring pronounced variations in the quantities of fur pelts received by the Hudson's Bay company 33,34 . As originally conceived, the Lotka-Volterra solution was expressed as the following pair of coupled linear first-order differential equations 27,28 : where x 1 is the prey population, x 2 is the predator population, α and β are coefficients describing the birth rate and predation rate of the prey, while γ and δ describe the death rate and birth rate of the predator. The model has oscillating solutions in which both populations increase and decrease periodically and out of step. Predator-prey modelling has been used to extract understanding from complexity in many contexts: chemical reactions 35 , astronomy 36 , economics 37 , and evolutionary game theory 38 . In climate modelling, the manner in which the structure and behaviour of atmospheric clouds are treated has a significant effect on overall climate sensitivity. This problem has been framed as a predator-prey situation where the rain acts as the predator, the cloud the prey, and the aerosol particles modulate the predator-prey response. Cloud droplets form on aerosol particles and then coalesce to form rain drops, which then consume the cloud and remove its cooling ability. The feedback mechanisms associated with the aerosol-cloud-precipitation system are many and complex, and had proven difficult to incorporate into climate models 39,40 . The predator-prey solution enabled the creation of a predictive model, which has proven itself reliable and accurate 41 .
Our model. We now describe dust-forming plasmas as a predator-prey system, in which the dust particulates are the predators and the electrons the prey. Figure 1 depicts the processes by which the species interact. The dust depletes the plasma by causing electrons and ions to recombine on its surfaces. The dust is therefore a predator of the plasma. Another process at work is the loss of electrons and ions from the plasma as has previously been studied by us 42 . A key finding using a simple feedback model was that a glow discharge plasma is not stable unless the losses of plasma are dependent on the concentration of electrons to a power greater than one; in other words, the rate of loss of electrons is more than linearly dependent on the electron concentration.
The problem of stabilising wild oscillations is a common theme across disciplines that use predator-prey modelling, and a variety of approaches have been taken to lead to population stability. Nonlinear stability theory has been considered as a means of stabilising predator-prey interactions 37,43 , as has the introduction of a firm population limit imposed on one or both of the populations 44 . In keeping with our previous findings, here we include a nonlinear loss term as a natural way of inducing stability.
The model is expressed as coupled equations resembling the original Lotka-Volterra equations by writing the "prey" population as the electron concentration x 1 and the "predator" population as the particle material concentration x 2 . An additional term is added to include the nonlinear loss of the prey. Our resultant equations are as follows: Figure 1. An overview of the processes occurring in a dusty plasma. (a) Dust particles effectively 'eat' or consume plasma as a cumulative effect of the following processes that occur in a dusty plasma: (b) Electrons collide with particles of the precursor gas, such as silane, methane, or acetylene, splitting them into reactive fragments. (c) The resulting fragments coalesce to form dust particles. (d) Once a dust particle grows too large, it is lost from the system. (e) Plasma is lost from the system, through recombination of ions and electrons or loss to the environment of the plasma. Visualisation based on Hollenstein 12 .
Scientific RepoRts | 6:24040 | DOI: 10.1038/srep24040 where n > 1. The model parameters are listed in Table 1 together with their relation to the dusty plasma environment. Equation set [2] represents an extension of the classic Lotka-Volterra equations, which are recovered when k = 0. Note that the dust population is in our definition a quantity that reflects the concentration of dust material, defined as the mass of dust per unit volume. This is done because the main measure of dust concentration is by light scattering which is greater for larger particles as well as for larger concentrations of particles, preventing a convenient separation into particle number and particle size. The predator loss term γ indicates the loss of grown particles from the plasma, and it incorporates the gravitational force and the ion drag force, both directly proportional to the total mass of dust per unit volume, as in the definition of x 2.
The electron and ion loss exponent n reflects the order of the reactions that cause plasma loss. One of these reactions is recombination between an electron and an ion, according to the reaction e − + Ar + → Ar where an argon plasma is used as an example. This is a second-order reaction, since its rate depends on the populations of two species: electrons and ions. The population of electrons is approximately equal to the population of ions in normal laboratory plasmas, as a result of the overall charge neutrality, giving a loss rate of electrons proportional to x 1 2 , so for recombination, n is 2. Another process that causes loss is the loss of electrons to the chamber wall, which is a reaction of order 1, since it depends directly on the population of electrons in the chamber. Hence the resultant n falls between 1 and 2, since it reflects both loss mechanisms.
The modified Lotka-Volterra equations [2] have a class of solutions in which the populations of prey (electrons in a plasma) and predators (dust particles in a plasma) oscillate out of phase and with diminishing amplitude. This is the same behaviour as reported in measurements of electron and dust populations in dust-forming plasmas 19,[21][22][23] .
It is useful to plot the predator population as a function of the prey population, to reveal the basic behaviour of the solutions. For the relevant class of solutions, this plot spirals in to a point of stability, as shown in Fig. 2. In this section we show that solutions of the modified Lotka-Volterra equations give a great deal of insight into the behaviour of dust-forming plasmas. They contain within them an explanation of how the period of the oscillations varies with precursor gas concentration and with electrical power input to the plasma.
The choice of parameters can result in a steady state, a perpetual oscillation, a convergence to equilibrium, or the demise of one or the other, or both, of the species. Analysis of the modified Lotka-Volterra equations [2] allows an assessment of when a discharge might stabilise, if at all. Stability primarily depends on n and k, because the nonlinear loss is necessary for stability 42 . The modified Lotka-Volterra equations comprise a system of two nonlinear ordinary differential equations. Equilibrium points can be calculated for both the prey population (x 1 ) and the predator population (x 2 ) separately, by setting the opposite population's derivative to zero. Simple algebra yields: as the equilibrium points. Stability analysis requires linearising the modified Lotka-Volterra equations about the equilibrium points to produce the Jacobian, or community, matrix J 45 :

Term Quantity
Method of evaluation x 1 Electron density (equivalent to ion density) Resonance frequency of the chamber used as a microwave cavity n n n 1 2 2 2 2 1 Figure 3 shows eigenvalues for the solution to the modified Lotka-Volterra equations calculated for α = β = 1.5, and γ = δ = 4, as a function of n and k. This set of parameters allows a wide range of eigenvalues to be seen in a relatively small n-k parameter space. For much of the n-k parameter space shown, the eigenvalues are complex. However, where n = 1, the eigenvalues are purely imaginary, and where n and k are large, they are purely real.
The eigenvalues of the Jacobian indicate the stability behaviour of the system 45 . When the eigenvalues are purely imaginary, the dynamics of the system are on the cusp between stability and instability. This case is known as neutral stability, in which the system oscillates around equilibrium, without moving closer or farther away. The standard Lotka-Volterra equations show this kind of stability, which gives rise to perpetual oscillations of constant amplitude. Eigenvalues that are complex indicate oscillatory systems, where the predator-prey population plot is spiral in shape. Stability is determined by the sign of the real part of the complex eigenvalue: if it is negative, the system spirals in to an equilibrium point; if it is positive, the system spirals away from the equilibrium point. Population oscillations in time are damped for the stable spiral, while they increase without bound for the unstable spiral. Purely real eigenvalues that are negative indicate the presence of a stable node, where the system quickly equilibrates to a stable point without oscillating. Purely real positive eigenvalues yield an unstable node, where the system neither oscillates nor equilibrates. Figure 4 shows the stability properties of the system corresponding to the eigenvalues shown in Fig. 3. Much of the n-k parameter space shows stable behaviour, with stable spirals yielding to stable node behaviour where n and k are large. Small n and k produce unstable spiralling, and the special case of neutral stability is seen along most of the line where n = 1. It is interesting to note that the line of neutral stability does not go all the way to k = 0; rather, it turns upwards where k is small, resulting in a strip of parameter space where small k produces unstable spirals, regardless of n. It is at this point that the behaviour does not correspond to the eigenvalues' prediction.
It is important to observe that if either n or k is too small, the system is unstable and does not converge. In fact, convergence is only found where n > 1. In other words, the magnitude of the loss must not be too small, and the exponent of the loss term must be greater than 1. The latter result is in agreement with our previous findings on the relevance of the loss exponent n to plasma stability 42 and with a value of n > 1 reflecting the order of the reaction leading to electron loss. Comparison with observations. The phenomenon we seek to explain is manifested in several ways: the density of the plasma 19 , defined as the number of electrons in unit volume, a quantity equal to the ion charge density per unit volume because of plasma quasi-neutrality; the intensity of scattering of a laser beam directed through the plasma 19,22 and of the absorption of a laser beam 22,46 , both of which are dependent on the dust particle size and density; the intensity of spontaneous light emission from the plasma 21 , dependent on both the plasma density and the dust particle size and density; the unprocessed precursor gas partial pressure or concentration 22 where the rate of consumption is dependent on the concentration of precursor gas molecules; the mass-spectrometer signal of the active gas molecule 21 ; and the absorbance at a particular IR frequency 24 as well Figure 3. Eigenvalues of the Jacobian (eq. 5) as a function of n and k, as given by equation 6, with real and imaginary parts plotted on the same z axis. The data shown here were calculated with α = β = 1.5, and γ = δ = 4; the eigenvalue plot is different for each combination of α , β , γ , and δ . The imaginary parts of each pair of eigenvalues mirror each other across the n-k plane, which is indicated where the z axis is zero. The real part of the eigenvalues reaches a minimum where the imaginary parts reach zero, where n and k are both large. Beyond this line, the imaginary parts remain zero, and the real part remains negative but approaches zero. The behaviour plot is different for each combination of α , β , γ , and δ . This set of parameters allows a wide range of behaviours to be observed in a relatively small n-k parameter space. The stability of the system falls into four categories, based on the predator-prey population, or limit-cycle, graph behaviour. An 'unstable spiral' spirals away from the initial condition, yielding population oscillations that grow without bound. A 'neutral centre' indicates a closed loop, where population oscillations stay steady without damping. A 'stable spiral' spirals in to reach a point of stability, and population oscillations are damped until both populations remain steady. A 'stable node' indicates movement towards a stable point without oscillation, with both populations equilibrating quickly.
as IR absorption spectra that show the chemical signature of dust particles and precursors 20,24 . Table 1 shows the relation of the populations to the measured quantities.
Our model predicts trends in dusty plasma behaviour in good agreement with published observational data. The most direct measurements of the electron density x 1 , via resonance frequencies in a microwave cavity, have been shown to oscillate with time, and the presence or absence of damping in the oscillations has been shown to depend on plasma parameters 19 , consistent with our model's solutions in the stability behaviour classes of the 'neutral centre' and the 'stable spiral' .
An interesting feature of the solution for the predator population as a function of time is its tendency to approach toward zero during each oscillation, for some combinations of parameters, while the prey population remains significantly above zero. This scenario is consistent with the observation of Cavarroc et al. of periodic void (dust-free region) formation at the centre of a dust-forming plasma 26 . Figure 5a shows laser light intensity scattered and transmitted from an Ar/C 2 H 2 discharge, while it is producing dust particles. The intensity of the scattered laser light is a measure of x 2 , the number density and size of dust particles 22 , as set out in Table 1. Therefore, when the scattered signal oscillations die away, we can infer that the oscillations in dust population are also dying away, in good agreement with the model's 'stable spiral' solution, shown in Fig. 4 as well as in Fig. 2b. Figure 5b shows the luminous intensity l of an Ar/C 2 H 2 discharge (open circles), as well as the mass spectrometer signal of C 2 H 2 (filled circles). Both signals oscillate over tens of seconds, decreasing in amplitude over time. As set out in Table 1, light emission is proportional to the number of atoms and ions capturing or being excited by electrons and thus to the electron population x 1 . However, the presence of dust particles in the plasma has the effect of attenuating the optical signals used for characterisation 47 . Therefore, we define the luminous intensity l as w 1 x 1 − w 2 x 2 , where w 1 and w 2 are fitting parameters that determine the relative weighting of the electron density x 1 and the dust particle mass per unit volume x 2 .
The mass spectrometer signal arising from C 2 H 2 derived species indicates the concentration of C 2 H 2 molecules present in the plasma. These species are precursors to the carbonaceous particulates, or dust, which form in the plasma, as set out in Fig. 1, where C 2 H 2 molecules are dissociated by collision with electrons and then agglomerate to form more massive particles. Hence the C 2 H 2 derived signal reflects only the population of the dust particulate precursors, and the rate of consumption of C 2 H 2 is proportional to the electron population in the plasma. Therefore, the C 2 H 2 signal corresponds to the negative integrated electron population, such that the maxima in electron population align with the fastest decrease in C 2 H 2 population. In Fig. 5, the light emission and mass-spectrometer data shown in b can be compared with the calculations shown in c, which shows the corresponding quantities of l and the negative integral of the prey population curve, as calculated from equation pair [2]. The simulation is consistent with experiment, as both show damped oscillations that appear to approach steady state, and the calculated phase delay is such that the maxima in electron population align with the fastest decrease in C 2 H 2 concentration, as seen experimentally.
Oscillation frequency. The dust creation rate parameter δ is primarily determined by the precursor gas concentration, or, more specifically, the concentration of precursor elements. δ increases from a CH 4 plasma to a C 2 H 2 plasma, because C 2 H 2 has twice as much carbon per molecule and thus is likely to be more effective at making dust. One event of an electron encountering a carbonaceous particle produces twice as much carbon which can then be used to agglomerate to form larger dust particles.
Figures 5b and d show light emission and mass-spectrometer signals for the active gas in plasmas of two different compositions, Ar/C 2 H 2 and Ar/CH 4 , measured on the same system. Measurements of the Ar/C 2 H 2 plasma show much faster oscillation than those for the Ar/CH 4 plasma, as can be expected from the higher concentration of carbon atoms. A similar result has been reported experimentally, namely, that the period of oscillation of the electron density x 1 depends inversely on the fraction of the reactive gas in the gas mixture 19 . (No pressure dependence has been reported, other than a threshold pressure below which oscillations are not seen; above the threshold, the period of oscillation does not vary with pressure 19 ).
The calculated periodicity is affected by the dust creation rate parameter δ according to Fig. 6a, where a change in δ whilst holding the other parameters constant results in a change in period of oscillation in a quasi-parabolic shape. The data points shown here are calculated by averaging the first five periods of oscillation for each δ . Where the period decreases with increasing δ , the period of oscillation is relatively constant over the first five calculated periods of oscillation, as is seen experimentally. This left-hand side of the paraboloid can be thought of as the physically realistic side, because it corresponds to what is seen experimentally. The right-hand side of the paraboloid, where the period increases with increasing δ is effectively a mathematical artefact, and can be thought of as non-physical. On this side, for each δ , the period decreases considerably over the first five oscillations, by up to 50% of the duration of the first oscillation. The period of electron density fluctuation has been shown experimentally to depend inversely on the gas flow rate (Fig. 6c) as well as inversely on the silane fraction 19 , and it has been noted elsewhere that a lower flow rate of SiH 4 results in slower oscillations 46 . The gas flow rate influences δ , because it is the presence of the active gas that allows the creation of dust, according to the model described in Fig. 1. Hence the experimental curve of Fig. 6c corresponds to the left-hand side of the paraboloid in Fig. 6a.
The electron birth rate parameter α depends on both the rf power supplied to the plasma and the pressure, as these determine the rate at which electrons can create new electrons. The power determines how much energy the electrons are given, which in turn determines the rate at which electrons create new electrons. Pressure determines the mean free path and thus the frequency of collisions of an electron with a given energy. Both of these effects increase the electron concentration x 1 , which is reflected in an increased rate of growth of the particles. The increased rate of growth increases the frequency of the oscillations. Figure 6b shows the effect of α on the period of oscillation, as calculated from equation pair [2]. For comparison, experimental data show that rf power affects the period of electron density fluctuation as shown in Fig. 6d, where period decreases with increasing rf power until it reaches a minimum and then increases with increasing rf power. The calculated inverse relation of b corresponds to the initial decrease seen experimentally in d.
The predation rate β corresponds to the particles' affinity for electrons, which would vary for different species of predator. For example, a silicon-based predator such as silane would have a different predation rate of electrons from a carbon-based predator such as methane or acetylene. A particularly interesting case to consider is that of candle flames, where the particles (predators) create electrons (prey) directly, via thermionic emission from the soot. It is a further twist on the classic Lotka-Volterra relation, a complement to the situation set out by Koren and Feingold in the context of aerosol-cloud-precipitation systems, where the prey (cloud) spawns the predator (rain), as raindrops coalesce from cloud droplets 41 . In a candle flame where the electron emission is stronger than the electron predation because of thermionic emission, then β would have a negative sign, indicating creating rather than destroying. In this case, it is still possible to attain stability, as n and k determine the death rate of the electrons, rather than β . In fact, the eigenvalues as set out in equation [6] happen not to include β , so its sign or indeed even its magnitude does not appear to affect the stability behaviour at all.
The essence of the predator-prey problem, namely, that a population does something in their own interest, short-term, that has the longer-term effect of being detrimental to them, is a phenomenon that appears in a wide variety of fields. It is manifested in this work by the formation of predator-free regions (voids) in dusty plasmas. We have demonstrated a natural way of inducing stability in the predator-prey problem by nonlinear losses, which we commend for wider application. The phenomenon of the electron (prey) population reaching its maximum just before the particle (predator) population is elegantly captured by the predator-prey model in a way no other model does. The success of ecological population dynamics in describing the population interactions of dust-forming plasmas may serve as an inspiration for others who seek insight via interdisciplinary integration.