# Josh Testing 1

###### This wiki is incomplete.

PS: Here are examples of great wiki pages — Fractions, Chain Rule, and Vieta Root Jumping

x=2

this is inline \(x=2\) and we

\[\overbrace{1+1}^\textrm{thing}\]

the box is \(\SI{5}{\watt\per\square\meter}\) long

\[\ce{c6h12o6 + o2 -> co2 + h2o}\]

\[x=2\]

Gravitational waves are all the rage these days, signature of events like binary stars, mergers, supernovae, etc. In general a very big signal dominated by the source's bulk distribution of matter and a much smaller time-dependent component reflecting the local motion.

One of the defining features of Newtonian gravity is that it adjusts instantaneously everywhere in space (there is no delay in field strength). We can write the field \(\varphi\) in terms of the distribution of matter \(\rho\) throughout space.

\[\varphi = -G\int\frac{\rho{r, t}}{\lvert r\rvert} d^3 r\]

Over a small region, a wave is a variation in the field like the one below (a wiggle, not a uniform rise). Is a field profile like this possible under Newtonian gravity?

These wiggles are impossible under NG because every source contributes a term proportional to 1/r, and no wiggles are possible from adding a bunch of monotonically decreasing curves. Suppose it took finite time for the gravitational field strength at a distance to reflect the arrangement of mass at the source. Could any of the arrangements below lead to wiggles in the field?

<sphere of mass rotating, binary star spinning, black hole merger, a star moving uniformly away>

The crucial ingredients for generating gravitational waves are thus

- a time-varying distribution of mass at the source and
- finite propagation speed of gravitational field strength.

The second of these is a natural outcome of Einstein's general relativity, but there is nothing stopping us from modifying Newton’s formulation on the fly.

Which of the following expressions reflects the idea that the field takes a finite time to propagate through space?

\[\varphi = -G\int\frac{\rho{r, t - r/c}}{\lvert r\rvert} d^3 r\] \[\varphi = -G\int\frac{\rho{r + ct, t - r/c}}{\lvert r\rvert} d^3 r\] \[\varphi = -G\int\frac{\rho{r - ct, t}}{\lvert r\rvert} d^3 r\]

Thus, in concept, gravitational waves are fairly simple. The gravitational field is equal to the expression from the last question and it is our task to find its time dependent components. We’ll take the approach of treating these wiggles as small perturbations on top of the background field. We can do this because the length scale of the variations (roughly the size of the source) are much smaller than the source’s distance from us, in other words \(\langle r\rangle \gg \ell_\textrm{source}.\)

We can see this from an expansion of \(r^{-1},\)

\[r^{-1} = \lvert x\rvert^{-1} + 2 \frac{\mathbf{x}}{\lvert x\rvert}\cdot\mathbf{y}\frac{1}{\lvert x\rvert} + O\left(\lvert x\rvert^{-2}\right) \]

Since \(\lvert x\rvert \gg \lvert y \rvert\) small variations in the source simply do not contribute substantial terms to \(r^{-1}\) and we can ignore any terms higher order in \(\lvert x\rvert{-1}.\)

Under this approximation, all of the variation is captured in the integral of \(\rho\) — but this expression isn’t very transparent. The source may have complex motion, and the closed form integral doesn’t tell us how significant ripples due to source motions will be relative to the background field.

Put simply, the source is a distribution located roughly \(\lvert\mathbf{x}\rvert\) away and contributions to gravitational field strength are delayed, \(\Delta t = \lvert r\rvert/c\) so that they arrive at time \(t - \Delta t.\) The delay is slightly different for each point \(\mathbf{y}\) and we can write \(\Delta t = \lvert r\rvert/c = \overbrace{\lvert x \rvert/c}^{\Delta t_0} + \overbrace{\left(\lvert x-y\rvert - \lvert x\rvert\right)/c}^{\Delta t_y}\)

Because \(\Delta t_y\) is small relative to \(\Delta t_0,\) we can expand \(\rho\) in \(\Delta t_y.\) Writing \(t_0 = t - \Delta t_0,\) we have

\[\begin{align} \varphi &= -\frac{G}{\lvert x\rvert}\int\rho(r, t_0 - \Delta t_y) d^3 y \\ &= -\frac{G}{\lvert x\rvert}\int \left(\rho(\mathbf{y}, t_0) - \dot{\rho}(\mathbf{y}, t_0)\Delta t_y + \frac12 \ddot{\rho}(\mathbf{y}, t_0) \left(\Delta t_y\right)^2\right) d^3\mathbf{y} \end{align}\]

This expansion will simplify our task immensely.

Find \(\Delta t_y\) in terms of \(\mathbf{x}\) and \(\mathbf{y}.\)

