# Predicting the Rossby number in convective experiments

Evan H. Anders Dept. Astrophysical & Planetary Sciences, University of
Colorado – Boulder, Boulder, CO 80309, USA  Laboratory for Atmospheric and
Space Physics, Boulder, CO 80303, USA  Cathryn M. Manduca  Laboratory for
Atmospheric and Space Physics, Boulder, CO 80303, USA  Benjamin P. Brown Dept.
Astrophysical & Planetary Sciences, University of Colorado – Boulder, Boulder,
CO 80309, USA  Laboratory for Atmospheric and Space Physics, Boulder, CO
80303, USA  Jeffrey S. Oishi Department of Physics and Astronomy, Bates
College, Lewiston, ME 04240, USA  Geoffrey M. Vasil  University of Sydney
School of Mathematics and Statistics, Sydney, NSW 2006, Australia

(Received October 19, 2018; Revised January 14, 2019; Accepted January 15,
2019)

###### Abstract

The Rossby number is a crucial parameter describing the degree of rotational
constraint on the convective dynamics in stars and planets. However, it is not
an input to computational models of convection but must be measured ex post
facto. Here, we report the discovery of a new quantity, the Predictive Rossby
number, which is both tightly correlated with the Rossby number and specified
in terms of common inputs to numerical models. The Predictive Rossby number
can be specified independent of Rayleigh number, allowing suites of numerical
solutions to separate the degree of rotational constraint from the strength of
the driving of convection. We examine the scaling of convective transport in
terms of the Nusselt number and the degree of turbulence in terms of the
Reynolds number of the flow, and we find scaling laws nearly identical to
those in nonrotational convection at low Rossby number when the Predictive
Rossby number is held constant. Finally, we describe the boundary layers as a
function of increasing turbulence at constant Rossby number.

convection — hydrodynamics — turbulence — dynamo — Sun: rotation

††journal: ApJ

##  1 Introduction

Rotation influences the dynamics of convective flows in stellar and planetary
atmospheres. Many studies on the fundamental nature of rotating convection in
both laboratory and numerical settings have provided great insight into the
properties of convection in both the rapidly rotating regime and the
transition to the rotationally unconstrained regime (King et al., 2009; Zhong
et al., 2009; Schmitz & Tilgner, 2009; King et al., 2012; Julien et al., 2012;
King et al., 2013; Ecke & Niemela, 2014; Stellmach et al., 2014; Cheng et al.,
2015; Gastine et al., 2016) The scaling behavior of heat transport, the nature
of convective flow structures, and the importance of boundary layer-bulk
interactions in driving dynamics are well known. Yet, we do not know of any
simple procedure for predicting the magnitude of vortical flow gradients
purely from experimental control parameters, such as bulk rotation rate and
thermal input.

In the astrophysical context, many studies of rotating convection have
investigated questions inspired by the solar dynamo (Glatzmaier & Gilman,
1982; Busse, 2002; Brown et al., 2008, 2010, 2011; Augustson et al., 2012;
Guerrero et al., 2013; Käpylä et al., 2014). Even when these simulations
nominally rotate at the solar rate, they frequently produce distinctly
different behaviors than the true Sun, such as anti-solar differential
rotation profiles (Gastine et al., 2014; Brun et al., 2017). It seems that
these differences occur because the simulations produce less rotationally
constrained states than the Sun. The influence of rotation results from the
local shear gradients, and these are not direct input parameters. Recent
simulations predict significant rotational influence in the deep solar
interior, which can drastically affect flows throughout the solar convection
zone (Featherstone & Hindman, 2016; Greer et al., 2016). In the planetary
context, the balance between magnetic and rotational forces likely leads to
the observed differences between ice giant and gas giant dynamos in our solar
system (Soderlund et al., 2015). The work of Aurnou & King (2017) demonstrates
the importance of studying a dynamical regime with the proper balance between
Lorentz, Coriolis, and inertial forces when modeling astrophysical objects
such as planetary dynamos.

In short, simulations must achieve the proper rotational balance if they are
to explain the behavior of astrophysical objects. In Boussinesq studies,
rotational constraint is often measured by comparing dynamical and thermal
boundary layers or deviation in heat transport from the non-rotating state
(King et al., 2012; Julien et al., 2012; King et al., 2013). Such measurements
are not available for astrophysical objects, where the degree of rotational
influence is best assessed by the ratio between nonlinear advection magnitude
and the linear Coriolis accelerations. The Rossby number is the standard
measure of this ratio,

$\text{Ro}\ \equiv\ \frac{|\nabla\times\bm{u}|}{2|\bm{\Omega}|}\ \sim\
\frac{|(\nabla\times\bm{u})\times\bm{u}|}{|2\bm{\Omega}\times\bm{u}|},$ (1)

where $\bm{\Omega}$ denotes the bulk rotation vector. Many proxies for the
dynamical Rossby number exist that are based solely on input parameters, most
notably the convective Rossby number. However, all proxies produce imperfect
predictions for the true dynamically relevant quantity.

> _In this letter, we demonstrate an emperical method of predicting the output
> Rossby number of convection in a simple stratified system._

In Anders & Brown (2017) (hereafter AB17), we studied non-rotating
compressible convection without magnetic fields in polytropic atmospheres. In
this work, we extend AB17 to rotationally-influenced, $f$-plane atmospheres
(e.g., Brummell et al., 1996, 1998; Calkins et al., 2015). We determine how
the input parameters we studied previously, which controlled the Mach and
Reynolds numbers of the evolved flows, couple with the Taylor number (Ta,
Julien et al., 1996), which sets the magnitude of the rotational vector.

In section 2, we describe our experiment and paths through parameter space. In
section 3, we present the results of our experiments and in section 4 we offer
concluding remarks.

##  2 Experiment

We study fully compressible, stratified convection under precisely the same
atmospheric model as in AB17, but here we have included rotation. We study
polytropic atmospheres with $n_{\rho}=3$ density scale heights and a
superadiabatic excess of $\epsilon=10^{-4}$ such that flows are at low Mach
number. We study a domain in which the gravity, $\bm{g}=-g\hat{z}$, and
rotational vector, $\bm{\Omega}=\Omega\hat{z}$, are antiparallel (as in e.g.,
Julien et al., 1996; Brummell et al., 1996).

We evolve the velocity ($\bm{u}$), temperature ($T$), and log density
($\ln\rho$) according to the Fully Compressible Navier-Stokes equations in the
same form presented in AB17, with the addition of the Coriolis term,
$2\bm{\Omega}\times\bm{u}$, to the left-hand side of the momentum equation. We
impose impenetrable, stress-free, fixed-temperature boundary conditions at the
top and bottom of the domain.

We specify the kinematic viscosity ($\nu$), thermal diffusivity ($\chi$), and
strength of rotation ($\Omega$) at the top of the domain by choosing the
Rayleigh number (Ra), Prandtl number (Pr), and Taylor number (Ta),

$\text{Ra}=\frac{gL_{z}^{3}\Delta
S/c_{P}}{\nu\chi},\,\,\,\text{Pr}=\frac{\nu}{\chi},\,\,\,\text{Ta}=\left(\frac{2\Omega
L_{z}^{2}}{\nu}\right)^{2},$ (2)

where $L_{z}$ is the depth of the domain as defined in AB17, $\Delta
S\propto\epsilon n_{\rho}$ is the specific entropy difference between the top
and bottom of the atmosphere, and the specific heat at constant pressure is
$c_{P}=\gamma/(\gamma-1)$ with $\gamma=5/3$. Throughout this work we set Pr =
1. The Taylor number relates to the often-quoted Ekman number by the equality
$\text{Ek}\equiv\text{Ta}^{-1/2}$.

Due to stratification, Ra and Ta both grow with depth as
(Ra,Ta)$\,\propto\rho^{2}$ (see AB17). We nondimensionalize our atmospheres at
the top of the domain, and so all values of Ra and Ta quoted in this work are
the minimal value of Ra and Ta in the domain at $z=L_{z}$. For direct
comparison to Boussinesq studies, past work has found that the value of Ra at
the atmospheric midplane ($z=L_{z}/2$) varies minimally with increasing
stratification (Unno et al., 1960). For the atmospheres presented in this
work, midplane Ra and Ta values are larger than reported top-of-atmosphere
values by a factor of $\sim$70, and values at the bottom of the atmosphere are
larger by $\sim$400.

When Ta is large, the wavenumber of convective onset increases according to
$k_{\text{crit}}\propto\text{Ta}^{1/6}$ (Chandrasekhar, 1961; Calkins et al.,
2015). We study horizontally-periodic, 3D Cartesian domains with extents of
$x,y=[0,4(2\pi/k_{\text{crit}})]$ and $z=[0,L_{z}]$. At large values of Ta,
these domains are tall and skinny, as in Stellmach et al. (2014). We evolve
our simulations using the Dedalus111http://dedalus-project.org/ pseudospectral
framework, and our numerical methods are identical to those presented in AB17.
The supplemental materials of this paper include a .tar file which contains
the code used to perform the simulations in this work, and this tarball is
also published online in a Zenodo repository (Anders et al., 2019).

Figure 1: (a) The critical Rayleigh number, as a function of the Taylor
number, is plotted as a solid black line. The grey shaded region is
subcritical, and rotation supresses convection there. Paths of constant
Convective Rossby number ($\text{Ro}_{\text{c}}$, green dashed line), constant
supercriticality ($\mathcal{S}$, blue dash-dot line), and constant Predictive
Rossby number ($\text{Ro}_{\text{p}}$, orange solid lines) are shown. From
thickest to thinnest, paths with $\text{Ro}_{\text{p}}=[1.58,0.96,0.6]$ are
plotted, and the value of ($\text{Ta}_{\text{crit}},\text{Ra}_{\text{crit}}$)
for each path is denoted by a circular marker (see Table 1). (b) Evolved Ro is
plotted vs. Ra along paths of $\text{Ro}_{\text{p}}$= [1.58, 0.96, 0.6] for
[big, medium, small] orange circles. For comparison, paths of constant
$\mathcal{S}$ (blue squares, $\mathcal{S}=[3,2]$ for [big, small] squares) and
constant $\text{Ro}_{\text{c}}$ (green triangles, $\text{Ro}_{\text{c}}$ = [1,
0.3, 0.1] for [big, medium, small] triangles) are shown. (c) The evolved value
of Ro is shown as a function of $\text{Ro}_{\text{p}}$ and
$\text{Ro}_{\text{c}}$. Each of the experiments in (b) is outlined by a black
(circle, triangle, square) for points along constant ($\text{Ro}_{\text{p}}$,
$\text{Ro}_{\text{c}}$, $\mathcal{S}$) paths. The color inside of the marker
represents the exact measured Ro of that experiment, while the colormap
outside of markers is a linear interpolation of the data set.

The critical value of Ra at which rapidly rotating convection onsets also
depends on Ta (see the black line in figure 1a), roughly according to
$\text{Ra}_{\text{crit}}\sim\text{Ta}^{2/3}$ (Chandrasekhar, 1961; Calkins et
al., 2015). Even taking account of linear theory, the dependence of the
evolved nonlinear fluid flows on the input parameters makes predicting the
rotational constraint very challenging. We will explore three paths through
Ra-Ta space:

$\text{Ra}=\begin{cases}\mathcal{S}\,\text{Ra}_{\text{crit}}(\text{Ta}),&(\text{I})\\\
(\text{Ro}_{\text{c}})^{2}\,\text{Pr}\,\text{Ta},&(\text{II})\\\
(\text{Ro}_{\text{p}})^{2}\,\text{Pr}^{1/2}\,\text{Ta}^{3/4}&(\text{III}).\end{cases}$
(3)

Paths on constraint I are at constant supercriticality,
$\mathcal{S}\equiv\text{Ra}/\text{Ra}_{\text{crit}}(\text{Ta})$ (blue dash-dot
line in figure 1a). Paths on constraint II (green dashed line in figure 1a)
are at a constant value of the classic convective Rossby number,

$\text{Ro}_{\text{c}}\ =\ \sqrt{\frac{\text{Ra}}{\text{Pr}\,\text{Ta}}}\ =\
\frac{1}{2\Omega}\sqrt{\frac{g\,\Delta S}{c_{p}L_{z}}},$ (4)

