Parametric Resonance Production of Ultralight Vector Dark Matter

Vector bosons heavier than $10^{-22}$ eV can be viable dark matter candidates with distinctive experimental signatures. Ultralight dark matter generally requires a non-thermal origin to achieve the observed density, while still behaving like a pressureless fluid at late times. We show that such a production mechanism naturally occurs for vectors whose mass originates from a dark Higgs. If the dark Higgs has a large field value after inflation, the energy in the Higgs field can be efficiently transferred to vectors through parametric resonance. Computing the resulting abundance and spectra requires careful treatment of the transverse and longitudinal components, whose dynamics are governed by distinct differential equations. We study these equations in detail and find that the mass of the vector may be as low as $10 ^{ - 18 }$ eV, while making up the dominant dark matter abundance. This opens up a wide mass range of vector dark matter as cosmologically viable, further motivating their experimental search.


I. INTRODUCTION
The existence of dark matter (DM) is one of the observational evidences for physics beyond the Standard Model (SM). Recently, vectors (X µ ) have gained significant attention as an intriguing dark matter candidates with unique experimental signatures. Theoretically, light vectors arise as gauge bosons of dark U(1)s, a simple extension of the SM and a common prediction of high energy theories. The origin of the vector mass is modeldependent and can either be a fundamental parameter in the full theory via the Stueckelberg mechanism, or can be generated through its coupling to an additional field which spontaneously breaks the corresponding U(1) via the Higgs mechanism. In either scenario the mass of the vector is stable under quantum corrections, motivating the possibility of vectors with ultralight masses, m X MeV, limited only to having wavelengths small enough to form galaxies, m X 10 −22 eV.
Experimentally, light relic vectors present different opportunities depending on their coupling to the SM. The overarching challenge in experimental prospects is competing with the powerful limits from stellar cooling [1,2] and fifth forces [3] while restricting considerations to the (approximately) conserved currents of the SM (otherwise one generically expects dominant constraints from flavor changing neutral currents [4,5]). Nevertheless there exist many experimental proposals to search for vector DM in unexplored parameter space. Such states can be observed through their coupling to electrically charged particles that could be searched for in resonant cavities [6], LC circuits [7], dish antennas [8], absorption in direct detection experiments [9][10][11], and low-energy threshold detectors [12][13][14][15][16]. If the vector couples to an unscreened force such as B − L then its coupling to neutral matter can be searched for in torsion balances and atom interferometry [17], gravitational wave detectors [17,18], and pulsar binary systems [19]. With current and proposed experiments, searches for vector dark matter can be undertaken over almost the entire mass range 10 −22 eV m X MeV.
While ideas to detect vector DM are plentiful, the theoretical prospects for producing ultralight vector DM are much less explored. The underlying challenge is the difficulty of producing such light states by matter-radiation equality while being consistent with the requirement of cold dark matter. For light vectors there are three classes of production which have been studied in the literature: freeze-in [20], misalignment [1], and inflationary fluctuations [21]. Freeze-in production is generically constrained by the bound on warm dark matter. Particles "frozen-in" through an interaction with the SM bath are produced with momentum ∼ T . Without additional dynamics, the momentum of the relics will redshift with the expansion of the universe and hence track the SM temperature limiting the produced DM mass to be above a keV to be consistent with the observation of small scale structure.
Misalignment has long been a standard non-thermal production mechanism for light bosons, first proposed for axions [22][23][24], and later considered for light vectors [1,25]. Here, a zero-momentum condensate of particles is produced as a result of the coherent oscillations of the field initially displaced from its minimum. Unfortunately, misalignment production of a massive vector fails. The energy density of a vector field redshifts as ρ X ∝ a −2 when H m X , and hence any proposed dark matter abundance will be diluted by a period of inflation. This makes the misalignment mechanism for vector DM highly inefficient unless an O(1) non-minimal coupling to gravity is added to make the vector conformally invariant and hence impervious to the expansion of the universe [1]. Such a special coupling quantum-mechanically destabilizes the mass of X, thus destroying one of the primary motivations for considering vector DM.
Alternatively, vector DM can be produced by quantum fluctuations during inflation [21]. This is an interesting possibility, as such a production has no dangerous large-scale isocurvature perturbations and appears unique to vectors. Here, the observed DM abundance is saturated for m X 10 −5 eV(10 14 GeV/H inf ) 4 , and thus observational constraints on the Hubble scale during inflation [26] limit the production to masses greater than about 10 −5 eV.
In this paper, we propose a new production mechanism for vector DM that occurs naturally if it obtains mass through a dark Higgs field. Generically, the production relies on a scalar field being displaced far from its minimum by the end of inflation. As the field rolls down its potential and begins to oscillate, its coupling to a vector results in a rapidly oscillating mass for the vector. This leads to non-perturbative production of X particles through a parametric resonance (PR) instability (as is the case in theories of reheating [27,28], though the dynamics we consider take place solely during radiation domination). Crucially, the rate of production is much greater than that of any possible perturbative process. The produced particles then have more time to red-shift, significantly relaxing the coldness constraint and allowing for the production of ultralight DM. This is in analogy with earlier work [29] on the PR production of axion DM via dynamics of a Peccei-Quinn symmetry breaking field (for other work on non-perturbative production of relics, see [30][31][32][33]). In this paper we focus on the minimal case where the scalar field is a dark Higgs. The nature of the resonance and resulting abundance of vectors and dark Higgses is different depending on the strengths of the gauge coupling e and dark Higgs quartic coupling, λ. We examine both limits and find that vector DM can be produced with masses as light as m X 10 −18 eV, consistent with all constraints. This opens up most of the mass range for vector DM as cosmologically viable, and further motivates the experimental program searching for such particles.
The paper is organized as follows. In Sec. II we outline the model of interest and show the limitations of a perturbative Higgs decay in producing light vector DM. In Sec. III we review the relevant non-perturbative dynamics of PR, specifically as it applies to vector production. In Sec. IV we examine the PR production of ultralight vector DM, and in Sec. V we discuss additional cosmological consequences and constraints on the mechanism. Finally, in Sec. VI we conclude and discuss future directions.

