\section{Introduction}
Since the detection of the first molecules in the interstellar medium (ISM; \citealt{swings1937,douglas1941}), there has been an increasing interest in the study of chemistry in astronomical contexts, up to the point that, nowadays, the study of both chemical structure and evolution of any given system under scrutiny is widely accepted as crucial to completely understand it. Unfortunately, the inclusion of self-consistent chemical evolution in theoretical models in the aim at further constraining the processes happening in astronomical environments is very complicated. Chemical reactions work on very short time-scales compared to other dynamical processes, and, on top of this, they strongly depend on the local properties (density and temperature) of the medium and the impinging radiation flux.

From a numerical perspective, an accurate treatment of any system evolution as a whole would require on-the-fly (OTF) non-equilibrium chemistry calculations in hydrodynamic simulations, especially when the chemical state of the gas affects the thermodynamics of the system, for instance on galactic scales \citep[see, e.g.][]{2014MNRAS.440.3349R,2016A&A...590A..15B}. While most of the galaxy-scale studies to date are still based on the post-processing of the simulations, either via pre-computed \textsc{cloudy} \citep{ferland1998} tables \citep{pallottini2017b,olsen2017,pallottini2019,katz2019,arata2020}, or determined from chemical networks under the assumption of photo-ionisation equilibrium \citep{keating2020}, the last decade has seen an increasing effort by several groups to properly include non-equilibrium chemistry in simulations, from simple models tracing primordial species \citep{gnedin2009,christensen2012,tomassetti2015,katz2017,pallottini2017b,lupi2018,nickerson2018} to more accurate ones that also include metal elements \citep{Glover2010,Hu2016, richings2016,capelo2018,lupi2020a,lupi2020b}.

Conversely to galactic scales, the chemistry complexity exponentially grows when simulating the small scales typical of pre-stellar cores and proto-planetary disks. During the initial phases of the star formation process, pre-stellar cores are characterised by extremely high densities ($n_{\rm H_2} \gtrsim 10^4$ cm$^{-3}$) and low temperatures ($T<20$ K), and exhibit a well-defined chemical structure determined by characteristic chains of chemical reactions \citep{caselli2012,2021PhR...893....1O}. These conditions are ideal for the adsorption of heavy species on the surface of dust grains, a process known as freeze-out \citep{bacmann2002,bacmann2003,2006A&A...455..577T}, and thus a plethora of gas-grain interactions and reactions on the surfaces of such grains starts to take place, further increasing the chemical complexity of these objects \citep[see, e.g.,][]{2018IAUS..332....3V}. As a consequence, chemical networks exponentially grow, making on-the-fly calculations in already expensive 3D magneto-hydrodynamic (MHD) simulations computationally unfeasible. Relevant examples of still unresolved features due to the lack of robust and consistent chemical modelling are, for instance, the high-abundances of NH$_3$ observed towards the central regions of pre-stellar cores \citep{2004A&A...416..191T,crapsi2007,sipila2019}, or the convoluted processes governing the formation and evolution of interstellar complex organic molecules (iCOMs) \citep{jimenez-serra2016,punanova2018}.

The need to include gas-grain interaction and surface chemistry in the study of star-forming regions becomes then fundamental for the complete understanding of particular features and chemical complexities observed towards these objects. The common approach, in these cases, is to build low-dimensionality models, which are qualitatively useful to explore the different formation paths of a given iCOM or to perform sensitivity studies, by mapping a large parameter space. To date, there are no works where the formation of iCOMs is consistently followed in three-dimensional simulations, and the available studies rely instead upon the post-processing of the simulation outputs with a non-equilibrium chemistry solver \citep[see, e.g.][]{ruaud2018,coutens2020}. This approach is generally based on the selection of a sub-sample of resolution elements (particles) from standard pure MHD simulations, i.e. performed without OTF chemistry, and on the subsequent integration of the chemical equations on the already pre-determined dynamical history of said particles. On the other hand, \citet{priestley2018,priestley2019} suggested an even more convenient post-processing approach based on the parametrisation of the density structure of a collapsing pre-stellar core as a function of both time and radius, facilitating the inclusion of more complex dynamical processes, such as ambipolar diffusion, which become more important during the process of star formation. However, forcing the electron density to change as the solely function of dynamics without including any chemistry, as done by \citet{priestley2019}, can introduce large errors in the subsequent post-processing analysis.

Another important limitation is that the aforementioned techniques are usually applied to isothermal simulations. When the thermodynamics of the gas is taken into account, indeed, a minimal network of chemical reactions should always be included, to guarantee the correct thermal behaviour, which in turns depends on key species abundances (coolants), as for instance C$^+$ and CO. In these cases, it is unavoidable to couple OTF chemistry with dynamics. As an example, a minimal network able to follow CO formation/destruction and to give the right thermal behaviour in molecular cloud simulations, already requires $\sim$40 species and hundreds of reactions \citep[see e.g.][]{2017MNRAS.466.1259G}. Once the thermal behaviour is correctly assessed, post-processing can be still pursued to include more complex processes, like the formation of specific iCOMs, or isotopologues and isomers of specific tracers.

Overall, post-processing techniques have the ability to explore the impact of different sources of chemical degeneracy on the evolution of specific tracers, without incurring in computationally intensive problems, but, despite the different results obtained so far, a proper assessment of their reliability and limitations, performed against self-consistent MHD simulations with OTF chemistry, is still missing.  With this in mind, we assess here the reliability of post-processing chemistry in 3D MHD simulations in the context of star-forming collapsing cloud fragments, by directly comparing it to the results obtained by MHD simulations including coupled OTF non-equilibrium chemistry calculations. Thanks to the latter, we have the unique opportunity to evaluate the uncertainties associated with the post-processing. For the first time, we can constrain parameters like the number of resolution elements in the sub-sample, numerical approximations in the interpolation schemes, and the time-discretisation for the dynamical properties update.

The paper is organised as follows. In Section~\ref{sec:methodology}, we briefly summarise the setup used and describe how chemistry calculations are performed in the post-processing of the simulation. The discussion of the validity and limitations of the post-processing approach relative to OTF chemistry calculations is presented in Sections~\ref{sec:validation} and~\ref{sec:comparison}, respectively. Finally, in Section ~\ref{sec:conclusions} we draw our conclusions.

\section{Methodology}
\label{sec:methodology}
Here, we describe the setup of the simulation considered in this work, providing a summary of the initial conditions, and the details of the chemical network employed, and then introduce the post-processing method we apply.

\subsection{The reference MHD simulation}
As a reference case for this work, we consider the slow-collapse `M1' core from the simulation suite presented in \citet{Bovino2019}, a 3D MHD simulation of the isothermal collapse of a high-mass pre-stellar core. The assumption of an isothermal gas helps us further facilitate the treatment of chemistry in post-processing, not requiring any `feedback' of the temperature resulting from the chemical calculations on to the gas dynamics. The simulation was performed with the MHD code \textsc{gizmo} \citep{hopkins2015,Hopkins2016}, using its meshless finite-mass method and the standard cubic-spline kernel with an effective number of neighbours of 32.

The initial conditions for the simulation consist of a $M_{\rm BE}=20\rm\, M_{\odot}$ core following a Bonnor-Ebert density profile \citep{Ebert1955,Bonnor1956} with core radius $R_{\rm BE} = 0.17$ pc and an homogeneous temperature $T = 15$ K. The average number density of the core is $\langle n \rangle \equiv (3M_{\rm BE})/(4\upi R_{\rm BE}^3) = 2.21\times 10^{4}$ cm$^{-3}$, corresponding to an average free-fall time of $t_{\rm ff} = 260$ kyr. The number of particles (from both the core and the background gas) was set to $6\times 10^5$, corresponding to a mass resolution of $3.33\times10^{-5}\,{\rm M_\odot}$.

Alongside self-gravity, our simulation also includes ideal MHD, with the magnetic field aligned along the $z$-direction, and scaling with the perpendicular distance $R_\perp$ (in the $xy$ plane) as
\begin{equation}
    B_z(R_\perp)=B_0\left(1 + \frac{R_\perp}{R_{\rm BE}}\right)^{-\kappa},
\end{equation}
where B$_z$ is the component of the magnetic field parallel to the $z$-axis, $k=1.5$ an input scaling exponent, and $B_0$ the magnetic field strength at $R_\perp=0$. The latter is obtained by enforcing a user-defined mass-to-magnetic flux ratio $\mu\equiv M/\Phi = 2$, commonly used in simulations of pre-stellar cores \citep[see, e.g.][]{bastian2017,Kortgen2018,Goodson2016}. This choice yields an average magnetic field strength in our core of $\langle B \rangle = 46$ $\mu$G.

Finally, the core was assumed to be stirred by an initial turbulent velocity field following a Burgers-like power spectrum, with $E(k)\propto k^{-2}$ (see \citealt{bastian2017}). The turbulent field was normalised according to the assumed temperature $T$ and the desired Mach number $\mathcal{M} = 3$. The main parameters of the run are summarised on Table~\ref{tab:ICs1}.
\begin{table}
    \centering
    \caption{Fiducial values of important physical parameters used in the validation of the post-processing methodology.}
    \label{tab:ICs1}
    \begin{tabular}{ll}
    \hline Parameter & Value \\\hline
    \multicolumn{2}{c}{\textbf{Core properties}}\\
    Core mass $M_{\rm BE}$ & 20 M$_\odot$\\
    Core radius $r_{\rm BE}$ & 0.17 pc\\
    Average volume density $\langle n \rangle$ & 2.21$\times10^{4}$ cm$^{-3}$\\
    Mach number $\mathcal{M}$ & 3\\
    Virial parameter $\alpha_{\rm vir}$ & 4.32\\
    Average magnetic field $\langle B \rangle$ & 46 $\mu$G\\
    Free-fall time $t_{\rm ff}$ & 260 kyr\\
    \\
    \multicolumn{2}{c}{\textbf{Other parameters}}\\
    Temperature $T^{a}$ & 15 K \\
    CR ionisation rate per H$_2$ $\zeta_2$ & $2.5\times10^{-17}$ s$^{-1}$\\
    Average grain size $a$ & 0.035 $\mu$m \\
    Dust grain density $\rho_0$ & 3.0 g cm$^{-3}$ \\
    Dust-to-gas mass ratio $\mathcal{D}$ & $7.09\times10^{-3}$ \\
    Mean molecular weight $\mu$ & 2.4\\\hline
    \end{tabular}
    \\\flushleft $^a$ We assume $T_{\rm gas} = T_{\rm dust} = T$ throughout this work.
