# Modelling Sampling Distributions of Test Statistics with Autograd

Ali Al Kadhim  Harrison B. Prosper  Department of Physics, Florida State
University, 600 W College Ave, Tallahassee, 32306, Florida, USA

###### Abstract

Simulation-based inference methods that feature correct conditional coverage
of confidence sets based on observations that have been compressed to a scalar
test statistic require accurate modelling of either the p-value function or
the cumulative distribution function (cdf) of the test statistic. If the model
of the cdf, which is typically a deep neural network, is a function of the
test statistic then the derivative of the neural network with respect to the
test statistic furnishes an approximation of the sampling distribution of the
test statistic. We explore whether this approach to modelling conditional
1-dimensional sampling distributions is a viable alternative to the
probability density-ratio method, also known as the likelihood-ratio trick.
Relatively simple, yet effective, neural network models are used whose
predictive uncertainty is quantified through a variety of methods.

###### Contents

  1. 1 Introduction
  2. 2 Related work
    1. 2.1 Likelihood Ratio Trick
    2. 2.2 LF2I and ALFFI
    3. 2.3 Conformal Inference
    4. 2.4 Multistage Modelling of Neural Networks
    5. 2.5 Bootstrap NNs
  3. 3 Example 1: signal/background or ON/OFF model
    1. 3.1 Modelling the CDF with ALFFI
    2. 3.2 Directly modelling the empirical cdf
    3. 3.3 Multistage Modelling of the CDF
  4. 4 Example 2: SIR Model
    1. 4.1 Directly modelling the empirical cdf
    2. 4.2 Uncertainty quantification using Bayesian neural networks
    3. 4.3 Bootstrap NN
  5. 5 Discussion and conclusion
  6. A ALFFI algorithm for modelling CDF for ON/OFF Example
  7. B Residuals

##  1 Introduction

Automatic differentiation (see, for example, Ref.[1]) has revolutionized
machine learning, permitting the routine application of gradient descent
algorithms to fit to data models of essentially unlimited complexity. The same
technology can be used to take the derivative of these models with respect to
their inputs without the need to explicitly calculate the derivatives [2]. A
potentially useful application of this capability is approximating the
probability density function (pdf), $f(x\mid\boldsymbol{\theta})$, given an
accurate neural network model of the associated conditional cumulative
distribution function (cdf), $F(x\mid\boldsymbol{\theta})$, using the fact
that

$\displaystyle f(x\mid\boldsymbol{\theta})$ $\displaystyle=\frac{\partial
F(x\mid\boldsymbol{\theta})}{\partial x},$ (1)

where $\boldsymbol{\theta}$ are the parameters of the data-generation
mechanism, which we distinguish from the parameters $\mathbf{w}$ of the neural
network model. This paper explores this possibility in the context of
simulation-based frequentist inference [3, 4, 5, 6, 7].

Equation (1) furnishes an approximation of the pdf
$f(x\mid\boldsymbol{\theta})$ whether $x$ is a function of the underlying
observations $D$ only or if $x=\lambda(D;\boldsymbol{\theta})$ is a test
statistic that depends on $D$ as well as on the parameters
$\boldsymbol{\theta}$. Moreover, computing the derivative of the cdf using
autograd to obtain the pdf is exact; autograd does not use finite difference
approximations.

Strictly speaking, Eq. (1) applies only if $x$ is from a continuous set.
However, discrete distributions are frequently encountered in high-energy
physics and other fields, and are often approximated by continuous
distributions through suitable coarse-graining of $x$. Therefore, it is of
interest to explore the degree to which we can approximate the cdf of discrete
distributions with continuous ones and then apply Eq. (1) to obtain smooth
approximations of the probability mass functions (pmf).

Given the pdf $f(x\mid\boldsymbol{\theta})$, where the data $D$ have been
compressed to the 1-dimensional quantity $x$, further inferences can be
performed by treating the pdf as a statistical model for $x$. These include
constructing approximate confidence intervals by profiling [8]
$f(x\mid\boldsymbol{\theta})$ or constructing posterior densities for the
parameters $f(\boldsymbol{\theta}\mid x)$ given a suitable prior. However, if
$x=\lambda$ is a test statistic some care must be exercised because, in
general, a variation of $\boldsymbol{\theta}$ will induce a change in
$\lambda$, which should be taken into account. The key issue is whether a
sufficiently accurate model of the cdf $F(x\mid\boldsymbol{\theta})$ can be
constructed. This is the issue we explore in this paper where the focus is on
test statistics that arise in simulation-based inference. We consider the
classic ON/OFF problem of astronomy [9] and high-energy physics as a benchmark
example and apply the insights gained to the SIR problem [10, 11] in
epidemiology.

This paper is structured as follows. In Sec. 2 we review related work in the
areas of smiulation-based inference and uncertainty quantification in neural
networks. In Sec. 3 we model the pdf from the cdf of a test statistic for
problem where the likelihood is tractable: the prototypical signal/background
problem in high energy physics, which in astronomy is known as the On/Off
problem. In sec. 4 we model the pdf from the cdf of a test statistic for a
problem where the likelihood is not tractable: the SIR model in epidemiology.
In both examples, various techniques are applied to improve the accuracy and
quantify the uncertainty. The paper ends with a discussion and conclusion in
Sec. 5.

##  2 Related work

###  2.1 Likelihood Ratio Trick

It has been known since at least the early 1990s [12, 13, 14, 15] that
probabilities can be modeled with neural networks using the density ratio
method, also known as the likelihood ratio trick [3, 16, 17, 18, 19]. Suppose
one has two data samples,
$\\{(\boldsymbol{x}_{i},\boldsymbol{\theta}_{i})\\}_{A}$ and
$\\{(\boldsymbol{x}_{i},\boldsymbol{\theta}_{i})\\}_{B}$. Sample $A$ comprises
pairs $(\boldsymbol{x}_{i},\boldsymbol{\theta}_{i})$ in which the components
are sampled sequentially: $\boldsymbol{\theta}_{i}\sim\pi$ followed by
$\boldsymbol{x}_{i}\sim G(\boldsymbol{\theta}_{i})$, where $\pi$ is a known
prior and $G$ is a simulator. Sample $B$ differs from $A$ in that
$\boldsymbol{x}_{i}\sim g$, where $g(\boldsymbol{x}\mid\boldsymbol{\theta})$
is a _known_ density that may, or may not, depend on $\boldsymbol{\theta}$. If
sample $A$ is assigned target $t=1$ and $B$ is assigned $t=0$, then a
sufficiently flexible function — fitted by minimizing the cross-entropy loss
averaged over a large balanced training sample — yields an approximation to
the function

$\displaystyle D(\boldsymbol{x},\boldsymbol{\theta})$
$\displaystyle=\frac{f(\boldsymbol{x},\boldsymbol{\theta})}{f(\boldsymbol{x},\boldsymbol{\theta})+g(\boldsymbol{x}\mid\boldsymbol{\theta})\pi(\boldsymbol{\theta})},$
(2)

where $f(\boldsymbol{x},\boldsymbol{\theta})$ is the (generally unknown) joint
density associated with sample $A$. A rearrangement of Eq.(2) leads to the
result

$\displaystyle f(\boldsymbol{x}\mid\boldsymbol{\theta})$
$\displaystyle=\frac{f(\boldsymbol{x},\boldsymbol{\theta})}{\pi(\boldsymbol{\theta})},$
$\displaystyle=g(\boldsymbol{x}\mid\boldsymbol{\theta})\left(\frac{D}{1-D}\right).$
(3)

This method and related methods are available in the Madminer package [5].

###  2.2 LF2I and ALFFI

Our motivation for starting with the cdf (or p-value) is that constructing a
model of it is the key step in the likelihood-free frequentist (LF2I) approach
[6] and in an extension of it (ALFFI) [7] in which the test statistic
$x=\lambda(D;\boldsymbol{\theta})$ is included as an input to the neural
network-based model of the cdf. If the test statistic $\lambda$ is chosen so
that large values of $\lambda$ disfavor the hypothesis
$H:\boldsymbol{\theta}=\boldsymbol{\theta}_{0}$ it follows that a confidence
set, $R(D)$, can be constructed at confidence level $\tau=1-\alpha$ where
$\alpha$ is the miscoverage rate. By definition the set $R(D)$ is all values
of $\boldsymbol{\theta}_{0}$ for which
$F(\lambda\mid\boldsymbol{\theta}_{0})\leq\tau$ given data $D$ [6, 7]. The
fact that $\lambda$ is an input to the model in the ALFFI algorithm presents
an opportunity: taking the derivative of the approximate cdf $\hat{F}$ with
respect to $\lambda$ provides an approximation
$\hat{f}(\lambda\mid\boldsymbol{\theta})$ of the sampling distribution of the
test statistic. In this paper we explore the accuracy with which both the cdf
and pdf can be approximated.

###  2.3 Conformal Inference

Suppose we have $n$ training samples
$\left(\mathbf{x}_{i},y_{i}\right)\in\mathbb{R}^{d}\times\mathbb{R},i\in[1,n]$,
where $x_{i}$ is a $d$-dimensional feature vector and $y_{i}$ is the response
variable (target). Let $f(x)$ denote the regression function (such as a neural
network), which is often fitted by minimizing the average quadratic loss
between the target and the function, in which case
$\hat{f}(x)\approx\mathbb{E}[y\mid x]$. We are interested in predicting a new
response $y_{n+1}$ from a new feature vector $x_{n+1}$. Given a miscoverage
rate $\alpha\in[0,1]$ we wish to build a confidence set $\mathcal{C}_{\alpha}$
with the property that

$\mathbb{P}\left\\{y_{n+1}\in\mathcal{C}_{\alpha}\left(x_{n+1}\right)\right\\}\geq
1-\alpha=\tau.$ (4)

The confidence set is a measure of uncertainty in the prediction of the
function $f(x)$.

Conformal prediction (a.k.a. conformal inference) [20, 21, 22] is a general
procedure for constructing such confidence sets/intervals for any predictive
model (such as a neural network). These sets are valid (i.e., they satisfy Eq.
4) in finite samples without any assumptions about the distribution or the
data other than the latter are exchangeable [20]. Split conformal prediction
[23, 24] achieves this by splitting the $n$ points into a training set and a
calibration set. A regression model $\hat{f}$ is fitted on the training set
and then used to predict on the calibration set. Next a conformity score, for
example, $s_{i}=|\hat{f}(x_{i})-y_{i}|$ is used to assess the agreement
between the calibration’s response variable and the predicted value. Next,
define $\hat{q}_{1-\alpha}$ to be the $\lceil(n+1)(1-\alpha)\rceil/n$ quantile
of the scores $s_{1},...,s_{n}$, where $\lceil\cdot\rceil$ is the ceiling
function. Finally, for a new test feature vector $x_{n+1}$ construct the
conformal interval as

$\mathcal{C}_{\alpha}=[\hat{f}(x_{n+1})-\hat{q}_{1-\alpha},\
\hat{f}(x_{n+1})+\hat{q}_{1-\alpha}].$ (5)

More details on conformal prediction are provided in [25, 26]. We use
conformal inference to construct confidence sets for the cdf and pdf of
$\lambda$ as a way to quantify their accuracy and for potentially correcting
the cdf to arrive at a more accurate cdf model.

###  2.4 Multistage Modelling of Neural Networks

Recently, multi-stage neural networks (MSNN) [27] (see also [28]) has been
proposed as a strategy to approximate the target function of neural networks
with remarkable accuracy, with the prediction errors approaching machine
precision $\mathcal{O}(10^{-16})$ for double floating point numbers. The
method constructs a sequence of neural networks, each fitted to the residuals
from the previous stage. We follow the notation in [27].

Let $x$ be the input features and $u_{g}(x)$ be the target function. The data
$(x,u_{g}(x))$ are used to train a neural network $u_{0}(x)$ to regress
$u_{g}(x)$. The error, or residual, between the neural network and the target
function $e_{1}(x)=u_{g}(x)-u_{0}(x)$ is calculated. The residual data
$(x,e_{1}(x)/\epsilon_{1})$ are used to train a second neural network
$u_{1}(x)$ to regress $e_{1}(x)/\epsilon_{1}$, where $\epsilon_{1}$ is the
root-mean square,

$\epsilon_{1}=\sqrt{\frac{1}{N_{d}}\sum_{i=1}^{N_{d}}\left[e_{1}\left(x\right)\right]^{2}}=\sqrt{\frac{1}{N_{d}}\sum_{i=1}^{N_{d}}\left[u_{g}-u_{0}\left(x\right)\right]^{2}}.$
(6)

The target function is rescaled so that its range is of order one. The
corrected regression function is given by

$u_{c}^{(1)}(x)=u_{0}(x)+\epsilon_{1}u_{1}(x).$ (7)

By fitting a second neural network on the residuals and adding it to the
original network one improves the precision of the regression function. The
algorithm can be continued to train further neural networks on
$\left(x,e_{n}\left(x\right)/\epsilon_{n}\right)$ to reach higher accuracy and
all $(n+1)$ neural networks are combined as follows,