II. THE MODEL
We now present an outline of the model and detail the dynamics of an oscillating scalar field in the potential. As our starting point we consider a complex scalar field, ϕ, that will give a mass to the vector: where D µ = ∂ µ + ieX µ , e is the dark gauge coupling constant (we absorb the scalar charge into the definition of e), and X µν is the field strength tensor. We consider the simplest model of spontaneous symmetry breaking with a potential parameterized as Expanding ϕ around the vacuum expectation value (VEV), we obtain: where The vacuum masses of X and dark Higgs boson are m X = ev and m φ = √ 2λv, respectively. 1 Furthermore, we refrain from making any assumptions about the magnitude of the vector coupling to the SM, up to assuming the coupling is not so large that it efficiently thermalizes the two sectors (or is phenomenologically excluded in other ways).
We assume φ starts out displaced from its minimum after inflation with an initial field value, φ 0 . The classical equation of motion for φ is which is valid as long as the back-reaction due to any created particles is negligible (these effects are crucial in the termination of non-perturbative particle production and will be addressed later). The field is stuck until is the effective (field-dependent) mass, at which point φ begins oscillating about the minimum. As long as φ 0 M pl (regardless of the hierarchy between φ 0 and v) the universe is radiation-dominated at the onset of oscillations which begin at, where M pl = 2.4 × 10 18 GeV is the reduced Planck mass. We now consider the two limits for the initial field value, φ 0 v and φ 0 v. If φ 0 v, oscillations start at T osc and the solution is the well-known harmonic oscillations, φ(t) = Φ cos (m φ t). The amplitude of oscillations red-shifts with the scale factor a (we use the convention that a = 1 at the onset of oscillation) as Φ(t) = φ 0 a −3/2 , 1 Note that while the vector mass is radiatively stable, the scalar mass is not and naturalness would suggest a cut-off of order Λ v min{1, m φ /m X }. Ultimately we will be interested in VEVs much larger than the weak scale, so fine-tuning in the dark Higgs sector is not a serious constraint and we will not address it further. and the energy density in coherent oscillations acts as non-relativistic matter ρ φ ∝ a −3 .
Conversely, if the field value is large, φ 0 v, then the effective mass is m eff (φ 0 ) √ 3λφ 0 with T osc λφ 0 M pl . Due to the conformal invariance of the quartic potential, it is most convenient to switch to conformal coordinates. Furthermore, it is convenient to absorb the oscillation time into our definition suggesting the coordinate transformationφ ≡ aφ/φ 0 and dz ≡ λφ 0 dt/a. The equation of motion is then simply: where we use primes to denote derivatives with respect to z. The exact solution is an elliptic cosine function with elliptic modulus of 1/2, This function is usually well-approximated by the simple cosine function,φ cos(0.85z), the first term in its Lambert series expansion, but some features require keeping higher order terms and so we refrain from making this approximation. Here, the (original) field amplitude instead red-shifts as Φ(t) = φ 0 a −1 and the energy density in coherent oscillations acts like radiation ρ φ ∝ a −4 .
A. Perturbative Decay The dynamics of particle production depend critically on the initial field value, and we postpone a careful treatment to Section III. However, we generally expect nonperturbative effects are negligible if φ 0 v and we briefly review the physics in this limit. We first compute the production of vector DM from perturbative decay of the dark Higgs-this will eventually highlight the effectiveness of parametric resonance. Coherent oscillations of the φ field result in an yield of dark Higgs: where s is the entropy density. This population can decay into X if it is kinematically allowed, i.e., m φ > 2m X .
Since the co-moving number density in the dark sector is conserved, the dark Higgs condensate will fully convert into a co-moving number density of vectors Y X = 2Y φ . The timescale for this conversion is set by the decay rate Γ φ→XX , which is dominated by the decay into longitudinal modes of X: The underlying challenge with DM production via decays is that the X particles are initially highly boosted with momentum O(m φ ). In this case, the produced vectors begin red-shifting as non-relativistic matter when the universe is at a temperature As expected, the temperature at which vectors become non-relativistic increases with λ which corresponds to earlier decays. Observations of cosmological large-scale structure require that the DM be non-relativistic by around a keV and so we require T NR keV [34][35][36] (precise constraints range from ∼ 1−5 keV though suffer from astrophysical uncertainties). Based on (9) and (11), we find the vector abundance equals the relic density of DM for masses: where T eq 0.75 eV is roughly the temperature at matter-radiation equality. Note that production of light vector DM here favors large values of λ, which is ultimately limited by perturbativity λ < 2π (this is in fact a stronger condition than φ 0 < M pl ). Saturating the coldness and perturbativity constraints, we conclude that it is impossible to produce vector DM with mass less than a keV using perturbative decays of the scalar field.

