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

The 21-cm line of neutral hydrogen has been used to see clouds of gas in the local universe for decades \citep{Ewen:51}, but in recent years, it has been theorized that highly redshifted emission and absorption from this line could track the history of the hydrogen gas in the early universe \citep[see, e.g.][]{Furlanetto:06,Morales:10,Pritchard:12,Loeb:13}. The postulated observable that allows this is known as the 21-cm signal and is essentially a perturbation on top of the cosmic microwave background (CMB), which is positive (negative) when the line is in emission (absorption) with respect to the CMB.

The 21-cm signal can be observed in two different ways, an angular power spectrum that characterizes spatial variations as a function of redshift \citep[as measured by, e.g.~HERA;][]{Deboer:17} and a sky-averaged spectrum known as the global signal that statistically describes the spatial mean behavior of the universe over time, which is the subject of this paper.~Multiple single-antenna experiments that are currently observing or are under active development are attempting to measure the global signal, such as the Large-aperture Experiment to detect the Dark Age \citep[LEDA;][]{Price:18}, Experiment to Detect the Global Epoch of Reionization (EoR) Signature \citep[EDGES;][]{Bowman:18,Monsalve:19,Mahesh:21}, Shaped Antenna measurement of the background RAdio Spectrum \citep[SARAS;][]{Singh:17}, Radio Experiment for the Analysis of Cosmic Hydrogen \citep[REACH;][]{Eloy:19,Anstey:20,Shen:21}, the Probing Radio Intensity at high-$z$ from Marion \citep[PRIzM;][]{Philip:19}, the Cosmic Twilight Polarimeter \citep[CTP;][]{Nhan:17,Nhan:19}, and the space-based mission concept Dark Ages Polarimeter PathfindER  \citep[DAPPER;][]{Burns:17,Burns:21a,Burns:21b}.

As has been pointed out by many, the main problem that must be solved for the global 21-cm signal, which is on the order of a few hundred millikelvin, to be confidently quantified is its extraction from the large foreground emission, which is on the order of a few thousand kelvin \citep{Liu:13,Switzer:14,Vedantham:14,Anstey:20}.\footnote{For a more comprehensive review of techniques proposed to solve the foreground problem in the global signal context, see \cite{Tauscher:18}.}~This paper is the fourth and final installment in a series laying out a pipeline designed to solve this problem. In \cite{Tauscher:18} (hereafter referred to as Paper I), we put forth a method of using training sets of foreground and signal spectra and singular value decomposition to form models of the components instead of assuming analytical, a priori models like polynomials, which are unlikely to fit the true foreground in the presence of a typically chromatic beam used for global 21-cm experiments, as shown in \cite{Hibbard:20} and \cite{Tauscher:20b}. The methods of Paper I result in constraints on the signal in frequency-temperature space. In \cite{Rapetti:20} (hereafter referred to as Paper II), we extended the method to an efficient exploration of the full posterior parameter distribution using any chosen nonlinear signal model, resulting in parameter constraints and covariances. Paper II also presents an analytical marginalization technique that greatly reduces the number of parameters that must be explored in order to rigorously sample the posterior signal parameter distribution. In \cite{Tauscher:20a} (hereafter referred to as Paper III), we showed that, in order to achieve uncertainties allowing for confident detection of the global 21-cm signal, one must perform a fit that correlates and models many spectra simultaneously. This stands in contrast with the analysis methods of many current experiments, including EDGES \citep{Bowman:18}, which only fit individual time-averaged spectra.

In this paper, we will combine the methods of papers I-III as summarized above into a complete form that includes receiver uncertainties, which have been neglected until now.~We will lay out a pipeline that first obtains signal estimates in frequency-temperature space using a modeling technique similar to the one in Paper I and then uses those estimates to follow up with a fit using any chosen nonlinear signal model that employs an analytical marginalization technique very similar to the one described in Paper II. Throughout this paper, we will be fitting ten concatenated spectra that are simulated at different local sidereal times (LST) because, as found in Paper III, it is necessary to achieve uncertainties below the K level.

In Section~\ref{sec:AMLP}, we describe the general technique we use to explore posterior distributions with large numbers of parameters as efficiently as possible, which we term analytical marginalization of linear parameters (AMLP).
In Section~\ref{sec:fitting-strategy}, we lay out how we use the technique described in Section~\ref{sec:AMLP} to fit 21-cm global signal data. In Section~\ref{sec:models}, we specify the models used for the beam-weighted foreground, global 21-cm signal, and receiver gain and noise temperature. In Section~\ref{sec:results}, we show the results of our method in simulation. In Section~\ref{sec:discussion}, we discuss caveats and simplifications used in this proof-of-concept work and how we plan to handle them when fitting observed data. Finally, we conclude in Section~\ref{sec:conclusions}.

\section{Analytical marginalization of linear parameters (AMLP)} \label{sec:AMLP}

In this section, we describe our method of exploring certain classes of distributions with large numbers of parameters as efficiently as possible.~In particular, AMLP allows any distribution with one or more Gaussian conditional distributions to be explored much more efficiently.

\subsection{Form of joint posterior}

In the following, we assume a Gaussian likelihood function,
\begin{subequations}
\begin{align}
  {\mathcal{L}}({\boldsymbol{\theta}}) &= \Pr[{\boldsymbol{y}}|{\boldsymbol{\theta}}], \\
  &\propto \exp\left\{-\frac{1}{2}\left[{\boldsymbol{y}}-{\boldsymbol{m}}({\boldsymbol{\theta}})\right]^T{\boldsymbol{C}}^{-1}\left[{\boldsymbol{y}}-{\boldsymbol{m}}({\boldsymbol{\theta}})\right]\right\}, \label{eq:gaussian-likelihood}
\end{align}
\end{subequations}
where ${\boldsymbol{y}}$ is the data being fit, ${\boldsymbol{m}}({\boldsymbol{\theta}})$ is the model of the data evaluated at parameters ${\boldsymbol{\theta}}$, and ${\boldsymbol{C}}$ is the noise covariance of the data. Using Bayes' theorem we can find the posterior distribution on the parameters,
\begin{subequations}
\begin{align}
  p({\boldsymbol{\theta}}) &= \Pr[{\boldsymbol{\theta}}|{\boldsymbol{y}}] \\
  &= \pi({\boldsymbol{\theta}})\ {\mathcal{L}}({\boldsymbol{\theta}}),
  \label{eq:full-posterior}
\end{align}
\end{subequations}
where $\pi({\boldsymbol{\theta}})$ is the prior distribution that quantifies our a priori knowledge of the parameters.

\subsection{Conditions of AMLP} \label{sec:conditions-of-AMLP}

AMLP relies on splitting ${\boldsymbol{\theta}}$ into two groups, ${\boldsymbol{\theta}_{\text{NL}}}$ (nonlinear parameters) and ${\boldsymbol{\theta}_{\text{L}}}$ (linear parameters),\footnote{The use of NL short for nonlinear and L short for linear will be consistent throughout. In the following, we will sometimes write ${\boldsymbol{\theta}}$ dependence as dependence on ${\boldsymbol{\theta}_{\text{L}}}$ and ${\boldsymbol{\theta}_{\text{NL}}}$, e.g. ${\mathcal{L}}({\boldsymbol{\theta}})={\mathcal{L}}({\boldsymbol{\theta}_{\text{L}}},{\boldsymbol{\theta}_{\text{NL}}})$.} such that the conditional distribution of ${\boldsymbol{\theta}_{\text{L}}}$ at constant ${\boldsymbol{\theta}_{\text{NL}}}$, $\Pr[{\boldsymbol{\theta}_{\text{L}}}|{\boldsymbol{y}},{\boldsymbol{\theta}_{\text{NL}}}]$, is Gaussian. With a Gaussian likelihood as given in Equation~\ref{eq:gaussian-likelihood}, the following is a common set of sufficient conditions for this to be true:
\begin{enumerate}
  \item The model conditioned on ${\boldsymbol{\theta}_{\text{NL}}}$ is linear in ${\boldsymbol{\theta}_{\text{L}}}$, i.e. $\frac{\partial {\boldsymbol{m}}}{\partial{\boldsymbol{\theta}_{\text{L}}}}={\boldsymbol{A}}({\boldsymbol{\theta}_{\text{NL}}})$ where the gradient matrix ${\boldsymbol{A}}$ is independent of ${\boldsymbol{\theta}_{\text{L}}}$ but in general depends on ${\boldsymbol{\theta}_{\text{NL}}}$.
  \item The priors on the two sets of parameters are independent, i.e. $\pi({\boldsymbol{\theta}})=\pi_{\text{L}}({\boldsymbol{\theta}_{\text{L}}})\ \pi_{\text{NL}}({\boldsymbol{\theta}_{\text{NL}}})$.
  \item $\pi_{\text{L}}({\boldsymbol{\theta}_{\text{L}}})$ is Gaussian.\footnote{The distribution $\pi_{\text{L}}({\boldsymbol{\theta}_{\text{L}}})$ could also be an improper uniform prior, i.e. $\pi_{\text{L}}({\boldsymbol{\theta}_{\text{L}}})=1$, which would imply that any value of ${\boldsymbol{\theta}_{\text{L}}}$ is equally likely.}
\end{enumerate}

\subsection{Efficient calculation of marginal posterior}

Under the conditions described in Section~\ref{sec:conditions-of-AMLP}, it is possible to characterize the distribution only through numerical sampling of ${\boldsymbol{\theta}_{\text{NL}}}$ instead of ${\boldsymbol{\theta}}$ as a whole. This can greatly reduce the number of dimensions being explored by the sampling algorithm, which often drastically reduces the time necessary to effectively converge to the posterior distribution. To do this, we must find the posterior marginal distribution of ${\boldsymbol{\theta}_{\text{NL}}}$, which is given by
\begin{subequations}
\begin{align}
  p_{\text{NL}}({\boldsymbol{\theta}_{\text{NL}}}) &= \int p({\boldsymbol{\theta}_{\text{NL}}},{\boldsymbol{\theta}_{\text{L}}})\  d{\boldsymbol{\theta}_{\text{L}}}, \\
  &= \pi_{\text{NL}}({\boldsymbol{\theta}_{\text{NL}}})\ {\mathcal{L}_{\text{eff}}}({\boldsymbol{\theta}_{\text{NL}}}). \label{eq:marginal-posterior-intermediate}
