Fluid Flows and Infinite-Dimensional Manifolds (Part 2)

Or: ideal fluids—dry water?

guest post by Tim van Beek

Last time in this series, we set the stage by explaining infinite dimensional manifolds. Then we looked at a simple example: the inviscid Burgers equation. We saw this was the equation for geodesics in the diffeomorphism group of the circle.

Now let’s look at a more interesting example! It will still be a simplified model of fluid flow: it will describe an ideal fluid that is incompressible. I’ll start by explaining these concepts. We will then see how the equation of motion for ideal incompressible fluids can be interpreted as a geodesic equation.

En route I will also repeat some stuff from classical vector analysis, mostly for my own sake. The last time I seriously had to calculate with it was when I attended a class on “classical electrodynamics”, which was almost 15 years ago!

When we delve into differential geometry, it is always a good idea to look both at the “coordinate free” formulation using abstract concepts like differential forms, and also at the “classical vector analysis” part, that is best for calculating stuff once suitable coordinates have been chosen. Our fluid flows will take place in a smooth, orientable, compact, -dimensional Riemannian manifold possibly with a smooth boundary

I will frequently think of as an open set in or so I will use the globally defined coordinate chart of Euclidean coordinates on denoted by (and if needed) without further warning.

Before we continue: Last time our reader “nick” pointed out a blog post by Terence Tao about the same topic as ours, but—as could be expected—assuming a little bit more of a mathematical background: The Euler-Arnold equation. If you are into math, you might like to take a look at it.

So, let us start with the first important concept: the ‘ideal fluid’.

What is an ideal fluid?

When you are a small parcel in a fluid flow, you will feel two kinds of forces:

• external forces like gravity that are there whether or not your fellow fluid parcels surround you or are absent,

• internal forces that come from your interaction with the other fluid parcels.

If there is friction between you and other fluid parcels, for example, then there will be a force slowing down faster parcels and speeding up slower parcels. This is called viscosity. I already explained it back in the post Eddy Who? High viscosity means that there is a lot of friction: think of honey.

The presence of viscosity leads to shear stress whenever there are differences in the velocities of nearby fluid parcels. These lead to the formation of eddies and therefore to turbulence. This complicates matters considerably! For this reason, sometimes people like to simplify matters and to assume that the fluid flow that they consider has zero viscosity. This leads us to the physics definition of an ideal fluid:

An ideal fluid (as physicists say) is a fluid with zero viscosity.

As you can guess, I have also a mathematical definition in store for you:

An ideal fluid (as mathematicians say) is a fluid with the following property: For any motion of the fluid there is a (real valued) function called the pressure such that if is a surface in the fluid with a chosen unit normal the force of stress exerted across the surface per unit area at at time is

This implies that there is no force acting tangentially to the surface :

This picture is from

• Alexandre Chorin and Jerrold E. Marsden, A Mathematical Introduction to Fluid Mechanics, 3rd edition, Springer, New York 1993.

An ideal fluid cannot form eddies by itself without the help of external forces, nor can eddies vanish once they are present. So this simplification exclude a lot of very interesting phenomena, including everything that is usually associated with the term ‘turbulence’. But it is a necessary simplification for describing fluid flow using geodesic equations, because something moving along a geodesic doesn’t lose energy due to friction! So we will have to stick with it for now.

Historically, ideal fluids were almost exclusively studied during the 19th century, because the mathematics of viscous fluids seemed to be too hard—which it still is, although there has been a lot of progress. T his led to a schism of theoretical hydrodynamics and engineering hydrodynamics, because engineers had to handle effects like turbulence that ideal fluids cannot model. A very problematic aspect is that no body with a subsonic velocity feels any drag force in an ideal fluid. This is known as D’Alembert’s paradox. This means that one cannot find out anything about optimal design of ships or aircrafts or cars using ideal fluids as a model. This situation was overcome by the invention of ‘boundary layer techniques’ by the physicist Ludwig Prandtl at the beginning of the 20th century.

John von Neumann is cited by Richard Feynman in his physics lectures as having said that ideal fluids are like “dry water”, because they are so unlike real water. This is what the subtitle of this post alludes to. I don’t think this is quite fair to say. Along these lines one could say that quantum mechanics is the theory of stagnant light, because it does not include relativistic effects like quantum field theory does. Of course every mathematical model is always just an approximation to a Gedankenexperiment. And ideal fluids still have their role to play.

Maybe I will tell you more about this in a follow-up post, but before this one gets too long, let us move on to our second topic: incompressible fluids and ‘volume preserving’ diffeomorphisms.

What is an incompressible fluid flow?

If you are a parcel of an incompressible fluid, this means that your volume does not change over time. But your shape may, so if you start out as a sphere, after some time you may end up as an ellipsoid. Let’s make this mathematically precise.

But first note, that “incompressible” in the sense above means that the density of a given fluid parcel does not change over time. It does not mean that the density of the whole fluid is everywhere the same. A fluid like that is actually called homogeneous. So we have two different notions:

• incompressible means that the volume of an infinitesimal fluid parcel does not change as it moves along the fluid flow,

• homogeneous means that the density at a given time is everywhere the same, that is: constant in space.

This distinction is important, but for now we will study fluid flows that are both homogeneous and incompressible.

Let us see how we can make the notion of “incompressible” mathematically precise:

Remember from the last post: The flow of each fluid parcel is described by a path on parametrized by time, so that for every time there is a diffeomorphism

defined by the requirement that it maps the initial position of each fluid parcel to its position at time :

Now let’s assume our fluid flow is incompressible. What does that mean for the diffeomorphisms that describe the flow? Assuming that we have a volume form on these diffeomorphisms must conserve it:

For people who need a reminder of the concepts involved (which includes me), here it is:

Remember that is a smooth orientable Riemannian manifold of dimension A volume form is a -form that vanishes nowhere. In with Cartesian coordinates the canonical example would be

The dual basis of is denoted by in our example.

Given two manifolds and a differentiable map we can pull back a differential form on to one on via

For the übernerds out there: remember that we see the group of diffeomorpisms as a Fréchet Lie group modelled on the Fréchet space of vector fields on For those who would like to read more about this concept, try this:

• Karl-Hermann Neeb, Monastir Summer School: infinite-dimensional Lie groups.

is clearly a subgroup of It is less obvious, but true, that it is a closed subgroup and therefore itself a Lie group. What about its Lie algebra? For a vector field to give a flow that’s volume preserving, it must have zero divergence. So, the vector fields that form the tangent space consist of all smooth vector fields with zero divergence:

These vector fields form a vector space we denote by Remember stands for the tangent space at the identity element of the group which is the identity diffeomorphism of The tangent space at the identity of a Lie group is a Lie algebra, so is a Lie algebra.

I will need a little refresher about the definition of divergence. Then I will point you to a proof of the claim above, namely that zero-divergence vector fields form the Lie algebra of volume preserving diffeomorphism. This may seem obvious on an intuitive level, if you ever learned that the zero-divergence vector fields have ‘no sinks and no sources’, for example in a course on classical electromagnetism.

So, what is the divergence, again? You’ve probably seen it somewhere if you’ve survived reading this so far, but you may not have seen it in full generality.

The divergence of a vector field with respect to a volume form is the unique scalar function such that:

Here, is the contraction with Contraction means that you feed the vector in the first slot of the differential form and therefore reduce the function of vector fields to one of vector fields.

When we use our standard example we of course write a vector field as

where and are smooth real-valued functions. The divergence of is then

which we get if we plug in the expression for into the formula

So, how does one see that ‘zero divergence’ of a vector field is equivalent to ‘volume preserving’ for the flow it generates?

If we write

for the path of a fluid particle and $u$ for its velocity, then of course we have:

For a scalar function we get

Here is the inner product. The latter part is often written with the help of the nabla operator as

This is really just a handy short notation, there is no mystery behind it: it’s just like how we write the divergence as and the curl as

The operator

appears so often that it has its own name: it is called the material derivative.

Why ‘material’? Because if we follow a little bit of material—what we’re calling a parcel of fluid—something about it can change with time for two different reasons. First, this quantity can explicitly depend on time: that’s what the first term, is about. Second, this quantity can depend on where you are, so it changes as the parcel moves: that’s what is about.