The terms in this expansion might look foreboding, but the first two are quite simple.

The first, \(\int \rho\, d^3\mathbf{y}\) simply yields the source’s total mass, \(M.\) Stopping here, we’d have \(\varphi = -\frac{GM}{x},\) which is obviously independent of time.

The second, \(\int \dot{\rho} \left(y_1 + y_2 + y_3\right)\, d^3\mathbf{y},\) is a bit tricky. By the conservation of mass, the source density has to follow the continuity equation. E.g. for any volume, the change in mass contained per unit time \(\dot{\rho}\) is equal to the rate at which matter flows across its bounding surface. \(\)

\(\)

Thus we have \(\dot{\rho} = \frac{\partial}{\partial y_1}\left(\rho v_1\right) + \frac{\partial}{\partial y_2}\left(\rho v_2\right) + \frac{\partial}{\partial y_3}\left(\rho v_3\right).\) \(\)

\(\)

With this insight, the term becomes a sum of three integrals of the form \[\int \left(y_1 + y_2 + y_3\right) \frac{\partial}{\partial y_i}\left(\rho v_i\right)\, dy_1\, dy_2\, dy_3\] (for \(i = 1\,2\,3\)).

What does this equal for \(i=1\)?

**Hint**: integrate by parts!

This shows that the second term is nothing but \[\int\rho\left(V_1 + V_2 + V_3\right)\, d^3\mathbf{y} = M\mathbf{V},\] the total momentum of the source. Sadly, this is conserved and does not change with time.

But what’s that, the faint beat of a drum? That’s right, we’ve arrived at the third term in the expansion and the party is about to begin.

Suppose we decide to throw caution to the wind, and calculate the gravitational field of a spherical mass shell by hand, adding the points one by one.

We can start small and approximate the sphere by the cube of points (black points in the diagram) that intersect the surface \(\mathcal{C} = \{\pm1/\sqrt{3}, \pm1/\sqrt{3}, \pm1/\sqrt{3}\}.\) We could make things a bit easier on ourselves by calculating the field at a point that follows a symmetry of the cube, e.g. directly above the center of a face. But we're not interested in abiding symmetry of our simplest approximation of a sphere — we're interested in the field of a sphere!

So let's consider a general point like the yellow point in the diagram, \(\mathbf{r}=\hat{\mathbf{x}} + \hat{\mathbf{y}}+ 3\hat{\mathbf{z}}.\) What is the field at \(\mathbf{r}\) due to \(\mathcal{C}?\) As we noted in the introductory pane, the field at a point \(\mathbf{r}\) is given by \(\sum_i \mathbf{g}\left(\mathbf{d}_i - \mathbf{r}\right)\) where the \(\mathbf{g}_i\) are given by

\[\mathbf{g}_i = G\frac{m_i}{\left(\mathbf{r} - \mathbf{d}_i \right)}\Delta \hat{\mathbf{d}}_i\]

and the unit vector \(\Delta \hat{\mathbf{d}}_i\) is given by

\[\Delta \hat{\mathbf{d}}_i = \frac{\left(r_x-d_x^i\right)\hat{\mathbf{x}} + \left(r_y-d_y^i\right)\hat{\mathbf{y}} + \left(r_z-d_z^i\right)\hat{\mathbf{z}}}{\sqrt{\left(r_x-d_x^i\right)^2 + \left(r_y-d_y^i\right)^2 + \left(r_z-d_z^i\right)^2}}.\]

There's no way around this, we just have to add up the terms and pray for a pattern. As it turns out, there is none!

The sum over contributions from all the points is given by (brace yourself):

Each term has a unique denominator, so the terms pile up and we're left with a rather opaque, unsimplifiable expression for the gravitational field that leaves even the most seasoned physicist crying tears of blood.

The brute force analytic approach grows increasingly unmanageable as the number of particles increases, and if we lose symmetry it will be nearly impossible. In the real world, we can't afford to avoid computation.

As the recent fires in California have reminded us, forest fires can devastate communities in hours or days, propagating more quickly than we can put them out. But they are a fact of life on a planet covered in trees.

By understanding the physics of forest fires, we can help form policies that make them more manageable. It's an important application of computational science and the conclusions affect millions of households in the U.S. alone.

In this Left to the Reader, we will build up to a basic understanding of the factors that govern the life cycle of a forest fire, including the **structure of the landscape** and the **physics of the spread**.

Here's how we'll get there:

- We'll learn about some basic features of forest fires, look at data, and observe some simple behaviors.
- We'll design a minimal model that captures the basic structural and dynamic features of wildfire.
- We'll apply a technique from statistical physics, called mean-field theory, to make mathematical predictions about the model.
- We'll employ a computer simulation to evaluate the predictions of our model and see what it has to say about public policy regarding wildfire.

