Dark Matter in the Standard Model?

We critically reexamine two possible Dark Matter candidate within the Standard Model. First, we consider the $uuddss$ exa-quark. Its QCD binding energy could be large enough to make it (quasi) stable. We show that the cosmological Dark Matter abundance is reproduced thermally if its mass is 1.2 GeV. However, we also find that such mass is excluded by the stability of Oxygen nuclei. Second, we consider the possibility that the instability in the Higgs potential leads to the formation of primordial black holes while avoiding vacuum decay during inflation. We show that the non-minimal Higgs coupling to gravity must be as small as allowed by quantum corrections, $|\xi_H|<0.01$. Even so, one must assume that the Universe survived in $e^{120}$ independent regions to fluctuations that lead to vacuum decay with probability 1/2 each.


I. INTRODUCTION
In this work we critically reexamine two different intriguing possibilities that challenge the belief that the existence of dark matter (DM) implies new physics beyond the standard model (SM).

A. DM as the uuddss hexaquark
The binding energy of the hexaquark dibaryon uuddss is expected to be large, given that the presence of the strange quark s allows it to be a scalar, isospin singlet [1], called H or S, and sometimes named sexaquark. A large binding energy might make S light enough that it is stable or long lived. All possible decay modes of a free S are kinematically forbidden if S is lighter than about 1.87 GeV. Then S could be a dark matter candidate within the standard model [2][3][4].
In Sec. II A we use the recent theoretical and experimental progress about tetra-and pentaquarks to infer the mass of the S hexaquark. In Sec. II B we present the first cosmological computation of the relic S abundance, finding that the desired value is reproduced for M S ≈ 1.2 GeV. In Sec. II C we revisit the bound from nuclear stability (NN → SX production within nuclei) at the light of recent numerical computations of one key ingredient: the nuclear wave function [5], finding that S seems excluded.

B. DM as primordial black holes
Primordial black holes (PBH) are hypothetical relics that can originate from gravitational collapse of sufficiently large density fluctuations. The formation of PBH is not predicted by standard inflationary cosmology: the primordial inhomogeneities observed on large cosmological scales are too small. PBH can arise in models with large inhomogeneities on small scales, k ≫ Mpc −1 . PBH as DM candidates are subject to various constraints. BH lighter than 6 × 10 −17 M ⊙ are excluded because of Hawking radiation. BH heavier than 10 5 M ⊙ are safely excluded. In the intermediate region, a variety of bounds makes the possibility that Ω PBH ¼ Ω DM problematic but maybe not excluded-the issue is presently subject to an intense debate. According to [6] DM as PBH with mass M ∼ 10 −15 M ⊙ is not excluded, as previously believed. And the HSC/Subaru microlensing constraint on PBH [7] is partially in the wave optics region. This can invalidate its bound below ∼10 −11 M ⊙ .
Many ad hoc models that can produce PBH as DM have been proposed. Recently the authors of [8] claimed that a mechanism of this type is present within the Standard Model given that, for present best-fit values of the measured SM parameters, the SM Higgs potential is unstable at h > h max ∼ 10 10 GeV [9]. We here critically reexamine the viability of the proposed mechanism, which assumes that the Higgs, at some point during inflation, has a homogeneous vacuum expectation value mildly above the top of the barrier and starts rolling down. When inflation ends, reheating adds a large thermal mass to the effective Higgs potential, which, under certain conditions, brings the Higgs back to the origin, h ¼ 0 [10]. If falling stops very close to the disaster, this process generates inhomogeneities that lead to the formation of primordial black holes. In Sec. III we extend the computations of [8] adding a nonvanishing nonminimal coupling ξ H of the Higgs to gravity, which is unavoidably generated by quantum effects [11]. We find that ξ H must be as small as allowed by quantum effects. Under the assumptions made in [8] we reproduce their results; however, in Sec. III F we also find that such assumptions imply an extreme fine-tuning.
The first mechanism is affected by the observed baryon asymmetry, but does not depend on the unknown physics that generates the baryon asymmetry. The second possibility depends on inflation, but the mechanism only depends on the (unknown) value of the Hubble constant during inflation. In both cases the DM candidate is part of the SM. Conclusions are given in Sec. IV.

II. DM AS THE uuddss HEXAQUARK
The hexaquark S ¼ uuddss is stable if all its possible decay modes are kinematically closed: A stable S is a possible DM candidate. A too light S can make nuclei unstable. Scanning over all stable nuclei, we find that none of them gets destabilized by single S emission if M S > 1.874 GeV, with 6 Li giving the potentially highest sensitivity to M S .

A. Mass of the hexaquark from a diquark model
We estimate the mass of the hexaquark S viewing it as a neutral scalar dibaryon constituted by three spin zero diquarks S ¼ ϵ αβγ ½ud α;s¼0 ½us β;s¼0 ½ds γ;s¼0 ; ð2Þ where α, β, γ are color indices. This is possible thanks to the strange s quark, while spin zero diquarks of the kind ½uu; ½dd; ½ss are forbidden by Fermi statistics because of antisymmetry in color and spin. We assume the effective Hamiltonian for the hexaquark [12] H ¼ X i≠j¼fu;d;sg where the κ ij are effective couplings determined by the strong interactions at low energies, color factors, quark masses, and wave functions at the origin. The m ij are the masses of the diquarks in S made of i and j constituent quarks [13]. S i is the spin of the ith quark. Another important assumption, which is well motivated by studies on tetraquarks [12], is that spin-spin interactions are essentially within diquarks and zero outside, as if they were sufficiently separated in space.
Considering diquark masses to be additive in the constituent quark masses, and taking q and s constituent quark masses from the baryons, one finds The chromomagnetic couplings κ ij could be derived as well in the constituent quark model using data on baryons However, it is known that to reproduce the masses of light scalar mesons, interpreted as tetraquarks, σð500Þ; f 0 ð980Þ; a 0 ð980Þ; κ [14], we need Spin-spin couplings in tetraquarks are found to be about a factor of 4 larger compared to the spin-spin couplings among the same pairs of quarks in the baryons, which also make diquarks. It is difficult to assess if this would change within a hexaquark. At any rate we can attempt a simple mass formula for S, which in terms of light tetraquark masses means Using the determination of chromomagnetic couplings from baryons we would obtain whereas keeping the chromomagnetic couplings needed to fit light tetraquarks gives if the same values for the chromomagnetic couplings to fit light tetraquark masses are taken [or 1.4 GeV using (6)]. There is quite a lot of experimental information on tetraquarks [12], whereas hexaquarks, for the moment, are purely hypothetical objects. On purely qualitative grounds we might expect that the mass of S could be closer to the heavier value being a dibaryon and not a dimeson (tetraquark) like light scalar mesons.
In the absence of any other experimental information it is impossible to provide an estimate of the theoretical uncertainty on M S .
Lattice computations performed at unphysical values of quark masses find small values for the S binding energy, about 13, 75, 20 MeV [15][16][17]. Extrapolations to physical quark masses suggest that S does not have a sizable binding energy; see, e.g., [18]. Furthermore, the binding energy of the deuteron is small, indirectly disfavoring a very large binding energy for the (somehow similar) S, which might too be a moleculelike state. 1 Despite this, we overoptimistically treat M S as a free parameter in the following.
We also notice that the S particle could be much larger than what is envisaged in [2][3][4] and that its coupling to photons, in the case of 2-3 fm size (see the considerations on diquark-diquark repulsion at small distances in [19]), could be relevant for momentum transfers k as small as k ∼ 60 MeV, compared to k ≈ 0.5 GeV considered by Farrar.

