[1]A. Philip Dawid
2]Stephen Senn
Personalised Decision-Making without Counterfactuals
Personalised Decision-Making
[1]University of Cambridge<EMAIL_ADDRESS>[2]Statistical Consultant, Edinburgh<EMAIL_ADDRESS>
This article is a response to recent proposals by Pearl
and others for a new approach to personalised treatment decisions,
in contrast to the traditional one based on statistical decision
theory.  We argue that this approach is dangerously misguided and
should not be used in practice.

§ INTRODUCTION

In recent works
Judea Pearl and collaborators have set out an approach to personalised
treatment that is radically different from that based on traditional
statistical decision theory.  It is based on the conception that we
should care, not only about the outcome that actually materialises,
but also about the (necessarily unobserved, counterfactual) outcome
that, it is supposed, would have occurred under the treatment that was
not applied.  A similar conception forms the basis of other recent
work <cit.>.

We consider this approach to be dangerously misguided, and believe
that real harm will ensue if it is applied in practice.  We argue our
case from a number of different viewpoints, and explain why this
approach should not be regarded as a viable alternative to standard
statistical decision theory.

§.§ Basic set-up

The context is that of a “target” patient suffering from a disease,
for which a treatment is available.  The treatment is far from
perfect, so that not all treated patients recover, while some
untreated patients may recover anyway.  There is information available
on recovery rates for treated and untreated patients; both these rates
may depend on measured individual patient characteristics.  The basic
problem is to decide, on the basis of the target patient's own
characteristics, whether or not to treat him.  A variation is how to
prioritise patients for treatment when there are limited doses

We introduce notation as follows:

Treatment decision Binary decision variable $X$, coded 1 for
treat, 0 for don't treat

Response Binary stochastic variable $Y$, coded 1 for recovery,
0 for no recovery

Individual background characteristics Stochastic variable $L$,
potentially multivariate, unaffected by the treatment decision

We suppose that there are available substantial data on $(L,X,Y)$, from
either experimental or unconfounded observational studies on patients
we can regard as similar to the target[See prag for
further discussion of this point.], from which we can estimate,
essentially perfectly[In a Bayesian setting it is
straightforward to relax this condition, using predictive
distributions based on finite samples.  However the main issues are
most clearly expressed in the case of essentially known
probabilities.], the distribution of $Y$, conditional on $L$, under
either treatment intervention.  That is, we know the probability
$\Pr(Y=1 \mid L=l, X \setto x)$, for any value $l$ of $L$ and $x = 0$
or $1$.  (Here $X \setto x$ denotes an external intervention to set
$X$ to $x$.)

§.§ Outline

In DT we recall the straightforward decision-theoretic
analysis of this problem.  Then in approach we briefly
outline the approach proposed by Pearl , followed by some
critical comments in comments. analysis describes
this approach in more detail, following a logical path that relates it
to other problems, in particular the use of general covariate
information to strengthen conclusions, and the specific case of an
“intention to treat” covariate, whose properties can be identified
by combining experimental and observational data.  In ex we
give critical consideration to some examples from Mueller and
Pearl <cit.>.  prag notes some important
assumptions that are implicitly made in the analysis, and points out
that they are unlikely to hold in practice.  disc summarises
our analysis and conclusions.

§ DECISION-THEORETIC APPROACH

We first describe the standard decision-theoretic (DT) approach to
treatment selection.

§.§ Single patient decision problem

Consider first the case of the single target patient.  We have to
decide whether to offer this patient treatment, or not.

Having access only to the target patient's value $L=l$, our objective
is to choose the treatment that will maximise the probability of
recovery.  We should thus treat this patient if
$p:= \Pr(Y=1 \mid L=l, X \setto 1) > q:=\Pr(Y=1 \mid L=l, X \setto
0)$.  That is, we should treat just when $\cate(l) > 0$, where
$\cate(l) = p-q$ is the “conditional average treatment effect”.

If the outcome $Y$ is not necessarily binary, for example a survival
time, we need to associate a utility $U(y)$ with the outcome $y$, and
treat the patient just when
$\E\{U(Y) \mid L=l, X \setto 1)\} > \E\{U(Y) \mid L=l, X \setto 0)\}$.
Applied to the binary case this reduces to the prescription above (so
long as $U(1)> U(0)$).

In <cit.>, a companion paper to this one which treats
utilities explicitly, the above is termed the interventionist
utility and approach, and contrasted with the counterfactual
utility and approach implicit in <cit.> and
explicit in <cit.>.  Here we restrict to the binary
case and do not use utilities.

Faced with a large collection $G$ of patients to treat, and unlimited
supplies of the treatment, managing each patient (each with their own
value $l$ of $L$) according to the above rule will maximise the number
of recoveries.  That is, any other (deterministic or randomised)
decision rule that uses only the information on $L$ would lead to
fewer recoveries.

Consider a case where

\begin{eqnarray*}
\Pr(Y=1 \mid  L=l, X\setto 1) &=& 0.49\\
\Pr(Y=1 \mid  L=l, X\setto 0) &=& 0.21.
\end{eqnarray*}

The conditional average treatment effect is
$\cate(l) = 0.49 - 0.21 = 0.28$.  Since $\cate(l) > 0$, the optimal
action is to treat this patient.  If we have a large collection of
similar patients, with the same value $L=l$, they should all be
treated—in which case the overall proportion of recovered patients
will be $49\%$.  This is the best outcome that can be achieved by
any treatment strategy.

§.§.§ Missing information

It may happen that, while we have full information on $(L,X,Y)$ for
the study individuals, the value $l$ of $L$ for the target patient is
not available.  In that case we can not condition on $L=l$, and we
have no option but to base the management of the patient on the
unconditional probabilities $\Pr(Y=1 \mid X \setto x)$ ($x=1,0$).
Nothing is gained by, for example, trying to impute the unknown value
of $L$.  If this is not obvious (as it should be), suppose we tried to
do so.  The recovery probabilities, conditional on a hypothesised
value $l$ for $L$, are $\Pr(Y=1 \mid L=l, X \setto x)$ ($x=1,0$).  But
as we do not know $l$, we need to take the expectation of
$\Pr(Y=1 \mid L, X \setto x)$ over the distribution of $L$, when
setting $X\setto x$ (which is the known marginal distribution of $L$,
unaffected by the intervention).  And this is just the unconditional
recovery probability $\Pr(Y=1 \mid X \setto x)$.  [With
finite data, on taking account of known structure in the
interventional distributions of $(L,Y)$ it may be possible to
improve the estimation of $\Pr(Y=1 \mid X \setto x)$.  But this
still remains what we need to focus on.]

§.§ Unit selection

Again consider a large collection $G$ of patients $i=1,\ldots,N$, with
individual recovery probabilities
$p_i = \Pr(Y_i = 1 \mid L_i = l_i, X_i \setto 1)$,
$q_i = \Pr(Y_i = 1 \mid L_i = l_i, X_i \setto 0)$.  If we treat just
those in a subset $S$, the expected number of recoveries will be
$\sum_{i\in S} p_i + \sum_{i \in G \setminus S} q_i = \sum_{i\in G}
q_i + \sum_{i\in S} \cate_i$.  Consequently, to maximise this expected
number, we should choose $S$, subject to any constraints, to maximise
$\sum_{i\in S} \cate_i$.  If we have limited treatments available, we
should thus prioritise individuals in decreasing order of their (while of course not treating any one for whom $\cate<0$.)  Again, any
other policy (subject to the same constraints) will have a smaller
number of recoveries.

§.§ Potential outcomes?

The “potential outcome” approach to causal inference <cit.>
conceives of the existence, even prior to treatment choice, of the
pair of variables $\bY = (Y(1), Y(0))$, where $Y(x)$ denotes the value
that, it is supposed, $Y$ will take if intervention $X \setto x$ is
applied.  The pretreatment variables $(Y(1), Y(0), L)$ are supposed to
have a joint distribution, unaffected by treatment.  With this
interpretation, we have

\begin{equation}
\label{eq:interp}
\Pr(Y=y, L=l \mid X\leftarrow x) =   \Pr(Y(x)=y, L=l),
\end{equation}

and $p = \E\{Y(1) \mid L = l \}$, $q = \E\{Y(0) \mid L = l\}$,

In this approach, inference is ideally desired for the “individual
treatment effect”, $\ite := Y(1)-Y(0)$, which can take values $+1$
(treatment benefits the patient), $-1$ (treatment harms the patient)
or $0$ (treatment has no effect).  Then $\cate=\E(\ite \mid L = l)$.
However, typically $\ite$ is unobservable, since it is impossible
simultaneously both to treat and not to treat the same patient.  In
particular, no information can be gained about the dependence between
$Y(1)$ and $Y(0)$, nor about the distribution (marginal, or
conditional on $L$) of $\ite$.  All that can be inferred is the
(conditional) expectation, as above, of $\ite$, depending as this does
only on the individual distributions of $Y(1)$ and $Y(0)$, which can
be identified from experimental data.

In certain very special and atypical cases, essentially those where we
have a fully deterministic and completely understood mechanistic
system, it may be that the background knowledge $L$ is detailed enough
to support perfect prediction of the eventual response, under either
intervention.  Then we will know, in advance of treatment choice, both
potential outcome variables.  In this case $p = Y(1)$, $q = Y(0)$, and
$\cate = \ite$.  Clearly we should treat as many of those who will (we
know for sure) benefit from the treatment as we can.

However, in typical cases perfect prediction is impossible, and then
it is arguable whether the potential responses even have any
meaningful existence.  In any case, there is nothing to be gained by
trying to impute potential responses: as in missinf above
(taking $L = \bY$), we should again simply focus on $\cate = p-q$.

In summary, consideration of potential responses (even if regarded as
meaningful) does not add any value to the decision-theoretic approach.

§ THE APPROACH OF MUELLER AND PEARL <CIT.>

In contrast to the above decision-theoretic approach, Mueller amd
Pearl <cit.> (henceforth MP) opt to take potential
outcomes seriously, and focus attention on $\ite=Y(1)-Y(0)$.  They
argue that we should ideally aim to treat those patients having
$\ite=1$, for whom the treatment made a difference: they would not
have recovered without it.  It would be wasteful to treat a patient
with $\ite = 0$, for whome the treatment made no difference, and
positively harmful to treat a patient with $\ite = -1$, who would have
recovered if untreated, but not if treated.

However, this ideal is unattainable, as we will not know a patient's
$\ite$ before treatment.  Concern is therefore transferred to the
“probability of benefit”, $\pb = \Pr(Y(1)=1, Y(0)=0) = \Pr(\ite=1)$,
and the “probability of harm”,
$\ph = \Pr(Y(1)=0, Y(0)=1)= \Pr(\ite=-1)$, which are now regarded as
the criteria by which to assess any treatment strategy.

But not only can we not know a patient's  before the treatment
decision is made, we can not even know it later, when the outcome $Y$
is observed.  For if we treat the patient we will observe $Y = Y(1)$,
but can not then observe the counterfactual outcome $Y(0)$ relevant
when we don't treat; similarly, for an untreated patient we can
observe $Y(0)$, but not $Y(1)$.  So $\ite$ is always unobservable.
This means that, even with extensive data on other patients, it will
not be possible fully to identify $\pb$ and $\ph$.  Such data can,
however, be used to set interval bounds on these quantities.
MP <cit.> further show how combining experimental
and observational data can narrow these bounds.  In certain very
special cases the bounds narrow to a single point, leading to full
identification of  and .

§ COMMENTS ON THE APPROACH

Our comments on the MP programme are arranged along several

§.§ Philosophy

