%% Copyright 2007, 2008, 2009 Elsevier Ltd
%%
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%\documentclass[preprint,authoryear,11pt]{elsarticle}
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}
%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%%%%\documentclass[final,authoryear,1p,times,fleqn]{elsarticle}
%% \documentclass[final,authoryear,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,authoryear,3p,times]{elsarticle}
%\documentclass[final,authoryear,3p,times,twocolumn]{elsarticle}
%\documentclass[final,authoryear,5p,times,fleqn]{elsarticle}
\documentclass[final,numbers,5p,times,fleqn]{elsarticle}
%\documentclass[final,authoryear,5p,times,twocolumn,fleqn]{elsarticle}
%\documentclass[final,numbers,5p,times,twocolumn]{elsarticle}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%\usepackage{lineno}
%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:
%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon (default)
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   authoryear - selects author-year citations (default)
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%   longnamesfirst  -  makes first citation full author list
%%
%% \biboptions{longnamesfirst,comma}

%\biboptions{numbers}
\usepackage{amssymb}

\journal{Icarus}

\usepackage{epsfig}
\usepackage{psfrag}
\usepackage{amsfonts}          % |
\usepackage{amsmath}           % |
\usepackage{ulem}
\usepackage{xspace}
%\usepackage{widetext}
%\usepackage{slashbox}
\usepackage{rotating}
\usepackage[a4paper,ps2pdf,colorlinks=true,linkcolor=red,urlcolor=blue,breaklinks=true]{hyperref}
\usepackage[yyyymmdd]{datetime}
\usepackage{natbib}


%-------- Abbreviation
\providecommand{\ie}{i.e.\xspace}
\providecommand{\eg}{e.g.\xspace}
%-------- Symbolic names for complicated things
%\providecommand{\A}{\textbf{aDD}\xspace}
%\providecommand{\B}{\textbf{uTL}\xspace}
%\providecommand{\C}{\textbf{sTL}\xspace}
%\providecommand{\Done}{\ensuremath{I_1}\xspace}
%\providecommand{\Dtwo}{\ensuremath{I_2}\xspace}

\input{custom}
\usepackage{charter}


\usepackage{lipsum}
\usepackage{color}

