\section{Introduction}

\noindent
Dark matter (DM) constitutes 84\% of the matter content in the universe \cite{aghanim2020planck} and plays an important role in the evolution of the early universe. It has so far eluded detection in all channels other than gravitational interactions. DM annihilation or decay could inject energy in the form of Standard Model particles, modifying the temperature and ionization of the intergalactic medium (IGM) and the anisotropies of the cosmic microwave background (CMB); studies of these observables have placed strong constraints on such energy injections (e.g. \cite{adams1998cosmic, chen2004particle, padmanabhan2005detecting, zhang2006impacts, slatyer2009cmb, kanzaki2010effects, slatyer2013energy, galli2013systematic,Slatyer:2015jla,Slatyer:2015kla,liu2016contributions, slatyer2017general, poulin2017cosmological, liu2021lyman}).

\texttt{DarkHistory}\xspace \cite{liu2020darkhistory} is a \texttt{Python} package developed to calculate the evolution of the IGM temperature and ionization in the early universe in the presence of such exotic energy injections. For an injected spectrum of Standard Model (SM) particles, it calculates the particle cascade by computing
(1) the production of photons and electrons/positrons by the decay of the originally injected SM particles;
(2) the subsequent secondary particle cascade and energy deposition arising from this exotic injection of photons/electrons/positrons, due to interaction with the IGM and the photon bath;
(3) modifications to the IGM's temperature and ionization from the secondary particles and their energy deposition, using a simple Three-Level Atom (TLA) model.

These calculations are carried out in redshift steps starting prior to recombination (at redshift $1+z=3000$ by default) and ending well after reionization near the present day (redshift $1+z=4$). In particular, the particle cascade in step (2) is evaluated using precomputed transfer functions, which are matrices that take an input spectrum and output the spectrum of secondary particles (for a given redshift step). \texttt{DarkHistory}\xspace includes the backreaction effects of changes to the ionization level of matter, which means the transfer functions themselves are functions of the gas ionization levels, as well as redshift. In the previous version of \texttt{DarkHistory}\xspace, this dependence is realized by interpolating tables of transfer function matrices on a grid of values for the hydrogen and helium ionization fractions, as well as a grid of redshift values.

\begin{figure}[ht]
\includegraphics[clip, width=0.4\textwidth]{model.pdf}
\caption{\textbf{Model schematics of a Dense Neural Network transfer function.} A DNN takes in (the logarithm of) the input and output photon/electron energy and necessary physical parameters including the redshift $z$, ionized hydrogen fraction $x_\text{HII}$, and singly ionized helium fraction $x_\text{HeII}$, and outputs (the logarithm of) the transfer function value $P$. The transfer function matrix acting on a given discretized spectrum is then obtained by evaluating the DNN on the given energy abscissa.}
\label{fig:model}
\end{figure}

At around $1.5$~Gb per table with 12 tables around this size, the transfer functions take up significant storage space as well as memory during runtime, since they are all loaded in a standard run. They will also be difficult to scale up to include additional parametric dependence, as the expected size scales exponentially with the number of added parameters. Intuitively, storing the transfer functions as tables is an over-representation of the information content within, since the transfer functions, despite not being smooth globally, can be divided into multiple regions that are relatively smooth, each of which could plausibly be fitted with analytical functions. However, in practice, finding such a solution is quite non-trivial, and even if a solution was found by \emph{ad hoc} methods, it would be quite specific to individual transfer functions and potentially difficult to maintain under changes to the code.
In this work, we present a general solution to these issues by replacing the transfer function tables with trained Dense Neural Networks.

Recently, machine learning and especially Neural Networks (NNs) has seen many applications in high energy physics and astrophysics \cite{lecun2015deep, bourilkov2019machine, radovic2018machine, albertsson2018machine}. NNs are often used to capture highly nonlinear and complex relations between inputs and outputs, and in particular, Dense Neural Networks (DNNs), also known as fully connected Neural Networks, are a class of general function approximators given sufficient neuron numbers \cite{Goodfellow-et-al-2016}.

\begin{figure*}[ht]
\includegraphics[clip, width=0.7\textwidth]{dhtfs.pdf}
\caption{\textbf{\texttt{DarkHistory}\xspace flow chart and transfer functions.} Starting with injected photons and electrons, this flow chart illustrates the work flow of \texttt{DarkHistory}\xspace in the evolution of each redshift step. The boxed quantities are particle spectra and the arrows indicate transfer functions, which act on the particle spectra. The orange transfer functions are generated or reconstructed using DNNs, while the blue ones only use tabulated data. (Modified from Fig.~1 in Ref.~\cite{liu2020darkhistory}.) }
\label{fig:dhtfs}
\end{figure*}

In the updated version of \texttt{DarkHistory}\xspace we present in this work, we use lightweight DNNs to store and automatically interpolate \texttt{DarkHistory}\xspace's transfer functions. A transfer function is essentially a multi-dimensional table, with 2 dimensions corresponding to the input and output particle energy, and the rest corresponding to physical parameters which in this case will be a subset of redshift $z$, ionized hydrogen fraction $x_\text{HII}\equiv n_\text{HII}/n_\text{H}$, and singly ionized helium fraction $x_\text{HeII}\equiv n_\text{HeII}/n_\text{H}$. (As a matter of convenience, we define the singly ionized helium fraction with the \emph{hydrogen} number density as the denominator, so that $x_\text{HeII}$ and $x_\text{HII}$ can be easily summed.) The DNN-based transfer functions are shown in schematic form in Fig.~\ref{fig:model}. We let a DNN take in all of the relevant parameters on equal footing and predict (the natural logarithm of) the transfer function value $P$.

Using DNNs as transfer functions has several benefits:

\begin{itemize}
    \item The DNNs we use are smaller in size compared to the stored transfer function tables, by a factor of $\sim$400.
    \item The computed matter temperature history and ionization history
    match those calculated using the transfer function tables to within a few percent relative difference (with sub-percent relative differences in regions when the species in question are more than 10\% ionized), while the spectral distortion due to upscattered CMB photons (see Sec.~\ref{sec:tf} for precise definition) matches to below the 10 percent level. These errors are small compared to current experimental uncertainties.
    \item We expect the DNNs to have improved scaling in size when additional parameters are added compared to the original tables. (Including a smooth dependence on an additional parameter might result in a $\order1$ increase in the DNN neuron number to reach similar accuracy due to increased information content. On the other hand, adding an additional dimension to a data table would increase its size by a multiplicative factor of the number of bins in the new parameter.)
    \item The DNNs automatically interpolate to any the physical parameter values and input/output particle energies within the trained range. This allows the use of flexible binning in \texttt{DarkHistory}\xspace, and will also allow us to perform interpolation on sparse training data \cite{duchi2011adaptive, duchi2013estimation}. The latter may become necessary in future extensions of \texttt{DarkHistory}\xspace, when probing dependence on an increasing number of physical parameters and generating dense grids of training data becomes computationally expensive.
    \item The DNNs predict transfer functions quickly, taking up a similar amount of time to the rest of the evolution routine for injected particles. Thus compared to retrieving tabular data from memory on a personal computer, the use of DNNs results in only a $\order1$ increase in total runtime.
    \item Open source building and training tools for NNs and especially simple architectures like DNNs are readily available. (In this work we use Tensorflow 2.0 \cite{tensorflow2015-whitepaper} with Keras \cite{chollet2015keras}.)

\end{itemize}

In Sec.~\ref{sec:tf}, we introduce \texttt{DarkHistory}\xspace and the roles of transfer functions. In Sec.~\ref{sec:nntf}, we detail the training and implementation of the DNN transfer functions. In Sec.~\ref{sec:performance} we present test runs and discuss the performance of DNN transfer functions compared to baseline \texttt{DarkHistory}\xspace. Finally, in Sec.~\ref{sec:conclusion} we summarize our results, briefly discuss other possible approaches, and outline some ideas for future expansion with DNN transfer functions in \texttt{DarkHistory}\xspace.

\section{Transfer functions in DarkHistory}
\label{sec:tf}

To better illustrate the role of transfer functions, we briefly introduce the procedure followed in \texttt{DarkHistory}\xspace and sketched in Fig.~\ref{fig:dhtfs}, which is modified from a flow chart in Ref.~\cite{liu2020darkhistory}. In Fig.~\ref{fig:dhtfs}, boxed quantities represent particle spectra, and arrows represent transfer functions, which are functions acting on spectra. \texttt{DarkHistory}\xspace stores the free streaming photon spectrum, the IGM temperature, and the IGM's ionized hydrogen fraction and singly ionized helium fraction at each redshift step (these quantities are assumed to be homogeneous). For each redshift step, \texttt{DarkHistory}\xspace:
\begin{enumerate}
    \item Converts injected SM particles at that redshift to injected photons $\mathbf{N}^\gamma_\text{inj}$ and electron/positrons $\mathbf{N}^e_\text{inj}$ (hereafter referred collectively as electrons). In \texttt{DarkHistory}\xspace, high energy ($>3$~keV) positrons are treated as electrons since their dominant energy-loss process, Inverse Compton Scattering (ICS) on the CMB, does not depend on the particle charge.  Lower-energy positrons are assumed to annihilate and the resulting photon spectrum is tracked; their kinetic energy is approximated as following the same pattern of energy deposition as that of the electrons (see Ref.~\cite{liu2020darkhistory} for a more in-depth discussion).
    \item Computes any injected electron spectrum's energy deposition into ionization, excitation, or heating by applying the transfer function $\overline{\mathbf{R}}_c$. Computes secondary photon and electron spectra produced from injected electrons due to ICS, positronium formation and decay, and atomic processes by applying the ICS transfer function $\overline{\mathbf{T}}_\text{ICS}$ and secondary electron transfer function $\overline{\mathbf{T}}_e$, and evaluating the spectrum of gamma rays produced from positron annihilation  $\overline{\mathbf{N}}^\gamma_\text{pos}$. These secondary photon spectra are added to the spectrum of photons $\mathbf{N}^\gamma$ propagated from the previous step, plus any injected photon spectrum.
    \item Computes the secondary particles and energy deposition produced by a propagating photon spectrum $\mathbf{N}^\gamma$, due to a variety of processes, including photon-photon scattering, Compton scattering, pair production, photoionization, and redshifting. The production of secondary electrons/positrons in the same redshift step, and their subsequent production of photons via ICS on the CMB or positron annihilation, are also included. The propagating photon spectrum for the next redshift step is obtained by applying the high energy photon transfer function $\overline{\mathbf{P}}_\gamma$ to $\mathbf{N}^\gamma$. The low energy photon spectrum, which stores photons below 3 keV that either photoionize within the redshift step or lie below 13.6 eV, is obtained by applying the low energy photon transfer function $\overline{\mathbf{D}}_\gamma$.
    The low energy electron spectrum, which stores electrons with kinetic energy below 3 keV where atomic cooling dominates ICS and is treated separately in the electron cooling module, is obtained by applying the low energy electron transfer function $\overline{\mathbf{D}}_e$. Finally, the photon's energy deposition into ionization, excitation, and heating is obtained by applying the high energy deposition transfer function $\overline{\mathbf{D}}_c^\text{high}$.
    \item Computes the change to IGM temperature and ionization by first calculating the energy deposition fraction $f_c$'s, from the low-energy electron/photon spectra and direct energy deposition by higher-energy particles, and then performing the TLA integration (see Ref.~\cite{liu2020darkhistory} for details).
\end{enumerate}