\end{table}

The chemical network employed in the simulation (and in our post-processing) is the same of \citet{Bovino2019,Bovino2020}. The network was created to explore the evolution of simple, yet fundamental C-N-O bearing tracers and deuterated species during the collapse of massive star-forming cores and clumps, following, for the first time, time-dependent freeze-out processes of key species in a 3D MHD setup. The species were initialised with the abundances reported in the second column of Table~\ref{tab:ICs2}, originally taken from \citet{sipila2015a}, but adapted to the molecular environment considered in our setup. In detail, species heavier than helium were no longer initialised in their atomic form, with the only exception being nitrogen, which was initialised half-atomic, half-molecular, due to the uncertainties about its main chemical state in dense environments. Our network contains mainly gas-phase chemistry, but also incorporates depletion and non-thermal desorption reactions for O, N, CO, and N$_2$. The adopted binding energies for these species are taken from \citet{Wakelam2017}, and are shown in the third column of Table~\ref{tab:ICs2}.
\begin{table}
    \centering
    \caption{Fiducial initial abundances $n$ of the species included in our network relative to the atomic hydrogen abundance $n_{\rm H}$. For the species that go through time-dependent depletion (see text), we also report the corresponding binding energy in the third column \citep[see][]{Wakelam2017}. Compared to the H$_2$ ortho-to-para ratio (OPR) usually assumed ($\leq0.1$), in this work we start from a more conservative value of 3.}
    \label{tab:ICs2}
    \begin{tabular}{lll}
    \hline Species & $n/n_{\rm{H}}$ & $E_{\rm B}$ [K] \\\hline
    oH$_2$ & $3.75\times10^{-01}$ & -- \\
    pH$_2$ & $1.25\times10^{-01}$ & -- \\
    HD & $1.50\times10^{-05}$ & -- \\
    He & $1.00\times10^{-01}$ & -- \\
    oH$_3^+$ & $1.80\times10^{-10}$ & -- \\
    pH$_3^+$ & $3.00\times10^{-09}$ & -- \\
    N & $1.05\times10^{-05}$ & 720 \\
    N$_2$ & $5.25\times10^{-06}$ & 1100 \\
    CO & $1.20\times10^{-04}$ & 1300 \\
    O & $1.36\times10^{-04}$ & 1600 \\
    GRAIN0 & $5.27\times10^{-11}$ & -- \\\hline
    \end{tabular}
\end{table}

\subsection{Post-processing hydrodynamic simulations with non-equilibrium chemistry}
\label{sec:postproc_method}
The computational cost of non-equilibrium chemistry in simulations poses significant limitations to its applicability, and evolving chemistry in post-processing can provide a viable cheaper alternative. Obviously, if all the particles in the simulation are used, we do not expect a large gain in terms of computational time, hence a more efficient solution is to only evolve a fraction of the gas resolution elements and then, if desired, interpolate the other elements according to the closest neighbour element in density (and/or any other relevant parameter) space.

The post-processing method is based on the assumption that chemistry does not significantly affect the dynamical evolution of the system. While this is certainly true under the isothermal assumption made in this work, relaxing this condition would necessarily require on-the-fly chemistry with at least a minimal network to track the thermal evolution of the gas \citep[see, e.g.,][]{2017MNRAS.466.1259G}.  We thus start from a `standard' hydrodynamic simulation of our system, i.e. without chemistry, for which we collect the outputs in which the history of its dynamical properties is stored.
On these outputs, we then apply the non-equilibrium chemical evolution, following the procedure reported here, and also summarised in the scheme of Fig.~\ref{fig:flowchart}:
\begin{enumerate}
    \item we extract the density of every particle in the run from the simulation outputs, and build the corresponding density evolution history (by tracking the particle ID across multiple snapshots);
    \item we randomly select a subset $N'=f_{\rm part}N$ of the $N$ particles in the simulation initial conditions (i.e. $t=0$), with $f_{\rm part}\leq 1$ a user-defined parameter. In order to properly sample the entire dynamic range of the simulation, independent of the total number of particles used to map the pre-stellar core and the background, we also enforce a user-defined number fraction $\xi_{\rm core}$ to be extracted from the particles sampling the core and the remaining ones from the low-density background surrounding it;\footnote{In principle, $\xi_{\rm core}$ should be tuned to the actual number of particles employed to map the core and the background in the simulation, respectively. We notice that, neglecting $\xi_{\rm core}$, hence assuming an homogeneous sampling of the entire particle distribution, can result in the oversampling of the background, depending on the assumptions made in the simulations, as the minimum spatial resolution in the background region and the box size.}
    \item for the selected particles, we initialise the species abundances in the same way as the reference simulation with OTF chemistry, with the values reported in Table~\ref{tab:ICs2};
    \item we evolve the chemical abundances of the $N'$ subset with \textsc{krome}, assuming a coarse integration time-step (i.e. the time at which the chemical abundances are actually stored) $\Delta t = m \Delta t_0$, with $\Delta t_0$ the time separation between the snapshots, and $m\geq 1$; over each integration step, the density of the particle is kept constant at the value obtained\footnote{We note that, although we assume a constant density at each step, different choices could be made, for example linearly interpolating between subsequent values, and the procedure would still be valid. However, such interpolation would simply represent an effective higher temporal resolution, and would require shorter $\Delta t$ for the chemistry solver, since each integration step necessarily assumes a constant gas density.} from the evolution history of every particle in the simulation at the beginning of each time-step;
    \item we estimate the species abundances for the remaining ($1-f_{\rm part}$) particles in the simulation from the closest post-processed particle in density-position space.
    In detail, we bin the post-processed particle sub-sample as a function of density. Then, we iterate over all non-evolved particles, find the corresponding density bin, and identify the spatially closest evolved particle in the bin, finally associating its species abundances to the target particle.
    It is important to notice that, because of the intrinsic scatter in the particle density distribution resulting from the hydrodynamic evolution, the interpolation can be moderately affected by the particles chosen for the chemical evolution.
    In addition, the accuracy of the results will depend on the chosen interpolation scheme. For instance, a simple density-only interpolation (using the average abundances in each density bin) would be much cheaper, but at the expense of losing information about local effects, hence of worsening the accuracy of the interpolation scheme. Most importantly, both these approaches are valid as long as the gas is isothermal, whereas a temperature-dependent binning might be also necessary when the isothermality is not guaranteed (unless a 1:1 relation between density and temperature can be found).
    \item After each step, we create a new snapshot file of the simulation that also includes the species abundances, and we repeat the procedure from step (iv), using the resulting abundances as initial values for the new step.
\end{enumerate}
This procedure is applied up to the time when a sink particle is identified in the snapshot, in order to avoid any artificial effect that the removal of density can have on the chemistry calculations. We note, however, that the calculations could be performed even after sink particles formed, if the proper caveats are considered.
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{postproc_flowchart.pdf}
    \caption{Schematic flowchart of the presented post-processing method.}
    \label{fig:flowchart}
\end{figure}
\label{sec:subset}
\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{pp_maps_183kyr.pdf}
    \caption{Column density maps of CO (first row), N$_2$H$^+$ (second row), and N$_2$D$^+$ (last row) at $t\sim184$ kyr for the `M1' core in \citet{Bovino2019}. The first column corresponds to the reference simulation (with OTF chemistry), and the second, third and fourth one to PP3, PP2, and PP1, respectively. While the agreement is reasonable in all cases, with the global distribution and density interval being always reproduced, the post-processing smears out some important features found in the full simulation (see, as an example, the ``filamentary'' structures in the top and bottom part of the N$_2$D$^+$ maps, or the location of the local minima in the CO maps), with the PP3 case showing the largest differences, as expected.}
     \label{fig:postproc_maps}
\end{figure*}

\section{Post-processing validation}
\label{sec:validation}
The approach considered in this work has two main free parameters: the fraction of particles used for the actual chemical evolution $f_{\rm part}$ and the time interval for the density update $\Delta t$. In this section, we validate the post-processing technique by applying it to the aforementioned M1 slow-collapse core run presented in \citet{Bovino2019}, employing the same network they used for consistency, and show how the free parameters affect our results.
\begin{figure*}
     \centering
     \includegraphics[width=.9\textwidth]{pp_pfracs_RPs_183kyr.pdf}
     \caption{Radial profiles at $t\sim 180$ kyr of CO (top panels), N$_2$H$^+$ (middle panels), and N$_2$D$^+$ (bottom panels) column density (left column) and number density (right column), for the simulation with OTF chemistry (black solid lines), and post-processed chemistry, PP1 (using 10 percent of the particles;  green dashed lines), PP2 (using 1 percent; blue dot-dashed lines), and PP3 (using 0.1 percent; red dotted lines), respectively.}
     \label{fig:profiles_fpart}
\end{figure*}

\subsection{The sub-sample selection}
We start by considering how the fraction of particles taken from the simulation for the post-processing affects our results. We select three different post-processing (PP) cases, comparing them with our reference simulation: one taking 10 percent (PP1), another taking 1 percent (PP2), and a third taking 0.1 percent (PP3) of the simulation particles, respectively, and assume an integration time-step of 0.5~kyr, corresponding to the time interval between the snapshots in our hydrodynamic simulations, i.e. $\Delta t = \Delta t_0$. This means that, with the original snapshot having $6\times 10^5$ particles, the post-processing of these three validation cases is performed on 60000, 6000, and 600 particles, respectively. Furthermore, for all the analysis reported in the rest of this work, we assume $\xi_{\rm core}=0.5$.\footnote{We explored the effect of choosing different core-to-background particle ratios and, in this case, the differences were entirely negligible, both qualitatively and quantitatively.}

In Fig.~\ref{fig:postproc_maps}, we show a visual comparison of the CO, N$_2$H$^+$, and N$_2$D$^+$ column densities at $t\simeq180$~kyr. The MHD simulation with coupled on-the-fly (OTF) chemistry is shown in the first column, while the results of the post-processing are reported in the other columns: PP1 in the second one, PP2 in the third one, and PP3 one in the last one, respectively. Qualitatively, the global distribution is well recovered, but we notice that the post-processing is not always able to reproduce all the features in the simulation.

For example, for N$_2$D$^+$ and N$_2$H$^+$ some artificial features in the maps tend to appear/disappear in the post-processing (e.g. filaments, overdensities, clumps etc.). On the contrary, for CO, which goes through the strongly time- and density-dependent freeze-out process, we observe an overabundance in the central region where we expect high depletion, which is progressively corrected as $f_{\rm part}$ increases.

