Fingerprint matching of beyond-WIMP dark matter: neural network approach

TL;DR

This paper introduces a neural network approach to analyze and compare beyond-WIMP dark matter models by capturing complex suppression patterns in galactic structure formation.

Contribution

It proposes using neural networks to effectively characterize and communicate the suppression features of various beyond-WIMP dark matter models.

Findings

Neural networks can model complex suppression shapes in matter power spectra.

The approach facilitates comparison across different beyond-WIMP models.

Demonstrated on a simplified light feebly interacting massive particles model.

Abstract

Galactic-scale structure is of particular interest since it provides important clues to dark matter properties and its observation is improving. Weakly interacting massive particles (WIMPs) behave as cold dark matter on galactic scales, while beyond-WIMP candidates suppress galactic-scale structure formation. Suppression in the linear matter power spectrum has been conventionally characterized by a single parameter, the thermal warm dark matter mass. On the other hand, the shape of suppression depends on the underlying mechanism. It is necessary to introduce multiple parameters to cover a wide range of beyond-WIMP models. Once multiple parameters are introduced, it becomes harder to share results from one side to the other. In this work, we propose adopting neural network technique to facilitate the communication between the two sides. To demonstrate how to work out in a concrete manner, we consider a simplified model of light feebly interacting massive particles.

Equations

$$f_{\chi }(t_f,p)={\int }_{t_i}^{t_f}dt\frac{1}{E_{\chi }}C\left(t,\frac{a(t_f)}{a(t_i)}p\right),$$

$$q^{2}f(q)$$

$$T^2(k)\equiv \frac{P(k)}{P_{\mathrm{CDM}}(k)}.$$

$${\Omega }_{\mathrm{WDM}}{h}^2=\left(\frac{m_{\mathrm{WDM}}}{94\mathrm{eV}}\right){\left(\frac{T_{\mathrm{WDM}}}{T_{\nu }}\right)}^3=7.5\left(\frac{m_{\mathrm{WDM}}}{7\mathrm{keV}}\right)\left(\frac{106.75}{g_{\ast }^{\mathrm{WDM}}}\right).$$

$$\frac{d{N}_{\mathrm{sat}}}{d\ln M}=\frac{1}{C_n}\frac{1}{6{\pi }^2}\frac{M_0}{M}\frac{P(1/R)}{R^3\sqrt{2\pi (S-S_0)}},$$

$$S=\frac{1}{2{\pi }^2}{\int }_0^{1/R}dk{k}^{2}P(k).$$

$$M=\frac{4\pi }{3}{\rho }_m(cR{)}^3,$$

$$P_{1D}(k)$$

$$r(k)=\frac{P_{1D}(k)}{P_{1D}^{\mathrm{CDM}}(k)}.$$

$$A={\int }_{k_{\min }}^{k_{\max }}dkr(k),$$

$$\delta A=\frac{A_{\mathrm{CDM}}-A}{A_{\mathrm{CDM}}}.$$

$$T_{\mathrm{WDM}}^2(k)$$

$$\alpha =0.049\mathrm{Mpc}/h{\left(\frac{m_{\mathrm{WDM}}}{\mathrm{keV}}\right)}^{-1.11}{\left(\frac{{\Omega }_{\mathrm{WDM}}}{0.25}\right)}^{0.11}{\left(\frac{h}{0.7}\right)}^{1.22}$$

$$T^2(k)$$

$${\mathcal{L}}_{\chi }$$

$$T_{\mathrm{DM}}={\left(\frac{g_{\ast }(T)}{g_{\ast }(T_{\mathrm{dec}})}\right)}^{1/3}T,$$

$$T_{\mathrm{DM}}={\left(\frac{g_{\ast }(T)}{\Delta \times g_{\ast }^{\mathrm{SM}}}\right)}^{1/3}T.$$

$$\frac{m_2}{m_1},m_{\mathrm{DM}},\Delta .$$

$${\mathcal{L}}_{\varphi }$$

$$\frac{m_2}{m_1},\frac{Y_{\mathrm{scat}}}{Y_{\mathrm{total}}},m_{\mathrm{DM}}.$$

$${\vec{x}}_1$$

$${\vec{x}}_n$$

$${\vec{x}}_{\mathrm{out}}$$

$$\vec{f}(\vec{y})$$

$$f(y)$$

$$E$$

$${\vec{x}}_{\mathrm{out}}$$

$$\gamma$$

$${\vec{x}}_{\mathrm{in}}$$

Full text

1. 1 Introduction

2. 2 Procedure for FIMP DM as an example

3. 2.1 Model $\to$ DM phase space distribution $\to$ Linear matter power spectrum

4. 2.2 Linear matter power spectrum $\to$ Observables

5. 2.3 $\{\alpha,\beta,\gamma\}$ parametrization of the transfer function

6. 3 Simplified FIMP model

7. 3.1 Freeze-in production

8. 3.2 Constraints

9. 4 Neural network approach

10. 4.1 Neural network setup

11. 4.2 Model parameters $\to$ $\{\alpha,\beta,\gamma\}$

12. 4.3 $\{\alpha,\beta,\gamma\}$ $\to$ Observables

13. 4.4 Combined results

14. 5 Summary

15. A Comparison with an analytic map

16. B Further check: original data vs neural network

17. C How to use the neural network data

1 Introduction

Dark matter (DM) is an essential component for the Universe to form the current shape. Its existence and abundance are probed by gravitational observations such as galaxy rotation curves, bullet cluster collision, and cosmic microwave background (CMB) anisotropy. On the other hand, we have not seen any DM signal by any non-gravitational interactions, and thus we still do not know the identity of DM: what it is and how it is produced. One intriguing possibility is that DM consists of a new particle, which provides a clue to physics beyond the standard model (SM) (see Ref. [1] for a review).

One of the early attempts is a weakly interacting massive particle (WIMP) (see Refs. [2, 3] for recent reviews). In this direction, much efforts have been devoted at the large hadron collider (LHC) (for example, mono-jet searches [4, 5]) and at direct/indirect detection searches [6, 7, 8, 9]. However, no firm signals have been reported yet. It may motivate us to consider beyond-WIMP scenarios that can be probed by cosmological/astrophysical observations.♢♢\diamondsuit1♢♢\diamondsuit1$\diamondsuit$1 We refer readers to Ref. [10] for a recent review of gravitational probes of DM properties.