$u_{c}^{(n)}(x)=\sum_{j=0}^{n}\epsilon_{j}u_{j}(x),$ (8)

where $\epsilon_{i}$ refers to the root mean square for the $i$-th neural
network with $\epsilon_{0}=1$. We explore whether the cdf to be modeled is
smooth enough for this method to work (see Sec. 3.3).

###  2.5 Bootstrap NNs

The bootstrap [29] is a statistical method that treats an observed dataset as
if it were a population. This makes it possible, for example, to approximate
the sampling distributions of statistics. Given an estimator $\hat{\theta}$ of
$\theta$ the bootstrap quantifies the uncertainty in $\hat{\theta}$ by
repeatedly sampling the original dataset with replacement to create new
datasets. The uncertainty in $\hat{\theta}$ can be quantified with a measure
of the variability of predictions between different bootstrap datasets [30,
31].

In the context of supervised learning, the samples are given in pairs
$z_{i}=(\mathbf{x}_{i},y_{i})$ from a joint distribution
$p_{\theta}(\mathbf{x},y)$. Suppose we are given a set of features and targets
composed of $n$ examples
$\mathbf{z}=\\{\left(\mathbf{x}_{i},y_{i}\right)\\}_{i=1}^{n}$. One draws a
bootstrap dataset from the original training dataset of the same size as the
original dataset. This is done $K$ times and the same model is fitted to each
of the $K$ bootstrap datasets using the same training protocol. This yields
$K$ neural networks and, therefore, $K$ outputs for a given input. A measure
of the spread of the outputs quantifies the uncertainty in the model output.

##  3 Example 1: signal/background or ON/OFF model

For our first example, we choose a problem that is ubiquitous in high-energy
physics and astronomy: the signal/background problem, also known as the ON/OFF
problem in astronomy [9] and we consider its simplest realization. An
observation is made which consists of counting $N$ events (these are particle
collisions in particle physics or photon counts in astronomy). A second
independent observation is made where no signal is present by design, yielding
$M$ counts. Following [7] the statistical model is taken be a product of two
Poisson probability mass functions (pmf) for the observed data $D=\\{N,M\\}$,

$\displaystyle\mathbb{P}(N,M\mid\mu,\nu)=\mathcal{L}(D;\mu,\nu)$
$\displaystyle=\frac{(\mu+\nu)^{N}\exp(-(\mu+\nu))}{N!}\,\frac{\nu^{M}\exp(-\nu)}{M!},$
(9)

where $\mu$ and $\nu$ are the mean signal and background counts, respectively,
and we use the same likelihood ratio test statistic as in [7]

$\lambda(D;\mu,\nu)=-2\log\left[\frac{\mathcal{L}(D;\mu,\nu)}{\mathcal{L}(D;\hat{\mu},\hat{\nu})}\right],$
(10)

where the maximum likelihood estimate of the signal is $\hat{\mu}=N-M$, which
can be positive or negative. Negative signal estimates are avoided by using
the “non-maximum likelihood estimate”