#### Contents

## Details and mechanisms

To take a stab at modeling a growing wildfire, let's get familiar with the details and mechanisms of forest fires.

As forests age and hydration fluctuates, trees and other foliage go through cycles of dryness that make them more or less susceptible to incineration. This leads to seasonal and longer term windows where forests are susceptible, characterized by a combination of drought, peak solar radiation, low humidity, and warm ambient temperature, and can be triggered by random events like lightning strikes.

There are also dynamical factors. For example, wind speed varies in space and time, and changes the pace of fire spreading and bias its course toward more or less flammable regions. Whether a fire turns toward a field of old pine or into a rocky ravine with sparsely planted, young spruce can be the difference between a fire's end and its beginning anew. Harder to predict, a spark or ember can catch a gust of wind, and transport very long distances and start a new fire in a distant part of the forest.

Knowing the particular details of a given forest would be great, and might allow us to make more accurate models of that forest — but it is hard to collect timely, forest-wide information for one forest, let alone several. Before getting too far into the weeds we should ask: how important are the details?

## Looking at the data

To get an idea about the severity of forest fires, we can look at some data. One might think that forest fires are, in a sense, the balance between incendiary factors like fuel, heat, etc. and the things that stand in the way of combustion like hydration and rocks. If forest fires were strongly dependesee data like the following, where there is an average behavior and the particular features of a forest make it more or less severe.

Instead what is found are patterns like those shown below. Notice that these distributions don't behave in the ways we're accustomed with the normal or a Poisson distribution, where there is a clear average around which the population varies — the burn areas of span well across seven orders of magnitude! Nor do they behave like the exponential distribution, where events much bigger than the average are unobservably rare.

In fact, these probability distributions show the hallmark of a remarkable property known as **scale-invariance** — let's dive into what that means.
Because the data is plotted on logarithmic axes and exhibits a linear relationship, we can write down the math straight from the graph — the \(\log\) of the frequency, \(\log p,\) and the \(\log\) of the area, \(\log A,\) are related by a line, or

\[\log p = \alpha - \beta\log A,\]

where \(\alpha\) is the \(y-\)intercept of the data and \(-\beta\) is its slope. If we exponentiate both sides, this relationship becomes \(p\sim A^{-\beta}.\) In other words, the probability of observing a forest fire of size \(A\) is proportional (the \(\sim\) symbol) to \(A^{-\beta}.\)

Approximate

Scale-invarianceLooking at the probability distribution of forest fire areas, which has the form \(p(A) \sim A^{-\beta},\) we can see that the forest fires that burn out the greatest area of forest floor, \(A_\textrm{max},\) are the most rarely observed. What if we're interested in the likelihood of fires of area \(A < A_\textrm{max}?\)

Given the simple form for \(p(A)\) we can see that \[\boxed{p(A) = \displaystyle p(A_\textrm{max})\times\left(\frac{A}{A_\textrm{max}}\right)^{\beta}}.\]

When this property holds over all values of forest fire area, it is known as

scale-invariance. Though forest fires may not be truly scale-invariant (see Grassberger 2002), the forest fire data appears to show scale-invariant behavior over a wide range of forest fire areas.The revelation is this: scale-invariance allows us to estimate the frequency of the largest (rarest and hardest to observe) forest fires in terms of the more easily measured frequencies of small ones.

Suppose that the measured frequency of forest fires that burn down \(\SI{100}{\kilo\meter\squared}\) is found to be \(\SI{0.5}{fires\ per\ year},\) and that \(\beta\) is known to be approximately \(1.2.\)

What do you expect the frequency of forest fires of burn area \(\SI{10000}{\kilo\meter\squared}\) to be?

If this widely observed behavior were dependent on the fine details of any particular forest, then it would be quite a conspiracy between the global growth of trees, landscapes, wind patterns, seed dispersal, etc. What's more likely is that fundamental features, shared by all forests, are what produce the patterns we see (in fact that's one of the implications of scale-invariance).

Scale-invariance raises a troubling question: if forest fires exhibit an indifference to details, then how can we anticipate the severity of a given forest fire? As we alluded above (link to details and mechanisms), the prediction of individual forest fires is subject to all sorts of complexities, model chaos, and in general is not a very reliable tool. Thus, while a holistic understanding of any particular forest gives us more information, it is not necessarily useful in understanding the observable patterns in wildfire severity.

Happily, there are less particular, structural features that we can use to characterize forests, and use to model the course of nascent forest fires. We're going to design a general model to interrogate forests in terms of their basic structure (tree density) and lightning frequency. In the process, we'll distill the complexity of factors like prevailing winds, geography, and other pesky factors and extract general lessons about how all forest fires behave.

