\section{Introduction}
\label{sec:intro}

Rotating convection   occurs in the interiors of  many stars and planets. 
But there is no adequate theory  for it yet,
despite many decades of research. 
In order to determine the  structure of a non-rotating star or planet, one typically employs 
 mixing length theory 
 (e.g.~\citealt{BohmVitense1958}), which
 despite its crudeness, 
 accounts for the main structural feature of a convection zone:  a nearly uniform entropy.
 But it remains unclear whether mixing length theory applies quantitatively, and how it may be extended to treat more subtle effects such as rotation and overshooting.

Convection in astrophysical bodies---rotating or not---is  difficult theoretically  because the
flow is
 turbulent.  Compounding the difficulty,   microscopic 
 viscosities 
($\nu$) and thermal diffusivities ($\kappa$)
   are 
 typically extremely small in stars and planets, orders of magnitude smaller than  
 accessible by experiment or simulation \citep[e.g.,][]{Spiegel1971,Miesch2005}. 
  It is generally believed that  bulk properties
 of turbulence should be independent of microscopic diffusivities in the limit
 that these are extremely small.  But it remains unclear 
 if that belief is correct, and if it is,  whether
  experiments and simulations are adequately probing that limit.

  Rather than study the full problem in a star or planet, a common approach is to consider
 a simpler setup, Rayleigh-Benard convection (RBC), 
which has been extensively studied theoretically \citep[e.g.,][]{C1961,GrossmanLohse2000}, numerically, and experimentally (e.g.~\citealt{Ahlers2009}, and references therein). 
 In RBC, 
a  fluid layer is sandwiched
 between two horizontal plates. Convection is driven in the interior
 by holding the bottom plate  hotter  than the top.
 The fluid obeys the
 Boussinesq equations,  
  in which the density is constant, 
  the velocity field is incompressible, and temperature fluctuations give rise to vertical buoyancy
 forces \citep[e.g.][]{C1961}.  This is perhaps the simplest
 set of equations that self-consistently evolve turbulent convection.  

Because of its simplicity,  RBC is an ideal testbed for mixing length theory.  
The typical goal is to
 predict the heat flux through
the fluid given the temperature drop between the plates.  
But the mixing length prediction,  taken at face value, 
 wildly disagrees with experimental and numerical results.\footnote{In the absence of rotation, mixing length theory predicts 
 Nu $ = $ const$\times $Ra$^{1/2}$ \citep{Kraichnan1962}, where Ra (the Rayleigh number) quantifies the temperature
 drop and Nu (the Nusselt number) quantifies the  heat flux. 
  By contrast, experiments
 and simulations obtain an exponent 
 $\approx 1/3$ at high Ra, i.e., at small $\nu$ and $\kappa$
 \citep{Spiegel1971,ShraimanSiggia1990,Ahlers2009,2010NJPh...12k5002N}.
 The discrepancy increases as $\nu,\kappa\rightarrow 0$.
Note that the mixing length prediction  may 
be derived by dimensional analysis under the assumption that
 viscosity and thermal diffusivity play no
 role (\S \ref{mlt}).}
The reason for the discrepancy is that a naive application of mixing length theory
assumes  a constant temperature gradient throughout the fluid.  But
in reality,
most of the temperature drop 
occurs in extremely thin boundary layers near the top and bottom plates,
where fluid  velocities are nearly zero and heat is transported primarily by
conduction. 
In fact, 
the relationship between temperature 
drop and heat flux can be accounted for by considering only the behavior of
 the boundary layers \citep{Malkus1954}. 
For {\it rotating} RBC, the situation is  similar:  \cite{King2012}  
 extend the boundary layer analysis of \cite{Malkus1954} 
to include rotation, and thereby obtain
impressive agreement with experiments and simulations (see also \citealt{King2009,King2013a}). 

However, for the purpose of explaining convection in astrophysical bodies, 
it is the behavior of the bulk of the fluid that is primarily of interest---not the boundary 
layers. The boundary layers in RBC are very different from the boundary 
of a  convection zone in a star or planet. But one might expect that the turbulent dynamics
of the interior fluid will be similar in the two cases.  Therefore in this paper we focus on the
dynamics of the
 interior
fluid in an RBC-like system.  We shall show that not only do the properties of the interior fluid converge
in the limit of small diffusivities, but  they converge to the prediction
of the rotating mixing length theory  first proposed by \cite{Stevenson1979}.
 
The organization of this paper is as follows. We first set up the problem (\S \ref{Theory}), and then derive the predictions of
 mixing length theory (\S \ref{mlt}).  Next, we test the theory in detail with a suite of numerical
 simulations that are similar to---but slightly different from---standard rotating RBC  (\S \ref{Results}).  We also run some comparison simulations with standard
 RBC  (\S \ref{RBC}). We conclude with a summary and discussion (\S \ref{Conclusions}).
 
