Dynamics of the pseudo-FIMP in presence of a thermal dark matter

In a two-component dark matter (DM) set up, when DM 1 is in equilibrium with the thermal bath, the other DM 2 can be equilibrated just by sizeable interaction with DM 1 , even without any connection with the visible-sector particles. We show that such DM candidates (DM 2 ) have unique ‘freeze-out’ characteristics impacting the relic density, direct and collider search implications, and propose to classify them as ‘pseudo-FIMP ’ (p FIMP ). The dynamics of p FIMP is studied in a model-independent manner by solving generic coupled Boltzmann Equations (cBEQ), followed by a concrete model illustration.

Albeit tantalizing astrophysical and cosmological evidences, what constitutes the dark matter (DM) component of the universe remains a mystery with no direct experimental evidence in terrestrial experiments so far.Amongst DM geneses, thermal freeze-out of Weakly Interacting Massive Particle (WIMP) [1,2] and non-thermal freeze-in of Feebly Interacting Massive Particle (FIMP) [3] can both address the observed DM relic density Ω DM h 2 ∼ 0.12 ± .0012[4] with widely different interaction strengths with the Standard Model (SM) particles.Multipartite dark sectors offer interesting cosmological and phenomenological consequences, like modified freezeout [32], relieving direct search constraints [6], addressing small scale structure anomalies of the universe [7] etc.The distinguishability of a two-component WIMP scenario, from that of single-component one, in terms of a doublehump in the missing energy distribution at collider [8,9] or a kink in the recoil energy spectrum [10,11] of the direct search experiment have been studied.
When the interaction between a WIMP and a FIMP remains 'feeble', the FIMP continues to be out of equilibrium [12][13][14], but enhanced interaction (of weak strength) can make the FIMP equilibrate to thermal bath and undergo freeze-out with some unique features, which are essentially model-independent, but relies on the properties of its thermal DM partner.We propose to classify such a particle into a category called 'pseudo-FIMP' (pFIMP), as it is a thermal DM in spite of having feeble connection with the SM.Also, pFIMP can be realized in presence of any thermal DM like Strongly Interacting Massive Particle (SIMP), independent of the depletion mechanism.Some model-dependent studies which mimic a pFIMP-like situation, have been studied [15,17,32,34], but without elaborating the generic features that such particles possess.This letter illustrates the pFIMP characteristics in a model-independent way followed by a concrete model illustration and detection possibilities.Actually the idea of pFIMP has been conceptualised and demonstrated here for the first time, that occurs in a particular limit of the DM-DM and DM-SM interactions in a multicomponent set up.This leads to several model possibilities to be explored in the pFIMP limit as detailed in [18].
The evolution of DM number density in WIMP-FIMP framework is governed by a set of coupled Boltzmann equations (cBEQ) as shown in Eq. ( 4), in terms of yield (Y = n s ) as a function of x = µ 12 /T , where µ 12 = m 1 m 2 / (m 1 + m 2 ) denotes the reduced mass of the two DMs, and T denotes the temperature of the thermal bath.
Other notations in Eq. ( 4) are all standard and available in any DM text [6].
For solving Eq. ( 4) to obtain DM yields in a modelindependent way, we choose some benchmark values of the DM masses and ⟨σv⟩.As the temperature dependence in ⟨σv⟩ requires σ(s) (s is c.o.m energy) i.e. diagrams that contribute to annihilation, in absence of which we assume here ⟨σv⟩ ≈ σ T ; a temperature-independent threshold value representing (σv) at the freeze-out (freeze-in) temperature for WIMP/pFIMP (FIMP).Full ⟨σv⟩ is used for a model-specific analysis and the dynamical features remain the same.
Plots in fig. 1 summarise the main outcome of the model independent analysis.In figures 1a, 1b, Y 1,2 are plotted as a function of x, for different conversion rates.In figures 1c, 1d, Ω 1,2 h 2 are plotted as a function of 2 /Y eq 1 ) 2 .All the cross-sections are in the units of GeV −2 and in the following text we may omit writing it explicitly.Both the mass hierarchies (1) m 1 > m 2 (on left) and (2) m 1 < m 2 (on right) are shown.In hierarchy (1), conversion from WIMP to pFIMP is kinematically allowed, while the reverse process is Boltzmann suppressed (∼ e −2δm x , δm = m 1 − m 2 ).For hierarchy (2), its the other way round.The whole analysis can be divided into four different regions of varying conversion cross-section,  1a), as it attains a modified equilibrium distribution following Eq.( 2), after pFIMP decouples [32], When pFIMP is heavier (hierarchy 2), the conversion from WIMP to pFIMP is kinematically suppressed, and no modified freeze-out pattern for WIMP is observed (see fig. 1b).This provides another distinction of pFIMP-WIMP scenario from the WIMP-WIMP one.After FIMP equilibrates to thermal bath, larger production of WIMP from pFIMP enhances WIMP relic density, see the jump from red thick line (σ T 11→22 = 10 −24 ) to green thick line (σ T 11→22 = 10 −17 ) in figs.1a, 1b.This feature is also observed in the violet/blue/green dashed lines in figs 1c, 1d.For hierarchy (1), the departure of the WIMP from the original equilibrium distribution to modified equilibrium (Eq.2) depends on σ T 11→22 (compare green and orange thick lines in fig.1a), however the freeze-out from modified equilibrium is primarily governed by σ T 11→SM , yielding same Ω 1 h 2 for different σ T 11→22 .On the contrary, for hierarchy (2), when pFIMP is heavier, the WIMP production from pFIMP is substantial due to kinematic accessibility, as a result the WIMP freeze out is delayed, and the WIMP yield goes down with enhanced conversion rate (see fig. 1d), providing a crucial hierarchical distinction.
In region IV , with σ T 11→22 > σ T 11→SM , we see even more interesting consequences.For hierarchy (1), the WIMP relic density falls drastically with larger σ T 11→22 due to larger depletion (violet/green/blue dotted lines in fig.1c), while the pFIMP yield remains the same, and is dependent only on σ T 11→SM (horizontal brown/yellow/red dotted lines in fig.1c, or overlapping dashed lines in fig.1a for σ T 11→22 = 10 −6 , 10 −4 ).This is because pFIMP yield in hierarchy ( 1) is governed by two quantities, WIMP number density and conversion rate.Now, the WIMP number density falls with larger σ T 11→22 , while that is compensated by the enhanced production via σ T 11→22 which keeps the pFIMP yield almost constant in this region.For hierarchy (2), the situation is however different.With increasing conversion, pFIMP yield decreases significantly, while the WIMP yield remains almost constant (see fig. 1d).This is attributed again to the combined effect of decreasing pFIMP number density along with larger WIMP production from pFIMP with larger conversion rate.Since pFIMP has no SM interaction, it goes out of the original equilibrium distribution as WIMP freezes out.However, due to considerable conversion to WIMP, it achieves a modified equilibrium distribution as seen from the orange and blue dashed lines in fig.1b.
We now perform a model-independent scan of the WIMP-pFIMP parameter space, where we assume that annihilation cross-sections and masses have adequate freedom to be varied independently.In fig.2, we scan points that satisy relic density (Ω 1 h 2 + Ω 2 h 2 = 0.12 ± .0012) in σ T 11→SM − σ T 22→11 plane for hierarchy (2).WIMP mass is kept fixed at 100 GeV and the mass splitting between WIMP and pFIMP δm = |m 1 − m 2 | is varied as shown in the colour bar.When the conversion rate is small, σ T 11→SM requires to be large in order to achieve WIMP density within the limit.Understandably δm does not play a major role here, and a large range of δm is allowed (right side of fig.2).With larger conversion rate, smaller σ T 11→SM is required (left of fig.2).If we decrease σ T 11→SM keeping σ T 22→11 fixed, the WIMP tends to become overabundant, unless the conversion contribution from pFIMP to WIMP is decreased by increasing δm via Boltzmann suppression.However, for hierarchy (1), the situation is different and the whole σ T 11→SM − σ T 11→22 parameter space within δm = |m 1 −m 2 | ≲ 10 GeV becomes allowed.
In summary, the allowed points accommodate mostly small mass difference (δm ≲ 10 GeV) and the mass hierarchy between WIMP and pFIMP play a crucial role in the consequent DM phenomenology, absent in usual WIMP-WIMP set up.Here we have ensured that the conversion cross-section remains well below the bullet cluster observational data [19] that suggests the upper limit on the self-interaction cross-section of DM(s) per unit mass to be, σ/m < 0.7 cm 2 g −1 ∼ 3.2 × 10 3 GeV −3 .
In Table I, we depict a few examples of two com-   ponent DM scenarios, with their strength of DM-DM conversion, possibility of producing correct relic abundance, along with direct (indirect) detection prospects at present/future experiments.Generically, the pFIMP detection is heavily dependent on the detection prospects of the WIMP, as the pFIMP direct search or collider search signal proceeds via WIMP loop as shown in Fig. 3.The strength of the pFIMP interaction is less than the WIMP, but can be of similar order when the loop factor is compensated by the large pFIMP-WIMP interaction.
We would now like to validate our claims of pFIMP characteristics for a specific model.Here we use temperature dependent ⟨σv⟩ in terms of the model parameters.We choose the simplest model consisting of two real scalar singlets, where ϕ 1 behaves as WIMP and ϕ 2 as FIMP/(pFIMP), both of which are rendered stable under Z 2 ⊗ Z ′ 2 symmetry.The Lagrangian density is given by, ( This model in WIMP-WIMP limit has been explored in many texts including [6].Here we focus on the pFIMP regime, partly explored in [34].For that, we choose ϕ 2 to have negligible coupling with the SM, λ 2H ≈ 10 −12 .However, in presence of the WIMP-like ϕ 1 , the conversion governed by λ 12 plays a key role and with large λ 12 ∼ 1 it reproduces all the characteristics of pFIMP.Relic density as a function of λ 12 is plotted in the supplementary material, which shows similar features to that of figs.1c, 1d and validates the pFIMP characteristics.The WIMP-pFIMP limit of this model has also been verified with the code micrOMEGAs4.1 [24]. Direct search of pFIMP is possible via WIMP loop, absent a tree-level connection to SM.The Feynman graph is shown in fig. 3 interaction (λ 12 ) helps in the detection of pFIMP, the loop suppression factor ∼ 1/(16π 2 ) allows it to evade existing bound and make it accessible to future sensitivities.In fig.4, spin-independent (SI) direct detection cross-section for the WIMP (ϕ 1 ) and pFIMP (ϕ 2 ) is shown as a function of WIMP mass for relic density allowed points of the model.The sensitivity of existing direct search data have all been compared [25][26][27][28].A clear dip is seen in the resonance region m ϕ1 ∼ m h /2, where due to resonance enhancement λ 1H is relaxed and points satisfy direct search constraints.One should however remember the limitations of Lee-Weinberg mechanism in the vicinity of resonance as pointed out by [29,30].When m ϕ1 > m ϕ2 , small δm enhances the ϕ 1 ϕ 1 → ϕ 2 ϕ 2 conversion, consequently decreasing Ω ϕ1 and effective direct detection cross-section σ SI i eff = Ωi Ωtot σ SI i , where i = {ϕ 1 , ϕ 2 }.However, most of the parameter space of   3)) in WIMP-pFIMP limit, where the loop factor is approximated as ∼ (λ 1H λ 12 )/(16π 2 ).The last two columns depict the WIMP and pFIMP production cross-section at LHC via Gluon fusion yielding mono-jet+ / E T in the final state at √ s = 14 TeV.this WIMP-pFIMP limit of the two component scalar DM model is on the verge of exclusion by the LZ (2022) data, leaving the Higgs resonance, m ϕ1 ∼ m h regions and some points close to the ±1σ exclusion limit up to 500 GeV survive to be explored.Unlike WIMP-WIMP case, where two DMs may have widely different masses, pFIMP limit necessitates the mass splitting to be small and therefore serves as a guiding principle for searching both DM components.Similarly pFIMPs can also be produced in collider via one-loop WIMP graphs yielding mono-X plus missing energy signals as estimated in the Table II for some benchmark points where the model satisfies relic density and direct search bounds.In this case, detection of pFIMP seems only possible in the resonance region, owing to a resonant enhancement in the production.Interesting to note that, even in this region, the effective coupling between Higgs and pFIMP is rather small due to loop suppression, thereby making it safe from the Higgs invisible decay constraints.One can have plethora of other possibilities, like a fermion WIMP instead of a scalar WIMP, which enhances the detection possibility of a scalar pFIMP, see for example, [18].We again note here that the model chosen for illustration is not the best for detection, rather it is the simplest one.Furthermore, pFIMP having a similar mass to that of the WIMP provides a challenge in distinguishing the DM components in direct and collider searches, as the distinguishability (a kink in the recoil spectrum for direct search [10,11] or two bumps in missing energy at colliders [8,9]) often depends on the mass splitting of the DM components.However mono-X signal can be useful, as the position of the peak here primarily depends on the operator structure, angular momentum conservation etc.We are exploring such possibilities, however the two component scalar case will be difficult to disentangle.
Finally we summarise here.We have shown that the FIMP can reach equilibrium in presence of a thermal DM component and sizeable conversion cross-section.We call it pFIMP.The freeze-out of pFIMP precedes that of its thermal partner and requires a small mass splitting with it.The numerical solution of Eq. 4 presented here and approximate analytical solution provided in the supplementary material validates all the pFIMP properties in a model independent way.The pFIMP realization in presence of SIMP is also provided in the supplementary material.The freeze-out properties of pFIMP is verified for the simplest two-component DM model with two scalar singlets.The detection prospect of pFIMP via WIMP loop in future direct, indirect and collider search experiments opens up new possibilities.Such investigations in context of direct DM search for two-component WIMP-pFIMP model has been pursued in [18].
Acknowledgments:-SB and JL acknowledge the grant CRG/2019/004078 from SERB, Govt. of India.DP thanks University Grants Commission for senior research fellowship and Heptagon, IITG for useful discussions.

SUPPLEMENTAL MATERIAL
Semi-analytic solution of the cBEQ for WIMP-pFIMP scenario A two-component WIMP-pFIMP framework is governed by a coupled Boltzmann equations (cBEQ), We elaborate upon the semi analytical solution of the cBEQ as in Eq. 4 in WIMP-pFIMP limit, for both the hierarchies: • Case-I: WIMP mass > pFIMP mass (i.e.m 1 > m 2 ).
The methodology follows the procedure in [32].
Case-I (m1 > m2) In WIMP-pFIMP limit, we neglect the pFIMP production from the visible sector as the pFIMP coupling with SM is very small, compared to the conversion cross-section.So, in the following we consider, Γ SM→22 = ⟨σv⟩ SM→22 = 0.For notational ease, let us also use, 19 GeV.Recall that σ 12 = σ 21 (n eq 2 /n eq 1 ) 2 .With above, the cBEQ in terms of bath temperature (T ) becomes, ( The equilibrium yield is given by Maxwell-Boltzmann distribution, Y eq 2 2 = y eq 2 1 , Equation 5 becomes, In order to solve the cBEQ analytically, it would be convenient to divide this whole scenario in three regions on the bath temperature: Region A: T ≫ T fi , Region B: T ≃ T fi and Region C: T ≪ T fi where T fi denotes the freeze-out temperature of both species i = 1, 2. Let use consider the difference of DM yield from equilibrium yield as ∆ 1 = Y 1 − y eq 1 and ∆ 2 = Y 2 − Y eq 2 .Using these relations, Equation 6 becomes, • Region A: For T ≫ T fi , d∆i dT is negligible and Equation 7 becomes, • Region B: Let us suppose T f2 > T f1 , then both Equation 8 and Equation 9 hold at T f2 , but only Equation 8 holds at T f1 as pFIMP already freezes out at T f1 .In the vicinity of T ≃ T f , we may further assume ∆ 1 (T f1 ) = c 1 y eq 1 (T f1 ) and ∆ 2 (T f2 ) = c 2 Y eq 2 (T f2 ) [31], where c i 's are unknown constants whose values are determined by matching the analytical result with the numerical one.Substituting in Equation 8and Equation 9, we get, In the above, Equation 10 and Equation 11 can be written further as, where, Equation 17 and Equation 18 are transcendental equations, analytical solutions are difficult to obtain, but using numerical method it is possible to extract the values of freeze out temperature T f1 and T f2 , which we did.
• Region C: are exponentially suppressed, so, ∆ i ≫ Y eq i and Equation 7 becomes, After solving the differential equations between T fi to T with ∆ i (T ) ≫ ∆ i (T fi ), we get, Following, ∆ 1 (T f1 ) ≈ c 1 y eq 1 (T f1 ), ∆ 2 (T f2 ) ≈ c 2 Y eq 2 (T f2 ), at T = T CMB = 2.35 × 10 −13 GeV, the DM yields are, Using the analytical approximate solutions, relic density as a function of conversion cross-section for both WIMP In the similar way it is easy to evaluate the present DM relic for the inverse mass hierarchy (m 1 < m 2 ) where one necessary input is T f1 > T f2 and ∆ 2 (T f2 ) = c 2 y eq 2 (T f2 ) with 1 .The transition of FIMP into a pFIMP depends not only on the temperature but also on the interaction rate.This can be understood, just by analysing the relation between interaction rate (Γ) and Hubble expansion rate (H). Figure 6 depicts it for both mass hierarchies: m 1 > m 2 (left), m 2 > m 1 (right).Γ ≃ H has been represented by the green curve.The red star points correspond to Γ(T ∼ T RH ) = H(T ∼ T RH ) in both the figures.In our model-independent analysis, we have neglected the temperature dependence of the conversion cross-section, which may alter the high temperature behaviour of the plot.

Dynamics of the pFIMP in presence of a SIMP
The pFIMP behaviour could also be achieved in presence of SIMP [33] when DM-DM conversion is sufficiently large (∼ weak scale).The cBEQ for pFIMP-SIMP scenario can be written as: In the above subscripts s denote SIMP and w denote pFIMP, x .For simplicity we only have taken 3DM → 2DM depletion process for SIMP, but (nDM → 2DM with n > 3) can also be taken to show the same effect.We solve the cBEQ with increasing order of conversion cross-section to show the change from FIMP to pFIMP case in Figure 7 for both mass hierarchies m w > m s (left panel) and m w < m s (right panel).The features remain very similar to pFIMP-WIMP scenario.In figure 8, we show the relic density of the DM components ϕ 1 , ϕ 2 for a two component scalar DM scenario, as a function of the DM-DM interaction coupling λ 12 .Both the hierarchies m ϕ1 > m ϕ2 (a) and m ϕ1 < m ϕ2 (b) are shown.We see a spectacularly similar feature to that of model independent analysis presented in fig. 1 of the main text, showing that the generic properties assigned to pFIMP are generically valid.Direct search prospect of pFIMP via WIMP loop The Feynman diagrams for pFIMP (ϕ 2 ) to interact with SM via Higgs.Sum of these tree (left), 1-loop (middle) and counter (right) diagram gives the total renormalized λ hϕ2ϕ2 coupling at 1-loop level.
We will now provide the details of direct search prospect of the pFIMP (ϕ 2 ) in the two component scalar singlet model.The renormalized amplitude is given by (see Figure 9), Here, The transfer momentum q h is momentum transfer from initial to final state particles and 1/2 is the symmetry factor due to the scalar loop.There is an ambiguity of choosing the renormalisation scale, its unknown essentially.The direct search cross-section depends on the choice.For, q 2 h → 0, the cross-section is vanishingly small [34], however if we choose q 2 h = 4m 2 ϕ2 , a scale relevant for the pFIMP production, the contribution is substantial [35].So we choose q 2 h = 4m 2 ϕ2 as the renormalization condition setting point and calculate the counter term which removes the divergence from 1-loop amplitude Γ 1−loop hϕ2ϕ2 to yield, Therefore, using Equation 32 the effective spin-independent DM-nucleon direct-detection cross-section in zero transfer momentum (q 2 h = t → 0) limit turns out to be, = 0.037(0.037)[35].We use Equation 33 for pFIMP direct search cross-section while scanning the parameter space as shown in the main text.

FIG. 4 :
FIG. 4: Relic density allowed points of the model described by Eq. 3 in σ SI ϕ1 eff − m ϕ1 plane.The parameters varied are {m ϕ2 , λ 1H , λ 12 }.σ SI ϕ2 eff is shown in colour bar.See heading and figure inset for other parameters and experimental limits.

FIG. 5 :FIG. 6 :
FIG. 5: Analytic solution of the cBEQ (represented by thick line) and numerical solution (represented by points) of Eq. (4) for pFIMP-WIMP case is shown as a function of conversion cross-section.The unknown constants are chosen as c 1 = 2 and c 2 = 7. Vertical dot-dashed line corresponding to σ 12 = σ 1 .

FIG. 8 :
FIG. 8: Variation of relic density with λ 12 for the the two component scalar DM model, where ϕ 1 is WIMP and ϕ 2 is pFIMP for m ϕ1 > m ϕ2 (a) and m ϕ1 < m ϕ2 (b).See figure insets and headings for other parameters kept fixed for the plot.
In the above,µ n = mnm ϕ 2 mn+m ϕ 2where m n is the nucleon mass, f n = 2

TABLE I :
Generic two-component DM scenarios, order of the interaction with SM as well as between DMs, possibilities of satisfying relic abundance and being probed in direct (indirect) search experiments.

TABLE II :
Some benchmarks for the two-component scalar singlet model (Eq.(