\end{align}
\end{subequations}
Here, we have implicitly defined the effective likelihood function ${\mathcal{L}_{\text{eff}}}({\boldsymbol{\theta}_{\text{NL}}})$ as
\begin{equation}
  {\mathcal{L}_{\text{eff}}}({\boldsymbol{\theta}_{\text{NL}}}) = \int \pi_{\text{L}}({\boldsymbol{\theta}_{\text{L}}})\ {\mathcal{L}}({\boldsymbol{\theta}_{\text{NL}}},{\boldsymbol{\theta}_{\text{L}}})\ d{\boldsymbol{\theta}_{\text{L}}}.
\end{equation}
This integral can be computed directly by noting that, under the conditions mentioned above, the integrand is proportional to a multivariate Gaussian distribution\footnote{In particular, the integrand is proportional to the conditional distribution $\Pr[{\boldsymbol{\theta}_{\text{L}}}|{\boldsymbol{y}},{\boldsymbol{\theta}_{\text{NL}}}]$, which should be read as the probability density of the linear parameters, ${\boldsymbol{\theta}_{\text{L}}}$, conditioned on the data, ${\boldsymbol{y}}$, and the nonlinear parameters, ${\boldsymbol{\theta}_{\text{NL}}}$.} in ${\boldsymbol{\theta}_{\text{L}}}$ with mean ${\boldsymbol{\mu}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})={\text{E}}[{\boldsymbol{\theta}_{\text{L}}}|{\boldsymbol{\theta}_{\text{NL}}}]$ and covariance ${\boldsymbol{\Sigma}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})={\text{Cov}}[{\boldsymbol{\theta}_{\text{L}}}|{\boldsymbol{\theta}_{\text{NL}}}]$. The conditional mean ${\boldsymbol{\mu}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})$ and covariance ${\boldsymbol{\Sigma}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})$ can be computed quickly because the model is linear when conditioned on ${\boldsymbol{\theta}_{\text{NL}}}$. We find
\begin{multline}
  {\mathcal{L}_{\text{eff}}}({\boldsymbol{\theta}_{\text{NL}}}) = \left\Vert 2\pi {\boldsymbol{\Sigma}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})\right\Vert^{1/2}\\ \times \pi_{\text{L}}\big({\boldsymbol{\mu}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})\big)\ {\mathcal{L}}\big({\boldsymbol{\theta}_{\text{NL}}},{\boldsymbol{\mu}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})\big),
\end{multline}
where $\left\Vert\cdot\right\Vert$ indicates the determinant. Plugging this into Equation~\ref{eq:marginal-posterior-intermediate}, we can write
\begin{multline}
  p_{\text{NL}}({\boldsymbol{\theta}_{\text{NL}}}) = \left\Vert 2\pi{\boldsymbol{\Sigma}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})\right\Vert^{1/2}\ \pi_{\text{NL}}({\boldsymbol{\theta}_{\text{NL}}})\\\times \pi_{\text{L}}\big({\boldsymbol{\mu}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})\big)\ {\mathcal{L}}\big({\boldsymbol{\theta}_{\text{NL}}},{\boldsymbol{\mu}_{\text{L}}}({\boldsymbol{\theta}_{\text{NL}}})\big). \label{eq:marginal-posterior}
\end{multline}

\subsection{Recreating sample from joint posterior}

The distribution $p_{\text{NL}}$ from Equation~\ref{eq:marginal-posterior} can be sampled numerically with a method such as Markov Chain Monte Carlo (MCMC) or nested sampling to yield a sequence of values $\{{\boldsymbol{\theta}}_{\text{NL}}^{(1)},{\boldsymbol{\theta}}_{\text{NL}}^{(2)},\ldots,{\boldsymbol{\theta}}_{\text{NL}}^{(N)}\}$.\footnote{For this work, we use the \texttt{emcee} code described in \cite{DFM:13}, although we also introduced a custom Metropolis Hastings MCMC sampler in \cite{Rapetti:20} that can be found in the \texttt{pylinex} code \citep{pylinex:21}.} From this sequence, we can create a sample of ${\boldsymbol{\theta}_{\text{L}}}$ by sampling the Gaussian conditional distributions. For each integer $k$ satisfying $1\le k\le N$, we sample\footnote{In this paper, ${\mathcal{N}}({\boldsymbol{v}},{\boldsymbol{Z}})$ and ${\mathcal{CN}}({\boldsymbol{v}},{\boldsymbol{Z}})$ denote normal and complex normal distributions with mean ${\boldsymbol{v}}$ and covariance ${\boldsymbol{Z}}$.}
\begin{equation}
  {\boldsymbol{\theta}}_{\text{L}}^{(k)}\sim {\mathcal{N}}\bigg({\boldsymbol{\mu}_{\text{L}}}\big({\boldsymbol{\theta}}_{\text{NL}}^{(k)}\big),{\boldsymbol{\Sigma}_{\text{L}}}\big({\boldsymbol{\theta}}_{\text{NL}}^{(k)}\big)\bigg). \label{eq:sampling-conditional-distribution}
\end{equation}
Then, if we define ${\boldsymbol{\theta}}^{(k)}$ as $\begin{bmatrix} {\boldsymbol{\theta}}_{\text{NL}}^{(k)} \\ {\boldsymbol{\theta}}_{\text{L}}^{(k)} \end{bmatrix}$, $\left\{ {\boldsymbol{\theta}}^{(1)}, {\boldsymbol{\theta}}^{(2)}, \ldots, {\boldsymbol{\theta}}^{(N)} \right\}$ is a sample from $p({\boldsymbol{\theta}})$.\footnote{Note that if we only want a sample of ${\boldsymbol{\theta}_{\text{L}}}$, then instead of sampling the distribution in Equation~\ref{eq:sampling-conditional-distribution} only once per $k$-value, we can sample it $M$ times, leading to a total sample of $N\times M$ realizations of ${\boldsymbol{\theta}_{\text{L}}}$.}

\subsection{Choosing how to split parameters}

One remaining question about the AMLP technique is how to choose which parameters to include in ${\boldsymbol{\theta}_{\text{NL}}}$ and which to include in ${\boldsymbol{\theta}_{\text{L}}}$. For example, if the data vector consists of the product of two components, $A$ and $B$, which each have linear models with parameter vectors ${\boldsymbol{\varepsilon}}$ and ${\boldsymbol{\rho}}$, respectively, then there are multiple ways of applying AMLP. For instance, $\begin{bmatrix} {\boldsymbol{\theta}_{\text{NL}}} \\ {\boldsymbol{\theta}_{\text{L}}} \end{bmatrix}$ can be set to either $\begin{bmatrix} {\boldsymbol{\varepsilon}} \\ {\boldsymbol{\rho}} \end{bmatrix}$ or $\begin{bmatrix} {\boldsymbol{\rho}} \\ {\boldsymbol{\varepsilon}} \end{bmatrix}$. To decide which of these choices is the best, we examine the dimensions of ${\boldsymbol{\varepsilon}}$ and ${\boldsymbol{\rho}}$. The longer of the two vectors should be marginalized over (i.e. set to ${\boldsymbol{\theta}_{\text{L}}}$), whereas the other should be explored numerically (i.e. set to ${\boldsymbol{\theta}_{\text{NL}}}$).

\subsection{Implementation in the \texttt{pylinex} code} \label{sec:pylinex-implementation}

Through the \texttt{MultiConditionalFitModel} class in the \texttt{pylinex} code \citep{pylinex:21}, we can perform AMLP with any model that can be written using a combination of products and sums of constituent submodels (such as the one described in Section~\ref{sec:full-model} and shown in Figure~\ref{fig:specific-model-tree}) and any split of parameters into ${\boldsymbol{\theta}_{\text{NL}}}$ and ${\boldsymbol{\theta}_{\text{L}}}$ that satisfies the condition that the model at constant ${\boldsymbol{\theta}_{\text{NL}}}$ is linear in ${\boldsymbol{\theta}_{\text{L}}}$.

\section{Fitting data with foreground, signal, and receiver} \label{sec:fitting-strategy}

In this section, we lay out the general framework we use to apply AMLP to 21-cm observations with receiver effects. Here, it is important to note that when we refer to linear and nonlinear, we are referring to parameters and not the receiver, which is assumed to be operated in a linear regime.

\subsection{Modeling observations} \label{sec:full-model}

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.45\textwidth]{model_tree.pdf}
  \caption{Visualization of the model tree corresponding to Equation~\ref{eq:block-model}. The ellipses (nodes with no children) represent the four components of the model, while the orange squares (branching nodes) represent their connections through either addition or multiplication. Blue ellipses indicate components (foreground and receiver offset) whose parameters will be marginalized at each step of an MCMC, i.e.~models whose parameters will be part of ${\boldsymbol{\theta}_{\text{L}}}$ as described in Section~\ref{sec:AMLP}. The green ellipse indicates that the gain parameters are included in ${\boldsymbol{\theta}_{\text{NL}}}$ and are explored via MCMC. The ellipse in magenta indicates a component (signal) whose parameters are marginalized over if they are linear, but explored via MCMC if they are nonlinear.} \label{fig:specific-model-tree}
\end{figure}

The full observation equation of the temperature data $T_{\nu,t}$ as a function of frequency $\nu$ and time $t$ is given by
\begin{equation}
  T_{\nu,t} = G_\nu\times(f_{\nu,t}+s_\nu+o_\nu) + n_{\nu,t},  \label{eq:observation-equation}
