Hidden Sector Monopole Dark Matter with Matter Domination

The thermal freeze-out mechanism for relic dark matter heavier than $O(10-100 $ TeV$)$ requires cross-sections that violate perturbative unitarity. Yet the existence of dark matter heavier than these scales is certainly plausible from a particle physics perspective, pointing to the need for a non-thermal cosmological history for such theories. Topological dark matter is a well-motivated scenario of this kind. Here the hidden-sector dark matter can be produced in abundance through the Kibble-Zurek mechanism describing the non-equilibrium dynamics of defects produced in a second order phase transition. We revisit the original topological dark matter scenario, focusing on hidden-sector magnetic monopoles, and consider more general cosmological histories. We find that a monopole mass of order ($1-10^5$) PeV is generic for the thermal histories considered here, if monopoles are to entirely reproduce the current abundance of dark matter. In particular, in a scenario involving an early era of matter domination, the monopole number density is always less than or equal to that in a pure radiation dominated equivalent provided a certain condition on critical exponents is satisfied. This results in a larger monopole mass needed to account for a fixed relic abundance in such cosmologies.


Introduction
The period between the end of inflation and the beginning of big bang nucleosynthesis (BBN) is a natural period for the production of dark matter (DM), though it is currently inaccessible to observations. The most popular dark matter candidate has traditionally been a weakly interacting massive particle (WIMP), produced in the right abundance by thermal freeze-out in the standard thermal history of radiation domination (RD) between inflation and BBN. This standard picture is now increasingly strained, with certain models excluded by indirect searches over much of the cosmologically interesting range for the WIMP mass [1,2]. Nonthermal production mechanisms, which depart from the assumptions of local thermal and chemical equilibrium of dark matter with Standard Model particles in the early Universe, and/or radiation domination, have become more widespread [3].
Spontaneous symmetry breaking in the early Universe prior to BBN provides a natural mechanism to produce interesting objects through an out-of-equilibrium process. Specifically, symmetry breaking via a second order phase transition can produce a large density of topological defects via the Kibble-Zurek mechanism (KZM) [4][5][6], and their density, can "leave an immediate imprint on the Universe and will be critically important" [6]. While the KZM theory was developed some time ago, it is only recently that the theory has received firm experimental support, at least for describing classical second order phase transitions, as certain key predictions of the theory have been confirmed in laboratory settings. In particular, the scaling of the density of topological defects with respect to the quenching rate has been verified in a number of two-and three-dimensional materials [7][8][9][10][11].
What is of focus here, is that the KZM is a plausible nonthermal mechanism for the production of an interesting class of dark matter candidates dubbed topological dark matter [12]. A key finding of [12] is that in this scenario, the dark matter mass must be of O(PeV) scale to obtain the correct relic abundance. Our main motivation for the present work is to explore the robustness of this finding, when other cosmological histories in the early Universe are considered. Topological dark matter is studied by [12] in the context of a standard thermal history, in which the phase transition that produces topological defects occurs during a radiation dominated era, and where the temperature of the symmetry breaking and visible sectors are assumed for simplicity to be equal. We explore this scenario in several different directions. We allow for an intervening phase of matter domination (MD) in the early Universe, during which the symmetry breaking occurs. We also allow the symmetry breaking sector to have a temperature different than that of the visible sector (VS) of Standard Model particles. For since the two sectors interact only very weakly, if at all, there is no reason to expect them to have the same temperature.
Phases of early matter domination (EMD) in the period between inflation and BBN are a generic prediction of early Universe string constructions and are commonly achieved via moduli which acquire a pressureless equation of state and drive the Universe toward matter domination before their eventual decay [13][14][15][16][17][18][19]; for a review see [20]. An early matter dominated era can also easily happen when a decoupled massive particle comes to dominate the energy density for some time before decaying and subsequently reheating the Universe. We will consider an era of early matter domination to be caused by either a modulus or a decoupled particle, and allow the phase transition to occur anywhere before, during, or after this era. Our cosmological scenario actually consists of two hidden sectors: a sector driving an early matter domination phase; and a second sector with the symmetry breaking by a second order phase transition. Couplings between these two sectors would be interesting to explore -leading to a more complicated cosmological historybut we do not do so here, simply to avoid over complicating the narrative.
While the original work on topological dark matter [12] considered the production of domain walls, strings, monopoles, or skyrmions, here we focus for simplicity on the case where the produced defects are magnetic monopoles, charged under an unbroken U (1) left over after the phase transition. 1 The abundance of magnetic monopoles charged under the U (1) of electromagnetism is constrained by observations, such as the Parker limit, to be less than that required for it to account for all of the DM [21,22]. We will therefore avoid such constraints altogether in this work by considering the simplest scenario in which the monopoles are not charged under electromagnetism, but instead charged under a hidden sector U (1), and further, that the hidden sector U (1) does not kinetically mix with electromagnetism, so that monopoles of the hidden sector do not couple to (visible sector) electromagnetic fields. 2 Our scenario begins in a radiation dominated phase after inflation, where we allow for the dominant energy component to be radiation in either the visible or hidden sector. As the Universe expands, each sector cools independently of the other, and we enter an early matter dominated phase caused by a modulus or by a heavy particle which has decoupled from either sector. As this phase proceeds, the dominating field continually decays into radiation in the visible sector, until the decay completes (at reheating) and we transition back to a radiation dominated phase of Standard Model particles, leading to the standard cosmology at the onset of BBN. We suppose that a secondorder phase transition occurs in the hidden sector as the temperature in the hidden sector drops below some critical temperature T (hid) C , resulting in a significant production of magnetic monopoles in the hidden sector due to the Kibble-Zurek mechanism. We allow the phase transition to occur at any time in the pre-BBN thermal history of our scenario. We a posteriori neglect any subsequent annihilations of monopoles due to their high mass (PeV and above) and consequently low number density. As mentioned above, we also do not consider any non-gravitational interactions between the sectors, other than that which provides the decay that reheats our Universe.
Our main results are shown in Figures 4, 7, and 8. We generally find that hidden sector monopoles in the mass range O(1-10 5 ) PeV can be dark matter candidates, with values for the monopole mass giving rise to the current dark matter relic abundance correlated with other particle and cosmological parameters. Furthermore, a long intervening era of matter domination in the early Universe significantly increases the hidden sector monopole mass needed to obtain the current relic abundance, compared to a purely radiation-dominated history, provided that the critical exponents, defined below, satisfy 2ν ≤ 1+µ. An analytic argument for this observation is presented in Section 4, which is also confirmed by our numerical results given in subsequent sections.
We begin with an overview of monopole production via the Kibble-Zurek mechanism in Section 2, followed by an overview of a cosmological history involving EMD in Section 3. In Section 4, we present analytical forms for the monopole abundance in the presence of an EMD phase, including monopole production before, during, and after EMD. We then present numerical results for the cases of EMD by a modulus or a heavy decoupled particle in Sections 5 and 6 respectively. Section 7 shows the monopole mass and cosmological parameters that give the correct present-day relic abundance for dark matter, using an analytic approximation that we show well-describes the relic abundance obtained using numerical methods. We conclude with a brief discussion, including a summary of important caveats to our work, in Section 8.
A number of detailed results are summarized in several Appendices. We include a table of notation in Appendix A. Appendix B describes the relation of a key cosmological parameter in our work -the length of the matter-dominated phase -to other defined cosmological parameters. Appendices C and D gather usual formulae for the decoupling of a relativistic particle, and Appendix E gives the constraint on cosmological parameters from requiring that a matter-dominated phase caused by a decoupled particle lasts at all.

