Electroweak Symmetric Dark Matter Balls

In the simple Higgs-portal dark matter model with a conserved dark matter number, we show that there exists a non-topological soliton state of dark matter. This state has smaller energy per dark matter number than a free particle state and has its interior in the electroweak symmetric vacuum. It could be produced in the early universe from first-order electroweak phase transition and contribute most of dark matter. This electroweak symmetric dark matter ball is a novel macroscopic dark matter candidate with an energy density of the electroweak scale and a mass of 1 gram or above. Because of its electroweak-symmetric interior, the dark matter ball has a large geometric scattering cross section off a nucleon or a nucleus. Dark matter and neutrino experiments with a large-size detector like Xenon1T, BOREXINO and JUNO have great potential to discover electroweak symmetric dark matter balls. We also discuss the formation of bound states of a dark matter ball and ordinary matter.

Dark matter (DM) is one of the remaining mysteries in particle physics after the discovery of Higgs boson in 2012. After a few decades of searching for electroweak-sector-related dark matter particles with a mass around 100 GeV and with a null result [1], we have started to enlarge the scope of dark matter masses from both the theoretic model and the experimental search sides. For our visible sector, we have many interesting states of ordinary matter ranging from diluted gas to a dense neutron star. Analogously, it will not be surprising that there are many types of states for dark matter. Under certain circumstance, the majority of dark matter could be in a form of macroscopic state instead of free particle states. The well-known example is the primordial black hole dark matter [2], which has the Schwarzschild radius as its macroscopic size. Another established example the so-called "quark nugget" [3] with around nuclear energy density, which has the constituents of dark matter to be fermionic quarks and the geometrical size of 0.01-10 cm. In this paper, we will focus on another type of macroscopic dark matter with a bosonic constituent or "non-topological soliton" as named in the literature. For a scalar field with non-linear interactions, it has long been pointed out that there exists a spacially-localized state that can be a solution to the scalar classical equation [4]. The existence and properties of the non-topological soliton as a field-theory object have been studied extensively by T. D. Lee [5] and S. R. Coleman [6] and their collaborators (see Ref. [7] for a review), while its primordial production from early universe physics has also been worked out in Refs. [8][9][10]. In supersymmetrical models, Q-balls (the soliton states and named by Coleman), built of squarks and sleptons have also been proposed as a potential dark matter candidate [11][12][13]. With a conserved global internal symmetry, the non-topological soliton is an extended object with the lowest value of the energy for a fixed conserved charge, and therefore is stable at quantum level. The non-topological soliton is simply different from topological solitons, which has a quantized charge related to the algebraic geometry. For instance, a nucleon can be regarded as a topological soliton state of pions or Skyrmion because of π 3 [SU (2)] = Z [14].
After some preparation of soliton basics, we want to point out the main observation of this paper: in the simple Higgs-portal complex scalar dark matter model, a non-topological soliton state exists for dark matter and could be the lowest energy state per dark matter number. For such a simple dark matter model, the dark matter could be in the macroscopic soliton state with a very large dark matter number, which we will refer as dark matter balls (DMBs). One possible mechanism to produce dark matter soliton states from early-universe dynamics could come from the first-order phase transition of electroweak (EW) symmetry, which can be naturally realized based on the quantum-corrected Higgs potential from the complex dark matter particle loop. Below the EW phase transition temperature, the EW symmetry-breaking bubbles grow and push the dark matter number to be in front of the bubble wall. After a few bubbles meet each other and coalesce, the dark matter number is enclosed in a small region and is still in the hightemperature EW-unbroken phase. Based on our later estimation and assuming some initial dark matter-antimatter asymmetry, we have found that the dark matter number is mainly stored in the soliton or DMB state instead of free dark matter particle state.
One interesting feature of DMBs based on the Higgs-portal dark matter model is the interplay between the Higgs potential and dark matter field strength. For a positive quartic interaction of the two fields, a large complex scalar dark matter field inside the DMB provides an effective positive mass for the Higgs field and prefers the Higgs field to sit at zero or a negligible value. So, the dark matter state could be a Electroweak Symmetric DM Ball (EWS-DMB) in sense that the electroweak symmetry is unbroken inside DMBs. This particular feature of DMB also means that the interaction of DMB with ordinary matter is relative large. From our later detailed calculations, we have found that the scattering cross section of DMBs off a nucleon or nucleus saturates to the geometrical cross section, when a DMB has a large radius. Our observation could dramatically change the experimental search strategy for dark matter: instead of searching for the single-hit scattering event with a small recoil energy in a location deep underground, one could search for multi-hit scattering events at a location not necessarily underground. In consequence, neutrino-oriented experiments with a large size detector become suitable for this type of macroscopic dark matter. Or, searching for tracks in an ancient mineral like Mica may also discover this type of DMBs because of its very long, billion-year, exposure time. We will also discuss various search strategies for EWS-DMB.
The paper is organized as follows. We first work out the properties of soliton states with and without the dark matter bare mass and self-quartic interaction in Section 2. In Section 3, we study the early-universe productions of DMBs based on the first-order electroweak phase transition and obtain the characteristic charge, mass and radius for DMBs. We then calculate the scattering cross sections of DMBs with Standard Model (SM) particles in Section 4. The detection of DMBs in various experiments will be discussed in Section 5. We summarize our results in Section 6. Furthermore, we have also included four Appendixes: the calculation of the number of DMB nucleation sites in Appendix A, the calculation of the binding energy of bound states of EWS-DMBs and ordinary matter in Appendix B, the calculation of the bound states in a Higgs potential well in Appendix C and a simple example of scattering against a heavy object in Appendix D.

Soliton States in a Higgs-Portal Dark Matter Scenario
In the Higgs-portal dark matter scenario with a complex scalar particle Φ, 1 the most general renormalizable Lagrangian preserving a U (1) Φ symmetry is The U (1) Φ symmetry ensures that the elementary Φ quanta are stable, and therefore a DM candidate. This is one of the simplest extension of the SM to include dark matter. For reasons that will become clear in the following, we will focus on the region of parameter space with λ φh > 0 and m 2 φ,0 ≥ 0, so that the physical Φ mass squared is never negative, even in the absence of a vacuum expectation value (VEV) for H. We will also take λ φ > 0. 2 In this case, the global minimum of the tree-level potential breaks the EW symmetry spontaneously: H T = (0, v/ √ 2) with v = 246 GeV, and Φ = 0. The quartic coupling λ h is related to the Higgs boson mass m h ≈ 125 GeV [15] by λ h = (m h /v) 2 /2 ≈ 0.13. After electroweak symmetry breaking (EWSB), the free Φ particle mass is When the bare dark matter mass m φ,0 = 0, the Φ particle obtains all of its mass from EWSB and m φ = λ φh v/ √ 2. We are interested here in non-vacuum field configurations that are nevertheless stable due to the conservation of the charge associated with the global U (1) Φ symmetry. In the theory given by Eq. (1), the existence and properties of such solutions were worked out in [5] (assuming m φ,0 = 0 and λ φ = 0), thus providing an example of a "non-topological soliton" (for a review, see [7]). We will briefly review how these solutions arise and their salient features. We start with the case m φ,0 = 0 and λ φ = 0, to establish that such DM solitons exist even in this minimal case, which depends on a single free parameter, λ φh . This will also highlight the crucial role played by this coupling. In a second stage we will include the effects of the remaining two free parameters, m φ,0 and λ φ , which can affect the qualitative properties of the soliton solutions. We will describe the relevant features in Section 2.2.
The DM solitons are characterized by a non-vanishing charge which is obtained from the time-dependence Φ(x) = e −iωt φ( x)/ √ 2, with φ( x) real. We will focus on spherically symmetric solitons (that have the lowest energy) with φ( x) = φ(r) and H T = 0, h(r)/ √ 2 , obeying the classical equations of motion and subject to the boundary conditions In order to develop an intuition it is useful to write down an approximate description by neglecting the Higgs derivatives in Eq. (5). The motivation is that often the Higgs profile is nearly vanishing inside the DM soliton and takes the (almost) constant value v outside, approximating a step function. Thus, apart from the relatively small transition region, the neglect of the spatial derivatives can be justified a posteriori, thus permitting an effective description in terms of a single degree of freedom. 3 Eq. (5) then shows that one can have configurations obeying Inserting Eq. (6) into Eq. (4) one gets where the effective potential is obtained by using Eq. (6) in (minus) the potential terms of Eq. (1), but including the terms coming from the time derivatives: giving For later use, we reintroduced here the pure φ-dependent terms although for the time being we are setting them to zero. We then see that at large φ values, U eff increases quadratically with φ. Importantly, the origin is unstable provided where we again wrote the m 2 φ,0 dependence for later reference. As we will see, the small ω limit corresponds to large charge Q. Assuming Eq. (11), one can see that there is a solution that starting "at rest" (using an effective 1-particle mechanics in 1D language, with time evolution parametrized by r) at φ(r = 0) = φ 0 , rolls down the effective potential towards the hill at φ = 0, loosing in the process energy due to the effective friction term in Eq. (7). It is clear that by adjusting φ 0 , it is always possible to arrange for this motion to stop at φ(r = ∞) = 0. One can also see that since U eff (φ = 0) = 0, one must have U eff (φ 0 ) > 0. At the saturation point of Eq. (11), the term in braces in U eff (φ) evaluated at the matching point λ φh φ 2 = m 2 h takes It describes the 1-particle mechanics in 1D analogue, with a friction term, as given in Eq. (7). The particle starts at rest at φ = φ 0 (for r = 0) and comes to rest at φ = 0 (for r → ∞). For the ω = 200 GeV case one has φ 0 ≈ 179 GeV.
the positive value m 4 h /(16λ 2 h ). Thus, for ω 2 λ φh m 2 h /(4λ h ) it is possible to find solutions fully contained in the region λ φh φ 2 < m 2 h . Such solutions can have a non-negligible h inside the core of the DM soliton, and typically require a more careful analysis that takes into account the h derivatives that have been neglected in the effective description. However, when ω is very small, φ 0 must be such that λ φh φ 2 0 > m 2 h to satisfy U eff (φ 0 ) > 0. Translated into the behavior for h this corresponds to situations with (nearly) vanishing h inside the core of the DM soliton. We will therefore sometimes refer to such solutions as Electroweak Symmetric DM Balls (EWS-DMBs), or DMBs for short. The associated h-profile typically resembles the step-like profile captured by the effective description.
We show in Fig. 1 the effective potential, U eff , taking λ φh = 3 and m φ,0 = λ φ = 0, for several values of ω. The threshold value defined by the saturation of the inequality (11) is about 301 GeV, and we show an example in its vicinity. Together with the ω = 200 GeV case, it gives rise to DM solitons with a sizable Higgs VEV inside the core. (As we will see, for ω = 100 GeV one obtains solutions displaying a core with a small Higgs VEV, i.e. an EWS-DMB.) Finally, the ω = 400 GeV case does not satisfy Eq. (11) and leads to oscillating solutions that tend slowly to zero as r → ∞, and are therefore not localized. These do not belong in the class of solitonic solutions.
From the previous discussion, we also see that for DMB solutions one must have the scaling The size of the DMBs can be easily estimated as follows: setting h = 0 in Eq. (4), as is appropriate inside the soliton, leads to This function has an infinite number of zeros, each of which corresponds to a solution. We will focus on the solution associated with the first zero, which has the lowest energy. Near this first zero, h turns on leading quickly to the asymptotic value φ → 0 as r → ∞ (the excited solutions arise in a similar manner, but with additional nodes). The size of the transition, i.e. the thickness of the surface boundary separating the EW breaking and EW preserving phases is of order π/v. We therefore see that the size of the lightest DMB, denoted by R Φ , is about Inserting the approximate solution (13) in Eq. (3), together with Eq. (12), one finds As stated earlier, small ω maps into large Q. We also see that in this limit, we have the scaling With this qualitative understanding, let us now consider some examples of the full solutions to Eqs. (4) and (5).