B. Cosmological relic density of the hexaquark
We here compute the relic density of S dark matter, studying if it can match the measured value Ω DM h 2 ¼ 0.1186, i.e., Ω DM ≈ 5.3 Ω b [20]. A key ingredient of the computation is the baryon asymmetry. Its value measured in cosmic microwave background (CMB) and Big Bang Thereby the baryon asymmetry before S decoupling must be One needs to evolve a network of Boltzmann equations for the main hadrons: p, n, Λ, Σ 0 , Σ þ , Σ − , Ξ 0 , Ξ − , and S. Strange baryons undergo weak decays with lifetimes τ ∼ 10 −10 s, a few orders of magnitude faster than the Hubble time. This means that such baryons stay in thermal equilibrium. We thereby first compute the thermal equilibrium values taking into account the baryon asymmetry. Thermal equilibrium implies that the chemical potentials μðTÞ satisfy Their overall values are determined imposing that the total baryon asymmetry equals The equilibrium values can be analytically computed in Boltzmann approximation (which becomes exact in the nonrelativistic limit) where the þ (−) holds for (anti)particles. We then obtain the abundances in thermal equilibrium plotted in Fig. 1, assuming M S ¼ 1.2 GeV. We see that the desired S abundance is reproduced if the interactions that form/destroy S decouple at T dec ≈ 25 MeV. This temperature is so low that baryon antiparticles have negligible abundances, and computations can more simply be done neglecting antiparticles. 2 Then, the desired decoupling temperature is simply estimated imposing n eq S =n eq p ∼ Y BS e ð2M p −M S Þ=T ∼ 1, and it decreases if M S is heavier: To compute the decoupling temperature, we consider the three different kinds of processes that can lead to the formation of S: (1) Strong interactions of two heavier QCD hadrons that contain the needed two s quarks. One example is ΛΛ ↔ SX, where X denotes pions. These are doubly Boltzmann suppressed by e −2m s =T at temperatures T < m s . (2) Strong interactions of one heavier strange hadron and weak ΔS ¼ 1 interactions that form the other s (as dū → sū) from lighter hadrons. One example is pΛ ↔ SX. These are singly Boltzmann suppressed by e −m s =T and by G 2 F Λ 4 QCD ∼ 10 −10 . (3) Double-weak interactions that form two s quarks starting from lighter hadrons. One example is pp → SX. These are doubly suppressed by ðG 2 F Λ 4 QCD Þ 2 . At T ∼ 25 MeV the abundance of strange hadrons is still large enough that QCD processes dominate over EW processes: interactions that form and destroy S proceed dominantly through QCD collisions of strange hadrons, where X can be a π 0 or a γ, as preferred by approximate isospin conservation. The Λ can be substituted by the Σ 0 . Defining z ¼ M S =T and Y p ¼ n p =s, the Boltzmann equation for the S abundance is where the superscript "eqb" denotes thermal equilibrium at fixed baryon asymmetry and Y B is summed over all baryons but the dibaryon S. A second equation for Y B is not needed, given that the baryon number is conserved: Furthermore, YS is negligible, and YB is negligible around decoupling. The S production rate is obtained after summing over all bb 0 ↔ SX processes of Eq. (18). In the nonrelativistic limit the interaction rate gets approximated as The opposite process is more conveniently written in terms of the S breaking width defined by γ eqb S ðTÞ ¼ n eqb S Γ eqb S and given by This gets Boltzmann suppressed at T ≲ M Λ − M p ∼ m s , when hyperons disappear from the thermal plasma. Assuming Y S ≲ Y B , the Boltzmann equation is approximatively solved by Y S ∼ Y eqb S evaluated at the decoupling epoch where Γ eqb S ∼ H, which corresponds to T dec ∼ m s = ln Y B . This leads to the estimated final abundance The fact that S is in thermal equilibrium down to a few tens of MeV means that whatever happens at higher temperatures gets washed out. Notice the unusual dependence on the cross section for S formation: increasing it delays the decoupling, increasing the S abundance. Figure 2 shows the numerical result for the relic S abundance, computed by inserting in the Boltzmann equation a s-wave σ bb 0 v rel ¼ σ 0 , varied around 1=GeV 2 . The cosmological DM abundance is reproduced for M S ≈ 1.2 GeV, while a large M S gives a smaller relic abundance. Bound-state effects at BBN negligibly affect the result and, in particular, do not allow one to reproduce the DM abundance with a heavier M S ≈ 1.8 GeV.
We conclude this section with some sparse comments. Possible troubles with bounds from direct detection have 2 Let us consider, e.g., the process Λ þ Λ ↔ S þ X where X denotes other SM particles that do not carry the baryon asymmetry, such as pions. Thermal equilibrium of the above process implies Inserting n Λ ∼ n p e ðM p −M Λ Þ=T with n p ≈ Y BS s gives Namely, n S ≪ n p at large T; n S ≫ n p at low T. A DM abundance comparable to the baryon abundance is obtained only if reactions that form S decouple at the T in Eq. (17).
been pointed out in [4,21,22]: a DM velocity somehow smaller than the expected one can avoid such bounds reducing the kinetic energy available for direct detection.
Using a target made of antimatter (possibly in the upper atmosphere) would give a sharp annihilation signal, although with small rates. The magnetic dipole interaction of S does not allow one to explain the recent 21 cm anomaly along the lines of [23] (an electric dipole would be needed). The interactions of DM with the baryon/photon fluid may alter the evolution of cosmological perturbations leaving an imprint on the matter power spectrum and the CMB. However, they are not strong enough to produce significant effects. The S particle is electrically neutral and has spin zero, such that its coupling to photons is therefore suppressed by powers of the QCD scale [3]. So elastic scattering of S with photons is not cosmologically relevant. A light S would affect neutron stars, as they are expected to contain Λ particles, made stable by the large Fermi surface energy of neutrons. Then, ΛΛ → S would give a loss of pressure, possibly incompatible with the observed existence of neutron stars with mass 2.0M sun [1]. However, we cannot exclude S on this basis, because production of Λ hyperons poses a similar puzzle. S as DM could interact with cosmic ray p giving and photon and other signals [24] and would be geometrically captured in the Sun, possibly affecting helioseismology. 3 In the next section we discuss the main problem, which seems to exclude S as DM.