Brief review of Kibble-Zurek mechanism theory
We now summarize the theory of the Kibble-Zurek mechanism describing the non-equilibrium dynamics of topological defects produced in a second order phase transition. We refer the reader to the original references [4][5][6] and recent review [26], which give several reasons for why (2.2) shown below gives the typical distance scale between topological defects.
In the KZM theory, a system is assumed to be driven through a second-order phase transition at temperature T C by a quench that importantly, is assumed to be of a finite timescale; it is neither instantaneous, nor extremely long. In a cosmological context, the quench is driven by the cosmological expansion of the Universe itself, a point we return to below.
If the quench is slow enough, the system has time to quasi-equilibrate and therefore as t → t C the correlation length continues to grow with some critical scaling, namely for some critical exponent ν.
The key point is that there is a time scale t * prior to the phase transition, such that for times t > t * , the correlation length exceeds the sound horizon. Subsequent to that time, the quench is fast compared to the timescale over which the system can respond. According to the KZM theory, after this cross-over time t * , fluctuations become frozen, and therefore ξ(t * ) sets the scale of the topological defects, namely [4][5][6], where u(t) = u 0 (t) µ−ν , for a critical exponent µ and typical velocity u 0 , is the characteristic velocity of perturbations in the system. 3 The characteristic correlation time scale τ (t) is then for a typical timescale τ 0 = ξ 0 /u 0 . We now arrive at the main prediction of the KZM theory. For this finite speed quench, the frozen correlation length is then predicted to be with approximately one topological defect (monopole) produced per correlation volume ξ(t * ) −3 [4][5][6]. The size of the frozen length scale is set by physical properties in ξ 0 , τ 0 , and the critical exponents, and by the timescale of the quench τ Q set by either the laboratory conditions or by the Hubble expansion rate, depending on the context. It follows that the number density of point-like defects in D = 3 spatial dimensions is 4 This scaling of defect density has been experimentally confirmed in a number of two and three dimensional condensed matter systems, such as 3-D ferroelectric crystals [7], 2-and 3-D Bose-Einstein condensate gases [8][9][10], and multiferroic hexagonal manganite crystals [11]. A critical dynamical assumption leading to these predictions is that fluctuations in spatial regions separated by more than this correlation length are randomly oriented and, subsequent to the above cross-over time, independent of each other. While this is a reasonable expectation for a classical phase transition, Zurek raises a caveat for systems such as the normal-to-superfluid transition in 4 He in which quantum mechanical effects are all important [4]. Namely, correlations between regions separated by several correlation lengths may only appear to be random and independent, but in fact could be secretly strongly correlated due to conservation laws (for the vortices studied in [4], notably angular momentum), in analogy to spin correlations in EPR experiments. Should this situation occur, the predicted topological number density would be smaller and these estimates for the cosmological relic density would need to be revisited [4].
But recent experimental results do suggest that -at least in the case of vortex formationdefects are indeed random and independent, reaffirming the KZM expectations. For the KZM theory also makes some statements about this randomness, as it specifies how the net winding number of vortices W in a fixed spatial region of circumference C should scale with the correlation length. Specifically, the typical absolute value |W| and dispersion W 2 are both predicted to have the same scaling at large |W| 1, namely W 2 ∼ |W| ∼ C/ξ, whereas at small winding number the KZM predicts different scaling laws for the absolute value and dispersion of W [4,27]. In both limits the KZM predictions for these two quantities have been dramatically confirmed in 3-dimensional ferroelectric crystals [7].
In a laboratory setting, in the non-relativistic mean field approximation (i.e., Landau-Ginzburg theory), the potential part of the free-energy of a system described by an order parameter φ is approximated by the Landau-Ginzburg potential, with the time-evolution of φ approximately described by the Gross-Pitaevskii equation, which is first order in time. This leads to the critical exponents µ = 1 and ν = 1/2, predicting ξ(t * ) ≈ ξ 0 (τ Q /τ 0 ) 1/4 . But in a relativistic quantum field theory context the scaling laws are different because the equation of motion for φ is second-order in time. For example, in a cosmological context the equation of motion for φ leads to the critical exponents µ = ν = 1/2. Here then, [12]. As noted above, when the phase transition occurs in an expanding Universe, the quench time can be re-expressed in terms of the Hubble rate at the critical time as H −1 C . To see that, first note that the quench is characterized by where T (t) is the time-dependent temperature of the system. Close to the time of the phase transition t C , this quantity scales linearly with time, which also defines the quenching time-scale τ Q . For example, in a cosmological context where the scale factor a increases as a(t) = (t/t C ) p , p = 2/3 (1/2) for MD (RD), then with t ≡ t C − ∆t, |∆t| t C , (t) = p∆t/t C and τ Q = t C /p, or in other words, That is, the characteristic time-scale τ Q of the quench is always given by the Hubble parameter at the time of the phase transition, generalizing from the pure RD scenario given in [6] to more general equations of state. We take the initial correlation sizes to be set by the mass m σ of the σ particle, which for a pure scalar φ 4 theory at weak coupling is given by m σ λT C /4 [28]. That is, [12]. Although µ = ν = 1/2 is the prediction for the critical exponents in the approximation that the second-order phase transition is described by a weakly coupled scalar field, for our analysis we consider more general values for the critical exponents. 5 In terms of cosmological quantities, the frozen correlation length is then regardless of the type of dominant energy density (matter or radiation), with the understanding that the temperature dependence of the Hubble parameter when the system is at the critical temperature, H C , does depend on the form of the dominant energy density component.
After the phase transition is complete, the monopole number density is n M ≈ ξ(t * ) −3 and the comoving number density is fixed as their abundance simply redshifts through the remaining history of the Universe. We will neglect any subsequent annihilations of monopoles because the masses needed to account for the entire current DM abundance will turn out to be quite high, with correspondingly low number densities. 6 For a general second order phase transition, quantum or classical, in the KZM theory the frozen correlation length setting the density of topological defects depends only on the critical temperature of the phase transition, the typical timescale of the quench, and the critical exponents. For a classical Landau-Ginzburg second order phase transition, however, the mass of the defect -here the monopole mass m M -is not independent of the critical temperature. For a 't Hooft-Polyakov monopole, m M = hT C , with h the magnetic coupling 2π/e h , and recall that φ ∼ T C . Thus for a classical phase transition, the monopole mass and critical temperature are parametrically at the same scale. Throughout this work we will assume the monopoles are produced in the early Universe by a classical second order phase transition, so the implied relation between the critical temperature and monopole mass is an important caveat to many of our results.
But such a mass-temperature (m − T ) relation is not expected to be true in general. On the contrary, one expects the monopole mass and critical temperature to be unrelated. The N = 2 Seiberg-Witten theory [32,33] is a prominent example of this kind, where near certain points on the moduli space the low-energy theory contains nearly massless composite particles charged under a magnetic U (1). Here one would like to know whether the theory ends up near these points as the theory is cooled through the phase transition, and what the order of the transition is. For the former question, the answer is affirmative, at least in the pure N = 2 SU (2) theory [34]. The latter remains an open question.
Because of this expectation, we will indicate which of our results are independent of any assumption about a m − T relation. The most important of these is the ratio of the monopole number density to photon entropy density, such as (4.2), (4.4), and (4.5) given below. In the low-density limit where monopole annihilations are negligible, these depend only on the critical temperature but not the monopole mass.
As previously mentioned, we will also vary the critical exponents µ and ν away from the Landau-Ginzburg value of 1/2, as a guide to future work.

Summary of the cosmological history with an early matter-dominated era
In order to proceed, we must address the relationship between the Hubble expansion rate and the temperatures of the different radiation components of the Universe. In this section we therefore introduce the general expansion history we will be considering, define terminology, and obtain relations between the Hubble parameter and key parameters during the different eras prior to reheating.
First, we begin with radiation domination (RD) by either the hidden or visible sector (or any combination) some time after inflation, with other energy densities comparatively negligible. In this era, the Hubble expansion rate is given by where the second equation implicitly defines the factor f ≡ ρ as the ratio of the radiation energy densities of the visible and hidden sectors. Also, T (hid) is the temperature of the HS, g (hid) * is the number of relativistic degrees of freedom in the HS at temperature T (hid) , and M P ≈ 2.4 × 10 18 GeV is the reduced Planck mass. In this period, the factor (1 + f ) is well approximated by its initial value (1 + f i ) regardless of the distribution of initial radiation among the two sectors, and we will make this substitution when using (3.1) below.
We consider the visible and hidden sectors to have independent temperatures, each with their own g * factors depending on the specific particle content (Standard Model for the visible sector), and we could have equivalently expressed (3.1) in terms of visible sector quantities. The g * factors of course depend on the temperature of their respective sector, but we will treat g (hid) * as roughly constant at high temperatures in order to avoid overly specifying the details of the HS.
We achieve early matter domination (EMD) through the presence of a scalar modulus, or by the decoupling of a heavy particle from either the hidden or visible sectors during this initial RD phase. In both cases we refer to the modulus and the heavy particle as Φ, and based on the context, there should not be any confusion. We assume that Φ couples to lighter particles through higher dimension operators suppressed by the Planck scale, with a decay rate where m Φ is the Φ mass. We have also included a possible loop factor α in the case that Φ decay occurs predominantly through a loop, but we will set α = 1 throughout unless otherwise noted. The decay is complete when H ≈ H RH ≡ Γ Φ , which marks the approximate time of reheating, and we avoid having significant amounts of left over hidden radiation by requiring Φ to decay predominantly to the Standard Model particles, where T