Solutions to the Classical Equations of Motion
It is possible to obtain numerical solutions to the system (4) and (5) and the specified boundary conditions by the "shooting method". This depends on two variables: φ(0) = φ 0 and h(0) = h 0 . The first derivatives vanish, which provides the four initial conditions to uniquely specify the solution. In practice, one starts at a small r 0 to avoid the singular point at the origin. One can then adjust φ 0 and h 0 to obtain the solution that obeys φ(∞) = 0 and h(∞) = v. In practice one takes r = ∞ to mean an r max large enough that the neglected part can be seen to be numerically close to the desired solution. This procedure can be followed for any fixed set of Lagrangian parameters, and fixed ω. For a given model, one is interested in scanning over ω, i.e. in obtaining soliton solutions of different charge Q. We show in the left panel of Fig. 2 the φ (blue) and h (orange) profiles for three different charges, in the model defined by λ φh = 3 and m φ,0 = λ φ = 0. The choice of λ φh = 3 is motivated by the mechanism of formation of such DM solitons, to be described later, but the features discussed in this section are similar for any λ φh of order one. The Q = 20.3 case (flatter, dashed  The Q ≈ 20.3 case (dashed lines) corresponds to ω = 300 GeV, close to the threshold value ω th = m φ ≈ 301 GeV that allows such types of configurations, as determined from Eq. (11). This solution is quantum mechanically unstable against decay into Q free particle states. The Q = 21 case (ω = 280 GeV), although having a similar charge, is stable. The Q ≈ 505 case (solid lines) corresponds to ω = 110 GeV. The size of this DMB, as estimated from Eq. (14), is R Φ ≈ 0.03 GeV −1 , which is reasonable from the figure. Right panel: difference between the DM soliton mass, M Φ , and the energy of Q free Φ particles of mass m φ ≈ 301 GeV, for λ φh = 3, as a function of the charge Q (low Q region). The orange branch corresponds to soliton solutions that are unstable against decay into such non-bound free Q-particle states. The blue branch is stable. For the model shown, the boundary between the two branches is at Q S ≈ 19.5, corresponding to ω S ≈ 286 GeV. This is the smallest charge for a stable DM soliton.
profiles) corresponds to a choice of ω close to the threshold value ω th = m φ defined by the saturation of the inequality (11). One can see that it is very close to the vacuum solution. There is a second solution with a similar charge with Q = 21 that displays a better defined core. As we will explain next, the former solution is unstable against decay into Q free elementary Φ quanta, while the latter is a stable DM soliton. The third example has a larger charge Q ≈ 505, corresponding to a smaller ω = 110 GeV. It shows a clear core with a small Higgs VEV and a large value for φ. This DM soliton would fall in the category of EWS-DMBs defined above. Although the transition in the Higgs profile from zero to v is comparable to the core, one can see that the φ profile is reasonably well described by the approximate solution (13) (for r R Φ ). Indeed, for ω = 110 GeV, Eq. (14) gives R Φ ≈ 0.03 GeV −1 , in good agreement with what is seen in the figure. Obtaining full numerical solutions with larger cores is challenging, as the solutions are sensitive to an exponentially small h 0 . Such solutions can nevertheless be easily obtained in the framework of the effective description. We will also use the effective description to discuss the effects of the two additional parameters, m φ,0 and λ φ .
Before turning to the general case, we consider the mass of the DM soliton. In the mean field approximation we are using, this can be obtained by computing the classical energy of the configuration: where V H (h) is the SM Higgs potential and V Φ (φ) was defined in Eq. (10). Here, φ and h are derivatives with respect to r. Note that localized field configurations cannot have a welldefined energy: although the mean energy in the rest frame is given by M Φ , there are quantum mechanical fluctuations in the 3-momentum. Such effects can be taken into account by a proper separation between the collective center of mass coordinates and the vibrational modes. This can be achieved by defining appropriate coherent states followed by a projection onto zeromomentum eigenstates [16]. We will ignore such corrections and use M Φ above as a proxy for the soliton mass, since the above precision is sufficient for our purposes. We show in the right panel of Fig. 2 the mass of the lightest DM soliton as a function of Q. 4 One can distinguish two branches as one increases ω from small values up to ω th = m φ where the DM soliton solutions cease to exist. The charge decreases monotonically down to a minimum value (∼ 17.9 in the figure), then increases rapidly again and diverges as ω → m φ . It can be shown that all such DM soliton solutions are stable against small classical fluctuations [5]. It is, however, important to compare the DM soliton mass against the energy of Q free elementary Φ quanta. We plot the difference M Φ − Qm φ , which shows that there are two branches to be distinguished. The first branch (orange) is unstable against quantum mechanical decay into Q free particle states. The second branch (blue) is forbidden from decaying by a combination of energetic considerations and the conservation of the Q charge. In fact, they correspond to stable quantum mechanical states. This defines a Q S that separates the two types of solutions. For the model parameters used in the figure, one finds Q S ≈ 19.5 (corresponding to ω S = 286 GeV). Thus, the Q ≈ 20.3 profiles in the left panel of the figure correspond to an unstable soliton, while the Q ≈ 21 and Q ≈ 505 cases correspond to stable DM solitons.
In order to establish the scaling of M Φ with Q, let us focus on DMBs by assuming that the h field vanishes inside the core. Then φ is given by Eq. (13) for r < R Φ ≈ π/ω (and zero for r > R Φ ). Neglecting the surface tension contributions from h (i.e. setting h = 0 everywhere), one has from Eq. (17): where we exchanged φ 0 for Q using Eq. (15). We are interested in the lowest energy solution for fixed Q. Minimizing against the dynamical variable ω, we find that the minimum is at 4 We show here the low Q region of a scan over ω, obtained by solving the EOM numerically, as described above. However, for the (almost) horizontal (orange) part of the curve we use instead the first order approximation (for ω ≈ ω c ) derived in [5]: where the moment M 2 is calculable and gives M 2 ≈ 0.75. The reason is that in this region the excited states are split by small energy differences (that tend to zero as Q → ∞ on this branch), and it is difficult to isolate the ground state numerically. Thus, together with the results of the previous section, we have that for DMBs in the case that λ φ = 0, the following scaling laws between the charge, the size, and the mass of the DMB apply 5 These scaling laws hold for DM solitons with a large Q. In the low Q region displayed in the right panel of Fig. 2 somewhat different relations are obeyed, that can only be found by a more detailed numerical analysis. We focus on the stability point associated with the charge Q S that delimits the stable/unstable soliton configurations. In the left panel of This is a parametric relation across models as we vary λ φh . The parametric dependence for Q S (λ φh ) is shown in the right panel of Fig. 3, and can be reproduced by a broken power law in this range, as shown in the figure. From the information in both panels one can get M S Φ (λ φh ). One can similarly consider the radius for such a minimum charge DM soliton, which is well described by Thus, the typical radii for such charges are of order R ∼ few 10 −2 GeV −1 ∼ 10 −3 fm, while the masses are in the tens of TeV and above region. These scales arise from the weak scale due 5 We will see in the next section that m 2 φ,0 by itself does not change these scaling relations. 6 If one uses Eq. (19) to determine Q S , even though it is not meant to hold for the lowest charges corresponding to DM solitons that are not DMBs, one can estimate Q S ∼ 512π 4 m 2 h /(81λ 2 φh v 2 ). This gives a reasonable order of magnitude estimate for Q S , working better for larger λ φh . to charge enhancements, qualitatively similar to the scaling laws discussed above, but not as simple. Based on these plots we can estimate the energy density associated with the DM soliton configuration to be of order To the extent that the scaling laws given in Eq. (20) connect the low Q and high Q cases, we expect the same estimate to hold for very large Q DMBs.

