Self-interacting dark matter from late decays and the $H_0$ tension

We study a dark matter production mechanism based on decays of a messenger WIMP-like state into a pair of dark matter particles that are self-interacting via exchange of a light mediator. Its distinctive thermal history allows the mediator to be stable and therefore avoid strong limits from the cosmic microwave background and indirect detection. A natural by-product of this mechanism is a possibility of a late time, i.e., after recombination, transition to subdominant dark radiation component through three-body and one-loop decays to states containing the light mediator. We examine to what extent such a process can help to alleviate the $H_0$ tension. Additionally, the mechanism can provide a natural way of constructing dark matter models with ultra-strong self-interactions that may positively affect the supermassive black hole formation rate. We provide a simple realization of the mechanism in a Higgs portal dark matter model and find a significant region of the parameter space that leads to a mild relaxation of the Hubble tension while simultaneously having the potential of addressing small-scale structure problems of $\Lambda$CDM.


I. INTRODUCTION
The standard ΛCDM cosmological model incorporates dark matter (DM) in the simplest way possible, i.e., a non-interacting cold matter component with constant equation of state throughout its cosmological evolution. A scenario of this type is not only simple and remarkably successful in explaining the Universe at large scales but also well motivated in many theories beyond the Standard Model (SM) of particle physics. However, the shortcomings of ΛCDM at small scales, e.g., the diversity [1,2], too big to fail [3], missing satellites [4][5][6] and core-cusp [7][8][9] problems, as well as tensions between parameters inferred from local and global cosmological measurements, most notably the Hubble parameter H 0 [10][11][12] (see, e.g., [13] for a review), may be viewed as a hint that the CDM paradigm is in fact too simple. Indeed, it is well known that at least some of the small scale problems can be simultaneously addressed if DM possesses significant self-interactions preferably with velocitydependent cross section (see, e.g., [14] for a review). Additionally, varying equation of state, e.g., due to late time conversion of a small fraction of DM into radiation, has been shown to have the potential of reducing the H 0 tension [15,16] (but see also [17,18]).
It is an intriguing question if both small scale problems and ΛCDM tensions can be simultaneously resolved through a modification of only the DM component. This point has been addressed in thermally produced self-interacting dark matter models featuring strong Sommerfeld enhancement in [19,20], where it has been demonstrated that late time annihilations can indeed be efficient enough to sufficiently modify the cosmological evolution.
However, models predicting thermally produced DM self-interacting via light mediator often run into problems with observations (see e.g., [21]). The DM annihilation to the mediator pair is greatly enhanced during the recombination epoch by the Sommerfeld effect [22,23] leading to too large energy injection into the plasma, if the mediator decays to visible states. On the other hand, for stable light mediators the overclosure bound is greatly constraining due to their large thermal population.
Several possibilities of how to avoid such limits are known, e.g., by having the mediator decay only to neutrinos or dark radiation (subject to much weaker bounds) or by assuming that the dark sector (DS) is effectively secluded and has much lower temperature than the one of the photon bath. In this paper we propose to utilize a mechanism for DM production akin to the one used in the superWIMP scenario [24] and show that it introduces alternative way of constructing models with velocity dependent self-interactions. In such a setting the DM component arises from decays of an intermediate weakly interacting massive state, which in turn is thermally produced via the usual freeze-out process. This production mode allows the mediator of the interactions in the DS to be absolutely stable while at the same time not overclosing the Universe. Additionally, if the decays of the WIMP-like particle happen at very late times, it is exactly the framework needed for the conversion of dark matter to radiation that might help alleviate the H 0 tension. What is more, only small fraction of the WIMPs decay to light mediators quite naturally, as it is a higher order process compared to the tree-level decay to the DM particles. This is a feature that is difficult to obtain without resorting to multi-component DM models.
This paper is organized as follows. In Sec. II we introduce the mechanism and the example from a generic class of Higgs portal models. Sec. III describes the thermal history, lays out calculations of DM self interactions and late time decay impact on cosmology. In Sec. IV we show and discuss the results of a numerical analysis. Finally, we conclude in Sec. V.

