## Abstract

The work extends the linear fields’ solution of compressible nonlinear magnetohydrodynamics (MHD) to the case where the magnetic field depends on superlinear powers of position vector, usually, but not always, expressed in Cartesian components. Implications of the resulting Lie–Taylor series expansion for physical applicability of the Dolzhansky–Kirchhoff (D–K) equations are found to be positive. It is demonstrated how resistivity may be included in the D–K model. Arguments are put forward that the D–K equations may be regarded as illustrating properties of nonlinear MHD in the same sense that the Lorenz equations inform about the onset of convective turbulence. It is suggested that the Lie–Taylor series approach may lead to valuable insights into other fluid models.

## 1. Introduction

Magnetohydrodynamics (MHD) is the study of electrically conducting fluids, where the flow is affected by the presence of a magnetic field. Its importance lies in the fact that ionized plasmas conduct electricity and may often be treated in a fluid approximation. The plasma state, as well as constituting over 99% of the matter in the observable universe, is also the state of the fuel in magnetic confinement fusion experiments that offer the potential for almost unlimited nuclear energy production. ‘Ideal’ in the MHD context implies that the fluid has no electrical resistance and normally that the flow affects the field, so that nonlinear effects are important. These must be treated because ideal MHD has important applications, namely to terrestrial confinement devices such as tokamaks, and to the generation of the solar magnetic field and its interaction with the Earth. (There is usually a small resistivity in both these applications, however, see discussion in §4.) The ‘linear’ fields approach reduces ideal MHD to a low-order set of ordinary differential equations (ODEs) capable of further simplification, so has the potential to enrich understanding of this difficult subject.

The recent work [1] showed how the equations of ideal, compressible magnetohydrodynamics may be elegantly formulated in terms of Lie derivatives, building on the work of Helmholtz, Walen and Arnold. The magnetic induction equation for compressible flow may be formulated in terms of a Lie derivative of a vector by introducing the field **B** divided by the mass density
**u**, *ρ* is mass density. The dynamical, potential vorticity equation [1] may also be put into the Lie derivative form
*p*(*ρ*) or sometimes in the isentropic approximation. Observe that the vectors which are evolved by equations (1.1) and (1.2) satisfy ∇⋅(*ρ***F**)=**0**, provided there is mass conservation and *ρ***F** is solenoidal initially.

Now, it is known since Dungey [2] that if the velocity field depends linearly on the Cartesian position vector, then compressible MHD is reducible exactly to a set of ODEs in the coefficients of the proportionality constants. (There is a much longer and complicated history regarding classical hydrodynamics which will not be discussed herein, other than to indicate that two largely complementary recent references for this material are reference [3] discussed below, and McKiver & Dritschel [4].) This ‘linear fields’ theory has been developed further as described in Arnold & Khesin [5], §I.10.C to three-dimensional MHD. The most comprehensive published treatment in the MHD context is that of Holm [6], wherein ‘linear fields’ are described using the perhaps less confusing term ‘affine motions’. Distinct from the present approach, reference [6] uses Lagrangian variables rather than Eulerian, and begins with the Lagrangian formulation of MHD rather than with equations (1.1) and (1.2). Dolzhansky [3] explains clearly how a special choice of three-dimensional vector basis for both velocity and magnetic field leads to the Kirchhoff equations, a sixth-order ODE describing the motion of an ellipsoid immersed in fluid. The difficulty with the three-dimensional vector basis is that it requires initial current distributions corresponding to the linear magnetic field that are not easy to realize in practice. It is possible to conceive that an ellipsoidal blob of uniform vorticity might somehow appear; indeed, an *elliptical* blob, embedded in a two-dimensional potential flow as was proposed by Helmholtz in 1889 as a model for a tornado [7], §159, indeed, Holm [6] refers to his solutions as ‘magnetic tornadoes’. However, it strains the imagination as to how an isolated *ellipsoidal* current distribution might be spontaneously produced.