which has provided (e.g., Julien et al., 1996; Brummell et al., 1996) a common
proxy for the degree of rotational constraint. This parameter measures the
importance of buoyancy relative to rotation without involving dissipation.
Paths on constraint III (e.g., orange solid lines in figure 1a) set constant a
ratio which we call the “Predictive Rossby number,”

$\text{Ro}_{\text{p}}=\sqrt{\frac{\text{Ra}}{\text{Pr}^{1/2}\,\text{Ta}^{3/4}}}\
=\ \frac{1}{(2\Omega)^{3/4}}\sqrt{\frac{g\,\Delta S}{c_{p}\chi^{1/2}}}$ (5)

Unlike paths through parameter space which hold $\text{Ro}_{\text{c}}$
constant, paths with constant $\text{Ro}_{\text{p}}$ feel changes in
diffusivities but not the depth of the domain. To our knowledge, these paths
have not been reported in the literature, although the importance of
$\text{Ra}/\text{Ta}^{3/4}=\text{Ra Ek}^{3/2}$ has been independently found by
King et al. (2012) using a boundary layer analysis. We compare our results to
their theory in Section 4.

In this work, we primarily study three values of $\text{Ro}_{\text{p}}$. These
values are shown in Fig. 1a and Table 1. Table 1 lists the values of
(Ra${}_{\text{crit}}$, Ta${}_{\text{crit}}$) for each value of
$\text{Ro}_{\text{p}}$, and also the maximum value of (Ra, Ta) studied in this
work for each path. We additionally walked two pathways at constant
supercriticality (constraint I, $\mathcal{S}=\\{2,3\\}$) and three pathways at
constant convective Rossby number (constraint II,
$\text{Ro}_{\text{c}}=\\{1,0.3,0.1\\}$). Full details on all cases are
provided in Appendix A and the supplemental materials.

Table 1: Parameter space

$\text{Ro}_{\text{p}}$ | (Ra${}_{\text{crit}}$, Ta${}_{\text{crit}}$) | (Ra${}_{\text{max}}$, Ta${}_{\text{max}}$)
---|---|---
0.60 | ($10^{4.88},10^{7.10}$) | ($10^{9.09},10^{12.72}$)
0.96 | ($10^{2.44},10^{3.30}$) | ($10^{8.58},10^{11.49}$)
1.58 | ($10^{1.39},10^{1.33}$) | ($10^{7.14},10^{8.99}$)

Note. — Values of the critical Ra and Ta for each $\text{Ro}_{\text{p}}$ track
are reported, as well as the maximal values of Ra and Ta studied on each
track. All values reported are for the top of the atmosphere. A fuller set of
simulations are reported in Table 2 with midplane Ra and Ta values as well.

##  3 Results

In our stratified domains, for $\text{Ta}\geq 10^{5}$, a best-fit to results
from a linear stability analysis provides
$\text{Ra}_{\text{crit}}(\text{Ta})=1.459\text{Ta}^{2/3}$ and
$k_{\text{crit}}(\text{Ta})=0.414\text{Ta}^{1/6}$ for direct onset of
convection. In figure 1a, the value of Ra${}_{\text{crit}}(\text{Ta})$ is
shown. Sample paths for each criterion in equation 3 through this parameter
space are also shown. In this work, we often find it instructive to use one
critical Ra for an entire $\text{Ro}_{\text{p}}$ path. This
$\text{Ra}_{\text{crit}}$ is determined by the intersection of the onset curve
and $\text{Ro}_{\text{p}}$ path (indicated by the orange circles in figure 1a,
and quoted in Table 1). In the high Ta regime, we find that
$\text{Ra}_{\text{crit}}=18.5\text{Ro}_{\text{p}}^{-16}$.

In figure 1b, we display the evolution of Ro with increasing Ra along various
paths through parameter space. We find that Ro increases on constant
$\text{Ro}_{\text{c}}$ paths, decreases on constant $\mathcal{S}$ paths, and
remains roughly constant along constant $\text{Ro}_{\text{p}}$ paths. In
figure 1c, the value of Ro is shown simultaneously as a function of
$\text{Ro}_{\text{p}}$ and $\text{Ro}_{\text{c}}$ for all experiments
conducted in this study. We find a general power-law of the form
$\text{Ro}=\mathcal{C}\text{Ro}_{\text{c}}^{\alpha}\,\text{Ro}_{\text{p}}^{\beta}$.
In the rotationally-dominated regime where $\text{Ro}<0.2$ and
$\text{Re}_{\perp}>5$ (see Eqn. 6), we find $\alpha=-0.02$, and Ro can be said
to be a function of $\text{Ro}_{\text{p}}$ alone. Under this assumption, we
report a scaling of $\text{Ro}=(0.148\pm 0.003)\text{Ro}_{\text{p}}^{3.34\pm
0.07}$. In the less rotationally dominated regime of $\text{Ro}>0.2$ and
$\text{Re}_{\perp}>5$, we find $\\{C,\alpha,\beta\\}=\\{0.2,-0.19,1.5\\}$.

Figure 2:  Horizontal slices of the evolved entropy deviation from the mean at
$z=0.95L_{z}$ (top row) and $z=0.5L_{z}$ (bottom row) are shown for select
simulations. All runs displayed here have an evolved volume-averaged
$\text{Re}_{\perp}\approx 32$. As Ro decreases from O(1) on the left to
O(0.03) on the right, and thus the rotational constraint on the flow
increases, significant changes in flow morphology are observed. As Ro
decreases, Coriolis forces more effectively deflect the convective flows, and
the classic granular convective pattern gives way to vortical columns that are
quasi-two-dimensional.  Figure 3:  Scaling laws for paths at
$\text{Ro}_{\text{p}}=1.58$ ($\text{Ro}\approx 0.4$),
$\text{Ro}_{\text{p}}=0.96$ ($\text{Ro}\approx 0.1$), and
$\text{Ro}_{\text{p}}=0.6$ ($\text{Ro}\approx 0.03$) are shown. Numbers are
plotted vs. Ra/Ra${}_{\text{crit}}$, where Ra${}_{\text{crit}}$ is given in
Table 1. (a) Nu, as defined in AB17, is shown. (b) $\text{Re}_{\parallel}$ and
$\text{Re}_{\perp}$, as defined in equation 6, are shown. All values of
$\text{Ro}_{\text{p}}$ trace out similar Nu and $\text{Re}_{\perp}$ tracks,
whereas $\text{Re}_{\parallel}$ tracks shift upwards as Ro decreases.  Figure
4: Horizontally-averaged profiles of the standard deviation of entropy
($\sigma_{\text{s}}$, a) and Rossby number (Ro, b) are shown vs. height for
$\text{Ro}_{\text{p}}=1.58$ ($\text{Ro}\approx 0.4$). Similar profiles are
shown in (c) and (d) for $\text{Ro}_{\text{p}}=0.96$ ($\text{Ro}\approx 0.1$).
The color of the profiles denotes the value of Ra/Ra${}_{\text{crit}}$, where
Ra${}_{\text{crit}}$ is given in Table 1. (e) The ratio of the thicknesses of
the dynamical boundary layers ($\delta_{\text{Ro}}$) and thermal boundary
layers ($\delta_{\text{s}}$) is shown vs. $\text{Ra}/\text{Ra}_{\text{crit}}$
for fixed $\text{Ro}_{\text{p}}$. (f) $\delta_{\text{s}}$ is plotted vs.
$\text{Ra}/\text{Ra}_{\text{crit}}$ in units of the density scale height at
the top of the atmosphere ($H_{\rho}$). Vertical lines denote when
$\delta_{\text{s}}/H_{\rho}=1$ for each value of $\text{Ro}_{\text{p}}$.

In figure 2, sample snapshots of the evolved entropy field in the $x$-$y$
plane near the top and at the middle of the domain are shown. In the left
column, flows are at $\text{Ro}\sim 1$ and resemble the classic granular
structure of nonrotating convection (see e.g., figure 2 in AB17), where strong
narrow downflow lanes punctuate broad upwellings. The narrow downflows at the
top organize themselves into intense coherent structures at the midplane, and
at the midplane the downflows have much stronger entropy fluctuations than the
broad and slower upflows.

Ro decreases from left to right into the rotationally constrained regime. As
Ro decreases, the narrow downflow lanes begin to disappear and the flows at
midplane become more symmetric. In the rotationally constrained regime
($\text{Ro}\sim 0.03$), the convective structures are distinctly different.
Here we observe dynamically persistent, warm upflow columns surrounded by bulk
weak downflow regions. At the midplane, the upflow columns have substantially
higher entropy perturbations than the surrounding weak downflows which sheathe
them, and the locations of the columns are tightly correlated with their
positions at the top of the domain. These quasi-two-dimensional dynamics are
similar to those seen in rapidly rotating Rayleigh-Bénard convection (e.g.,
Stellmach et al., 2014). The select cases displayed in figure 2 each have an
evolved volume-averaged $\text{Re}_{\perp}\approx 32$ (defined below in
equation 6).

We measure the Nusselt number (Nu), which quantifies heat transport in a
convective solution, as defined in AB17. In figure 3a, we plot Nu as a
function of $\text{Ra}/\text{Ra}_{\text{crit}}$ at fixed
$\text{Ro}_{\text{p}}$. We find that $\text{Nu}\propto\\{\text{Ra}^{0.29\pm
0.01},\text{Ra}^{0.29\pm 0.01},\text{Ra}^{0.24}\\}$ for
$\text{Ro}_{\text{p}}=\\{0.6,0.957,1.58\\}$. In the regime of
$\text{Ro}\lesssim 0.1$, these scaling laws are indistinguishable from a
classic Ra2/7 power law scaling, which is observed in nonrotating Rayleigh-
Bénard and stratified convection (Ahlers et al., 2009, AB17). Our results seem
consistent with the stress-free, rotating Rayleigh-Bénard convection results
of Schmitz & Tilgner (2009), whose re-arranged Eqn. 7 returns a best-fit of
$\text{Nu}\propto\text{Ra}^{0.26}$ at fixed $\text{Ro}_{\text{p}}$ . Their
work primarily spans the transition regime between rotationally constrained
and unconstrained convection, and so it is perhaps not surprising that their
power law is a blend of our rotationally-constrained $\text{Ra}^{2/7}$ power
law and the fairly rotationally unconstrained $\text{Ra}^{0.24}$ at
$\text{Ro}_{\text{p}}=1.58$.

Flows are distinctly different parallel to and perpendicular from the rotation
vector, which aligns with gravity and stratification. We measure two forms of
the RMS Reynolds number,

$\text{Re}_{\parallel}=\frac{|\bm{u}|L_{z}}{\nu},\qquad\text{Re}_{\perp}=\frac{|\bm{u}|}{\nu}\frac{2\pi}{k_{\text{crit}}},$
(6)

where the length scale in $\text{Re}_{\perp}$ is the wavelength of convective
onset, and is related to the horizontal extent of our domain (see section 2).
From our work in AB17, we expect the RMS velocity to scale as
$|u|\propto\sqrt{\Delta S}$. By definition,
$\nu\propto\sqrt{\text{Ra}/(\text{Pr }\Delta S)}$, and $L_{z}$ is a constant
set by the stratification while $k_{\text{crit}}\propto\text{Ta}^{1/6}$. Along
paths of constant $\text{Ro}_{\text{p}}$, we thus expect
$\text{Re}_{\parallel}\propto\text{Ra}^{1/2}$ and
$\text{Re}_{\perp}\propto\text{Ra}^{5/18}$ when Pr is held constant.

