Dark matter freeze-in from semi-production

We study a novel dark matter production mechanism based on the freeze-in through semi-production, i.e. the inverse semi-annihilation processes. A peculiar feature of this scenario is that the production rate is suppressed by a small initial abundance of dark matter and consequently creating the observed abundance requires much larger coupling values than for the usual freeze-in. We provide a concrete example model exhibiting such production mechanism and study it in detail, extending the standard formalism to include the evolution of dark matter temperature alongside its number density and discuss the importance of this improved treatment. Finally, we confront the relic density constraint with the limits and prospects for the dark matter indirect detection searches. We show that, even if it was never in full thermal equilibrium in the early Universe, dark matter could, nevertheless, have strong enough present-day annihilation cross section to lead to observable signals.


Introduction
Although the nature of dark matter (DM) remains an open question, the measurement of its present-day energy density constitutes a robust positive signal of physics beyond the Standard Model (SM). A successful explanation of the properties of DM, thus, by necessity needs to provide a production mechanism that can fit the observed value of Ω DM h 2 = 0.12 ± 0.0012 [1]. Even though this is quite a constrictive requirement within any given particle physics model, numerous ways to explain the current abundance are known when being completely agnostic about the nature of DM. Among these of special interest are mechanisms that rely solely on interactions of the dark sector states with the SM ones from the thermal plasma and do not carry any dependence on the details of the initial conditions at the time of reheating. In particular, the thermal freeze-out from the primordial plasma [2], where the DM at some time in its evolution was in chemical equilibrium, and freeze-in [3,4], where its abundance was always much lower than in equilibrium and DM was gradually produced from negligible starting population.
While the freeze-out paradigm often leads to DM candidates that have appreciable signals in direct, indirect and collider searches and is subject to increasingly far-reaching constraints (for a review see e.g. [5]) the freeze-in naturally predicts the DM particles to be much more weakly interacting and more challenging to detect (see e.g. [6]). In particular, across the board expectation for stable DM produced through freeze-in is to provide null signal in the indirect detection searches (see [7,8] for exceptions).
The specific processes that lead to the freeze-in production vary depending on the model, but the most well-studied realizations are coming from the decays or 2 → 2 pairproduction processes from the states that are in equilibrium with the SM plasma (see again [6] for a comprehensive review). Analogous processes play a vital role also in the majority of models based on the freeze-out mechanism. Interestingly, in the latter context a different process has become of growing interest in the recent years, namely the semi-annihilation [9], i.e. the interaction of the type χχ → χφ, where χ represents the DM particle, while φ is either a SM state or a mediator. Semi-annihilations are present in a number of realistic models (see e.g. [10][11][12][13]) and not only are they subject to weakened direct detection limits, but also lead to peculiar phenomenology concerning the indirect detection [14] and the core-cusp cosmological problem [15][16][17].
In this work we study a novel freeze-in production mechanism based on the inverse semi-annihilation processes, which we will refer to as semi-production. A peculiar feature of this mechanism is that the production rate is suppressed not only by small couplings, but also by a small initial abundance of DM such that it vanishes in the limit when no DM particles are present in the plasma. Hence, creating the observed abundance requires much larger coupling values than in the usual freeze-in scenario. Consequently, even if DM was never in full thermal equilibrium (with the SM plasma or within the dark sector) in the early Universe, it can still have strong enough present-day annihilation cross section to be a source of observable signals in the indirect searches.
An interesting effect that arises in the study of the proposed mechanism is that the semi-production rate can strongly depend not only on the number density of DM particles at any given time, but also on their energy. Especially when the semi-annihilating partner is lighter than the DM particle, which causes a kinematical suppression of the interactions involving low energy tails of the momentum distribution. Therefore, in order to provide a reliable calculation of the resulting abundance we solve the coupled system of Boltzmann equations for both the number density and the temperature of DM.
This paper is organized as follows. We start in Sec. 2 with the basic illustration of the mechanism. In Sec. 3 we describe the specifics of the model chosen to demonstrate the properties of the proposed scenario. Sec. 4 discusses the formalism suitable for determining the DM relic density in this setup. Finally, in Sec. 5 we present results of the scan in the parameter space of the model and discuss phenomenology, while Sec. 6 concludes.