II. THE MECHANISM
The main idea behind the production mechanism studied here is that, if the dark sector is populated by decays taking place late enough that it never reaches chemical equilibrium with the visible sector, then the light mediator is effectively absent from the plasma while still carrying a long range force between DM particles. Therefore, it can be absolutely stable and completely naturally evade all the limits from CMB observations and indirect detection. 1

A. The SM-DS coupling through a portal
A very generic framework naturally encompassing the above mechanism is the scenario when the dark sector is connected to the visible sector only through a weak portal. Here, for concreteness, let us concentrate on a so-called Higgs portal connection between the SM and the DS, which is one of the most simple and natural choices. Multitude of examples of such DM models can be found in the literature (see, e.g., [25] for a review). The scenario is illustrated in Fig. 1 where the connecting SM singlet scalar S is assumed to mix weakly with the Higgs and also have a weak or very weak coupling to the states in the dark sector. 2 A natural choice for S is to be a pseudo-WIMP, i.e., particle undergoing thermal freeze-out with nearstability guaranteed by imposed spontaneously or explicitly broken Z 2 symmetry S ↔ −S.
The perspective of DM portal framework highlights an alternative angle on the studied mechanism: it can be viewed as an extension of the usual Higgs portal freeze-out or freeze-in models to even weaker couplings to the dark sector. Indeed, parameterizing the breaking by a small parameter , one can quite generally distinguish four regimes: 0) weak : the DS reaches chemical equilibrium with the SM independently of the reheating details leading to a thermal population of the dark matter and light mediator -one recovers usual thermal self-interacting model subject to strong limits A) very weak weak: the DS is produced through decay of S and never reaches chemical equilibrium with the SM; the light interaction mediator can be stable and avoid overclosure and CMB limits; viable regime for selfinteracting DM B) ultra weak very weak: the same as A, but leading to S having lifetime on cosmological scales; regime for self-interacting DM with an impact on the H 0 tension C) ultra weak: S is quasi-stable with onset of its decays reaching times of order of Gyr; one 1 Although, for simplicity, we will limit ourselves to stable mediators, we remark that introducing a small decay width provides still a viable model with additional potential phenomenology and detection possibilities. 2 The simplest realization of such setup would assume only one state in the dark sector, which would be stable and provide the DM candidate. Such case, however, cannot accommodate any significant velocity dependent self-interactions between DM particles. ends up with two component DM with only a fraction being self-interacting which can play a role of ultra-strong self-interacting dark matter (uSIDM) [26]; regime potentially addressing the H 0 tension and providing an uSIDM candidate.

SM DS
Regimes 0 and A point to the Z 2 breaking at relatively low energy scales, not much larger than the DM particle mass. Smaller values of leading to scenarios B and C naturally emerge when the breaking comes from some new physics at very high scale, e.g., GUT or even Planck scale.

B. Toy model example
For concreteness, let us consider a dark sector comprised of a Dirac fermion χ charged under new gauged U (1) X broken spontaneously at some higher scale resulting in massive vector A µ . 3 This choice is not crucial in what follows, but exemplifies a very simple and natural realization within a renormalizable model.
The dark sector part of the Lagrangian after the U (1) X breaking reads while the connection with the visible sector is given by the portal where H denotes the SM Higgs boson doublet and in the trilinear term we explicitly pulled out the factor to emphasize that this term is allowed only due to Z 2 breaking. This is a crucial observation because it ensures that S decays predominantly to DS states, if only µ HS is small enough or S light enough that the resulting branching ratio to the SM particles is strongly suppressed compared to BR(S →χχ). Phenomenologically interesting interactions are given in the second lines of both Eq. (1) and (2).

III. PHENOMENOLOGY
Having introduced the framework and defined concrete realization we describe in this section the main properties of such a scenario.

