\section*{Acknowledgements}
MAF is supported by a National Science Foundation Graduate Research Fellowship under grant No. DGE-1326120.
MFH is supported by a National Aeronautics and Space Administration FINESST under grant No. ASTRO20-0022.
SB is supported by NSF grant AST-1817256.

Computing resources were provided by Frontera LRAC AST21005.
The authors acknowledge the Frontera computing project at the Texas Advanced Computing Center (TACC) for providing HPC and storage resources that have contributed to the research results reported within this paper.
Frontera is made possible by National Science Foundation award OAC-1818253.
URL: \url{http://www.tacc.utexas.edu}

\section*{Data Availability}
Flux power spectra generated from the low resolution, high resolution, and testing sets are available at \url{https://github.com/mafern/MFEmulatorLyaData}.
HDF5 and plain text (appropriate for multi-fidelity emulation) formats are available.
Select single- and multi-fidelity emulator predictions for the $10$ testing simulations are also available from the same repository.
The spectra underlying the flux power are available upon request.

\bibliographystyle{mnras}

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

The modern and future testbed for cosmology lies in small scales and non-linear structures.
Cosmological analyses exploit observations of these scales to explore questions such as the nature of dark matter, the total neutrino mass, and the thermal history of the intergalactic medium (IGM).
One of the most powerful probes of small scale structure is the Lyman-$\alpha$ forest, a series of absorption features in the spectrum of quasars \citep{1965ApJ...142.1633G, 1998MNRAS.301..478T, 2000ApJ...543....1M, 2001ApJ...552...15H, 2002MNRAS.329..848V, 2006AJ....132..117F, 2006MNRAS.365..231V, 2006ApJS..163...80M}.
Numerical simulations are required to analyze these observations as they probe the non-linear regime.
As these simulations are expensive, cosmologists build emulators \citep{Heitmann:2006,Habib:2007,Heitmann:2009}, which interpolate a summary statistic (in this case the 1D Lyman-$\alpha$~forest flux power spectrum) between simulation outputs at different cosmological parameters.

A recent development in cosmology is the application of multi-fidelity emulators, which allow simulations with different particle loads, and thus costs, to be combined together \citep{2022MNRAS.509.2551H}.
Here we adapt the multi-fidelity emulation technique to the Lyman-$\alpha$~forest 1D flux power spectrum.
Multi-fidelity emulation is especially useful in this context because the Lyman-$\alpha$~forest probes a range of redshifts, and is sensitive to smaller scales, which require higher resolution simulations, at higher redshifts \citep{2009MNRAS.398L..26B}.
The models developed here will allow a single emulator to target the wide range of scales probed by the Lyman-$\alpha$ forest, which would otherwise require a computationally infeasible number of very large simulations \citep{2014JCAP...07..005B}.

The Lyman-$\alpha$ forest is the result of overlapping neutral hydrogen absorption profiles in the spectra from distant luminous quasars, processed through the expansion of the universe \citep{1965ApJ...142.1633G}.
As light travels from the quasar, it passes through neutral hydrogen gas of varying densities.
In the rest frame of those neutral hydrogen islands, light that has been redshifted close to the Lyman-$\alpha$ transition at $1215.67${\AA} will be absorbed and the rest transmitted.
This is repeated as the light continues to intersect more neutral hydrogen islands on its path towards us, the observers.
The result is a quasar transmission spectra containing an overlapping field of absorption features that provides a proxy to the dark matter density along that sightline \citep{1998ApJ...495...44C}.

The densities probed by the Lyman-$\alpha$ forest, from redshift $2-5$, are $\sim 1-100 \ \times$ the cosmological mean density.
At these densities stellar winds and star formation effects are negligible, although black hole feedback is important \citep{2013MNRAS.429.1734V, 2020MNRAS.495.1825C}.
The densities along with the range of scales accessed has made the Lyman-$\alpha$ forest popular in cosmological studies, including: constraining the thermal history of the IGM and thus reionization \citep{2008MNRAS.386.1131B,2014MNRAS.438.2499B, 2016MNRAS.463.2335N, 2019ApJ...872..101B, 2019MNRAS.490.3177W,2021MNRAS.506.4389G, 2021arXiv211100019V}, constraining cosmological parameters including the neutrino mass \citep{2004MNRAS.354..684V, 2005ApJ...635..761M, 2006MNRAS.370L..51V, 2005PhRvD..71j3515S, 2006JCAP...10..014S, 2020JCAP...04..038P, 2021JCAP...03..049G}, and testing alternatives to cold dark matter \citep{2005PhRvD..71f3534V,  2013PhRvD..88d3502V, 2017PhRvD..96b3522I, 2020JCAP...04..038P, 2021MNRAS.502.2356G, 2021PhRvL.126g1302R}.

The Lyman-$\alpha$ forest 1D flux power spectrum is the most commonly used summary statistic for Lyman-$\alpha$~forest spectra.
It probes small scale structure by measuring the two-point Fourier-space correlation between neutral hydrogen absorption within a sightline \citep{1998ApJ...495...44C}.

Current observational measurements of the Lyman-$\alpha$ forest flux power spectrum come from either a lower resolution, larger sample survey (SDSS, \cite{2019JCAP...07..017C}), or various higher resolution, smaller sample surveys \citep{2017MNRAS.466.4332I, 2022MNRAS.509.2842K, 2019MNRAS.489.2536D}.
In \cite{2019JCAP...07..017C}, the flux power spectrum constructed from BOSS and eBOSS spectra accesses redshifts from $z=2.2-4.6$ ($6\%$ \& $18\%$ average uncertainty, respectively) and scales from $k\approx0.001-0.02$ km$^{-1}$ s ($6\%$ \& $14\%$ average uncertainty, respectively).
The small sample, higher resolution surveys generally access a similar redshift range, but shift both the largest and smallest scales to higher $k$.
For example, in \cite{2022MNRAS.509.2842K} (their conservative results), using spectra from multiple surveys (XQ-100, KODIAQ, and SQUAD) they access redshifts from $z=2-4.6$ ($7\%$ \& $27\%$ average uncertainty, respectively) and scales from $k\approx0.005-0.1$ km$^{-1}$ s ($12\%$ \& $9\%$ average uncertainty, respectively).

The Dark Energy Spectroscopic Instrument (DESI) will soon report its first year results.
Ultimately it will increase the number of Lyman-$\alpha$~quasar spectra by a factor of four over SDSS.
This corresponds to  $\sim 50$ quasars per square degree and a total of $7\times10^5$ quasars over the $14,000$ square degree survey footprint \citep{2016arXiv161100036D}.
In addition, DESI is expected to measure the 1D flux power spectrum at smaller scales ($k < 0.035$ km$^{-1}$ s) and higher redshifts ($z>4.6$) than SDSS, achieving order of a few percent accuracy \citep{2022arXiv220307491V}.

Extracting cosmological information from these observations will require simulations which follow the distribution of gas at relevant densities and on relevant scales.
For the Lyman-$\alpha$~forest, box sizes of at least $100$ Mpc h$^{-1}$ and mean particle spacing of $100/3072 \approx 0.03$ Mpc h$^{-1}$ are necessary \citep{2014JCAP...07..005B}.
Earlier work has focused on methods which can reduce the cost of such simulations.
\cite{2014JCAP...07..005B} used a splicing technique to produce high resolution, large volume outputs from three sets of less computationally intensive simulations: low resolution, large volume; high resolution, small volume; and low resolution, small volume.
\cite{2015MNRAS.446.3697L} explored the use of Richardson extrapolation to enhance output resolution, in addition to testing the splicing technique.

Parameter inference tasks, such as a direct Markov Chain Monte Carlo analyses, require $\sim10^5-10^6$ model evaluations, indicating the number of simulations required by a naive approach.
Even using techniques such as splicing, this is computationally infeasible.
However, using a significantly reduced number of simulations ($\sim30$), an emulator can be constructed that effectively interpolates between this smaller set of simulations.
In addition to the Lyman-$\alpha$~forest, emulators have been used extensively in cosmology for studying: the matter power spectrum \citep{Heitmann:2009, Heitmann:2014, Lawrence:2017, Giblin:2019, Euclid:2021, Arico:2021, Giri:2021}, weak lensing \citep{Harnois:2019, Davies:2021}, the halo mass function \citep{McClintock:2019, Nishimichi:2019, Bocquet:2022}, and the 21-cm signal \citep{Kern:2017, Cohen:2020, Bevins:2021, Bye:2022}.
Still, the computational resources required to run $\sim30$ simulations with the requisite volume and resolution is highly restrictive, especially for the Lyman-$\alpha$~forest.

Here, we use modern machine learning techniques to alleviate the computational resource cost associated with constructing an emulator.
Specifically, we are concerned with using machine learning to predict high resolution simulation outputs to a high degree of accuracy, while running a minimal number of high resolution simulations.
One machine learning method that is suited to this task is a Gaussian Process (GP) emulator.
Gaussian processes \citep{2006gpml.book.....R} are a means of interpolating between the simulation outputs, providing function prediction in a Bayesian framework.
Essentially, a distribution of functions is learned through training on simulations, and the mean (best estimate) and variance (interpolation error) of the output can be returned for arbitrary simulation inputs.

Previous uses of GP emulators for the Lyman-$\alpha$ forest have been shown to be effective at predicting summary statistics \citep{2019JCAP...02..050B, Rogers:2019, 2021JCAP...05..033P, 2021JCAP...04..059W, Rogers:2021a,2021PhRvL.126g1302R}.
\cite{2019JCAP...02..050B} self-consistently showed that the predicted flux power from their GP emulator (trained with $21$ simulations) agreed to within $1-2\%$ of the corresponding simulation flux power spectrum.
These GP emulators still require a substantial computational cost, as the full simulation suite must be run with sufficient volumes and resolutions for the Lyman-$\alpha$ forest.

Recently, \cite{2022MNRAS.509.2551H} implemented a multi-fidelity GP emulator \citep{10.1093/biomet/87.1.1} for the matter power spectrum.
Here, we combine and expand on the methods outlined in \cite{2019JCAP...02..050B} and \cite{2022MNRAS.509.2551H}, to produce a multi-fidelity GP emulator for the Lyman-$\alpha$ forest flux power spectrum.
In our multi-fidelity model, the training simulations are split into two fidelities; a large sample of low resolution simulations (low fidelity, LF), and a small subset of these simulations run at higher resolution (high fidelity, HF).
Note that we use fidelity and resolution interchangeably throughout this work.
Using these two training sets, the multi-fidelity emulator is trained to predict the 1D flux power spectrum that would be output by a \textit{high} resolution simulation for arbitrary cosmological and astrophysical parameters.

A multi-fidelity emulator allows us to replace some of the HF simulations that would be needed in a single-fidelity emulator with LF simulations.
This can dramatically reduce the computational cost of constructing an emulator, while retaining predictive power across parameter space.
Using this method, emulators can be constructed that make use of the full range of scales and redshifts probed by Lyman-$\alpha$ forest observations.
This enables analyses which can jointly constrain thermal, astrophysical, and cosmological parameters.

In our high resolution simulations, the Lyman-$\alpha$ forest flux power spectrum is converged to $\approx5\%$.
While the scales we probe in this work are resolved, the box size we use is smaller than that required to analyze the full range of scales available in Lyman-$\alpha$ forest data, i.e. we cannot compute a likelihood function using all the real data without larger boxes.
We therefore defer a full cosmological likelihood analysis to future work.
Our goal is to quantitatively test the accuracy of the emulator output, as compared with the output from a set of testing simulations run at the same resolution as the HF training set, and demonstrate the validity and utility of the multi-fidelity technique.
Specifically, we quantify the accuracy of the single- and multi-fidelity emulators with respect to true values from the testing simulations, thus determining how effective the multi-fidelity model is at producing high resolution outputs at minimal computational cost.

\section{Simulations}\label{sec:sims}

\begin{figure*}
    \centering
	\includegraphics[width=2\columnwidth]{figures/spectra_simulation.png}
    \caption{Example Lyman-$\alpha$ forest spectra and corresponding gas density from simulations at redshift $4$.
    The top and bottom panel show the simulated gas surrounding the skewer which produced the spectra shown in the middle panel.
    Examples are shown for simulations run at high (top panel, yellow line) and low resolution (bottom panel, red line).}
    \label{fig:spec_sim}
\end{figure*}

Simulations were performed using MP-Gadget\footnotemark, an N-body and smoothed particle hydrodynamics (SPH) simulation code.\footnotetext{\url{https://github.com/MP-Gadget/MP-Gadget}}
MP-Gadget is adapted from Gadget-3, which is an update of Gadget-2 \citep{2005MNRAS.364.1105S}.
The initial power spectrum and transfer functions are generated with the Boltzmann code CLASS \citep{2011arXiv1104.2932L}.
Species-specific initial conditions are generated for baryons and dark matter \citep{2020JCAP...06..002B, 2021MNRAS.503.1668F}.
We include radiation in the cosmological background model and assume massless neutrinos.

The physics models largely follow those used for the \texttt{ASTRID} simulation, which are described in \cite{2022MNRAS.512.3703B}.
Primary sources for these models, as well as changes from \cite{2022MNRAS.512.3703B} are described below.

We use a cubic kernel for our density estimator (rather than a quintic kernel).
For these simulations, the cubic kernel, in addition to running faster, produced a neutral hydrogen column density distribution that was more consistent with observations for column densities between $10^{20}$ and $10^{22}$ cm$^{-2}$ \citep{2021MNRAS.507..704H}.
We continue to use the pressure-entropy formulation of SPH.
The smaller SPH kernel increases the noise within galaxies, but has minimal effect on the Lyman-$\alpha$ forest \citep{2013MNRAS.429.3341B}.

Star formation follows the model of \cite{2003MNRAS.339..289S}, with our specific implementation as described in \cite{2016MNRAS.455.2778F}.
We lower the number of stars produced per gas particle from $4$ (used in \texttt{ASTRID}) to $1$ (as in Illutris-TNG), which speeds up the simulation without having an effect on the Lyman-$\alpha$ forest.

Black holes follow the model of \cite{2022MNRAS.513..670N}.
We found that for the resolutions used here, the dynamic friction from gas led to a few black holes escaping from their dark matter halo, so we only use dynamic friction from dark matter and stars.
The black hole feedback factor, which controls the fraction of luminosity that is converted to thermal energy, is an emulator parameter (BHF), with the associated parameter limits in Figure~\ref{fig:samples}.
The black hole feedback radius is fixed to $3$ kpc h$^{-1}$, selected to be the average black hole feedback radius at the highest tested resolution when using a nearest neighbour distance.
To accommodate a lower mass resolution than \texttt{ASTRID}, the minimum stellar mass needed in a halo to seed a black hole was increased to $2\times10^{8} M_{\odot}$, and all black hole seeds start with a mass of $5\times10^{4} M_{\odot}$.

Stellar winds are modeled following \cite{2010MNRAS.406..208O}.
The decoupling distance for the winds is increased from $20$ kpc h$^{-1}$ to $1$ Mpc h$^{-1}$, which allows the winds to recouple due to density changes rather than travel distance.
The density threshold for wind recoupling is set to $10\%$ of the star formation density threshold (which is $57.7$ times the critical density).
The minimum wind velocity is set to $100$ km/s.
Finally, metal return (gas enrichment) is disabled as it is not important for the Lyman-$\alpha$ forest and can be computationally expensive.

Gas is assumed to be in ionization equilibrium with a uniform ultraviolet background using the model of \cite{2020MNRAS.493.1614F}.
We boost the temperature of the gas to $15000$ K the timestep after the gas is reionized, to model impulsive heating during hydrogen reionization from ionization fronts \citep{2019ApJ...874..154D}.

We implement He~{\sc ii} reionization using the model of \cite{2020MNRAS.496.4372U}.
The input parameters for this model are: quasar mean bubble size and variance, redshifts for the start and completion of He~{\sc ii} reionization z$^{\text{He~{\sc ii}}}_i$, z$^{\text{He~{\sc ii}}}_f$, and the quasar spectral index $\alpha_q$ (which effectively scales the peak temperature during He~{\sc ii} reionization).
The quasar bubble size is reduced from the default of $\sim30$ Mpc, motivated by radiative transfer simulations \cite{2009ApJ...694..842M}, to $5$ Mpc, due to our small box size.

Simulations are initialised at $z=99$ and finish at $z=2$, and use periodic boundaries.
Box volume, particle number, and gas particle mass resolution are reported in Table~\ref{table:simulations}.
The range given for the gas resolution is due to the varying value of $h$ in our simulation suite.
The gas particle mass resolution for our HF simulations does not meet the required resolution to resolve the forest at all redshifts of interest \citep{2009MNRAS.398L..26B}.
However, the Lyman-$\alpha$ forest flux power spectrum from our HF simulations is converged to within $\approx5\%$ of a simulation that does meet the required resolution of \cite{2009MNRAS.398L..26B}.
We are interested in the performance of the multi-fidelity GP emulator in learning the mapping from low to high resolution, thus this slight lack of numerical convergence does not affect our results.
Examples of the gas density (at $z=3.6$) for the two resolutions are shown in the top and bottom panels of Figure~\ref{fig:spec_sim}.

\begin{table}
	\centering
	\caption{Table of simulation sets}
	\label{table:simulations}
	\begin{tabular}{lcccc}
		\hline
		Simulation & Box Volume & N$_{\text{part}}$ & M$_{\text{gas}}$ (M$_{\odot}$ h$^{-1}$)\\
		\hline
		LF & $(30$ Mpc h$^{-1})^3$ & $2\times256^3$ & $[1.78, 2.37]\times10^7$\\
		HF & $(30$ Mpc h$^{-1})^3$ & $2\times512^3$ & $[2.22, 2.96]\times10^6$\\
		Test & $(30$ Mpc h$^{-1})^3$ & $2\times512^3$ & $[2.22, 2.96]\times10^6$\\
		\hline
	\end{tabular}
\end{table}

Lyman-$\alpha$ forest absorption spectra are generated using Fake Spectra Flux Extractor \citep{2017ascl.soft10012B}\footnotemark, described in \citet{2015MNRAS.447.1834B}.
\footnotetext{\url{https://github.com/sbird/fake_spectra}}
We generate $32,000$ (seeded) randomly placed skewers for each snapshot, from $z=5.4$ to $z=2.0$ in increments of $\Delta z=0.2$.
The pixel resolution is set to $10$ km s$^{-1}$.
An optical depth threshold of $\tau < 10^6$ is set to eliminate damped Lyman-$\alpha$ systems.
Example spectra from one LF and one HF simulation (at $z=4$) are shown in the middle panel of Figure~\ref{fig:spec_sim}.
Note that the LF simulation is not simply a smoothed version of the HF simulation, as the fine velocity structure of the gas moves the location of the absorption peaks.

These sets of neutral hydrogen absorption spectra are used to construct the Lyman-$\alpha$ forest flux power spectrum for each simulation, at each redshift.
The flux power spectrum is defined as $P_F(k) = |L^{-1}\tilde{\delta}^2_F(k)|$, where $\tilde{\delta}^2_F(k)$ is the Fourier transform of the flux excess, $\delta_F(k) = F(k)/\langle F(k) \rangle - 1$, and $L$ is the length of the sightline.
The reported flux power spectrum is averaged over all $32,000$ spectra (see Figure~\ref{fig:low_high_comp} in Section~\ref{sec:accuracy} for example flux power spectra).

Figure~\ref{fig:samples} lists the input parameters that are varied across our suite of simulations, as well as their limits.
Two parameters control the primordial power spectrum: $n_s$ is the scalar spectral index (slope) and $A_p$ is the amplitude (see \cite{2019JCAP...02..050B} for more details).
Three of the parameters relate to the He~{\sc ii} reionization model: z$^{\text{He~{\sc ii}}}_i$ and z$^{\text{He~{\sc ii}}}_f$ are the redshifts for the start and end of He~{\sc ii} reionization, and $\alpha_q$ is the quasar spectral index.
We vary the Hubble constant through $h$, and the total matter density through $\Omega_M h^2$.
One parameter is varied for H~{\sc i} reionization: z$^{\text{H~{\sc i}}}$ is the midpoint redshift of H~{\sc i} reionization.
Finally, we vary one parameter for the black hole model: $BHF$ is the black hole feedback factor, which controls the fraction of luminosity that is converted to thermal energy.
Note that our simulations do not contain a kinetic feedback model.
However, at $z > 2$ it is expected that the thermal mode dominates.
Also shown in Figure~\ref{fig:samples} are the LF and HF training samples.
The selection of these samples is described in Section~\ref{subsubsec:sampling_strategy}.

\begin{figure}
    \centering
	\includegraphics[width=\columnwidth]{figures/30box_40sample_params.pdf}
    \caption{Simulation parameter limits, and samples used to construct the emulators presented.
    Parameters for the low resolution simulations (crosses) were determined by filling a Latin hypercube.
    From the low resolution set, the optimal high resolution set is determined, and these are shown as red circles.
    Initially, $30$ low resolution samples were generated, then an additional $10$ were added while maintaining the Latin hypercube method, hence the non-uniform spacing for the low resolution samples.}
    \label{fig:samples}
\end{figure}

\section{Emulators}\label{sec:emulators}

In Section~\ref{subsec:gp_emu}, we will briefly review emulation using a Gaussian process.
In Section~\ref{subsec:mfemu_emu}, we will review how the Gaussian process emulator can be extended to model simulations with different qualities using a multi-fidelity emulator, \mfemu.
The multi-fidelity emulation technique of \cite{10.1093/biomet/87.1.1} will be reviewed in Section~\ref{subsubsec:mfemu_emu_ko}.
Section~\ref{subsubsection:differences} will discuss the differences in multi-fidelity emulator design between this paper and \cite{2022MNRAS.509.2551H}.
Finally, Section~\ref{subsubsec:sampling_strategy} will outline how we select our training simulations in the parameter space.

\subsection{Gaussian process emulator}\label{subsec:gp_emu}

Gaussian process (\textsc{gp}) regression models \citep{2006gpml.book.....R} have been widely used to build cosmological emulators \citep{Heitmann:2006, Habib:2007, Heitmann:2009}.
A {\textsc{gp}} provides closed-form expressions for predictions.
In addition, a {\textsc{gp}} naturally comes with uncertainty quantification, which is useful when building an inference framework.
In the context of emulation, a {\textsc{gp}} can be seen as a Bayesian prior for the simulation response. It is a prior because the emulator model is chosen to ensure smoothness and monotonicity features of the simulation response \textit{before} data are collected \citep{Santner:2003}.

Let $\thetavec \in \Theta \subseteq \mathbb{R}^{d}$ be the input cosmologies for the simulator, where $d$ is the dimension of the parameters ($d = 9$ for our emulator).
$f(\thetavec)$ is the corresponding output summary statistic.
In this work, the summary statistic, $f(\thetavec)$, is the Lyman-$\alpha$ forest flux power spectrum.
A {\textsc{gp}} regression model can be viewed as a prior on the response surface of our simulated Lyman-$\alpha$ forest flux power spectrum:
\begin{equation}
    f(\thetavec) \sim \mathcal{GP}(\mu(\thetavec), k(\thetavec, \thetavec')),
\end{equation}
where $\mu(\thetavec) = \mathbb{E}[f(\thetavec)]$ is the mean function, and $k(\thetavec, \thetavec') = \mathrm{Cov}[f(\thetavec), f(\thetavec')]$ is the covariance kernel function.
In this work, we assume a zero mean function, and we used the same covariance function as \cite{2019JCAP...02..050B}, which will be defined later in this section.

Suppose we run the simulations at $n$ carefully chosen input cosmologies, $\mathcal{D} = \{\thetavec_1, \cdots, \thetavec_n\}$, and we generate the corresponding Lyman-$\alpha$ forest flux power spectrum for each simulation, $\outputVector = \{ f(\thetavec_1), \cdots, f(\thetavec_n) \}$.
Conditioning on this training data, we can get the predictive distribution of $f$ at a new input cosmology $\thetavec$ through the closed-form expression:
\begin{equation}
    f(\thetavec) \mid \outputVector, \mathcal{D} \sim \mathcal{N}(\mu_n(\thetavec), \sigma_n^2(\thetavec)),
\end{equation}
where the mean and variance are:
\begin{equation}
    \begin{split}
        \mu_n(\thetavec) &= \kvec(\thetavec, \mathcal{D})^\intercal \Kvec(\mathcal{D})^{-1} \outputVector;\\
        \sigma_n^2(\thetavec) &= k(\thetavec, \thetavec) - \kvec(\thetavec, \mathcal{D})^\intercal \Kvec(\mathcal{D})^{-1} \kvec(\thetavec, \mathcal{D}).
    \end{split}
\end{equation}

The vector $\kvec(\thetavec, \mathcal{D}) = [ k(\thetavec, \thetavec_1), \cdots, k(\thetavec, \thetavec_n) ]$ represents the covariance between the new input cosmology, $\thetavec$, and the training data.
The matrix $\Kvec(\mathcal{D})$ is the covariance for the training data.

For the covariance kernel function, we choose the same kernel as in \cite{2019JCAP...02..050B}, which is a combination of a linear kernel and a radial basis kernel (RBF):
\begin{equation}
    \begin{split}
        k(\thetavec,\thetavec'; \sigma_0, \boldsymbol{l}, \boldsymbol{\sigma}) &= k_\mathrm{RBF}(\thetavec, \thetavec'; \sigma_0, \boldsymbol{l}) + k_\mathrm{LIN}(\thetavec, \thetavec'; \boldsymbol{\sigma})\\
        &= \sigma_0^2 \exp{\left( \sum_{i=1}^{d} -\frac{(\thetavec_i - {\thetavec_i}')^2}{2 l_i^2} \right)} +  \sum_{i=1}^{d} \sigma_i^2 \thetavec_i {\thetavec_i}',
    \end{split}
\end{equation}
where $\sigma_0^2$ and $\boldsymbol{\sigma}^2$ are the variance hyperparameters for the RBF kernel and the linear kernel, respectively.
$\boldsymbol{l}$ is the lengthscale parameter that controls the smoothness of the Gaussian process function.
We applied Automatic Relevance Determination (ARD) for both linear and RBF kernels.
That is, we assign one lengthscale $l_i$ (variance $\sigma_i$) hyperparameter for each input dimension $i$ for the RBF and linear kernels.
This allows the GP to dynamically learn the scale over which each input dimension varies, which corresponds to the degree of sensitivity of the flux power spectrum to the input parameter.

Although we do not explicitly write in the notation, $f(\thetavec)$ is a single-valued output.
Since our target summary statistic is a vector, we model each $k$-bin of the flux power spectrum with a separate {\textsc{gp}}.
The primary reason for this choice is that the correlation between the low-fidelity and high-fidelity flux power spectrum changes depending on the scale considered.
The multi-fidelity method can only capture this scale dependence if we model each scale separately.

\subsection{Multi-Fidelity Emulation}\label{subsec:mfemu_emu}

We first introduce the Kennedy-O'Hagan model (KO model) \citep{10.1093/biomet/87.1.1} in Section~\ref{subsubsec:mfemu_emu_ko}.
Section~\ref{subsubsection:differences} describes the changes we have made to adapt the model from \cite{2022MNRAS.509.2551H} to the Lyman-$\alpha$~forest.
Finally, the strategy we employ for choosing parameters at which to generate high-fidelity training simulations is described in Section~\ref{subsubsec:sampling_strategy}.

\subsubsection{Kennedy O'Hagan Method}
\label{subsubsec:mfemu_emu_ko}

The KO model \citep{10.1093/biomet/87.1.1} was first introduced to model a sequence of computer codes with increasing fidelity.
For simplicity, we assume there are only two fidelities:
low-fidelity (LF) simulations with low resolution,
and high-fidelity (HF) simulations with high resolution.

We define $\{\outputVector_\mathrm{LF}, \outputVector_\mathrm{HF}\}$ as the {Lyman-$\alpha$} forest flux power spectra in the training set. $\outputVector_\mathrm{LF} = \{f_\mathrm{LF}(\thetavec^\mathrm{LF}_i) \}_{i=1}^{n_\mathrm{LF}}$  and $\outputVector_\mathrm{HF} = \{f_\mathrm{HF}(\thetavec^\mathrm{HF}_i) \}_{i=1}^{n_\mathrm{HF}}$, where $n_\mathrm{LF}$ and $n_\mathrm{HF}$ are the number of simulations in the low- and high-fidelity training sets.
We use the KO method to model Lyman-$\alpha$~forest flux power spectra from different fidelities:
\begin{equation}
    f_\mathrm{HF}(\thetavec) = \rho \cdot f_\mathrm{LF}(\thetavec)
    + \delta(\thetavec),
    \label{eq:ko_model}
\end{equation}
where $\rho$ is a trainable parameter describing a multiplicative correction between the low- and high-fidelity {Lyman-$\alpha$} forest flux power spectra.
$\delta(\thetavec)$ is a GP independent of $f_\mathrm{LF}(\thetavec)$, describing an additive correction between fidelities.
In other words, Equation~\ref{eq:ko_model} assumes the high-fidelity {Lyman-$\alpha$} forest flux power can be decomposed as the low-fidelity flux power multiplied by a correction parameter, $\rho$, and an additive bias function $\delta(\thetavec)$.

As mentioned in \cite{2022MNRAS.509.2551H}, the $\rho$ parameter has to be scale-dependent (a function of $k$) to model the well-known fact that small scales are less well resolved in smaller simulations.
Here we use the same method as \cite{2022MNRAS.509.2551H} and assume Equation~\ref{eq:ko_model} is a single-output GP model. We assign a KO model to each $k$ bin of the data.\footnote{We can easily get the same set of $k$ bins for low- and high-fidelity by using the same spectral resolution for both simulations.}
In this way, we can model $\rho$ as a function of $k$, as shown in Figure~\ref{fig:ko_rho}.

\begin{figure}
    \includegraphics[width=\columnwidth]{figures/rho_ko_mf_linear.pdf}
    \caption{The scale parameter, $\rho$, of the KO model (Equation~\ref{eq:ko_model}) as a function of $k$.
    Different colors represent different redshifts.
    We separate the redshifts into four bins, $2 \leq z < 3$, $3 \leq z <4$, $4 \leq z < 5$, and $5 \leq z \leq 5.4$.
    Within each redshift bin, we average over fixed $k$ modes, which are linearly spaced between the maximum $k$ and minimum $k$ at the given redshift range.
    This emulator used 40 LF and 6 HF for training.
    }
    \label{fig:ko_rho}
\end{figure}

We also assign KO models for each redshift.
As shown in Figure~\ref{fig:ko_rho}, $\rho$ is a non-trivial function of both $k$ and $z$, so we cannot simply use an emulator trained on one redshift to apply on another redshift.\footnote{See \cite{2021JCAP...05..033P} for a Lyman-$\alpha$ forest emulator that uses a single GP for all redshifts, and achieves sub-percent accuracy, albeit with some ambiguity between model parameters and redshift.}
We note that it is possible to assume a smooth function to model $\rho(k, z)$.
However, validating $\rho(k, z)$ is out-of-scope for this paper.
In practice, observational data are conditioned on a specific redshift, so training emulators on separate redshifts is sufficient for cosmology inference.

Figure~\ref{fig:ko_rho} shows that $\rho$ stays close to unity at large scales for most of the redshifts.
At small scales, however, different redshifts require different values of $\rho$.
At the middle redshifts ($3 \leq z < 5$), $\rho$ has a large positive deviation from unity.
At the low redshifts ($2 \leq z < 3$), $\rho$ has a moderate deviation toward $\rho > 1$.
The only exception is at the high redshifts ($5 \leq z \leq 5.4$), which stays close to $\rho = 1$ for all scales.
This indicates the resolution correction in {Lyman-$\alpha$} forest flux power varies with redshift, depending on the over-density probed by the forest.

\subsubsection{Model differences from Ho et al. (2022)}\label{subsubsection:differences}

Here we highlight the ways in which we have adapted the model from \cite{2022MNRAS.509.2551H}, which emulated the non-linear matter power spectrum at $z= 0 - 2$, to the Lyman-$\alpha$~forest at $z=2 - 5.4$.
Since the redshift range is larger in this work, we employed a new strategy to select the optimal HF training set by averaging over the interpolation loss for all redshift bins.
We will describe the strategy in detail in Section~\ref{subsubsec:sampling_strategy}.

In \cite{2022MNRAS.509.2551H}, we outlined two multi-fidelity methods: a linear multi-fidelity emulator (the KO model, or the autoregression model (AR1)), and a non-linear multi-fidelity emulator (non-linear autoregressive GP, or NARGP, \cite{2017RSPSA.47360751P}).
However, we found that NARGP requires more HF training simulations for the {Lyman-$\alpha$} forest flux power than AR1, perhaps due to the wide range of redshifts used.
We use the KO model for our main results, and describe the NARGP results in Appendix~\ref{sec:nonlinear}.

In this work, instead of emulating logarithm scaled powers, we adopted the mean-flux normalization strategy proposed in \cite{2019JCAP...02..050B}.
We normalize all flux power spectra in the training set by the median spectrum:
\begin{equation}
    \begin{split}
        \boldsymbol{y}_\mathrm{LF} &\leftarrow \frac{\boldsymbol{y}_\mathrm{LF}}{\mathrm{median}_i{(\boldsymbol{y}_\mathrm{LF})}} - 1;\\
        \boldsymbol{y}_\mathrm{HF} &\leftarrow \frac{\boldsymbol{y}_\mathrm{HF}}{\mathrm{median}_i{(\boldsymbol{y}_\mathrm{LF})}} - 1.
    \end{split}
    \label{eq:normalization}
\end{equation}
The index $i$ refers to one of the spectra in the training set, $\outputVector_\mathrm{LF} = \{f_\mathrm{LF}(\thetavec^\mathrm{LF}_i) \}_{i=1}^{n_\mathrm{LF}}$.
Equation~\ref{eq:normalization} ensures the training sample distribution is close to having a zero mean, matching the prior of the GP emulator.
We found that in practice this normalisation makes training the emulator substantially easier. Note that we normalize the HF training set using the same LF median spectrum. As the HF training set is small, the median spectrum estimate for HF is noisy, and so using it for normalization may introduce some unwanted training bias.

\subsubsection{Sampling strategy for high-fidelity simulations}\label{subsubsec:sampling_strategy}

The KO model approach can be seen as a Bayesian way to correct an emulator from low-fidelity to high-fidelity.
Thus, if $\outputVector$ is the high-fidelity {Lyman-$\alpha$} forest flux power spectrum used for training, and $\thetavec$ is the corresponding input parameters:
\begin{equation}
    \begin{split}
        \outputVector &= f_\mathrm{LF}(\thetavec) + (f_\mathrm{HF}(\thetavec) - f_\mathrm{LF}(\thetavec))\\
              &= f_\mathrm{LF}(\thetavec) + \mathrm{error}(\thetavec).
    \end{split}
\end{equation}
The emulation accuracy will be directly affected by how well an autoregressive construction can model $\mathrm{error}(\thetavec)$.
Usually, a large set of low-fidelity simulations are used as training data for $f_\mathrm{LF}(\thetavec)$ because they can be obtained cheaply.
The quality of training data for $\mathrm{error}(\thetavec) = f_\mathrm{HF}(\thetavec) - f_\mathrm{LF}(\thetavec)$ thus relies on the choice of high-fidelity simulations.

In \cite{2022MNRAS.509.2551H}, we proposed an optimization strategy to select high-fidelity training simulations.
A low-fidelity only emulator (LFEmu)\footnote{In a similar way, we call a high-fidelity only emulator, HFEmu.} is trained on a subset of low-fidelity training simulations.
The posterior means of the trained LFEmu are used to calculate the emulation errors from the remaining LF samples in the Latin Hypercube Sampling (LHS).
By minimizing the emulation errors of LFEmu, we can grid search for the optimal set of cosmologies that best interpolates the parameter space using a small number of training simulations.
Assuming LFEmu is correlated with HFEmu, we can use the selected optimal cosmologies as inputs for the HF training set.
By ensuring the HF training set achieves a good interpolation, we mitigate emulation errors for the multi-fidelity emulator.

In practice, we employ a three-stage procedure for building a multi-fidelity emulator:
\begin{enumerate}
    \item Prepare LF simulation suite.
    \item Prepare HF simulation suite. This is done by using LFEmu to find the set of cosmologies that minimizes the interpolation loss.
    \item Build {\mfemu}. If the accuracy is not enough, go back to stage 1 or 2 to run more training simulations.
\end{enumerate}
For stage (ii), to avoid wasting computational resources running more LF simulations, we directly use the LF simulation suite in stage 1 to build and validate the LFEmu.
Thus, the chosen cosmologies are a subset of the LF simulation LHS, which fulfills the nested training dataset design suggested in \cite{10.1093/biomet/87.1.1}.
The benefit of using a nested data structure, $\thetavec_\mathrm{HF} \subseteq \thetavec_\mathrm{LF}$, is that we can directly compute posterior means from the LF training set for cosmologies $\thetavec_\mathrm{HF}$, without any interpolation in LF.

We note that it is possible to train a {\mfemu} without using the LF simulations to optimize the HF points.
However, if the selection of HF points are suboptimal (i.e., can barely interpolate in the prior volume), then the {\mfemu} accuracy will be suboptimal.
This is because the $\mathrm{error}(\theta)$ cannot be decomposed into an autoregressive structure easily.

To find the optimal HF training set across the full redshift range, $z = 2 - 5.4$, we train a LFEmu for each redshift and get the validation loss (we used mean squared errors).
We sum up the validation loss for all redshifts and find a subset of cosmologies that minimizes the summed validation loss.

In Figure~\ref{fig:hf_optimize_flowchart}, we summarize the above described procedure in a flowchart.
For a formal description, we refer to \cite{2022MNRAS.509.2551H}.
The only difference between the proposed procedure in Figure~\ref{fig:hf_optimize_flowchart} and the procedure in \cite{2022MNRAS.509.2551H} is that we optimize $\thetavec_\mathrm{HF}$ across redshifts $z = 2 - 5.4$ in this work.
Thus, we have an additional step to sum the LFEmu validation loss across redshifts.

\begin{figure}
    \includegraphics[width=\columnwidth]{figures/mf_lya_hf_optimize_flowchart.pdf}
    \caption{A flowchart for training a multi-fidelity emulator.
    We start with a low-fidelity (LF) set of simulations.
    We then select a subset of low-fidelity points to train an emulator, LFEmu, at each redshift.
    The validation loss for each LFEmu is computed using the rest of the LF simulations as a validation set.
    Finally, we sum up the validation losses at each redshift, and use the summed loss to propose a set of cosmologies $\thetavec_\mathrm{HF}$ which can best minimize the summed loss.}
    \label{fig:hf_optimize_flowchart}
\end{figure}

Training an LFemu on all possible low-fidelity subsets is computationally intensive.
To reduce costs, we employed the greedy optimization strategy from \cite{2022MNRAS.509.2551H}.
We first explored all possible subsets for $3$ design points within the LF LHS.
For the optimal $4$ design points, instead of exploring all possible subsets, we grew the subset one point at a time, fixing the previously chosen optimal $3$ HF points.
In the same line of thought, we grew the subset to $6$ optimal design points for HF training cosmologies.
Our final simulation suite of $40$ LF and $6$ HF samples, along with parameter limits, is shown in Figure~\ref{fig:samples}.

\section{Results \& Discussion}\label{sec:results}

Using the flux power spectra from our LF and HF simulations, we train single-fidelity (one LF only, and one HF only) and multi-fidelity emulators.
These trained emulators are used to predict the flux power spectrum output for a set of $10$ simulation input parameters.
We then compare these predictions to the corresponding testing simulations which were run at the same resolution as the HF simulations (see Table~\ref{table:simulations}).

\subsection{Emulator Accuracy}\label{sec:accuracy}

\begin{figure}
    \centering
	\includegraphics[width=\columnwidth]{figures/fluxpower_comp.pdf}
    \caption{Lyman-$\alpha$ forest flux power spectrum from LF and HF simulations.
    The top panel shows the flux power spectrum at redshift $3$ from an LF simulation (blue), and its HF counterpart (yellow).
    The lower panels show the ratio of the flux power for these two resolutions, for all LF-HF simulation pairs, at $z=3.8, 4.6, \& \ 5.4$.
    Dashed lines shown the highest $k$ probed by BOSS/eBOSS \citep{2019JCAP...07..017C}, and the estimated reach for DESI \citep{2022arXiv220307491V}.}
    \label{fig:low_high_comp}
\end{figure}

In the following, we only show results for the emulators that use the full available set of training simulations ($40$ LF and $6$ HF).
We have verified that using all available training simulations leads to the most accurate emulator.
Section~\ref{sec:runtime} shows how emulator accuracy degrades when a smaller subset of the available simulations is used.

Using the full set of available simulations, the mean prediction error for the multi-fidelity emulator is $\langle| P_F^{\text{pred}}/P_F^{\text{true}} - 1 |\rangle \approx0.8\%$ (averaging across all scales, redshifts, and testing simulation outputs).
For the LF single-fidelity emulator, the mean prediction error is $\approx4\%$.
This is not unexpected; there are real differences between the flux power spectra output by the low and high resolution simulations that are not being captured with this method.
Figure~\ref{fig:low_high_comp} shows flux power spectra from a LF simulation and its HF counterpart, and the ratio of these at several redshifts.
While the exact difference between the LF and HF flux power spectra depends on simulation input parameters and redshift, in general the LF differs most from the HF at small scales and high redshifts.
The $4\%$ error may seem quantitatively quite good, considering the simpler methodology and reduced resource cost.
However, there is no indication that the error could be reduced further with additional simulations (see Section~\ref{sec:runtime}).

For the HF single-fidelity emulator, the mean prediction error is $\approx3\%$.
This is likely limited by the sample size of the training set ($6$ simulations), leading to increased errors when making predictions for inputs that are far away from the training samples.
It is important to note that the HF samples are selected to optimize the multi-fidelity emulator, rather than as an independent emulator (i.e. as a Latin hypercube sample).
There is some indication that prior information from the LF training samples provides useful information about the best areas of parameter space to sample.
To test this, we split our testing set ($10$ simulations, same resolution as the HF samples) into training and testing sets, then train all $210$ combinations of $6$ samples, and predict the outputs for the remaining $4$ samples.
The error range from this exercise is $2.5-11.5\%$ ($5.5\%$ mean error, $1.5\%$ standard deviation).
Though not a direct comparison, the $3\%$ error we obtain from the HF single-fidelity emulator compares favorably with this, indicating that the HF samples selected are an improvement over using a Latin hypercube sampling scheme.

\begin{figure}
    \centering
	\includegraphics[width=\columnwidth]{figures/kf_error.pdf}
    \caption{Comparing the prediction error as a function of (linearly binned) wavevector for multi- and single-fidelity emulators.
    This is the mean error across all redshifts and $10$ test simulations.
    The shaded regions are the variance in the prediction error.
    Dashed lines show the highest $k$ probed by BOSS/eBOSS \citep{2019JCAP...07..017C}, and the estimated reach for DESI \citep{2022arXiv220307491V}.}
    \label{fig:k_error}
\end{figure}

Figure~\ref{fig:k_error} shows the mean prediction error, averaged over all redshifts and $10$ test simulation outputs, as a function of wavevector $k$.
In this, and the following figures, the shaded region around the curves is the variance in prediction error ($| P_F^{\text{pred}}/P_F^{\text{true}} - 1 |$), to give a sense of how much the error varies beyond the mean.
The multi-fidelity emulator outperforms the single-fidelity emulators at all scales, with an error between $0.5-1.5\%$.
The LF (HF) single-fidelity emulator has error between $2-7\%$ ($2-6\%$).

Both single-fidelity emulators and the multi-fidelity emulator trend towards higher error for small scales.
The LF emulator dips $1-2\%$ around $k\approx0.02$ s/km.
The dip occurs on scales at which the low resolution flux power spectra goes from overestimating to underestimating the high resolution power (see Figure~\ref{fig:low_high_comp}, for $z\le4.6$).
The uptick in the multi-fidelity emulator error for the largest k-bins is also present in the HF emulator, indicating that there is a scarcity of large scale modes available in the emulator training.

\begin{figure*}
    \centering
	\includegraphics[width=2\columnwidth]{figures/z_scales_error.pdf}
    \caption{Comparing the prediction error as a function of redshift and scale for multi- and single-fidelity emulators.
    This is mean error across all $10$ test simulations.
    The blue shaded region shows the extent of helium reionization in our simulations (see parameter limits in Figure~\ref{fig:samples}).}
    \label{fig:z_scales_error}
\end{figure*}

Figure~\ref{fig:z_scales_error} shows how the emulators perform as a function of redshift and scale.
In the following, we define small scales as $k > 0.02$ s/km and large scales as $k \le 0.02$ s/km (divided at the smallest scale accessible by BOSS/eBOSS data \citep{2019JCAP...07..017C}).
On large scales (Figure~\ref{fig:z_scales_error}, left panel), the LF emulator error decreases with time until $z\approx4$, then slowly increases.
This trend is because, as can be seen in Figure~\ref{fig:low_high_comp} (left of the BOSS/eBOSS dashed line), the low resolution flux power comes into better agreement with the high resolution flux power as it nears $z\approx4$.
We also found that the LF simulations do not cool as efficiently as the HF simulations after He~{\sc ii} reionization, likely leading to the rise in error for $z\le3.2$.

On small scales (Figure~\ref{fig:z_scales_error}, right panel), the LF emulator error is more variable.
The dip in error seen around $z\approx4.6$ occurs as the low resolution flux power crosses from overestimating to underestimating the high resolution flux power (Figure~\ref{fig:low_high_comp} right of the BOSS/eBOSS dashed line).
The subsequent rise in error is due to the loss of small-scale power and consequent under-estimation of the flux power spectrum in the low-fidelity simulations.

The trends in redshift and scale seen in the LF emulator performance are due not to interpolation error, but to the different numerical resolution of the two simulation fidelities, since this emulator is not predicting the flux power at the higher resolution.
Some differences are connected to temperature differences between the LF and HF simulations.
For low densities ($\sim1$ times the mean density), the LF simulations are colder than the HF simulations at high redshift, but come into better agreement leading up to He~{\sc ii} reionization, after which they diverge from the HF simulations again.
For higher densities ($1-100$ times the mean density), the LF simulations are once again colder than the HF simulations at high redshift, but at lower redshift they are too hot (with a crossover at $z\approx4.6$).
As higher redshifts probe lower densities, the error initially decreases with redshift, before rising again towards the lowest redshifts.

On both large and small scales the HF emulator errors are around $3.5\%$, dominated by sampling variance.
During He~{\sc ii} reionization the HF emulator has more variation in error, which is probably exacerbated by our small box sizes.
The increased variation during He~{\sc ii} reionization further indicates that the primary source of error for the HF emulator is the sample size of the training set.

On small scales, the multi-fidelity emulator error is insensitive to redshift and small ($0.9\%$).
On large scales, the multi-fidelity emulator error slightly decreases until $z=4.2$, then increases with the onset of He~{\sc ii} reionization, before flattening again.
The trend is also more variable during He~{\sc ii} reionization, indicating that emulator finds it more difficult to learn the mapping during this process.
However, the multi-fidelity emulator still outperforms the single-fidelity emulators, with an error between $\approx0.4-1\%$.

\subsection{Emulator Runtime}\label{sec:runtime}

While we have shown that the multi-fidelity emulator outperforms the single-fidelity emulators presented here, it still remains to show that it is more computationally cost efficient.
We could, for example, add more training simulations to our single-fidelity HF emulator and get a similarly accurate high resolution emulator.
However, the computational cost would increase significantly.
By comparing the total emulator runtime to prediction error, we can determine the choice that balances computational cost and accuracy.
In practice, the important question is to determine the computational cost at which a given emulation technique can achieve a desired accuracy.
The computational cost of training the emulators is subdominant ($\mathcal{O}(1)$ cpu-hours) to running the training simulations, so in the following we only consider the runtime for the simulation suites.

\begin{figure}
    \centering
	\includegraphics[width=\columnwidth]{figures/nsim_scales_error.pdf}
    \caption{Emulator prediction error as a function of the number of simulations used in training the emulator.
    This is the mean error across all redshifts and $10$ test simulation outputs.
    The prediction error is broken down into large ($k < 0.02$ s/km) and small ($k > 0.02$ s/km) scales, as in Figure~\ref{fig:z_scales_error}.
    The solid lines show how the average error depends on the number of LF simulations, while the dashed lines show how the average error depends on the number of HF simulations once all LF simulations are included.
    }
    \label{fig:nsim_scales}
\end{figure}

Figure~\ref{fig:nsim_scales} shows the mean prediction error (averaged over all redshifts and test outputs) as a function of the number of simulations used in the training set, for small and large scales (as defined in Section~\ref{sec:results} and Figure~\ref{fig:z_scales_error}).
The solid lines show prediction errors for multi-fidelity emulators trained using $6$ HF and a varying number of LF simulations.
The small and large scale errors flatten out after $\approx 30$ LF simulations are used in the training set.
The LF simulations allow the emulator to determine how the flux power spectrum depends on the cosmological input parameters, and so this indicates that $30$ LF simulations are needed to explore our $9$ parameter space.
Other emulators range from using $\approx6$ simulations per parameter (e.g. \cite{McClintock:2019, 2022MNRAS.509.2551H}) to $30$ per parameter (e.g. \cite{Euclid:2021}).
The $3-4$ simulations per parameter required here is unusually low, perhaps because the input parameters affect the flux power spectrum close to linearly in much of parameter space.

Dashed lines show prediction errors for multi-fidelity emulators trained using $40$ LF and a varying number of HF simulations.
Adding extra HF simulations to the training set has a larger impact than adding LF simulations.
The addition of each HF simulation generally improves the emulator accuracy for small scales more than for large scales.
This is as expected, since the main purpose of the HF simulations is to learn the mapping from low to high resolution output, with small scales being more resolved in the HF simulations.

Figure~\ref{fig:runtime} shows the emulator prediction errors as a function of the total runtime (cost of running the training simulations).
All simulations were run on the Frontera supercomputer at the Texas Advanced Computing Center.
The cost is divided between the LF training simulations, which cost $\approx10$ node hours each, and the HF training simulations, which cost $\approx150$ node hours each.

\begin{figure}
    \centering
	\includegraphics[width=\columnwidth]{figures/runtime_error.pdf}
    \caption{Prediction error as a function of total training simulation computational cost.
    This is the mean error for $10$ test simulations over all redshifts and scales.
    The solid line shows the prediction error when changing the number of HF simulations used in training the multi-fidelity emulator.
    Also shown are the single-fidelity emulators (red squares), and the multi-fidelity emulator trend with a varying number of LF training simulations (yellow, dashed).}
    \label{fig:runtime}
\end{figure}

The dashed trend shows the same emulators as Figure~\ref{fig:nsim_scales}, but no longer divided by scale.
Qualitatively it looks the same as both the large scale and small scale results from the previous figure.
The error is flat after $\approx30$ LF training simulations, indicating that a similar accuracy can be achieved using the multi-fidelity emulator with $\approx30$ rather than $40$ LF simulations.
The most efficient $1\%$ error emulator in this study is a multi-fidelity emulator using $30$ LF and $5$ HF training simulations (the cost for this was $\approx1050$ node hours).
The most accurate emulator is the $40$ LF, $6$ HF multi-fidelity emulator, with error $0.8\%$, and cost $\approx1300$ node hours.

The dotted line (squares) shows the error and runtime for the single-fidelity emulators.
Following from the $6$ HF single-fidelity emulator result to the dashed (yellow) line, it can be seen that the addition of just a few LF training simulations quickly improves the accuracy.
It can also be seen that in terms of computational cost, the multi-fidelity emulator is more efficient\footnote{At least for errors less than $4\%$, the approximate amount by which the LF simulations fail to be converged.}.
Note that the HF training simulations are not selected to optimize a single-fidelity emulator, but instead are selected to optimize the multi-fidelity emulator.
They thus use prior information provided by the LF training simulations and so perform better than a naive Latin hypercube construction of a HF emulator using $6$ training samples.
Our multi-fidelity scheme is thus an even larger improvement on a single-fidelity model than Figure~\ref{fig:runtime} suggests.

The solid line shows the error and runtime for the multi-fidelity emulator trained using $40$ LF simulations, and $2-6$ HF simulations.
The point on the solid line corresponding to $40$ LF, $2$ HF has a similar cost, but slightly worse performance than the $5$ HF single-fidelity result.
Adding a third HF training sample decreases the error more for the multi-fidelity emulator (error for $40$ LF, $3$ HF emulator) than it does for the single-fidelity emulator (error marked $6$ HF).
Adding a $6$th HF simulation to the $40$ LF, $5$ HF multi-fidelity emulator produces a relatively small improvement in error, perhaps indicating that stochasticity in the simulations due to our relatively small box size, is beginning to dominate over interpolation error.
\section{Conclusions}\label{sec:conclusions}

In this work we developed and tested a multi-fidelity emulator for the simulated Lyman-$\alpha$ forest flux power spectrum.
Emulators address the growing computational demands of simulations, which must be run at increasingly high resolutions to allow analysis of the increasing quality and quantity of observational data.
Here, we use a Gaussian process based emulator that addresses this demand by, in a Bayesian framework, training an interpolating function to predict the output (Lyman-$\alpha$ forest flux power spectrum) for a given input (simulation input parameters).
Relatively few simulations are required to accurately predict across the span of input parameter space, making emulators especially useful for parameter inference problems.

The multi-fidelity framework allows a further reduction in computational cost by dividing the emulator training samples into multiple (in our case two) fidelities.
The low-fidelity (low-resolution) training samples allow the emulator to learn how the the outputs depend on input parameters.
The high-fidelity (high-resolution) training samples correct numerical errors in the low-fidelity emulator with a (parameter-dependent) mapping from low- to high-fidelity.
Thus, the emulator can be trained with a large sample of low-fidelity training simulations and a small subset of high-fidelity training simulations.

Our training suite included $40$ low resolution hydrodynamical simulations ($30$ Mpc/h simulation box length, $256^3$ particles) and $6$ high resolution hydrodynamical simulations ($30$ Mpc/h simulation box length, $512^3$ particles).
Using the Lyman-$\alpha$ forest flux power spectrum extracted from these simulations, we trained single- and multi-fidelity emulators to predict the high resolution flux power spectrum.
Ten independent simulations were run to test the prediction accuracy of the trained emulators.

In summary, the multi-fidelity emulator:
\begin{itemize}
    \item Modelled a redshift range $5.4 - 2$ (in $18$ redshift bins with $\Delta z=0.2$), on scales ranging from $k=1.4\times10^{-3}$ to $k=5.7\times10^{-2}$ s/km (in $25$ bins).
    \item Achieved sub-$1\%$ error on most scales and redshifts when averaged over $10$ test simulations.
    \item Reached an average error of $0.8\%$.
    \item Achieved $1\%$ average error most cost efficiently using a training set with $30$ low resolution and $5$ high resolution simulations.
\end{itemize}

The low resolution single-fidelity emulator ($4\%$ average error) predicts the low resolution flux power, so it is limited by real differences between the output of the two resolutions.
The high resolution single-fidelity emulator ($3\%$ average error) is limited by the small number of training samples.
It is likely that the average error for the high resolution single-fidelity emulator could be improved to match the multi-fidelity emulator performance with the addition of more training simulations.
However, the high resolution single-fidelity emulator quickly increases in computational cost with additional samples, and we expect it would thus be more expensive than our multi-fidelity emulator.

Some important caveats to our results are that the Lyman-$\alpha$ forest is converged at the $\approx5\%$ level in our high resolution simulations, and the box size is small.
In future work, the effect box size has on the multi-fidelity model will be explored, and a model where the fidelities have different box sizes will be tested.
While there is no direct evidence to suggest that changing the resolution or box size would significantly enhance or diminish the accuracy of the emulators presented here, it still remains to be tested on simulations with higher resolution and larger box sizes.
In a forthcoming work, we test the multi-fidelity framework on larger box size, higher resolution simulations, as well as using such an emulator for cosmological inference.

\section{Non-linear Multi-Fidelity Emulator}\label{sec:nonlinear}

\begin{figure*}
    \centering
	\includegraphics[width=2\columnwidth]{figures/non-linear_error.pdf}
    \caption{Comparing the prediction error as a function of redshift and scale for linear and non-linear multi-fidelity emulators.
    This is mean error across all $10$ test simulations.
    The shaded region shows the extent of helium reionization in our simulations (see parameter limits in Figure~\ref{fig:samples}).}
    \label{fig:nonlinear}
\end{figure*}

In the main text we have explored the effectiveness of a linear multi-fidelity emulator (the KO model, or AR1).
In the linear model, the mapping from LF to HF is $f_{\text{HF}}(\thetavec) = \rho f_{\text{LF}}(\thetavec) + \delta(\thetavec),$ where $f_{\text{HF}}$ and $f_{\text{LF}}$ are the emulator predictions at those resolutions, and $\rho$ is independent of the input parameters $\thetavec$.

Here, we compare those results with the results using a non-linear multi-fidelity emulator (non-linear autoregressive GP, or NARGP).
In the non-linear multi-fidelity model, proposed by \cite{2017RSPSA.47360751P}, the mapping is a function of both the LF output and the input parameters.
We model this as:
$$f_{\text{HF}}(\thetavec) = \rho(\thetavec, \tilde{f}_{\text{LF}}(\thetavec)) + \delta(\thetavec),$$
such that $\rho$ depends on both the input parameters and LF posterior output.
The LF outputs, as is the case with the linear model, are median normalized such that the assumption on the Gaussian process of zero mean is more reasonable, $\tilde{f}_{\text{LF}}(\thetavec) = f_{\text{LF}}(\thetavec)/\mu_{\text{LF}} - 1$.

Following \cite{2017RSPSA.47360751P}, $\rho$ is modelled as a Gaussian process with input from both input cosmologies for HF, $\thetavec_\mathrm{HF}$, and the output from LF, $\tilde{f}_{\text{LF}}(\thetavec)$.
The NARGP construction results in a deep Gaussian process model \citep{pmlr-v31-damianou13a}.
We follow the approximation in \cite{2017RSPSA.47360751P} and replace $\tilde{f}_{\text{LF}}(\thetavec)$ with its posterior distribution.
Thus, the training reduces to training two regular GPs recursively.

In Figure~\ref{fig:nonlinear} we show the prediction errors separated into small and large scales, as a function of redshift for both the linear and non-linear multi-fidelity emulators.
They perform similarly, with the linear emulator being more accurate at low redshifts on all scales, and high redshifts for small scales.
The difference between the average error for the linear and non-linear models (over all scales and redshifts) is $0.08\%$.
This is in contrast to emulation of the matter power spectrum in \cite{2022MNRAS.509.2551H}, where the non-linear model outperformed the linear model.

It is worth noting that the non-linear model agrees closely with the linear model when using the full suite of training simulations, but lags behind the linear model when using fewer HF training simulations.
For example, the difference in the average error between linear and non-linear models using $40$ LF and $3$ HF is $\approx2\%$ ($1.9\%$ error for linear, $3.7\%$ error for non-linear).
When using $4$ HF the difference is $\approx1\%$ ($1.5\%$, $2.5\%$), and when using $5$ HF the difference is $\approx0.4\%$ ($0.9\%$, $1.3\%$).
While differences in the effectiveness of the non-linear model may be due to the quantity being emulated (matter power versus flux power), one likely reason for the difference is the number of input parameters.
In \cite{2022MNRAS.509.2551H}, five input parameters are used, while in this work we use nine.
The non-linear model uses the posterior of the LF output, which requires Monte-Carlo sampling.
It is possible that the additional dimensions degrade the performance of the Monte-Carlo integration, and thus the performance of the non-linear model.
One other reason may be the larger number of hyperparameters that need to be optimized in the training.


\section{Introduction}

  Two of the most fundamental goals of planet formation studies are to understand how and when planets gain mass. The planetary mass assembly process can be directly constrained by characterizing actively accreting protoplanets. So far, the discoveries of protoplanets are limited to PDS~70 b and c, two gas giants directly imaged within the same planetary system \citep{Keppler2018,Haffert2019}. A plethora of follow-up studies on this system confirmed the presence of circumplanetary disks \citep{Isella2019,Benisty2021}, estimated the mass accretion rates \citep[e.g.,][]{Wagner2018, Haffert2019, Hashimoto2020,Zhou2021}, investigated the planet-disk interactions \citep[e.g.,][]{Bae2019}, and constrained the planetary orbits \citep{Wang2021}. Despite this steady progress, these studies are confined to one protoplanetary system and only represent a single outcome of the planet formation process observed at an instant in time.

  To expand the sample of accreting planets, we launched the Hubble Accreting Luminous Protoplanets in H-Alpha (HALPHA) Survey (Program ID: 16651\footnote{Program information can be found at \url{https://www.stsci.edu/cgi-bin/get-proposal-info?id=16651&observatory=HST}}, PI: Zhou) to search for accreting planets in ten transition disk systems. We are exploiting the recently demonstrated optical and ultraviolet (UV) high-contrast imaging capibility of Wide Field Camera 3 (WFC3) UVIS \citep{Zhou2021,Sanghi2022} and using the narrow-band F656N (H\ensuremath{\alpha}\xspace) filter to look for H\ensuremath{\alpha}\xspace-emitting planets. Accretion onto planets produces strong hydrogen line emission, reduces planet-star contrasts, and can in principle improve search efficiency. Any detections result in H\ensuremath{\alpha}\xspace luminosities (\ensuremath{L_\mathrm{H\alpha}}\xspace) and accretion rate estimates of candidate planets.

  The strategy of observing transition disk systems is motivated by models showing that the disk gaps and cavities are consistent with being sculpted by multiple giant planets \citep{Dodson2011,Close2020}. These models are broadly supported by the discoveries and characteristics of PDS 70 b and c \citep{Bae2019}. We identified ten systems that harbor giant cavities demonstrated in ALMA dust continuum images \citep{Francis2020}. The host stars are generally faint or located in the northern hemisphere and, consequently, are not easily accessible by ground-based visible-light adaptive-optics (AO) systems.

  \begin{figure*}[!th]
    \centering
    \includegraphics[width=\textwidth]{Figure1.pdf}
    \caption{Primary-subtracted images (upper) and the SNR maps (lower) of AB Aur observed by HST/WFC3/UVIS in the F656N filter. The left column shows the primary-subtracted natively sampled images (pixel scale\,=\,40~mas) and the right column shows the  reconstructed images (pixel scale\,=\,20~mas). The upper panels show the intensity maps (in unit of $\mathrm{e^{-}s^{-1}}$). These images, as well as the upper panels of Figure~\ref{fig:epoch}, the left panel of Figure~\ref{fig:injection}, and Figure~\ref{fig:detection} are smoothed using a Gaussian kernel with $\mbox{FWHM}=0.07''$ ($1.0\times$ PSF FWHM). Smoothing improves visualization by reducing high-frequency noise but is not applied in the sensitivity or contrast curve estimates. The lower panels are the SNR maps. For all four images, the color maps have a linear stretch. The central $r=0.3''$ radius circle is masked out, because in this region the detection sensitivity is poor due to high-contrasts and saturation. In the upper-left panel, an ellipse marks the transition disk cavity identified by \citet{Francis2020}. In the lower-left panel, gray dashed lines indicate the spiral disk structures (S1 to S4) detected in our images. The point-like source at $0.6''$ south of AB Aur is visible in all panels. }
    \label{fig:image}
  \end{figure*}

  \begin{figure*}[!th]
    \centering
    \includegraphics[width=\textwidth]{Figure2.pdf}
    \caption{Primary-subtracted images and SNR maps in Epochs 1 (left four panels) and 2 (right four panels). The intensity and SNR maps obtained from processing the natively sampled and reconstructed images are organized in the same way as Figure~\ref{fig:image}. The point source is consistently recovered in both epochs, although the detection SNRs are below 5 in the second epoch due to shorter integration time. Systematic residuals near the diffraction spikes are also prominent in the second epoch.}
    \label{fig:epoch}
    \end{figure*}

  The first target of the survey is AB Aurigae (AB Aur), a  Herbig Ae/Be star \citep{DeWarf2003} with a spectral type of A0, a mass of $2.4\pm0.2M_{\odot}$, and an isochronal age of $6.0^{+2.5}_{-1.0}$~Myr\footnote{The age is estimated with the PARSEC isochrones \citep{Bressan2012} and the observationally constrained $L_{\mathrm{bol}}$ and spectral type of AB Aur \citep{Herczeg2014}.}. It is actively accreting at a rate of $\dot{M}\approx10^{-7}M_{\odot}\mathrm{yr^{-1}}$ \citep{Garcia-Lopez2006} and exhibits significant variability (${>}10\%$) in both broadband photometry and the H\ensuremath{\alpha}\xspace line \citep{Harrington2007,Cody2013}. AB Aur hosts a transition disk that includes a highly structured inner disk with multiple spirals \citep[e.g.,][]{Boccaletti2020,Jorquera2022} and a cavity with an outer edge at 156~AU \citep{Francis2020}. Several spiral structures have been identified in the inner disk and their presence has been attributed to ongoing planet formation \citep{Fukagawa2004, Oppenheimer2008, Boccaletti2020}. Specifically, \citet{Tang2017} identified two spirals in high-resolution ALMA  $^{12}$CO $J{=}2{-}1$ emission observations and suggested that the spiral morphology could be explained by tidal disturbance caused by a companion located at 60 to 80 AU ($0.4''$ to $0.6''$ angular separation) and a ${\sim}180^{\circ}$ position angle (P.A.). Recently, \citet{Currie2022} presented evidence for an embedded protoplanet, referred to as AB Aur b, in visible and infrared images near the location predicted by \citet{Tang2017}. In this \emph{Letter}, we report results from two epochs of HST/WFC3 F656N (H\ensuremath{\alpha}\xspace) high-contrast imaging observations of AB Aur.

  \section{Observations}

  \subsection{HST Observations}
  AB~Aur\xspace was observed by HST/WFC3/UVIS on UT 2022-02-07 (Epoch 1) and 2022-03-28 (Epoch 2) for three (Orbits 1 to 3) and two orbits (Orbits 4 and 5), respectively. A 0.5-pixel (20~mas), four-point-box dithering strategy was adopted in all five orbits to enable the reconstruction of Nyquist-sampled images. In every orbit, HST sequentially pointed at each position and took eleven identical 2.7~s exposures in the F656N filter (H\ensuremath{\alpha}\xspace narrow-band; $\lambda_{\mathrm{c}}=6561.5\,\text{\normalfont\AA}$, $\Delta\lambda=17.9\,\text{\normalfont\AA}$). The exposures were captured by the \texttt{c512c} subarray, which has a native pixel scale of 40~mas and a field of view of $20''\times20''$.  Due to guide star acquisition failures, the first frame of Epoch 1 and the first two frames of Epoch 2 did not acquire the target. We discarded these images in our analysis. In total, the HST observations comprised 217 raw frames, amounting to 586~s of on-source exposure time.

  HST's roll angles (the V3 axis orientation) were $246.3^{\circ}$, $271.3^{\circ}$, $246.3^{\circ}$, $260.6^{\circ}$, and $260.6^{\circ}$  in Orbits 1 to 5, respectively. The angular differentials between any two distinct roll angles are $10.7^{\circ}$, $14.3^{\circ}$, and $25^{\circ}$. At a separation of $0.6''$, these correspond to spatial displacements of 108, 149, and 262~mas, or 1.6, 2.2, and 3.7 times the full width at half maximum (FWHM) of the F656N point spread function (PSF), respectively. These large azimuthal displacements enable the use of angular differential imaging (ADI, e.g., \citealt{Marois2006}) to subtract the PSF of the host star between any pair of images with different telescope rolls.

  The 2022-03-28 observations were a repeat of the first and third orbits of the 2022-02-07 observations that suffered from a fine-guiding sensor failure resulting in 0.1 to 0.3 pixel pointing errors. This issue impairs image reconstruction accuracy at small scales (${<}15\,\mbox{mas}$) but is not detrimental to this study. Therefore, we present results from all five orbits.

  \subsection{High-resolution Spectroscopy of AB Aur}

  The HST images do not constrain the continuum or resolve the H\ensuremath{\alpha}\xspace line. To precisely characterize the H\ensuremath{\alpha}\xspace emission from the host star, we obtained an optical \'echelle spectrum of AB Aur using the Tull Coud\'e spectrograph \citep{Tull1995} on the 2.7-m Harlan J. Smith Telescope at McDonald Observatory on UT 2019-09-26. The exposure time was 600~s and conducted with a $1.2''$ slit, resulting in a resolving power of $R=60,000$ and a wavelength span from $3870\,\text{\normalfont\AA}$ to $10,450\,\text{\normalfont\AA}$ in a total of 56 orders. The wavelength calibration was carried out with a ThAr emission lamp spectrum.

  The Tull spectrum is reduced with a custom pipeline. We adopt the order that spans from $6547.7\,\text{\normalfont\AA}$ to $6656.7\,\text{\normalfont\AA}$ and contains the H\ensuremath{\alpha}\xspace line. The spectrum is normalized by the continuum determined by the best-fitting fifth order polynomial of the line-excluded region ($\lambda<6550\text{\normalfont\AA}$ or $\lambda>6580\text{\normalfont\AA}$). We compute the synthetic photometry in the F656N filter using the \texttt{pysynphot} package and find that the band-average flux is 2.58 times of the continuum.

  \section{Data Reduction}

  \begin{figure*}
    \centering

    \includegraphics[width=\textwidth]{Figure3.pdf}
    \caption{The injection-and-recovery image (left) and the $5\sigma$ contrast curve (right).
      The injection-and-recovery test is conducted on the natively sampled images. Injected artificial planets with the same brightness as the one measured in the point-like source are consistently recovered and their recovered photometry is used in correcting companion flux loss in primary subtraction and calibrating the contrast curve. The inset in the left panel is the inner $r=0.7''$ circular region of the injection-and-recovery image, where the artificial sources have the same angular separations as AB Aur b.
      In the right panel, the $5\sigma$ contrast measurements are shown in squares, which are connected by a spline interpolation (solid line).
      The loss of sensitivity at $3.6''$ to $4.0''$ is due to a known filter ghost.}
    \label{fig:injection}
  \end{figure*}

  Our data reduction starts with the \texttt{flc} files downloaded from the HST archive. The nominal bias, dark, flat field calibrations, as well as charge transfer efficiency correction are performed by the \texttt{CalWFC3} pipeline. We first identify cosmic-ray (CR) affected pixels using the data quality flag (DQ=1024) and replace these pixels with two-dimensional linear interpolations of the respective neighboring pixels. Every set of four-point dithered and CR-corrected frames are processed by a Fourier-interlacing pipeline \citep{Lauer1999a,Zhou2021} to construct one $2\times$ finer sampled image. This pipeline is optimized for point-source PSFs and distinct from the drizzling method that is part of the regular data reduction performed at STScI.

  We then independently perform a complete set of primary subtraction procedures on the natively sampled (pixel scale${=}40$~mas) and reconstructed (pixel scale${=}20$~mas) image sets. First, the image cube is registered by the primary star centroids determined by the \texttt{photutils} package’s \texttt{centroid\_2dg} function \citep{Robitaille2013}. Then, the images are split into two annuli that have inner/outer radii of $0.3''/0.7''$ and $0.7''/2.5''$, respectively. We subtract the host star PSF in each annulus using the Karhunen-Lo\`eve Image Projection (KLIP, \citealt{Soummer2012}) algorithm with a custom pipeline. For each image, all frames obtained at a different roll angle are used as references and the number of KLIP components is equal to the number of reference frames. We have also experimented with limiting the KLIP components to 5, 10, and 15 and find that the variation has negligible effect on the subtraction results. After primary subtraction, we rotate all frames to align their $y$-axes to the true north and co-add the aligned frames to form the final image. The co-adding step is also performed for the two epochs individually. In total, we create six primary-subtracted images for two pixel scales each with three sets of co-adds (the combined set: Figure~\ref{fig:image}; individual epochs: Figure~\ref{fig:epoch}).

  We convert the primary-subtracted images into signal-to-noise ratio (SNR) maps for point-source detections. The SNR is defined as the ratio between the flux integrated in a 1~FWHM radius aperture and the standard deviation of fluxes along the remaining azimuthal region sampled by non-overlapping apertures. Biases induced by small-number statistics were corrected based on \citet{Mawet2014}. The SNR maps are shown in Figures~\ref{fig:image} and \ref{fig:epoch} under the respective primary-subtracted images. Substantial residuals appear near the diffraction spikes in the reconstructed images, which are related to the telescope pointing errors.

  A point-like source emerges in every PSF-subtracted image at an identical location ${\sim}600$~mas south of the primary star. We apply the KLIP forward modeling method \citep{Pueyo2016} on the natively sampled images to precisely determine the astrometry and photometry of this source and correct measurement biases introduced by primary subtraction. Synthetic WFC3 UVIS2 F656N PSFs are generated with the TinyTim package \citep{Krist2011} as the point-source model. We inject a negative-flux PSF to the original frames and optimize the P.A., separation, and flux of the injection so that the residual sum of squares in a 3-pixel radius circular aperture minimizes.

    To estimate the astrometric and photometric uncertainties, as well as the KLIP throughput, we conduct injection-and-recovery tests on the natively sampled images. Positive-flux TinyTim PSFs are injected at P.A.s of $0^{\circ}$, $\pm45^{\circ}$, $\pm90^{\circ}$, and $\pm135^{\circ}$ spanning eight equal-interval separations between $0.62''$ and $4.26''$ (Figure~\ref{fig:injection}). These PSFs have the same flux as the point-like source. We subtract the primary PSF  and then measure the positions and fluxes of the injected sources. At a specified separation, the average ratio between the injected and recovered fluxes is adopted as the throughput, and the standard deviations of the recovered fluxes, P.A.s, and separations are adopted to be the photometric and astrometric uncertainties.

    Figure~\ref{fig:injection} includes the $5\sigma$ contrast curve that characterizes our detection sensitivity at angular separations spanning $0.4''$ to $4.5''$. In deriving the contrast curve,  we calculate the mean and standard deviation of the fluxes in non-overlapping 1-FWHM-radius apertures centered at a given separation. At 600~mas, the resolution element containing the detected point source is excluded. We follow the SNR definition in \citet{Mawet2014} to convert the mean and standard deviation into a flux corresponding to $\mbox{SNR}=5$ and normalize this value by the flux of AB Aur.  Finally, we divide the contrast by the KLIP throughput determined by positive PSF injections (Figure~\ref{fig:injection}) to correct for flux loss in primary subtraction.

  \section{Results}

  We detect a point-like source near AB Aur in the primary-subtracted image (Figure~\ref{fig:detection}). In the natively sampled frames, the detection SNRs are 5.4, 3.8, and 6.4 in Epoch 1, Epoch 2, and the combined image, respectively; in the reconstructed images, the SNRs are 5.6, 4.4, and 7.8. The consistent results in natively sampled and reconstructed images, as well as in individual epochs, further strengthen the confidence in the detection. The source is located at a separation of $600\pm22$~mas ($93\pm3$~AU in projected physical distance) away from AB~Aur with a P.A. of $182.5^{\circ}\pm1.4^{\circ}$. This position is consistent with the protoplanet candidate AB Aur b reported in \citet{Currie2022}. Photometry on the combined image yields $64\pm16$ $\mathrm{e}^{-}$/s, corresponding to a band-averaged flux density of $F_{\nu}=1.5\pm0.4$~mJy or $F_{\lambda}=1.0\pm0.3\times10^{-15}$~erg/s/cm$^{2}$/$\text{\normalfont\AA}$. Fitting a 2D Gaussian to the source yields a FWHM of 150~mas in the tangential direction and 127~mas in the radial direction, which are 2.1 and 1.8 times the FWHM of a F656N PSF. The FWHMs of recovered artificial companions are between 85 to 120 mas. The larger size of the detected source suggests it may trace slightly extended emission.

  The H\ensuremath{\alpha}\xspace emission from the companion is consistent between the two epochs separated by 50 days, but our constraint on variability is weak. Photometric results of the individual epochs are $55\pm18$ $\mathrm{e}^{-}$/s and $76\pm30$ $\mathrm{e}^{-}$/s. The corresponding flux densities are $F_{\nu,1}=1.3\pm0.4$~mJy and $F_{\nu,2}=1.8\pm0.7$~mJy or $F_{\lambda,1}=0.9\pm0.3\times10^{-15}$~erg/s/cm$^{2}$/$\text{\normalfont\AA}$ and $F_{\lambda,2}=1.2\pm0.5\times10^{-15}$~erg/s/cm$^{2}$/$\text{\normalfont\AA}$. Based on the joint flux uncertainty between the two epochs, we reject variability with amplitude greater than 50\%. Using the same HST images, we precisely determine the H\ensuremath{\alpha}\xspace variability of the star AB~Aur to be $13.1\pm0.2\%$, consistent with the known variability of the host star \citep{Harrington2007,Cody2013}.

  Four spiral-like scattered-light features are detected at $\mathrm{SNR}$ values between 2 to 3. In Figure~\ref{fig:image}, we mark these structures as S1 to S4. These spiral features have all been previously identified in high-contrast imaging observations of AB Aur \citep[e.g.,][]{Fukagawa2004, Boccaletti2020}. Notably, the leading knots of S1 and S2 point towards the position of the point-like source.

  \begin{figure}[!t]
    \centering
    \includegraphics[trim=1cm 0cm .5cm 0cm,clip,width=\columnwidth]{Figure4.pdf}
    \caption{A zoomed-in ($2''\times2''$) natively sampled primary-subtracted image highlighting the detection of the point-like source. The source is at a $\mbox{P.A.}$ of $182^\circ$ and a separation of $600$~mas relative to the host star. This detection is at the location of AB Aur b, a candidate embedded planet recently reported by \citet{Currie2022}.}
    \label{fig:detection}
  \end{figure}

  We do not detect any sources at $2.75''$ or $3.72''$ reported in \citet{Currie2022}. The $5\sigma$ upper limits in our combined WFC3 frames are $0.098$~mJy and $0.091$~mJy, respectively, and are 4.8 and 22 times the flux measured by HST/STIS \citep{Currie2022}. Therefore, our observations do not have sufficient sensitivity to determine the nature of the two sources.

  \begin{figure*} [!ht]
    \centering
    \includegraphics[width=\columnwidth]{Figure5a}
    \includegraphics[width=\columnwidth]{Figure5b}
    \includegraphics[width=\columnwidth]{Figure5c}
    \includegraphics[width=\columnwidth]{Figure5d}
    \caption{The SED of AB Aur b (upper left), PDS 70 b (upper right), a comparison of the two SEDs (lower left), and a comparison to the stellar H\ensuremath{\alpha}\xspace{} emission of AB Aur (lower right). The SED of AB Aur b includes optical to NIR measurements from \citet[the blue squares and solid line]{Currie2022} and the H$\alpha$ flux from the combined HST data in this work (orange square). All AB Aur b fluxes are de-reddened based on $A_{V}=0.5$~mag and $R_{V}=3.1$. The best-fitting blackbody (3800~K) is shown in a green-dashed line. In contrast, the NIR SED of PDS~70~b (gray circles, data source: \citealt{Muller2018,Mesa2019a,Wang2020,Stolker2020}) is best-fit by a 1200~K blackbody (black dashed line). In the bottom-left panel, we normalize the two SEDs by their respective H$\alpha$ flux, so the H\ensuremath{\alpha}\xspace points (orange square for AB Aur b, black circle for PDS 70 b) overlap. AB Aur b is bluer and has a stronger optical continuum. In the lower right panel, flux is normalized to the continuum (for AB Aur b, it is the HST/STIS flux). The H\ensuremath{\alpha}\xspace of AB Aur b is consistent with the stellar value within $1\sigma$. }
    \label{fig:SED}
  \end{figure*}

  \section{Discussion}

  Rings and spirals in transition disks may introduce ambiguous protoplanet signals in high-contrast imaging observations. Emission from these structures may appear clumpy and even point-like after PSF-subtraction \citep{Follette2017,Rameau2017}. For example, the protoplanet candidates LkCa 15 bcd \citep{Kraus2012,Sallum2015} and HD 100546 bc \citep{Currie2015} were later found to be consistent with disk features (LKCa 15: \citealt{Currie2019}; HD 100546: \citealt{Follette2017} and \citealt{Rameau2017}). Because H\ensuremath{\alpha}\xspace emission from an accreting planet distinguishes the planet from the disk in their spectral energy distributions (SED), H\ensuremath{\alpha}\xspace images can provide an independent way of identifying protoplanets. For PDS 70 b and c, their H\ensuremath{\alpha}\xspace flux densities are over 100 times that of the estimated optical continuum, which help confirm these objects as accreting protoplanets. In contrast, in the HD~100546 case, non-detections of H\ensuremath{\alpha}\xspace-emitting point sources call into question the candidates as accreting protoplanets \citep{Follette2017,Rameau2017}. For AB Aur b, \citet{Currie2022} presented several lines of evidence to support the protoplanet identification: orbital motion has been detected, and the source's spectrum and polarimetric intensities differ from the rest of the disk. However, the companion resides in a region of extended scattered light and photospheric features appear to be absent in the SED \citep{Currie2022}. Below, we discuss whether the H\ensuremath{\alpha}\xspace detection adds support to the protoplanet interpretation.

  The H\ensuremath{\alpha}\xspace emission at the position of AB Aur b has been detected by Subaru/VAMPIRES \citep{Currie2022} and HST/WFC3 (this work). However, it is unclear whether the VAMPIRES detection is physically associated with AB~Aur~b due to uncertainties in the instrument's astrometric calibration \citep{Currie2022}. The astrometry of the WFC3 detection is well-calibrated and consistent with measurements in other bands. Therefore, we only adopt the WFC3 measurement for this discussion. The WFC3 band-averaged H\ensuremath{\alpha}\xspace flux density of $f_{\nu}{=}1.5\pm0.4$~mJy is higher than the optical continuum measured by HST/STIS ($f_{\nu}{=}0.47\pm0.05$~mJy; $\lambda{=}0.57 \mu$m, $\Delta\lambda{=}0.43\,\mu$m) and Subaru/VAMPIRES ($f_{\nu}{=}1.13\pm0.37$~mJy; $\lambda{=}0.647 \mu$m), but the difference is moderate to marginal, respectively. The line-to-continuum flux ratio (defined as the band-averaged F656N flux divided by the continuum flux) is $3.2\pm0.9$ or $1.3\pm0.6$, depending on whether the STIS or the VAMPIRES flux is adopted for the continuum, much lower than the flux ratios ${>}100$ estimated for the PDS 70 planets, {as well as several accreting planetary mass companions \citep{Zhou2014,Eriksson2020,Stolker2021}}. AB Aur b's low line-to-continuum flux ratio is inconsistent with the accretion shock models that explains the observations of PDS 70 b and c \citep[e.g.,][]{Aoyama2018}, which produce pronounced H\ensuremath{\alpha}\xspace lines.  As a result, the H\ensuremath{\alpha}\xspace detection does not strengthen the interpretation that AB Aur b is an accreting protoplanet.

If the observed H\ensuremath{\alpha}\xspace emission is indeed from an accretion shock, the shock also produces strong optical continuum emission that reduces the line-to-continuum flux ratio. Follow-up characterization of this source to measure ultraviolet and optical excess emission would provide useful comparisons to model SEDs and enable a diagnosis of the accretion shock interpretation \citep[e.g.,][]{Zhu2015}. Assuming that the observed F656N flux is entirely produced by the accretion shock of AB Aur b, we can determine the instantaneous mass accretion rate of the protoplanet based on the observed H\ensuremath{\alpha}\xspace luminosity.  Adopting a line-of-sight extinction of $A_{V}=0.5$~mag (the same as the star, \citealt{Garcia-Lopez2006}) and assuming no additional extinction opacity, we find a H\ensuremath{\alpha}\xspace line luminosity of $\ensuremath{L_\mathrm{H\alpha}}\xspace=2.2\pm0.7\times10^{-5}\ensuremath{L_{\odot}}\xspace$. At this \ensuremath{L_\mathrm{H\alpha}}\xspace, gas in the accretion flow should not lead to significant absorption \citep{Marleau2022}. The \ensuremath{L_\mathrm{H\alpha}}\xspace value corresponds to a total accretion luminosity of $\log(\ensuremath{L_\mathrm{acc}}\xspace/\ensuremath{L_{\odot}}\xspace)=-2.8\pm0.3$ by assuming the planetary surface shock model \citep[][]{Aoyama2018,Aoyama2021} or $\log(\ensuremath{L_\mathrm{acc}}\xspace/\ensuremath{L_{\odot}}\xspace)=-3.5\pm0.3$ using the empirical relation for classic T Tauri stars \citep{Alcala2017}, respectively. Adopting a planetary mass of $9M_{\mathrm{Jup}}$ and radius of $2.7R_{\mathrm{Jup}}$ \citep{Currie2022}, these accretion luminosities can be translated to mass accretion rates of $\log(\dot{M}/(M_{\mathrm{Jup}}\mathrm{yr^{-1}}))$ between $-5.8$ and $-6.6$. These estimates roughly agree with the result $\log(\dot{M}/(M_{\mathrm{Jup}}\mathrm{yr^{-1}}))=-6.0$ obtained from SED fitting \citep{Currie2022}.

Alternatively, the observed H\ensuremath{\alpha}\xspace flux can also be stellar light scattered by a compact disk structure or an envelope surrounding the protoplanet. In this case, the optical SED of the candidate protoplanet, which requires additional multi-band photometry to precisely determine, should mimic the SED of the star.  Remarkably, the line-to-continuum ratios of AB Aur (2.58) and AB Aur b ($3.2\pm0.9$, adopting HST/STIS photometry as the continuum) are consistent within $1\sigma$ (Figure~\ref{fig:SED}, bottom right panel). We note that this agreement might be serendipitous because AB Aur's H\ensuremath{\alpha}\xspace emission is variable \citep{Harrington2007} and our comparison is based on asynchronous observations (the stellar spectrum is obtained 15 months earlier than the HST images). Our HST observations do not measure continuum flux and thus cannot simultaneously constrain the stellar line-to-continuum flux ratios.  Continuous and simultaneous H\ensuremath{\alpha}\xspace monitoring of AB Aur and AB Aur b can help test this scattering scenario. \citet{Currie2022} concluded that the scattered light could not account for the infrared emission, because AB Aur b was undetected in infrared polarized images. If the H\ensuremath{\alpha}\xspace flux is scattered light, the disk structure or dust envelope should predominantly scatter at short wavelengths such that scattered stellar emission is only significant in optical bands.

Finally, we emphasize the striking difference in the SED shape between AB Aur b and the PDS 70 planets. In Figure~\ref{fig:SED}, we compare the SEDs of AB Aur b and PDS 70 b and find that AB Aur b is significantly brighter and bluer. Fitting single blackbody models to the SEDs of AB Aur b and PDS 70 b (excluding the H\ensuremath{\alpha}\xspace points) yield effective temperatures/radii of 3800~K/1.6~$R_{\mathrm{Jup}}$ and 1200~K/2.7~$R_{\mathrm{Jup}}$, respectively. These results correspond to over 1~dex difference in their bolometric luminosities ($\log(L_{\mathrm{bol}}/L_{\odot})=-2.3$ for AB Aur b \footnote{This value is higher than the one ($\log(L_{\mathrm{bol}}/L_{\odot})=-2.695\pm0.095$) reported by \citet{Currie2022}, because we adopted a blackbody rather than a photospheric model SED to calculate $L_{\mathrm{bol}}$.} and $\log(L_{\mathrm{bol}}/L_{\odot})=-3.8$ for PDS 70 b). We can compare the bolometric luminosity of AB Aur b to evolutionary track predictions while safely ignoring contributions from an accreting circumplanetary disk, because the expected disk luminosity is at least one order of magnitude lower\footnote{$\log(L_{\mathrm{CPD}}/L_{\odot})\approx-3.5$. This is estimated using Equation (5) of \citet{Zhu2015}, assuming $M=9M_{\mathrm{Jup}}$, $R_{\mathrm{in}}=7.5R_{\mathrm{Jup}}$ \citep{Currie2022}, and $\log(\dot{M}/(M_{\mathrm{Jup}}\mathrm{yr^{-1}})=-5.8$}. Assuming a hot-start model \citep[e.g.,][]{Burrows1997,Mordasini2017} and that AB Aur b has a planetary mass ($M{<}13M_{\mathrm{Jup}}$), AB Aur b has to be younger than 3~Myr to have such a high $L_{\mathrm{bol}}$. In contrast, $L_{\mathrm{bol}}$ of the PDS 70 planets are consistent with $1$ to $4$~$M_{\mathrm{Jup}}$ planets at an age of ${\sim}$5~Myr \citep{Wang2020,Stolker2020}. Because the isochronal age of AB Aur ($6.0^{+2.5}_{-1.0}$~Myr) is older than 3~Myr, the planetary interpretation of AB Aur b and the inferred young age might imply delayed planet formation in the AB Aur disk.

  \section{Summary}

  1. We observed the Herbig Ae/Be star AB Aur using HST/WFC3/UVIS in the F656N (H\ensuremath{\alpha}\xspace) band in two epochs separated by 50 days. After subtracting the primary-star PSF, we detected a point-like source $0.6''$ away from the star in both epochs. In the most optimal reduction, the detection SNRs are 5.6, 4.4, and 7.8 in Epoch 1, Epoch 2, and the combined data-sets, respectively.

  2. The P.A. and separation of the companion are $182.5^{\circ}\pm1.4^{\circ}$  and $600\pm22$~mas relative to the host star. This location is consistent with the astrometry of AB Aur b, a recently reported protoplanet candidate. The candidate companions at separations of $2.75''$ and $3.72''$ are below our detection limits.

  3. The band-averaged H\ensuremath{\alpha}\xspace flux densities are $0.9\pm0.3$, $1.2\pm0.5$, and $1.0\pm0.3$~$\times10^{-15}$erg/s/cm$^{2}$/$\text{\normalfont\AA}$ ($F_{\lambda}$) or $1.3\pm0.4$, $1.8\pm0.7$, $1.5\pm0.4$~mJy ($F_{\nu}$) in Epoch 1, Epoch 2, and the combined images, respectively.

  4. In comparison to PDS 70 b and c, the H\ensuremath{\alpha}\xspace-to-continuum flux ratio of AB Aur b is significantly lower, suggesting that an accretion shock that produces strong H\ensuremath{\alpha}\xspace line that has flux density over 100 times the continuum is not visible from AB Aur b. Both planetary accretion and scattered stellar emission from either a compact disk structure or an envelope surrounding AB Aur b can contribute to the H\ensuremath{\alpha}\xspace flux, and the origin of the H\ensuremath{\alpha}\xspace emission remains unconstrained.

  5. Assuming that the observed H\ensuremath{\alpha}\xspace emission is entirely powered by accretion onto the protoplanet and there is no gas or dust extinction from the disk nor from the accretion flow itself, we estimate an mass accretion rate of $\log(\dot{M}/M_{\mathrm{Jup}}\mathrm{yr^{-1}})=-5.8$ to $-6.6$, depending on the assumed $\ensuremath{L_\mathrm{H\alpha}}\xspace-\ensuremath{L_\mathrm{acc}}\xspace$ relations. The accretion rate estimate suffers from significant systematic uncertainty due to the unknown origin of the H\ensuremath{\alpha}\xspace emission and unconstrained extinction.

  The detection of H\ensuremath{\alpha}\xspace emission from the candidate protoplanet AB~Aur b offers exciting opportunities to  understand the mass assembly processes of giant planets forming in a highly structured disk. These interpretations are conditioned on definitively resolving the origin of the H\ensuremath{\alpha}\xspace emission. Disentangling the planetary accretion versus stellar scattering of H\ensuremath{\alpha}\xspace emission possibilities could not only prove to be crucial in understanding the origin of AB Aur b, but also for the interpretation of future discoveries of protoplanets embedded in the disks of their host stars.

  \vspace{2em}
  We thank the referee for a constructive report that improves the rigor of the manuscript. We acknowledge excellent observing support from STScI staffs, especially J. Green and A. Armstrong. We thank Thayne Currie for the CHARIS spectrum and Gabriel-Dominique Marleau for helpful comments. Y.Z. acknowledges support from the Heising-Simons Foundation 51 Pegasi b Fellowship. B.P.B. acknowledges support from the National Science Foundation grant 450 AST-1909209, NASA Exoplanet Research Program grant 451 20-XRP20\_2-0119, and the Alfred P. Sloan Foundation. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.  The observations and data analysis works were supported by program HST-GO-16651. Supports  for  Program  numbers  HST-GO-16651 were provided by NASA through a grantfrom the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

  \software{Astropy \citep{Robitaille2013}, Numpy \citep{Walt2011}, Scipy \citep{Virtanen2020}, Matplotlib \citep{Hunter2007}, pyKLIP \citep{Wang2015}}