Illustration of the mechanism
To introduce the mechanism let us begin with a simple example. Consider a model extending the SM by a real scalar φ, which we will refer to as the mediator, and a complex scalar χ charged under a Z 3 symmetry, constituting the DM. Assume that the interaction Lagrangian is given by for now neglecting all the other possible interaction terms. Furthermore, let us assume that the mediator is tightly coupled to the SM plasma keeping it both in thermal and chemical equilibrium for all the time that number changing processes for χ are active. Finally, let us focus on the region where m φ < 3m χ , such that the φ → χχχ decay is kinematically forbidden and that λ 1 such that χ never thermalizes. In such an example toy model, the DM relic density is build up by the freeze-in mechanism due to the φχ → χχ semi-production process. However, in contrast to the usual freeze-in, the rate for this process vanishes in the limit of vanishing number density of χ. In other words, in order for it to take place there needs to be some initial population of χ's present. This could be e.g. a direct remnant of the reheating process [18], effect of the gravitational production [19,20], ultraviolet freeze-in [21,22], forbidden freeze-in [23] (see also [24]) or decay of the false vacua [25]. In fact, a general renormalizable Lagrangian that accommodates semi-annihilation also typically allows for the pair-production process (see Sec. 3). Such a process can, thus, naturally accompany semi-annihilation and explain the genesis of the initial dark matter population.
Irrespectively of the details of the mechanism giving rise to this initial population of χ, as long as it happens at some early time corresponding to the temperature of the Universe being T in , the yield Y χ (T in ) ≡ Y in χ = 0, the semi-production will be active. However, its effect is going to be suppressed by a very low number density of χ's, leading to the requirement for the coupling λ to be higher than for the usual 2 → 2 freeze-in.
The evolution of the DM yield Y χ = n χ /s is then given by the Boltzmann equation which can be written as (see Sec. 4.1 for more details): where the function γ φχ→χχ in general depends not only on the SM temperature parameter x = m χ /T , but also on the form of the distribution function of χ. In particular, making a simplifying assumption that its shape traces the equilibrium one, but allowing for a different temperature T χ , introducing x χ = m χ /T χ one has γ φχ→χχ (x, x χ ) = 1 xsH whereH ≡ H/ [1 + 1/3d(log h eff )/d(log T )], with H being the Hubble rate, h eff is the number of entropy degrees of freedom and s is the entropy density. If the function x χ (x) is known a priori, the number density can be calculated simply by integrating (2.2) giving a formal solution: (

2.4)
Inverting this relation one gets the prediction for the coupling λ that results in the observed relic abundance to be where Y 0 = 0.12(ρ c /(m χ s 0 )) with ρ c and s 0 being the critical density and entropy density today respectively. 10 -30 The value of the coupling λ (left) and the present day annihilation cross section (right) giving Ωh 2 = Ω DM h 2 through the semi-production mechanism for different choices of masses and fixed T χ /T . The solid lines give the corresponding values if the process was φφ → χχ instead.
To illustrate the difference with respect to the usual freeze-in and highlight the importance of the evolution of T χ alongside its number density in the left panel of Fig. 1 we show the value of the coupling (2.5) resulting in the observed relic abundance as a function of Y in (x in ) for three different choices of the masses. The shaded regions encompass the nearly an order of magnitude variation of the result when assuming the temperature ratio to be constant and a factor of 2 smaller and larger than T χ /T = 1. In the bottom part of the plot the solid lines show the corresponding values of λ if the process is instead φφ → χχ. The huge shift in its value highlights the effect of the suppression of the semi-production channel.
The right panel of Fig. 1 shows the value of the corresponding present-day annihilation cross section as a function of the DM mass for two choices of Y in (x in ) and for fixed m χ /m φ = 10. Even though in this simple example the resulting cross sections are out of the observational reach, it is clear that the semi-production freeze-in points to the parameter regions that have much greater potential to be probed by indirect detection observations. It is also worth mentioning that if the process was instead φφ → χχ, then, as expected, there would be no visible dependence on Y in , as can be seen from the solid lines coinciding for both choices.
In reality both the temperatures of the DM and the mediator are affected by the decay and annihilation processes in a non-trivial way, e.g. the semi-production leads to the selfcooling due to a portion of the kinetic energy being used to produce the heavier state. Therefore, in what follows we study the presented mechanism in more detail, dynamically solving for the temperature evolution. In order to do so we first define a complete example model which also includes self-consistently the initial production process, followed by the period of semi-production.