A. Thermal history
The underlying assumption in the discussion of the thermal history of χ is the one of the freeze-in models, i.e., that only SM sector is populated during reheating, while the dark sector has negligible initial number and energy densities.
The connector S undergoes usual WIMP-like evolution where it thermalizes with the SM plasma at early times due to mixing and, typically more importantly, the quartic λ HS coupling. When its annihilation rate drops below the Hubble rate it goes through the freeze-out process. At later times, possibly even after recombination, it decays via S →χχ and also sub-dominantly by construction to SM through the Higgs mixing. In Fig. 2 an illustration of example evolution of mass densities of S, χ and A µ is shown for decay regimes A (solid lines), B (dashed) and C (dotted). In all the cases the χ and A µ undergo a freeze-in type production, which is very inefficient due to smallness of the coupling to S. It follows that their number densities are extremely small until the onset of S decay.
The transitions between the regimes are only indicative and not sharply defined. In particular, the chosen redshift z ∼ 7 line separating cases B and C corresponds to times of oldest observed quasars with supermassive black holes (SMBHs) [27][28][29]. Decays of S around that time can impact the formation rate of the SMBHs, see Sec. IV C. The onset of S decays can also happen later until and beyond the present day, meaning that regime C extends to cover all the possible lifetimes of S.
As can be seen in Fig. 2 in case A the connector S typically needs to chemically decouple with larger number density than would give the correct thermal abundance, since during the decay some of its energy is transferred to the kinetic energy of the χ, which gets redshifted. Note also that annihilation of χχ → AA can have some effect, even if the number densities do not reach equilibrium values, as seen in the small drop of χ density at early times. For later decays in case B and C the χ particles need to be produced with very small kinetic energy, as discussed in Sec. III D 2 below, otherwise will negatively affect the structure formation. It follows that S needs to have the number density just a bit over the observed one which is then nearly completely transferred to the DM.

B. Dark matter self-interactions
In calculating the strength of the elastic scattering between two DM particles at present day velocities v ∼ 10 −3 we follow standard numerical procedure of solving Schrödinger equation described in [14,30]. We use natural units c = = 1.
The relevant quantity with respect to selfinteractions is transfer cross section which is defined as a weighted average of the differential cross section with respect to the fractional longitudinal momentum transfer (1 − cos θ): The differential cross section is given by series expansion into Legendre polynomials corresponding to orthogonal partial waves: The phase shift δ for a partial wave is obtained by solving Schrödinger equation for the radial wavefunction R (r), which describes reduced χ-χ system, given by where v is relative velocity of χ's, µ = m χ /2 is reduced mass of the system and k = µv. Potential term comes from the gauge interactions in Eq. (1). Multiple exchanges of A µ coupled to χ with coupling strength α = g 2 /(4π), result in a Yukawa-type potential: Since we took A µ to be a vector, the interactions are attractive (−) for χχ scattering and repulsive (+) for χχ orχχ scattering. The interaction cross section is then taken as the average of repulsive and attractive interactions.
Far away from Yukawa potential range Eq. (5) has well known solution in terms of spherical Bessel functions j (r) and n (r) (for definitions and properties of spherical Bessel functions see, e.g., Sec. 10.47 in [31]): Therefore, one needs to numerically solve Eq. (5) for a ≤ r ≤ b and match numerical solution at b to the analytical one. We use Numerov method [32,33] which is fourth-order linear method in step size h = (b − a)/n, where n is number of points in the grid. The limit points a and b are determined by demanding that at a Eq. (5) is dominated by the centrifugal term, which means a µv . The upper bound, b, is determined by demanding that the potential term is much smaller than the kinetic term: The illustration of the thermal history of S (blue), χ (black) and A µ (orange) with example parameter choices leading to early (regime A, solid lines), late (regime B, dashed) and very late (regime C, dotted) decays of S. The borders of the regimes are indicative and not sharply defined. In particular, the redshift z ∼ 7 line corresponds to times of oldest observed quasars with SMBHs -see text and Sec. IV C for details.
The resulting phase shift is determined by gluing the numerical solution with asymptotic one at the endpoint of the grid [34]: where R is a wave-function obtained numerically and j , n are spherical Bessel functions. We calculate phase shifts until convergence of Eq. (3) where we consider σ to be converged if successive values obtained for max and max → max + 1 differ by less than 0.1%.
The numerical solution is strictly needed only in the resonant regime, which occurs when αmχ m A 1. In other regions of parameter space one can use analytical formulae to speed up the numerical scan. These can be obtained either from perturbative expansion in α (Born regime [35]; applicable when αmχ m A 1) or from classical calculations of charged particles moving in plasma (classical regime [35][36][37][38]; applicable when vmχ m A 1). We find agreement between numerical results and analytical formulae whenever they are applicable.
Both the coupling g and light mediator mass m A governing the scattering cross section are free, essentially unconstrained parameters of the model. It follows that a very wide range of possible self-interaction strengths can be obtained. Two regions are of particular phenomenological interest, on which we will focus: • The first is when σ/m χ ∈ (10 −1 , 10 1 ) g/cm 2 leading to momentum transfer rates in the correct ballpark to address the small-scale structure problems of ΛCDM. Theories giving rise to cross sections in this range are often referred to as the strongly interacting dark matter (SIDM) models.
• The second is the so-called ultra SIDM (or uSIDM) regime with σ/m χ 10 3 g/cm 2 which could resolve the puzzle of supermassive black holes formation.
One possible solution is that a small uSIDM component can, through a gravothermal collapse, form an initial seed which is what is needed for accelerating growth rate of SMBHs at their early stages of evolution [26,39].