Effects of the Dark Matter Bare Mass and Self-Quartic Interaction
We now comment on the effects of the remaining two parameters of the model, m 2 φ,0 and λ φ . Within the context of the effective potential description defined in Eq. (9), one can see that 1. The bare mass m 2 φ,0 can be easily included by defining an effective ω 2 ≡ ω 2 − m 2 φ,0 in the effective potential and associated EOM. One must only remember that when computing observables such as the mass of the DM soliton via Eq. (17), it is the orthogonal combination ω 2 + m 2 φ,0 that appears. Similarly, the charge Q is proportional to ω, and the combination ω enters only through φ. With the solutions for the m 2 φ,0 = 0 case at hand this can be easily taken into account.
2. The quartic self-interaction λ φ has a more significant effect: it changes the large φ behavior of the effective potential from the quadratic one used in the previous section, turning it down to reach an asymptotic behavior − 1 4 λ φ φ 4 (for λ φ > 0) (see the left panel of Fig. 4). This creates a maximum in the potential at some φ max . The soliton solutions must therefore satisfy φ 0 < φ max , since for φ 0 > φ max the solutions would run down the hill in the wrong direction and not be bounded. This is the scenario considered for Q-balls in Ref. [6], and we know that stable solitonic configurations exist in this case.
The first point could mean that even for ultraheavy elementary Φ particles that receive only a small contribution to their mass from EWSB, it could be possible to have solitonic configurations related to the weak scale, i.e. sustaining an EW symmetric "vacuum" in a finite region of space, inside the normal EW breaking vacuum.
Let us now describe some of the consequences of the quartic coupling λ φ , assuming for simplicity that we are interested in DM solitons with a large charge Q, such that they fall in the class of EWS-DMBs. In this case, the maximum of the effective potential described in point 2 above lies in the region λ φh φ 2 > m 2 h , where according to Eq. (9),  Fig. 1 for the effective potential as a function of φ, but including the Φ bare mass and quartic self-interactions. For λ φh = 3 and λ φ = 1, only the range 147.5 GeV < ω < 301 GeV allows for DM solitons. The particles marked as DMB (2) and DMB (4) are meant to illustrate the two classes of DMBs that can appear in this system (see footnote 7). Right panel: DM soliton profiles for several values of ω near ω c ≈ 147.5 GeV, taking λ φh = 3, λ φ = 1. The plateau approaches φ ≈ φ max ≈ 147.5 GeV as ω approaches ω c . The transition in the approximate Higgs profile, Eq. (6) is expected to also occur near r = R Φ , where R Φ is the size of the DM soliton, as read from the figure. This such that soliton solutions must obey ω > ω c . The conditon (11) must also be imposed, so that the origin be a maximum as opposed to a minimum, as discussed in the previous section. Thus, in the presence of λ φ , ω is bounded by non-zero values both from below and above. In the left panel of Fig. 4, we show the effective potential as a function of φ for several choices of ω and fixed λ φh = 3 and λ φ = 1. Only when ω ∈ (147.5, 301) GeV one finds trajectories where the effective particle, starting at an appropriate φ 0 , comes to rest at φ = 0 at r = ∞, with an exponential approach, so that it effectively reaches the second hill in a finite r. These are the finite energy, finite Q, DM solitons. We also indicate a categorization of two distinct classes of DM solitons in terms of the initial conditions in the particle mechanics analogue. The "quadratic DM solitons", discussed in the previous section, are denoted by DMB (2) in the figure.
Those for which the quartic Φ self-interaction plays an essential role, as we are discussing in this section, are denoted by DMB (4) , i.e. we will refer to them sometimes as "quartic DMBs". 7 Although the allowed range of ω is limited, soliton solutions with arbitrarily large charges exist. These occur for ω close to ω c , and are obtained by making the volume large, as they display a uniform charge density. Thus, such balls behave like aggregates of Φ matter [6]. We show in the right panel of Fig. 4 the numerical solutions obtained within the effective potential approach, for λ φh = 3, λ φ = 1, and for several values of ω near ω c ≈ 147.5 GeV, as obtained from Eq. (25). One can see that as ω approaches ω c , the size of the soliton increases, and the profiles resemble a step-function, much more than when the quartic coupling is absent or negligible. This can be easily understood from the particle mechanics analogy: one starts at rest near the top of the potential maximum at φ max , and slowly picks up speed for a long "time" r, generating a nearly constant inner core. At some point enough speed is attained and the particle falls down the potential in a short time, decelerates rapidly as it approaches the local maximum at the origin, and eventually comes to rest there. We can therefore use a simple step function profile for φ to estimate quantities of interest, where for concreteness we can take the size of the φ profile as the R Φ such that φ(R Φ ) = φ max /2. In the limit of ω − ω c ω c , the DMB radius shows a simple scaling as R Φ ∝ 1/(ω − ω c ), which can be seen in the right panel of Fig. 4. For λ φh = 3 and λ φ = 1, the overall coefficient can be fitted from the numerical results: The Higgs profile also takes a step-like form, with The Higgs profile "size" is determined by the R h such that φ(R h ) = m h / λ φh , which we have assumed is smaller than φ max . Taking the ratio of the two φ values that define these radii, we have where in the second relation we assumed ω ≈ ω c . For order one couplings λ φh and λ φ , the ratio is of order one, so that the two radii can be identified: When λ φ is small, there can be some difference between the two radii. However, one expects at least a 1-loop size of order λ 1−loop φ ∼ λ 2 φh /(16π 2 ), so that the above ratio of VEVs is not expected to be greater than (2π 2 ) 1/4 ∼ 2, and therefore the two radii should be close enough to be identified as in the case of order one couplings.
The charge of the DM soliton, Eq. (3), in the case ω ≈ ω c , is then approximately given by The soliton mass, Eq. (17), neglecting the h surface tension contributions, is given here by Compared to the energy, Q m φ , of Q free quanta, each of mass m φ as given in Eq. (2), we have We see that the above ratio is less than one when , this is always satisfied. The DM soliton is then the lowest energy per dark number state, and stable.
We also see that for large Q, large R Φ DMBs in the presence of a λ φ = 0, which have ω ≈ ω c , one has the following scaling laws between the DMB's charge, size and mass: Thus while in both types of DMBs we have discussed, M Φ scales with volume, they can carry very different charges. The formation mechanism for such aggregates of charges will determine the type of DMBs one would expect. We discuss these issues next.

Early Universe Production of DMBs
Depending on the early-universe history, there could be several possible ways to produce DMBs, in this section we concentrate a simple mechanism based on first-order phase transition. Especially, with the extension of the singlet scalar of the SM Higgs potential, the electroweak symmetry breaking is naturally a first-order one. Furthermore, the typical dark matter number in one DMB depends whether there is an asymmetry in dark matter and antimatter or not. For simplicity, we just assume that the dark matter asymmetry is given by some ultra-violent physics and have the production mechanism similar to the "quark nugget" in Ref. [3].

First-Order Electroweak Phase Transition
The tree level scalar potential is given by V ω=0 (h, φ) in Eq. (8), setting ω = 0. The form of this potential is equivalent to one obtained by addition to the SM Higgs potential of two real scalar singlets, corresponding to the real and imaginary parts of Φ. In the early universe, at very high temperature, the global minimum occurs at the electroweak symmetry preserving point ( h , φ ) = (0, 0). As the universe cools down, the global minimum happens at an EWSB vacua with ( h , φ ) = (v, 0). Depending on the coupling λ φh , one can have a "one-step" phase transition where the phase transition occurs purely along the Higgs direction.
To study the electroweak phase transition (EWPT), we consider the effective finite-temperature potential V eff (h, T ) where h is the real component of the SM Higgs doublet, H T = (0, h/ √ 2) and T is temperature [17][18][19][20][21][22][23][24][25] where Π i is thermal masses (or Debye masses) (see its formulas in [24] for instance). The first is the tree-level SM Higgs potential. The second term V CW is the one-loop contribution to the zero-temperature effective potential, also known as Coleman-Weinberg potential [26]. Using the on-shell renormalization scheme in the Landau gauge, it is given by [21] where g i is the degree of freedom for each particle, F i = 1(0) for fermions(bosons), m i (h) are masses in the presence of a background Higgs field with i = t, W, Z, h, Φ and ignoring lighter fermions. The finite-temperature correction term has where the integral with "−/+" sign denotes the thermal bosonic/fermionic function. Before we provide the parameter space for first-order phase transition, we want to note that requiring the ordinary electroweak vacuum with h = v = 246 GeV as the global vacuum at T = 0 or V eff (v, 0) < V eff (0, 0) sets a constraint on the coupling λ φh and the bare mass m φ,0 [27]. When m φ,0 = 0, this requires The two-loop effective potential could slightly change this numerical number. For the range of 0 ≤ m φ,0 ≤ 200 GeV, the upper bound on λ φh varies from 9.0 to 10.0. Therefore, in our numerical calculation for the parameter space of phase transition, we will restrict ourselves to this allowed range.
In the left panel of Fig. 5, we show the first-order phase transition temperature as a function of λ φh for different bare mass m φ,0 = {0, 100, 200} GeV. In each curve, we also separate it into two regions with the strong first-order phase transition region in blue with v(T c )/T c ≥ 0.6 [28] and the weak first-order phase transition region in red with v(T c )/T c < 0.6. For m φ,0 = 0 GeV, the first-order phase transition happens for λ φh 2, while the strong first-order phase transition happens for λ φh 2.6. As λ φh increases but below the upper bound in (36), the phase transition temperature decreases. For the benchmark point with λ φh = 3, the phase transition temperature has T c ≈ 134 GeV. In the right panel of Fig. 5, we show the ratio of the T c -dependent EWSB VEV v(T c ) over T c as a function of λ φh . Again, the strong(weak) first-order phase transition region is denoted in blue(red) color. As the coupling λ φh increases, the ratio of v(T c )/T c increases. We also note that for both plots, the Φ self-interaction quartic coupling λ φ does not play a role for the one-step phase transition evaluated at the one-loop level.