Potential responses such as $Y(0)$ and $Y(1)$, first introduced by
Neyman <cit.>, have been considered as fundamental to the
conduct of causal inference ever since reintroduced by
Rubin <cit.>.  However this conception was challenged by
Dawid <cit.>, who pointed out that, so far from
being fundamental, they are entirely unnecessary, and that a fully
satisfactory theory can be based on standard decision-theoretic
elements.  Indeed, there are serious philosophical objections to
regarding potential responses as having real existence.  Only if we
take a fully Laplacean view of the universe, in which the future of
the universe is entirely determined by its present state and the laws
of Physics, does this make any sense at all—and even then, it is
difficult to incorporate the whims of an unconstrained external agent
who decides whether or not to give treatment, or to account for the
effect of external conditions arising after treatment. Even under
Laplacean determinism, our ignorance of the information needed to
predict the future means that we are unable to make use of it.
Whether or not we believe in a deep-down deterministic universe, our
predictions of the future can only be based on the limited information
we do have at our disposal, and must necessarily be
probabilistic.[See Dawid <cit.> for an approach
to understanding non-extreme probabilities based on imperfect
information about a deterministic world.] Imagining what we could
know or do, if only we had more information than we actually do have,
is just pointless.

§.§ Applicability

Another important dimension of criticism is that the strong conditions
needed for application of the MP theory will almost never obtain in
practice.  See prag below for details.

§.§ Helpfulness

The output of an MP analysis will, at very best, be point estimates of
the probabilities of benefit and of harm—in most cases, we won't
even get these, but can only bound these quantities within an
interval.  But even when we have these quantities, it is far from
clear how they help to inform treatment decisions.

§.§ Ethics

Our final criticism is the simplest, but most incisive.  The treatment
decisions made using the DT approach are guaranteed to be better than
those made by any other decision rule, in the sense that they will
maximise the number of recoveries in the population.  So whenever the
MP approach leads to different decisions, it will produce a decrease
in the number of recoveries.  We find it hard to construe this as

§ ANALYSIS

Here we provide a deconstruction of the analysis of
MP <cit.>—which should not, however, be taken as
agreement with their arguments and interpretations.  There are a
number of crucial assumptions required, but to avoid cluttering the
argument we leave these implicit, postponing specification and
discussion of them to prag.

We develop the story-line in a number of stages.

In simple we consider the case where we have access to
experimental data on treatment $X$ and response $Y$, and show how this
can be used to bound the probabilities of benefit and of harm.  We
also discuss the special circumstances in which these interval bounds
shrink to a point.

In cov we further suppose that we can measure additional
covariate information $L$ on individuals.  If we have these values in
the experimental data, this additional information can lead to a
narrowing of the bounds for $\pb$ and $\ph$ for the target case, even
when $L$ for that case is unobserved.

itt introduces a particular, potentially useful, covariate,
“intention to treat”, $X^*$—the treatment that a patient (or their
doctor) would like to choose, if unconstrained.  This may well be
informative about their state of health, and thus their outcome.  In
some experiments it may be possible to obtain information about $X^*$,
and this can then be used as $L$ in cov.  However it will
often not be possible to observe $X^*$ in the experiment.
comb considers how this problem can be overcome by the
incorporation of observational data, if we can assume that, in such
data, the desired treatment was the one actually applied, so that
$X^* = X$ becomes observable.  The combination of experimental and
observational data allows us to identify the distribution of $X^*$
(together with the other variables), and so once again allows us to
apply the theory of cov to obtain improved bounds for $\pb$
and $\ph$, which are detailed in app.

§.§ Simplest case

We start by presenting the basis of the approach in the simplest case,
where the data are experimental, and there is no additional covariate
information.  We thus have access to the interventional response
probabilities $\Pr(Y=y \mid X\leftarrow x) = \Pr(Y(x)=y)$,
$(x,y=0,1)$.  What can be inferred, from these, about the
probabilities of benefit and of harm?

As described by Dawid and Musio <cit.>, it is helpful to
express the interventional probabilities in terms of parameters $\tau$
and $\rho$, where

\begin{eqnarray}
\label{eq:tau}
\tau &\defeq& \Pr(Y=1\mid X\setto 1) - \Pr(Y=1\mid X\setto 0)\\
\label{eq:rho}
\rho &\defeq& \Pr(Y=1\mid X\setto 1) - \Pr(Y=0\mid X\setto 0).
\end{eqnarray}

Then $\tau$ is the average treatment effect, , of $X$ on
$Y$, while $\rho=\Pr(Y=1\mid X\setto 1) + \Pr(Y=1\mid X\setto 0)-1$ is
a measure of how common the outcome is.

The transition matrix $(\Pr(Y=y \mid X\setto x))$ from $X$ to $Y$ is

\begin{equation}
\label{eq:P}
P = P(\tau,\rho) \defeq\left(
\begin{array}[c]{cc}
\half(1+\tau+\rho) & \half(1-\tau-\rho)\\
\half(1-\tau+\rho) & \half(1+\tau-\rho)
\end{array}
\right),
\end{equation}

where the row and column labels are implicitly $1$ and $0$ in that
order.  The necessary and sufficient condition for all the transition
probabilities to be non-negative is

\begin{equation}
\label{eq:rhotau}
|\tau| + |\rho| \leq 1.
\end{equation}

We have equality in (<ref>) only in the degenerate case that one of the entries of $P$ is 0.

We can express the joint distribution for $\bY=(Y(1),Y(0))$ as in

1c     $Y(0)=1$     1c$Y(0)=0$
$Y(1) = 1$      $\half(1+\rho-\xi)$      $\half(\xi+\tau)$     $\half(1+\tau+\rho)$

$Y(1) = 0$      $\half (\xi-\tau)$     $\half(1-\rho - \xi)$      $\half(1-\tau-\rho)$

$\half(1-\tau+\rho)$       $\half(1+\tau-\rho)$    $1$