C. Late time S decay
Due to the breaking of the stabilizing Z 2 symmetry, the S decays both to DS and SM states. We will assume that the latter are negligible compared to the former, which is the case if only the trilinear coupling µ HS is small enough leading to small mixing with the Higgs. At tree-level the only decay is then S →χχ with width taking the form: is the parameter governing the mass splitting and we introduced exemplary parameter values that lead to late decays. However, at higher order the three-body S →χχA and loop decay S → AA are present and parametrically Γ S→χχA /Γ S→χχ ∼ g 2 and Γ S→AA /Γ S→χχ ∼ g 4 where the former is also potentially significantly affected by the available phase space, especially if δ 1. One can see that S decay naturally results in few % of energy being transferred to radiation and therefore one obtains a complete one-component DM model with the property desired for alleviating the H 0 tension.
In more detail, final decay products will be either non-relativistic (in tree decay S →χχ, acts as dark matter), relativistic (in loop S → AA, acts as dark radiation) or mixed (in three body S →χχA). In the latter case we adopt a prescription that χ will always act as matter (very good approximation as long as δ is small, as assumed), while A µ will be counted as matter if its kinetic E A < m A , otherwise as radiation. 4 The differential three-body decay rate reads where the amplitude M S→χχA is given by: where p 1 is momentum of A and p 2 , p 3 are momenta ofχ and χ, respectively, in the rest frame of S. * r (p 1 ) is polarization vector coming from external A µ .
Integrating over the whole kinematically allowed region, we get total Γ S→χχA . However, to calculate the fraction of energy transferred to radiation we need to separate the region where A µ is relativistic at decay. We will approximate this fraction by the quantity: where is the fraction of the decay width resulting in A µ having kinetic energy equal or larger to its mass.
The one loop decay S → AA is of higher order in perturbation theory, but does not suffer from phase space suppression and transfers all the energy of S to radiation. For calculations we used Mathematica packages FeynCalc [40][41][42] and Package-X [43] to symbolically calculate the amplitude and evaluate the numerical expressions: The amplitude M S→AA is given by: where B 0 and C 0 are two and three-points Passarino-Veltman [44] scalar functions, respectively, while C 00 is coefficient of three-point tensor function proportional to the metric. We follow conventions of [43], where, in particular, the 1/(16π 2 ) is factored-out in their expressions, hence it reappears in Eq. (14). Note that B 0 and C 00 are UV divergent, however their divergent parts actually cancel out in Eq. (16), which renders the whole expression finite.
Before concluding this subsection a comment is in order. If m S ≈ 2m χ , which as we discuss later is expected to be necessary not to spoil large structure formation, then the χs produced in S decay will have small velocities. Since they interact via light mediator creating long range force, there can be a substantial threshold correction. If present, it would mainly result in a shift of the coupling which is not consequential for what follows. The reason is that such a threshold effect would appear in all three decay processes and while one would expect some change in their relative size the inclusion of this effect would be necessary only when high precision is called for and goes beyond the scope of our work.