III. PARAMETRIC RESONANCE
If φ 0 v, the production rate of vectors can be much larger than the perturbative rate. Such a large initial field value is a generic expectation unless the coupling with an inflaton strongly fixes φ to the origin. In the classical background of an oscillating φ field, the field X feels a large, oscillating, mass. This may lead to a period of non-perturbative, exponential production of vectors though parametric resonance (PR). 2 Particle production by PR is a well-studied phenomena, particularly in the context of reheating after inflation (so-called preheating) [27,28]. However, vector production by PR has not been studied nearly as extensively as for scalars. PR production of vectors has been previously considered for the transverse modes in [37] to explain the large observed magnetic fields in the galaxy [38] as well as in [39] to reheat the SM through the electroweak sector after inflation. In this section we review the theory of PR for a vector, and show that the dynamics in our case depend delicately on the hierarchy between couplings e and λ and require a careful treatment of the transverse and longitudinal modes. We present a differential equation which governs the production of longitudinal modes that is distinct from the well-studied Mathieu and Lamé equations (the typical differential equations studied in the context of PR). The classes of solutions are presented in an instability chart of the exponentially growing momenta modes as a function of e/λ, and we compare the PR production of longitudinal and transverse modes in the different limits of interest. Ultimately we show that, as a consequence of both initial conditions set by inflation and a longitudinal mode enhancement in the coupling, the longitudinal mode dominates production for a wide range of couplings.