The example model
For a particular realization of the mechanism described in the previous section we consider a Higgs-portal type model with a scalar mediator φ. Such models are extensively studied in the literature and provide a vast phenomenology for collider searches and astroparticle physics (see e.g. [26] for a recent review). We take the DM candidate χ to be a scalar particle, protected by the Z 3 symmetry. In the absence of any stabilizing symmetry the φ is expected to mix with the SM Higgs and be unstable, decaying into visible states. 1 The general and renormalizable Lagrangian for such a model is where H stands here for the SM Higgs doublet. To preserve the Z 3 symmetry of the vacuum, we require that χ does not acquire a vacuum expectation value (VEV). For light φ mediator the branching ratio of the SM Higgs decay to invisible particles is constrained to be [27] smaller than 0.19, which translates to λ hφ 10 −2 . Furthermore, we will focus on the region where couplings λ hφ , λ 1 and λ 2 are extremely small, as required by the freeze-in framework. Finally, λ φ and λ χ are free as they are inconsequential for the DM production, as long as they are small enough not to lead to self-thermalization through 2 ↔ 3 or 2 ↔ 4 processes. Before the electroweak phase transition (EWPT) the mass of χ is simply determined by the parameter µ χ , while the four massless degrees of freedom of the Higgs doublet have the effective masses arising from the temperature corrections. We adopt [28] where g 1 , g 2 and y t are the gauge and top Yukawa coupling constants. The mass of the φ can in principle receive thermal corrections as well, coming from both three-and four-point self-interaction vertices, but in what follows these can be neglected as φ is significantly underpopulated compared to its equilibrium number density. Below T = T EW ∼ 160 GeV both H and φ acquire their VEVs, such that H = (0, v + h H )/ √ 2 (in the unitary gauge) and φ = v φ + ϕ. The physical scalars appear as the combinations of h H and ϕ

5)
1 As before we focus on the regime where m φ < 3mχ to avoid the φ decay to be the main production channel during the freeze-in process.
where we restore the previously used notation φ for the mass eigenstate of the new scalar.
Assuming that the minimum of the scalar potential corresponds to H † H = v 2 /2, where v = 246 GeV, and that the parameters A, µ 3 and λ φ are very small, we obtain the following minimization conditions where we have used the expression m 2 φ ≈ µ 2 φ + λ hφ v 2 /2. We are interested in the case when h is essentially the SM Higgs with the mass m h ≈ 125 GeV and φ is weakly coupled to it, so that the mixing angle θ 1. Nominally, the mass of χ also gets a correction and would become m χ = µ 2 χ + λ 2 v 2 φ /2, but as in our setup both λ 2 1 and v φ v, this shift is negligible.
At very early times, for T > T EW the system is in an unbroken phase where the production of φ is dominated by the 2 → 2 pair-production processes hh → φφ. After the EWPT, for T < T EW , the system moves to a broken phase where h → φφ opens up, as long as m φ < m h /2, and takes over the role of the most efficient production mode. 2 For a small mixing angle the mediator φ is very long-lived, such that it can be considered stable during the freeze-in period. Later it will decay with the lifetime given by the inverse of the decay width Γ φ ≈ θ 2 Γ h→SM (m φ ), where Γ h→SM (m φ ) is the total width of the SMlike Higgs boson with mass m φ . It follows that the mediator is extremely long-lived at low mass, where only its decays into lepton pairs and photons are kinematically allowed. As we will discuss in the Sec. 5, such long lifetime is severely constrained by the astrophysical limits and the beam dump limits. Therefore, we will restrict ourselves to m φ > 100 MeV in the following.

Relic density calculation
In this section we introduce the formalism we adopt for determining the amount of produced DM. As we will show, in order to determine the DM relic density in this model one needs to trace the number densities, as well as temperatures, of both φ and χ. 3 This in general leads to a large set of coupled differential equations that, though tractable, are numerically rather expensive. However, in our setup, due to a very weak coupling of φ and χ one can simplify the calculation by neglecting the backreaction of the χ evolution on the temperature of φ, which in turn is then determined solely by the coupling to the SM plasma. Consequently, one can first solve for (n φ , T φ ) and then treat the resulting T φ (T ) as an input for the remaining evolution of (n χ , n φ , T χ ). 2 During the EWPT itself there is a short period when the Higgs bosons can convert to φ through oscillations [29], which however does not give rise to an appreciable contribution in the parameter regions we study here. 3 For concreteness we will assume that all the interactions conserve CP and that there is no asymmetry between χ and χ * , while relaxing this assumption may lead to interesting phenomenological consequences.