Formation of DMBs from First-Order Phase Transition
We discuss now how DMBs might be formed during the EWPT in the early universe. As we will see, the formation of the DMBs requires the transition to be a strong first-order. We will also discuss their expected average properties such as charge, mass and size.
For the purpose of this section, we assume that some high-scale physics, analogous to leptogenesis, has already generated a Q asymmetry, that we will call DM number, 8 with a yield Y Φ ≡ n Φ /s, which we treat as a UV-dependent free parameter. Here, the entropy density is s = (2π 2 /45)g * S T 3 , with g * S being the effective number of relativistic degrees of freedom. As a reference point, the SM baryon number asymmetry is measured to be Y B 10 −10 [29]. It would be interesting if there was a common origin for Y Φ and Y B , in which case one would expect Y Φ ∼ Y B , at least if the generation occurs at the same time. Realizing such a scenario would require additional model assumptions. However, one should note that the presence of the complex scalar can already lead a strong first-order EWPT, which is one of the conditions for EW baryogenesis. Thus, one may be able to build a model to also generate the DM number asymmetry within the framework of EW baryogenesis, which we will not explore in this paper.
We organize the analysis in three stages for conceptual clarity: 1. The "snowplow" stage, taking place around T c , when the EWSB nucleation process happens. We will argue that a large fraction of the DM number ends up in the unbroken phase, as opposed to the true vacuum (broken) phase.
2. The second stage is delimited by the formation of DMBs from the DM number stored inside regions of unbroken phase.
3. Subsequent to the DMB formation, the free DM number gets rid of its symmetric component, leaving behind the asymmetric yield Y Φ .
We will argue that up until the freeze-out temperature T F of the free Φ particles in the broken phase, DM number continues being accumulated inside the DMBs. The end result is that the amount of DM number stored in elementary Φ quanta is exponentially suppressed. We start with the snowplow stage. Just below the EWPT temperature T c ∼ 130 GeV, the EWSB (true vacuum) bubbles start to pop up, and grow when they surpass a critical size. During the bubble nucleation process, one immediate question is whether the DM number stays mainly in the unbroken or broken phases. To address this, we first give a simple kinetic argument, assuming that m 2 φ,0 = 0 (or that it can be neglected). 9 At leading order, the answer involves the Φ particle mass m φ (T ) ≈ λ φh /2 v(T ), the phase transition temperature, T c , and the bubble wall speed β w (or the corresponding boost factor γ w = 1/ 1 − β 2 w ). It is convenient to work in the bubble wall's rest frame, which sees a stream of Φ particles moving in theẑ direction (this is just the direction of expansion of the bubble wall in the plasma frame). From energy conservation, the condition for a Φ particle to remain in the unbroken phase, where it is massless, Here the hat denotes thatp z is the momentum of the particle in the wall's rest frame. Boosting this condition back to the plasma frame, we arrive at For a non-relativistic wall speed, β w 1, this condition simplifies to p 2 z ≤ m 2 φ,c . Using the Bose-Einstein statistics distribution, the average momentum is p 2 For large m 2 φ,0 , such that the Φ particles are non-relativistic already at T c in the unbroken phase, taking into account the conserved DM number is more involved. Considerations analogous to the ones detailed in this section would allow to determine how much of the DM number ends up in the broken versus unbroken phases. However, this would be relevant only in the presence of additional physics that would account for the first-order phase transition, since the small abundance of Φ particles would have a negligible effect on the finite-temperature Higgs effective potential. Thus, we do not consider this case, and focus on 0 ≤ m 2 φ,0 ≤ (200 GeV) 2 , as discussed in Section 3.1. Note however, that Eqs. (37) and (38) remain unchanged in the presence of an arbitrary m 2 φ,0 .
3.5 T 2 c . So for the bubble wall to "snowplow" the DM number into the unbroken phase one needs From the relation between λ φh and v(T c )/T c shown in the right panel of Fig. 5, one can infer that one needs a modestly large value of λ φh 4 so that most of the DM number stays in the unbroken phase. Instead of kinematic arguments, one can also provide an estimation based on chemical equilibrium considerations. Here the condition of chemical equilibrium, µ Φ , allows to estimate the ratio of DM number in the high-temperature, "h", and low-temperature, "l", phases. In both phases and not far below T c , one has µ Φ /T 1 (small asymmetry, see footnote 9). For a relativistic gas of elementary Φ particles, one has at T where Thus, when in chemical equilibrium, For a heavy elementary Φ particle, r is suppressed. In the case that m 2 φ,0 = 0, and assuming that the inequality (38) is saturated, one has m φ,c /T c ≈ 1.86, and r ≈ 0.15. However, the chemical equilibrium between inside and outside of the DMB could be kept until a lower temperature, T F . This is because the free Φ and Φ † can be absorbed by the DMBs or a large binding energy can be released when free Φ and Φ † particles enter DMBs.
The relevant process is Φ Q + Φ ↔ Φ Q+1 + X with X denoting SM particles. We first note that when the temperature is above the "binding energy", E bind (T ) ≡ m φ (T ) − ω(T ) [with ω(T ) as the temperature-dependent energy per charge for the soliton state], both forward and backward processes are efficient. The chemical equilibrium between DMB and free Φ state is reached. As T < E bind (T ), the free Φ can be absorbed by the DMBs, but not the other way. The freeze-out temperature, T F , for Φ Q + Φ → Φ Q+1 + X, is anticipated to be satisfy T F < E bind (T F ). The free Φ particle absorbing rate by DMBs is estimated to be with the radius of DMB as a function of T . Just below the temperature of the ending of nucleations T f , the number of DMB, N Hubble DMB , within one Hubble patch is estimated in Eq. (92). Using it, the averaged radius of DMB is around R Φ (T f ) ∼ (N Hubble DMB ) 1/3 d H . The inverse Hubble distance is 1/d H = H(T f ) = π 2 g * /90 T 2 f /M P , where the reduced Planck mass M P = 2.43 × 10 18 GeV. As the Universe cools, the radius of DMB also reduces, which requires a more detailed understanding of how DMB evolves at a non-zero temperature and non-zero vacuum pressure. As a simplistic estimation, we assume its radius shrinking velocity is β w from T f to the freeze-out temperature T F . Using the relation t = 1/(2H), we have the radius as a function of temperature as Substituting R Φ (T ) into Eq. (42) and requiring Γ Q+Φ→Q+1 H(T ), we have the approximate freeze-out temperature as For sure, our estimation of the freeze-out temperature and the ratio r is a naive one, but the results do suggest that the fraction of dark number in the free Φ particle state can be neglected for the phenomenological purpose.
As the universe cools down, Φ and Φ † states annihilate into SM particles and leave behind the asymmetric component. If one is allowed to ignore additional DM number shuffling processes from one phase to the other or below the freeze-out temperature T f , DM in our Universe could be composed of both macroscopic DMBs and microscopic Φ-particle states, if both states are stable. Focusing on quartic DMBs, we saw in Eq. (29) and the subsequent discussion, that the energy per charge of the DMB is given by ω c , as given in Eq. (25), which is less than the free particle mass, ω c < m φ , for the generic parameter space defined by λ φ /λ φh < 1.4 [see Eq. (31)]. Given an asymmetric yield Y Φ , one can calculate the ratio of the DM and ordinary matter energy densities as To fit the measured value of Ω DM /Ω B ≈ 5.4 [29], one needs ω c Y Φ ≈ 5.4 × 10 −10 GeV. If the yield Y Φ is comparable to the ordinary baryon one, the model parameters are then required to satisfy ω c ∼ 5.4 GeV, which is the well-known situation of asymmetric DM models. We note that, to achieve such a small value of ω c , we need λ φ ≈ 1.8 × 10 −6 , which is much smaller than its natural lower-limit value of O(λ 2 φh /16π 2 ) for λ φh = O(1). We therefore take λ φ ∼ 10 −2 and choose ω c ∼ 50 GeV as a benchmark model point, for which one has Y Φ ∼ 10 −11 .
Having discussed the DM abundance and its rough composition in terms of free elementary Φ particles versus such trapped inside DMBs, we can now estimate the average DM number in a DMB. Since, as we have argued, we expect most of the DM number to stay in the false EWS vacuum, we will neglect the small contribution inside the EWSB bubbles in the following estimates. Given the DM number density around T c , the average DM number inside a DMB can be estimated as the ratio of the total number within one Hubble patch over the number of DMBs in one Hubble patch. The total DM number within one Hubble patch has Here, we also took g * S ≈ g * ≈ 108.75. The number of DMBs in one Hubble patch is of the same order of magnitude as the number of EWSB nucleation sites, which is sensitive to the detailed properties of the EW bubbles or the model parameter λ φh . In Appendix A, we have estimated the number of nucleation sites in terms of the model parameter λ φh [see Eq. (92)]. After some numerical fit, the number of DMBs within one Hubble patch reads which captures the dominant dependence on λ φh (see also Table 1 for numbers for λ φh from 3 to 7). When λ φh varies from 3 to 7, we find that N Hubble DMB decreases from 1.1 × 10 13 to 4.0 × 10 7 . Finally, the average DM number in one DMB is estimated to be In the range λ φh ∈ [2,9], the average DMB mass ranges from 1.1×10 24 GeV to 9.2×10 33 GeV or from 1.9 g to 1.6 × 10 10 g. In our subsequent phenomenological considerations we will allow for a wider range of DMB masses, and use the first-order phase transition values as guidance. From Eq. (28), which applies to quartic DMBs, we see that the DMB radius scales like ω −1 c λ 1/3 φ Q 1/3 , in the limit that m 2 φ,0 = 0 (i.e. ω c = ω c ). Using also Eq. (25), we therefore have For the range of λ φh ∈ [2,9], the DMB radius varies from 8.1 × 10 4 GeV −1 to 1.7 × 10 8 GeV −1 or from 1.6 × 10 −9 cm to 3.3 × 10 −6 cm.

Scattering of DMBs with SM Particles
In this section we discuss a number of issues related to the scattering of DMBs. We start with a discussion of bound states of SM particles inside the DMB, as this bears on its scattering cross section, and other possible effects.