\section{Problem Setup}
\label{Theory}

We consider a convecting fluid  in a region of depth $H$ that carries thermal flux $F$ and rotates at rate $\Omega$.
Our goal is to determine 
  the    properties of the turbulently convecting fluid, 
 focusing in particular on the mean temperature gradient, the rms velocity
and temperature fluctuations, and the characteristic size  of the eddies that dominate the flux. 
We eschew the commonly used dimensionless numbers (Ra and Nu) because
  they depend on the temperature
drop across a  region that includes both convecting fluid and boundary layers, and
are thus a better diagnostic of the boundary layers than of the convective fluid.
(But see Table \ref{Table}, where we catalog these numbers in our simulations to 
facilitate comparison with prior work).

 Our setup is perhaps the simplest one possible:
 Boussinesq convection in which 
the rotation vector is aligned with gravity.
 But the physical argument  will be
  straightforward to extend to more complicated situations.
The Boussinesq equations are 
 \begin{eqnarray}
\label{eq1aa}
&& \nabla \cdot \boldsymbol{v} = 0,
\label{eq1b}
\\
\label{eq1a}
&&\left(\partial_{t} + \boldsymbol{v}\cdot \nabla\right)\boldsymbol{v}  +2\Omega
\boldsymbol{e}_{z}
\times\boldsymbol{v}= -\nabla p +  T\boldsymbol{e}_{z} ,
 \\
&& \left(\partial_{t} + \boldsymbol{v}\cdot \nabla \right)T = 0,
\label{eq1c}
\end{eqnarray}
(e.g.~\citealt{Spiegel1960})\footnote{
If one decomposes the temperature into a horizontally 
uniform part and a remainder, 
$T=T_0(z)+T_1(\bld{x},t)$, then
Eq. \ref{eq1c}
becomes
$ \left(\partial_{t} + \boldsymbol{v}\cdot \nabla \right)T_1 = -v_z{dT_0/ dz} $
while Eq. \ref{eq1a} remains unchanged, aside from replacing $T\rightarrow T_1$ (since
 $T_0$ is absorbed by $p$).}
which 
 express, respectively,  incompressibility, force balance, and heat balance.
Here, $\bld{v}$ is the velocity field;
$p$ is a pressure variable that  is determined
by incompressibility; 
 the gravitational
acceleration points downwards (in the $-\bld{e}_z$ direction);
the rotation vector (with magnitude $\Omega$) points upwards; 
 and $T$ is a scaled temperature.  
More precisely, $T=g\alpha (\delta T_{\rm true}/T_{\rm true})$, where $g$ is gravity, 
 $\alpha$ is a thermal expansion coefficient, and $\delta T_{\rm true}/T_{\rm true}$ is the fractional difference in the true temperature relative to a reference value.
 One may also  think of $T$ as the buoyancy variable, 
 since in the Boussinesq approximation  temperature only affects the equation of motion through buoyancy; as such,  $T$ has dimensions of acceleration. 

Equations \ref{eq1a} and \ref{eq1c} must be supplemented with  diffusive terms
(e.g.,  
 $\nu\nabla^2\bld{v}$ and $\kappa\nabla^2T$, where $\nu$ is viscosity and $\kappa$ is thermal
 diffusivity).
In most astrophysical applications, $\nu$ and $\kappa$ are extremely small.
 Thus, while they can affect dynamics on extremely small lengthscales, 
 that presumably has little effect on large-scale properties.
 Our results below will confirm this presumption.  


\section{Mixing Length Theory of Rotating Convection}
 \label{mlt}
 
We assume (subject to verification below) that
the equations of motion without diffusivities are responsible for the large-scale dynamics
in the bulk of the convecting region.
We decompose the temperature  into a constant-gradient background and fluctuation
\begin{eqnarray}
T={\rm const}-N_*^2z+T_1 \ ,
\end{eqnarray}
where $N_*$ is a positive constant with dimensions of frequency. It is related
to the usual Brunt-Vaisala frequency $N$ via $N_*^2=-N^2$.
Linearizing
the Boussinesq equations  in $T_1$ and $\bld{v}$ yields the dispersion relation
for the growth rate $\gamma$:
\begin{equation}
\gamma^2= N_*^2{k_\perp^2\over k^2}-(2\Omega)^2{k_z^2\over k^2} \ ,
\label{eq:dr}
\end{equation}
where $k_\perp$ and $k_z$ are the components of the wavevector in the horizontal
and vertical directions.
Modes with $k_\perp> k_z(2\Omega/N_*)$ are exponentially unstable ($\gamma^2>0$), 
as they would be in a fluid without rotation. Although the Coriolis force stabilizes some modes, 
it cannot stabilize those whose phasefronts are sufficiently vertical.
But such modes obviously cannot continue to grow indefinitely.  Rather,  they 
 grow until nonlinearities cause  them to cascade. That is the basis for the mixing length
theory, described  in what follows. 