$\hat{\mu}=\left\\{\begin{array}[]{ll}N-M&\text{ if }\quad N>M\\\
0&\quad\textrm{ otherwise.}\end{array}\right.$ (11)

and

$\hat{\nu}=\left\\{\begin{array}[]{ll}M&\text{ if }\quad\hat{\mu}=N-M\\\
(M+N)/2&\quad\textrm{ otherwise.}\end{array}\right.$ (12)

Since $\mu$ is the _parameter of interest_ and $\nu$ is a _nuisance parameter_
one is generally interested in confidence intervals for $\mu$ regardless of
the true value of $\nu$. However, the sampling distribution of $\lambda$,
$f(\lambda\mid\mu,\nu)$, in general depends on both the parameter of interest
and the nuisance parameter. This remains true even if we were to replace Eq.
10 with the profile likelihood ratio in which the nuisance parameter $\nu$ is
replaced by its conditional estimate. This is because of the exclusion of
negative estimates of the signal. Therefore, we still must contend with the
nuisance parameter.

There are at least two plausible ways one might proceed given an approximation
to the sampling distribution, $f(\lambda\mid\mu,\nu)$: either replace the
nuisance parameter by a plug-in estimate or replace it through a procedure
analogous to profiling but applied to $f(\lambda\mid\mu,\nu)$. However, while
profiling one has to account for the fact that $\lambda(D;\mu,\nu)$ depends on
the nuisance parameter. The viability of using a profiled sampling
distribution for creating approximate confidence intervals for $\mu$ alone
remains to be explored.

###  3.1 Modelling the CDF with ALFFI

ALFFI is an algorithm for approximating the conditional cdf, $F(x\mid\theta)$,
of a scalar random variable $x$ where $\theta$ denotes the parameters of the
statistical model. If $x=\lambda(D;\theta)$ is a test statistic then the cdf
can be used to construct confidence sets for all parameters simultaneously.
Here we explore the accuracy of the ALFFI algorithm for modelling the cdf with
a view to deriving the sampling distribution, $f(\lambda\mid\mu,\nu)$ of the
ON/OFF test statistic by differentiating the cdf with respect to $\lambda$.

Again following [7] we start by sampling $\mu$ and $\nu$ from uniform priors
and sample $n_{i}\sim\text{Poiss}(\mu_{i}+\nu_{i})$ and
$m_{i}\sim\text{Poiss}(\nu_{i})$. At each parameter point, we calculate
$\lambda_{i}=\lambda(n_{i},m_{i}\mid\mu_{i},\nu_{i})$ according to Eq. 10 and
the procedure is repeated to sample $N_{i},M_{i}$ and
$\lambda_{D}=\lambda(N_{i},M_{i}\mid\mu_{i},\nu_{i})$. Finally, the indicator
$Z$, which is unity if $\lambda_{i}\leq\lambda_{D}$ and zero otherwise, is
computed. This results in a training set of size $B$,
$\mathcal{T}=\\{(\mu_{i},\nu_{i},\lambda_{i},Z_{i})\\}_{i=1}^{B}$. The
observation in ALFFI and LF2I is that critical value functions, such as the
cdf $F(\lambda\mid\mu,\nu)=\mathbb{P}(\lambda\leq\lambda_{D}\mid\mu,\nu)$, are
the expectation value $\mathbb{E}(Z\mid\lambda_{D},\mu,\nu)$ of the discrete
random variable $Z$, of which a smooth approximation can be created with a
deep neural network $f(\mathbf{x}_{i},\mathbf{w})$. Critically, this network
must be trained to minimize the mean square error (MSE) loss111In LF2I and
ALFFI, the predicted value is $\mathbb{E}[Z\mid x]=\mathbb{P}(Z=1\mid
x)=\frac{\mathbb{P}(x\mid Z=1)\mathbb{P}(Z=1)}{\mathbb{P}(x\mid
Z=1)+\mathbb{P}(x\mid Z=0)}$. In other words, it behaves like a classifier,
and because of that the cross entropy loss can also be used.,

$L(\mathbf{w})=\frac{1}{B}\sum_{i=1}^{B}\left(y_{i}-f(\mathbf{x}_{i},\mathbf{w})\right)^{2},$
(13)

where $\mathbf{w}$ are the parameters (weights) of the neural network and
$\mathbf{x}_{i}=\\{\mu_{i},\nu_{i},\lambda_{i}\\}$ is a batch of training
data, with targets $y_{i}=Z_{i}$. Details of this algorithm are provided in
Appendix A.

The cdf, using the ALFFI algorithm, was approximated using a fully connected
neural network with 3 input features $\mathbf{x}=\\{\mu,\nu,\lambda\\}$, 12
hidden layers with 12 nodes each, and a single output. A sigmoid was used in
the output layer to constrain the output to the unit interval. The activation
function at each hidden node is a SiLU [32]. The SiLU nonlinearity was used in
all the models in this study, since the SiLU allows the model to be smooth and
differentiable multiple times. The network was trained with the NAdam
optimizer [33] with a fixed learning rate of $6\times 10^{-4}$. The training
set is composed of $10^{7}$ examples, which were used in batches of size $512$
for $10^{5}$ iterations, that is, for 5 epochs. Longer training runs were also
performed, but did not yield improvements in the results.

Figure 1: ON/OFF cdfs modelled with ALFFI at different $(\mu,\nu)$ points.
Figure 2: ON/off pdfs obtained by differentiating cdfs modelled with ALFFI at
different $(\mu,\nu)$ points.

Although the ALFFI algorithm yields valid confidence sets for the ON/OFF
problem, it does not yield a sufficiently accurate approximation of the cdf,
and hence of the pdf. This is shown in Fig. 1, where the empirical
cdf222Suppose that we have a sample of observations
$\lambda_{1},...,\lambda_{n}$. The empirical cdf is defined by
$F(\lambda)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}\left(\lambda_{i}\leq\lambda\right)$,
$F$, is displayed together with the approximated cdf
$\hat{F}_{\text{ALFFI}}(\lambda\mid\mu,\nu)$. Taking the derivative of the
approximated cdf yields the ALFFI approximation of the pdf,
$\hat{f}_{\text{ALFFI}}(\lambda)=\frac{\partial\hat{F}_{\text{ALFFI}}}{\partial\lambda}$,
displayed in Fig. 2. We observe that $\hat{f}_{\text{ALFFI}}(\lambda)$
displays sharp fluctuations and turning points, which was also confirmed by
taking a numerical derivative of $\hat{F}_{\text{ALFFI}}$. This is because the
actual slope of the cdf changes dramatically from point to nearby point,
especially at low values of $\lambda$. It can also be seen that the
histogrammed pmf of $\lambda$ at a particular $(\mu,\nu)$ point is a high
frequency function, or one that is noisy and features very sharp peaks. The
fluctuations seen in the pmf also depend on the choice of binning: finer
binning shows larger fluctuations because this is a discrete distribution.
This high-frequency nature of the pmf complicates our attempt to create a
smooth modelling of the pmf with ALFFI.

###  3.2 Directly modelling the empirical cdf

Figure 3: ON/OFF cdfs from modelling the empirical cdf with the latter as the
targets at different $(\mu,\nu)$ points, with the associated $68\%$ confidence
band attained from conformal inference. Figure 4: ON/OFF pdfs obtained by
differentiating the models of the empirical cdf at different $(\mu,\nu)$
points, with the associated $68\%$ confidence band attained from conformal
inference..

Trying to model an intrinsically discrete distribution with a smooth
approximation may be beyond an algorithm such as ALFFI and one might argue it
is in any case a hopeless task. Nevertheless it is useful to check if an
alternative, more direct, approach works better: directly modelling the
empirical cdf. At each $\\{\mu_{i},\nu_{i}\\}$ point, we generate $K=100$
experiments: $\\{N_{i}^{j}\\}_{j=1}^{K}\sim\text{Poiss}(\mu_{i}+\nu_{i})$ and
$\\{M_{i}^{j}\\}_{j=1}^{K}\sim\text{Poiss}(\nu_{i})$, for which we generate
$K$ instances of the test statistic $\\{\lambda_{i}^{j}\\}_{j=1}^{K}$
according to Eq. 10. These data are used to calculate the empirical cdf
$F_{i}(\lambda_{i})$ at every parameter point. The data were then flattened to
result in a training set of size $B=10^{7}$,
$\mathcal{T}=\\{(\mu_{i},\nu_{i},\lambda_{i},F_{i})\\}_{i=1}^{B}$.

The cdf was again approximated with a neural network with 3 input features
$\mathbf{x}=\\{\mu,\nu,\lambda\\}$ but this time with the targets set to
values of the empirical cdf $F$. The architecture consisted of 6 hidden layers
with 12 nodes each and a single output. The activation function at each hidden
node is a SiLU. The network was trained to minimize the mean square error with
the NAdam optimizer with a fixed learning rate of $3\times 10^{-5}$. The
training set is composed of $10^{7}$ examples, which were used in batches of
size $60$ for $10^{6}$ iterations, that is, for 6 epochs.

A much better approximation of the cdf is obtained, as shown in Fig. 3 and
more importantly in the associated pdf. Fig. 4 shows that the jagged
fluctuations in the pdf disappear and we obtain a smooth approximation of the
discrete distribution. Modelling the empirical cdf presents yet another
advantage: the possibility of using the empirical cdf as the true response
variable in the conformal inference algorithm (see Sec. 2.3) and thereby
providing a quantification of the uncertainty in the modelling of the cdf and
pdf. We use conformal inference to calculate a $68\%$ confidence interval in
the cdf space for each value of the test statistic. For the associated pdf, a
coarse-graining of $\lambda$ is performed by histogramming $\lambda$ and the
coarse-grained $f(\lambda)$ was used as the true response variable (that is,
the true pdf). The excellent agreement between the predicted pdf and the
histogrammed pdf is shown in Fig. 4 along with the $68\%$ conformal confidence
interval at each value of $\lambda$.

###  3.3 Multistage Modelling of the CDF

Another approach to improve the accuracy of regression functions is described
in [27], which proposes a multi-stage neural network (MSNN) strategy in which
the network training is divided into multiple stages, where each stage fits a
separate neural network to the residuals from the previous stage (see Sec.
2.4).

Quantile residuals ($\Delta C$): Let $p=\\{F_{1},...,F_{k}\\}$ be a set of
evenly spaced probabilities and let
$\hat{F}=\\{\hat{F}_{1},...,\hat{F}_{n}\\}$ be a set of estimated cdf values
at a given $(\mu,\nu)$ parameter point with $n>>k$. Define $q$ as the
quantiles of $\hat{F}$ associated with the probabilities $p$. If we can model
the correction function $p=C(q)$ that maps the estimated cdfs $q$ to the exact
cdf $p$ with sufficient accuracy then, in principle, the correction function
can be used to improve the estimated cdf. Writing $p=q+\Delta C(q)$, we define
the quantile residual function as $\Delta C=q-p$ and try to model it as a
function of $q$ following the multistage method.

In order to follow the multistage method for constructing a model of the cdf
one needs an accurate model of the residuals. We anticipate that this will be
challenging because the distribution of $\lambda$ is discrete and their values
are not evenly spaced. Indeed, this is borne out in Fig. 9 which shows that
$\Delta C$ is a noisy high frequency function. Furthermore, as can be seen in
Fig. 9, the functions do not seem to vary smoothly with the parameters
$(\mu,\nu)$ (small changes in the parameter values yield large changes in the
residual function).

Nevertheless, we wish to see whether a neural network can produce a smooth
interpolation of this noisy function by implementing the suggestions described
in [27]: 1. Use a sine nonlinearity in the first layer and tanh nonlinearity
in subsequent layers. 2. Use Glorot weight initialization [34] with a constant
factor $\kappa=60$ multiplying the weights in the first layer and no such
factor in subsequent layers. It was observed that the $\kappa$ weight factor
helps to stabilize the training and the sine nonlinearity helps capture the
high frequency. 3. Normalize the target function by the RMS of the outputs of
the target function $\Delta C\rightarrow\Delta C/RMS(\Delta C)$, as in Eq. 6.
Further experiments were performed such as adding or omitting dropout, batch
normalization, weight regularization, and experimenting with different
nonlinear functions such as ReLU, SiLU, tanh, sine, and groupsort [35]. None
of these experiments were able to capture the high frequency behavior of the
residual target function when trained on the residuals of all parameters in
the training data simultaneously. However, we were able to use a simple neural
network to model the residual function when modeled point-by-point in the
parameter space.

The frequency principle [36] states that neural networks first fit low
frequencies and later fit high frequencies during the training process.
However, prior studies [37, 38] have demonstrated that in practice standard
neural networks struggle to fit high frequency functions. We confirm these
results and observe that the neural networks we used struggle to fit the
residual functions even when trained for very long training runs. LSTMs, [39]
on the other, hand are known to outperform standard neural networks in fitting
high frequency and sequential data, where the order of the data is significant
[40]. They have been successfully applied in many areas in sequence learning
and time series forecasting such as financial forecasting [41], handwriting
recognition [42] and speech recognition [43]. In our study, LSTMs achieve
reasonably good predictions of the residuals when trained over all $(\mu,\nu)$
values. Although our LSTM models outperform our multi-layer perceptron (MLP)
models in this task the results are unsuitable for use in the MSNN approach
due to inadequate precision in the resulting predictions.

Cdf residuals ($\Delta r$): An alternative to modelling the quantile residuals
is to model the cdf residuals as a function of $\lambda$, which is closer to
what is done in the MSNN approach. At a particular $(\mu,\nu)$ point, the cdf
residual is simply the difference between the predicted and the empirical cdf:
$\Delta r=\hat{F}-F$. Figure 10 shows that the cdf residuals as a function of
$\lambda$ are a slightly smoother function of the parameter points than is
$\Delta C$. Nevertheless, the MLP and LSTM models still struggle to capture
the high-frequency nature of this residual function. A coarse-graining of
$\lambda$ was also performed to reduced the high-frequency of the residual
function but to no avail.

##  4 Example 2: SIR Model

The lessons learned in the ON/OFF example are applied to modelling the cdf in
a different context: epidemiology. We use the SIR model [10, 11] which treats
an epidemic as a system in which transitions occur between three states or
compartments: susceptible (S), infected (I), and recovered (R). This model has
both a stochastic realization as well as a description in terms of ordinary
differential equations

$\displaystyle\frac{dS}{dt}$ $\displaystyle=-\beta\,SI,$
$\displaystyle\frac{dI}{dt}$ $\displaystyle=-\alpha I+\beta\,SI,$
$\displaystyle\frac{R}{dt}$ $\displaystyle=\alpha I,$ (14)

whose solutions approximate the mean number of susceptible $(S(t))$, infected
($I(t)$) and recovered ($R(t)$) individuals as a function of time, $t$. The
model depends on two parameters $\alpha$, the rate of recovery, and $\beta$,
the rate of transmission per infected individual. We fit the model to the data
$D=\\{x_{1},...,x_{n}\\}$ by minimizing the following test statistic [7]

$\displaystyle\lambda(D;\theta)$
$\displaystyle=\frac{1}{50}\sqrt{\frac{1}{N}\sum_{n=1}^{N}\frac{[x_{n}-I_{n}]^{2}}{I_{n}}},$
(15)

where $\theta=\\{\alpha,\beta\\}$, $x_{n}$ is the observed number of infected
individuals at the observation time $t_{n}$ and $I_{n}=I(t_{n},\theta)$ is the
predicted mean infected count obtained by solving the system of ordinary
differential equations, Eq. 14. This example is chosen because it exemplifies
the utility of the methods discussed in this paper for inference with
intractable statistical models. All of our major findings from the ON/OFF
example apply to this example, but here we wish to improve the modelling of
the cdf and employ more techniques for uncertainty quantification.

###  4.1 Directly modelling the empirical cdf

One of the lessons from the ON/OFF problem is that for discrete distributions
the ALFFI algorithm does not yield a sufficiently accurate smooth model of the
cdf, which was confirmed in the SIR example. We therefore directly model the
empirical cdf as a function of the model parameters $(\alpha,\beta)$ and the
test statistic. The training dataset comprises $250$ uniformly sampled
$(\alpha,\beta)$ parameter points, where at each $(\alpha_{i},\beta_{i})$
point, $K=400$ epidemics were simulated yielding a set of test statistics
$\\{\lambda_{i}^{j}\\}_{j=1}^{K}$ and a set of empirical cdfs
$\\{F_{i}^{j}\\}_{j=1}^{K}$. The data were then flattened to result in a
training set of size $B=10^{7}$,
$\mathcal{T}=\\{(\alpha_{i},\beta_{i},\lambda_{i},F_{i})\\}_{i=1}^{B}$.

The MSE loss is known to be sensitive to outliers in the data. Therefore, we
chose to use the more robust Huber loss [44], which is an MSE loss for
relatively small errors and an absolute loss for larger errors. The Huber loss
is given by

$L_{\delta}(\boldsymbol{\omega})=\left\\{\begin{array}[]{ll}\frac{1}{2}(y-f(\mathbf{x},\boldsymbol{\omega}))^{2}&\text{for
}|y-f(\mathbf{x},\boldsymbol{\omega})|\leq\delta\\\
\delta\left(|y-f(\mathbf{x},\boldsymbol{\omega})|-\frac{1}{2}\delta\right),&\text{otherwise.}\end{array}\right.$
(16)

where $\delta$ is a tunable parameter. The Huber loss is both differentiable
everywhere and robust to outliers. We employ models trained with MSE and
models trained with the Huber loss and observe that for relatively short
training runs ($\approx 10^{5}$ iterations) the Huber loss with $\delta=0.7$
significantly outperforms the MSE loss when all other hyperparameters of the
model are held fixed. For larger training runs ($>10^{6}$ iterations),
however, using the Huber loss and MSE loss yields approximately the same
results. In order to reduce the number training iterations, we use the model
with Huber loss for all of our subsequent models.

All of the models for this example approximate the cdf using a fully-connected
neural network with 3 input features $\mathbf{x}=\\{\alpha,\beta,\lambda\\}$,
5 hidden layers with 10 nodes each, and a single output using a sigmoid. The
activation function at each hidden node is a SiLU. The network was trained
with the NAdam optimizer with a fixed learning rate of $3\times 10^{-4}$. The
training set, which comprises $10^{7}$ examples, was used in batches of size
$60$ for $3\times 10^{6}$ iterations, that is, for 18 epochs. As the training
proceeds the network with the lowest validation loss is saved. Figures 5 and 6
show the resulting estimated cdfs and pdfs with the 68% conformal confidence
intervals at different values of $\lambda$ and for different $(\alpha,\beta)$
points. A coarse-graining of $\lambda$ was performed as a way to estimate the
true pdf in order to apply the conformal inference algorithm.

Figure 5: Cdfs for the SIR model obtained from modelling the empirical cdf at
different $(\alpha,\beta)$ points with the associated $68\%$ conformal
confidence band. Figure 6: Pdfs for the SIR model obtained from modelling the
empirical cdf at different $(\alpha,\beta)$ points with the associated $68\%$
conformal confidence band. The spikes around $\lambda\approx 0.15$ arise from
simulated epidemics that die off within a day or so. In a considerably larger
sample of simulated epidemics the spikes would not be visible.

###  4.2 Uncertainty quantification using Bayesian neural networks

Suppose we have a training dataset
$\mathcal{D}=\\{\mathbf{x}_{i},y_{i}\\}_{i=1}^{n}$ (where $\mathbf{x}$ and $y$
are assumed to be drawn from a joint distribution $p(\mathbf{x},y)$). Our goal
is to predict the target variable $y$ given a new input value $\mathbf{x}$.
From a Bayesian perspective, the goal is to approximate the posterior
predictive distribution $p(y\mid\mathbf{x},\mathcal{D})$ given new input
$\mathbf{x}$ and training data $\mathcal{D}$,

$p(y\mid\mathbf{x},\mathcal{D})=\int
p(y\mid\mathbf{x},\mathbf{w})\,p(\mathbf{w}\mid\mathcal{D})\,\mathrm{d}\mathbf{w},$
(17)

which entails a marginalization over the network parameters. By Bayes’
theorem, the posterior density over the neural network parameter space is
given by

$p(\mathbf{w}\mid\mathcal{D})=\frac{p(\mathcal{D}\mid\mathbf{w})p(\mathbf{w})}{p(\mathcal{D})},$
(18)

which requires the model evidence,

$p(\mathcal{D})=\int p(\mathcal{D}\mid\mathbf{w})p(\mathbf{w})d\mathbf{w}.$
(19)

Unlike the best-fit approach to neural networks, a Bayesian neural network
(BNN) is the posterior density $p(\mathbf{w}\mid\mathcal{D})$ over the network
parameter space.

The severe bottleneck with BNNs is computing the high-dimensional integrals
Eq. 17 and Eq. 19. In practice, these integrals must be approximated using
various sampling schemes including Hamiltonian Monte Carlo [45] or variational
inference [46]. Furthermore, a prior density $p(\mathbf{w})$ over the neural
network parameter space must be specified. Given a collection of $K$ neural
networks with parameters sampled from the posterior density, the predictive
distribution can be approximated as follows

$\displaystyle p(y\mid\mathbf{x},\mathcal{D})$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}p(y|\mathbf{x},\mathbf{w}_{k}),$ (20)

where, typically, one models the density $p(y|\mathbf{x},\mathbf{w}_{k})$ as

$\displaystyle
p(y|\mathbf{x},\mathbf{w}_{k})\approx\mathcal{N}(y;f(\mathbf{x},\mathbf{w}),\sigma),$
(21)

with $f(\mathbf{x},\mathbf{w})$ the neural network.

Figure 7: Fluctuating NN parameters by sampling from a prior modeled as a
multivariate diagonal normal distribution with a single standard deviation
$\sigma$. The spread of results is seen to be particularly sensitive to the
choice of $\sigma$. We speculate that this could be due to the lack of
adaptation of $\sigma$ to the different directions in the neural network
parameter space.

In practice, given the computational burden of sampling from the true
posterior density, we adopt a simpler approach. Let
$\mathbf{w^{*}}\in\mathbb{R}^{k}$ be the best-fit neural network parameters
obtained by our training protocol. We fluctuate the neural network parameters
by sampling
$\mathbf{w}_{1},...,\mathbf{w}_{n}\sim\mathcal{N}(\mathbf{w^{*}},\sigma^{2}\mathbf{I}_{k})$
for different choices of $\sigma$. With $n=100$, the resulting ensemble of
models evaluated at the sampled parameters
$\hat{f}(\mathbf{w}_{1}),...,\hat{f}(\mathbf{w}_{n})$ are shown in Fig. 7.

###  4.3 Bootstrap NN

We follow the bootstrap idea (see Sec. 2.5) to draw datasets with replacement
from the training data, such that each bootstrap sample has the same size as
the original training set. This is done $K=200$ times, resulting in $K=200$
bootstrap training datasets. The same model with the architecture described in
Sec. 4.1 is re-trained for $10^{5}$ iterations, with each model fitted to a
different bootstrap sample. The spread of the resulting models over the
$K=200$ neural networks is displayed in Fig. 8. The envelope of the spread is
observed to be much smaller than for the previous methods. This could be due
to insufficient training of the models. Other possible reasons for this are
outlined below.

Figure 8: Responses of the 200 neural networks each trained on a different
bootstrap training sample. See Sec. 4.3 for details.

Although the bootstrap has been hailed as a satisfactory method for
uncertainty quantification due to its theoretical guarantees [47] it has
limitations. For example, it is known to fail in situations where the data
have a non-trivial dependency structure or when the functionals of interest
are not smooth [48].

Another question that was recently explored is the efficacy of bootstrap in
the high-dimensional regime $n/d<1$ relevant to machine learning, where $d$ is
the dimension of the machine learning model parameter space and $n$ the
training sample size. In [49] it is shown that even for the simple case of
inference in the linear model with unregularized least squares, bootstrap
techniques perform poorly in only moderately high dimensions even with $n>d$.
It was also shown in [50] that resampling methods such as the bootstrap yield
reliable error estimates only in the very low-dimensional regime $n>>d$, and
is fraught with problems in high-dimensional regimes. Therefore, perhaps one
shouldn’t be surprised by the results described above.

##  5 Discussion and conclusion

In this paper we try to answer the question whether one can accurately model
the sampling distribution of a test statistic starting from a neural network
model of the cdf. 333Code to reproduce all of our results is available at
https://github.com/hbprosper/cdf2pdf.. Our motivation for modelling the cdf
$F(\lambda\mid\boldsymbol{\theta})$ as opposed to an _ab initio_ modelling of
$f(\lambda\mid\boldsymbol{\theta})$ is that in general a better precision is
reached approximating integrals of functions rather than the functions
themselves. Furthermore, since the step of taking the derivative with autograd
is exact, the accuracy of the approximated cdf directly translates to the
accuracy of the desired pdf.

It was found that a recently introduced method (ALFFI) that approximates the
cdf as the mean of a certain discrete random variable appears not to be
accurate enough in modelling the cdf for the purpose of deriving its
associated pdf. We caution however that this may be the inevitable consequence
of trying to model the cdf of an intrinsically discrete distribution with a
smooth approximation.

We employ relatively simple yet effective neural network models such that the
amount of fine-tuning that is involved in going to more complex data settings
is minimal. Several neural network uncertainty quantification techniques were
reviewed and implemented and we studied possible corrections of the cdf. Of
the uncertainty quantification methods considered only one is calibrated by
construction, namely, conformal inference. We recommend that it be used as a
simple benchmark with respect to which other uncertainty quantification
methods can be compared and calibrated.

As can be seen in Eq. 4 conformal inference results in confidence sets that
guarantee marginal coverage,
$\mathbb{P}\left\\{y_{n+1}\in\mathcal{C}_{\alpha}\left(x_{n+1}\right)\right\\}\geq
1-\alpha$, that is, the coverage marginalized (integrated) over the population
of which a dataset is a sample. That property is not at all controversial as
coverage is a property of the population from which the dataset is presumed to
have been sampled. But if each element $x$ of the dataset contains a subset
that are parameters $\theta$ sampled from a prior $\pi_{\theta}$ then the
marginalization is also with respect to this prior. When a prior is involved
this form of coverage is weaker than the more sought after conditional
coverage guarantee
$\mathbb{P}\left\\{y_{n+1}\in\mathcal{C}_{\alpha}\left(x_{n+1}\right)\mid\theta\right\\}\geq
1-\alpha$. It was shown in [23] that such distribution-free conditional
coverage is impossible to achieve with a finite sample, which has inspired
studies that try to bridge the gap between marginal and conditional coverage
[51]. In our case, $x=\\{\theta,\lambda\\}$ and we generate calibration sets
at each $\theta$ point, thereby making our conformal confidence sets adaptive,
that is, change with $\theta$. However, they are still not adaptive with
respect to $\lambda$; the conformal confidence intervals are the same width at
each $\lambda$ value.

Another issue is that our conformal prediction sets are adaptive only because
we have calibration data at the desired parameter points. It would be helpful
to have a smooth interpolation of the width of the conformal confidence
interval so that an interval can be computed at points where no calibration
data exist.

Acknowledgements This work is supported in part by the US Department of Energy
under Award No. DE-SC0010102 and in part by the University Research
Association (URA) Visiting Scholars Program Award No. 23-S-01. We are deeply
thankful to Pushpalatha Bhat and Stephen Mrenna for their insightful feedback
on this project and their patience. Special thanks to the Fermilab LHC Physics
Center (LPC) for providing access to resources and facilities that were
instrumental to the completion of this work.

##  Appendix A ALFFI algorithm for modelling CDF for ON/OFF Example

ALFFI algorithm for modelling CDF in the On/Off problem, inspired by [7].

Algorithm 1 Estimate the CDF $\mathbb{C}(\lambda\mid\mu,\nu)$, given the
observed data $\\{N,M\\}$ and the observed test statistic $\lambda$.

1:estimated CDF $\widehat{\mathbb{C}}(\lambda\mid\mu,\nu)$ for all
$\theta=\theta_{0}\in\Theta$

2:set $\mathcal{T^{\prime}}\leftarrow\emptyset$

3:for $i$ in $\\{1,...,B^{\prime}\\}$ do

4: Draw parameter $\mu_{i}\sim\pi_{\mu}=\text{Unif}(0,20)$

5: Draw parameter $\nu_{i}\sim\pi_{\nu}=\text{Unif}(0,20)$

6: Draw $n_{i}\sim\text{Poiss}(\mu_{i}+\nu_{i})$

7: Draw $m_{i}\sim\text{Poiss}(\nu_{i})$

8: Draw parameter $\mu_{i}^{\prime}\sim\pi_{\mu^{\prime}}=\text{Unif}(0,20)$

9: Draw parameter $\nu_{i}^{\prime}\sim\pi_{\nu^{\prime}}=\text{Unif}(0,20)$

10: Draw $N_{i}\sim\text{Poiss}(\mu_{i}^{\prime}+\nu_{i}^{\prime})$

11: Draw $M_{i}\sim\text{Poiss}(\nu_{i}^{\prime})$

12: Compute test statistic under the null
$\lambda_{i}\leftarrow\lambda(n_{i},m_{i}\mid\mu_{i},\nu_{i})$

13: Compute indicator
$Z_{i}\leftarrow\mathbbm{1}(\lambda_{i}\leq\lambda(N_{i},M_{i}\mid\mu_{i},\nu_{i}))$

14:
$\mathcal{T^{\prime}}\leftarrow\mathcal{T^{\prime}}\cup\\{(\mu_{i},\nu_{i},\lambda_{i},Z_{i})\\}$

15:end for

16:Use $\mathcal{T^{\prime}}$ to learn the parameterized function
$\widehat{\mathbb{C}}(\lambda\mid\mu,\nu):=\mathbb{E}[Z\mid\mu,\nu]$ via
regression of $Z_{i}$ on $\\{\mu_{i},\nu_{i},\lambda_{i}\\}$ using mean square
error as the loss function

17:return $\widehat{\mathbb{C}}(\lambda\mid\mu,\nu)$, which estimates
$\mathbb{P}(\lambda_{i}\leq\lambda(N,M\mid\mu,\nu))$

##  Appendix B Residuals

Figure 9: $\Delta C=\hat{q}-p$ residuals as a function of $\hat{q}$ at
different $(\mu,\nu)$ points. Figure 10: $\Delta r=\hat{F}-F$ a a function of
$\lambda$ at different $(\mu,\nu)$ points.

## References

  *   * Baydin et al. [2017] Baydin, A.G., Pearlmutter, B.A., Radul, A.A., Siskind, J.M.: Automatic differentiation in machine learning: a survey. J. Mach. Learn. Res. 18(1), 5595–5637 (2017)
  * Chen et al. [2020] Chen, C.-H., Song, F., Hwang, F.-J., Wu, L.: A probability density function generator based on neural networks. Physica A: Statistical Mechanics and its Applications 541, 123344 (2020) https://doi.org/10.1016/j.physa.2019.123344
  * Cranmer et al. [2016] Cranmer, K., Pavez, J., Louppe, G.: Approximating Likelihood Ratios with Calibrated Discriminative Classifiers (2016)
  * Cranmer et al. [2020] Cranmer, K., Brehmer, J., Louppe, G.: The frontier of simulation-based inference. Proc. Nat. Acad. Sci. 117(48), 30055–30062 (2020) https://doi.org/10.1073/pnas.1912789117 arXiv:1911.01429 [stat.ML]
  * Brehmer et al. [2020] Brehmer, J., Kling, F., Espejo, I., Cranmer, K.: MadMiner: Machine learning-based inference for particle physics. Comput. Softw. Big Sci. 4(1), 3 (2020) https://doi.org/10.1007/s41781-020-0035-2 arXiv:1907.10621 [hep-ph]
  * Dalmasso et al. [2023] Dalmasso, N., Masserano, L., Zhao, D., Izbicki, R., Lee, A.B.: Likelihood-Free Frequentist Inference: Confidence Sets with Correct Conditional Coverage (2023)
  * Al Kadhim et al. [2024] Al Kadhim, A., Prosper, H.B., Prosper, O.F.: Amortized simulation-based frequentist inference for tractable and intractable likelihoods. Mach. Learn. Sci. Tech. 5(1), 015020 (2024) https://doi.org/10.1088/2632-2153/ad218e arXiv:2306.07769 [stat.ME]
  * Algeri et al. [2020] Algeri, S., Aalbers, J., Dundas Morå, K., Conrad, J.: Searching for new phenomena with profile likelihood ratio tests. Nature Rev. Phys. 2(5), 245–252 (2020) https://doi.org/10.1038/s42254-020-0169-5 arXiv:1911.10237 [physics.data-an]
  * Li and Ma [1983] Li, T.-P., Ma, Y.-Q.: Analysis methods for results in gamma-ray astronomy. Astrophys. J. 272, 317–324 (1983) https://doi.org/10.1086/161295
  * Kermack and McKendrick [1927] Kermack, W.O., McKendrick, A.G.: A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 115(772), 700–721 (1927)
  * Allen [2017] Allen, L.J.S.: A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis. Infectious Disease Modelling 2(2), 128–142 (2017) https://doi.org/10.1016/j.idm.2017.03.001
  * Ruck et al. [1990] Ruck, D.W., Rogers, S.K., Kabrisky, M., Oxley, M.E., Suter, B.W.: The multilayer perceptron as an approximation to a bayes optimal discriminant function. IEEE Transactions on Neural Networks 1(4), 296–298 (1990) https://doi.org/10.1109/72.80266
  * Wan [1990] Wan, E.A.: Neural network classification: a bayesian interpretation. IEEE Transactions on Neural Networks 1(4), 303–305 (1990) https://doi.org/10.1109/72.80269
  * Richard and Lippmann [1991] Richard, M.D., Lippmann, R.P.: Neural network classifiers estimate bayesian a posteriori probabilities. Neural Computation 3(4), 461–483 (1991) https://doi.org/10.1162/neco.1991.3.4.461
  * Prosper [1993] Prosper, H.B.: Some Mathematical Comments on Feed-Forward Neural Networks: D0-Note-1606 (1993)
  * Mohamed and Lakshminarayanan [2017] Mohamed, S., Lakshminarayanan, B.: Learning in Implicit Generative Models (2017)
  * Tran et al. [2017] Tran, D., Ranganath, R., Blei, D.M.: Hierarchical implicit models and likelihood-free variational inference. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. NIPS’17, pp. 5529–5539. Curran Associates Inc., Red Hook, NY, USA (2017)
  * Goodfellow et al. [2014] Goodfellow, I.J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative Adversarial Networks (2014)
  * Baldi et al. [2016] Baldi, P., Cranmer, K., Faucett, T., Sadowski, P., Whiteson, D.: Parameterized neural networks for high-energy physics. The European Physical Journal C 76(5), 235 (2016) https://doi.org/10.1140/epjc/s10052-016-4099-4
  * Vovk et al. [1999] Vovk, V., Gammerman, A., Saunders, C.: Machine-learning applications of algorithmic randomness. In: International Conference on Machine Learning (1999). https://api.semanticscholar.org/CorpusID:1636783
  * Vovk et al. [2005] Vovk, V., Gammerman, A., Shafer, G.: Algorithmic Learning in a Random World, (2005). https://doi.org/10.1007/b106715
  * Papadopoulos et al. [2002] Papadopoulos, H., Proedrou, K., Vovk, V., Gammerman, A.: Inductive confidence machines for regression. In: Elomaa, T., Mannila, H., Toivonen, H. (eds.) Machine Learning: ECML 2002, pp. 345–356. Springer, Berlin, Heidelberg (2002)
  * Lei and Wasserman [2014] Lei, J., Wasserman, L.: Distribution-free prediction bands for non-parametric regression. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76 (2014) https://doi.org/10.1111/rssb.12021
  * Lei et al. [2016] Lei, J., G’Sell, M.G., Rinaldo, A., Tibshirani, R.J., Wasserman, L.A.: Distribution-free predictive inference for regression. Journal of the American Statistical Association 113, 1094–1111 (2016)
  * Zaffran et al. [2022] Zaffran, M., Féron, O., Goude, Y., Josse, J., Dieuleveut, A.: Adaptive conformal predictions for time series. In: International Conference on Machine Learning, pp. 25834–25866 (2022). PMLR
  * Angelopoulos and Bates [2021] Angelopoulos, A.N., Bates, S.: A gentle introduction to conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511 (2021)
  * Wang and Lai [2024] Wang, Y., Lai, C.-Y.: Multi-stage neural networks: Function approximator of machine precision. Journal of Computational Physics, 112865 (2024)
  * Aldirany et al. [2023] Aldirany, Z., Cottereau, R., Laforest, M., Prudhomme, S.: Multi-level Neural Networks for Accurate Solutions of Boundary-Value Problems (2023)
  * Efron [1979] Efron, B.: Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics 7(1), 1–26 (1979) https://doi.org/10.1214/aos/1176344552
  * Bishop [2006] Bishop, C.M.: Pattern Recognition and Machine Learning. Springer, ??? (2006)
  * Clarté et al. [2024] Clarté, L., Vandenbroucque, A., Dalle, G., Loureiro, B., Krzakala, F., Zdeborová, L.: Analysis of bootstrap and subsampling in high-dimensional regularized regression. arXiv preprint arXiv:2402.13622 (2024)
  * Elfwing et al. [2018] Elfwing, S., Uchibe, E., Doya, K.: Sigmoid-weighted linear units for neural network function approximation. Neural Networks 107, 3–11 (2018) https://doi.org/10.1016/j.neunet.2017.12.012 . Special issue on deep reinforcement learning
  * Dozat [2016] Dozat, T.: Incorporating nesterov momentum into adam. (2016). https://api.semanticscholar.org/CorpusID:70293087
  * Glorot and Bengio [2010] Glorot, X., Bengio, Y.: Understanding the difficulty of training deep feedforward neural networks. In: Teh, Y.W., Titterington, M. (eds.) Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research, vol. 9, pp. 249–256. PMLR, Chia Laguna Resort, Sardinia, Italy (2010). https://proceedings.mlr.press/v9/glorot10a.html
  * Tanielian and Biau [2021] Tanielian, U., Biau, G.: Approximating lipschitz continuous functions with groupsort neural networks. In: International Conference on Artificial Intelligence and Statistics, pp. 442–450 (2021). PMLR
  * Xu et al. [2019] Xu, Z.-Q.J., Zhang, Y., Luo, T., Xiao, Y., Ma, Z.: Frequency principle: Fourier analysis sheds light on deep neural networks. arXiv preprint arXiv:1901.06523 (2019)
  * Tancik et al. [2020] Tancik, M., Srinivasan, P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoorthi, R., Barron, J., Ng, R.: Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems 33, 7537–7547 (2020)
  * Wang et al. [2021] Wang, S., Wang, H., Perdikaris, P.: On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering 384, 113938 (2021)
  * Hochreiter and Schmidhuber [1997] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural Computation 9(8), 1735–1780 (1997) https://doi.org/10.1162/neco.1997.9.8.1735
  * Lindemann et al. [2021] Lindemann, B., Müller, T., Vietz, H., Jazdi, N.: A survey on long short-term memory networks for time series prediction. Procedia CIRP 99, 650–655 (2021) https://doi.org/10.1016/j.procir.2021.03.088 . 14th CIRP Conference on Intelligent Computation in Manufacturing Engine
  * Bhandari et al. [2022] Bhandari, H.N., Rimal, B., Pokhrel, N.R., Rimal, R., Dahal, K.R., Khatri, R.K.C.: Predicting stock market index using lstm. Machine Learning with Applications 9, 100320 (2022) https://doi.org/10.1016/j.mlwa.2022.100320
  * Graves et al. [2009] Graves, A., Liwicki, M., Fernández, S., Bertolami, R., Bunke, H., Schmidhuber, J.: A novel connectionist system for unconstrained handwriting recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence 31(5), 855–868 (2009) https://doi.org/10.1109/TPAMI.2008.137
  * Graves et al. [2013] Graves, A., Mohamed, A.-r., Hinton, G.: Speech recognition with deep recurrent neural networks. In: 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 6645–6649 (2013). https://doi.org/10.1109/ICASSP.2013.6638947
  * Huber [1964] Huber, P.J.: Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics 35(1), 73–101 (1964) https://doi.org/10.1214/aoms/1177703732
  * Neal [1996] Neal, R.M.: Bayesian Learning for Neural Networks. Springer, Berlin, Heidelberg (1996)
  * Jordan et al. [1999] Jordan, M., Ghahramani, Z., Jaakkola, T., Saul, L.: An introduction to variational methods for graphical models. Machine Learning 37, 183–233 (1999) https://doi.org/10.1023/A:1007665907178
  * Rayner [1988] Rayner, R.K.: Some asymptotic theory for the bootstrap in econometric models. Economics Letters 26(1), 43–47 (1988) https://doi.org/10.1016/0165-1765(88)90049-3
  * Bickel and Freedman [1981] Bickel, P.J., Freedman, D.A.: Some Asymptotic Theory for the Bootstrap. The Annals of Statistics 9(6), 1196–1217 (1981) https://doi.org/10.1214/aos/1176345637
  * Karoui and Purdom [2016] Karoui, N.E., Purdom, E.: Can we trust the bootstrap in high-dimension? (2016)
  * Clarté et al. [2024] Clarté, L., Vandenbroucque, A., Dalle, G., Loureiro, B., Krzakala, F., Zdeborová, L.: Analysis of Bootstrap and Subsampling in High-dimensional Regularized Regression (2024)
  * Gibbs et al. [2023] Gibbs, I., Cherian, J.J., Candès, E.J.: Conformal prediction with conditional guarantees. arXiv preprint arXiv:2305.12616 (2023)


# Simulations of high-redshift [OIII] emitters: Chemical evolution and multi-
line diagnostics

Yurina Nakazato Department of Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo, Tokyo 113-0033, Japan  Naoki Yoshida Department of Physics, The
University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), UT Institute
for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Research Center for the Early Universe, School of Science, The University of
Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Daniel Ceverino Universidad
Autonoma de Madrid, Ciudad Universitaria de Cantoblanco, E-28049 Madrid, Spain
CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, E-28049 Madrid,
Spain