Boltzmann equations
The Boltzmann equation for the evolution of the distribution function f i={χ,φ} has the form where the collision term C [f i ] is described in the next section. Typically the only relevant quantity in the freeze-in scenarios is the produced number density , since the dark sector is never populated densely enough for the backreaction processes, nor for the actual momentum distribution of DM, to have any impact. Therefore, only the 0-th moment of the equation (4.1) is usually considered (but see [30] for an exception). In the case of semi-production, however, χ is also present in the initial state, leading to a dependence of the production rate on its distribution function.
The most direct approach to tackle this problem is to solve Eq. (4.1) completely numerically. This has been done recently in the context of DM freeze-out in [31,32], where it has been also shown that solving the system of equations for coupled 0-th and 2-nd moment (cBE) of (4.1) instead often gives from fairly to very good estimate of the resulting relic density. 4 At the same time the structure of the semi-annihilation and semi-production rates is technically more challenging due to the necessity of performing thermal averages on pairs of particles with different temperatures [14,15,33]. Having all this in mind we therefore leave full numerical solution of (4.1) to future work, while here we focus on the system of cBE.
The derivation of the system of cBE follows exactly [31] (see also [34] and [35] for a different formulation). Integrating (4.1) over g i d 3 p/(2π) 3 /E i and g i d 3 p/(2π) 3 p 2 /E 2 i respectively and introducing x = m χ /T , the yield Y i ≡ n i /s and the 'temperature' parameter y i ≡ m i T i /s 2/3 , one gets Here we also introduce the first two non-vanishing moments of the collision term As it can be seen from (4.3), however, additional assumptions are needed to close the Boltzmann hierarchy. In particular, it contains a higher moment for which, as in [31], we will make an ansatz for both i = {χ, φ} that where µ i is the chemical potential. This ansatz would be exact in the limit of very efficient self-interactions keeping the shape of the distribution functions close to thermal, albeit with a different temperature. Note that we approximate the equilibrium distribution by the Maxwell-Boltzmann one even though the freeze-in occurs dominantly when χ and φ are relativistic. This is justified as both are extremely diluted throughout the whole evolution. 5 However, in the numerical computations at high temperatures we treat the SM Higgs in equilibrium as having relativistic distribution, cf. [36].

The structure of the collision term
The collision term C[f i ] can be expressed as a sum of integrals that correspond to different processes in which the particle i participates, i.e. φφ → χχ pair-production, φχ → χχ semi-production, production from the interactions with the SM particles (mainly the Higgs boson) and elastic scatterings: In principle, each of these terms is different for χ than for φ, although the structure of the expressions remains similar, while numerical factors and signs can vary. For example, the general expression for the semi-annihilation collision term for the evolution of χ's distribution function is The corresponding term for the evolution of f φ has the same structure, but a reverse order of momenta, a different sign and an additional symmetry factor of 1/2. Assuming that the populations of χ and φ are very diluted in comparison to their equilibrium values, we can safely neglect the backreaction of φ on the density of SM plasma and omit all of the contributions to the collision term that are O(f 2 χ,φ ). Thus, for calculation purposes is convenient to further rewrite all the expressions like Eq. (4.9) in terms of specific combinations of distribution functions. For the system of cBE we do not require the full unintegrated collision term as above, but rather only its 0-th and 2-nd moments. The 0-th moment, Eq. (4.4), governs the rate of the number density evolution and is basically an integral of the collision term over the phase-space of the particle under consideration. The elastic scattering part contribution to the collision term vanishes after this integration and does not affect the density evolution, as expected. Considering Eq. (4.9) as an example and using the ansatz from Eq. (4.7) and the approach described in the previous paragraph, one can notice that the 0-th moment can be formulated simply in terms of the velocity-averaged cross sections for various 2 → 2 processes (or in terms of the width for decays), e.g.
The 2-nd moment, Eq. (4.5), governs the rate of the temperature evolution and is considerably more complicated. The additional p 2 /E factor makes the integration more involved and the 2-nd moment terms cannot be in general expressed as simply as in Eq. (4.10). In particular in the case of semi-annihilations such terms require additional numerical integrations and hence the solution of the cBE in full generality takes noticeably more time. Similarly, the 2-nd moment of the elastic scatterings term also connects particles with different temperatures leading to the same complications. Note, that unlike in the freeze-out case here one cannot justify the expansion in the Fokker-Planck type collision term [37,38]. However, we have checked with explicit integration of the scattering term that for the interaction strengths considered, the momentum transfer rate is too small to have any impact.
Therefore, after including all the decay, pair-annihilation and semi-annihilation processes, in the absence of an asymmetry between χ and χ * one arrives at the system of cBE that we solve numerically: where Y DM = Y χ + Y χ * and the complete set of expressions entering the moments of the collision term can be found in the Appendix A.
As mentioned above, due to a very weak coupling between the states in the dark sector, for the temperature evolution of φ we can consider only the interactions with the SM plasma, while neglecting the backreaction of χ. That is, we first solve Eq. (4.11) coupled only with (4.12) to determine the relation T φ (T ). Next, we use it as an input for the system of coupled Eqns. (4.11), (4.13) and (4.14).