(vis)
RH is the visible sector temperature at reheating, and g (vis) * RH is the number of relativistic degrees of freedom in the visible sector at this temperature. In order to preserve standard Big Bang Nucleosynthesis (BBN), the visible sector reheat temperature must be larger than O(10 MeV).
The ratio of the visible sector radiation energy density to that of the HS at reheating, denoted by f RH , depends on the duration of the EMD phase as well as the initial factor f i , but is typically large due to our visible sector reheating requirement, and thus always satisfies f RH > 1 and f RH > f i (this statement is demonstrated in Appendix B). This conclusion, together with our assumption that Φ predominantly decays to SM particles, ensures that the temperature of the HS at reheating, T (hid) RH , is correspondingly always smaller than that of the visible sector. We also point out that this ratio remains fixed after reheating due to the absence of any further decays. From (3.2) and (3.3), we additionally see that a given choice for the visible sector reheat temperature and α determines a corresponding Φ mass.
In order to have a well defined EMD phase, we assume the energy density of Φ is large enough to dominate well-before reheating. During EMD, the scaling of the Hubble rate with the visible sector temperature is altered from a typical MD redshift relation because the visible sector is fed by the decay of Φ; however, from entropy conservation, the scaling of H with the HS temperature remains unaffected: H 2 ∝ h * (hid) T (hid) 3 . Based on the initial energy density of VS radiation, there can be a phase of ordinary redshift for the VS temperature even during EMD, but once the effect of the decay wins over this dilution, the relation becomes (see (20) of [35] for a derivation): This relation is always true just before reheating, but may not start until deep within the EMD phase if the initial VS radiation energy density is large. 7 At the end of the EMD phase, once reheating completes, we enter the RD era with the Hubble rate given by where the factor f RH is large such that the visible sector is dominant, thus recovering the standard thermal history leading up to BBN.

Monopole production with an era of early matter domination
Recall that we are interested in producing monoples during a second order phase transition occurring in a hidden sector, so the critical temperature appearing in (2.10) refers to the temperature of the hidden sector at the critical time. In this section we address monopole production in the context of the thermal history presented in the previous section. The effects of EMD on the monopole abundance can be understood regardless of the mechanism for establishing MD in this early period, and we obtain analytical expressions below that do not depend on the identity of the field Φ. In addition to the start time of EMD, what matters is that the dominant energy density component decays to visible sector radiation at a rate Γ Φ , thus setting the end time of EMD. The overall effect is to slow the redshift of visible sector radiation relative to the HS such that only the visible sector is dominant after EMD even if it was not initially. Because we only consider HS magnetic monopoles, this offset in the visible sector and HS temperatures generally results in a lower number density of monopoles of a given mass, where the magnitude of the offset is determined by the duration of EMD and the initial abundances of visible and hidden radiation. We label the start of EMD by H = H MD , with visible and HS temperatures T respectively, and the end of the EMD phase occurs when H ≈ Γ Φ . Recall that the visible sector reheat temperature, which we restrict to be larger than O(10 MeV) such that reheating occurs before BBN, is the primary parameter that determines the end of EMD.

Case I: phase transition occurs before EMD
We will start with the case where the HS phase transition occurs in the RD period before EMD, resulting in a frozen monopole number density that is redshifted through the remainder of the RD phase as well as the full EMD period. This results in considerable dilution and a need for higher monopole masses in order to maintain a fixed contribution to the energy density of the Universe. Using (2.10) and recalling that the number density of monopoles produced in the phase transition is approximately one per correlation volume, we have (see Appendix A for a table of notation) where the first factor in parentheses on the right-side accounts for the redshift of the monopole number density from the critical time to the start of EMD, and the second factor gives the redshift from the start of EMD to reheating. We have also defined a MD and a RH to be the scale factors at the onset of matter domination and at reheating, respectively. At this point we do not need to redshift any further, and can obtain a fixed comoving abundance by normalizing by the visible sector entropy density at reheating, as both number density and entropy density dilute as the cube of the scale factor once the significant entropy production from reheating stops. This leads to The factor h (vis) * tracks the visible sector relativistic degrees of freedom for entropy and is nearly equal to g (vis) * for the high temperatures in our scenario as well as the low temperature today [36,38] (it is evaluated at reheating in the expression above, as indicated by the subscript). Note that the Hubble rate at the critical time is given by (3.1).

Case II: phase transition occurs during EMD
If the phase transition occurs during the EMD phase, the frozen monopole number density only redshifts through the remaining duration of EMD, and we have Again normalizing to the visible sector entropy density at reheating, one has The dependence of H C on the HS temperature is that of ordinary MD redshift, while the relation to the visible sector temperature is more complicated, for it depends on how much visible sector radiation was present at the onset of EMD. If the visible sector energy density at H = H MD is greater than the subsequent contribution from the decay of Φ at H = H MD , then to evaluate H C one will need to include the effect of a period of ordinary MD redshift for the visible sector temperature as well. Once the decay contribution takes over well within the EMD phase, we have the relation (3.4). We note that this modified scaling can begin much earlier, even before EMD, if the initial visible sector radiation energy density is small.

Case III: phase transition occurs after EMD
Finally, if the phase transition occurs in the RD period after reheating but still before BBN, so as to leave the later evolution of the Universe unchanged, the abundance can be evaluated directly at the critical time, without need of redshifting: This expression is also valid for a thermal history that does not involve EMD at all, where the HS radiation energy density is lower than or equal to that of the visible sector by a constant factor, as both energy densities simply redshift with time. The Hubble rate at the critical time is given by (3.5) in terms of visible sector quantities, but is easily related to the corresponding HS quantities by multiplying by the square root of the constant factor. Finally, we note that all of the results in these three subsections are independent of any possible relation between the monopole mass and the critical temperature.

Monopole production: analytic approximation at boundaries
In this subsection we obtain analytical expressions to better understand the effect of EMD in more detail. The three cases of monopole production described above are separated by production at the start and end of EMD, and we can easily obtain expressions below for the monopole abundance corresponding to these boundaries.
For production at the start of EMD, the HS temperature at the critical point is T with corresponding H C = H MD . From (3.1) and (4.4), we obtain the frozen abundance of monopoles at reheating: Aside from the parameters of the phase transition, the final abundance is determined by the visible sector reheat temperature, the initial ratio of visible sector to HS radiation, and the monopole mass.
Monopole production at the end of EMD corresponds to a HS critical temperature of T with the implicit relation between Γ Φ and T (vis) RH given by (3.3). Note that this expression does not depend on the initial ratio of radiation energy densities as it only involves the time of reheating.
Requiring EMD to start before reheating, these two expressions for production at the boundaries of EMD significantly constrain the allowed parameter space. For a realistic scenario, even the shortest EMD period will have a finite duration such that EMD is well defined, ensuring that we never quite access the limiting case where the start and end of EMD are coincident. This case, rather, corresponds to the absence of EMD altogether.

Present-day hidden sector monopole abundance
We will now obtain the present day relic abundance of monopoles. In the three main cases of monopole production -before, during, or after EMD -as well as the two boundary cases of production at the start and end of EMD, the parameters µ, ν, and λ, are determined by the details of the phase transition, as is the ratio The ratio x M is the magnetic coupling, and typically has a value of O(10) [12] -we will assume x M = 50 in our numerical results below. The current abundance of monopoles, expressed as a fractional energy density Ω M h 2 , is related to the frozen abundance provided in the previous sections by where Ω γ h 2 = 2.47 × 10 −5 corresponds to the current photon energy density, ρ (vis)