###### Abstract

Recent observations by James Webb Space Telescope discovered a number of high-
redshift galaxies with strong emission lines from doubly ionized oxygen.
Combined with ALMA observations of far-infrared lines, multi-line diagnostics
can be applied to the high-redshift galaxies in order to probe the physical
conditions of the inter-stellar medium. We study the formation and evolution
of galaxies using the FirstLight simulation suite, which provides outputs of
62 high-resolution, zoom-in galaxy simulations. We devise a physical model of
Hii regions and calculate spatially resolved [Oiii] line emission. We show
that massive galaxies with stellar masses of $M_{*}>10^{9}M_{\odot}$
chemically evolve rapidly to $z=9$. Young stellar populations in the star-
forming galaxies boost the [Oiii] line emission, rendering the ratio of line
luminosity to star formation rate larger than that for low-redshift galaxies,
which is consistent with recent observations. Measuring the flux ratios of
rest-frame optical and far-infrared lines allows us to estimate the physical
conditions such as density and metallicity of the star-forming gas in high-
redshift [Oiii] emitters.

##  1 Introduction

Understanding the formation and evolution of the first galaxies is one of the
key scientific goals of new generation telescopes including James Webb Space
Telescope (JWST) and Atacama Large Millimetre/Submillimetre Array (ALMA).
High-redshift galaxies can be detected and identified using strong emission
lines, among which [Oiii] 88$\rm\mu\mathrm{m}~{}$line is thought to be
promising (Inoue et al., 2014). A number of galaxies have been found at $z>7$
by ALMA observations targeting the [Oiii] 88$\mu\mathrm{m}$ line (e.g. Inoue
et al., 2016; Hashimoto et al., 2018), including the most distant galaxy
candidate at $z=13.27$ with a $4\sigma$ [Oiii] 88$\mu\mathrm{m}$ detection
(Harikane et al., 2022). Since the [Oiii] line emission originates from Hii
regions around young massive stars, it can be used to trace the star formation
activities and also the physical properties of the inter-stellar medium (ISM)
in the early galaxies.