Joint probability distribution of $(Y(1),Y(0)$

The margins are determined by (<ref>) (with $L$ absent) and
(<ref>); but the internal entries are indeterminate, having one
degree of freedom crystallised in the unspecified “slack variable”
$\xi$, which is not identified by the experimental data.  The only
constraint on $\xi$ is the logical one that all internal entries of
r0r1 be non-negative. This holds if and only if

\begin{equation}
\label{eq:basicineq}
|\tau| \leq \xi \leq  1-|\rho|.
\end{equation}

This interval information is all that can be concluded about the joint
distribution for $\bY$ when we have data on the behaviour of $Y$ under
intervention on $X$, and no additional information.

The interval (<ref>) shrinks to a point, so that the joint
distribution of $\bY$ is fully determined by the experimental data,
if and only if we have equality in (<ref>), , just when
$P$ is degenerate, so that, for some $x,y =0,1$,
$\Pr(Y=y \mid X \setto x) =0$.  That is to say, for at least one of
the interventions, the resulting outcome $Y$ can be predicted with
certainty—a most unusual state of affairs.  In this case
$\Pr(Y(x) = y) = 0$, so that both joint events
$(Y(x) = y, Y({\overline x})=0)$ and $(Y(x) = y, Y(\overline x)=1)$
(where $\overline x = 1-x)$ have probability $0$.

§.§.§ Benefit and harm

The probability of benefit $\pb$ is the upper right entry of
r0r1, $\pb = \Pr(Y(1)=1,Y(0)=0) = \half(\xi+\tau)$, which by
(<ref>) is bounded between
$\pb_- := \half(|\tau|+\tau) = \max\{\tau,0\}$ and
$\pb_+:=\half(1-|\rho|+\tau) = \min\{\Pr(Y=1 \mid X \setto 1), \Pr(Y=0
\mid X \setto 0)\}$.  The probability of harm is the lower left entry
of r0r1,
$\ph=\Pr(Y(1)=0,Y(0)=1) = \half(\xi-\tau) = \pb-\tau$.

For the case of simple, we have $\tau= 0.28$, $\rho= -0.3$.
Without any further information, we can only infer
$0.28 \leq \pb \leq 0.49$, and correspondingly $0 \leq \ph \leq 0.21$.

§.§ Covariate information

Now suppose that, again with experimental data, we can obtain
additional information on some pre-treatment covariate information $L$
(for simplicity assumed discrete), unaffected by intervention (so
$\Pr(L=l \mid X\leftarrow x) = \Pr(L=l)$, assumed known and positive).
We thus have access to the conditional interventional probabilities
$\Pr(Y=y \mid L=l, X \leftarrow x)$.

Let $\tau(l), \rho(l)$ be defined as in (<ref>) and (<ref>),
but with probabilities further conditioned on $L=l$.  If, for the
target case, we observe $L=l$, then we simply apply the above
analysis, conditional on $L=l$.  In particular, the joint distribution
for $\bY$, given $L=l$, will be as in r0r1, with
$\rho,\tau,\xi$ replaced, respectively, by $\rho(l),\tau(l),\xi(l)$,
where $\xi(l)$ is subject only to

\begin{equation}
\label{eq:sineq}
|\tau(l)| \leq \xi(l) \leq  1-|\rho(l)|.
\end{equation}

Finally, suppose that, while having access, from the experimental
data, to the probabilities $\Pr(Y=y \mid L=l, X \leftarrow x)$, we do
not observe $L$ for the target patient.  In this case (and unlike the
situation for decision theory) the additional background knowledge can
make a difference.  In r0r1 we now have
$\xi = \sum_s \xi(l)\times \Pr(L=l)$, and we get the new interval

\begin{equation}
\label{eq:impb}
L := \sum_s |\tau(l)|\times \Pr(L=l) \leq \xi \leq 1-\sum_s |\rho(l)| \times \Pr(L=l)=:U.
\end{equation}

Since $\tau= \sum_s \tau(l)\times \Pr(L=l)$,
$\rho = 1-\sum_s \rho(l) \times \Pr(L=l)$, this interval will be
strictly contained in that of (<ref>) so long as not all the
$(\tau(l))$, or not all the $(\rho(l))$, have the same sign.

The probability of benefit is now bounded below by
$\sum_l \pb_-(l)\Pr(L=l)$ and above by $\sum_l \pb_+(l)\Pr(L=l)$,
where $\pb_-(l)$ and $\pb_+(l)$ can be computed as in
benharm1 with $\tau$ and $\rho$ replaced by $\tau(l)$ and
$\rho(l)$, respectively.

Applying pointsimp, and noting $|\tau(l)|\leq 1-|\rho(l)|$,
all $l$, we see that the interval (<ref>) will reduce to a
point, yielding full identification of the joint distribution of
$\bY$, if and only if $|\tau(l)| = 1-|\rho(l)|$, all $l$, so that,
for each $l$, at least one of $\Pr(Y=y \mid L=l, X\setto x)$, for
$x,y=0,1$, is zero.  In this case, both
$\Pr(Y(x) = y, Y({\overline x})=0 \mid L=l)$ and
$\Pr(Y(x) = y, Y(\overline x)=1\mid L=l)$ will be $0$.  Knowing the
value of $L$ will then always allow us to predict at least one of
the interventional outcomes with certainty.  However, the relevant
$x$ and $y$ may vary with $l$, in which case such certainty will not
be possible in the absence of knowledge of $L$.

§.§.§ Observational data

Consider now the case that our data are observational, rather than
experimental.  Suppose we can observe a “sufficient covariate”: a
covariate $L$ such that, conditional on $L$, we can assume there is no
residual confounding.  That is to say, the observational probability
$\Pr(Y =y \mid L=l, X = x)$ can be equated with the interventional
probability $\Pr(Y =y \mid L=l, X \setto x)$.  To ensure meaningful
conditioning, we further need the positivity condition: in the
observational setting,

\begin{equation}
\label{eq:pos}
\Pr(L=l, X= x) >0\qquad\mbox{all $l$, and $x=0$ or 1}.
\end{equation}

We can then proceed exactly as in cov above.

§.§ Intention to treat

Allocation of treatment to patients can be usefully decomposed into
two steps:

Intention The patient, or their doctor, decides on which
treatment they would ideally want.  This decision will typically be
related to their health status and other background information that
could be predictive of recovery, so that we cannot regard those who
desire, and those who reject, active treatment as comparing like
with like.  This is the genesis of confounding.

We introduce a binary stochastic “intention to treat” (ITT)
variable $X^*$ to denote the treatment desired.

Application A treatment $X$ is imposed on the patient.

It is important to distinguish $X$ and $X^*$.[We should
further distinguish between imposed treatment and received
treatment, as in Dawid <cit.>.  Here we notate both as $X$,
hoping this will cause no confusion.  We write $X\setto x$ when $X$
refers to the imposed treatment, and $X = x$ when $X$ refers to the
received treatment.]  The ITT variable $X^*$ exists prior to
application of treatment, and can thus be regarded as independent of

\begin{equation}
\Pr(X^* = x^* \mid X \leftarrow x) = \Pr(X^*=x^*).
\label{eq:xint}
\end{equation}

This expresses the covariate nature of $X^*$.

We assume that, in an observational setting, the desired treatment is
the one that is actually administered (there being no reason to do
otherwise).  Thus the received treatment $X$ will be the same as the
desired treatment $X^*$.  In particular, since we observe $X$, we can
infer the value of $X^*$.

In an experiment, however, the treatment $X$ will be imposed (, by
randomization), in a way that will typically take no account of $X^*$.
Even though we can still conceive of the ITT variable $X^*$ as
existing, it may or—more usually—may not be possible to observe
it.  When $X^*$ is observable, it can be used, just like any other
covariate, to improve decision-making, as in DT (when $X^*$
is observed for the target patient), or, in the approach of MP, to
narrow the bounds on $\pb$ and $\ph$, as in cov.

§.§.§ ITT as a sufficient covariate

In an observational setting, where $X^*=X$ is observed, it is natural
to assume “distributional consistency” <cit.>: the
distribution of $Y$ given intended treatment $X^*=x$—and so, also,
given received treatment $X=x$—is the same as that of $Y$, given
$X^*=x$, under an imposed intervention $X \leftarrow x$ that happens
to coincide with the treatment that would have been chosen anyway:

\begin{equation}
\label{eq:distcons}
\Pr(Y=y \mid X^*=x,X = x) = \Pr(Y=y \mid X^*=x,X  \leftarrow x).
\end{equation}

For $x^*\neq x$, the event $(X^*=x^*,X = x)$ does not occur in the
observational regime, so we can interpret
$\Pr(Y=y \mid X^*=x^*,X = x)$ however we want, in particular as

\begin{equation}
\label{eq:distcons2}
\Pr(Y=y \mid X^*=x^*,X = x) = \Pr(Y=y \mid X^*=x^*,X  \leftarrow x),
\end{equation}

and then (<ref>) implies that (<ref>) holds for all

Properties (<ref>) and (<ref>) imply that $X^*$, which
is observed in the observational setting, behaves as a sufficient

§.§ Combination of data

It would be nice if, with observational data, we could profit from the
fact that $X^*$ is a sufficient covariate, as in obs.
However, this is not straightforward, since the positivity condition
(<ref>) fails: for $x^*\neq x$, even though we may assume
$\Pr(Y=y \mid X^*=x^*, X\setto x) = \Pr(Y=y \mid X^*=x^*, X = x)$, we
have no data to estimate the latter term.  Again, when our data are
experimental but we can not directly observe $X^*$, we can not
identify $\Pr(Y=y \mid X^*=x^*, X\setto x)$.  However, it turns out
that we can do so if we can also obtain observational data: the
combination of both types of data allows us, after all, to identify
$\Pr(Y=y \mid X^*=x^*, X\setto x)$, even for $x\neq x^*$.  This we
show in the following theorem.

Suppose we can identify the joint distribution of $X$ and $Y$ in the
observational context, where $0<\Pr(X=1)<1$, and can also identify the
distribution of $Y$ under either intervention $X \leftarrow x$
$(x=0,1)$.  Then, under conditions (<ref>) and (<ref>),
all the probabilities $\Pr(Y=y \mid X^*=x^*, X \leftarrow x)$
$(x,x^* = 0,1)$ are identified.  Specifically,

\begin{eqnarray}
\label{eq:11}
\Pr(Y=y \mid X^* = 1, X \leftarrow 1) &=& \Pr(Y=y \mid X=1)\\
\label{eq:00}
\Pr(Y=y \mid X^* = 0, X \leftarrow 0) &=& \Pr(Y=y \mid X=0)\\
\label{eq:10}
\Pr(Y=y \mid X^*=1, X \leftarrow 0) &=& \frac{\Pr(Y=y \mid  X \leftarrow 0) - \Pr(Y=y, X=0)}{\Pr(X=1)}\\
\label{eq:01}
\Pr(Y=y \mid X^*=0, X \leftarrow 1) &=& \frac{\Pr(Y=y \mid  X \leftarrow 1) - \Pr(Y=y, X=1)}{\Pr(X=0)}.
\end{eqnarray}

(<ref>)  and (<ref>) follow from (<ref>).

To identify $\Pr(Y=y \mid X^* = 0, X \leftarrow 1)$, we argue as
follows.  We have

\begin{eqnarray}
\nonumber  \Pr(Y=y \mid X \leftarrow 0) &=& \Pr(Y=y \mid X^*=0,  X\leftarrow 0)\times\Pr(X^*=0 \mid  X \leftarrow 0)\\
\nonumber                               &&{}+\Pr(Y=y \mid X^*=1,  X\leftarrow 0)\times\Pr(X^*=1 \mid  X \leftarrow 0)\\
\nonumber                               &=& \Pr(Y=y \mid X=0)\times\Pr(X=0)\\
\label{eq:follow}                               &&{}+\Pr(Y=y \mid X^*=1,  X\leftarrow 0)\times\Pr(X=1),
\end{eqnarray}

on using (<ref>) and (<ref>), and the fact that $X^*=X$
in the observational setting.  Since all the other terms in
(<ref>) are identifiable in either the observational or the
experimental context, and $\Pr(X=1)\neq 0$, we can solve for
$\Pr(Y=y \mid X^*=1, X\leftarrow 0)$, obtaining (<ref>).  Then
(<ref>) follows similarly.

The above proof relies on $X$ (and so $X^*$) being binary, but $Y$
need not be.  Versions of this argument have appeared in

The joint distribution of $(X^*, Y)$ under the intervention
$X \leftarrow x$ is then identified.

Follows since, by (<ref>),
$\Pr(X^* = x^* \mid X \leftarrow x) = \Pr(X=x^*)$ is identified in
the observational context.

Since $\Pr(Y=y \mid X^*=1, X\leftarrow 0)\geq 0$, , we deduce
from (<ref>) and (<ref>) the consistency constraint
$\Pr(Y=y \mid X\leftarrow x) \geq \Pr(Y=y, X=x)$, all $x, y$.  When
this fails, and that failure can not be ascribed to sampling
variation or bias, that is evidence of violation of the conditions
of prag below, that have, implicitly, been used to justify
the above argument.

dawid and joint express just what the combination of
observational and experimental data is doing for us: it allows us to
identify distributions involving the ITT variable $X^*$.

§.§ Benefit and harm

Taking now $X^*$ as our sufficient covariate $L$, we can apply the
formulae of (<ref>)–(<ref>) to compute the quantities
$\tau(x^*)$, $\rho(x^*)$ required for the analysis of cov.
Noting that $X^*=X$ in the observational regime, so that
$\Pr(X=x)=\Pr(X^*=x)$, we obtain

\begin{eqnarray*}
\label{eq:tau1}
{\Pr(X^*=1)} \times \tau(1) &=& {\Pr(Y=1) - \Pr(Y=1 \mid X \setto 0)}\\
{\Pr(X^*=0)}\times\tau(0) &=& {\Pr(Y=1 \mid X \setto 1)-\Pr(Y=1)}\\
{\Pr(X^*=1)} \times\rho(1) &=& K - \Pr(Y=0 \mid X \setto 0)\\
{\Pr(X^*=0)}\times\rho(0) &=& \Pr(Y=1 \mid X \setto 1)-K
\end{eqnarray*}

where $K= \Pr(Y=1, X=1) + \Pr(Y=0, X=0)$.  Then from (<ref>) we
bound $\xi$ within $(L,U)$, where

\begin{eqnarray*}
L &=& |\Pr(Y=1) - \Pr(Y=1 \mid X \setto 0)| + |\Pr(Y=1) - \Pr(Y=1 \mid X \setto 1)|\\
1-U &=& |\Pr(Y=0 \mid X \setto 0) - K| + |\Pr(Y=1 \mid X \setto 1) - K|
\end{eqnarray*}

Then $\pb$ lies in $(\half(L+\tau), \half(U+\tau))$, and
$\ph = \pb-\tau$ lies in $(\half(L-\tau), \half(U-\tau))$.  Although
expressed differently, these results agree with those of
MP <cit.>.

By points, the joint distribution of $\bY$, and in particular
$\pb$, $\ph$, will be point identified just when, for both $x^* = 0$
and $x^* = 1$, there exist $x$, $y$ such that
$\Pr(Y=y \mid X^* = x^*, X \setto x)= 0$.  In non-trivial cases we
will have $\Pr(Y = y \mid X=x) \neq 0$, in which case, by (<ref>)
and (<ref>), this would need to happen with $x \neq x^*$.  For
that, by (<ref>) and (<ref>), we require

\begin{equation}
\label{eq:a10}
\Pr(Y=y \mid X\setto 0) = \Pr(Y=y,  X=0)
\end{equation}

for either $y=1$ or $y=0$; as well as

\begin{equation}
\label{eq:a01}
\Pr(Y=y \mid X\setto 1) = \Pr(Y=y,  X=1)
\end{equation}

for either $y=1$ or $y=0$.

§ EXAMPLES

MP <cit.> consider two cases, in both of
which the interventional probabilities of recovery are as in our
simple, having $\Pr(Y=1 \mid X\setto 1) = 0.49$,
$\Pr(Y=1 \mid X\setto 0) = 0.21$, and so $\ate = 0.28$.  However, they
have different observational data.  We now analyse these in detail.

This example relates to females, for whom the observational joint
probabilities are as in f.

1c     $Y=1$     1c$Y=0$
$X =1$     $0.19$     $0.51$     $0.70$

$X = 0$     $0.21$     $0.09$     $0.30$
$0.40$     $0.60$     $1$

Joint observational distribution of $(X,Y)$ for females

Applying the formulae of app we find:

\begin{eqnarray*}
0.7 \times\tau(1) &=& 0.19\\
0.3 \times\tau(0) &=& 0.09\\
0.7 \times\rho(1) &=& -0.51\\
0.3 \times\rho(0) &=& 0.11
\end{eqnarray*}

It follows that $\pb_-(1) = \tau(1) = 19/70$.  Also,
$\pb_+(1) = \Pr(Y=1 \mid X^* =1, X\setto 1) = \Pr(Y=1 \mid X=1) =
19/70$.  Hence, given $X^* = 1$, we have exact identification:
$\pb(1) = 19/70$.  This occurs because
$\Pr(Y=1, X=0) = 0.21 = \Pr(Y=1 \mid X\setto 0)$, implying the
deterministic property $\Pr(Y=1 \mid X^* = 1, X \setto 0) = 0$: a
female who desires treatment will never recover if untreated.
Consequently such a female should be treated.

Also, $\pb_-(0) = \tau(0) = 0.3$, while
$\pb_+(0) = \Pr(Y=0 \mid X^* =0, X\setto 0) = \Pr(Y=0 \mid X=0) =
0.3$.  Given $X^* = 0$, we again have exact identification:
$\pb_-(0) = 0.3$.  This occurs because
$\Pr(Y=0, X =1) = 0.51 = \Pr(Y=0 \mid X \setto 1)$, so that
$\Pr(Y=0 \mid X^*=0, X\setto 1) = 0$: a female who does not desire
treatment will always recover if treated.  Again, such a female should
be treated.

Finally we obtain exact identification marginally: $\pb= 0.28$.
Correspondingly, $\ph = \pb-\tau = 0$.  As there is no possibility of
harm, any female should be treated.

All the above conclusions agree with the DT prescription, based on the
experimental data alone: since $\ate > 0$, a female should be treated.

For males, the observational joint probabilities are as in

1c     $Y=1$     1c$Y=0$
$X =1$     $0.49$     $0.21$     $0.70$

$X = 0$     $0.21$     $0.09$     $0.30$
$0.70$     $0.30$     $1$

Joint observational distribution of $(X,Y)$ for males

Proceeding similarly to f, we find that $\pb$ and $\ph$ are
again identified exactly: $\pb = 0.49$, $\ph = 0.21$.  Indeed,
$\Pr(Y=1, X=1) = 0.49 = \Pr(Y=1 \mid X\setto 1)$, implying the
deterministic property $\Pr(Y=1 \mid X^* = 0, X \setto 1) = 0$: a male
who does not desire treatment will never recover if treated.
Consequently such a male should not be treated.  Also,
$\Pr(Y=1, X=0)= 0.21 - \Pr(Y=1 \mid X \setto 0)$, so that
$\Pr(Y=1 \mid X^* = 1, X \setto 0) = 0$: a male who desires treatment
will never recover if untreated, so that such a male should be

However, if we do not observe $X^*$ for the target male patient, the
above does not tell us how to proceed.  We might try to balance $\ph$
($= 0.49$) and $\pb$ ($=0.21$) somehow: for example, treat just when
$\pb > \lambda \ph$ for some chosen value of $\lambda$.  In the light
of the clinical maxim primum non nocere, a value $\lambda = 3$
might be chosen—in which case the target male would not be
treated.[This argument parallels one in
MP <cit.>, having different numbers, and
$\lambda = 2$.]

By contrast, in the absence of knowledge of $X^*$ for the target male,
the DT approach would take no account of the observational data, again
focusing simply on $\ate = 0.28$—and so decide to treat.  In a large
population of similar cases, this would lead to an overall recovery
rate of $49\%$, the maximum possible; whereas the above strategy based
on balancing $\pb$ and $\ph$ would only have a $21\%$ recovery rate.
It is difficult to see how this could be regarded as ethical.

Consider another case.  Again, the interventional probabilities are
$\Pr(Y=1 \mid X\setto 1) = 0.49$, $\Pr(Y=1 \mid X\setto 0) = 0.21$,
with $\tau = 0.28$.  Now the observational joint probabilities are
as in x.

1c     $Y=1$     1c$Y=0$
$X =1$     $0.2$     $0.5$     $0.7$

$X = 0$     $0.1$     $0.2$     $0.3$
$0.3$     $0.7$     $1$

Another joint observational distribution of $(X,Y)$

We compute

\begin{eqnarray*}
0.7 \times \tau(1) &=& 0.09\\
0.3 \times \tau(0) &=& 0.19\\
0.7 \times \rho(1) &=& -0.39\\
0.3 \times \rho(0) &=& 0.09.
\end{eqnarray*}

Using (<ref>) we find $0.28 \leq \xi \leq 0.52$, whence
$\pb = \half(\xi+\tau)$ is bounded between $0.28$ and $0.40$—and so
$\ph = \pb-\tau$ lies between $0$ and $0.12$.  So, even with the aid
of the additional observational data, we have not been able to
identify these probabilities exactly.  And even if we were to resolve
the ambiguity somehow, for example by taking the midpoints of these
intervals as suggested by Li and Pearl <cit.>, we
would be no better off than we were in m, where trying to
balance $\pb$ against $\ph$ could lead to a decision opposite to the
rcommendation of the simple DT analysis, so leading to fewer

§ ASSUMPTIONS AND CRITICAL COMMENTS

Here we identify and discuss some of the assumptions underlying the
foregoing analyses.

§.§ Representative data

A fundamental assumption underlying both the decision-theoretic
analysis of DT and the alternative approach of
approach is that the data available for estimating the
interventional probabilities $\Pr(Y=y, L=l \mid X\setto x)$ are on
individuals who can be regarded as “similar to” (“exchangeable
with”) the target case, so that these estimated probabilities are
applicable to the target.[For application to the MP arguments
of approach, the representativeness assumption should
apparently be extended to the (typically unidentifiable) bivariate
distribution, along with the other variables, of the pair of
potential responses $(Y(1),Y(0))$.  For the interval-valued
inferences made, however, this is not crucial, since these allow for
arbitrary dependence in this bivariate distribution.]  In reality
this is highly implausible.  For example, a clinical trial will have
entry criteria and processes that make its subjects quite untypical of
the population from which they are drawn, or indeed of the individuals
recruited into another such trial.  In any case, despite the name,
entry criteria govern who does not get into a trial: they
cannot guarantee that those who enter are representative even of a
target individual meeting the same criteria.  A clinical trial gains
its value, not from representativeness, but from the internal
randomisation that ensures that a comparison between its treated and
untreated groups is indeed a comparision of like with like, and that
valid probability statements can be made about likely differences, so
enforcing internal validity.  Because of unrepresentativeness it would
not be appropriate to regard $\Pr(Y=y, L=l \mid X\setto x)$, estimated
from the data, as being directly relevant to the target case—the
problem of external validity.  (One cheating way round this is to
focus on a hypothetical target individual who can be regarded
as exchangeable with those in the study.)  Nevertheless, it may still
be reasonable to regard the estimated  or  as applying to
the target—if not in its exact numerical value, at least in its
sign, which is what is required, for DT application, to solve the
single patient treatment problem; or in its ordering of the $\cate_i$,
as required to solve the DT unit selection problem.

To underline how unreasonable the representative assumption is, it
should be noted that even when clinical trials with similar protocols
are compared this assumption is not made.  A striking example of its
failure for nearly identical protocols is given by the TARGET study
<cit.>, in which osteoarthritis patients in some centres
were randomised to receive either lumiracoxib or naproxen, and
patients in other centres either lumiracoxib or ibuprofen.  The degree
of comparability in design of the two sub-studies thus defined was
greater than one would typically expect between two randomised
controlled trials (RCTs), and a fortiori than between an RCT
and an observational study, such as MP consider.  Nevertheless, very
important differences at baseline were seen between the two
sub-studies, even though within-sub-study treatment arms were
comparable.  Furthermore, it was possible to demonstrate differences
at outcome between the two studies using lumiracoxib data only, a
striking illustration of a study effect.  It is generally accepted by
sponsors and regulators that as soon as concurrent control is
abandoned the greatest of care must be taken in drawing inferences.
Modern work on using data on historical controls to try and improve
the efficiency of clinical trials takes such study-to-study variation
as a given that must be allowed for <cit.>.

§.§ Combination of data

An essential requirement for the application of dawid is that
the observational and experimental datasets comprise similar
individuals, so that the same probabilities for $X, X^*, Y$ apply to
both groups.  This is even more implausible than the
representativeness of either group.  In particular, the assumption of
a common distribution for the desired treatment $X^*$, in both
datasets and in the target patient, is vital but highly questionable.
Even if we were to accept the arguments of
MP <cit.> based on combining observational and
experimental data, without this property they are simply irrelevant.

§.§ What do clinical trialists do in practice?

The key to using RCTs is to identify reasonable assumptions, and use
theory to transfer results from trial to practice.  A striking example
is given by bioequivalence studies.  The subjects are usually young
healthy volunteers, frequently male.  However, the results will be
used to decide on appropriate treatments for elderly frail patients,
some female.  There is no pretence of representativeness.  Instead,
tight control and sensible scales of analysis are used.  The purpose
of such studies is to compare two formulations in terms of
bioavailability, and this is typically done using a cross-over trial
in which each subject is their own control, the order of
administration being randomised.  On separate days, concentration of
the test and reference pharmaceuticals are measured, and the ratio of
the areas under the two concentration time curves (AUCs) are
calculated for each subject, then analysed over all subjects,
typically after log-transformation.  What is relevant for treating an
individual patient is their own AUC: too low and efficacy may be
disappointing, too high and the drug may not be tolerated.  However,
no inference is made from a bioequivalence study in terms of AUCs
alone, since they would be quite different in healthy volunteers and
patients.  Instead, the idea is that the ratio between test and
reference ought to be the same in volunteers and patients, and this
ratio can be used to make predictions as to how the test drug will
behave in clinical practice.  An interesting example of such a study
is reported by Shumaker and Metzler <cit.>.  They used a more
elaborate design in which test and reference drugs were given in a
double cross-over, thus permitting them to analyse the
formulation-by-subject interaction.  They were able to demonstrate
that there was no evidence of an individual bioequivalence effect:
although you could estimate the individual relative bioavailability,
using the average over all subjects would be superior than any such
naïve estimate.  This raises a further issue with MP, who assume that
individual causal effects are stable over time.  Moreover, typical
causal analysis assumes an infinite sample size, but no infinities are
available for individual subjects, and estimating individual causal
effects requires close attention to components of variance.
Bioequivalence studies are an extreme example, but the general idea of
transferring results using a suitable scale for analysis, and
back-transforming to a scale suitable for decision analysis, is
commonplace: see <cit.> for a general discussion and
<cit.> in the specific context of vaccine efficacy.  Of
course, as the COVID-19 pandemic has reminded us, there are no
guarantees.  Things that work at one time may not do so at another.
It behoves all those proposing solutions to be cautious and humble.

§ SUMMARY

We have given careful accounts of the DT and MP approaches to
individualised treatment choice.  The DT approach is simple in the
extreme, and selects the treatment strategy that maximises the number
of recoveries.  In contrast, the MP approach fixates on
philosophically questionable and unknowable counterfactual concepts,
and when its recommendations differ from those of DT will lead to
fewer recoveries.  This has been illustrated in a number of examples.

One feature of the MP approach is the combination of experimental and
observational data.  When some very strong, and practically
implausible, conditions are satisfied, this permits identification of
the distribution of a special covariate, the intention to treat (ITT).
As with any other covariate whose distribution is known, this can then
feed back to tighten the MP inferences.  But it would be better to
observe this—or any other—covariate in the target patient, which
would then lead to better results from the DT point of view.  In
particular we have shown that, in just those very special cases that
use of ITT leads to point identification of the MP probabilities of
benefit and of harm, knowledge of the target patient's ITT value
allows perfect prediction of the outcome under at least one of the
treatment interventions, and so to a trivial solution to the decision

The DT approach has a long history of fruitful application to an
enormous variety of fields, from clinical trials to rocket science.
Attempts to replace it with another approach, based on
counterfactuals, are totally unnecessary and dangerously misguided.
This approach should not be used in practice.

§ ACKNOWLEDGMENTS

We have benefited greatly from discussions with Mats Stensrud and
Aaron Sarvet.

Conflict of interest: Prof. Philip Dawid is a member of
the Editorial Board in the Journal of Causal Inference but was not
involved in the review process of this article.


# Large Scale Dynamo in a Primordial Accretion Flow – An Interpretation from
Hydrodynamic Simulation

Wei-Ting Liao  Department of Astronomy, University of Illinois, 1002 West
Green Street, Urbana, IL 61801, USA<EMAIL_ADDRESS>Matthew Turk  School
of Information Sciences, University of Illinois, 501 East Daniel Street,
Champaign, IL 61820, USA  Department of Astronomy, University of Illinois,
1002 West Green Street, Urbana, IL 61801, USA  Hsi-Yu Schive  Institute of
Astrophysics, National Taiwan University, No.1 Sec.4 Roosevelt Road Taipei
10617, Taiwan  Center for Theoretical Physics, National Taiwan University,
No.1 Sec.4 Roosevelt Road Taipei 10617, Taiwan

###### Abstract

Without an existing large scale coherent magnetic field in the early Universe,
Population III (PopIII) stars would likely rotate at or near break-up speed.
In this work, focusing on the accretion phase of PopIII stars, we investigate
the possibility of generating a coherent magnetic field through large scale
dynamo processes, as well as the corresponding field saturation level. Using
results from hydrodynamic simulations, we demonstrate that primordial
accretion disks are turbulent with a Shakura-Sunyaev disk parameter
$\alpha_{ss}\gtrsim 10^{-3}$, and evidence helical turbulence with a dynamo
number $|D_{\alpha\Omega}|\gg 10$. The presence of helical turbulence at these
levels allows large scale dynamo modes to grow, and the saturation level is
determined by the amount of net helicity remaining in the dynamo-active
regions (a.k.a. the quenching problem). We demonstrate that, if the accretion
could successfully alleviate the quenching problem, the magnetic field can
reach approximate equipartition with $B/B_{\rm eq}\sim 3$.

##  1 Introduction

Although the formation of PopIII stars is one of the rarely well posed
problems in Astronomy, with a properly understood initial condition from the
cosmic microwave background measurement and a relative simple chemistry
network that could be followed numerically, the currently unknown magnetic
field in the early Universe presents itself as one of the biggest
uncertainties.

Magnetic fields in the present day star forming regions are known to play the
major role in regulating the star formation processes, including braking the
rotation of stars through magnetic braking and/or through star-disk
interactions (Mouschovias, 1991; Matt & Pudritz, 2005; Matt et al., 2010). It
could also cause the accretion disk around the protostar to be turbulent
through magneto-rotational instability (MRI) (Velikhov, 1959; Chandrasekhar,
1960; Balbus & Hawley, 1991). Additionally, the presence of magnetic field is
known to suppress the fragmentation (Hennebelle & Teyssier, 2008; Price &
Bate, 2008), a coherent strong magnetic field could launch collimated jets
(Lee et al., 2018), and, in some conditions, magnetic field could even forbid
the disk formation – the braking catastrophe (Mellon & Li, 2008, 2009; Li et
al., 2011).

In the early Universe, magnetic field is largely unknown, and is difficult, if
not impossible, to measure. Theoretically, several magnetogenesis models have
been proposed, yet it is hard to produce a strong and large scale coherent
field (Biermann, 1950; Weibel, 1959; Harrison, 1970). Consequently, the role
of magnetic fields in PopIII star formation processes falls back to the
turbulent dynamo problem, i.e., the amplification of magnetic field through
turbulence (Tan & Blackman, 2004; Schleicher et al., 2010; Schober et al.,
2012; Turk et al., 2012; Latif et al., 2013; Schober et al., 2015).

Magnetic field impacts PopIII stars in two major ways – the fragmentation
processes (and thus the mass and the initial mass function, IMF), and the
rotataion of PopIII stars. Both of the mass and the rotation have major roles
in affecting the stellar evolution path, and thus the subsequent chemical
enrichment processes and the observables that imprinted on the immediate
descendants – low mass extremely metal poor (EMP) stars (Maeder & Meynet,
2012; Yoon et al., 2012).

Several efforts, from both observational constraint and the results of
numerical simulations, have been devoted to understand the mass and IMF of
PopIII stars (observation: (Aoki et al., 2014), Ishigaki et al. (2018);
simulation: Clark et al. (2011), Greif et al. (2011), Hosokawa et al. (2011),
Hirano et al. (2014), Susa et al. (2014), Stacy et al. (2016), Hosokawa et al.
(2016)). Although fragmentation is a common product in the PopIII star forming
environments in numerical simulations and could potential build up the IMF for
PopIII stars, currently only few simulations have explored the role of
magnetic field (Turk et al., 2012; Latif et al., 2013), which is partly
because the current computer power is not yet possible to fully resolve the
dynamo processes and to investigate the saturation level. Observationally, the
peak of IMF is believed to be at $25M_{\odot}$, which is determined through
the chemical abundance of EMP stars, although the full shape of PopIII stars’
IMF is still unclear. The unclear IMF at both the low and high mass end is due
to the fact that the chemical abundance from these two mass ranges might not
be imprinted in the EMP stars.

The other important component, the rotation of PopIII stars, is relatively
unexplored (Stacy et al., 2011; Choplin et al., 2019). Rotation at the end of
the pre-main sequence phase plays an important role in determining the stellar
evolution path, since it governs the mass loss rate and the internal mixing
(Maeder & Meynet, 2012). In present day star formation, the rotation of a star
is in part regulated by magnetic field, although gravitational torque would
also be important for massive stars (Lin et al., 2011; Rosen et al., 2012).
Without an efficient braking mechanism, the star would gain angular momentum
through accretion and ultimately reach nearly its break-up speed. Efficient
magnetic braking requires a large scale, coherent magnetic field. In addition,
the efficiency of magnetic braking depends on the strength of the magnetic
field. As illustrated in Hirano & Bromm (2018), a nearly equipartition large
scale field is preferred for braking PopIII stars.

In the small-scale dynamo process, the growth of magnetic field is mediated
through energy exchange with local turbulent kinetic energy. The small-scale
dynamo proceeds quickly, with a rate comparable to the eddy turn overtime at
the dissipation scale, and is able to reach a saturation level of
$\mathcal{O}(1\%)-\mathcal{O}(10\%)$ compared to the kinetic turbulent energy
(Schober et al., 2015). However, the produced field geometry is like to be
entangled (Schekochihin et al., 2004). The large scale dynamo, on the other
hand, is able to produce a magnetic field beyond the turbulence driving scale
with a super-equipartition saturation level (Brandenburg, 2001). An efficient
large scale dynamo process would be a good candidate for braking the rotation
of PopIII stars, although its efficiency is known to be limited by the
subsequently developed magnetic current helicity.

Although magnetic field might have an interesting role in the PopIII star
forming regions, current numerical simulations are not able to fully resolve
the dynamo processes and thus are not able to study the role of magnetic field
ab initio. The potential magnetic field geometry and magnitude in the early
Universe thus relies on our understandings of dynamo processes from analytic
modelings and from the turbulent statistics with a guide using results of
turbulent box simulations. In this paper, using the results from pure
hydrodynamic simulations, we aim to understand the possibility for a large
scale dynamo to operate in a primoridal accretion flow, and to estimate the
corresponding saturation level. We outline the numerical setup in S2, discuss
the large scale dynamo process in S3, and discuss the results in S4.

##  2 Simulation Setup

The simulation is carried out using GAMER (Schive et al., 2018), which we have
augmented with a cylindrical fluid solver and a cylindrical self-gravity
solver111https://github.com/wt-liao/gamer_popIII_accretion (see Appendix A for
detail). The Navier-Stokes equation is solved on a uniform cylindrical grid
using a Godnov scheme, applying piecewise-linear reconstruction and using an
HLLC Riemann solver for the flux prediction. In order to consistently follow
the chemistry evolution, chemical abundance is evolved through GRACKLE with 9
species – electron, HI, HII, HeI, HeII, HeIII, H-, H2, $H_{2}^{+}$ (Smith et
al., 2017). For molecular hydrogen formation/dissociation at high density, we
use the three body rate from Glover (2008). The chemical abundances are then
used as input to the thermodynamic calculations used in the hydro solver.
Analysis and visualization of the simulations was conducted with Turk et al.
(2011).