In figure 3b, we plot $\text{Re}_{\parallel}$ and $\text{Re}_{\perp}$ as a
function of Ra$/\text{Ra}_{\text{crit}}$ at fixed $\text{Ro}_{\text{p}}$. We
find that $\text{Re}_{\parallel}\propto\\{\text{Ra}^{0.44\pm
0.01},\text{Ra}^{0.45\pm 0.01},\text{Ra}^{0.44}\\}$ and
$\text{Re}_{\perp}\propto\\{\text{Ra}^{0.22\pm 0.01},\\\ \text{Ra}^{0.23\pm
0.01},\text{Ra}^{0.21}\\}$ for $\text{Ro}_{\text{p}}=\\{0.6,0.957,1.58\\}$.
These scalings are similar to but slightly weaker than our predictions in all
cases. However, the scaling of $\text{Re}_{\parallel}\propto\text{Ra}^{0.45}$,
is once again a power law observed frequently in nonrotating convection
(Ahlers et al., 2009, AB17). We also observe that $\text{Re}_{\perp}$
collapses for each $\text{Ro}_{\text{p}}\,$ track, while
$\text{Re}_{\parallel}$ experiences an offset to larger values as
$\text{Ro}_{\text{p}}$ shrinks. The offset in $\text{Re}_{\parallel}$ is
unsurprising, because more rotationally constrained flows result in smaller
boundary layers relative to the vertical extent of our stratified domain. The
horizontal extent of our domain scales with the strength of rotation, and so
regardless of $\text{Ro}_{\text{p}}$, flows perpendicular to the rotational
and buoyant direction are comparably turbulent at the same
$\text{Ra}/\text{Ra}_{\text{crit}}$. We find $\text{Re}_{\perp}$ and
$\text{Re}_{\parallel}$ are, respectively, good proxies for the horizontal and
perpendicular resolution required to resolve an experiment.

Figure 4 shows time- and horizontally-averaged profiles of Ro and the standard
deviation of the entropy, $\sigma_{\text{s}}$. Figures 4a&b show these
profiles for $\text{Ro}_{\text{p}}=1.58$ ($\text{Ro}\approx 0.4$), while
Figures 4c&d show these profiles for $\text{Ro}_{\text{p}}=0.96$
($\text{Ro}\approx 0.1$). The transition in profile behavior from low Ra
(yellow) to high Ra (purple) is denoted by the color of the profile. As Ra
increases at a constant value of $\text{Ro}_{\text{p}}$, both the thermal
($\sigma_{\text{s}}$) and dynamical (Ro) boundary layers become thinner. We
measure the thickness of the thermal boundary layer ($\delta_{\text{s}}$) at
the top of the domain by finding the location of the first maxima of
$\sigma_{\text{s}}$ away from the boundary. We measure the thickness of the Ro
boundary layer ($\delta_{\text{Ro}}$) in the same manner. In figure 4e, we
plot $\delta_{\text{Ro}}/\delta_{\text{s}}$, the ratio of the sizes of these
two boundary layers. As anticipated, the dynamical boundary layer
($\delta_{\text{Ro}}$) becomes relatively thinner with respect to the thermal
boundary layer ($\delta_{\text{s}}$) as Ro and $\text{Ro}_{\text{p}}$
decrease. However, the precise scaling of this boundary layer ratio with
$\text{Ro}_{\text{p}}$ and Ra is unclear, and we cannot immediately compare
these ratios to similar measures from the Rayleigh-Bénard convection
literature, such as Fig. 5 of King et al. (2013). They measure the dynamical
boundary layer thickness as the peak location of the horizontal velocities,
but our horizontal velocities are subject to stress-free boundary conditions,
and we find that the maxima of horizontal velocities occur precisely at the
boundaries. In figure 4f, we plot $\delta_{\text{s}}$ in units of the density
scale height at the top of the atmosphere, and we plot vertical lines when
this crosses 1. We find no systematic change in behavior when
$\delta_{\text{s}}$ is smaller than the local density scale height.

##  4 Discussion

In this letter, we studied low-Mach-number, stratified, compressible
convection under the influence of rotation. We examined three paths through
Ra-Ta space, and showed that the newly-defined Predictive Rossby number,
$\text{Ro}_{\text{p}}=\text{Ra}/(\text{Pr }^{1/2}\text{Ta}^{3/4})$, determines
the value of the evolved Rossby number.

Shockingly, along these constant $\text{Ro}_{\text{p}}$ pathways, particularly
when $\text{Ro}\lesssim 0.1$, we find $\text{Nu}\propto\text{Ra}^{2/7}$ and
$\text{Re}_{\parallel}\propto\text{Ra}^{0.45}$. These scalings are
indistinguishable from the scalings of Re and Nu with Ra in non-rotating
Boussinesq convection (Ahlers et al., 2009). Julien et al. (2012) theorized
that in the rapidly rotating asymptotic limit,
$(\text{Nu}-1)\propto(\text{Ra}^{3/2}/\text{Ta})=(\text{Ra}/\text{Ra}_{\text{crit}}(\text{Ta}))^{3/2}$.
Thus, at fixed Ta, a very sharp $\text{Ra}^{3/2}$ scaling law is expected. At
a fixed Ta = $10^{14}$, Stellmach et al. (2014) found that the Ra3/2 scaling
described the results of stress-free DNS in Boussinesq cylinders very well.
Gastine et al. (2016) studied Boussinesq convection in spherical shells with
no-slip boundaries, and also found good agreement with the theory of Julien et
al. (2012) for various Ra at $\text{Ta}\geq 10^{10}$.

Here, when we run simulations at fixed $\text{Ro}_{\text{p}}$, the value of Ta
is coupled to the value of Ra, and both increase simultaneously. Recasting the
scaling of Julien et al. (2012) into this perspective, we find
$(\text{Nu}-1)\propto\text{Ra}^{3/2}/\text{Ta}=\text{Ro}_{\text{p}}^{8/3}\text{Ra}^{1/6}\propto(\text{Ra}/\text{Ra}_{\text{crit}})^{1/6}$,
where in this final result we use the $\text{Ra}_{\text{crit}}$ value of the
whole $\text{Ro}_{\text{p}}$ path, such as those specified in Table 1. This
$\text{Ra}^{1/6}$ scaling is much weaker than the $\text{Ra}^{2/7}$ law we
find here. We leave it to future work to explain this discrepancy between
Boussinesq theory and our observed Nu vs. Ra scaling.

In this work, we experimentally arrived at the
$\text{Ra}/\text{Ta}^{3/4}=\text{Ra}\,\text{Ek}^{3/2}$ scaling in
$\text{Ro}_{\text{p}}$, but this relationship was independently discovered by
King et al. (2012). Arguing that the thermal boundary layers should scale as
$\delta_{S}\propto\text{Ra}^{-1/3}$ and rotational Ekman boundary layers
should scale as $\delta_{\text{Ro}}\propto\text{Ta}^{-1/4}=\text{Ek}^{1/2}$,
they expect these boundary layers to be equal in size when
$\text{Ra}/\text{Ta}^{3/4}\sim 1$. They demonstrate that when
$2\lesssim\text{Ra}/\text{Ta}^{3/4}\lesssim 20$ flows are in the transitional
regime, and for $\text{Ra}/\text{Ta}^{3/4}\lesssim 2$, flows are rotationally
constrained. We remind the reader that Boussinesq values of Ra and Ta are not
the same as their values in our stratified domains here, as diffusivities
change with depth (see section 2). Taking into account this change with depth,
our simulations fall in King et al. (2012)’s rotationally constrained
($\text{Ro}_{\text{p}}=0.6$) and near-constrained transitional regime
($\text{Ro}_{\text{p}}=\\{0.957,1.58\\}$). The measured values of Ro in Fig.
1b and the observed dynamics in Fig. 2 agree with this interpretation.

We note briefly that the scaling $\text{Ra}\propto\text{Ta}^{3/4}$ is very
similar to another theorized boundary between fully rotationally constrained
convection and partially constrained convection predicted in Boussinesq
theory, of $\text{Ra}\propto\text{Ta}^{4/5}$ (Julien et al., 2012; Gastine et
al., 2016). This Ta4/5 scaling also arises through arguments of geostrophic
balance in the boundary layers, and is a steeper scaling than the Ta3/4
scaling present in $\text{Ro}_{\text{p}}$. This suggests that at sufficiently
low $\text{Ro}_{\text{p}}$, a suite of simulations across many orders of
magnitude of Ra will not only have the same volume-averaged value of Ro (as in
Fig. 1b), but will also maintain proper force balances within the boundary
layers.

Our results suggest that by choosing the desired value of
$\text{Ro}_{\text{p}}$, experimenters can select the degree of rotational
constraint present in their simulations. We find that
$\text{Ro}\propto\text{Ro}_{\text{p}}^{3.34\pm 0.07}$, which is within
2$\sigma$ of the estimate in King et al. (2013), who although defining Ro very
differently from our vorticity-based definition here, find
$\text{Ro}\propto\text{Ro}_{\text{p}}^{3.84\pm 0.28}$. We note briefly that
they claim that the value of Ro is strongly dependent upon the Prandtl number
studied, and that low Ro can be achieved at high Pr without achieving a
rotationally constrained flow. We studied only $\text{Pr}=1$ here, and leave
it to future work to determine if the scaling of
$\text{Ro}_{\text{p}}\propto\text{Pr}^{-1/4}$ is the correct scaling to
predict the evolved Rossby number.

Despite the added complexity of stratification and despite our using stress-
free rather than no-slip boundaries, the boundary layer scaling arguments put
forth in King et al. (2012) seem to hold up in our systems. This is
reminiscent of what we found in AB17, in which convection in stratified
domains, regardless of Mach number, produced boundary-layer dominated scaling
laws of Nu that were nearly identical to the scaling laws found in Boussinesq
Rayleigh-Bénard convection.

We close by noting that once $\text{Ro}_{\text{p}}$ is chosen such that a
convective system has the same Rossby number as an astrophysical object of
choice, it is straightforward to increase the turbulent nature of simulations
by increasing Ra, just as in the non-rotating case. Although all the results
reported here are for a Cartesian geometry with antiparallel gravity and
rotation, preliminary 3D spherical simulations suggest that
$\text{Ro}_{\text{p}}$ also specifies Ro in more complex geometries (Brown et
al. 2019 in prep).

We thank Jon Aurnou and the anonymous referee who independently directed us to
the work of King et al. (2012) during the review process, which greatly helped
us understand the results of the experiments here. We further thank the
anonymous referee for other helpful and clarifying suggestions. This work was
supported by NASA Headquarters under the NASA Earth and Space Science
Fellowship Program – Grant 80NSSC18K1199. EHA further acknowledges the
University of Colorado’s George Ellery Hale Graduate Student Fellowship. This
work was additionally supported by NASA LWS grant number NNX16AC92G.
Computations were conducted with support by the NASA High End Computing (HEC)
Program through the NASA Advanced Supercomputing (NAS) Division at Ames
Research Center on Pleiades with allocation GID s1647.

##  Appendix A Table of Simulations

Information for select simulations in this work are shown in Table 2. The
simulation at minimum (Ra, Ta) and maximum (Ra, Ta) for each of the
$\text{Ro}_{\text{p}}$, $\text{Ro}_{\text{c}}$, and $\mathcal{S}$ paths in
Fig. 1b are shown. This information for the displayed simulations and all
other simulations in this work is included as a .csv file in the supplemental
materials and is published online in a Zenodo repository (Anders et al.,
2019).

Table 2: Table of simulation information