D. H0 tension and structure formation
In recent years, cosmological probes become increasingly more precise which further constraints alternatives to the standard ΛCDM model. One of the persistent tensions, which actually became more severe with more data, is determination of Hubble parameter. Early Universe observations such as CMB or baryonic acoustic oscillations (BAO) prefer significantly lower value H 0 ∼ 67 km/s/Mpc in comparison to the local Universe observations which determine H 0 ∼ 74 km/s/Mpc. The uncertainties of the measurements are ∼ 1-2% and the resulting discrepancy reaches ∼ 4σ. Currently, no universally accepted solution is known [13], however it is believed that systematic errors in both measurements are unlikely to completely relieve the difference, as they probe the history of the Universe billions of years apart from each other and they would have to skew the results in the opposite directions. One of the possibilities is decaying dark matter (DCDM) where dark matter particle decays partly into dark radiation. As radiation redshifts faster than dark matter, it results in reduced expansion rate at late times as compared to the early times. Therefore, in DCDM model, the Hubble parameter at z = 0, H 0 , can be put in agreement with the evolution of H(z) at higher redshifts, as measured in, e.g., the CMB.
From the point of view of the impact on cosmology, the scenario under consideration has significant similarities with the DCDM model. Therefore, in this section we describe the details of the analysis for the latter and later we use the obtained results to constrain our model.
We used publicly available Boltzmann solver code CLASS [45] in combination with MCMC code MontePython [46,47] to constrain DCDM model and compare with standard ΛCDM cosmology.
In addition to 6 standard cosmological parameters {ω b , ω cdm , ln 10 10 A s , n s , 100θ s , τ reio } [10], we scan over two additional ones: Γ and F. They denote decay width and fraction of DCDM that decays into dark radiation, respectively. Note that in the context of our model, the latter parameter was already introduced in Eq. (12), while Γ is the total decay width of S. We use thus obtained cosmological limits on Γ and F to find the parameter space regions of our model that is preferred from the perspective of cosmological data.

Cosmological scan
We performed three separate scans using in each case the same likelihoods. They correspond to ΛCDM, DCDM with broad prior on decay lifetime (later called short) and DCDM with prior on decay lifetime constrained to be comparable to current age of the Universe (later called long). The last scan is motivated by [16] which found late DCDM model with lifetime ∼ 20 Gyr to relieve the Hubble tension. In this last case we fixed the reionization time, initial perturbation amplitude A s and its spectral index n s to ΛCDM best fit value, similar to what was done in [16].
We generated chains until the Gelman-Rubin criterion R − 1 < 0.2 is satisfied. The results of the scans are presented on Fig. 3 and 4.
We find two disconnected regions that improve the fit by mildly increasing H 0 , relatively to ΛCDM. They correspond to early decay lifetime (∼ 4 Myr) with small (∼ 1%) fraction going into dark radiation and to late decay lifetime (∼ 5 Gyr) with significant fraction (∼ 10%) going into dark radiation. Such anticorrelation between F and Γ is expected and was previously noted in e.g. [53,54]. In the first case, all of S decayed into χs by the onset of structure formation, therefore χ's self-interactions can improve the structure formation at small scales relative to the ΛCDM. In the second case, a potentially large fraction of final DM component is still in the non-interacting form of S particles that did not yet managed to decay until the present day. In this case the scattering cross section can be even larger and therefore even tiny fraction of ultra-SIDM can serve as seeds of SMBHs.   Comparison with ΛCDM in H 0 -σ 8 plane is shown in Fig. 3. Mean values of the parameters are presented in Table I. We see mild reduction in tension between CMB and low-redshift observations of H 0 and σ 8 in DCDM model.

Structure formation
Late time decays can affect not only the Hubble parameter, but also structure formation as the product of the decay can obtain sufficient energy to freestream.
We impose the bounds coming from halo massconcentration, galaxy-cluster mass function and Lyman-α power spectrum [55][56][57][58][59] as an upper bound on mass splitting between decaying (mother) particle and the resulting massive (daughter) particle. It provides the so-called kick velocity to the daughter particle, which at time of decay is v kick ∼ δ.
We can estimate the free-streaming length of daugh- ter particle using formula from [59]: where, here, τ is the conformal time, integration limits are conformal times corresponding to the time of decay and to the present, Γ is the decay width and a d is the scale factor at the time of decay. Lifetimes considered herein, correspond to Γ −1 10 Gyr for which mass splitting is constrained [57][58][59] to be: δ 10 −2 for short lifetime regime and δ 10 −3.5 for long lifetime regime (note Fig. 11 of [59]). It is worth noting that in short lifetime regime, virtually all of S will decay into self-interacting DM, and the elastic scatterings between the DM particles additionally should suppress free-streaming and somewhat relax the bound on mass splitting. For longer lifetime regime, the limits are stronger because the daughter particle had less time to redshift.