Results for benchmark points
In the next section we present the results for our model's parameter space scan, but first let us start with a discussion of the thermal history of χ and φ using as examples two representative benchmark points. These are shown in Figs. 2 and 3 for two choices of the parameters that lead to the DM creation mode, which is dominated by semi-production and pair-production respectively. In both figures the left panel gives the evolution of the yields Y χ and Y φ , while the right panel shows the departure of the y parameters from the  Figure 2. The evolution of χ and φ number densities (left) and temperatures (right) for a benchmark point with m χ = 100 GeV, µ φ = 1 GeV, λ 1 = 1.1 × 10 −2 , λ 2 = 10 −8 , λ hφ = 6 × 10 −11 and θ = 10 −5 , which makes semi-production the dominant process. Left: the solid lines show the evolution of Y φ (black), Y χ (orange) and Y χ (T χ = T ) (green), while the dotted correspond to their equilibrium values. The process specified below the orange line highlights the dominant production channel in a given regime. Right: in the main plot the solid lines show the evolution of y φ (black) and y χ (orange), while in the inset plot their ratio to the equilibrium value is shown. The part of the orange line with lighter shading shows the regime where no significant population has yet been produced. In all of the plots the light gray shaded region depicts the time window in which the production of χ is most efficient.
corresponding equilibrium ones illustrating the magnitude of the departure from kinetic equilibrium with the SM plasma.
Looking at the left panels of the figures one can distinguish four phases of the evolution: I) At very high temperatures one is in an unbroken SU(2) phase and φ is being produced from 2 → 2 processes, while the production of χ, being dependent upon Y φ , is still negligible.
II) After some non-negligible population of φ is built up, the φφ → χχ process starts to slowly create DM particles. This phase is essential to give the initial seed for the semi-production to occur.
III) After the EWPT the h → φφ decay is switched on, which leads to an accelerated production of φ and indirectly χ as well. The growing populations of both particles exponentially boost the increase of Y χ via semi-production, given that the corresponding coupling is considerably large.  Figure 3. The same as in Fig. 2, but for m χ = 100 GeV, µ φ = 1 GeV, λ 1 = 10 −10 , λ 2 = 10 −3 , λ hφ = 2 × 10 −11 and θ = 10 −5 , which makes pair-production the dominant process.

IV) When the temperature drops further all the interaction processes freeze in and both
number densities reach a plateau. Not shown in the plot, but worth mentioning is that later on φ will eventually decay through the mixing with the SM Higgs and thus will not contribute to the DM relic abundance.
Turning now to the right panels of Figs. 2 and 3 showing the temperature change one can observe that for most of the evolution, both φ and χ have temperatures of the same order of magnitude as T . This is of course expected, as at early times the 2 → 2 processes tend to produce φ and χ with energies of order ∼ T and since T m χ , m φ , they redshift in the same way as the SM plasma. Note that in reality the T φ and T χ are somewhat lower than T , which is related to the fact that in the relativistic regime particles in the low energy part of the distribution have a higher probability of interacting. After the EWPT point is reached, the h → φφ introduces a visible drop in the φ temperature, and then indirectly T χ as well, followed by a slow rise later on. This is because at the onset of EWPT, when h is still relativistic, the decay leads to E φ ∼ T /2, while later the same process gives T . Finally, for the evolution of T χ at x ∼ few the production freezes-in and χ completely decouples from the rest of the states, both chemically and kinetically.
The light gray region on both plots is shown to guide the eye to the region of x in which the most of the DM production takes place, to help visualise how much the actual temperatures are different then from the assumption of kinetic equilibrium. The part of the orange line with lighter shading shows the regime where no significant population has yet been built up and the value of y there is dependent on the initial condition, which later becomes washed out after more DM is produced.
Additionally, the green solid line shows the result for Y χ under the simplifying assumption that both χ and φ remain in kinetic equilibrium throughout the whole production period. For both benchmark points this overshoots the more accurate result based on solving the cBE system from factor of a few, when pair-production dominates, to more than on order of magnitude, when it is the semi-production that is more efficient. Two effects contribute to such a significant overestimate when tracing only the number density. Firstly, neglecting the change in the temperature of φ leads to an incorrect thermal average for both the φφ → χχ and χφ → χχ processes, impact of which is the greater the smaller m φ is compared to m χ . In particular, since the actual T φ is smaller than T , this thermal average overestimates the production rate. Secondly, with strong semi-production comes the dependence on T χ as well. Together with an exponential growth period it leads to a much larger change with respect to the naive expectation.
Therefore, one can appreciate the necessity of solving for not only the number density, but also for the temperature in order to arrive at an accurate estimate of the relic abundance. It should be mentioned, however, that from the phenomenological perspective, at least for the model at hand, this exponential dependence means only logarithmic change of the coupling λ 1 is sufficient to mitigate the error made by assuming the kinetic equilibrium. Nevertheless, it is worth stressing that the departure from the kinetic equilibrium can influence the production process even more strongly if, unlike in the case of our example model, the interactions are substantially velocity-dependent.

Scan results and phenomenology
When φ is in the chemical equilibrium with the SM plasma the semi-production mechanism leads to the enhanced values of the coupling compared to the typical freeze-in, but as illustrated in Sec. 2, the resulting cross sections are still somewhat below the sensitivity of the near-future indirect searches. However, combined with the idea of the sequential freeze-in [30], i.e. simply when φ is coupled very weakly to the SM states such that it also undergoes the freeze-in, the resulting signals can be much stronger.
In order to study the phenomenological implications of the proposed mechanism we performed a numerical scan within the specified model to find the points that satisfy the relic density constraint. We scanned over the following ranges with the logarithmic sampling for the five parameters: m χ ∈ 10 −1 , 10 3 GeV, m φ ∈ 10 −1 , 50 GeV, λ 1 ∈ 10 −10 , 1 , λ 2 ∈ 10 −10 , 1 , λ hφ ∈ 10 −15 , 10 −8 and θ ∈ 10 −8 , 10 −2 . The ranges were chosen to cover a wide range of masses for χ and φ, with the conditions that m χ > m φ and m φ < m h /2. As for the couplings, we restricted λ hφ to be concentrated in the freeze-in region of φ and the mixing angle to be small enough (θ 10 −2 ) not to be in conflict with the collider measurements [40].
The results of the scan are shown in Figs. 4 and 5. First in Fig. 4 all the points found are projected onto the mass vs. lifetime of φ plane and confronted with the limits from the light boson searches. The points are represented as blue (green) diamonds for the regime of semi-annihilation domination σv χχ→χφ > σv χχ→ φφ (pair-annihilation σv χχ→χφ < σv χχ→φφ ). Additionally, some of the blue diamonds are filled, signifying the points that have the present-day semi-annihilation cross section large enough to potentially affect the core formation in the dwarf galaxies, following the results of [17]. , E949 [43] and LHCb [44,45] (gray); the limits from BBN [46] (orange) and SN1987a [47] (yellow). Also we show the prospects for the future experiments with dashed lines: FASER [48] (violet), SHiP [49] (blue) and MATHUSLA [50] (red) and LHCb-Run3 (black) adopted from [51].
The shaded regions show the recast of the existing limits, while the dashed lines show the prospects for the upcoming upgrades and planned experiments. There is a visible expected correlation between the m φ and the lifetime t φ , but apart from that, one can see that the most of the allowed parameter region in this plane is accessible for the studied model.
In Fig. 5 the same points are separated into ones with dominant semi-annihilation (left panel, blue points) and pair-annihilation (right panel, green points) and projected onto the plane of dominant cross section vs. m χ . 7 Filled diamonds indicate the points that are within the reach of the future searches for the mediator φ, cf. Fig. 4, while the empty ones are beyond the prospects. The points that are in conflict with the above measurements are discarded.
In both panels we show the constraints from the CMB measurements by the Planck 7 A small subset of the points will have both processes contributing at the same level. For these one should actually combine them both and only then compare with the constraints. However, this does not change the overall results of the scan and would unnecessarily complicate the discussion.   Fermi-LAT+MAGIC combined analysis of the gamma-ray observations of dSphs [52] (orange) are shown as shaded regions, while the dashed red line gives the projected sensitivity of the CTA for the bb channel [53]. Left: points with dominant semi-annihilation σv χχ→χφ > σv χχ→ φφ projected onto the (m χ , σv χχ→χφ ) plane. The black ellipse shows the best fit region for the GC excess recasted from [54]. The dot-dashed gray line indicates the semi-annihilation cross section below which there is no appreciable core formation in the dSph by the self-heating of DM [17]. Right: points with dominant pair-annihilation σv χχ→χφ < σv χχ→ φφ projected onto the (m χ , σv χχ→φφ ) plane. collaboration [1], recasted using [55], the limits from the joint analysis of the gammaray flux from dwarf spheroidal galaxies (dSphs) performed by Fermi-LAT and MAGIC collaborations [52] and the prospects for CTA [53]. All these limits need to be taken as an illustration of the actual constraints that by necessity apply differently to different points, e.g. due to varying m φ that changes the final photon spectra. For the purpose of demonstration we took the constraints on the DM annihilations to bb and shifted them accordingly such that the σv limit for a given value of m χ corresponds to the aforementioned limits for the DM particle that has the same mass as half of the energy of φ. We take into account the corresponding normalization for the yield and the DM number density and we assume that the difference in the gamma-ray spectra between these two cases can be neglected. For pair-annihilation E φ = m χ and for semi-annihilation E φ = (m 2 φ +3m 2 χ )/4m χ (assuming both processes happen at rest). In the broad mass range of interest φ decays predominantly to quarks through the mixing with the SM Higgs, so the bb mode is a satisfactory approximation.
In the same fashion we recast the 2σ best-fit parameter region for the explanation of the GC gamma-ray excess with a quite similar scalar-singlet model [54]. These regions are marked with black ellipses in Fig. 5. This is done in order to emphasize that the freeze-in production can indeed be used to construct models explaining the GC excess (see also [8]). Note, that although we use a rather simplified approach to obtain the best-fit region of the parameters, our estimations are in a reasonable agreement with the results of a dedicated fit to the semi-annihilation signal described in [56].
Finally, in the semi-annihilation case we show the line separating the cross sections that might lead to a sufficient heat transfer in the inner parts of the dwarf galaxies to contribute to the formation of the core [17]. Above this thin gray line, depending on the efficiency of the energy redistribution among the DM particles, one might affect the core formation. 8 From these results one can infer that the best search strategy for the model at hand is through the searches of the light mediator: a significant portion of the studied parameter space is going to be covered by the beam dump experiments. On the indirect detection frontier typically the predicted (semi-)annihilation cross sections are well below the values probed by current and near-future experiments. However, some of the model realizations do give the signals that are potentially detectable. Moreover, one should stress that even though there are not so many points found in the scan that are within the current limits, the high cross-section region requires more fine-tuning and is harder to find through a random scan. Therefore, the existence of even relatively few points like this in the scan results does prove that these regions exist even in the simple example model at hand.

Conclusions
Dark matter produced through freeze-in is an attractive alternative to its more studied thermal freeze-out cousin which is typically significantly more challenging to detect, especially in the indirect searches. In this work we introduced a novel thermal freeze-in production mechanism which we called semi-production, that leads to DM candidates having much larger coupling values than in the usual freeze-in, and consequently, present day annihilation cross section potentially giving observable signals in the indirect searches.
In order to calculate the relic density in an accurate way we employed a set of coupled Boltzmann equation for the number density and temperature of DM interacting with the out-of-equilibrium mediator. We show that neglecting the temperature evolution can lead to a result for the relic abundance that is off by more than an order of magnitude, especially in the parameter space regions with dominant semi-production. Furthermore, we have performed a numerical scan of the parameter space of an example model which shows on a quantitative level that one can indeed explain the current abundance of DM through the freeze-in mechanism and at the same time predict large enough (semi-)annihilation cross section to be phenomenologically relevant. In particular, it can reach the values needed for fitting the GC excess, as well as semi-annihilation cross sections that can help to alleviate the core-cusp problem in dSphs.
Even though most of the parameter space of the model studied here predicts the signals that are too weak to be observed even in a not-so-near future, we would like to emphasize that the proposed mechanism is quite general and can be used in other models as well. The fact that the semi-production freeze-in relies on the only prerequisite that the φχ → χχ process dominates over all other interaction channels and that it leads to a quite intriguing phenomenology constitutes a promising basis for further investigations into model building employing this mechanism.
Note added: During the final stages of completion of this paper, an article [58] appeared being first to propose in essence the same mechanism as studied here. Although, as far as the production mechanism is concerned, there are no crucial qualitative differences between our works, nevertheless, our study goes quite beyond the discussion of [58] by calculating the temperature evolution, considering the out-of-equilibrium mediator and investigating indirect detection signals for a specific model realization.