Ra${}_{\text{top}}$ | Ta${}_{\text{top}}$ | Ro${}_{\text{p, top}}$ | Ra${}_{\text{mid}}$ | Ta${}_{\text{mid}}$ | Ro${}_{\text{p, mid}}$ | $\text{Ro}_{\text{c}}$ | $\mathcal{S}$ | $L_{x}/L_{z}$ | (nz, nx, ny) | Ro | Re∥ | Re⟂ | Nu
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Constant $\text{Ro}_{\text{p}}$ , path III |  |  |  |  |  |  |  |  |  |  |
1.8 $\times 10^{5}$ | 4.1 $\times 10^{7}$ | 0.60 | 1.4 $\times 10^{7}$ | 3.0 $\times 10^{9}$ | 1.03 | 0.067 | 1.1 | 0.51 | (256, 32, 32) | 0.015 | 19.4 | 2.5 | 1.2
1.2 $\times 10^{9}$ | 5.2 $\times 10^{12}$ | 0.60 | 9.2 $\times 10^{10}$ | 3.8 $\times 10^{14}$ | 1.03 | 0.015 | 3.0 | 0.07 | (2048, 64, 64) | 0.026 | 1771 | 32.0 | 15.4
5.2 $\times 10^{2}$ | 4.6 $\times 10^{3}$ | 0.96 | 3.8 $\times 10^{4}$ | 3.4 $\times 10^{5}$ | 1.64 | 0.333 | 1.13 | 2.28 | (64, 64, 64) | 0.074 | 4.2 | 2.5 | 1.2
3.8 $\times 10^{8}$ | 3.1 $\times 10^{11}$ | 0.96 | 2.8 $\times 10^{10}$ | 2.3 $\times 10^{13}$ | 1.64 | 0.035 | 6.0 | 0.12 | (2048, 64, 64) | 0.129 | 3906 | 113 | 66.9
7.9 $\times 10^{1}$ | 1.0 $\times 10^{2}$ | 1.58 | 5.8 $\times 10^{3}$ | 7.4 $\times 10^{3}$ | 2.70 | 0.888 | 1.56 | 4.44 | (64, 64, 64) | 0.303 | 4.4 | 4.9 | 1.7
1.4 $\times 10^{7}$ | 9.7 $\times 10^{8}$ | 1.58 | 1.0 $\times 10^{9}$ | 7.2 $\times 10^{10}$ | 2.70 | 0.119 | 10.0 | 0.30 | (512, 128, 128) | 0.376 | 1257 | 94.9 | 40.1
Constant $\text{Ro}_{\text{c}}$ , path II |  |  |  |  |  |  |  |  |  |  |
8.6 $\times 10^{4}$ | 8.6 $\times 10^{6}$ | 0.74 | 6.3 $\times 10^{6}$ | 6.3 $\times 10^{8}$ | 1.26 | 0.1 | 1.47 | 0.68 | (128, 128, 128) | 0.051 | 40.2 | 6.7 | 2.2
2.6 $\times 10^{6}$ | 2.6 $\times 10^{8}$ | 1.13 | 1.9 $\times 10^{8}$ | 1.9 $\times 10^{10}$ | 1.93 | 0.1 | 4.64 | 0.39 | (256, 512, 512) | 0.27 | 565 | 53.2 | 33.3
1.4 $\times 10^{3}$ | 1.6 $\times 10^{4}$ | 1.01 | 1.1 $\times 10^{5}$ | 1.2 $\times 10^{6}$ | 1.72 | 0.3 | 1.47 | 1.90 | (64, 128, 128) | 0.124 | 11.3 | 5.4 | 1.8
1.1 $\times 10^{6}$ | 1.2 $\times 10^{7}$ | 2.29 | 7.8 $\times 10^{7}$ | 8.6 $\times 10^{8}$ | 3.93 | 0.3 | 14.7 | 0.65 | (192, 384, 384) | 0.808 | 529 | 83.5 | 27.3
5.5 $\times 10^{1}$ | 5.5 $\times 10^{1}$ | 1.65 | 4.0 $\times 10^{3}$ | 4.0 $\times 10^{3}$ | 2.82 | 1.0 | 1.47 | 4.84 | (64, 128, 128) | 0.303 | 3.6 | 4.4 | 1.5
2.8 $\times 10^{6}$ | 2.8 $\times 10^{6}$ | 6.39 | 2.0 $\times 10^{8}$ | 2.0 $\times 10^{8}$ | 10.93 | 1.0 | 100 | 0.82 | (256, 512, 512) | 3.357 | 1099 | 220 | 46.6
Constant $\mathcal{S}$, path I |  |  |  |  |  |  |  |  |  |  |
1.9 $\times 10^{1}$ | 1.1 $\times 10^{-1}$ | 10.00 | 1.4 $\times 10^{3}$ | 7.9 $\times 10$ | 17.12 | 13.2 | 2.0 | 9.91 | (64, 64, 64) | 3.668 | 3.0 | 10.2 | 1.7
3.0 $\times 10^{6}$ | 1.1 $\times 10^{9}$ | 0.70 | 2.2 $\times 10^{8}$ | 8.1 $\times 10^{10}$ | 1.20 | 0.052 | 2.0 | 0.30 | (512, 128, 128) | 0.053 | 242 | 17.9 | 6.0
3.0 $\times 10^{1}$ | 2.0 $\times 10^{-1}$ | 10.00 | 2.2 $\times 10^{3}$ | 1.5 $\times 10^{1}$ | 17.12 | 12.2 | 3.0 | 9.48 | (64, 64, 64) | 4.418 | 5.0 | 15.5 | 2.1
1.3 $\times 10^{7}$ | 5.6 $\times 10^{9}$ | 0.80 | 9.7 $\times 10^{8}$ | 4.2 $\times 10^{11}$ | 1.37 | 0.048 | 3.0 | 0.23 | (512, 128, 128) | 0.08 | 592 | 33.4 | 13.6

Note. — Input parameters and output parameters for select simulations are
shown. For each of the eight paths in Fig. 1b, we show information for the
lowest and highest (Ra, Ta) point on that path. The first six rows show
information for constant $\text{Ro}_{\text{p}}\,$ paths, the next six for
constant $\text{Ro}_{\text{c}}$ paths, and the last four for constant
$\mathcal{S}$ paths. We show the input Ra, Ta, and $\text{Ro}_{\text{p}}$ at
the top of the atmosphere, as well as their stratification-weighted values at
the midplane of the atmosphere, which provide a more direct comparison to
Boussinesq values (Unno et al., 1960). We also provide the input
$\text{Ro}_{\text{c}}$ at the top of the atmosphere, $\mathcal{S}$, aspect
ratio ($L_{x}/L_{z}$), and coefficient resolution (nz, nx, ny). Each dimension
of the physical grid is 3/2 the size of the coefficient grid for adequate
dealiasing of quadratic nonlinear terms. Output values of Ro, Re∥, Re⟂, and Nu
are also provided. This table in its entirety is published as a supplemental
.csv file with this manuscript and also online in a Zenodo repository (Anders
et al., 2019).

## References

  * Ahlers et al. (2009) Ahlers, G., Grossmann, S., & Lohse, D. 2009, Rev. Mod. Phys., 81, 503
  * Anders & Brown (2017) Anders, E. H., & Brown, B. P. 2017, Physical Review Fluids, 2, 083501
  * Anders et al. (2019) Anders, E. H., Manduca, C. M., Brown, B. P., Oishi, J. S., & Vasil, G. M. 2019, Supplemental Materials: Predicting the Rossby number in convective experiments, v.1.0, Zenodo, doi:10.5281/zenodo.2539500
  * Augustson et al. (2012) Augustson, K. C., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2012, ApJ, 756, 169
  * Aurnou & King (2017) Aurnou, J. M., & King, E. M. 2017, Proceedings of the Royal Society of London Series A, 473, 20160731
  * Brown et al. (2008) Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2008, ApJ, 689, 1354
  * Brown et al. (2010) —. 2010, ApJ, 711, 424
  * Brown et al. (2011) Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69
  * Brummell et al. (1996) Brummell, N. H., Hurlburt, N. E., & Toomre, J. 1996, ApJ, 473, 494
  * Brummell et al. (1998) —. 1998, ApJ, 493, 955
  * Brun et al. (2017) Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192
  * Busse (2002) Busse, F. H. 2002, Physics of Fluids, 14, 1301
  * Calkins et al. (2015) Calkins, M. A., Julien, K., & Marti, P. 2015, Geophysical and Astrophysical Fluid Dynamics, 109, 422
  * Chandrasekhar (1961) Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability
  * Cheng et al. (2015) Cheng, J. S., Stellmach, S., Ribeiro, A., et al. 2015, Geophysical Journal International, 201, 1
  * Ecke & Niemela (2014) Ecke, R. E., & Niemela, J. J. 2014, Phys. Rev. Lett., 113, 114301
  * Featherstone & Hindman (2016) Featherstone, N. A., & Hindman, B. W. 2016, ApJ, 830, L15
  * Gastine et al. (2016) Gastine, T., Wicht, J., & Aubert, J. 2016, Journal of Fluid Mechanics, 808, 690
  * Gastine et al. (2014) Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76
  * Glatzmaier & Gilman (1982) Glatzmaier, G. A., & Gilman, P. A. 1982, ApJ, 256, 316
  * Greer et al. (2016) Greer, B. J., Hindman, B. W., & Toomre, J. 2016, ApJ, 824, 4
  * Guerrero et al. (2013) Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176
  * Julien et al. (2012) Julien, K., Knobloch, E., Rubio, A. M., & Vasil, G. M. 2012, Physical Review Letters, 109, 254503
  * Julien et al. (1996) Julien, K., Legg, S., McWilliams, J., & Werne, J. 1996, Journal of Fluid Mechanics, 322, 243
  * Käpylä et al. (2014) Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43
  * King et al. (2012) King, E. M., Stellmach, S., & Aurnou, J. M. 2012, Journal of Fluid Mechanics, 691, 568
  * King et al. (2013) King, E. M., Stellmach, S., & Buffett, B. 2013, Journal of Fluid Mechanics, 717, 449
  * King et al. (2009) King, E. M., Stellmach, S., Noir, J., Hansen, U., & Aurnou, J. M. 2009, Nature, 457, 301
  * Schmitz & Tilgner (2009) Schmitz, S., & Tilgner, A. 2009, Phys. Rev. E, 80, 015305
  * Soderlund et al. (2015) Soderlund, K. M., Sheyko, A., King, E. M., & Aurnou, J. M. 2015, Progress in Earth and Planetary Science, 2, 24
  * Stellmach et al. (2014) Stellmach, S., Lischper, M., Julien, K., et al. 2014, Phys. Rev. Lett., 113, 254501
  * Unno et al. (1960) Unno, W., Kato, S., & Makita, M. 1960, Publications of the Astronomical Society of Japan, 12, 192
  * Zhong et al. (2009) Zhong, J.-Q., Stevens, R. J. A. M., Clercx, H. J. H., et al. 2009, Physical Review Letters, 102, 044502


$B_{j}$ for $i\neq j$. To eliminate this problem, we will add an explicit
identifier for the index $j$ to the functor. To view the index $i$ as an
element of the input structure $A_{i}$, we use the first element in the
universe of the corresponding list element $A_{i}$. Here, when we refer to the
first element in the universe, we mean the natural linear order on the
universe in a structure from a graded list type, which arises from the ordered
nature of lists and pairs. Therefore, instead of $H_{\ell}(A)$, to represent
item 2* we use the product $A\times H_{\ell}(A)$, with the $A$ part
representing the index $i$. Summing up, the functor $G_{\ell}$ that describes
elements in each $B_{i}$ is

$\displaystyle
G_{\ell}(A)=A+\underbrace{1+\cdots+1}_{\mathclap{\text{\scriptsize\begin{tabular}[]{c}number
of elements\\\ in $B_{0}$ that have grade $\ell$\end{tabular} }}}+A\times
H_{\ell}(A).$

In the rest of this proof, we will view the universe of $B_{n}$ as being a
subset of

$\displaystyle G(A)=G_{0}(A)+\cdots+G_{m}(A),$

with $G_{\ell}(A)$ representing the elements of grade $\ell$. The polynomial
functor $G(A)$ will be the polynomial functor for the mso interpretation in
the conclusion of the lemma. To conclude the proof of the lemma, we need to
show that in mso we can define which elements of $G(A)$ belong to the universe
of $B_{n}$, and what relations from the output vocabulary are satisfied by
tuples of such elements. In other words, we need to define in mso the
quantifier-free theory of tuples from $G(A)$ in the output structure. This is
done in the following claim, which completes the proof of the lemma.

###### Claim C.3.

For every $\ell,k\in\\{0,1,\ldots\\}$ the following function is mso definable:

  * •

Input. A structure $A\in\Sigma^{*}$ with elements $\bar{a}\in A^{k}$.

  * •

Output. The quantifier-free theory of $G(\bar{a})$ in $B_{n}|\ell$.

Furthermore, the output depends only on $A$ and $\bar{a}$ restricted to
elements of grade at least $\ell$.

Proof
Fix some $\ell$ and $k$ as in the statement of the claim. The claim is proved
by induction on $\ell$, in reverse order $m,\ldots,0$. Suppose that we want to
prove the claim for some grade $\ell$, and assume that it has already been
proved for strictly bigger grades.

We use the same idea as in the proof of Claim 3.4. Consider a finite
automaton, in which the states are all possible theories that arise by taking
some $k$-tuple $\bar{a}$, and returning the quantifier-free theory of
$G(\bar{a})$ in some structure from $\Gamma$. This set of states is finite,
since the length of the tuple and the vocabulary are fixed.

We will design an automaton with this set of states, together with an input
string (which will be called the _advice string_), so that it satisfies the
following invariant: after reading the first $i$ letters of the advice string,
the state of the automaton is the quantifier-free theory of $G(\bar{a})$ in
$B_{i}|\ell$.