Also, h
(vis) * 0 = 43/11 = 3.91 is the present-day era total entropy density pre-factor, assuming three massless species of neutrinos. The subscript '0' labels the current era, and the final term labeled by '(EMD)' refers to any one of the five above cases. The subscript 'RH' on the final term means this quantity is evaluated at reheating if the phase transition occurs before reheating, whereas in the circumstance that the phase transition occurs after reheating, 'C' means the quantity is simply evaluated at the time of the phase transition. In order for monopoles to constitute all of dark matter, the value of Ω M h 2 must reach the observed value of 0.12 [39].
For comparison with our numerical results in subsequent sections, analytical expressions for Ω M h 2 can be obtained in the three main periods of our scenario by noting that where the cases refer to monopole production before, during, or after the EMD phase. In the period before EMD, we have the RD relation (3.1), while in the period after EMD we have this same functional form, but with a different constant factor offsetting the visible sector and HS radiation energy densities. The expression for H C during EMD is obtained by using entropy conservation in the hidden-sector radiation, together with redshifting during the EMD era between the start of EMD to when the temperature of the hidden sector reaches T T (hid) ∝ H 2/3 . Next, using (4.2), (4.4), (4.5), (4.8), and (4.9), one obtains analytical estimates for the monopole abundance produced in the three periods by direct substitution 8 In the model-independent discussion of this section, the Hubble rate at the onset of EMD has been an independent parameter. In Sections 5 and 6 below, where we address two examples for establishing a period of EMD, we provide expressions for H MD in terms of the underlying model parameters.
It is useful to extract the functional dependence of the energy density of monopoles on the monopole mass, produced during any of the three periods of before, during, or after EMD. From (4.10)-(4.12) above, we have  Here we have factored the dependence of the energy density on the mass into an explicit factor arising from the mass itself, and an implicit factor due to the number density. The RD case applies to monopole production both before and after EMD, and we have again assumed a constant factor, x M , between the monopole mass and T (hid) C . Note that in general, the type of cosmology in which the phase transition occurs -here either an EMD or RD era -affects the monopole energy and number densities through a different power-law dependence on the critical exponents.
Before moving on to consider specific scenarios for establishing EMD, we can see that, depending on the relative sizes of the critical exponents, the presence of an intervening EMD phase in the period before BBN can push the preferred monopole mass for DM higher than in a purely RD equivalent. For the two prefactors in (4.13) are not the same in each case. Fixing the phase transition parameters (µ, ν, λ, and x M ) as well as the monopole mass, m M , we must first identify the equivalent RD scenario, which comes down to specifying the constant factor f (RD) between the VS and HS radiation energy densities in the RD scenario. We obtain this by decreasing the duration of EMD until we arrive at the limiting RD scenario to use for comparison. If EMD is preceded by a period of RD by the VS, the limiting scenario is one which preserves the initial ratio of VS-to-HS radiation: f (RD) = f i . However, if HS radiation is dominant before EMD, the limiting case is one of f (RD) = 1 because we wish to avoid RD by the HS at the onset of BBN. In short, (4.14) and consequently, To proceed, for all three cases we define the ratio of the scale factors at reheating and the onset of the EMD phase to be which we show in Appendix B (see (B.6)) to be equivalent to f (1 + f i )e f , and because of (B.8), is always larger than f (RD) by the factor e f > 1, so long as Φ preferentially decays to the VS. The factor e f is fixed for a given EMD phase, regardless of the value of f i or the timing of the phase transition.
Using (4.2)-(4.5), and recalling that the HS temperature redshifts as T (hid) ∝ a −1 in all periods of our scenarios, be they EMD or RD, we arrive at the ratio of the current monopole abundance between an EMD and a pure RD scenario: (case III : after) (4.17) As with the previous expressions (4.2), (4.4), and (4.5) given above for the ratio of monopole number density to visible sector entropy density, in deriving these equations we have not made use of any relationship between the monopole mass and the temperature of the phase transition.
In all three cases, the products involving f 's and the critical exponents are the ratios of the monopole number densities produced at the critical time between the EMD and RD scenarios. We note that since T (hid) C and λ appearing in the correlation length (2.10) are fixed between the two scenarios, this ratio is simply given by the ratio of Hubble parameters H C . In the first two cases, we normalize the monopole number densities by the VS entropy density at the time of reheating (when the VS temperature is equal to the reheat temperature), accounting for the redshift factors, while in the third case, because monopole production occurs in RD after EMD, there is no need for redshifting, and we normalize by the VS entropy densities at the critical time. The factor of 1/e 3/4 f , in the first two cases, is the ratio of the redshift factors from the time of monopole production to the time when T (vis) = T (vis) RH between the EMD and RD scenarios respectively, while in the third case, it, along with the terms involving the relativistic degrees of freedom, comes from the ratio of entropy densities at the critical time between the two scenarios. Note that the relativistic degrees of freedom in the VS can be different at the critical time between the EMD and RD scenarios because it is the HS critical temperature, not the visible, that is the same across the scenarios.
We note in the limit of no EMD phase, the above expressions for the three cases smoothly go over to Ω . For cases I and III this statement is readily apparent, since in this limit . To see that for case II requires one additional remark. By definition of this scenario, (T , and also (T • For case I, of monopole production before EMD, the right-side of (4.17) is always less than one. To see that, first focus on the ratio of f factors appearing in (4.17). Recall that (1, f i ), and therefore Thus the number density of monopoles just after their production is smaller than, or at most equal to, the number density in a RD equivalent scenario. Furthermore, the factor e f > 1, and therefore the number density experiences more redshift due to the EMD phase than the RD equivalent number density, resulting in a smaller frozen abundance.
For the other two cases, whether the monopole relic abundance is larger or smaller in the EMD scenario compared to the RD-equivalent scenario depends on the relative sizes of the critical exponents, and for case II, additionally on the ratio of the temperature of the hidden sector at reheating to the critical temperature. A sufficient condition for the right-side of (4.17) to be less than or equal to one is 2ν ≤ 1 + µ . This condition can be verified by considering the relative sizes of the numerical factors involved: • For case II, note that 1 2 so that this fraction of f 's is bracketed by e f . Thus for critical exponents satisfying (4.19), the factor of e 3/4 f in the denominator of (4.17) due to the redshift is always larger than the ratio of number densities at the production time, irrespective of the relative size of (T . But for critical exponents violating (4.19), then the right-side of (4.17) can in principal be larger than 1, but whether that occurs depends on the relative of size of e f and the ratio of temperatures (T (hid) • In the last case, of monopole production after EMD, the ratio of f 's is the same as for case II, because f . To further simplify the analysis, assume that the visible sector degrees of freedom are the same in the two scenarios when the phase transition occurs in the hidden sector (which may occur at different visible sector temperatures). Then if the critical exponents satisfy (4.19) the ratio on the right-side of (4.17) is always less than one.
We therefore conclude that provided the critical exponents satisfy 2ν ≤ 1 + µ, the current frozen monopole abundance in a scenario involving EMD is always less than or equal to that in a pure RD equivalent, for a fixed monopole mass. This, along with the mass-dependence of (4.13), results in a larger monopole mass needed to account for a fixed Ω M h 2 when EMD is involved.

EMD by a modulus: numerical results
We now move to consider specific mechanisms for establishing a period of EMD, beginning with the case where the matter-dominating field Φ is a scalar modulus with mass m Φ and initial amplitude Φ i M P [20]. The modulus begins to oscillate, acquiring a matter equation of state, when H ≈ m Φ , at which time its energy density is given by ρ Φ (t i ) = (1/2)m 2 Φ Φ 2 i . This initial energy density, along with the matter-like redshift relation ρ Φ ∼ a −3 , determines how quickly Φ can dominate over the background radiation energy density, be it of the hidden or visible sectors. The initial ratio of the VS radiation energy density to that of the hidden sector is given by the factor f i . The Hubble factor during the period before EMD by Φ is given by (3.1).
The modulus amplitude, initially fixed at Φ i , starts to oscillate once H m Φ , and an EMD phase begins shortly after the energy densities of Φ and radiation become comparable. Solving for H m Φ and redshifting to this first era of matter-radiation equality, one finds the expansion at this time approximately corresponds to In calculating this, we have assumed the energy density of Φ is dominant over, as opposed to equal to, that of radiation, which results in a better agreement between our analytical calculations and numerical results shown below. For a modulus with maximal amplitude, we note that the modulus essentially dominates the energy density of the Universe as it begins to oscillate, while a smaller amplitude results in a delay. In order to successfully establish EMD, Φ must also be sufficiently long lived such that its decay completes well after the start of EMD. The minimum value of the initial amplitude, corresponding to decay at the onset of EMD, can be estimated from (3.2) and (5.1) to be For tree-level decays, a given visible sector reheat temperature determines not only the end of EMD, but also the mass of Φ and thus the minimum amplitude to have an EMD era at all. A choice of Φ i , within the allowed limits, then determines how early the EMD phase starts. We parenthetically note that for a given visible sector reheat temperature, the inclusion of a loop factor in Γ Φ shifts the values of m Φ and Φ i which correspond to a particular EMD duration. There is however, some degeneracy in the corresponding cosmologies. For instance, a change in initial amplitude of 10 −1 can be compensated by a change in mass of 10 4 and a loop factor α of 10 −6 , such that the resulting EMD phase is unchanged, having the same H MD , Γ Φ , and boundary condition (5.2). As mentioned previously, we will set α = 1 throughout unless otherwise specified.
The evolution of the three background energy density components (that of Φ and the radiation from the hidden and visible sectors) is governed by the following usual set of Boltzmann equations: . We emphasize that, for simplicity, in the Boltzmann equations above we have taken Φ to decay only to the visible sector, though it is straightforward to include branching fractions for decay to both sectors. We numerically solve this set of equations beginning in a period of RD by any combination of visible sector and HS radiation, and track the evolution sufficiently beyond reheating such that RD in the visible sector is well-established. In our numerical calculations, we use a smooth function to estimate the temperature dependence of the relativistic degrees of freedom for energy density in the VS, g (vis) * , shown in Figure  1. At temperatures greater than ∼100 GeV, when all SM species are relativistic, g (vis) * takes its maximum value of 106.75. As the temperature decreases, the value smoothly drops as the various particle species become nonrelativistic. We only show temperatures greater than 1 GeV because the VS reheat temperature in our scenarios is typically larger. The minimum value of g (vis) * , corresponding to the present era, is 3.36 assuming 3 massless neutrino species. For the HS we assume a constant g (hid) * = 100. Figure 2 shows the energy density evolution in the two cases of initial RD by the HS (f i << 1) and VS (f i >> 1) respectively, for an example set of parameters.
We allow the phase transition of the HS to occur at any time in the background evolution, and obtain the resultant current monopole abundance from the numerical solution. This is done by evaluating (2.10), the equation for the correlation length at the phase transition, when the temperature of the hidden sector reaches T (hid) C , and then approximating the number density of monopoles at that time as n M ∼ ξ(t * ) −3 . Subsequently, the number density is simply redshifted numerically through the EMD era and then normalized to the VS entropy density at reheating. We now turn to our numerical results. In Figure 3 we plot the present-day relic monopole abundance, Ω M h 2 , as a function of monopole mass, m M , where we have taken x M = 50 to be fixed, as well as λ = 1. In what follows we will set x M = 50 and λ = 1 throughout unless otherwise noted. The other parameter values match those of Figure 2. We show both numerical results, obtained from numerically solving the Boltzmann equations, and the three analytical approximations of Section 4, (4.10), (4.11), and (4.12). The numerical curve, shown in dark blue, has three distinct segments corresponding to the three regimes of production time: in the top right, monopoles are produced in the RD period before EMD -the slope of the curve in this region is the same as that of a pure RD monopole production scenario; the central segment of the curve corresponds to production during EMD, with a slope given by (4.11); and in the bottom left section, production after EMD recovers the RD slope. As can be seen by inspection, the analytic approximations, (4.10), (4.11), and (4.12), have extremely good agreement with the numerical results -the analytic results correspond to the light-blue dotted line "lying inside" the numerical curve. Figure 3 also shows colored regions depicting the three regimes of monopole production time. A given parameter set {m M , T (vis) RH , Φ i , f i , α, x M , λ, µ, ν} corresponds to a single point on Figure  3, so that as m M is varied, a single (blue) curve is traced out, passing through the colored regions that correspond to production after, during, or before the time of the phase transition. In this way only a subset of the colored regions are accessed. However, other points in the colored regions can be accessed by varying m M together with one or more of these other parameters. This behavior can be seen in Figure 4, which we discuss in more detail below. Figure 3 also shows as black dashed lines the two analytical expressions for production at the beginning (4.6) and end (4.7) of EMD, separating these three regimes. One way to interpret the boundary curves is the following. These two lines give analytic predictions for monopole production if, for a given monopole mass, production occurs at the end of initial RD and start of EMD (upper), or end of EMD and start of second RD (lower). The intersection of either of these dashed lines and the solid blue (numerical) line gives the mass for which production did occur at cross-over, for the parameters assumed for the solid line. These intersection points therefore mark the transitions between the three behaviors of the numerical line discussed in the previous paragraphs.
Lastly, we note that the entire numerical curve sits at higher monopole masses when compared to a pure RD production scenario (shown by the red dashed line) because of the offset of the hidden and visible sector energy densities. This is consistent with the behavior of (4.13) and (4.17), specifically that the right-side of (4.17) is always less than one when 2ν ≤ 1 + µ.
In Figure 4 we show how the curves of Figure 3 change for a variety of parameter values. As the beginning of EMD is placed earlier (by increasing the initial modulus amplitude Φ i ) while keeping the VS reheat temperature T (vis) RH fixed, the numerical curves (along with their analytical counterparts) shift farther away from the RD line toward larger monopole masses due to the increased amount of dilution from a progressively longer EMD period. If instead the end time of EMD is placed later (by decreasing T (vis) RH ) while holding the start time fixed, the curves again shift toward higher monopole masses due to the longer EMD period, but now the corresponding dashed boundary lines shift downward due to their dependence on the reheat temperature. Finally, as the critical exponents, µ and ν, are varied, the slopes of the curves change as expected. The upper-right shaded region (orange) corresponds to monopole production having occurred during the initial RD phase prior to EMD; the large central/lower-right region (magenta) corresponds to production during the EMD phase; and the small lower-left region (green) corresponds to production in the RD epoch after EMD has ended. Where the blue lines overlap with these three regions specifies the period in which monopole production occurred. For reference across the two panels, the dotted horizontal and vertical lines in both panels mark Ω M h 2 = 0.12 and m M = 1 PeV respectively. The entire set of curves and region boundaries in the right panel is shifted downward and to the left relative to the left panel, along the RD equivalent line due to the larger final offset between the visible and hidden radiation energy densities after reheating (see Figure 2).
In all panels of Figure 4, all of the numerical curves retain the three-region slope behavior displayed in Figure 3, with the regions separated by the two dashed boundary lines regardless of the specific parameter values, as expected. We note that the change in slope between the three regimes of production time is most noticeable in the bottom blue curve of the bottom two panels, for which µ = ν = 1. As in Figure 3, the left panels correspond to initial RD by HS radiation (with f i < 1), while the right panels correspond to initial VS domination (f i > 1). The full set of lines shown in each right panel is shifted downward and to the left as f i is increased above 1 relative to the corresponding left panels. Otherwise, the scale and orientation is the same between the left and right panels.

EMD by a decoupled particle: numerical results
Rather than being a modulus, the field Φ that drives EMD can instead be a heavy particle which decouples from either the hidden or visible sector at a very early time and subsequently dominates the energy density of the Universe as a non-relativistic matter component before eventually decaying (see Figure 5). We will parameterize the interaction rate of Φ with the sector from which it is decoupling (the "host" sector) by the thermally averaged annihilation cross-section times relative velocity, σ Φ v . 9 The Boltzmann equation for the number density of Φ is then where Γ Φ is the decay rate given in (3.2), and the Hubble parameter H is again given by the sum of all energy density components. In our numerical calculations, we use the integral expression for the equilibrium number density, where + is for fermions, − is for bosons, E(p) 2 = m 2 Φ + |p| 2 , g Φ is the number of internal degrees of freedom for Φ, and the temperature T is of the host sector. If Φ decouples from the HS, the remaining two Boltzmann equations for the radiation components are The energy density of Φ is given by ρ Φ = E Φ f Φ , which we have approximated as E Φ n Φ , Figure 6. Numerical evolution of the background energy density components with scale factor in the case of EMD by a decoupled particle Φ. EMD begins once ρ Φ dominates over both radiation components, and lasts until Φ decays. Left panels: initial RD by the hidden sector. Right panels: initial RD by the visible sector. Top panels: Φ decoupling from the dominant sector. Bottom panels: Φ decoupling from the subdominant sector. The values of σ Φ v in each panel are chosen to correspond to relativistic freeze-out, thus yielding the longest possible EMD phase for the chosen background parameters. Larger values of σ Φ v will result in nonrelativistic freeze-out of Φ while smaller values lead to freeze-in, both of which reduce the duration of EMD by lowering the frozen Φ abundance and hence delaying the start time. Note that in the bottom two panels, relativistic freeze-out of Φ essentially results in the limiting EMD case where the start and end are nearly coincident. The mass of Φ in all panels is m Φ ≈ 10 9 GeV, due primarily to the value of T with the average energy per particle given approximately as E Φ ≈ m 2 Φ + 9T 2 [35,40]. The temperature T is that of the host sector. Note that we retain the decay of Φ predominantly to the visible sector in order to preserve the standard history from BBN onward. 10 We numerically solve the Boltzmann equations, in both decoupling cases, for the background energy densities, as shown in Figure 6. As before, we use a smooth function for the temperature dependence of the relativistic degrees of freedom in the VS, g (vis) * , shown in Figure 1. To obtain the energy density evolution, we start in RD at some initial early time, with the HS and VS radiation related by the factor f i , and with negligible Φ energy density. 11 As the Universe cools, Φ decouples from its host sector via freeze-out or freeze-in, leaving a frozen energy density that redshifts like matter once Φ becomes non-relativistic. This matter energy density can then dominate over radiation, provided that the frozen energy density is high enough for domination to occur before the eventual decay of Φ. The decay completes near H ≈ Γ Φ , and we are subsequently left with the standard phase of domination by visible sector radiation.
The evolution of the equilibrium number density for Φ transitions from relativistic to nonrelativistic when the temperature of the host sector drops below m Φ . Because of this transition, there is a maximum frozen number density for a given m Φ , which is achieved through the decoupling of Φ while it is relativistic and in chemical equilibrium with its host sector. This is relativistic freeze-out. If Φ were to start with a number density larger than equilibrium, annihilations would drive it down to the equilibrium density, unless the annihilation rate was too small, which is not a scenario we will consider here because we assume RD at the initial time in order to justify an origin for the intervening EMD phase. Decoupling through relativistic freeze-out results in the earliest possible start time for the EMD phase caused by Φ of a given mass, and requires the annihilation rate to be large enough such that Φ reaches equilibrium while still relativistic, but not too large such that it remains in equilibrium after becoming non-relativistic. The largest value of σ Φ v that corresponds to relativistic freeze-out (which is the transition between relativistic and non-relativistic freeze-out) can be approximated by for relativistic decoupling, and where ζ(s) is the Riemann zeta function of s. 12 If instead the annihilation rate of Φ is large enough to maintain equilibrium with its host sector below T ≈ m Φ , then decoupling will occur via non-relativistic freeze-out, resulting in a smaller frozen number density and thus a later start time for EMD. As the annihilation rate increases further, the frozen Φ energy density decreases and the start of EMD approaches the time of reheating, resulting in a shorter duration for the EMD phase. This gives an upper limit, corresponding to H MD Γ Φ , on the value of σ Φ v , for a given mass and decay rate (or equivalently visible sector reheat temperature) for EMD to happen at all: where H F is the expansion rate at freeze-out and given in Appendix C, and we have used (C.2) for the expansion rate H MD at the time of matter domination. Now going in the other direction, if the annihilation rate is smaller than that needed for relativistic freeze-out, Φ will never reach local chemical and thermal equilibrium, which may possibly lead to a freeze-in process [41]. If freeze-in does occur, lowering σ Φ v further reduces the out-of-equilibrium number density, and thus the duration of EMD, down to a minimum value corresponding to the absence of EMD altogether. The value of σ Φ v corresponding to the transition between freeze-in and relativistic freeze-out (which defines the lower limit of the range of values leading to relativistic freeze-out) is approximately 9) and the minimum value corresponding to H MD Γ Φ is (see Appendix D) (6.10) We summarize these three different regimes of the annihilation rate. Starting with small annihilation rates, the decoupling of Φ proceeds as follows. For σ Φ v less than the right-side of (6.10), Φ decouples via freeze-in at such low energy densities that it will never dominate over radiation before decaying. For rates that satisfy (6.10) but are less than (6.9), the frozen-in energy density of Φ is large enough to dominate, leading to longer EMD durations as σ Φ v , and thus the frozen-in energy density, is increased. Between (6.9) and (6.7), decoupling occurs via relativistic freeze-out, which yields the largest frozen Φ energy density and the longest possible EMD duration, indepen-dent of σ Φ v . We note that essentially the only difference in (6.9) and (6.7) is the presence of the initial host sector temperature or the Φ mass in the denominator. Because the initial temperature can in general be quite large compared to m Φ , the regime of σ Φ v corresponding to relativistic freeze-out can extend for many orders of magnitude. For σ Φ v larger than (6.7) but satisfying (6.8), Φ decouples via nonrelativistic freeze-out, resulting in smaller frozen-out energy densities, and thus shorter EMD durations, as σ Φ v is increased. Finally, for rates larger than the right-side of (6.8), the frozen-out energy density is again too small to establish EMD before Φ decays.
Other than defining the range of annihilation rates that can yield an EMD phase 13 , the significance of these regimes of σ Φ v is that a particular EMD phase, with a fixed start time and end time, can be established by two different values of σ Φ v , one corresponding to freeze-out and the other to freeze-in.
The abundance of monopoles produced by the HS phase transition is determined by using (4.10)-(4.12), which are given in Section 4. These expressions were obtained in a modelindependent context and are valid in the cases presented in this section, provided that we use the appropriate expressions for quantities such as H MD .
The present-day relic monopole abundance is shown in Figure 7 as a function of monopole mass for some example parameter values, and we have again taken x M ≡ m M /T (hid) C = 50 and α = λ = 1. We in particular consider several values for σ Φ v , and we have checked that these values are well-below the perturbativity limit for the Φ mass inferred from (3.2), (3.3), and the assumed reheat temperature. As in the modulus case, there are three regions corresponding to monopole production before, during, and after EMD, and the curves have the same behavior as before. The main feature that sets the decoupled-particle case apart from the modulus case is that any particular curve can be obtained be either non-relativistic freeze-out or freeze-in, meaning the value of the annihilation rate of Φ can be quite different while still reproducing the same curve. Otherwise, the same regions are generally accessible to a modulus or decoupled-particle scenario, where the maximum extent toward larger monopole masses is set by either the maximum initial modulus amplitude or by relativistic freeze-out in the two cases respectively.
We finally note that the case of freeze-in depends on the initial host-sector temperature because freeze-in of Φ occurs in RD, such that the time of peak Φ production from the background occurs at the initial time (see [37] for details of freeze-in during RD before EMD). In our numerical calculations, we chose the initial time arbitrarily, with an initial energy density configuration consisting of dominant radiation and negligible Φ. For a given initial time, there is a unique annihilation rate that results in a particular freeze-in Φ energy density, provided that we remain within the freeze-in regime of the annihilation rate. The important thing to note is that the accessible region in Ω M h 2 vs m M is generally independent of the initial time because it is determined by the start and end of EMD, which can be obtained by multiple values of the initial time and annihilation rate. 13 We include an additional constraint in Appendix E on the parameter values that must hold for an EMD phase to have nonzero duration. Figure 7. Dependence of the present-day monopole relic abundance on the monopole mass in the case of EMD driven by a decoupled particle. As in Figures 3 and 4, the solid curves (purple and green) are obtained from a numerical evolution of the background, while the dotted lines (light purple and light green) on top of the numerical curves are the analytical expressions (4.10), (4.11), and (4.12). The purple color denotes Φ decoupling from the HS, while the green color corresponds to decoupling from the VS. All other lines have the same meaning as in Figures 3 and 4, which we repeat here. The red dashed line in all panels marks the purely RD equivalent scenario. The two black dashed lines in all panels indicate monopole production occurring at the start or end of EMD. The dotted horizontal and vertical lines in all panels mark Ω M h 2 = 0.12 and m M = 1 PeV respectively. Left panels: initial RD occurring in the hidden sector. Right panels: initial RD occurring in the visible sector. Top panels: Φ decoupling from the dominant sector. Bottom panels: Φ decoupling from the subdominant sector. In each panel, the curves which sit farthest to the right correspond to relativistic freeze-out of Φ from its host sector and thus mark the largest monopole masses accessible for the chosen parameters. The dependence on T (vis) RH and the critical exponents µ and ν is the same as in Figure 4.

Parameter values giving observed dark matter relic abundance
In this section, we will consider the values of our various parameters that result in the observed present-day DM abundance of Ω M h 2 = 0.12. As we've seen in the two previous sections, our analytical and numerical results agree very well, and we will therefore present an analytical analysis of the main parameters of our scenario, rather than a full numerical parameter scan.
We will primarily use (4.10)-(4.12) as well as (B.6) which gives f RH ∝ e f , requiring that the observed DM relic abundance is achieved. For clarity in the analysis below, we will not specify the identity of the field Φ, taking the beginning and end of EMD as the more fundamental parameters. We will use the VS reheat temperature T (vis) RH to set the end of EMD, and the factor e f = a RH /a MD to fix the duration of EMD.
Recall that e f can be expressed as (see Appendix B): The remaining parameters are the initial ratio of the VS to HS radiation energy density f i , the monopole mass m M , as well as the various parameters associated with the details of the phase transition, x M , λ, µ, and ν. Four of these eight parameters can vary by many orders of magnitude in the cosmological histories we have been considering: m M , T RH , e f , and f i , so here we will focus on those as they lead to a more direct effect on the resulting cosmology. The others have much narrower ranges, and for these we will consider a discrete set of possibilities. Also, we will not vary parameters such as α, m Φ , Φ i (in the case of the modulus), or σ Φ v (in the case of the decoupled particle), as including variations in these parameters is degenerate, in the sense that they lead to the same cosmology, as discussed in Section 5. Figure 8 shows contours of T (vis) RH in the m M −e f plane, with the monopole abundance held fixed at Ω M h 2 = 0.12. The region above each contour results in overproduction of DM, while the region below results in underproduction. What can immediately be seen from the figure is that most lines shown have positive slopes in this plane, meaning that a longer EMD duration (i.e. a larger value of e f ) requires a larger monopole mass in order to achieve the same monopole abundance. This is consistent with the behavior in Figures 4 and 7, where the curves corresponding to longer EMD periods cross the Ω M h 2 = 0.12 line at larger monopole masses. Furthermore, for fixed monopole mass, a longer EMD duration results in too much dilution and thus underproduction of DM, while a shorter duration doesn't dilute the monopole abundance enough, leading to overproduction.
In each panel of the figure, the region accessible to the T RH indicate that EMD ends at an earlier time, while larger values of e f correspond to longer EMD durations. Each solid contour has two segments with different slopes (a few of which occur beyond the range shown in the figure). Contours above the dashed contour overlap in their steeper segments, which follow the upper dotted black boundary line corresponding to monopole production before EMD, while those below overlap in their shallower segments, and follow the lower dotted black boundary line corresponding to production after EMD. Segments that are parallel to the dashed contour indicate monopole production during EMD. The region above a given contour results in overproduction of DM, while the region below results in underproduction. The slight differences in the overlap of the upper contours are due to changes in g (vis) * RH . We include a horizontal line at m M = 100 PeV for reference across panels, as well as a red 'star' which marks the pure RD scenario at e f = 1. The green circles located at e f ≈ 1.2 × 10 9 along the EeV contour, and e f ≈ 8.7 × 10 11 along the 7.2 PeV contour, correspond to the bound on H MD from (7.4). Please see the text for more details.
Each panel additionally shows a special, dashed blue-green contour which separates two regimes of T (vis) RH , and passes through the RD point mentioned above without changing slope. Relative to Figures 4 and 7, this contour corresponds to the special value of T (vis) RH which places the intersection of the two black dashed lines (representing the start and end of EMD) at Ω M h 2 = 0.12 (this is most easily seen in the middle panels of Figure 4, where the intersection point of the two black dashed lines shifts along the RD line as T (vis) RH is changed). As the duration of EMD is increased along this dashed contour, the contour rises away from e f = 1 with a slope given by RH ) has two segments with different slopes: beginning on the left side at e f 1, the contours rise along the upper boundary line, corresponding to monopole production before EMD, until they reach a point which corresponds to production at the start of EMD -beyond this point, the contours deviate from the upper boundary with a slope parallel to the dashed contour -this segment corresponds to monopole production during EMD.
The contours located below the dashed contour (with higher values of T (vis) RH ) have a similar two-segment behavior: beginning again at e f 1, the contours rise at a shallow slope along the lower boundary line (monopole production after EMD), until they reach a point corresponding to production at the end of EMD -from here on the contours leave the lower boundary and continue with the same slope as the dashed contour -production in this region occurs during EMD. The region above the dashed contour can therefore only access monopole production before and during EMD, while the region below only accesses production during and after EMD. Additionally, we note that in the lower right panel, with µ = ν = 1, the slope of the "after EMD" segment is essentially independent of e f , consistent with the lower panels of Figure 4 where the segments of the numerical curves corresponding to monopole production after EMD coincide with the pure RD scenario, thus erasing any dependence on the prior EMD history.
The boundaries of the accessible region in the m M − e f plane, which correspond to monopole production before and after EMD, are given by (4.10) and (4.12), and are independent of T (vis) RH . This can be trivially understood for production after EMD, while in the case of production before, the monopole abundance experiences dilution from the full EMD phase, regardless of it's specific timing. However, as e f increases, a given contour turns away from the boundary at a point that corresponds to the start (upper contours) or the end (lower contours) of EMD, which does depend on T where we have additionally made use of f RH 1. This expression can then be used along with (4.10) to locate the monopole mass and EMD duration which result in the observed DM abundance for monopole production at the start of EMD.
For monopole production at the end of EMD, (3.3) can be expressed in terms of VS quantities and set equal to itself in terms of HS quantities to obtain This expression, together with (4.12), then yields the monopole mass and EMD duration which result in the observed DM abundance for monopole production at the end of EMD. The value of T (vis) RH for the dashed contour shown in Figure 8, which separates the two sets of contours, can similarly be obtained by first eliminating e f in (7.2) and (7.3). This corresponds to the RD point at e f 1, marked by the red star, where the before and after boundaries meet (as do production at the start and end of EMD). Then, using either (4.10) or (4.12) gives the monopole mass required for Ω M h 2 = 0.12, which then yields the special value of T  (4.13). Note that because the monopole abundance produced during the EMD and RD periods displays different power-law dependence on the monopole mass, the EMD and RD segments of the curves will shift by different amounts, resulting in movement of the turn-off points along the before and after boundaries.
As is evident from (7.1), a long duration for the EMD phase requires a large separation between H MD and Γ Φ . This is rather easy to achieve, even for high reheat temperatures. However in inflationary models, the Hubble parameter at the start of EMD, H MD , is bounded from above by the value of the Hubble parameter H I at the end of inflation. This would correspond to an interesting scenario in which after inflation the early Universe directly enters the EMD phase, with some reheating in the hidden sector so that the initial temperature in that sector is above the critical temperature. However, from the non-detection of tensor modes, PLANCK data gives an upper limit to H I [42] H I < 2.5 × 10 −5 M P . Lastly, we comment on some interesting effects when the critical exponents, µ and ν, satisfy µ = ν > 1. Though we have specifically considered µ = ν = {1/2, 1} in our figures, the expressions presented throughout the text are applicable to more general values of the critical exponents. 15 In particular, we recall that as µ and ν approach 1, the case of monopole production after EMD (case III) approaches a purely RD scenario, so that when µ = ν = 1, the dependence on the prior EMD history is completely removed. This suggests that for µ = ν > 1, or more generally 2ν > 1 + µ, the monopole mass required for the observed dark matter abundance can actually be smaller than the RD case, at least for monopole production after EMD or shortly before its end. We have checked that this is indeed the case, however, the RD curve itself gets shifted to higher monopole masses when µ = ν > 1 such that case III actually results in heavier masses as compared to µ = ν < 1 (keeping the relic abundance fixed). This can be seen from expressions such as (4.5) and (4.12), where increasing the critical exponents above 1 results in an increase in the required monopole mass for both EMD and RD scenarios, but the increase is larger in the RD case. We also note that, from (2.10), the correlation length gets larger as the critical exponents are increased, resulting in less correlation volumes per Hubble volume, which in turn results in a smaller monopole number density at production.
In this work we have broadened the scale for hidden sector monopoles masses to O(1-10 5 ) PeV. One may wonder how robust the lower limit of 1 PeV actually is. The effect of lowering the monopole mass relative to a RD scenario when µ = ν > 1 is greater for a longer EMD duration, as the lower boundary line in Figure 8 acquires a negative slope. Additionally, the visible-sector reheat temperature needs to be larger than that for µ = ν < 1 in order for contours of the observed dark matter relic abundance to access the lower boundary line -note the different positioning of the PeV contour in the upper-left and lower-right panels of Figure 8. Because of these two effects, an extended EMD period occurring very early will have the greatest effect in producing enough lower-mass monopoles to reproduce the observed DM abundance. Perhaps if the phase transition occurs toward the end of (or after) a period of EMD caused by inflationary reheating at very high temperatures, the monopole mass may be able to be brought below the PeV scale and still result in the full DM relic abundance. Furthermore, having the HS temperature be extremely suppressed below the VS actually helps lower the needed monopole mass significantly, as long as the VS reheat temperature is large enough to bring up the abundance. This suppression effect also applies to a purely RD scenario.
In passing, we finally note that like µ = ν = 1, setting µ = ν = 2 is another special case in which the monopole abundance produced during EMD is now independent of the Hubble rate at the time of production, and only depends on the critical temperature. This can easily be seen in (4.4), where the factor of H 2 C in the denominator due to redshift cancels the dependence on the critical exponents. If the altered phase of expansion is instead caused by a form of energy density other than matter, this effect would occur for a different value of the critical exponents.
Overall, with the exception of the effect of the critical exponents discussed above, as we vary the parameters of our scenarios, the accessible regions which reproduce the observed DM relic abundance do not change drastically. As we saw in Figure 4, the largest shifts occur when the critical exponents are changed. Our main finding that for hidden sector monopoles to be dark matter candidates, their masses must be larger than O(PeV) scale appears generic, with longer EMD periods leading to larger monopole masses when 2ν ≤ 1 + µ.