C. Super-Kamiokande bound on nuclear stability
Two nucleons N ¼ fn; pg inside a nucleus N can make a double weak decay into S, emitting π, γ, or e [2]. This is best probed by Super-Kamiokande (SK), which contains ∼8 × 10 32 oxygen nuclei. No dedicated search for 16 O 8 → N 0 SX (where X can be one or two π, γ, e and N 0 can be 14 O 8 , 14 N 7 , 14 C 6 , depending on the charge of X) has been performed, 4 but a very conservative limit is obtained by requiring the rate of such transitions to be smaller than the rate of triggered background events in SK, which is about 10 Hz [27]. A more careful analysis would likely improve this bound by 3 orders of magnitude [2]. The amplitude for the formation of S is reasonably dominated by the sample diagram in Fig. 3: doubly weak production of two virtual strange Λ Ã baryons (e.g., through p → π þ Λ Ã and n → π 0 Λ Ã ; at quark level u → sdu and d → sūu), followed by the strong process Λ Ã Λ Ã → S, The predicted lifetime is then obtained as where the smaller value holds if S is so light that the decay can proceed through real π or ππ emission, while the longer lifetime is obtained if instead only lighter e þ ν or γ can be emitted.
The key factor is the dimensionless matrix element FIG. 3. Sample diagram that dominates nucleon decay into S inside nuclei. The initial state can also be nn or pp. 3 We thank M. Pospelov for suggesting these ideas. 4 SK searched for dinucleon decays into pions [25] and leptons [26] and obtained bounds on the lifetime around ∼10 32 years. However, these bounds are not directly applicable to 16 O 8 → N 0 SX where the invisible S takes away most of the energy reducing the energy of the visible pions and charged leptons, in contrast to what is assumed in [25,26]. 5 A numerical factor of 1440 due to spin and flavor effects has already been factored out from jMj 2 Λ Ã Λ Ã →S here and in the following. Note also that the threshold 2M N − M π ¼ 1.74 GeV neglects the small difference in binding energy between N and N 0 .
which we now discuss. Following [2], we assume that the initial state wave function can be factorized into wave functions of the two Λ Ã baryons and a relative wave function ψ nuc ð⃗ aÞ for the separation ⃗ a between the center of mass of the Λ Ã 's. The matrix element is given by the wave-function overlap Here, ⃗ ρ a;b , ⃗ λ a;b are center-of-mass coordinates that parametrize the relative positions of the quarks within each Λ Ã . Using the Isgur-Karl (IK) model [28] the wave functions for the quarks inside the Λ Ã and inside the S are approximated by where r N and r S are the radii of the nucleons, respectively, of S. 6 Performing all integrals except the final integral over a ≡ j⃗ aj gives As shown in Fig. 5 below (and as discussed in [2]), if r S is not much smaller than r N , the overlap integral is not very much suppressed and τð 16 O 8 → N 0 SXÞ is tens of orders of magnitude below the experimental limit and is clearly excluded. This conclusion is independent of the form of ψ nuc . However, if r S were a few times smaller than r N -a possibility that seems unlikely due to diquark repulsions (see, e.g., [19]) but cannot firmly be excluded-then τð 16 O 8 → N 0 SXÞ is extremely sensitive to the probability of the overlap of two nucleons inside the oxygen core at very small distances (less than, say, 0.5 fm). The wave function of nucleon pairs ψ nuc at such small distances has not been probed experimentally. In fact, at such small distances nucleons are not the appropriate degrees of freedom. 7 Thus, for a very small S one can only make an educated guess of τð 16 O 8 → N 0 SXÞ, since the form of ψ nuc is uncertain. Nevertheless, we will show in the following that for a reasonable form of ψ nuc a stable S is excluded even if it were very small.
Numerical computations of the ground-state wave functions of nuclei, including 16 O 8 , have been performed, e.g., in [5]. The quantity that determines ψ nuc is the two-nucleon point density ρ NN , defined in Eq. (58) of [5]. We obtain ρ NN ðaÞ by interpolating the data given in [5] and adding the constraint ρ NN ð0Þ ¼ 0, which is a conservative assumption for our purposes since ρ NN ð0Þ ≠ 0 would lead to a larger matrix element. There are 28 neutron-neutron pairs and 64 proton-neutron pairs in 16 O 8 so one has R da4πa 2 ρ nn ðaÞ ¼ 28 and R da4πa 2 ρ pn ðaÞ ¼ 64. We therefore define the wave functions These wave functions are plotted in Fig. 4, together with the Miller-Spencer (MS) and the Brueckner-Bethe-Goldstone (BBG) wave function used in [2]. The BBG wave functions assume a hard repulsive core between nucleons such that ψ nuc vanishes at a < r core . We take r core ¼ 0.5 fm for illustration. This is not realistic but allows one to see what kind of nuclear wave function would sufficiently suppress the rate of S formation in nuclei, if S is small enough. The resulting jMj 2 Λ Ã Λ Ã →S is plotted in Fig. 5, again compared to that obtained using the Miller-Spencer and BBG wave functions. 6 One should be aware that the IK model has serious shortcomings. One issue is that it is a nonrelativistic model-an assumption that is problematic in particular for small S. Another problem is that the value of r N that gives a good fit to the lowest lying 1 2 þ and 3 2 þ baryons-r N ¼ 0.49 fm-is smaller than the charge radius of the proton: r N ¼ 0.87 fm. Therefore we consider both r N ¼ 0.49 fm and r N ¼ 0.87 fm, as done in [2]. 7 Data indicate that about 20% of the nucleons form pn pairs so close (about 1 fm) that the local density reaches the nucleon density (about 2.5 times larger than the nuclear density) and thus that the quark structure of nucleons starts becoming relevant already at a ∼ 1 fm [29].
The resulting matrix elements from the MS wave function qualitatively agree to what is obtained using the wave functions extracted from [5]. By contrast, the matrix element using the BBG wave function with hard core radius r core ¼ 0.5 fm is very much suppressed, especially if S is small. The reason is that, according to the assumption of a hard core repulsive potential, the nucleons cannot get close enough to form the small state S. Since we do not consider a ψ nuc ðaÞ that vanishes for a ≲ 0.5 fm realistic, we conclude that a stable S is excluded.
Weaker bounds on S production are obtained considering baryons containing Λ's.

