1 Introduction
Throughout the universe electrically conducting fluid flows interact with magnetic fields. By the stretching and folding of magnetic field lines, initially weak fields can grow, a process known as dynamo action (Moffatt Reference Moffatt1978). As the magnetic field increases in strength it then resists deformation through Lorentz forces exerted on the flow. Eventually a state of fully developed MHD turbulence ensues, composed of a superposition of coherent structures and random eddies interacting with magnetic fields.
In order to understand the complex interactions of flows and magnetic fields over a wide range of spatial and temporal scales it is valuable to investigate simplified models. In geophysical and astrophysical flows, coherent structures are believed to play a key role in generating large-scale magnetic fields (Tobias & Cattaneo Reference Tobias and Cattaneo2008). Since these structures have correlation times much greater than the turnover time of the turbulence, one of the most illuminating models of flow-field interaction is the effect of a steady fluid flow with closed streamlines – a single closed eddy – on an initially uniform magnetic field. Motivation for this approach comes from the study of the coupling of magnetic fields and convection (Weiss & Proctor Reference Weiss and Proctor2014); furthermore, in rapidly rotating convection the dynamics may be dominated by long-lived coherent structures taking the form of vortices (Julien et al. Reference Julien, Rubio, Grooms and Knobloch2012). The problem has been addressed in a kinematic regime in which the magnetic field is presumed to be so weak as not to affect the flow, which may then be specified and ceases to have any dynamical attributes. This classic problem was first studied for a smooth flow in the pioneering numerical study of Weiss (Reference Weiss1966), and by Parker (Reference Parker1966) for the case of a piecewise smooth flow. These works led to identification of the fundamental process of flux expulsion, whereby the field is destroyed within regions of closed streamlines; in cellular flows the resulting magnetic fluxes are then concentrated on bounding separatrices.
The mathematical theory of flux expulsion is elucidated for linear and axisymmetric, smooth flow fields by Moffatt & Kamkar (Reference Moffatt, Kamkar and Soward1983) and for more general streamline geometry by Rhines & Young (Reference Rhines and Young1983). The effect of a closed eddy on a weak imposed field at high magnetic Reynolds number
$R_{m}$
(that is,
$UL/{\it\eta}$
where
$U$
is the characteristic flow speed,
$L$
is a characteristic length and
${\it\eta}$
is the magnetic diffusivity) is to expel the field towards the cell boundaries on a time scale of order
$TR_{m}^{1/3}$
, where
$T=L/U$
is the turnover time scale. The key mechanism is the effect of shear or differential rotation in reducing length scales and so accelerating diffusion, be it of magnetic vector potential, passive scalar or vorticity: the useful and general term ‘shear–diffuse mechanism’ was coined by Bernoff & Lingevitch (Reference Bernoff and Lingevitch1994). These studies are elaborated in Bajer (Reference Bajer1998) and Bajer, Bassom & Gilbert (Reference Bajer, Bassom and Gilbert2001), referred to as BBG in what follows. In BBG a further time scale is identified at the eddy centre: here any differential rotation must vanish for a smooth flow, and the shear–diffuse process is weaker. The flux expulsion time scale here increases to order
$TR_{m}^{1/2}$
and a weak remnant, which we will call a ‘magnetic core’, is created and decays exponentially on this time scale.
Our goal in the present paper is to extend these kinematic studies into the dynamical regime, in which the field affects the flow via the Lorentz force. While it is clear that for sufficiently weak magnetic fields kinematic results are recovered, we address the problem of determining the threshold for the field to have a dynamical effect for the classic problem of flux expulsion in an axisymmetric flow. From another viewpoint, the question becomes: for what field strengths (in a two-dimensional flow) will a magnetic field have an impact on the material conservation of vorticity? Many interesting dynamical effects in quasi-two-dimensional hydrodynamics in rotating systems, such as zonal flow generation, have been ascribed to the material conservation properties of potential vorticity (Dritschel & McIntyre Reference Dritschel and McIntyre2008). It is therefore critical for MHD studies to determine the effect of magnetic fields in modifying these conservation properties. It is known that such dynamical effects of the magnetic field can be subtle and depend sensitively on molecular values of the transport coefficients. Perhaps the key study that highlights this is Cattaneo & Vainshtein (Reference Cattaneo and Vainshtein1991): the authors consider transport in turbulent two-dimensional fluid flows with a mean magnetic field
$b_{0}$
across the system. For field strengths
$b_{0}$
(in Alfvén or velocity units) of order
$UR_{m}^{-1/2}$
, the mean field becomes dynamically important in suppressing the stretching of Lagrangian parcels and so transport; see also Vainshtein & Cattaneo (Reference Vainshtein and Cattaneo1992), Cattaneo (Reference Cattaneo1994) and more recently Kim (Reference Kim2006), Keating & Diamond (Reference Keating and Diamond2008) and Keating, Silvers & Diamond (Reference Keating, Silvers and Diamond2008). The reason is that the stretching of the field leading to small scales (transverse to field lines) is accompanied by increasing magnetic energy, which is only limited by the effect of molecular (not turbulent or effective) diffusion. To put this more bluntly, in the magnetic context a ‘cascade’ (turbulent or otherwise) of field to small scales is far less passive than appears to be the case for energy or enstrophy in three- or two-dimensional turbulence. Similar effects arise in magnetoconvection (Galloway, Proctor & Weiss Reference Galloway, Proctor and Weiss1978; Weiss & Proctor Reference Weiss and Proctor2014), in which flux expulsion drives magnetic flux to the boundary of convective cells. In dynamical regimes, the resulting peak fields are limited not by equipartition values
$b_{0}=O(U)$
but can be substantially larger, with dependence on the magnetic diffusivity and fluid viscosity.
To address the problem of dynamical effects and thresholds for flux expulsion, in this paper we will work in the most straightforward setting of a flow with initially circular streamlines permeated by a uniform magnetic field of strength
$b_{0}$
. The problem is set up in § 2, parameterised by
$b_{0}$
and values of the magnetic Reynolds number
$R_{m}$
and (fluid) Reynolds number
$R_{e}$
. We will work in the quasi-linear approximation in which we keep the fluid flow axisymmetric and truncate the Lorentz force feedback by retaining only the mean azimuthal component. This is reasonable as we are examining the onset of the importance of the Lorentz force in the dynamics. The approach can be justified in some contexts such as dynamos in high-Reynolds-number rotating flow (as we have here) (Bassom & Gilbert Reference Bassom and Gilbert1997) and transport and jet formation in geophysical systems, for example see Tobias, Dagon & Marston (Reference Tobias, Dagon and Marston2011) and Srinivasan & Young (Reference Srinivasan and Young2012). Note that in a rotating fluid flow, the axisymmetric component of the Lorentz force can only be balanced against viscous terms, whereas a non-axisymmetric component can be balanced against weak inflows and outflows (cf. the Taylor (Reference Taylor1963) constraint), and for this reason the response to the axisymmetric component dominates at large Reynolds number, as in the quasi-linear model.
In § 3 of the paper we present numerical simulations of the evolution of field and flow (within the quasi-linear model) for a range of initial magnetic field strengths. Here we observe the key competition or race between the processes of flux expulsion and Lorentz force feedback: will the Lorentz force act early enough to halt the stretching in the flow and so defuse the dramatic effect of flux expulsion in destroying the field? Or, will flux expulsion act first and cut elastic field lines so as to remove the Lorentz force feedback and leave a magnetic core behind? We work in a regime in which
$R_{e}\gg R_{m}\gg 1$
to highlight the transition between these effects without consideration of viscosity, and this gives us our first threshold
$b_{core}({\it\eta})$
for
$b_{0}$
, as discussed below. However we also obtain a second, lower threshold
$b_{dynam}({\it\eta})$
, below which the magnetic field has no discernible effect and the kinematic picture holds. In §§ 4 and 5 we develop the theory for the two thresholds identified in § 3, based both on the classic asymptotics in Moffatt & Kamkar (Reference Moffatt, Kamkar and Soward1983) and Rhines & Young (Reference Rhines and Young1983) and on the more elaborate picture in BBG, which we need to identify the lower threshold. Finally § 6 offers some concluding comments and avenues for further research.
2 Governing equations
Our starting point is the equations for MHD, written in the form
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}=\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{b}-\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}\boldsymbol{b}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{b}=\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+{\it\eta}{\rm\nabla}^{2}\boldsymbol{b}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{b}=0. & \displaystyle\end{eqnarray}$$
The two-dimensional flow and field (measured in velocity units) are confined to the plane
$z=0$
in cylindrical polar coordinates
$(r,{\it\theta},z)$
, and setting
$\boldsymbol{u}=\boldsymbol{{\rm\nabla}}\times ({\it\psi}\hat{\boldsymbol{z}})$
,
$\boldsymbol{b}=\boldsymbol{{\rm\nabla}}\times (a\hat{\boldsymbol{z}})$
we obtain the equations in terms of the stream function
${\it\psi}$
, vorticity
${\it\omega}$
, flux function
$a$
and current
$j$
, all functions of
$(r,{\it\theta},t)$
, as
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}{\it\omega}=J({\it\psi},{\it\omega})-J(a,j)+{\it\nu}{\rm\nabla}^{2}{\it\omega}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}a=J({\it\psi},a)+{\it\eta}{\rm\nabla}^{2}a, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\omega}=-{\rm\nabla}^{2}{\it\psi},\quad j=-{\rm\nabla}^{2}a. & & \displaystyle\end{eqnarray}$$
Here the Jacobian is given by
$J({\it\psi},{\it\omega})=r^{-1}[(\partial _{r}{\it\psi})(\partial _{{\it\theta}}{\it\omega})-(\partial _{{\it\theta}}{\it\psi})(\partial _{r}{\it\omega})]$
.
Following the discussion in Moffatt & Kamkar (Reference Moffatt, Kamkar and Soward1983), we commence with an initial uniform magnetic field in the
$x$
-direction and an axisymmetric flow field. Thus
${\it\psi}$
and
${\it\omega}$
are taken to be independent of
${\it\theta}$
and the fields
$a$
and
$j$
are represented using an
$\text{e}^{\text{i}m{\it\theta}}$
dependence with
$m=1$
. Although our focus is always on
$m=1$
, it is helpful in the analytical development to leave a general, integer value of
$m>0$
. We therefore set for the flow
$$\begin{eqnarray}\displaystyle {\it\omega}={\it\omega}(r,t)+\cdots ,\quad {\it\psi}={\it\psi}(r,t)+\cdots , & & \displaystyle\end{eqnarray}$$
and for the field
$$\begin{eqnarray}\displaystyle \hspace{-10.00002pt}a=\tilde{a}(r,t)\text{e}^{\text{i}m{\it\theta}}+\text{c.c.}+\cdots ,\quad j=\tilde{j}(r,t)\text{e}^{\text{i}m{\it\theta}}+\text{c.c.}+\cdots ,\quad \boldsymbol{b}=\tilde{\boldsymbol{b}}(r,t)\text{e}^{\text{i}m{\it\theta}}+\text{c.c.}+\cdots .\hspace{-30.00005pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The tildes denote the harmonic
$m>0$
in
${\it\theta}$
, but for readability we drop these in what follows.
Now the Lorentz force feedback from the field to the flow will incorporate a mean part, independent of
${\it\theta}$
, and harmonics
$\text{e}^{2\text{i}m{\it\theta}}$
, which will then proliferate, giving the trailing terms not written down explicitly in (2.7a,b
) and (2.8a−c
). We employ a truncation by neglecting these higher-order harmonics, and retain only the terms shown. This leaves the quasi-linear system for the field harmonic
$m$
and the mean flow, written compactly as
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}a+\text{i}m{\it\alpha}\,a={\it\eta}{\it\Delta}_{m}a, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}{\it\omega}=r^{-1}\partial _{r}G+{\it\nu}{\it\Delta}_{0}{\it\omega}. & \displaystyle\end{eqnarray}$$
The angular velocity is
${\it\alpha}(r,t)=-r^{-1}\partial _{r}{\it\psi}$
, and the current and vorticity satisfy
$$\begin{eqnarray}\displaystyle j=-{\it\Delta}_{m}a,\quad {\it\omega}=-{\it\Delta}_{0}{\it\psi}, & & \displaystyle\end{eqnarray}$$
where
${\it\Delta}_{m}=\partial _{r}^{2}+r^{-1}\partial _{r}-m^{2}r^{-2}$
. The Lorentz force term in (2.10) is
$G(r,t)$
given by
$$\begin{eqnarray}\displaystyle G=\text{i}m(aj^{\ast }-a^{\ast }j). & & \displaystyle\end{eqnarray}$$
We will use this quasi-linear approximation to gain an understanding of the essential processes in the competition between flux expulsion and the Lorentz force, both numerically and analytically.
Figure 1. The angular velocity profile
${\it\alpha}$
(solid) in (2.16a,b
) and
${\it\alpha}^{\prime }$
(dashed), with the location of maximal differential rotation
$|{\it\alpha}^{\prime }|$
identified as
$r^{\dagger }\simeq 1.8$
.
For initial conditions, our focus is on the case of a uniform field in the
$x$
-direction of strength
$b_{0}$
, which has current
$j=0$
. However we generalise to an arbitrary value of
$m$
, and consider an initial multipole, current-free field given by the vector potential
$$\begin{eqnarray}\displaystyle a(r,0)=-{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m} & & \displaystyle\end{eqnarray}$$
(with
$b_{0}$
real). We also need an initial axisymmetric fluid flow
$\boldsymbol{u}$
and have chosen a Gaussian vortex with
$$\begin{eqnarray}\displaystyle \hspace{-10.00002pt}{\it\omega}(r,0)=T^{-1}(4{\rm\pi})^{-1}\text{e}^{-r^{2}/4L^{2}},\quad {\it\alpha}(r,0)=T^{-1}(2{\rm\pi}r^{2}/L^{2})^{-1}(1-\text{e}^{-r^{2}/4L^{2}}), & & \displaystyle\end{eqnarray}$$
where
$L$
is a length scale,
$T$
a time scale and below
$U=L/T$
. We will non-dimensionalise the system using these scales and defining, for example,
$$\begin{eqnarray}\displaystyle r=L\hat{r},\quad t=T\hat{t},\quad \boldsymbol{u}=U\hat{\boldsymbol{u}},\quad \boldsymbol{b}=U\hat{\boldsymbol{b}},\quad {\it\nu}=UL\hat{{\it\nu}},\quad {\it\eta}=UL\hat{{\it\eta}}, & & \displaystyle\end{eqnarray}$$
with hats denoting non-dimensional quantities. With this we can also identify
$R_{m}\equiv \hat{{\it\eta}}^{-1}$
as a magnetic Reynolds number,
$R_{e}\equiv \hat{{\it\nu}}^{-1}$
as a Reynolds number and
$M\equiv \hat{b}_{0}^{-1}$
as a magnetic Mach number. We will in what follows drop the hats and work with dimensionless quantities, except when we refer to our results in the final discussion section.
Our goal then is to solve the PDEs specified in (2.9)–(2.12), for the initial conditions (2.13) and now
$$\begin{eqnarray}\displaystyle {\it\omega}(r,0)=(4{\rm\pi})^{-1}\text{e}^{-r^{2}/4},\quad {\it\alpha}(r,0)=(2{\rm\pi}r^{2})^{-1}(1-\text{e}^{-r^{2}/4}). & & \displaystyle\end{eqnarray}$$
The angular velocity
${\it\alpha}$
and differential rotation
${\it\alpha}^{\prime }\equiv \partial _{r}{\it\alpha}$
are depicted in figure 1, where for convenience we often use a prime to denote a radial derivative. The radius
$r^{\dagger }$
marks the location of maximal differential rotation
$|{\it\alpha}^{\prime }|$
, where we will see that flux expulsion commences in a kinematic regime.
The parameter set comprises the three non-dimensional quantities
$\{{\it\eta},{\it\nu},b_{0}\}$
, and we use this form as it is more convenient to place ‘
${\it\eta}$
’ rather than the bulky term ‘
$R_{m}^{-1}$
’ in our calculations. We are interested in the regimes that are realised depending on the strength of the initial field for different diffusive parameters. In this study our primary interest is in the interaction of flux expulsion (depending on
${\it\eta}$
) and the Lorentz force (linked to
$b_{0}$
), rather than viscous effects. We will thus take
${\it\nu}\ll {\it\eta}$
and so work at a low value of the magnetic Prandtl number
$P_{m}\equiv R_{m}/R_{e}={\it\nu}/{\it\eta}$
. (Note that viscous damping would emerge on a time scale of order
$R_{e}$
whereas the effects we consider occur on time scales of order
$R_{m}^{1/3}$
and
$R_{m}^{1/2}$
, and thus our results are likely to be correct over a wider range of values of
$P_{m}$
including
$P_{m}=O(1)$
.) In all the simulations shown, we have simply taken
$P_{m}=0.01$
, and tests confirm that our results are insensitive to this precise value. Our parameter set is thus reduced to
$\{{\it\eta},b_{0}\}$
and we are interested in thresholds for different types of behaviour, giving
$b_{0}$
as a function of
${\it\eta}$
with a power-law scaling: how strong must the initial field be for the Lorentz force to modify the classic picture of flux expulsion?
Theory for these thresholds is developed in § 4. Before giving numerical results in the next section it is worth outlining the origin of the classic approximation for kinematic flux expulsion (for more detail see Moffatt & Kamkar Reference Moffatt, Kamkar and Soward1983 and Rhines & Young Reference Rhines and Young1983). Consider when there is no dissipation; then the vector potential and current are given by
$$\begin{eqnarray}\displaystyle a=-{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m}\text{e}^{-\text{i}m{\it\alpha}t},\quad j\simeq -{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m}m^{2}{\it\alpha}^{\prime 2}t^{2}\text{e}^{-\text{i}m{\it\alpha}t}, & & \displaystyle\end{eqnarray}$$
where we have retained in
$j$
just the terms that grow fastest, i.e. quadratically with time
$t$
. For small
${\it\eta}$
, the right-hand side
${\it\eta}{\it\Delta}_{m}a=-{\it\eta}j$
in (2.9) grows quadratically from small values and the accumulated effect of the dissipation gives a term cubic in
$t$
. Incorporating this damping in the evolution of
$a$
then yields the approximate solution
$$\begin{eqnarray}\displaystyle a\simeq -{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m}\text{e}^{-\text{i}m{\it\alpha}t-(1/3){\it\eta}m^{2}{\it\alpha}^{\prime 2}t^{3}},\quad j\simeq -{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m}m^{2}{\it\alpha}^{\prime 2}t^{2}\,\text{e}^{-\text{i}m{\it\alpha}t-(1/3)m^{2}{\it\eta}{\it\alpha}^{\prime 2}t^{3}}. & & \displaystyle\end{eqnarray}$$
This gives dramatic suppression of the vector potential and so of the magnetic field, which commences at the radius
$r^{\dagger }$
where
$|{\it\alpha}^{\prime }|$
is maximised, indicated on figure 1 for the Gaussian profile (2.16a,b
), on a time scale
$t^{\dagger }=O({\it\eta}^{-1/3})$
. Naturally, in the kinematic problem the flow profile is taken to be fixed for all
$t$
and not viscously damped; only the magnetic field evolves via (2.9) in the given flow field with
${\it\alpha}(r,t)={\it\alpha}(r,0)$
. Specifically
$r^{\dagger }$
and
$t^{\dagger }$
(up to an order-one constant) are given by
$$\begin{eqnarray}\displaystyle r^{\dagger }=\arg \max _{r}|{\it\alpha}^{\prime }(r)|^{2},\quad t^{\dagger }=\left[{\textstyle \frac{1}{3}}m^{2}{\it\eta}{\it\alpha}^{\prime }(r^{\dagger })^{2}\right]^{-1/3}. & & \displaystyle\end{eqnarray}$$
The dagger helpfully denotes the cutting of magnetic field lines caused by flux expulsion.
For our purposes there are two problems with this approximation. The first is that the Lorentz force term (2.12) vanishes identically to this order. The next-order terms need to be taken into account, and we undertake this in a systematic development in § 4.2. The second issue is discussed in BBG and is that this approximation breaks down near the origin where
${\it\alpha}^{\prime }\rightarrow 0$
. In other words the approximation is predicated on a shear
${\it\alpha}^{\prime }=O(1)$
that makes the quadratic multiplier
$t^{2}$
in
$j$
in (2.18a,b
) dominant; towards the centre of a vortex another, inner approximation has to be used. Although this may seem a technical issue, in fact it is key to understanding the dynamical problem, and we develop this in § 4.1 following BBG. In preparation for this we note here that the angular velocity in (2.16a,b
) behaves as
$r\rightarrow 0$
according to
$$\begin{eqnarray}\displaystyle {\it\alpha}={\it\alpha}_{0}+{\it\alpha}_{1}r^{2}+{\it\alpha}_{2}r^{4}+\cdots , & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\alpha}_{0}=1/8{\rm\pi},\quad {\it\alpha}_{1}=-1/64{\rm\pi}. & & \displaystyle\end{eqnarray}$$
From the point of view of the kinematic problem, with the axisymmetric flow specified independent of the field, the form (2.20) gives the behaviour near the origin. Smoothness considerations eliminate any odd powers of
$r$
and generically the constant
${\it\alpha}_{1}$
is non-zero. This constant plays an important role in kinematic theory as it controls how the differential rotation
${\it\alpha}^{\prime }$
responsible for flux expulsion switches off near the origin. Dynamically the magnetic field can change the form of the flow near the origin, as we shall see. We also note that the numerical values of
${\it\alpha}_{0}$
and particularly
${\it\alpha}_{1}\simeq -1/200$
are rather low, meaning that the turnover time of the flow is large in our non-dimensionalisation and making the corresponding time scales appear long on our plots. For best comparison with simulations we will retain factors of
${\it\alpha}_{1}$
in our theoretical development.
3 Numerical simulations
3.1 Illustrative runs
Our goal in this section is to present numerical simulations of the model (2.9)–(2.12) with the initial conditions (2.13) and (2.16a,b
). This will motivate the theory in the next section, but we will also refer ahead to results in that section. Our parameters are only
$\{{\it\eta},b_{0}\}$
and we first show the various phenomena that occur when we fix
${\it\eta}=10^{-7}$
(with
$m=1$
,
${\it\nu}=0.01{\it\eta}$
as always) and allow a range of values of
$b_{0}$
.
Figure 2. Kinematic magnetic field evolution. In (a,b)
$t=10^{4}$
and in (c)
$t=4\times 10^{4}$
. In (a)
$-10\leqslant x,y\leqslant 10$
, while in (b,c) only the central region
$-2\leqslant x,y\leqslant 2$
is depicted.
Our starting point is the kinematic problem when the field
$b_{0}$
is sufficiently weak that dynamical effects may be neglected. Although the Lorentz force feedback is easily switched off in our simulations, one of our goals is to quantify just how weak the initial field
$b_{0}$
needs to be for kinematic theory to apply, for
${\it\eta}\ll 1$
. In any case we integrate (2.9) in isolation. Figure 2 shows colour-scale plots of the vector potential reconstructed from (2.8a−c
), so that lines of constant colour give magnetic field lines, from blue (low values of
$a$
in (2.8a−c
)) to red (high values); zero or weak values of
$a$
are green. The panels show different times, and in each panel the colour scale is normalised separately on the maximum/minimum vector potential for that panel. Panel (a) gives a wide view, showing spiral wind-up of field lines by the flow and suppression of magnetic field by diffusion for moderate
$r$
. Panels (b,c) show a zoom into the central region, where the onset of flux expulsion around
$r^{\dagger }$
can be seen in (b), followed by further destruction of field leading to the final phase in (c), of a decaying, remnant field structure localised at the origin.
Figure 3. (a) Schematic of kinematic evolution of magnetic field for
${\it\eta}\ll 1$
, and (b) log–log plot of
$E_{M}$
(upper solid) and
$E_{A}$
(lower solid) as functions of time
$t$
. In (b) the scaling laws
$t^{-7}$
(4.10) and
$t^{-6}$
(4.11) are depicted (dashed), and the formulae (4.13) and (4.14) are shown (dotted).
Figure 3 shows a schematic of kinematic field evolution (adapted from BBG): at
$t^{\dagger }=O({\it\eta}^{-1/3})$
there is the onset of flux expulsion at a radius
$r^{\dagger }$
. This then spreads outwards and inwards, as indicated in (2.18a,b
). At late times
$t=O(\mathscr{T})=O({\it\eta}^{-1/2})$
(see (4.1a−c
)) this wave of destruction reaches the centre of the axisymmetric flow and all that is left behind is a rotating, exponentially decaying remnant, an eigenfunction of the scalar advection–diffusion equation, seen in figure 2(c).
For times greater than the flux expulsion onset time
$t^{\dagger }$
, we can identify the field inside
$r^{\dagger }$
as the magnetic core and it is convenient to define the corresponding magnetic energy and (half) the integrated square vector potential by
$$\begin{eqnarray}\displaystyle E_{M}(t)=\int _{0}^{r^{\dagger }}2{\rm\pi}r\,|\boldsymbol{b}|^{2}\,\text{d}r,\quad E_{A}(t)=\int _{0}^{r^{\dagger }}2{\rm\pi}r\,|a|^{2}\,\text{d}r. & & \displaystyle\end{eqnarray}$$
These are shown as functions of time in the log–log plot figure 3(b): commencing at
$t^{\dagger }\simeq 1.2\times 10^{4}$
they rapidly adopt a
$t^{-7}$
(4.10) and
$t^{-6}$
(4.11) decay in time. At time of order
$\mathscr{T}\simeq 3.2\times 10^{4}$
(4.1a−c
) the power-law decay is replaced by exponential decay (4.13), (4.14). Note that there is no absolute definition of
$t^{\dagger }$
: instead for practical purposes we check when
$|a(r,t)|<{\it\delta}|a(0,t)|$
with a small number
${\it\delta}$
: the earliest time at which this occurs determines
$t^{\dagger }$
, with
$r^{\dagger }$
as the corresponding radius. We have chosen
${\it\delta}=0.001$
for results shown here: other values make minor changes to the values of
$t^{\dagger }$
but do not affect
$r^{\dagger }\simeq 1.8$
and make no visible difference to the curves for
$E_{M}$
or
$E_{A}$
once a magnetic core is so defined.
We now leave the kinematic problem and report on runs for varying values of the initial magnetic field
$b_{0}$
. First figure 4 shows a run for a relatively strong field
$b_{0}=5\times 10^{-4}$
. Here we see the process of spiral wind-up commence in (a), but stretching of magnetic field lines saps energy from the flow field, reducing the shear and stopping the flux expulsion process. In panels (b,c) we see the field in the simulation start to unwind, reversing the direction of the original vortex. This disturbance then propagates outwards to hit the numerical boundary and, unphysically, bounces back and forth (not shown). In reality the original flow would turn into Alfvén waves propagating to infinity along the field lines. For this reason we do not spend further effort on the cases where the field is above the threshold for flux expulsion to occur, even though one of our primary goals is to establish this threshold, working from below.
Figure 4. Magnetic field evolution for
$b_{0}=5\times 10^{-4}$
. The region
$-10\leqslant x,y\leqslant 10$
is shown for (a)
$t=10^{4}$
, (b)
$t=2\times 10^{4}$
and (c)
$t=2.5\times 10^{4}$
.
Figure 5. Magnetic field evolution for
$b_{0}=2\times 10^{-4}$
. The central region
$-2\leqslant x,y\leqslant 2$
is shown for (a)
$t=10^{4}$
, (b)
$t=2\times 10^{4}$
and (c)
$t=5\times 10^{4}$
.
We now reduce the initial field
$b_{0}$
in subsequent plots, starting with figure 5 which shows magnetic field evolution for
$b_{0}=2\times 10^{-4}$
. In this case flux expulsion occurs to cut the field lines at approximately
$r^{\dagger }\simeq 1.8$
; see panel (a). However, once the field lines are cut to leave a magnetic core, the Lorentz force now acts to reduce field line curvature and magnetic energy within the core, in panel (b). What remains is a core consisting of two lobes of field relaxed to a state of low energy, and the flow field is modified so that there is solid body rotation in the region occupied by the field, shown in figure 6(a). We say that a dynamical core has been formed, dynamical as the Lorentz force has acted to control the flow. This core will then decay, but on a longer Ohmic time scale (presumably
$O({\it\eta}^{-1})$
, though we will not try to verify this).
Figure 7. Magnetic field evolution for
$b_{0}=5\times 10^{-5}$
. The central region
$-2\leqslant x,y\leqslant 2$
is shown for (a)
$t=10^{4}$
, (b)
$t=2\times 10^{4}$
and (c)
$t=5\times 10^{4}$
.
Finally figure 7 shows evolution of a yet weaker field
$b_{0}=5\times 10^{-5}$
. Here again there is a process of flux expulsion leaving a core which shrinks as the field is diffusively destroyed. The Lorentz force again acts to leave a flattened region in the flow field, constant angular velocity seen in figure 6(b), and a dynamical core with two lobes, albeit now smaller than in the previous case. For fields yet weaker than this, we soon find results that become indistinguishable from the kinematic regime in figure 2 and no dynamical core forms.
Figure 8. Plotted are (a) the magnetic energy
$E_{M}$
and (b)
$E_{A}$
as functions of time
$t$
for varying field strength
$b_{0}$
. The dashed curves are kinematic traces, while
$b_{0}=2\times 10^{-4}$
,
$10^{-4}$
,
$5\times 10^{-5}$
,
$2\times 10^{-5}$
,
$10^{-5}$
,
$5\times 10^{-6}$
and
$2\times 10^{-6}$
, going down the curves on the right of the plots.
Another viewpoint is given in figure 8, which plots (a)
$E_{M}$
and (b)
$E_{A}$
against time for a range of
$b_{0}$
. In each panel the lowest, dashed curve gives the kinematic time trace (from figure 2). For strong fields (highest, outermost curves), the onset
$t^{\dagger }$
of flux expulsion is delayed and the field decays very slowly, noting that
$E_{M}$
shows oscillations (torsional oscillations) while the dynamical magnetic core relaxes. On the other hand
$E_{A}$
behaves monotonically, as it must since the quantity
$a$
obeys a scalar advection–diffusion equation. As the initial field is reduced the traces follow the kinematic curves for a period and then depart: this marks the formation of a dynamically controlled core sitting in a region of solid body motion. Finally very weak field follows the kinematic traces.
3.2 Thresholds and scaling laws
Figure 9. (a) The onset of flux expulsion
$t^{\dagger }$
as a function of initial field strength
$b_{0}$
using different thresholds
${\it\delta}=10^{-4}$
,
$10^{-3}$
and
$10^{-2}$
, reading down the curves. (b)
$E_{A}$
plotted against
$b_{0}$
at times
$\mathscr{T}/2$
,
$\mathscr{T}$
and
$2\mathscr{T}$
, reading down the curves.
We can summarise the results of the previous section as the presence of two thresholds. For a given flow field, the first is the threshold
$b_{0}=b_{core}({\it\eta})$
above which the Lorentz force is so strong that flux expulsion does not occur and so
$t^{\dagger }$
,
$E_{M}$
and
$E_{A}$
cannot be defined. Figure 9(a) shows
$t^{\dagger }$
as measured numerically using several values of
${\it\delta}$
(see below (3.1a,b
)). There is a sharp transition from the onset of flux expulsion at a time independent of initial field strength, to one that diverges rapidly with
$b_{0}$
; at the same time the radius
$r^{\dagger }$
increases from 1.8 to approximately 2.2 (not shown). This sharp transition makes the threshold
$b_{0}=b_{core}({\it\eta})$
easy to measure, at least approximately.
The second, lower threshold is
$b_{0}=b_{dynam}({\it\eta})$
above which the field is sufficiently strong (i.e. sufficiently dynamical) to modify the flow field and suppress diffusive decay. In this case the Lorentz force opposes differential rotation and results in solid body motion in a region near the origin. In terms of the flow field, this has the effect of turning off
${\it\alpha}_{1}$
in (2.20) and so removes the weak differential rotation at the origin that otherwise controls the exponential decay of the core in the kinematic regime; see figure 3(a) and § 4.1. Above the threshold we are left with a slowly decaying dynamical core. To measure this threshold we have considered the quantity
$E_{A}$
(unlike
$E_{M}$
this is free from torsional oscillations – see figure 8) at three times
$\mathscr{T}/2$
,
$\mathscr{T}$
and
$2\mathscr{T}$
and plot this against
$b_{0}$
in figure 9(b). Where the curves cluster together on the right-hand side, the energy decays negligibly on the
$\mathscr{T}$
time scale and we are in the dynamical core regime. Where the curves are spread out on the left-hand side, the field is so weak that the evolution is kinematic with a core decaying rapidly, that is, on the
$\mathscr{T}$
time scale. The somewhat broad transition between these marks the threshold
$b_{dynam}({\it\eta})$
, for this value of
${\it\eta}$
.
Figure 10. Thresholds
$b_{core}({\it\eta})$
(upper dataset) and
$b_{dynam}({\it\eta})$
(lower dataset) plotted against
${\it\eta}$
. In each case the data points come from a series of runs with varying
$b_{0}$
. The solid lines give the scalings
${\it\eta}^{1/3}$
(upper) and
${\it\eta}^{3/4}$
(lower). The dotted lines give the scalings from (5.4) and (5.7).
Neither of these thresholds is precisely defined, but in any case we are primarily interested simply in how they scale with
${\it\eta}$
. We select a representative transitional field strength
$b_{dynam}$
to be the initial value of
$b_{0}$
for which
$\log [E_{A}(2\mathscr{T})/E_{A}(\mathscr{T})]$
is half of its kinematic value. We estimated
$b_{dynam}$
using a series of runs with 10 values of
$b_{0}$
per decade on the logarithmic scale. For example in figure 9 we obtain
$b_{dynam}=7.9\times 10^{-6}$
for
${\it\eta}=10^{-7}$
. Figure 10 shows
$b_{dynam}$
estimated this way, in the lower dataset, showing good agreement with the scaling law
$b_{dynam}\sim {\it\eta}^{3/4}$
(solid) from (5.5).
To estimate the threshold
$b_{core}({\it\eta})$
for core formation we adopt two methods and these are shown in the upper dataset in figure 10. The first is obtained by selecting the minimum field strength
$b_{0}$
for which the ratio
$E_{A}(\mathscr{T})/E_{A}(2\mathscr{T})<1.08$
(asterisks), and this shows a good fit to the scaling law
$b_{core}\sim {\it\eta}^{1/3}$
in (5.3), with
$b_{core}\simeq 4\times 10^{-4}$
at
${\it\eta}=10^{-7}$
. For another method, we chose the largest
$b_{0}$
for which core formation, i.e. the existence of a value of
$t^{\dagger }$
, was detected in our code (plus signs), giving
$b_{core}\simeq 5\times 10^{-4}$
at
${\it\eta}=10^{-7}$
. This also confirms the scaling law, albeit with more scatter. Note that alternatively we can use the scalings to rescale the vertical and the horizontal axes in figure 9(a,b) for a set of values of
${\it\eta}$
and so collapse the curves, which works well but which we do not show here.
Finally we comment on the numerical methods used. Equations (2.9)–(2.12) were written as a system of six first-order PDEs for the real and imaginary parts of
$a$
,
$\partial _{r}a$
and
$v$
,
$\partial _{r}v$
. This system was then passed to the NAG solver d03pef, which employs a Keller box method and integrates in time, with
$r_{lo}\leqslant r\leqslant r_{hi}$
for
$r_{lo}=0.05$
,
$r_{hi}=20$
and up to
$2\times 10^{4}$
radial grid points in typical runs. At the inner and outer boundaries, the condition of behaviour as
$r^{m}$
was imposed on the vector potential, namely
$r\partial _{r}a-ma=0$
. For the flow component
$v$
, the code similarly imposed behaviour proportional to
$r$
at the inner boundary, and to
$r^{-1}$
at the outer. In the cases where flux expulsion occurs and a magnetic core is formed, the outer boundary does not play a role, as discussed above, and
$r_{hi}=20$
is sufficiently large (at least over the time scales shown here).
4 Calculation of the Lorentz force feedback
Our goal in this section is to derive the various scaling laws above, in particular those in figure 10. It is evident from figure 8 that the kinematic evolution is key, in that the dynamical curves typically follow the kinematic evolution, at least for a time. We consider the kinematic regime, corresponding to the limit
$b_{0}\rightarrow 0$
and a fixed flow field. We then assess the feedback on the flow through the Lorentz force.
4.1 Flux expulsion: inner solution
We begin by recalling the classic flux expulsion calculation yielding (2.18a,b
) as an approximate outer solution to (2.9). This gives
$a$
and
$j$
at radii where
${\it\alpha}^{\prime }(r)\neq 0$
(fixed as
${\it\eta}\rightarrow 0$
) and is valid for a general smooth profile
${\it\alpha}(r)$
. However, the approximation breaks down at the centre of the flow field where
${\it\alpha}^{\prime }(0)=0$
necessarily, and so another approximation, the inner solution, is needed there. Generally we assume that
${\it\alpha}(r)$
expands as (2.20) with
${\it\alpha}_{1}\neq 0$
, working in a kinematic regime, but we note that in a dynamical regime the Lorentz force acts to suppress differential rotation given by
${\it\alpha}^{\prime }(r)$
. Note that the theoretical development in § 4 below is written as if
${\it\alpha}_{1}$
is positive to follow BBG but the final results are dependent only on
$|{\it\alpha}_{1}|$
in any case. There are also several inessential notational differences between the theory in the earlier paper and in the self-contained development here. The kinematic theory is provided in BBG and is most easily set out by defining a length scale
$\mathscr{L}$
, time scale
$\mathscr{T}$
, their inverses and a velocity scale
$\mathscr{V}$
by
$$\begin{eqnarray}\displaystyle k\equiv \mathscr{L}^{-1}=(m{\it\alpha}_{1}/2{\it\eta})^{1/4},\quad p\equiv \mathscr{T}^{-1}=(2{\it\eta}m{\it\alpha}_{1})^{1/2},\quad \mathscr{V}=\mathscr{L}/\mathscr{T}=(8{\it\eta}^{3}m{\it\alpha}_{1})^{1/4} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
and new variables by
$$\begin{eqnarray}\displaystyle {\it\tau}=pt,\quad {\it\rho}=kr. & & \displaystyle\end{eqnarray}$$
The exact solution to (2.9) may then be obtained with only the quadratic term retained in
${\it\alpha}(r)$
in (2.20), that is,
${\it\alpha}\simeq {\it\alpha}_{1}r^{2}$
(where we also set
${\it\alpha}_{0}=0$
as solid body rotation is irrelevant). The solution is
$$\begin{eqnarray}\displaystyle & \displaystyle a=-{\textstyle \frac{1}{2}}\text{i}b_{0}k^{-m}{\it\rho}^{m}g\text{e}^{-\text{i}f{\it\rho}^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle j=-{\textstyle \frac{1}{2}}\text{i}b_{0}k^{2-m}{\it\rho}^{m}g\,[4\text{i}(m+1)f+4f^{2}{\it\rho}^{2}]\text{e}^{-\text{i}f{\it\rho}^{2}}, & \displaystyle\end{eqnarray}$$
where
$f({\it\tau})$
and
$g({\it\tau})$
satisfy
$$\begin{eqnarray}\displaystyle \partial _{{\it\tau}}g=-2\text{i}(m+1)fg,\quad \partial _{{\it\tau}}f-1=-2\text{i}f^{2}, & & \displaystyle\end{eqnarray}$$
and so
$$\begin{eqnarray}\displaystyle f({\it\tau})=(1+\text{i})^{-1}\tanh [(1+\text{i}){\it\tau}],\quad g({\it\tau})=\{\cosh [(1+\text{i}){\it\tau}]\}^{-m-1}. & & \displaystyle\end{eqnarray}$$
Now once a core has formed, for
$t>t^{\dagger }$
, undertaking the integrals in (3.1a,b
) gives straightforwardly
$$\begin{eqnarray}\displaystyle & \displaystyle E_{M}={\rm\pi}(m+1)!\,b_{0}^{2}\,k^{-2m}|g|^{2}|f|^{2}\,(-2f_{\text{i}})^{-m-2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{A}={\textstyle \frac{1}{4}}{\rm\pi}\,m!\,b_{0}^{2}\,k^{-2m-2}|g|^{2}\,(-2f_{\text{i}})^{-m-1}, & \displaystyle\end{eqnarray}$$
with
$f=f_{\text{r}}+\text{i}f_{\text{i}}$
for brevity.
This is all exact if
${\it\alpha}(r)$
is given by only the leading terms, in
${\it\alpha}_{0}$
and
${\it\alpha}_{1}$
in (2.20). The solution combines two different processes, with a cross-over time
$\mathscr{T}$
defined in (4.1a−c
). For
$t\ll \mathscr{T}$
, the approximation captures a wave of flux expulsion, destruction of
$a$
, in the locally quadratic angular velocity profile (2.20). In this regime, we have
$$\begin{eqnarray}\displaystyle f={\it\tau}-{\textstyle \frac{2}{3}}\text{i}{\it\tau}^{3}+\cdots ,\quad g\simeq 1,\quad -2f_{\text{i}}\simeq {\textstyle \frac{4}{3}}{\it\tau}^{3}, & & \displaystyle\end{eqnarray}$$
giving
$$\begin{eqnarray}\displaystyle & \displaystyle E_{M}\simeq {\rm\pi}\left({\textstyle \frac{4}{3}}\right)^{-m-2}(m+1)!\,b_{0}^{2}\,k^{-2m}\,{\it\tau}^{-3m-4}=O(t^{-3m-4}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{A}\simeq {\textstyle \frac{1}{4}}{\rm\pi}\left({\textstyle \frac{4}{3}}\right)^{-m-1}\,m!\,b_{0}^{2}\,k^{-2m-2}\,{\it\tau}^{-3m-3}=O(t^{-3m-3}) & \displaystyle\end{eqnarray}$$
and so algebraic behaviour of these quantities, with dependence as
$t^{-7}$
and
$t^{-6}$
in the important case
$m=1$
of an initial uniform field, as seen in figure 4.
On the other hand for
$t\gg \mathscr{T}$
the solution describes an exponentially decaying core taking a Gaussian form near the origin, with
$$\begin{eqnarray}\displaystyle f\simeq (1+\text{i})^{-1},\quad g\simeq 2^{m+1}\text{e}^{-(m+1)(1+\text{i}){\it\tau}},\quad -2f_{\text{i}}\simeq 1 & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle & \displaystyle E_{M}\simeq 2^{2m+1}\,{\rm\pi}\,(m+1)!\,b_{0}^{2}\,k^{-2m}\text{e}^{-2(m+1){\it\tau}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{A}\simeq 2^{2m}\,{\rm\pi}\,m!\,b_{0}^{2}\,k^{-2m-2}\text{e}^{-2(m+1){\it\tau}}. & \displaystyle\end{eqnarray}$$
Again this is confirmed by the results shown in figure 4 for
$m=1$
.
4.2 Lorentz force from the outer solution
We need to evaluate the Lorentz torque
$G$
in (2.12) from the outer solution; however, there is a cancellation at leading order when we substitute
$a$
and
$j$
from (2.18a,b
). To compute
$G$
we need to expand the flux expulsion solution systematically. This is most easily done by setting
$$\begin{eqnarray}\displaystyle a=-{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m}\text{e}^{-\text{i}m{\it\alpha}t+{\it\chi}}, & & \displaystyle\end{eqnarray}$$
where the complex function
${\it\chi}$
, of space and time, gives the effect of flux expulsion. The current and Lorentz force are then
$$\begin{eqnarray}\displaystyle j & = & \displaystyle -{\textstyle \frac{1}{2}}\text{i}b_{0}r^{m}\{\!m^{2}{\it\alpha}^{\prime 2}t^{2}+\text{i}m[{\it\alpha}^{\prime \prime }+(2m+1){\it\alpha}^{\prime }r^{-1}+2{\it\alpha}^{\prime }{\it\chi}^{\prime }]t\nonumber\\ \displaystyle & & \displaystyle -\,[{\it\chi}^{\prime \prime }+{\it\chi}^{\prime 2}+(2m+1){\it\chi}^{\prime }r^{-1}]\!\}\!\text{e}^{-\text{i}m{\it\alpha}t+{\it\chi}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle G & = & \displaystyle {\textstyle \frac{1}{4}}mb_{0}^{2}r^{2m}\{\!2m[{\it\alpha}^{\prime \prime }+(2m+1){\it\alpha}^{\prime }r^{-1}+{\it\alpha}^{\prime }({\it\chi}^{\prime }+{\it\chi}^{\prime \ast })]t\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}[({\it\chi}^{\prime \prime }-{\it\chi}^{\prime \prime \ast })+({\it\chi}^{\prime 2}-{\it\chi}^{\prime \ast 2})+(2m+1)\,({\it\chi}^{\prime }-{\it\chi}^{\prime \ast })r^{-1}]\!\}\!\text{e}^{{\it\chi}+{\it\chi}^{\ast }},\end{eqnarray}$$
without approximation.
Now to calculate
${\it\chi}$
we introduce
$T={\it\eta}^{1/3}t$
as the time scale on which flux expulsion occurs, and set
${\it\chi}={\it\chi}(r,T)$
. The exact equation for
${\it\chi}$
follows from (2.9) and is
$$\begin{eqnarray}\displaystyle \partial _{T}{\it\chi} & = & \displaystyle -m^{2}{\it\alpha}^{\prime 2}T^{2}-\text{i}m{\it\eta}^{1/3}[{\it\alpha}^{\prime \prime }+(2m+1){\it\alpha}^{\prime }r^{-1}+2{\it\alpha}^{\prime }{\it\chi}^{\prime }]T\nonumber\\ \displaystyle & & \displaystyle +\,{\it\eta}^{2/3}[{\it\chi}^{\prime \prime }+{\it\chi}^{\prime 2}+(2m+1){\it\chi}^{\prime }r^{-1}],\end{eqnarray}$$
where as usual the prime denotes an
$r$
-derivative at constant
$T$
(or
$t$
). A series approximation for
${\it\chi}$
can now be developed, with expansion parameter
${\it\eta}^{1/3}\ll 1$
. Explicitly
$$\begin{eqnarray}\displaystyle {\it\chi}(r,t)={\it\chi}_{0}(r,T)+{\it\eta}^{1/3}{\it\chi}_{1}(r,T)+\cdots , & & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\chi}_{0}=-{\textstyle \frac{1}{3}}m^{2}{\it\alpha}^{\prime 2}T^{3}=-{\textstyle \frac{1}{3}}m^{2}{\it\eta}{\it\alpha}^{\prime 2}t^{3}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\chi}_{1} & = & \displaystyle -{\textstyle \frac{1}{2}}\text{i}m[{\it\alpha}^{\prime \prime }+(2m+1){\it\alpha}^{\prime }r^{-1}]T^{2}+{\textstyle \frac{4}{15}}\text{i}m^{3}{\it\alpha}^{\prime 2}{\it\alpha}^{\prime \prime }T^{5}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle -{\textstyle \frac{1}{2}}\text{i}m[{\it\alpha}^{\prime \prime }+(2m+1){\it\alpha}^{\prime }r^{-1}]{\it\eta}^{2/3}t^{2}+{\textstyle \frac{4}{15}}\text{i}m^{3}{\it\eta}^{5/3}{\it\alpha}^{\prime 2}{\it\alpha}^{\prime \prime }t^{5}.\end{eqnarray}$$
To obtain the Lorentz torque at leading order (first line of terms in
$G$
in (4.17)) we need only
${\it\chi}_{0}$
, which then gives
$$\begin{eqnarray}\displaystyle G\simeq {\textstyle \frac{1}{2}}m^{2}b_{0}^{2}r^{2m}\left\{[{\it\alpha}^{\prime \prime }+(2m+1){\it\alpha}^{\prime }r^{-1}]t-{\textstyle \frac{4}{3}}m^{2}{\it\eta}{\it\alpha}^{\prime 2}{\it\alpha}^{\prime \prime }t^{4}\right\}\!\text{e}^{-(2/3)m^{2}{\it\eta}{\it\alpha}^{\prime 2}t^{3}}. & & \displaystyle\end{eqnarray}$$
Here we see that despite the initial quadratic growth with time
$t$
of the current
$j$
in (2.18a,b
), the Lorentz torque
$G$
grows only linearly via the first two terms (the latter term being small for
$t=O(1),{\it\eta}\ll 1$
). However, at the time of flux expulsion
$t=O({\it\eta}^{-1/3})$
all the terms in
$G$
in (4.23) are of the same order, and so all are important in computing the Lorentz force up to and during the destruction of field through flux expulsion.
4.3 Feedback through the Lorentz force: outer solution
We now discuss the effects of the Lorentz force on the flow: this features in the evolution of azimuthal velocity
$v$
,
$$\begin{eqnarray}\displaystyle \partial _{t}v=r^{-1}G+{\it\nu}({\it\Delta}_{0}v-r^{-2}v), & & \displaystyle\end{eqnarray}$$
which in turn gives the equation for the angular velocity gradient
${\it\alpha}^{\prime }$
,
$$\begin{eqnarray}\displaystyle \partial _{t}{\it\alpha}^{\prime }=(r^{-2}G)^{\prime }+{\it\nu}[r^{-1}({\it\Delta}_{0}v-r^{-2}v)]^{\prime }. & & \displaystyle\end{eqnarray}$$
It is this key quantity we shall use. Our working hypothesis is that the Lorentz force has the effect of flattening the angular velocity at radius
$r$
if the term
$(r^{-2}G)^{\prime }$
, evaluated kinematically and integrated over all time, is comparable to or bigger than
${\it\alpha}^{\prime }$
at that radius. The threshold to consider is then when
$$\begin{eqnarray}\displaystyle {\it\alpha}^{\prime }(r)\sim \int _{0}^{\infty }[r^{-2}G(r,t)]^{\prime }\,\text{d}t. & & \displaystyle\end{eqnarray}$$
Bearing in mind that the right-hand side is quadratic in
$b_{0}$
and we are interested in thresholds in terms of
$b_{0}$
, we set the key function we need as
$$\begin{eqnarray}\displaystyle h(r)=-{\it\alpha}^{\prime }(r)^{-1}\,b_{0}^{-2}\int _{0}^{\infty }[r^{-2}G(r,t)]^{\prime }\,\text{d}t. & & \displaystyle\end{eqnarray}$$
The minus sign here arises as in (4.25) it is the integral on the right-hand side which is being compared with reducing the angular velocity gradient from
${\it\alpha}^{\prime }(r)$
at
$t=0$
to zero at large times
$t$
. We will take up the discussion of the physics in § 5, but for the moment just focus on calculating
$h(r)$
for the outer solution and then the inner one.
First we take
$G$
as given from (4.23) with
$m=1$
now (to avoid unnecessary complexity) and obtain
$$\begin{eqnarray}\displaystyle [r^{-2}G]^{\prime }={\textstyle \frac{1}{2}}b_{0}^{2}[G_{0}(r)t+G_{1}(r){\it\eta}t^{4}+G_{2}(r){\it\eta}^{2}t^{7}]\text{e}^{-(2/3){\it\eta}{\it\alpha}^{\prime 2}t^{3}}, & & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle G_{0}(r)={\it\alpha}^{\prime \prime \prime }+3{\it\alpha}^{\prime \prime }r^{-1}-3{\it\alpha}^{\prime }r^{-2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{1}(r)=-4{\it\alpha}^{\prime }{\it\alpha}^{\prime \prime }({\it\alpha}^{\prime \prime }+{\it\alpha}^{\prime }r^{-1})-{\textstyle \frac{4}{3}}{\it\alpha}^{\prime 2}{\it\alpha}^{\prime \prime \prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{2}(r)={\textstyle \frac{16}{9}}{\it\alpha}^{\prime 3}{\it\alpha}^{\prime \prime 2}. & \displaystyle\end{eqnarray}$$
Next we integrate (4.28) from time zero to infinity using integration by parts and by setting the constant
$c_{0}$
defined by
$$\begin{eqnarray}\displaystyle c_{0}=\int _{0}^{\infty }t\exp (-t^{3})\,\text{d}t=3^{-2/3}{\rm\pi}\text{Hi}^{\prime }(0)\simeq 0.45137 & & \displaystyle\end{eqnarray}$$
(see Olver et al. Reference Olver, Lozier, Boisvert and Clark2010, § 9.12) to obtain
$$\begin{eqnarray}\displaystyle h(r)=-c_{1}{\it\eta}^{-2/3}{\it\alpha}^{\prime -7/3}\,\left[G_{0}(r)+{\it\alpha}^{\prime -2}G_{1}(r)+{\textstyle \frac{5}{2}}{\it\alpha}^{\prime -4}G_{2}(r)\right]\!,\quad c_{1}\equiv 2^{-5/3}3^{2/3}c_{0}, & & \displaystyle\end{eqnarray}$$
or, with the formulae for the
$G_{j}(r)$
in (4.29) we obtain
$$\begin{eqnarray}\displaystyle h(r)=c_{1}{\it\eta}^{-2/3}{\it\alpha}^{\prime -7/3}\left[{\textstyle \frac{1}{3}}{\it\alpha}^{\prime \prime \prime }+{\it\alpha}^{\prime \prime }r^{-1}+3{\it\alpha}^{\prime }r^{-2}-{\textstyle \frac{4}{9}}{\it\alpha}^{\prime \prime 2}{\it\alpha}^{\prime -1}\right]\!. & & \displaystyle\end{eqnarray}$$
This (implausible looking) expression is correct for any
${\it\alpha}(r)$
. However, note that when we approach the origin,
${\it\alpha}(r)$
expands as in (2.20) and then
$$\begin{eqnarray}\displaystyle h(r)\simeq c_{2}{\it\eta}^{-2/3}{\it\alpha}_{1}^{-4/3}r^{-10/3},\quad c_{2}\equiv 2^{2}3^{-4/3}c_{0}. & & \displaystyle\end{eqnarray}$$
This is valid in the overlap region
$\mathscr{L}\ll r\ll 1$
where both inner and outer solutions are valid (referring to figure 3). Before using this, we proceed to look at the feedback in the inner solution.
4.4 Feedback through the Lorentz force: inner solution
Near the origin where the above formulation breaks down, we need to use the inner solution (4.1a−c )–(4.6a,b ) instead, which yields
$$\begin{eqnarray}\displaystyle G=mb_{0}^{2}|g|^{2}k^{2-2m}{\it\rho}^{2m}(f+f^{\ast })[(m+1)-\text{i}(f-f^{\ast }){\it\rho}^{2}]\text{e}^{-\text{i}(f-f^{\ast }){\it\rho}^{2}} & & \displaystyle\end{eqnarray}$$
and hence (with
$m=1$
) we find
$$\begin{eqnarray}\displaystyle (r^{-2}G)^{\prime }=2b_{0}^{2}|g|^{2}k^{3}(f^{2}-f^{\ast 2})[-3\text{i}{\it\rho}-(f-f^{\ast }){\it\rho}^{3}]\text{e}^{-\text{i}(f-f^{\ast }){\it\rho}^{2}}. & & \displaystyle\end{eqnarray}$$
Now to obtain
$h(r)$
in (4.27) we divide by
$b_{0}^{2}$
and by
${\it\alpha}^{\prime }=2{\it\alpha}_{1}k^{-1}{\it\rho}$
and integrate over all time. With the use of (4.1a−c
) we obtain
$$\begin{eqnarray}\displaystyle h({\it\rho})=\mathscr{V}^{-2}\int _{0}^{\infty }|g|^{2}(f^{2}-f^{\ast 2})[3\text{i}+(f-f^{\ast }){\it\rho}^{2}]\text{e}^{-\text{i}(f-f^{\ast }){\it\rho}^{2}}\,\text{d}{\it\tau}. & & \displaystyle\end{eqnarray}$$
We have taken the liberty of thinking of
$h$
now as a function of
${\it\rho}=kr$
and taken the integral over
${\it\tau}=pt$
. We cannot evaluate this analytically, except for large
${\it\rho}$
when the approximation (4.9a−c
) is valid throughout the time range giving the leading contribution to the integral, with
$$\begin{eqnarray}\displaystyle h({\it\rho})\simeq c_{3}\mathscr{V}^{-2}{\it\rho}^{-10/3},\quad c_{3}\equiv 2^{8/3}3^{-4/3}c_{0}. & & \displaystyle\end{eqnarray}$$
This is valid in the overlap region
$1\ll {\it\rho}\ll \mathscr{L}^{-1}$
and making use of the definitions of
$p$
,
$k$
and
${\it\rho}$
in (4.1a−c
) we recover (4.33a,b
) as we must. More generally we can plot
$\mathscr{V}^{2}h({\it\rho})$
in (4.36) against
${\it\rho}$
as in figure 11 to give a universal curve for the Lorentz feedback at the centre of a vortex with a general, smooth angular velocity profile (that is, with
${\it\alpha}_{1}\neq 0$
).
5 Scaling laws and related information
We can use the results in the last section to make a number of predictions for scaling laws. We begin with crude estimates and then give more precise versions. Our aim is to evaluate the accumulated effect of the Lorentz force at a given radius
$r$
in a kinematic regime and compare this with
${\it\alpha}^{\prime }(r)$
at that radius. We use
${\gtrsim}$
and
${\lesssim}$
to denote inequality up to a constant of order unity. Recall that we are comparing
$$\begin{eqnarray}\displaystyle {\it\alpha}^{\prime }(r)\quad \text{and}\quad \int _{0}^{\infty }[r^{-2}G(r,t)]^{\prime }\,\text{d}t & & \displaystyle\end{eqnarray}$$
at a given radius
$r$
or range of radii. From the definition of
$h(r)$
in (4.27) the integrated effect of the Lorentz force at radius
$r$
is sufficient to suppress the differential rotation if
$$\begin{eqnarray}\displaystyle b_{0}\sim h(r)^{-1/2}. & & \displaystyle\end{eqnarray}$$
Also to fix ideas we plot
$h(r)$
in figure 12 for the Gaussian vortex (2.16a,b
) for a range of
${\it\eta}$
values using the outer expansion (4.32) and the inner expansion (4.36).
Figure 12. Plot of
${\it\eta}^{2/3}h(r)$
(solid) given in (4.32) for the Gaussian vortex against
$r$
with (a) linear scales, (b) log–log scales (taking the absolute value). In (a) the approximation (4.33a,b
) is shown dashed: it is also present in (b) but is somewhat overlapped with dotted curves showing
${\it\eta}^{2/3}h({\it\rho})$
from (4.36) plotted against
$r$
for
${\it\eta}=10^{6}$
to
$10^{9}$
, reading up the curves.
Core formation threshold: first consider
$b_{0}\gtrsim b_{core}$
; then at all radii
$r$
the left-hand side in (5.2) dominates. The Lorentz force (estimated by kinematic evolution) is strong enough to modify the flow field and suppress differential rotation and so flux expulsion. This represents the upper threshold for core formation. For greater fields, elastic forces dominate and prevent the onset of flux expulsion.
First of all, consider a basic estimate. The location where flux expulsion would commence is at the radius
$r^{\dagger }$
of order unity where the differential rotation
${\it\alpha}^{\prime }$
is maximised, and it is here or nearby where the maximum in
$b_{core}$
will be realised. At order-one values of radius we can use (4.32), which gives
$h(r)=O({\it\eta}^{-2/3})$
and results in the threshold
$$\begin{eqnarray}\displaystyle b_{core}=O({\it\eta}^{1/3}), & & \displaystyle\end{eqnarray}$$
as observed in figure 10 (upper dataset).
To improve on this, it would make sense to minimise
$h(r)$
over all radii in view of (5.2) so as to maximise the field for flux expulsion to take place. However,
$h(r)$
has a zero crossing, and we should note that the threshold for core formation is linked to Alfvén wave generation, a process very far from the kinematic model we are using here. We do suggest instead a more precise estimate which brings in factors of
${\it\alpha}_{1}$
, taken by substituting
$r^{\dagger }$
in (4.33a,b
) to give
$$\begin{eqnarray}\displaystyle b_{core}({\it\eta})=c_{2}^{-1/2}{\it\eta}^{1/3}{\it\alpha}_{1}^{2/3}(r^{\dagger })^{5/3}. & & \displaystyle\end{eqnarray}$$
This is shown on figure 10 with good agreement, fortuitously good given that it is only correct up to a constant of order unity.
Dynamical core threshold: now let us go to the other extreme. If
$b_{0}\lesssim b_{dynam}$
then for all radii the right-hand side of (5.2) dominates and the field will remain kinematic at leading order. A core will form through flux expulsion and will shrink to an exponentially decaying remnant at the origin following BBG. We thus need to look at where the maximum of
$h(r)$
is realised. From (4.33a,b
) (see also figures 11, 12), from the point of view of the outer solution, we can increase
$h(r)$
by reducing
$r$
. This represents the physical fact that the flow field is ‘naturally’ fairly flat near the origin, the field is long-lived, and the Lorentz force required to flatten it completely becomes vanishing small. However, this cannot continue indefinitely and we can argue this in two ways. First of all it is clear that this increasing
$h(r)$
in (4.33a,b
) (in the overlap region) must be cut off at radii
$r\sim \mathscr{L}=O({\it\eta}^{1/4})$
in (4.1a−c
) where the overlap region ends and the inner solution really takes over. Substituting this in (4.33a,b
) gives the maximum of
$h(r)$
and the threshold estimated as
$$\begin{eqnarray}\displaystyle b_{dynam}=O({\it\eta}^{3/4}), & & \displaystyle\end{eqnarray}$$
in agreement with figure 10 (lower dataset).
With more precision and elegance, we can move to the inner solution and figure 11(a), which shows that
$\mathscr{V}^{2}h$
is maximised at
$$\begin{eqnarray}\displaystyle c_{4}^{-2}\equiv h(0)=\int _{0}^{\infty }3\text{i}|g|^{2}(f^{2}-f^{\ast 2})\,\text{d}{\it\tau}\simeq 0.91732 & & \displaystyle\end{eqnarray}$$
(obtained numerically) to give
$$\begin{eqnarray}\displaystyle b_{dynam}=c_{4}\mathscr{V}=c_{4}(8{\it\eta}^{3}m{\it\alpha}_{1})^{1/4}, & & \displaystyle\end{eqnarray}$$
from (4.1a−c
). Thus we link the threshold field to the velocity scale
$\mathscr{V}$
based on the inner solution, which has the crucial
${\it\eta}^{3/4}$
dependence. This is plotted in figure 10 (lower dotted line), showing good agreement up to a modest constant.
Onset of Lorentz force: finally suppose that
$b_{0}$
lies in the range
$b_{dynam}\lesssim b_{0}\lesssim b_{core}$
. Then flux expulsion occurs at a radius
$r^{\dagger }$
, a core is formed and a wave of flux expulsion spreads inwards. The Lorentz force does, however, become important at a radius
$r_{\ast }(b_{0})$
and time
$t_{\ast }(b_{0})$
given by
$$\begin{eqnarray}\displaystyle b_{0}=h(r_{\ast })^{-1/2},\quad t_{\ast }=\left[{\textstyle \frac{1}{3}}{\it\eta}{\it\alpha}^{\prime }(r_{\ast })^{2}\right]^{-1/3}. & & \displaystyle\end{eqnarray}$$
These functions would depend on the original flow profile, but when
$r_{\ast }\ll 1$
the approximation (4.33a,b
) in the overlap region comes into play (also
${\it\alpha}^{\prime }(r_{\ast })\simeq 2{\it\alpha}_{1}r_{\ast }$
) and we can estimate that
$$\begin{eqnarray}\displaystyle r_{\ast }\simeq c_{2}^{3/10}b_{0}^{3/5}{\it\eta}^{-1/5}{\it\alpha}_{1}^{-2/5},\quad t_{\ast }\simeq c_{5}{\it\eta}^{-1/5}b_{0}^{-2/5}{\it\alpha}_{1}^{-2/5},\quad c_{5}=2^{-2/3}3^{1/3}c_{2}^{-1/5}.\quad & & \displaystyle\end{eqnarray}$$
For example for
${\it\eta}=10^{-7}$
and
$b_{0}=10^{-5}$
we obtain
$r_{\ast }\simeq 0.16$
and
$t_{\ast }\simeq 2.3\times 10^{4}$
, in agreement with figure 8. Note that at the lower and upper limits of the range
$b_{dynam}\lesssim b_{0}\lesssim b_{core}$
, the core size is of order
$\mathscr{L}$
in (4.1a−c
) and
$r^{\dagger }$
, respectively.
6 Discussion
We have studied some of the effects of dynamical feedback on flux expulsion using a quasi-linear model, using both numerical simulation and theory based on the kinematic picture. We have identified two thresholds: returning to our original dimensional formulation before (2.15a−f
) and measuring magnetic field in velocity units, the first threshold is
$b_{core}\sim UR_{m}^{-1/3}$
below which flux expulsion still operates, cutting field lines at the radius
$r^{\dagger }$
and leaving a magnetic core within. The second is
$b_{dynam}\sim UR_{m}^{-3/4}$
below which the field evolves as in the kinematic regime. Between the two thresholds a magnetic core is formed in which the flow field near the origin is modified to be solid body rotation at leading order, so halting the diffusive decay of the core. In the range between
$b_{core}$
and
$b_{dynam}$
the core radius scales as
$L(b_{0}/U)^{3/5}R_{m}^{1/5}$
. Recall that
$L$
is the spatial scale of the flow and
$U$
the velocity scale of the flow. However, the actual key parameters at the lower threshold are the scale
$L$
of the flow and the magnitude of
${\it\alpha}_{1}$
in dimensional units of
$(\text{length}^{2}\times \text{time})^{-1}$
; in other words it is appropriate to take
$U=L^{3}{\it\alpha}_{1}$
.
In each case diffusive processes are key and the results depend sensitively on the magnetic diffusivity – only the cutting of field lines by diffusion can halt the increase of Lorentz force as the field reduces in scale. With this in mind, we have made careful estimates based on the kinematic solutions, inner and outer, and we note that cruder arguments could easily lead to incorrect conclusions. For example, although the magnetic field grows linearly with time, the Lorentz force term also grows linearly in (4.23) and not quadratically as might be suggested by a simplistic estimate
$\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{b}\sim b^{2}/L$
,
$L$
being the scale of the eddy, this overestimating the variation of
$\boldsymbol{b}$
along field lines. Even worse would be to estimate the Lorentz force as
$jb=O(t^{3})$
. Likewise in (4.26) it is necessary to calculate the effect of
$G$
integrated from time zero up to the time when flux expulsion occurs (similar remarks apply in the kinematic regime, as discussed in Moffatt & Kamkar Reference Moffatt, Kamkar and Soward1983). It is also important to look at the effect on the differential rotation, not the velocity or angular velocity. Finally to pick up the lower threshold
$b_{dynam}$
the scaling structure of the inner solution from BBG is needed, in particular relating the field magnitude to
$\mathscr{V}=O(UR_{m}^{-3/4})$
. We also note that any replacement of the true magnetic diffusion term involving the Laplacian, by some hyperdiffusion or similar cutoff, would change these scaling laws. Any use of hyperdiffusion in the induction equation must be treated with caution: although the change may have a minor impact on small-scale fields at any moment, here it would have a significant impact on the large-scale, long-time evolution.
Future directions of research could include generalising the geometry, for example to an axisymmetric eddy localised in three dimensions and a magnetic field with initially an arbitrary orientation with respect to the eddy. It would also be interesting to study other regimes of the Reynolds numbers. We have taken only
$R_{e}\gg R_{m}\gg 1$
, although we expect our results to have wider applicability, we think at least up to
$R_{e}\sim R_{m}$
. The ordering
$R_{e}\gg R_{m}\gg 1$
is relevant in typical astrophysical and geophysical contexts, and so the modelling could be broadly relevant to the formation of magnetic fields in the early life of some astrophysical objects, for example the relict magnetic field likely to lie in the radiative zone of the Sun (e.g. Mestel & Weiss Reference Mestel and Weiss1986). At least it indicates the importance of taking into account small-scale reconnection processes that depend on molecular transport coefficients, even in the formation of large-scale fields, this also being the point originally stressed by Vainshtein & Cattaneo (Reference Vainshtein and Cattaneo1992) in the context of dynamo theory.
A related problem would be to consider an initial two-dimensional turbulent flow containing eddies or vortices, on a range of length scales: which of these would develop dynamical magnetic cores, and what would be their distribution? For what field threshold would the evolution over all length scales be entirely kinematic? Finally, it would be valuable to study dynamical flux expulsion numerically by means of full simulations in unbounded geometry and explore the limitations of the quasi-linear approximation as set out here. Simulations with
$R_{e}\gg R_{m}\gg 1$
can be undertaken efficiently using contour–spectral methods (Dritschel & Tobias Reference Dritschel and Tobias2012) while regimes with
$R_{e}\sim R_{m}\gg 1$
would be best simulated with standard pseudo-spectral codes.
Acknowledgements
This paper is dedicated to K. Bajer (Warsaw University), an inspiring colleague who contributed greatly to the field of fluid mechanics and who sadly passed away in August 2014. We thank H. K. Moffatt, A. M. Soward and N. O. Weiss for useful discussions and the referees for helpful comments leading to improvements in the manuscript.
