Freeze-in bino dark matter in high scale supersymmetry

We explore a scenario of high scale supersymmetry where all supersymmetric particles except gauginos stay at a high energy scale $M_{\rm SUSY}$ which is much larger than the reheating temperature $T_\text{RH}$. The dark matter is dominated by bino component with mass around the electroweak scale and the observed relic abundance is mainly generated by the freeze-in process during the early universe. Considering the various constraints, we identify two available scenarios in which the supersymmetric sector at an energy scale below $T_\text{RH}$ consists of: a) bino; b) bino and wino. Typically, for a bino mass around 0.1-1 TeV and a wino mass around 2 TeV, we find that $M_{\rm SUSY}$ should be around $10^{12-14}$ GeV with $T_\text{RH}$ around $10^{4-6}$ GeV.


INTRODUCTION
Supersymmetry (SUSY) [1][2][3][4][5][6] is a significant theoretical framework aiming at extending the Standard Model (SM), drawing inspiration from the pursuit of a quantum gravity theory, particularly within the context of superstring theory.In the field of phenomenology, SUSY not only provides a viable candidate for dark matter (DM) which plays a crucial role in the formation of large-scale structures in the universe, but also contributes to the renormalization group running of gauge couplings through the inclusion of additional particles near the electroweak scale.This property of SUSY facilitates the potential unification of the three fundamental forces at high energy scales.It has long been postulated that SUSY DM takes the form of Weakly Interacting Massive Particles (WIMPs) that can be probed through diverse experiments [7][8][9][10][11][12][13][14][15][16][17].However, the absence of confirmed DM signals poses significant challenges to the standing of SUSY DM.The current LHC search results indicate that SUSY particles seem to be heavier than the electroweak (EW) scale [18,19], thus challenging the WIMP paradigm of SUSY (for recent reviews on SUSY in light of current experiments, see, e.g., [20][21][22]).
Given the current situation, in this study we consider an alternative scenario of SUSY DM in which gauginos are located at a low energy scale while all other SUSY partners exist at a significantly higher scale M SUSY .This scenario is a special case of the Split SUSY [23][24][25][26] where higgsinos are also taken to be a similar scale as sfermions.One should note that the Higgs sector in this scenario is fine-tuned [27][28][29][30][31][32][33][34] and it might be a consequence of the anthropic principle.However, in this work we will assume that SUSY still provides a candidate of DM and we will specifically consider the Minimal Supersymmetric Standard Model (MSSM).Since the measurement of gamma-ray from the MAGIC [35] has strongly constrained the possibility of wino DM1 , the only viable DM candidate in the MSSM is bino.However, it is widely known that pure bino DM is typically overabundant from the freezeout mechanism [36] due to its weak coupling with the visible sector [37,38].Alternatively, a bino particle with a rather weak coupling may serve as a suitable candidate for Feebly Interacting Massive Particle (FIMP) DM with a correct relic abundance generated via the freeze-in mechanism [39], with assumptions that the reheating process solely occurs in the Standard Model (SM) sector and the reheating temperature T RH is lower than the SUSY scale M SUSY .
In this work we study the possibility that the bino DM in MSSM is generated via the freeze-in process during the early universe.We assume that all MSSM particles except gauginos share similar mass M SUSY which is much higher than the reheating temperature T RH of the universe.To generate enough relic abundance of bino dark matter, we always require the bino mass lower than the reheating temperature.While for the mass of wino or gluino, they could be either higher or lower than the reheating temperate T RH depending on the different scenarios we consider.
The paper is organized as follows.In Section 2 we present the model set up.In Section 3 we first overview the physics related to dark matter and then study the dominate channels for bino freeze-in production.In Section 4 we give the numerical results and discuss the experimental limits on the model parameter space relevant for our scenarios.We draw the conclusions in Section 5 and leave the calculation details in Appendices.
Since we are considering a scenario of high scale supersymmetry in which only gauginos are at low energy scale, the relevant Lagrangian terms are where A = 1, 2 For the Higgs sector, we need a SM-like Higgs boson H SM near the electroweak scale [40,41].This is obtained from the mixing between the two Higgs doublets H u and H d in the MSSM: where σ2 is the second Pauli matrix, and tan β = ⟨H 0 u ⟩/⟨H 0 d ⟩ with ⟨H 0 u ⟩ and ⟨H 0 d ⟩ being the vacuum expectation values (VEVs).Such mixings can be realized by properly choosing Higgs mass parameters µ, M Hu , M H d , and b.The subscription "NP" in H NP denotes the new physics (NP) Higgs doublet in the MSSM accompanying the SM one 2 .Since the mass parameters M Hu , M H d , b, µ are all much larger than the electroweak scale, a tuning of these parameters are needed to get a light Higgs at electroweak scale [27,28,[30][31][32].We need also match the Higgs self-coupling to be the SUSY value at the scale of M SUSY , Note that the Higgs self-coupling λ becomes very small at high energy scale due to the RGE running, and thus the β value should get close to π/4 and tan β ≈ 1.We will fix tan β = 1 as the benchmark parameter throughout this work for simplicity.Generally, when considering physical processes at temperature T ≪ M SUSY , we can integrate out the heavy mediators with mass µ, M f ∼ M SUSY ≫ T RH and get the following of the complex but electrically neutral scalars H 0 SM , H 0 NP into real and imaginary parts.However, one needs to beware that G ± SM , H 0 SM contain the Goldstone boson modes to be absorbed into vector gauge bosons W ± , Z 0 after the electroweak symmetry breaking (EWSB).
effective operators at the level of dimension 5 and 6, respectively, dimension-5: (2.7) Since we assume the mass parameters of higgsinos µ and sfermions M f around M SUSY , the dominant process would be from the dimension-5 (dim-5) operators.Nevertheless we also present the processes related to dim-6 operators for completeness.We acknowledge that a majority of the significant processes are evaluated at energy scales considerably beneath M SUSY .The recommended approach entails initiating the integration procedure for the massive particle to derive the effective operators of dimension 5 and 6, along with their corresponding Wilson coefficients, within the realm of M SUSY .Subsequently, the computation of these Wilson coefficients at the pertinent scale is achieved by employing the Renormalization Group Equations to track the evolution of the operators.
Notably, there exists a potential correction to the primary outcome, potentially on the order of O( 1), yet the fundamental framework remains robust.We leave the investigation of this effect for future study.

Particle spectrum
Despite the existence of new Higgs bosons and many supersymmetric partners of the SM particles, the MSSM particle spectrum we consider in this work consist of two sectors distinguished by their characteristic mass scales.Although not making significant difference for the mass spectrum structure before and after EWSB, we take the pre-EWSB case as an illustration.
• Heavy sector, inactive after cosmological reheating In the above we utilized gauge eigenstates for description, since B, W do not mix with higgsinos Hu , Hd before EWSB when the SM Higgs H SM has not acquired the VEV.

Bino production from freeze-in mechanism
In the early stage of universe before EWSB when the gaugino states B, W do not mix with higgsinos Hu , Hd , pure B acting as DM can only interact with SM via mediators with heavy mass near the scale M SUSY , as shown in Fig. 1.Due to the suppressed interacting strength, the cosmological production of bino DM in our scenario proceed via the freezein mechanism.In the follows we consider the contributions to bino DM production from several typical processes 3 .Hermitian conjugated processes also exist when the amplitudes are complex.See more discussions in the main texts.

Case I: bino freeze-in from HH * → B B
This case corresponds to the left panel of Fig. 1 but without winos W .After integrating out the heavy higgsinos, the relevant dim-5 effective interaction is given by (the details are given in Appendix A) where In the subscription HH * → B B on the left side (and hereafter when not causing any confusion), we denote H SM as H to simplify the notation, and all fields in the initial and final states of the process should be understood in the sense of physical particles 4 .With more details given in Appendix B, Eq.(3.1) would induce the Boltzmann equation of the bino number density: The above equation can be modified to a differential equation about bino yield Y B = n B /S (S is the entropy density) and temperature T : where M Pl ≈ 1.22 × 10 19 GeV is the Planck mass, S = 2π 2 g * T 3 /45, and Hubble expansion rate H ≈ 1.66 √ g * T 2 /M Pl with g * = 106.75before EWSB.Performing a simple integral from reheating temperature, it can be found that the final yield of B depends on the reheating temperature T RH which corresponds to the Ultraviolet (UV) freeze-in scenario [39,43]: and the corresponding current relic abundance is given by (3.4)

Case II: fermion scattering process f f → B B
After integrating out sfermions with heavy mass M q, l ∼ M SUSY in the right panel of Fig. 1, the effective interactions between SM fermion pair and B pair have the following form at dimention 6 (for more details, see Appendix C): where for simplicity, we consider an universal mass for all the fermions, i.e.
Thus the Boltzmann equation is and correspondingly, )