###  2.1 Initial and Boundary Condition

We initialize our simulation with a force balanced, axially-symmetric disk
with circular rotation. A 0.2 $M_{\odot}$ central star is sitting at $r=0$,
that locates outside the simulation domain and provides a static gravitational
potential. The initial surface density $\Sigma$ and disk midplane temperature
$T_{0}$ are

$\displaystyle\displaystyle\Sigma(r)=\frac{0.01M_{\odot}}{\pi(100{\rm\
AU})^{2}}\left(\frac{r}{1{\rm\ AU}}\right)^{-1}\,,$ (1)
$\displaystyle\displaystyle T_{0}(r)=2000{\rm\ K}\left(\frac{r}{1{\rm\
AU}}\right)^{-0.15}\,,$ (2)

which are chosen according to previous cosmological simulation results. The
initial density field is then seeded with a random perturbation of a magnitude
between zero to $10\%$ of the local density, and relaxed for three outer
orbits with a constant cooling time $\tau_{\rm cool}=15\times\tau_{\rm orbit}$
to allow for the development of turbulence and to avoid artificial
fragmentation caused by the initial condition (Deng et al., 2017). Following
this initial relaxation, chemistry is enabled, with an initial chemical
abundance determined by local chemical and ionization equilibrium.
Additionally, the equilibrium is chosen such that the local H2 destruction
timescale is equal to the local orbital time. Following chemical
initialization, the system is then relaxed for another outer orbit, during
which the cooling rate is increased linearly from zero to the appropriate
value. To ensure the cooling curve is well-sampled, we also require the
hydrodynamic time step to be smaller than $5\%$ of the cooling time.