Bound States
We have seen how DMBs sustain a core where the EW symmetry is essentially unbroken, while outside the soliton the usual EWSB vacuum is quickly reached. Such Higgs profiles display a sharp transition of size ∼ 1/v, separating a much wider EW preserving region from the symmtrybreaking vacuum outside. It acts as a potential well, seen by any SM particle or bound state with a strength dictated by its coupling strength to the Higgs boson. Typically, the Higgs well is invisible only to massless particles like photons (at tree level), as well as neutrinos due to their lightness. We can model this physics as a 3D potential well (i.e. with a sharp transition), and use the intuition from the quantum mechanical treatment of such a problem. However, since elementary particles trapped inside the DMB see essentially no Higgs VEV, they behave like trapped massless states. We will therefore include the kinematic relativistic effects. For the SM fermions and gauge bosons, spin effects are also expected to be important in determining the spectrum of bound states. Although our methods can be generalized in a straightforward manner to include such effects, we will neglect them for simplicity, and aim at getting only a qualitative understanding. Hence, we will be thinking of appropriate scalar particles as proxies for the SM fermions and gauge bosons. Similarly, we will neglect corrections from the possible creation of particle-antiparticle pairs, which would require a significantly more complex quantum field theory treatment. In summary, we are interested here in the Klein-Gordon equation in the presence of a 3D well, in its one-particle interpretation. We will also be interested in particles, such as hadrons or nuclei, that get most of their mass from sources (QCD dynamics) other than EWSB. Such bound states can be described by a non-relativistic analysis, which can be obtained by taking the non-relativistic limit of the case above.
We describe the appropriate treatment in Appendix C, and invoke here only the main results, which are sufficiently intuitive. In Appendix B we discuss the backreaction of such bound states on the DMB, which is expected to be small due to the large dark matter or Φ number composing a DM soliton with a large Q. This is in spite of the large number of particles that can get bound by the DMB, as we shall discuss.
The structure of the Klein-Gordon equation is identical to the Schrödinger equation in a 3D potential well, with the replacement E → E 2 /2m χ (where m χ is the mass of the particle in the normal EWSB vacuum), together with an appropriate mapping of the potential (by a constant rescaling). Thus, the solutions are given by spherical Bessel functions inside and outside the DMB, whose matching at the boundary lead to the condition for the spectrum. For example, for particles that get their mass solely from EWSB, one gets for the s-wave states for a DMB with a radius R which determines the bound state spectrum E n . The threshold radius to have a bound state can be estimated to be Numerically, one has R t th ≈ 0.009 GeV −1 , R e th ≈ 3142 GeV −1 , R ν th 1.6 × 10 10 GeV −1 for m ν < 0.1 eV. So, for the DMB radius from first-order phase transition in Eq. (50), neutrinos can not be bounded, but other fermions do have bound states.
For particles that get only part of their mass from EWSB, the s-wave spectrum is determined by where we denote the mass in the EWSB vacuum by m N (as for nucleons), and the coupling of the N particle to a single Higgs boson by y hN N (e.g. the nucleon Yukawa coupling). Note that Eq. (51) can be obtained from Eq. (53) by setting m N = y hN N v. Also, if E ≡ m N + ∆E with |∆E| m N , one can check that Eq. (53) reduces to the non-relativistic result determining the binding energies |∆E n | in the presence of a potential well of depth V 0 = (y hN N v) 2 /2m N . The threshold radius to have a bound sate for nucleons is or numerically R N th ≈ 5.8 GeV −1 . The corresponding threshold radius for a nucleus is even smaller by a factor of 1/A. As a first application, consider Eq. (51) in the case of a large DMB with m χ R 1. The eigenvalues are then close to the poles of the cotangent. In particular, keeping the first order correction, the lowest energy solution is at As shown in Appendix B, when backreaction effects can be neglected, the mass of the DMB with a bound state χ as above is given simply by M is the DMB mass in the absence of χ, and E χ n is the bound spectrum as described above. So, for the particles satisfying m χ R 1, the bound state masses for different species are still ordered in their SM relations or E t 0 > E b 0 > E c 0 for instance. On the other hand, the mass spectrum is very degenerate, which may kinematically forbid some SM decaying channels. The energy balance for a process A → B + C + · · · can be described simply in terms of the bound state spectrum. For example, for the SM process t → W + b, assuming all the particles are bound to the DMB in the ground state, and using Eq. (55), one can easily see that this decay channel is kinematically forbidden in the limit of m b R 1. However, the top quark in DMB could still decay into light fermions via three-body processes like t → e + ν e d for m e R, m d R 1 and no bound states for electron, neutrinos and down quark.
Consider now a nucleon. Here, for y hN N vR 1, the binding energy for the lowest energy state is determined by the vanishing of the denominator on the r.h.s. of Eq. (53), E bind N,0 ≈ m 2 N − y 2 hN N v 2 − m N ≈ −0.04m N and around 38 MeV numerically. 11 Hence the nuclei can be safely treated in the nonrelativistic limit. Compared to the typical electron binding energies which are of order m e (at least if R is large enough for the potential well to sustain several bound states), the typical nucleon binding energies are larger, but only by a factor of order 0.04m N /m e ≈ 80. Let us compare the proton versus neutron cases. Using m p = 938.27 MeV and m n = 939.56 MeV [30], together with y hpp = 1.12 × 10 −3 and y hnn = 1.14 × 10 −3 [31], we get Thus, the neutron is slightly more deeply bound and "lighter" than the proton inside a DMB. This reflects the fact that inside the DMB the quarks are massless, and the fact that m d > m u in the normal vacuum gives a positive contribution to m n − m p (that dominates over the negative electromagnetic contribution) is absent inside the soliton. Hence, the proton is "heavier" than the neutron when bound to a DMB. This means, in particular, that for such bound states, we can have the process with the neutrino escaping the DMB and the positron bounded inside DMB. 12 Thus, trapped protons tend to decay into neutrons, in analogy to a neutron star. One can reason along similar lines to address other processes involving particles bound to a DMB. One should keep in mind that, in general, a DMB can be expected to be a complicated object carrying with it a cloud of different types of particles. Finally, we also note that there are many bound states in a DMB with a large R. The maximum number of angular momentum for the bound states has l max ∼ R/R th , i.e. roughly when the order of the spherical Bessel function passes the number of cycles that the l = 0 Bessel function can fit in the size R. Thus, we can put an upper bound on the number of orbital states as (R/R th ) 3 or more explicitly (58) 11 Although the expression for E bind N,0 we have quoted is only approximate, one can check that it reproduces the correct difference in binding energies between the proton and neutron bound states by solving Eq. (53) for the ground states numerically. 12 We are assuming here that QCD is in its standard chirally broken phase, an issue that requires analysis. Notice also that in this case the W gauge boson receives a contribution to its mass, unrelated to the Higgs, which for simplicity we have ignored in Eq. (56).
One should keep in mind that the above estimates do not take into account the strong and electromagnetic interactions between these particles. Using a benchmark R = 2 × 10 4 GeV −1 , the number of states for various states have N t orbital ∼ 10 18 , N N orbital ∼ 10 10 and N e orbital ∼ 10 2 . If all the nucleon states are occupied, we estimate a nucleon density of about 10 10 /(4πR 3 /3) ∼ (100 MeV) 3 , which corresponds to an inter-nucleon distance of about 2 fm, a density of the order of the nuclear density [as seen from Eq. (58), this is independent of R]. The total contribution to the DMB mass due to these nucleons is about 10 10 m N ≈ 10 10 GeV, where the small binding energy per nucleon is neglected. The captured nuclear matter gives a negligible contribution to the DMB mass, which is dominated by the φ and h contributions.

Scattering Off a Nucleon or Nucleus
After the previous cursory description of bound states of SM particles in a DMB, we turn to the question of DMB scattering from normal matter such as nucleons or nuclei.
The first observation is that DMBs are expected to be heavy compared to the target particles. For instance, Figs. 2 and 3 show examples with masses in the multi-TeV range and above, 13 and we will actually discuss much heavier objects like in (49), as expected from our discussion on how they can be produced. This means that one can analyze the scattering problem treating the DMB as an infinitely heavy object generating a fixed Higgs background field, while treating the nucleon or nucleus as a light particle scattering against such a fixed potential-a quantum mechanics scattering problem. We describe the procedure in Appendix D, in the context of a simple toy model. Since the typical velocities involved if the DMBs are gravitationally bound to our galaxy are of order 10 −3 c, one is safely in the non-relativistic regime. As discussed previously, we can further model the Higgs background as a 3D potential well.
As is familiar from the non-relativistic quantum mechanical treatment of scattering processes, the scattering cross section can be affected significantly when bound states are available. In the present context, the presence of bound states depends on the size of the DMB, with a threshold radius R th ∼ (y hN N v) −1 that depends on the scattering particle (normal vacuum) mass m N (as discussed in Appendix C). When the DMB is small enough so as to not allow for bound states, one can treat the problem in the Born approximation (and safely in the q = 0 limit). For DMBs with a large charge, on the other hand, they are large enough to contain a very large number of bound states. In such a situation a partial wave analysis is more appropriate. We will describe the partial wave analysis in some generality first, as the Born approximation case can be understood as an appropriate limit of the partial wave result.

Partial Wave Analysis
As explained above and in Appendix C our problem is mapped into the corresponding nonrelativistic problem by a simple reinterpretation of certain quantities. In particular, for a partial wave l in the 3D potential well we are using, the spectrum is determined by matching the logarithmic derivatives at r = R, and takes exactly the form of the non-relativistic result: K j l (KR) j l (KR) = k j l (kR) cos δ l + n l (kR) sin δ l j l (kR) cos δ l + n l (kR) sin δ l , where j l and n l are spherical Bessel functions of the first and second kind respectively, and δ l is the scattering phase shift for the l-th partial wave. The only difference is in the interpretation of the wavevectors k (outer region) and K (inner region). For particles that get their mass only from EWSB, such as an electron, they are given by where m is the mass of the particle in the normal EWSB vacuum. For particles like nucleons that get contributions to their mass from sources other than EWSB, they are given by where m N is the mass in the normal vacuum, and y hN N is their coupling to a single Higgs (e.g. a Yukawa interaction). The cross section is given by For given parameters, it is straightforward to obtain numerically sin δ l from Eq. (59), and do it up to large enough l that the sum in Eq. (64) is observed to converge.
In the left panel of Fig. 6 we show the result for scattering against a nucleon, as a function of the DMB radius R. We plot the cross section, as computed from Eqs. (62), (63) and (64), summing up to l = 100, where we have checked that the sum has been saturated in the range of R shown. We use m N = 938.9 MeV and y hN N = 1.1 × 10 −3 , and assume a typical DMB scattering momentum of k = 10 −3 m N . We see that, as a function of R, the cross section displays a complicated resonant structure. For small DMBs with R < 1 GeV −1 , as those shown in Figs. 2 and 4, the cross section is suppressed. This is the regime where it can be computed in the Born approximation, as will be discussed in the next subsection. Well above the most prominent resonances, it displays a "hard ball" behavior that drops from σ = 4πR 2 to σ = 2πR 2 , although with some additional small scale resonant structures. These are large DMBs that are of phenomenological interest given our previous considerations on their production.
In the right panel of Fig. 6 and summing up to l = 100(900) for nucleon(oxygen), we show the cross section for DMB scattering against a nucleon or an oxygen nucleus, 14 for a benchmark 14 Here we use m O ≈ A m N and y hOO ≈ A y hN N , with the atomic number A = 16. DMB of radius R = 2×10 4 GeV −1 ≈ 0.04Å, which corresponds to a geometrical cross section of πR 2 ≈ 4.9×10 −19 cm 2 . One can see that the cross sections follow roughly a "hard ball" behavior. For a smaller radius (for instance R = 6000 GeV −1 ), there exists also a superimposed resonant structure. We note, however, that the detailed structure is rather sensitive to the coupling of nucleons to the Higgs. One should keep this in mind, as we have made approximations and neglected physical effects such as spin or other relativistic effects that can affect the spectrum of bound states. Nevertheless, we learn that for such large DMBs, the cross section is expected to be between 2 and 4 times the geometric cross section. This holds, in particular, whether the scattering is off a nucleon or a much heavier nucleus. Our considerations in this section depend only on the DMB radius. However, the charge and mass of the DMB for this radius depend on the type of DMB. For instance, for DMBs sustained with λ φ = 0 (quadratic DMBs), which obey the scaling laws of Eq. (20), one has for R = 2 × 10 4 GeV −1 For a DMB sustained by the quartic self interaction λ φ = 10 −2 [see paragraph after Eq. (45)], one finds from Eqs. (25) and (28), We see that their masses are comparable, as expected from the fact that in both cases the mass scales with volume, and the underlying scales involved are just the weak scale times functions of couplings taken to be of order one. Only the associated charges are significantly different. From our discussion in Section 2.2 we expect that the second case is more likely (i.e. the quartic instead of the quadratic DMBs).