A. Parametric Resonance for a Higgsed Vector
Using the conventional diag(1, −a 2 , −a 2 , −a 2 ) metric in an expanding universe, we can write the kinetic and mass term of X explicitly in terms of its temporal and spatial components: wherem ≡ m X (1 + φ/v). Since X t does not contain a kinetic term, it is an auxiliary field and can explicitly be integrated out using its equation of motion. Switching to k-space, we separate X into its longitudinal and transverse components such that k · X = kX L and k · X T = 0. As a result, the action for the vector field separates for the transverse and longitudinal components S = S T +S L : Note for here and throughout we use k to denote comoving momentum.
We now study the PR production of the transverse and longitudinal modes. Exponentially growing modes naturally occur in specific resonance bands, and it is conventional to map out the regions of unstable momenta as a function of the couplings. The dominantly produced modes lie in the widest resonance band and generally have a large exponential instability, resulting in a rapid conversion of the energy density in the oscillating scalar field into these modes. We present the instability charts for the transverse and longitudinal modes in Fig. 1 and refer to it throughout.
a. Transverse Modes It is convenient to switch to dimensionless, conformally invariant quantities: a time variable, dz ≡ λφ 0 dt/a as introduced in (7), a momenta κ ≡ k/λφ 0 , and a conformal oscillating mass, µ ≡ am/λφ 0 = (e/λ)cn(z). Doing so, the equation of motion for transverse modes becomes: This equation is known as the Lamé equation and has been extensively studied in the literature (see e.g. [40]). 3 Solutions to this equation are exponentially growing for certain momentum modes X T ∝ e µκz . The characteristic exponents, µ κ , are a non-trivial function of momenta as well as the ratio of couplings. PR is often classified as either broad or narrow, based on the width of resonance bands and the size of the characteristic exponents. For the mode equation of (17), the resonance is broad if e λ and narrow if λ e. We will be interested in both these limits, which have previously been solved analytically.
In the case of λ e the first resonance band around κ 2 1 dominates production while subsequent resonance bands (at larger κ) become increasingly narrow. Such narrow resonances are known to have suppressed production with a small range of produced momenta and characteristic exponents, ∆κ 2 , µ κ ∝ e 2 /λ 2 . Thus we conclude production of transverse modes is not efficient in this regime. For e λ the resonance is instead broad and can achieve much more efficient production. Inspection of Fig. 1 shows that the structure of the resonance, in particular the size of κ in the first resonance band, depends critically on the value of e/λ. Interestingly, it is still the case that for e λ there is an upper bound on the produced momenta which can be estimated analytically [40]. A necessary condition for exponential instability in the regime of broad PR is nonadiabatic change in the frequency of fluctuations. The (dimensionless) frequency felt by the transverse modes is ω(t) κ 2 + (e/λ) 2 cn 2 (z). If one defines an adiabatic parameter R ≡ |ω |/ω 2 , then for most of the oscillation period this is close to zero and the frequency changes adiabatically. The only time R > 1 is when the background field oscillates toward zero, cn(z) → 0 and κ 2 e/λ, which is an estimate of the upper bound on the dominantly produced momenta. In fact this is bound is evident in Fig. 1, where the widest resonance band (red) always lies below the line κ = e/λ. We thus find that the typical physical momenta produced by PR here is much less than the time-averaged mass of the vector, ∼ eφ 0 , and vectors are produced non-relativistically. We compute the maximum characteristic exponent numerically for e λ and find µ κ 0.2, in agreement with the previous literature [40].
b. Longitudinal Mode PR for the transverse modes reduce to equations that have been solved extensively in the literature. We now move to the longitudinal mode which, as we will show, dominates the production of vectors in a wide range of parameters. Starting from (16) and making the transformations to conformal fields, we find the equation of motion: To our knowledge the dynamics governed by this differential equation have not been studied in the literature.
Here we present a brief analysis, and leave a detailed study for future work. First we note that in the limit κ → 0, the equation of motion (18) reduces to precisely that of the transverse modes: This is expected, since at low energies the longitudinal mode can no longer be distinguished from the transverse modes and should obey the same dynamics. We also see this directly in Fig. 1, where the resonance bands of the two modes roughly coincide (except for very particular values of the couplings) in the limit of small momentum κ 1. The high energy limit is more challenging to analyze since the physics is obscured by a divergence in the friction term as the oscillating field passes through the origin. While it is in principle possible to solve the equation as is, it is simpler to introduce a field redefinition, Being a linear transformation, this does not mix the different momentum modes and hence does not obscure the structure of the resonance. The resulting equation of motion is: (21) If we then take the high-energy limits κ µ and µ 1, we recover a familiar form: This is analogous to the equation of motion for transverse modes (17), though crucially the amplitude of the oscillations is enhanced by a factor λ 2 /e 2 . As a result, the PR dominantly produces longitudinal modes with κ 2 1, which can also be seen directly in the inset of Fig. 1 where there is a wide instability band (blue) for the longitudinal mode in the limit e/λ → 0. This result in the high-energy limit can also be derived directly from the action of the dark Higgs ϕ using the Goldstone boson equivalence theorem. Expanding ϕ = (φ + v + iχ)/ √ 2 and switching to conformal fields, we find the same equation of motion for χ as found for the longitudinal mode in this limit (22).
We are now in a position to complete the discussion of PR for the longitudinal mode as a function of the couplings. Firstly we consider the limit of λ e (where we found the transverse modes are not efficiently produced). In this case the longitudinal mode is produced strictly in the high-energy regime κ µ, and the results follow the approximate form of the mode equation (22). Here we find the resonance is efficient for κ 2 ∆κ 2 1 and we again have a large characteristic exponent µ κ 0.1. We emphasize that, in contrast to the transverse modes, the longitudinal mode has a marginally narrow resonance allowing it to be produced efficiently. In addition, the longitudinal modes are produced highly boosted with relativistic momenta. We now turn to the limit e λ which is much more interesting. Although the mode equations become identical in the limit κ 1, PR production only occurs when the vector mass is rapidly varying (i.e. when adiabaticity is violated). At this point κ is of order the oscillating mass µ, and the longitudinal mode equation (18) does not approximately reduce to any well-known forms (due to the non-negligible friction term). Indeed, as is evident from Fig. 1, there are substantial differences between the resonance structures of the longitudinal and transverse modes in this regime. While the solutions for longitudinal mode similarly suggest an upper bound on the dominantly produced momenta, the bound may be larger than that of the transverse modes depending on the coupling. This is an intriguing feature that opens up the possibility of producing relatively boosted longitudinal modes, although we still expect that the momenta in the first resonance bands satisfy κ 2 < e/λ such that produced modes are not relativistic. Finally, we compute the typical characteristic exponent for longitudinal mode production in this limit to be µ κ 0.2.
c. φ Fluctuations In addition to vector production, an oscillating φ field will inevitably also resonantly pro-duce φ fluctuations with non-zero momentum (denoted as δφ to differentiate from the zero-momentum condensate which we continue to denote by φ) from the self-coupling, λ. We emphasize that these excitations are in addition to the zero-mode condensate that results from coherent oscillations and carry a particle interpretation similar to the vector fluctuations. The mode equation can be derived from (5) by restoring the momentum term and expanding the field as φ + δφ, keeping order linear terms in the fluctuations. The resulting equation of motion is identical to that of the transverse modes (17) but with the replacement e 2 → 3λ 2 : The PR is qualitatively similar to that of the longitudinal mode in the λ e case (22). Fluctuations of φ are dominantly produced at momenta κ 2 1 with a width ∆κ 2 1 and (slightly smaller) characteristic exponent µ κ 10 −2 .
d. Initial Conditions We have seen that due to parametric resonance, there is an exponential amplification of fluctuations in the fields X and φ for certain momentum modes. However, an important effect we have yet to address are the initial conditions for the fields. Assuming a period of inflation, we can estimate the initial conditions for each field. Transverse components of the vector are conformally invariant and do not experience the expansion suggesting that they should have an initial field value given by the Bunch-Davies vacuum with a power spectrum, P T (k) ∼ k 2 . The initial conditions of the longitudinal mode are more dramatic. These are created by coupling to the metric during inflation and can far exceed their transverse counterparts [21] with a power spectrum, P L (k) ∼ (H inf k/eφ 0 ) 2 (this applies for both λ e and e λ and assumes the vector mass during inflation is eφ 0 ). This gives a ratio of the longitudinal to transverse mode amplitudes at the end of inflation as, which is independent of k. Since we do not consider parameter space such that the vector mass is above the scale of inflation, the longitudinal mode will dominate the transverse mode production as long as they can both be produced efficiently. The scalar fluctuations during inflation behave similarly to the longitudinal mode and will have comparable initial conditions.