The simulation domain covers $[0.5{\rm\ au},50{\rm\ au}]$ and $[-50{\rm\
au},50{\rm\ au}]$ in the radial and vertical direction respectively, and a
full $2\pi$ in the azimuthal direction. The boundary conditions are periodic
in the azimuthal boundary and outflow in the radial boundary. For the vertical
boundary condition, the density and pressure are determined using an
isothermal and force-balanced boundary. An outflow boundary condition is
applied to the velocity fields. The initial resolution is chosen to be
$(N_{r},N_{\theta},N_{z})=(N,4N,2N)$, with $N=256$, resulting in a roughly
cubical cell at $r=25{\rm\ AU}$. The resolution corresponds to at least $8$
cells per Jeans length for a $\rho=10^{-9}{\rm\ g\ cm^{-3}}$ clump at
$10^{3}{\rm\ K}$, which is above the Truelove condition, $4$ cells per Jeans
length, to avoid artificial fragmentation (Truelove et al., 1997). Past work
has suggested that higher resolution is necessary to achieve numerical
convergence and to allow for the operation of dynamo processes (Sur et al.,
2010; Federrath et al., 2011; Turk et al., 2012). The high resolution run
would be possible through the adaptive mesh refinement technique, and will be
performed in a future study.

###  2.2 Optical Depth Approximation

Molecular hydrogen has two primary cooling channels: one is through ro-
vibrational line cooling and the other is through continue cooling triggered
by collision, which is commonly known as collision induced emission (CIE). CIE
cooling is more efficient at high density, and is typically optically thin
when $n\lesssim 10^{16}{\rm\ cm^{-3}}$. In contrast, H2 line cooling becomes
optically thick at $n\sim 10^{10}{\rm\ cm^{-3}}$, and thus requires a model
for the radiative transfer processes to account for the corresponding cooling.

To model for the optical depth of H2 line cooling, we apply a fitting function
adopted for disk geometry (Liao & Turk, 2019). We assume a vertically
isentropic profile within a disk scale height. The optical depth is:

$\displaystyle\tau(\eta)=76\,f_{\rm H_{2}}\left(\frac{\rho_{0}}{10^{-10}{\rm\
g\ cm^{-3}}}\right)\left(\frac{H}{10{\rm\
AU}}\right)\,\widehat{\tau}\left(T_{0},\eta\right)\,,$ (3)

with

$\displaystyle\widehat{\tau}\left(T_{0},\eta\right)=\widehat{\tau}_{0}\,\left(\frac{T_{0}}{10^{3}{\rm\
K}}\right)^{c1}(1-\eta)^{c2}\,.$ (4)