Arter [8] pointed out that the Cartesian linear fields could be regarded as the truncation of a Taylor series expansion solution in position to first order. Thus, a problem with boundary conditions at a finite distance from the origin might be formulated, by allowing higher-order Taylor terms to help say fix the current on a flat surface at distance *λ* rather than on the problematic ellipsoidal surface. Substituting the higher-order terms in the governing partial differential equations (PDEs) leads to complicated sets of ODEs, but fortunately the use of the above Lie derivative form for MHD is especially convenient for such analysis; see §2. Section 3 first discusses implications of the results in §2d for Dolzhansky’s model and then explores important mathematical features of the Dolzhansky–Kirchhoff (D–K) equations. As might be expected, nearly all progress beyond the initial derivations relies on the assumption of solenoidal velocity field. Section 4 discusses the introduction of resistivity into the D–K model and conclusions are drawn in §5.

## 2. Lie–Taylor expansion and implications

### (a) Lie derivative expansions

Suppose that the vectors **u** and **q** have components labelled *j* and moreover that each component may be separately expressed as a Taylor series in coordinates *x*^{i}, i.e.
**q** is set equal to **q** from equation (1.1) may be written as
**v**,**w**] for general vectors. Hence, in component form

It follows that it is necessary to calculate the derivative Taylor series:

If the Taylor series representation, equations (2.2) and (2.8) and equivalents, is now substituted in equation (2.6), the term which is independent of position vector gives ODEs
*U*^{i}=0, *Q*^{j} depends on

Following this recipe produces

At *N*th order, it is apparent that coefficients with *N* suffixes evolve according to sums of nonlinear terms, each containing a total of *N*+1 suffixes. It follows that if *U*^{i}=*Q*^{i}=0, then there is no closure problem, each order varies in time depending only on itself and lower-order contributions.

Equations for the evolution of *Q*^{j}*Q*^{j}, *U*^{i}=*Q*^{i}=0,
*i* and *j* shows that the terms in *QUQ* cancel, and hence
*L*>0
*Q*^{2} (as well as *Q*) is constant in time (provided *U*^{i}=*Q*^{i}=0). The more powerful result would be that *tr*(*QQ*^{T}) is conserved, instead *tr*(*Q*^{2}) is the sum of squares of elements of *S* less the sum of squares of elements of *A*, where *S* and *A* are the symmetric and skew-symmetric parts of *Q*, respectively. Considering the separate cases *Q*=*S* and *Q*=*A* together with need for the time evolution to maintain the (skew-) symmetry, equation (2.18) implies that solutions of equation (2.12) are bounded if *Q* and *U* are both skew-symmetric. When they are both symmetric, their commutator is skew-symmetric and there is no consistent dynamic. A stronger result can be deduced by contracting equation (2.12) with *U* is skew-symmetric.

Because they are based purely on Taylor series’ manipulations, and the Lie derivative has the same form in any non-degenerate coordinate system if the vector components are ‘raised’, i.e. treated as contravariant, it should be clear that all the above results apply in an arbitrary coordinate system. The induction equation may also [9] be expressed as the vanishing of a four-dimensional Lie derivative, with implications for deriving solutions which are polynomial in time.

### (b) Scalar transport equation

Although not strictly needed in the current work, the compressible MHD equations are completed by scalar transport equations for quantities such as internal energy, and at minimum the mass density *ρ*. For consistency with the development in the previous section, it is necessary to introduce the point mass *ϱ* evolves as
*ρ***F** is also frame independent in this same sense, becoming

Taylor-expanding the point mass
*U*^{j}=0, or for all higher derivatives of *ϱ* to vanish, and for higher-order derivatives of *U* also to be zero.

In the cases considered, *ρ* replacing *ϱ*. Moreover, the model mainly considered is inherently incompressible, so *ρ* is also time- and position-independent.

### (c) Vorticity and current

In a general, non-orthogonal coordinate system with coordinates *x*^{j}, the curl operator relating velocity to vorticity is such that
*e*^{ijk}=*e*_{ijk} is the permutation symbol (*e*_{123}=−*e*_{132}=1, *e*_{112}=0 etc.), *g*_{ij} is the metric tensor and

Supposing *g*_{sq} to be constant, equivalently assuming an affine transformation, the second two terms vanish; then if the linear fields’ assumption is made for the velocity

Assuming incompressible fluid, write *ω*^{j} for *Q*^{j} in equation (2.9), so that this represents the vorticity evolution equation (1.2) without the forcing terms. It follows that the background flow is immaterial for a purely linear velocity field since then