B. Final Relic Abundance and Momenta
The exponential production from PR does not last indefinitely. Thus far we have neglected the non-linear back reaction of these fluctuations on PR itself. There are three kinds of back reactions. 1) The vector and scalar fluctuations grow large enough and give large mass contributions to both φ and X that subsequently red-shift as ∝ 1/a, and can lead to other interesting cosmological effects that will be discussed later. Here, we see that a changing mass acts to shift the resonance bands and can thus ruin the important Bose enhancement in final states that leads to continued exponential production for growing modes. 2) Scattering of fluctuations with the zero-mode condensate as well as fluctuations shift the particle momenta out of resonance bands. Again, this destroys the Bose enhancement in produced fluctuations and can also shut down exponential production.
3) The scattering also depletes the zero-mode condensate and terminates PR. (A fourth effect, due to the expansion, is not present in this theory due to its conformal nature.) In practice these effects occur simultaneously and act to cease particle production when the energy density of the fluctuations becomes comparable to the original energy density in the condensate. While these effects are highly non-linear and challenging to compute, if particle production lasts long enough the condensate will completely convert into the produced φ and X particles, regardless of the detailed processes involved. We assume that the zero-momentum field is completely depleted and does not make up any of the dark matter today (we expect this is a reasonable approximation due to significant scattering with produced fluctuations at the end of PR). In this sense, a full solution to the equations of motion, including back reactions, gives us the relative fraction in these two populations. We can parameterize the yields after the conclusion of PR production as: .
Here ρ φ,osc = 1 4 λ 2 φ 4 0 , f is the relative fraction of the condensate co-moving energy density dumped into vectors, and E 2 i = k 2 * ,i + m 2 i are the co-moving energies of the particle species i (m i denotes the time-varying mass of the particle). For simplicity, we assume particles are produced with a co-moving momenta peaked at k * ,i though in practice there will be small corrections associated with an O(1) spread around this typical value.
Once PR stops being efficient the produced particles are in a highly non-equilibrium state with peaked momenta. These particles can still undergo collisions within the sector scattering their momenta and changing their number densities. This includes simple 2 → 2 elastic scattering as well inelastic processes inducing cannibalization. Tracking the dynamics rigorously throughout this "post-scattering" phase requires a dedicated computation putting vectors and scalars on lattices and is beyond the scope of this work. However, we can still qualitatively estimate the behavior in each limit.
For λ e, the longitudinal mode and the fluctuations only differ by their parity and hence have comparable energy densities after decay (f 1/2) as well as momenta k * ∼ λφ 0 . In this limit, the symmetry between the longitudinal mode and the fluctuations of φ allows us to treat them as a single fluid regardless of the details of the postscattering phase. Furthermore, we do not expect these processes to be active even if either species becomes nonrelativistic. We thus expect that the vectors and scalars should have comparable number densities and momenta at late times: For e λ, the situation is more subtle. As shown, vectors are primarily produced non-relativistically from PR (a detailed spectrum will depend on the coupling e/λ), while fluctuations of φ are produced mildly relativistically. At this point, elastic and inelastic processes are efficient in driving the sector toward a state equating the momenta and number densities of X and φ (f 1/2). Effective scattering during this time relies on a Bose enhancement of the final state which is spoiled at large enough momenta. We estimate that such processes cannot produce vectors with momenta larger than their mass eφ 0 , and as a result both species should eventually be upscattered to a co-moving momenta as large as k * ∼ eφ 0 . We thus expect yields of X and φ at late times of order: We have confirmed this expectation employing a lattice computation using LATTICEEASY [41] and approximating the vector interaction by that of a scalar field with a quartic coupling to φ.