In general, transfer functions from an input spectrum to another spectrum (such as $\overline{\mathbf{T}}_\text{ICS}$, $\overline{\mathbf{P}}_\gamma$) take up much more space than those outputting an energy value (such as $\overline{\mathbf{D}}_c^\text{high}$) or those that are diagonal (such as $\overline{\mathbf{D}}_\gamma$). For this version of \texttt{DarkHistory}\xspace, we focus on replacing the largest spectral transfer functions with DNNs, but similar procedures can be applied to the lower-dimension transfer functions in the future. In the following we introduce the two major types of transfer functions we will replace in more detail.

\subsection{ICS transfer functions}

As discussed above, the ICS transfer functions describe the spectrum of scattered photons from the complete cooling of injected electrons. Unlike the transfer functions applied to the photon spectrum, the ICS transfer functions $\overline{\mathbf{T}}_\text{ICS}$, $\overline{\mathbf{T}}_e$ and $\overline{\mathbf{R}}_c$ are not directly interpolated from tables, but reconstructed from reference ICS transfer function tables. Sec.~III.D and Appendix A of Ref.~\cite{liu2020darkhistory} describe in detail how this is achieved, and we only provide a simple summary here: since an electron quickly deposits all of its energy in one of our redshift steps, in order to obtain the total secondary spectrum or energy deposition by an electron, we need to consider multiple interaction events (via ICS and atomic processes). This can be done recursively: one can reconstruct the full ICS secondary spectrum and energy deposition for an electron of energy $E$ knowing the same information for all electrons with $E'<E$. With discretized energy abscissa, the full ICS-induced photon spectrum and energy output can be solved recursively from the lowest energy bin. As a result, at each redshift, one can solve for the full electron ICS transfer functions using transfer functions describing a single ICS scattering event (as well as functions describing the interaction rates due to atomic processes).

The transfer functions for a single ICS scattering event on the CMB have simple (approximate) scaling relations with respect to the temperature of the CMB $T$, as described in detail in Appendix A of Ref.~\cite{liu2020darkhistory}. As such one can derive ICS transfer functions at different redshifts from a single transfer function at a fixed redshift, assuming the CMB is the dominant radiation background (in \texttt{DarkHistory}\xspace, $1+z=400$ is used). These reference transfer functions are interpolated from tables \texttt{ics\_thomson}, \texttt{ics\_rel}, and \texttt{ics\_engloss}, corresponding respectively to transfer functions for the secondary photon spectra of nonrelativistic electrons and relativistic electrons, and the relativistic electron energy loss in a single ICS scattering event on the CMB. It is these tables that we will fit with DNNs.

\subsection{Photon transfer function}
\label{subsec:photontf}