III. DM AS BLACK HOLES TRIGGERED BY HIGGS FLUCTUATIONS
We here present the technical computations relative to the mechanism anticipated in the Introduction. The SM potential is summarized in Sec. III A. In Sec. III B we outline the mechanism that generates black holes. Section III C studies the generation of Higgs inhomogeneities. Postinflationary dynamics is studied in Sec. III D. The formation of black holes is considered in Sec. III E. The viability of a critical assumption is discussed in Sec. III F.

A. The Higgs effective potential
The effective potential of the canonically normalized Higgs field during inflation with Hubble constant H 0 is at h ≫ 174 GeV. Here λ eff is the effective quartic coupling computed including quantum corrections. The second mass term in V eff ðhÞ can be generated by various different sources [8]. We consider the minimal source: a Higgs coupling to gravity, during inflation. Finally, during inflation the effective potential in Eq. (31) is augmented by the vacuum energy associated with the inflaton sector, 18 GeV is the reduced Planck mass.
We implement the renormalization group (RG) improvement of the effective potential at next-to-next-to-leading order (NNLO) precision: running the SM parameters at three loops and including two-loop quantum corrections to the effective potential. We consider fixed values of α 3 ðM Z Þ ¼ 0.1184 and M h ¼ 125.09 GeV, and we vary the main uncertain parameter, the top mass, in the interval M t ¼ 172.5 AE 0.5 GeV [30]. In Fig. 6(a) we show the resulting λ eff ðhÞ as a function of h.  The darker (lighter) green curve shows ψ nn nuc , obtained from [5], using the AV18 (AV18+UIX) potential; the darker (lighter) red curve shows ψ np nuc using the AV18 (AV18 +UIX) potential; the blue curve is the MS wave function used in [2] and the dashed black curve is the BBG wave function with hard core radius r core ¼ 0.5 fm. Left: Wave functions over the entire range. Right: Zoom-in plot of the region relevant for our calculation.
The nonminimal coupling to gravity ξ H receives SM quantum corrections encoded in its renormalization group equation (RGE), which induce ξ H ≠ 0 even starting from ξ H ¼ 0 at some energy scale. The RGE running of small values of ξ H ðμÞ is shown in Fig. 6(b). As mentioned before, a nonzero ξ H can be considered as a proxy for an effective mass term during inflation. The latter, for instance, can be generated by a quartic interaction λ hϕ jHj 2 ϕ 2 between the Higgs and the inflaton field or by the inflaton decay into SM particles during inflation. For these reasons, it makes sense to include ξ H as a free parameter in the analysis of the Higgs dynamics during inflation, at most with the theoretical bias that its size could be loop suppressed.

Analytic approximation
We will show precise numerical results for the SM case. However, the discussion is clarified by introducing a simple approximation that encodes the main features of the SM effective potential in Eq. (31), where h cr is the position of the maximum of the potential with no extra mass term, ξ H ¼ 0. The parameters b and h cr depend on the low-energy SM parameters such as the top mass: they can be computed by matching the numerical value of the Higgs effective potential at the gauge-invariant position of the maximum, V eff ðh cr Þ ¼ bh 4 cr =8. The result is shown in the right panel of Fig. 7.
Results will be better understood when presented in terms of the dimensionless parameters b, ξ H , h cr =H 0 , and T=H 0 , where T is the temperature, as they directly control the dynamics that we are going to study. The parameter b controls the flatness of the potential beyond the potential barrier at h cr , with smaller b corresponding to a flatter potential. The nonminimal coupling ξ H controls the effective Higgs mass during inflation. Finally,M Pl =H 0 will set the reheating temperature in Eq. (35) and thus the position and size of the thermal barrier. The position of the potential barrier-defined by the field value where the effective potential has its maximumstrongly depends on the value of the top mass, on the nonminimal coupling to gravity, and, after inflation, on the temperature of the thermal bath that provides an extra mass term. For ξ H ≠ 0, the maximum of the Higgs potential gets shifted from h cr to where WðzÞ is the product-log function defined by must be satisfied; otherwise, the effective mass is too negative, and it erases the potential barrier, thus leading to a classical instability.