IV. VECTOR DARK MATTER FROM PARAMETRIC RESONANCE
The above results apply for any Higgsed vector in the early universe, and we now consider the implications of PR for the production of ultralight vector DM. For the rest of this section, we assume a large initial field value φ 0 v. In practice, PR production is not instantaneous but requires sufficiently long exponential growth so we in fact have the condition φ 0 /v 10 − 100. 4 There are four (a priori) independent parameters in the model: {φ 0 , v, m X , m φ } or alternatively {φ 0 , v, e, λ}. As we have seen, the nature of PR depends on the relative strengths of the couplings in a non-trivial way. This is also true for the resulting constraints on safely obtaining the correct relic abundance of vector DM. Thus we look at the two simplifying limits separately: λ e and e λ. Here we focus on the fundamental challenge of being consistent with constraints on warm dark matter while producing the entire dark matter abundance. Additional constraints and phenomenological consequences of the vector production are examined in Section V, which we refer to in the results of Figures 2 and 3.
We begin with the case where the gauge coupling is small with respect to the quartic (and hence also m X m φ ). In Section III we estimated a yield for the X and φ fluctuations to be roughly equal at late times: In the absence of any additional interactions these yields are conserved until today. Once X and φ become nonrelativistic, X constitutes a small fraction of the energy density of the dark sector: Furthermore, the typical co-moving momenta of each species is of order k * ∼ λφ 0 . This is related to the physical momenta by red-shifting from the time of production. Crucially, particles are produced from PR at very early times near the start of oscillations. Due to the conformal invariance, we can effectively treat the yields (35) as being produced with a physical momenta λφ 0 at a temperature T osc even if the particles are dominantly created somewhat later (PR results in rapid, through not instantaneous, particle creation). Given the relative energy densities between the scalar and vector, it is most natural that the dark Higgs constitutes nearly all of the DM today, with the vector being a subdominant component. In order for this scenario to be consistent with observations we require the dark Higgs be both non-relativistic by around a keV and satisfy the relic density condition. Requiring the dark Higgs yield to be the right relic abundance fixes the required initial field value: The temperature at which the Higgs becomes nonrelativistic is given by If the fraction of vector DM is greater than a few percent, it must also be sufficiently cold; otherwise, the relic vectors will be a hot DM subcomponent which is ruled out by the cosmic microwave background [36]. We show the viable parameter space for a vector subcomponent of DM in Fig. 2 (left), fixing the fraction of vector DM to be 10 −2 (with the rest made up by the dark Higgs). Even if vectors make up a small fraction of the DM abundance, they may still be detectable. Nevertheless, it is interesting to consider the possibility that vectors make up all of the DM due to some dynamics which Viable parameter space for ultralight vector production using parametric resonance in the limit that λ e. We fix the initial condensate amplitude φ0 to obtain the correct DM relic abundance at every point in the figures. Shown are constraints from coldness (yellow), cosmic strings (red), isocurvature (purple), late time dark Higgs decays (brown) and sufficiently long PR (blue) as described in the text. Left: We do not incorporate any additional interactions, and the dark Higgs makes up the dominant source of DM. We fix the fraction of vector DM to be 10 −2 . Right: Vectors make up all of the DM, and the dark Higgs is eliminated at late times. We fix the scalar mass to be 10 GeV and show the corresponding constraint from thermalization requirements (green) as described in Appendix A . eliminated the dark Higgs yield at later times. In particular, it is possible to introduce additional couplings to the model that drastically alter the expected relic abundance of dark Higgses, without affecting the abundance of vectors produced from PR. (Although we might have naively suspected that the large initial yield of φ particles could simply decay away to vectors through the perturbative process (10), such population of vectors constitutes an O(1) hot DM component.) For now we will take it as a given that Y X reproduces the entire observed DM density. This fixes the initial field value: The X population becomes non-relativistic when the universe is at a temperature: As before, we require T NR keV. Note that the initial field amplitude has a maximum value consistent with the vector DM abundance and coldness constraints: This makes the condition on oscillation during radiationdomination (φ 0 < M pl ) trivially satisfied.
We now return to the elimination of the dark Higgs yield in the above scenario. If we assume the vector constitutes all of DM as per (32), then to avoid the dark Higgs dominating the energy density of the universe at an intermediate time the yield should have been destroyed by the temperature ∼ T eq m φ /m X . This is not a constraint, though a necessary condition for the above formula to hold as they assume radiation domination throughout. If, on the other hand, the universe has gone through a period of dark Higgs domination that later gets dumped into the SM this could have profound implications on small scale structure [42][43][44] and changes the predicted relic abundances. If we assume this matterdominated era lasts until the dark Higgs reheats the universe to a temperature T R , the resulting entropy production dilutes the yield of relic vectors Y X ∼ T R /m φ . A concrete example of the such a cosmology occurs if the dark Higgs is able to thermalize with the SM. This generically requires the dark Higgs to have a substantial coupling to the SM, and as a result the allowed mass range of φ will be subject to experimental constraints. The simplest interaction of the dark Higgs with the SM is a Higgsportal coupling. As we show in Appendix A, this has severe constraints from star cooling and rare meson decays below around 5 GeV. We show the viable parameter space for vector DM production in Fig. 2 (right), assuming the large dark Higgs yield is eliminated at late times before dominating the energy density of the universe.
Here we fix the dark Higgs mass to be m φ = 10 GeV and show the requirements on dark Higgs thermalization through the Higgs portal interaction, leaving a detailed examination of the necessary conditions to Appendix A. We do not explicitly show the parameter space for vector DM production in the case of dark Higgs domination though we have checked the lower reach in m X is ultimately the same as that in the case of no entropy production.