Case III: gluino/wino scattering or decay processes
As indicated by blue colored arrows in Fig. 1, the 2 → 2 scattering processes consist of two ways of generating bino DM when combining U(1) Y with SU(2) L or SU(3) C interactions, related by the cross symmetry.Moreover, we can also have the red colored arrow indicating decay processes generating binos before (after) EWSB when the cosmological temperature drops below the scale of M 2 or M 3 (equivalently, when the age of the universe reach the lifetime of W and G).Similar to the previous two cases, integrating out heavy higgsino and sfermions would generate the following dim-5 and dim-6 effective operators: Note that the index f in the second line includes only SU(2) L doublets, while the index f in the third line includes only quarks.To highlight the difference, we use index a and b to denote generators of SU(3) C and SU(2) L interactions, respectively.Correspondingly, λ a and σ b are Gell-Mann and Pauli matries, respectively.
In the following, we consider the contributions to the bino DM production from 2 → 2 scattering and 1 → 3 decay separately, while leaving the effects of 1 → 2 decay appearing after EWSB in Section 4.1.

Case III A: 2 → 2 scattering involving gluino/wino
With more detailed given in Appendix D, the collision terms in the Boltzmann equation for dim-5 and dim-6 operators are approximated as (ignoring the masses of all external particles) where N conj = 2 denotes the effects of conjugated process.