The thermal potential
After the end of inflation, the Higgs effective potential receives large thermal corrections from the SM bath at generic temperature T. The initial temperature of the thermal bath is fixed by the dynamics of reheating after inflation. We assume instantaneous reheating, as this is most efficient for rescuing the falling Higgs field. The reheating temperature is then given by where g Ã ¼ 106.75 is the number of SM degrees of freedom. After reheating, the Universe becomes radiation dominated, the Ricci scalar vanishes, and so does the contribution to the effective potential from the nonminimal Higgs coupling to gravity. The effective Higgs potential at finite temperature is obtained adding an extra thermal contribution V T , which can be approximated as an effective thermal mass for the Higgs field, M 2 T ≃ 0.12T 2 (see, e.g., [10]) At h ≲ T we can neglect the exponential suppression in the thermal mass, and the maximum of the effective potential in Eq. (36) is given by B. Outline of the mechanism During inflation, the Higgs field is subject to quantum fluctuations. Depending on the value of H 0 , these quantum fluctuations could lead the Higgs beyond the barrier and make it roll toward Planckian values. If T RH is high enough and h is not too far, thermal corrections can "rescue" the Higgs, bringing it back to the origin [10]. The mechanism relies on a tuning such that the following situation occurs [8]: (i) At N in ∼ 20 e-folds before the end of inflation, the Higgs background value h is brought by quantum fluctuation to some h in ≠ 0. This configuration must be spatially homogeneous on an inflating local patch large enough to encompass our observable Universe today. We consider the de Sitter metric in flat slicing coordinates, ds 2 ¼ −dt 2 þ a 2 ðtÞd⃗ x 2 . We will discuss later how precisely this assumption must be satisfied, as well as its plausibility. (ii) When the classical evolution prevails over the quantum corrections, the Higgs field, starting from the initial position h in , begins to slow roll down the negative potential. This condition reads where c is a constant of order 1, fixed to c ¼ 1 in [8].
We will explore what happens choosing c ¼ 0.9 or c ¼ 1.1. From this starting point t in on, the classical evolution of the background Higgs value is described bÿ where the subscript cl indicates that this is a classical motion. Dots indicate derivatives with respect to time t. (iii) At the end of inflation, N end ¼ 0, the Higgs is rescued by thermal effects. This happens if the value of the Higgs field h end at the end of inflation is smaller than the position of the thermal potential barrier at reheating, h T RH max . A significant amount of PBH arises only if this condition is barely satisfied in all the Universe. This is why the homogeneity assumption in (i) is needed. To compute condition (iii) we fix the initial value of the classical motion h in such that Eq. (38) is satisfied with c ¼ 1; next, we maximize the h end obtained solving Eq. (39) by tuning the amount of inflation where the fall happens, as parametrized by N in . The left panel of Fig. 7 shows the initial value N in obtained following this procedure as a function of h cr in units of H 0 . Smaller values of b (i.e., smaller values of M t ) imply a flattening of the potential, and the classical dynamics during inflation is slower. The right side of the curves is limited by the classicality condition in Eq. (38). A ξ H < 0 shifts the position of the potential barrier toward the limiting value h T max in Eq. (37)-which does not depend on ξ H -above which the rescue mechanism due to thermal effects is no longer effective: its net effect is to reduce the number of e-folds during which classical motion can happen (for fixed b).
We anticipate here the feature of PBH formation, which implies the restriction on the parameter space mentioned at point (i): the Higgs fall must start at least N in ∼ 20 e-folds before the end of inflation. The collapse of the mass inside the horizon N e-folds before the inflation end forms a PBH with mass (see also Sec. III E), where F → 1 at M PBH > 10 17 g. BH heavier than M PBH > 10 15 g are cosmologically stable, and BH heavier than M PBH > 10 16.5 g are allowed by bounds on Hawking radiation as a (significant fraction of) DM. Since N < N in , imposing M PBH > 10 16.5 g implies a conservative lower limit on N in ,

C. Higgs fluctuations during inflation
We now consider the evolution of Higgs perturbations during inflation. Expanding h in Fourier space with comoving wave number k, 8 the equation for the mode δh k takes the form where we neglected metric fluctuations. In terms of the number of e-folds N and of the Mukhanov-Sasaki variable u k ≡ aδh k , it becomes It is convenient to consider the evolution of the perturbation making reference to a specific moment before the end of inflation: at the initial value N in defined in Sec. III B. We recall that in our convention N end ¼ 0 at the end of inflation. Equation (44) becomes In this form, the Mukhanov-Sasaki equation is particularly illustrative. Consider the evolution of the perturbation for a mode of interest k that we fix compared to the reference value a in H 0 at t ¼ t in . In particular, we consider the case of a mode k that is subhorizon at the beginning of the classical evolution, that is, k ≫ a in H 0 . From Eq. (45), we see that in the subsequent evolution with N < N in the exponential suppression will turn the mode from subhorizon to superhorizon.
We are now in the position to solve Eq. (45). To this end, we need boundary conditions for u k and its time derivative. We use the Bunch-Davies conditions at N ¼ N in for modes that are subhorizon at the beginning of the classical evolution, k=a in H 0 ≥ 1, and we treat the real and imaginary parts of u k separately since they behave as two independent harmonic oscillators for each comoving wave number k. At generic e-fold time N, the perturbation δh k is related to the Mukhanov-Sasaki variable u k by We show in the left panel of Fig. 8 our results for the time evolution of the classical background h cl and the perturbation δh k (both real and imaginary parts) during the last 20 e-folds of inflation. As a benchmark value, we consider an initial subhorizon mode with k=a in H 0 ¼ 100. After a few 8 The comoving wave number k ¼ j ⃗ kj is time independent, and it is related to the physical momentum via k phys ¼ k=aðtÞ, which decreases as the space expands. e-folds of inflation, such a mode exits the horizon: oscillations stop, and from this point on, further evolution is driven by the time derivative of the classical background. This is a trivial consequence of the equations of motion on superhorizon scales. Differentiating Eq. (39) with respect to the cosmic time shows that _ h cl and δh k satisfy the same equation on superhorizon scales, and, therefore, they must be proportional, δh k ¼ CðkÞ _ h cl for k ≪ aH 0 [8]. The proportionality function CðkÞ can be obtained by a matching procedure. Deep inside the horizon, in the limit k ≫ aH 0 , the Mukhanov-Sasaki variable u k ≡ aδh k reproduces the preferred vacuum of a harmonic oscillator in flat Minkowski space, and we have, after introducing the conformal time τ as dt ¼ adτ, u k ¼ e −ikτ = ffiffiffiffiffi 2k p . Roughly matching the absolute value of the solutions at horizon crossing, we determine the absolute value of CðkÞ as where we indicate with t k the time of the horizon exit for the mode k-the time at which k ¼ aðt k ÞH 0 ≡ a k H 0 (equivalently, −kτ k ¼ 1). The number of e-folds at the horizon exit is given by

Primordial curvature perturbations
The primordial curvature perturbation ζð⃗ x; tÞ on uniform energy density slices ρ is defined (at the linear order) by the perturbed line element [31,32] and it is related to the (total) energy density perturbation δρ and to the curvature perturbation on a generic slicing Ψ by the gauge invariant formula where we introduced the conformal time by means of dt ¼ adτ and H ¼ aH where H ¼ a 0 =a and prime 0 indicates the derivative with respect to τ. The virtue of this definition is that, on superhorizon scales, ζ is practically identical to the comoving curvature perturbation R defined on hypersurfaces of constant comoving time. Furthermore, on superhorizon scales ζ is conserved provided that the pressure perturbation is adiabatic. In conventional singlefield inflation models, purely adiabatic perturbations are generated due to quantum fluctuations of the single field driving inflation. The considered setup differs from a conventional scenario because of the presence, in addition to the inflaton, of the Higgs field [8]. During inflation, the curvature perturbation on uniform energy density slices is given, in a spatially flat gauge, by where we separated the inflaton component ζ ϕ and the Higgs component ζ h ≡ H 0 δρ h =_ ρ h [8]. We assume that there is no energy transfer between the Higgs and the inflaton sectors, and that the latter generates, on large scales, the perturbations responsible for the CMB anisotropies and large scale structures. Given these assumptions, ζ ϕ and ζ h are separately conserved on superhorizon scales [33]. As customary, we can decompose ζ h in Fourier modes. For a given comoving wave number k we have, in a spatially flat gauge, where the last analytical approximation is valid for the absolute value of ζ ðkÞ h in the superhorizon limit, and where we used