B. Case 2: e λ
We now turn to the limit where the gauge coupling is much larger than the quartic (and so m X m φ ). Due to the effects of post-scattering, the co-moving number densities of φ and X at late again become comparable and are given by, This difference in mass of φ and X leads to vectors dominating the energy density at late times, and the dark Higgs becomes a subdominant component with a fractional abundance λ/e. In addition, vectors are produced non-relativistically with typical co-moving momentum k * eφ 0 , while the dark Higgses are dominantly produced from vector fluctuations with a similar spectrum and are thus highly relativistic.
The observed DM abundance is reproduced for the initial field amplitude of The temperature at which the vectors become nonrelativistic is given by: Note that by T NR , the vector mass (initially dominated by fluctuations of φ after PR) assumes the vacuum value.
Since the vector makes up most of the DM, we require the coldness constraint T NR keV. On the other hand, if the fraction of produced dark Higgses is roughly greater than 10 −2 , then this subdominant component must also be sufficiently cold. We show the viable parameter space for vector DM production in Figure 3 for e/λ = 10 (left) and e/λ = 10 3 (right), with the value of φ 0 fixed at every point to achieve the correct relic abundance. For e/λ = 10 the dark Higgs is a non-negligible subcomponent and in addition to the vectors being sufficiently cold we also require the dark Higgs is non-relativistic by a keV, while for e/λ = 10 3 we only require that the vector population satisfies the coldness constraint (37). The lowest possible vector masses can be obtained by saturating e → λ where we find we can produce cold DM for m X 10 −18 eV (though saturating this limit results in the vectors being accompanied by non-negligible dark Higgs abundance).

V. PHENOMENOLOGY
In this section, we summarize some distinctive features of vector production through parametric resonance which could be used to differentiate it from other non-thermal cosmologies.
a. Dark Higgs Perhaps the most prominent prediction of PR production would be searching directly for the accompanying light scalar. The detectability of the scalar depends on its model-dependent coupling (if any) to the SM, and in general no such coupling is required to produce vectors. However, if λ e than the scalar is either a large fraction of the DM abundance today, or the scalar is destroyed by some additional mechanism (e.g. thermalization with the SM) so that vectors make up all of the observed DM abundance. If the scalar is a non-negligible relic today then it could be searched for directly through experiments sensitive to light scalars. Furthermore, if it dominates the DM density then it could be observed as a (cosmologically slow) dark decay into the vectors from anomalous changes of equation of state of the universe [45]. The current consistency with the ΛCDM picture allows us to set a bound on this decay rate as given in 2.
Alternatively, if the dark Higgs is assumed to thermalize with the SM then the minimal required coupling to achieve thermalization sets a convenient target for experimental searches. We study these specific requirements in the context of a Higgs portal coupling in Appendix A.
b. Cosmic Strings As we have seen, PR produces large quantum fluctuations in the X and φ fields. These fluctuations can lead to a large positive effective mass for φ resulting in the symmetry being temporarily restored once PR terminates and a subsequent non-thermal phase transition once the mass of φ becomes negative [46]. This has an intriguing prediction of the formation of cosmic strings [47]. Cosmic strings are one-dimensional topological defects, characterized by a string tension µ ∼ v 2 . After formation, it is expected that the string network quickly approaches a scaling regime, i.e., energy density in strings scales with the energy density of the universe but roughly suppressed by the factor Gµ, where G = 1/8πM 2 pl is Newton's constant. Such strings have several characteristic predictions owning to their induced large energy gradients in the universe. Perhaps the most robust detection of cosmic strings can be extracted from the cosmic microwave background, whose gravitational interaction would induce small temperature distortions leading to inhomogeneities in the temperature map [48]. Using the WMAP data with a combination of cosmological observations such strings have yet to be observed, putting a constraint Gµ 10 −7 [49]. Viable parameter space for ultralight vector production using parametric resonance in the limit that e λ. We fix the initial condensate amplitude, φ0, to obtain the correct relic abundance at every point in the figures. Shown are constraints from coldness (yellow), cosmic strings (red), and sufficiently long PR (blue) as described in the text. Both plots have vectors making up nearly all the DM. Left: We take e = 10λ. Right: We take e = 10 3 λ.
An additional prediction of cosmic strings comes from gravitational radiation emitted by the string oscillations. The evolution of a scaling cosmic string network is expected to contribute to the stochastic gravitational wave background [50] as well as induce gravitational wave bursts [51]. This is contrast to global strings, which predominantly radiate massless Goldstone bosons (e.g., axion strings). The gravitational wave spectrum from a cosmic string network can be computed, under basic assumptions. The, thus far, null observation of a stochastic gravitational wave background by LIGO and pulsar timing arrays constrain Gµ 10 −11 [52][53][54], which roughly translates to a bound on the VEV v 10 14 GeV. Future pulsar timing array measurements are expected to have improved sensitivity with the upcoming future Square Kilometer Array [55] and provide an opportunity to probe these non-thermal phase transitions.
c. Isocurvature Perturbations Another prediction of this production mechanisms is due to the lightness of the dark Higgs, inducing isocurvature perturbations in the cosmic microwave background (CMB). During inflation, we presume φ is stuck with an initial field amplitude obeying λφ 0 H inf , and fluctuations, δφ ∼ H inf /2π. During PR the energy density of the φ condensate is transferred to the observed DM abundance and instills these isocurvature perturbations in the dark matter spectrum. These perturbations can be looked for in the CMB though they have yet to be seen [26]. This can be interpreted as a bound on the Hubble scale during inflation H inf 3 × 10 −5 φ 0 , which in the simplest picture suggests a bound λ 3 × 10 −5 . This puts a relevant constraint for λ e if the dark Higgs is required to thermalize with the SM but turns out to be negligible when we do not enforce this requirement. We note that, in principle, this isocurvature perturbation can be suppressed if the Hubble induced mass of φ is larger than H inf .