The high energy photon transfer function $\overline{\mathbf{P}}_\gamma$, low energy photon transfer function $\overline{\mathbf{D}}_\gamma$, and low energy electron transfer function $\overline{\mathbf{D}}_e$ are interpolated from corresponding tables \texttt{highengphot}, \texttt{lowengphot}, and \texttt{lowengelec}. They in general depend on the CMB temperature through redshift and matter ionization levels (ionized hydrogen fraction and singly ionized helium fraction). For each combination of these physical parameters, \texttt{lowengphot} can be represented as 1-D arrays (with the one dimension being input/output energy) and is a factor of 500 smaller than the other transfer functions, so it is at present not replaced with a DNN. We also found that the numerical calculation of Compton scattering used in the previous version of \texttt{DarkHistory}\xspace was inaccurate in some parts of parameter space (in particular populating kinematically forbidden regions), and so we have updated the relevant tabular transfer functions to ensure sufficient accuracy.\footnote{This issue primarily affected secondary electrons in the 10 eV to 3 keV range produced by Compton scattering of photons with energies between 100 eV and 10s of keV; the effects on observable quantities were very small for all cases we checked. The updated tables on \href{https://doi.org/10.5281/zenodo.6819281}{Zenodo} incorporate this correction.} We now describe some special features of the photon transfer functions:

\begin{figure*}[ht]
\includegraphics[clip, width=\textwidth]{comptf_hep.pdf}
\caption{\textbf{Comparison of DNN and tabulated transfer functions.} The first two panels show the high energy photon transfer functions generated at redshift $1+z=300$, ionized hydrogen fraction $x_\text{HII}=0.6$ (and singly ionized helium fraction $x_\text{HeII}=0.0$). The transfer function value in the region encircled by the black dashed line is negative and the $\log_{10}$ of the absolute value is shown. The third panel shows the error in $\log_{10}{\mathbf{\overline{P}}_\gamma}$ (defined as $\Delta \log_{10} {\mathbf{\overline{P}}_\gamma} \equiv\log_{10}|{\mathbf{\overline{P}}_\gamma}_\text{DNN}|-\log_{10}|{\mathbf{\overline{P}}_\gamma}_\text{table}|$), which is concentrated on the physical features of the transfer functions. The fourth panel shows the relative error in the energy transfer function (see more details in Sec.~\ref{subsec:tfvaluepred}). The errors are concentrated on the diagonal that corresponds to photon propagation and redshifting.}
\label{fig:comptf_hep}
\end{figure*}

\paragraph{Redshift regimes and matter ionization dependence.}
All photon transfer functions depend on the redshift, as the photon cooling processes involve interactions with the redshift-dependent photon background and/or interstellar medium.
For late redshifts ($z<40$) encompassing the epoch of reionization, the transfer functions are allowed to vary with both the ionized hydrogen fraction $x_\text{HII}$ and the singly ionized helium fraction $x_\text{HeII}$, which can be altered by exotic energy injections.
For redshifts between helium recombination and reionization ($40<z<1600$), the ionized helium fraction can be safely approximated as zero \cite{liu2020darkhistory}; the transfer functions are pre-computed assuming no helium ionization, but can depend on the hydrogen ionized fraction.
Before helium recombination ($z>1600$), exotic energy injections consistent with current experimental bounds have little impact on the thermal equilibrium determining hydrogen and helium ionization levels \cite{liu2020darkhistory}, so \texttt{DarkHistory}\xspace uses \texttt{RECFAST} \cite{wong2008well} ionization fractions as a baseline to pre-compute the transfer functions, which only depends on redshift.

For the DNN implementation, the flexibility of the network allows one DNN to be trained on the entire redshift range $4<z<3000$ with a learned dependence on the ionization levels, for each transfer function. However, using different networks for different redshift regimes gives slightly better accuracy, and the latter approach is chosen in this version of \texttt{DarkHistory}\xspace.

\paragraph{Redshift step coarsening and energy conservation.}
In the previous version of \texttt{DarkHistory}\xspace, photon transfer functions are computed with redshift step $\Delta\log(1+z)=0.001$. One can choose to increase the (log) step size to multiples of 0.001 to speed up computation. To combine multiple redshift steps, \texttt{DarkHistory}\xspace pre-composes multiple propagation transfer functions $\overline{\mathbf{P}}_\gamma$ and applies them appropriately to the deposition transfer functions $\overline{\mathbf{D}}_\gamma$, $\overline{\mathbf{D}}_e$, and $\overline{\mathbf{D}}_c^\text{high}$. Since the DNN implementation of transfer functions introduces a small amount of error that can accumulate over many redshift steps, in order to increase numerical stability,
we train the DNNs to reproduce the pre-composed transfer functions, and require the use of a fixed log redshift step of $\Delta\log(1+z)=0.012$.
\footnote{It was shown in \cite{liu2020darkhistory} that an increased redshift step size at this level will lead to errors smaller than $10^{-3}$, and so will not be the primary contributor to numerical error in our context. Note also that using a $n$-fold coarsened log redshift step is not equivalent to using every $n^\text{th}$ entry in the transfer function table. The redshift abscissa used for the table are not directly related to (and in general are coarser than)
the redshift steps used during a run; the table is interpolated to obtain the transfer functions at each redshift step.} To further decrease numerical error, we store the total fraction of the injected energy entering each type of secondary particle spectrum, for an injected photon of any given energy, and use these data (which are of $\sim1\%$ the size of the transfer functions) to ensure energy conservation while accounting for energy loss to redshifting. We discuss the precise procedure of transfer function pre-composing and imposing energy conservation in Appendix~\ref{appd:coarsening}.

\begin{figure*}[ht]
\includegraphics[clip, width=\textwidth]{tfinterp.pdf}
\caption{\textbf{The DNN high energy photon transfer function network's fitting and interpolation of the table transfer function values.} The panels show the DNN transfer function's (redshift regime $4<z<40$) dependence on each of the five input variables, in a limited region around a central value indicated on the plot. The shaded region is where both the table value and the DNN prediction are negative (and the log of the absolute value is shown). A fairly consistent relative error below $10\%$ is observed.}
\label{fig:interp}
\end{figure*}

\section{Transfer Functions from Neural Networks}
\label{sec:nntf}

As described earlier and as indicated in Fig.~\ref{fig:model}, we update the largest transfer function tables (\texttt{ics\_thomson}, \texttt{ics\_rel}, \texttt{ics\_engloss}, \texttt{highengphot}, and \texttt{lowengelec}) to DNN networks that take in input and output particle energies, redshift, and possibly (depending on redshift) the ionized hydrogen fraction and singly ionized helium fraction, to produce the transfer function value $P$. In general, $P$ can vary greatly across many orders of magnitude. (For example, after recombination the probability for a 10 keV photon to free stream, losing energy only through redshifting, is substantial, while the chance of it directly producing secondary photons of 1 keV is very close to 0, since this outcome is not kinematically allowed in a single Compton scattering event, nor is the scattered electron produced by Compton scattering or photoionization able to up-scatter other photons to 1 keV.) As such, we train the networks to output the natural logarithm of the transfer function value $\log P$.
Similarly, since the input/output energy abscissa and our redshift steps are also binned in log space by default, the networks also take the log value of these as inputs. The ionized hydrogen fraction and singly ionized helium fraction are linearly scaled to match the spread of other parameters before being fed into the DNN.

For high energy photons, the redshift-coarsened transfer functions (see Sec.~\ref{subsec:photontf}) are fitted. Note also that the transfer function value $P$ can be negative since by convention the CMB spectrum is subtracted from the transfer function and the negative values (and the positive values at higher energies) represent CMB photons being upscattered. In this case $\log|P|$ is predicted by the DNNs, and the negative value region is recovered in a post-processing step that identifies local minima of $\log|P|$ (which takes up an negligible time compared to the rest of the \texttt{DarkHistory}\xspace routine). Additionally, the transfer functions values near the diagonal in the input/output energy dimensions (corresponding to the free-streaming photons) are adjusted to enforce energy conservation up to redshifting. For details please refer to Appendix \ref{appd:coarsening}.

After adjusting hyperparameters,
we find that for all transfer functions in question, it is sufficient to use DNNs with 7 hidden layers with 400 neurons each, making the number of parameters per DNN $\sim 400^2 \times (7 - 1)\sim9.7\times10^5$,
and $2.9\times10^6$ for each transfer function built from 3 such networks.

Training is done with TensorFlow 2.0~\cite{tensorflow2015-whitepaper} and Keras~\cite{chollet2015keras}; Adagrad~\cite{lydia2019adagrad} is used as the optimizer, with the mean squared error of $\log P$ as the loss function. Each DNN is trained on 2 V100 GPUs for $\order{10}$ hours or equivalent. For each epoch, training data is generated by interpolating the multi-dimensional transfer function table on uniformly random sampled inputs. Since training data is not reused across epoch, there is no concern of overfitting to a fixed subset of full data set. Trainings are terminated with the evaluation after each iteration stops improving significantly. To check for systematic offsets of the table transfer functions and the NNs, we trained multiple NNs with random initial values and randomly sampled training data. We found no obvious systematic offsets between the ensemble of NNs and the tables; e.g. at any given point the different NNs both underpredicted and overpredicted the table data.

The codes related to the DNN transfer functions are stored in the \texttt{nntf} module under \texttt{DarkHistory}\xspace. A new example file ``Example 12: Using Neural Network transfer functions.ipynb"\href{https://github.com/hongwanliu/DarkHistory/blob/nntf_alt_hep_clean/examples/Example_12_Using_Neural_Network_Transfer_Functions.ipynb}{\faGithub} is provided to demonstrate using the DNN transfer functions and comparison with the baseline \texttt{DarkHistory}\xspace (which will be available if the appropriate data tables are present).

\section{Performance}
\label{sec:performance}

In this section, we describe the accuracy and speed of generation of the DNN transfer functions, as well as the accuracy of full runs over a range of simple DM injection scenarios using the DNN transfer functions.

\subsection{Transfer function value prediction}
\label{subsec:tfvaluepred}

In Fig.~\ref{fig:comptf_hep}, a high energy photon transfer function generated by a DNN is compared against one interpolated from tables. As one can see, the errors in the raw output i.e. logarithm value of the transfer functions are concentrated near the distinct physical features in the transfer function, such as at the output photon energy of $\sim0.5$~MeV corresponding to positronium decay. The absolute values in logarithm errors $\Delta\log_{10} |P|$ can be interpreted as relative errors in $|P|$ (up to a $\ln 10$ factor). In this particular slice through the high energy photon transfer function, the average $\Delta\log_{10}|P|$ when $|P|>10^{-20}$ is 0.017, corresponding to a relative error of $\sim4$\%. (The range of $\log_{10}|P|$ is about $\sim-45$ to 6). Overall $\Delta\log_{10} |P|$ is comparable to this value, for all DNN transfer functions. A summary of the errors can be found in Tab.~\ref{tab:tf_err_time}.

\begin{figure*}[ht]
\includegraphics[clip, width=\textwidth]{Txx_errors_decay_elec.pdf}
\caption{\textbf{Example evolution of matter temperature and ionization levels.} The evolution of matter temperature and ionization levels, in an example scenario where 0.1~GeV DM is decaying into $e^+e^-$ pairs, as predicted using the original tabulation ({\it black lines}) and using the DNN transfer functions ({\it red dashed lines}). The blue line shows the difference between the two, and the gray dashed line shows the standard evolution without any external injections.}
\label{fig:TmxHxHe}
\end{figure*}

To see the impact of these errors in a \texttt{DarkHistory}\xspace evolution run, it is also useful to look at errors in energy (per bin) transition rates, besides the particle number transition rates. The energy transfer function de-emphasize errors at low energies where many particles can be produced with only a small fraction of the original particle's energy, and thus such errors have a small effect on heating and ionization. Since the standard photon transfer functions are maps between particle number spectra $N_i$, the transfer function values have the physical meaning of number transition rates. The particle energy spectrum is related to the number spectrum by
\begin{equation}
    E_i = \mathcal E_i N_i=\mathcal E_i^2\left(\frac{\text{d} N}{\text{d}\mathcal E}\right)_i,
\end{equation}
where $N_i$ is the number of particles in the $i$-th bin, and $\mathcal E_i$ its central energy. (Note that the energy bins are log-spaced.) For a particle number transfer function $P$ with
\begin{equation}
    N^\text{out}_i = \sum_j N^\text{in}_j P_{ji},
\end{equation}
the corresponding energy transfer function $P^E$ is defined by
\begin{equation}
    E^\text{out}_i = \sum_j E^\text{in}_j P_{ji}^E
    \implies P_{ji}^E = P_{ji}\mathcal E_i/\mathcal E_j.
\end{equation}
In the last panel of Fig.~\ref{fig:comptf_hep}, we show the relative error in the high energy photon transfer function. As expected, the errors are concentrated on the highest output energy for a given input energy. The relative errors are generally sub-percent.

Fig.~\ref{fig:interp} shows that the DNN networks interpolate sensibly between the fixed abscissa values in the transfer function tables, for the high energy photon transfer function in the lowest redshift regime ($4<z<40$), as an example. The errors in the transfer function interpolation are consistent with average values shown in Tab.~\ref{tab:tf_err_time}.

The time it takes to generate a photon transfer function is about or less than 1 second on an 8-CPU personal computer. ICS transfer functions that are only generated once per run takes only slightly longer at $\lesssim3$~s. The accuracy and prediction time for all other transfer functions can be seen in Tab.~\ref{tab:tf_err_time}.

\begin{figure*}[ht]
\includegraphics[clip, width=0.8\textwidth]{scan_errors_decay_elec.pdf}
\caption{\textbf{Relative errors in temperature, ionization levels, and low energy photon distortion across a range of DM mass in a scenario with DM decaying to $e^+e^-$.} The left panel shows the maximum relative error at any point in the evolution in $4<z<3000$ of matter temperature, ionized hydrogen fraction, and singly ionized helium fraction across a range of DM exotic electron injection energies*. (*Due to its very small absolute value, the relative error for singly ionized helium fraction $x_\text{HeII}$ when $n_\text{HeII}/n_\text{He}<10^{-3}$ is not included in this plot. Since $n_\text{HeII}/n_\text{He}$ changes rapidly between $\order{1}$ and $<10^{-3}$, we are essentially only counting its relative error when its value is order unity.) The relative errors are generally below 5\%. The right panel shows the maximum error over the maximum value in the low energy photon spectral distortion at $z=0$. The spectral distortion errors at the two dashed-line values of $m_\text{DM}=0.1$~GeV and $10$~GeV are shown in Fig.~\ref{fig:spec}. }
\label{fig:scan_errors_decay_elec}
\end{figure*}

\begin{figure*}[ht]
\includegraphics[clip, width=0.8\textwidth]{spec_decay_elec_2.pdf}
\caption{\textbf{Example partial low-energy photon spectral distortion in the present day.} The two panels show two examples of low energy photon spectral distortion outputs (see definition in Sec.~\ref{sec:tf}) from two different runs: 0.1~GeV and 10~GeV DM decaying to $e^-e^+$ pairs. They represent the extremes of large and small relative errors for runs over the full range of DM masses we consider, as shown in Fig.~\ref{fig:scan_errors_decay_elec}. The photon number density (per bin) is normalized against the baryon number density.
The black lines represent outputs generated using tabulated transfer functions while the red dashed lines represent that using the DNN transfer functions. The blue line shows the different between the two. Note that the relatively large errors shown in the right panel are in part due to the error in the location of the zero.}
\label{fig:spec}
\end{figure*}

\subsection{Performance over a range of scenarios}

Fig.~\ref{fig:TmxHxHe} shows the evolution of integrated variables in one particular setting: 0.1~GeV DM particles decaying into electron-positron pairs (this is the same as the example used in Fig. 4 of Ref.~\cite{liu2020darkhistory}). In this example, the error introduced by the DNN in the matter temperature history and hydrogen ionization history is consistently sub-percent, while the error in singly ionized helium fraction is sub-percent when $n_\text{HeII}/n_\text{He}>10^{-3}$. For a scenario where DM decays into electron-positron pairs, over a range of DM rest masses, the maximum relative error for temperature and ionized hydrogen fraction over the entire redshift range is consistently below 2\%, as shown in Fig.~\ref{fig:scan_errors_decay_elec}. Taking into account injection scenarios with DM decaying to photons, and also undergoing $s$-wave annihilation into $e^+e^-$ or photons, the relative error is always below 8\%. \texttt{DarkHistory}\xspace also computes the partial photon spectral distortion from high energy photon and electron processes, mostly from ICS of electrons or positrons on the CMB. This distortion is stored in the low energy photon spectrum. Note that spectral distortions arising from atomic transitions, due to photons and electrons below 3~keV, are not included in this spectrum (which is why we label it as ``partial'' or ``incomplete''). The DNN transfer functions introduce a small amount of error in this spectral distortion, as shown in Fig.~\ref{fig:spec} and Fig.~\ref{fig:scan_errors_decay_elec}. While the shape of the spectral distortion is generally correct, the error in the location of the distortion zero can cause errors with a magnitude up to 10\% of the distortion magnitude due to the errors in the distortion zero location. In a future update of \texttt{DarkHistory}\xspace, we are anticipating an update to include the correct treatment of the complete photon spectral distortion including contributions from atomic transitions. A small photon spectral error would allow the DNN functions to be used simultaneously with these updates.

In Appendix \ref{appd:errors}, we include the errors in some other exotic injection scenarios, including DM decaying to photon pairs, and DM annihilating to photon or $e^+e^-$ pairs. Although these examples cover only a few simple injection scenarios, they can serve to test the whole range of \texttt{DarkHistory}\xspace's dependence on transfer functions, since all exotic energy injections are converted to either $e^+e^-$ or photon injections. We expect exotic energy injections with a more complicated injection spectrum (such as e.g. annihilation to quarks) to have similar amounts of error associated with using the DNN transfer function.

\begin{figure*}[ht]
\includegraphics[clip, width=0.8\textwidth]{scan_errors_decay_phot.pdf}
\caption{\textbf{Relative errors in physical quantities across a range of DM mass, in a scenario with DM decaying to a pair of photons.}
The panels have similar construction to Fig.~\ref{fig:scan_errors_decay_elec}. Note that the matter temperature and ionization fractions have max relative errors below 2\%, and the errors in the low-energy photon spectral distortion are also below 10\%.}
\label{fig:scan_errors_decay_phot}
\end{figure*}

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

In this work, we have made use of simple Dense Neural Networks to approximate complex and multi-dimensional transfer functions in \texttt{DarkHistory}\xspace to reduce storage and memory usage, as well as to enable the possibility of adding more parameter dependence to the transfer function. The DNN transfer functions achieve good accuracy in computing the evolution history of matter temperature and ionization, as well as the partial CMB spectral distortion evaluated by the current version of \texttt{DarkHistory}\xspace; typical errors are at the few percent level, and are comparable to or smaller than estimates of systematic uncertainties in previous studies of constraints on energy injection \cite{galli2013systematic,Slatyer:2015jla,liu2020darkhistory}.

The DNN-based functionality is available in the \texttt{DarkHistory}\xspace Github repository at \url{https://github.com/hongwanliu/DarkHistory}, and the necessary data files (one can choose to download the large tables, or DNN and auxiliary files, or both) are hosted on \href{https://doi.org/10.5281/zenodo.6819281}{Zenodo} (see the Github repository for details).

While the use of DNNs offers one solution to this challenge, there may well be other viable solutions. The information in the DNNs still seems likely to be an over-representation of the piecewise-smooth transfer functions. We have briefly explored some alternative methods, including fitting to conventional functions directly, and with the assistance of symbolic regression techniques \cite{schmidt2009distilling}; however, DNN networks stand out as the best solution (so far) in terms of fitting accuracy and ease of implementation. There is work ongoing to expand the capabilities of \texttt{DarkHistory}\xspace, and we look forward to exploring NN-based and alternative techniques in this context.

\section*{Appendix}

\subsection{Redshift step coarsening and energy conservation}
\label{appd:coarsening}

In this Appendix, we describe how the photon transfer function change with coarsened redshift step, (expanding on Sec.~III.E.3 of Ref.~\cite{liu2020darkhistory}), and how energy conservation while correctly accounting for photon energy loss due to redshift is implemented.

\subsubsection{Transfer functions without coarsening}
In \texttt{DarkHistory}\xspace's \texttt{main.evolve} function, evolution is discretized into log-normal redshift steps with fixed spacing $d$ (\texttt{dlnz} in code) where the next redshift $z'$ is expressed in terms of $z$ such that:
\begin{equation}
    \log(1+z')=\log(1+z)-d.
\end{equation}

Following \texttt{DarkHistory}\xspace's flow described in Fig.~\ref{fig:dhtfs}, to obtain the propagating photon spectrum ${\mathbf{N}_\text{prop}^\gamma}$, low energy photon spectrum ${\mathbf{N}_\text{low}^\gamma}$, low energy electron spectrum ${\mathbf{N}_\text{low}^e}$ and energy deposition array ${E^\text{high}_c}$, we retrieve the corresponding transfer functions ${\mathbf{\overline{P}}_\gamma}$, ${\mathbf{\overline{D}}_\gamma}$, ${\mathbf{\overline{D}}_e}$, and ${\mathbf{\overline{D}}^\text{high}_c}$ at a middle redshift $z_\text{mid}$ given by
\begin{equation}
    \log(1+z_\text{mid})=\log(1+z)-d/2.
\end{equation}
The corresponding transfer functions ${\mathbf{\overline{P}}_\gamma}$, ${\mathbf{\overline{D}}_\gamma}$, and ${\mathbf{\overline{D}}_e}$ all take in the propagating photon spectrum ${\mathbf{N}_\text{prop}^\gamma}$ at redshift $z$ and produce secondary spectra at $z'$. (${\mathbf{\overline{D}}^\text{high}_c}$ produces energy deposition in this redshift step.)
They are implemented as:
\begin{equation}
\label{eq:applytf}
\begin{aligned}
    {\mathbf{N}_\text{prop}^\gamma}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\cdot{\mathbf{\overline{P}}_\gamma}(z_\text{mid})=\sum_j{\mathbf{N}_\text{prop}^\gamma}_j{\mathbf{\overline{P}}_\gamma}_{ji} \\
    {\mathbf{N}_\text{low}^\gamma}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\cdot{\mathbf{\overline{D}}_\gamma}(z_\text{mid}) \\
    {\mathbf{N}_\text{low}^e}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\cdot{\mathbf{\overline{D}}_e}(z_\text{mid}) \\
    {E^\text{high}_c}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\cdot{\mathbf{\overline{D}}^\text{high}_c}(z_\text{mid})=\sum_j{\mathbf{N}_\text{prop}^\gamma}_j{\mathbf{\overline{D}}^\text{high}_c}_j,
\end{aligned}
\end{equation}
where $i$ and $j$ are indices of discretized energy abscissa.

\begin{figure*}[ht]
\includegraphics[clip, width=0.8\textwidth]{scan_errors_swave_elec.pdf}
\caption{\textbf{Relative errors in physical quantities across a range of DM mass in a scenario where DM undergoes $s$-wave annihilation to a electron/positron pair.}
The panels have similar construction to Fig.~\ref{fig:scan_errors_decay_elec}.
}
\label{fig:scan_errors_swave_elec}
\end{figure*}

\begin{figure*}[ht]
\includegraphics[clip, width=0.8\textwidth]{scan_errors_swave_phot.pdf}
\caption{\textbf{Relative errors in physical quantities across a range of DM mass in a scenario where DM undergoes $s$-wave annihilation to a photon pair.}
The panels have similar construction to Fig.~\ref{fig:scan_errors_decay_elec}.
}
\label{fig:scan_errors_swave_phot}
\end{figure*}

Energy conservation can be enforced straightforwardly: for an injection in any photon energy bin with central energy ${\mathcal{E}_i}$, the total output energy on the right hand side of the above equations should add up to ${\mathcal{E}_i}$ minus the loss to redshifting of the propagating photons. Photons below a certain energy $E_\text{relevant}$ do not contribute to redshift energy loss because they either dump all of their energies efficiently within one redshift step or free stream and no longer interact, in which case they are stored as an array history of low energy photons ${\mathbf{N}_\text{low}^\gamma}(z')$) and not immediately redshifted.

For the propagating photons that are redshifted by ${\mathbf{\overline{P}}_\gamma}$, the energy lost is approximately:
\begin{equation}
\begin{aligned}
    |E_{\text{redshift},i}|=\left(1-\frac{1+z'}{1+z}\right){\mathcal{E}_i}&=\frac{z-z'}{1+z}\,{\mathcal{E}_i}\approx d\,{\mathcal{E}_i}\\
    &\text{for }{\mathcal{E}_i}>E_\text{relevant}.
\end{aligned}
\end{equation}
Let the energy abscissa (log-central energies of each bin) for photons and electron be ${\mathcal{E}^\gamma}_i$ and ${\mathcal{E}^e}_i$ respectively. We can express the above equation as
\begin{equation}
    |E_{\text{redshift},i}|=d_i{\mathcal{E}_i}
\end{equation}
where $d_i=d$ when the photon with energy ${\mathcal{E}_i}$ should be redshifted and $d_i=0$ otherwise. (This renders $d_i$ a function of redshift and hydrogen and helium ionization levels in general.)

Then the energy conservation constraint can be written as
\begin{equation}
    {\mathcal{E}_i}=\sum_j\left({\mathbf{\overline{P}}_\gamma}_{ij}{\mathcal{E}^\gamma}_j+{\mathbf{\overline{D}}_\gamma}_{ij}{\mathcal{E}^\gamma}_j+{\mathbf{\overline{D}}_e}_{ij}{\mathcal{E}^e}_j+{\mathbf{\overline{D}}^\text{high}_c}_i+d_i{\mathcal{E}_i}\right)
\end{equation}
This relation is imposed by shifting the near diagonal (propagating) part of the high energy photon transfer function ${\mathbf{\overline{P}}_\gamma}$. Any energy non-conservation due to numerical errors or errors from approximating the transfer function as DNNs can be absorbed into this shift.

\begin{table*}[ht]
\centering
\begin{tabular} {|c||c|c||c|c|}
\hline\hline
DNN transfer function & $|\Delta\log_{10}|P||$ & Prediction time & Table size & NN size \\
\hline
high energy photon (compounded) regime 0 & 0.029 & 1.01 s & & \\
high energy photon (compounded) regime 1 & 0.043 & 1.04 s & 4.7Gb & 11.4Mb \\
high energy photon (compounded) regime 2 & 0.016 & 1.02 s & & \\ \cline{4-5}
high energy photon (propagator) regime 0 & 0.029 & 1.02 s & & \\
high energy photon (propagator) regime 1 & 0.043 & 1.05 s & 4.7Gb & 11.4Mb \\
high energy photon (propagator) regime 2 & 0.062 & 1.07 s & & \\ \cline{4-5}
low energy electron regime 0 & 0.047 & 0.382 s & & \\
low energy electron regime 1 & 0.049 & 0.376 s & 4.7Gb & 11.4Mb \\
low energy electron regime 2 & 0.040 & 0.370 s & & \\ \cline{4-5}
ICS Thomson              & 0.00199 & 2.82 s & 0.93Gb & 3.8Mb \\
ICS relativistic         & 0.00410 & 2.12 s & 0.93Gb & 3.8Mb \\
ICS electron energy loss & 0.00250 & 2.30 s & 0.93Gb & 3.8Mb \\
\hline
\end{tabular}
\caption{\textbf{Error, prediction time, and size of DNN transfer functions used in \texttt{DarkHistory}.} The error of the DNN transfer functions in comparison to the tabular transfer functions is defined as $|\Delta\log_{10}|P||\equiv|\log_{10}|P_\text{DNN}|-\log_{10}|P_\text{table}||$ when $|P|>10^{-20}$. Both the errors and prediction times are evaluated from random draws in the relevant domains for each DNN until the values stabilize. The prediction times are evaluated on a 8-CPU personal computer.
All DNN transfer functions use 6 hidden layers of 400 neurons each, but the evaluation time differs due to the length of the input array corresponding to the size of the evaluated transfer function matrix ($\sim10^5$ for high energy photon, $\sim4\times10^4$ for low energy electron, $\sim2.5\times10^5$ for ICS).}
\label{tab:tf_err_time}
\end{table*}

\subsubsection{Coarsening}
\texttt{DarkHistory}\xspace can enlarge the redshift step $d$ to an multiple of a preset step $d_0$. Let the coarsening multiple be $c$, then the next redshift step $z'$ to $z$ is such that
\begin{equation}
    \log(1+z')=\log(1+z)-c\cdot d_0.
\end{equation}
The photon and electron transfer functions are built with log redshift step $d_0$, so we have to reconstruct the transfer functions for $d$. We first obtain single-step transfer functions (and $d_i$) with current ionization levels and a redshift value $z_\text{mid}$ in the middle of the large redshift step $d$:
\begin{equation}
    \log(1+z_\text{mid})=\log(1+z)-c\cdot d_0/2.
\end{equation}
Then, we approximate the large redshift step as consisting of $c$ single redshift steps, each with the same transfer function applied. The compounded transfer functions can be expressed as
\begin{equation}
\begin{aligned}
    {\mathbf{N}_\text{prop}^\gamma}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)~{\mathbf{\overline{P}}_\gamma}^c \\
    {\mathbf{N}_\text{low}^\gamma}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\left(1+{\mathbf{\overline{P}}_\gamma}+{\mathbf{\overline{P}}_\gamma}^2+\cdots+{\mathbf{\overline{P}}_\gamma}^{c-1}\right){\mathbf{\overline{D}}_\gamma} \\
    {\mathbf{N}_\text{low}^e}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\left(1+{\mathbf{\overline{P}}_\gamma}+{\mathbf{\overline{P}}_\gamma}^2+\cdots+{\mathbf{\overline{P}}_\gamma}^{c-1}\right){\mathbf{\overline{D}}_e} \\
    {E^\text{high}_c}(z')&={\mathbf{N}_\text{prop}^\gamma}(z)\left(1+{\mathbf{\overline{P}}_\gamma}+{\mathbf{\overline{P}}_\gamma}^2+\cdots+{\mathbf{\overline{P}}_\gamma}^{c-1}\right){\mathbf{\overline{D}}^\text{high}_c}, \\
\end{aligned}
\end{equation}
with all transfer function evaluated at $z_\text{mid}$. In the DNN implementation, the compounded transfer functions ${\mathbf{\overline{P}}_\gamma}^c$, $\left(1+{\mathbf{\overline{P}}_\gamma}+{\mathbf{\overline{P}}_\gamma}^2+\cdots+{\mathbf{\overline{P}}_\gamma}^{c-1}\right)$, and ${\mathbf{\overline{D}}_e}$ are learned as DNN networks, and this step can be carried out without using the value of ${\mathbf{\overline{P}}_\gamma}$ itself. (The low energy photon transfer function ${\mathbf{\overline{D}}_\gamma}$ can be quickly reconstructed using the CMB energy loss information.)

Imposing energy conservation is similar: At each $d_0$ step, the redshift energy loss is ${\mathbf{N}_\text{prop}^\gamma}(z+nd_0)\cdot d_{0i}{\mathcal{E}_i}$, where ${\mathbf{N}_\text{prop}^\gamma}(z+nd_0)={\mathbf{N}_\text{prop}^\gamma}(z)\cdot{\mathbf{\overline{P}}_\gamma}^n$. So the total redshift energy loss for a ${\mathcal{E}_i}$ photon is
\begin{equation}
     |E_{\text{redshift},i}|=\left[\left(1+{\mathbf{\overline{P}}_\gamma}+{\mathbf{\overline{P}}_\gamma}^2+\cdots+{\mathbf{\overline{P}}_\gamma}^{c-1}\right) d_0\mathcal E\right]_i.
\end{equation}
Let ${\mathbf{\overline{S}}_{\gamma,c}}=1+{\mathbf{\overline{P}}_\gamma}+{\mathbf{\overline{P}}_\gamma}^2+\cdots+{\mathbf{\overline{P}}_\gamma}^{c-1}$. Then the energy conservation condition is
\begin{equation}
\begin{aligned}
    {\mathcal{E}_i}={\mathbf{\overline{P}}_\gamma}^c_{ij}{\mathcal{E}^\gamma}_j+{\mathbf{\overline{S}}_{\gamma,c}} \Big( {\mathbf{\overline{D}}_\gamma}_{ij}{\mathcal{E}^\gamma}_j+{\mathbf{\overline{D}}_e}_{ij}{\mathcal{E}^e}_j&+{\mathbf{\overline{D}}^\text{high}_c}_i+d_i{\mathcal{E}_i} \Big) \\
    &\forall i\text{, sum on }j.
\end{aligned}
\end{equation}
Again, the propagating photon spectrum can be adjusted to account for energy non-conservation from numerical and DNN prediction errors.

\subsection{Performances in other DM injection scenarios}
\label{appd:errors}

Performance metrics of \texttt{DarkHistory}\xspace using DNN transfer functions are further demonstrated in Fig.~\ref{fig:scan_errors_decay_phot} (DM decaying to photon pairs), Fig.~\ref{fig:scan_errors_swave_elec}, (DM annihilating to electron/positron pairs with $s$-wave cross section), and Fig.~\ref{fig:scan_errors_swave_phot} (DM annihilating to photon pairs with $s$-wave cross section).

\subsection{Performances and prediction time of individual DNN transfer functions.}
\label{appd:tferrors}
The prediction accuracy and prediction time of individual DNNs are summarized in Tab.~\ref{tab:tf_err_time}.

\section*{Acknowledgements}
It is a pleasure to thank Hongwan Liu, Wenzer Qin, Gregory Ridgway, Jesse Thaler, Ken van Tilburg, Kathrin Nippel, Nils Schoeneberg, Anastasia Fialkov, Stefan Heimersheim, Julian Mu\~noz, Laura Lopez Honorez, Vivian Poulin, and Yacine Ali-Haimoud for useful conversations pertaining to this work. YS and TRS were supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics of U.S. Department of Energy under grant Contract Number DE-SC0012567 through the Center for Theoretical Physics at MIT, and the National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, \url{http://iaifi.org/}).
The computations in this paper were run on MIT LNS's Erebus machine and the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.
This research made use of the
\texttt{IPython}~\cite{PER-GRA:2007}, \texttt{Jupyter}~\cite{Kluyver2016JupyterN}, \texttt{matplotlib}~\cite{Hunter:2007},  \texttt{NumPy}~\cite{harris2020array}, \texttt{Tensorflow}~\cite{tensorflow2015-whitepaper}, \texttt{Keras}~\cite{chollet2015keras}, and \texttt{SciPy}~\cite{2020SciPy-NMeth} software packages.


\section{Introduction}

	The doping of superconductors with non-magnetic and magnetic impurities paved the way for understanding the physics of isotropic and anisotropic superconductivity \cite{Balatsky:2006_RMP,Mockli:2020_JAP,Gorkov:2008_Book}. The effects of impurities on the electronic properties of both conventional and unconventional superconductors were explained in terms of electron lifetime and pairing symmetries due to scattering by an ensemble of impurities \cite{Anderson:1959_JPCS,Abrikosov:1961_JETP,Abrikosov:1969_JETP,Balatsky:2006_RMP,Zittartz:1970_ZPhys_1,Mueller-Hartmann:1970_ZPhys,Zittartz:1970_ZPhys_2,Woolf:1965_PhysRev,Skalski:1964_PhysRev,Fulde:1966_PhysRev,Fulde:1973_Adv.Phys.,Maki:1967_PhysRev,Anderson:1959_JPCS}. According to the Cooper's pairing model \cite{Cooper:1956_PhysRev} and to the Bardeen-Cooper-Schreiffer's (BCS) theory, \cite{BCS:1957_PhysRev_1,BCS:1957_PhysRev_2}, superconductivity is a result of the instability of the Fermi surface against the pairing of time-reversed quasiparticle states. Any perturbation that is unable to lift the Kramers degeneracy of these states does not affect the superconducting transition temperature. While non-magnetic or scalar impurities are not known to affect the isotropic singlet $s$-wave order parameter in conventional BCS superconductors, magnetic impurities efficiently break the Cooper pairs, leading to the suppression of any long range superconducting order \cite{Balatsky:2006_RMP} or to the emergence of gapless superconductivity \cite{Woolf:1965_PhysRev,Abrikosov:1969_JETP,Balatsky:2006_RMP}. Most theoretical approaches \cite{Gorkov:2008_Book,Fulde:1966_PhysRev,Pint:1989_PhysC} based on the BCS theory \cite{BCS:1957_PhysRev_1,BCS:1957_PhysRev_2} treat the interaction of the impurity spin $\vec{S}$ with the spin $\vec{\sigma}$ of the conduction electrons $via$ an exchage interaction $\vec{\sigma}\cdot{\vec{S}}$. For a magnetic impurity embedded in the superconducting matrix, there is a coupling between the local spin on the impurity site and the conduction electrons of the matrix represented by the Hamiltonian:

	\begin{equation}
		H_\mathrm{imp}=\sum_{\alpha\beta}^{}\int_{}^{}d\vec{r} J(\vec{r})\Psi^{\dagger}_{\alpha}(\vec{r})\hat{U}^\dagger_{\mathrm
			imp}\Psi_{\beta}(\vec{r})
	\end{equation}
	\noindent

	where $\hat{U}^\dagger_{\mathrm{imp}}=J(\vec{r})\vec{\sigma}\cdot\vec{S}$ \cite{Maki:1967_PhysRev} and $J(\vec{r})$ is the exchange interaction between the impurity ion and the conduction electrons while $\alpha$ and $\beta$ represent the particle indices. Furthermore, it was also predicted that single magnetic impurities lead to pair breaking and result in the formation of quasiparticle Yu-Shiba-Rusinov (YSR) states within the energy gap and localized in the vicinity of the impurity atom \cite{Heinrich:2018_PSS}. However, for dilute magnetic impurities with atomic concentrations $\lesssim{10^{-4}}\%$ the spins can be considered randomly oriented and uncorrelated \cite{Pint:1989_PhysC}.
	Another approach to study the effect of magnetic impurities on superconductivity is to map the nature of superconductivity with a ferromagnetic background resulting from the proximity of a superconducting thin film to magnetic layers, promoting the onset of emergent phenomena, including spin-triplet Cooper pairing, Majorana fermions and spin superfluids \cite{Linder:2015_NatPhys,Mueller:2021_PRL}. The nature of superconductivity in strong spin-exchange fields was extensively discussed by Fulde-Ferrel \cite{Fulde:1964_PhysRev} and Larkin-Ovchinnikov \cite{Larkin:1964_JETP}. The existence of the Fulde-Ferrel-Larkin-Ovchinnnikov (FFLO) pairing was reported for heavy fermion systems, organic superconductors and recently for combinations of superconductors (S) and ferromagnets (F) in F/S/F and S/F/S heterostructures \cite{Lenk:2016_PRB,Buzdin:2005_RMP}.

	Within the family of the conventional superconductors, NbN with a bulk superconducting transition temperature of $\sim$16 K, has been widely studied both in the bulk crystal phase and as thin film \cite{Chockalingam:2008_PRB,Koushik:2013_PRL,Ganguly:2015_PRB,Chand:2009_PRB,Chand:2012_PRB,Mondal:2011_PRL,Chockalingam:2009_PRB,Destraz:2017_PRB}. Magnetron sputtered NbN thin films are intensively investigated \cite{Mondal:2011_PRL,Chand:2009_PRB,Chand:2012_PRB,Chockalingam:2009_PRB,Ganguly:2015_PRB} in view of diverse relevant applications of superconducting NbN in Josephson junctions, hot electron bolometers, single photon detectors \cite{Nikzad:2016_Sensors,Polakovic:2020_Nanomater} and in devices for quantum information and circuit quantum electrodynamics \cite{Blais:2021_RMP}. Disordered NbN thin films grown by sputtering are used as the workbench to investigate the Berezinskii-Kosterlitz-Thouless (BKT) phase transition \cite{Mondal:2011_PRL,Yong:2013_PRB}, superconductor-insulator transitions \cite{Mondal:2011_JSNM,Chand:2012_PRB}, conductance fluctuation \cite{Koushik:2013_PRL}, Andreev reflection and the Higgs-Anderson mechanism of superconductivity \cite{Sherman:2015_NatPhys,Tsuji:2020_PRRes}. Doping of conventional superconductors like NbN with magnetic impurities is expected to broaden the perspectives for hybrid structures-based applications in superconducting spintronics \cite{Linder:2015_NatPhys}, spin-orbitronics, in dark matter detectors \cite{Hochberg:2016_PRL}, integrated resonators and superconducting qubit processors for quantum computation.

	While superconducting layers of NbN are traditionally grown on conventional substrates such as Si and MgO, the use of GaN and Al$_{1-x}$Ga$_{x}$N templates and substrates has emerged as an alternative \cite{Krause:2014_SST,Sam-Giao:2014_AIPAdv,Kobayashi:2020_APE,Kobayashi:2020_APL} and can serve as basis for all nitride integrated superconductor/semiconductor devices \cite{Yan:2018_Nature}. The Si substrates are commonly employed to grow superconducting NbN films intended for single photon detection and hot bolometer applications, due to the advantages in device processing and relatively low losses at THz frequencies. The significant lattice mismatch between Si and NbN generally results in polycrystalline films with $T_\mathrm{c}\le{10}$~K for layers $\sim(5-10)$ nm thick. The lattice matched substrate MgO presents challenges in device processing, due to its hydrophobic nature and sensitivity to alkaline solutions that are used during the fabrication processes \cite{Krause:2014_SST}. Since the crystallographic orientation of the NbN films does not affect their superconducting properties, the $c$-plane of hexagonal wurtzite template/buffer layers is also suitable for NbN growth, provided that the lattice parameter $a^{\mathrm{w}}_\mathrm{hkl}$ of the wurtzite template matches the one of NbN along the (111) plane $i.e.$ $a^{\mathrm{NbN}}_{111}=\frac{a^{\mathrm{w}}_\mathrm{100}}{\sqrt{2}}$. Thus, the epitaxial relation $\left [ 111 \right ] \left ( \mathrm{NbN} \right )\parallel \left [ 100\right ] \left ( \mathrm{GaN} \right )$ is established between NbN and $c$-plane wurtzite template\cite{Krause:2014_SST,Kobayashi:2020_APL,Yan:2018_Nature}.

	The effect of dilute magnetic impurities on the superconducting properties of NbN thin films is scarcely reported in literature. The ferromagnetic proximity in the NbN/FeN bilayer system has been shown to lower the superconducting transition temperature \cite{Hwang:2017_PhysC}, and Gd ion implantation was reported to have similar effects on the superconducting properties of sputtered NbN films \cite{Jha:2013_JSNM}. In the following, the emergence of a bosonic insulator (BI) state in Fe implanted sputtered NbN thin film using low-$T$/high-$\mu_{0}H$ magnetotransport studies is analysed.

	\section{Samples and Experimental Details}

	\begin{figure*}[htbp]
		\centering
		\includegraphics[scale=1.25]{Fig_1_XRD.pdf}

		\caption{\label{fig:XRD} High resolution ($2{\theta}-{\omega}$) scan of \#A and \#B.}
		\label{fig:fig1}
	\end{figure*}

	\subsection{Growth and structure}

	The NbN films considered in this work are grown by reactive magnetron sputtering with high purity Nb as target. Epitaxial single crystalline $\sim$ 1 µm GaN grown on one-side polished $c$-plane Al$_{2}$O$_{3}$ (0001) substrates by means of metal organic vapor phase epitaxy is taken as template for the deposition of the NbN thin films. The GaN templates are cut into $\left(1\times1\right)$ cm$^{2}$ specimens for the sputtering of the NbN layers. The polycrystalline thin films are grown in an ultrahigh vacuum (UHV) chamber with a base pressure of $(2\times{10^{-9}})$ mbar. A high purity (99.99\%) Nb target is employed for the reactive magnetron sputtering process under a plasma consisting of Ar:N$_{2}$ in the ratio 10:5 standard cubic centimeters per minute (sccm) with a power $P=40$ W and at a constant substrate temperature $T_\mathrm{sub}=500^{\circ}$C.

	\begin{table*}

		\centering

		\caption{Relevant parameters of the representative samples considered in this work.}

		\renewcommand{\arraystretch}{2.0}

		\resizebox{\columnwidth}{!}{\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
				\hline
				Sample & Material & Template & Nominal $d_\mathrm{NbN}$ (nm) &  Ar:N$_{2}$ ratio  & P (W) & T$_\mathrm{sub}$ (°C) & Implanted ion & $E_\mathrm{ion}$ (keV) & Dose (at/cm$^{3}$) \\ \hline
				A          & NbN      & wz-GaN   & 100    & 10:5 & 40    & 500       & -                 & -                  & -            \\ \hline
				B          & Fe:NbN   & wz-GaN   & 100    & 10:5 & 40    & 500       & Fe                & 35                 & $1 \times 10^{14}$    \\ \hline
		\end{tabular}}
		\label{tab:tab1}
	\end{table*}

	The NbN films are implanted with $^{56}{\mathrm{Fe}}$ using powdered Fe$_{3}$O$_{4}$ as the source for Fe ions. The implantation is carried out at room temperature at a base pressure of $(1\times{10^{-6}})$ mbar. Half of the surface of the $\left(1\times1\right)$ cm$^{2}$ samples is covered with a brass mask and the uncovered part of the NbN layers is exposed to the Fe ion beam for implantation. The energy and ion dose are selected after simulating the stoppage and range of ions in matter (SRIM). The simulation of the ion distribution is considered for ion energies of 35 keV and 50 keV, respectively. A low ion current $\sim$ 100 nA is required to ignite the plasma and for the subsequent ion extraction. The ion current is kept constant by tuning the filament current and ion beam focus during the entire implantation process. The Fe implantation dose used in this work is kept constant at $\left(1\times10^{14}\right)$ atoms/cm$^{2}$. After implantation, the samples are not subjected to any thermal annealing, in order to avoid precipitation of secondary Fe or FeN phases and also to preserve the defects generated by the implanted ions. The $\left(1\times1\right)$ cm$^{2}$ samples are then cut into $\left(5\times5\right)$ mm$^{2}$ specimens for further measurements and characterization. In the following, a representative as-grown NbN layer sample A (\#A), is considered, while a 35 keV Fe implanted NbN layer is taken as sample B (\#B). Both the as-grown NbN and Fe implanted NbN are studied using x-ray diffraction (XRD) and the high crystallinity of the sputtered films is confirmed \cite{Vorhauer:2021_Thesis}. An overview of the relevant growth parameters for \#A and \#B is provided in Table I.

	Radial scans are acquired using a PANalytical X'Pert PRO Materials Research Diffractometer (MRD) equipped with a hybrid monochromator with a  divergence slit of 1/4$^\circ$. The diffracted beam is collected by a solid-state PixCel detector provided with a 9.1 mm anti-scatter slit. The $\left(2\theta-\omega\right)$ scans for \#A and \#B are shown in  Fig.~\ref{fig:fig1}. The Bragg peaks of NbN(111) and NbN(222) are observed alongside the GaN(002) and GaN(004) peaks, as indicated in Fig.~\ref{fig:fig1}. The obtained result also points at an epitaxial relation $a^{\mathrm{NbN}}_{111}=\frac{a^{\mathrm{w}}_\mathrm{100}}{\sqrt{2}}$ of the NbN (111) on $c$-Al$_2$O$_3$ (0006) as reported in literature \cite{Krause:2014_SST}. No secondary Bragg peaks of precipitated or clustered Fe and Fe$_{4}$N are detected, pointing out the homogeneous incorporation of the implanted Fe in the cubic NbN crystal lattice.

	\subsection{Magnetotransport}

	The low-$T$/high-$\mu_{0}H$ magnetotransport studies are performed on \#A and \#B in a four probe van der Pauw geometry. Ohmic contacts to the $(5\times{5})$ mm$^{2}$ specimens are achieved using electrically conducting Ag epoxy and are bonded with Au wires of diameter 25 $\mu$m. The magnetotransport experiments are carried out in a Janis Super Variable Temperature 7TM-SVM cryostat (Janis Cryogenics, Westerville, OH, USA) equipped with a 7 T superconducting magnet. A lock-in amplifier (LIA) ac technique is employed for measuring the magnetotransport properties of the NbN and of the Fe:NbN thin films. The constant current $I_\mathrm{ac}$ is sourced from a Stanford Research SR830 LIA $via$ a resistance decade box, while the longitudinal voltage $V_\mathrm{xx}$ is measured in a phase-locked mode. The lock-in expand function is employed to enhance the sensitivity of the LIA. All measurements have been performed at a frequency of 127 Hz, while $I_\mathrm{ac}=1$~$\mu$A for all measurements. The low input current minimizes the thermal drift due to Joule heating of the samples. The magnetic fields are varied between $-7$ T and $+7$ T for the magnetoresistance measurements.

	\begin{figure*}[htbp]
		\centering
		\includegraphics[scale=2.50]{Fig_2_ZFC.pdf}

		\caption{\label{fig:ZFC} ZFC $\rho_{xx}$ as a function of $T$ for samples \#A and \#B.}
		\label{fig:fig2}
	\end{figure*}

	All experimental set-ups including the LIA, the magnet power supply, and the temperature controller	are regulated by means of an indigenously developed software. Prior to the measurements, the Ohmic characteristics of the Ag epoxy contacts to the samples are tested with a high resolution Keithley 4200 SCS source-measure unit (SMU). The longitudinal resistivity $\rho_\mathrm{xx}$ is estimated as a function of temperature $T$ and magnetic field $\mu_{0}H$ from the measured longitudinal voltage $V_\mathrm{xx}$. The evolution of thermoresistivity $\rho_{\mathrm{xx}}(T)$ as a function of $T$ is estimated by cooling the sample both in the absence and in the presence of a transverse magnetic field $\mu_{0}H_\perp$ and of a longitudinal magnetic field $\mu_{0}H_\parallel$. While $\mu_{0}H_\perp$ is applied parallel to the surface normal of the sample, $\mu_{0}H_\parallel$ is applied perpendicular to it. The zero field cooled (ZFC) $\rho_\mathrm{xx}-T$ behavior is studied while the samples are cooled down for $\mu_{0}H_{\perp}=0$ ($\mu_{0}H_{\parallel}=0$) and the field cooled (FC) for $\mu_{0}H_{\perp}\ne0$ ($\mu_{0}H_{\parallel}\ne0$).

	\section{Results}

	The evolution of the ZFC $\rho_\mathrm{xx}$, measured as a function of $T$ in the interval $10~\mathrm{K}\le{T}\le16.5~\mathrm{K}$ for samples \#A and \#B is reported in Fig.~\ref{fig:fig2}. The normal state resistivity $\rho_{\mathrm{xx}}$ of \#B at $T\ge{T_\mathrm{loc}^\mathrm{onset}}$ is recorded to be $\sim$40 $\Omega\cdot{\mathrm{cm}}$, while for \#A, $\rho_{\mathrm{xx}}$ $\sim$22 $\Omega\cdot{\mathrm{cm}}$. The presence of Fe dopants in the NbN matrix is identified to be the reason for the increase of $\rho_{\mathrm{xx}}$. As evidenced in Fig.~\ref{fig:fig2}, the ZFC $\rho_\mathrm{xx}(T)$ of \#B presents five characteristic temperatures labeled by 1, 2, 3, 4 and 5. The characteristic temperatures 1, 4 and 5 are found also for the ZFC $\rho_\mathrm{xx}(T)$ of \#A. These characteristic temperatures are marked as follows:	1: $T_\mathrm{loc}^\mathrm{onset}$;	2: $T_\mathrm{d}$;	3: $T_\mathrm{peak}$; 4: $T = T^{*}$ and 5: $T_\mathrm{glo}^\mathrm{offset}$.

	It is also evidenced in  Fig.~\ref{fig:fig2}, that the thermodynamic transition from the normal state N to the superconducting state S occurs over a temperature range $\Delta{T_\mathrm{c}}$ defined as:
	\begin{equation}
		\Delta{T_\mathrm{c}}=\left(T_\mathrm{loc}^\mathrm{onset}-T_\mathrm{glo}^\mathrm{offset}\right)
	\end{equation}
	\noindent

	\begin{figure*}[htbp]
		\centering
		\includegraphics[scale=1.25]{Fig_3_RTMR.pdf}

		\caption{\label{fig:RTMR}  ZFC and FC $\rho_{xx}$ as a function of $T$ under the application of $\mu_{0}H_\perp$ = 0 T, 1 T, 3 T, 5 T and 7 T measured for samples (a) \#A and (b) \#B. Transverse $\rho_\mathrm{xx}$ as a function of $\mu_{0}H_\perp$ recorded at different $T$ for samples (c) \#A and (d) \#B.}
		\label{fig:fig3}
	\end{figure*}

	In the as-grown NbN film \#A, the dominant microstructural disorder originates from the grain boundaries of the polycrystalline film, as reported for similar polycrystalline and amorphous systems \cite{Zhang:2013_PRL,Zhang:2016_PRApp,Postolova:2017_SciRep,Mondal:2011_JSNM,Ganguly:2015_PRB}. Such a 3-dimensional (3D) system with only a structural disorder is expected to undergo a N/S transition with delayed onset of the global superconductivity. These disordered superconductors are characterized by a finite non-zero value of $\Delta{T_\mathrm{c}}$ \cite{Zhang:2016_PRApp}. The physical mechanism of this transition is described in terms of the percolation of  bosonic clusters - composed of Cooper pairs - to form bosonic islands at the $T=T_\mathrm{loc}^\mathrm{onset}$ indicated in Fig.~\ref{fig:fig2}. These bosonic islands percolate to produce a network of bosonic conduction channels, which in turn replaces the single particle fermionic conduction channels of the N state of NbN. The grain boundaries serve as scattering centers for the bosonic islands, leading to the observed kink at $T=T^{*}$, as indicated in the ZFC $\rho_{\mathrm{xx}}(T)$ curve in Fig.~\ref{fig:fig2}. At $T=T_\mathrm{glo}^\mathrm{offset}$, the global superconducting phase sets in, as evidenced in Fig.~\ref{fig:fig2}.

	However, in \#B two more characteristic temperatures indicated by $T_\mathrm{d}$ and $T_\mathrm{peak}$ are observed in the evolution of $\rho_{\mathrm{xx}}(T)$, as presented in Fig.~\ref{fig:fig2}. The anomalous peak at $T=T_\mathrm{peak}$ in the $\rho_{\mathrm{xx}}(T)$ points at the existence of disorder induced electronic phase transitions for the N/S thermodynamic phase transition in \#B, due to the presence of the implanted Fe ions. In order to understand the mechanism of the electronic processes involved in the observed behavior of \#B, FC $\rho_{\mathrm{xx}}(T)$ measurements are performed on both \#A and \#B for $\mu_{0}H_\perp$ and $\mu_{0}H_\parallel$. The behavior of $\rho_{\mathrm{xx}}(T)$ for \#A and \#B under $\mu_{0}H_\perp$ = 0 T, 1 T, 3 T, 5 T and 7 T is provided in Figs.~\ref{fig:fig3}(a) and \ref{fig:fig3}(b), respectively. It is observed, that with the increase of $\mu_{0}H_\perp$, for both samples \#A and \#B  $T_\mathrm{loc}^\mathrm{onset}$ and $T_\mathrm{glo}^\mathrm{offset}$ are shifted to lower values compared to the ones in the ZFC case. The kink in the $\rho_{\mathrm{xx}}(T)$ detected in \#A vanishes upon the application of $\mu_{0}H_\perp$. For \#B, the anomalous peak in $\rho_{\mathrm{xx}}(T)$ diminishes in amplitude, $T_\mathrm{peak}$ shifts to lower $T$ and vanishes for $\mu_{0}H_\perp=3$~T. In addition, the resistivity minima at $T = T_\mathrm{d}$ and the kink at $T = T^{*}$ are also found to quench for $\mu_{0}H_\perp\ge3$~T. The resistivity $\rho_{\mathrm{xx}}(\mu_{0}H)$ of \#A and \#B is recorded at different $T$ for both $\mu_{0}H_\perp$ and $\mu_{0}H_\parallel$. In Figs.~\ref{fig:fig3}(c) and \ref{fig:fig3}(d), the $\rho_{\mathrm{xx}}(\mu_{0}H)$ of \#A and \#B measured for an applied $\mu_{0}H_\perp$ are presented. In agreement with the results of $\rho_{\mathrm{xx}}(T)$ measurements, the anomalous peak is also observed in the measured $\rho_{\mathrm{xx}}(\mu_{0}H)$ recorded for $\mu_{0}H_\perp$. The five characteristic temperatures observed in the ZFC $\rho_{\mathrm{xx}}(T)$ behavior of \#A and \#B along with $\Delta{T_\mathrm{c}}$ are summarized in Table 2.

	\begin{table}[]
		\caption{Characteristic temperatures of samples \#A and \#B recorded for ZFC $\rho_{\mathrm{xx}}(T)$.}
		\begin{tabular}{|l|l|l|l|l|l|l|}
			\hline
			Sample & $T_\mathrm{loc,ZFC}^\mathrm{onset}$ (K)    & $T_\mathrm{d,ZFC}$ (K)    & $T_\mathrm{peak,ZFC}$  (K)  & $T_\mathrm{ZFC}^{*}$   (K)  & $T_\mathrm{glo,ZFC}^\mathrm{offset}$ (K)  & $\Delta{\mathrm{T}}$ (K) \\ \hline
			A      & 15.72 & -     & -     & 15.27  & 15.05 & 0.67 \\ \hline
			B      & 15.1  & 14.71 & 14.60 & 14.305 & 13.5  & 1.6  \\ \hline
		\end{tabular}
		\label{tab:tab2}
	\end{table}

	The $\rho_{\mathrm{xx}}(T)$ and $\rho_{\mathrm{xx}}(\mu_{0}H)$ have been also measured for an applied $\mu_{0}H_\perp$. From the evolution of $\rho_{\mathrm{xx}}(T)$ and $\rho_{\mathrm{xx}}(\mu_{0}H)$ measured at different $T$ and as a function of $\mu_{0}H_\perp$ and $\mu_{0}H_\parallel$, the characteristic temperatures namely -  $T_\mathrm{loc,\perp/\parallel}^\mathrm{onset}$, $T_\mathrm{d,\perp/\parallel}$, $T_\mathrm{peak,\perp/\parallel}$, $T_{\perp/\parallel}^{*}$ and $T_\mathrm{glo,\perp/\parallel}^\mathrm{offset}$ - are estimated. The subscripts $\perp$ and $\parallel$ refer to $\mu_{0}H_\perp$ and $\mu_{0}H_\parallel$, respectively.

	\begin{figure*}[htbp]
		\centering
		\includegraphics[scale=1.20]{Fig_4_HT.pdf}

		\caption{\label{fig:HT} The behavior of the characteristic temperatures  $T_\mathrm{loc}^\mathrm{onset}$, $T_\mathrm{d}$, $T_\mathrm{peak}$, $T = T^{*}$ and $T_\mathrm{glo}^\mathrm{offset}$ as a function of (a) $\mu_{0}H_\perp$. (b) $\mu_{0}H_\perp$. (c) Evolution of $T_\mathrm{loc,\perp/\parallel}^\mathrm{onset}$ recorded for \#B.}
		\label{fig:fig4}
	\end{figure*}

	The behavior of the characteristic temperatures of \#B as a function of  $\mu_{0}H_\perp$ and $\mu_{0}H_\parallel$ is shown in
	Figs.~\ref{fig:fig4}(a) and \ref{fig:fig4}(b), respectively. While $T_\mathrm{d,\parallel}$ and $T_\mathrm{peak,\parallel}$ are found to persist over the whole range of the applied $\mu_{0}H_\parallel$ ($i.e.$ $0\,\mathrm{T}\le{\mu_{0}H_\parallel}\le7\,\mathrm{T}$), $T_\mathrm{d,\perp}$ and $T_\mathrm{peak,\perp}$ are suppressed for $\mu_{0}H_\perp>{3\,\mathrm{T}}$. This anisotropy in the behavior of $T_\mathrm{d}$ and $T_\mathrm{peak}$ for magnetic fields applied parallel and perpendicular to the surface normal is not observed in similar disordered systems such as B-doped diamond \cite{Zhang:2013_PRL,Zhang:2016_PRApp} and $a$-InO \cite{Sacepe:2015_PRB} thin films. A comparison between $T_\mathrm{loc,\perp/\parallel}^\mathrm{onset}$ and $T_\mathrm{glo,\perp/\parallel}^\mathrm{offset}$ as a function of $\mu_{0}H$ for \#B is shown in Fig.~\ref{fig:fig4}(c). The similarity of these two characteristic temperatures for \#B indicates that the system is in a 3D regime and a 3D to 2D crossover of dimensionality, as the one reported for ultrathin TiN layers, can be ruled out \cite{Postolova:2017_SciRep}. The calculated full-width-and-half-maxima of $T_\mathrm{Peak}$ and $\Delta{T_\mathrm{c}}$ estimated for samples \#A and \#B recorded for $\mu_{0}H_{\parallel}$ and $\mu_{0}H_{\perp}$ are summarized in Table 3.

	Highly disordered superconducting thin films can show a reentrant resistive behavior, either due to a dimensionality crossover \cite{Postolova:2017_SciRep} or due to the application of a strong $\mu_{0}H$ \cite{Sacepe:2015_PRB}. The implantation of Fe into the NbN matrix enhances the disorder of the system compared to that of the as-grown NbN. The reentrant $T$ dependence of the normal state is referred to as a $N$-shaped $T$ dependence \cite{Postolova:2017_SciRep}, as the one recorded for the ZFC $\rho_{\mathrm{xx}}(T)$ of \#B. Such $N$-shaped $T$ dependence was reported in high $T_\mathrm{c}$ superconductors (HTS) \cite{Peng:2013_NatCom,Tagaki:1992_PRL,Ono:2000_PRL,Semba:2001_PRL,Komiya:2005_PRL,Oh:2006_PRL,Moshchalkov:2001_PRB} and ultrathin films of conventional superconductors \cite{Postolova:2017_SciRep,Sacepe:2015_PRB,Zhang:2013_PRL,Zhang:2016_PRApp}. Even though this phenomenon is ubiquitious, it still remains a subject of debate. In HTSs, this effect is described in terms of scaling functions \cite{Oh:2006_PRL,Moshchalkov:2001_PRB} or attributed to the emergence of the pseudogap phase \cite{Daou:2009_NatPhys}.

	\begin{table}[]
		\caption{The full width and half maxima of the $T_\mathrm{Peak}$ and $\Delta{T_\mathrm{c}}$ estimated for samples \#A and \#B recorded for $\mu_{0}H_{\parallel}$ and $\mu_{0}H_{\perp}$.}
		\begin{tabular}{|l|ll|ll|ll|ll|}
			\hline
			$\mu_{0}H_{\parallel/\perp}$ (T) & \multicolumn{2}{l|}{$\Delta{\mathrm{T}_\mathrm{Peak}^\mathrm{A}} (\mathrm{K})$}    & \multicolumn{2}{l|}{$\Delta{\mathrm{T}_\mathrm{Peak}^\mathrm{B}} (\mathrm{K})$}         & \multicolumn{2}{l|}{$\Delta{\mathrm{T}_\mathrm{c}^\mathrm{A}} (\mathrm{K})$}      & \multicolumn{2}{l|}{$\Delta{\mathrm{T}_\mathrm{c}^\mathrm{B}} (\mathrm{K})$}      \\ \hline
			& \multicolumn{1}{l|}{$\mu_{0}H_{\perp}$} & {$\mu_{0}H_{\parallel}$} & \multicolumn{1}{l|}{$\mu_{0}H_{\perp}$}   & {$\mu_{0}H_{\parallel}$}    & \multicolumn{1}{l|}{$\mu_{0}H_{\perp}$} & {$\mu_{0}H_{\parallel}$}  & \multicolumn{1}{l|}{{$\mu_{0}H_{\parallel}$}} & {$\mu_{0}H_{\parallel}$}  \\ \hline
			0     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{0.0653} & 0.0653 & \multicolumn{1}{l|}{0.67} & 0.67 & \multicolumn{1}{l|}{1.6}  & 1.6  \\ \hline
			1     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{0.1376} & 0.1337 & \multicolumn{1}{l|}{0.72} & 0.94 & \multicolumn{1}{l|}{1.64} & 1.71 \\ \hline
			2     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{0.0858} & 0.1514 & \multicolumn{1}{l|}{-}    & -    & \multicolumn{1}{l|}{1.5}  & 1.34 \\ \hline
			3     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{-}      & 0.1392 & \multicolumn{1}{l|}{0.74} & 0.94 & \multicolumn{1}{l|}{1.25} & 1.54 \\ \hline
			4     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{-}      & 0.0841 & \multicolumn{1}{l|}{-}    & -    & \multicolumn{1}{l|}{1.23} & 1.25 \\ \hline
			5     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{-}      & -      & \multicolumn{1}{l|}{0.8}  & 0.96 & \multicolumn{1}{l|}{1.12} & 1.37 \\ \hline
			6     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{-}      & -      & \multicolumn{1}{l|}{-}    & -    & \multicolumn{1}{l|}{1.31} & 1.05 \\ \hline
			7     & \multicolumn{1}{l|}{-}    & -   & \multicolumn{1}{l|}{-}      & -      & \multicolumn{1}{l|}{0.78} & 0.91 & \multicolumn{1}{l|}{0.96} & 1.18 \\ \hline
		\end{tabular}
		\label{tab:tab3}
	\end{table}

	\begin{figure*}[htbp]
		\centering
		\includegraphics[scale=1.75]{Fig_5_ZFCM.pdf}

		\caption{\label{fig:ZFCM} ZFC $\rho_{\mathrm{xx}}$ as a function of $T$. The six electonic phases and the phase boundaries of \#B. are evidenced}
		\label{fig:fig5}
	\end{figure*}

	In the case of ultrathin $<15$ nm layers of conventional superconductors, quantum contributions due to a 3D to 2D dimensionality crossover have been identified as the mechanism behind the observed $N$-shaped $T$ dependence of $\rho_{\mathrm{xx}}$. In 3D films of disordered granular superconductors such as B-doped diamond \cite{Zhang:2013_PRL,Zhang:2016_PRApp}, $a$-InO \cite{Sacepe:2015_PRB} or AlGe \cite{Zaken:1994_JPCM} this behavior is interpreted within the framework of a N-BI-S transition originating from granularity induced disorder.

	In the ZFC $\rho_{\mathrm{xx}}(T)$ of \#B a series of electronic phases can be identified. These electronic phases are marked by EP-1, EP-2, EP-3. EP-4, EP-5 and EP-6 in Fig.~\ref{fig:fig5}. The phase boundaries of the adjacent electronic phases of \#B correspond to the characteristic temperatures $T_\mathrm{loc}^\mathrm{onset}$, $T_\mathrm{d}$, $T_\mathrm{peak}$, $T^{*}$ and $T_\mathrm{glo}^\mathrm{offset}$, as described earlier.

	The evolution of the six electronic phases, $i.e.$ from EP-1 to EP-6, across the phase boundaries is explained using an empirical model of competing impurity scattering and a conduction network of fermionic conduction channel and bosonic conduction channels  across a N-BI-S phase transition. A schematic of this empirical model is presented in Fig.~\ref{fig:fig6}. The Fe implantation dose of $\left(1\times10^{14}\right)$ at/cm$^3$ results in a highly dilute Fe doped system, restricting the Fe-Fe interaction to the paramagnetic limit. Since after implantation \#B is not thermally treated, the Fe ions are incorporated into the NbN lattice randomly. Moreover, the incorporation efficiency of implanted Fe ions is different for NbN grains with different crystallographic orientations. The disorder of the Fe implanted NbN system stems from an extrinsic granularity and an intrinsic one. The disorder due to the polycrystalline texture and grain boundaries of the Fe:NbN lattice is the extrinsic granularity. On the other hand, the electronic disorder due to the randomly implanted Fe ions results in modulations of the chemical potential and thus in intrinsic electronic granularity. The intrinsic granularity plays a significant role in defining the observed electronic properties of \#B.

	The mechanisms leading to the specific electronic phases are:

	(i) Phase EP-1: The normal electronic phase in the ZFC $\rho_{\mathrm{xx}}(T)$ of Fe:NbN sets in for $T>T_\mathrm{loc}^\mathrm{onset}$, as shown in Fig.~\ref{fig:fig5}. This electronic phase is represented in Fig.~\ref{fig:fig6}(a) and is referred to as EP-1. For $T>T_\mathrm{loc}^\mathrm{onset}$, the electrical transport takes place $via$ the percolation of thermally activated fermionic conduction channels with normal electrons being the medium of transport. These fermionic conduction channels are marked in the schematic diagram as solid lines and the fermions are visualized by solid ellipses. The Fe impurities in the matrix are represented by solid circles. The resistivity behavior of this phase is the same as the one of the N phase of Fe:NbN.

	(ii) Phase EP-2: As $T$ decreases, a superconducting gap opens at $T=T_\mathrm{loc}^\mathrm{onset}$. A spatial inhomogeneity associated with the superconducting gap due to the intrinsic granularity leads to the formation of bosonic clusters. With a gradual decrease of $T$ below $T=T_\mathrm{loc}^\mathrm{onset}$, the density of the bosonic clusters increases and becomes locally phase locked, resulting in the initial decrease of $\rho_{\mathrm{xx}}$ with decreasing $T$ through the percolation of the phase locked bosonic clusters merging into bosonic islands. These islands are responsible for the bosonic conduction channels for electronic transport in the temperature range $T_\mathrm{d}<T<T_\mathrm{loc}^\mathrm{onset}$. The bosonic conduction channels and the bosonic islands are represented in Fig.~\ref{fig:fig6}(b) by connectors and dumb bells, respectively. In this regime, a competition of the bosonic conduction channels with the fermionic conduction channels leads to a decrease in $\rho_{\mathrm{xx}}$.

	\begin{figure*}[htbp]
		\centering
		\includegraphics[scale=1.50]{Fig_6_SCHM.pdf}

		\caption{\label{fig:SCHM} Schematic of the empirical model to describe the N-BI-S transition in \#B. The connector lines connecting the ellipses and dumbbells represent fermionic conduction channels and bosonic conduction channels, respectively. The normal state fermions and the Cooper pairs are shown by ellipses and dumbbells, respectively, while the Fe impurities are depicted by solid circles. The electronic phase transitions taking place as the sample \#B is cooled from $T>T_\mathrm{loc}^\mathrm{onset}$ down to $T\le{T_\mathrm{loc}^\mathrm{onset}}$ are: (a) Phase EP-1 for $T>T_\mathrm{loc}^\mathrm{onset}$. (b) Phase EP-2 for $T_\mathrm{d}<T<T_\mathrm{loc}^\mathrm{onset}$. (c) Phase EP-3 for $T^\mathrm{peak}<T<T_\mathrm{d}$. (d) Phase EP-4 for $T_\mathrm{Peak}\le{T}\le{T^{*}}$. (e) Phase EP-5 for $T^{*}\le{T}\le{T_\mathrm{glo}^\mathrm{offset}}$ and (f) Phase EP-6 for $T\le{T_\mathrm{glo}^\mathrm{offset}}$.}
		\label{fig:fig6}
	\end{figure*}

	(iii) Phase EP-3: As the $T$ is further reduced, $i.e.$ for $T<T_\mathrm{d}$, in the presence of the intrinsic granularity due to the paramagnetic Fe impurities, the bosonic conduction channels dominate over the fermionic conduction channels $via$ the removal of single fermions or normal electrons from the fermionic conduction channels to fuel the bosonic conduction channels. This process takes place by forming groups of phase-locked bosonic clusters and producing isolated BIs. The consequent anomalous increase of $\rho_{\mathrm{xx}}$ for $T^\mathrm{peak}<T<T_\mathrm{d}$ is likely due to two mechanisms, $i.e.$

	- removal of the single fermions from the fermionic conduction channels to estblish the percolation path for the bosonic conduction channels;

	- scattering of the bosonic clusters from the isolated and randomly distributed Fe ions in the matrix.

	The anomalous increase of $\rho_{\mathrm{xx}}$ continues up to $T=T_\mathrm{peak}$ The electronic phase EP-3 is referred to as a bosonic insulator phase and is represented in Fig.~\ref{fig:fig6}(c).

	(iv) Phase EP-4: As $T$ is further lowered, the density of the bosonic conduction channels increases $via$ the percolation of more bosonic islands, leading to a decrease in $\rho_{\mathrm{xx}}$, as more and more phase locked coherent bosonic clusters and islands percolate in the bulk of the system. The dissipationless transport due to the coherent bosonic islands dominates over the scattering from the paramagnetic Fe impurities, leading to the gradual decrease of $\rho_{\mathrm{xx}}$ until $T = T^{*}$. The EP-4 phase is shown in Fig.~\ref{fig:fig6}(d).

	(v) Phase EP-5: At $T = T^{*}$, a kink in $\rho_{\mathrm{xx}}(T)$ is observed, which is attributed to the scattering of the percolating bosonic islands due the extrinsic granularity. This is also detected in the as-grown NbN \#A, as shown in Fig.~\ref{fig:fig6}(e).

	(vi) Phase EP-6: At $T\le{T_\mathrm{loc}^\mathrm{onset}}$, the global phase coherence sets in between all the phase-locked bosonic islands and the electrical transport is dominated solely by the bosonic conduction channels, leading to a superconducting state of the Fe implanted NbN layer, as sketched in Fig.~\ref{fig:fig6}(f).

	The empirical model also describes the behavior of $\rho_{\mathrm{xx}}(T)$ for applied $\mu_{0}H_\perp$ and $\mu_{0}H_\parallel$. The presence of $\mu_{0}H_\perp$ breaks the Cooper pairing and shifts $T_\mathrm{loc}^\mathrm{onset}$ to a lower $T$. This also suppresses the onset of the bosonic insulator phase. In addition, the formation of the Abrikosov vortices for $\mu_{0}H_\perp$ also plays a role in the suppression of the BI phase and of the intrinsic granularity on the electronic transport of \#B.

	For $\mu_{0}H_\parallel$ however, the fields required to break the Cooper pairs are generally orders of magnitude stronger than the ones for $\mu_{0}H_\perp$ and do not influence the formation of the BI state. As a result both $T^\mathrm{peak,\parallel}$and $T_\mathrm{d,\parallel}$ persist up to the highest applied field  $\mu_{0}H_\parallel=7\,\mathrm{T}$.

	\section{Conclusions}

	In conclusion, 100 nm thick crystalline films of NbN are deposited by reactive magnetron sputtering on single crystalline GaN templates grown on epi ready $c$-plane sapphire substrates. Paramagnetic Fe doping of the NbN layers is achieved by ion implantation using Fe$_{3}$O$_{4}$ as the source for Fe ions. An implantation energy of 35 keV and a dose of $\left(1\times10^{14}\right)$ at/cm$^3$ result in a highly dilute Fe doped system, restricting the Fe-Fe interaction to the paramagnetic limit. Low-$T$/high-$\mu_{0}H$ magnetotransport measurements confirm, that the Fe doping does not suppress the superconductivity of the sputtered NbN thin films, but decreases the superconducting transition temperature. A reentrant resistive BI phase is observed in the Fe doped NbN samples, which is explained by an empirical model of a competition between the percolation of bosonic conduction channels at the expense of fermionic conduction channels and by the scattering of the BI as a result of intrinsic granularity due to the random Fe dopants. The observation of a robust superconductivity in the dilute magnetic conventional superconductor Fe:NbN mediated $via$ percolation of bosonic insulator states is foreseen to violate the symmetry of electron-like and hole-like excitations due to the formation of subgap bound Andreev states in the vicinity of magnetic impurities, leading to giant thermoelectric effects \cite{Kalenkov:2012_PRL}. A system like the one reported in this work is expected to find applications in zero biased thermoelectric bolometers with reduced power dissipation in large scale multi-pixel arrays and in hybrid quantum interference devices (HyQUID) \cite{Shelley:2021_SciAdv}. Moreover, these systems are the workbench for understanding quantum emergent phenomena, including gapless superconductivity, triplet Cooper pairings, YSR states and odd frequency superconductivity.

	\section*{Acknowledgements}

	The work was funded by the Austrian Science Fund (FWF) through Projects No. P26830 and No. P31423.

	\bibliographystyle{apsrev4-2}