D. Higgs dynamics after inflation
Assuming instantaneous reheating, the energy density of the inflaton is instantaneously converted into radiation at the end of inflation. The inflaton energy density ρ R is related to the temperature of the thermal bath by The dynamics after reheating is described by the following system of coupled Higgs-radiation equations: where the energy density of the Higgs field is given by and the Hubble parameter is related to the total energy density by The damping factor Γ takes into account the Higgs decays at temperature T and represents the energy transportation from the Higgs field to radiation. We use the expression quoted in [8], Γ ≃ 10 −3 T. In the following, we shall adopt the assumption of sudden decay approximation. In this approximation-corresponding to Γ ¼ 0 in Eqs. (54)-Higgs and radiation evolve separately, and the Higgs decay occurs instantaneously at H ¼ Γ. The system in Eqs. (54) can be solved using the following boundary conditions. As far as the classical Higgs field is concerned, we use the field values at the end of inflation computed in Sec. III C. The evolution of ρ R , in contrast, starts from the temperature in Eq. (35). 9 Using the solution for ρ R , Eq. (53) gives the time evolution of the temperature. After instantaneous reheating, the Higgs potential suddenly changes from the expression in Eq. (32) to the one in Eq. (36), and interactions with the SM bath generate a large thermal mass. If h end < h T max , the fall of the Higgs field is rescued by thermal corrections, and its background value starts oscillating around the minimum at zero until the decay becomes efficient at H ¼ Γ. At this stage, the Higgs component ζ h of the curvature perturbation is not conserved compared to the value computed in Eq. (52) during inflation. This is because the interactions with the SM plasma that are responsible for the appearance of the thermal mass introduce a nonadiabatic component in the pressure perturbation.
We compute-in a spatially flat gauge, and for a given comoving wave number k-the total curvature perturbation after reheating according to where radiation perturbations are set to zero. We compute ζ ðkÞ numerically at the time of Higgs decay, H ¼ Γ. After Higgs decay, radiation remains as the only component of the Universe and ζ ðkÞ is, therefore, conserved.

E. The power spectrum and the PBH abundance
The curvature power spectrum is given by A numerical example is shown in Fig. 9. Its left panel shows that a small ξ H ¼ −10 −3 has only a minor effect. The cut in the power spectrum at small k arises because of assumption (i): before that classical rolling starts, the Higgs 9 More precisely, conservation of total energy determines the reheating temperature as The approximation in Eq. (35) is valid because the inflaton energy density dominates over the Higgs contribution.
field is away from its minimum and exactly homogeneous. Relaxing this assumption would increase the power spectrum at small k.
The right panel shows the significant effect of a small change in the arbitrary order one parameter c that defines the classicality condition in Eq. (38). Reducing c anticipates the initial moment t in (or equivalently N in ) where the Higgs starts to roll down the potential. As a consequence, earlier quantum fluctuations get taken into account by our computation of Sec. III C.
Finally, we can now compute the mass and amount of PBH generated by Higgs fluctuations. The radius of the Hubble horizon or the wavelength of the modes determines the typical mass of the PBHs [34], where N is the number of e-folds when the k mode leaves the horizon; γ ≈ 0.2 is a correction factor [35]. For example, H 0 ¼ 10 12 GeV and N ¼ 20 To compute the fraction of the Hubble volume collapsing to PBHs, we need the variance of the smoothed density perturbation over a radius R, RÞ is the density fluctuation smoothed by a window function Wðx; RÞ, assumed to be Gaussian, The variance can be computed in terms of density power spectrum P δ ðkÞ, which is related to the curvature perturbation ζ power spectrum P ζ ðkÞ as [36] whereWðkÞ ¼ exp ð−k 2 =2Þ is the Fourier transform of the window function. Assuming that PBH are formed when the density perturbation exceeds δ th ≈ 0.5 [37], the fraction of the Universe ending up in PBH is given by the tail of the assumed Gaussian distribution, The latter approximation is relevant, given that the obtained power spectra P ζ are of order 10 −2 . We convert σðRÞ to a function σðM PBH Þ, taking into account that the size R is related to the mass M PBH as R ≈ 2GM PBH =γa. The fraction of PBH relative to the DM abundance at given mass, f PBH ðM PBH Þ, is given by [36] f PBH ≈ 2.7 × 10 8 10.75 The distribution of PBH as a function of their mass M PBH is strongly peaked at the value that maximizes σðMÞ, in view of the exponential factor in Eq. (63). In terms of the power spectrum P ζ this means that the PBH abundance roughly scales as e −1=P ζ and is dominated by the top of the peak of P ζ ðkÞ. The PBH mass distribution is peaked at the value corresponding to the k that maximizes P ζ ðkÞ. The desired abundance is reproduced for max P ζ ∼ 10 −2 , and slightly larger (smaller) values produce way too many (too few) PBH. This means that, due to the fine-tuned nature of the mechanism (as in all models that can generate PBH), the amount and mass of PBH depend in an extremely sensitive way on the uncertain SM and cosmological parameters, mainly the top mass and the Hubble constant H 0 . Furthermore, uncertainties in the computation of black hole formation (such as the value of δ th ) imply uncertainties of many orders of magnitude in the PBH density.
For example, a 10% variation in the order one arbitrary constant c that parametrizes the classicality condition has an order one impact on the power spectrum (Fig. 9), and consequently a huge impact on the PBH abundance.
In addition to SM parameters and to H 0 , the PBH abundance also depends on two extra parameters that characterize the assumed initial condition: constant hðxÞ ¼ h in at N in e-folds before the end of inflation, when the approximated classicality condition is satisfied. The next section reconsiders if this assumption is justified.