To better quantify how important these differences are, we report in Fig.~\ref{fig:profiles_fpart} the column (left-hand panels) and number density (right-hand panels) radial profiles of the same three species shown in Fig.~\ref{fig:postproc_maps}. The column density profiles are obtained by radially averaging at set radii the correspondent column density map, while for the number density profiles we obtain the average over spherical shells of the actual 3D particle distribution. Evaluating the error on the number densities helps understanding how much uncertainty the integration along the line-of-sight is bringing, and to better assess the real error on local quantities, rather than averaged ones. Our reference run is shown as a black solid line, PP1 as a green dashed one, PP2 as a blue dot-dashed one, and PP3 as a red dotted one. The profiles clearly show that species depending mostly on gas-phase reactions (N$_2$H$^+$ in this case) are well reproduced, with low relative errors, even when only 1 percent of the particles is used. Species that heavily depend on gas-grain interactions, on the other hand, like N$_2$D$^+$ and CO, are more sensitive to the number of particles chosen for the post-processing, although the errors remain small overall.
This might be especially relevant in regions where the density rapidly changes in time, on time-scales much shorter than $\Delta t_0$. Nevertheless, considering the huge dynamic range involved, we can conclude that the error is in general very small ($\lesssim 10\%$), hence making our post-processing method a viable and reasonably cheap alternative to on-the-fly calculations.

By looking more in detail at the differences between the number density and column density profiles, we can notice a similar behaviour in both number and column densities, that can be attributed to the incomplete tracing of the core density structure when using a low $f_{\rm part}$, that favours higher density regions, sampled with more particles. These errors are then consistently propagated to the column densities.

This, in general, means that the $f_{\rm part}$ choice does not directly affect the calculations of the chemical abundances per se, but using a low value will still amplify the errors because of the poor sampling of the structure under analysis. In addition, as already mentioned, the value of $f_{\rm part}$ does have a more important role on the task of reproducing ``observational'' features, with the species relying on gas-grain processes being slightly more affected in general that those controlled by mostly gas-phase reactions.

\subsection{Integration time effect}
\label{sec:timesteps}
\begin{figure*}
     \centering
     \includegraphics[width=.9\textwidth]{pp_Deltat_RPs_183kyr.pdf}
     \caption{Same as Fig.~\ref{fig:profiles_fpart}, but now comparing the integration time-step choices. Here, the cyan dot-dashed line corresponds to the PP4 case, while the magenta dotted line corresponds to the PP5 case.}
     \label{fig:profiles_times}
\end{figure*}
The other important free parameter in our method is the time interval for the density update $\Delta t$. In the previous section, we assumed $\Delta t\equiv\Delta t_0 = 0.5$~kyr for simplicity, even though, in principle, $\Delta t$ might be arbitrarily set to larger values. However, we should keep in mind that, because of the dependence of chemical reactions on gas density, a different choice for the time-interval between density updates could have a significant impact on the evolution. In order to assess how sensitive the post-processing technique is to this choice, and thus its reliability, we analyse two additional cases, where we fix $f_{\rm part}=0.1$: a case with a five times longer integration step ($\Delta t=5\Delta t_0$; PP4), and a second one using an eight times longer one ($\Delta t=8\Delta t_0$; PP5). We stress that the density used at each time-step corresponds to the value obtained from the particle evolutionary history at the beginning of the time-step, in the aim at mimicking a coarser time sampling of the simulation outputs.

In Fig.~\ref{fig:profiles_times}, we show radial profiles of the number and column densities of the same species discussed before, comparing our reference run with OTF chemistry (black solid line) with PP1 (green dashed line), PP4 (cyan dot-dashed line), and PP5 (magenta dotted line). The right-hand panels, where we report the number density profiles, clearly show that the time-step choice heavily affects the high-density gas abundances of two of the three species shown (CO and N$_2$H$^+$), while at large radii (lower densities) the effect becomes negligible. In particular, when longer integration time-steps at constant density are used, the solution obtained for high-density gas, that is expected to evolve more rapidly, tends to approach the equilibrium solution, which does not necessarily match the solution in the reference simulation, leading to the observed overestimation. Interestingly, this effect is smeared out in the column densities  shown in the left-hand panels, because of the averaging effect resulting from the integration along the line-of-sight (LOS), and makes the time-step choice almost irrelevant. This result highlights how crucial a careful choice of the integration time-step is to improve the agreement with the reference simulation, and that it should be calibrated on the local quantities rather than on integrated ones.

\section{Error analysis and computational efficiency of the post-processing}
\label{sec:comparison}
\begin{figure*}
     \centering
     \includegraphics[width=\columnwidth]{Fig_A1_183kyr.png}
     \includegraphics[width=\columnwidth]{Fig_A2_183kyr.png}
     \caption{Fractional abundance of N$_2$D$^+$ as a function of density for our reference and PP runs at $t=183$~kyr. In the left panels, we show the impact of the $f_{\rm part}$ parameter, whereas in the right ones we show the effect of a longer integration time-step. The average abundance of N$_2$D$^+$ for the different runs is reported in the top panels, and the relative error between each PP run and our reference run in the bottom ones, following the same colour scheme of the previous figures. Finally, as a reference, in the top panels we also show the actual chemical abundance distribution from the reference run as a grey shaded area.}
     \label{fig:species}
 \end{figure*}
To better assess the observed differences in the validation processes, we now look more in detail at the error associated with the species abundance distribution. In Fig.~\ref{fig:species}, we focus in particular on the mass fraction of N$_2$D$^+$ ($f_{\rm N_2D^+}$) as a function of density, comparing the reference simulation with the post-processing at $t=183$~kyr.
It is important to mention that for all the subsequent analysis, we focus on  particles above a density threshold of $10^3$ cm$^{-3}$ only, in the aim at consistently tracing particles belonging to the core, removing any possible influence of the low density background particles on the results.
In the top panels, we report the average mass fractions directly obtained from the reference run (in black) and those obtained by the different post-processing cases, as described in Sections~\ref{sec:subset}~and~\ref{sec:timesteps}, in bins 0.4~dex wide. In the background, we also show the full particle distribution from the reference run as a grey shaded area. PP1 is shown in red, PP2 in blue, PP3 in green, PP4 in cyan, and PP5 in magenta. From the figure, it is clear that the agreement is extremely good in all cases, with deviations happening mostly in the density intervals where the abundance evolution is fast, hence the scatter large (i.e. $n\sim10^4-10^7\rm\, cm^{-3}$), where the post-processing is not always able to accurately follow the chemical evolution.
In the bottom panels, we show the average absolute value of the relative error, defined as
\begin{equation}
    \varepsilon\equiv \frac{f_{\rm N_2D^+,PP}-f_{\rm N_2D^+,OTF}}{f_{\rm N_2D^+,OTF}},
\end{equation}
where $f_{\rm N_2D^+,PP}$ is the mass fraction from the considered PP run and $f_{\rm N_2D^+,OTF}$ the one from the reference simulation with chemistry solved OTF.
The line styles and colours are the same as in Figures~\ref{fig:profiles_fpart}~and~\ref{fig:profiles_times}, and correspond to results obtained with the post-processing for the evolved particles only (not the interpolated ones).
The impact of $f_{\rm part}$, reported in the left-hand panel, is almost negligible across the entire density range sampled, and the abundance errors never exceed the ten percent level. On the other hand, the integration time-interval has a stronger impact, as can be inferred from the larger variations at the highest densities, where the errors can reach higher values. In general, the discrepancies between the post-processing and the OTF calculations can be easily explained with the density evolution time-scale. In fact, when the density changes slowly, the constant density assumption over the time-step between simulation outputs is reasonable, with the abundances approaching the chemical equilibrium solution. When the density changes quickly, instead, the gas is not able to approach the equilibrium conditions, and the post-processing, that does not keep up with this evolution, gives inaccurate results.

To summarise, we report in Table~\ref{tab:errors} an estimate of the typical error of the post-processing technique, defined by the 10-th ($\varepsilon_{10}$) and 90-th ($\varepsilon_{90}$) percentiles of the relative error distribution for each post-processing case, resulting either from the post-processed particles only (fourth column), or from the whole particle set (last column). Also in this case, we only consider N$_2$D$^+$, for consistency with Fig.~\ref{fig:species}. The extremely small variation among PP1, PP2, and PP3 confirms that the post-processing is only moderately affected by the sampled fraction of particles chosen, with the subsequent interpolation scheme to recover the whole particle distribution being the main source of error, especially when a small fraction of particles is used, raising it to about a factor of $\sim2$ at most. On the other hand, for the runs with longer time-steps between chemistry calculations (PP4 and PP5), the typical error gets progressively worse both on the evolved particle subset and the whole distribution, as expected from the key dependence of the calculation of chemical abundances on the time evolution of the system.
\begin{table}
	\centering
	\caption{Typical error of the PP runs for N$_2$D$^+$ for the directly post-processed particles only ($\varepsilon_x^{\rm e}$) and the entire particle distribution (evolved and interpolated; $\varepsilon_x^{\rm t}$), with $x=10$ or 90 the error distribution percentiles (see text). Consistently with the previous analysis, the impact of $f_{\rm part}$ is very moderate, whereas the time interval choice has a much stronger impact.}
	\label{tab:errors}
	\begin{tabular}{lccccc}
		\hline
		Model & $f_{\rm part}$ & $\Delta t$ & $[-\varepsilon_{10}^{\rm e};\varepsilon_{90}^{\rm e}]$ & $[-\varepsilon_{10}^{\rm t};\varepsilon_{90}^{\rm t}]$\\
		\hline
		PP1 & 0.1 & $\Delta t_0$ & $[<0.001;0.013]$ & [0.259;0.384] \\
		\\
		PP2 & 0.01 & $\Delta t_0$ & $[<0.001;0.014]$ & [0.476;1.139] \\
		PP3 & 0.001 & $\Delta t_0$ & $[<0.001;0.013]$ & [0.667;1.888] \\
        \\
        PP4 & 0.1 & $5\Delta t_0$ & [0.012;0.059] & [0.254;0.489] \\
		PP5 & 0.1 & $8\Delta t_0$ & [0.023;0.100] & [0.248;0.569] \\
		\hline
	\end{tabular}
\end{table}