where the coefficients $(\widehat{\tau},c1,c2)=(1.34,\,-0.79,\,2.18)$. In the
above equations, we introduce a self-similar disk coordinate $\eta$. At the
disk midplane, $\eta\rightarrow 0$; $\eta\rightarrow 1$ corresponds to disk
surface. The disk scale height $H=c_{s}/\Omega$, where $c_{s}$ is the sound
speed at the midplane and $\Omega$ is the angular momentum. $T_{0}$ is the
temperature at midplane. To determine the midplane flow state variables, we
utilize the vertically isentropic condition. In addition, in the disk
atmosphere that $z\gg H$, the density barely contributes to the column
density. Photon thus travels freely and optical depth becomes much smaller
than unity. Accordingly, we set optical depth to a small number when $|z|>2r$.

We show surface density, midplane temperature and midplane photon escape
fraction as a function of radius in Fig. 1, and the midplane density-
temperature phase diagram in Fig. 2. Note that photon escape fraction is
overall smaller than $10^{-2}$, which is caused by large absorption
coefficient around $T=10^{3}{\rm\ K}$. We would also like to caution that our
optical depth approximation is similar to MH-Ray-M2 method in Greif (2014),
but including the thermal broadening effect. Since the effect of Doppler shift
is neglected, photon escape fraction is under-estimated. According to Fig 9 in
Greif (2014), the introduced error would be around a factor of two.

Figure 1: Surface density, midplane temperature, and midplane ro-vibrational
photon escape fraction at $t=0$. Solid line shows the azimuthally averaged
values. The color denotes the probability for a cell to be in the
corresponding states, with dark/light indicating high/low probability.  Figure
2: Phase diagram (density-temperature) at mid-plane when $t=0$.

##  3 Large Scale Dynamo

The most important ingredient for a large scale dynamo process is the
$\alpha_{D}$ component, with $\alpha_{D}=\alpha_{K}+\alpha_{M}$ (Pouquet et
al., 1976; Blackman & Field, 2002). In the above expression,
$\alpha_{K}=-\frac{1}{3}\tau_{\rm
eddy}\,\overline{\vec{v^{\prime}}\cdot(\vec{\nabla}\times\vec{v^{\prime}})}$
is due to the kinetic heliecity, and $\alpha_{M}=\frac{1}{3}\rho^{-1}\tau_{\rm
eddy}\,\overline{\vec{j}\cdot\vec{b}}$ is originated from the magnetic current
helicity, $\vec{j}\cdot\vec{b}$. In the above expressions, $\tau_{\rm eddy}$
is the eddy turn over time, $\vec{v^{\prime}}$ and $\vec{b}$ represent the
turbulent velocity and magnetic fields, and
$\vec{j}=\vec{\nabla}\times\vec{b}$. As a result, at the linear stage,
$\alpha_{D}\approx\alpha_{K}$ and kinetic helicity is the key ingredient to
trigger a large scale dynamo. However, the saturation level depends on the
whole $\alpha_{D}$ expression and thus on the amount of magnetic current
helicity retained in the dynamo active region.

In the mean field approximation, a dynamo can be activated from two
eigenmodes. In the absence of a background mean flow, a pure helical
turbulence is able to drive a large scale dynamo; using forced turbulent box
simulations, Brandenburg (2001) demonstrates that helical turbulence is able
to generate a field beyond the turbulence driving scale, and reach a super-
equipartition level. This is known as the $\alpha^{2}$-dynamo, with a growth
rate

$\displaystyle\gamma_{\alpha^{2}}=|\alpha_{D}^{2}k^{2}|^{1/2}-\left(\eta_{t}+\eta\right)k^{2}\,,$
(5)

where $\eta_{t}$ is the turbulent diffusion coefficient, $\eta$ is the
magnetic diffusivity, and $k$ is the corresponding wave number.

In the presence of a background large scale rotating mean flow, rotation tends
to wind up the poloidal field into toroidal components. Subsequently, the
feedback from the toroidal field back to the poloidal field is through the
$\alpha$-effect. This is the well-known $\alpha\Omega$ process. The growth
rate of $\alpha\Omega$ process is

$\displaystyle\gamma_{\alpha\Omega}=\frac{1}{2}|\alpha_{D}q_{d}\Omega
k_{z}|^{1/2}-\left(\eta_{t}+\eta\right)k^{2}\,,$ (6)

where $\Omega$ is the disk angular frequency, $q_{d}=d\ln\Omega/d\ln r$ is the
disk shear parameter, and $k_{z}$ is the wave number along the rotation axis –
z-axis.

###  3.1 Dynamo Number

Large scale dynamo requires a helical turbulent velocity field. In an
accretion disk, it is commonly believed that turbulence is generated through
MRI. However, MRI requires the growth rate to be, at least, comparable to the
non-ideal diffusion processes. In a primordial accretion flow, the required
minimum field strength is at the order of $10^{-4}{\rm\ G}$, which is several
orders of magnitude higher than the field froze in the gas with a primordial
origin from, e.g., the Biermann process (Tan & Blackman, 2004). Numerically,
it has also been found that with an initial weak (but finite) and random
field, MRI cannot be sustained (Bhat et al., 2017).

However, several other pure hydrodynamic instabilities are known to operate in
an accretion disk, e.g., the vertical shear instability (VSI) (Urpin &
Brandenburg, 1998; Nelson et al., 2013; Richard et al., 2016), and the
parametric instability due to the existence of spiral wave (Bae et al., 2016;
Riols et al., 2017). We suspect that both instabilities are active in the
primordial accretion flow, based on the cooling condition, disk geometry and
the massive nature of primordial disks.

Fig 3 demonstrates the turbulent field by showing the $\alpha_{ss}=\rho
v^{\prime}_{r}v^{\prime}_{\theta}/q_{d}\overline{P}$ from the Shakura-Sunyaev
disk model (Shakura & Sunyaev, 1973), at the disk midplane after an evolution
of one outer orbital time. And, $\overline{P}$ is the vertically averaged
pressure weighted by density. We also show the kinetic helical through
$\alpha_{K}$ component for a given azimuthal angle in Fig 4. In deriving the
turbulent field, we assume that the background velocity has reached a qusai-
steady state and use the azimuthally averaged velocity as the mean flow. For
simplicity, we use isotropic helicity in $\alpha_{K}$, and, in Fig 4,
normalize $\alpha_{K}$ with the unit of $1{\rm\ au}$ and the corresponding
orbital time. Since the mean flow approximation is more appropriate in the
body of a disk, we only show the turbulence related fields at $|z|\leqq 2r$.

Figure 3: Shakura–Sunyaev parameter $\alpha_{ss}$ at the disk midplane after
an evolution of one outer orbital time. Primordial accretion flow is turbulent
with $\alpha_{ss}\gtrsim 10^{-3}$.  Figure 4: Turbulent kinetic helicity
$\alpha_{K}$, in the unit of $1{\rm\ au}$ and the corresponding orbital time.
Helical turbulence allows large scale dynamo modes to grow.

To activate an $\alpha^{2}$-dynamo at a given wavenumber $k$, it requires
$D_{\alpha^{2}}=|\alpha_{D}|/\eta_{t}\,k>1$ (cf. Eq 5). In Fig 5, we show
$D_{\alpha^{2}}(r)$ at $k(r)$, the wavenumber corresponding to the local
radius. For $D_{\alpha^{2}}(r)>1$, $\alpha^{2}$-dynamo is possible at a scale
larger than the local radius. In addition, by requiring the growth rate from
Eq. 6 to be positive, one can derive the expression of dynamo number
$D_{\alpha\Omega}=\alpha\Omega/\eta_{t}^{2}k^{3}$. With an assumed disk
profile, the critical dynamo number is usually found to be
$D_{\alpha\Omega,\rm crit}\sim\mathcal{O}(1)-\mathcal{O}(10)$ (Ruzmaikin et
al., 1988). We show the azimuthally averaged absolute dynamo number,
$|D_{\alpha\Omega}|$, in Fig 6. Since $|D_{\alpha\Omega}|\gg 10$, it confirms
the existence of mean field dynamo in a primordial accretion disk.

Figure 5: Azimuthally averaged $|D_{\alpha^{2}}|$, dynamo number of
$\alpha^{2}$-process. For $|D_{\alpha^{2}}|>1$, $\alpha^{2}$-dynamo is
possible at the scales larger than local radius. Figure 6: Azimuthally
averaged $|D_{\alpha\Omega}|$. Since $|D_{\alpha\Omega}|\gg 1$, mean field
dynamo would be efficient.

###  3.2 $Re_{M}$ and $Pr_{M}$

Unlike small scale dynamo that requires a magnetic Reynolds number
$Re_{M}\gtrsim 10^{2}$, large scale dynamo can be activated with
$Re_{M}\gtrsim 10$ (Brandenburg, 2001; Iskakov et al., 2007; Brandenburg et
al., 2012). $Re_{M}$ of a primordial composition has been calculated by, e.g.,
Schober et al. (2012); Nakauchi et al. (2019), and is known to be larger than
$\mathcal{O}(10^{2})$. In the above cited works, the chemical composition and
temperature field are derived from the one-zone model. Here, for completeness,
we also calculate $Re_{M}$ using the result from our simulation.

$Re_{M}$ is a measure between advection and diffusion dissipation of the
magnetic field. It is defined as

$\displaystyle Re_{M}\equiv\frac{c_{s}H_{d}}{\eta}\,,$ (7)

where $c_{s}$ is the local sound speed, $H_{d}$ is the disk scale height with
$H_{d}\sim r$ for a geometrically thick disk, and $\eta$ is the magnetic
diffusivity. In principle, the non-ideal effects, including Ohmic dissipation,
ambipolar diffusion, and Hall effect, could all contribute to increasing
$\eta$. However, ambipolar diffusion and the Hall effect depend on magnetic
field strength and are therefore functions of the dynamo processes. For
simplicity, we consider the contribution from Ohmic dissipation only. We show
$Re_{M}$ in Fig 7, and leave the detail calculation in Appendix B. Indeed, we
found $Re_{M}>10^{2}$ in the midplane of the disk, and beyond
$\mathcal{O}(10^{5})$ at the disk atmosphere. In addition, the fact that
$Re_{M}>10^{4}$ in most part of the disk indicates that both small and large
scale dynamo are active. Also note that the estimation of $Re_{M}$ here should
be a conservative estimation, since we do not include lithium, a primary
electron donor at $T\gtrsim 2000{\rm\ K}$ (Nakauchi et al., 2019), in our
chemistry network.

Figure 7: Magnetic Reynolds number, $Re_{M}$ at disk midplane (upper) and at a
given azimuthal angle (lower). $Re_{M}\gg 10$ allows large scale dynamo to
activate.

Another important parameter is the magnetic Prandtl number, $Pr_{M}=\nu/\eta$,
where $\nu$ is the kinetic viscosity. For $Pr_{M}\gg 1$, small scale dynamo is
easier to be activated and has a higher saturated level, since the stretch and
folding is at the scale above the resistive dissipation scale. The critical
$Re_{M}$ also decreases when $Pr_{M}$ increases (Iskakov et al., 2007).
However, large scale dynamo does not seem to be sensitive to the variation of
$Pr_{M}$ (Brandenburg, 2009). We show $Pr_{M}$ in Fig 8. The small $Pr_{M}$ at
the disk midplane is because that a high density leads to a small kinetic
viscosity and a large magnetic diffusivity. Since $Pr_{M}\ll 1$, magnetic
field is likely saturate at a level of $\mathcal{O}(10^{-2})$ if the
amplification is through small scale dynamo only (Schober et al., 2015).

Figure 8: Magnetic Prandtl number, $Pr_{M}$ at disk midplane (upper) and at a
given azimuthal angle (lower). $Pr_{M}\ll 1$ indicates that, through purely
small scale dynamo, magnetic field is likely saturated at a level of
$\mathcal{O}(10^{-2})$.