VI. DISCUSSION
In this work we present a new production mechanism for vector DM in the early universe through its (possible) coupling to a dark Higgs. The mechanism relies on the non-perturbative dynamics associated with parametric resonance, thus allowing the produced vectors to be ultralight while still being consistent with the stringent constraints on warm DM.
Vector production from parametric resonance has qualitative differences from the well-studied theory of scalar production. We study the equations governing the PR production of transverse and longitudinal modes and present an instability chart. For λ e the transverse mode production is highly inefficient while the longitudinal mode is rapidly produced (this can be understood as a consequence of the Goldstone boson equivalence theorem). Fluctuations of dark Higgses are also produced which results in a Higgs-dominated dark sector. In order for vectors to make up the entire DM abundance, additional interactions can be considered to thermalize the dark Higgs with the visible sector. We find produced vectors can be as light as 10 −20 eV if they form 1% of the energy density in DM (with the dark Higgs making up the rest). If, on the other hand, we require the dark Higgs to thermalize with the SM it is difficult to foresee a viable model without making the dark Higgs heavier than around 10 GeV (otherwise there are tight constraints on its coupling). This restricts the produced vector DM to having masses above around 10 −4 eV. In the case where e λ both the transverse and longitudinal mode can be efficiently produced, though as a consequence of initial conditions set by inflation we still expect the longitudinal mode to dominate for a wide range of parameters. As in the previous case fluctuations of φ are rapidly produced resulting in comparable number densities between vectors and dark Higgs, although due to the ratio of masses the DM energy density today is dominated by vectors. Ultimately, the coldness constraint restricts the viable vector DM mass to be above 10 −18 eV.
Our study of PR production of ultralight vectors was not meant to be exhaustive, and there is still much work to do to better understand the viable parameter space and implications of such a production mechanism. We conclude by commenting on directions we feel merit further attention. Firstly, the focus of this work was entirely on vectors which get their mass from a dark Higgs. In principle, this could easily be generalized to other types of scalars which obtain a large field value. Secondly, in this work we did not attempt a complete lattice simulation of the non-linear effects. This would be particularly important in the limit of e λ since in this case it is possible that the coldness constraint is significantly weakened if the vectors do not get boosted to their maximum possible momenta, eφ 0 . Furthermore, it is important to note that a general feature of this framework is the necessity for tiny couplings (for the mass range in the e λ case, gauge couplings in the viable parameter space go down to as low as ∼ 10 −40 ). While such couplings are technically natural, it would be interesting to see how viable these are in a UV model. Lastly, in this work we briefly explored the prominent phenomenological signatures of this production mechanism though it may be fruitful to consider these in more detail as well as others to differentiate PR production from other possible production mechanisms. mass. Therefore, the Higgs portal coupling gives a mass correction δm 2 φ ∼ y 2 T 2 osc to the dark Higgs at the onset of oscillations. If this is larger than the assumed effective mass-squared m 2 eff ∼ λ 2 φ 2 0 , the dark Higgs would have oscillated earlier with a large frequency (set by the size of y 2 ), thus rendering the PR production of vector modes narrow and relatively inefficient. Requiring this thermal contribution to be sufficiently small δm 2 φ λ 2 φ 2 0 translates into an upper bound on the mixing angle, For Fig. 2 where we fix m φ 10 GeV, this is the dominant upper bound on the mixing angle. A second possible effect of the Higgs portal coupling is PR production of the SM Higgs. For thermalization we generically require κ λ which would suggest that PR may be efficient in producing the SM Higgs; this is ultimately not the case due to the large thermal mass of the Higgs. We can verify this by computing the adiabatic parameter R = |ω|/ω 2 for the SM Higgs with a timedependent frequency ω(t) k 2 + m H (T ) 2 + 1 2 y 2 φ(t) 2 . Since the thermal mass is strictly larger than the mass contribution from the oscillating field at the time of PR (y 2 φ 2 0 T 2 osc ), the adiabatic parameter is always less than unity. We conclude that there is no significant nonperturbative production of the SM Higgs.
There are a number of phenomenological constraints on a Higgs-portal scalar. These constraints depend on whether the scalar decays visibly or appears invisible (at least on detector scales.) For detailed review see, e.g., [60]. For low masses m φ MeV, we find stellar cooling constraints [61] are in tension with the requirements of thermalization. For intermediate masses, MeV m φ 300 MeV constraints from Supernova 1987A and rare Kaon decays are powerful, but it may be possible to have a scalar in this mass range with a sufficiently large θ consistent with thermalization if it sits in the small gap between these constraints. For masses GeV m φ 5 GeV, rare B meson decays roughly constrain θ 10 −3 . Above this scale, scalar production at LEP constrains θ 10 −1 although this is weaker than the condition (A3). Interestingly, if φ decays visibly then lower values of θ in this mass range could also be probed by future experiments designed to look for long-lived particles [62][63][64][65][66].
Lastly one may wonder if the dark Higgs ever dominates the energy density of the universe. Dark Higgs domination will take place if the temperature of thermalization is less than T eq (m φ /m X ), or in terms of the mixing angle: We see that for m φ 10 GeV and the relevant vector mass range m X 10 −4 in Fig. 2, this is never the case as long as we satisfy the condition of thermalization.