The initial state of the automaton is determined by the invariant, it must be
the quantifier-free theory of $G(\bar{a})$ in $B_{0}$. Since the universe of
$B_{0}$ is equal to $G(\emptyset)$, it follows that the initial state does not
depend on the tuple $\bar{a}$ or the input structure $A$.

We now describe the transition function of the automaton, as well as the
advice string. By unfolding the definition of the graded mso interpretation
$\delta$, there is some quantifier rank $s$ such that the state of the
automaton after reading $i$ letters is uniquely determined by the following
four pieces of information:

  1. (1)

the quantifier-free theory of $G(\bar{a})$ in $B_{i-1}$,

  2. (2)

the quantifier-free theory of $G(\bar{a})$ in $A_{i}$,

  3. (3)

the rank $s$ mso theory of $G(\bar{a})$ in $B_{i-1}|\ell+1$,

  4. (4)

the rank $s$ mso theory of $G(\bar{a})$ in $A_{i}|\ell+1$.

The first piece of information is the previous state of the automaton. The
remaining infomration will be the stored in the advice string; i.e. the $i$-th
letter of the advice string will contain the information described the last
three items above. Note that the advice string can be computed in mso, by the
induction assumption. Therefore, since the automaton can be simulated in mso,
it follows that the last state of this automaton can be defined in mso, thus
proving the claim. $\Box$ $\Box$

##  Appendix D Proof of Theorem 6.1

In this section, we prove that the system in Theorem 5.3 is sound and complete
with respect to linear regular functions.

##### Soundness.

The soundness proof follows the same lines as the soundness proof in Theorem
5.3. The general idea is that we use graded mso interpretations where all
components have dimension at most one. This, however, on its own is not going
to be enough. To see why, let us compare the two absorption functions

$\displaystyle\underbrace{!\Sigma\to\Sigma\times!\Sigma}_{\mathclap{\text{\scriptsize\begin{tabular}[]{c}not
allowed\end{tabular}
}}}\quad\underbrace{!\Sigma\to\Sigma\times\Sigma}_{\mathclap{\text{\scriptsize\begin{tabular}[]{c}allowed\end{tabular}
}}}.$

Both of them have linear size increase – each element of the input structure
contributes two copies to the output structure. What is wrong with the
function that is not allowed? The problem is that one of the copies has the
same grade, and the other has lower grade. In the presence of folding, we can
get an unbounded number of copies, by spawning a new lower grade copy in each
iteration. This phenomenon will not occur in the allowed function, since both
copies have lower grade. The phenomenon discussed above is formalised in the
following definition:

###### Definition D.1.

A _linear graded mso interpretation_ is a graded mso interpretation in which
the underlying functor is linear, i.e. all components have dimension one, and
which furthermore satisfies the following _downgrading_ condition: if an
element of the input structure has at least two copies in the output
structure, then all of the copies have strictly lower grade.

In the definition above, the copies of an element in the output structure are
defined in the natural way; this definition makes sense when the functor is
linear. For example, if the functor is

$\displaystyle A+A+A+1+1$

then each input element spawns at most three copies. The components of
dimension zero, of which there are two in the above example, are not counted
as copies of any input elment.

To prove completeness of the system from Theorem 6.1, we show that all
functions that are strongly derived in it are linear graded mso
interpretations. The proof is a simple inducton on the derivation. The most
interesting cases are composition and folding. For composition, we simply
observe that the condition on lower grades from Definition D.1 is preserved
under composition.

We are left with folding. where we use the following lemma, which is the same
as Lemma C.2 except that the functions in the assumption and conclusion are
required to be linear. In the assumption, we use linearity as defined in
Definition D.1, in particular the downgrading condition is assumed; in the
conclusion we have an ungraded function, and therefore only the linearity of
the functor and not the downgrading condition are assumed.

###### Lemma D.2.

For every linear graded mso interpretation

$\displaystyle\delta:\Gamma\times\Sigma\to\Gamma,$

between graded list types, and every $B_{0}\in\Gamma$, the following function
is an (ungraded) linear mso interpretation

$\displaystyle
A=\underbrace{[A_{1},\ldots,A_{n}]}_{\mathclap{\text{\scriptsize\begin{tabular}[]{c}list
of structures in $\Sigma$,\\\ with the grades forgotten\end{tabular}
}}}\qquad\mapsto\qquad\underbrace{B_{n}}_{\mathclap{\text{\scriptsize\begin{tabular}[]{c}defined
based on $A$\\\ as in the proof of Calim~{}\ref{claim:qf-folder}\end{tabular}
}}}.$

Proof
We use the same proof as in Lemma C.2. However, there is one difficulty, which
is that the functor $G$ defined in that proof is not linear, even if $\delta$
is linear. This is because of the product $A\times H_{\ell}(A)$ which is used
to code indexes. In fact, the functor $G$ can have arbitrarily high dimension.
However, thanks to the downgrading condition on $\delta$, one show by
induction that for every grade $\ell$ there is some constant
$c_{\ell}\in\\{0,1,\ldots\\}$ such that for every grade $\ell$ element $a$ in
the input structure, there are at most $c_{\ell}$ elements in the output
structure which use $a$. Here, we say that an element uses $a$ if it belongs
to $G(A)$ but not to $G(A\setminus\\{a\\})$. Using this property, we can turn
$G$ into a linear functor. $\Box$

This finishes the soundness proof. Below, we give two completeness proofs.

Proof (First completeness proof.)
This proof uses the sst model from Example 3.4, which is complete for linear
regular functions, in the case where the input and output types are strings
over finite alphabets [1, Theorem 3]. In Example 3.4, we show how to weakly
derive every sst that uses each input letter at most once. To get the general
form of sst, where an input letter can be used a constant number of times, it
is enough to generalize the model from Example 3.4 so that the initial
function is weakly derivable, and the transition function can be derived in
type

$\displaystyle\Delta\times!^{k}\Sigma\to\Delta$

for some $k$. With these relaxations, we get all copyless sst, and retain weak
derivability. This proof works only for functions of string-to-string type
(admittedly, this is the case that we really care about), and for this reason
we also present a second proof, which can also handle types such as strings of
strings or pairs of strings. $\Box$

Proof (Second completeness proof.)
In this proof, similarly to the completeness proof from Theorem 5.3, we reduce
to a known complete system. In the case of linear mso interpretations, the
corresponding known system is from [7]. It is the same as in Theorem B.1,
except that the split function is removed. In the completeness proof of
Theorem 5.3, only the proof for split used general absorption (as opposed to
linear absorption). Therefore, the system with linear absorption is complete
for the linear regular functions. $\Box$


# Data-Efficient Reinforcement Learning for Malaria Control