Now suppose we have a little parcel of fluid. We’ve been talking about it intuitively, but mathematically we can describe it at time zero as an open set in our manifold. After a time it will be mapped by the fluid flow to

This describes how our parcel moves. We define the fluid to be incompressible if the volume of for all choices of is constant, that is:

If we write for the Jacobian determinant of then we have

So in a first step we get that a fluid flow is incompressible iff the Jacobian determinant is for all times, which is true iff is volume preserving.

It is not that hard to show by a direct calculation that

If you don’t want to do it yourself, you can look it up in a book that I already mentioned:

• Alexandre Chorin and Jerrold E. Marsden, A Mathematical Introduction to Fluid Mechanics, 3rd edition, Springer-Verlag, New York 1993.

This is the connection between ‘volume preserving’ and ‘zero divergence’! Inserting this into our equation of incompressibility, we finally get:

which is true for all open sets iff The equation of continuity for a fluid flow is:

This says that mass is conserved. Written with the material derivative it is:

So, since we’re assuming we get

which is what we intuitively expect, namely that the density is constant for a fluid parcel following the fluid flow.

Euler’s equation for the ideal incompressible fluid

The equation of motion for an ideal incompressible fluid is Euler’s equation:

is the pressure function mentioned in the mathematical definition of an ideal fluid above. As I already mentioned, to be precise I should say that we also assume that the fluid is homogeneous. This means that the density is constant both in space and time and therefore can be cancelled from the equation of motion.

If has a nonempty (smooth) boundary the equation is supplemented by the boundary condition that is tangential to

How can we turn this equation into a geodesic equation on ? Our strategy will be the same as last time when we handled the diffeomorphism group of the circle. We will define the necessary gadgets of differential geometry on using the already existing ones on First we define them on Then, for any diffeomorphism we use right translation by to define them on After that, we can use the version of the abstract version of the geodesic equation for right invariant metrics to calculate the explicit differential equation behind it.

Let us start with defining right invariant vector fields on A right invariant vector field is a vector field such that there is a such that In the following, we restrict ourselves to right invariant vector fields only.

We define the usual inner product of vector fields on just as last time:

The inner product used on the right is of course the one on

For two right invariant vector fields with and we define the inner product on by

This definition induces a right invariant metric on Note that it is right invariant because we are only considering volume preserving diffeomorphisms. It is not right invariant on the larger group of all diffeomorphims !

For an incompressible ideal fluid without external fields the only kind of energy one has to consider is the kinetic energy. The inner product that we use is actually proportional to the kinetic energy of the whole fluid flow at a fixed time. So geodesics with respect to the induced metric will correspond to Hamilton’s extremal principle. In fact it is possible to formulate all this in the language of Hamiltonian systems, but I will stop here and return to the quest of calculating the geodesic equation.

Last but not least, we define the following right invariant connection:

Here on the right is the connection on —sorry, this is not quite the same as the we’d been using earlier! But in or Euclidean space of any other dimension, is just another name for so don’t get scared.

Remember from last time that the geodesic equation says

where is the velocity vector of our geodesic, say

where is the curve describing our geodesic. We saw that for a right-invariant metric on a Lie group, this equation says

where the coadjoint operator is defined by

For simplicity, let us specialize to or an open set in there. What can we say about the right hand side of the above equation in this case? First, we have the vector identity

Since we are talking about divergence-free vector fields, we actually have

Also note that for a scalar function and the divergence-free vector field we have

The last term is zero because of our boundary condition that says that the velocity field is tangent to

So, now I am ready to formulate my claim that

for some yet undetermined scalar function This can be verified by a direct calculation:

What next? We can use the following 3 dimensional version of Green’s theorem for the curl operator:

That is, the curl operator is symmetric when acting on vector fields that have no component that is tangent to Note that I deliberately forgot to talk about function spaces that our vector fields need to belong to and the regularity assumptions on the domain and its boundary, because this is a blog post and not a math lecture. But the operators we use on vector fields obviously depend on such assumptions.

If you are interested in how to extend the symmetric curl operator to a self-adjoint operator, for example, you could look it up here:

• R. Hiptmair, P. R. Kotiuga, S. Tordeux, Self-adjoint curl operators.

Since our vector fields are supposed to be tangent to we have that the boundary term in our case is

because is normal, and therefore is tangent to so its inner product with the normal vector is zero.

So we can shift the curl operator from right to left like this:

In the last step we used the cyclicity of the relation of the vector product and the volume spanned by three vectors:

This verifies the claim, since the part does not contribute, as stated above.

And now, yet another vector identity comes to our rescue:

So, we finally end up with this:

for some function Why? Since the middle term is actually a gradient, we can absorb this summand and into one summand with a new function,

Thanks to this formula we derived, the abstract and elegant equation for a geodesic on any Lie group

becomes, in this special case

If we can convince ourselves that is the pressure of our fluid, we get Euler’s equation:

Wow! Starting with abstract stuff about infinite-dimensional Lie groups, we’ve almost managed to derive Euler’s equation as the geodesic equation on ! We’re not quite done: we still need to talk about the role of the function and why it’s minus the pressure. But that will have to wait for another post.

Metallic Hydrogen

 

See that gray stuff inside Jupiter? It’s metallic hydrogen—according to theory, anyway.

But how much do you need to squeeze hydrogen before the H₂ molecules break down, the individual atoms form a crystal, and electrons start roaming freely so the stuff starts conducting electricity and becomes shiny—in short, becomes a metal?

In 1935 the famous physicist Eugene Wigner and a coauthor predicted that 250,000 times normal Earth atmospheric pressure would do the job. But now they’ve squeezed it to 3.6 million atmospheres and it’s still not conducting electricity! Here’s a news report:

• Probing hydrogen under extreme conditions, PhysOrg, 13 April 2012.

and here’s the original article, which unfortunately ain’t free:

• Chang-Sheng Zha, Zhenxian Liu, and Russell J. Hemley, Synchrotron infrared measurements of dense hydrogen to 360 GPa, Phys. Rev. Lett. 108 (2012), 146402.

Three phases of highly compressed solid hydrogen are known, with phase I starting at 1 million atmospheres and phase III kicking in around 1.5 million. I would love to know more about these! Do you know where to find out? Some people also think there’s a liquid metallic phase, and a superconducting liquid metallic phase. In fact there are claims that liquid metallic hydrogen has already been seen:

• W.J. Nellis, S.T. Weir and A.C. Mitchell, Metallization of fluid hydrogen at 140 GPa (1.4 Mbar) by shock compression, Shock Waves 9 (1999), 301–305.

1.4 Mbar, or megabar, is about 1.4 million atmospheres of pressure. Here’s the abstract:

Abstract. Shock compression was used to produce the first observation of a metallic state of condensed hydrogen. The conditions of metallization are a pressure of 140 GPa (1.4 Mbar), 0.6 g/cm (ninefold compression of initial liquid-H density), and 3000 K. The relatively modest temperature generated by a reverberating shock wave produced the metallic state in a warm fluid at a lower pressure than expected previously for the crystallographically ordered solid at low temperatures. The relatively large sample diameter of 25 mm permitted measurement of electrical conductivity. The experimental technique and data analysis are described.

Apprently the electric resistivity of fluid metallic hydrogen is about the same as the fluid alkali metals cesium and rubidium at 2000 kelvin, right when they undergo the same sort of nonmetal-metal transition. Wow! So does that mean that at 2000 kelvin but at lower pressures, these elements don’t act like metals? I hadn’t known that!

Another reason this is interesting is that if you look at hydrogen on the periodic table, you’ll see it can’t make up its mind whether it’s an alkali metal—since its outer shell has just one electron in it—or a halogen—since its outer shell is just one electron short from being full! You could say compressing hydrogen until it becomes metallic is like trying to get it to break down and admit its an alkali metal.

Apparently the metal-nonmetal transition for for liquid cesium, rubidium and hydrogen all happen when the stuff gets squashed so much that the distance between atoms goes down to about 0.3 times the size of these atoms in vacuum… where by ‘size’ I mean the Bohr radius of the outermost shell.

How did Huntington and Wigner get their original calculation so wrong? I don’t know! Their original paper is here:

• Eugene Wigner and Hillard Bell Huntington, On the possibility of a metallic modification of hydrogen J. Chem. Phys. 3 (1935), 764-771.

It’s not free; I guess the American Institute of Physics is still trying to milk it for all it’s worth. One interesting thing is that they assumed the crystal stucture of metallic hydogen would be a ‘body-centered cubic’… it’s rather hard to compute these things from scratch without computers. But this more recent paper claims that a diamond cubic is energetically favored at 3 million atmospheres:

• Crystal structure of atomic hydrogen, Phys. Rev. Lett. 70 (1993), 1952–-1955.

I explained the diamond cubic crystal structure in my recent post about ice. Remember, it looks like this:

Since the body-centered cubic is one of the crystal lattices I didn’t talk about in that post, let me tell you about it now. It’s built of cells that look like this:

… which explains its name. In the same style of drawing, the face-centered cubic looks like this:

In my post about ice, I mentioned that if you pack equal-sized spheres with centers at points in the face-centered cubic lattice, you get the maximum density possible, namely about 74%. The body-centered cubic does slightly worse, about 68%.

So, I always thought of it as a kind of a second-best thing. But apparently it’s the best in some ‘sampling’ problems where you’re trying to estimate a function on space by measuring it only at points in your lattice! That’s because its dual is the face-centered cubic.

Eh? Well, the dual of a lattice in a vector space consists of all vectors in the dual vector space such that is an integer for all points in The dual of a lattice is again a lattice, and taking the dual twice gets you back where you started. Since the Fourier transform of a function on a vector space is a function on the dual vector space, I don’t find it surprising that this is related to sampling problems. I don’t understand the details, but I bet I could find them here:

• Alireza Entezari, Ramsay Dyer and Torsten Möller, From sphere packing to the theory of optimal lattice sampling.

So: from the core of Jupiter to Fourier transforms! Yet again, everything seems to be related.

Khumbu Icefall and the Valley of Silence

National Geographic has a blog written by people who are now climbing Mount Everest. Here’s Sam Elias training in the Khumbu Icefall near the Everest Base Camp:

As usual, it’s the Sherpas who impress me most:

Years of experience, or maybe the mountain itself, had told the Sherpas that passing through the Ballroom on this day was not a good idea, something would happen. “Big ice will fall.” Panuru’s words echoed in my head. “How do they know?” I wondered.

I was sitting in my tent fitting my crampons onto my boots when I heard it. I know the sound now. Before, when the loud rumbling began I instinctively thought of a giant semi barreling down a highway. But there are no vehicles here.

Also:

Every year, the route through the Khumbu is set by the “ice doctors,” a small team of Sherpas who take mortal risks to navigate the safest passage through the Icefall, putting up ropes in the steep sections and stretching ladders across the abyss-like crevasses.

Crossing the ladders is an adventure for some. For the Sherpas, setting them up is a job.

Khumbu Icefall

Suppose you take the southeast route to Mount Everest, on the Nepal side. When you climb up from Base Camp, the first thing you’ll hit is the Khumbu Icefall, a crazy and ever-changing mass of ice at the bottom of the Khumbu Glacier:

As the National Geographic blog put it:

Like a gargantuan bulldozer, the Khumbu glacier plows down off the Lhotse Face between Mounts Everest and Nuptse. Dropping over a cliff just above Base Camp, this mile-wide river of ice shatters into building-size blocks and steeple-size spires called seracs. It’s riven with cracks called crevasses that can be hundreds of feet deep. To reach our expedition’s two goals — the Southeast Ridge and the West Ridge, which both begin atop the Khumbu glacier in the Western Cwm — we must travel up through this labyrinth of raging ice.

To cross the crevasses, you use bridges that the Sherpas have made by lashing ladders together with rope. Here’s Nima Dorje Tamang crossing one. The clouds are like a ceiling… but there’s no floor:

The picture above is again from National Geographic.

The glacier advances about a meter each day around here. Most climbers try to cross before the sun rises, when the cold keeps things frozen. As the intense sunlight warms things, the icefall becomes more dangerous. Blocks of ice tumble down the glacier from time to time, ranging in size from cars to houses… and sometimes entire large towers of ice collapse. They say bodies of people who die in here sometimes show up at the base of the icefall years later.

Here’s Kenton Cool talking about the Khumbu Icefall. “It can implode underneath you, it can drop on you above – or god forbid, you can fall into its inner depths, never to be seen again.”

And this is photographer Leo Dickinson speaking about the dangers of this place. Look at the fellow poking at snow with a pick around 0:58, revealing that it would be deadly to step there!

 

The Valley of Silence

Suppose you succeed in crossing the Khumbu Icefall—including the last crevasse, shown in this photo by Olaf Rieck. Then you have reached the Western Cwm, also known as the Valley of Silence:

In the middle background is Lhotse. At far right you see a bit of Nuptse. And at left there’s Sāgārmatha, also known in Tibetan as Chomolungma… or in English, Mount Everest.

‘Cwm’, pronounced ‘coom’, is Welsh for a bowl shaped valley, also known as a ‘cirque’. This one is a 4-kilometer-long valley carved out by the Khumbu Glacier, which starts at the base of Lhotse. It’s the easiest way to approach Everest from the southeast. However, it’s cut by massive crevasses that bar entrance to the upper part: here you must cross to the far right, over to the base of Nuptse, and through a narrow passageway known as the Nuptse corner.

It’s called the Valley of Silence because it’s often windless and deathly quiet. On days like that, the surrounding snow-covered slopes surrounding are so bright that the valley becomes a kind of solar oven, with temperatures soaring to 35 °C (95 °F) despite an elevation of 6000 to 6800 metres (19,600-22,300 feet). But when sun turns to shade, the temperature can plummet to below freezing in minutes!

The photo above was taken by the Moving Mountains Trust. See the people? You may need to click for a bigger version! For more, see:

• Alan Arnette, Life in the Western Cwm.

Want to go further? When you’ve reached Base Camp II near the top of the Western Cwm, you still have 2300 meters to climb… and now it gets steep! I’m sorry, I’m quitting here and heading back down—it’s my bedtime. Good luck!

For more

Cut your carbon footprint. Travel virtually:

• Mount Everest summit—interactive 360 degree panorama.

• Reality Maps viewer for Everest.

Michael Murphy writes:

I had become intrigued by the story of Marco Siffredi, a French snowboarder who was the first to successfully descend Everest on a snowboard via the Norton Couloir. His second attempt to descend a far more serious route, the Hornbein Couloir ended in his demise.

Here’s the video of him leaving the summit. I used Reality Maps to trace his route. It is no wonder he did not make it.

Enormous Integers

Sometimes math problems have unexpected answers. For example, consider the integral of this infinite product:

The answer matches up to its 43rd digit. But it’s not actually Weird, huh? But it’s not a coincidence; there’s an explanation.

Here’s a puzzle due to the logician Harvey Friedman. It too has an unexpected answer.

Say you have a finite alphabet to write with. How long can a word be if no block of letters in this word, from the nth letter to the 2nth, is allowed to appear as a subsequence of a bigger block from the mth letter to the 2mth?

If you have just one letter, this is the longest it can be:

AAA

If you have two, this is the longest it can be:

ABBBAAAAAAA

Puzzle: How long can the word be if you have three letters in your alphabet?

Friedman showed there’s still a finite upper bound on how long it can be. But, he showed it’s incomprehensibly huge!

Now Friedman is one of the world’s experts on large cardinals—large infinite numbers. So when he says a finite number is incomprehensibly huge, you sit up and listen. It’s like seeing a seasoned tiger hunter running through the jungle with his shotgun, yelling “Help! It’s a giant ant!”

The answer to the puzzle involves the 7th Ackermann function. The first Ackermann function is just doubling:

The second Ackermann function of n is 2 to the nth power:

To calculate the third Ackermann function of n, you write down a tower of n 2’s. For example:

And so on: to calculate the (k+1)st Ackermann function applied to n, you apply the kth Ackermann function n times to the number 1:

where there are n of these ‘s.

For example, with a little work you can show

where the number of 2’s in the tower is 65,536.

But this is still puny, compared to the answer to the puzzle!