We focus here and in the remainder of this paper on the rapidly 
rotating limit, which we take to mean
\begin{equation}
2\Omega\gg N_* \ .
\end{equation}
The dispersion relation   implies that modes with $k_\perp/k_z> {\rm few \ }\times (2\Omega/N_*)\gg 1$ grow at the rate $\gamma\approx N_*$. 
Although all such modes are unstable, one expects that the modes
with the largest wavelengths are responsible for the bulk of the heat transport (or equivalently
that
   convection is dominated by the mode that transports the most heat -- e.g.~\citealt{Malkus1954,Stevenson1979}).
The largest wavelength in the vertical direction is the one that fits in a convection zone
of depth $H$:
\begin{equation}
k_z \sim {1\over H} \ .
\label{eq:kz}
\end{equation}
In the horizontal direction, the largest unstable mode then has
\begin{equation}
k_\perp \sim {1\over H}{2\Omega\over N_*} \ .
\label{eq:kperp}
\end{equation}
In  convective steady state the growth rate of such modes  equals their nonlinear
cascade rate, $\sim k_\perp v_\perp$,\footnote{
When employing the ``$\sim$'' symbol, we take  $v_\perp$, $v_z$, and $T_1$ as referring to R.M.S. values.} since that is the rate at which fluid is mixed across horizontal distance $1/k_\perp$.
We thus infer that
\begin{equation}
N_*\sim k_\perp v_\perp.
\label{eq:relation1}
\end{equation}
We wish to relate the temperature gradient to the thermal flux $F=\left<v_z T \right>$, where
angled brackets denote averaging over a horizontal plane; i.e.,
\begin{equation}
F \sim v_z T_1 \  ,
\end{equation}
since hot fluid typically rises and cold fluid sinks.
To complete the calculation requires an estimate for $T_1$.  As hot plumes rise, they carry the background temperature upwards for a time $1/\gamma\sim 1/N_*$ before cascading, 
i.e., 
$T_1\sim {d\left<  T\right>/ dz}\times {v_z/ N_*}$. 
 Since $d\left<  T\right>/ dz=-N_*^2$, we have
\begin{equation}
T_1
\sim v_zN_* \ .
\label{eq:relation3}
\end{equation}

Combining Equations (\ref{eq:kperp})--(\ref{eq:relation3}) and setting $v_\perp\sim v_z$\footnote{The linear dispersion relation yields $v_\perp\approx v_z(2\Omega/N)(k_z/k_\perp)$ for modes
with $\gamma\approx N_*$, and hence $v_\perp\sim v_z$ for wavenumbers given by Equations (\ref{eq:kz}) and (\ref{eq:kperp}).}
 yields the principal results:
the temperature gradient  
 is 
\begin{eqnarray}
\label{eq:ml1}
-{d\left<  T\right>\over dz}=N_*^2\sim {F^{2/5}\Omega^{4/5}\over H^{4/5}} \ ;
\end{eqnarray}
 the velocity and temperature fluctuations are
\begin{eqnarray}
v_z&\sim& {H^{1/5}F^{2/5}\over\Omega^{1/5}} \ ,
{\label{eq:ml2}}
 \\
T_1&\sim& {F^{3/5}\Omega^{1/5}\over H^{1/5}} \ ;
\label{eq:ml3}
\end{eqnarray}
and
the lengthscale of the modes that dominate heat transport is $1/k_\perp$ where
\begin{eqnarray}
k_\perp\sim {\Omega^{3/5}\over H^{3/5}F^{1/5}} \ .
\label{eq:ml4}
\end{eqnarray}

 Eqs. \ref{eq:ml1}--\ref{eq:ml4} are applicable 
for fast rotation ($\Omega\gg N_*$). 
For negligibly slow rotation,  one should
 replace $\Omega\rightarrow N_*$, yielding
   $N_*^2\sim F^{2/3}/H^{4/3}$, $v_z\sim (FH)^{1/3}$, and $k_\perp\sim 1/H$.
  These relations agree with 
the predictions of standard (non-rotating) mixing length theory. 
One may alternatively ``derive'' them by dimensional analysis.  However, 
they have proven more difficult to verify in simulations than their rotating counterparts. 
We defer further investigation of slowly rotating convection to future work.


The theory above was  proposed by \cite{Stevenson1979} 
(see his Eq. 43)
who additionally considered the case when rotation and gravity are misaligned,
 as well as magnetic fields.  