Lixin Zou1 Long Xia1 Linfang Hou2 Xiangyu Zhao3&Dawei Yin1 1Baidu Inc., 2
JD.com
3Michigan State University
{zoulixin15, long.phil.xia<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Sequential decision-making under cost-sensitive tasks is prohibitively
daunting, especially for the problem that has a significant impact on people’s
daily lives, such as malaria control, treatment recommendation. The main
challenge faced by policymakers is to learn a policy from scratch by
interacting with a complex environment in a few trials. This work introduces a
practical, data-efficient policy learning method, named Variance-Bonus Monte
Carlo Tree Search (VB-MCTS), which can copy with very little data and
facilitate learning from scratch in only a few trials. Specifically, the
solution is a model-based reinforcement learning method. To avoid model bias,
we apply Gaussian Process (GP) regression to estimate the transitions
explicitly. With the GP world model, we propose a variance-bonus reward to
measure the uncertainty about the world. Adding the reward to the planning
with MCTS can result in more efficient and effective exploration. Furthermore,
the derived polynomial sample complexity indicates that VB-MCTS is sample
efficient. Finally, outstanding performance on a competitive world-level RL
competition and extensive experimental results verify its advantage over the
state-of-the-art on the challenging malaria control task.

##  1 Introduction

Malaria is a mosquito-borne disease that continues to pose a heavy burden on
South Sahara Africa (SSA) Moran (2007). Recently, there has been significant
progress in improving treatment efficiency and reducing the mortality rate of
malaria. Unfortunately, due to financial constraints, the policymakers face
the challenge of ensuring continued success in disease control with
insufficient resources. To make intelligent decisions, learning control
policies over the years have been formulated as Reinforcement Learning (RL)
problems Bent et al. (2018). Nevertheless, applying RL to malaria control
seems to be a tricky issue since RL usually requires numerous trial-and-error
searches to learn from scratch. Unlike simulation-based games, _e.g.,_ Atari
games Mnih et al. (2013) and game GO Silver et al. (2017), the endless
intervention trial is unacceptable to regions over the years since the actual
cost of life and money is enormous. Hence, as in many human-in-loop systems
Zou et al. (2019, 2020b, 2020a), it is too expensive to apply RL to learn
malaria intervention policy from scratch directly.

Therefore, to reduce the heavy burden of malaria in SSA, it is urgent to
improve the data efficiency of policy learning. In Bent et al. (2018), novel
exploration techniques, such as Genetic Algorithm Holland (1992), Batch Policy
Gradient Sutton et al. (2000) and Upper/Lower Confidence Bound Auer and Ortner
(2010), have been firstly applied to learn malaria control policies from
scratch. However, these solutions are introduced under the Stochastic Multi-
Armed Bandit (SMAB) setting, which myopically ignores the delayed impact of
the interventions in the future and might result in serious problems. For
example, the large-scale use of spraying may lead to mosquito resistance and
bring about the uncontrolled spread of malaria in the coming years. Hence, it
requires us to optimize disease control policies in the long run, which is far
more challenging than the exploration in SMAB.

Considering the long-term effects, the finite horizon continuous-space Markov
Decision Process is employed to model the disease control in this work. Under
this setting, we propose a framework named Variance-Bonus Monte Carlo Tree
Search (VB-MCTS) for data-efficient policy searching, illustrated in Figure 1.
Particularly, it is a model-based training framework, which iterates between
updating the world model and collecting data. In model training, Gaussian
Process (GP) is used to approximate the state transition function with
collected rollouts. As a non-parametric probabilistic model, GP can avoid the
model bias and explicitly model the uncertainty about the transitions, _i.e.,_
the variance of the state. In data collection, we propose to employ MCTS for
generating the policy with the mean MDP plus variance-bonus reward. The
variance-bonus reward can decrease the uncertainty at the state-action pairs
with high potential reward by explicitly motivating the agent to sample the
state-actions with the highest upper-bounded reward. Furthermore, the sample
complexity of the proposed method indicates that it is a PAC optimal
exploration solution for malaria control. Finally, to verify the effectiveness
of our policy search solution, extensive experiments are conducted on the
malaria control simulators111https://github.com/IBM/ushiriki-policy-engine-
library Bent et al. (2018), which are gym-like222https://gym.openai.com
environments for KDD Cup 2019 Zhou et al. (2020). The outstanding performance
on the competition and extensive experimental results demonstrated that our
approach could achieve unprecedented data efficiency on malaria control
compared to the state-of-the-art methods.

Figure 1: Overview of the proposed learning system on Malaria Control. The
system alternates model training and data collection.

Our main contributions are: (1) We propose a highly data-efficient learning
framework for malaria control. Under the framework, the policymakers can
successfully learn control policies from scratch within 20 rollouts. (2) We
derive the sample complexity of the proposed method and verify that VB-MCTS is
an efficient PAC-MDP algorithm. (3) Extensive experiments conducted on malaria
control demonstrate that our solution can outperform the state-of-the-art
methods.

##  2 Related Work

As a highly pathogenic disease, malaria has been widely studied from the
perspective of predicting the disease spread, diagnosis, and personalized care
planning. However, rarely work focuses on applying RL to learn the cost-
effectiveness intervention strategies, which plays a crucial role in
controlling the spread of malaria Moran (2007). In Bent et al. (2018), malaria
control has firstly been formulated as a stochastic multi-armed bandit (SMAB)
problem and solved with novel exploration techniques. Nevertheless, SMAB based
solutions only myopically maximize instant rewards, and the ignorance of
delayed influences might result in disease outbreaks in the future. Therefore,
in this work, a comprehensive solution has been proposed to facilitate policy
learning in a few trials under the setting of finite-horizon MDP.

Another topic is data-efficient RL. To increase the data efficiency, we are
required to extract more information from available trials Deisenroth and
Rasmussen (2011), which involves utilizing the samples in the most efficient
way (_e.g.,_ exploitation) and choosing the samples with more information
(_e.g.,_ exploration). Generally, for exploitation, model-based methods Ha and
Schmidhuber (2018); Kamthe and Deisenroth (2017) are more sample efficient but
require more computation time for the planning. Model-free Szita and Lörincz
(2006); Krause et al. (2016); Van Seijen et al. (2009) methods are generally
computationally light and can be applied without a planner, but need
(sometimes exponentially) more samples, and are usually not efficient PAC-MDP
algorithm Strehl et al. (2009). For exploration, there are two options: (1)
Bayesian approaches, considering a distribution over possible models and
acting to maximize expected reward; unluckily, it is intractable for all but
very restricted cases, such as the linear policy assumption in PILCO
Deisenroth and Rasmussen (2011). (2) intrinsically motivated exploration,
implicitly negotiating the exploration/exploitation dilemma by always
exploiting a modified reward for directly accomplishing exploration. However,
on the one hand, the vast majority of papers only address the discrete state
case, providing incremental improvements on the complexity bounds, such as
MMDP-RB Sorg et al. (2012), metric-E3 Kakade et al. (2003), and BED Kolter and
Ng (2009). On the other hand, for more realistic continuous state space MDP,
over-exploration has been introduced for achieving polynomial sample
complexity in many work, such as KWIK Li et al. (2011), and GP-Rmax Grande et
al. (2014). These methods will explore all regions equally until the reward
function is highly accurate everywhere. By drawing on the strength of existing
methods, our solution is a model-based RL framework, which efficiently plans
with MCTS and trade-off exploitation and exploration by exploiting a variance-
bonus reward.

##  3 Proposed Method: VB-MCTS

###  3.1 Malaria Control as MDP

Finding an optimal malaria control policy can be posted as a reinforcement
learning task by illustrating it as a Markov Decision Process (MDP).
Specifically, we formulate the task as a finite-horizon MDP, defined by the
tuple $\langle\mathcal{S},\mathcal{A},P,R,\gamma\rangle$ with $\mathcal{S}$ as
the potential infinite state space, $\mathcal{A}$ as the finite set of
actions, $P:\mathcal{S}\times\mathcal{A}\rightarrow\mathcal{S}$ as the
deterministic transition function, $R:\mathcal{S}\rightarrow\mathbb{R}$ as the
reward function, and $\gamma\in(0,1]$ as the discount factor. In this case, we
face the challenge of developing an efficient policy for a population over a 5
year intervention time frame. As shown in Figure 2, the corresponding
components in malaria control are defined as,

Figure 2: Malaria intervention as MDP.

#### Action

The actions are the available means of interventions, including the mass-
distribution of long-lasting insecticide-treated nets (ITNs) and indoor
residual spraying (IRS) with pyrethroids in SSA Stuckey et al. (2014). In this
work, the action space $\mathcal{A}$ is constructed through
$a_{i}\in\mathcal{A}=\\{(a_{ITN},a_{IRS})\\}$ with $a_{ITN},a_{IRS}\in(0,1]$,
which represent the population coverage for ITNs and IRS of a specific area.
Without significantly affecting performance, we discrete the action space with
an accuracy of 0.1 for simplicity.

#### Reward

The reward is a scalar $r_{t}=R(s_{t+1})\in[R_{\min},R_{\max}]$, associated
with next state $s_{t+1}$. In malaria control, it is determined through an
economic cost-effectiveness analysis. In Bent et al. (2018), an overview of
the reward calculation is specified. Without loss of generality, the reward
function $R(\ast)$ is assumed known to us since an MDP with unknown rewards
and unknown transitions can be represented as an MDP with known rewards and
unknown transitions by adding additional states to the system.

#### State

The state contains important observations for decision making in every time
step. In malaria control, it includes the number of life with disability, life
expectancy, life lost, and treatment expenses Bent et al. (2018). We set the
state in the form $s_{t}=(r_{t-1},a_{t-1},t)$, including current reward
$r_{t-1}$, previous action $a_{t-1}$ and the current intervention timestamp,
which covers the most crucial and useful observations for malaria control, as
shown in Figure 2. For the start state $s_{1}$, the reward and action are
initialized with 0.

Let $\pi:\mathcal{S}\rightarrow\mathcal{A}$ denote a deterministic mapping
from states to actions, and let
$V_{\pi}(s_{t})=\mathbb{E}_{a_{i}\sim\pi(s_{i})}\left[\sum_{i=t}^{T}\gamma^{i-t}r_{i}\right]$
denote the expected discounted reward by following policy $\pi$ in state
$s_{t}$. The objective is to find a deterministic policy $\pi^{\ast}$ that
maximizes the expected return at state $s_{1}$ as

$\displaystyle\pi^{\ast}={\arg\max}_{\pi\in\Pi}V_{\pi}(s_{1}).$

###  3.2 Model-based Indirect Policy Search

In the following, we detail the key components of the proposed framework VB-
MCTS, including the world model, the planner, and variance-bonus reward with
its sample complexity.

#### World Model Learning

The probabilistic world model is implemented as a GP, where we use a
predefined feature mapping $x_{t}=\phi(s_{t},a_{t})\in\mathbb{R}^{m}$ as
training input and the target state $s_{t+1}$ as the training target. The GP
yields one-step predictions

$\displaystyle p(s_{t+1}|s_{t},a_{t})$ $\displaystyle=$
$\displaystyle\mathcal{N}(s_{t+1}|\mu_{t},\sigma^{2}_{t})$
$\displaystyle\mu_{t}$ $\displaystyle=$
$\displaystyle\bm{k}_{\ast}^{\top}(K+\omega^{2}_{n}I)^{-1}\bm{y},$
$\displaystyle\sigma^{2}_{t}$ $\displaystyle=$ $\displaystyle
k(x_{t},x_{t})-\bm{k}_{\ast}^{\top}(K+\omega^{2}_{n}I)^{-1}\bm{k}_{\ast},$

where $k$ is the kernel function, $\bm{k}_{\ast}\equiv k(\bm{X},x_{t})$
denotes the vector of covariances between the test point and all training
points with $\bm{X}=[x_{1},\cdots,x_{n}]$, and
$\bm{y}=[s_{1},\dots,s_{n}]^{\top}$ is the corresponding training targets.
$\omega_{n}$ is the noise variance. $K$ is the Gram matrix with entries
$K_{ij}=k(x_{i},x_{j})$.

Throughout this paper, we consider a prior mean function $m\equiv 0$ and a
squared exponential (SE) kernel with automatic relevance determination. The SE
covariance function is defined as

$\displaystyle
k(x,x^{\prime})=\alpha^{2}\exp(-\frac{1}{2}(x-x^{\prime})^{\top}\mathbf{\Lambda}^{-1}(x-x^{\prime})),$

where $\alpha^{2}$ is the variance of state transition and
$\bm{\Lambda}\equiv\operatorname{diag}\left(\left[l_{1}^{2},\ldots,l_{m}^{2}\right]\right)$.
The characteristic length-scale $l_{i}$ controls the importance of $i$-th
feature. Given $n$ training inputs $\bm{X}$ and the corresponding targets
$\bm{y}$, the posterior hyper-parameters of GP (length-scales $l_{i}$ and
signal variance $\alpha^{2}$) are determined through evidence maximization
technique Williams and Rasmussen (2006).

#### Exploration with Variance-Bonus Reward

The algorithm we propose is itself very straightforward and similar to many
previously proposed exploration heuristics Kolter and Ng (2009); Srinivas et
al. (2009); Sorg et al. (2012); Grande (2014). We call the algorithm Variance-
Bonus Reward, since it chooses action according to the current mean estimation
of the reward plus an additional variance-based reward bonus for state-actions
that have not been well explored as

$\displaystyle\tilde{R}(s_{t+1}|s_{t},a_{t})=R(s_{t+1})|_{s_{t+1}=\mathbb{E}_{s_{t+1}\sim
GP(s_{t},a_{t})}[s_{t+1}]}+$ $\displaystyle\beta_{1}\text{Var}_{s_{t+1}\sim
GP(s_{t},a_{t})}[R(s_{t+1})]+\beta_{2}\text{Var}_{s_{t+1}\sim
GP(s_{t},a_{t})}[s_{t+1}],$

where $\beta_{1}$ and $\beta_{2}$ are the parameters that trade-off the
balance of exploitation and exploration. $\text{Var}_{s_{t+1}\sim
GP(s_{t},a_{t})}[R(s_{t+1})]$ and $\text{Var}_{s_{t+1}\sim
GP(s_{t},a_{t})}[s_{t+1}]$ are the predicted variances for state and reward.
The variance of reward can be exactly computed following the law of iterated
variances.

#### Planning with Mean MDP + Reward Bonus

MCTS is a strikingly successful planning algorithm Silver et al. (2017), which
can find out the optimal solution with enough computation resources. Since
disease control is not a real-time task, we propose to apply MCTS (Figure 1)
as the planner for generating the policy to maximize the variance-bonus
reward. In the executing process, the MCTS planner incrementally builds an
asymmetric search tree guided to the most promising direction by a tree
policy. This process usually consists of four phases — selection, expansion,
evaluation, and backup (as shown in Figure 3).

Specifically, each edge $(s,a)$ of the search tree stores an average action
value $\tilde{Q}(s,a)$ and visit count $N(s,a)$. In the selection phase,
starting from the root state, the tree is traversed by simulation (that is,
descending the tree with the mean prediction of states without backup). At
each time step $t$ of each simulation, an action $a_{t}$ is selected from
state $s_{t}$

$\displaystyle
a_{t}={\arg\max}_{a\in\mathcal{A}}(\tilde{Q}(s_{t},a)+\frac{c_{puct}}{|\mathcal{A}|}\frac{\sqrt{\sum_{b\in\mathcal{A}}N(s_{t},b)}}{1+N(s_{t},a)}),$

so as maximize action value plus a bonus that decays with repeated visits to
encourage exploration for tree search. Here, $c_{puct}$ is the constant
determining the level of exploration. When the traversal reaches a leaf node
$s_{L}$ at step $L$, the leaf node may be expanded with each edge initialized
as $\tilde{Q}(s_{L},a)=0$, $N(s_{L},a)=0$, and the corresponding leaf nodes
are initialized with the mean GP prediction
$s_{L+1}=\mathbb{E}[s_{L+1}|s_{L},a]$. Then, the leaf node is evaluated by the
average outcome $v(s_{L})$ of rollouts, which played out until terminal step
$T$ using the fast rollout policies, such as random policy and greedy policy.
At the end of the simulation, the action values and visit counts of all
traversed edges are updated, _i.e.,_ backup. Each edge accumulates the visit
counts and means the evaluation of all simulations passing through that edge
as

$\displaystyle\tilde{Q}(s_{j},a_{j})$ $\displaystyle\leftarrow$
$\displaystyle\frac{N(s_{j},a_{j})\times\tilde{Q}({s}_{j},{a}_{j})+v({{s}_{j}})}{N(s_{j},a_{j})+1},$
$\displaystyle v({s}_{j})$ $\displaystyle\leftarrow$ $\displaystyle
v({{s}_{j+1}})+\tilde{R}(s_{j+1}|s_{j},a_{j}),$ $\displaystyle N(s_{j},a_{j})$
$\displaystyle\leftarrow$ $\displaystyle N(s_{j},a_{j})+1,$

where $(s_{j},a_{j})$ with $j<L$ is the edge in the forward trace.

Figure 3: An illustration of MCTS for policy generation.

This cycle of selection, evaluation, expansion and backup is repeated until
the maximum iteration number has been reached. At this point, the best action
has been chosen by selecting the action that leads to the highest reward (max
child) as follows

$\displaystyle\pi(s_{t})=\arg\max_{a\in\mathcal{A}}\tilde{Q}(s_{t},a),$ (1)

where $\pi(s_{t})$ is the policy generated by MCTS. To be noticed, the
variance-bonus reward $\tilde{R}$ can be replaced with the $R$ for generating
the best-known policy since it maximizes the expected reward based on the
posterior so far.

Input: Hype-parameters $\beta$, maximum trials number $N$.

Output: The malaria control policy $\bm{\pi}$.

# Initialization.

Generate random policy $\bm{\pi}=\left[a_{1},\dots,a_{T}\right]$.

Apply policy $\bm{\pi}$ in real world and collect training samples
$\mathcal{D}=\\{(s_{i},a_{i},s_{i+1})\\}_{i=1}^{T}$.

Update GP world model with $\mathcal{D}$.

Set number of trials $k=1$.

for  _$k <=N$_ do

Initialize time step $t=1$, and state $s_{t}$.

for  _$t <=T$_ do

# Data Collection.

Generate action $a_{t}=\arg\max_{a\in\mathcal{A}}Q(s_{t},a)$ by MCTS with
variance-bonus reward $\tilde{R}$.

Apply $a_{t}$ in real world and collect $s_{t+1}$.

Update training samples
$\mathcal{D}=\mathcal{D}\cup\\{(s_{t},a_{t},s_{t+1})\\}$.

# Model Training.

Update GP world model with $\mathcal{D}$.

Update $s_{t}=s_{t+1}$, and $t=t+1$.

end for

Increase the training trials $k=k+1$.

end for

# Generate final decision.

Initialize time step $t=1$, state $s_{t}$, and policy
$\bm{\pi}=\left[\right]$.

for  _$t <=T$_ do

Generate action $a_{t}=\arg\max_{a\in\mathcal{A}}Q(s_{t},a)$ by MCTS with mean
reward $R$.

Predict $s_{t+1}=\mathbb{E}[s_{t+1}|s_{t},a_{t}]$ with Gaussian Process.

Update $s_{t}=s_{t+1}$, and $t=t+1$.

Append action to final decision $\bm{\pi}=\bm{\pi}\oplus a_{t}$.

end for

Return the malaria control policy $\bm{\pi}$.

ALGORITHM 1 Variance-Bonus Monte Carlo Tree Search

Finally, we implement an iterative training procedure, as shown in Algorithm
1, where we specify the order in which each component occurs within the
iteration.

###  3.3 Sample Complexity

In the following, we derive sample complexity (_e.g.,_ the required samples to
learn near-optimal performances) for our proposed solution. In the worst case,
when the reward function is equal everywhere, and all state-action pairs will
be equally explored, VB-MCTS will have the same complexity bounds as over-
exploration methods, such as GP-Rmax Grande et al. (2014). In practice, VB-
MCTS will learn the optimal policy in fewer steps than the over-exploration
methods because the high uncertainty region with low reward is not worth
exploration. Following the general PAC-MDP theorem, Theorem 10 in Strehl et
al. (2009), we derive the polynomial sample complexity of VB-MCTS.

######  Theorem 1.

Assume that the feature space $\mathcal{X}\subset\mathbb{R}^{m}$ of state-
action pairs is a compact domain and the target value is bounded
$y\in[0,V_{m}]$. The reward and action value are Lipschitz continuous w.r.t
the state-action pairs and let $L_{q}$ and $L_{r}$ be the Lipschitz constants
for reward and action value respectively. If VB-MCTS is executed with
$\beta_{1}=\frac{L_{r}}{2\sigma_{n}^{2}}$ and
$\beta_{2}=\frac{L_{p}}{2\sigma_{n}^{2}}$ for any MDP $M$, then with
probability $1-\delta$, VB-MCTS will follow a $4\epsilon$-optimal policy from
its current state on all but

$\displaystyle
O\left(\frac{V_{m}\zeta(\epsilon,\delta)}{\epsilon(1-\gamma)}\log\frac{1}{\delta}\log\frac{1}{\epsilon(1-\gamma)}\right)$
(2)

timesteps, with probability at least $1-2\delta$, where
$\zeta(\epsilon,\delta)=\left(\frac{4V_{m}^{2}}{\epsilon^{2}(1-\gamma)^{2}}\log\left(\frac{2\left|\mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\text{tol
}}^{2}\right)\right)\right|}{\delta}\right)\right)\mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\text{tol
}}^{2}\right)\right)$,
$\sigma_{\text{tol}}^{2}=\frac{2\omega_{n}^{2}\epsilon_{1}^{2}}{V_{m}^{2}\log\left(\frac{2}{\delta_{1}}\right)}$.
$\mathcal{N}_{\mathcal{X}}(\tau)$ is the covering number of domain
$\mathcal{X}$, which defines the cardinality of the minimal set
$C=\left\\{c_{i},\ldots,c_{N_{c}}\right\\}\text{ s.t. }\forall
x\in\mathcal{X},\exists c_{j}\in C\text{ s.t. }d\left(x,c_{j}\right)\leq\tau,$
where $d(\cdot,\cdot)$ is a distance measure.