Born Limit
Let us now turn to the Born limit, applicable for small enough DMBs that do not sustain any bound states (for the given target particle). We will see that in this case one can compute in more detail the scattering cross section, without assuming a 3D well model. It also addresses directly the constraints on DMBs with masses not exceedingly above the weak scale.
In the presence of the nontrivial Higgs background induced by the DMB, a particle approaching this region sees a potential where we choose that V (r) → 0 as r → ∞, and y is the coupling of the scattering particle to the Higgs (e.g. y hN N for nucleons). In the first Born approximation the scattering amplitude is given by where m is the reduced mass of the system and q = | q| is the momentum transfer. For the case at hand, the range of the Higgs potential well is much shorter than the length scales that can be probed by the typical q ∼ 10 −3 m. We can therefore set q = 0, and compute Eq. (68) numerically for the Higgs profiles found as described in Section 2.
We have computed the DMB-nucleon scattering cross sections for a number of models with different λ φh and different DMB charges (see Section 2.1). Depending on the value of λ φh , these span DMB masses from about 5 − 500 TeV. 15 We find that the cross sections for this range of masses are well described by the relation where we have considered all the models together as they all fall reasonably close to the above parametrization. Comparing to the Xenon1T bound [1], which in this region can be parametrized as we see i) that at the lowest masses of order several TeV, σ ΦN is already excluded by direct detection searches, and ii) that σ ΦN increases with mass faster than the linear growth in Eq. (70). We note that the models considered above all have a vanishing bare mass m φ,0 = 0 (and λ φ = 0). Increasing the value of m φ,0 does not qualitatively relax the constraint. For quartic DMBs, where λ φ is important, the DMBs have a large radius and a large geometric nucleon scattering cross section, much larger than those captured by the Born approximation above.
Hence, all the models at low masses are excluded by direct detection constraints up to the experimental reach (about 2.8 × 10 18 GeV for Xenon1T and 1.4 × 10 21 GeV for BOREXINO), and we do not consider them in further detail. Our focus is instead in the large DMB mass limit.

DMB Detection
Before we discuss the detection of DMB, we want to briefly discuss the collider search for Higgsportal dark matter. For the dark matter particle mass is below one half of the Higgs boson, or m φ < m h /2, the Higgs boson can decay into two dark matter particles and has additional invisible decay branching ratio [32]. For the DMB formation scenario from first-order phase transition, the coupling λ φh is needed to be above around 2, such that the dark matter particle mass m φ is heavier than the Higgs boson mass. For this case, the collider constraints from the LHC are dramatically reduced because of the three-body production phase space for producing off-shell Higgs-mediated two dark matter particles and one additional jet. For instance, the parton-level production cross section for λ φh = 2 or m φ = 246 GeV with m φ,0 = 0 is around 0.1 fb with a missing transverse energy cut above 250 GeV at the 13 TeV LHC. With 36.1 fb −1 , the number of signal events is three orders of magnitude smaller than the uncertainty of measurement [33]. So, there is no collider constraint on the model parameter space with λ φh 2 considered in this paper.
Another constraints on the model parameter space come from direct detection of the free Φ particle. As discussed around Eq. (44), the free dark matter particle is subdominant of the total dark matter energy density. However, given the stringent direct detection constraints on a dark matter with a mass around 100 GeV, a non-trivial constraint on the coupling λ φh may apply here. Using the coupling λ φh v h ΦΦ † and y hN N h N N with N = n, p and y hN N ≈ 1.1 × 10 −3 , the spin-independent scattering cross section has the formula of where we have used the reduced mass to be µ Φ−N ≈ m p in the limit of m p m φ and the small bare mass limit with m φ,0 m φ . The latest constraints from Xenon1T have set an upper bound for dark matter scattering cross sections in Eq. (70), which can be translated into a constraint on our model parameter space as which is well satisfied for a freeze-out temperature below 1 GeV in Eq. (41).

Multiple Scattering Signals for a DMB with a Large Q
For DMBs with a small Q and as we discussed around Eq. (70), the scattering cross section is small for DMB to reach the underground detectors, but not small enough to satisfy the direct detection constraints. On the other hand, the DMB could have a very heavy mass and a large scattering cross section (see Fig. 6). Although we did discuss the potential early universe productions of DMB from a first-order phase transition and obtained some benchmark radii for DMB in Eq. (50), we will keep the DMB radius and mass as free parameters to discuss its detection potential. As a simple recap what we have learned for the properties of DMB in Section 2: for m φ,0 = 0 with ω c = ω c , the mass DMB has M Φ = Q ω c , while the radius has R Φ = (3λ φ /4π) 1/3 ω −1 c Q 1/3 . For a large value of R Φ , the scattering cross section of DMB with a nucleus reaches a geometrical one with (see Fig. 6) which will be used as a prediction based on DMB properties. We emphasize here that once the radius of DMB is large enough, the scattering cross section with a nucleus is approximately independent of the nucleus mass number. For the heavy DMB with a large scattering cross section, the ordinary underground dark matter direct detection experiments become less sensitive. The ideal experiment would be the one with a large product of the exposure time and the effective detector area. For a long exposure time, one could consider Mica type experiment [34] that looks for tracks generated by DMBs passing by. For a large effective detector area, one could adopt the neutrino-oriented detectors with a small enough energy-trigger threshold. In the following, we discuss the search potential for a few experiments.

Mica Constraints
For our DMB with or without quark or neutron matter inside, the scattering cross section is in the range of 2-4 times the geometric cross section (see Fig. 6). Since the scattering off a nucleon and a nucleus have similar cross sections, we can ignore the detailed chemical components of Mica for constraining the scattering cross section or the geometric size of DMB. To have a dark mattergenerated track in Mica, one at least requires one encounter event or ρ DM /m DM v A det t exp ∼ 1. Based on Ref. [34], the experiment has A det ≈ 595 cm 2 and t exp ≈ 0.6 × 10 9 yr. For the averaged dark matter velocity with v ∼ 10 −3 c, one has So, DMB with a mass 1.0 × 10 26 GeV may leave a track in Mica. Following the treatment in the paper by De Rújula and Glashow [35], the energy loss rate is Here, ρ A i is the individual element energy density and ρ mica ≈ 2.88 g/cm 3 . The velocity v should be the speed at a depth of around L ∼ 3 km underground. 16 For the velocity, we simple take it to be the averaged dark matter velocity around our solar system. To leave a track in the Mica experiment [34], the Mica stopping power has to be above So, the constraint on the DMB scattering cross section has

Xenon-1T Sensitivity
For the Xenon1T experiment, one could search for multi-hit signals as discussed in Ref. [36].
Since the liquid Xenon TPC has 1 m in diameter and 1 m in height, we simply take the detector area as A det ≈ 10 4 cm 2 . Taking the observation time of one year or t exp = 1 yr, the upper limit on dark matter mass is 2.8 × 10 18 GeV. Different from the situation in a neutrino detector, one could count the numbers of hits from both the scintillation and electroluminescence of electrons that have drifted into the gas above the target liquid. The number of hits is estimated to be N hit ∼ σ ΦA n det L det . The liquid Xenon has a density of 3.1 g/ml with the main abundant atomic mass number of around A = 131. Taking L det ≈ 1 m and N hit = 5, the projected probing cross section has σ ΦA 3.5 × 10 −24 cm 2 . (78)

BOREXINO Sensitivity
For the neutrino experiment, BOREXINO, located in the Gran Sasso underground laboratory with the average rock cover of about 1.4 km and based on organic scintillator with a low energy trigger threshold [37], it can also constrain some of the DMB parameter space. Using Eq. (74) and A det = 5 × 10 5 cm 2 and t exp = 10 yr, the dark matter mass is required to be below around 1.4 × 10 21 GeV to have a few encounter events. Following the experimental paper [38], the trigger of BOREXINO requires at least 25 to 30 hits of PMTs or 50 to 60 keV energy threshold during a selected time window around t select ≈ 100 ns. However, there is some efficiency factor κ for the nuclear recoil energy converting to detectable photons. Following the experimental measurement [39] and the semi-empirical calculation [40], the factor κ ≈ 10% for the recoil energy below around 20 keV. In the experimental paper [39], the possible threshold energy around 2.8 keV for Carbon is mentioned, but 16 One may worry about the overburden effects on reducing the DMB velocity when it reaches Mica. The relative change of the velocities from ground with v i to underground with L and v f has (v f − v i )/v i ≈ σ ΦA ρ ⊕ L/M Φ , which is around 10 −11 and negligible for L = 3 km, ρ ⊕ = 2.9 g/cm 3 (for crust) and the benchmark model point in Eq. (73). Beyond the DMB model, we will consider the parameter region with σ 2 × 10 −14 cm 2 × (m DM /10 16 GeV) in order to ignore the overburden effects.
was not crystal clear. The BOREXINO uses the organic scintillator pseudocumene (C 9 H 12 ) with a density of 0.88 g/cm 3 . Since Carbon occupies around 90% of the total density and has a larger recoil energy, we mainly keep Carbon when we derive a constraint on scattering cross sections.
Since the kinematics has M Φ m C , the reduced mass of the two-body scattering system has m r ≈ m C . In terms of the scattering angle θ * defined in the center-of-mass frame, the energy deposited in the detector is For a simple Maxwellian halo and ignoring the motion of the Sun and Earth and the escaping velocity, the differential rate per recoil energy has a simple dependence on E R via [41] T with v 0 = 220 km/s. The averaged recoil energy can be estimated to be For one incident DMB hitting the detector, the interaction rate is Γ = n C σ ΦA v rel . The constraints on the cross section after satisfying the trigger requirement is For the multi-hit signal events from DMB passing by, one could use the event shape to isolate them from backgrounds. Before dedicated experimental searches, we take the above equation as the potential reach from BOREXINO.