Discussion
In this work we have considered a scenario for dark matter production via a second order phase transition in the early Universe, where the dark matter (DM) candidate is a hidden-sector magnetic monopole. Such a topological dark matter scenario has been studied before, with the entire relic DM abundance being produced in the standard radiation-dominated (RD) era before BBN [4][5][6]12]. We have expanded the parameter space region of viability to allow the different sectors to have different temperatures, and by generalizing the cosmological history to include a period of early matter domination (EMD). By allowing the phase transition to occur at any time before, during, or after EMD, we have shown that histories involving EMD generally require heavier monopole masses in order to produce the entire DM relic abundance. Along with this general result, we have considered two specific examples of how a period of EMD may be generated: by a modulus, or by a heavy decoupled particle. These examples illustrate how one can embed our scenario in a specific model, and how the underlying model parameters influence the monopole abundance. Our main results are summarized in Figures 4, 7, and 8. We generally find that hidden sector monopoles in the mass range O(1-10 5 ) PeV can be dark matter candidates.
We now summarize our main caveats, address some ways our scenario can be changed for future work, and what we expect that will do.
Throughout this work we have assumed the number density for PeV scale monopoles is small enough to ignore the effects of monopole-anti-monopole annihilation, as shown in [12] following [31]. But because the scattering cross-section between fermions and monopoles is a strongly coupled problem, it is possible that the final monopole abundance is depleted more than the diffusion approximation studied in [31], due to the interaction with the hidden sector plasma (if present). The interaction of the monopole with the plasma may be more critical to understand if the monopole is a dyon, a possibility not considered here. Of course, if the number density decreases further due to annihilation, a higher monopole mass will be needed to get the same DM abundance.
Another key assumption pervading this work is that the second order phase transition is classical, although we have strayed from that strict assumption by allowing the critical exponents to have generic values. But a consequence of assuming the monopole to be a classical topological object is that the monopole mass and the temperature of the phase transition are at similar mass scales, m M ∼ T (hid) C . Our conclusions will change substantially in theories for which this relation no longer holds. A prominent counter-example is provided by the N = 2 Seiberg-Witten theory [32,33] near the massless monopole or massless dyon points of the moduli space, in which the effective theory below the symmetry breaking scale contains nearly massless composites -'mesons' and 'baryons' of a magnetic U (1). Additionally, here the effect of annihilations at energies near the scale of the transition are expected to be important.
Another fundamental assumption in our work is the set-up of our sectors, where we have assumed the sector which hosts the phase transition to interact very weakly, if at all with the visible sector of standard model particles. This can in general be different, and can result in changes to the monopole abundance after their production. For example, kinetic mixing between the visible and hidden sectors can lead to a long-range force which can then deplete the monopole abundance via annihilation. We expect this to have a similar effect to the scattering of monopoles with a HS plasma followed by annihilating, but in this case the monopole abundance can depend more strongly on visible-sector as well as hidden-sector properties. See for example [23].
Along these lines, we have also assumed that the energy density component driving the EMD period decays almost entirely to visible sector radiation. With additional interactions between the sectors, the EMD driving field may decay to hidden sector radiation as well. This can easily be incorporated into our analysis by generalizing the decay rate Γ Φ to include branching fractions to both visible and hidden radiation. One must then be careful to not produce too much hidden (or "dark") radiation by restricting the branching fractions with current limits on dark radiation [43].
Aside from the set-up of our sectors, another important generalization of our work is to allow for early domination by a component with a generic equation of state, rather than focusing on EMD alone. The redshift relation for the dominating energy density is then ρ ∝ a −3(1+w) , with the parameter w determining the behavior, which modifies subsequent calculations.
A specific alternative to EMD is a period of kination, where the kinetic energy of a scalar field dominates the energy density of the universe for a time. In such a period, the dominant form of energy density redshifts faster than radiation, with w = 1 and ρ ∝ a −6 , which can have interesting consequences for the monopole abundance if the phase transition occurs during or before such a period. In fact, the phase transition occurring after a period of kination can also affect the resultant monopole abundance, for example by flipping the radiation energy densities of the two sectors. Kination would typically not last very long because it dilutes as a 6 , but if other components are suppressed, it can last longer -perhaps the same EMD driving field can have an early period of kination which later transitions to EMD before decaying. One should track the behavior of radiation in the two sectors during such a history to see how it affects the temperatures and thus the final monopole abundance.
Lastly, in our decoupled particle example, the mechanism of Φ decoupling need not be velocity independent. This can lead to temperature dependence in the interaction rate of Φ with its host sector and can alter the details of the decoupling. Such effects, however, shouldn't change our main results, just the specifics of the particle decoupling models (what values of Φ mass and decoupling parameter lead to an EMD phase of a given start and end).
We hope this work stimulates further research into topological dark matter scenarios. and HS radiation energy densities at reheating is where we have redshifted hidden sector quantities back to the start of EMD, and where f i is defined as the ratio of the visible sector to hidden sector radiation energy densities at some time t i prior to the onset of the EMD phase, To facilitate our comparison between scenarios which include a phase of EMD and those which remain purely RD, we make use of the double ratio where in the second line we have made use of (B.6) and where we have included superscripts on the two f RH 's on the left-side for clarity (whenever f appears without a superscript label, it refers to the EMD case). We note that since in any given RD-equivalent scenario f (RD) RH is just a number, to simplify our notation we will often drop the subscript and just write this term as f (RD) . The energy density ratio in a purely RD scenario corresponding to an EMD scenario with initial domination by visible sector radiation is given by f (RD) RH = f i , while in the case of an EMD scenario with initial domination by HS radiation, it is f (RD) RH = 1. We have additionally numerically verified the value of e f as the ratio of the scale factors at the end and beginning of the EMD period, as well as the double ratio of radiation energy densities.