\newcommand{\mypsfrag}[2]{\psfrag{#1}{\small{#2}}}

\allowdisplaybreaks

% Alter some LaTeX defaults for better treatment of figures:
    % See p.105 of ``TeX Unbound'' for suggested values.
    % See pp. 199-200 of Lamport's ``LaTeX'' book for details.
    %   General parameters, for ALL pages:
\renewcommand{\topfraction}{0.99}% max fraction of floats at top
\renewcommand{\bottomfraction}{0.99}% max fraction of floats at bottom
%   Parameters for TEXT pages (not float pages):
\setcounter{topnumber}{2}
\setcounter{bottomnumber}{2}
\setcounter{totalnumber}{2}     % 2 may work better
\setcounter{dbltopnumber}{2}    % for 2-column pages
\renewcommand{\dbltopfraction}{0.99}% fit big float above     2-col. text
\renewcommand{\textfraction}{0.01}% allow minimal text w. figs
%   Parameters for FLOAT pages (not text pages):
\renewcommand{\floatpagefraction}{0.99}% require fuller float pages
% N.B.: floatpagefraction MUST be less than topfraction !!
\renewcommand{\dblfloatpagefraction}{0.99}% require fuller float    pages

\newlength{\shortwidth}
\setlength{\shortwidth}{\textwidth}
\addtolength{\shortwidth}{-20mm}

\newlength{\textwidtha}
\setlength{\textwidtha}{\textwidth}
\addtolength{\textwidtha}{-1mm}

\makeatletter
\setlength{\@mathmargin}{0cm}
\makeatother

\raggedbottom

% labels of the comments
%\makeatletter
%\input{Response.aux}
%\makeatother

\begin{document}

\begin{frontmatter}
%% Title, authors and addresses

\title{Remarks on anelastic spherical dynamos and their transition to a buoyancy-dominated regime}

\author[label1,label2,label3]{Radostin D.~Simitev}
\ead{Radostin.Simitev@glasgow.ac.uk}
\author[label2,label4]{Friedrich H.~Busse}
\author[label3]{Alexander G.~Kosovichev}
\address[label1]{School of Mathematics and Statistics, University of
  Glasgow -- Glasgow G12 8QW, UK, EU}
\address[label2]{Department of Earth and Space Sciences, University of
  California, Los Angeles -- Los Angeles, CA 90095, USA}
\address[label3]{W.W. Hansen Experimental Physics Laboratory,
    Stanford University -- Stanford CA 94305, USA}
\address[label4]{Institute of Physics, University of Bayreuth --
  Bayreuth 95440,  Germany}

\begin{abstract}
We present a set of convective dynamo simulations in rotating
spherical fluid shells based on an anelastic approximation of
compressible fluids. 
The simulations extend into a ``buoyancy-dominated'' regime where the
buoyancy forcing is dominant while the Coriolis force is no longer
balanced by pressure gradients and strong retrograde differential
rotation develops as a result. Dynamos in this regime are strongly
dominated by dipole components but at the same time their magnetic
energies are relatively small compared to the corresponding kinetic
energies of the flow. Despite being relatively weak the self-sustained
magnetic fields are able to reverse the direction of differential
rotation to prograde and give rise to some similarities with Solar
convection. 
\end{abstract}

\begin{keyword}
Rotating spherical shells \sep
Convection \sep
Dynamo \sep
Anelastic approximation \sep
Differential rotation \sep
Solar convection zone
\end{keyword}
\end{frontmatter}
%\linenumbers


\section{Introduction}

\paragraph{Motivation and hypotheses}
The following facts about the Solar structure and dynamics are
established from observations or are well-accepted assumptions in the
research community.
(a) Solar convection is very vigorous and turbulent with
characteristic surface velocities of the order ... \cite{xx}; (b) Solar
rotation is nonuniform but relatively slow at an average of 27 days
per revolution, e.g. \cite{}; (c) Molecular viscosity of Solar plasma is nearly
negligible with a characteristic value of the order of ... \cite{xx};
(d) The energy of the global magnetic field of the Sun is small
compared to the kinetic energy of plasma motions and differential
rotation does not vary much with the Solar cycle \cite{xx}.
This suggests that Solar convection is in a ``buoyancy-dominated
regime'' where buoyancy forces dominate over viscous, Lorentz and
Coriolis forces in the momentum balance. 
Thus, it is of interest to study convection and dynamo action in the
buoyancy-dominated regime. 

In fact, non-magnetic convection in the buoyancy- dominated regime has
already attracted attention. The regime was first identified by Fokal
\cite{xx} and Gilman \cite{Gilman76} and more recently it was studied
by Aurnou et al.~\cite{Aurnou2007} and Gastine et
al.~\cite{Gastine2013a}. Bousinesq models were used in \cite{Gilman76}
and \cite{Aurnou2007} and an anelastic model was used in
\cite{Gastine2013a}. These studies consistently found that due to
vigorous mixing angular momentum is homogenized within the whole
volume of the shell and this leads to strong retrograde zonal flow in
the equatorial region. Retrograde differential rotation in the
equatorial region is antisolar. So the question arises how can the
Solar convection zone be in a buoyancy-dominated regime and at the
same time sustain prograde equatorial rotation. This is the question
we attempt to address in this paper. To this end, we demonstrate that 
demonstrate that relatively weak self-sustained field will act
to reverse the retrograte differential rotation thus achieving a
buoyancy-dominated regime with a solar-like zonal flow. 

It is well-established \cite{Bu03,Simitev2005a} that the main effect
of self-sustained magnetic field on convection is to suppress
differential rotation. Dynamo action beyond the transition from
rotation-dominated regime to buoyancy-dominated regime has not been
previously studied.  

Buoyancy-dominated regime of convection may be relevant to Neptune and
Uranus. There is evidence that the magnetic field of Uranus is
non-dipolar (provide references)   and that the differential rotation
is retrograde at the  equator. Retrograde differential rotation
naturally develops in the  buoyancy-dominated regime and it is of
interest whether the magnetic  fields that are generated resemble that
of Neptune and Uranus. 


\paragraph{Aims and plan of paper}
In section \ref{sec.2} we introduce the mathematical model based on
the anelastic approximation and briefly touch on some of the numerical
implementation details. In this paper we present first results from a
new anelastic code. For this reason, we include in section \ref{sec.3}
a brief validation test against the benchmark solutions recently
presented in \cite{JonesBench}. In secion \ref{sec.4} we comment in
passing on some similarities between Boussinesq and anelastic
convection. Section \ref{sec.5} is the focus of the paper. Here we
discuss the properties of dynamos in the buoyancy-dominated regime of
convection. We conclude in section \ref{sec.6} with a summary of our
main results and discussion of their implications. We also outline
questions for future research. 


\section{Mathematical model and numerical method}
\label{sec.2}

\subsection{Anelastic equations}

We consider an electrically conducting, self-gravitating, perfect gas
confined to a spherical shell. 
The shell rotates with a fixed angular velocity $\Omega \k$ about the
vertical axis and an entropy contrast $\Delta S$ is imposed between
its inner and outer surfaces. 

Assuming a gravity field proportional to $1/r^2$, a hydrostatic
polytropic reference state exists of the form  
\begin{gather}
\rhobar = \rho_c\zeta^n, \quad \Tbar=T_c\zeta, \quad \Pbar = P_c
\zeta^{n+1}, \quad
\zeta= c_0+c_1 d/r,
\end{gather}
with parameters $c_0=(2\zeta_o-\eta-1)/(1-\eta)$,
$c_1=(1+\eta)(1-\zeta_o)/(1-\eta)^2$, $\zeta_o=(\eta+1)/(\eta
\exp(\Nrho/n)+1)$.
%$\zeta_i=(1+\eta-\zeta_o)/\eta$. 
The parameters $\rho_c$, $P_c$ and $T_c$ are reference values of
density, pressure and temperature at the middle of the shell, and the
gas polytropic index $n$, the density scale height 
$N_\rho$ and the shell thickness ratio $\eta$ are defined further below. 
%and $N_\rho=\ln(\rhobar(r_i)/\rhobar(r_o))$. 
%Here,  $n$ is the polytropic index, $\eta = r_i/r_o$ is the radius
%ratio, 
Convection and magnetic
field generation set in for sufficiently large values of the
entropy contrast, $\Delta S$, and can be described by the 
equations of continuity, momentum, energy and magnetic flux. In the
anelastic approximation \citep{Gough1969,JonesBench} these equations
take the form   
\bs
\label{gov}
\begin{gather}
\label{gov.01}
\n\.\rhobar\u =0, \quad \n\.\B =0, \\
\label{gov.02}
\partial_t \u + (\n\x\u)\x\u \\
\hspace*{6.5mm}=-\n\Pi -\tau(\k\x\u)+\f{\R}{\Pr}\f{S}{r^{2}}\ru
+ \v F_\nu + \f{1}{\rhobar} (\n\x\B)\x\B, \nonumber \\
%
\label{gov.03}
\partial_t S + \u\.\n S \\
\hspace*{6.5mm}= \f{1}{\Pr \rhobar\Tbar} \n\.\kapbar\rhobar\Tbar \n S
  + \f{c_1 \Pr}{\R \Tbar}\left(Q_\nu+ \f{1}{\Pm\rhobar} Q_j\right) 
%+ H, 
\nonumber\\
%
\label{gov.04}
\partial_t \B = \n\x(\u\x\B)+\Pm^{-1} \n^2 \B,
\end{gather}
\es
where $\u$ is the velocity, $\B$ is the magnetic flux density, $S$
is the entropy and $\nabla \Pi$ includes all terms that can be
written as gradients. The viscous force, and the viscous and Joule
heating, 
\begin{gather}
\label{not}
\v F_\nu = \frac{\rho_c}{\rhobar}\n\.\v {\hat S}, \quad Q_\nu=\v{\hat S}:\v e,  \quad Q_j=(\n\x\B)^2, 
\end{gather}
are defined in terms of the deviatoric stress tensor  
\begin{gather}
\label{tens}
\hat S_{ij}=2\nubar\rhobar(e_{ij}-e_{kk}\delta_{ij}/3), \quad e_{ij}=(\partial_i u_j +\partial_j u_i)/2,
\nonumber
\end{gather}
where double-dots (\textbf{:}) denotes the Frobenius inner product.
%and the driving term 
%\begin{gather}
%\label{drive}
%H=-\f{c_1 n d}{r^2
%  \rhobar\Tbar}\f{\rhobar(r_i)\rhobar(r_o)}{\rhobar(r_o)-\rhobar(r_i)} \ru\.\u.
%\end{gather}
We assume that the viscosity and the entropy diffusivity vary in
radius as $\nubar(r)= \nu_c  \rhobar^{k}$ and $\kapbar(r) = \kappa_c
\rhobar^{k}$ with some negative power $k\le 0$, where $\nu_c$ and
$\kappa_c$ are their reference values at midshell. The governing
equations have been nondimensionalised using the thickness of the  shell $d=r_o-r_i$
as a unit of length, $d^2/\nu_c$ as a unit of time, $\Delta S$ as a
unit of entropy, $\nu_c \sqrt{\mu_0\rho_c}/d$ as a unit of magnetic
induction, $\rho_c$ as a unit of density and $T_c$ as a unit of
temperature. Here, $r_i$ and $r_o$ are the inner and the outer radius,
$\lambda$ and $\mu_0$ are the magnetic diffusivity and permeability,
respectively. The system is then characterized by 
eight dimensionless parameters: the radius ratio, the polytropic
index of the gas, the density scale number, the radial dependence power, 
the Rayleigh number, the ordinary and the magnetic Prandtl numbers and
the Coriolis number, defined as   
\begin{gather}
\eta=r_i/r_o, \quad
n, \quad
N_\rho=\ln\Big(\rhobar(r_i)/\rhobar(r_o)\Big), \quad
k,
\nonumber \\
%\R=\f{GMd\Delta S}{c_p \nu_c \kappa_c}, \quad 
\R=\f{c_1 T_c d^2 \Delta S}{\nu_c \kappa_c}, \quad 
\Pr=\f{\nu_c}{\kappa_c}, \quad
\Pm=\f{\nu_c}{\lambda},\quad 
\tau = \f{2\Omega d^2}{\nu_c},
\label{numbers}
\end{gather}
respectively.
%respectively, where $G$ is the gravitational constant, $M$ is the mass
%within the core of the shell, $c_p$ is the specific heat at constant
%pressure (these can be avoided if $R$ is defined as $R=c_1T_cd^2\Delta 
%S/\nu_c \kappa$),

Since the mass flux $\rhobar \u$, and the magnetic flux
density $\vec B$ are solenoidal vector fields, it is advantageous to
employ a decomposition in poloidal and toroidal components,  
\bs
\begin{gather}
\rhobar \vec u = \nabla \times ( \nabla \times \ru r\pol) + \nabla \times
\ru r^2 \tor, \\
\vec B = \nabla \times  ( \nabla \times \ru \polB) + \nabla \times
\ru \torB,
\end{gather}
\es
where $\ru$ is the radial unit vector, $r$ is its length, $\pol$,
$\tor$, $\polB$ and $\torB$ are the poloidal and toroidal scalars of
the momentum and magnetic field, respectively.
Equations~\eqref{gov.01} are then satisfied by construction. Scalar
equations for $v$ and $w$ are obtained, and effective pressure
gradients are eliminated by taking $\ru\cdot\n\x\n\x$ and
$\ru\cdot\n\x$ of equation \eqref{gov.02}. Similarly, equations for
$h$ and $g$ are obtained by taking $\ru\cdot \n\x$ and $\ru\cdot$ of
equation \eqref{gov.04}. The resulting poloidal-toroidal equations are
somewhat lengthy and will not be listed here \red{(or may be listed in
an appendix)}. A minor disadvantage of this representation is that a
fourth-order poloidal equation is obtained.  


\subsection{Boundary conditions}

We explore various assumptions for the boundary conditions imposed on
velocity, entropy and magnetic flux at the surface of the shell. The
alternatives are listed below and the particular choice will be
specified in each individual case. At a no-slip, impenetrable boundary
we impose
\begin{gather}
\pol = 0, 
\quad 
\dd{r}{} \pol = 0, 
\quad 
\tor=0.
\end{gather}
At a stress-free, impenetrable boundary we require
\begin{gather}
\pol = 0,
\quad
\partial_r^2 \pol - \f{\rhobar'}{\rhobar r}\dd{r}{} (r\pol) = 0, 
\quad
\partial_r \tor - \f{\rhobar'}{\rhobar} \tor = 0.
\end{gather}
Values of the entropy may be fixed at the boundaries, then 
\begin{gather}
S=1 \text{ at } r=r_i, \qquad S=0 \text{ at } r=r_o.
\end{gather}
Alternatively, the entropy flux may be specified at the top,
\begin{gather}
\partial_r S=0 \text{ at } r=r_o.
\end{gather}
Boundary conditions for the magnetic field may be derived from the
assumption of an electrically insulating external region. The
poloidal function $h$ is then matched to a function $h^{(e)}$, which
describes an external  potential field,
\begin{gather}
\torB =0, \quad  \polB-\polB^{(e)} = 0, \quad \partial_r ( \polB-\polB^{(e)})=0.
\end{gather}
Alternatively, a perfectly conducting external region may be assumed,
requiring 
\begin{gather}
\partial_r \torB =0, \quad  \polB=0.
\end{gather}
Finally, the, so-called, ``pseudo-vacuum'' condition offers another choice,
\begin{gather}
\torB =0, \quad  \partial_r \polB=0.
\end{gather}


\subsection{Numerical method}
For the numerical solution of the problem we have adapted the
pseudo-spectral method described by \cite{Tilgner1999}. 
The scalar unknowns $\pol$, $\tor$, $\polB$, $\torB$ and $S$, are
expanded in Chebychev polynomials $T_p$ in the radial direction $r$, and
in spherical harmonics in the angular directions $(\theta,\varphi)$ e.g.,
\begin{equation}
v = \sum\limits_{l=1}^{N_l} \sum\limits_{m=-l}^{l}
\sum\limits_{p=0}^{N_r} V_{lp}^m (t)\, T_p\Big(x(r)\Big)\, P_l^m ( \cos \theta )\, \exp(im \varphi),
\end{equation}
where $P_l^m$ denotes the associated Legendre functions, 
$x(r)=2(r-r_i)-1$, and $N_l$ and $N_r$ are truncation parameters.
A system of equations for the coefficients in these expansions is
obtained by a combination of a Galerkin spectral projection of the
governing equations in the angular directions and a collocation
constraint in radius. Computation of nonlinear terms in spectral space is expensive,
so nonlinear products and the Coriolis term are computed in physical
space and then projected to spectral space at every time
step. A standard 3/2-dealiasing in $\theta$ and $\varphi$ is used
at this stage. 
A hybrid of a Crank-Nicolson scheme for the diffusion terms and a
second order Adams-Bashforth scheme for the nonlinear terms is used
for integration in time.

A range of numerical resolutions has been used in this study varying
from ($N_r=61$, $N_l=96$)  in less demanding cases to ($N_r=121$,
$N_l=144$) in more strongly stratified or turbulent runs. Correspondingly, the
physical gridpoints on which non-linear terms are evaluated have been varied up to
$N_r=121$, $N_\theta=216$, $N_\varphi=437$.
\input{tab01}
\input{fig01}    
\begin{figure*}
\begin{center}
\psfrag{om}{$\omega_c$}
\psfrag{Rc}{$R_c$}
\psfrag{mc}{$m_c$}
\psfrag{P}{$P$}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/01Calcs/2013.07.16.linear/p1t2
\epsfig{file=../01Calcs/2013.07.16.linear/p1t2/P_RaMin.eps,height=0.32\textwidth,clip=}
\epsfig{file=../01Calcs/2013.07.16.linear/p1t2/P_mCrit.eps,height=0.32\textwidth,clip=}
\epsfig{file=../01Calcs/2013.07.16.linear/p1t2/P_woCrit.eps,height=0.32\textwidth,clip=}
\end{center}
\vspace*{-6mm}
\caption{
(a) The critical Rayleigh number $R_c$, (b) the critical wave number
$m_\text{crit}$ and (c) the critical frequency for onset of convection
as functions of the Prandtl number $P$ for $\eta=0.65$, $\tau=100,
200, 300, 500, 1000, 2000, 3000$, $n=2$, (red dashed lines) and
$N_\rho=0.1, 1, 2, 3, 4, 5, 6$ (black solid lines) both increasing from
bottom to top. The values are calculated using the quasi-geostrophic
approximation of Busse and Simitev (2013).
}
\label{f00100}
\end{figure*}

\subsection{Diagnostic output quantities} 

Our numerical convection and dynamo solutions are characterized
by their kinetic and magnetic energy and heat transport given by a
Nusselt number.  
The energies can be conveniently split into mean and fluctuating
components, into poloidal and toroidal components and further into
equatorially-symmetric and -antisymmetric components, thus giving a
rather complete description of the scales of the convective flow and
the multipole structure of dynamos. The mean and fluctuating toroidal
and poloidal components of the kinetic energy are defined as
\begin{subequations}
\begin{eqnarray}
\label{7}
\bar E_p = \langle \big(\nabla \times ( \nabla \bar \pol \times \vec r )
\big)^2/(2\rhobar)  \rangle,
& 
\bar E_t = \langle \big(\nabla r \bar\tor \times \vec r \big)^2
/(2\rhobar)  \rangle, \\
\label{8}
\check E_p =  \langle \big(\nabla \times ( \nabla \check \pol
\times \vec r)  \big)^2/(2\rhobar) \rangle, &
\check E_t = \langle \big( \nabla r \check \tor \times\vec r \big)^2
/(2\rhobar) \rangle,
\end{eqnarray}
\end{subequations}
where  angular brackets $\langle \rangle$ denote averages over the
spherical volume of the shell. Magnetic energy components are defined
analogously with $\polB$ and $\torB$ replacing $\pol$ and $\tor$ and
without the factor $\rhobar^{-1}$ within the angular brackets. The
total energies are, of course, the sum of all components. The Nusselt
number is defined as the ratio between the values of the luminosity of
the convective state and of the basic conduction state, 
$$
Nu=-\frac{\exp(\Nrho)-1}{4\pi n c_1 \rhobar(r_i)^n}\int_{\partial V}
\kappa\rhobar\Tbar (\partial_r S) \r^2 \sin\theta d\theta d\varphi,
$$
with the integral taken over the top surface $\partial V$. Apart from
quantifying the 
heat transport of convection, the value of the Nusselt serves as a
convenient proxy for the supercriticallity of the solution.

Other diagnostic quantities that are sometimes used to quantify
convective and dynamo solutions can be derived from these
quantities. For example, a non-dimensional magnetic Reynolds number,
Rossby number and Lorentz number are given by
$$
Rm= \Pm \sqrt{2 E_\text{kin}}, \quad 
Ro= \frac{2}{\tau}\sqrt{2 E_\text{kin}}, \quad
Lo= \frac{2}{\tau}\sqrt{2 E_\text{mag}},
$$
respectively.

\section{Benchmarking and validation}
\label{sec.3}
To perform the numerical simulations of this study, we have extended
our mature Boussinesq code
\citep{TilgnerBusse1997,Bu03,Simitev2005a,Simitev2006b,Simitev2008a,Simitev2009,Simitev2012a}
to solve the anelastic problem described in section \ref{sec.2}.
Despite similarities with the Boussinesq code, this is a major
modification both in terms of the mathematical model and the numerical
code. In order to validate 
the new code, we wish to report here a comparison with the anelastic
dynamo benchmarks recently published by \cite{JonesBench}. To
aid comparison with the latter paper in this section only we employ
the alternative nondimensionalisation used in \citep{JonesBench} where
the magnetic diffusion timescale rather that the viscous diffusion
scale is adopted. Our output results from the three benchmark cases
defined in \citep{JonesBench} are summarized in table \ref{tab01}, and
selected components of the solution are plotted in figure \ref{fig01}. We
achieve near exact agreement with the results reported in
\citep{JonesBench} for the hydrodynamic case  and the steady dynamo
case, labeled \textit{Benchmark 1} and \textit{Benchmark 2} in table
\ref{tab01}. Our results for the unsteady dynamo case labeled
\textit{Benchmark 3} in table \ref{tab01} show certain
insignificant differences from the values reported in
\citep{JonesBench}. The mean zonal flow we obtain in this case is 33\%
larger than that reported in the benchmark paper even though the
angular momentum of our run remains smaller than $10^{-3}$. This
somewhat larger differential rotation gives rise to a magnetic field
with energy 6\% larger than that reported in the benchmark paper. The
reasons for the discrepancy are that in this case we have imposed a two-fold
azimuthal symmetry and used markedly lower resolution that the codes
participating in \citep{JonesBench}. This has been done to reduce
computing time in this challenging case.

\section{Similarities of anelastic and Boussinesq spherical convection and dynamos}
\label{sec.4}

\subsection*{Some linear results}



These results are not intended to be comprehensive. Since the
anelastic and the Boussinesq models are essentially different it is no
obvious how to compare results directly. We will introduce a
``supercriticallity'' criterion -- two cases with the same values of
all basic parameters \eqref{numbers} apart from the Rayleigh number
$R$ are directly comparable if they have equal values of the ratio
$R/R_c$, where $R_c$ is the corresponding critical Rayleigh number for
the onset of convection.

For this purpose, here we present $R_c$ curves of typical cases we explore.

\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/01Calcs/2013.07.16.linear/Nr_xmin.v1.agr
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/01Calcs/2013.07.16.linear/Nr_xmin.v1.fig
\psfrag{Nr}{$N_\rho$}
\psfrag{x}{$x$}
\psfrag{inner boundary}{inner boundary}
\psfrag{outer boundary}{outer boundary}
\psfrag{0}{0}
\psfrag{2}{2}
\psfrag{4}{4}
\psfrag{6}{6}
\psfrag{2.0}{2.0}
\psfrag{2.5}{2.5}
\psfrag{3.0}{3.0}
\epsfig{file=../01Calcs/2013.07.16.linear/xmin/Nr_xmin.v1.eps,width=\columnwidth,clip=}
\end{center}
\vspace*{-6mm}
\caption{
}
\label{f00101}
\end{figure}

The main difference with Boussinesq convection is that onset of
convection moves closer to the outer boundary - make new plot of
position. Illustrate with case: 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/a.mv.e065p1t2r2500000m1p4N4/pol.48.pdf}{a.mv.e065p1t2r2500000m1p4N4}.

\subsection*{Relaxation oscillations in anelastic convection}

It is interesting that relaxation oscillations that have been found in
Boussinesq simulations also appear in anelastic convection as seen in
figures \ref{f001} and \ref{f002} for the case 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r20000000m4n2N3}{ae035p1t40r20000000m4n2N3}.

Another, case is 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r15000000m4n2N3}{ae035p1t40r15000000m4n2N3}. 

\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p1t40r20000000m4n2N3/version2/3mmm.tex
\epsfig{file=../figs/ae035p1t40r20000000m4n2N3/version2/3mmm.epsi,width=\columnwidth,clip=}
\end{center}
\caption{(\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r20000000m4n2N3}{ae035p1t40r20000000m4n2N3}) 
  A period of relaxation oscillations of convection in the case
  $\eta=0.35$, $P=1$, $\tau=4\times10^4$, $R=2\times10^7$, $n=2$,
  $N_\rho=3$, and four-fold azimuthal symmetry. The left halves of each
  plot show lines of constant $\overline{u}_\varphi$ in a meridional
  cut and the right halves streamlines, $r\partial_\varphi v = $
  constant, in the equatorial plane. The separation in time is 0.02.}