We have rederived it here in a slightly simpler way in order to highlight  
the main assumptions so that they may be tested in simulations. In addition, Stevenson's work
appears to be insufficiently
 well-known or  accepted in the literature.
 For example, \cite{Julien2012} derive Eq. \ref{eq:ml1} in a different way by scaling
 to the state at marginal stability---similar to the approach used by \cite{Kraichnan1962} for non-rotating convection. (Note that Eq. 2 in \cite{Julien2012} is equivalent to our Eq. \ref{eq:ml1} if one replaces $\Delta T/H\rightarrow N_*^2$ inside their dimensionless numbers, as is appropriate if most of the temperature drop is not within a boundary layer).
 And \cite{Aubert2001}  rederive Eq.~\ref{eq:ml2} (see their Eq. 10); note that \cite{Hide1974} Eq. A32 also gives the same prediction for the rms velocity.
  We suspect the main reason for this neglect is that Eqs. \ref{eq:ml1}--\ref{eq:ml4} have been 
insufficiently tested by numerical simulations---although the  simulations by \cite{Julien2012}
using a simplified set of equations provide  support, and
the prediction for the convective velocity has been compared with numerical simulations by \cite{Gillet2006} using a  quasi-geostrophic model. 
Note that
  \cite{2011Icar..211.1258S}
     suggest that a relation different than Eq. \ref{eq:ml2}
(specifically, $v_z\sim (F/\Omega)^{1/2}$) might be applicable.
The main goal of this paper is to verify Equations (\ref{eq:ml1})--(\ref{eq:ml4}) definitively with
high resolution simulations, which we turn to now.


\section{Simulations With Heating and Cooling Zones}
\label{Results}

As described in \S \ref{sec:intro}, RBC produces
thin boundary layers adjacent to the top and bottom plates through which heat enters and exits the fluid
via conduction.  As one decreases the diffusivities
$\nu$ and $\kappa$ towards  astrophysically tiny values, the boundary layers
become increasingly thin, and hence difficult to resolve in numerical simulations.  
But the  boundary layers themselves are of little interest if one's goal is to explain
convection zones in stars and planets.  
  To avoid this difficulty, we 
  directly heat a   zone of fluid at the bottom of the simulation domain
($z\in[0,0.2]$)
and cool a zone at the top ($z\in[1.2,1.4]$).
  That allows us to pass a flux of our choosing through the ``convection zone''
  ($z\in[0.2, 1.2]$) without suffering from overly thin boundary layers \footnote{ 
  Another possible route to avoiding boundary layers is to employ periodic boundary conditions  in the vertical
  direction  (e.g.~\citealt{Lohse2003,Garaud2010}).  However, such
  simulations develop large-scale up-flows that exit at
  the top of the box and enter at the bottom, as well as oppositely directed down-flows.
  These flows typically dominate the heat transport. 
    By contrast, in our setup  fluid is forced to turn around before reaching the top or bottom
  of the box, as in realistic convection zones.}.
  In the heating and cooling zones, the dynamics will depend 
  on the heating/cooling prescription.
  But  in the convection zone, the dynamics should be universal when 
  the diffusivities are small enough.  That will be tested both with convergence studies as
 $\nu$ and $\kappa$ are decreased, and with a few direct comparisons with
   standard Rayleigh-Benard simulations (in the next section).
  
  For most of our runs we integrate the equations of motion    (Eqs. \ref{eq:eoma1}--\ref{eq:eoma2} in the Appendix) with
   the pseudospectral code SNOOPY 
  \citep{Lesur2005}, and also rerun some of the simulations 
  with a
  very different code, the spectral element code  Nek5000 \citep{nek5000}, as an independent check.
  We set $\nu=\kappa$ throughout, and 
  boundary conditions are horizontally periodic and vertically impenetrable.
  Further numerical details are described in   \S \ref{methods}, including
  the implementation of heating/cooling, as well as a table of 
  simulations
  (Table~\ref{Table}). The choice $\nu=\kappa$ is made for simplicity. However, if the bulk properties of the convection are independent of the diffusivities (as predicted by Eqs.~\ref{eq:ml1}--\ref{eq:ml4} and verified
  below in our simulations), then they should not depend on this ratio.
   Nonetheless, it would be desirable to verify this directly with future simulations.


  We choose the amount of  heating/cooling so that the total heat flux
  (conductive + convective) through the 
   convection zone is unity in steady state.  The heat flux is defined as
   \begin{eqnarray}
   F&\equiv& {-\kappa}{d\left<T\right>\over dz}+\left< v_zT\right>,
   \end{eqnarray}
  where  angled brackets represent averaging in the horizontal direction at fixed height.  
  In steady state, the time-averaged $F$ is independent of height in the convection zone, 
  and we choose its value to be \textbf{ $F=1$}.
    That, together with our choice for the depth of the convection 
  zone ($H=1$), sets the simulation units\footnote{
  Heat flux
   has  units $[F]$= [{\it length}$]^2$/[{\it time}]$^3$ in the Boussinesq equations, and
   temperature has units $[T]=[{\it length}]/[{\it time}]^2$.}.
   Equivalently, lengths are measured in units of $H$ and time in units of $H^{2/3}/F^{1/3}$.
  The predictions of mixing length theory (Eqs. \ref{eq:ml1}--\ref{eq:ml4}) become particularly
  simple in our chosen units, and depend solely on $\Omega$.
 