C Decoupling of Φ from either sector via freeze-out
In order to analytically estimate the relic abundance of topological DM from (4.10)-(4.12), we need to obtain an expression for the Hubble rate at the onset of EMD, H MD . We do so by redshifting the frozen number density of Φ at the time of freeze-out, given by n Φ,F , to the start of EMD: What remains is to specify H F , which we do below for a number of cases.
C.1 Non-relativistic freeze-out from hidden sector Using the usual freeze-out condition of n Φ,eq σ Φ v = H F , with the non-relativistic form of the equilibrium number density for a boson Φ, we have where we have used H F ≈ (1 + f i ) . Rearranging yields an expression that can be solved for x F : If Φ is instead a fermion, the left-side of (C.3) is multiplied by a factor of 3/4, with a corresponding change in the expression for x F . The solution to this can then be used in the expression for H F above to complete its specification in terms of the parameters of our scenario.

C.2 Non-relativistic freeze-out from visible sector
Here we define x F ≡ m Φ /T (vis) F , resulting in (C.5) and x F ≈ ln Otherwise, this case is the same as above.

C.3 Relativistic freeze-out from hidden sector
In this case, we use the relativistic expression for the equilibrium number density, giving (1 + f i ) 1/2 . (C.8)

C.4 Relativistic freeze-out from visible sector
In this case, we have