IV. RESULTS
In this section we present and discuss the results of the numerical scans for the three regimes A, B and C. In all the cases we implicitly assume that the correct observed relic abundance of DM is set by adjusting the details of the freeze-out and decay process of S.

A. The SIDM regime
For the values small enough that the dark sector does not thermalize with the SM, but at the same time large enough that S decays happen before recombination the scenario effectively boils down to a self-interacting ΛCDM model. Phenomenologically it has the same properties as many well studied SIDM models (again we refer to, e.g., [14] for a review), with two important distinctions. First, the self-interaction strength is governed by a different coupling that the one giving rise to the relic abundance, opening much wider parameter space. And second, the light mediator can be completely stable rendering the most constraining limits ineffective. In this regime the whole phenomenology is governed by m χ , m A and α.
In Fig. 5 we present the cross sections of the strength needed for solving small-scale structure problems of ΛCDM with rainbow-like palette. The left panel shows the case of small (α = 10 −4 ) while the right panel large (α = 10 −1 ) values of the coupling. One can notice well-known resonant behaviour in lower right part of the plot, which gets more pronounced as α increases. For fixed m A , correct σ/m χ is inversely proportional to m χ and directly proportional to α, as expected. In gray region, parameter space is excluded due to too strong self-interactions [60][61][62]. The light-green region predicts too week selfinteractions to affect cosmology at the small scales in any visible way. The existence of color bands in between, spanning more than an order of magnitude in both masses when taking into account varying α is a demonstration that the proposed mechanism can successfully give rise to the viable SIDM candidate.
Before ending this section let us mention that even in the regime where the S decays happen well before recombination the resulting DM component can help alleviate the cosmological tensions. This was observed and studied in detail in [20] where it was found that if annihilation happens very close to the peak of one of the Sommerfeld effect resonances, the DM can undergo a second period of annihilations at late times [63], leading to conversion of some fraction of matter to radiation. The same effect can appear in our setup, with the modification due to different thermal histories of the DM component. In particular, in [20] the time of kinetic decoupling from the SM thermal bath plays a significant role. However, if χs came from decays of S and were never in equilibrium, then the evolution of their velocity distribution, and consequently the impact of possible late time Sommerfeld enhanced annihilations, would require a separate study.

B. The SIDM from late decays regime
Lowering the values, the lifetime of S extends beyond the recombination and the following decays modify the cosmological model. In this regime the resulting dark matter phenomenology is still governed by m χ , m A and α, but m S (or equivalently δ) and start to have important consequences as well by affecting the kinematics of the decay and the lifetime, respectively.
The main results for this regime are given in Fig. 6. It shows the results of the cosmological scan with priors set to short S lifetime projected onto fixed 5 The gray area on the bottom left is excluded as it leads to too strong DM self-interactions, while the pale green region above is allowed, but does not affect structures at small scales.

FIG. 6:
The results for the SIDM regime B originating from late S decays. The points shown in color coding the value of coupling g predict σ/mχ in the range σ/mχ ∼ (1 ± 10%)cm 2 /g. On top of that the dark green shade denotes the region at the 1σ (68%) level around the mean values of DCDM parameters, which relax Hubble tension in the short lifetime scenario. Gray pluses overlay points that have δ > 0.01 which are in this model in tension with structure formation.
σ/m χ ∼ (1 ± 10%) cm 2 /g in the m A − m χ plane, with colour bar indicating coupling strength g. The dark semi-transparent green region shows 1σ range around the best fit values relaxing the Hubble tension, i.e. the DR fraction of F = 10 −2.41 ≈ 0.004. The light green line denotes the best fit parameters, which depend on m S , hence it is a continuum and not a point, with small width due to numerical resolution. The numerical scan was performed in a grid over four parameters uniquely specifying this fraction: m S , m A , m χ and g, with the condition that the mass splitting, Eq. (9), is small, δ ∈ [10 −6 , 10 −1 ]. The only remaining relevant cosmological parameter, the decay width Γ, can always be brought to correct value by re-scaling the coupling constant.
hance readability of this particular figure, while we emphasize that allowing larger range for the cross section enlarges the allowed parameter space.
The lower right region starting roughly at the right tip of best fit and going along right diagonal, represents the resonant regime. One sees smaller density of points here, compared to Born and classical regimes, and higher values of g are allowed. For largest m A , points are very sparse which comes from the irregular pattern of consecutive resonances which have very small width for large value of α. Roughly half of resonant parameter space is also marked by gray pluses, which denote that those points require large δ, which is in tension with structure formation limits.
The 1σ region is bounded from above by the condition on F . The points above this bound are giving too efficient conversion to DR and manifest in two regimes. The resonant and α ∼ 1 regimes are dominated by loop decay into two As. This region is, partially, also constrained by the limit on δ. For the rest of the parameter space, three body decay of S is dominant.
It is worth stressing that a large parameter space of the model allows for both the self-interactions to be at the right range to potentially solve small scale problems and to decay to correct amount of radiation to help relieving the H 0 tension.