JWST is opening a new window into the early universe through its superb
observational capability in near-infrared. For example, JWST Early Research
Observation (ERO) in the lensing field SMACS 0723 already reported three
galaxies confirmed spectroscopically by NIRSpec (Schaerer et al., 2022; Curti
et al., 2022; Heintz et al., 2022). NIRSpec instrument is capable of detecting
and identifying various rest-frame optical lines such as [Oii] 3727Å, [Oiii]
4959Å and [Oiii] 5007Å. The relatively weak [Oiii] 4363Å line has been
detected for all the three galaxies, enabling us to estimate the ISM
metallicity in a direct manner.

Detailed numerical simulations are indispensable to study the physical
conditions of the ISM. There have been several studies focusing on [Oiii]
emission lines from high-z galaxies (Hirschmann et al., 2017; Olsen et al.,
2017; Moriwaki et al., 2018; Katz et al., 2019; Arata et al., 2020; Ceverino
et al., 2021; Pallottini et al., 2022). Moriwaki et al. (2018) use a
cosmological simulation with a large boxsize of 50 Mpc (Shimizu et al., 2016)
to calculate the [Oiii] 88$\mu\mathrm{m}$ line intensities for a few hundred
galaxies with stellar masses of $\sim 10^{8}~{}M_{\odot}$. High-resolution,
zoom-in simulations have also been performed to study the internal structure
of early galaxies (Katz et al., 2019; Arata et al., 2020). For the upcoming
observations conducted by JWST, it is urgently needed to study the population
of high-redshift galaxies with high resolution in a fully cosmological
context. In this Letter, we use the outputs of FirstLight simulation (Ceverino
et al., 2017). The simulation suite is motivated to produce a statistically
significant number of galaxies with very high resolution at the epoch of
reionization. Thanks to the mass and volume complete sample of more than 60
massive galaxies and to the high-resolution of $\sim 20$ pc, we can
investigate the internal structure as well as statistics of the high-redshift
galaxies.

Throughout this Letter, we assume $Z_{\odot}=0.02$ as the solar metallicity
(Anders & Grevesse, 1989).

##  2 method

###  2.1 Cosmological Simulation

We use mass-limited galaxy samples selected from the FirstLight simulation
suite (Ceverino et al., 2017). The simulations are performed with ART code
(Kravtsov et al., 1997; Kravtsov, 2003; Ceverino & Klypin, 2009; Ceverino et
al., 2014), which follows gravitational $N$-body dynamics and Eulerian
hydrodynamics using an adaptive mesh refinement method. Besides the two
processes, the code incorporates astrophysical processes relevant for galaxy
formation. The so-called subgrid physics includes atomic and molecular cooling
of hydrogen and helium, photoionization heating by a cosmological UV
background with partial self-shielding, and star formation and the associated
stellar feedback. Details are described in Ceverino et al. (2017). The
simulations track metals released from SNe-Ia and from SNe-II, using
supernovae yields from Woosley & Weaver (1995).

Our simulated galaxies are hosted by dark matter haloes with maximum circular
velocity ($V_{\mathrm{max}}$) higher than 178 km/s at $z=5$ in a cosmological
volume of 40 $h^{-1}$Mpc on a side. The host haloes are selected in a low-
resolution $N$-body only simulation, for which refined initial conditions are
generated using a standard zoom-in technique (Klypin et al., 2011). The
refinement achieves the dark matter particle mass of $m_{\mathrm{DM}}=8\times
10^{4}~{}M_{\odot}$, the minimum star particle mass of $10^{3}~{}M_{\odot}$,
and the maximum spatial resolution is a few tens proper parsec depending on
the refinement level.

We calculate the stellar mass distribution for the selected 62 massive
galaxies at $z=9,8,7,6$. The maximum stellar mass is 9.5, 9.7, 10.1,
10.7$\times 10^{9}~{}M_{\odot}$, respectively. The sample allows us to study
the evolution of more massive galaxies than in previous simulations, e.g.,
Moriwaki et al. (2018), SERRA simulation (Pallottini et al., 2022) and SÍGAME
simulation (Olsen et al., 2017), and thus is well-suited to compare with
observed massive galaxies by HST, ALMA, and JWST (e.g. Tacchella et al., 2022;
Graziani et al., 2020; Topping et al., 2022; Trussler et al., 2022; Barrufet
et al., 2022; Leethochawalit et al., 2022).

###  2.2 Line emissivity calculation