A The moments of the collision term
Here we derive the expressions for the moments of the collision terms introduced in Eqs. (4.4) and (4.5). Note that presented results are not completely general, but follow the approximations as described in Sec. 4.2.
It is convenient to express these moments for each type of particle as a sum of the integrals N i C A→B that correspond to different processes A → B, in which a particle i is involved and N denotes the order of the moment of the collision term. The particle i under consideration always corresponds to the momentum p and its position in the reaction is underlined. We are going to use the following notation dΠ (n) i = dΠ p · ... · dΠ (n) (2π) 4 δ (4) where dΠ p ≡ d 3 p (2π) 3 2Ep , n represents the number of states in the corresponding reaction A → B and the delta-function determines the 4-momentum conservation. For simplicity we use the same notation for each amplitude squared |M| 2 , but one has to note that it always corresponds to the reaction specified in the notation of the respective contribution to the collision term. Moreover, in the model that we study |M| 2 is constant for all of the processes, so in the explicit expressions we are going to put it as a factor in front of the integrals. For the distribution functions of χ and φ everywhere below we substitute the ansatz from Eq. (4.7) with the corresponding temperatures.

The 0-th moment
Here we list all the relevant contributions needed for both evolution of the abundance of χ and φ in a form suitable for numerical computations.
Pair-production. These contributions contain the same type of distribution functions among the particles in both initial and final states 0 χ C φφ→χχ * = 1 2m χ n χ dΠ (4) Thus, they can be expressed in the usual form using the thermally averaged cross section 0 χ C φφ→χχ * = 1 m χ n χ n 2 φ σv φφ→χχ * (T φ ) − n 2 χ σv χχ * →φφ (T χ ) . (A. 3) The factor 1/2 takes into account the symmetry between the identical particles from the perspective of the Boltzmann equation. 9 The pair-annihilation term for φ can be expressed through the respective term for the evolution of χ Note that the factor of 2 here compensates the symmetry factor from the respective term for χ.
Semi-production. These contributions can still be formulated in terms of the velocityaveraged cross sections for the corresponding processes, however, due to different distribution functions cannot be cast in a form of one dimensional integral. In general we have The inverse, semi-annihilation, term can be calculated exactly as in the case of pairannihilation above, while the more complicated semi-production one can be expressed as wherep(s) = s − 4m 2 χ /2 √ s is the momentum of χ in the center-of-mass (CM) frame, s max/min (E q , E r ) = m 2 φ + m 2 χ + 2E q E r ± 2qr denote the maximal and the minimal values of the s-invariant for the states with momenta q and r and E min r (q) is the energy limit such that the reaction with the given energy of E q is kinematically possible. The integral over s can be done analytically, so the computation of such terms for the processes with a constant squared amplitudes require double numerical integration. 9 The Boltzmann equation keeps track of the states with a certain momentum, while the other momenta are integrated over. Integrating over a set of n identical states yields a factor of 1/n!. For example, in the process φφ → χχ * the BE keeps track of the φ(p), while φ(k) is integrated over (χ and χ * are not identical). In the process φφ → χχ * the BE keeps track of χ(p), while both φ states are integrated over, which yields a symmetry factor of 1/2. A useful discussion on these symmetry factors can be found in the notes inside [59].

(A.7)
For φ the semi-annihilation contributions can be also expressed through the terms above Production from the SM plasma. Interactions between φ and the SM Higgs boson h give the following contributions to the 0-th moment of the collision term (A.10) The first term here can be included in full detail, however in practice it is sufficient to use a simpler, analytical expression, where one adapts the Maxwell-Boltzmann approximation for the distribution function of h (note that this process is active only after EWPT where T is of order m h or lower): where K 1 is the modified Bessel function of the second kind. The second term containing f φ (p) is more complicated where we assume for the simplicity of the expression that m h m φ . Finally, the hh pair-production gives the following contribution 0 φ C hh→φφ = 1 2m φ n φ dΠ (4) φ f eq h (q)f eq h (r) − f φ (p)f φ (k) |M| 2 , (A. 13) which again boils down to the same structure of integration as for the pair-production above. Note, however, that for the SM Higgs we do keep the Fermi-Dirac form of the distribution function. The other contributions from the Higgs-mediated interactions with the SM plasma are strongly suppressed by the small mixing angle θ.

The 2-nd moment
As we state in the text, the calculation of the 2-nd moment terms is more elaborate. If the factor p 2 /E corresponds to the momentum of a product particle, the integration over s (see Eq. (A.6)) cannot be performed analytically, so several terms require triple numerical integrations. 10 10 Including the terms with a distribution function in the integration over the final states (such as the terms of O(f 2 χ,φ )) would lead to an even larger number of numerical integrations.