In 1998 Friedman show that the longest word you can write with 3 letters, following the rule I described, is at least A₇(184) letters long.

As he notes, when we reach numbers this large:

… a profound level of incomprehensibility emerges. The definitions are not incomprehensible, but the largeness is incomprehensible. These higher levels of largeness blur, where one is unable to sense one level of largeness from another.

I don’t understand his proof, but you can see it here:

• Harvey M. Friedman, Long finite sequences, 8 October 1998.

When I posted about this on Google+, Andrés Caicedo helped me notice that later in the paper, Friedman goes further. In fact, A₇(184) is a ridiculous underestimate of the true answer! And apparently now there’s an upper bound on the answer, too.

When Andrés Caicaido heard me say the word could be at least A₇(184) letters long, he wrote:

Much much longer!

I wrote a bit on this:

• A lower bound for n(3), 14 February 2012.

Using some computer assistance from Ramdall Dougherty, Friedman shows that a lower bound is A₇₁₉₈(158386). I asked Dougherty, and got a copy of his files. I am hoping to find some decent time in the not too distant future to examine them in detail and see how far this can be pushed.

Friedman also found an upper bound, A_n(n), where n = A₅(5). He mentions this in a note from 2001, but the note is a series of announcements with no proofs. Actually, it is exciting stuff, in spite of the lack of proofs. He discusses some other simply defined combinatorial constructions that give numbers much much longer than this one:

• Harvey M. Friedman, Lecture notes on enormous integers, 22 November 2011.

So, we have a well-defined enormous integer, and we know—or at least Friedman claims—that it’s somewhere between and That’s an impressively large spread. I wonder how accurately we will ever know this number.

I should add that beside the sheer shock value of seeing such a huge number show up as the answer to a simply stated puzzle, there’s some deep math lurking in these games. It gets rather scary. Near the end of his lecture notes, Friedman mentions a number so big that any proof of its existence—in a certain system of mathematics—has an incomprehensibly large number of symbols.

But if all this talk of enormous integers isn’t scary enough for you, just wait until you read about dark integers:

• Greg Egan, Dark integers, 2007.

Dark matter, dark energy . . . dark integers. They’re all around us, but we don’t usually see them, because they don’t quite play by the rules.

Personal Rapid Transportation

guest post by Todd McKissick

We all can’t wait for High Speed Rail to come to our town. Whether we’re referring to fast traditional trains on wheels (HSR) or those that float down the track on magnetic fields (maglev), this is the 21st Century, so what most people desire is the full featured deal. Anything less is just another compromise. They have been touted now for 40 years that I know of.

But, as usual, it begs the question: is this really the best solution? I’ve found a lot of different solutions and reviewed everyone of them, but only two stand head and shoulders above the rest.

First, let’s look at some of the specifics of what we’re asking for. It runs really fast so there’s lots of possibility of crossing the country in a couple hours. It gets its efficiency mostly from packing lots of cargo into a very efficient vehicle as most trains do. It’s clearly better than getting 25, 60 or 85 passenger-kilometers per liter in a fully loaded airplane, Suburban or Prius. (That’s 65, 140 or 200 passenger-miles per gallon.) As long as it’s comfortable, this is all good stuff.

Unfortunately, to accomplish this, maglev takes a fairly standard sized train and hovers it over a massive rail with thousands of high-powered electromagnets to float this 80+ tonne piece of machinery from one town to another. It accelerates slowly and brakes slowly, unless you want to double the incredible amount of power it uses already. Its rail system consists of hundreds of tons of concrete per 30 meter segment to make the rail and the support beams, and we all know that concrete is horrible for the environment. And lastly, it costs a billion dollars to build each couple miles of the system.

I’m thinking there’s a better way.

To truly combat the automobile, the airplane and other forms of transportation that use lots of fossil fuels, let’s first look at the last mile segment. This is from your door to the store, work or school lobby or even to your friend’s door. Check out this picture from a company called SkyTran.

It’s called Personal Rapid Transit (PRT for short), and it’s courting numerous locations around the world right now. Basically, it’s a small carbon fiber pod that holds two seats and an iPad. This pod hangs from a small rail which then hangs from arms on regular telephone-style poles. At certain locations, a set of steps to a platform is placed to allow people to call for one and hop in. These terminals are cheap enough that they can be placed so that you’re never more than about 2 blocks from one in town or a couple miles in the country. In fact, they can be incorporated into lobbies, shopping malls, schools, sports arenas and even the higher floors of high rises because they are really just a landing to stand on (or ramp for those using wheels). Up to 350 kilogram pallets can be shipped autonomously. The one-way rails pass each other at different elevations so collisions are avoided while off-ramps to terminals simply drop down to a separate rail to stop on. This allows following traffic to continue at full speed while you merge on or off.

Now for the fun stuff. Once you get in and select your destination, it takes you straight to your end destination. Initial speeds are advertised at 70 kilometers per hour (45 mph) for town and 240 km/h (150 mph) for country but the top speed is well over 320 km/h (200 mph). It’s really only limited by wind resistance. This means you can board one at your sidewalk and step directly onto the upper level of a sports stadium 80 miles away in 25 minutes. Need a 2 hour nap before reaching the kids’ house? You’ll have to be farther than 450 kilometers away. When you figure in all the time needed for the various legs required to travel, this is the fastest way to travel any distance between 6 and 600 kilometers. Ya just gotta love avoiding parking at the airport to switch to a plane, because you can get there faster by avoiding the plane altogether—not to mention airport security!

The ride is perfectly smooth and quiet and offers iPad access with wifi support for your own devices. The pods can even be ‘ganged’ for when you have more than two riders (or kids) so that you load your kids in one car, connect it to your car virtually and then follow them to the destination. Along the way, it can switch the order so you can get out first for safety. Each rail has the capacity of a 3-lane freeway. Since the rail is upside down, balancing suspension is not needed because there’s no chance of the pod falling off. In other words, the curves are all banked for the set speed so the passengers feel no side force.

The energy required to operate it is equivalent to getting 85 kilometers per liter (200 mpg) in a loaded car because it is lifting small enough loads to take advantage of ‘passive levitation’. This is a type of maglev that uses drag to levitate the car. This kicks in at about 1 kilometer per hour, raises the car off the wheels, and diminishes to negligible drag once you pass 22 kilometers per hour (14 mph). Coupling this with regenerative braking means that you really only need to provide the energy to push against the wind. In fact a canopy of solar cells over it could power the entire system during rush hour for free. The rail is designed to also incorporate transmission and/or distribution power lines, 3G / WiFi internet connectivity (including backbone and user distribution) and possibly other utility services. A nationwide network installation could reduce US oil imports from 12 mbbl/day to around 3 mbbl, cutting the cost we pay to OPEC from $700B/yr to $175B/yr.

The cost of building such a system is still fairly high at €4.7/km ($10M/mile), but even that’s 1/19th of the cost of most HSR, and it’s and expected to come way down. The cost to the riders for a privately funded system (before profit, of course) is about €0.02/km, or 4 US cents per mile. Compare this to the cost of a personal vehicle which comes in at 14 times more. When you envision the scope this could be implemented, you can see that many targeted communities could do away with roads altogether, and opt for wider bike paths (to accommodate the occasional moving truck) and more nature. Parking lots could be located in cheap real estate areas or eliminated altogether. Delivery trucks could be replaced with individual pallet deliveries directly inside the factory. In short, all deliveries would eliminate the return trip. Cargo sharing could be implemented along with interstate passenger ‘opportunity’ trips for low cost travel for those who can’t afford travel. If there were a public outcry for this system and we decided to install it nationwide, each community could fund a significant chunk of it from existing road funds with no change in taxes. Or private individuals could invest to install it on a for profit basis.

I mapped out my small 12,000 person town, hit all the major points directly and put a terminal within 2 1/2 blocks of every house, for a total price of $160 million. The annual cost to pay it off completely in 10 years would be around $5-6,000 per family before considering the added savings like providing transportation to those in our population that can’t legally drive now. That’s 60% of what people spend on cars today. What better way to help young adults get their lives started without debt? All infrastructure additions can work in parallel with existing roads and utilities and installation time is roughly 2-3 miles of rail per week per crew. You do the math. If there were a global push behind PRT, we could cut our energy dependence, our environmental impact (not to mention the impact on people’s lives) and bring nature back into our communities.