C. The uSIDM regime
Finally, for even longer S lifetimes we enter the two-component DM regime where the χ can be much more strongly interacting. As was noticed in [26] and followed by, e.g., [39], such uSIDM could provide a mechanism of formation of supermassive black holes with masses of order 10 9 M Sun which formed by z ∼ 7. Such SMBHs were observed recently [27][28][29] and provide a challenge for standard formation mechanisms because of their large masses forming at such an early time. The proposed mechanism of [26] is similar to ordinary gravothermal collapse which is believed to be responsible for formation of globular clusters [64] and takes place by ejection of most energetic stars, allowing the rest of system to contract. Concerning black holes formation, uSIDM causes similar process in DM halo and as there is no inhibitor to the pro-cess, SMBH forms. Unfortunately, if uSIDM constitutes the whole of DM, self-interaction rate necessary for gravothermal collapse exceeds the bound set by, e.g., the Bullet cluster. However, a small fraction of even ultra strongly interacting DM is allowed by observations and as showed by detailed simulations in a framework of multi-component DM models in [26,39], can be responsible for boosting the formation rate of SMBHs.
In Fig. 7 we show the results of the long lifetime scan in the F -σ/m χ plane, where F denotes the fraction of uSIDM that existed by z = 7. In light blue we show a region at the 2σ (95%) level around the mean values of DCDM parameters which relax Hubble tension in the long lifetime scenario. Vertical dashed lines denote resulting fractions of uSIDM component at present day. These are significantly larger than the values of F on the x-axis, because of the decays that take place between z = 7 and z = 0. It follows that the whole light blue region leads to a scenario where ultra strongly interacting component constitutes unacceptably large fraction ( 0.4) of DM at late times. Therefore, we find that if the uSIDM arises from decays of an intermediate unstable state the requirement of significant fraction of uSIDM to be already present at z ∼ 7 implies very long lifetimes 40 Gyr giving small fraction F and large scattering cross sections. This is not the parameter region that is preferred for the requirement of relaxing the Hubble tension.
The light green region denotes the parameter space where decay of S happens too late to significantly influence the H 0 tension, but with large enough σ/m χ and F to be relevant for accelerating SMBHs formation. Therefore, significant part of the parameter space corresponds to a scenario of two-component DM which provides a viable mechanism of production of subdominant uSIDM. Note that although some parts of this region lead to a substantial present day uSIDM component as well, the exact limits on F (t 0 ) are rather uncertain and do not exclude the whole parameter space of the model.
The red lines are the results of numerical simulations performed in [39] (Fig. 5, Model A for elastic scatterings) and denote redshifts z = 7 (solid) and z = 15 (dashed). In that work, two component DM scenario was assumed, with constant fraction of uSIDM, F . In our case, F depends both on time of the decay 1/Γ and the fraction 1 − F going into DM component. Hence, the limits presented here should be taken as exemplary and further study conducting numerical simulation would be needed.
To summarize, we find that production of uSIDM via late decay is strongly constrained if one restricts the decay lifetime to be 40 Gyr, which would at the same time relax the Hubble tension. The difficulty lies at the very early time of SMBH formation as z ∼ 7 corresponds to ∼ 0.77 Gyr, while we find that the decay times relevant to Hubble tension correspond to either earlier (∼ 4 Myr) or longer (∼ 5 Gyr) times. However, if the decay is assumed to happen even later than 40 Gyr, it could be a viable mechanism of accelerating the formation of early SMBHs.