We generate emission-line maps for our galaxy samples by choosing a region
enclosed by 0.3 times the virial radius of the host halo as same as Mandelker
et al. (2014, 2017). We configure a uniform 3D grid with a side length of 100
pc. We locate the star particles and gas elements within each grid, and store
the mass of stars younger than 10 Myr, the average density of the gas with
$n_{\mathrm{H}}>0.1~{}\mathrm{cm^{-3}}$, and the average metallicity of the
cold/warm gas with $T<5\times 10^{4}~{}\mathrm{K}$. These physical quantities
assigned to the individual grids are then used to compute the line
emissivities in a similar manner to those in Hirschmann et al. (2017);
Moriwaki et al. (2018); Ceverino et al. (2021). We generate a library of
emission lines using CLOUDY (Ferland et al., 2013).The library covers a wide
range of gas metallicity $Z$ and ionization parameter $U$ as given in Table 1.

The library lists the individual line luminosity, $L_{\mathrm{line}}$,
normalized by the H$\beta$ line luminosity calculated with the case-B
approximation (Dopita & Sutherland, 2003), $L^{\mathrm{caseB}}_{{\rm
H}\beta}$, as

$\displaystyle L_{\mathrm{line}}$
$\displaystyle=(1-f_{\mathrm{esc}})\,C_{\mathrm{line}}(Z_{\mathrm{gas}},U,n_{\textsc{Hii}})\,L^{\mathrm{caseB}}_{\mathrm{H\beta}},$
(1) $\displaystyle L^{\mathrm{caseB}}_{\mathrm{H\beta}}$ $\displaystyle=4\pi
j_{\mathrm{H\beta}}V=h\nu_{\mathrm{H\beta}}\left(\frac{\alpha_{\mathrm{H\beta}}^{\mathrm{eff}}}{\alpha_{\mathrm{B}}}\right)Q,$
(2)

where $f_{\mathrm{esc}}$ is the Lyman continuum escape fraction and
$C_{\mathrm{line}}$ is the line luminosity ratio. The $\mathrm{H\beta}$
emission rate per unit volume per unit time per unit solid angle is denoted as
$j_{\mathrm{H\beta}}$, and $\alpha_{\mathrm{H\beta}}^{\mathrm{eff}}$ is an
effective recombination coefficient, $Q$ is the production rate of ionizing
photons from each star particle, and $\alpha_{\mathrm{B}}$ is the case-B
hydrogen recombination coefficient given by

$\alpha_{\mathrm{B}}=2.6\times
10^{-13}\left(\frac{T_{\mathrm{e}}}{10^{4}~{}\mathrm{K}}\right)^{-0.85}~{}\mathrm{cm^{3}s^{-1}}$
(3)

with a constant electron temperature $T_{\mathrm{e}}=10^{4}~{}\mathrm{K}$.

We set $f_{\mathrm{esc}}=0.1$, which is consistent with previous radiative
transfer simulations for massive galaxies with
$M_{\mathrm{halo}}>10^{10-11}M_{\odot}$ (Yajima et al., 2011; Kimm & Cen,
2014; Wise et al., 2014; Paardekooper et al., 2015; Xu et al., 2016). It is
also consistent with recent observational estimates at $z\sim 6-8$ (Castellano
et al., 2017; Robertson et al., 2013). We note that some galaxies have been
reported to have an even higher escape fraction of over 20 percent (e.g.
Marques-Chaves et al., 2022; Vanzella et al., 2016; Fletcher et al., 2019;
Bian & Fan, 2020; Flury et al., 2022) at $z<4$.

Since individual Hii regions are not resolved in our simulations, we resort to
a physical model of the ISM structure to calculate the line emissivities of
Hii regions. We characterize the ISM by the local gas density $n$ and
metallicity $Z$, and also by a volume-averaged ionization parameter

$\langle
U\rangle=\frac{3\alpha_{\mathrm{B}}^{2/3}}{4c}\left(\frac{3Q\epsilon^{2}n_{\textsc{Hii}}}{4\pi}\right)^{1/3}.$
(4)

Our fiducial model assumes a constant gas density $n_{\textsc{Hii}}$ in a
spherical Hii region surrounding a star particle (see, e.g. Panuzzo et al.,
2003; Gutkin et al., 2016). We set the Hii region density
$n_{\textsc{Hii}}=100~{}\mathrm{cm}^{-3}$ (e.g. Osterbrock & Ferland, 2006;
Hirschmann et al., 2017, 2022). We define the volume-filling factor of the gas
as

$\epsilon=\frac{n_{\mathrm{gas,grid}}}{n_{\textsc{Hii}}},$ (5)

where $n_{\mathrm{gas,grid}}$ is the gas number density in each grid. In rare
cases where the volume-averaged gas density exceeds $n_{\textsc{Hii}}$
($\epsilon>1$), we set the filling factor to unity. Note that a larger
$n_{\mathrm{gas,grid}}$ for a fixed $n_{\textsc{Hii}}$ yields a larger filling
factor $\epsilon$. Hence the resulting line emissivity depends only weakly on
the assumed $n_{\textsc{Hii}}$ in our model. We have tested with a few
variations with $n_{\textsc{Hii}}=50,~{}300~{}\mathrm{cm}^{-3}$, and
explicitly checked that our main findings in the following sections are not
sensitively affected by this choice.

$\log_{10}~{}(Z_{\mathrm{gas}}/Z_{\odot})$ | -1.30, -0.70, -0.40, 0., 0.30
---|---
$\log_{10}~{}U$ | -4.0, -3.9, …, -1.1, -1.0
$\log_{10}~{}(n_{\textsc{Hii}}/\mathrm{cm}^{-3}$) | 2.0 (fixed)
Table 1: The parameters used to calculate the line luminosities with CLOUDY.

We compute the production rate of ionizing photons $Q$ of the simulated
galaxies using publicly available tables from the Binary Population and
Spectral Synthesis (BPASS) model (Byrne et al., 2022). Our simulations adopt a
stellar initial mass function represented by broken power laws as

$\displaystyle\begin{aligned} &N(M<M_{\mathrm{max}})\propto\\\
&\int^{M_{1}}_{0.1}\left(\frac{M}{M_{\odot}}\right)^{\alpha_{1}}dM+M_{1}^{\alpha_{1}}\int^{M_{\mathrm{max}}}_{M_{1}}\left(\frac{M}{M_{\odot}}\right)^{\alpha_{2}}dM\end{aligned}$
(6)

with
$\alpha_{1}=-1.3,~{}\alpha_{2}=-2.35,~{}M_{1}=0.5,~{}M_{\mathrm{max}}=300~{}M_{\odot}$
as in Ceverino et al. (2019). We use a grid of 13 values of metallicity, from
$Z=10^{-5}$ to $0.04$, and 50 logarithmic bins in stellar population ages
between 1 Myr and 100 Gyr.

|
---|---
|
Figure 1: The [Oiii] 88 $\mu\mathrm{m}$ luminosity versus SFR for our 62
simulated galaxies at $z=9$ (top left), $z=8$ (top right), $z=7$ (bottom
left), and $z=6$ (bottom right). The solid circles are colored with the gas
metallicity (see the colorbar on the right). For comparison, we show the
[Oiii] -SFR relation derived from observations of local galaxies by De Looze
et al. (2014). Gray points are the observational results of high-z $(z>6)$
galaxies from Hashimoto et al. (2018); Laporte et al. (2017); Tamura et al.
(2019), Inoue et al. (2016)(I16), Hashimoto et al. (2019)(H19), Carniani et
al. (2017)(C17), Wong et al. (2022)(WG22), Witstok et al. (2022)(WT22) and
Harikane et al. (2020).

We re-assign ”fine” ages to star particles in order to mitigate the
discreteness effect caused by our simulation set up. Our simulations produce
new star particles with a fixed time step of $\Delta
t_{\mathrm{SF}}=5~{}\mathrm{Myr}$, and the simulation output timings are not
synchronized with $\Delta t_{\mathrm{SF}}$. In a snapshot, young stars
typically have discretized ages such like $t_{\rm age}=2$ Myr, 7 Myr, etc. The
apparently minor gap in stellar ages causes a large impact when we compute the
line emissivities because the ionization photon production rate quickly
decreases with age. For instance, in the BPASS SED of a single stellar
population that we use, the number of ionizing photons decreases over a factor
of 100 from age 1 Myr to 10 Myr (Xiao et al., 2018). We thus re-assign the
stellar age as follows. We consider star particles younger than 15 Myr, with
stamped ages at $T_{1},~{}T_{2},~{}T_{3}~{}(T_{1}<T_{2}<T_{3})$ Myr. We do
random sampling within each age interval. For instance, to a star with
$T_{1}$, we randomly draw a new age within $[1,T_{1}]$ Myr and re-assign to
it. Finally, we select star particles younger than 10 Myr for our emission
line calculation. We calculate the ionizing photon production rate $Q$ for
each stellar particle using the BPASS table.

We consider stellar atmosphere models with different elemental compositions,
i.e., different values of $[\alpha/\mathrm{Fe}]$. In the BPASS_v2.3 (Byrne et
al., 2022), there are five models with the mass fractions in $\alpha$-elements
relative to iron of $\Delta(\log(\alpha/\mathrm{Fe}))=-0.2,+0.0,+0.2,+0.4$ and
+0.6. For the calculation of [$\alpha$/Fe], the $\alpha$-element abundance is
approximated by the oxygen abundance ($\log N_{\mathrm{O}}$) assuming that a
half of the mass in metals produced by SNII are in the form of oxygen atoms;

$\log N_{\mathrm{O}}=\log(f_{\mathrm{O}}z_{\mathrm{SNII}}/A_{\mathrm{O}}),$
(7)

where $f_{\mathrm{O}},\,z_{\mathrm{SNII}}$ are the fraction of oxygen released
by Type-II SNe, and the mass fraction of metals released from Type-II SNe,
respectively. Here, the atomic weight of oxygen is $A_{\mathrm{O}}=16$ and we
assume $f_{\mathrm{O}}=0.5$ (Woosley & Weaver, 1995). We calculate the iron
abundance ratio considering both contributions from Type-Ia and II SNe as

$N_{\mathrm{Fe}}=\frac{(f_{\mathrm{Fe,Ia}}\,z_{\mathrm{SNIa}}+f_{\mathrm{Fe,II}}\,z_{\mathrm{SNIa}})}{A_{\mathrm{Fe}}},$
(8)

where $z_{\mathrm{SNIa}}$ is the mass fraction of metals released from Type-Ia
SNe and $A_{\mathrm{Fe}}=56$. We set the fractions $f_{\mathrm{Fe,Ia}}=0.5$
(Thielemann et al., 1986) and $f_{\mathrm{Fe,II}}=(0.026,0.033)$ for metal
mass ratio between zero and solar metallicity (Nomoto et al., 2006; Ceverino
et al., 2019), respectively. Finally, [$\alpha$/Fe] is obtained from

$[\alpha/\mathrm{Fe}]=\log N_{\mathrm{O}}-\log
N_{\mathrm{Fe}}-\log(N_{\mathrm{O}}/N_{\mathrm{Fe}})_{\odot},$ (9)

where $(N_{\mathrm{O}}/N_{\mathrm{Fe}})_{\odot}=1.17$ is the solar value of
O/Fe abundance ratio.

##  3 results

We focus on rest-frame sub-millimeter and optical [Oiii] lines from high-
redshift galaxies, which are detected by ALMA and JWST.

###  3.1 $L_{[\textsc{Oiii}]\,}$ vs SFR

Figure 1 shows the [Oiii] 88$\mu\mathrm{m}$ luminosity against star formation
rate (SFR) for our galaxy samples. The color-bar indicates the nebular
metallicity $Z_{\mathrm{neb}}$, which is the line luminosity-weighted gas
metallicity. We compare with the observed local galaxies (De Looze et al.,
2014) and with the observed [Oiii] 88$\mu\mathrm{m}$ luminosities of high-
redshift galaxies (see the caption). At $z=9$ to $z=7$, most of our simulated
galaxies are located above the local galaxy relation (solid line), similar to
the results of Moriwaki et al. (2018); Arata et al. (2020); Pallottini et al.
(2022).

At $z=7-9$, our galaxy samples are distributed around the observed galaxies.
It is interesting that luminous galaxies are already chemically enriched with
$\log(Z/Z_{\odot})\sim-0.5$ at the early epochs. Our simulations predict a
slightly steeper relation at $z=7-9$ than the local relation:

$\displaystyle
L_{[\textsc{Oiii}]\,,{88}}\propto\left(\frac{\mathrm{SFR}}{M_{\odot}\,{\rm
yr}^{-1}}\right)^{0.9-1.2}.$ (10)

Figure 2: Projected gas density (left), averaged ionization parameter
(middle), and [Oiii] 88$\rm\mu\mathrm{m}~{}$distribution (right) for a galaxy
sample FL964 at $z=7$. Each panel shows a region with a side length and depth
of 0.3$R_{\mathrm{vir}}(=7.4~{}\mathrm{kpc})$.