\label{f001}
%\end{figure}
%\begin{figure}
\psfrag{E}{$E$}
\psfrag{t}{$t$}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p1t40r20000000m4n2N3/Engs.agr
\epsfig{file=../figs/ae035p1t40r20000000m4n2N3/Engs.eps,width=\columnwidth,clip=}
\end{center}
\caption{(\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r20000000m4n2N3}{ae035p1t40r20000000m4n2N3}) 
 Time series of energy densities
  of convection in the case plotted in \ref{f001}. Red, blue and green
  lines indicate $\overline{E}_t$, $\overline{E}_t$, $\widetilde{E}_p$,
  respectively.  $\overline{E}_p$ smaller by more than an order of
  magnitude than the other energy densities and has not been plotted
  for this reason. The three dotted vertical lines indicate the time
  instances where the snapshots in figure \ref{f001} are taken.
}
\label{f002} 
\end{figure}


In the presence of magnetic field relaxation oscillations
disappear. However, the dynamic balance between diff rotation and
magnetic fields produces oscillating dynamos. Here is a typical
oscillation hemispherical dynamo, as shown in figures \ref{f003} and
\ref{f004}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r10000000m4p2n2N3}{(ae035p1t40r10000000m4p2n2N3)}. 

This case is also interesting because 
the mean poloidal magnetic energy is dominant (both dipole and
quadrupole) so it is a MD dynamo but the dynamo is oscillating as if
it was and FD dynamo. The oscillations of the magnetic field are very
regular with similar period twice as large as the period of the
relaxation oscillations that have died out.
\begin{figure}
\begin{center}
\psfrag{Ek}{$E$}
\psfrag{Md,q}{$M_\text{dip,quad}$}
\psfrag{Mquad}{$M_\text{quad}$}
\psfrag{coeff}{$H^0_{1,2}$}
\psfrag{t}{$t$}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p1t40r10000000m4p2n2N3/Engs.v2..agr
\epsfig{file=../figs/ae035p1t40r10000000m4p2n2N3/Engs.v2.eps,width=\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r10000000m4p2n2N3/}{(ae035p1t40r10000000m4p2n2N3)} 
Time series of energy densities of a hemispherical dynamo in the case
$\eta=0.35$, $P=1$, $\tau=4\times10^4$, $R=10^7$, $P_m=2$, $n=2$,
$N_\rho=3$, and four-fold azimuthal symmetry.  
The upper panel shows kinetic energy densities of the velocity field. 
The second panel shows energy densities of dipolar components (thin
solid lines) and quadrupolar components (thick brocken lines) of the
magnetic field. The mean toroidal components are represented by red
lines, the fluctuating toroidal -- by blue lines, the mean poloidal --
by black lines and the fluctuating poloidal -- by green lines. The
lowest panel shows the coefficient of the axial dipole component
$H_1^0$ (red line) and the coefficient of the axial quadrupole
component (blue line) at $r=r_i+0.5$. 
}
\label{f003} 
%\end{figure}
%\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p1t40r10000000m4p2n2N3/3mmm.tex
\epsfig{file=../figs/ae035p1t40r10000000m4p2n2N3/3mmm.eps,width=\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r10000000m4p2n2N3/}{(ae035p1t40r10000000m4p2n2N3)}
A period of hemispherical dynamo oscillations in the case shown in
figure \ref{f003}. The left half of each plot shows lines of constant
$\overline{B}_\varphi$ and the right half shows meridional field lines
$r\sin \theta \partial_\theta \overline{h} = $ const. Plots are
arranged row-by-row and the separation in time is 0.08.}
\label{f004}
\end{figure}