## Modeling parameters

For the rest of our time here, we're going to model forest fires using the model outlined below.

Model rulesWe can model the forest as a 2D lattice of dimensions \(L\times L.\) The spatial density of trees \(\rho_\textrm{tree},\) measured in trees per unit area, captures the concentration of trees per unit area of the forest floor. We therefore also have a density of empty sites \(\rho_\textrm{empty} = 1 - \rho_\textrm{tree}.\) At time zero the forest contains no trees and the lattice sites are updated according to the rules below.

At each time step,

- Trees are replenished at empty lattice sites with probability \(r.\)
- Trees are lit on fire (e.g. by lightning) with probability \(f.\)
- If any of a non-burning tree's neighbors are on fire, then the tree is set on fire.

- In the setup shown above, the trees to the left and top of the flaming site would be on fire in the next time step.
- In addition, the trees in the four corners each have probability \(f\) of incinerating due to lightning strike.
- The two gray sites each have probability \(r\) of containing a tree in the next time step.

**Is that all?**

It might be surprising that a model meant to describe phenomena as rich and variegated as the range of forest fires can be stated so simply. It is important to stop here and ask whether we've missed something important.

On the surface, this model seems to contain no information about forest aging, hydration, or variations in landscape and it certainly does not incorporate prevailing winds, the succession of different plant species, or climate effects.
In fact, the model actually does incorporate the first three effects, and the contribution of the last three **do not appear necessary to explain the data**.

## Basic Model

Before we simulate a universe of forests, we can do some simple mathematics to see the basic relationships at play.

For example, how does the frequency of lightning strikes \(f\) affect the size of the average forest fire?

Let's introduce one last quantity—we'll call the number of trees that burn down in a lightning strike \(s\) and the average number of trees that burn down in any given lightning strike \(\langle s\rangle = \sum_s p(s)\times s.\) If we calculate this on the lattice, taking into account the possible correlations between neighboring sites and the exact dynamic of the lattice, we'd be swiftly overwhelmed by the complexities, at least as rookie forest fire modelers.

To facilitate this calculation, we're going to use the idea of an ensemble of forests. This is an idea from statistical physics where we imagine there to be a large number of forests over which we take an average. In the ensemble, we can forget about the behavior of any one system and instead focus on how things work on average, across all of the systems. This is far easier than the detailed work of finding the trajectory of one forest in particular.

Concretely, the average forest has tree density \(\bar{\rho}_\textrm{tree},\) and sees \(\langle s\rangle\) trees burn down per lightning strike, on average.