We find three galaxies with $L_{[\textsc{Oiii}]\,}>10^{9}~{}L_{\odot}$ at
$z=7$, which are as bright as several observed galaxies. We study the
structure of one of them (sample FL964) in detail. It has
$M_{\mathrm{gas}}=6.41\times
10^{9}~{}M_{\odot},~{}M_{\mathrm{\star}}=9.96\times 10^{9}~{}M_{\odot}$, and a
specific SFR of $11~{}\mathrm{Gyr}$ at $z=7$. Figure 2 shows the projected
maps of number density of gas, ionization parameter, and [Oiii] 88$\mu$m.
Clearly, regions with high ionization parameters of $\log U\sim-2$ cause high
emissivities, consistent with the observation by Harikane et al. (2020) and
also with recent simulations by Kohandel et al. (2022). The total luminosity
of [Oiii] 5007Å of FL964 is $7.60\times 10^{9}L_{\odot}$, which is about 5
times larger than $L_{[\textsc{Oiii}]\,,88}$.

###  3.2 The mass-metallicity relation

It is important to examine the metallicity evolution of our simulated
galaxies. We study the so-called mass-metallicity relation (MZR) by
calculating the gas-phase metallicity for individual galaxies. Figure 3 shows
the stellar mass-gas phase oxygen abundance relation. We calculate the gas
phase oxygen abundance by adopting the conversion equation of Mandelker et al.
(2014);

$\frac{\mathrm{O}}{\mathrm{H}}=\frac{f_{\mathrm{O}}}{z_{{\rm
SNII}}}{XA_{\mathrm{O}}}.$ (11)

We set the hydrogen mass fraction $X=0.755$ and the other values of
$f_{\mathrm{O}}$ and $A_{\mathrm{O}}$ are the same as those of eq.(7), which
adopts the solar oxygen abundance $12+\log(\mathrm{O/H})=8.9$. We then
calculate the averaged $z_{\mathrm{SNII}}$, weighted by the [Oiii] luminosity
of each grid. This weighting is compatible with observational methods such as
direct method or strong line method, which use oxygen emission lines (e.g.
Bian et al., 2018; Izotov et al., 2019).

We calculate the mass of stars within the region of 0.3 $R_{\mathrm{vir}}$. In
Figure 3, we also plot the MZR for local galaxies from Curti et al. (2020)
(dashed line) and recent JWST observation results of high-redshift galaxies
(Sun et al., 2022; Curti et al., 2022; Langeroodi et al., 2022; Williams et
al., 2022). Curti et al. (2022) estimated metallicities of SMACS field
galaxies by direct method, Sun et al. (2022) adopt strong line calibration by
Bian et al. (2018) using O32, and Langeroodi et al. (2022) and Williams et al.
(2022) adopt strong line method by Izotov et al. (2019).

Our simulated galaxies have similar metallicities (oxygen abundance) and
stellar masses to the observed ones. Note that Figure 3 shows the evolution
for a fixed sample of simulated galaxies, rather than for all the galaxies at
respective epochs. Namely, we select the galaxies at $z=5$ by mass and plot
their progenitors at $z=6-9$. Hence we likely miss low-mass, low-metallicity
galaxies at $z=9$ (see Langan et al. (2020) for the mass-metallicity of low-
mass galaxies in FirstLight). Some galaxies with
$M_{\star}>10^{9}~{}M_{\odot}$ have gas-phase metallicities of
$12+\log\mathrm{(O/H)}\sim 8.5$ even at $z=9$, suggesting that metal-
enrichment can proceed rapidly in the early galaxies.

Figure 3: Gas-phase metallicity versus stellar mass for our galaxy samples
from $z=9$ to $z=6$. The solid lines show the median and the colored bands
indicate the sample dispersion in the range of 5-95%. The dashed line is the
local mass-metallicity relation from Curti et al. (2020). Red, blue, and
purple symbols show the mass and metallicity of the $z>7$ galaxies observed in
SMACS J0723 field (Curti et al., 2022), $z\sim 6$ galaxies observed by JWST/
NIRCam WFSS mode (Sun et al., 2022), and $z=8.1-9.5$ galaxies observed in the
cluster RX J2129.4+0009 field (two galaxies at $z\sim 8.15$ from Langeroodi et
al. (2022) and one at $z=9.51$ from Williams et al. (2022)), respectively.

###  3.3 Far-IR/optical line ratios

It is interesting and timely to explore line-ratio diagnostics using three
[Oiii] lines; 88 $\mu$m, 52 $\rm\mu\mathrm{m}~{}$and 5007Å. The former two
fine-structure lines are observed by ALMA whereas the latter is to be observed
by JWST. Hereafter we denote the line luminosity ratios using the wavelength
such as $R_{5007/88}=L_{5007\mathrm{\mathring{A}}}/L_{88\mu\mathrm{m}}$.
Figure 4 shows $R_{5007/88}$ against $R_{52/88}$ for our simulated galaxies.
We also show the model line ratios obtained by our set of CLOUDY calculations
(Table 1).

The ratio $R_{5007/88}$ is commonly thought to be a sensitive temperature
indicator (e.g. Fujimoto et al., 2022). Interestingly, Figure 4 shows that
$R_{5007/88}$ may also trace the mean gas metallicity of a galaxy. We argue
that it is a model-dependent, indirect indicator because of the complex
dependence of the line emissivities on the relevant physical quantities.
Typically, the oxygen line emissivity increases with increasing oxygen
abundance (metallicity), but there is a critical abundance beyond which the
emissivity decreases because of the temperature decrease of Hii regions owing
to metal line cooling. The critical ”peak” abundance is different for
different lines and thus line ratios vary non-trivially as metallicity
increases.

In Figure 4, we plot local metal-rich galaxies observed with both FIR (Brauher
et al., 2008) and optical emission lines (Moustakas et al., 2006). Most of the
plotted local galaxies have high metallicities with $Z>1Z_{\odot}$ and are
located in the lower portion (low $R_{52/88}$) in the figure. Only NGC 1569,
the left most symbol with $R_{5007/88}=4.9$, has a sub-solar metallicity of
$\log(Z/Z_{\odot})=-0.6$ (Israel, 1988), which is located near the same
metallicity line as our high-redshift galaxy samples. The local planetary
nebulae data from Dinerstein et al. (1985) are also plotted as red stars. It
can be easily estimated that the planetary nebulae have electron densities of
$n_{e}([\textsc{Oiii}]\,)\gtrsim 10^{3}\mathrm{cm^{-3}}$, which are consistent
with those derived from [Oii] line ratios.

The line emissivities and hence the ratios have implicit dependence on
ionization parameter through other quantities such as electron temperature,
but the dependence is weak at $\log U\sim(-3,-2)$. Our simulated galaxies have
generally high ionization parameter with $\log U\simeq-2$ (Figure 2), and thus
we may use $R_{5007/88}$ as a metallicity indicator as well.

In our emission line model (Section 2.2), the Hii regions have a fixed density
of $n_{\textsc{Hii}}=100~{}\mathrm{cm^{-3}}$. Hence our galaxy samples are
populated in the left-upper portion with $R_{52/88}\lesssim 1$. Since
$R_{52/88}$ varies weakly with $Z$ and $U$ (Yang & Lidz, 2020), galaxies with
high $Z$ and high $U$ are distributed toward bottom/right in Figure 4.

Figure 4: Line luminosity ratio $R_{5007/52}$ against $R_{52/88}$. Our
simulated galaxies at $z=7$ are represented by solid circles colored with gas
metallicity. Gray star symbols show the local galaxies from Brauher et al.
(2008); Moustakas et al. (2006) and red ones show the local planetary nebulae
from Dinerstein et al. (1985). The results of CLOUDY calculations are
represented by lines colored with metallicity
($\log(Z/Z_{\odot})=-1.30,-0.70,-0.40,0.0$). Solid, dashed, and dotted lines
are the case of $\log U=-1.5,~{}-2,~{}-3$ respectively. The number densities
of Hii region $\log n_{\mathrm{HII}}[\mathrm{cm^{-3}}]=1,~{}2,~{}3$ are also
marked by ticks from left to right on each CLOUDY line.

##  4 Discussion

In this Letter, we have studied the chemical evolution of early star-forming
galaxies from $z=9$ to $z=6$ by using zoom-in hydrodynamics simulations. We
find that oxygen line emission galaxies with stellar masses of
$M_{\star}=10^{9-9.5}~{}M_{\odot}$ have large ionization parameter of $\log
U=-2$ and metallicity of $\log(Z/Z_{\odot})\sim(-1,-0.5)$. In these galaxies,
metal-enrichment occurs early and quickly over a few hundred million years.

We have examined line diagnostics using [Oiii] 5007Å, 88 $\mu$m, and 52
$\rm\mu\mathrm{m}~{}$for future observation synergies of JWST and ALMA. There
have already been a few interesting observations of high-redshift galaxies.
Killi et al. (2022) use ALMA and detect [Oiii] 52 $\rm\mu\mathrm{m}~{}$line
from a galaxy at $z=7$ for the first time. The derived value of $R_{52/88}\sim
0.7$ is close to our galaxy samples (Figure 4), and indicates a relatively low
electron density of $n_{e}\sim 50-260~{}{\rm cm}^{-3}$. Observations of
SMACS0723-4590 at $z=8.5$ by Fujimoto et al. (2022) show a large line ratio
$R_{5007/88}=15.8$, which is slightly larger than our galaxy samples,
suggesting a low-metallicity of $Z\sim 0.04~{}Z_{\odot}$. Combining
$R_{52/88}$ from future observation will constrain the values of metallicity
and ionization parameter at the same time according to Figure 4. Planned
observations using JWST NIRSpec are targetted to several [Oiii] 88
$\rm\mu\mathrm{m}~{}$emitters (e.g., GO-1740, PI: Harikane, and GO-1840,
PI:Álvarez-Márquez & Hashimoto). Multi-line diagnostics such as those
presented in this Letter holds promise to reveal the physical conditions of
the ISM in the high-redshift galaxies.

Our simulations show rapid chemical evolution at high redshift. The resulting
MZR relation is consistent with up-to-date JWST observations (Figure 3).
Langan et al. (2020) use 300 less massive galaxy samples with $M_{\star}\leq
10^{8.5}~{}M_{\odot}$ at $z=8$ and derive MZR from $z=8$ to $z=5$ (see also a
similar study by Noel et al. (2022)). Our galaxy samples with larger stellar
masses with $M_{\star}=10^{8.5-10.0}~{}M_{\odot}$ show a steeper MZR, which
indicates rapid chemical evolution at the early epoch. It would be highly
interesting to study relatively massive galaxies with JWST observations such
as B14–65666 (Roberts-Borsani et al., 2020), A2744–YD4 (Morishita et al.,
2022), and MACS1149–JD1 (Hashimoto et al., 2018).

There are a few caveats in our emission line model. Most notably we do not
account for dust extinction. Recent ALMA surveys report the existence of a
substantial amount of dust in star-forming galaxies at $z\sim 6-8$ (Fudamoto
et al., 2020; Burgarella et al., 2022; Bakx et al., 2021; Schouws et al.,
2022; Tamura et al., 2019; Inami et al., 2022). Given the importance of
emission line ratios including [Oiii] 5007Å, accurate modeling of dust
extinction may be needed for future studies.

We have studied the statistics of early emission-line galaxies and compared
with recent observations. It will be possible and important to study the
internal structure of galaxies using both JWST observations and numerical
simulations. We have shown in Figure 2 that there are large
variations/fluctuations of line emissivities, metal and density distributions
within a galaxy. Cameron et al. (2022) argue that unresolved variations of the
electron temperature within a galaxy results in a biased estimate when the so-
called $T_{e}$-method is applied. JWST’s NIRSpec IFU can resolve with a pixel
scale of 0.1 [arcsec/pixel]111https://jwst-docs.stsci.edu/jwst-near-infrared-
spectrograph. For our configuration shown in Figure 2, the 7.4 kpc region at
$z=7$ can be resolved with $13\times 13$ pixels. Gravitational lensing
magnification will greatly help resolving further the structure of individual
galaxies. In our future work, we will generate mock two-dimensional maps for
our simulated galaxies with the same resolution of NIRSpec IFU, and will
address how well the physical quantities such as gas density and temperature
distribution can be reconstructed.

##  5 Acknowledgements

We thank Kana Moriwaki and Yuichi Harikane for fruitful discussions. This work
made use of v2.3 of the Binary Population and Spectral Synthesis (BPASS)
models as described in Byrne et al. (2022) and Stanway & Eldridge (2018). The
authors thankfully acknowledges the computer resources at MareNostrum and the
technical support provided by the Barcelona Supercomputing Center (RES-
AECT-2020-3-0019). Numerical analyses were carried out on the analysis servers
at Center for Computational Astrophysics, National Astronomical Observatory of
Japan. YN has been supported by International Graduate Program for Excellence
in Earth-Space Science (IGPEES) of the University of Tokyo. DC is a Ramon-
Cajal Researcher and is supported by the Ministerio de Ciencia, Innovación y
Universidades (MICIU/FEDER) under research grant PGC2018-094975-C21.