###  3.3 Quenching and Saturation

The saturation level of large scale dynamo processes depends on the amount of
magnetic current helicity retained in the dynamo active region. During the
large scale dynamo processes, $\alpha_{M}$ develops together with the
amplification of magnetic field. Ultimately, it brings the
$\alpha_{D}\rightarrow 0$, and quenches the large scale dynamo processes. A
resolution of this is through magnetic helicity flux that brings away small
scale $\alpha_{M}$ and rejuvenates $\alpha_{D}$ (Blackman & Field, 2000).

Subramanian & Brandenburg (2006) provides an expression for the evolution of
magnetic helicity density ($h_{M}$) due to the helicity flux
($\mathcal{\vec{F}}_{h}$):

$\displaystyle\frac{\partial h_{M}}{\partial
t}+\vec{\nabla}\cdot\mathcal{\vec{F}}_{h}=-2\mathcal{\vec{E}}\cdot\vec{B}-2\eta\overline{\vec{j}\cdot\vec{b}}\,.$
(8)

In the above equation,
$\mathcal{\vec{E}}=\overline{\vec{v}\times\vec{b}}=\alpha\vec{B}-\eta_{t}\vec{J}$
is the turbulent electromotive force (EMF), and the last equality is valid
under the framework of a mean field dynamo and neglects the higher order
terms. Following Shukurov et al. (2006), Eq 8 can be re-written as

$\displaystyle\frac{\partial\alpha_{M}}{\partial
t}\approx-\frac{2\eta_{t}}{l_{0}^{2}}\left(\alpha_{D}\hat{B}^{2}+\alpha_{M}\left(\hat{k}_{f}^{2}+\hat{\eta}_{t}^{-1}\right)\right)-\vec{\nabla}\cdot\left(\alpha_{M}\vec{U}\right)\,,$
(9)

where $\vec{U}$ represents the large scale flow that brings $\alpha_{M}$ away
from the dynamo active zone, $l_{0}$ is the integral scale of the turbulence,
$\hat{B}^{2}=B^{2}/B_{eq}^{2}$, $B_{eq}=\rho\overline{v^{\prime 2}}$,
$\hat{k_{f}}=k_{f}/k_{0}$, $k_{f}$ is the wavenumber corresponding to the
driving force (i.e., at the scale of the disk scale height), and
$\hat{\eta}_{t}=\eta_{t}/\eta$.

In a primordial accretion disk, the large scale flow could be from the
accretion, which brings some amount of the $\alpha_{M}$ away from the bulk of
the disk – dynamo active zone – to the protostar where strong differential
rotation would create a different dynamo condition. At the linear stage, the
dynamo is governed by the kinetic helicity, and is gradually quenched due to
the development of current helicity. At steady state, Eq 9 becomes

$\displaystyle\hat{B}_{\rm
sat}^{2}\sim|\frac{\hat{\alpha}_{K}}{D_{\alpha\Omega}/Re_{M,t}^{2}}|\left(\hat{k}_{f}^{2}+\hat{k}_{f}^{2}\hat{U}Re_{M,t}+\hat{\eta}_{t}^{-1}\right)\,,$
(10)

where $\hat{\alpha}_{i}=\alpha_{i}/c_{s}$, $\hat{U}=U/c_{s}$, and
$Re_{M,t}=c_{s}H_{d}/\eta_{t}=Re_{M}/\hat{\eta_{t}}$. We also assume that, at
saturation, $|\alpha_{D}|\ll|\alpha_{K}|$ and $|\alpha_{M}|\sim|\alpha_{K}|$.
Using the mean value from the simulation result at the midplane with
$10<r{\rm\ (au)}<20$, we have $|\hat{\alpha}_{K}|\rightarrow 0.01$,
$\hat{\eta}_{t}\rightarrow 250$, $Re_{M}\rightarrow 4\times 10^{5}$, so
$Re_{M,t}\rightarrow 1.6\times 10^{3}$. Typically,
$\hat{k}_{f}\sim\mathcal{O}(10^{-1})$, and we choose
$(\hat{k}_{f}^{2}+\hat{k}_{f}^{2}\hat{U}Re_{M,t}+\hat{\eta}_{t}^{-1})\rightarrow\hat{\eta}_{t}^{-1}$
for a conservative estimation of $\hat{B}_{\rm sat}$. The above approximation
also assumes a slow accretion with $\hat{U}\ll 1$. With $D_{\alpha\Omega}\sim
D_{\rm crit}=10$ at the saturation, we have $\hat{B}_{\rm sat}\sim 3>1$. The
large scale dynamo is able to amplify the magnetic field to an equipartition
level in a primordial accretion flow. And, the timescale of a large scale
dynamo is the dissipation time, $\tau_{\alpha\Omega}\sim
H_{d}^{2}/\eta_{t}\sim Re_{M,t}\Omega^{-1}$.

##  4 Discussion & Conclusion

Based on the turbulent velocity field from pure hydrodynamic simulaion
results, we investigate the likelihood of observing large scale dynamo
processes in a primoridal accretion flow. We demonstrate that, as the
accretion flow is turbulent with $\alpha_{ss}\gtrsim 10^{-3}$, and has
$D_{\alpha\Omega}\gg 10$, a fast, large scale dynamo can result during the
linear stage of evolution.

In the presence of a large-scale dynamo, the saturation level depends on the
quantity of current helicity retained in the dynamo-active regions. Under a
slow accretion, with $\hat{U}\ll 1$ that may be able to prevent the dynamo
from being completely quenched, we estimate that the saturation level would
have $B/B_{\rm eq}\sim 3$. This estimation provides a point of view at the
saturation when $|\hat{\alpha}_{D}|\ll 1$ but $\hat{\alpha}_{D}\neq 0$. We
note that for an accretion rate of $\dot{M}\sim 10^{-3}M_{\odot}{\rm
yr^{-1}}$, $\hat{U}\sim\mathcal{O}(10^{-1})$; as a result, the helicity flux
term in Eq. 8 may quickly re-distribute $\alpha_{M}$. A self-consistent
calculation that includes the dynamo equations is necessary to determine the
evolution of the magnetic field, although our calculation is able to place
constraints.

Our analysis also relies on the assumption that the turbulent statistics
remain the same when magnetic fields are present in the accretion disk. While
this may be true during the kinetic stage of dynamo processes, as the magnetic
field is (by definition) too weak to be dynamically important, if a coherent
field with $B\gtrsim 10^{-4}{\rm\ G}$ is attained, the MRI would also be
active in the accretion disk (Tan & Blackman, 2004). This could then become
the primary source of turbulence, and furthermore, a strong magnetic field
would suppress both MRI and VSI (Urpin & Brandenburg, 1998). If there is no
other active instability in the strong magnetic field limit, the saturation
level is unlikely to exceed equipartition.

If the field saturation level reaches equipartition, it may be able to
efficiently brake the rotation of PopIII stars (Hirano & Bromm, 2018). It
affects the efficiency of rotational mixing inside the star, and thus the
elements that would be produced during the stellar evolution. In addition, the
presence of a strong, coherent field would favor asymmetric supernovae
(Ezzeddine et al., 2019). If the parent dark matter halo ultimately serves as
the host of a supermassive black hole, the strong coherent field would also be
able to produce sychrontron radiation that could explain for the detected
excess radio background (Feng & Holder, 2018), and thus the global
21-centimeter absorption depth at $78{\rm\ MHz}$ (Bowman et al., 2018).
Although the dynamo problem is numerically difficult to solve since the most
efficient dynamo mode locates at the dissipation scale that is beyond the
resolution limit, future numerical studies that include magnetic field would
still be beneficial to understand the linear phase, the impact from field
topology, the helicity flux and the turbulence statistics.

W.-T. L. and M.J.T. are supported in part by the Gordon and Betty Moore
Foundation’s Data-Driven Discovery Initiative through Grant GBMF4561 to
Matthew Turk. W.-T. L. acknowledges government scholarship to study aboard
from the ministry of education of Taiwan. H.S. acknowledges the funding
support from the Jade Mountain Young Scholar Award no. NTU-108V0201, MOST of
Taiwan under the grant no. 108-2112-M-002-023-MY3, and the NTU Core Consortium
project under the grant no. NTU-CC-108L893401. This research is part of the
Blue Waters sustained-petascale computing project, which is supported by the
National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state
of Illinois. Blue Waters is a joint effort of the University of Illinois at
Urbana-Champaign and its National Center for Supercomputing Applications.

##  Appendix A Cylindrical Solver in GAMER

In this section, we demonstrate the accuracy of the cylindrical fluid and
self-gravity solvers, and leave the detail implementations and scaling tests
to a future discussion. The cylindrical fluid solver aims to conserve the
angular momentum, which is similar to the implementation in, e.g., Athena
(Skinner & Ostriker, 2010) and PLUTO (Mignone et al., 2007). All the test
problems presented here use piecewise-linear reconstruction and HLLC Riemann
solver for the flux prediction – the same numerical schemes used in the
production run.

The first test problem demonstrates the capability of angular momentum
conservation by examining the Rayleigh’s criterion. Rayleigh’s criterion
requires the angular momentum to increase outward, i.e.,
$\partial\left(r^{2}\Omega(r)\right)/\partial r>0$. For $\Omega(r)\propto
r^{-n_{\Omega}}$, $n_{\Omega}<2$ is the stability condition. For
$n_{\Omega}>2$, a perturbation of $\Omega$ would trigger angular momentum
exchange between different annuli. In Fig 9, we show the angular momentum flux
$\mathcal{F}_{AM}=\int\rho rv_{r}\delta v_{\theta}\,dV/\int RP\,dV$, where
$\delta v_{\theta}=v_{\theta}-r\Omega_{0}(r)$ and $\Omega_{0}(r)$ is the
background equilibrium angular velocity (Skinner & Ostriker, 2010). In the
test problem, $\rho=100$, $P=1$, and $v_{\theta}(r=1)=1$, and we perturb the
angular velocity $\Omega(r)$. The simulation box has a size of
$[1,5]\times[0,2\pi]\times[0,1]$, and has a resolution of $(64,128,16)$ in
$(r,\theta,z)$ directions. For $n_{\Omega}>2$, the unstable velocity profile
initiates angular momentum exchange, whereas, for $n_{\Omega}<2$, the angular
momentum flux is fluctuated around the perturbation level.

Figure 9: Angular momentum flux, $\mathcal{F}_{AM}$, for a rotating flow with
different equilibrium velocity profiles.

As another test problem, in Fig 10, we show the solutions of the blast wave
problem solved in both Cartesian (shown in gray solid line) and cylindrical
coordinate (shown in blue circles). The values shown in Fig 10 are averaged
over the same radial distance from the blast center. In the test problem, the
total energy and the initial radius of the blast wave are: $E_{\rm blast}=10$
and $r_{\rm blast}=0.1$. The background density is set to be unity and
background energy density is $2.5\times 10^{-5}$. The simulation box has a
size of $[0.5,1.5]\times[0,1]\times[0,1]$ in the $(r,\theta,z)$ or $(x,y,z)$
directions. The resolution of both run are $128^{3}$.

For the self-gravity solver in cylindrical coordinate, we solve the integral
form of the Poisson equation, and use fast Fourier transform in the $\theta$
and $z$ directions. As a test problem, we solve the gravitational potential of
four uniform spheres. Each of them has a density $\rho_{0}=10$ and radius
$R_{0}=3.5$. The centers of the spheres are located at $(4.5,\,0.0,\,4.0)$,
$(4.5,\,0.5\pi,\,4.0)$, $(4.5,\,1.0\pi,\,4.0)$, and $(4.5,\,1.5\pi,\,4.0)$.
The simulation domain is $[0.5,\,8.5]\times[0,\,2\pi]\times[0,\,8]$. The
gravitational potential of each sphere is