That covers the short range travel but we still have long distance air, rail and international travel to address. Enter Evacuated Tube Transport.

This system suspends a long vacuum tube overhead or under water to guide mini-trains of 6 passengers on extremely high speed, long range trips. By evacuating the entire length of tube, most of the wind drag is removed, allowing it to travel at speeds up to 6,500 kilometers per hour (4,000 mph) with a maximum of 1 g of acceleration in any direction. The ET3 website suggests that intra-state travel will run at around 550 kilometers per hour (allowing for 2.5 kilometer radius u-turns) while the higher speed legs across the country or under the ocean surface can do a loop in a 320 kilometer radius.

A sample trip from L.A. to N.Y. would take 3 minutes to accelerate over the first 160 kilometers, 42 minutes to cruise the middle and 3 more minutes to slow down while capturing the remaining momentum as electrical regeneration.

Of course, riding in a vacuum requires a pod capable of safely withstanding dangerous pressures, but even transoceanic underwater travel poses no problems we don’t already deal with for other causes. It would be worth the ride for just the scenery if there were transparent sides on the tube, but one has to wonder what you could actually see at that speed below sea level. And since a 6 hour trip across the Pacific would not include any stops, there are some obvious human considerations which would need to be dealt with. Even considering these issues, the economics are sound in dollars, resources and energy. As you can see on their page comparing to standard trains, the ET3 system far surpasses even the Skytran for efficient long distance transport. One could only hope for the merger of the two, where you hop in a Skytran car on your corner, zoom straight into a ET3 loading system, jet up to high speed, cross half the globe, reverse the process at some foreign destination and charge $100 on your card, all in a couple hours.

If you think this is is all hype and fairy land, you might want to search around for some of the projects around the world that are reviewing these two little gems. It’s only a lack of popular opinion that’s holding these two back. Let’s make this happen with a little viral support!

Ice

Water is complicated and fascinating stuff! There are at least sixteen known crystal phases of ice, many shown in this diagram based on the work of Martin Chaplin:

(Click for a larger version.)

For example, ordinary ice is called ice I_h, because it has hexagonal crystals. But if you cool it below -37 °C, scientists believe it will gradually turn into a cubic form, called ice I_c. In this form of ice the water molecules are arranged like carbons in a diamond! And apparently some of it can be found in ice crystals in the stratosphere.

But if you wait longer, ice below -37 °C will turn into ice XI. The transformation process is slow, but ice XI has been found in Antarctic ice that’s 100 to 10,000 years old. This ice is ferroelectric, meaning that it spontaneously become electrically polarized, just like a ferromagnet spontaneously magnetizes.

That’s the usual textbook story, anyway. The true story may be even more complicated:

• Simon Hadlington, A question mark over cubic ice’s existence, Phys.org, 9 January 2012.

But now a team led by Benjamin Murray at the University of Leeds has carried out work that suggests that cubic ice may not in fact exist. The researchers searched for cubic ice by suspending water droplets in oil and gradually cooling them to -40 °C while observing the X-ray diffraction pattern of the resulting crystals. “We modelled the diffraction pattern we obtained and compared it to perfect cubic and perfect hexagonal ice, and it was clearly neither of them,” Murray says. “Nor is it merely a mixture of the two. Rather it is something quite distinct.”

Analysis of the diffraction data shows that in the ice crystals the stacking of the atomic layers is disordered. “The crystals that form have randomly stacked layers of cubic and hexagonal sequences,” Murray says. “As each new layer is added, there is a 50% probability of it being either hexagonal or cubic.” The result is a novel, metastable form of ice with a stacking-disordered structure.

Re-examination of what had previously been identified as cubic ice suggests that this was stacking-disordered structures too, Murray says. “Cubic ice may not exist.”

Even plain old ice I_h is surprisingly tricky stuff. The crystal structure is subtle, first worked out by Linus Pauling in 1935. Here’s a nice picture by Steve Lower:

There’s a lot of fun math here. First focus on the oxygen atoms, in red. What sort of pattern do they form?

Close-packings

To warm up for that question, it really helps to start by watching this video:

It’s not flashy, but it helped me understand something that had been bothering me for decades! This is truly beautiful stuff: part of the geometry of nature.

But I’ll explain ice I_h from scratch—and while I’m at it, ice I_c.

First we need to think about packing balls. Say you’re trying to pack equal-sized balls as densely as possible. You put down a layer in a hexagonal way, like the balls labelled ‘a’ here:

Then you put down a second hexagonal layer of balls. You might as well put them in the spots labelled ‘b’. There’s another choice, but the symmetry of the situation means it doesn’t matter which you pick!

The fun starts at the third layer. If you want, you can put these balls in the positions labelled ‘a’, directly over the balls in the first layer. And you can go on like this forever, like this:

-a-b-a-b-a-b-a-b-a-b-a-

This gives you the hexagonal close packing.

But there’s another choice for where to put the balls in the third layer… and now it does matter which choice you pick! The other choice is to put these balls in the spots labelled ‘c’:

And then you can go on forever like this:

-a-b-c-a-b-c-a-b-c-a-b-c-

This gives you the cubic close packing.

There’s also an uncountable infinity of other patterns that all give you equally dense packings. For example:

-a-b-a-c-b-a-c-b-c-b-a-c-

You just can’t repeat the same letter twice in a row!

To review, here’s the difference between the hexagonal and cubic close packings:

But what’s ‘cubic’ about the cubic close packing? Well, look at this. At left we see a bit of a hexagonal close packing, while at right we see a bit of a cubic close packing:

The dashed white circle at right shows what you shouldn’t do for a cubic close packing. But the red lines at right show why it’s called ‘cubic’! If you grab that cube and pull it out, it looks like this:

As you can see, there’s a ball at each corner of a cube, but also one at the middle of each face. That’s why this pattern is also called a face-centered cubic.

If we shrink the balls a bit, the face-centered cubic looks like this:

What used to confuse me so much is that the cubic close packing looks very different from different angles. For example, when you see it like this, the cubical symmetry is hard to spot:

But it’s there! If you stack some balls in the hexagonal close packing, they look different:

And to my eye, they look uglier—they’re less symmetrical.

Anyway, perhaps now we’re in a better position to understand what Benjamin Murray meant when he said:

As each new layer is added, there is a 50% probability of it being either hexagonal or cubic.

I’m guessing that means that at each layer we randomly make either of the two choices I mentioned. I’m not completely sure. My guess makes the most sense if the ice is growing by accretion one layer at a time.

Diamonds

But there’s more to ice than close-packed spheres. If we ignore the hydrogen atoms and focus only on the oxygens, the two kinds of ice we’re talking about—the hexagonal ice I_h and the cubic ice I_c, assuming it exists—are just like two kinds of diamond!

The oxygens in ice I_c are arranged in the same pattern as carbons in an ordinary diamond:

This pattern is called a diamond cubic. We can get it by taking the cubic close packing, making the balls smaller, and then putting in new ones, each at the center of a tetrahedron formed by the old ones. I hope you can see this. There’s a ball at each corner of the cube and one at the center of each face, just as we want for a face-centered cubic. But then there are are 4 more!

But what about good old ice I_h? Here the oxygens are arranged in the same pattern as carbons in a hexagonal diamond.

You’ve heard about diamonds, but you might not have heard about hexagonal diamond, also known as lonsdaelite. It forms when graphite is crushed… typically by a meteor impact! It’s also been made in the lab.

Its crystal structure looks like this:

We get this pattern by taking a hexagonal close packing, making the balls smaller, and then putting in new ones, each at the center of a tetrahedron formed by the old ones. Again, I hope you can see this!

Entropy

Now we understand how the oxygens are arranged in good old ice I_h. They’re the red balls here, and the pattern is exactly like lonsdaelite, though viewed from a confusingly different angle:

But what about the hydrogens? They’re very interesting: they’re arranged in a somewhat random way!

Pick any oxygen near the middle of the picture above. It has 4 bonds to other oxygens, shown in green. But only 2 of these bonds have a hydrogen near your oxygen! The other 2 bonds have hydrogens far away, at the other end.