###### Sketch of Proof.

The proof is organized by showing the required key properties (optimism,
accuracy, and learning complexity) for general PAC-MDP Strehl et al. (2009)
are satisfied.

#### Optimism

Assume that MCTS can find the optimal decision under every state $s$ with
enough computation resources. Then, $\forall(s_{t},a_{t})$, the optimal action
value $Q^{\ast}(s_{t},a_{t})$ satisfies

$\displaystyle
Q^{*}\left(s_{t},a_{t}\right)=\arg\max_{a_{t+1}\in\mathcal{A}}R\left(s_{t+1}\right)+\gamma
Q^{*}\left(s_{t+1},a_{t+1}\right)$ $\displaystyle\leq$ $\displaystyle
R\left(s_{t+1}^{\prime}\right)+\frac{L_{r}}{\omega_{n}^{2}}\operatorname{Var}_{s_{t+1}^{\prime}\sim
GP\left(s_{t},a_{t}\right)}\left[R\left(s_{t+1}^{\prime}\right)\right]$
$\displaystyle+$
$\displaystyle\gamma\tilde{Q}^{*}\left(s_{t+1}^{\prime},a_{t+1}\right)+\gamma\frac{L_{q}}{\omega_{n}^{2}}\operatorname{Var}_{s_{t+1}^{\prime}\sim
GP\left(s_{t},a_{t}\right)}\left[s_{t+1}^{\prime}\right]$ $\displaystyle\leq$
$\displaystyle\arg\max_{a_{t+1}^{\prime}\in\mathcal{A}}R\left(s_{t+1}^{\prime}\right)+\frac{L_{r}}{\omega_{n}^{2}}\operatorname{Var}_{s_{t+1}^{\prime}\sim
GP\left(s_{t},a_{t}\right)}\left[R\left(s_{t+1}^{\prime}\right)\right]$
$\displaystyle+$
$\displaystyle\gamma\tilde{Q}^{*}\left(s_{t+1}^{\prime},a_{t+1}^{\prime}\right)+\gamma\frac{L_{q}}{\omega_{n}^{2}}\operatorname{Var}_{s_{t+1}^{\prime}\sim
GP\left(s_{t},a_{t}\right)}\left[s_{t+1}^{\prime}\right]$ $\displaystyle=$
$\displaystyle\tilde{Q}\left(s_{t},a_{t}\right),$

where $s_{t+1}{=}P(s_{t},a_{t})$, $s^{\prime}_{t+1}{=}\mathbb{E}_{s_{t+1}\sim
GP(s_{t},a_{t})}[s_{t+1}|s_{t},a_{t}]$. The maximum error of propagating
$Q(s_{t+1},a_{t+1})$ instead of $Q(s^{\prime}_{t+1},a_{t+1})$ is given by
$\|s^{\prime}_{t+1}-s_{t+1}\|L_{q}$. The norm $\|s^{\prime}_{t+1}-s_{t+1}\|$
is upper bounded by the regularization error
$\frac{L_{q}}{\omega_{n}^{2}}\operatorname{Var}_{s_{t+1}^{\prime}\sim
GP\left(s_{t},a_{t}\right)}\left[s_{t+1}^{\prime}\right]$ with probability of
$(1-\delta)$ (Section A.2 in Grande (2014)). In the same way, we have
$R(s^{\prime}_{t+1})+\frac{L_{r}}{\omega^{2}_{n}}\text{Var}_{s^{\prime}_{t+1}\sim
GP(s_{t},a_{t})}[R(s^{\prime}_{t+1})]-R(s_{t+1})\geq 0$ with probability
$1-\delta$. It follows that Algorithm 1 remains optimistic with probability at
least 1-2$\delta$.

#### Accuracy

Following the Lemma 1 from Grande et al. (2014), the prediction error at
$x_{i}\in\mathcal{X}$ is bounded in probability
$Pr\left\\{\left|\hat{\mu}\left(x_{i}\right)-f\left(x_{i}\right)\right|\geq\epsilon_{1}\right\\}\leq\delta_{1}$
if the predictive variance of the GP at
$\sigma_{n}^{2}\left(x_{i}\right)\leq\sigma_{\text{tol}}^{2}=\frac{2\omega_{n}^{2}\epsilon_{1}^{2}}{V_{m}^{2}\log\left(\frac{2}{\delta_{1}}\right)}$.

#### Learning complexity

(a) For $\forall x\in\mathcal{X}$ with $\sigma_{0}^{2}=1$, let
$\rho=\min\limits_{j=1,\dots,n}k(x,x_{j})$,
$d_{\max}=\max\limits_{j=1,\dots,n}d(x,x_{j})$, then $d_{\max}$ can be
rewritten as $d_{\max}=\sqrt{2\theta^{2}\log{\frac{1}{\rho}}}$ with
$\rho=e^{-\frac{d_{\max}^{2}}{2\theta^{2}}}$. If
$n\geq\frac{\omega_{n}^{2}}{\sigma_{\text{tol}}^{2}-1+\rho^{2}}$, the
posterior variance of $x$ satisfies $\sigma_{n}^{2}(x)\leq
1-\frac{n\rho^{2}}{n+\omega^{2}_{n}}\leq\frac{n(1-p^{2})+\omega^{2}_{n}}{n}\leq\sigma^{2}_{\text{tol}}$.
(b) Since $\mathcal{X}$ is a compact domain with the length of $i$-th
$\left(i\in[1,m]\right)$ dimension as $L_{i}$, we have
$\mathcal{N}_{\mathcal{X}}\left(d_{\max}\right)\leq\mathcal{N}_{\mathcal{X}}\left(\frac{d_{\max}}{2}\right)$
and $\mathcal{X}$ can be covered by $2^{m}\prod_{j=1}^{m}L_{j}/d_{\max}^{m}$
balls $\left\\{B\left(c_{j},\frac{d_{\max}}{2}\right)\right\\}_{j}$, which
centers at $c_{j}$ with a radius of $\frac{d_{\max}}{2}$.

Given (a) and (b), $\forall x\in\mathcal{X}$, $\exists c_{j}\in\mathcal{X}$
s.t. $B(x,d_{\max})\supset B(x,d_{\max}/2)$. If there are at least
$\frac{\omega_{n}^{2}}{\sigma_{\text{tol}}^{2}-1+\rho^{2}}$ observations in
$B(x,d_{\max}/2)$, then (2) satisfied. Combining the Lemma 8 from Strehl et
al. (2009) and denoting $\tau(\sigma^{2}_{\text{tol}})={d_{\max}}/2$, the
total number of updates occurs will be bounded by
$\zeta(\epsilon,\delta)=\left(\frac{4V_{m}^{2}}{\epsilon^{2}(1-\gamma)^{2}}\log\left(\frac{2\left|\mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\mathrm{tol}}^{2}\right)\right)\right|}{\delta}\right)\right)\mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\mathrm{tol}}^{2}\right)\right)$
with probability of $1-\delta$.

Now that the key properties of optimism, accuracy, and learning complexity
have been established, the general PAC-MDP Theorem 10 of Strehl et al. (2009)
is invoked. ∎

Theorem 1 allows us to guarantee that the number of steps in which the
performance of VB-MCTS is significantly worse than that of an optimal policy
starting from the current state is at most log-linear in the covering number
of the state-action space with probability $1-\delta$.

##  4 Experiments

In this section, we conduct extensive experiments on two different OpenMalaria
Smith et al. (2008) based simulators: SeqDecChallenge and ProveChallenge,
which are the testing environments used in KDD Cup 2019. In these simulators,
the parameters of the simulator are hidden from the RL agents since the
“simulation parameters” for SSA are unknown for policymakers. Additionally, to
simulate the real disease control problem, only 20 trials were allowed before
generating the final decision, which is much more challenging than traditional
RL tasks. These simulating environments are available at
https://github.com/IBM/ushiriki-policy-engine-library.

#### Agents for Comparison