Increasing Rayleigh number we get slightly more dipolar oscillations
as in the case 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r15000000m4p2n2N3/}{(ae035p1t40r15000000m4p2n2N3)}.

This case is not
plotted as it is not significantly different from the previous case
shown in \ref{f003} and \ref{f004}.

Increasing both the Rayleigh number and the magnetic Prandtl number
leads to much more dipolar oscillating case with APERIODIC reversals

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r15000000m4p4n2N3}{(ae035p1t40r15000000m4p4n2N3)}

as shown in figures \ref{f005} and \ref{f006}. Aperiodic reversals are  interesting as a model of the
geodynamo polarity reversals. 
\begin{figure}
\begin{center}
\psfrag{H12}{$H^0_{1,2}$}
\psfrag{t}{$t$}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p1t40r15000000m4p4n2N3/coeff.agr
\epsfig{file=../figs/ae035p1t40r15000000m4p4n2N3/coeff.eps,width=\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r15000000m4p4n2N3/}{(ae035p1t40r15000000m4p4n2N3)} 
Aperiodicaly reversing dynamo. Shown are time series of coefficient
of the axial dipole component $H_1^0$ (red line) and the coefficient
of the axial quadrupole component (blue line) at $r=r_i+0.5$ in the case 
$\eta=0.35$, $P=1$, $\tau=4\times10^4$, $R=1.5\times10^7$, $P_m=4$, $n=2$,
$N_\rho=3$, and four-fold azimuthal symmetry. One reversal happens
within the shaded area when the dipole component changes sign.
}
\label{f005} 
%\end{figure}
%\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p1t40r15000000m4p4n2N3/3mmm.tex
\epsfig{file=../figs/ae035p1t40r15000000m4p4n2N3/3mmm.eps,width=\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r15000000m4p4n2N3/}{(ae035p1t40r15000000m4p4n2N3)}
An apperiodic reversal of the dynamo shown in figure \ref{f005}.  The
left half of each plot shows lines of constant $\overline{B}_\varphi$
and the right half shows meridional field lines $r\sin \theta
\partial_\theta \overline{h} = $ const. Plots are arranged row-by-row
and the separation in time is 0.04 falling in the shaded area of
figure \ref{f005}.} 
\label{f006}
\end{figure}