$\displaystyle\Phi_{\rm
sphere}(R)=\left\\{\begin{array}[]{ll}\displaystyle-\frac{2\pi
G\rho_{0}}{3}\left(3R_{0}^{2}-R^{2}\right)&{\rm\ ,if\ }R\leq R_{0}\\\
\displaystyle-\frac{4\pi G\rho_{0}}{3}\frac{R_{0}^{3}}{R}&{\rm\ ,if\
}R>R_{0}\end{array}\right.\,,$ (A3)

and the final resulting potential is the superposition of all spheres. In Fig
11, we show the relative error of the four sphere problem.

Figure 10: The density ($\rho$), pressure ($P$) and kinetic energy ($K.E.$) of
the blast wave test problem solved in both Cartesian coordinate (gray solid
line) and cylindrical coordinate (blue circles).

Figure 11: Relative error (lower is better) of the gravitational potential in
the four sphere test problem with a resolution of $N=256^{3}$. The left panel
shows the relative error map at $z=4$ plane; The right panel shows the map at
$\theta=0$ plane. Both maps cut through the center of the sphere(s).

##  Appendix B Magnetic Reynolds Number

Magnetic Reynolds number, to re-iterate,

$\displaystyle Re_{M}\equiv\frac{c_{s}H_{d}}{\eta}\,,$ (B1)

where $c_{s}$ is the sound speed, $H_{d}$ is the disk scale height,
$\eta\equiv c^{2}/4\pi\sigma$ is the magnetic diffusivity, and $\sigma$ is the
conductivity. By definition (Nicholson, 1983),

$\displaystyle\sigma\equiv\frac{n_{e}c^{2}}{m_{e}\nu_{\rm coll}}\,,$ (B2)

where we have assumed electron is the primary current carrier because
$m_{e}/m_{p}\ll 1$, and $n_{e}$ is the electron number density, $m_{e}$
($m_{p}$) is the electron (proton) mass, $\nu_{\rm coll}$ is the electron
collision frequency. Electron colliding with either another electron/ion and
with a neutral particle would both impede its mobility and ultimately
contribute to the conductivity. Therefore,

$\displaystyle\nu_{\rm coll}=\nu_{\rm coll,\,ei}+\nu_{\rm coll,\,en}\,,$ (B3)

where $\nu_{\rm coll,\,ei}$ and $\nu_{\rm coll,\,en}$ represent the collision
frequency between electron-ion and electron-neutral respectively. We also note
that the collision frequency between electron-electron and electron-ion is at
the same order, so it is sufficient to discuss one of them.

For collision between electron and neutral,

$\displaystyle\nu_{\rm coll,\,en}=n_{n}\left\langle\sigma v\right\rangle\,,$
(B4)

where $\left\langle\sigma v\right\rangle$ is the coefficient of momentum
transfer rate, and is taken from Pinto & Galli (2008). For collision between
electron and ion, we have (Nicholson, 1983)

$\displaystyle\nu_{\rm coll,\,ei}=\frac{8\pi
n_{i}e^{4}}{m_{e}^{1/2}\left(k_{B}T\right)^{3/2}}\ln\Lambda\,,$ (B5)

where $n_{i}$ is ion number density, $e$ is the electron charge, $k_{B}$ is
the Boltzmann constant, $\Lambda=n_{e}\lambda_{D}^{3}$ is the plsmaa
parameter, and $\lambda_{D}=\sqrt{k_{B}T/8\pi n_{e}e^{2}}$ is the Debye
length. We have also assumed that the integral of small angle scattering
effect between charged particles is more significant than the large angle
scattering, which is true when $\Lambda\gg 1$.

In Fig 12, we show the collision frequency between different species for a
given ionization fraction ($f_{i}$) and total particle number density ($n$).
For simplicity, we set $n_{e}=n_{i}=f_{i}n$, and set the number density of
each neutral species $n_{s}=n$. For $f_{i}\lesssim 10^{-8}$, corresponding to
the accretion phase, electron-neutral collision dominates. Also note that for
$f_{i}\gtrsim 10^{-5}$, corresponding to the collapsing phase and within the
protostar, electron-ion collision becomes important.

Figure 12: Collision frequency between different species. Left panel shows the
collision frequency at $f_{i}=10^{-8}$, and the right panel is the condition
when $f_{i}=10^{-5}$.

## References

  * Aoki et al. (2014) Aoki, W., Tominaga, N., Beers, T. C., Honda, S., & Lee, Y. S. 2014, Science, 345, 912
  * Bae et al. (2016) Bae, J., Nelson, R. P., Hartmann, L., & Richard, S. 2016, ApJ, 829, 13
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Bhat et al. (2017) Bhat, P., Ebrahimi, F., Blackman, E. G., & Subramanian, K. 2017, MNRAS, 472, 2569
  * Biermann (1950) Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Blackman & Field (2000) Blackman, E. G., & Field, G. B. 2000, ApJ, 534, 984
  * Blackman & Field (2002) —. 2002, Phys. Rev. Lett., 89, 265007
  * Bowman et al. (2018) Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67
  * Brandenburg (2001) Brandenburg, A. 2001, ApJ, 550, 824
  * Brandenburg (2009) —. 2009, ApJ, 697, 1206
  * Brandenburg et al. (2012) Brandenburg, A., Sokoloff, D., & Subramanian, K. 2012, Space Sci. Rev., 169, 123
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
  * Choplin et al. (2019) Choplin, A., Tominaga, N., & Ishigaki, M. N. 2019, arXiv e-prints, arXiv:1910.01366
  * Clark et al. (2011) Clark, P. C., Glover, S. C. O., Smith, R. J., et al. 2011, Science, 331, 1040
  * Deng et al. (2017) Deng, H., Mayer, L., & Meru, F. 2017, ApJ, 847, 43
  * Ezzeddine et al. (2019) Ezzeddine, R., Frebel, A., Roederer, I. U., et al. 2019, ApJ, 876, 97
  * Federrath et al. (2011) Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62
  * Feng & Holder (2018) Feng, C., & Holder, G. 2018, ApJ, 858, L17
  * Glover (2008) Glover, S. 2008, in American Institute of Physics Conference Series, Vol. 990, First Stars III, ed. B. W. O’Shea & A. Heger, 25–29
  * Greif (2014) Greif, T. H. 2014, MNRAS, 444, 1566
  * Greif et al. (2011) Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75
  * Harrison (1970) Harrison, E. R. 1970, MNRAS, 147, 279
  * Hennebelle & Teyssier (2008) Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25
  * Hirano & Bromm (2018) Hirano, S., & Bromm, V. 2018, MNRAS, 476, 3964
  * Hirano et al. (2014) Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60
  * Hosokawa et al. (2016) Hosokawa, T., Hirano, S., Kuiper, R., et al. 2016, ApJ, 824, 119
  * Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250
  * Ishigaki et al. (2018) Ishigaki, M. N., Tominaga, N., Kobayashi, C., & Nomoto, K. 2018, ApJ, 857, 46
  * Iskakov et al. (2007) Iskakov, A. B., Schekochihin, A. A., Cowley, S. C., McWilliams, J. C., & Proctor, M. R. E. 2007, Physical Review Letters, 98, 208501
  * Latif et al. (2013) Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013, MNRAS, 432, 668
  * Lee et al. (2018) Lee, C.-F., Hwang, H.-C., Ching, T.-C., et al. 2018, Nature Communications, 9, 4636
  * Li et al. (2011) Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180
  * Liao & Turk (2019) Liao, W.-T., & Turk, M. 2019, arXiv e-prints, arXiv:1911.00610
  * Lin et al. (2011) Lin, M.-K., Krumholz, M. R., & Kratter, K. M. 2011, MNRAS, 416, 580
  * Maeder & Meynet (2012) Maeder, A., & Meynet, G. 2012, Reviews of Modern Physics, 84, 25
  * Matt & Pudritz (2005) Matt, S., & Pudritz, R. E. 2005, MNRAS, 356, 167
  * Matt et al. (2010) Matt, S. P., Pinzón, G., de la Reza, R., & Greene, T. P. 2010, ApJ, 714, 989
  * Mellon & Li (2008) Mellon, R. R., & Li, Z.-Y. 2008, ApJ, 681, 1356
  * Mellon & Li (2009) —. 2009, ApJ, 698, 922
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mouschovias (1991) Mouschovias, T. C. 1991, ApJ, 373, 169
  * Nakauchi et al. (2019) Nakauchi, D., Omukai, K., & Susa, H. 2019, MNRAS, 488, 1846
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Nicholson (1983) Nicholson, D. R. 1983, Introduction to Plasma Theory
  * Pinto & Galli (2008) Pinto, C., & Galli, D. 2008, A&A, 484, 17
  * Pouquet et al. (1976) Pouquet, A., Frisch, U., & Leorat, J. 1976, Journal of Fluid Mechanics, 77, 321
  * Price & Bate (2008) Price, D. J., & Bate, M. R. 2008, MNRAS, 385, 1820
  * Richard et al. (2016) Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571
  * Riols et al. (2017) Riols, A., Latter, H., & Paardekooper, S.-J. 2017, MNRAS, 471, 317
  * Rosen et al. (2012) Rosen, A. L., Krumholz, M. R., & Ramirez-Ruiz, E. 2012, ApJ, 748, 97
  * Ruzmaikin et al. (1988) Ruzmaikin, A. A., Sokolov, D. D., & Shukurov, A. M. 1988, Magnetic Fields of Galaxies, Vol. 133, doi:10.1007/978-94-009-2835-0
  * Schekochihin et al. (2004) Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276
  * Schive et al. (2018) Schive, H.-Y., ZuHone, J. A., Goldbaum, N. J., et al. 2018, MNRAS, 481, 4815
  * Schleicher et al. (2010) Schleicher, D. R. G., Banerjee, R., Sur, S., et al. 2010, A&A, 522, A115
  * Schober et al. (2012) Schober, J., Schleicher, D., Federrath, C., et al. 2012, ApJ, 754, 99
  * Schober et al. (2015) Schober, J., Schleicher, D. R. G., Federrath, C., Bovino, S., & Klessen, R. S. 2015, Phys. Rev. E, 92, 023010
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shukurov et al. (2006) Shukurov, A., Sokoloff, D., Subramanian, K., & Brand enburg, A. 2006, A&A, 448, L33
  * Skinner & Ostriker (2010) Skinner, M. A., & Ostriker, E. C. 2010, ApJS, 188, 290
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217
  * Stacy et al. (2016) Stacy, A., Bromm, V., & Lee, A. T. 2016, MNRAS, 462, 1307
  * Stacy et al. (2011) Stacy, A., Bromm, V., & Loeb, A. 2011, MNRAS, 413, 543
  * Subramanian & Brandenburg (2006) Subramanian, K., & Brandenburg, A. 2006, ApJ, 648, L71
  * Sur et al. (2010) Sur, S., Schleicher, D. R. G., Banerjee, R., Federrath, C., & Klessen, R. S. 2010, ApJ, 721, L134
  * Susa et al. (2014) Susa, H., Hasegawa, K., & Tominaga, N. 2014, ApJ, 792, 32
  * Tan & Blackman (2004) Tan, J. C., & Blackman, E. G. 2004, ApJ, 603, 401
  * Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179
  * Turk et al. (2012) Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ, 745, 154
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399
  * Velikhov (1959) Velikhov, E. P. 1959, Sov. Phys. JETP, 36, 1398
  * Weibel (1959) Weibel, E. S. 1959, Phys. Rev. Lett., 2, 83
  * Yoon et al. (2012) Yoon, S. C., Dierks, A., & Langer, N. 2012, A&A, 542, A113