\end{equation}
where $G_\nu$ is the receiver's multiplicative bias (gain),\footnote{$G_\nu$ is the power gain, which is equal to $|g_\nu|^2$, where $g_\nu$ is the complex voltage gain. See Appendix~\ref{app:noise-level}.} $o_\nu$ is the receiver's additive bias (offset or noise temperature), $s_\nu$ is the global 21-cm signal, $f_{\nu,t}$ is the frequency- and time-dependent beam-weighted foreground emission, and $n_{\nu,t}$ is a random realization of a Gaussian noise vector. In order to fit the data, we form the model
\begin{multline}
  {\mathcal{M}}_{T_{\nu,t}}({\boldsymbol{\zeta}}_G,{\boldsymbol{\zeta}}_f,{\boldsymbol{\zeta}}_s,{\boldsymbol{\zeta}}_o) =\\ {\mathcal{M}}_{G_\nu}({\boldsymbol{\zeta}}_G)\times\left[{\mathcal{M}}_{f_{\nu,t}}({\boldsymbol{\zeta}}_f) + {\mathcal{M}}_{s_\nu}({\boldsymbol{\zeta}}_s) + {\mathcal{M}}_{o_\nu}({\boldsymbol{\zeta}}_o)\right], \label{eq:block-model}
\end{multline}
where ${\mathcal{M}}_X({\boldsymbol{\zeta}}_X)$ is the model of the component $X$ and ${\boldsymbol{\zeta}}_X$ is the vector of parameters of that model. Figure~\ref{fig:specific-model-tree} shows a tree representation of this model. As described in Section~\ref{sec:pylinex-implementation}, AMLP can be applied to any model that can be represented by this kind of sum-product tree using the \texttt{pylinex} code.

\subsection{Fitting strategy}

The main goal of this fit is a set of constraints on the parameters of a nonlinear signal model. While one may attempt to fit the desired nonlinear signal model directly, the posterior distribution is often difficult to numerically explore if the sampling algorithm is not started in a narrow region around the maximum. Therefore, as is laid out in the following sections, we fit the data with a linear signal model first to home in on a specific set of spectral shapes and follow up with a nonlinear signal model fit initialized near this region.

\subsubsection{Form models} \label{sec:form-models}

  Before we fit the data, we must form models for each of the four components.

  \begin{itemize}
    \item Generate a model to fit the receiver gain ${\mathcal{M}}_{G_\nu}({\boldsymbol{\zeta}}_G)$, which can be either linear or nonlinear.
    \item Create linear model of the receiver noise temperature, ${\mathcal{M}}_{o_\nu}({\boldsymbol{\zeta}}_o)$, either from a priori knowledge or from a training set of simulated noise temperature curves.\footnote{In this paper, we assume the model of the noise temperature is linear, although modifications to this method for the case of a nonlinear noise temperature model are discussed in Section~\ref{sec:nonlinear-offset-model}.}
    \item Create linear model of beam-weighted foreground spectra, ${\mathcal{M}}_{f_{\nu,t}}({\boldsymbol{\zeta}}_f)$, from a training set simulated using reasonably varied beams and foreground maps.
    \item Choose a nonlinear signal model, ${\mathcal{M}}_{s_\nu}^{\text{NL}}({\boldsymbol{\zeta}}_s^{\text{NL}})$.
    \item Create linear model, ${\mathcal{M}}_{s_\nu}^{\text{L}}({\boldsymbol{\zeta}}_s^{\text{L}})$, of the signal using a training set generated from the chosen nonlinear signal model. \label{item:linear-signal-model}
  \end{itemize}

  \subsubsection{Fit using the linear signal model} \label{sec:fit-using-linear-signal-model}

  Perform AMLP using the model of Equation~\ref{eq:block-model} with the linear signal model, i.e. ${\mathcal{M}}_{s_\nu}({\boldsymbol{\zeta}}_s)={\mathcal{M}}_{s_\nu}^{\text{L}}({\boldsymbol{\zeta}}_s^{\text{L}})$. Split the parameters through ${\boldsymbol{\theta}_{\text{NL}}}={\boldsymbol{\zeta}}_G$ and ${\boldsymbol{\theta}_{\text{L}}}=\begin{bmatrix} {\boldsymbol{\zeta}}_s^{\text{L}} \\ {\boldsymbol{\zeta}}_f \\ {\boldsymbol{\zeta}}_o \end{bmatrix}$. In this fit, only the parameters of the gain model must be explored numerically, so it should be completed very quickly and will produce a fast estimate of the signal in frequency-temperature space. The walkers of the MCMC should be initialized through a sample of the prior distribution on ${\boldsymbol{\zeta}}_G$.

  \subsubsection{Prepare to fit using the nonlinear signal model} \label{sec:prepare-to-fit-using-nonlinear-signal-model}

  The final fit will need to numerically explore both ${\boldsymbol{\zeta}}_G$ and ${\boldsymbol{\zeta}}_s^{\text{NL}}$. An estimate of the distribution of ${\boldsymbol{\zeta}}_G$ should be available from the fit done in the previous step; but, we must find an estimate of the distribution of ${\boldsymbol{\zeta}}_s^{\text{NL}}$.\footnote{Item~\ref{item:linear-signal-model} from Section~\ref{sec:form-models} and the steps laid out in Sections~\ref{sec:fit-using-linear-signal-model}~and~\ref{sec:prepare-to-fit-using-nonlinear-signal-model} can be skipped if the nonlinear signal model is very fast (e.g. if many signals can be evaluated in one second). These steps are meant to provide a starting point for the final fit so that the nonlinear signal model parameter distribution can be achieved in the fewest possible steps.} To do so, we perform a least square fit to the linear signal model parameter mean vector from the first AMLP fit described in Section~\ref{sec:fit-using-linear-signal-model} in mode coefficient space, i.e. we define the first guess nonlinear parameter vector ${\boldsymbol{\zeta}}_{s,\text{guess}}^{\text{NL}}$ through
  \begin{equation*}
    {\boldsymbol{\zeta}}_{s,\text{guess}}^{\text{NL}}=\underset{{\boldsymbol{\zeta}}_{s}^{\text{NL}}}{\text{argmin}}\left|{\boldsymbol{\Lambda}}_{\text{L}}^{-1/2} \left[{\boldsymbol{\mathcal{M}}}_{s}^{\text{L},-1}({\boldsymbol{\mathcal{M}}}_{s}^{\text{NL}}({\boldsymbol{\zeta}}_s^{\text{NL}}))-\overline{{\boldsymbol{\zeta}}}_s^{\text{L}}\right]\right|^2,
  \end{equation*}
  where $\overline{{\boldsymbol{\zeta}}}_s^{\text{L}}={\text{E}}[{\boldsymbol{\zeta}}_s^{\text{L}}]$, ${\boldsymbol{\Lambda}}_{\text{L}}={\text{Cov}}[{\boldsymbol{\zeta}}_s^{\text{L}}]$, and ${\boldsymbol{\mathcal{M}}}_s^{\text{L},-1}$ is the pseudo-inverse of ${\boldsymbol{\mathcal{M}}}_s^{\text{L}}$.\footnote{Here, ${\boldsymbol{\mathcal{M}}}_{s}$ is a vectorized form of ${\mathcal{M}}_{s_\nu}$. By pseudo-inverse, we mean the function that takes in a spectrum ${\boldsymbol{s}}$ and outputs the parameters ${\boldsymbol{\zeta}}_s^{\text{L}}$ that minimize the difference between ${\boldsymbol{s}}$ and ${\boldsymbol{\mathcal{M}}}_{s}^{\text{L}}({\boldsymbol{\zeta}}_s^{\text{L}})$, as measured by the noise covariance, ${\boldsymbol{C}}$. If ${\boldsymbol{\mathcal{M}}}_{s}^{\text{L}}({\boldsymbol{\zeta}}_s^{\text{L}})={\boldsymbol{a}}+{\boldsymbol{F}}{\boldsymbol{\zeta}}_s^{\text{L}}$, then ${\boldsymbol{\mathcal{M}}}_s^{L,-1}({\boldsymbol{s}})=({\boldsymbol{F}}^T{\boldsymbol{C}}^{-1}{\boldsymbol{F}})^{-1}{\boldsymbol{F}}^T{\boldsymbol{C}}^{-1}({\boldsymbol{s}}-{\boldsymbol{a}})$.} We also derive a covariance of the guess distribution from the Fisher information, i.e.
  \begin{equation*}
    {\boldsymbol{\Lambda}}_{\text{NL},\text{guess}} = ({\boldsymbol{Z}}^T{\boldsymbol{\Lambda}}_{\text{L}}^{-1}{\boldsymbol{Z}})^{-1},
  \end{equation*}
  where ${\boldsymbol{Z}}=\left[\frac{\partial}{\partial{\boldsymbol{\zeta}}_s^{\text{NL}}}{\boldsymbol{\mathcal{M}}}_s^{\text{L},-1}({\boldsymbol{\mathcal{M}}}_s^{\text{NL}}({\boldsymbol{\zeta}}_s^{\text{NL}}))\right]_{{\boldsymbol{\zeta}}_{s}^{\text{NL}}={\boldsymbol{\zeta}}_{s,\text{guess}}^{\text{NL}}}$. To initialize walkers for the MCMC of the final nonlinear signal model fit, we draw a sample of gain parameters from the chains of the MCMC fit from Section~\ref{sec:fit-using-linear-signal-model} and we draw the nonlinear signal parameters from the distribution defined by the mean and covariance given in this section, i.e.
  \begin{equation*}
    {\boldsymbol{\zeta}}_{s}^{\text{NL}}\sim{\mathcal{N}}({\boldsymbol{\zeta}}_{s,\text{guess}}^{\text{NL}},{\boldsymbol{\Lambda}}_{\text{NL},\text{guess}}).
  \end{equation*}

  \subsubsection{Fit using the nonlinear signal model} \label{sec:fit-using-nonlinear-signal-model}

  Finally, we perform AMLP using Equation~\ref{eq:block-model} with the nonlinear signal model, i.e. ${\mathcal{M}}_{s_\nu}({\boldsymbol{\zeta}}_s)={\mathcal{M}}_{s_\nu}^{\text{NL}}({\boldsymbol{\zeta}}_s^{\text{NL}})$, and splitting the parameters into ${\boldsymbol{\theta}_{\text{NL}}}=\begin{bmatrix} {\boldsymbol{\zeta}}_G \\ {\boldsymbol{\zeta}}_s^{\text{NL}} \end{bmatrix}$ and ${\boldsymbol{\theta}_{\text{L}}}=\begin{bmatrix} {\boldsymbol{\zeta}}_f \\ {\boldsymbol{\zeta}}_o \end{bmatrix}$. Because more parameters are explored numerically, this fit will naturally take longer than the fit with the linear signal model; but, it will ultimately output the desired constraints on the signal parameters.