\subsection*{Dipolar Oscillations and benchmarks}

The benchmark case is oscillatory and shows nice dipolar oscillations
as illustrated in figure \ref{f007}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ben3.ae035p2t40r25000000m2p2n2N3/}{(ben3.ae035p2t40r25000000m2p2n2N3)}

\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ben3.ae035p2t40r25000000m2p2n2N3/3mmm.tex
\epsfig{file=../figs/ben3.ae035p2t40r25000000m2p2n2N3/3mmm.eps,width=\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ben3.ae035p2t40r25000000m2p2n2N3/}{(ben3.ae035p2t40r25000000m2p2n2N3)}
A period of dipolar dynamo oscillations in the benchmark case
$\eta=0.35$, $P=2$, $\tau=4\times10^4$, $R=2.5\times10^7$, $P_m=2$, $n=2$,
$N_\rho=3$, and two-fold azimuthal symmetry.  The
left half of each plot shows lines of constant $\overline{B}_\varphi$
and the right half shows meridional field lines $r\sin \theta
\partial_\theta \overline{h} = $ const. Plots are arranged row-by-row
and the separation in time is 0.015.}
\label{f007}
\end{figure}

Dipolar oscillations in thin shells

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/cases/a.mv.e065p1t2r1200000m1p5N2}{a.mv.e065p1t2r1200000m1p5N2}
--
Plot a sequence near t=85.