## References

  * Anders & Grevesse (1989) Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197, doi: 10.1016/0016-7037(89)90286-X
  * Arata et al. (2020) Arata, S., Yajima, H., Nagamine, K., Abe, M., & Khochfar, S. 2020, MNRAS, 498, 5541, doi: 10.1093/mnras/staa2809
  * Bakx et al. (2021) Bakx, T. J. L. C., Sommovigo, L., Carniani, S., et al. 2021, MNRAS, 508, L58, doi: 10.1093/mnrasl/slab104
  * Barrufet et al. (2022) Barrufet, L., Oesch, P. A., Weibel, A., et al. 2022, arXiv e-prints, arXiv:2207.14733.  https://arxiv.org/abs/2207.14733
  * Bian & Fan (2020) Bian, F., & Fan, X. 2020, MNRAS, 493, L65, doi: 10.1093/mnrasl/slaa007
  * Bian et al. (2018) Bian, F., Kewley, L. J., & Dopita, M. A. 2018, ApJ, 859, 175, doi: 10.3847/1538-4357/aabd74
  * Brauher et al. (2008) Brauher, J. R., Dale, D. A., & Helou, G. 2008, ApJS, 178, 280, doi: 10.1086/590249
  * Burgarella et al. (2022) Burgarella, D., Bogdanoska, J., Nanni, A., et al. 2022, A&A, 664, A73, doi: 10.1051/0004-6361/202142554
  * Byrne et al. (2022) Byrne, C. M., Stanway, E. R., Eldridge, J. J., McSwiney, L., & Townsend, O. T. 2022, MNRAS, 512, 5329, doi: 10.1093/mnras/stac807
  * Cameron et al. (2022) Cameron, A. J., Katz, H., & Rey, M. P. 2022, arXiv e-prints, arXiv:2210.14234.  https://arxiv.org/abs/2210.14234
  * Carniani et al. (2017) Carniani, S., Maiolino, R., Pallottini, A., et al. 2017, A&A, 605, A42, doi: 10.1051/0004-6361/201630366
  * Castellano et al. (2017) Castellano, M., Pentericci, L., Fontana, A., et al. 2017, ApJ, 839, 73, doi: 10.3847/1538-4357/aa696e
  * Ceverino et al. (2017) Ceverino, D., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 470, 2791, doi: 10.1093/mnras/stx1386
  * Ceverino et al. (2021) Ceverino, D., Hirschmann, M., Klessen, R. S., et al. 2021, MNRAS, 504, 4472, doi: 10.1093/mnras/stab1206
  * Ceverino et al. (2019) Ceverino, D., Klessen, R. S., & Glover, S. C. O. 2019, MNRAS, 484, 1366, doi: 10.1093/mnras/stz079
  * Ceverino & Klypin (2009) Ceverino, D., & Klypin, A. 2009, ApJ, 695, 292, doi: 10.1088/0004-637X/695/1/292
  * Ceverino et al. (2014) Ceverino, D., Klypin, A., Klimek, E. S., et al. 2014, MNRAS, 442, 1545, doi: 10.1093/mnras/stu956
  * Curti et al. (2020) Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944, doi: 10.1093/mnras/stz2910
  * Curti et al. (2022) Curti, M., D’Eugenio, F., Carniani, S., et al. 2022, arXiv e-prints, arXiv:2207.12375.  https://arxiv.org/abs/2207.12375
  * De Looze et al. (2014) De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62, doi: 10.1051/0004-6361/201322489
  * Dinerstein et al. (1985) Dinerstein, H. L., Lester, D. F., & Werner, M. W. 1985, ApJ, 291, 561, doi: 10.1086/163096
  * Dopita & Sutherland (2003) Dopita, M. A., & Sutherland, R. S. 2003, Astrophysics of the diffuse universe
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137.  https://arxiv.org/abs/1302.4485
  * Fletcher et al. (2019) Fletcher, T. J., Tang, M., Robertson, B. E., et al. 2019, ApJ, 878, 87, doi: 10.3847/1538-4357/ab2045
  * Flury et al. (2022) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, ApJS, 260, 1, doi: 10.3847/1538-4365/ac5331
  * Fudamoto et al. (2020) Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4, doi: 10.1051/0004-6361/202038163
  * Fujimoto et al. (2022) Fujimoto, S., Ouchi, M., Nakajima, K., et al. 2022, arXiv e-prints, arXiv:2212.06863.  https://arxiv.org/abs/2212.06863
  * Graziani et al. (2020) Graziani, L., Schneider, R., Ginolfi, M., et al. 2020, MNRAS, 494, 1071, doi: 10.1093/mnras/staa796
  * Gutkin et al. (2016) Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757, doi: 10.1093/mnras/stw1716
  * Harikane et al. (2020) Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93, doi: 10.3847/1538-4357/ab94bd
  * Harikane et al. (2022) Harikane, Y., Inoue, A. K., Mawatari, K., et al. 2022, ApJ, 929, 1, doi: 10.3847/1538-4357/ac53a9
  * Hashimoto et al. (2018) Hashimoto, T., Laporte, N., Mawatari, K., et al. 2018, Nature, 557, 392, doi: 10.1038/s41586-018-0117-z
  * Hashimoto et al. (2019) Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71, doi: 10.1093/pasj/psz049
  * Heintz et al. (2022) Heintz, K. E., Giménez-Arteaga, C., Fujimoto, S., et al. 2022, arXiv e-prints, arXiv:2212.06877.  https://arxiv.org/abs/2212.06877
  * Hirschmann et al. (2017) Hirschmann, M., Charlot, S., Feltre, A., et al. 2017, MNRAS, 472, 2468, doi: 10.1093/mnras/stx2180
  * Hirschmann et al. (2022) —. 2022, arXiv e-prints, arXiv:2212.02522.  https://arxiv.org/abs/2212.02522
  * Inami et al. (2022) Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126, doi: 10.1093/mnras/stac1779
  * Inoue et al. (2014) Inoue, A. K., Shimizu, I., Tamura, Y., et al. 2014, ApJ, 780, L18, doi: 10.1088/2041-8205/780/2/L18
  * Inoue et al. (2016) Inoue, A. K., Tamura, Y., Matsuo, H., et al. 2016, Science, 352, 1559, doi: 10.1126/science.aaf0714
  * Israel (1988) Israel, F. P. 1988, A&A, 194, 24
  * Izotov et al. (2019) Izotov, Y. I., Guseva, N. G., Fricke, K. J., & Henkel, C. 2019, A&A, 623, A40, doi: 10.1051/0004-6361/201834768
  * Katz et al. (2019) Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902, doi: 10.1093/mnras/stz1672
  * Killi et al. (2022) Killi, M., Watson, D., Fujimoto, S., et al. 2022, arXiv e-prints, arXiv:2211.01424.  https://arxiv.org/abs/2211.01424
  * Kimm & Cen (2014) Kimm, T., & Cen, R. 2014, ApJ, 788, 121, doi: 10.1088/0004-637X/788/2/121
  * Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102, doi: 10.1088/0004-637X/740/2/102
  * Kohandel et al. (2022) Kohandel, M., Ferrara, A., Pallottini, A., et al. 2022, arXiv e-prints, arXiv:2212.02519.  https://arxiv.org/abs/2212.02519
  * Kravtsov (2003) Kravtsov, A. V. 2003, ApJ, 590, L1, doi: 10.1086/376674
  * Kravtsov et al. (1997) Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73, doi: 10.1086/313015
  * Langan et al. (2020) Langan, I., Ceverino, D., & Finlator, K. 2020, MNRAS, 494, 1988, doi: 10.1093/mnras/staa880
  * Langeroodi et al. (2022) Langeroodi, D., Hjorth, J., Chen, W., et al. 2022, arXiv e-prints, arXiv:2212.02491.  https://arxiv.org/abs/2212.02491
  * Laporte et al. (2017) Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ, 837, L21, doi: 10.3847/2041-8213/aa62aa
  * Leethochawalit et al. (2022) Leethochawalit, N., Trenti, M., Santini, P., et al. 2022, arXiv e-prints, arXiv:2207.11135.  https://arxiv.org/abs/2207.11135
  * Mandelker et al. (2017) Mandelker, N., Dekel, A., Ceverino, D., et al. 2017, MNRAS, 464, 635, doi: 10.1093/mnras/stw2358
  * Mandelker et al. (2014) —. 2014, MNRAS, 443, 3675, doi: 10.1093/mnras/stu1340
  * Marques-Chaves et al. (2022) Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2022, MNRAS, 517, 2972, doi: 10.1093/mnras/stac2893
  * Morishita et al. (2022) Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2022, arXiv e-prints, arXiv:2211.09097.  https://arxiv.org/abs/2211.09097
  * Moriwaki et al. (2018) Moriwaki, K., Yoshida, N., Shimizu, I., et al. 2018, MNRAS, 481, L84, doi: 10.1093/mnrasl/sly167
  * Moustakas et al. (2006) Moustakas, J., Kennicutt, Robert C., J., & Tremonti, C. A. 2006, ApJ, 642, 775, doi: 10.1086/500964
  * Noel et al. (2022) Noel, I., Zhu, H., & Gnedin, N. 2022, arXiv e-prints, arXiv:2210.16750.  https://arxiv.org/abs/2210.16750
  * Nomoto et al. (2006) Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424, doi: 10.1016/j.nuclphysa.2006.05.008
  * Olsen et al. (2017) Olsen, K., Greve, T. R., Narayanan, D., et al. 2017, ApJ, 846, 105, doi: 10.3847/1538-4357/aa86b4
  * Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Paardekooper et al. (2015) Paardekooper, J.-P., Khochfar, S., & Dalla Vecchia, C. 2015, MNRAS, 451, 2544, doi: 10.1093/mnras/stv1114
  * Pallottini et al. (2022) Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621, doi: 10.1093/mnras/stac1281
  * Panuzzo et al. (2003) Panuzzo, P., Bressan, A., Granato, G. L., Silva, L., & Danese, L. 2003, A&A, 409, 99, doi: 10.1051/0004-6361:20031094
  * Roberts-Borsani et al. (2020) Roberts-Borsani, G. W., Ellis, R. S., & Laporte, N. 2020, MNRAS, 497, 3440, doi: 10.1093/mnras/staa2085
  * Robertson et al. (2013) Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71, doi: 10.1088/0004-637X/768/1/71
  * Schaerer et al. (2022) Schaerer, D., Marques-Chaves, R., Oesch, P., et al. 2022, arXiv e-prints, arXiv:2207.10034.  https://arxiv.org/abs/2207.10034
  * Schouws et al. (2022) Schouws, S., Stefanon, M., Bouwens, R., et al. 2022, ApJ, 928, 31, doi: 10.3847/1538-4357/ac4605
  * Shimizu et al. (2016) Shimizu, I., Inoue, A. K., Okamoto, T., & Yoshida, N. 2016, MNRAS, 461, 3563, doi: 10.1093/mnras/stw1423
  * Stanway & Eldridge (2018) Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75, doi: 10.1093/mnras/sty1353
  * Sun et al. (2022) Sun, F., Egami, E., Pirzkal, N., et al. 2022, arXiv e-prints, arXiv:2209.03374.  https://arxiv.org/abs/2209.03374
  * Tacchella et al. (2022) Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022, ApJ, 927, 170, doi: 10.3847/1538-4357/ac4cad
  * Tamura et al. (2019) Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27, doi: 10.3847/1538-4357/ab0374
  * Thielemann et al. (1986) Thielemann, F. K., Nomoto, K., & Yokoi, K. 1986, A&A, 158, 17
  * Topping et al. (2022) Topping, M. W., Stark, D. P., Endsley, R., et al. 2022, MNRAS, 516, 975, doi: 10.1093/mnras/stac2291
  * Trussler et al. (2022) Trussler, J. A. A., Adams, N. J., Conselice, C. J., et al. 2022, arXiv e-prints, arXiv:2207.14265.  https://arxiv.org/abs/2207.14265
  * Vanzella et al. (2016) Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41, doi: 10.3847/0004-637X/825/1/41
  * Williams et al. (2022) Williams, H., Kelly, P. L., Chen, W., et al. 2022, arXiv e-prints, arXiv:2210.15699.  https://arxiv.org/abs/2210.15699
  * Wise et al. (2014) Wise, J. H., Demchenko, V. G., Halicek, M. T., et al. 2014, MNRAS, 442, 2560, doi: 10.1093/mnras/stu979
  * Witstok et al. (2022) Witstok, J., Smit, R., Maiolino, R., et al. 2022, MNRAS, doi: 10.1093/mnras/stac1905
  * Wong et al. (2022) Wong, Y. H. V., Wang, P., Hashimoto, T., et al. 2022, ApJ, 929, 161, doi: 10.3847/1538-4357/ac5cc7
  * Woosley & Weaver (1995) Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181, doi: 10.1086/192237
  * Xiao et al. (2018) Xiao, L., Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 477, 904, doi: 10.1093/mnras/sty646
  * Xu et al. (2016) Xu, H., Wise, J. H., Norman, M. L., Ahn, K., & O’Shea, B. W. 2016, ApJ, 833, 84, doi: 10.3847/1538-4357/833/1/84
  * Yajima et al. (2011) Yajima, H., Choi, J.-H., & Nagamine, K. 2011, MNRAS, 412, 411, doi: 10.1111/j.1365-2966.2010.17920.x
  * Yang & Lidz (2020) Yang, S., & Lidz, A. 2020, MNRAS, 499, 3417, doi: 10.1093/mnras/staa3000