In conclusion, our analysis suggests that, in order to get a good accuracy in the chemical evolution, a high output frequency of the underlying hydrodynamic simulation (without chemistry) is needed to avoid large errors in the recovered abundances. In particular, our choice for $\Delta t_0$ was dictated by the simultaneous need of a high output frequency (able to accurately track the typical dynamical evolution time-scales in the simulation) and a low enough number of snapshots (given the large storage overhead due to the chemical species abundances and the long timescales we wanted the evolution to last). Although our choice has been shown to be adequate, an even higher output rate should be preferred when chemistry is not solved OTF, as is the case for the runs aimed at be post-processed. Note that, however,  this would make the post-processing approach more computationally expensive, hence less convenient against a fully coupled simulation. When chemistry is solved on-the-fly, instead, such a high output frequency is not needed, since chemistry is consistently evolved on time-scales that are typically much shorter than those between outputs.

Finally, a last aspect of our post-processing method needs to be addressed, i.e. its computational efficiency, that has been improved in our analysis via OpenMP parallelisation. In Fig.~\ref{fig:postproc_cputime}, we report the computational cost of the different cases as the total time of the run in CPU hours, to highlight how much time is gained via the post-processing technique relative to OTF calculations. Here, our reference MHD simulation with OTF chemistry is shown as a black/cyan bar, whereas the post-processing models are split in a black part, the same for all models, corresponding to the cost of a MHD run without OTF chemistry, and an cyan part resulting from the individual cost of the post-processing step, i.e. evolution+interpolation. As expected, the cost increases in a sub-linear fashion with the number of particles in the post-processed sub-sample, because of the competing effect of the actual post-processing and the interpolation step.
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{postproc_bartime.pdf}
    \caption{CPU time in hours used for the reference simulation with OTF chemistry (Sim), and the different post-processing cases, PP1 ($f_{\rm part}=0.1$), PP2 ($f_{\rm part}=0.01$), PP3 ($f_{\rm part}=0.001$), PP4 ($\Delta t = 5 \Delta t_0$), and PP5 ($\Delta t = 8 \Delta t_0$).}
    \label{fig:postproc_cputime}
\end{figure}
Obviously, if we extrapolate the CPU time to all the particles of the simulation, we expect to obtain a cost comparable to that of the fully coupled simulation, hence the post-processing stops being convenient. If we instead increase the time interval between density updates, the cost mildly decreases, resulting in no significant gain. This is due to the fact that, although the number of density updates is reduced, the chemical calculations must be performed over a longer time-scale per step, resulting in only moderately different computational times. Moreover, as discussed above, longer time-steps can significantly overestimate the chemical abundances, especially at the highest densities, making this choice generally less convenient.

Concluding, our analysis suggests that the post-processing technique applied to small-scale simulations yields good results, with relative errors on average at the percent level being achieved by using a subset of just 10 percent of the total particles, and assuming short enough, but still reasonable, time-steps between density updates. The method appears reliable also when the number of particles used is significantly smaller than that of the simulation ($f_{\rm part}\lesssim 0.1$), and its accuracy is not affected by this choice, as long as the the time interval between density updates is not increased too much. Nevertheless, in our experiments, the actual gain of the post-processing compared to coupled OTF chemistry calculations is not more than a factor of a few, hence it appears convenient only when the same simulation is expected to be used multiple times to explore the relative importance and the effect of different chemical processes, rates, or when including huge networks, which are challenging to solve on the fly. As a final note, based on our findings, and on the average free-fall time $t_{\rm ff}$ of the considered core, we can speculate that an output frequency of $\Delta t\ll 10^{-2}t_{\rm ff}$ is sufficient to maintain the chemical abundances of the majority (80\%) of the post-processed particles below an average error of 1\%.

\section{Conclusions}
\label{sec:conclusions}
In this study, we assess the reliability and the uncertainties of the chemical post-processing technique applied on 3D MHD simulations, with the aim at facilitating the implementation of complex chemical networks in the context of star and planet formation processes. For our analysis, we employed one of the simulations of collapsing massive cloud fragments introduced by \citet{Bovino2019}. After having extracted the dynamical evolution of a subset of particles, we evolved the chemistry via the non-equilibrium chemical code \textsc{krome}, and interpolated the abundance for the remaining particles in each simulation output to recover an approximate version of the chemical abundances for the full particle distribution. To validate such an approach, we compared the results obtained with different combinations of key input parameters, namely, the fraction of particles extracted for the post-processing and the coarse time-step used for the chemical calculations, against the simulations performed by \citet{Bovino2019}, in which the chemistry calculations were performed on-the-fly fully coupled within the MHD code \textsc{gizmo}. This is the first time that the effective errors related to post-processing procedures is comprehensively explored, made it possible from our powerful 3D MHD simulations with on-the-fly non-equilibrium chemistry.

Overall, we found that the post-processing technique is quantitatively reliable, provided the particle sub-sample to be directly post-processed is, at least, above 1\% of the total simulation ensemble, in order to limit the impact of the interpolation scheme on the abundance error, and significantly cheaper than simulations with fully coupled OTF chemistry. However, if the aim is to recover specific features as those obtained in the MHD simulations with OTF chemistry, a larger amount of post-processed particles ($f_{\rm part} \geq 0.1$) must be used, resulting in a less significant gain of computational time. Although our results suggest that the exact choice of the parameter $f_{\rm part}$ only moderately affects the results, we stress that, when the particle sampling is not good enough (below 1\%), the chemical evolution will more coarsely sample the fiducial solution, independent of the chosen time-step, and the results will not always trace the chemical structure of the system under study.

A few caveats in the post-processing approach must be mentioned:
\begin{itemize}
\item First, a crucial issue is the interpolation of the non-post-processed particles, which is based on the density-binned abundances obtained from the post-processed particles. While the coupling of the density-binned interpolation with a distance-based neighbour selection to assign the species abundances allows to take into account the evolutionary history of the gas around the target particle, it is still not sufficient to perfectly recover the abundances, thus introducing some artefacts in the species abundances. However, we have to keep in mind that the intrinsic scatter can only be recovered with OTF chemistry, because of its strong dependence on particle-by-particle variations in the dynamic properties, typically occurring on timescales that cannot be easily tracked with the post-processing unless an extremely short output time-step is assumed (removing any gain of the post-processing). In addition, compared to a simpler density-based-only interpolation, in which the spatial localisation of the gas particles is completely neglected, and the scatter in the chemical abundances completely washed out, our density+proximity approach is able to improve the accuracy of the post-processing by a factor of 2--3.
\item An additional issue, which is connected to the one just discussed, is the poor sampling of some density intervals, that would produce large errors also in the average chemical properties. This limitation can be overcome by appropriately tuning the $\xi_{\rm core}$ parameter, and increasing $f_{\rm part}$, but at the cost of a more demanding analysis.
\item It is also important to notice that a low number of particles does not allow to recover possibly important spatial features in the system, and should be avoided, in particular considering the moderate increase in computational cost when a larger fraction of particles is employed.
\item Finally, a fundamental limitation of the post-processing approach is in general the time-step over which the density is kept fixed for the chemistry calculations. While a higher output frequency can ideally improve the accuracy of the results, it also makes the post-processing more expensive, hence closer to a proper OTF chemical evolution.
\end{itemize}

Apart from these caveats, the post-processing of simulations with the method described in the present work can be extremely useful when exploring the relative effect of key input parameters, saving a considerable amount of time and (super-)computing resources, or when it is necessary to employ large networks that would be otherwise impossible to include on-the-fly. However, when a very accurate model as close to reality as possible is required (within the limits of the underlying chemical model), on-the-fly chemical calculations within (magneto-)hydrodynamic simulations should always be preferred, especially when there is a direct feedback between the evolution of chemistry and dynamics, e.g. when the dynamical evolution of the gas is also determined by the temperature variations, that in turn are properly traced by the balance between chemical reactions.

\section*{Acknowledgements}
We thank the anonymous referee for the useful suggestions that improved the quality of the manuscript.
SFC thanks funding from CONICYT Programa de Astronom\'ia Fondo ALMA-CONICYT 2017 Project \#31170002.
AL acknowledges support from MIUR under the grant PRIN 2017-MB8AEZ.
SB acknowledges financial support from Millenium Nucleus NCN19\_058 (TITANs), and BASAL Centro de Astrofisica y Tecnologias Afines (CATA) AFB-17002. The simulations were performed with resources provided by the KULTRUN Astronomy Hybrid Cluster at Universidad de Concepción.

\section*{Data Availability Statement}
The data underlying this article will be shared on reasonable request to the corresponding author.

\bibliographystyle{mnras}


\section{Introduction}

Galaxy surveys have provided observations of  the growth of large-scale structure  with steadily increasing precision and range. This trend is set  to accelerate in the coming years with next-generation spectroscopic surveys such as those with {\em Euclid} \cite{Blanchard:2019oqi}, the Square Kilometre Array (SKA) \cite{Bacon:2018dui}, the Dark Energy Spectroscopic Instrument (DESI) \cite{Levi:2019ggs}, and the Nancy Grace Roman Space Telescope \cite{Dore:2019pld}. The evolution of matter perturbations gives direct insight into the underlying theory of gravity. Whether the correct description is  the standard model of cosmology $\Lambda$CDM, a dynamical dark energy model or modified gravity  (see e.g. \cite{Capozziello:2011et,Clifton:2011jh} for reviews), is still an open question.