\subsection*{Magnetic field produces Solar-like differential rotation}

A sufficiently strong magnetic field may counterbalance the
geostrophic constraint and produce conical profile of the differential
rotation with contours that are aligned in radial direction as shown
in figure \ref{f008}. This is reminiscent of the profile of the Solar
differential rotation measured by helioseismology e.g. \cite{Schou1998}.

\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ben3.ae035p2t40r25000000m2p2n2N3/12.um.82.ps
\epsfig{file=../figs/ben3.ae035p2t40r25000000m2p2n2N3/12.um.82.epsi,width=0.5\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ben3.ae035p2t40r25000000m2p2n2N3/}{(ben3.ae035p2t40r25000000m2p2n2N3)}
Solar-like differential rotation $\overline{u}_\varphi=$ const.~in the left
halves of the plot and azymuthally-averaged meridional circulation in
the right halves in the case shown in figure \ref{f007}. The plots are
time averaged over a period of 0.005.}
\label{f008}
\end{figure}


Comment that this is why Brandenburg and Jorg observe solar like diff
rot in their coronal mass ejections - Jorg's thesis.
\url{http://arxiv.org/pdf/1301.2248.pdf}


\subsection*{Banded differential rotation like Jupiter}

See figure \ref{f009}.

\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/ae035p2t40r15000000m4p2n2Nrho3a/12.um.11.ps
\epsfig{file=../figs/ae035p2t40r15000000m4p2n2Nrho3a/12.um.11.epsi,width=0.5\columnwidth,clip=}
\end{center}
\caption{\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p2t40r15000000m4p2n2Nrho3a/}{(ae035p2t40r15000000m4p2n2Nrho3a)}
Jupiter-like banded differential rotation $\overline{u}_\varphi=$
const.~in the left halves of the plot and azymuthally-averaged
meridional circulation in the right halves in the case $\eta=0.35$,
$P=2$, $\tau=4\times10^4$, $R=1.5\times10^7$, $P_m=2$, $n=2$,
$N_\rho=3$, and four-fold azimuthal symmetry. } 
\label{f009}
\end{figure}

 -\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/ae035p1t40r10000000m4p2n2N3/}{ae035p1t40r10000000m4p2n2N3}/um.59.pdf

 -Magnetic field is important
- see the effect of hemispherical dynamo - this case is  good to use
 as figs of magn  field already included \ref{f003} \& \ref{f004}. --- 
Perhaps no evidence for this.

At stronger magnetic fields band structure does not persist and
radially-constant conical Solar profiles of diff rot emerge as in
figure \ref{f007}.


\subsection*{Subsurface decrease of differential rotation}

This is very interesting, see cases 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/a.e065p1t2r2000000m1p6N3.mv.vc}{a.e065p1t2r2000000m1p6N3.mv.vc}
(time 18 to 20)

and 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/a.e065p1t2r1800000m1p8N3.mv.vc/}{(a.e065p1t2r1800000m1p8N3.mv.vc)} and 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/a.e065p1t2r1800000m1p8N3.mv.vc.per/}{(a.e065p1t2r1800000m1p8N3.mv.vc.per)}.

 --- Do these oscillate?

 --- Can we find equatorward propagation?

Oscillating cases:

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/a.mv.e065p1t2r2000000m1p4N3}{a.mv.e065p1t2r2000000m1p4N3}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/a.e065p1t2r2000000m1p6N3.mv.vc}{a.e065p1t2r2000000m1p6N3.mv.vc}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/cases/a.mv.e065p1t2r2000000m1p4N3}{a.mv.e065p1t2r2000000m1p6N3}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/Fri.2012.11.16/um.22.pdf}{a.mv.e065p1t2r2000000m1p4N3-sequence}

a.mv.e065p1t2r2500000m1p4N3

\newpage

\begin{figure*}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/figs/a.mv.e065p1t2r2500000m1p4N3/fig01/3mmm.tex
\epsfig{file=../figs/a.mv.e065p1t2r2500000m1p4N3/fig02/7.epsi,width=\textwidth,clip=}
\end{center}
\caption{(a.mv.e065p1t2r2500000m1p4N3)
One period of dipolar oscillations.
}
\label{f0401}
\end{figure*}
\begin{figure}
\psfrag{t}{$t$}
\psfrag{H012}{$H_0^{1,2}$}
\psfrag{G012}{$G_0^{1,2}$}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/figs/a.mv.e065p1t2r2500000m1p4N3/fig01/coeffs.agr
\epsfig{file=../figs/a.mv.e065p1t2r2500000m1p4N3/fig01/coeffs.eps,width=\columnwidth,clip=}
\end{center}
\caption{(a.mv.e065p1t2r2500000m1p4N3)
Dipolar oscillations shown in the timeseries of the the axisymmetric
toroidal coefficients $G_0^{1}$ (red) and $G_0^{2}$ 
(blue) and the axisymmetric poloidal coefficients $H_0^{1}$ (red) and
$H_0^{2}$ (blue) describing the main dipolar and poloidal
contributions in the spherical harmonics expansion of the magnetic
field in the case $\eta=0.65$, $P=1$,
$\tau=2\times10^3$, $R=2.5\times10^6$, $P_m=4$, $n=2$, $N_\rho=3$
shown in figure \ref{f0401}. 
}
\label{f0402}
\end{figure}