Case III B: decay of gluino/wino
Following the method in [39] with f G and f W approximated by e −E G/T and e −E W /T , the Boltzmann equation of freeze-in production for the 1 → 3 decay processes is where g G = 16 and g W = 6 are the internal d.o.f. of G and W , respectively.The expressions of decay width involved in the above results are listed in Appendix E. Changing variables to yield Y B and temperature T , we then integrate over temperature evolution to obtain the final yield.If reheating temperature T RH is much larger than M 2 and M 3 , then the final yield from 1 → 3 decay can be approximated by It is worth pointing out that the above result is not sensitive to T RH .Taking a low reheating temperature T RH = 1.1 M 3 as an example, increasing the value of T RH does no modify the result significantly.
In addition to the 1 → 3 decay, we should also note that wino W with mass M 2 < T RH keeps staying in the thermal bath until reaching its freeze-out moment yielding a relic wino number density, which would later convert to the equal amount of bino number density n B via 1 → 2 decay W → B + h after EWSB occurs.Depending on the bino mass M 1 , this freeze-out component would also contribute to the total bino DM abundance in today's epoch.We checked that with wino mass M 2 = 2 TeV, the 1 → 2 decay contribution of Y 1→2 B to final bino yield is around 25% (1%) on the percentage level for M 1 = 1(0.1)TeV [44], thus not affecting the freeze-in domination scenario of this work.We properly include the wino freeze-out contribution in our results.There is also contribution from gluino late time decay.However, to avoid the constraints from BBN, we have to set the gluino mass higher than the T RH , thus we do not include its contribution here.

NUMERICAL RESULTS AND DISCUSSIONS
In Fig. 2 we show the required scales of µ (M f ) for dim-5 (6) operators with various T RH to produce the observed bino DM relic abundance.The upper (lower) two lines correspond to dim-5 (6) operators.We can see that due to the more suppression of dim-6 operators, the needed M f are generally O(10 −4 ) smaller than µ in the dim-5 case.If we assume O(µ) ≈ O(M f ), in order not to overclose the Universe, the dim-6 contributions would be completely negligible.
From Fig. 2, we can see that for the case M B < T RH ≪ M W , the dominant production of bino dark matter is from the process HH * → B B from the dim-5 operator.Generally, M SUSY should be around 10 13−14 GeV for T RH < 10 6 GeV.Since the final relic abundance is proportional to T RH /µ 2 , the M SUSY could continue increasing if the reheating temperature T RH becomes higher.Note that this is similar to the model of Higgs portal to fermion dark matter which are studied in [45], with which we find our result are consistent.We emphasize that our model is motivated by a more complete framework and [45] falls into one of cases we consider.Moreover, For the case M B , M W < T RH , we find the wino-included process can largely enhance the annihilate rate and a higher scale is needed to satisfy the relic abundance.In this case, M SUSY should be around 10 14−15 GeV for T RH < 10 6 GeV.
Notice that if the gluino is in the thermal equilibrium with SM in the early universe and the sfermions mediating the gluino decay are heavier than 10 9 GeV, the lifetime of the gluino could be longer than the age of the Universe when the big bang nucleosynthesis (BBN) happens, leading to energy injection into the cosmic plasma and altering the BBN profile.In all cases considered in this work we find M SUSY is much larger than 10 9 GeV, therefore we always need M G ≫ T RH to avoid the limit from BBN [46].More discussions on BBN limits are given in 4.1.In Fig. 3 we show the comparison of final contributions and intermediate profile of UV and IR freeze-in processes to the bino DM relic abundance.It can be clearly seen that the IR freeze-in final yields from wino 3-body decays are negligible compared to that of UV freeze-in processes generated by 2 → 2 annihilation.Moreover, the critical production moment determining the final yield of UV freeze-in locates in a much smaller x (and thus much higher temperature) than the IR freeze-in case.