Studying gravity on cosmological scales has been pursued via two paths in recent years. The model-dependent approach constrains a particular model of dark energy or modified gravity with observations and compares it to $\Lambda$CDM. This approach has shown hints that alternative theories may be favoured via Bayesian evidence considerations \cite{Peirone:2019aua,Sola:2019jek,Aoki:2020oqc,DiValentino:2020kpf}. The `agnostic' approach avoids choosing a particular model and develops predictions as model independent as possible, which can be used to evaluate how they differ from those of a given model.

Machine learning offers several model-agnostic methods, including Gaussian Process (GP) regressions, which enable data-driven reconstructions of trends \cite{Rasmussen}. The characteristic input of GP is statistical and corresponds to a kernel to build a covariance matrix and a mean prior. The use of GP is now common in cosmology and has been focused mainly on the background evolution of the universe \cite{Holsclaw:2010nb,Holsclaw:2010sk,Holsclaw:2011wi,Shafieloo:2012ht,Seikel:2012uu,Yahya:2013xma,Bester:2013fya,Busti:2014dua,Bester:2015gla,Bester:2016fbs,Joudaki:2017zhq,Gomez-Valent:2018hwc,Haridasu:2018gqm,Gerardi:2019obr,Keeley:2019hmw,Bengaly:2019oxx,Bengaly:2019ibu,Aljaf:2020eqh,Keeley:2019hmw,Colgain:2021ngq,Dhawan:2021mel,Mukherjee:2021kcu,2021arXiv210407356V,Canas-Herrera:2021qxs}. GP has also proved useful to determine the value of the Hubble constant using strong lensing data and supernova data \cite{Liao:2019qoc,Liu:2019ddm,Pandey:2019yic,Liao:2020zko,Renzi:2020fnx,Bengaly:2020neu,Zhang:2021tmg}, to improve photometric redshift estimates \cite{2016MNRAS.462..726A}, and to speed up the processing of weak lensing data \cite{Mootoovaloo:2020ott}. Less work has been done with GP on measurements of perturbations. Examples are reconstructions of the linear anisotropic stress  \cite{Pinho:2018unz} and of the growth rate of large-scale structure \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt}.

In this paper we use GP to reconstruct the growth rate of large-scale structure. We produce forecasts for a nominal Stage IV spectroscopic survey. Our focus is not only on the growth variable $f\sigma_8$, but also on the separated growth rate $f$ and variance of matter fluctuations $\sigma_8$. We apply the concept of multi-task GP \cite{Caruana1998,Rasmussen,866abfcac5f9481d97628a546255878b,21ae58e5ad934d7494d64e5bfd4a1a52,Melkumyan,Vasudevan}, previously used in cosmology to improve reconstructions of the background acceleration \cite{Haridasu:2018gqm}.

The paper is structured as follows. We review the basics of the growth of large-scale structure and the methodology of our analysis in \autoref{sec:context}. We then analyse and compare the single- and multi-task reconstructions of mock data in \autoref{sec:perfo}, where we display the improvement from using multi-task GP. Previous single-task GP studies of growth rate in \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt} do not give much detail about the assumptions specific to their reconstruction of $f\sigma_8$. In \autoref{sec:current}, we discuss difficulties in reproducing such results and how we find current growth data to be ill-suited for robust GP reconstructions. We conclude in \autoref{sec:discussion}.

\section{Context}\label{sec:context}

We briefly review the growth of large-scale structure before presenting the GP reconstructions from the single-task to multi-task approach. Then we show how our mock data are derived and the numerical tools that we employ.

\subsection{Growth of large-scale structure}\label{sec:growthgen}

The galaxy distribution in redshift space is distorted  by galaxy peculiar velocities that are induced by the growth of cosmological structures, an effect known as redshift-space distortions (RSD) \cite{1977ApJ...212L...3S,Kaiser:1987qv,Guzzo:2008ac,Song:2008qt,Percival:2008sh}. RSD measurements determine the product of  the growth rate $f$ and variance of matter fluctuations $\sigma_8^2$, where
\begin{eqnarray}
\label{eq:defD}
    f(z)&=&-\frac{{\rm d}\ln{D(z)}}{{\rm d}\ln(1+z)}\, ,\\
\sigma^2_8(z) &=&\frac{1}{2\pi^2}\int_0^{\infty} P_m(k,z) \, \hat{W}^2(k)\, k^2 \, {\rm d} k\, .
\end{eqnarray}
Here $D(z)=\delta_m (\bm k,z)/\delta_m(\bm k,0)$ is the linear growth factor of the matter density contrast $\delta_m$, $P_m$ is the matter power spectrum and $\hat{W}$ is the Fourier transform of the top-hat window function that  smooths over spheres of radius $8\,{\rm Mpc}/h$.

Both these quantities depend on the evolution of $\delta_m$, whose redshift dependence and sometimes scale dependence  will differ between models of gravity. The reconstruction of $f$ and $\sigma_8$ is then an ideal tool to distinguish between  different theories. In general relativity with standard dark energy, the linear growth rate does not depend on the scale $k$, but this does not necessarily hold in alternative theories of gravity.

RSD are not sensitive to the two functions $f$ and $\sigma_8$ separately, but to their product $f\sigma_8$. Disentangling $f$ and $\sigma_8$ given a measurement of $f\sigma_8$, requires their degeneracy to be broken.  It is possible to break this degeneracy through two different techniques, both requiring additional observations: (i)~by combining RSD measurements in the power spectrum and bispectrum \cite{Gil-Marin:2016wya}; (ii)~by combining RSD data with galaxy-galaxy lensing data \cite{delaTorre:2016rxm,Shi:2017qpr,Jullo:2019lgq}.

The use of GP to reconstruct trends from current $f\sigma_8$ data has been studied previously in \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt}. Given the low precision and relative randomness of current data, we find that the reconstructions depend too strongly on the GP assumptions made. We expand on this  in \autoref{sec:current}.

\subsection{Basics of Gaussian processes}\label{sec:stgp}

Reconstructing a latent function from data with GP essentially corresponds to assuming that the data is a random realisation of the function with Gaussian noise given by the covariance matrix of the data. The reconstruction of $\bm \mathfrak{f}^*$ as a function of $\bm X^*$, given data {$\bm y$} as a function of  $\bm X$, with a covariance matrix $C$, is characterised by the joint probability distribution \cite{Rasmussen}
\begin{equation}\label{eq:jointdistr}
\begin{bmatrix} \bm y \\ \bm \mathfrak{f}^* \end{bmatrix}\! \sim \!
\mathcal{N}\! \left(\!
\begin{bmatrix}
\bm \mu  \\
\bm \mu^*
\end{bmatrix}
\!,
\begin{bmatrix}
K(\bm X,\bm X) +C & K(\bm X,\bm X^*)\\
K(\bm X^*,\bm X)    & K(\bm X^*,\bm X^*)
\end{bmatrix}\!
\right)\!.
\end{equation}
Here $\bm \mu$ and $\bm \mu^*$ correspond to  priors on the means of the reconstructed functions  of $\bm X^*$ and $\bm X$ respectively. $K$ is the kernel of the GP, i.e. the covariance function that relates the values of the function to be reconstructed over the points $x$ in $\bm X$ and $x^*$ in $\bm X^*$. From the joint distribution of \autoref{eq:jointdistr}, the mean and covariance of the reconstructed function are
\begin{align}\label{eq:gpmean}
{\rm mean}(\mathfrak{f}^*)  &=  \bm \mu^* \\\nonumber
& + K(\bm X^*, \bm X)\!\left[\!K(\bm X,\bm X) + C\right]^{-1}\!\! (\bm y -\bm \mu),  \\[2mm]
\label{eq:gpcov}
{\rm cov}(\mathfrak{f}^*)  & = K(\bm X^*,\bm X^*) \\\nonumber
& - K(\bm X^*,\bm X)\! \left[\!K(\bm X,\bm X) + C\right]^{-1} K(\bm X,\bm X^*) .
\end{align}
The marginal log likelihood of the reconstructed function can then be obtained as
\begin{align}\label{eq:gp_log_marginal}
\ln \mathcal{L} = & -\dfrac{1}{2}(\bm y - \bm \mu)^T \left[\! K(\bm X,\bm X) + C \right]^{-1} (\bm y - \bm \mu) \nonumber\\
& - \dfrac{1}{2}\ln\left| K(\bm X,\bm X) + C \right| - \frac{n}{2}\ln 2\pi  \;,
\end{align}
where $n$ is the number of data points.

It should  now be clear that a GP allows us, given a set of training data points, to reconstruct a cosmological function without the need to assume its trend or a parametrisation, but only relying on the assumption that the data are a Gaussian realisation of the underlying function. For this reason, GP reconstructions fall in the category of non-parametric methods.  However,  this does not imply that GP are independent of  free quantities to be fixed. A kernel must be chosen, which depends on a set of `hyperparameters' that need to be evaluated during the reconstruction process. In addition, the reconstruction also requires a mean prior function.

The kernel most often used in cosmology so far is the squared-exponential
\begin{align}\label{eq:SE}
k(x,x') & = \sigma^2 \exp\Big(-\dfrac{\vert x-x' \vert ^{2}}{2 \xi^{2}}\Big).
\end{align}
Two hyperparameters are related to the characteristic scales in the data:  the correlation length  $\xi$ and signal variance $\sigma$ \cite{Seikel:2013fda}. (See Chapter 4 of \cite{Rasmussen} for more details on covariance functions.)

\subsection{Multi-tasking the growth of structures}\label{sec:mtgp}

The previous section describes how to reconstruct one function from a data set, i.e. a single-task GP. In order to generalise this for several data sets and latent functions which describe the same underlying phenomena, we use a multi-task GP \cite{Caruana1998,Rasmussen}. This extension  properly accommodates possible correlations between the data sets for each function. As described in \autoref{sec:growthgen}, it is straightforward to obtain $f\sigma_8$ measurements, while a more involved procedure is then needed to obtain the disentangled $f$ and $\sigma_8$ contributions. If these data sets originate from different surveys and different methods, the data would in principle not be correlated between components. In general there will be correlations between the data sets. The multi-task GP is able to capture such correlations.

In addition to a joint likelihood that incorporates correlations between the three data sets, multi-task GP also allows us to takes into account that the growth functions are inter-dependent, since all three are probes of the evolution of matter perturbations. Furthermore, one function is a product of the other two.

In the general case, a combined data set including $f$, $\sigma_8$ and $f\sigma_8$ data has a covariance matrix in block form,
\begin{equation}
C=\!
  \left[\! \begin{array}{lll}
   {{\rm cov}(f,f)}     & {\rm cov}(f,\sigma_8)     & {\rm cov}(f,f\sigma_8)    \\
   {\rm cov}(f,\sigma_8)  & {{\rm cov}(\sigma_8,\sigma_8)}  & {\rm cov}(\sigma_8,f\sigma_8) \\
   {\rm cov}(f,f\sigma_8) & {\rm cov}(\sigma_8,f\sigma_8) & {{\rm cov}(f\sigma_8,f\sigma_8)}
  \end{array} \! \right],
\end{equation}
where diagonal entries are  auto-covariances  and off-diagonal are cross-covariances.  Similarly, the multi-task GP covariance function encodes not only the three auto-covariance functions, but also the cross-covariances arising from their  inter-dependence as theoretical probes of growth. It has the form
 \begin{equation}\label{eq:Ktilde}
\tilde{K} =
  \left[ {\begin{array}{lll}
   K_{f,f}      & K_{f,\sigma_8}    & K_{f,f\sigma_8}     \\
   K_{f,\sigma_8}  & K_{\sigma_8,\sigma_8}  & K_{\sigma_8,f\sigma_8}  \\
   K_{f,f\sigma_8} & K_{\sigma_8,f\sigma_8} & K_{f\sigma_8,f\sigma_8}       \\
  \end{array} } \right] + C .