JUNO Sensitivity
The JUNO neutrino experiment aiming to determine the neutrino mass hierarchy is located in Jiangmen of South China and around 700 m underground. It uses liquid scintillators with 3 g/L 2,5-diphenyloxazole as the fluor and 15 mg/L p-bis-(o-methylstyryl)-benzene as the wavelength shifter. It has the density of 0.859 g/ml with around 88% of Carbon. The radius of the JUNO detector is R JUNO = 17.7 m. So, we take the detector area to be approximately A det ≈ πR 2 JUNO = 9.8 × 10 6 cm 2 . Taking an experimental observation time of t det ≈ 10 yr, the reach on the DMB mass is estimated to be 2.7 × 10 22 GeV.
From the experimental studies in Refs. [42,43], one can have the selection time of t select = 300 ns with the trigger energy around E PMT R = 70 keV and around 80% trigger efficiency. Similar to the estimation in Eq. (82), we have the projected limit from the JUNO experiment as In Fig. 7, we summarize the experimental reaches from Mica, Xenon1T, BOREXINO and JUNO for fixed model parameters m φ,0 = 0, λ φ = 0.013 and ω c = 50 GeV. Varying these model  (73)]. Based on the formation from a first-order EWPT and for the range λ φh ∈ [2,9], the average DMB masses vary from 1.1 × 10 24 GeV to 9.2 × 10 33 GeV.
parameter values does not change the signal curve much, especially in the log plot here. From this figure, one can see that a wide range of DMBs with a mass below around 10 22 GeV could be discovered by experiments with a large exposure. For DMB with a mass above 10 26 GeV or the Mica reach, there is no experimental probe at the current moment. New ideas to probe a heavy DMB are worth of exploring. For instance, one could use seismic data to search for heavy DMBs with a large scattering cross section [44,45].
Finally, we also comment on other bigger-size neutrino experiments. For IceCube with A det ≈ 10 10 cm 2 and t det ≈ 10 yr, the typical trigger energy is 100 GeV with the deep-core part as low as 10 GeV [46]. Only relativistic incident particles can generate Cherenkov lights, which is not the case for the non-relativistic DMB at hand. For the DUNE experiment with A det ≈ 1.0×10 8 cm 2 and t det ≈ 10 yr [47], it potentially can probe DMB mass up to 3×10 23 GeV. However, the energy threshold is a few MeV and still too high to measure the summed recoil energy for the DMB multi-hit events. Beyond the recoil energy from the DMB elastic scattering, the bound state formation from DMB capturing nucleons and nuclei could deposit much larger energy, 38 MeV for a nucleon and 38 × A MeV for a nucleus, which requires additional dedicated studies to estimate the experimental reach.

Discussion and Conclusions
The EWS-DMB is a special dark matter candidate because of its close relation to the Higgs potential. Different from the collider studies for the Higgs boson properties with the Higgs boson as the quantum field around the global vacuum with h = v, the restoration of electroweak symmetry inside a DMB provides us an opportunity to probe a wide range of the Higgs VEV from 0 to v. When the DMB interacts with the ordinary matter, the vacuum energy difference inside and outside the DMB can generate an effective "Higgs force" on nucleons or nuclei. As we discuss before, this force could also bind nucleons or nuclei inside the DMB.
As discussed around Eq. (58), the number of allowed bound states for heavy quarks inside a DMB is very large. After the primordial formation of DMBs from the early-universe dynamics, many of those states could be filled. As a result, the baryon number density inside a DMB could be very high and even to have the QCD in the deconfined phase, when the energy per baryon is below the proton mass outside DMB. If this is indeed the case, the electroweak symmetry inside DMB is really unbroken, not even corrected by the QCD confinement related chiral symmetry breaking. The EW sphaleron process is therefore active and can change baryons to leptons, which similar to the monopole in Grand Unified Theory that can induce catastrophic proton decays [48,49]. Although the Earth is not dense enough to stop a DMB, a neutron star may have a DMB stuck inside. The subsequent sphaleron-induced nucleon decays may change the properties of neutron stars and even evaporate them away.
For the DMB considered in this paper, the constituent dark matter particle is a boson. A similar study can be performed for a fermionic dark matter particle. The equilibrium of the low-temperature DMB state is then reached by the vacuum pressure and the degenerate fermion pressure. The situation is similar to the quark nugget in Ref. [3], except that the energy density of EWS-DMB has the electroweak scale and higher than the QCD scale in the QCD quark nugget. For the early-universe formation of DMBs, we have simply used the firstorder electroweak phase transition, which is very natural given the Higgs-portal coupling to the complex scalar dark matter particle. The stochastic gravitational waves could be another correlated signatures to cross check the scenario in this paper.
In summary, starting from the simple Higgs-portal dark matter model, we have shown that the non-topological soliton state could be the main appearance of dark matter. The electroweaksymmetric DMB is another type of macroscopical dark matter models and has its mass of 1-10 10 g, depending on the portal quartic coupling strength and if formed from the first-order electroweak symmetry phase transition. The radius of a spherically-symmetric DMB varies from 10 −9 to 10 −6 cm. The energy density of DMB is at the electroweak scale such that it can penetrate the Earth without stopping. The experiments with a large detector size or a long exposure time are ideal to search for DMBs. Indeed, the existing Xenon1T and BOREXINO or the future JUNO experiments may discover a DMB with mass from 10 −14 g to 0.1 g, based on multi-hit events.

A Number of DMB Nucleation Sites
In this appendix we derive part of the information necessary to estimate the expected DM number carried by DMBs produced during a first-order phase transition. Assuming that the net Q charge in a Hubble volume has been set by some earlier DM number-genesis mechanism, the average DM number for such DMBs is controlled by how many DMBs are produced within that Hubble patch, which we call N Hubble DMB . Here we simply take N Hubble DMB ∼ N nucl , where N nucl is the number of EWSB nucleation sites in one Hubble patch. To estimate N nucl , we follow a standard calculation (see Ref. [50] for example).
The bubble nucleation rate is controlled by the SO(3)-symmetric bounce action, S 3 (T )/T , where S 3 (T ) is the energy of the static SO(3)-symmetric critical bubble, and V eff (h, T ) is the 1-loop, finite-temperature effective potential given in Eq. (33). The bounce solution can be obtained by solving the temperature-dependent Euclidean equation of motion with boundary conditions h(∞) = 0 and h (0) = 0. For different values of λ φh and in Fig. 8, we show the bounce profiles at T f (see Table 1 for numerical numbers), approximately the temperature at the end of nucleation process. In general, the bounce profiles have a "thickwall" feature. As one increases the coupling λ φh , the size of the profile decreases. For the electroweak phase transition, we are interested in times only slightly after the critical temperature is reached. We parametrize the bounce action, in the vicinity of T c , as where the two parameters S 3 and α are functions of the model-parameter λ φh , 17 and we wrote the second equality in terms of the relative difference between T c and T : η(T ) ≡ (T c − T )/T c .  Table 1 for numerical values).
For the range 3 λ φh 6, we find that the index α ≈ 1.7 is relatively insensitive to λ φh . 18 The coefficient S 3 /T c , however, has a strong dependence on λ φh . We show the numerical values for α and S 3 /T c , as well as T c , for a few values of λ φh in Table 1. Performing a numerical fit, we have a power-law dependence for S 3 /T c as a function of λ φh Given S 3 (T ), one can estimate the nucleation rate per unit volume as [51] γ ≈ ζ T 4 where the T -independent ζ is expected to be of order one. Once formed, the nucleated bubble expands due to the vacuum pressure difference between the broken and unbroken phases. Due to its interaction with the particles in the plasma, a non-relativistic terminal velocity β w is reached [52]. The bubble wall is also preceded by a shock front that moves, for a strong firstorder phase transition, at a speed close to the speed of sound. Here, we use the speed of sound also for the bubble wall velocity, β w ≈ 1/ √ 3, to estimate the temperature at the end of the nucleation process, and then the total number of nucleation sites produced. transition, where φ = 0, m 2 φ,0 enters only at 1-loop order through the h-dependent m φ , while λ φ enters only at 2-loop order. We have also investigated the impact of other values on T c . We don't expect the estimates derived here to change significantly, as long as a strong first-order phase transition can be obtained. In particular, m 2 φ,0 cannot be too large, or else Φ would decouple from the plasma and one goes back to the SM result. 18 Ref. [50] obtains α = 2 based on an analytic high-temperature-expansion. Here, we have used the full one-loop finite-temperature potential which leads to a slightly different index.  Starting at time t c = 1.509g 19 when the plasma temperature equals T c , the fraction of space that remains in the EW unbroken phase is given by [51] f where the exponentiation accounts for the overlap of the bubbles [53]. Defining T f as the temperature (corresponding to time t f ) when the nucleation process is essentially complete by f unbroken (t f ) = e −1 , we calculate T f from the bounce action as parametrized in Eq. (86) as follows. Using the steepest descent or saddle-point approximation to evaluate the integral in Eq. (89) [50,54], one gets the following approximate relation where η f ≡ η(T f ), and we wrote γ(T f ) simply as γ(η f ) without changing the notaton [see Eqs. (86) and (88)]. Eq. (90) holds in the limit that S 3 | T f 1 and 20 where we used that we are interested in times (from the phase transition temperature, T c ) up to T f , and that η(T f ) 1. In the left panel of Fig. 9, we show η f as a function of S 3 /T c , for ζ = 1 and α = 1.7, obtained numerically from Eq. (90). We see that it can be described by 19 In the radiation dominated era. Also, close to the EW phase transition, T ≈ T c , the effective number of relativistic degrees of freedom, g * , does not suffer any abrupt change. M P is the reduced Planck Mass. 20 For the Hubble scale we have H = 0.331g 1/2 * T 2 /M P , where the effective number of relativistic degrees of freedom is between g * = 106.75 (the SM value) and g * = 108.75, depending on the bare mass parameter m 2 φ,0 . For concreteness, we use g * = 108.75. The number of nucleation sites within a Hubble patch at t f is approximately given by [50,54] where in the third equality we have used the fitted relation between S 3 /T c and λ φh in Eq. (87). This captures the dominant dependence on λ φh . Because of the large value in the power index, the number of nucleation sites can change five orders of magnitude even for λ φh changes by a factor of two. We also note that taking α = 2, the power index −1.7 in Eq. (92) changes slightly from −1.7 to −1.5 [54]. Finally, for our calculation of DMBs charges, it is convenient to know the λ φh dependence of phase transition temperatures, which has