\subsection{Noise level}

The only aspect of the fitting strategy left to define is how to determine the noise covariance ${\boldsymbol{C}}$. As is described in Appendix~\ref{app:noise-level}, the noise in each frequency bin is independent and nearly equal to the data divided by the dynamic range factor $\sqrt{\Delta\nu\ \Delta t}$, where $\Delta\nu$ is the channel width in Hz and $\Delta t$ is the integration time per spectrum in $s$, i.e.
\begin{equation}
  C_{jk} = \frac{(y_j)^2}{\Delta\nu\ \Delta t} \times \begin{cases} 1 & j=k \\ 0 & j\ne k \end{cases}. \label{eq:noise-covariance}
\end{equation}
For the fits in this paper, in addition to being used to define ${\boldsymbol{C}}$, this covariance is used to simulate the random noise $n_{\nu,t}$ from Equation~\ref{eq:observation-equation}. We use $\Delta\nu=1$ MHz and $\Delta t = 800$ hr.

\section{Models} \label{sec:models}

While Section~\ref{sec:fitting-strategy} laid out the general strategy for fitting observations while including receiver models, in this section, we put forth the specific models used in this work, following the tasks from Section~\ref{sec:form-models}.

\begin{table}[t!]
  \centering
  \caption{Gain model parameters} \label{tab:gain-model-parameters}
  \begin{tabular*}{0.47\textwidth}{c @{\extracolsep{\fill}} ccc}
    \hline
    \hline
    Parameter & Unit & Input & Prior \\
    \hline
    $10\log_{10}{\left[1 + \left(\varepsilon^{\text{HP}}\right)^2\right]}$ & dB & $0.25$ & $\text{Unif}(0, 1)$ \\
    $10\log_{10}{\left[1 + \left(\varepsilon^{\text{LP}}\right)^2\right]}$ & dB & $0.25$ & $\text{Unif}(0, 1)$ \\
    $\nu_0^{\text{HP}}$ & MHz & 41.1 & $\text{Unif}(35, 50)$ \\
    $\nu_0^{\text{LP}}$ & MHz & 116.8 & $\text{Unif}(110, 125)$ \\
    \hline
  \end{tabular*}
    \vspace{0.8ex}

    {\raggedright \textbf{Notes.} The transformations of the $\varepsilon$ parameters in the first two rows are the levels of the in-band ripples in dB. \par}
\end{table}

\subsection{Gain model} \label{sec:gain-model}

In this paper, we use a gain model that is the product of a low-pass (LP) filter and a high-pass (HP) filter, i.e.
\begin{multline}
  {\mathcal{M}}_{G_\nu}\left(\nu_0^{\text{HP}}, \varepsilon^{\text{HP}}, \nu_0^{\text{LP}}, \varepsilon^{\text{LP}}\right) =\\ {\mathcal{M}}_{G^{\text{HP}}_\nu}\left(\nu_0^{\text{HP}}, \varepsilon^{\text{HP}}\right)\times {\mathcal{M}}_{G^{\text{LP}}_\nu}\left(\nu_0^{\text{LP}}, \varepsilon^{\text{LP}}\right).
\end{multline}
The individual filters are taken to have Chebyshev type I transfer functions, meaning that their gains are given by
\begin{subequations}
\begin{align}
  {\mathcal{M}}_{G_\nu^{\text{HP}}}\left(\nu_0^{\text{HP}},\varepsilon^{\text{HP}}\right) &= \frac{1}{1+\left[\varepsilon^{\text{HP}}\ T_n\left(\frac{\nu_0^{\text{HP}}}{\nu}\right)\right]^2}, \\
  {\mathcal{M}}_{G_\nu^{\text{LP}}}\left(\nu_0^{\text{LP}},\varepsilon^{\text{LP}}\right) &= \frac{1}{1+\left[\varepsilon^{\text{LP}}\ T_n\left(\frac{\nu}{\nu_0^{\text{LP}}}\right)\right]^2},
\end{align}
\end{subequations}
where $T_n$ is the $n^{\text{th}}$ order Chebyshev polynomial.\footnote{Chebyshev polynomials are defined by $T_n(\cos{\theta})=\cos{n\theta}$. We use $n=9$.} The input gain curve in the simulations (see the red curve in Figure~\ref{fig:nonlinear-uncertainties}) is generated using the parameter values shown in Table~\ref{tab:gain-model-parameters}. The $\varepsilon$ parameters determine the level of ripple inside the band. We choose values that generate a 0.25 dB ripple in the pass-band from both the high-pass and low-pass filters, leading to a total ripple of roughly 0.5 dB. The reference frequencies $\nu_0^{\text{HP}}$ and $\nu_0^{\text{LP}}$ were chosen so that the gain is close to $\frac{1}{2}$ ($-3$ dB) at the edges of the observed band (40-120 MHz).

\subsection{Nonlinear signal model} \label{sec:nonlinear-signal-model}

Here, we use the so-called turning point model of the signal (used also in Paper II), which is an interpolation between turning points (i.e. local extrema) of the 21-cm signal. These turning points are labeled A-E and are illustrated and described in Figure~\ref{fig:turning-point-explanation}. The free parameters of the model are the frequencies and brightness temperatures of the turning points B-D, whereas $(\nu_A, T_A)$ is fixed to the $\Lambda$CDM value of $(18\text{ MHz}, -40\text{ mK})$ and $(\nu_E,T_E)$ is fixed to $(180\text{ MHz}, 0\text{ mK})$.\footnote{We justify fixing $\nu_E$ by noting that it will not be constrained by observations between 40-120 MHz.} The model is a cubic spline between the turning points. In order to force the turning points to be extrema (i.e. have derivative zero), each turning point uses two spline knots placed at the same temperature and 20 kHz apart symmetrically around the turning point frequency. In addition to turning points A-E, there are two knots placed at 0 K and $10\pm10$ kHz to force the signal to approach zero smoothly at very large redshifts.~The model always evaluates to 0 K at frequencies above that of turning point E. The red rectangles in Figure~\ref{fig:turning-point-explanation} and the last column of Table~\ref{tab:turning-point-parameters} indicate the priors we place on the six varying parameters.

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.47\textwidth]{turning_point_explanation.pdf}
  \caption{Illustration of the turning point signal model. The signal is a cubic spline using knots at the labeled turning points to enforce that they are local extrema. The red rectangles around turning points B-D indicate the prior distribution of their frequencies and temperatures (see Table~\ref{tab:turning-point-parameters}). Turning points A and E are fixed because they will not be constrained by points within the observed band (40-120 MHz) marked by the vertical dashed lines.} \label{fig:turning-point-explanation}
\end{figure}

\begin{table}[t!]
    \centering
    \caption{Turning point signal model parameters} \label{tab:turning-point-parameters}
    \begin{tabular*}{0.47\textwidth}{c @{\extracolsep{\fill}} ccc}
        \hline
        \hline
        Parameter & Units & Input & Prior \\
        \hline
        $\nu_A$ & MHz & 18 & Fixed \\
        $T_A$ & mK & -40 & Fixed \\
        $\nu_B$ & MHz & 45 & Unif(30, 80) \\
        $T_B$ & mK & -3 & Unif(-5, 0) \\
        $\nu_C$ & MHz & 80 & Unif(60, 120) \\
        $T_C$ & mK & -250 & Unif(-350, -100) \\
        $\nu_D$ & MHz & 115 & Unif(100, 150) \\
        $T_D$ & mK & 15 & Unif(0, 25) \\
        $\nu_E$ & MHz & 180 & Fixed \\
        \hline
    \end{tabular*}
    \vspace{0.8ex}

    {\raggedright \textbf{Notes.} The frequencies of adjacent turning points are also constrained to differ by at least 10 MHz. \par}
\end{table}

\subsection{Linear models} \label{sec:linear-models}

As in Paper I, we form models through decomposition of training sets. For doing so in this paper, though, we apply a slightly modified procedure that employs a principal component analysis-like formalism described in Appendix~\ref{app:pca}. Here, we summarize how we form the training sets in the case of each component. In the following, we denote an individual training set curve of component $X$ as $b_X$.

\begin{itemize}
  \item \textit{Noise temperature training set}: ${\mathcal{M}}_{o_\nu}$ is formed from a training set where each curve is a line parameterized by its low- and high-frequency endpoints, $o_L$ and $o_H$, i.e.
  \begin{equation}
    b_{o_\nu}=o_L\left[\frac{(120\text{ MHz})-\nu}{80\text{ MHz}}\right]+o_H\left[\frac{\nu-(40\text{ MHz})}{80\text{ MHz}}\right].
  \end{equation}
  To generate many curves, we assume that $\frac{o_L}{1\text{ K}}\sim\text{Unif}(115, 120)$ and $\frac{o_H}{1\text{ K}}\sim\text{Unif}(110, 115)$.

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.47\textwidth]{fwhm_training_set.pdf}
  \caption{Ten percent of the training set of 500 beam full widths at half maximum (FWHM). The curves are quadratic Legendre polynomials with coefficient distributions given in Table~\ref{tab:legendre-coefficients}.} \label{fig:foreground-training-set}