FIG. 7:
The results for the regime C. The blue and green regions feature self-interactions strong enough to accelerate SMBHs formation rates, while on top of that the blue region is in the 2σ region around the best fit for the H0 parameter. The dotted vertical lines show contours of uSIDM fraction at present day. See text for more details.

V. DISCUSSION AND CONCLUSIONS
Motivated by the question of how far in solving or alleviating the tensions of ΛCDM one can go by modification of only the dark matter component in a complete particle physics model, we study in this paper the implications of the self-interacting dark matter production mechanism based on (late time) decays of an intermediate thermally produced WIMPlike state. The decay is at the tree-level only into pair of DM particles, while at higher order three-body and loop processes introduce small branching ratio to final states containing the light mediator. This leads to a very natural explanation of why only several percent of the dark matter energy was transferred into radiation, which is a necessary condition for improving the fit to the H 0 parameter. At the same time, if only the lifetime of the intermediate state is smaller than the age of the Universe, the whole non-interacting dark matter is converted into strongly interacting component capable of addressing as well the ΛCDM tensions at small scales. Moreover, this mechanism allows the mediator to be stable and therefore avoid strong limits from the observations of cosmic microwave background and indirect detection.
From a particle physics perspective such scenario is a natural extension of the very well studied models connecting the dark sector with the visible sector by a weak portal. We provide and study a simple example model of this kind, where for concreteness we focus on the Higgs portal. Within this model we perform numerical analysis with the emphasis on the dark matter self-interaction properties and fits to local and global cosmological measurements. We find that the proposed mechanism allows for a perfectly viable selfinteracting dark matter with large parameter space resulting in the elastic cross section of the correct range to address the small scale cosmological problems. Additionally, there exists a significant overlap with the parameter regions required to impact the Hubble tension. However, the resulting improvement of the fit is relatively mild, not offering any improvement over alternative methods to reduce the H 0 tension.
We also consider a scenario when the decays happen much later, with lifetimes of order O(1 Gyr) or larger. We find that in that regime the resulting dark matter consists of two components: dominant non-interacting one and a subdominant component of SIDM or uSIDM type. Although the former has typically too small number density to address the small scale problems of ΛCDM in that scenario, the latter provides a viable model of ultra strongly interacting DM that can help accelerate the formation of the SMBHs and by doing that explain how could they have been formed at times as early as z ∼ 7. Unfortunately, within the studied model we find that the region that could simultaneously alleviate the Hubble tension and provide a mechanism for speeding up the SMBHs formation is not allowed by the observations due to unacceptably large uSIDM component at the times of the Bullet cluster at z ∼ 0.5.
It is worth adding that although all the explicit results given in this work are for a DM being a Dirac fermion interacting via a vector mediator, we have also analyzed a scenario in which the mediator is a scalar leading to purely attractive interactions. This does not introduce any qualitative change and also quantitatively the results are similar to those presented on Fig. 6.
Last but not least, let us comment on the recently reported unaccounted excess of events over the background in electronic recoils around 1-7 keV in the XENON1T experiment [65]. One of the potentially most promising explanations of this excess in terms of new physics involves the existence of a dark photon coupled to SM via kinetic mixing term − κ 2 F µν F µν [66]. What was noticed in that paper, is that both XENON1T excess and observations of cooling anomalies in horizontal branch stars [67][68][69] could be explained by light ∼ keV dark photon with kinetic mixing parameter κ ∼ 10 −15 .
It is interesting to note that the light mediator A µ studied in this work for completely independent reasons, is in fact also the same as the aforementioned dark photon. Although, in our work for simplicity we considered no kinetic mixing in the interaction Lagrangian, one can naturally incorporate it as long as κ 10 −12 , i.e., when the resulting interactions will not significantly affect the thermal history of neither the DM nor the light mediator. 6 Therefore, it is intriguing to note that allowing κ ∼ 10 −15 could be relevant to XENON1T excess, in addition to production of SIDM while simultaneously mildly relaxing the Hubble tension. This could serve as a compelling motivation to perform further dedicated studies of models featuring the production mechanism put forward in this work.