Limits from BBN
After EWSB, the SM-like Higgs doublet needs to be replaced by: where v = 246 GeV is the VEV of SM Higgs5 and h is the observed SM-like Higgs scalar.G ± (G − = (G + ) * ) and G 0 are Goldstone bosons that form the longitudinal modes of SM gauge bosons W ± and Z.As mentioned earlier, the SM-like Higgs VEV will generate mixings among the gauge states B, W , Hu , Hd and form mass eigenstates of charge-neutral neutralinos χ1,2,3,4 and charged charginos χ± 1,2 (with ascending mass order inside sectors of neutralinos and charginos, respectively).For the scenario considered in this work, the

component of neutralino χ0
1 ( χ0 2 ) is dominated by bino B (wino W 3 ), and component of chargino χ± 1 is dominated by winos . More details of the approximated masses and couplings can be found in [47][48][49].In the following, we would utilize the language of gauge states (bino B, wino W , higgsinos Hu , Hd ) and mass eigenstates (neutralino χ0 , chargino χ± ) interchangeably before and after EWSB.Now we study the limit of BBN on our scenario from lifetimes of neutralinos, charginos.In our scenario, only neutralino χ0 1 ≈ B, χ0 2 ≈ W 3 and chargino χ± ) existed in the primordial thermal bath.Due to the loop induced mass-splitting between χ± 1 and χ0 2 , chargino χ± 1 can have the 2-body decay χ± 1 → χ0 2 π ± [50][51][52][53].It makes the lifetime of χ± 1 much shorter than 1 sec, and thus not affecting the BBN profile.However, we need to scrutinize the lifetime of χ0 2 more carefully.If χ0 2 decays after the onset of BBN, then the highly energetic decay products will cause the photodissociation or hadrodissociation and thus change the final abundances of light elements.So a bound from BBN can be put on the model parameters, especially on the SUSY scale M SUSY [54,55].
It is easy to see that Fig. 1 implies the 2-body decay mode of χ0 2 → χ0 1 h at the level of dim-5 after EWSB, in which case we will have: where the first term containing Goldstone boson G ∓ can be understood in the context of Goldstone equivalence theorem (GET) for χ± 1 → χ0 1 W ± .It should be noticed that Eq.(4.2) does not contain the three-particle coupling G 0 W B and thus would not provide a way of inferring the 2-body decay mode χ0 2 → χ0 1 Z via the GET.In fact, χ0 2 → χ0 1 Z comes from the gauge covariant kinetic terms of gauginos and higgsinos combined with gaugino mixings after EWSB.However, the decay width of χ0 2 → χ0 1 Z suffers from an extra suppression of 1 µ 2 embedded in the mass mixings compared to χ0 2 → χ0 1 h and thus can be ignored [56] .Therefore, we have the following dominant 2-body decay (see Appendix F for more details): Using the GET we would obtain the same results for Γ χ± 1 → χ0 1 W ± when neglecting the gauge boson masses.
In this work, we apply the limit of BBN to the requirement that lifetime of χ0 2 must be less than 0.3 second [39].In Fig. 4 , we show the interplay between BBN constraints and freeze-in production, where regions below black lines are allowed while region above blue lines are allowed.We can see that for bino mass around 0.1-1 TeV, an upper bound of M SUSY ∼ 10 14 TeV is needed to satisfy both phenomenological requirements.