There is a problem for the linear fields’ approach, unless the momentum equation is explicitly introduced, because the vorticity equation (2.28) represents an evolution equation for the three components *ω*^{j}, yet there *ϖ*_{k}. If ℓ_{ij} is the transformation matrix of coordinates, then identifying *U*_{qr}, the definition is completed as
^{T}=*I* the identity).

At this juncture, note that if analogously *B*=ℓ*K*ℓ^{−1}, then the matrix commutator
*Q*=*B* evolves as the commutator equation (2.19), then

### (d) Finite domain considerations

The difficulty in physically interpreting the outcome of a ‘linear fields’ model is that because the fields increase linearly in an unbounded domain, they have infinite energy at all times. (Beware that this also implies the Helmholtz decomposition of each field into gradient and curl is not unique.) Thus, it is not unexpected that linear fields solutions in general have finite time singularities [10], particularly when the linear fields represent inflow boundary conditions. Imshennik & Syrovatskii [10] interpret these singularities as implying current sheet formation, and go on to discuss how they might occur in a more realistic situation where the fields are linear only in a bounded region, inferring that the singularities require input of significant external energy. Hence singular linear fields cannot be regarded as self-consistent *local* models, leading to the emphasis of this work on ensuring bounded solutions. However, as discussed in §1, there are physical difficulties regarding the assumption of the existence of current blobs needed to ensure a bounded MHD problem in general.

The value of the Taylor series approach is that the higher-order terms allow for a more realistic current distribution. The question is to what extent is their presence consistent with the simple ‘linear fields’ model.