\end{figure}

\renewcommand{\tabcolsep}{8mm}
\begin{table}[t!]
    \centering
    \caption{Parameters of Legendre coefficient distributions} \label{tab:legendre-coefficients}
    \begin{tabular}{ccc}
        \hline
        \hline
        $k$ & $\mu_k$ & $\sigma_k$ \\
        & $[^\circ]$ & $[^\circ]$ \\
        \hline
        0 & 70 & 10 \\
        1 & -20 & 5 \\
        2 & 0 & 5 \\
        \hline
        \multicolumn{3}{l}{\textbf{Note}: See Equations~\ref{eq:legendre-polynomials},~\ref{eq:legendre-polynomial-definitions},~and~\ref{eq:coefficient-distribution}.}
    \end{tabular}
\end{table}

  \item \textit{Beam-weighted foreground training set}: A general beam-weighted foreground training set should be created from two sources: variations in the antenna beam and in foreground emission. In this paper, as in the rest in the series, however, we use many beams and one foreground map,\footnote{See \cite{Hibbard:20} for an initial attempt at varying the foreground in this formalism.} with the latter given by \cite{Remazeilles:15} \citep[which used the data of][]{Haslam:82} scaled with a spectral index of -2.5.~The beams are angular Gaussians, i.e.~they satisfy $B(\nu,\theta,\phi)\propto \exp{\left\{-\frac{\theta^2}{2[\alpha(\nu)]^2}\right\}}$, where $\theta$ and $\phi$ are the polar and azimuthal spherical coordinate angles, respectively. The scale $\alpha$ is a function of frequency, $\nu$, so that beam chromaticity can be included in the analysis robustly. The full width at half maximum (FWHM), given by $\text{FWHM}(\nu) = \sqrt{8\ln{2}}\ \alpha(\nu)$, is varied between training set curves and is generated by quadratic polynomials in frequency. To control magnitude variations in each order simply, we use second-order Legendre polynomials:
\begin{equation}
  \text{FWHM}(\nu) = \sum_{k=0}^2 a_kL_k\left(\frac{\nu-\nu_0}{\delta\nu}\right), \label{eq:legendre-polynomials}
\end{equation}
where $\nu_0=(\nu_{\text{max}} + \nu_{\text{min}})/2$ is the average frequency, $\delta\nu=(\nu_{\text{max}}-\nu_{\text{min}})/2$ is half the width of the frequency band, and
\begin{equation}
    L_0(x) = 1,\ \
    L_1(x) = x,\ \
    L_2(x) = \frac{3x^2-1}{2}. \label{eq:legendre-polynomial-definitions}
\end{equation}
In our case, $\nu_{\text{min}}=40$ MHz and $\nu_{\text{max}}=120$ MHz, so $\nu_0=80$ MHz and $\delta\nu=40$ MHz. To seed the beam variations in our training set, we draw $a_0$, $a_1$, and $a_2$ from independent normal distributions,
\begin{equation}
  a_k \sim {\mathcal{N}}(\mu_k, \sigma_k^2), \label{eq:coefficient-distribution}
\end{equation}
with the means and standard deviations $\mu_k$ and $\sigma_k$ given in Table~\ref{tab:legendre-coefficients}. An extra constraint is applied to exclude $\text{FWHM}(\nu)$ curves which dip below $15^\circ$ in the $40-120$ MHz band. The resulting training set of $\text{FWHM}$ curves is shown in Figure~\ref{fig:foreground-training-set}.

  We simulate the beam-weighted foreground temperature for ten different bins in local sidereal time (LST) by smearing the map through LST before computing the spectra.\footnote{By smearing, we refer to averaging many rotations of the foreground map corresponding to LSTs between the edges of each LST bin.}~The spectra from the ten LST bins are then concatenated into the final data vector.
  \item \textit{Signal training set}: ${\mathcal{M}}^{\text{L}}_{s_\nu}$ is formed from a training set where each curve is created by applying the turning point nonlinear signal model ${\mathcal{M}}^{\text{NL}}_{s_\nu}$ to a sample from the prior distribution of its parameters (see Section~\ref{sec:nonlinear-signal-model}, Figure~\ref{fig:turning-point-explanation}, and Table~\ref{tab:turning-point-parameters}).
\end{itemize}

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

This section lays out the results of applying the general AMLP method described in Section~\ref{sec:AMLP} to the procedure put forth in Section~\ref{sec:fitting-strategy} for fitting sky-averaged radio data measured with a non-ideal receiver using the models laid out in Section~\ref{sec:models}. The initial fit using the linear signal model (following Section~\ref{sec:fit-using-linear-signal-model}) is shown in Section~\ref{sec:results-fit-using-linear-signal-model} and the final fit using the nonlinear signal model (following Section~\ref{sec:fit-using-nonlinear-signal-model}) is shown in Section~\ref{sec:results-fit-using-nonlinear-signal-model}.

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.43\textwidth]{linear_signal_posterior_sample.pdf}
  \caption{\textit{Top}: Signal realizations from the first MCMC fit with the linear signal model, i.e. ${\mathcal{M}}_{s_\nu}^{\text{L}}({\boldsymbol{\zeta}}_s^{\text{L}})$, where ${\boldsymbol{\zeta}}_s^{\text{L}}$ is sampled from the posterior distribution, are shown in black. The red line shows the input signal, which is clearly inside the uncertainty interval. \textit{Bottom}: Same as top with the mean of the black lines subtracted from each curve. The main source of variance in the signal estimate comes from modes similar to lines, which are nearly degenerate with the linear modes describing the noise temperature.} \label{fig:linear-signal-model-uncertainties}
\end{figure}

\subsection{Fit with linear signal model} \label{sec:results-fit-using-linear-signal-model}

The signal uncertainties implied by the first MCMC fit with the linear signal model are shown in Figure~\ref{fig:linear-signal-model-uncertainties}. The main mode of uncertainty is line-like, with a width of approximately 50 mK. This width comes from the similarity, or overlap, between the linear signal model and the models of the other components. For simplicity, here we will discuss the effect of a similarity between the signal and noise temperature models, but it is important to note that similarities between the signal model and the gain or beam-weighted foreground models would produce the same effect. The uncertainties implied by the posterior are designed to answer a fundamental question: what size shift in signal parameters $\delta{\boldsymbol{\zeta}}_s^{\text{L}}$ can lead to a spectral change that is compatible with the Gaussian noise distribution of the data when accounting for offsetting changes $\delta{\boldsymbol{\zeta}}_o$ in parameters of the noise temperature model? If all combinations of $\delta{\boldsymbol{\zeta}}_s^{\text{L}}$ and $\delta{\boldsymbol{\zeta}}_o$ produce orthogonal effects on the data, then the allowed size of $\delta{\boldsymbol{\zeta}}_s^{\text{L}}$ is determined by the noise alone, and the resulting signal uncertainties should essentially match the noise level of the data. On the other hand, if the effects of some $\delta{\boldsymbol{\zeta}}_s^{\text{L}}$ can be offset very closely by a corresponding change $\delta{\boldsymbol{\zeta}}_o^{\text{NL}}$, then the uncertainties will be larger to account for this overlap. Moreover, if both the signal and noise temperature models are linear (as they are in this first MCMC fit) and they have overlapping gradients, then the uncertainties can grow greatly or even diverge because if, for example, $\delta{\boldsymbol{\zeta}}_s^{\text{L}}$ can be exactly offset by $\delta{\boldsymbol{\zeta}}_o$, then $2\delta{\boldsymbol{\zeta}}_s^{\text{L}}$ can be exactly offset by $2\delta{\boldsymbol{\zeta}}_o$. The uncertainties shown in the bottom panel of Figure~\ref{fig:linear-signal-model-uncertainties} are caused by the fact that there is a direction in the ${\boldsymbol{\zeta}}_s^{\text{L}}$ parameter space that causes a line-like spectral feature. Since the noise temperature model used in this paper is a linear function of frequency, it can closely offset a line-like signal change, leading to large signal errors with a line-like shape.

The signals from this first MCMC fit are used to initialize the MCMC fit with the nonlinear signal model as described in Sections~\ref{sec:prepare-to-fit-using-nonlinear-signal-model}~and~\ref{sec:fit-using-nonlinear-signal-model}.

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.47\textwidth]{parameter_covariance.pdf}
  \caption{The full parameter covariance matrix arising from the final posterior distribution in the fit with the nonlinear signal model.~Each row and column is normalized so that diagonal elements of one correspond to a difference of one noise level in a sum of squares sense.~The blocks of the matrix corresponding to the parameters of different components (i.e. gain, signal, foreground, and noise temperature) are separated by black lines. The noise temperature is the component that covaries most with the signal because it appears similarly in the observation equation (Equation~\ref{eq:observation-equation}).} \label{fig:nonlinear-fit-parameter-covariance}
\end{figure}

\begin{figure*}[h!]
  \centering
  \includegraphics[width=0.423\textwidth]{nonlinear_signal_posterior_sample.pdf}
  \hspace{1cm}
  \includegraphics[width=0.423\textwidth]{foreground_posterior_sample.pdf}
  \includegraphics[width=0.423\textwidth]{gain_posterior_sample.pdf}
  \hspace{1cm}
  \includegraphics[width=0.423\textwidth]{offset_posterior_sample.pdf}
  \caption{Each plot shows 100 curves (black lines) from the final MCMC-explored posterior distribution with the nonlinear signal model alongside the input (red line) for the given component. The lower of each pair of plots shows the uncertainties on that component by subtracting the mean of the black curves from the plot above it. The two plots in the top left corner (A) show the signal uncertainties, which are much improved from their counterparts in Figure~\ref{fig:linear-signal-model-uncertainties}. The beam-weighted foreground sample is shown in the top right pair of plots (B), where the 10 spectra included in the data vector are concatenated. The receiver gain and noise temperature samples are shown in the bottom left (C) and bottom right (D) plots, respectively.} \label{fig:nonlinear-uncertainties}