\begin{figure}
  \begin{center}
   \hspace*{\fill}
    \subfigure{\includegraphics[trim=0.5cm 6cm 8cm 1cm, clip=true,width=0.48\textwidth]{Fig1a}} \hfill
   \hspace*{\fill}
     \subfigure{\includegraphics[trim=0.5cm 6cm 8cm 1cm, clip=true,width=0.48\textwidth]{Fig1b}} \hfill
    \end{center}
  \caption{
   Top: 
  Steady-state temperature profiles in four simulations 
  (labelled 3M, 10M, 30L and 100L in Table~\ref{Table})
  with different values of $\Omega$, averaged over horizontal planes and in time. 
   The shading around each curve  represents the corresponding RMS fluctuation. 
Temporal averaging is over at least 50 time units.
   Bottom: 
   The black lines show the slopes of the four profiles in the top panel.
   Also shown as colored points are comparison simulations with 
   higher diffusivities (blue for ones labelled M in Table~\ref{Table} and light blue
   if labelled H).
  Our simulation units throughout this paper are determined by setting $F=H=1$.
  }
  \label{2}
\end{figure}


\begin{figure*}
  \begin{center}
    \hspace*{\fill}
        \subfigure{\includegraphics[trim=8cm 6.5cm 8cm 4cm, clip=true,width=0.46\textwidth]{Om10Key} }
      \hspace*{\fill}
     \subfigure{\includegraphics[trim=8cm 6.5cm 8cm 5cm, clip=true,width=0.46\textwidth]{Om1000Key} }
    \end{center}
  \caption{Vertical velocity in a simulation with $\Omega=10$ (left panel; simulation 10M) and
  $\Omega=1000$ (right panel; 1000L), both at $t=50$.
  In the right panel, the horizontal scale of the box has been stretched by a factor of $100^{3/5}$ relative to the left panel, and the colour scale has similarly been scaled by $100^{-1/5}$. These scalings correspond with the predictions in \S \ref{mlt}.}
  \label{1}
\end{figure*}




The top panel of Fig.~\ref{2} shows the  averaged 
temperature profiles
in steady state for four  simulations with different rotation rates.
 The temperature gradient in the convection zone
   steepens as $\Omega$ is increased. That is because rotation
   hinders convection: to achieve a given flux ($F=1$) requires a steeper gradient
   in a more rapidly rotating simulation. 
 The bottom panel of Fig.~\ref{2}  shows the slopes of the profiles in the top panel (black curves), 
 along with those from  simulations with higher diffusivities (colored points). 
The fact that the points agree with the  curves shows that
we are  probing the regime in which the bulk properties are independent (or at worst  weakly dependent) on the microscopic diffusion coefficients.\footnote{
The  importance of thermal diffusivity in the convection zone can be quantified by the ratio 
of conductive to total flux:
 $F_{\rm cond}/F=\kappa d\left<T\right>/dz$ (since $F=1$).  For example, in the 
 $\Omega=10$
 simulation, we find at the midplane $F_{\rm cond}/F \approx 10^{-3.3} \times 11\approx 0.006$. 
 The smallness of this ratio suggests that diffusivities play little role in the convection zone.
}

Fig.~\ref{1} shows a snapshot of $v_{z}$ in two simulations, 
one with $\Omega=10$ and the other with $\Omega=1000$.  
For the right panel, we  scaled the horizontal length scale, as well as the colour scale, by the amount predicted by the arguments in \S \ref{mlt}, relative to the left panel. 
 The similarity of the flow in both panels provides  support for the 
mixing length theory.
The dominant convective modes occur on smaller horizontal length-scales for more rapid rotation (Eq.~\ref{eq:ml4}), and the corresponding vertical velocities decrease (Eq.~\ref{eq:ml2}).

 Fig.~\ref{3} shows the spectrum of the heat flux 
  in three $\Omega=10$ simulations that have different diffusivities. 
 The heat flux in these simulations is dominated by wavenumbers near $k_\perp\sim 20-30$. 
 When $\nu$ and $\kappa$ are decreased, the spectrum does not change near those lengthscales, indicating
 that the modes that dominate the heat flux are well-resolved
 and hence  bulk properties are independent of $\nu$ and $\kappa$.
The convection is anisotropic, since the depth of the convection zone is considerably
larger than the scale of modes that dominate the heat flux.  This anisotropy
is evident in Fig. \ref{1}, and is 
  expected from linear theory (Eq. \ref{eq:kperp}).
We also deduce from Fig.~\ref{3} that
 when $\nu$ is decreased, 
  the inertial range is extended to smaller scales.
   But those small scales  have little influence on the
larger scales that carry the bulk of the heat transport. 
These deductions conform  with the expectation from \S \ref{mlt}.
 