B EWS-DMBs and Bound States
Consider adding a new (complex scalar) degree of freedom to the Higgs-Φ sector discussed in the main text. This will serve as a proxy for other "matter" and we will denote it by χ. Thus, the system to be considered in this appendix is given by where L is given by Eq. (1). Like the SM fermions and gauge bosons, χ interacts with the scalar sector only through the Higgs field. 21 We look for spherically symmetric solutions of the form where φ(r) and h(r) are real, while it is convenient to treat ψ(r) as complex. 22 They obey the boundary conditions We also assume 0 < E < m χ ≡ yv/ √ 2, so that ψ corresponds to a bound state in the Higgs well induced by φ, as well as the analogous condition ω 2 < λ φh v 2 /2 [see Eq. (11)], so that there is a soliton in the first place. This generalizes the case studied in the main text, and allows to discuss the issue of χ particles that get bound to the DMB.
The EOM are where V H (h) is the SM Higgs potential and V Φ (φ) was defined in Eq. (10). We see that, at the level of the EOM, the distinction between φ and ψ is mainly in λ φh vs y 2 (and potentially in their self-interactions, which play a secondary role in the present discussion). We make a notational distinction as a reminder that the y interaction is meant to mimic a Yukawa interaction. The 21 At loop-level, local interactions with Φ will be induced, but since our interest here is only to use χ as a proxy for the SM sector, we do not consider such terms. We also choose not to include a mass term or self-interactions for χ. This means that when H = 0, as happens inside a DMB, χ behaves as a massless particle and relativistic effects are important. We are including here only the kinematic ones. 22 We can chose to treat ψ as real, but then it is more convenient to include a factor of 1/ √ 2 in its definition to ensure canonical normalization. Treating it as complex makes the relation to the non-relativistic wavefunction more standard. However, note that here ψ is not exactly the radial wavefunction as usually defined in quantum mechanical central problems. In particular, ψ includes the factor Y 0 0 = 1/ √ 4π, and has mass dimension one, as for relativistic canonical normalization. important physical difference between φ and χ arises instead from the normalizations we impose. For φ we impose that the charge Q as defined by Eq. (3) be much larger than one: For ψ, we impose instead the relativistic normalization for "one particle".
These normalizations do not affect directly the EOM for φ or ψ, but they affect the Higgs EOM, saying that, even if λ φh ∼ y 2 , the Higgs well is sustained mainly by φ (there are many more particles in φ than in ψ). Then ψ can be thought as a bound state in this well. We are interested in the total energy of the system, denoted by and, in particular, on how it compares to the energy of an "empty" DMB, M , wih the same charge: We want to compare two solutions with the same Q: one where the soliton and a "free" χ are far away, so that the soliton mass is given by M Φ , and the solution where χ is bound to the soliton. Let us call φ 0 , h 0 the solutions to the EOM with ψ = 0, i.e. the DM soliton configurations we considered in Section 2. Since we want to keep fixed, we can consider the Legendre transform The EOM for φ and h, Eqs. (99) and (100) can be written as A bound state, i.e. ψ = 0, induces a perturbed φ-ball solution φ = φ 0 + δφ, h = h 0 + δh, and therefore where the EOM (108) and (101) have been used. We also have from Eq. (106) so we can finally write where we used Eq. (103). The second term corresponds to the backreaction of χ on the DM soliton when it gets bounded to it. One can get an estimate for δφ by considering the linearized EOM satisfied by the perturbations, but in essence one expects the backreaction of a single particle that gets bound to the Q 1 particle system to be small. Thus, We see that the change in the total mass is positive and approximately set by the energy of the bound state, as determined by the eigenvalue problem (101) in the unperturbed background for h. Note that to obtain the binding energy, we must compare the bound state mass to M (0) Φ +m χ , i.e. the total energy of the empty DM soliton plus a χ particle at rest at infinity (where h = v). This gives We can find the solution to the full EOM (99)-(101) numerically by proceeding iteratively as follows. We start from the 0-th order solutions to Eqs. (99) and (100) (with y = 0), that we are calling here φ 0 and h 0 . We then solve Eq. (101) in the fixed background h 0 using the "shooting method" to obtain a bound state wavefunction, ψ (n) 0 , obeying the desired boundary conditions, and corresponding to the n-th bound state. We then go back to Eqs. (99)-(100), inserting (the properly normalized) ψ (n) 0 as a fixed background, and solving the two-variable system as described in Section 2.1. Here we need to readjust ω = ω 0 + δω, where ω 0 corresponds to the y = 0 soliton, so as to keep Q fixed. This produces new solutions φ 1 = φ 0 + δφ and h 1 = h 0 + δh. The procedure is iterated to obtain the corrected n-th bound state in the h 1  Figure 10: Self-consistent solution with Q ≈ 8213 and containing a single bound particle of "vacuum mass" m χ = 174 GeV. The total energy is 49 GeV larger than the energy in the absence of the bound state (but with the same charge Q). The bound state energy (which corresponds to the ground state) is E 0 = 47 GeV. The small mismatch of 2 GeV is due to the backreaction. Note that the bound state wavefunction has been multiplied by 50. The model parameters are λ φh = 1, m 2 φ,0 = λ φ = 0.
background, which is then used to obtain (φ 2 , h 2 ), etc. We find that this procedure typically converges fast.
As an example, we consider the model defined by λ φh = 1, m 2 φ,0 = λ φ = 0, and the 0-th order soliton with ω = 50 GeV, which has Q ≈ 8213 and a radius of about R = 0.06 GeV −1 (this is an example of a DMB). Taking a particle with vacuum mass m χ = 174 GeV (corresponding to y = 1/ √ 2) one finds a ground state with energy E 0 ≈ 47 GeV, while the total mass M Φ [ψ] is about 49 GeV heavier than for the empty DMB. We see that the backreaction amounts to about 2 GeV. ∆M Φ = 49 GeV is the available energy if χ "disappeared". The corresponding profiles are shown in Fig. 10.
Consider a second, lighter particle with m χ ≈ 49 GeV (y = 0.2). One finds ∆M Φ = 35.93 GeV and E 0 = 32.48 GeV (not as deeply bound as χ above). The decay χ → 2χ is allowed outside the DMB. But if initially χ is bound inside the DMB in the ground state, and if we assume that a state with two bound χ particles has ∆M Φ = 2 × 35.93 ≈ 72 GeV (i.e. that the effects are approximately additive), we see that energy conservation would not allow the decay to proceed inside such a soliton.

C Bound States in a Higgs Potential Well
If one neglects the backreaction effects discussed in Appendix B, we can get a handle on the spectrum of bound states in a DM soliton by solving Eq. (101) in the fixed h background of the unperturbed DM soliton. With some abuse of notation, it will be useful to state the starting point here as follows: where the source term J is added for convenience in taking various limits below. Furthermore, we will assume that the H background corresponds to a (spherically symmetric) 3D potential well of size R: We shall treat the EOM that follows from L χ above as a single particle wave equation. When looking for solutions of the form specified in Eq. (95), it takes the form where ψ is regarded as complex. The structure of this equation is the same as the corresponding eigenvalue Schrödinger problem, using a non-relativistic Hamiltonian with a potential, V , defined by 2M V = y 2 |H| 2 + J, and taking E 2 /2M → E NR , where E NR stands for the non-relativistic energy eigenvalue, and M would correspond to the mass of the particle in the corresponding nonrelativistic problem. For potentials which are piecewise constant, the solutions can be expressed in terms of spherical Bessel functions. We are interested in bound states, which means that E 2 must be below m 2 + J and above J (for particles, as opposed to antiparticles, we also take E > 0). Here we defined m ≡ yv/ √ 2, the mass of the χ field in the normal EWSB vacuum, when J = 0.
For the s-wave case, we have ∇ 2 ψ = r −2 d dr (r 2 ψ ) and 23 ψ(r) = A r sin(Kr) for r < R B r e −κr for r > R , where are both real. Matching the functions and their first derivatives at r = R and dividing the two relations leads to the equation for the spectrum − cot(KR) = κ K .
Setting J = 0 leads to Eq. (51), describing the s-wave spectrum of a particle that gets its mass only from the background H. A similar analysis can be applied to higher partial waves. Let us establish the condition for the existence of at least one bound state. Writing E ≡ m + ∆E, the first possible bound state appears when ∆E < 0 with |∆E| → 0, i.e. just at threshold. In this case, Eq. (51) reads If mR = π 2 + for 0 < 1, the cotangent is negative and there is a solution for ∆E ≈ −( 2 /2) m, consistent with the approximation |∆E| m made above. Thus, we see that the (minimum) threshold radius for the existence of a bound state is R th = π/(2m). Note that this problem is not exactly the same as a non-relativistic 3D well, since we are taking into account here the fact that the mass vanishes inside R (recall that m above simply stands for yv/ √ 2). 24 So, even though the χ particle inside the well is "massless", it can be bound. For larger R, more bound states appear with an asymptotic spacing In the large R limit, the ground state has energy of order E ≈ 0, which corresponds to a binding energy of order −∆E ≈ m.
Let us now consider nucleons, whose mass is essentially independent of v. The relevant terms in the (physical) nucleon Lagrangian read For a "scalar nucleon" we can then set J = m 2 N and |H| 2 → 2|H| 2 − v 2 in Eq. (114), so that 25 where we also relabelled y → y hN N . This leads to Eq. (53) in the s-wave case. Setting E ≡ m N + ∆E, one can see that the first bound state with ∆E ≈ 0 appears for 26 where the second equation is specific to a nucleon. If we considered a nucleus with A nucleons one gets an additional factor of 1/A. Again, as R increases one gets additional bound states, and in the large R limit, their (s-wave) spacing is of order π/(2R). 24 We note that the threshold R is mapped into the nonrelativistic result R th = 1 √ 2mV0 π 2 with the identification V 0 = m/2, as expected from the discussion after Eq. (116), taking M = m. 25 Alternatively, one can make H → √ 2H − v, adding an explicit minus sign in front of the y 2 hN N ( √ 2H − v) 2 term in Eq. (114) to model a well instead of a barrier. 26 The threshold R coincides with the one found in the nonrelativistic analysis of a particle of mass m N in a potential well of depth 2m N V 0 = (y hN N v) 2 . See footnote 24.

D Simple Example of Scattering Against a Heavy Object
In this appendix we illustrate the method used in the main text to compute scattering cross sections with a simple toy model. We take the interaction terms where Φ and h are scalars (complex and real, respectively), and χ is a Dirac fermion. We assume that Φ is much heavier than χ. In this situation, the approach is to first find the h field induced by Φ, then consider the scattering of χ in such an h background. Consider an "elementary" Φ configuration in its rest frame given by where M Φ is the Φ mass. In the h EOM this enters as a source which can be interpreted as the superposition of point-like sources The Φ charge contained in volume ∆V in the configuration of Eq. (127) is so we can write Eq. (129) as The static and spherically symmetric h field induced by the point-like source is We can now turn to the fermion, which sees a potential V (r) = g h(r). In the Born approximation, the differential scattering cross section is given by where q ≡ | q|. The total cross section for q = 0 is Let us now repeat the computation following a textbook quantum field theory approach. The tree-level invariant amplitude for Φχ → Φχ scattering is where for elastic scattering q 2 = q µ q µ = −| q| 2 . In the heavy M Φ limit, i.e. taking the Mandelstam s ≈ M 2 Φ , we have Using uu = 2m χ , and taking q = 0, we get which reproduces Eq. (134) with Q = 1. This corresponds to the fact that in the QFT computation we are scattering a single Φ particle, and establishes how to implement this concept in the classical language, or generalize it for aggregates of particles. Note also that the only real assumption we have made is that M Φ m χ . While we assumed that q = 0 to compute the total cross section, this was done only for simplicity. Indeed, one can see from Eqs. (136) and (135) that with f (E, θ) given by Eq. (133), for any q. Although clearly the QFT approach is significantly more efficient in this example, the "fixed background" approach is better suited for the scattering against heavy extended objects, which is more naturally discussed in configuration space.