\end{figure*}

\begin{figure*}[t!]
  \centering
  \includegraphics[width=0.95\textwidth]{signal_triangle_plot.pdf}
  \caption{Triangle plot showing marginal distributions on the turning point signal parameters. The univariate marginal distributions of each parameter are shown on the diagonal, while the 95\% confidence intervals of the bivariate marginal distributions are shown in the off-diagonal plots. The dashed lines show the input parameters, which are inside the confidence intervals in each plot. While the temperatures of turning points B and D are not constrained (in a univariate sense) at the 95\% level beyond the narrow priors, the other parameters are tightly constrained.} \label{fig:triangle-plot}
\end{figure*}

\subsection{Fit with nonlinear signal model} \label{sec:results-fit-using-nonlinear-signal-model}

Figure~\ref{fig:nonlinear-fit-parameter-covariance} shows the posterior parameter covariance of the MCMC fit with the nonlinear signal model.~Each row and column is scaled in such a way that a value of unity in the diagonal elements corresponds to a variation of one noise level. The signal covaries most with the noise temperature because the two components appear similarly in the data, i.e. the same in each spectrum, and must be separated based on their spectral shape alone.

The final MCMC fit yields our best mean and uncertainty estimates of the signal, foreground, receiver gain, and receiver noise temperature, which are shown in Figure~\ref{fig:nonlinear-uncertainties}. The signal sample shown in Figure~\ref{fig:nonlinear-uncertainties}A is greatly improved over the signal sample found using the linear signal model (Figure~\ref{fig:linear-signal-model-uncertainties}) for multiple reasons:
\begin{enumerate}
  \item The turning point signal model has fewer freely varying parameters (6) than the linear signal model (33).
  \item Due to it being nonlinear, the gradient of the turning point signal model is not constant, meaning that a given variation $\delta{\boldsymbol{\zeta}}_s^{\text{NL}}$ in the parameter vector affects the signal differently based on the point ${\boldsymbol{\zeta}}_s^{\text{NL}}$ around which the parameters are varying. Using the terminology of Section~\ref{sec:results-fit-using-linear-signal-model}, this means that even if a signal parameter change $\delta{\boldsymbol{\zeta}}_s^{\text{NL}}$ can be exactly offset by a noise temperature parameter change of $\delta{\boldsymbol{\zeta}}_o$, it is generally not true that a signal parameter change of $2\delta{\boldsymbol{\zeta}}_s^{\text{NL}}$ can be offset by a noise temperature parameter change of $2\delta{\boldsymbol{\zeta}}_o$.
\end{enumerate}

Figure~\ref{fig:nonlinear-uncertainties}B shows the foreground sample from the final MCMC fit, with all of the spectra concatenated. The uncertainties encompass the input foreground, shown with the red line, and have widths of 20-60 mK at low frequencies and a few mK at high frequencies. The beam-weighted foreground must be fit to this level to extract the signal without bias.\footnote{In this section, we use ``bias'' to refer to a disparity between a) the difference between the input curve and fitted curve of a particular component and b) the posterior uncertainty band of that component. This means that fits where the mean of the posterior sample is different than the input curve, which is almost always the case, can be unbiased as long as the uncertainties accurately characterize the difference.}

Figure~\ref{fig:nonlinear-uncertainties}C~and~\ref{fig:nonlinear-uncertainties}D show the posterior samples of the receiver gain and noise temperature, respectively. The gain is fit down to the sub-parts per million (ppm) level due to the fact that the dynamic range of the measurement is $\frac{\sqrt{C_{kk}}}{y_k}=\frac{1}{\sqrt{\Delta\nu\ \Delta t}}\approx 0.6$ ppm (see Equation~\ref{eq:noise-covariance} and Appendix~\ref{app:noise-level}) and the noise temperature is fit down to the 10 mK level, which is necessary to fit the signal down to the same level.

A triangle plot showing the univariate and bivariate marginal distributions of the turning point signal parameters is shown in Figure~\ref{fig:triangle-plot}. The temperatures of turning points B and D are not constrained tighter than their narrow priors (see Table~\ref{tab:turning-point-parameters} and Figure~\ref{fig:turning-point-explanation}), although $T_D$ has interesting nontrivial correlations with the frequency and temperature of turning point C. The other four parameters are well constrained. The input values of the parameters are within the posterior uncertainties.

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

In this section, we discuss complexities beyond those accounted for in the fits in Section~\ref{sec:results} for the sake of simplicity and how the method can be modified to deal with them in future work.

\subsection{Receiver time-dependence} \label{sec:receiver-time-dependence}

In this paper, we have assumed that the noise temperature and gain are functions of frequency; but, in reality, both the gain and noise temperature of the receiver will also vary with time.

The main source of variation in the receiver is the ambient temperature of the components in the analog radio-frequency signal chain.\footnote{Both gain and noise temperature are functions of the input antenna impedance, which can vary with ambient temperature due to thermal expansion of the antenna components. However, for most cases, this change is relatively negligible comparing to the actual thermal-coupled variations of the electronic components.}~To a large extent, the effects of ambient temperature on the gain and noise temperature can be characterized well in the laboratory. For example, both may rise uniformly in frequency with increasing temperature. In this case, the time dependence of both the gain and noise temperature can be encapsulated in a simple temperature parameter for each LST bin, which would be explored via MCMC when applying the techniques demonstrated in this work.\footnote{As the ambient temperature can be measured alongside the temperatures in the field and saved as metadata, these temperatures would have tight priors and should therefore be sampled efficiently by the MCMC algorithm.}

A common method of calibrating variations due to temperature drift is to model the receiver's frequency response within a certain operating temperature range. A full receiver circuit model can be carefully constructed and constrained by laboratory measurements to determine the instrument state when operating at different temperatures.\footnote{The gain and the noise temperature of a two-port network can be quantified by two sets of network parameters, namely the complex scattering $S$-parameters and noise parameters, which are constrained by fitting circuit models to laboratory measurements, similarly to the technique that was demonstrated in \cite{Nhan:19}.} In addition, the instrument can have internal sources, such as broad-band noise sources, that aid in calibration.~However, even after calibration, uncertainty in the model must be accounted for in order to produce rigorous constraints.~Therefore, in this case, we must form a model for $G_{\nu,t}/\overline{G}_{\nu,t}$ where $\overline{G}_{\nu,t}$ is the gain curve assumed for calibration.~Similarly a model for $O_{\nu,t}-\overline{O}_{\nu,t}$ must be formed to allow for noise temperature variations in a fit to calibrated data.

While the receiver gain was considered to be a deterministic quantity in the fits presented in this paper, there are also stochastic perturbations that generally introduce $1/f$ noise which in turn creates a noise floor, below which the data cannot be averaged down regardless of integration time. In order to reach the ppm-level dynamic range necessary to extract the global signal, the instrument must be dynamically stabilized through state-dependent corrections so that the noise floor is not reached until after 1000 hours of integration.

\subsection{Pass-band ripples}

While the gain and noise temperature models given in Section~\ref{sec:models} are good approximations to the gain and noise temperature when the receiver is designed to fit them, the in-band ripples will not exactly match any analytical model such as that given by the type-I Chebyshev filter transfer function, due to circuit components that vary from their nominal parameters within some manufacturing tolerance and exhibit parasitic behavior and impedance mismatches between the antenna and receiver. The gain model will need to be modified in order to prevent a modeling error from biasing the fitted signal. One way of including these non-ideal ripples is by forming both the gain and noise temperature models via training sets as performed in Section~\ref{sec:linear-models} instead of as an a priori known function like the Chebyshev filter gain model. To create such training sets, the same receiver circuit model mentioned in Section~\ref{sec:receiver-time-dependence} can be evaluated at a large number of physically reasonable variations in the characteristics of individual electronic components. CTP and DAPPER will adopt a similar scheme to create their receiver gain and noise temperature models.\footnote{The CTP and DAPPER receivers will be designed using a type-II (inverse) Chebyshev filter gain model, which has flatter spectral behavior in the pass-band.}

\subsection{Correlations between gain and noise temperature} \label{sec:nonlinear-offset-model}

In this paper, we have implicitly assumed that variations in the gain are independent from variations in the noise temperature; but, this is not true in reality. For instance, as the physical temperature of the receiver increases, the noise temperature and gain should both vary across the band.~Taking advantage of these correlations in the fit should improve it because it would force the model to assign a smaller probability to variations in the gain and noise temperature that do not correspond to each other. To do this, we could enforce that the gain and noise temperature models share some parameters.~For instance, the gain model and noise temperature model could be written ${\mathcal{M}}_{g_\nu}({\boldsymbol{\zeta}}_{G+o},{\boldsymbol{\zeta}}_G)$ and ${\mathcal{M}}_{o_\nu}({\boldsymbol{\zeta}}_{G+o},{\boldsymbol{\zeta}}_o)$, respectively, where ${\boldsymbol{\zeta}}_G$ (${\boldsymbol{\zeta}}_o$) represents the parameters that are only relevant to gain (noise temperature) variations and ${\boldsymbol{\zeta}}_{G+o}$ represents parameters (such as temperature) that cause correlated variations in the gain and noise temperature. For both uses of AMLP in this paper, we could then numerically explore ${\boldsymbol{\zeta}}_G$, ${\boldsymbol{\zeta}}_o$, and ${\boldsymbol{\zeta}}_{G+o}$ instead of marginalizing over any of them since the majority of the benefit of AMLP comes from marginalizing over beam-weighted foreground parameters.