F. Is a homogeneous Higgs background a sensible assumption?
The computation was based on assumption (i): at N in ≈ 20 e-folds before the end of inflation, the Higgs field must be away from its minimum and constant within the presently observable horizon.

Approximate Higgs homogeneity?
We here show that approximate homogeneity is a natural product of inflation, provided that the SM Higgs potential satisfies certain conditions. Quantum corrections in inflationary (de Sitter) spacetime have been studied in [38], which showed that long wavelength fluctuations can be described by a Fokker-Planck equation for ρðh; NÞ, the probability of finding the Higgs field at the value h at N e-folds of inflation, The first term on the right-hand side is a diffusion term due to quantum fluctuations while the second term is a drift (or transport) term due to the potential. After some e-folds, the distribution converges to its equilibrium value ρ eq ðhÞ ∝ expð−8π 2 V eff ðhÞ=3H 4 Þ [38]. We are interested in the probability of having a roughly constant Higgs away from its minimum h ¼ 0 within the presently observable horizon. This configuration is a natural outcome of inflationary dynamics provided that the correlation length of Higgs fluctuations is larger than the present horizon. Following [38], the computation of the correlation length is simplified, observing that correlation functions depend on space and time separations respecting the O(4,1) invariance of de Sitter. Thereby we can compute the evolution of Higgs fluctuations at a fixed point in space and study the correlation as a function of time t, or equivalently as a function of the number of e-folds N. At a large time separation the correlation is well approximated by its dominant exponential and parametrized in terms of a correlation time τ corr or equivalently in terms of number of e-folds N corr ¼ H 0 τ corr as [38] O(4,1) invariance implies that the spatial correlation length is exponentially large e HN corr , namely that space is exponentially inflated. We demand that N corr is larger than about 40 in order to produce a smooth region as large as our Universe at N in ≈ 20 e-folds before the end of inflation. Before computing numerically N corr for the SM potential, it is useful to consider some simple limits: (0) A massless free scalar h is a simple but special case, because it does not "thermalize" to the equilibrium distribution. Rather, it undergoes random walk diffusion. Starting from h ¼ 0, after N e-folds one has hh 2 i 1=2 ¼ ffiffiffiffi N p H 0 =2π and the correlation length in a region with given hhi is N corr ¼ 2πhhi=H. Imposing N corr > 40 implies H 0 < 0.16hhi. (m) A free scalar with squared mass m 2 > 0 fluctuates as [38]. For m 2 ¼ −12ξ H H 2 one gets N corr ¼ −1=4ξ H > 40 for −0.006 < ξ H < 0. Roughly, this will be our final result.
036. This is satisfied in the SM at large energy. The above two results agree, up to order one factors, approximating a generic VðhÞ as a massive field with m 2 ¼ V 00 . The precise value of the correlation length can be obtained as proposed in [38]: the dominant exponential that solves Eq. (65) can be computed as the eigenvalue of a Schrödinger-like equation. In the SM λ eff ðhÞ runs to negative values, making the curvature V 00 negative at large field values; the patch with a large correlation length is created while h is climbing the potential in its region with positive curvature. The result is shown in Fig. 10, and basically agrees with the result anticipated at point (m): ξ H must be negative and small. Such a small ξ H is roughly compatible with the size of quantum corrections to ξ H . The running Higgs quartic coupling is small enough that it does not spoil the mechanism, as qualitatively understood at point (λ).
In conclusion, a roughly homogeneous Higgs field (up to fluctuation of order H 0 ) encompassing our whole horizon is a natural outcome of inflation, provided that ξ H is small enough. In a multiverse context, its fine-tuned value that leads to DM as PBH can be justified on an anthropic basis. Figure 11(a) shows that the initial homogeneous Higgs fields hðxÞ ¼ h in at N ¼ N in e-folds before the end of inflation must be tuned to a part in about 10 −3 in order to produce a final h end close to the maximal value that can be rescued by reheating, as needed to produce a substantial amount of primordial black holes. An interesting PBH abundance is obtained within the narrow strip in Fig. 11(a). Its boundaries have been computed as follows:

Exact Higgs homogeneity?
(i) A slightly smaller h in leads to a negligible PBH amount. In Fig. 11(a) we show the peak value of the power spectrum and shade in gray regions where it is below 10 −2 . (ii) A slightly larger h in leads to a too large h end not rescued by reheating. In Fig. 11(a) we shade in red regions where h end is above the maximal value that can be rescued by instantaneous reheating. In Fig. 11(b) we show that a variation in _ h in has a smaller effect as it gets redshifted away.
Within the assumption that h in is homogeneous, its tuned value can be justified on an anthropic basis provided that PBH make all DM [8,39]. However, one single fluctuation δh in ≳ 10 −3 h in away from the assumed perfect homogeneity can lead to one vacuum decay bubble that, after inflation, expands engulfing the observable Universe.
As discussed previously, inflation can produce an approximate homogeneity within a large patch, but up to fluctuations of order δh in ∼ H 0 =2π. 10 Additional fluctuations in _ h in are less significant and we ignore them. The Hubble constant H 0 cannot be significantly reduced, for two reasons: First, H 0 ≳ h cr is needed to keep the Higgs fluctuating around the top of its potential barrier, until it starts to classically roll down. Second, it would suppress the unwanted prefall fluctuations δh in at the price of (a) (b) FIG. 11. We fix the SM to M t ¼ 172 GeV [which corresponds to h cr ¼ 4 × 10 12 GeV, b ¼ 0.09=ð4πÞ 2 ] and H 0 ¼ 10 12 GeV and explore the dependence on the parameters that define the assumed initial condition. Left: Final outcome as a function of the homogeneous initial value of the Higgs field h in at N in e-folds before the end of inflation and of N in . The desired PBH abundance is produced inside the white band. Inflationary fluctuations in h in are much larger than the needed tuning. Right: We fix N in and introduce dh in =dN as an additional parameter. 10 Technically, such fluctuations appear in our equations when (motivated, e.g., by the arbitrariness in the classicality condition) we start from a different initial time: this shifts the position and the height of the presumed peak of the power spectrum, and thereby the abundance and mass of PBHs. suppressing postfall fluctuations that generate PBHs; see Eq. (52).
Thereby we must take into account the effect of fluctuations in the initial hðxÞ. Pictorially, the status of the Universe should now be represented in Fig. 10 not by a point (which can lie in the desired region), but by a dot, much thicker than the desired region.
This is a problem because, at N in ≈ 20 e-folds before the end of inflation (which needs at least 60 e-folds), the present Universe is composed by about e 120 causally disconnected regions. Within each region, the probability that the Higgs fluctuates to the desired tuned h in is about 10 −2−3 .
Underfluctuations produce a negligible PBH abundance. On small scales this is not a problem: only the average matters. Accretion of black holes after their formation would increase the average. On large scales of the order of the present horizon (those probed by observations that find a DM density more homogeneous than what could be justified anthropically), fluctuations of hðxÞ would produce a PBH density that is not homogeneous.
The effect of overfluctuations is much worse. As in a cosmic Russian roulette, a too large h in in one of the e 120 regions can form a vacuum decay bubble that, after the end of inflation, engulfs the whole Universe (a general relativity computation finds that innocuous bubbles that shrink and/ or expand behind a black hole horizon can form, but together with dangerous ones [8]). The probability of avoiding vacuum decay is roughly estimated as ℘ ∼ 2 −e 120 .
As we now discuss, this unlikely possibility cannot be justified on an anthropic basis.
To explain why, let us start from an analogous anthropic argument considered by Weinberg [40]: if observers exist only where the cosmological constant Λ is small enough to allow their development (at the price of a tuning with probability ℘ Λ ∼ 10 −120 , which is possible in multiverse with more than 10 120 different vacua), they should expect to see a cosmological constant around the anthropic bound. Once the desired Universe with small Λ is formed, it is relatively safe. Technically, the timescale for the variation of the ℘ Λ is the Hubble scale, presently 10 10 yr.
Within the "rescued Higgs fall" mechanism, the probability ℘ to form DM but no vacuum decay is much smaller. This might be by itself a problem, unless one relies on eternal inflation and argues that 1=℘ is smaller than infinity. The same Weinberg-like argument leads to expectations incompatible with experience. Indeed, the timescale for the growth of 1=℘ is very short. An observer that justifies its lucky survival to vacuum decay as needed for its existence should expect to be immediately executed by an expanding vacuum decay bubble, given that it arises with statistical certainty in the extra regions which continuously enter in causal contact with the observer due to the Universe expansion.
Needless to say, multiverse probability is a shaky concept, plagued by infinities. Nevertheless, the problem seems worrying enough that one wonders if it can be avoided or alleviated.
One possibility is devising a mechanism that suppresses vacuum decay by rescuing the Higgs more efficiently than the thermal barrier considered in [8,10]. For example, a nonthermal distribution 11 or an inflation that initially decays to the SM particles more coupled to h increases its thermal mass. The green dashed line in Fig. 11 shows the extra rescued region imposing what we believe is the most optimistic possibility: h end < h max end with V eff ðh max end Þ ¼ 0. Namely, we demand that the negative Higgs potential energy remains smaller than the inflaton potential V 0 ¼ 3M 2 Pl H 2 0 ; otherwise, nothing can stop the Higgs fall. Figure 11 shows that a more efficient rescue mechanism would not qualitatively change the picture. Alternatively, some mechanism beyond the SM could prepare the Higgs in the homogeneous state needed for the rescued Higgs fall mechanism. A study of this possibility goes beyond the scope declared in the title of this paper, where we wanted to see if some mechanism can generate dark matter within the Standard Model.

IV. CONCLUSIONS
We critically reexamined two different dark matter candidates that do not require new physics beyond the Standard Model.
The first is a hexaquark S ¼ uuddss dibaryon, which (being a spin 0 isospin singlet) might have a QCD binding energy large enough to make it lighter than M S < 1.876 GeV and thereby stable because all possible decay modes are kinematically closed. In Sec. II A we estimated its mass at the light of recently measured tetraquarks. We found that S could be light enough, possibly as light as 1.2 GeV, although we cannot provide a precise mass. In Sec. II B we computed its cosmological relic density, finding that it can reproduce the desired DM density if M S ≈ 1.2 GeV, while larger masses lead to smaller abundances. The dominant process that keeps S in thermal equilibrium is the scattering of two strange baryons, whose abundance gets Boltzmann suppressed at temperatures smaller than the strange quark mass, leading to S decoupling at the temperature T dec ≈ 25 MeV where S has the desired abundance. However, in Sec. II C we find, following the strategy of [2], that such a light S is excluded, because nucleons inside nuclei would bind in S faster than what is allowed by Super-Kamiokande bounds on the stability of oxygen. We reached this conclusion at the light of recent global fits of nuclear potentials used to compute the nuclear wave function of oxygen, which indicate that nucleons are close enough to make S production too fast. Both S and nuclei can be stable for M S ≈ 1.87 GeV; however, this mass leads to a relic S abundance much 11 We thank the authors of [8] for this suggestion. smaller than the DM abundance. Among sparse comments, we mention that direct detection of S on an antimatter target gives an annihilation signal.
In the second part of the paper, we considered the proposal of [8]: given that the SM Higgs potential can be unstable at large field values, during inflation the Higgs might fall from an assumed homogeneous vacuum expectation value h beyond the potential barrier toward the Planck scale. If the fall is tuned such that h almost reaches the maximal value that can be rescued by thermal effects, this process generates small-scale inhomogeneities that form primordial black holes. While DM as PBH seems excluded in the proposed mass range (just above the bound on Hawking radiation), the proposal is interesting enough to deserve further scrutiny. We confirmed the computations of [8] and extended them, adding a nonminimal coupling ξ H of the Higgs to gravity. We find that inflationary fluctuations can generate a quasihomogenous h only if ξ H is as small as allowed by quantum corrections that unavoidably generate it. Furthermore, we find in Sec. III F that the amount and mass of PBH depend on an extremely sensitive way on the uncertain SM and cosmological parameters, including two extra parameters introduced as assumptions: the Higgs starts falling from a homogeneous value at a given moment during inflation. The assumed homogeneity is, however, not the typical state present during inflation, where the Higgs has fluctuations of order Hubble. This is important for the present mechanism because it has a Russian roulette feature: the Universe is eaten by a vacuum decay bubble if the Higgs fluctuates to a value too high to be rescued by thermal effects in one of the e 120 causally disconnected patches in which the present horizon is divided while the Higgs starts falling. The probability of obtaining the observed Universe is so small, about ð1=2Þ e 120 , that trying to justify it through anthropic considerations leads to the issues discussed in Sec. III F.
We conclude that tentative searches and interpretations of dark matter as a phenomenon beyond the Standard Model remains a justified field.