WIMPs behave as cold dark matter (CDM) on galactic scales. They are in good agreement with many independent observations such as CMB anisotropy [11] and galaxy clustering [12]. On the other hand, their predictions of galactic-scale structure are in debate. On galactic scales, there have been issues that are difficult to explain in CDM (small-scale issues).♢♢\diamondsuit2♢♢\diamondsuit2$\diamondsuit$2 Prominent examples are the missing satellite problem [13, 14, 15, 16, 17], core-cusp problem [18, 19, 20, 21], and too-big-to-fail problem [22, 23, 24, 25, 26, 27]. We refer readers to Ref. [28] for a recent review and further details. State-of-the-art hydrodynamical simulations have been demonstrating that astrophysical processes also play an important role [29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. There have also been reports that small-scale issues persist even in state-of-the-art hydrodynamical simulations [47, 48, 49, 50, 51, 52, 53, 54, 55]. To our best knowledge, it is still controversial if astrophysical processes fully resolve the small-scale issues.

Alternatives to CDM may explain small-scale issues: warm dark matter (WDM) [56, 57, 58, 15, 59, 16, 60, 17, 52]; fuzzy DM [61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71]; and long-lasting DM interaction with primordial plasma or free-streaming light particles [72, 73, 74, 75, 76, 60, 77, 78, 79, 80, 81, 82, 83, 84, 85].

On the other hand, impacts on galactic-scale structure formation depend on beyond-WIMP scenarios. Free-streaming of light WDM particles smears out the primordial density contrast. Quantum pressure of fuzzy DM prevents DM from gravitational clustering. Pressure of radiation to which DM couples involves DM in acoustic oscillation rather than gravitational clustering. Such effects are reflected in the linear matter power spectrum, which one can obtain by following evolution of the primordial density contrast. Generally by performing a suit of simulations with the resulting linear matter power spectrum, one can obtain observable quantities, which can be directly compared with cosmological/astrophysical observations. In summary, we need to work out the following procedure on a model-by-model basis:

• Model $\to$ Linear matter power spectrum $\to$ Observables.

See the blue flow in Fig. 1. The whole procedure requires interdisciplinary expertise from particle phenomenology to (computational) astrophysics. Moreover, each step often requires a dedicated calculation. In particular, simulations in the last step are often too time-consuming to repeat.

One can work out each step independently by parametrizing the linear matter power spectrum. See the red flow in Fig. 1. A single parameter has been adopted conventionally: the thermal WDM mass $m_{\rm WDM}$.♢♢\diamondsuit3♢♢\diamondsuit3$\diamondsuit$3 An underlying model may be light gravitino [86, 87]. WDM particles are thermalized in the early Universe and decouple from thermal plasma at some point.

On the other hand, a single parameter is not enough to cover a wide range of beyond-WIMP scenarios. For this purpose, Ref. [88] introduces the 3-parameter ($\{\alpha,\beta,\gamma\}$) characterization of the linear matter power spectrum. On one side, one (likely particle physicist) can construct a map of model parameters onto $\{\alpha,\beta,\gamma\}$. On the other side, one (likely astrophysicist) can provide observational constraints on $\{\alpha,\beta,\gamma\}$, as indeed done for the Lyman-$\alpha$ forest data in Ref. [89]. By combining results from the two sides, one can obtain observational constraints on a given beyond-WIMP scenario. Nevertheless, once multiple parameters are introduced, it becomes hard to share results from one side to the other.

In this respect, we propose building ready-to-use networks: one maps model parameters onto $\{\alpha,\beta,\gamma\}$; and another maps $\{\alpha,\beta,\gamma\}$ onto observables. One can use these networks to examine models without repeating the aforementioned time-consuming procedure. Ideally, it would be the most efficient if one obtained analytic maps, but in reality, it is hard to establish such analytic maps. Thus, a numerical method is helpful to develop such effective maps. For this purpose, we adopt neural network technique.

To be concrete, in this paper, we consider a feebly interacting massive particle (FIMP) [90] (see Ref. [91] for a recent review). Light (keV-scale) FIMPs, which are produced through the freeze-in mechanism, are a compelling example of WDM. Even in FIMP models, the shape of suppression in the linear matter power spectrum depends on production processes such as 2-body decay, 3-body decay, and 2-to-2 scattering [92, 93, 94, 95, 96, 97, 98] (see Ref. [99] for a comprehensive discussion).♢♢\diamondsuit4♢♢\diamondsuit4$\diamondsuit$4 We refer readers to Refs. [100, 101, 102, 103] for sterile neutrino DM. Sterile neutrinos are produced through mixing with active neutrinos. We also refer readers to Ref. [104] for superWIMPs. SuperWIMPs are produced by the decay of WIMPs long after the WIMP freeze-out. If the WIMP decay occurs close after the WIMP freeze-out, one may need to take into account the momentum distribution function of WIMPs [105, 106, 107, 108, 109]. In this paper, we do not consider these possibilities, although they may be FIMPs in a broad sense.

Thus 3-parameter characterization rather than conventional single-parameter characterization is required to cover a wide range of FIMP models. By taking a simplified FIMP model, we demonstrate how one can work out the simplified procedure. We also provide the obtained neural networks through the arXiv website: one is a map of “model parameters $\to$ $\{\alpha,\beta,\gamma\}$” and the others are “$\{\alpha,\beta,\gamma\}$ $\to$ observables”.

The organization of this paper is following. In Sec. 2, we overview the conventional procedure to place constraints on FIMPs and describe the simplified procedure with the $\{\alpha,\beta,\gamma\}$ parametrization. In Sec. 3, we introduce a simplified FIMP model. Our FIMP model shares many common aspects with a broad class of FIMP models. The basic production process is $2$-body decay. We take into account late-time entropy production after freeze-in (case A) and also freeze-in production through $2$-to-$2$ scattering (case B). In Sec. 4, we introduce neural network technique and work out the simplified procedure. We compare the constraints from the simplified procedure and those from the conventional procedure. Sec. 5 is devoted to the summary. In Appendix A, we compare our constraints to those obtained through an analytic map from the conventional thermal WDM mass. In Appendix B, we examine precision of the neural networks in detail. In Appendix C, we explain how to use the neural networks we provide.

2 Procedure for FIMP DM as an example

As we described in introduction, to study galactic-scale structure formation of beyond-WIMP scenarios, generically one has to take a 2-step procedure on a model-by-model basis:

• Model $\to$ Linear matter power spectrum $\to$ Observables

(corresponding to the blue flow in Fig. 1). In the case of FIMP, the first step of “Model $\to$ Linear matter power spectrum” actually consists of two steps:

• Model $\to$ DM phase space distribution $\to$ Linear matter power spectrum.

To follow the two steps, one first needs to construct the collision term of the Boltzmann equation and integrate it to obtain the phase space distribution of the DM species. Then one has to follow evolution of the primordial density contrast with the obtained phase space distribution, possibly by using public cosmological Boltzmann solvers such as CLASS [110, 111]. In the following we overview this conventional procedure more specifically.

2.1 Model $\to$ DM phase space distribution $\to$ Linear matter power spectrum

We define the DM phase space distribution $f_{\chi}(t,p)$ as a function of the cosmic time $t$ and the physical momentum $p$, such that the DM number density is given by $n_{\chi}=g_{\chi}\int d^{3}p/(2\pi)^{3}f_{\chi}(t,p)$, where $g_{\chi}$ is the spin degrees of freedom. We assume that the DM phase space distribution $f_{\chi}$ is much smaller than unity. We then obtain the phase space distribution at a late cosmic time $t_{f}$ by integrating the collision term as

[TABLE]

where $t_{i}$ is the reheating time and $a(t)$ is the cosmic scale factor. Given a squared matrix element of a specific production process, one obtains a semi-analytic expression of the corresponding collision term (see Ref. [99] for expressions).

FIMP production is most efficient when the heaviest particle in the process becomes non-relativistic (freeze-in mechanism). After that, FIMP particles free-stream and the phase space distribution is invariant as a function of the comoving momentum $q\equiv p/T_{\chi}$, where $T_{\chi}$ is the effective DM temperature (see Sec. 3.1 for a specific expression of $T_{\chi}$). Thus we use $f(q)=f(t,p)$ to characterize the distribution. Practically, we fit the obtained phase space distribution of DM by

[TABLE]

where $(a_{i},b_{i},c_{i})$ are fitting parameters and $i$ runs for different production processes.♢♢\diamondsuit5♢♢\diamondsuit5$\diamondsuit$5 One may wonder if we can work out “Model $\to$ DM phase space distribution” and “DM phase space distribution $\to$ Linear matter power spectrum” separately by using $(a,b,c)$. On one side, one can report constraints on $(a,b,c)$. On the other side, one can calculate $(a,b,c)$ as a function of model parameters. It is worth investigating this possibility somewhere else.

We plug the fitting function into the public cosmological Boltzmann solver CLASS [110, 111] to obtain the linear matter power spectrum $P(k)$ as a function of the wavenumber $k$. We use the cosmological parameters from “Planck 2015 TT, TE, EE+lowP” in Ref. [112]. Practically, we use the CLASS fluid approximation of non-cold DM.

2.2 Linear matter power spectrum $\to$ Observables

Galactic-scale structure places constraints on the linear matter power spectrum $P(k)$, or, the transfer function that is defined by

[TABLE]

It generically requires a suit of time-consuming simulations to obtain constraints on FIMP DM. We may simplify this step by using semi-analytic models and/or somehow converting the conventional thermal WDM mass $m_{\rm WDM}$.

In the conventional thermal WDM model, WDM particles follow the Fermi-Dirac distribution with two spin degrees of freedom with temperature $T_{\mathrm{WDM}}$. The relic abundance is expressed by $m_{\mathrm{WDM}}$ and $T_{\mathrm{WDM}}$ as

[TABLE]

For a given WDM mass, the temperature is determined such that the relic abundance reproduces the observed DM density. Note that for a keV-scale mass, somewhat large entropy production after decoupling is required for $\Omega_{\mathrm{WDM}}h^{2}\simeq 0.12$. On the other hand, FIMP DM has a different thermal history and thus different temperature and does not follow the Fermi-Dirac distribution. Thus reported lower bounds on $m_{\rm WDM}$ is not directly applicable to FIMP DM.

In this paper we consider the number of satellite galaxies [113, 114, 115, 116, 117, 118] and Lyman-$\alpha$ forest [119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132] as observables.♢♢\diamondsuit6♢♢\diamondsuit6$\diamondsuit$6 Other used probes include the delay of the reionization [133, 134, 135, 136, 137, 138, 139], the counts of high-$z$ gamma-ray bursts [140, 141], the faint end of luminosity function of high-$z$ galaxies [142, 135, 136, 143, 144, 145, 139, 146], the flux anomaly of quadrupole lens systems [147, 148, 149, 150, 151, 152, 153, 154], and the redshifted 21 cm signal [155, 156, 157, 158, 159, 160, 161, 162]. The counts of lensed distant supernovae [163] and direct collapse black holes [164] are suggested for a future use. We also refer readers to Ref. [165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179] for hydrodynamical simulation results differentiating WDM and CDM in galaxy formation.

Our analysis, which follows Refs. [180, 88], uses a semi-analytic model for the number of satellite galaxies and converts the reported lower bound on $m_{\rm WDM}$ for the Lyman-$\alpha$ forest.

Number of satellite galaxies

One compares the predicted number of satellite galaxies $N_{\rm sat}$ in simulated Milky Way-size or M31-size haloes with the observed one. If the predicted number is smaller than the observed one, such FIMPs are excluded. This constraint may be conservative when one counts all the subhalos above a certain mass, since some of subhalos may not host galaxies bright enough to be detected.

We evaluate the number of satellite galaxies $N_{\rm sat}$ from the linear matter power spectrum in our FIMP model as follows. Ref. [118] develops a semi-analytic formula of the subhalo mass function in the conventional thermal WDM model. The formula uses the conditional mass function [181] based on the extended Press-Schechter approach [182] and the halo model (see Ref. [183] for a review). The formula adopts the top-hat filter function in the Fourier space (sharp-$k$ filter) to reproduce results of $N$-body simulations in the conventional thermal WDM model:

[TABLE]

where quantities with and without the subscript “0” denote those of the host halo and subhalo, respectively. For example, $M$ ($M_{0}$) is the subhalo (host halo) mass. The variance $S$ is given by the linear matter power spectrum as

[TABLE]

The filter scale $R$ is related with the mass as

[TABLE]

with the matter mass density at present $\rho_{m}$. Following Ref. [118], we adopt $c=2.5$ and $C_{n}=44.5$. We use $M_{0}=1.7\times 10^{12}\,h^{-1}\mathrm{M_{\odot}}$ as the Milky-Way mass, where $h$ is the dimensionless Hubble constant. With these values, the number of satellite galaxies above $M=10^{8}\,h^{-1}\mathrm{M_{\odot}}$ is $N_{\rm sat}=159$, which is consistent with the result of the Aquarius simulation [184]. $M>10^{8}\,h^{-1}\mathrm{M_{\odot}}$ roughly corresponds to the lower bound on the maximal circular velocity of $V_{\rm max}>10$ km/s.

We estimate the observed number of satellites above $M=10^{8}\,h^{-1}\mathrm{M_{\odot}}$ as $N_{\rm sat}^{\rm obs}=63$ (11 classical dwarf galaxies and $3.5\times 15$ ultra-faint dwarf galaxies).♢♢\diamondsuit7♢♢\diamondsuit7$\diamondsuit$7 Classical dwarfs: Sagittarius, LMC, SMC, Ursa Minor, Sculptor, Draco, Sextans, Carina, Fornax, LeoII, and LeoI. Ultra-faint dwarfs: Segue I, Ursa Major II, Segue II, Willman I, Coma Berenics, Bootes II, Bootes I, Pisces I, Ursa Major I, Hercules, Canes Venatici II, Leo IV, Leo V, Pisces II, Canes Venatici I. We refer readers to Refs. [185, 186] for dynamical properties. Note that $V_{\rm max}\geq V_{\rm 1/2}\simeq\sqrt{3}\sigma_{\rm l.o.s}$, where $V_{\rm 1/2}$ is the circular velocity at the half light radius and $\sigma_{\rm l.o.s}$ is the line-of-sight velocity dispersion [186].

We multiply 3.5 by the number of ultra-faint satellites found in SDSS to take account of the SDSS limited sky coverage as in Refs. [114, 60, 115, 116, 118, 180, 88]. $N_{\rm sat}>N_{\rm sat}^{\rm obs}$ places a lower bound on the conventional thermal WDM mass as $m_{\rm WDM}>2.9$ keV.♢♢\diamondsuit8♢♢\diamondsuit8$\diamondsuit$8 We remark that we do not use the fitting function given by Eq. (2.12), but directly compute the linear matter power spectrum by using CLASS [110, 111].

As we see, $N_{\rm sat}$ implicitly depends on the lower bound on the satellite mass. For example, once a number of smaller-size satellite galaxies are discovered in future, one has to repeat the above procedure by adjusting the lower bound on the satellite mass and scan model parameters on a model-by-model basis again. This drives us to use the $\{\alpha,\beta,\gamma\}$ parametrization. Once the observational constraint on $\{\alpha,\beta,\gamma\}$ is updated, one can easily update the constraint on models parameters by using a constructed map between model parameters and $\{\alpha,\beta,\gamma\}$.

Lyman-$\alpha$ forest

Another observable is the Lyman-$\alpha$ forest in high-resolution quasar spectra. The flux power spectrum is a powerful probe of underlying galactic-scale structure, while the thermal history of the intergalactic medium has uncertainties. The most stringent constraint seems to exclude the WDM solution to small-scale issues [187].

The procedure for the Lyman-$\alpha$ forest constraint is an example of mapping the reported lower bound on the conventional thermal WDM mass onto a given model. We evaluate the impact of the WDM model on the Lyman-$\alpha$ forest data as follows. This approach follows Ref. [88], which extends the approach of Ref. [180]. First, given a 3-dimensional linear matter power spectrum $P(k)$, we calculate the 1-dimensional power spectrum as

[TABLE]

Second, we normalize the 1D power spectrum by that in the CDM model:

[TABLE]

Third, we integrate $r(k)$ over the typical range of $k$ that a given Lyman-$\alpha$ forest spectrum probes:

[TABLE]

The dimensionless deviation of $A$ represents net suppression in the Lyman-$\alpha$ forest spectrum:

[TABLE]

Finally, we compare $\delta A$ between our FIMP model and the conventional thermal WDM model with the reported lower bound on $m_{\rm WDM}$. Note that one should use the typical range of $k$ for $A$ and the lower bound on $m_{\rm WDM}$ consistently from the same dataset or analysis. If $\delta A>\delta A_{\rm WDM}$, then we regard our FIMP model is excluded.

Ref. [88] suggests $k_{\rm min}=0.5\,h/{\rm Mpc}$ and $k_{\rm max}=20\,h/{\rm Mpc}$ for the MIKE/HIRES+XQ-100 combined dataset used in Ref. [126]. The dataset places the lower bound of $m_{\rm WDM}>3.5$ keV in the conventional thermal WDM model. We find that $\delta A_{\rm WDM}=0.46$ for $m_{\rm WDM}=3.5$ keV,♢♢\diamondsuit9♢♢\diamondsuit9$\diamondsuit$9 We again remark that we do not use the fitting function given by Eq. (2.12), but directly compute the linear matter power spectrum by using CLASS [110, 111]. This may be partially why our $\delta A_{\rm WDM}=0.46$ is different from $\delta A_{\rm WDM}=0.38$ in Ref. [88].

so we use $\delta A_{3.5\,{\rm keV}}\equiv 0.46$ as an upper bound of $\delta A$ of a given model. As we see, $\delta A$ needs a data-dependent input $k_{\rm min}$ and $k_{\rm max}$ and thus one has to repeat the procedure for different dataset. A more extendable procedure is presumable. Our proposal is the $\{\alpha,\beta,\gamma\}$ parameterization. For a given new dataset, while one has to update constraints in terms of in terms of $\{\alpha,\beta,\gamma\}$, one can use the constructed map between model parameters and $\{\alpha,\beta,\gamma\}$ as it is.

2.3 $\{\alpha,\beta,\gamma\}$ parametrization of the transfer function

As we described above, the thermal WDM model has been conventionally used to report observation constraints on the transfer function $T^{2}(k)$. The single-parameter fitting function of $T^{2}(k)$ in the thermal WDM model is given by [57, 188, 119]♢♢\diamondsuit10♢♢\diamondsuit10$\diamondsuit$10 We refer readers to Ref. [102] for a fitting function of $T^{2}(k)$ in the resonantly produced sterile neutrino DM.

[TABLE]

Here $\nu=1.12$ and thus only $\alpha$ is a parameter related with the thermal WDM mass:

[TABLE]

from Ref. [119].

However, the single-parameter ($m_{\mathrm{WDM}}$) characterization does not cover a wide range of beyond-WIMP models. Ref. [88] proposes characterizing the transfer function as

[TABLE]

This parametrization allow us to divide the procedure to place constraints on FIMPs into two with $\{\alpha,\beta,\gamma\}$ being a “common language”. On one side, one calculates $\{\alpha,\beta,\gamma\}$ as a function of model parameters in a given model (corresponding to the left red flow in Fig. 1). On the other side, one reports a likelihood function from observations as a function of $\{\alpha,\beta,\gamma\}$ (corresponding to the right red flow in Fig. 1). By combining these two, one can obtain the constraints on model parameters more easily. This procedure is also very extendable. Once a new observation date becomes available, what one has to do is just to update the latter, namely, constraints on $\{\alpha,\beta,\gamma\}$. One does not need to repeat the former. One can use a constructed map between model parameters and $\{\alpha,\beta,\gamma\}$ as it is.

A remaining challenge is how to share results from the two sides. It is not apparent how to share 3-parameter results efficiently. In this paper, we propose using neural network technique. In the context of the paper, advantages of using a neural network are:

• –

It expresses nonlinear relations quite efficiently.

• –

It learns nonlinearity without being explicitly taught.

• –

It provides us with a unified format in presenting results.

We indeed see these advantages in Sec. 4.

3 Simplified FIMP model

In this work, we consider a simple setup. The model contains a seemingly renormalizable interaction of Majorana DM $\chi$ with a heavy Dirac fermion $\Psi$ and a heavy scalar $\phi$:

[TABLE]

with the Yukawa coupling $y_{\chi}$. We assume the mass hierarchy of $m_{\Psi}>m_{\phi}\gg m_{\chi}$.♢♢\diamondsuit11♢♢\diamondsuit11$\diamondsuit$11 The result will change only slightly for $m_{\phi}>m_{\Psi}\gg m_{\chi}$ and for different quantum statistics of particles [99].

This simplified model virtually corresponds to a light axino FIMP model considered in Refs. [189, 99]. The axino FIMP model is based on a supersymmetric version of Dine-Fischler-Srednicki-Zhitnitsky axion model [190, 191]. Axino is a fermionic supersymmetric partner of axion that dynamically explains why the strong interaction preserves $CP$ very precizely [192, 193, 194, 195]. One can identify $\chi$, $\Psi$, and $\phi$ as light axino, Higgsino (supersymmetric partner of Higgs), and Higgs in the axino FIMP model.

3.1 Freeze-in production

We assume that $\Psi$ is equilibrated in thermal plasma. Freeze-in production of DM $\chi$ proceeds mainly through $2$-body decay of $\Psi\to\phi+\chi$. The production process ceases (decouples) when the plasma temperature $T$ gets comparable with the mother particle mass; i.e., the decoupling temperature is $T_{\rm dec}\sim m_{\Psi}$.

It is convenient to define a DM “temperature” as

[TABLE]

with the effective number of massless degrees of freedom $g_{*}(T)$ and the decoupling temperature $T_{\rm dec}$. This temperature scales as $T_{\rm DM}\propto 1/a$ with the cosmic scale factor $a(t)$ and thus the dimensionless momentum $q=p/T_{\rm DM}$ is conserved after the decoupling. In the following, we take $g_{*}(T_{\rm dec})=g_{*}^{\rm SM}=106.75$ (all the SM particles) as a baseline value.

Case A: Decay with entropy production

Meanwhile, we incorporate a different value of $g_{*}(T_{\rm dec})$ or entropy production after the decoupling, by introducing $\Delta$ as

[TABLE]

$\Delta>1$ takes account of entropy production after the decoupling, or lager degrees of freedom at the decoupling (e.g., minimal supersymmetric standard model, where $g_{*}^{\rm MSSM}=226.75$). $\Delta<1$ is applied to the case of late decoupling, i.e., $g_{*}(T_{\rm dec})<g_{*}^{\rm SM}$.

We take into account only relevant model parameters to “warmness” of FIMP DM. Note that warmness of FIMP DM depends on the phase space distribution $f(p)$ (equivalently, $f(q)$ and $T_{\rm DM}$) and the FIMP mass $m_{\chi}$. The phase space distribution does not depend on an absolute scale of $m_{\phi}$ and $m_{\Psi}$, but is sensitive to the ratio $m_{\phi}/m_{\Psi}$ since the ratio determines the kinematic phase space of decay product, i.e., $\chi$ in this case. If the two masses are degenerate, the energy of $\chi$ in $\Psi\to\phi+\chi$ is suppressed and thus the resultant $\chi$’s are colder [98, 189, 99].

In this class of models, therefore, the relevant parameters are

[TABLE]

Hereafter we use the notation of $m_{1}=m_{\Psi}$, $m_{2}=m_{\phi}$, and $m_{\mathrm{DM}}=m_{\chi}$ for the sake of notational simplicity. The Yukawa coupling $y_{\chi}$ is fixed by the observed DM abundance $\Omega_{\mathrm{DM}}=m_{\chi}s_{0}Y_{\chi}/\rho_{c}$. While the colder phase space distribution is realized for a more degenerate mass spectrum, the larger Yukawa coupling or lighter $\Psi$ is necessary to obtain the observed DM abundance.

Case B: Decay with scattering

Generally a daughter particle $\phi$ has another interaction with a light Dirac fermion $f$:

[TABLE]

with the Yukawa coupling $y_{f}$. One can identify $f$ as top quark (again $\phi$ as Higgs) in the axino FIMP model [189, 99]. We assume the mass hierarchy of $m_{\Psi}>m_{\phi}\gg m_{f}$. We also assume that $f$ is equilibrated in thermal plasma. In this case, freeze-in production of $\chi$ occurs through $s$-channel scattering of $f{\bar{f}}\to\Psi\chi$ and $t$-channel scattering of $\Psi f\to\chi f$ as well as through 2-body decay of $\Psi\to\phi+\chi$. The decoupling temperature is again $T_{\rm dec}\sim m_{\Psi}$. $y_{f}$ determines the scattering contribution to the yield, $Y_{\rm scat}$. Freeze-in production through scattering becomes more important for more degenerate $\phi$ and $\Psi$, since the partial decay width becomes smaller.

In summary, in this case, the relevant parameters are

[TABLE]

Again hereafter we use the notation of $m_{1}=m_{\Psi}$, $m_{2}=m_{\phi}$, and $m_{\mathrm{DM}}=m_{\chi}$ for the sake of notational simplicity. $y_{\chi}$ is fixed by the observed DM abundance: $Y_{\rm dec}+Y_{\rm scat}=Y_{\rm total}$. In this case, we do not vary $\Delta$ but take several values such as $\Delta=0.3,1$, and $3$.

3.2 Constraints

We derive constraints from $N_{\rm sat}$ and from $\delta A$ through the conventional procedure described in Sec. 2 (corresponding to the blue flow in Fig. 1).

First we present constrains from $N_{\rm sat}$ in Fig. 2. The top-left panel is for Case A (Decay with entropy production), while the other panels are for Case B (Decay with scattering). For Case A, bluer regions satisfy the condition $N_{\rm sat}>N_{\rm sat}^{\rm obs}=63$ for each value of $\Delta$. For Case B, the three panels correspond to $\Delta=0.3$ (top-right), $\Delta=1$ (bottom-left), and $\Delta=3$ (bottom-right), respectively. As in Case A, bluer regions satisfy $N_{\rm sat}>N_{\rm sat}^{\rm obs}$ for each value of $m_{\rm DM}$. We also display two lines corresponding to $y_{f}=\sqrt{\pi/3}$ (red-dashed) and $y_{f}=\sqrt{1/3}$ (red-dotted), to depict a perturbative Unitarity limit.

Next we show constraints from $\delta A$ in Fig. 3. The four panels are for Case A (top-left) and for Case B with $\Delta=0.3$ (top-right), $\Delta=1$ (bottom-left), and $\Delta=3$ (bottom-right), respectively. For each parameter, bluer regions satisfy the condition $\delta A<\delta A_{3.5\,{\rm keV}}$. The red lines are the same as Fig. 2. We see that $\delta A$ gives stronger constraints than $N_{\rm sat}$.

As repeatedly stated, constraints on the transfer function are often provided in terms of the conventional thermal WDM mass $m_{\rm WDM}$. In Appendix A we convert $m_{\rm WDM}>2.9$ keV corresponding to $N_{\rm sat}>N_{\rm sat}^{\rm obs}$ and $m_{\rm WDM}>3.5$ keV corresponding to $\delta A<\delta A_{3.5\,{\rm keV}}$ into constraints on our FIMP parameters. We see that the constraints are qualitatively similar but quantitatively slightly different ($\sim 10\%$ in $m_{\rm DM}$) from those derived in this section.

4 Neural network approach

As stressed in Sec. 1, one of the main purposes of this paper is to provide ready-to-use maps for “Model parameters $\to$ $\{\alpha,\beta,\gamma\}$” and also for “$\{\alpha,\beta,\gamma\}$ $\to$ Observables” (see the red flow in Fig. 1). Our proposal is to use a neural network for this purpose. In the following we first explain our neural network setup in Sec. 4.1, and then construct concrete neural networks for “Model parameters $\to$ $\{\alpha,\beta,\gamma\}$” and “$\{\alpha,\beta,\gamma\}$ $\to$ Observables” in Sec. 4.2 and Sec. 4.3, respectively. Finally we combine the two neural networks to reproduce the constraints presented in Sec. 3 to demonstrate the precision of the neural networks.

4.1 Neural network setup

The setup of our neural network is summarized in Fig. 4. We identify the input vector $\vec{x}_{\rm in}$ as the three model parameters for each of Case A and B. As the layer proceeds, the original layer is operated by linear algebra and then multiplied by a nonlinear function $\vec{f}$. More concretely, the connections among the layers are given by

[TABLE]

where $N$ is the number of hidden layers and $W$’s and $b$’s are called weight matrices and biases, respectively. The nonlinear function $\vec{f}$ is understood as acting on each component:

[TABLE]

and we adopt a Rectified Linear Unit (ReLU) [197] for the $f$ function:

[TABLE]

We train the neural network with supervised learning. As we explain in the next subsections, we collect $\mathop{\mathcal{O}}(10,000)$ combinations of the input $\vec{x}_{\rm in}$ and the true value (from direct calculations) of the output $\vec{x}_{\rm out}^{\rm(true)}$. Note that, with such a large number of data points, it is much more efficient to recast the obtained data onto the neural network and share the neural network parameters than to provide the data itself. Training of the neural network is performed through the updates of the weight matrices and biases so that the output of the neural network $\vec{x}_{\rm out}$ gets closer to the true value $\vec{x}_{\rm out}^{\rm(true)}$. The closeness is measured by the loss function $E$, which we take as

[TABLE]

where $(\vec{x})_{i}$ denotes the $i$-th component of $\vec{x}$.

For the number of hidden layers, we use $N=2$ in this paper. Then the relation between the input and output reduces to

[TABLE]

We construct the neural network using the public code TensorFlow [198],♢♢\diamondsuit12♢♢\diamondsuit12$\diamondsuit$12 We use the version r1.1.7.

and train it for $\mathop{\mathcal{O}}(10^{5})$ epochs. The whole dataset is split into training ($90\%$) and test ($10\%$) subsets, and the former is used to train the neural network, while the latter is used to monitor the training process and avoid possible overfitting. We also apply a 10% dropout [199] to avoid overfitting. We use Adam Optimizer [200] with a learning rate of 0.001.

4.2 Model parameters $\to$ $\{\alpha,\beta,\gamma\}$

We first construct a neural network connecting model parameters and the transfer function parameters $\{\alpha,\beta,\gamma\}$. Before moving on, however, we remark that parameter degeneracy often appears when we fit $\{\alpha,\beta,\gamma\}$ to the resulting power spectrum in the simplified FIMP model in Sec. 3. Indeed Ref. [88] also notices this parameter degeneracy (see Appendix. A of Ref. [88]). Meanwhile, Ref. [89] reports that the combination of $|\beta\times\gamma|$ is well constrained by observational data (see Fig. 4 of Ref. [89]), while the orthogonal direction $\beta/\gamma$ is not very sensitive. Therefore, in this paper, we fix this orthogonal direction by the relation

[TABLE]

As a result, the output $\vec{x}_{\rm out}$ becomes a two-component vector.

Case A: Decay with entropy production

Let us first take Case A (see Sec. 3.1). We identify the output $\vec{x}_{\rm out}$ and input $\vec{x}_{\rm in}$ as

[TABLE]

Here $\vec{x}_{{\rm in},0}$ and $\vec{x}_{{\rm out},0}$ are the means of the input and output data, respectively, while $\vec{\sigma}_{\rm in}$ and $\vec{\sigma}_{\rm out}$ are the standard deviations. These are constant vectors introduced to normalize the data and make learning more efficient.

For the dataset, we sample about $20,000$ data points from $0<m_{2}/m_{1}<1$, $1\leq m_{\rm DM}\,{\rm[keV]}\leq 20$, and $0.1\leq\Delta\leq 10$. We exclude data points in the gray-shaded regions of Figs. 5 and 6, and thus the resulting neural networks cannot be used for the input parameters in these regions.♢♢\diamondsuit13♢♢\diamondsuit13$\diamondsuit$13 The reason for excluding the gray-shaded regions is as follows. For Case A, the right-top corner of the parameter space corresponds to the CDM limit. Since the transfer function approaches unity in this region, the parameter set $\{\alpha,\beta,\gamma\}$ are not uniquely determined by fitting even after $\gamma=-\beta$ is imposed. For Case B, the left-top corner corresponds to the large Yukawa coupling limit and thus the perturbative Unitarity violation problem arises.

Case B: Decay with scattering

Next let us take Case B (see Sec. 3.1). We identify the output $\vec{x}_{\rm out}$ and input $\vec{x}_{\rm in}$ as

[TABLE]

For the dataset, we sample about $20,000$ data points from $0.5<m_{2}/m_{1}<0.9$, $0\leq Y_{\rm scat}/Y_{\rm total}\leq 0.7$, and $1\leq m_{\rm DM}\,{\rm[keV]}\leq 15$. We again exclude data points in the gray-shaded regions of Figs. 5 and 6.

4.3 $\{\alpha,\beta,\gamma\}$ $\to$ Observables

We next construct a neural network that maps $\{\alpha,\beta,\gamma\}$ onto the observables, more specifically, $N_{\rm sat}$ and $\delta A$ introduced in Sec. 2. We identify the input and output as

[TABLE]

[TABLE]

Note that we do not assume $\gamma=-\beta$ in contrast to the previous subsection, and thus $\vec{x}_{\rm in}$ is a three-component vector. This is to accommodate broader class of models than the models we adopt in this paper. Also note that $\vec{x}_{\rm out}$ is a one-component vector, which means that we construct neural networks for “$\{\alpha,\beta,\gamma\}\to N_{\rm sat}$” and for “$\{\alpha,\beta,\gamma\}\to\delta A$” separately.

For the dataset, we sample about $70,000$ points from $0.001\leq\alpha\leq 0.1$, $0.1\leq\beta\leq 10$, and $0.1\leq\gamma\leq 10$.

4.4 Combined results

Before combining the two neural networks constructed in the previous subsections, we remark that we discuss details about the precision of the neural network in Appendix B. We provide the resultant neural network parameters through the arXiv website. See Appendix C for further explanation of the data files. We also provide a Mathematica file (freeze-in.nb) for illustration.

Now let us check the precision of the neural network by combining the two neural networks. The results should coincide with the constraints obtained in Sec. 3 as long as the neural networks work well. Figs. 5 and 6 are the constraints from $N_{\rm sat}>N_{\rm sat}^{\rm obs}$ and $\delta A<\delta A_{3.5\,{\rm keV}}$ derived through the combination of the two neural networks and thus should be compared with Figs. 2 and 3, respectively. We see that the neural networks nicely reproduce the original constraints.

We again stress constructing nonlinear maps for “Model parameters $\to$ Linear matter power spectrum” and for “Linear matter power spectrum$\to$ Observables” separately is very useful and time-saving: given the common language of $\{\alpha,\beta,\gamma\}$, those interested in particle physics models can provide $\{\alpha,\beta,\gamma\}$ as functions of model parameters, while those who reports observational constraints can update the constraints in terms of $\{\alpha,\beta,\gamma\}$. Neural network technique provides us with a ready-to-use format for this procedure.

5 Summary

Galactic-scale structure formation of the Universe is of particular interest in DM research. Beyond-WIMP scenarios alter galactic-scale structure formation, while conventional WIMP DM behaves as CDM. Precise measurement of galactic-scale structure in near-future observations may hint beyond-WIMP scenarios. On the other hand, here is a practical bottleneck. Impacts of beyond-WIMP scenarios on galactic-scale structure vary model by model. In principle, one has to repeat the two-step procedure on a model-by-model basis:

• Model $\to$ Linear matter power spectrum $\to$ Observables,

which is sketched by the blue flow in Fig. 1. Each step requires different disciplines and dedicated computations. Following this procedure in the model-by-model basis is very time-consuming.

We may improve the situation by characterizing the transfer function (i.e., the linear matter power spectrum) with some parameter. One (likely particle physicist) calculates the transfer function parameter as a function of model parameters. Another reports observational constraints in terms of the transfer function parameter. Now we can get constraints on the model parameters very easily by combining the two results. Although a single-parameter characterization (the thermal WDM mass $m_{\rm WDM}$) has been conventionally used, 3-parameter characterization is proposed to cover a wide range of beyond-WIMP scenarios. Our main stress is that neural network technique facilitates sharing results from one side to another by providing the results in a ready-to-use format.

We devoted this paper to demonstrating how we can actually work out with $\{\alpha,\beta,\gamma\}$ and a neural network. To be specific, we considered a simplified model of light (keV-scale) FIMP DM Freeze-in production from $2$-body decay gives a main contribution to the relic abundance. We also took into account entropy production after the decoupling and freeze-in production from scattering. We constructed first a map between the FIMP model parameters and $\{\alpha,\beta,\gamma\}$ and next a map between $\{\alpha,\beta,\gamma\}$ and the observables, i.e., the number of satellite galaxies and Lyman-$\alpha$ forest, by adopting neural network technique. We provided the constructed maps in a ready-to-use format through the arXiv website. Meanwhile, we performed the conventional procedure to derive the direct constraints on the FIMP model parameters. The constraints derived through $\{\alpha,\beta,\gamma\}$ and a neural network are in good agreement with those derived through the conventional procedure.

Although we focused on a simplified model of FIMP DM in this paper, it is worth performing a similar study in other FIMP models such as sterile neutrino DM and superWIMP DM and also in other alternatives to CDM such as Fuzzy DM and late kinetic decoupling of DM. Our suggestion will facilitate comparison between beyond-WIMP models and future updates of constraints on galactic-scale structure formation, e.g., from redshifted 21cm surveys.

Acknowledgments

The work of KJB, RJ, and AK was supported by IBS under the project code, IBS-R018-D1. The work of RJ was supported by Grants-in-Aid for JSPS Overseas Research Fellow (No. 201960698). The work of RJ was supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC 2121 ,,Quantum Universe“ – 390833306. The work of KY was supported by JSPS KAKENHI Grant Number JP18J10202.

Appendix A Comparison with an analytic map

In this appendix we derive constraints on our FIMP model parameters by converting the thermal WDM mass $m_{\rm WDM}$. Proposed ways of converting $m_{\rm WDM}$ onto a given model are as follows:

• –

One compares the characteristic quantity such as the free-streaming length [201, 202] and Jeans length [203, 204, 205] between a given WDM model and the conventional thermal WDM model. If the free-streaming length in the given model is larger than that in the conventional thermal WDM model with an observational lower bound on $m_{\rm WDM}$, the given model is regarded as disfavored by the same observation. See Ref. [60] for comparison of the transfer function in different WDM models with the Jeans length fixed.

• –

One compares the transfer function $T^{2}(k)$ below some critical wavenumber between a given model and the conventional thermal WDM model. If $T^{2}(k)$ in the given model is smaller in amplitude than that in the conventional thermal WDM model with a n observational lower bound on $m_{\rm WDM}$, the given model is regarded as disfavored by the same observation. A suggested choice of the critical wavenumber is the half mode $k_{1/2}$ where $T^{2}_{\rm WDM}(k_{1/2})=1/2$ [109].

In this appendix, we adopt a “warmness” quantity (equivalently, the Jeans length) calculated from a DM phase space distribution [60]:

[TABLE]

where $\langle{p^{2}}\rangle$ is the 2nd moment of the DM phase space distribution and thus

[TABLE]

$\tilde{\sigma}$ depends on the shape of the phase space distribution. For a given observational lower bound on $m_{\mathrm{WDM}}$, a WDM model is regarded as disfavored by the same observation, if $\sigma>\sigma_{m_{\mathrm{WDM}}}$. Using the definition of DM temperature given by Eq. (3.2), we obtain the constraint on a FIMP as

[TABLE]

Note that in the conventional thermal WDM model, WDM particles follows the Fermi-Dirac distribution, and thus $\tilde{\sigma}_{\mathrm{WDM}}\simeq 3.6$.

In our simplified FIMP model, the phase space distribution can be expressed analytically [206], and thus $\tilde{\sigma}$ is also analytically derivable. As a result, we can construct an analytic map between $m_{\rm WDM}$ and the model parameters. The total $\tilde{\sigma}$ is calculated from each production process as

[TABLE]

where each $\tilde{\sigma}^{2}$ is calculated analytically as

[TABLE]

and each FIMP yield is also obtained as

[TABLE]

Here $r\equiv m_{2}/m_{1}$, prefactors count a number of particle spieces ($\Psi/{\bar{\Psi}}$ and $f/\bar{f}$), and $M_{0}$ is a dimensionful constant whose expression is not relevant in this appendix.

In this way, we derive the constraints on our simplified FIMP model from $m_{\rm WDM}$ through warmness. First, Fig. 8 shows constraints from the observed number of Milky Way satellites, $\sigma>\sigma_{2.9\,{\rm keV}}$. $m_{\rm WDM}>2.9$ keV corresponds to $N_{\rm sat}>N_{\rm sat}^{\rm obs}$. The left and right panels are for Case A (Decay with entropy production) and Case B (Decay with scattering) with $\Delta=1$ and should be compared with the top-left and top-right panels of Fig. 2, respectively. We see the results are $\sim 10\%$ different with each other, while are qualitatively equivalent.

Next, Fig. 8 shows constraints from the Lyman-$\alpha$ forest data, $\sigma>\sigma_{3.5\,{\rm keV}}$. $m_{\rm WDM}>3.5$ keV corresponds to $\delta A<\delta A_{3.5\,{\rm keV}}$. The left and right panels are for Case A and Case B with $\Delta=1$ and should be compared with the top-left and top-right panels of Fig. 3, respectively. We see the derived results are $\sim 10\%$ different from those from the direct modeling in Sec. 3.2, as for $N_{\rm sat}>N_{\rm sat}^{\rm obs}$.♢♢\diamondsuit14♢♢\diamondsuit14$\diamondsuit$14 We also derive the constraints through $k_{1/2}$. The derived constraints are again $\sim 10\%$ different from those from the direct modeling.

Appendix B Further check: original data vs neural network

In this appendix, we take a closer look at the difference between the original data and neural network.

First we check validity of the $\{\alpha,\beta,\gamma\}$ parametrization itself (which is irrelevant to the precision of the neural network). Fig. 9 shows the transfer function for Case A with $m_{2}/m_{1}=0.499$, $m_{\rm DM}=10.256$ keV, and $\Delta=1$ (left panel) and Case B with $m_{2}/m_{1}=0.701$, $Y_{\rm scat}/Y_{\rm tot}=0.367$, $m_{\rm DM}=4$ keV, and $\Delta=1$ (right panel). The red points are data points, while the blue lines are $T^{2}(k)$ given by Eq. (2.14) with the fitted values of $\{\alpha,\beta,\gamma\}$ ($\gamma=-\beta$ as explained in Sec. 4.2). We see that the $\{\alpha,\beta,\gamma\}$ parametrization nicely reproduces the original data.

Next we examine the precision of the neural network. Fig. 11 compares the original values of $\alpha$ and $\beta$ (upper panels) and the fit from the neural network (lower panels) for Case A with $\Delta=1$. Similarly, Fig. 11 is for Case B with $m_{\rm DM}=4$ keV. We see that the neural network not only reproduces the original data quite well, but also somewhat smoothens artificial fluctuations in the original data.

Figs. 12 and 13 are color plots for the relative error between the original data and neural network for $\alpha$ (left columns) and $\beta$ (right columns), respectively. Fig. 12 is for Case A with $\Delta=0.1$, $1$, and $10$ from top to bottom, while Fig. 13 is for Case B with $m_{\rm DM}=2$, $4$, and $6$ keV for $\Delta=1$ from top to bottom. The relative error for Case A is at most 1$\%$ in $\alpha$ and $\beta$, while for Case B the error is at most $2\%$ and $0.2\%$ in $\alpha$ and $\beta$, respectively.

We finally comment that the error of the neural network for “$\{\alpha,\beta,\gamma\}\to$ Observables” is much smaller than that for “Model parameters $\to\{\alpha,\beta,\gamma\}$”.

Appendix C How to use the neural network data

In this appendix we explain how to use the data provided through the arXiv website. The datafile we provide are

• •

mean.tsv, std.tsv,

• •

b1.tsv, b2.tsv, bout.tsv,

• •

w1.tsv, w2.tsv, wout.tsv.

The first items are the means, $\vec{x}_{0}\equiv(\vec{x}_{{\rm in},0}^{\rm T},\vec{x}_{{\rm out},0}^{\rm T})^{\rm T}$, and standard deviations, $\vec{\sigma}\equiv(\vec{\sigma}_{\rm in}^{\rm T},\vec{\sigma}_{\rm out}^{\rm T})^{\rm T}$, which shift and normalize the neural network input and output. The second items are the biases, $\vec{b}_{1}$, $\vec{b}_{2}$, and $\vec{b}_{\rm out}$, while the last items are the weight matrices, $W_{1}$, $W_{2}$, and $W_{\rm out}$.

The data files for “Model parameters $\to\{\alpha,\beta,\gamma\}$” are in the directory of freeze-in/CaseA for Case A, and in freeze-in/CaseB/Delta=... for Case B, respectively. The data files for “$\{\alpha,\beta,\gamma\}\to$ Observables” are in freeze-in/NSat and freeze-in/deltaA for $N_{\rm sat}$ and $\delta A$, respectively.