\subsection{Polarization measurements}

In addition to laying out that using multiple correlated spectra (as done in this paper) vastly improves global signal fits, Paper III showed that posterior signal uncertainties can be significantly improved by including full-Stokes polarization measurements, which were not included in this paper for simplicity. To fit them, the observation equation (\ref{eq:observation-equation}) must be modified. For example, in the case of a dual-dipole antenna, it becomes
\begin{widetext}
\begin{equation}
  \underbrace{\begin{bmatrix} I \\ Q \\ U \\ V \end{bmatrix}}_{T} = \underbrace{\begin{bmatrix} \frac{1}{2}(|g_X|^2+|g_Y|^2) & \frac{1}{2}(|g_X|^2-|g_Y|^2) & 0 & 0 \\ \frac{1}{2}(|g_X|^2-|g_Y|^2) & \frac{1}{2}(|g_X|^2+|g_Y|^2) & 0 & 0 \\ 0 & 0 & {\text{Re}}(g_X^\ast g_Y) & -{\text{Im}}(g_X^\ast g_Y) \\ 0 & 0 & {\text{Im}}(g_X^\ast g_Y) & {\text{Re}}(g_X^\ast g_Y) \end{bmatrix}}_{G} \left( \underbrace{\begin{bmatrix} f_I \\ f_Q \\ f_U \\ f_V \end{bmatrix}}_{f} + \underbrace{\begin{bmatrix} s \\ 0 \\ 0 \\ 0 \end{bmatrix}}_{s} + \underbrace{\begin{bmatrix} o_X+o_Y \\ o_X-o_Y \\ 0 \\ 0 \end{bmatrix}}_{o} \right) + \underbrace{\begin{bmatrix} n_I \\ n_Q \\ n_U \\ n_V \end{bmatrix}}_{n},
\end{equation}
\end{widetext}
where $I$, $Q$, $U$, and $V$ are the measured Stokes parameters to be fit, $o_Z$ and $g_Z$ are the noise temperature and complex voltage gain of the $Z$ dipole feed (where $Z$ is either $X$ or $Y$), $f_P$ and $n_P$ are the beam-weighted foreground and Gaussian noise components in the Stokes parameter $P$, and the $\nu$ and $t$ dependence of all quantities involved has been left off for clarity. As shown by the brackets beneath the individual terms, this is Equation~\ref{eq:observation-equation} with $T$, $f$, $s$, $o$, and $n$ generalized to be vectors and $G$ generalized to be a matrix.~It is important to note that even under this generalization, $T$ is still linear in $f$, $s$, and $o$ at constant $G$; so, the methods of this paper still apply.

CTP and DAPPER are planned to have a four-channel correlation receiver that treats the $+X$, $-X$, $+Y$, and $-Y$ antennas as monopoles. This can be described by another straightforward generalization of the observation equation with the Stokes parameters replaced by the auto- and cross-correlations of the monopole voltage signals, each with its own complex voltage gain and noise. Under this formalism, the Stokes parameters can be formed from cross-correlations alone, meaning that the uncorrelated noise of the monopoles can be avoided. Sources of noise that are correlated between different channels of the receiver will still produce Stokes parameter profiles. But, the vast majority of the noise should be uncorrelated between the channels; so, this four channel correlation receiver can avoid most additive biases.

\section{Conclusions} \label{sec:conclusions}

This paper concludes the series on our data analysis pipeline for global 21-cm signal experiments, which was designed as a rigorous alternative to the most common methods used by existing global signal experiments, such as the use of polynomial models.~We have introduced several key ideas throughout the series:
\begin{itemize}
  \item The beam-weighted foreground must be fit down to the noise level by the chosen model. The best way to achieve this is to build a model specific to the given experimental situation (e.g. antenna design, pointing direction, location, etc.).~We do so through matrix decompositions of training sets such as singular value decomposition or principal component analysis.
  \item Signal parameter constraints can be explored rigorously and efficiently by analytically marginalizing over as many parameters as possible while numerically sampling only the essentially nonlinear parameters, i.e. those that leave the model linear if they are fixed.
  \item Uncertainties are vastly improved when including more than one spectrum of data when fitting. For example, modeling the correlations between foreground spectra at different times can decrease uncertainties from the $\sim1$ K level to the $\sim10$ mK level, required for detection of the expected signal, for total integration times in the hundreds of hours.~Including measurements of all four Stokes parameters decreases the uncertainties down nearer to the few mK noise level, allowing for precision cosmology.
  \item While the receiver systematic effects introduce nonlinearity to the model and new potential for confusion when extracting the global signal, they can be rigorously included in the pipeline using the AMLP methodology described in this work. When is well characterized via lab measurements and simulated models, the receiver does not significantly impact the precision of signal extraction.
\end{itemize}