\begin{figure}
  \begin{center}
      \subfigure{\includegraphics[trim=0cm 0cm 0cm 0cm, clip=true,width=0.48\textwidth]{spectrum.pdf} }
  \end{center}
  \caption{
   Spectrum of the heat flux for three simulations with 
   $\Omega=10$ and
   different viscosities.
   The ordinate is  $k_\perp E(k_\perp)= k_\perp^2
   \mathrm{Re}\left<
    \hat{v}_{z}(\bld{k}_{\perp},z,t)\hat{T}(\bld{k}_{\perp},z,t)^{*}\right>$ (within a constant), where hats denote Fourier transforms in the horizontal directions with horizontal wavevector $\bld{k_\perp}$, and
      the averaging is threefold: (i) over $z$  
      inside the convection zone;
      (ii) over the orientation of $\bld{k}_\perp$; and
       (iii)  over time until convergence.
       Note that the spatially averaged flux is $F={\rm const}\times\int E(k_\perp)dk_\perp$, 
       so that the peak of the spectrum plotted indicates the wavenumbers that dominate
       the flux. 
       The $-2/3$ line corresponds to the Kolmogorov scaling, i.e., 
        the assumption that 
         fluctuations across lengthscale $\lambda_\perp$ scale as
       $\delta v_z, \delta T\propto\lambda_\perp^{1/3}$, and that
       the two are strongly correlated. 
     The jaggedness is due to the discrete nature of the simulation grid, not temporal
     fluctuations which have been averaged away. }
 \label{3}
\end{figure}


\begin{figure*}
 \begin{tabular}{l r}
 \subfigure{\includegraphics[trim=0cm 0cm 0cm 0cm, clip=true,width=0.48\textwidth]{Fig4a} } &
  \subfigure{\includegraphics[trim=0cm 0cm 0cm 0cm, clip=true,width=0.48\textwidth]{Fig4b} } \\
  \subfigure{\includegraphics[trim=2cm 0cm 0cm 0cm, clip=true,width=0.47\textwidth]{Fig4c} } &
   \subfigure{\includegraphics[trim=0cm 0cm 0cm 0cm, clip=true,width=0.48\textwidth]{Fig4d} }
  \end{tabular}
  \caption{Dependence of  bulk properties   on  rotation rate.
   The predictions from mixing length theory  are shown as lines.
   The points labelled
    L, M, and H are SNOOPY simulations with viscosities from low to high; 
    red points are Nek5000 results; and error bars are RMS fluctuations. Points with the same $\Omega$ are slightly offset horizontally for clarity. Our simulation units are determined by setting $F=H=1$.}
  \label{5}
\end{figure*}

Figure~\ref{5}, the main result of this paper, displays the bulk properties 
in all of the simulations listed in Table \ref{Table}.
The top-left panel shows  the average temperature gradient in the convection zone, 
which 
 we extract  by fitting  the temporally and horizontally averaged temperature profile in the central $\frac{1}{3}$ of the simulation domain with a straight line. The gradient is shown as a point and the  corresponding RMS fluctuation as error bars, after averaging over at least 50 time units.
 At each $\Omega$, the results from the different simulations---SNOOPY with
 various viscosities and Nek---agree quite well  with the prediction, shown as a solid line.
 We fit the points and RMS errors 
 for all simulations
 with a linear least squares fit in log-log space, keeping points with
  $\Omega>6$ only (i.e.~the rapidly rotating limit).
 We find 
 \begin{equation}
  -{d\left<T\right>\over dz}=1.9^{+0.3}_{-0.3}\; \Omega^{0.75\pm0.06}.
  \end{equation}
   The exponent agrees with the theoretical prediction  of 0.8 (Eq. \ref{eq:ml1}) within 
   the error bars.  The lower-left panel of Fig.~\ref{5} plots the same data after removing the predicted scaling. 
 The prediction  works remarkably well
   over more than two orders of magnitude; exponents that differ from 0.8 by more than $\sim 10\%$ are definitively ruled out.

The top-right panel of Fig.~\ref{5} shows  the RMS fluctuations in $v_z$ and $T$ at the 
midplane of the box $(z=L_z/2)$. These also agree very well with the predicted
scalings, shown as lines.
A least-squares fit to the simulation points gives
\begin{eqnarray}
\delta v_{z}&=&0.9^{+0.5}_{-0.3} \; \Omega^{-0.15\pm0.08}\ ,
\\
\delta T &=& 1.5^{+0.9}_{-0.6} \; \Omega^{0.21\pm0.13} \ ,
\end{eqnarray}
the exponents of which may be compared with the predictions 
of -0.2 for $\delta v_z$ (Eq. \ref{eq:ml2}), and +0.2 for $\delta T$ (Eq. \ref{eq:ml3}).

In the fourth panel of Fig.~\ref{5}, we plot the horizontal wavenumber that dominates the heat flux, $\hat{k}$.  We define it via
\begin{equation}
\hat{k}={\int k_\perp E(k_\perp) d\ln k_\perp\over \int E(k_\perp)d\ln k_\perp} \ ,
\end{equation}
where the flux spectrum $E(k_\perp)$ is defined in the caption of Fig. \ref{3}.  Note that both the
numerator and denominator of this expression are dominated by large scales. 
A least-squares fit to the plotted values yields
\begin{equation}
\hat{k}=7.3^{+7.2}_{-3.6}\; \Omega^{0.58\pm0.19}  \ ,
\end{equation}
which may be compared with the mixing length prediction $k_\perp\propto \Omega^{0.6}$
(Eq.~\ref{eq:ml4}).


The results of this section, and in particular Fig.~\ref{5}, provides strong support for the
mixing length theory  presented in \S \ref{Theory}.




\section{Simulations Without Internal Heating and Cooling}
\label{RBC}

For the simulations presented in \S \ref{Results}, we directly heated and cooled the fluid 
inside the simulation domain to avoid thin boundary layers.
That approach was 
predicated on the assumption that the details of how
heat enters and leaves the convection zone is only of  minor importance for determining the bulk properties. In this section, 
we test that assumption.  To do so, we run comparison simulations 
without any internal heating or cooling.
The temperature
is held fixed at the bottom boundary and the flux is fixed at the top, which is essentially
the setup for standard Rayleigh-Benard convection. 
 Since the fluid is driven by heating at the boundaries, it is essential to correctly resolve the thermal boundary layers
  -- failing to do so results in an incorrect heat flux through the domain \citep{Groetzbach1983,Shishkina2010}. When the diffusivities are small, 
   the boundary layers become very thin, and hence  computationally costly to resolve.

\begin{figure}
  \begin{center}
 \begin{tabular}{l r}
 \subfigure{\includegraphics[trim=1cm 8cm 8cm 1cm, clip=true,width=0.48\textwidth]{Fig5a} } \\
 \hspace{-0.4cm} \subfigure{\includegraphics[trim=0.6cm 8cm 8cm 2cm, clip=true,width=0.49\textwidth]{Fig5b} }
  \end{tabular}
     \end{center}
  \caption{
  Temperature profiles in simulations with (black) and without (red) internal heating at two 
  rotation rates.  Those with internal heating are 10M and 30M.  The ones with
  no internal heating were done with Nek5000. 
  The temperature was fixed 
  at the bottom and the flux fixed  at the top to $F=1$. We set $L_z=1.2$ in order to 
   crudely account for the fact that in the internally heated simulations the convection zone 
  extends somewhat into the heating/cooling zones. 
   The simulation without internal heating that has $\Omega=10$
   also has $L_{\perp}=1$,  $\nu=\kappa=10^{-3}$, and
  uses $200^3$ points split into 20 elements and 10th order polynomials within each element per dimension; the simulation with $\Omega=30$ has
  $L_{\perp}=0.5$, $\nu=\kappa=10^{-3.6}$, and similarly uses $200^3$ points.
 }
  \label{8}
\end{figure}

The equations of motion (Eqs. \ref{eq:eoma1}--\ref{eq:eoma2} with  $q=0$) were integrated
  with Nek5000, which is better suited than SNOOPY for resolving  thin boundary layers.   This is because grid points  in Nek5000 are clustered towards the boundaries, 
  whereas in SNOOPY sharp boundary layers produce unwanted Gibbs oscillations.
    Fig.~\ref{8} shows the temperature profiles and gradients from two such simulations,  with $\Omega=10$
and 30, and compares them with the corresponding heating/cooling zone simulations. 
The agreement
is  reasonable throughout the convection zone,  thus confirming our assumption\footnote{The agreement is not 
perfect primarily because of the difference in the depth of the convecting layer.}. The thin boundary layers are also evident in these figures.


\section{Discussion}
\label{Conclusions}

We  presented a simple derivation of mixing length theory 
in rapidly rotating convection, and then verified it with simulations.
The theory, postulated
by \cite{Stevenson1979},
 predicts the  properties of the convecting fluid under the assumption
that they are independent of microscopic diffusivities ($\nu$ and $\kappa$).
Equations \ref{eq:ml1}--\ref{eq:ml4} list the predictions for the mean temperature 
gradient,  the velocity and temperature fluctuations, and  the lengthscale of the modes that
dominate heat transport. 
Our simulation results, summarized in  
Figure \ref{5},  agree remarkably well with the theory, 
 across more than two orders of magnitude in rotation rate. 

We chose to focus on a very simple setup: Boussinesq convection in a box.
But despite its simplicity, and despite the vast literature already devoted to the topic,
 the result remains
under debate  \citep[e.g.,][]{King2012,Julien2012}, largely because of the complicating effect of
boundary layers. 
We circumvented this complication  by focusing 
on the properties of the convecting fluid---i.e., between the boundary 
layers.  
We did this by fixing the flux, 
and examining the  interior fluid's properties 
 for increasingly small  diffusivities.
 We thereby showed that the  convecting fluid's properties converged to the prediction 