Limits from direct/indirect detection
Our scenario can easily escape from the current limits from the direct and indirect detection.In the case of direct detection, Eq.(3.1) after EWSB would generate the t-channel scattering of χ0 1 with quarks and gluons in SM neucleons mediated by SM Higgs, of which the event rate is suppressed by 1/µ 2 and thus negligibly small.In the case of indirect detection, which is basically the inversed process of the freeze-in DM production, would generate cosmic rays via DM pair annihilations χ0 1 χ0 1 → h * → SM and χ0 1 χ0 1 → hh → SM, of which the flux is again suppressed by 1/µ 2 and thus not violating the current experimental bounds.

Limits from the LHC
The collider signals of our scenario mainly come from pp → χ± 1 h > O(10 −2 ) s would make χ0 2 traverse through the whole detector before decaying without leaving any energy deposit in the calorimeters, thus can easily evade the current ATLAS [57] and CMS [58] searches for displaced vertex signals at √ s = 13 TeV.As for the disappearing track signature of χ± 1 → χ0 2 π ± , ATLAS [59] and CMS [60] also performed dedicated searches using dataset at √ s = 13 TeV and imply that χ± 1 , χ0 2 should be heavier than 500-600 GeV, therefore our benchmark points with M 2 = 2 TeV are still available.

Discussions
Before ending this section, we discuss some details concerning the SUSY mass spectrum.Firstly, our findings indicate that, to achieve the correct dark matter relic abundance through the UV freeze-in mechanism, the typical mass scales of SUSY particles (excluding gauginos) should be in the range of 10 13−14 GeV.An intriguing question arises concerning whether the SM Higgs with mass around 125 GeV can be accommodated within this framework.In heavy SUSY scenarios, as discussed in [28,31], sparticle masses around 10 13 GeV are still viable, particularly when considering tan β = 1 and allowing for the uncertainty in SM parameters within 1σ range.Expanding the range of uncertainty in SM parameters, particularly the top Yukawa coupling, to 2σ range allows for a significant upward adjustment of the SUSY mass scale.Notably, the work of [32] delves into high-scale SUSY within 3σ uncertainty for Standard Model parameters, with findings indicating that for tan β = 1 a SUSY scale as high as 10 16 GeV remains consistent with the observed SM Higgs mass.This underscores the importance of considering a reasonable range of uncertainty in SM parameters when assessing SUSY scenarios.However, it is crucial to note that future precision measurements of SM parameters hold the potential to rigorously scrutinize SUSY scenarios.Therefore, our model stands poised for being tested against these precise measurements, providing an avenue for further validation and refinement.
Secondly, in our study we adopted the assumption of a gluino mass greater than the reheating temperature T RH to avoid potential conflict with BBN constraints.Concurrently, we presented typical mass ranges for the bino (and wino) falling within the span of 0.1-1 TeV, with T RH estimated at around 10 4 − 10 6 GeV.This naturally entails the requirement for a substantial hierarchy between the gluino mass and the bino (as well as the wino) mass.Achieving such a hierarchy within the domain of supersymmetry calls for a meticulous consideration of the scenarios associated with SUSY breaking and mediation.One plausible avenue involves postulating non-universal gaugino masses.This can be accomplished by ascribing distinct representations to the SUSY breaking superfield Φ with non-vanishing F-terms (see e.g.[61][62][63]).While the above framework provides a well-recognized means of introducing a phenomenologically oriented hierarchy among gaugino masses, attaining the desired mass ratio between the gluino and the bino/wino may necessitate a fine-tuning of the contributions arising from these different representations.
It is crucial to emphasize that our present research predominantly focuses on delving into the phenomenological aspects, especially within the realm of dark matter.Acknowledging that a comprehensive model incorporating precise calculations of the Higgs mass and the requisite mass hierarchy for gauginos is undoubtedly imperative, we intend to actively explore the feasibility of incorporating these elements in our future work.

CONCLUSION
We studied a scenario of dark matter generated from UV freeze-in mechanism, realized in the framework of high scale MSSM.The bino is the dark matter candidate and its relic abundance is generated by the freeze-in processes via the dim-5 or dim-6 operators.We found that the SUSY scale M SUSY should be around 10 13−15 GeV for reheating temperature in the range of 10 4−6 GeV.We also illustrated the interplay between BBN constraints from neutral wino decay and the experimentally observed dark matter relic abundance, implying an upper bound of M SUSY around 10 14 GeV for wino mass around 2 TeV and bino mass of 0.1 ∼ 1 TeV.
After integrating out higgsinos with mass µ, we obtain dim-5 operator between SM Higgs H SM and B DM: where