In addition to the points above focused on the global signal, the analytical marginalization of linear parameters (AMLP) technique introduced in this paper can be used to explore any posterior distribution that has Gaussian conditional distributions more efficiently. AMLP is also included in \texttt{pylinex},\footnote{Download available at \url{https://bitbucket.org/ktausch/pylinex}.} the general, publicly available Python code that implements the fitting procedures laid out throughout the series \citep{pylinex:21}.

The pipeline is built into the design of the newest version of the Cosmic Twilight Polarimeter (CTP) and the proposed Dark Ages Polarimeter PathfindER (DAPPER), allowing for training sets to be developed through lab measurements (receiver), theory (signal), external observations (foreground), and simulations (antenna beam and receiver) before observations begin.

\acknowledgments{\noindent We thank Eric Switzer for the early suggestion of numerically exploring only signal parameters, an idea which has evolved into AMLP. D.R.~was supported in the early stages of this work by a NASA Postdoctoral Program Senior Fellowship at the NASA Ames Research Center, administered by the Universities Space Research Association under contract with the National Aeronautics and Space Administration (NASA). This work was also supported by NASA under award number NNA16BD14C for NASA Academic Mission Services.~B.N.~is a Jansky Fellow of the National Radio Astronomy Observatory.~This work is directly supported by the NASA Solar System Exploration Virtual Institute cooperative agreement 80ARC017M0006.}

\software{numpy, scipy, matplotlib, distpy+pylinex}

\newpage

\bibliographystyle{aasjournal}


\section{Introduction}
\label{sec1}

After isolation of a single layer of graphite (Graphene) \cite{novoselov2004electric}, obtained by mechanical exfoliation, there is a renewed interest in obtaining new 2D materials \cite{zhang2005experimental,balandin2008superior}. From a topological point of view, graphene layers can be considered as the basis to generate other structures, such as carbon nanotubes \cite{iijima1993single} and nanoscrolls \cite{bacon1960growth}. Carbon nanotubes are understood as graphene layers rolled up into cylinders, while nanoscrolls are graphene layers rolled up into a papyrus-like form. There are many theoretical and experimental works investigating nanoscrolls of different materials \cite{tomanek2002mesoscopic,lisa2003viculis,viculis2003chemical,braga2004structure,coluci2007prediction,shi2010tunable,uhm2020structural,perim2013controlled,junior2021reactive,tromer_PE,perim_Frontiers}.

In contrast to carbon nanotubes (CNTs), carbon nanoscrolls (CNSs) have open ends, which allows an easily radial expansion (tunable interlayer distances) explored in different applications \cite{zhang2010carbon,shi2010mechanics}. The most significant barrier in their usage has been the limited synthesis of high-quality structures. It leads to a controlling inability of some crucial aspects, among them, the number of scrolled layers \cite{liu2016layered,savin2016symmetric,savin2017graphene}. To overcome this barrier, new alternative routes in obtaining CNSs that allow higher control over the final structure to broaden their applications are desirable \cite{kuroda2011one,tojo2013controlled,shin2014ice}. CNSs exhibit distinct electronic and mechanical behaviors when compared with graphene and CNTs \cite{chen2007structural,song2013atomic}. The CNSs nucleus presents vibrational properties that can be useful in developing nano-actuators and electromechanical devices with excellent performance as cathodes in lithium-sulfur batteries \cite{shi2009gigahertz,shi2010translational,zuo2016novel}. It is worth mentioning that Density Functional Theory (DFT) calculations were used to address the synthetic growth concept of carbon-based materials \cite{freitas_JPCC}.

Recently, molecular dynamics simulations have been used to investigate other possible synthetic routes \cite{zhang2010carbon,perim2013controlled,zhang2016formation}. Results showed that CNT triggers the self-scrolling process of graphene nanoribbons. As a result, the layer spontaneously rolls itself around the CNT to lower the surface energy, a process governed by van der Waals interactions. The advantage of this method is that it imposes dry, non-chemical, and room-temperature conditions \cite{zhang2010carbon,perim2013controlled,zhang2016formation}. An experimental study was conducted by \textit{Xie et. al.} to obtain stable CNSs focused on offering higher control over the final product \cite{xie2009controlled}. In their method, graphite is mechanically exfoliated, deposited over SiO$_2$ substrate, and then immersed in a solution of water and isopropyl alcohol. After resting for a few minutes, the system spontaneously formed CNSs \cite{xie2009controlled}. Another possibility is to use functionalized structures such as graphene oxides and hydrogenated ones \cite{zuo2016novel,amadei_nanoscale}.

Herein, we used fully-atomistic reactive (ReaxFF) molecular dynamics (MD) simulations to investigate the effects of hydrogenation on the self-scrolling mechanism of graphene and graphane nanoribbons when triggered by CNTs. Our results show that it is possible to obtain stable CNT-wrapped CNSs from pure and hydrogenated graphane nanoribbons.

\section{Methodology}
\label{sec2}

In the present work, we used fully atomistic (ReaxFF) molecular dynamics simulations for C/H/O \cite{ashraf2017extension}, as implemented in LAMMPS code \cite{plimpton1995fast}. The initial configurations of our model systems are presented in Figure \ref{fig:complexes}, where we shown three atomic models identified as Complex-A (graphene/graphene/CNT, Complex-B (graphene/graphane/CNT), and Complex-C (graphane/graphane/ CNT). The two graphene-based nanoribbons have a length and width of 150 \r{A} and 21 \r{A}, respectively. The CNT used to trigger the self-scrolling mechanism is a conventional zigzag CNT $(0,22)$ with a diameter of 17.21 \r{A} and length of 100 \r{A}. The graphene nanoribbon and CNT have 1296 and 2112 carbon atoms, respectively. The initial distance between the nanoribbons and the CNT wall is 5.41 \r{A}.In the initial assembly, the nanoribbons are placed into contact with the nanotubes as shown in Figures \ref{fig:complexes}(d-f). The nanoribbons were placed on opposite sides regarding the CNT center.

\begin{figure*}[pos=t]
\centering
\includegraphics[width=0.8\linewidth]{fig1}
\caption{Model systems used to simulate the formation of CNT-wrapped CNSs. Illustration of the atomic configuration for (a) graphene and (b) graphane nanoribbons. The dashed lines highlight the structural units. (c) Diagrammatic representation for the initial arrangement of the system involving two nanoribbons and a CNT, where $d_1=150$ \AA~(nanoribbon length), $d_2=17.21$ \r{A} (CNT diameter), and $d_3=28.03$ \r{A} (initial distance between the nanoribbons). Initial configurations for the systems studied here: (d) graphene/graphene/CNT (Complex--A), (e) graphene/graphane/CNT (Complex--B), (f) and graphane/graphane/CNT (Complex--C). In the color scheme, the gray and yellow spheres represent carbon and hydrogen atoms, respectively.}
\label{fig:complexes}
\end{figure*}

The composite systems (Complexes --A, --B, and --C) were equilibrated/thermalized within an NVT ensemble at $300$, $600$, and $1000$ K using the Nos\'e-Hoover thermostat \cite{hoover1985canonical}. We considered 1000 K as the uppermost temperature limit, envisioned carbon-based optoelectronic applications. Scroll structures are resilient to high temperatures (up to 4000K), as recently demonstrated in references \cite{junior2021reactive,lee_Polymers}. The system equilibration was initially performed to eliminate any residual stress during 100 ps. The time step for the simulations was $0.1$ fs.  The MD snapshots and trajectories were obtained by using the free visualization and analysis program VMD \cite{humphrey1996vmd}.

In the present work, we performed 10 MD runs (with different seeds) for each temperature. The random seeds defined distinct initial sets of atomic velocities in the equilibration of the systems. All the MD runs yielded the same products for each complex. In this sense, the results reported in the manuscript are the most representative outcomes obtained for all the simulations.

\section{Results}

We begin our discussion by presenting the dynamical behavior of the self-scrolling process for all complexes studied here (see Figure \ref{fig:complexes}). In Figure \ref{fig:MD} we show a sequence of representative MD snapshots (in lines) for Complex-A, h-B, and -C from top to bottom, respectively, for simulations at $300$ K. In the very first stage of the simulation (about $2.5$ ps), the nanoribbons are in contact with the CNT surface, and the self-scrolling process takes place. This trend is realized in the simulations for all complexes, as shown in Figure \ref{fig:MD}. Such behavior also occurred for the self-scrolling mechanism of a single graphene \cite{xia2010fabrication} and graphane (see Supplementary Material) nanoribbons in contact with a CNT.

\begin{figure*}[pos=ht]
\centering
\includegraphics[width=\linewidth]{fig2}
\caption{Representative MD snapshots for the spontaneous self-scrolling process at 300 K, for the complexes presented in Figure \ref{fig:complexes}.}
\label{fig:MD}
\end{figure*}

At $12.5$ ps of simulation, the nanoribbons start to curl around the CNT. Subsequently, for Complexes A and B, the overlap between the nanoribbons occurs. The graphane nanoribbons in Complex-C do not form a scroll. Instead, it simply wraps the CNT by covering its surface and is temperature independent (for the temperature values considered here, up to $1000$ K), creating a lock-like mechanism. This behavior is attributed to stereo hydrogen-hydrogen interactions when the graphane layers move to create an opposing force to the sliding. This trend is similar to observed for other hydrocarbon structures, as the lander molecule, where this lock-like mechanism is present \cite{otero2004lock}.

Differences in the complexes dynamics become clear at $37.5$ ps. At this stage, in Complexes A and B, the nanoribbons are almost completely wrapped, as a consequence of the interlayer interaction presented by graphene and graphane \cite{savin2016symmetric,perez2015pi,savin_PRB}. Graphene/graphane bilayer (Complex-B) presents an interlayer interaction energy smaller than graphene/graphene one (Complex-A) \cite{umadevi_PCCP}. Even so, the self-scrolling process for both complexes takes the same time yielding a completely scrolled configuration at $50$ ps. For Complex-C, the graphane nanoribbons wrap the entire CNT wall at $130$ ps. In all cases, the nanoribbons showed discontinuous wrinkles and corrugations during the spontaneous self-scrolling process. Other reactive MD simulations performed within the scope of the AIREBO interatomic potential \cite{huang2015mechanical,song2013atomic} also pointed for such instabilities. The Supplemental Material includes the videos for the dynamics of the complexes A-C. It is worth mentioning that the degree of instability for the nanoribbons was more pronounced for temperatures higher than $300$ K.

The Supplementary Material presents the case in which a single graphane nanoribbon is interacting with a CNT. The main conclusions on the scrolling process hold, and the system evolved to a single graphane-based nanoscroll. Isolated graphane nanoscrolls were already investigated by molecular dynamics simulations \cite{savin_PRB}, and the results are consistent with the ones presented here.

It is important to stress that depending on the area of the nanoribbon, the most stable configuration is not planar but a scrolled one. This trend is related to the size of the overlapped region that enables van der Waals gains to overcome the elastic cost to scroll the structures \cite{braga2004structure}. Even for small graphene nanobelts, the self-scrolling can be a self-sustained process \cite{martins_nanotechlogy}.

\begin{figure*}[pos=t]
\centering
\includegraphics[width=\linewidth]{fig3}
\caption{Total energy time evolution as function of simulation time for the Complex-A at 10 K. The insets are MD snapshots extracted from the simulation trajectory at $t_1=2.5$, $t_2=20$, $t_3=40$, $t_4=65$, and $t_2=80$ ps, respectively.}
\label{fig:energy01}
\end{figure*}

In Figure \ref{fig:energy01} we present the time evolution for the total energy for the Complex-A at $10$ K at five different stages. During the dynamical process for the CNSs formation, there is a linear downhill trend that ends at $65$ ps, which represents a completely formed scroll. From that moment on, the total energy stabilizes, and a flat region persists (neglecting thermal fluctuations, scroll ``breathing'' - see video) until the end of the simulation. For all the temperature regimes simulated here ($300$, $600$, and $1000$ K), the time evolution of the total energy shows a similar trend. An overall picture of the CNSs formation is presented in Figure \ref{fig:energy02}. The existence of an energy plateau is the signature of the scroll formation. At high temperatures, the CNSs total energies are larger due to the higher degree of atomic thermal vibrations, which tends to difficult the spontaneous scrolling process. Even so, one can realize a well-defined plateau in all the cases, standing for the scroll's formation.

\begin{figure*}[pos=t]
\centering
\includegraphics[width=\linewidth]{fig4}
\caption{Total energy time evolution as a function of the simulation time for all complexes at 300, 600, and 1000 K. The red, blue, and green data points refer to the total energy values during the simulation time for the Complex-A, Complex-B, and Complex-C, respectively. The lines were obtained from the data interpolation.}
\label{fig:energy02}
\end{figure*}

\section{Conclusion}
\label{sec4}

Through fully-atomistic reactive (ReaxFF) MD simulations, we have studied the hydrogenation effects on the self-scrolling process of graphene and graphane (hydrogenated graphene) nanoribbons interacting with a carbon nanotube (CNT). We have considered different configurations: (1) graphene/graphene/CNT; (2) graphene/graphane/CNT, (3) graphane/graphane/CNT; and different temperatures, $300$, $600$, and $1000$ K. Our results show that the scroll formation is possible for configuration (1) and (2), for all temperatures. Configuration (3) does not yield scrolls. The scroll formation is related to the existence of an energy plateau in the total energy as a function of the simulation time. The approach used here to assemble the nanoribbons and CNT is prototypical and extensible to other species of 2D materials beyond graphene. The usage of nanotubes for triggering the scrolling processes of sheets and nanoribbons was investigated for carbon and boron nitrides systems \cite{li_PCCP,perim_CPC,siahlo2018structure}. The tribology of carbon-based nanostructures has been an object of research in the last years \cite{guo2005coupled,lebedeva2009dissipation}. In this sense, the present study can stimulate further investigations along these lines.

\section*{Acknowledgements}

The authors gratefully acknowledge the financial support from Brazilian research agencies CNPq, FAPESP, and FAP-DF. M.L.P.J gratefully acknowledges the financial support from CAPES grant 88882.383674/2019-01. J.M.S. and D.S.G. thanks the Center for Computing in Engineering and Sciences at Unicamp for financial support through the FAPESP/CEPID Grants \#2013/08293-7 and \#2018/11352-7. L.A.R.J acknowledges the financial support from a Brazilian Research Council FAP-DF and CNPq grants $00193.0000248/2019-32$ and $302236/2018-0$, respectively. L.A.R.J and J.M.S acknowledge CENAPAD-SP for providing the computational facilities. L.A.R.J. gratefully acknowledges the financial support from IFD/UnB (Edital $01/2020$) grant $23106.090790/2020-86$. The authors acknowledge the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer, which have contributed to the research results reported within this paper. URL: \url{http://sdumont.lncc.br}.

\printcredits
\bibliographystyle{unsrt}