\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/2013-06-13/../figs/a.e065p1t2r1800000m1p8N3.mv.vc/12.um.18.ps
\epsfig{file=../figs/a.e065p1t2r1800000m1p8N3.mv.vc/12.um.18.epsi,width=0.5\columnwidth,clip=}
\end{center}
\caption{(a.e065p1t2r1800000m1p8N3.mv.vc)
Solar like (sub-surface decrease) differential rotation $\overline{u}_\varphi=$
const.~in the left halves of the plot and azymuthally-averaged
meridional circulation in the right halves in the case $\eta=0.65$,
$P=1$, $\tau=2\times10^3$, $R=1.8\times10^6$, $P_m=8$, $n=2$,
$N_\rho=3$. } 
\label{f010}
\end{figure}


\subsection*{High Prandtl number convection looks very solar like}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/Fri.2012.11.16/a.e075p5t05r800000m1N3.sf.vc.Di0.plt.22.pdf}{a.e075p5t05r800000m1N3.sf.vc.Di0.plt.22.pdf}

Note the double-cell meridional circulation and the tendency towards
conical profiles of diff rotation.

\subsection*{Variable viscosity}
Thin downwelling and thick upwelling 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/2012-12-10/a.e065p04t04r300000m1p2.sf.vr4.plt.13.pdf}{a.e065p04t04r300000m1p2.sf.vr4}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/2012-12-10/a.e065p06t06r600000m1p2.sf.vr4.plt.11.pdf}{a.e065p06t06r600000m1p2.sf.vr4.plt.11.pdf}

Compare cases with otherwise fixed parameters:

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/2013.07.19/a.e075p5t05r800000m1N3.sf.vr.plt.28.pdf}{a.e075p5t05r800000m1N3.sf.vr.plt.28.pdf}

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/2013.07.19/a.e075p5t05r800000m1N3.sf.vc.plt.16.pdf}{a.e075p5t05r800000m1N3.sf.vc.plt.16.pdf}


Cases with fixed params:

a.e065p05t03r300000m1N3.sf.vc,

a.e065p05t03r300000m1N3.sf.vr,

a.e065p05t03r300000m1N3.sf.vrkr

a.e065p1t1r800000m1N3.sf.vc,

a.e065p1t1r800000m1N3.sf.vr,

a.e065p1t1r800000m1N3.sf.vrkr


\section{Buoyancy-dominated convection}
\label{sec.5}


\begin{figure}
\begin{center}
% laptop /home/rds9k/ownCloud/046_RemarksAnelastic/engs_seqs/cases.engs03.agr
\psfrag{up}{\hspace{-6mm}$\overline{u}_\varphi(r_o,\pi/2)$}
\psfrag{R-6}{\hspace{-5mm}$R\times 10^{-6}$}
\psfrag{R-5}{\hspace{-5mm}$R\times 10^{-5}$}
\epsfig{file=../engs_seqs/cases.engs03.eps,width=\columnwidth,clip=}
\end{center}
\vspace*{-6mm}
\caption{
(a,b) Equatorial zonal flow at the equator $\overline{u}_\varphi(r_o,\pi/2)$.
(c,d) Average kinetic energy densities and (e,f) average magnetic energy
densities as functions of the Rayleigh number $R$ in the cases
(a,c,e) $\eta=0.65$, $P=1$, $\tau=2000$, $n=2$, $N_\rho=3$, and
(b,d,f), $\eta=0.65$, $P=0.5$, $\tau=300$, $n=2$, $N_\rho=3$
Nonmagnetic convection cases are denoted by thin symbols in
(a,b,c,d). Dynamo cases are denoted by thick symbols in all panels and
have $P_m=2$ in (a,c,e) and $P_m=6$ in (b,d,f). Black circles, red
squares, green pluses and blue crosses denote $\overline{X}_p$,
$\overline{X}_t$, $\tilde{X}_p$, $\tilde{X}_t$, with $X=E,
M$. Vertical dash-dotted lines denote the transition to buoyancy-dominated regime.
}
\label{f00101}
\end{figure}


\paragraph{Discussion of case (a.e065p1t2r10000000m1p2N3n2.sf.vc) shown in
figures \ref{f040} and \ref{f050}}

The influence of the magnetic field is quite remarkable. Consider the
following observations illustrated in figures \ref{f040} and
\ref{f050}:
\begin{itemize}
\item pure convection - diff rot is monotonously increasing/decreasing
  towards the outer surface in the prograde/retrograde case.
\item pure convection in the buoyancy-dominated regime is retrograde.
\item pure convection - does not show columnar structure at depth.
\item dynamo - the energy of the self-sustained magnetic field is an
  order of magnitude smaller that the kinetic energy.
\item dynamo - diff rotation is reversed from retrograde to prograde
  by even by such relatively  weak magnetic field.
\item dynamo - differential rotation has a maximum inside the shell
and then shows a subsurface decrease.
\item dynamo - convective columns are visible in depth.
\item dynamo - has a strong dipolar symmetry because of the strong
  polar convection. Unfortunately this dipole does not oscillate.
\end{itemize}

Question - Is there a transition from multipolar to dipolar dynamos
with the increase of $R$? There must be - in the prograde regime at
low $R$ thin shell dynamos are mixed multipolar; in the retrograde
regime at high $R$ dipolar dynamos are found as figure \ref{f050}
demonstrates. There must be a transition in between - look at newly
started sequences.



\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/figs/a.e065p1t2r10000000m1p2N3n2.sf.vc
\epsfig{file=../figs/a.e065p1t2r10000000m1p2N3n2.sf.vc/f01/3mmm.epsi,width=\columnwidth,clip=}
\end{center}
\caption{(a.e065p1t2r10000000m1p2N3n2.sf.vc and
  a.e065p1t2r10000000m1N3n2.sf.vc)
Left column: dynamo case a.e065p1t2r10000000m1p2N3n2.sf.vc; Right
column: convection only case a.e065p1t2r10000000m1N3n2.sf.vc.
Row one: $u_r$ at $r=0.95+r_i$; Row two: $u_r$ at $r=0.3+r_i$; Row 3:
pol streamlines in the equatorial plane; Row 4: $\overline{u}_\varphi$
and mer streamlines average over a period of 0.025.
}
\label{f040}
\end{figure}

\begin{figure}
\begin{center}
% laptop:/home/rds9k/ownCloud/046_RemarksAnelastic/figs/a.e065p1t2r10000000m1p2N3n2.sf.vc
\epsfig{file=../figs/a.e065p1t2r10000000m1p2N3n2.sf.vc/f01/5mag.epsi,width=\columnwidth,clip=}
\end{center}
\caption{(a.e065p1t2r10000000m1p2N3n2.sf.vc)
Magnetic components of the same case as in the left column of figure
\ref{f040}.
Left plot:  $B_r$ at $r=1.13+r_i$; Right plot:: $\overline{B}_\varphi$
and mer fieldlines average over a period of 0.025.
}
\label{f050}
\end{figure}







\subsection{Interesting Boussinesq cases}

a.mv.e065p1.2t2r500000m1p4.5N001  - maybe we should start cases with
slowly increasing $N_\rho$ to follow transitions. What will
quadrupolar dynamos become?




\subsection*{Effects of radial variations of viscosity and diffusivity}

\subsection*{Evidence of localized convection or spotty convection} 




\begin{figure*}
\epsfig{file=../plots/a.e065p03t02r...m1p10N3.sf.vr/3mmm.epsi,height=0.7\textheight}
\caption{A figure like this but from the sequences in figure \ref{f00101}.}
\end{figure*}


\bi
\item I have also thinking about your paper. It could be interesting
  to plot wavenumber m (number of equatorial cells) as function of tau
  and Pr, 0.1<Pr<2 for convenient value of r and N3 because there is
  this   remarkable transition from low m to high m as function of
  mainly Pr. I am   wondering whether this transition may be hysteretic. Yours,F.

Thank you for the suggestion. Few questions, that will help a lot:
 - Why do you think this is remarkablle? In Boussinesq linear
analysis it is not uncommon to see a jump in wavenumber, even a  
large jumps say fig 1 from our paper in GAFD-2005.

 Here we have a much more dramatic phenomenon. I GAFD m changed from  2 
 at Pr=0.2 to 3 for Pr=2 and similar changes with higher m at at
  higher tau. In the a-case it changes from m=7 at Pr=0.3 to m=20 at Pr=1.

- How will we select comparable cases when we vary Pr? Clearly
      this
will be only a rough comparison, but we need to achive approx the   
 same supercriticallity as Pr is varied. This question is actually
worrying me.

 At fixed t the critical R should not vary that much. It seems to me
 that there could be substantial hysteresis.

 It would be nice to have a liner code so we can find critical R.

 Apparently Wicht used one supplied by Chris Jones. It will be
 important to use a case of nu(r) which has a reasonable deviation of
s
 from the S-reference throughout the layer such the influence of
 convection filling only part of the layer is minimized.


\item Neat separation between convection in the tangent cylinder and
  convection in the polar regions
\item Asymmetry of descending and ascending flows - this is visible in most
cases at small $\Pr$ and $\tau$.  - include a characteristic plot,
say, 

\href{http://www.maths.gla.ac.uk/~rs/res/B/Anel/2012-11-27/a.e065p05t03r300000m1N3.sf.vr.plt.26.pdf}{a.e065p05t03r300000m1N3.sf.vr}

and similar.
\item Transition to retrograde differential rotation. - Many cases,
  most accessible for low values of $\tau$. Can we reach beyond this
  regime?  What are dynamos like in the regime of retrograde diff
  rotation. In this regime diff rotation is substantially decreasing
  towards the surface, so according to Yoshimura dynamo waves should
  travel towards the equator.
\item Spread of convection throughout the layer - in models with variable viscosity (when $\nu=\rhobar^{-1}$ is
  constant DYNAMIC viscosity).
\item In models with variable viscosity (when $\nu=\rhobar^{-1}$ is
  constant DYNAMIC viscosity) convection spreads throughout the layer
  and not only near the surface as in the simulations of Whicht and of
  Jones. So each time when they claim that a particular mechanism is
  due to this effect they may not be right.
\item Double-cell meridional circulation as per Junwei Zhao, poleward
  flow near surface.

\ei



\section{Conclusion}
\label{sec.6}

\paragraph{Summary of results}

\begin{itemize}
\item
We have shown that even a rather small increase of the value of the
Prandtl number above unity can produce profiles of the differential
rotation that deviate significantly from geostrophy even in rather
laminar solutions. When such differential rotation profiles are
further affected by self-sustained magnetic field differential
rotation not unlike the one inferred from helioseismological
observations can be found.
\end{itemize}

\paragraph{Discussion and implications}

\begin{itemize}
\item The results may be relevant as a model of  Solar magnetic field
  oscillations, structure of the Solar convection zone,  and
  differential rotation. 
\begin{itemize}
\item
  \textbf{Solar cycle --} The direction of dynamo waves in
  self-consistent simulations with prograde differential rotation is
  always from equator to   poles in contrast to observations. We wish
  to check whether retrograde   differential rotation generated in the
  buoyancy-dominated regime may   reverse the direction of dynamo
  waves so that they   propagate from   poles to equator in agreement
  with observations.  
\item 
\textbf{Structure of convection zone --} It will be of interest to find whether different regimes of
  convection may develop in the inner and at the outer part of the
  shell simultaneously  so that organised geostrophic convection is
  hidden below a near-surface layer of well-mixed (``turbulent'')
  convection.   We believe that such a configuration is similar to the
  structure of the Solar convection zone. Indeed, it is well established by observation
  of the surface  velocities and magnetic fields that buoyancy effects
  dominate near   the surface. As a result the near-surface flows and
  magnetic fields   are strongly turbulent and lack any global
  organization. On the   other hand, the regularity of the 11-year
  Solar cycle   is also well   established and the regular reversals
  cannot be   achieved without   large-scale organized field and
  flow. We would   like to argue that   the organized columnar
  structures that are   characteristic in the   rotation-dominated
  regime will play this role.  

  Variation in density is a necessary ingredient in this scenario.
  Indeed, in the Boussinesq  approximation, the rotation parameter
  $\tau$ and the critical  Rayleigh number $\R_c$ are constant
  throughout the shell. This means  that the entire volume of the shell will be either in the
  rotation-dominated or the buoyancy-dominated regime. In contrast,
  the anelastic approximation allows for radial variation in density,
  viscosity and entropy diffusivity so that $\tau$ and $\R_c$ vary
  with radius. In this way, one may hope to achieve a configuration
  where the inner part of the shell is in the rotation-dominated
  regime and the outer part is in the buoyancy-dominated regime. 
\item 
\textbf{Solar differential rotation --}
Finally, differential rotation may also assume more Solar-like
profile. 
\end{itemize}
\item In fact, variation in density only may not be sufficient to
produce Solar-like differential rotation as indicated by Gastine et
al. In this case convection flows in the rotation-dominated region are
relatively weaker by comparison with flows in the buoyancy-dominated
near-surface layer. A stronger organized flow is likely needed to
produce the solar-lice cyclic oscillations and differential rotation
profile. 
\end{itemize}


\paragraph{Future directions}

-- Include radial variation of viscosity and entropy diffusivity. 

-- Go beyond the buoyancy-dominated regime into regime III.

-- Explore parameter regime so as to find periodic reversals. Will
they propagate equator-wise?



\section*{Acknowledgment}
This research has been supported by NASA Grant NNX-09AJ85G.
RDS and FHB acknowledge the hospitality of Stanford University and
UCLA. RDS enjoyed a period of study leave granted by the University of 
Glasgow.

\bibliographystyle{plain}
\bibliography{biblio}

\end{document}