\end{equation}
Provided the data covariance matrix $C$ is block diagonal, setting all non-diagonal blocks of $\tilde{K}$ to zero leads back to the standard approach where each latent function for each data set would be reconstructed as a `single-task', i.e. independently of the other sets. In practice, this is indeed equivalent to considering separately each single-task $i$, where the relevant covariance matrix for the GP is $K_i+C_i$, and doing the respective reconstructions. By contrast, writing the full $\tilde{K}$ as in \autoref{eq:Ktilde} allows multiple functions to be reconstructed simultaneously, including the covariance between the different tasks \cite{Haridasu:2018gqm,866abfcac5f9481d97628a546255878b,21ae58e5ad934d7494d64e5bfd4a1a52,Melkumyan,Vasudevan}.

The $\tilde{K}$ kernel of the multi-task GP is constructed as follows. First, the sub-kernels for each diagonal  latent function must be chosen. Then, the corresponding cross-kernels are obtained from the convolution of the basis functions of each  sub-kernel (see for instance \cite{Melkumyan} and Section 2.2 of \cite{Haridasu:2018gqm} for details). The multi-task GP construction implies that the final reconstruction is characterised by a set of hyperparameters brought by the kernels describing the auto-covariance of each latent function, as in the single-task approach, while the cross-terms do not contribute any additional parameter to the reconstruction. For the case of the squared exponential kernel \autoref{eq:SE}, the basis function is defined as:
\begin{equation}
    g(x) = \sigma\left(\frac{2}{\pi \xi^2}\right)^{1/4}\exp{\left(-\frac{x}{\xi}\right)^2}\;,
\end{equation}
and the convolution of two squared exponential kernels is
\begin{equation}
    \label{eqn:ksese}
    k_{1 \times 2}(x,x') = \sigma_1 \sigma_2 \left(\frac{2\, \xi_1 \xi_2 }{\xi_1^2+\xi_2^2}\right)^{\!1/2}\!\! \exp{\left(-\frac{\vert x-x' \vert^2}{\xi_1^2+\xi_2^2}\right)}\!.
\end{equation}

An important issue for multi-task GP is {`negative transfer'} and {`over-fitting'}. As discussed in \cite{Haridasu:2018gqm}, these occur when one data set is much more constraining than the others and so dominates the reconstruction. It is particularly problematic in regions where data are sparse. This is something which will be of importance with future data. For instance, higher quality and quantity data on $f\sigma_8$ is likely to be released sooner than new data on $f$ or $\sigma_8$. We leave the study of these caveats for future work.

\subsection{Mock data and analysis method}\label{sec:models}