To show the advantage of VB-MCTS, many benchmarking reinforcement learning
methods and open source solutions have been deployed to verify its
effectiveness: Random Policy: The random policy is executed in 20 trials and
chooses the generated policy with the maximum reward as the final decision.
SMAB: This kind of policy treats the problem as a Stochastic Multi-Armed
Bandit problem and independently optimizes the policy every year with Thompson
sampling Chapelle and Li (2011). CEM: Cross-Entropy Method is a simple
gradient-free policy searching method Szita and Lörincz (2006). CMA-ES: CMA-ES
is a gradient-free evolutionary approach to optimizing non-convex objective
functions Krause et al. (2016). Q-learning-GA: It learns the malaria control
policy by combining Q-learning and Genetic Algorithm. Expected-Sarsa: It
collects 13 random episodes and runs expected value SARSA Van Seijen et al.
(2009) for 7 episodes to improve the best policy using the collected
statistics. GP-Rmax: It uses GP learners to model $T$ and $R$ and replaces the
value of any $Q^{\ast}(s,a)$ where $T(s,a)$ is ”unknown” with a value of
$\frac{R_{\max}}{1-\gamma}$ Li et al. (2011); Grande et al. (2014). GP-MC: It
employs Gaussian Process to regress the world model. The policy is generated
by sampling from the posterior and choosing the max rewarded action. VB-MCTS:
Our proposed method.

#### Implementation Details

We build a feature map of 14-dimension for this task, which includes the
periodic feature and cross term feature. Specifically, the feature map is set
as,

$\displaystyle\phi(s_{t},a_{t})=[t,t\%2,t\%3,r_{t-1},a_{(t-1,{ITN})},a_{(t-1,{IRS})},$
$\displaystyle a_{(t,{ITN})},a_{(t,{IRS})},a_{(t,{ITN})}\times a_{(t,{IRS})},$
$\displaystyle a_{(t-1,{ITN})}\times a_{(t-1,{IRS})},a_{(t,{ITN})}\times
a_{(t-1,{ITN})},$ $\displaystyle a_{(t,{IRS})}\times
a_{(t-1,{IRS})},a_{(t,{ITN})}\times\left(1-a_{(t-1,{ITN})}\right),$
$\displaystyle a_{(t,{IRS})}\times\left(1-a_{(t-1,{IRS})}\right)],$

where $t,t\%2,t\%3$ are the periodic features and $a_{(\ast,{\ast})}\times
a_{(\ast,{\ast})}$ is the cross term feature. Since the predicted variance of
state and reward are the same in our setting, we empirically set the sum of
exploration/exploitation parameters as $\beta_{1}+\beta_{2}=3.5$ in the
experiments. In MCTS, $c_{puct}$ is 5, and only the top 50 rewarded child
nodes are expanded. The number of iterations does not exceed 100,000. For
Gaussian Process, to avoid overfitting problems, 5-fold cross-validation is
performed during the updates of the GP world model. Particularly, we use
1-fold for training and 4-fold for validation, which ensures the
generalizability of the GP world model over different state-action pairs. Our
implementation and all baseline codes are available at
https://github.com/zoulixin93/VB_MCTS.

Table 1: Performance comparisons between different agents.
---
Agents | SeqDecChallenge
Med. Reward | Max Reward | Min Reward
Random Policy | 167.79 | 193.24 | 135.06
SMAB | 209.05 | 386.28 | -6.44
CEM | 179.30 | 214.87 | 120.92
CMA-ES | 185.34 | 246.12 | 108.18
Q-Learning-GA | 247.75 | 332.40 | 171.33
Expected-Sarsa | 462.76 | 495.03 | 423.93
GP-Rmax | 233.95 | 292.99 | 200.35
GP-MC | 475.99 | 499.60 | 435.51
VB-MCTS | 533.38 | 552.78 | 519.61
Agents | ProveChallenge
Med. Reward | Max Reward | Min Reward
Random Policy | 248.25 | 464.92 | 55.24
SMAB | 18.02 | 135.37 | -56.86
CEM | 229.61 | 373.83 | 20.09
CMA-ES | 289.03 | 314.57 | 92.95
Q-Learning-GA | 242.97 | 325.24 | 88.70
Expected-Sarsa | 190.08 | 296.16 | 140.86
GP-Rmax | 287.45 | 371.49 | 153.98
GP-MC | 300.37 | 447.15 | 263.96
VB-MCTS | 352.17 | 492.23 | 259.97

###  4.1 Results

#### Main Results

In Table 1, we report the median reward, the maximum reward, and the minimal
reward over 10 independent repeat runs. The results are quite consistent with
our intuition. We have the following observations: (1) For the finite-horizon
decision-making problem, treating it as SMAB does not work, and the delayed
influence of actions can not be ignored in malaria control. As presented in
Table 1, SMAB’s performances have large variance and are even worse than
random policy in ProveChallenge. (2) Overall, two model-based methods (GP-MC
and VB-MCTS) consistently outperform the model-free methods (CEM, CMA-ES,
Q-learning-GA, and Expected-Sarsa) in SeqDecChallenge and ProveChallenge,
which indicates that empirically model-based solutions are generally more
data-efficient than model-free solutions. From the results in Table 1,
performances of model-free methods are defeated by model-based methods with a
large margin in SeqDecChallenge, and their performances are almost the same as
the random policy in ProveChallenge. (3) The proposed VB-MCTS can outperform
all the baselines in SeqDecChallenge and ProveChallenge. Compared with GP-MC
and GP-Rmax, the advantage from the efficient MCTS with variance-bonus reward
leads to the success on SeqDecChallenge and ProveChallenge.

Figure 4: Malaria policy learning curves in five different agents.

#### Data Efficiency

This paragraph compares data-efficiency (required trials) of VB-MCTS with
other RL methods that learn malaria policies from scratch. In Figure 4(a) and
4(b), we report agents’ performances after collecting every trial episode in
SeqDecChallenge and ProveChallenge. The horizontal axis indicates the number
of trials. The vertical axis shows the average performance after collecting
every episode. Figure 4(a) and 4(b) highlighted that our proposed VB-MCTS
approach (brown) requires on average only 8 trials to achieve the best
performances in SeqDecChallenge and ProveChallenge, including the first random
trial. The results indicate that VB-MCTS can outperform the state-of-the-art
method on both data efficiency and performance.

##  5 Conclusion and Future Work

We proposed a model-based approach employing the Gaussian Process to regress
the state transition for data-efficient RL in malaria control. By planning
with variance-bonus reward, our method can naturally deal with the dilemma of
exploration and exploitation by efficiently MCTS planning. Extensive
experiments conducted on the challenging malaria control task have
demonstrated the advantage of VB-MCTS over state-of-the-arts both on
performance and efficiency. However, the stationary setting of MDP may be
unrealistic due to the development of disease control tools and the evolution
of the disease. Therefore, data-efficient reinforcement learning under
nonstationary settings will be more realistic and more challenging task.

## References

  * Auer and Ortner [2010] Peter Auer and Ronald Ortner.  Ucb revisited: Improved regret bounds for the stochastic multi-armed bandit problem.  Periodica Mathematica Hungarica, 61(1-2):55–65, 2010.
  * Bent et al. [2018] Oliver Bent, Sekou L Remy, Stephen Roberts, and Aisha Walcott-Bryant.  Novel exploration techniques (nets) for malaria policy interventions.  In AAAI, 2018.
  * Chapelle and Li [2011] Olivier Chapelle and Lihong Li.  An empirical evaluation of thompson sampling.  In NIPS, 2011.
  * Deisenroth and Rasmussen [2011] Marc Deisenroth and Carl E Rasmussen.  Pilco: A model-based and data-efficient approach to policy search.  In ICML, 2011.
  * Grande et al. [2014] Robert Grande, Thomas Walsh, and Jonathan How.  Sample efficient reinforcement learning with gaussian processes.  In ICML, 2014.
  * Grande [2014] Robert Conlin Grande.  Computationally efficient gaussian process changepoint detection and regression.  PhD thesis, Massachusetts Institute of Technology, 2014.
  * Ha and Schmidhuber [2018] David Ha and Jürgen Schmidhuber.  Recurrent world models facilitate policy evolution.  In NeurIPS, 2018.
  * Holland [1992] John H Holland.  Genetic algorithms.  Scientific american, 267(1):66–73, 1992.
  * Kakade et al. [2003] Sham Kakade, Michael J Kearns, and John Langford.  Exploration in metric state spaces.  In ICML, 2003.
  * Kamthe and Deisenroth [2017] Sanket Kamthe and Marc Peter Deisenroth.  Data-efficient reinforcement learning with probabilistic model predictive control.  arXiv:1706.06491, 2017.
  * Kolter and Ng [2009] J Zico Kolter and Andrew Y Ng.  Near-bayesian exploration in polynomial time.  In ICML, 2009.
  * Krause et al. [2016] Oswin Krause, Dídac Rodríguez Arbonès, and Christian Igel.  Cma-es with optimal covariance update and storage complexity.  In NIPS, 2016.
  * Li et al. [2011] Lihong Li, Michael L Littman, Thomas J Walsh, and Alexander L Strehl.  Knows what it knows: a framework for self-aware learning.  Machine learning, 82(3):399–443, 2011.
  * Mnih et al. [2013] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  arXiv:1312.5602, 2013.
  * Moran [2007] Mary Moran.  The malaria product pipeline: planning for the future.  Health Policy Division, The George Institute for International Health, 2007.
  * Silver et al. [2017] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al.  Mastering the game of go without human knowledge.  Nature, 550(7676):354, 2017.
  * Smith et al. [2008] T Smith, N Maire, A Ross, M Penny, N Chitnis, A Schapira, A Studer, B Genton, C Lengeler, Fabrizio Tediosi, et al.  Towards a comprehensive simulation model of malaria epidemiology and control.  Parasitology, 135(13):1507–1516, 2008.
  * Sorg et al. [2012] Jonathan Sorg, Satinder Singh, and Richard L Lewis.  Variance-based rewards for approximate bayesian reinforcement learning.  arXiv:1203.3518, 2012.
  * Srinivas et al. [2009] Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger.  Gaussian process optimization in the bandit setting: No regret and experimental design.  arXiv:0912.3995, 2009.
  * Strehl et al. [2009] Alexander L Strehl, Lihong Li, and Michael L Littman.  Reinforcement learning in finite mdps: Pac analysis.  JMRL, 2009.
  * Stuckey et al. [2014] Erin M Stuckey, Jennifer Stevenson, Katya Galactionova, Amrish Y Baidjoe, Teun Bousema, Wycliffe Odongo, Simon Kariuki, Chris Drakeley, Thomas A Smith, Jonathan Cox, et al.  Modeling the cost effectiveness of malaria control interventions in the highlands of western kenya.  PloS one, 9(10):e107700, 2014.
  * Sutton et al. [2000] Richard S Sutton, David A McAllester, Satinder P Singh, and Yishay Mansour.  Policy gradient methods for reinforcement learning with function approximation.  In NIPS, 2000.
  * Szita and Lörincz [2006] István Szita and András Lörincz.  Learning tetris using the noisy cross-entropy method.  Neural computation, 18(12):2936–2941, 2006.
  * Van Seijen et al. [2009] Harm Van Seijen, Hado Van Hasselt, Shimon Whiteson, and Marco Wiering.  A theoretical and empirical analysis of expected sarsa.  In DPRL, 2009.
  * Williams and Rasmussen [2006] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * Zhou et al. [2020] Wenjun Zhou, Taposh Dutta Roy, and Iryna Skrypnyk.  The kdd cup 2019 report.  ACM SIGKDD Explorations Newsletter, 22(1):8–17, 2020.
  * Zou et al. [2019] Lixin Zou, Long Xia, Zhuoye Ding, Jiaxing Song, Weidong Liu, and Dawei Yin.  Reinforcement learning to optimize long-term user engagement in recommender systems.  In SIGKDD, 2019.
  * Zou et al. [2020a] Lixin Zou, Long Xia, Pan Du, Zhuo Zhang, Ting Bai, Weidong Liu, Jian-Yun Nie, and Dawei Yin.  Pseudo dyna-q: A reinforcement learning framework for interactive recommendation.  In WSDM, 2020.
  * Zou et al. [2020b] Lixin Zou, Long Xia, Yulong Gu, Xiangyu Zhao, Weidong Liu, Jimmy Xiangji Huang, and Dawei Yin.  Neural interactive collaborative filtering.  In SIGIR, 2020.