The first point to note is that the equations for evolution of the hierarchy of derivatives in §2a may be formally, consistently ordered if the second- and higher-order derivatives are supposed to be smaller than first-order derivatives by a factor of order *ϵ*≪1. Inspection of equation (2.1) shows that the second-order term rises up to equal the first when *r*=∥*x*∥/*a*_{0} measures distance from the origin scaled by lengthscale *a*_{0}. At *x*≈*a*_{0}, inflow or indeed outflow boundary conditions with *u*∝±*x* could be imposed.

The simplest example of a function *f* satisfying *a*_{0} needs to be smaller than the domain size *λ*, and in practice, it will be set by the initial ratio of function to first derivative *B*^{j}≠0.

There remains the question as to whether the ordering remains consistent under time evolution, that is to say whether equation (2.15) implies that second-order field derivatives also do not grow. In the case where *U*^{i}=0 contracting *U* is contracted with a symmetric tensor consisting of the product of *Q* with itself, hence each vanishes separately. This statement does not depend on the coordinate system used hence taking *Q*=*B*, the sum of the squares of the *V* was antisymmetric). Given the same result at the end of §2a for first-order *Q* derivatives, it should be evident that a similar analysis could be conducted at any higher order.

## 3. Case study: Dolzhansky–Kirchhoff equations

### (a) Derivation

The Kirchhoff equations for *incompressible* MHD after Dolzhansky [3] follow upon assuming that the transformation of coordinates applied to an antisymmetric matrix representation of the velocity gradient matrix *U* is a simple anisotropic scaling
_{ij} introduced in §2c is diagonal. The convention is adopted that suffixes on *x*, as distinct from those on fields *u*^{i} and *ϱ*, denote Cartesian coordinates.

Without the anisotropy induced by this scaling, it helps to note that Dolzhansky has introduced a basis for the fields which takes the form
*i*, **e**_{i}⋅**x**=0 and the linear fields flow expressed in terms of this basis, viz.,
*ϖ*′′_{k} vary only with time, have streamlines and indeed streaklines which are confined to spherical surfaces. Under the transformation equation (3.1), and indeed, any affine scaling, spheres become ellipsoids and the basis becomes non-orthogonal, but the same local confinement property applies. Rather less satisfactorily from the physical point-of-view as mentioned in §1, although justification was attempted in §2d, the magnetic field has to have the same basis and therefore is forced to have this local confinement property, viz.,

The derivation of the Kirchhoff equations from the Lie–Taylor series approach of §2 proceeds by introducing the transformed vector ** Ω**=ℓ

^{−1}

**, in components**

*ω**Ω*

^{i}=ℓ

_{ji}

*ω*

^{j}, so that equation (2.9) for

*Q*

^{i}as the vorticity (i.e. with

*Q*

^{i}=

*ω*

^{i}) becomes

*Ω*

_{i}is expressed in terms of the entries in

*ϖ*

_{k}in

*V*. Equations (2.28) and (2.30) when combined give

_{ij}or equivalently

*G*

_{ij}to be diagonal, equation (3.7) yields the relation

*I*

_{i}=

*g*′

_{ii}denotes, e.g. that the term

*g*

_{11}is omitted from the trace of

*g*in the definition of

*ϖ*

_{1}, etc.

Analogous to equation (2.30), write
*r*_{i}|≤1 each *i*, and that the *r*_{i} are not independent, but related by
*r*_{3}=−(*r*_{1}+*r*_{2})/(1+*r*_{1}*r*_{2}), and the inverse relations are, for example,
*ϖ*′′ and *ϖ*, *ι*′′ and *ι* may, respectively, be identified. Further, although the above approach may seem to enable the generalization of Dolzhansky’s approach to a non-orthogonal coordinate system, this does not, in fact, constitute a different physical situation, because the affine transformation of an ellipsoid is another ellipsoid.

The Kirchhoff equations have a conserved Hamiltonian
*H*_{0}=*I*_{1}*C*_{1}+*I*_{2}*C*_{2}+*I*_{3}*C*_{3} and *C*_{0}=*C*_{1}+*C*_{2}+*C*_{3}.

### (b) Catastrophic behaviour

It is of interest for comparison with ideal MHD to understand the transient behaviour of Dolzhansky’s equations. Supposing that the magnetic field has negligible dynamical effect, then it evolves kinematically in a flow, described by vorticity variables that obey Euler’s equations for the motion of a massless spinning top in classical mechanics. For such a body, it is a classical result that if the *I*_{i}, regarded now as moments of inertia of the top, satisfy *I*_{3}>*I*_{1}>*I*_{2} then motion about the 1-axis is unstable. Thus a significant transient is expected when a slow variation of *a*_{i} is arranged such that *I*_{3} approaches and then drops below *I*_{1}, so that rotation about the 3-axis is destabilized.

This transient corresponds to the near disappearance (see appendix Ab) of the effective potential when *r*_{2}=0, so that the system moves ballistically on the timescale of *ω*_{3}≈1 away from a now unstable equilibrium, i.e. the rate at which the instability threshold is crossed is of little importance. This is illustrated in figure 2, where the initial conditions ensure in fact that *ι*_{i}=0 for all time. This behaviour may also be deduced from the properties of Jacobi elliptic functions under parameter variation.

Energy remains bounded in the case of top dynamics, because the system ends in a different well with different rotation and translation direction. The implications for field kinematics is the possibility of a sudden transient which redirects not only the vortex, but also the electric current direction. This is reminiscent of the ideal MHD kink where components of field and current in new directions appear. It has, however, to be established that this behaviour extends to the case of a dynamically active magnetic field.

Figures 3 and 4 support this contention. All the field components now have initial conditions 0.01, except *ω*_{3}=0.1 and *ι*_{3}=1, so that this is a study of the full Dolzhansky equations where the magnetic field is dominant. Figure 3 exhibits a similar transient to the ‘spinning top’ or field-free case when the stability boundary *r*_{1}=0 (see appendix Aa) is again crossed at the rate 0.001. The exact nature of this transient requires further examination, but it will be seen from figure 4 that the timescale for, say, *ι*_{3} to reverse changes only by a factor of three (from approx. 40 to 110) when the crossing rate drops by a factor of 100.

### (c) Lie algebra

The scaling equation (3.1) applied to the basis functions equation (3.2) leads to Dolzhansky’s basis functions **e**_{i} (**W**_{i} in the notation of reference [3])
**e**_{i}=0, hence, the Lie derivative of one basis vector with another may be efficiently evaluated as
**e**_{i} and its curl.

The latter quantity is
**c**_{i},**c**_{j}]=**0**. Given these new Lie brackets, the ODE, equation (3.11), may be written down immediately from the PDE, equation 1.2.

## 4. Resistivity

Although the electrical resistance of very hot plasma is very low, plasma flow may concentrate magnetic variation into thin layers or current sheets, where field diffusion becomes important locally. The resulting local reconnection of fieldlines matters because it can cause a global change in topology, allowing energy to escape both from laboratory devices and the solar vicinity. Resistivity, either in the classical isotropic case, or possibly as the result of a renormalization approach to turbulence, leads to an additional term in the induction equation, which, assuming *ρ*=const., becomes
**B** has a linear field representation in Cartesian coordinates, then ∇×*η*∇×**B**=∇*η*×∇×**B**. Considering first the spherical, or the unscaled case as discussed at the start of §3a, the resistive term contains contributions of the form **e**_{i} if ∇*η*∝**x**. It is then plausible, as may be verified by direct substitution, that, in scaled coordinates, if the resistivity is written
*η*_{0} and *η*_{1} do not directly affect the model evolution of linear field **B**).

It follows that, for quadratic spatial dependence of resistivity of the form equation (4.2), the equations for electric current variation acquire terms
*η* must be maximal at the origin (*η*_{2}<0) otherwise the new resistivity terms imply exponential growth of the square current *et al.* [13], who find that only for locally maximum resistivity is the Petschek mechanism for reconnection structurally stable.

## 5. Conclusion

This work has reinforced the contentions of ref [3] that the D–K equations exhibit mathematical properties important for understanding nonlinear MHD in the limit of small or vanishing resistivity. Section 2 illustrates the general mathematical framework into which the equations fit; then §3a shows how the D–K equations emerge when a solution for ideal MHD is sought as a Taylor series in Cartesian coordinates. Section 3b shows that catastrophism is natural in the system, in a consistent sense, namely that model variables remain bounded, despite the dynamic timescale. These results enable optimism regarding the answer to the long-standing question about the structural stability of all ‘linear field’ solutions to the growth of higher-order field terms [6], because fast timescales have now been seen for linear fields.

The derivation of the model also illustrates important features of the Lie derivative, specifically its anti- or skew-symmetric nature as demonstrated by its replacement by the matrix-commutator in §2. In §2 also, important conservation relations, extending to arbitrary order of power of the unknowns, are illustrated. Continuing the mathematical note, §3c shows the importance of the concept of Lie algebra when seeking time-dependent solutions of nonlinear MHD. All these properties imply that the D–K equations should also be an aid to understanding the reduction of PDEs to non-canonical Hamiltonian systems and subsequent analysis [5,14]; cf. also reference [6].

As already incidentally demonstrated by the catastrophic simulations, the sixth-order D–K equations, with four invariants, are capable of exhibiting oscillation with two different timescales and amplitudes over a wide range of frequencies. The reduction of nonlinear MHD to a low-order Hamiltonian system highlights the likely prevalence of oscillation in ideal MHD, because such Hamiltonian systems do not generically possess attracting steady solutions and their stability can only be established in a Lyapunov sense. It is plausible that these considerations extend to the case of small resistivity when current sheets do not form.

Other important physical behaviour, such as the existence and behaviour of nonlinear Alfvenic solutions, corresponding to *ι*_{k}=*c*_{(k)}*ϖ*_{k}, for some constants *c*_{k}, may be deduced from known results for the Kirchhoff equations [12].

Section 4 shows how resistivity may be included in the model so that the physics of reconnection, believed important in many laboratory and astrophysical contexts, may be studied. In the context of laboratory plasmas, specifically the central region of tokamaks, the structural stability question raised by Forbes *et al.* requires further analysis.

Dolzhansky [3] describes how other important physics such as rotation and gravity (with buoyancy force owing to a density evolving due to a thermal evolution equation) may be included. Another application of linear fields, which will require development of the Lie–Taylor expansion for the momentum conservation equation along the lines of the mass conservation equation in §2b, is to partially ionized plasma where mass and momentum sources allow *out*flow boundary conditions, e.g. in one-dimensional models of the tokamak edge.

Further mathematical and physical insights into all these more complicated inviscid or almost-inviscid situations may be anticipated. As touched upon in reference [9], the vorticity evolution equation presents the problem that the algebra must involve not only basis functions but their curls. Progress may be made using Beltrami or ‘screw’ fields because these offer the potential to explore non-orthogonal geometries without the explicit appearance of the metric tensor, and indeed work on these lines is represented by the helical mode approach involving a triad of three-dimensional Fourier modes introduced by Waleffe and extended to MHD in reference [15]. Details of this approach, which generates equations resembling D–K but in complex mode amplitudes, are however perhaps better explored elsewhere.

## Data accessibility

This work contains sufficient information to repeat the calculations described, using where necessary the publicly available package mentioned in the Acknowledgements.

## Competing interests

I have no competing interests.

## Funding

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement number 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## Acknowledgements

The ODE integrations were performed using Hindmarsh’s LSODE package.

## Appendix A

**(a) Stability analysis**

Note that *ω* is used as a synonym for *ϖ* in the following.

**(i) Aligned vorticity and current**

Without loss of generality, assume that the vorticity and current are aligned in the direction of the 3-axis, with *ω*_{3}=*W* and *ι*_{3}=*J* and all other components of *t*=0. With these assumptions equations (3.11) and (3.12) become, respectively,

Seeking solutions varying in time *m*_{11}=*r*_{1}*J*^{2}−*W*^{2}, *m*_{12}=*JW*(1−*r*_{1}) and *m*_{22}=*r*_{2}(*r*_{1}*W*^{2}−*J*^{2}). The resulting stability polynomial is
*x*=*s*^{2}) are negative, implying *β*<0 and *r*_{1}<0, where *β* is the coefficient of *s*^{2}. The former means
*r*_{1}=0 is indeed a stability boundary.

**(ii) Orthogonal vorticity and current**

Without loss of generality, assume that the vorticity and current are aligned in the directions of the 3-axis and the 1-axis, respectively, with *ω*_{3}=*X* and *ι*_{1}=*K* and all other components of *t*=0. With these assumptions equations (3.11) and (3.15) become, respectively,
*ι*_{2} grows on a *X*=*o*(1) or *K*=*o*(1).

In the former case *X*=*o*(1), differentiating equation (A 10) gives
*r*_{2}<0, but even so *ω*_{3} varies in time proportional to *ι*_{2}. Similarly, in the latter case *K*=*o*(1)
*ι*_{1}.

Thus it seems that there is no stable steady solution with orthogonal vorticity and current. This conclusion is supported by analysis with *ω*_{3} and *ι*_{1} are constant. For then,
*ω*_{3} oscillates with frequency *r*_{3}>0). It also follows that *ι*_{1} and *ι*_{2} oscillate with frequencies *X* and *ω*_{3} oscillation is slow.

For the *ω*_{1} and *ι*_{3} dynamic, there is a determinantal equation *m*′_{11}=*r*_{1}*r*_{2}*X*^{2}, *m*′_{12}=*r*_{2}*KX* and *m*′_{22}=*r*_{2}*K*^{2}. This leads to a quadratic in *x*=*s*^{2} with roots *x*=0 and
*r*_{2}<0 and *r*_{1}<−*K*^{2}/*X*^{2} or *r*_{2}>0 and *r*_{1}>−*K*^{2}/*X*^{2}.

**(b) ‘Disappearance’ of the potential**

The surprising, at first hearing, remark that the potential almost disappears when *I*_{3}=*I*_{1} is explained if each component of vorticity is treated as evolving in its own separate potential. The relevant equations may be derived by differentiating each of the equations for *ϖ*_{i}=*ω*_{i} in equation (3.11) separately with respect to time, then eliminating first derivatives in terms of the *ω*_{i}. This gives equations specifiable, given
*C*_{1},*C*_{2},*C*_{3} in equation (A 17) and permutations then gives

Each equation (A 18) represents the motion of a particle in a potential *V* _{i}, e.g. given by
*ω*_{3}, and obvious permutations for the other *ω*_{i}. When *I*_{3}=*I*_{1}, *r*_{2}=0 from equation (3.12), it follows that the quadratic potential contributions to *V* _{i} vanish for each *i*, and only the quartic term in *V* _{2} survives. The corresponding variable is *ω*_{2}, which is small at least initially, hence the particle representing vector ** ω** is able to move ballistically, as though there were no potential forces present.

Table 1 gives values to two decimal places for the coefficients of 2*V* _{i} at the start of the simulation shown in figure 2, all of which change sign at *r*_{2}=0 except the one marked with a dagger. The corresponding dependence of the potential *V* _{i} on *ω*_{i} is denoted by ‘*U*’,‘*W*’ or ‘*M*’, with the letter corresponding to the approximate shape of the potential. Thus *ω*_{3} sits stably in the well at the bottom of the right-hand of the ‘*W*’ potential, which flips to become an ‘*M*’ shape in *r*_{2}>0 whereupon *ω*_{3} oscillates in the well in the centre of the ‘*M*’.

- Received June 28, 2016.
- Accepted November 25, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.