The reference model we consider is flat $\Lambda$CDM, using the parameter values from the \emph{Planck} collaboration 2018 release \cite{Aghanim:2018eyx}. We set parameters to the mean values obtained via the robust combination of constraints from\footnote{See \emph{Planck} Legacy Archive: \url{https://pla.esac.esa.int/\#cosmology}.}: cosmic microwave background (CMB) base temperature, polarisation and lensing; baryon acoustic oscillations (BAO);  supernovae type Ia (SNIa). We then use \texttt{CLASS} \cite{2011JCAP...07..034B} to compute the redshift evolution of $f$, $\sigma_8$ and $f\sigma_8$ fiducials, in order to build our mock data sets. We generate simplified mocks as follows:
\begin{itemize}
\item
Each mock contains 20 measurements uniformly distributed between redshift 0 and 2.
\item
We assume a constant relative error of 1\%, consistent with Stage IV expectations.
\item
We define a multivariate Gaussian distribution with mean equal to the fiducial values of each growth function, and diagonal covariance  matrix defined by the fiducial values times the relative error. We do not incorporate cross-correlations in the mock data since these are too survey-specific. Mock data are then generated as random realisations of this distribution.
\end{itemize}
\vskip1mm
We optimise the hyperparameters, i.e. find their best fit, by minimising the log marginal likelihood in \autoref{eq:gp_log_marginal} applied to the data set. We use the stochastic minimiser `differential evolution' from the python scipy package\footnote{\url{https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html}} in our GP code\footnote{ \url{https://gitlab.com/perenonlouis/GPreconstructions}}. The mean and  covariance of the final reconstruction are then obtained using this set of parameters in \autoref{eq:gpmean} and \autoref{eq:gpcov}. For the sampling of the posteriors of the hyperparameters via Markov chain Monte Carlo (MCMC) methods, we interface our GP code with \texttt{Cobaya\footnote{\url{https://cobaya.readthedocs.io/en/latest}} }\cite{Torrado:2020dgo}. This allows us to use a Metropolis Hasting sampler \cite{Lewis:2002ah,Lewis:2013hha} and the nested sampler \texttt{Polychord} \cite{Handley:2015fda,2015MNRAS.453.4384H}. We also interfaced the Affine-Invariant Ensembler Sampler \texttt{emcee\footnote{\url{https://emcee.readthedocs.io/en/stable/}} } \cite{ForemanMackey:2012ig}. All our statistical analyses and post processing of our results are done with \texttt{GetDist\footnote{\url{https://getdist.readthedocs.io/en/latest/}}} \cite{Lewis:2019xzd}.

\section{Performance and robustness}\label{sec:perfo}

Here we investigate the performance of multi-task GP compared to single-task and make some tests of the results.

\subsection{`Sharing is caring'}\label{sec:sharing}

\begin{figure}[!]
\begin{center}
\includegraphics[trim=0 77 0 0,clip,width=0.92\linewidth]{Figure_1_zplot_fsigma8}\\
\hskip3mm\includegraphics[trim=0 77 0 0,clip,width=0.88\linewidth]{Figure_1_zplot_f}
\includegraphics[width=0.92\linewidth]{Figure_1_zplot_sigma8}
\caption{Optimised single- (grey) and multi-task (red) reconstructions of $f\sigma_8$, $f$ and $\sigma_8$ with 1\% relative error mocks, using $\Lambda$CDM fiducial constrained by CMB T+P+L, BAO, SNIa (see \autoref{sec:models}), for squared-exponential kernel and null mean prior. Shading indicates 68\% confidence intervals about the mean (solid curve). The multi-task GP reconstructs the three growth functions simultaneously by minimising a common likelihood (see \autoref{sec:mtgp}) -- thus `sharing' information between the three reconstructions. As a result the single- and multi-task predictions do not match exactly: multi-task reconstruction for a given function benefits from the constraints on the other functions, thereby decreasing the error on all the reconstructions. Quantitatively, the multi-task GP provides the best fit to the mock data.}
\label{fig:stgp_vs_mtgp_plotz}
\end{center}
\end{figure}

We use a squared-exponential kernel \autoref{eq:SE} and a null mean prior as the base assumptions for the single- and multi-task GP. Note that the reconstructions are robust under the change of these assumptions for reasonable choices. This is studied in \autoref{sec:assumptions}.

An overlay of the single- and multi-task reconstructions is shown in \autoref{fig:stgp_vs_mtgp_plotz}, for the generation of a single mock. At this level of precision, the reconstructions are very sensitive to the distribution of the data around the fiducial. We discuss  the deviations of the reconstructions in \autoref{sec:deviations}. We concentrate for now on the differences between single- and multi-task reconstructions.

\autoref{fig:stgp_vs_mtgp_plotz} further shows that multi-task GP increases the precision of reconstructions compared to single-task GP. It also changes the redshift trends relative to the single-task case. Although the error for the multi-task GP is in general reduced, this is not guaranteed at all redshifts. The changes in the redshift trends and error sizes are a clear manifestation of the sharing of information between the reconstructions of functions in the multi-task case. Quantitatively, we find that the multi-task log marginal likelihood is larger than the sum of the single-task ones.

\begin{figure*}[!]
\begin{center}
\includegraphics[scale=0.58]{Figure_2_contours_stgp_fsigma8}
\includegraphics[scale=0.58]{Figure_2_contours_stgp_f}
\includegraphics[scale=0.58]{Figure_2_contours_stgp_sigma8}\vskip4mm
\includegraphics[scale=0.45]{Figure_2_contours_mtgp}
\caption{Posterior distributions from the MCMC sampling of the hyperparameters for the single-task (grey contours) and multi-task (red contours) reconstructions shown in \autoref{fig:stgp_vs_mtgp_plotz}. Dashed lines indicate optimised values. The transfer of information between reconstructions in the multi-task case adds correlations between hyperparameters of each kernel and the hyperparameter posteriors become much wider. It also increases values of all hyperparameters relative to single-task, so that the reconstructions will be more `wiggly', as seen in \autoref{fig:stgp_vs_mtgp_plotz}.}
\label{fig:stgp_vs_mtgp_contours}
\end{center}
\end{figure*}

A complementary approach to optimising the GP is to reconstruct the full posterior distributions of the hyperparameters. This allows in particular the possibility to study the degeneracies between them. Sampling the space of hyperparameters from the log marginal likelihood in \autoref{eq:gp_log_marginal} can be achieved using Monte Carlo Markov Chains (MCMC). The multi-task case is more difficult to sample because of non-global minima. We found that a nested sampler is more robust than a Metropolis-Hasting sampler. For the same reason, optimising the hyperparameters with a stochastic minimiser is safer than any method based on gradient descents or similar. The results of the MCMC sampling corresponding to the GP of \autoref{fig:stgp_vs_mtgp_plotz} are shown in \autoref{fig:stgp_vs_mtgp_contours}. Sampling the hyperparameters also allows us to derive the marginalised reconstruction. We verify that the Gaussianity of the mocks renders each best-fit value of the reconstructions to be very close to the maximum of their marginalised posterior. This implies that the optimised reconstructions are virtually the same as the marginalised ones \cite{Seikel:2012uu}.

\autoref{fig:stgp_vs_mtgp_contours} shows that the 2-dimensional posteriors between two hyperparameters of the same kernel display correlations whatever the approach. The transfer of information between the reconstructions in the multi-task case results in correlations between hyperparameters of each kernel. The hyperparameter posteriors are therefore much wider. This means that previously ruled out values of hyperparameters for one of the functions are now allowed, as long as the others change accordingly. This transfer of information also has the effect of increasing the values of all hyperparameters relative to the single-task case. Interestingly, this implies that one effect of the convolutions in the multi-task kernel is to increase the correlation of the data as seen by the GP. The hyperparameters for each reconstruction are larger and therefore they will be more `wiggly'. This is consistent with what seen in \autoref{fig:stgp_vs_mtgp_plotz}.

A Monte Carlo generation of $10^4$ mocks confirms that the above conclusions are robust. We find that the distribution of best-fit multi-task log marginal likelihood is higher than the sum of the three single-task ones, as shown in \autoref{fig:stgp_vs_mtgp_chi2s_x}. The multi-task thus performs a better fit to the union of the three growth data sets. This highlights an important point: the multi-task approach reduces the width of the confidence interval around the mean, but it also changes its redshift evolution. The single-task approach reconstructs more faithfully the redshift dependence of each data set -- simply because it is oblivious to the data on the other functions. Since the log marginal likelihood remains better for the multi-task case, this means that the reduction of the confidence interval width of the reconstruction dominates over the change of redshift evolution in the goodness-of-fit.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.92\linewidth]{Figure_3_lnlike}
\caption{Distribution of best-fit values for the multi-task log marginal likelihood (red) and the sum of single-task log marginal likelihoods (grey) from MC generation of {$10^4$} mocks. The sum of the single-task log marginal likelihoods is smaller than the multi-task log marginal likelihood, indicating that the multi-task approach fits the data better.}
\label{fig:stgp_vs_mtgp_chi2s_x}
\end{center}
\end{figure}

\subsection{Assumptions}\label{sec:assumptions}

The first step to evaluate the robustness of GP reconstructions  is to compare predictions with different kernels. The more dependent the results are on this choice, the less  model-independent is the reconstruction. For the growth of structure, we reconstruct underlying smooth functions. We therefore consider only stationary kernels. For these tests, we also concentrate on the single-task reconstruction of $f\sigma_8$.

We test all the common stationary and non-periodic kernels: exponential, squared exponential, gamma exponential, the Matern class and rational quadratic (see \cite{Rasmussen}). We find that the exponential, gamma exponential, and Matern with $\nu < 5/2$ produce jagged reconstructions, as shown in \autoref{fig:assumptions}. These are discarded  on the basis of the smoothness requirement. Reconstructions performed with the other kernels all give very similar results. The rational quadratic kernel produces virtually the same reconstruction and goodness of fit as the squared exponential. This kernel can be understood as a scale mixture, i.e. an infinite sum, of squared-exponential kernels with different correlation lengths \cite{Rasmussen}. It also has an extra  hyperparameter. The absence of difference between  reconstructions with these kernels thus highlights that adding hyperparameter freedom does not produce new features.

The Matern class  reconstructions produce errors that increase slightly as $\nu$ decreases. In the limit $\nu \rightarrow \infty$ the Matern kernel recovers the squared exponential. Not surprisingly, the case $\nu=9/2$ is very similar to the squared exponential, as illustrated in \autoref{fig:assumptions}. On the other hand, the smaller is $\nu$ the more sensitive will this kernel be to the noise of the data. This is why a reconstruction with $\nu=3/2$ shows short wiggles on top of the evolution reconstructed by the squared exponential (see \autoref{fig:assumptions}).

In conclusion, from a GP point of view the simplest kernel, the squared exponential, captures all the information about the smooth evolution underlying the data. We  do not find any need  to consider combinations of kernels.

The next step is to look at the other input defining a GP, i.e. its mean prior. The influence of the mean prior becomes important when the latter is either far away from the data or when the data displays gaps in redshift. The former case should not occur here, given reasonable user choices. The latter could be relevant in the future, depending on survey choices. To see how sensitive the input is with our mock set-up, we compare two common choices of mean prior, i.e. equal to zero or to the mean of the data, and an unreasonable choice, equal to 20. The results are shown in \autoref{fig:assumptions}.

The unreasonable choice displays large and noisy errors. Nonetheless, the precision of the data is such that the reconstruction does not go to the mean prior in between any neighbouring measurements,  even though the mean prior is chosen significantly far away from the data. For mean prior equal to the mean of the data, the maximum log marginal likelihood is slightly enhanced relative to the null mean prior. However, the reconstruction $\chi^2$ is slightly increased. For real future data, it could be more meaningful to use the mean of the data as mean prior.

Alternatively, one could input the standard model prediction as mean prior. Here, this is exactly the fiducial, and it produces a flat reconstruction around the fiducial, i.e. with correlation length hyperparameter $\xi$ diverging to infinity. Overall, we avoid a mean prior that is dependent on each mock, to ensure that each mock is treated on the same ground during the MC generations.
We conclude that the null mean prior is a sound benchmark.
\begin{figure}[!h]
\begin{center}
\includegraphics[trim=0 77 0 0,clip,width=0.9\linewidth]{Figure_7_test_kernel_fsigma8}\\
\includegraphics[trim=0 0 0 0,clip,width=0.9\linewidth]{Figure_7_test_mean_fsigma8}
\caption{Optimised single-task reconstructions of $f\sigma_8$ with 1\% relative error mock, $\Lambda$CDM fiducial constrained by CMB T+P+L, BAO, SNIa. Shading indicates 68\% confidence intervals about the mean (solid curves). Reconstructions are shown with different kernels (top) and different mean priors (bottom). The Matern kernel with $\nu =3/2$ produces jagged reconstructions, while squared exponential and Matern $\nu =9/2$ give very similar results.
The unreasonable choice of mean prior shows large and noisy errors, while sensible mean priors  produce similar reconstructions.}
\label{fig:assumptions}
\end{center}
\end{figure}

Regarding the priors on the range of the hyperparameters, we choose them large enough that they have no effect on the reconstructions. We thus ensure that the optimised value of a hyperparameter is not  its maximum or minimum allowed. Note that these considerations do not apply when using currently available growth data (see \autoref{sec:current}), highlighting how the precision of observations is a crucial ingredient to apply machine learning tools to cosmological analyses.

We also verified that marginalising over the hyperparameters produces virtually the same reconstructions as when optimised. This is expected, since the precision of the mock data and its Gaussian distribution around the fiducial produces posteriors that are very close to Gaussian, as shown in \autoref{fig:stgp_vs_mtgp_contours}.

In conclusion, this analysis shows that at this level of precision of the mock data, the reconstructions are fully driven by data, indicating that our comparisons of the single- and multi-task approaches are robust.

\subsection{Bias and deviations}\label{sec:deviations}

\begin{figure}[!]
\begin{center}
\includegraphics[trim=0 0 0 00,clip,width=0.9\linewidth]{Figure_4_average_fsigma8}
\caption{Average distribution over $10^4$ mocks of the $f\sigma_8$ reconstructions, computed by drawing 1000 realisations of each optimised GP, from the multivariate Gaussian built with their mean and covariance. We concatenate all 1000$\times${$10^4$} draws. The average fluctuates around the fiducial function, which is well recovered within the error.}
\label{fig:bias}
\end{center}
\end{figure}

Reconstructions can deviate quite significantly from their fiducials given the precision and distribution of the mock data. How likely are these departures and are any biases hidden? We run a series of tests to answer these questions.

First, we check for any overall bias in our reconstruction by performing the equivalent of an `ensemble average' of the reconstructed functions, over the $10^4$ mocks. To estimate the average distribution analytically would require the computation of a mixture of multivariate Gaussian distributions (i.e.\ the outputs from each GP optimisation) which has no closed form. Thus, we compute this average distribution numerically by drawing 1000 realisations of each optimised GP -- from the multivariate Gaussian built with their mean and covariance -- and we concatenate all 1000$\times${$10^4$} draws. We then derive the overall mean and the 68\% limits of the resulting distribution. The results for $f\sigma_8$ in are displayed \autoref{fig:bias}. We note that the average fluctuates around the fiducial function -- but nevertheless the fiducial is well recovered within the error and no significant biases can be seen. We also verified that the distributions at a given redshift are Gaussian. The same conclusions hold for the other two functions.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.68\linewidth]{Figure_5_chi2}
\includegraphics[width=0.68\linewidth]{Figure_5_chi2_50_data}
\includegraphics[width=0.68\linewidth]{Figure_5_chi2_5_percent}
\caption{Distributions of the reduced $\chi^2/n_{\rm dof}$ (see \autoref{eq:chi2}) of the combined $f\sigma_8\;, f\;, \sigma_8$ reconstructions for the $10^4$ mocks. We vary some characteristics of the mocks: 20 points uniformly distributed between redshift 0 and 2 with 1\% relative error (top); increased number of data points to 50 (middle); decreased precision of the data to 5\% for 20 points (bottom). The distributions of the reduced $\chi^2/n_{\rm dof}$ change depending on the configurations. This highlights that the sensitivity of GP's to the probability distribution of the data depends on the method -- single-task or multi-task -- and on the characteristics of the mock data itself.}
\label{fig:stgp_vs_mtgp_chi2s}
\end{center}
\end{figure}

The second check is to quantify how significantly  the reconstructions deviate from the fiducial. In order to do so, we exploit the fact that we know the true underlying model in our case study, i.e. the fiducial. We define a hybrid $\chi^2$ between the optimised reconstructions and the fiducial as
\begin{equation}\label{eq:chi2}
    \chi^2 = \sum_{g}\sum_{ij} \left(Y^i_{\rm fid}-Y^i_{g}\right) C^{-1}_{ij,{g}}\left(Y^j_{\rm fid}-Y^j_{g}\right)\, ,
\end{equation}
summed over the three growth reconstructions $g=\lbrace{f,\, \sigma_8,\, f\sigma_8 \rbrace}$. For each reconstruction, the $\chi^2$ is computed for a vector $\bm{Y}$ of 100 redshifts uniformly distributed between 0 and 2. $\bm{C}_{g}^{-1}$ is the corresponding inverse diagonal covariance matrix of the reconstruction. The number of degrees of freedom is therefore $n_{\rm dof}=3 \times 100$. The results are displayed in \autoref{fig:stgp_vs_mtgp_chi2s}.

We find that the reduced $\chi^2$ of the GP reconstructions is not distributed around unity (see top plot). This shows that given the characteristics of the mocks -- relative error, Gaussian probability distribution, number of data -- the reconstructions are too sensitive to the probability distribution of the data and do not recover the fiducial efficiently. In addition, the multi-task distribution is shifted to higher values and is flattened relative to the single-task. This is the result of the tighter confidence intervals produced by the former and not to any particular bias, since we found the reconstructions to be Gaussian distributed on average around the fiducial.

The third check is to test how modifying the characteristics of the mocks influences the previous conclusions. On the one hand, adding data points in the same redshift range (\autoref{fig:stgp_vs_mtgp_chi2s}, middle) does not affect significantly the distribution of reduced $\chi^2$ for the multi-task case. Interestingly, the reduced $\chi^2$ of the single-case is pulled towards reproducing that of multi-task. On the other hand, lowering the relative error of the data shifts the distributions toward a better fiducial recovery.

The fourth test we conduct is to analyse quantitatively how much  the reconstructions deviate from the fiducial. To this end, we compute for each mock the fraction of the redshift range over which the fiducial falls outside of the 95\% confidence interval of the reconstructions \cite{Seikel:2012uu}. We display the distribution of this fraction over the $10^4$ mocks for our three reconstructed functions in \autoref{fig:stgp_vs_mtgp_counts}. The results show a trend similar to our previous $\chi^2$ test -- i.e., the tighter multi-task errors exclude the fiducial function more often. Another important result highlighted by these plots is that the distributions depend significantly on the function considered, and hence on the characteristics of the mocks. For instance, the redshift evolution of the function has an importance since we consider relative errors to construct the mocks.

The tests in this section highlight some of the subtleties associated with GP regressions. Even when the reconstruction is not biased on average, the fact that it performs `too well' in some cases (i.e. has stringent error bars) can complicate the assessment of whether it excludes a given cosmological model or not. Moreover, as we showed above, such effects can depend on the specific choice of mocks, hence ultimately on the intrinsic properties of the data (noise amplitude, redshift dependence, etc.) to which we apply the GP reconstruction. A full diagnostic of such effects will be the topic of an in-depth future work.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.7\linewidth]{Figure_6_counts_fsigma8}
\includegraphics[width=0.7\linewidth]{Figure_6_counts_f}
\includegraphics[width=0.7\linewidth]{Figure_6_counts_sigma8}
\caption{Distribution of the fraction of the redshift range over which the fiducial falls outside of the 95\% confidence interval of the reconstructions. The distribution for each  reconstruction is computed over the $10^4$ mocks for 100 points uniformly distributed between redshift 0 and 2. The distributions depend significantly the function considered.}
\label{fig:stgp_vs_mtgp_counts}
\end{center}
\end{figure}

\section{Applications to current data}\label{sec:current}

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.9\linewidth]{Figure_8_current_zplot_fsigma8}
\includegraphics[width=0.9\linewidth]{Figure_8_current_contours_fsigma8}
\caption{{\it Top:} Optimised single-task reconstructions of current $f\sigma_8$ data with squared-exponential kernel and null mean prior. Upper prior bounds on correlation length hyperparameter $\xi_{f\sigma_8}$ are either restricted the maximum redshift range of the data (grey) or unbounded (red). Shading indicates 68\% confidence intervals. Dashed curve is the $\Lambda$CDM fiducial. {\it Bottom:} The corresponding MCMC sampling of posterior distributions on the hyperparameters. Shading shows 68\% and 95\% confidence intervals. Dashed lines indicate optimised values of hyperparameters. If $\xi_{f\sigma_8}$ is unbounded from above then its best fit wants to be infinite which produces a flat reconstruction; in practice it will be very close to its maximum range prior.}
\label{fig:reco_current}
\end{center}
\end{figure}

To conclude our analysis we investigate reconstructions of current growth data. We consider first the compilation of measurements of $f\sigma_8$ in \cite{Perenon:2019dpc} (updated using \cite{Nesseris:2017vor}). By using this data from different and non-overlapping galaxy surveys, we avoid possible double counting due to overlapping areas or different releases of the same survey. (This differs from the approach of \cite{Yin:2018mvu,Benisty:2020kdt,Li:2019nux}, where an extended set of RSD measurements is considered.)

The result (with squared-exponential kernel) is shown in \autoref{fig:reco_current}. It is apparent that current RSD data does not allow for a robust GP reconstruction, which depends on the allowed range for the correlation hyperparameter $\xi_{f\sigma_8}$ (red and grey shadings). This dependence persists for different choices of mean prior, including the standard model, and for different kernels.

Allowing the range of $\xi_{f\sigma_8}$ to be very large leads to a posterior that is unbounded from above when explored via MCMC, as shown in the bottom panel. The best-fit value for $\xi_{f\sigma_8}$ tends to infinity and in practice will always be close to the maximum range allowed. This is why the reconstruction is flat (red shading). The error on the reconstruction is then carried solely by the signal variance $\sigma_{f\sigma_8}$, whose 1-dimensional posterior distribution is on the contrary well defined.

This illustrates that current RSD data is simply too noisy for a GP reconstruction independent of assumptions. The data is effectively seen as random noise by the GP.

If we use the extended data set of \cite{Yin:2018mvu,Benisty:2020kdt,Li:2019nux}, the same problem arises. The only way we can reproduce a growth reconstruction similar to \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt} is to require the range of $\xi_{f\sigma_8}$ to be smaller than the redshift range of the data. This forces the best-fit value of $\xi_{f\sigma_8}$ to be that of the chosen maximum range, as shown in the grey shadings in \autoref{fig:reco_current}.

The ratio of the average size of data errors over the total variation of the function is a probe of how faithful the reconstruction is to the fine features of the data. We find that the smaller this ratio is, the more the reconstruction is faithful to data features, independent of priors  and with well-defined posteriors for the hyperparameters. For current data this ratio is much larger than for the mock Stage IV data. If we multiply the current RSD measurements by a factor, or equivalently reduce their covariance by the same factor squared, the reconstruction improves.

In conclusion, we find that the distribution and precision of current data produces GP reconstructions that are strongly dominated by the choices of priors, either the range of the hyperparameters or the mean. We tested whether a multi-task GP adding the few measurements of $f$ and $\sigma_8$ ameliorates the situation. The low number of data points and their precision does not allow us to evade the prior dominance. This scenario is also a typical example of an over-fitting configuration and should not be given much significance.

\section{Conclusion}\label{sec:discussion}

In this paper we explored Gaussian processes as a crucial tool for `agnostic' probes of gravity via the growth of large-scale structure. In particular, we focused on how GP can reconstruct the growth functions, namely the growth rate $f(z)$, the rms of matter fluctuations $\sigma_8(z)$ and their product $(f\sigma_8)(z)$.

We generated simulated measurements for the growth functions using a nominal sensitivity consistent with upcoming Stage IV surveys, and applied the GP pipeline to assess its ability to reconstruct the functions. Crucially, we reconstructed these functions not only separately, using single-task GP, but also via the multi-task GR approach. In multi-task GP, the three functions are reconstructed simultaneously, with a single likelihood. The kernel for the joint reconstruction contains sub-kernels for each function but also cross-convolution kernels that encode the links between the reconstructions.

A multi-task approach is more faithful than single-task, since a common likelihood allows us to use a data covariance matrix which describes fully the data, i.e. it can take into account the possible correlations in the three data sets. At a more fundamental level, multi-task GP also incorporates the inter-dependence of the growth functions -- all three are probes of the evolution of matter perturbations. Furthermore, one function is the product of the other two.

We showed how the sharing of information between the reconstructions of each growth function in the multi-task approach is not only necessary but also enhances their precision. By performing a Monte Carlo generation of a large number of mocks, we show statistically that, relative to the single-task, the multi-task approach:

\begin{itemize}
    \item
fits the data better;
    \item
is more precise;
\item
is more sensitive to the probability distribution of the data and thus can show more deviations from the fiducial.
\end{itemize}

We verified the robustness of our findings against the change of the assumptions defining GP reconstructions.

Robustness however applies only to the simulated data sets of a nominal future survey. Indeed, we found that when applying our pipeline to currently available data, the reconstruction is very sensitive to the assumptions on the GP settings: priors on the hyperparameter ranges and the mean choices have a significant impact on the final reconstruction. We found that using our baseline settings (a zero prior mean and a very extended allowed range for the hyperparameters), the GP reconstructs a flat function. In order to reconstruct functions that behave as expected, one would instead need to fine tune the hyperparameters. This behaviour is due to the noise of current data being too random.

Next-generation surveys will reach a sensitivity that will enable the robust reconstruction of the growth functions $f(z)$, $\sigma_8(z)$ and $(f\sigma_8)(z)$, using model-agnostic methods such as GP. High-precision measurements will allow us to avoid fine tuning the parameters of the reconstruction, using only limited assumptions. The data-driven model-agnostic reconstruction of the three growth functions is further improved by a multi-task approach that incorporates the inter-dependence of these functions and possible correlations between the three data sets. A multi-task GP avoids the theoretical bias inherent in a single-task approach, which neglects inter-dependence and correlations.

The ability of the GP to distinguish between a `true' deviation from a given model or an artefact from the noise of the data must be assessed with care. We showed that the accuracy of the reconstructions is highly dependent on the characteristics of the mock data, the function reconstructed and the GP method. Detecting a possible departure from the standard model with a GP reconstruction will require a more refined GP modelling in the future.

\section*{CRediT authorship contribution statement}

\textbf{Louis Perenon}: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing - Original Draft, Writing - Review \& Editing, Visualization, Supervision, Project administration. \textbf{Matteo Martinelli}: Conceptualization, Methodology, Software, Validation, Formal analysis, Writing - Original Draft, Writing - Review \& Editing. \textbf{St\'ephane Ili\'c}: Methodology, Validation, Writing - Review \& Editing. \textbf{Roy Maartens}: Conceptualization, Methodology, Writing - Original Draft, Writing - Review \& Editing, Supervision. \textbf{Michelle Lochner}: Methodology, Validation, Writing - Review \& Editing. \textbf{Chris Clarkson}: Conceptualization, Validation, Writing - Review \& Editing.

\section*{Declaration of competing interest}

The authors declare that they have no known competing financial interests or personal relationships that
could have appeared to influence the work reported
in this paper.

\section*{Acknowledgements}

We thank Sambatra Andrianomena, Mario Ballardini, Stefano Camera, Ed Elson, Julien Larena, Marco Raveri and Mario Santos for useful discussions. We thank the authors of \cite{Shi:2017qpr,Jullo:2019lgq} for making their data covariance matrix available to us. We acknowledge the Sciama High Performance Compute Cluster, which is supported by the ICG at the University of Portsmouth, and the Centre for High Performance Computing, South Africa (under the project Cosmology with Radio Telescopes, ASTRO-0945), for providing computational resources for this research project. LP and RM are supported by the South African Radio Astronomy Observatory and the National Research Foundation (Grant No. 75415). MM is supported by a fellowship (code LCF/BQ/PI19/11690015) from `la Caixa' Foundation (ID 100010434) and by the Spanish Agencia Estatal de Investigacion through the grant IFT Centro de Excelencia Severo Ochoa SEV-2016-0597. RM is also supported by the United Kingdom Science \& Technology Facilities Council (STFC) Consolidated Grant ST/S000550/1. CC is supported by STFC Consolidated Grant ST/P000592/1.

\bibliographystyle{elsarticle-num-names}