of mixing length theory as $\nu$, $\kappa\rightarrow 0$.
Moreover, 
the numerical resolutions required to demonstrate convergence were relatively modest,
after artificially thickening the boundary layers with heating/cooling zones. 
  For example,
 our SNOOPY simulations had 256$^3$ gridpoints or fewer. 
Our numerical results  provide strong  support for  those   of  \cite{Julien2012}, who simulate
 a set of reduced equations  valid in the limit of rapid rotation.
 
Our work 
  lends confidence to mixing length theory's ability to accurately
 model highly turbulent convection. 
 We hope 
to extend it to include a variety of more complicated---and  realistic---effects, some
 of the most important of which are as follows:
\begin{itemize}
 	\item Including a background density gradient. One must then distinguish
	between entropy and  temperature.  The argument presented in \S \ref{mlt} should
	remain largely unchanged,  after replacing temperature with entropy.  
	But a possible complication is the asymmetry between upflows and downflows in the
	presence of 
	a density gradient (e.g.~\citealt{Hurlburt1984,Cattaneo1991,Miesch2005}).
	\item Allowing for a more realistic geometry, i.e., quasi-spherical rather than cubical.  
	The work presented in this paper strictly applies only to a small patch of the convective region near the poles of the star or planet. 
	An intermediate step before considering the full spherical problem would be
	 to allow rotation and gravity to be misaligned. \cite{Stevenson1979} predicts that in that case 
	Eqs. \ref{eq:ml1}--\ref{eq:ml4} should be altered by replacing $\Omega\rightarrow\Omega\cos\theta$, where $\theta$ is co-latitude. But that has yet to be confirmed by simulations.
	\item  Including the boundaries of a convection zone and the possibilities of 
	penetration and
	overshooting  into neighbouring stable layers (e.g.~\citealt{Hurlburt1994,Brummell2002,Rogers2005})
	\item  Allowing for the interaction between convection and differential rotation, and the generation of secondary flows.
	\item  Including  magnetic fields.
\end{itemize}


A large body of work has already been devoted to simulations of convection in rotating stars and planets, from Boussinesq (e.g.~\citealt{Hathaway1983,Tilgner2009,King2012}) to fully compressible \citep{Brummell1996,Brummell1998,Kapyla2005} Cartesian box simulations to Boussinesq \citep{Christensen2002,ChristensenAubert2006}, anelastic \citep{Glatz1984,Miesch2000,Kaspi2009,Jones2009,Gastine2012} and fully compressible simulations \citep{Kapyla2011} in spherical shell geometry.  Given the evident complexities of some of these simulations, it is our view that a more complete understanding of simpler models is required to enable us to understand these simulation results. Our work complements the literature by definitively verifying the rotating mixing length theory described in \S \ref{mlt} for the case of Boussinesq convection in the polar regions of a planet or star. We anticipate that the theory described in this paper, as well as the extensions discussed above, will help provide a theoretical basis for the simulation results.


Turning  to astrophysical applications of the theory, we note first that 
rotation
changes the entropy gradient relative to that
predicted by  standard (non-rotating) mixing length theory
by an order-unity factor---at least for the Sun, where
the rotation rate is comparable to the convective turnover time.
Thus the inclusion of rotation
 will not substantially change static
structure calculations, 
since 
it hardly affects the conclusion
that convection zones have a near-constant entropy throughout \citep{Stevenson1979}. But a potentially important application is explaining the differential rotation profile of the Sun and other stars. 
In particular, small latitudinal entropy gradients drive differential rotation via the thermal wind equation
 (e.g.~\citealt{Thompson2003,Miesch2006,Balbus2009a,Balbus2009b}). 
 Therefore to predict the differential rotation profile from first principles requires 
 one to understand how the entropy gradient depends on latitude.  It appears likely that 
 rotating mixing length theory (at least when extended to the case in which rotation and gravity are misaligned) will provide an important piece towards solving this puzzle.
 
 Another potential application is to tidal dissipation in a convective star or planet that has
  an orbiting companion.\footnote{We thank Jeremy Goodman for pointing out this application to us.}
This is important for understanding, for example, the tidal circularization of  solar-type binary stars out to approximately ten day orbits.
Previous  work estimates the turbulent viscosity due to convection by employing non-rotating
mixing length theory
\citep{Zahn1966,GN1977}.  
It would be of interest to see how the predictions are affected by employing rotating mixing
length theory.
 
 
\section*{Acknowledgments}
We  thank Jonathan Aurnou, Keith Julien and the referee for suggestions which have improved the manuscript.
YL acknowledges the support of NSF grant AST-1109776 and NASA grant NNX14AD21G. The computations in this paper were performed on Northwestern University's HPC cluster Quest.