There are many ways for this to happen, and they’re all allowed—so ice is like a jigsaw puzzle that you can put together in lots of ways!

So, even though ice is a crystal, it’s disordered. This gives it entropy. Figuring out how much entropy is a nice math puzzle: it’s about 3/2 times Boltzmann’s constant per water molecule. Why?

Rahul Siddharthan explained this to me over on Google+. Here’s the story.

The oxygen atoms form a bipartite lattice: in other words, they can be divided into two sets, with all the neighbors of an oxygen atom from one set lying in the other set. You can see this if you look.

Focus attention on the oxygen atoms in one set: there are of them. Each has 4 hydrogen bonds, with two hydrogens close to it and two far away. This means there are

allowed configurations of hydrogens for this oxygen atom. Thus there are configurations that satisfy these atoms.

But now consider the remaining oxygen atoms: in general they won’t be satisfied: they won’t have precisely two hydrogen atoms near them). For each of those, there are

possible placements of the hydrogen atoms along their hydrogen bonds, of which 6 are allowed. So, naively, we would expect the total number of configurations to be

Using Boltzmann’s ideas on entropy, we conclude that

where is Boltzmann’s constant. This gives an entropy of 3.37 joules per mole per kelvin, a value close to the measured value. But this estimate is ‘naive’ because it assumes the 6 out of 16 hydrogen configurations for oxygen atoms in the second set can be independently chosen, which is false. More complex methods can be employed to better approximate the exact number of possible configurations, and achieve results closer to measured values.

By the way: I’ve been trying to avoid unnecessary jargon, but this randomness in ice I_h has such a cool name I can’t resist mentioning it. It’s called proton disorder, since the hydrogens I’ve been talking about are really just hydrogen nuclei, or protons. The electrons, which form the bonds, are smeared all over thanks to the wonders of quantum mechanics.

Higher pressures

If I were going to live forever, I’d definitely enjoy studying and explaining all sixteen known forms of ice. For example, I’d tell you the story of what happens to ordinary ice as we slowly increase the pressure, moving up this diagram until we hit ice XI:

Heck, I’d like to know what happens at every stage as we crush water down to neutronium! Unfortunately I’m mortal, with lots to do. So I recommend these:

• Martin Chaplin, Water phase diagram.

• Norman Anderson, The many phases of ice.

But there’s a bit of news I can’t resist mentioning…

Researchers at Sandia Labs in Albuquerque, New Mexico have been using the Z Machine to study ice. Here it is in operation, shooting out sparks:

Click for an even more impressive image of the whole thing.

How does it work? It fires a very powerful electrical current—about 20 million amps—into an array of thin, parallel tungsten wires. For a very short time, it uses a power of 290 terawatts. That’s 80 times the world’s total electrical power output! The current vaporizes the wires, and they turn into tubes of plasma. At the same time, the current creates a powerful magnetic field. This radially compresses the plasma tubes at speeds of nearly 100,000 kilometers per hour. And this lets the mad scientists who run this machine study materials at extremely high pressures.

Back in 2007, the Z Machine made ice VII, a cubic crystalline form of ice that coexists with liquid water and ice VI at 82 Celsius and about 20,000 atmospheres. Its crystal structure is theorized to look like this:

But now the Z Machine is making far more compressed forms of ice. Above 1 million atmospheres, ice X is stable… and above 8 million atmospheres, a hexagonal ferroelectric form called ice XI comes into play. This has a density over 2.5 times that of ordinary water!

Recent experiments with the Z Machine seem to show water is 30% less compressible than had been thought under conditions like those at the center of Neptune. These experiments seem to match new theoretical calculations, too:

• Experiments may force revision of astrophysical models: ice giant planets have more water volume than believed, Science Daily, 19 March 2012.

While it’s named after the Roman god of the seas, and it’s nice and blue, Neptune’s upper atmosphere is very dry. However, there seems to be water down below, and it has a heart of ice. A mix of water, ammonia and methane ice surrounding a rocky core, to be precise! But the pressure down there is about 8 million atmospheres, so there’s probably ice X and ice XI. And if that ice is less compressible than people thought, it’ll force some changes in our understanding of this planet.

Not that this matters much for most purposes. But it’s cool.

And if you don’t love math, this might be a good place to stop reading.

A little more math

Okay… if you’re still reading, you must want more! And having wasted the day on this post, I might as well explain a bit of the math of the crystal structures I described. I’ll want to someday; it might as well be now.

A lattice in n-dimensional Euclidean space is a subset consisting of all linear combinations of basis vectors. The centers of the balls in the cubic close packing form a lattice called the face-centered cubic or A₃ lattice. It’s most elegant to think of these as the vectors in 4-dimensional space having integer components and lying in the hyperplane

That way, you can see this lattice has the group of permutations of 4 letters as symmetries. Not coincidentally, this is the symmetry group of the cube.

In this description, you can take the three basis vectors to be

Note that they all have the same length and each lies at a 120° angle from the previous one, while the first and last are at right angles.

We can also get the A₃ lattice by taking the center of each red cube in a 3d red-and-black checkerboard. That means taking all vectors in 3-dimensional space with integer coefficients that sum to an even number. In this description, you can take the three basis vectors to be

This seems like an ugly choice, but note: they all have the same length and each lies at a 120° angle from the previous one, while the first and last are at right angles. So, we know it’s the A₃ lattice!

We see more symmetry if we look at all the shortest nonzero vectors in this description of the lattice:

These form the 12 midpoints of the edges of a cube, and thus the corners of a cuboctahedron:

So, in the cubic close packing each ball touches 12 others, centered at the vertices of a cuboctahedron, as shown at top here:

Below we see that in the hexagonal close packing each ball also touches twelve others, but centered at the vertices of a mutant cuboctahedron whose top has been twisted relative to its bottom. This shape is called a triangular orthobicupola.

I should talk more about the lattice associated to the hexagonal close packing, but I don’t understand it well enough. Instead I’ll wrap up by explaining the math of the diamond cubic:

The diamond cubic is not a lattice. We can get it by taking the union of two copies of the A₃ lattice: the original lattice and a translated copy. For example, we can start with all vectors with integer coefficients summing to an even number, and then throw in all vectors This is called the D₃⁺ pattern.

The A₃ lattice gives a dense-as-possible packing of balls, the cubic close packing, with density

The D₃⁺ pattern, on the other hand, gives a way to pack spheres with a density of merely

This is why ordinary ice actually becomes denser when it melts. It’s not packed in the diamond cubic pattern: that would be ice I_c. Ordinary ice is packed in a similar pattern built starting from the hexagonal close packing! But the hexagonal close packing looks just like the cubic close packing if you only look at two layers… so this similar pattern gives a packing of balls with the same density as the diamond cubic.

Finally, for a tiny taste of some more abstract math: in any dimension you can define a checkerboard lattice called D_n, consisting of all n-tuples of integers that sum to an even integer. Then you can define a set called D_n⁺ by taking the union of two copies of the D_n lattice: the original and another shifted by the vector

D_n⁺ is only a lattice when the dimension n is even! When the dimension is a multiple of 4, it’s an integral lattice, meaning that the dot product of any two vectors in the lattice is an integer. It’s also unimodular, meaning that the volume of the unit cell is 1. And when the dimension is a multiple of 8, it’s also even, meaning that the dot product of any vector with itself is even.

Now, D₈⁺ is very famous: it’s the only even unimodular lattice in 8 dimensions, and it’s usually called E₈. In week193, I showed you that in the packing of balls based on this lattice, each ball touches 240 others. It’s extremely beautiful.

But these days I’m trying to stay more focused on the so-called ‘real world’ of chemistry, biology, ecology and the like. This has a beauty of its own. In 3 dimensions, D₃ = A₃ is the face-centered cubic. D₃⁺ is the diamond lattice. So, in a way, diamonds come as close to E₈ as possible in our 3-dimensional world.

The diamond cubic may seem mathematically less thrilling than E₈: not even a lattice. Ordinary ice is even more messy. But it has a different charm—the charm of lying at the bordeline of simplicity and complexity—and the advantage of manifesting itself in the universe we easily see around us, with a rich network of relationships to everything else we see.

This Week’s Finds (Week 319)

This week I’m trying something new: including a climate model that runs on your web browser!

It’s not a realistic model; we’re just getting started. But some programmers in the Azimuth Project team are interested in making more such models—especially Allan Erskine (who made this one), Jim Stuttard (who helped me get it to work), Glyn Adgie and Staffan Liljgeren. It could be a fun way for us to learn and explain climate physics. With enough of these models, we’d have a whole online course! If you want to help us out, please say hi.

Allan will say more about the programming challenges later. But first, a big puzzle: how can small changes in the Earth’s orbit lead to big changes in the Earth’s climate? As I mentioned last time, it seems hard to understand the glacial cycles of the last few million years without answering this.

Are there feedback mechanisms that can amplify small changes in temperature? Yes. Here are a few obvious ones:

• Water vapor feedback. When it gets warmer, more water evaporates, and the air becomes more humid. But water vapor is a greenhouse gas, which causes additional warming. Conversely, when the Earth cools down, the air becomes drier, so the greenhouse effect becomes weaker, which tends to cool things down.

• Ice albedo feedback. Snow and ice reflect more light than liquid oceans or soil. When the Earth warms up, snow and ice melt, so the Earth becomes darker, absorbs more light, and tends to get get even warmer. Conversely, when the Earth cools down, more snow and ice form, so the Earth becomes lighter, absorbs less light, and tends to get even cooler.

• Carbon dioxide solubility feedback. Cold water can hold more carbon dioxide than warm water: that’s why opening a warm can of soda can be so explosive. So, when the Earth’s oceans warm up, they release carbon dioxide. But carbon dioxide is a greenhouse gas, which causes additional warming. Conversely, when the oceaans cool down, they absorb more carbon dioxide, so the greenhouse effect becomes weaker, which tends to cool things down.

Of course, there are also negative feedbacks: otherwise the climate would be utterly unstable! There are also complicated feedbacks whose overall effect is harder to evaluate:

• Planck feedback. A hotter world radiates more heat, which cools it down. This is the big negative feedback that keeps all the positive feedbacks from making the Earth insanely hot or insanely cold.

• Cloud feedback. A warmer Earth has more clouds, which reflect more light but also increase the greenhouse effect.

• Lapse rate feedback. An increased greenhouse effect changes the vertical temperature profile of the atmosphere, which has effects of its own—but this works differently near the poles and near the equator.

See week302 for more on feedbacks and how big they’re likely to be.

On top of all these subtleties, any proposed solution to the puzzle of glacial cycles needs to keep a few other things in mind, too:

• A really good theory will explain, not just why we have glacial cycles now, but why we didn’t have them earlier. As I explained in week317, they got started around 5 million years ago, became much colder around 2 million years ago, and switched from a roughly 41,000 year cycle to a roughly 100,000 year cycle around 1 million years ago.

• Say we dream up a whopping big positive feedback mechanism that does a great job of keeping the Earth warm when it’s warm and cold when it’s cold. If this effect is strong enough, the Earth may be bistable: it will have two stable states, a warm one and a cold one. Unfortunately, if the effect is too strong, it won’t be easy for the Earth to pop back and forth between these two states!

The classic example of a bistable system is a switch—say for an electric light. When the light is on it stays on; when the light is off it stays off. If you touch the switch very gently, nothing will happen. But if you push on it hard enough, it will suddenly pop from on to off, or vice versa.

If we’re trying to model the glacial cycles using this idea, we need the switch to have a fairly dramatic effect, yet still be responsive to a fairly gentle touch. For this to work we need enough positive feedback… but not too much.

(We could also try a different idea: maybe the Earth keeps itself in its icy glacial state, or its warm interglacial state, using some mechanism that gradually uses something up. Then, when the Earth runs out of this stuff, whatever it is, the climate can easily flip to the other state.)

We must always remember that to a good approximation, the total amount of sunlight hitting the Earth each year does not change as the Earth’s orbit changes in the so-called ‘Milankovich cycles’ that seem to be causing the ice ages. I explained why last time. What changes is not the total amount of sunlight, but something much subtler: the amount of sunlight at particular latitudes in particular seasons! In particular, Milankovitch claimed, and most scientists believe, that the Earth tends to get cold when there’s little sunlight hitting the far northern latitudes in summer.

For these and other reasons, any solution to the ice age puzzle is bound to be subtle. Instead of diving straight into this complicated morass, let’s try something much simpler. Let’s just think about how the ice albedo effect could, in theory, make the Earth bistable.

To do this, let’s look at the very simplest model in this great not-yet-published book:

• Gerald R. North, Simple Models of Global Climate.

This is a zero-dimensional energy balance model, meaning that it only involves the average temperature of the earth, the average solar radiation coming in, and the average infrared radiation going out.

The average temperature will be measured in Celsius. We’ll assume the Earth radiates power square meter equal to

where watts/meter² and watts/meter² per degree Celsius. This is a linear approximation taken from satellite data on our Earth. In reality, the power emitted grows faster than linearly with temperature.

We’ll assume the Earth absorbs solar energy power per square meter equal to

Here:

• is the average insolation: that is, the amount of solar power per square meter hitting the top of the Earth’s atmosphere, averaged over location and time of year. In reality is about 341.5 watts/meter². This is one quarter of the solar constant, meaning the solar power per square meter that would hit a panel hovering in space above the Earth’s atmosphere and facing directly at the Sun. (Why a quarter? That’s a nice geometry puzzle: we worked it out at the Azimuth Blog once.)

• is the coalbedo: the fraction of solar power that gets absorbed. The coalbedo depends on the temperature; we’ll have to say how.

Given all this, we get

where is Earth’s heat capacity in joules per degree per square meter. Of course this is a funny thing, because heat energy is stored not only at the surface but also in the air and/or water, and the details vary a lot depending on where we are. But if we consider a uniform planet with dry air and no ocean, North says we may roughly take equal to about half the heat capacity at constant pressure of the column of dry air over a square meter, namely 5 million joules per degree Celsius.

The easiest thing to do is find equilibrium solutions, meaning solutions where so that

Now doesn’t matter anymore! We’d like to solve for as a function of the insolation but it’s easier to solve for as a function of :

To go further, we need to guess some formula for the coalbedo The coalbedo, remember, is the fraction of sunlight that gets absorbed when it hits the Earth. It’s 1 minus the albedo, which is the fraction that gets reflected. Here’s a little chart of albedos:

If you get mixed up between albedo and coalbedo, just remember: coal has a high coalbedo.

Since we’re trying to keep things very simple right not, not model nature in all its glorious complexity, let’s just say the average albedo of the Earth is 0.65 when it’s very cold and there’s lots of snow. So, let

be the ‘icy’ coalbedo, good for very low temperatures. Similarly, let’s say the average albedo drops to 0.3 when its very hot and the Earth is darker. So, let

be the ‘ice-free’ coalbedo, good for high temperatures when the Earth is darker.

Then, we need a function of temperature that interpolates between and Let’s try this:

If you’re not a fan of the hyperbolic tangent function this may seem scary. But don’t be intimidated!

The function is just a function that goes smoothly from 0 at low temperatures to 1 at high temperatures. This ensures that the coalbedo is near its icy value at low temperatures, and near its ice-free value at high temperatures. But the fun part here is a parameter that says how rapidly the coalbedo rises as the Earth gets warmer. Depending on this, we’ll get different effects!

The function rises fastest at since that’s where has the biggest slope. We’re just lucky that in Celsius is the melting point of ice, so this makes a bit of sense.

Now Allan Erskine’s programming magic comes into play! Unfortunately it doesn’t work on this blog—yet!—so please hop over to the version of this article on my website to see it in action.

You can slide a slider to adjust the parameter to various values between 0 and 1.

You can then see how the coalbedo changes as a function of the temperature In this graph the temperature ranges from -50 °C and 50 °C; the graph depends on what value of you choose with slider.

You can also see how the insolation required to yield a given temperature between -50 °C and 50 °C:

It’s easiest to solve for in terms of But it’s more intuitive to flip this graph over and see what equilibrium temperatures are allowed for a given insolation between 200 and 500 watts per square mater.

The exciting thing is that when gets big enough, three different temperatures are compatible with the same amount of insolation! This means the Earth can be hot, cold or something intermediate even when the amount of sunlight hitting it is fixed. The intermediate state is unstable, it turns out. Only the hot and cold states are stable. So, we say the Earth is bistable in this simplified model.

Can you see how big needs to be for this bistability to kick in? It’s certainly there when since then we get a graph like this:

When the insolation is less than about 385 W/m² there’s only a cold state. When it hits 385 W/m², as shown by the green line, suddenly there are two possible temperatures: a cold one and a much hotter one. When the insolation is higher, as shown by the black line, there are three possible temperatures: a cold one, and unstable intermediate one, and a hot one. And when the insolation gets above 465 W/m², as shown by the red line, there’s only a hot state!

Why is the intermediate state unstable when it exists? Why are the other two equilibria stable? To answer these questions, we’d need to go back and study the original equation:

and see what happens when we push slightly away from one of its equilibrium values. That’s really fun, but we won’t do it today. Instead, let’s draw some conclusions from what we’ve just seen. There are at least three morals: a mathematical moral, a climate science model, and a software moral.

Mathematically, this model illustrates catastrophe theory. As we slowly turn up we get different curves showing how temperature is a function of insolation… until suddenly the curve isn’t the graph of a function anymore: it becomes infinitely steep at one point! After that, we get bistability:

This is called a cusp catastrophe, and you can visualize these curves as slices of a surface in 3d, which looks roughly like this picture:

from here:

• Wolfram Mathworld, Cusp catastrophe. (Includes Mathematica package.)

The cusp catastrophe is ‘structurally stable’, meaning that small perturbations don’t change its qualitative behavior. This concept is made precise in catastrophe theory. It’s a useful concept, because it focuses our attention on robust features of models: features that don’t go away if the model is slightly wrong, as it always is.

As far as climate science goes, one moral is that it pays to spend some time making sure we understand simple models before we dive into more complicated ones. Right now we’re looking at a very simple model, but we’re already seeing some interesting phenomena. The kind of model we’re looking at now is called a Budyko-Sellers model. These have been studied since the late 1960’s:

• M. I. Budyko, On the origin of glacial epochs (in Russian), Meteor. Gidrol. 2 (1968), 3-8.

• M. I. Budyko, The effect of solar radiation variations on the climate of the earth, Tellus 21 (1969), 611-619.

• William D. Sellers, A global climatic model based on the energy balance of the earth-atmosphere system, J. Appl. Meteor. 8 (1969), 392-400.

• Carl Crafoord and Erland Källén, A note on the condition for existence of more than one steady state solution in Budyko-Sellers type models, J. Atmos. Sci. 35 (1978), 1123-1125.

• Gerald R. North, David Pollard and Bruce Wielicki, Variational formulation of Budyko-Sellers climate models, J. Atmos. Sci. 36 (1979), 255-259.

It also pays to compare our models to reality! For example, the graphs we’ve seen show some remarkably hot and cold temperatures for the Earth. That’s a bit unnerving. Let’s investigate. Suppose we set on our slider. Then the coalbedo of the Earth becomes independent of temperature: it’s 0.525, halfway between its icy and ice-free values. Then, when the insolation takes its actual value of 342.5 watts per square meter, the model says the Earth’s temperature is very chilly: about -20 °C!

Does that mean the model is fundamentally flawed? Maybe not! After all, it’s based on very light-colored Earth. Suppose we use the actual albedo of the Earth. Of course that’s hard to define, much less determine. But let’s just look up some average value of the Earth’s albedo: supposedly it’s about 0.3. That gives a coalbedo of If we plug that in our formula:

we get 11 °C. That’s not too far from the Earth’s actual average temperature, namely about 15 ° C. So the chilly temperature of -20 °C seems to come from an Earth that’s a lot lighter in color than ours.

Our model includes the greenhouse effect, since the coeficients and were determined by satellite measurements of how much radiation actually escapes the Earth’s atmosphere and shoots out into space. As a further check to our model, we can look at an even simpler zero-dimensional energy balance model: a completely black Earth with no greenhouse effect. Another member of the Azimuth Project has written about this:

• Tim van Beek, Putting the Earth in a box, Azimuth, 19 June 2011.

• Tim van Beek, A quantum of warmth, Azimuth, 2 July 2011.

As he explains, this model gives the Earth a temperature of 6 °C. He also shows that in this model, lowering the albedo to a realistic value of 0.3 lowers the temperature to a chilly -18 ° C. To get from that to something like our Earth, we must take the greenhouse effect into account.

This sort of fiddling around is the sort of thing we must do to study the flaws and virtues of a climate model. Of course, any realistic climate model is vastly more sophisticated than the little toy we’ve been looking at, so the ‘fiddling around’ must also be more sophisticated. With a more sophisticated model, we can also be more demanding. For example, when I said 11 °C is “is not too far from the Earth’s actual average temperature, namely about 15 ° C”, I was being very blasé about what’s actually a big discrepancy. I only took that attitude because the calculations we’re doing now are very preliminary.

Finally, here’s what Allan has to say about the software you’ve just seen, and some fancier software you’ll see in forthcoming weeks:

Your original question in the Azimuth Forum was “What’s the easiest way to write simple programs of this sort that could be accessed and operated by clueless people online?” A “simple program” for the climate model you proposed needed two elements: a means to solve the ODE (ordinary differential equation) describing the model, and a means to interact with and visualize the results for the (clearly) “clueless people online”.

Some good suggestions were made by members of the forum:

• use a full-fledged numerical computing package such as Sage or Matlab which come loaded to the teeth with ODE solvers and interactive charting;

• use a full-featured programming language like Java which has libraries available for ode solving and charting, and which can be packaged as an applet for the web;

• do all the computation and visualization ourselves in Javascript.

While the first two suggestions were superior for computing the ODE solutions, I knew from bitter experience (as a software developer) that the truly clueless people were us bold forum members engaged in this new online enterprise: none of us were experts in this interactive/online math thing, and programming new software is almost always harder than you expect it to be.

Then actually releasing new software is harder still! Especially to an audience as large as your readership. To come up an interactive solution that would work on many different computers/browsers, the most mundane and pedestrian suggestion of “do it all ourselves in Javascript and have them run it in the browser” was also the most likely to be a success.

The issue with Javascript was that not many people use it for numerical computation, and I was down on our chances of success until Staffan pointed out the excellent JSXGraph software. JSXGraph has many examples available to get up and running, has an ODE solver, and after a copy/paste or two and some tweaking on my part we were all set.

The true vindication for going all-Javascript though was that you were subsequently able to do some copy/pasting of your own directly into TWF without any servers needing configured etc., or even any help from me! The graphs ought to be viewable by your readership for as long as browsers support Javascript (a sign of a good software release is that you don’t have to think about it afterwards).

There are some improvements I would make to how we handle future projects which we have discussed in the Forum. Foremost, using Javascript to do all our numerical work is not going to attract the best and brightest minds from the forum (or elsewhere) to help with subsequent models. My personal hope is that we allow all the numerical work to be done in whatever language people feel productive with, and that we come up with a slick way for you to embed and interact with just the data from these models in your webpages. Glyn Adgie and Jim Stuttard seem to have some great momentum in this direction.

Or perhaps creating and editing interactive math online will eventually become as easy as wiki pages are today—I know Staffan had said the Sage developers were looking to make their online workbooks more interactive. Also the bright folks behind the new Julia language are discussing ways to run (and presumably interact with) Julia in the cloud. So perhaps we should just have dragged our feet on this project for a few years for all this cool stuff to help us out! (And let’s wait for the Singularity while we’re at it.)

No, let’s not! I hope you programmers out there can help us find good solutions to the problems Allan faced. And I hope some of you actually join the team.

By the way, Allan has a somewhat spiffier version of the same Budyko-Sellers model here.

For discussions of this issue of This Week’s Finds visit my blog, Azimuth. And if you want to get involved in creating online climate models, contact me and/or join the Azimuth Forum.

---

Thus, the present thermal regime and glaciations of the Earth prove to be characterized by high instability. Comparatively small changes of radiation—only by 1.0-1.5%—are sufficient for the development of ice cover on the land and oceans that reaches temperate latitudes. – M. I. Budyko