D Decoupling of Φ from either sector via freeze-in
Because Φ is the source of the EMD period, at some point it decouples in the prior RD phase. If the annihilation rate to produce Φ is too tiny, Φ may never reach local, chemical and thermal equilibrium with the ambient radiation. However, the produced number density of Φ particles may be large enough to eventually dominate the energy density. This is known as freeze-in [41]. In this case, freeze-in in a RD period is dominated by the relativistic component and the abundance is set at the initial time. We begin with d(a 3 n Φ ) dt = a 3 σ Φ v (n 2 Φ,eq − n 2 Φ ) − a 3 Γ Φ n Φ , (D.1) We are interested in the early evolution of the Φ number density in a freeze-in scenario well-before it decays, as well as well-before it reaches equilibrium. Thus we may drop the decay term relative to the decoupling term above, as well as the actual number density relative to the thermal equilibrium value. With these approximations we have which for a = a i (t/t i ) 1/2 and H = 1/2t, appropriate for RD, one has To continue, we must express the temperature dependence of the equilibrium number density in terms of H, which is most easily done by specializing to the two decoupling cases.

D.1 Freeze-in from hidden sector
If Φ is produced from the HS, we have (1 + f i ) 3/2 .

(D.5)
Assuming this can be large enough to dominate the energy density at, by definition, the beginning of EMD, and using m Φ n Φ,MD ≈ 3H 2 MD M 2 P , setting H F = H MD gives (D.6)

D.2 Freeze-in from visible sector
If Φ is produced from the visible sector, we similarly have (D. 8) In sum, the equations in this Appendix give the number density n F of Φ particles in a freeze-in scenario, assuming it is produced in the early Universe from either the hidden or visible sectors, evaluated well-before it decays. And by definition of the freeze-in scenario, the number density n F is assumed to be well-below its equilibrium number density.

E Additional consistency constraint for the decoupled Φ scenario
We obtain another constrain that must be satisfied in order for the EMD phase caused by the decoupled Φ to have nonzero duration. If Φ decouples from the subdominant sector, the value of f i must be such that the decoupled number density is large enough to lead to EMD. Using (6.7) for an annihilation rate that achieves relativistic freeze-out (which corresponds to the maximum frozen number density and thus longest possible duration for EMD), we require H MD Γ Φ . Using (C.2) for H MD and (C.8) and (C.10) for x F in their respective cases, we have in the case of decoupling from the HS while the VS is dominant, and in the case of decoupling from the VS while the HS is dominant.