B BOLTZMANN EQUATION AND CALCULATION DETAILS OF FREEZE-IN DM IN CASE I
In the homogeneous and isotropic universe, the production of bino is described by following Boltzmann equation [36]: with n B denoting the number density of bino particle, and H is the Hubble expansion rate.
Taking HH * → B B ( B means the physical bino particle) in Case I of Section 3.2 as an example, we have [64] where f i,j,k,l are the phase space distribution functions.The number density, taking f i as example, is defined as in which g i is the internal degree of freedom (d.o.f.) of particle i.The factor N denotes the number of particles under consideration produced in the final state and the factor 1/S originates from the phase space suppression due to the identical particles in the initial and final states.For HH * → B B we have N = 2 and 1/S = 1/(N !) = 1/2.After some manipulations and neglecting the negligible backward process, we have [64] where p kl is similar to p ij .After summing over all bino spin states s 1 , s 2 and isospin states of the SM-like Higgs , we have the amplitude square (s is the square of the central energy): We modify the MSSM model file available in FeynRules [65,66] to highlight the gauge state interactions and then export to FeynArts [67] augmented with FeynCalc [68] to perform the calculation.
Since we are considering freeze-in production of B, f 1,2 in Eq. (B.2) can be ignored.We can further approximate f 3,4 by Maxwell-Boltzmann distribution, i.e. f 3,4 ≈ e −E 3,4 /T .Then the collision term can be rewritten as [43,64,69] Here K 1 is the Bessel function of the second kind, and we treat the SM-like Higgs in the initial state as being massless.In the case where M 1 ≪ T , the collision term can be approximated as (using After integrating out sfermions with mass M f in the right panel of Fig. 1, we obtain dim-6 operators between SM fermion pair and B pair: where for simplicity we consider an universal mass for all the fermions, i.e.M f ≡ M q = M l.
The amplitude squared terms in the collision term for f f → B B scattering process is given by 6 where N flavor = N color = 3.As in Eq. (B.2), if we neglect bino mass, then the collision term can be approximately given by (using 6 Again, fields in the initial and final states in the process should be understood in the sense of physical particles, where f denotes the physical anti-particle.Discussion on the naming convention of particles, states and filed can be found in, e.g.[42].

D THE CALCULATION DETAILS IN CASE III A
When neglecting all particle masses in the final state, we have where dm 2 12 , dm 2 23 are defined in [70].where Finally, we have [47] Γ χ0 2 → χ0 (F.4) Using the GET we would obtain the same results in the high energy limit for Γ χ±

Figure 2 .
Figure 2. Values of µ and M f to produce the observed DM abundance via the UV freeze-in processes.See more discussions in the main texts.

Figure 3 .
Figure 3.Comparison between UV freeze-in and IR freeze-in.Note the difference between temperatures indicated by x = M 2 /T producing the correct relic density of bino DM.

Figure 4 .
Figure 4. Interplay between BBN constraints and freeze-in production, where regions below black lines are allowed while region above blue lines are allowed.
denote the left-handed twocomponent Weyl spinor of SM quarks and leptons, where the bars are simply notations and do not mean the Dirac conjugation.Hypercharges are given by {Y

3
(3) correspond to the SM gauge group U(1) Y , SU(2) L , SU(3) C , respectively, and a denotes the corresponding indices in adjoint representation of group A. Fields Ṽ A,a , Hu , Hd , f are the superpartners of the SM vector gauge bosons V A,a = B, W 1∼3 , G 1∼8 , scalar doublets H u , H d and fermions f .The fields H u , H d , Hu , Hd are defined as Schematic plots for interactions of DM composed of pure B with SM after cosmological reheating considered in this work, which would induce dimension-5 (left) and dimension-6 (right) effective operators.The SM Higgs H SM originates from the mixing between MSSM Higgs doublets H u , H d .Colored lines indicate the direction of freeze-in production when applicable.Additional χ± THE CALCULATION DETAILS OF 2-BODY DECAY AFTER EWSBAs discussed in Section 4.1, we have the following 1 → 2 decay possibly affecting the cosmological BBN: