## Abstract

The Earth's main magnetic field is generally believed to be due to a self-exciting dynamo process in the Earth's fluid outer core. A variety of antidynamo theorems exist that set conditions under which a magnetic field cannot be indefinitely maintained by dynamo action against ohmic decay. One such theorem, the *Planar Velocity Antidynamo Theorem*, precludes field maintenance when the flow is everywhere parallel to some plane, e.g. the equatorial plane. This paper shows that the proof of the Planar Velocity Theorem fails when the flow is confined to a sphere, due to diffusive coupling at the boundary. Then, the theorem reverts to a conjecture. There is a need to either prove the conjecture, or find a functioning planar velocity dynamo. To the latter end, this paper formulates the toroidal–poloidal spectral form of the magnetic induction equation for planar flows, as a basis for a numerical investigation. We have thereby determined magnetic field growth rates associated with various planar flows in spheres. For most flows, the induced magnetic field decays with time, supporting a planar velocity antidynamo conjecture for a spherical conducting fluid. However, one flow is exceptional, indicating that magnetic field growth can occur. We also re-examine some classical kinematic dynamo models, converting the flows where possible to planar flows. For the flow of Pekeris *et al*. (Pekeris, C. L., Accad, Y. & Shkoller, B. 1973 Kinematic dynamos and the Earth's magnetic field. *Phil. Trans. R. Soc. A* **275**, 425–461), this conversion dramatically reduces the critical magnetic Reynolds number.

## 1. Introduction

The Earth's main magnetic field exhibits behaviour on time-scales ranging from seconds to millennia. Characteristic features over longer time-scales of years include secular variation, polar wandering, polarity reversals and long-term maintenance. The self-exciting dynamo hypothesis, that such large-scale behaviour is due to magnetohydrodynamic interactions in the Earth's fluid core, was advanced by Larmor (1919). An ionospheric disturbance dynamo hypothesis had been previously proposed by Stewart (1883) to explain daily magnetic variations. The feasibility of self-exciting dynamos has been well established by numerous kinematic and dynamic investigations, and by the establishment of dynamical benchmark computer codes (Dudley & James 1989; Glatzmaier 1997; Christensen *et al*. 2001).

However, from almost the beginning of dynamo theory, so-called antidynamo theorems have hindered the construction of simple models based on laminar flows. Such theorems rule out magnetic field maintenance when certain simplifying geometric constraints are imposed either on the flow or the magnetic field.

One such theorem is the *Toroidal Velocity Theorem* (Bullard & Gellman 1954) which rules out field maintenance by a purely toroidal fluid velocity , where *f* is some defining scalar and * r* is the radius vector. An analogous result exists in Cartesian coordinates. The

*Planar Velocity Theorem*(PVT) precludes maintenance of

*by incompressible flow with no*

**B***z*-component, i.e. a

*planar flow*(1.1)The conditions under which the PVT holds are central to this paper. In §2, it will be argued that a proof for the PVT has not been established when the conducting fluid volume is a sphere, the context of primary interest in self-gravitating, rapidly rotating bodies such as the Earth. In §3, we derive toroidal–poloidal spectral equations, to be used later in §5 for a numerical investigation of specific planar velocity models.

## 2. The Planar Velocity Theorem

We consider two cases: (i) where the conducting fluid occupies all space , and (ii) where the fluid occupies a spherical volume *V*, of radius *a* and surface , surrounded by free space . The fluid is assumed to have constant electric conductivity *σ*, and magnetic diffusivity where *μ* is the free space permeability. The magnetic field satisfies everywhere, and inside the conductor evolves in time *t* according to(2.1)where is the prescribed fluid velocity. Suppose firstly that the fluid occupies , and that both and decrease sufficiently fast as (as detailed below). For incompressible planar , and the *z*-component of (2.1) is(2.2)Multiplying by , integrating by parts and applying the divergence theorem yields(2.3)where is the outward directed surface differential. Letting gives(2.4)assuming , are as . From (2.4), we conclude that either , or decays strictly monotonically with *t*. This is the crux of the PVT.

To show decay of the component , introduce scalars , , such that(2.5)Letting , implies and hence(2.6)where the integral is over the plane . From (2.6), one infers that or decays as , following the behaviour of . So, asymptotically as , . For this , uncurling (2.1) yields(2.7)where . Let where the average is over the disc *A* of radius , at fixed *z*. Averaging (2.7), and applying the divergence theorem in the plane, yields(2.8)where the line integral is around the perimeter *C* of *A*. Since can be adjusted by adding a function of *z* without affecting , we may take . Letting in (2.8) then shows that . Assuming as , manipulations of (2.7) analogous to those leading to (2.3) implyThus, or decays strictly monotonically. Hence, from (2.5), and decay, completing the proof.

The above proof in essence parallels that of Zel'dovich & Ruzmaikin (1980). And the PVT is commonly thought to apply when the conducting fluid is confined to a finite volume (Moffatt 1978). However, here we argue the contrary. If the fluid occupies finite *V* surrounded by insulating then(2.9)with continuous across and as . Let represent values on the inside and outside of , and represent the jump across . If the boundary is fixed, at *r*=*a*. In this case, (2.3) reduces to(2.10)

Since the current sources are confined to *V*, as . Applying the divergence theorem, the *z*-component of the curl of (2.9) yields(2.11)

Combining (2.10) and (2.11), and using the continuity of across , leads to(2.12)In general, is discontinuous across . The first term on the right side of (2.12) therefore prevents a conclusion about the decay of . The initial-value problem defined by (2.1), (2.9) and , with initial value , has a unique solution if is continuous across . To impose continuity on would in general result in that initial-value problem being over-determined, with no solution.

Antidynamo theorems generally rely on the induction equation (2.1) decoupling somehow. In the PVT, the induction equation for has decoupled, so that there is no generation of from the horizontal component . When the conductor occupies , evolves to match the conditions as . So, , if not identically zero, decays. However, the decoupling of the induction equation depends on , and *η* being constant in . Introducing an insulating exterior introduces variable *η*, diffusive coupling between and at the boundary *r*=*a*, and the possibility of being maintained. Coupling between and can be seen explicitly, by noting that implies . HenceAlternatively, can be written in terms of the electric current density . Since , and ,So the coupling term on the right side of (2.12) can be expressed in various ways:(2.13)The presence of the boundary *r*=1 must promote the existence of tangential currents, generally with a component, and the coupling term (2.13) may be significant. However, on is an integrated contribution from throughout *V*. So, the sign of (2.13) is unknown. If the sign is negative, then (2.13) may act as a generation term in (2.12). Given the lack of a proof of decay in the present case, we therefore investigated numerical solutions to see if non-decaying magnetic fields can be found. Maximum principles applied to (2.2) imply that if does not decay then its maximum must occur at the boundary *r*=1 (Sperb 1981). One might expect then that a growing field will be more likely with a flow that has more shear at the boundary.

## 3. Numerical methods

The non-dimensionalized form of the induction equation (2.1) is(3.1)Here, is the magnetic Reynolds number, based on a characteristic speed *U* and the radius *a* of the spherical conducting fluid, with *t* measured in diffusion time units . For the angular coordinates , we employ a spherical harmonic spectral representation as described in §3*a*. The size of the resulting numerical problem is greatly reduced using the decoupling rules described in §3*b*. For solution of the resulting radial equations we used finite differences as described in §3*c*.

### (a) The spherical harmonic spectral equations

We represent and in poloidal–toroidal formalism, using scalars *S*, *T* and *s*, *t*, respectively:(3.2)(3.3)and expand these scalars in spherical harmonics,(3.4)where and . We denote the individual vector modes of by , , and similarly for . In the above equation,(3.5)where is the Schmidt-normalized Legendre function. The overbar denotes complex conjugate, and is the Kronecker delta. Since and are real, (3.5) implies(3.6)and similarly for and .

With *m*-superscripts suppressed, and , the spectral form of (3.1) is(3.7)

(3.8)There is no interaction in (3.7), since a poloidal magnetic field cannot be created by the interaction of a toroidal flow and toroidal field. The interaction terms are given explicitly, albeit in a different formalism, by Bullard & Gellman (1954). Equations (3.7) and (3.8) are solved subject to the conditions(3.9)

(3.10)The velocity satisfies(3.11)

(3.12)Equations (3.9) and (3.11) are required for differentiability of and with respect to at *r*=0. Equation (3.10) reflects the continuity of across and the current-free nature of including the vanishing of at . Equation (3.12) applies at a rigid boundary. Some flows also satisfy the no-slip conditions(3.13)

### (b) Selection rules and field decoupling

The interaction terms on the right side of (3.7) and (3.8) are proportional to the Adams *A* and Elsasser *E* integrals (James 1973), which impose the following selection rules.

, and depend on

*A*and are zero unless is even, and ., , and depend on

*E*and are zero unless is odd, and .are all zero unless .

These rules can be used to greatly reduce the size of the computational problem. For the present purposes, it is convenient to classify fields and as having dipole (D) or quadrupole (Q) parity if their poloidal and toroidal components belong to one of the following sets:(3.14)Here, denotes the presence of both and , and similarly for *T*. The D chain comprises *S*-harmonics with odd *n*, and *T*-harmonics with even *n*; the Q chain comprises *S*-harmonics with even *n*, and *T*-harmonics with odd *n*.

Rules SR1 and SR2 imply that if has Q-parity, then decouples into independent chains, a D chain, and a Q chain. This decoupling greatly reduces the size of the matrices in the computations described in §3*c*. The decoupling is verified by considering the interactions shown in table 1, where *bold* entries show the type of harmonic produced by the various interactions. For example, must satisfy SR1, so is zero unless is even. Thus, a Q-parity interacts with a D-parity from chain D to produce a D-parity also in chain D. However, a D-parity interacts with a D-parity from chain D to produce a Q-parity that is not in chain D. Similar considerations lead to the chain-preserving interactions shown asterisked in table 1, and show that Q-parity preserve chains D and Q, whereas D-parity destroy chains D and Q.

Independently of D–Q decouplings, the field can also be decoupled with respect to the order *m*, using SR3. Suppose only contains harmonics with orders chosen from the set for some integer *k*. Then, allowing for -pairs, as required by (3.6), SR3 implies that decouples into independent chains containing harmonics with orders chosen from one of the sets , . If *k*=1 then no *m*-decoupling occurs, i.e. mode M_{01} includes all such that . However, for example, if *k*=4, there are three decoupled -chains, namely M_{04}, M_{12} (equivalent to M_{14}) and M_{24}. Combining this *m*-decoupling with the D–Q decoupling in (3.14) further significantly reduces the dimensions of the computational problem. We denote such chains as and . The various chains used in the computations for this paper are given explicitly in §5.

### (c) The radial equations

Given that is stationary, we look for solutions of the formFor fixed *R* we let denote the with maximum real part. For non-decaying modes, for some *R* greater than a critical value. Having formed the spectral equations (3.7) and (3.8), we truncate the series (3.4) at some finite degree *N*. For most of our results, we discretized the radial direction using second-order finite differences with *J* subintervals. We also used some cross checks against several other independent programs using fourth-order finite differences and Chebyshev collocation.

For a given magnetic Reynolds number *R*, and truncation levels *J* and *N*, is the solution of an eigenvalue problem(3.15)Convergence was sought by increasing the truncation levels *J*, *N*. The solution of (3.15) was straightforward except that very large *J*, *N* were required in some cases. The matrix is very sparse and banded. Moreover, the matrix band itself has many zero elements due to the selection rules given in §3*b*. In most cases, (3.15) was solved by inverse iteration (Wilkinson 1965) without pivoting, which preserves the bandwidth and allows higher truncation levels.

## 4. Planar flows in poloidal–toroidal form

For the planar flow (1.1), the streamlines are the level surfaces of the stream function *f*, lying in planes perpendicular to the *z*-axis. To apply the numerical methods of §3 we need to express in the poloidal–toroidal form (3.3), with coefficients , as in (3.4). Expanding , and dotting (1.1) with gives(4.1)whereIn the following, we often suppress the superscript *m*, letting represent , etc. Thus, (4.1) implies(4.2)To obtain the toroidal coefficient of the flow, we use the curl of (1.1):(4.3)With , the spherical harmonics defined by (3.5) satisfy the recurrence relations (Chapman & Bartels 1962):Using these, and letting , the spectral form of (4.3) is(4.4)Thus, for each , there are two toroidal coefficients:(4.5a)(4.5b)as well as the poloidal coefficient given by (4.2). So, for a flow derived from a prescribed *f*, it is straightforward to find the corresponding poloidal and toroidal coefficients. This approach is the basis of our numerical investigation in §5*a*.

The inverse problem is also of interest. That is, is it possible to convert a given flow to a planar flow, by addition of components defined via (4.2), (4.5*a*) and (4.5*b*)? It is convenient to use the terminology ‘planarizing’ for such a construction. Clearly, if , then for given , one can find from (4.2), and construct from (4.5*a*) and (4.5*b*). The combination is then planar, and has been planarized. However, to planarize a given component, there are two paths one might follow. One has to either (i) derive from using (4.5*a*), or (ii) derive from using (4.5*b*). In (i), one generates planarizing and from (4.2) and (4.5*b*). In (ii), one generates planarizing and from (4.2) and (4.5*a*). For given , the *f*-solutions of (4.5*a*) and (4.5*b*) satisfying the rigid boundary condition (3.12) at *r*=1 are(4.6a)

(4.6b)In case (i), must also satisfy the consistency condition(4.7)in order for given by (4.6*a*) to satisfy the differentiability condition (3.11). If does satisfy (4.7), then the corresponding planar flow is ; otherwise cannot be planarized by method (i). For method (ii), given by (4.6*b*) satisfies (3.11), and for , a planarized flow of the form can be constructed. However, if , the flow cannot be planarized by method (ii) since then . Later, in §5*b*, we will consider the possibility of planarizing some classic flows, where is defined by a given set of , coefficients.

## 5. Results for various planar flows

In this section, we solve the kinematic dynamo problem comprising (3.7)–(3.10) for various planar flows. First, we consider flows with a single harmonic degree, then some flows of historical interest.

### (a) Single harmonic flows

Throughout the following we use for the flow, and , for the magnetic field, as done earlier in the spectral equations (3.7) and (3.8). With , we considered the stream function with(5.1)The associated planar flow is(5.2)where are derived from via (4.2), (4.5*a*) and (4.5*b*). This flow satisfies the differentiability condition (3.11), the fixed boundary condition (3.12) and, if *p*>1, the no-slip conditions (3.13).

Computer limitations and convergence difficulties restricted our investigation to relatively low . To allow decoupling as in (3.14), we mainly considered , where the flow based on (5.1) has Q-parity, and preserves D-parity magnetic chains as discussed in §3*b*. Allowing , , and D-parity magnetic field chains, led to 36 models overall. For brevity, we will refer to a model with flow parameters and magnetic field chain DM_{12} by the acronym Y42DM12, and if *p*=1 then by p1Y42DM12, and similarly for other models. Some were found to be complex, representing oscillatory magnetic field modes. In this report, we are primarily interested in and determining whether the field grows or not.

Streamlines for various flows at height are shown in figure 1 for *p*=1. Increasing *m* increases the horizontal cellular structure. The vertical cell structure depends on both *n* and *m*, and the direction of circulation depends on *z* and . In particular, the equator *z*=0 is stagnant if *n*−*m* is odd. Various contours of in the *xz*-plane are shown in figure 2. Other flows defined by (5.1) exhibit similar patterns, with flows for increasing *p* having decreasing shear near *r*=1.

Generally, we tracked out from its known free-decay value (where *ν* is a spherical Bessel function root, usually *π*), to at least *R*=200. For models defined by (1.1) and (5.1), , so reversed flows need not be considered. Qualitatively, our results fell into four categories: (i) axisymmetric flows with decaying magnetic fields, (ii) non-axisymmetric flows with decaying fields, (iii) a non-axisymmetric flow with a growing field and (iv) models where insufficient convergence with increasing *J*, *N* was obtained. The convergent results for (i)–(iii) are discussed in the following sections.

#### (i) Axisymmetric flows

For axisymmetric flows, *m*=0 and hence by (4.2). The streamlines are just circles in planes of constant *z* exemplified by frame 1 of figure 1. There is full decoupling with respect to the order of the magnetic field. So, we considered cases for the dipole magnetic field chainsSuch models are covered by the Toroidal Velocity Theorem of §1, the magnetic field having a theoretical decay rate upper bound from variational calculus. So, such flows cannot maintain a magnetic field, but they were used as one of our various computational checks.

#### (ii) Non-axisymmetric flows

Typical multicellular streamlines are shown in figure 1. For flow parameters , the various possible decoupled dipole chains are listed explicitly below, the list being shortened by the equivalences , .

The chain DM_{01} has no decoupling with respect to the magnetic field parameter , but all other chains shown above give significant computational saving due to -decoupling on top of D-parity decoupling. All feasible combinations of the above chains with flows derived from (5.1), for , and , were considered. Convergence was considerably more difficult to obtain for *p*=1, than for , especially for higher *n*.

The exceptional model p1Y22DM12 is discussed further on. For all other models where convergence was evident, the decay was found to be faster than free decay for the entire interval , and the decay rates were faster for larger *R*, and smaller *p*. While there are some differences in detail, the decay profiles are similar to those shown in figure 3*a* for model Y21DM01, and indicate support for a PVT. A mode crossing is evident in figure 3*b* for *p*=1 at , above which a steady magnetic field mode is dominant. Similar mode crossings occasionally occurred in the other models.

#### (iii) The exceptional case

The profiles for model Y22DM12 are shown in figure 3*b*. In all three cases , is real. For , good convergence was obtained and for all *R*. The decay profile for p2Y22DM12 does tend slightly upwards for larger *R*, but determining whether that trend persists and gives field growth for requires truncation levels well above our present capacity. The case *p*=1 stands out in various respects. For *p*=1, convergence was much harder to obtain, the growth rate estimates only starting to settle down above and . Figure 3*b* shows that the *p*=1 profile differs markedly from the reasonably typical decay profiles for in figure 3*a*. The *p*=1 profile consistently increases with *R* for the wide range of high truncation levels .

The magnetic field vector is more sensitive to truncation levels, so in figure 4 we show the projections of the components in the plane *ϕ*=0. There is good qualitative agreement between the different truncation levels. The differences in and between the truncations levels occur in regions where those components are small. It is notable that is confined to a band near the equator, to a band near the polar axis and to a band near the polar axis, with a weaker band near the boundary *r*=1. And is highly concentrated near the boundary, consistent with the maximum principle argument in §2.

While the level of convergence is not complete, the consistently growing profiles in figure 3, and the field projections in figure 4 are strong prima facie evidence that the *p*=1 flow supports magnetic field growth for , and that the PVT is therefore invalid for bounded conducting regions.

### (b) Planarizing some classic flows

In this section, we consider planarizing three classic flow models using the method described in §4. In all three cases, it is not possible to fully planarize, but for two different reasons. In one case, part planarization destroys the pre-existing field growth of the original model. In another case, case part planarization enhances the dynamo action. Details of the original non-planarized flows are given by Dudley & James (1989).

#### (i) Bullard & Gellman flow

In the present formalism, the classic flow of Bullard & Gellman (BG, 1954) is(5.3)whereThis motion satisfies the fixed boundary condition (3.12) and no-slip conditions (3.13), but not the differentiability condition (3.11). A planarized version of the BG flow is(5.4)From (4.2) and (4.5*b*),(5.5)where . We report here on the case *ϵ*=5 with DM_{02} magnetic decoupling. All numerical evidence suggests that the original BG flow fails to function as a dynamo, and our work suggests that the planarized version (5.4) also fails. Indeed, as figure 5*a* shows, the planar -profile is much flatter than the original, being closer to free-decay, and decaying less rapidly. Furthermore, the planar is purely real, whereas the BG is complex for .

#### (ii) Kumar & Roberts flow

The flow (5.6) was one of the first shown to function as a dynamo (Kumar & Roberts (KR) 1975). The differences from the BG flow (5.3) are the addition of the large-scale meridional circulation , and the dependence of on both and . It is a convenient model to consider because of its fast numerical convergence. The main flow considered by KR was, in the present formalism,(5.6)whereThis motion satisfies the fixed boundary condition (3.12) and no-slip conditions (3.13), but not the differentiability condition (3.11). We used , , *k*=3 and DM_{02} decoupling, a configuration shown by KR to produce dynamo action. This flow cannot be fully planarized using (4.2), since *m*=0 for the meridional circulation . However, the component is planarized by the addition of a component as done for the BG flow in §5*b*(i). The result is a meridional flow superimposed on a planar flow:(5.7)Here, is again given by (5.5) with . As shown by figure 5*b*, this part planarization greatly flattens the -profile, and removes the possibility of dynamo action. Although not shown on the graph, this flat decay profile persisted out to . Also, the -profile for the partly planarized flow is purely real, whereas the original KR profile shows that an oscillatory mode dominates for . Even without meridional circulation, i.e. with , the KR flow (5.6) still acts as a dynamo, although *R*>20 000 is required (Sarson & Gubbins 1996). However, for the fully planar flow (5.7) with , we found that differed only slightly from the free-decay value over the entire interval .

It is perhaps of some interest to note that a similar construction, i.e. addition of a component to the KR flow, has previously been used by Sarson (2003) to convert the KR flow to a flow interpretable as a thermal wind.

### (c) Pekeris, Accad & Shkoller flow

Lastly, we consider the possibility of planarizing the flow of Pekeris, Accad & Shkoller (PAS, 1973), another early successful kinematic dynamo. This is an example where it is not possible to planarize the flow even though it has no axisymmetric meridional component. The motion of PAS is(5.8)whereHere, , and is a positive root of the spherical Bessel functionThis motion satisfies the fixed boundary condition (3.12) and the differentiability conditon (3.11), but not the no-slip conditions (3.13). By (4.2), in (5.8) is easily planarized by usingobtaining from (4.5*b*), and using . By (4.5*a*), in (5.8) is planarized by finding fromHowever, since does not satisfy the consistency condition (4.7), the resulting cannot satisfy both the differentiability condition (3.11) and the fixed boundary condition (3.12). Hence, PAS cannot be fully planarized by this method. We thus left unchanged and considered the effect of just planarizing the component. The resulting flow is then the superposition of a purely planar component and a sectorial component :(5.9)Neither the original PAS flow (5.8), nor the partly planarized version (5.9) support D or Q magnetic decoupling. Here, we report on the magnetic chain M_{02}, shown by PAS to be maintainable by the flow (5.8). As shown by our numerical results for Λ=12.3229 in figure 5*c*, the partly planarized flow (5.9) reduces the original PAS critical *R* from about 0.37 to about 0.26, and greatly increases the rate of growth of the magnetic field for supercritical *R*.

## 6. Conclusions

In §2, the Planar Velocity Theorem was reviewed, and it was argued that in a spherical conductor the proof fails due to diffusive coupling of the horizontal and vertical fields and at the boundary. In §4, the relations between the stream function for a planar flow and the spherical harmonic poloidal and toroidal spectral components were derived. The spectral equations in §3*a* were used to numerically investigate possible field maintenance by planar flows. In §5, we reported on the numerical growth rates determined for a number of single harmonic planar flows , where is the stream function. Where convergence was evident, our results generally indicate field decay, typified by the -profiles in figure 3*a*. However, for the exceptional case , the -profiles in figure 3*b* behave remarkably differently. These growing profiles are strong prima facie evidence that planar flows can maintain magnetic fields. These results and the theoretical arguments given in §2 imply that a PVT theorem does not apply to bounded conductor regions. A more definite conclusion requires considerably higher truncation levels and computing capacity, or perhaps an alternative numerical method avoiding spherical harmonics.

At the outset of this project, and given formula (2.13), it was expected that lower *p* and hence higher velocity shear near the boundary would be more effective in acting as a self-sustaining dynamo. Superficially, this has proven true in that the successful model p1Y22DM12 has stream function with *p*=1. However, in other cases, flows with higher *p* have produced slower decay rates, preventing a conclusion about effectiveness based simply on *p*. For example, in figure 3*a*, the *p*=1 flow produces the most rapidly decaying magnetic field.

For the Earth, rotation is a very significant factor. Fully dynamical geodynamo models (Glatzmaier 1997) show concentrations of planar flow, indeed zonal flow, about the inner core boundary and extending into the tangent cylinder above and below the inner core. Thus, the study of the efficacy of planar flows is of considerable ongoing interest. The resolution of the convergence difficulties referred to earlier, and the establishment of the validity or otherwise of the Planar Velocity Theorem remain the subject of ongoing work.

## Acknowledgements

One of us (A.A.B.) wishes to acknowledge the financial support from an Australian Development Scholarship at the University of Sydney, while on leave from the University of Indonesia. We also wish to acknowledge the support of the Australian Research Council via the Sydney Regional Visualisation Laboratory where some of our results were generated. We also thank the referees for their helpful comments.

## Footnotes

- Received September 12, 2005.
- Accepted January 6, 2006.

- © 2006 The Royal Society