We can relate these quantities to each other by conservation equations. For example, because the average forest has tree density \(\bar{\rho}_\textrm{tree}\) (which we'll write as \(\rho_\textrm{tree}\) from now on) it cannot be losing or gaining any trees on average.

Therefore we can say

\[\langle\textrm{new trees planted}_t\rangle - \langle\textrm{trees burned down by lightning}_t\rangle = 0.\]

In other words, the number of trees that disappear due to fire in any time step \(t\) must be equal to the number of trees that are planted to replace them in that same time step, on average.

**The number of trees planted per time step**

This style of calculation is called **mean-field** theory, because we are assuming that we can ignore spatial variations and consider the interaction of the average values with themselves.
Since this is a strange and fascinating tool, we'll calculate the first part together.

Suppose we find the lattice in the state above and that \(p = 0.6\bar{6}.\) On the left we can see three empty sites that each have the potential to be filled by planting events. How many of these to we expect to be filled on average?

The answer is \(3\times0.6\bar{6} = 2.\) We mark the newly birthed trees on the right with asterisks to showcase their newcomer status.Looking case by case is a lot of fun, but it isn't going to get us anywhere. Instead, we can write the expected number of newly birthed trees as the density of empty sites \(\rho_\textrm{empty},\) times the probability of planting per lattice point per unit time \(p,\) times the number of lattice sites \(L^2.\)

Thus,

\[\boxed{\sum_\textrm{forests}\textrm{new trees planted}_t\times P(\textrm{forest}) = \rho_\textrm{empty}\times p\times L^2}.\]

With this under your belt, finding the number of trees that burn down per time step should be a breeze.

What is the average number of trees that are consumed by flames per unit time?

**Recall** that new fires are initiated by lightning striking a random site on the lattice; if it's occupied by a tree then it will ignite; if it's empty then nothing will happen.
The average fire rages until \(\langle s\rangle\) trees are burnt down.

With these two quantities in hand, we can equate them to find \(\langle s\rangle\) as a function of tree density and the frequency of lightning strikes.

Use the conservation equation to find \(\langle s\rangle\) in terms of \(\rho_\textrm{tree},\) \(f,\) and \(r.\)

What is \(\langle s\rangle\) when \(\rho_\textrm{tree} = 0.4,\) \(f = 0.01,\) and \(r = 0.05?\)

**Hint:** Recall that \(\rho_\textrm{empty} + \rho_\textrm{tree} + \rho_\textrm{burning} = 1.\)

**Assume** that \(\rho_\textrm{burning} \ll \{\rho_\textrm{empty}, \rho_\textrm{tree}\}.\)

After some algebra, and recalling that \(\rho_\textrm{empty} = 1 - \rho_\textrm{tree}\) we find\[\boxed{\displaystyle\langle s \rangle = \frac{p}{f}\frac{1-\rho_\textrm{tree}}{\rho_\textrm{tree}} \propto \frac{p}{f}}\]

## The tradeoff

The last result suggests an intriguing relationship between the number of trees burned down in the typical forest fire, and the frequency of lightning strikes that start new fires.
Naively, we might associate a great prevalence of fire with more fearsome infernos, but our **equation of state suggests the actual relationship is more subtle**.

Given the result we just obtained, what can we say about the relationship between the rate of fire setting events \(f\) and the size the average forest fire \(\langle s\rangle,\) all else being equal?

To really interrogate this question, we'll have to move beyond mathematical analysis and simulate the forest environment with code.

## Simulating the forest

Luckily, we can make some minor modifications to code we saw in a recent **Problem of the Week**.

In that problem we employed the same model we use here, with the exception that we looked at the course of fires resulting from single lightning strikes upon preformed forests. Averaging over thousands of instances of the forest, we found an intriguing relationship between tree density \(\rho_\textrm{tree}\) and the duration of the forest fire \(T_\textrm{burn},\) shown to the right. Counterintuitive though it may be, there is a simple explanation

- at low tree densities the forest is not connected enough for fires to spread far beyond their starting point
- at very high densities the forest is fully connected and fires spread at their fastest rate
- it is at the critical density \(\rho_c\approx 0.6\) where the forest first becomes dense enough to sustain large scale fires, but still has hard to access pockets, that fires burn the longest.

In statistical physics, this is known as a phase transition as we can think of the value \(\rho_c\) as separating two qualitatively different kinds of forest: sparsely connected forests that don't suffer large-scale wildfires, and thick forests that essentially burn to completion whenever lightning strikes. This result reflects a deep connection to the problems of percolation theory, the focus of a future Left to the Reader.

## Modeling the intervention strategy

Now that we've established our model, we return to our original question.

Wildfires are a problem society must face as long as people insist on living near forests — so how best do we manage them? The tradeoff we found between the size of the average fire and the frequency of lightning strikes, \(\langle s\rangle \sim 1/f\) presents an interesting question to those interested in controlling wildfires: what, if any, effect do common fire control tactics have upon the severity of wildfires?

We can simulate the impact of fire control policies by varying the value of \(f\) from low to high. To see this, realize that the impact of a lightning strike is to burn down some fraction of the forest. If we make it a policy to try to put fires out as soon as they start burning, so that they don't burn down any significant area of the forest, it is as if we've simply eliminated some fraction of the fires, so that \(f\) becomes \(f^\prime < f.\)

We'll call the high-\(f\) limit the "

natural regime" and the low-\(f\) limit the "Smokey-the-Bear regime", for the famous mascot of fire control who encouraged every one of us to put out even the smallest forest fire, lest it transform into a raging inferno.

To see how changing \(f\) affects the nature of forest fires, we can use the codex below where we have programmed the model we developed above. The simulation starts off with an empty of \(L^2\) sites, after which we apply the following routine for 4000 timesteps:

- Each site in the lattice is tested to see if it's empty. If it's empty, then we try to plant a tree there with probability \(r.\)
- If it isn't empty, we check to see whether any of its neighbors are on fire in which case we set it on fire.
- If none of its neighbors are on fire, but it is a tree, we allow lightning to set it on fire with probability \(f.\)
- If the tree is currently burning, we set the lattice site to \(\textrm{Empty}.\)

We can be as sophisticated as we want to in measuring the dynamics of forest fire, but a quick way to gain insight is to simply keep track of how many trees are burning as a function of time (inspect for yourself to see how the number of burning trees is kept track of in the code). Our lattice is set to \(L=25\) and \(r\) is set to \(0.02\) which means that in the absence of fire, it should take \(T\approx 1/r = 50\) time steps to plant the entire forest. We initially set \(f = 0.01\) so that the probability of lightning striking a lattice point is half the probability of a tree being planted at an empty site.

How does the nature of forest fire vary as you lower the incidence of lightning strikes?

Noteyou're encouraged to play with all the parameters (\(r,\) \(f,\) and \(L\)) but the codex has a runtime limit of \(\SI{10}{\second},\) so keep that in mind as you explore the forest.

[[codex-lttr-fire-time-series]]

For the default parameters of the system, we observe the time series shown in the top row of the figure to the right. The number of trees hovers around \(\langle s\rangle = 200\) so that \(\hat{\rho}_\textrm{tree}\approx s/L^2 = 200/625 = 0.32.\)

We can see that from moment to moment, the number of trees fluctuates, but the fluctuations are small relative to the mean number of trees. This implies that the number of trees that burn down per unit time is small relative to the number of trees. In other words, although there are many fires they are each relatively inconsequential and short lived.

As we lower \(f\) we notice several things.

- The average number of trees in the forest increases, rising up to nearly \(\langle s\rangle \approx 400\) when \(f = 0.00001.\)
- The fluctuations in the number of trees increases tremendously. When \(f = 0.01,\) the fires are not big enough to meaningfully affect the state of the forest, as we noted above. As \(f\) drops, the number of fires is reduced, but the swing in tree count explodes. This suggests that the severity of forest fires is increasing.
By the time \(f = 0.0001,\) the flucutations are almost equal to the size of the system \(\sqrt{\langle T^2\rangle - \langle T\rangle^2} \approx L^2.\) At this point the fires are so severe that the typical fire can essentially burn down the entire forest!

Notethat the drops in these plots slightlyunderestimatethe severity of fire, as new trees can start to grow in the fire's wake even as it continues to burn across the system.

To get a quantitative measure for the severity of fire, we should measure the average number of trees burned down in a fire, \(\langle s\rangle.\) It is not trivial to do this directly, as tracking the provenance of a single tree's burning down requires the use of recursive search algorithms over the lattice. Instead we can use the connection described below.

Note, that the average number of trees that grow between lightning strikes, \(\langle s\rangle_\textrm{strike to strike},\) is equal to the rate of tree growth divided by the rate of lightning strikes.

Thus we have

\[\begin{align} \langle s\rangle_\textrm{strike to strike} &= \frac{\left(1-\rho_\textrm{tree}\right)pL^2}{\rho_\textrm{tree}fL^2} \\ &= \frac{p}{f}\frac{1-\rho_\textrm{tree}}{\rho_\textrm{tree}} \end{align}\]

However, the right side is the expression for \(\langle s\rangle\) that we found before. Thus, we can

measure the size of the average clusterof trees burned down by a single lightning strike bycounting the number of trees that grow between strikes.

Measuring \(\langle s\rangle\) as a function of \(f\) is important, but we'll leave this for the first of the **Discussion Challenges** at the end of the post.
In the meantime, we can already explain a puzzling and painful lesson learned by the US Forest Service over a century of trying to contain wildfires.

## The Yellowstone Effect

Up until the early 1970s, the US Forest Service maintained a policy to suppress every fire discovered in the forests it managed within within 24 hrs of their discovery. This zeal for fireless forests knew few bounds, perhaps peaking with the invention of so-called "smokejumpers", firefighters who would parachute out of airplanes to fight fires located in hard-to-reach regions of the forest. The policy was well-intentioned, influenced by a string of devastating fires in the late 1800's that tore through midwest territories and resulted in significant casualties.

But based on our studies here, we can see a fatal flaw in the thinking of the US Forest Service. By stopping most fires dead in their tracks, before they have a chance to burn through any significant portion of the forest, the rangers were effectively setting \(f\) below its natural value. As we saw above, this can be expected to provide**short term stability at the expense of catastrophic, system wide forest fires**.

In June of 1988 this came to a head when lightning strikes set off a handful of small fires in the park. The service had recently come to appreciate the utility of natural fires and allowing small fires to burn out, and therefore decided to let them burn themselves out.

However, the forest was already tuned to \(f^\prime < f\) and was in a state ripe for inferno.
By the end of July, the fires had not burned themselves out and had destroyed \(\num{100000}\) acres of forest.
At this point, it proved too late to implement fire controls and the fires raged on, eventually consuming roughly \(\num{800000}\) acres of Yellowstone, or \(36\%\) of the park.
In executing the no-fire policy the service was myopic, and **tuned the forest to a state destined for a system-scale inferno**.

## What did we learn?

We've covered a lot of ground here, before getting into the discussions let's take a step back and recall the path that got us here. We started out by looking at data from real forest fires where we noticed that the probability of forest fires of different sizes is scale-invariant. This implies that the fine details of individual forests are not responsible for the global behavior. Inspired by this, we designed a simple, universal model that considered only two explicit characteristics of the forest: the rate of tree replenishment \(r\) and the rate of lightning strikes \(f.\)

We then used mean-field theory to calculate a simple prediction for the size of forest fires as a function of lightning frequency (\(\langle s\rangle \sim 1/f\)) and verified this by explicit simulation in the codex. Simple though it is, this model helped us to understand the devastating 1988 Yellowstone fires and the contribution of forest control policies to the severity of wildfire.

In the end, forest fires are a natural clearinghouse for old trees and flammable materials that work to keep the forest in a healthy state. Attempts to naively constrain them can bottle up the natural volatility of the forest, causing it to erupt violently all at once.

In the **Discussion Challenges** below we'll build on top of these lessons and see if we can take things to the next level.

## Discussion Challenges

**Discussion #1**

Modify the code we used to measure the number of trees as a function of time to measure \(\langle s\rangle.\) Remember, this is difficult to measure directly, and it's much easier to use the mathematical connection we derived in the post.

**Discussion #2**

Our simulation shows that a zero-tolerance policy increases the likelihood of large-scale fires, but it doesn't suggest that it's impossible to do fire control in a smarter way. Can you design a fire suppression policy that increases the density of trees in the forest while reducing the magnitude of the average forest fire? Can you show that it works using the simulation?

**Discussion #3**

Moving beyond measurement of the average fire size \(\langle s\rangle\) is non-trivial, as we mentioned in the post. Can you design an algorithm that measures the size of the fire that results from each individual lightning strike? How does your estimate of fire size counts \(N(s)\) change with \(s?\) Note that to do this with statistical significance, and to take the measurement into larger systems, you'll likely need to run your code offline.

Elementary problems in physics often yield smooth curves as trajectories. [GIF of particle doing a smooth motion on a surface, another of a fired cannonball.]

Smooth solutions are commonplace, though they aren’t dictated by the laws of physics. The same rules that describe the flight of cannonballs, or the path of a bicycle can accommodate random motions as well. The determining factor is the form of the forces or, equivalently, the potential energies that govern the system.

Concretely, in the case of a particle near the surface of Earth, we have \(a = g\) which leads to the smooth velocity \(v = gt.\) Alternatively, for a car with smooth shifts in horsepower, we have the smooth velocity \(v = \int a(t) dt.\) In each case, a smooth force gives rise to smooth velocities and positions.

This is fine for a large slice of common dynamics like the orbits of the planets or the flight of a rocket, but falls flat for other common motions like the diffusion of a gas, the folding of proteins, or the formation of river morphology. Such problems generate trajectories that are continuous, but not smooth. Smooth dynamics are also ineffective for applying the methods of physics to non-physical problems like the stock market, or traffic systems.

# Stochastic dynamics

The fundamental distinction is that unlike the reliable strength of the Earth’s gravitational field, these dynamics are dominated by **stochastic** events like

- the thermal jostling of biomolecules by their fluid environment
- the unpredictable confluence of random timing and decision making in the stock market
- or the interaction of stream water with the random surface structure of rock

This means we have to find ways to represent uncertain events and incorporate them into the laws of motion. In this post, we’ll explore stochastic motion through the particular case of the Brownian particle—a small molecule that experiences random collisions with the fluid it's in.

At the macroscopic level, we might liken this to the position of a soccer ball as its passed around on the field. Of course it has some non-random features, like the fact that people are trying to put it through the goal, but the timing and directionality of the passes and shots are largely unpredictable and up to in-the-moment decisions by the players.

Without writing down a particular model, we know that the soccer ball's motion has several prominent features:

- Its current position and velocity depend only on its previous position and velocity.
- The timing and directionality of kicks are unpredictable (from the modeler’s perspective).
- Its trajectory is continuous though its velocity isn’t.

Let’s try to capture some of these features with a simple model in one dimension.

# One dimensional model—the random walk

To get started, we'll first consider a model in which the position of the soccer ball is confined to move along one spatial dimension. In this initial attempt at describing random motion, we’ll relax some of the features we listed above and employ a simpler set of assumptions:

- We discretize space and time so that the ball can reside at positions \(r_t \in \mathbb{Z},\) and its position is updated at discrete times \(t \in \mathbb{N}.\)
- We replace velocity with an
**update rule**: at each moment the particle can make a move \(m_i=\pm1\) that shifts one unit up or down the field with equal probability.

We can therefore represent the trajectory, \(\gamma,\) of the ⚽️ by the sequence of moves it takes from its starting position: \[\gamma = \{m_i\} = \{\uparrow, \uparrow, \downarrow, \uparrow, \downarrow, \downarrow, \downarrow,\downarrow\}.\] In this case, the ball ends up two units downfield of its starting position after eight time steps, i.e. \(r_8 \sum_i m_i = -2.\) For this trajectory, our ball ends up at \(r_8 = -2,\) but where do the rest of the balls go, and how likely are they to end up there?

Let’s enumerate all the two step trajectories \[\begin{array}{c|c} \text{Sequence} & \sum_i m_i \\ \hline \{\uparrow, \uparrow, \uparrow\} & 3 \\ \{\uparrow, \downarrow, \uparrow\} & 1 \\ \{\downarrow, \uparrow, \uparrow\} & 1 \\ \{\downarrow, \downarrow, \uparrow\} & -1 \\ \{\uparrow, \uparrow, \downarrow\} & 1 \\ \{\uparrow, \downarrow, \downarrow\} & -1 \\ \{\downarrow, \uparrow, \downarrow\} & -1 \\ \{\downarrow, \downarrow, \downarrow\} & -3 \\ \end{array}\]

We can see that the ball has one way to end up 3 positions upfield, three ways to end up 1 positions upfield, three ways to end up 1 position downfield, and one way to end up three positions downfield. We can average the displacement \(\sum_i m_i\) over all the trajectories and we find

\[\begin{align} \langle r_3 \rangle &= \langle\sum_i m_i\rangle \\ &= \langle m_1\rangle + \langle m_2\rangle + \langle m_3\rangle \\ &= 3\times\frac{1-1}{2} \\ &= 0 \end{align}\]

This is a little unsatisfying as it suggests that, on average, our ball goes nowhere when we know that the average ball does (in fact all of them do) move away from the origin. It suggests that a straight sum loses important information about our trajectories.

A direct approach is to average the absolute values of the displacements, i.e. to use \(\langle r_3\rangle = \lvert \sum_i m_i\rvert = \frac{3 + 1 + 1 + 1 + 1 + 1 + 1 + 3}{8} = 1.5.\) This is great, because it effectively captures the ball’s displacement, but in practice absolute values are awkward to work with and make analysis difficult.

## Random walk - calculation

Instead we can use a more tractable measure of displacement, the square of the sum \(\left(\sum_i m_i\right)^2.\) Writing this out for \(r_3\), we have

\[\begin{align} \langle r_3\rangle &= \left(m_1 + m_2 + m_3\right)^2 \\ &= m_1^2 + m_2^2 + m_3^2 + 2\left(m_1m_2 + m_2m_3 + m_3m_1\right) \end{align}\]

It might look like things just got more complicated since we have 9 terms that we have to average over all 8 possible trajectories. The squared terms \(m_i^2\) are just 1 since \(m_i = \pm 1.\)

For the cross terms like \(m_1m_2,\) we can average over all possibilities: \(\langle m_1m_2\rangle = \frac14\left(1\times1 + 1\times -1 + -1 \times 1 + -1\times-1\right) = 0.\) This is true about all of the cross terms, so we end up with \(\langle r_3^2\rangle = 1 + 1 + 1 = 3.\) We can use the square root of this quantity, \(\sqrt{\langle r_3^2\rangle} = \sqrt{3}\approx 1.732,\) as our measure of displacement.

In general, this calculation shows that \(\langle r_n^2\rangle = n.\) The more common notation is to change the step number \(n\) to time \(t\) so that \(\sqrt{\langle r_t^2\rangle} = \sqrt{t}.\) Playing loose with notation, we could say that the position scales like the square root of time: \(r \sim \sqrt{t}.\)

[GIF of a cloud of random walkers with expanding circle and a circle representing the sample mean]

This is a neat result—despite the fact that each soccer ball has its own trajectory, we can capture the average motion of any given soccer ball. Moreover, the motion doesn’t follow a usual function of \(t\) like \(r \sim t\) for constant velocity motion, or \(r \sim t^2\) in a gravitational field, it follows a fractional exponent of time \(r \sim t^{\frac12}.\)

## Random walk - recap

As cool as it is, the road ends there. It is difficult to use this modeling approach to extract much more about the system. Moreover, we know it’s missing some essential ingredients of realistic random motions.

In the random walk, our ball has no sense of ballistic momentum in that its current velocity is completely independent of its velocity at previous times. In reality, we know that even if a series of impulses are random in nature, it takes time to change momentum. This leads us to expect some signature of prior dynamics in the current dynamics; instead we have perfect independence.

Another shortcoming is that every random move has the same magnitude—one unit up or downfield. In reality, we know that the soccer ball (or our Brownian particle) can experience impulses of different strengths, resulting in different jump sizes.

There’s no easy way to incorporate these features into the discrete modeling strategy of the random walk. Let’s see if we can find a way to keep moving.