Vector dark matter production from inflation with symmetry breaking

We present a scenario of vector dark matter production from symmetry breaking at the end of inflation. In this model, the accumulated energy density associated with the quantum fluctuations of the dark photon accounts for the present energy density of dark matter. The inflaton is a real scalar field while a heavy complex scalar field, such as the waterfall of hybrid inflation, is charged under the dark gauge field. After the heavy field becomes tachyonic at the end of inflation, rolling rapidly towards its global minimum, the dark photon acquires mass via Higgs mechanism. To prevent the decay of the vector field energy density during inflation, we introduce couplings between the inflaton and the gauge field such that the energy is pumped to the dark sector. The setup can generate the observed dark matter abundance for a wide range of the dark photon's mass and with the reheat temperature around $10^{12}$ GeV. The model predicts the formation of cosmic strings at the end of inflation with the tensions which are consistent with the CMB upper bounds.


Introduction
The nature of dark matter is still unknown and remains one of the most compelling evidence for the physics beyond the standard model of particle physics (SM). While the most well-known candidate are the so-called Weakly Interacting Massive Particles (WIMPs), in the absence of observation of such new particles in colliders the researchers tend to look for other possibilities as well (see Refs. [1][2][3][4][5] for attempts in this regard). Recently, dark photon became one of the candidate dubbed "vector dark matter" [6][7][8][9]. They are supposed to have negligible direct interactions with the SM particles while they can have significant direct coupling with other particles in the dark sector as well as the universal minimal coupling to gravity. It is then shown that they can be efficiently produced before the time of matter radiation equality and play the role of dark matter relic density [10][11][12][13]. Although dark photons interact only weakly with the SM particles, they can be probed indirectly through the gravitational waves observations and this possibility has been the subject of many studies in recent years [14][15][16][17][18][19].
On the other hand inflation has emerged as the leading paradigm to solve the horizon and the flatness problems associated with the big bang cosmology and also to explain the origin of the large scale structures in the Universe. Among the basic predictions of simple models of inflation are that the primordial perturbations are nearly scale invariant, nearly Gaussian and nearly adiabatic which are well consistent with the cosmological observations [20,21]. The inflaton field is usually taken to be a scalar field which slowly rolls on top of a nearly flat potential leading to a quasi-de Sitter background expansion. However, there is no primary reason which forbids considering other fields like gauge bosons and even fermions during inflation. Evidently, it is not an easy task to find a consistent inflationary dynamics completely based on the fermionic fields while inflationary models based on the gauge fields are widely studied in recent years. Indeed, it is possible to find an inflationary scenario completely based on vector/gauge fields [22,23] or considering scenarios in which gauge fields partially contribute to the dynamics of inflation. One problem associated to the presence of the gauge fields during inflation is the conformal symmetry which causes the exponential decay of the energy density of the gauge fields in a quasi-de Sitter background. One natural possibility to remedy this issue is to consider a direct coupling between the inflaton and the vector/gauge fields. In this respect the gauge field drags energy from the inflaton which can prevent the decay of its energy density. There are some models based on this picture among which are the so-called anisotropic inflation [24][25][26][27][28] and pseudoscalar inflaton [29][30][31][32]. This idea was already implemented to explain the origin of the primordial magnetic field [29,[33][34][35][36][37][38][39][40][41][42][43].
Dark matter and inflation are usually considered to be unrelated. However, inflation and reheating can be the frameworks for the productions of particles beyond the SM. As the dark matter particles may be from beyond the SM sector, it is interesting to look for the possibility of the production of dark matter particles during inflation. In this regard, it was shown that if one considers a massive dark gauge field during inflation without any direct coupling to inflaton, the associated longitudinal mode can be responsible for the dark matter energy density after the time of matter and radiation equality [8]. Very recently, this idea was further extended and it was shown that even the transverse modes can be responsible for the dark matter relic if a massive gauge boson directly interacts with the inflaton [44][45][46][47]. The difficulty with the dark matter models based on the transverse modes of the gauge fields produced during inflation is that in one hand, the gauge field should have a small mass compared to the energy scale of inflation to have efficient dark photon productions during inflation. While, on the other hand, the mass should be large enough to make the produced dark photons non-relativistic before the time of matter and radiation equality. These criteria make the appropriate mass range of the dark photon to be very model-dependent. For example in Ref. [44] with the parity violating interaction between the inflaton and the gauge field the mass of dark photons should be m A ≥ 10 −6 eV while in a model with non-minimal conformal coupling one finds m A ≥ 10 −21 eV [47].
In this paper, we consider a scenario in which the dark gauge field acquires mass dynamically through a symmetry breaking, i.e. the Higgs mechanism, at the end of inflation. Apart from the fact that the Higgs mechanism is a natural way of giving mass to a gauge boson in the SM or beyond SM, it opens up a new window to the scenarios of vector dark matter production during inflation. The reason is that the gauge field is massless during inflation and therefore dark photons can be produced very efficiently. The gauge field then acquires mass only toward the end of inflation and one can easily find appropriate range of parameters for which the produced dark photons can become non-relativistic before the time of matter-radiation equality and also to give the correct dark matter abundance. In this model the mass of dark photon can easily be so large that they become non-relativistic even right after inflation. In this respect, we will be able to find new novel mechanisms which were not possible in previous works due to the assumption of small mass for the dark photons during inflation.
The rest of the paper is organized as follows. In Section 2 we present a general inflationary model in which a dark gauge boson acquires mass through the symmetry breaking at the end of inflation. In Section 3, we show that dark photon can be efficiently produced during inflation and they can play the role of (vector) dark matter after inflation. We show that the symmetry breaking scenario opens a large parameter space for these models. In Section 4 we employ our general scenario to a particular case and show that it has a rich phenomenology. In Sections 5.1 and 5.2 we consider two subsets of this specified model which resemble the already well-known models of kinetic conformal coupling and the parity violating coupling respectively. In Section 5.3 we then look at the more general case which includes both of the conformal and parity violating couplings. Section 6 is devoted to the summary and conclusions.

Inflation with symmetry breaking
In this section we present an inflationary model in which a gauge boson in the beyond SM dark sector (dark photon) obtains mass through the dynamical symmetry breaking at the end of inflation. The symmetry breaking mechanism can be achieved via a waterfall phase transition as in hybrid inflation [48,49]. Note that the details of the inflationary model is not important for our computation of the dark photon evolution and its relic energy density. The only important aspect of the model is that the symmetry breaking occurs approximately at the end of inflation such that the gauge field is massless during inflation but massive afterwards. Having said this, it is possible to consider a model in which the symmetry breaking takes place during inflation but we leave this possibility to a future study.
The standard scenario where the symmetry is broken at the end of inflation is the hybrid inflation in which the inflaton is a real scalar field and there is another complex scalar field, the waterfall field, which is responsible for the symmetry breaking and the termination of inflation. We extend this standard picture by assuming that the waterfall field is charged under the dark photon field. The theory is described by the following action [50,51] where M Pl = 1/ √ 8πG is the reduced Planck mass with G being the Newton constant. As mentioned before, the inflaton is the real scalar field φ and the waterfall is the complex scalar field ψ which is charged under the U (1) gauge field A µ . The components of the covariant derivative of the waterfall field is given by where e is the dimensionless gauge coupling between the gauge field A µ and ψ. As usual, the components of the field strength tensor F associated to the gauge field is given by andF σρ = σρδη F δη /2 is the corresponding dual field strength tensor. We have added the conformal couplings f (φ) and h(φ) to break the conformal invariance and to pump energy from the inflaton sector to the gauge field sector. Indeed, it is well known that in the simple Maxwell theory with no such couplings, the electromagnetic fields energy density are exponentially diluted during the quasi-de Sitter inflationary expansion. To prevent this, one usually introduces couplings via appropriate forms of f (φ) and h(φ). This formalism was employed extensively in the context of anisotropic inflation where it is shown that with appropriate functional form for f (φ) with h(φ) = 0, the energy density of the gauge field survives the de Sitter expansion while its perturbations become nearly scale invariant [24][25][26][27][28]. For a review of anisotropic inflation see Ref. [52]. On the other hand, in the pseudo-scalar inflation with the parity violating coupling h(φ) 2 ∝ φ and f (φ) = 1, vector field perturbations are enhanced through a tachyonic instability [29][30][31][32]. Both of these scenarios are also implemented as a mechanism for generating primordial magnetic field during inflation [29,[33][34][35][36][37][38][39][40][41][42][43]. We will study the production of dark photon during inflation in the following section. Here, we describe the inflationary model.
We assume the axial symmetry so the potential depends on the modulus of the waterfall field |ψ| along with the inflaton field. We thus decompose the waterfall field into the radial and angular components as follows ψ ≡ |ψ|e iθ , (2.4) where the field θ determines the argument of the complex waterfall field. As usual, the action (2.1) is invariant under the local gauge transformation As a result, we can exploit this gauge symmetry to set θ = 0, i.e. the unitary gauge. In this gauge, we deal with a real waterfall field in practice and, therefore, from now on we do not write the absolute value sign anymore. The typical choice for the potential in hybrid inflation is [48,49] Figure 1: A schematic view of symmetry breaking based on hybrid inflation. The red line is the chaotic potential for φ while the blue curve is the potential of ψ which is initially very heavy but at φ c becomes tachyonic and soon settles at its global minimum at ψ min . Note that we always have ψ ≥ 0.
where λ and g are two dimensionless couplings. The potential is depicted in Fig. 1. Initially the inflaton starts in the region φ > φ c . Here, φ c ≡ M/g is a critical value beyond which the waterfall field ψ has positive mass squared. The assumption is that ψ is very heavy so during inflation and it is stuck in its local minimum ψ = 0 and one can safely ignore its quantum fluctuations. We assume that inflation is mainly driven by the vacuum term M 4 /4λ. In other words the potential during inflation is This is crucial because we need inflation to end after the symmetry breaking and it requires that the waterfall dominates the energy density during inflation. As the inflaton field rolls during the expansion of the Universe and reaches the critical value φ c , the waterfall field becomes tachyonic and rapidly rolls towards its global minimum at ψ min = M/ √ λ (note that ψ is the absolute value of the waterfall field so it can only be positive) terminating inflation efficiently. We define the following dimensionless parameters which measure the ratios of the mass scales of the inflaton and the waterfall field to the Hubble parameter H during inflation. The condition that φ is slow-rolling implies that α 1 while the assumption that ψ is very heavy during inflation such that we can ignore its quantum fluctuations, requires β 1. Further, we want the procedure of phase transition happens very rapidly within less than one e-fold of inflation which implies αβ 1 [53][54][55]. Since during inflation the waterfall is stuck in its local minimum ψ = 0, the induced mass of the dark photon field is zero during inflation, m A = 0. As a result, the dark photon mode function receives no corrections from the mass term and it is given by the conventional massless modes in a nearly fixed de Sitter background. As inflation ends and the waterfall field rapidly settles to its global minimum, the induced mass of the dark photon is given by (at the end of inflation) . (2.9) We also define the ratio of the gauge field mass to the Hubble scale at the end of inflation H e as follows Depending on the three parameters e, λ and β we can have dark photon heavier (R > 1) or lighter (R < 1) than the Hubble scale at the end of inflation H e . The above scenario was a working example of symmetry breaking at the end of inflation to yield mass to the dark photon via the Higgs mechanism. The specific setup discussed above is based on hybrid model but we emphasize that this picture is more general and not restricted to the particular setup of hybrid inflation. Indeed, the standard hybrid model is ruled out from the CMB observations as it generates a blue spectral index for the curvature perturbation power spectrum [56]. Instead, all we need is a phase of inflation with a nearly constant value of H terminated by a heavy waterfall field which is not excited during inflation. There are various ways that the standard hybrid setup can be modified such that the setup will be consistent with the CMB observations. One simple extension is to consider a multiple field models where more than two fields can drive inflation and generate the large scale perturbations. One can massage the inflationary potential such that on the observed CMB scales the spectrum is red-tilted while towards the final stage of inflation the model effectively looks like Eq. (2.6) with one field essentially playing the role of the inflaton field. Alternatively, one may embed our setup in the context of brane inflation [57,58] in string theory where the inflaton field is the geometric distance between a stack of branes and that of anti-branes in the string theory compactification. As the branes and anti-branes approach each other a tachyon develops and inflation terminates rapidly via the branes and anti-branes annihilation. One may think of the dark photon as living in some left over brane in these configuration which is different than our three-brane where the SM particles are located.
A generic prediction of our setup is that cosmic strings are produced at the end of inflation. This is a general outcome when a U (1) symmetry is broken in early universe [59,60]. Unlike monopoles and domain walls, cosmic strings are allowed to form in cosmic history with interesting cosmological implications such as lensing, CMB anisotropies or generating stochastic gravitational waves [60]. The tension of string µ is related to the scale of symmetry breaking via [60], The CMB observations impose the upper bound Gµ 10 −7 [61]. We will use this upper bound on µ to impose some bounds on our model parameters such as the value of R and the reheat temperature.
In the following sections we will focus on the production of dark photons during inflation and their evolution afterwards to form dark matter observed today.

Dark photons: Production and evolution
As we explained in the previous section, around the time of end of inflation, the dark gauge field acquires mass from the waterfall field. Thus, the dark photon sector of the action (2.1) takes the following form Note that the mass term in the above action is generated dynamically through the symmetry breaking scenario that we presented in the previous section. Therefore, during inflation dark photons are massless m A = 0 while at the end of inflation they almost instantaneously acquire mass m A = 0. This is the key feature of our model which makes it completely different from other vector dark matter models which are based on the inflationary scenarios. Taking variation of the action (3.1) with respect to A σ , we find the equation of motion for the gauge field as The energy momentum tensor of the gauge field can also be obtained by taking the variation of Eq. (3.1) with respect to the metric, yielding We consider a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric for the inflationary background geometry in which a(τ ) is the scale factor and τ is the conformal time related to the cosmic time via dτ = dt/a(t).
For the gauge field sector, our assumption is that the gauge field has no background values A σ = 0. Taking this fact into account, from the temporal component of Eq. (3.2) we find where a prime denotes derivative with respect to the conformal time. Therefore, as usual, the component A 0 is a constraint and its solution can be imposed at the level of the action. Considering the background symmetry we decompose the spatial components of the gauge field as where χ and A T i are longitudinal and transverse parts of the gauge field respectively. Substituting the above relation in Eq. (3.5), we can solve A 0 in favour of the longitudinal mode χ as follows The energy density of the gauge field is defined by ρ A ≡ −T 0 0 . We decompose it into ρ A = ρ T +ρ L +∂ i J i in which ρ T and ρ L are the energy densities for the transverse and longitudinal modes respectively. The term ∂ i J i is a total divergence term which after taking the vacuum expectation value does not contribute to the background [47]. Now, plugging Eq. (3.7) in Eq. (3.3) and after a couple of lines of algebra, we find Substituting Eq. (3.6) into the action (3.1) and then eliminating A 0 in favour of the longitudinal mode χ through Eq. (3.7), we can also decompose the gauge field quadratic action into the transverse and longitudinal parts as follows where we have defined (3.11) During inflation the gauge field is massless m A = 0 and therefore the longitudinal mode is absent.
After the symmetry breaking we have m A = 0 and the longitudinal mode is dynamical with a large mass. As the conformal and parity violating couplings are stabilized after inflation the massive longitudinal mode cannot be excited efficiently. As a result, from now on we neglect the longitudinal mode and focus only on the transverse modes.

Quantization
To study the production of dark photons, we need to look at the quantum fluctuations of the transverse modes of the dark gauge field during inflation. We define the canonical field V i ≡ f A T i and expand it in terms of the creation and annihilation operators a k and a † k as usual where the time dependence of the gauge field is encoded in the mode functions v k,λ (τ ). In the above relation, ε λ i (k) are the polarization vectors for λ = ± which satisfy ε λ During inflation, the field is massless m A = 0 so substituting from Eq. (3.12) in the action (3.10), and taking the variation, the equation of mode function is given by Note that for nonzero time-dependent coupling h the mode functions for different polarizations evolve differently.
The energy density of the gauge field during inflation can be computed from Eq. (3.8) (3.14) We will use this expression in computing the backreaction of the vector field on inflationary background dynamics. However, at the end of inflation the dark photon becomes massive through the symmetry breaking process. As a result, after that the dark photon is a massive vector field with a mode function denoted by u k,λ and its equation of motion reads Note that the symmetry breaking occurs at the end of inflation so the above equation is relevant only around the time τ = τ e . We want to compute the energy density of the dark photons after acquiring mass at the end of inflation which from Eq. (3.8) is given by This will be carried to the radiation and matter dominated eras as the observed relic abundance of the dark matter. To solve Eq. (3.15) exactly, we need to take into account the dynamics of the waterfall field during and after the phase transition (see Refs. [53,62]). However, since we have assumed that the waterfall field is heavy enough, the symmetry breaking happens almost instantaneously around the time of end of inflation τ e [48,53], and we can avoid the technicalities of considering the waterfall field dynamics.
Assuming an instantaneous phase transition, we only need to impose the junction conditions u k,λ (τ e ) = v k,λ (τ e ) and u k,λ (τ e ) = v k,λ (τ e ) which enables us to find the energy density of the dark photons at the end of inflation in terms of the massless mode function v k,λ as follows , (at the end of inflation) . (3.17) As we will see in the concrete models in the following, we can parametrize the energy density of the produced dark photons at the end of inflation as where H e is the Hubble parameter at the end of inflation and C is a dimensionless function which depends on the parameters of the system under consideration. We will find explicit form of C for some particular cases in the next sections.
In the following subsections we use expressions (3.14) and (3.17) for the backreaction constraint and the relic abundance of the dark matter respectively. We remind that Eq. (3.14) is the energy density accumulated from the gauge field fluctuations during inflation while ρ e in Eq. (3.17) is the energy accumulated at the end of inflation including the effects of mass term m 2 A . The energy density of the dark photon during the big bang expansion is related to ρ e via ρ A ∝ ρ e /a 4 or ρ A ∝ ρ e /a 3 , depending on whether the dark photon is relativistic or non-relativistic.

Constraints on the model
Our assumption is that the gauge field does not destroy the coupled equations of motion for the inflaton and the background geometry. In this respect, there are two sources of backreaction from dark photons to the dynamics of the inflaton and the geometry. Correspondingly, there are two conditions to be imposed. The first condition is that the energy density of dark photons must be small compared to the vacuum energy during inflation, i.e.
Second, the source term that appears in the Klein-Gordon equation due to the conformal couplings must be small. The equation of motion for the inflaton reads Note that the right hand side of the Klein-Gordon equation is nonzero because of the couplings f (φ) and h(φ) between the inflaton and the gauge field as given in the action of gauge field Eq. (3.1).
In the FLRW background (3.4), the equation of motion for the inflaton (3.20) takes the form are the corresponding electric and magnetic fields of the dark gauge field and we have used df /dφ = f /aφ with a dot denoting the derivative with respect to the cosmic time t. The second backreaction condition then is given by For the later purposes let us write S φ given in Eq. (3.21) in terms of the mode functions The conditions (3.19) and (3.22) should be satisfied in order to have a consistent inflationary background. Note that both of the above conditions must always be valid during inflation. However, since the gauge field energy density and also S φ are growing, the two conditions are strongest at the end of inflation. Thus we will evaluate them at the end of inflation.
Moreover, we should look at the effects of the gauge field perturbations on the curvature perturbations as well. We need to demand that the corrections in power spectrum, bispectrum and trispectrum induced from the gauge field perturbations are suppressed with respect to the standard curvature perturbations given by the perturbations of the inflaton field. Since the gauge field in our model does not have a vev, demanding that the energy density of the vector field is negligible in the inflationary background is sufficient for this purpose. However, if we assume that inflation has started much earlier than the observed number of e-folds on CMB then the accumulated IR effects of the vector field perturbations can effectively change the background and put stronger constraints than (3.19) and (3.22). This possibility has been studied for example in [63,64] for the case of the kinetic conformal coupling, i.e. h = 0. The result of these studies is that the total number of e-folds cannot be much different from the observed number of e-folds on CMB. Here, we take the conservative approach that the total number of inflation N is not very long in the past. In particular, to solve the flatness and the horizon problems we assume N ∼ 60. As a result, it is enough to only consider the backreaction constraints Eqs. (3.19) and (3.22).
Besides the above backreaction conditions, we also have to take into account the constraints from the primordial black hole (PBH) formations and the isocurvature perturbations. This is somewhat model-dependent and we impose the constraints from PBHs and isocurvature perturbations after presenting our models in some details.

Relic abundance
The dynamics after inflation is in principle complicated and needs numerical analysis. One reason is that the oscillating inflaton field can cause production of dark photons during (p)reheating through the couplings f and h [65,66]. To take the problem under theoretical control we assume that the conformal couplings are such that they are stabilized at the end of inflation and afterwards to the trivial values f | τe = 1 and ∂ τ h| τe = 0. In other words we assume that the dependence on φ in the couplings f and h becomes very weak or disappears after inflation. With these assumptions the dark photons may not be produced efficiently after inflation and the resultant gauge field energy density at the end of inflation is just redshifted to the late time universe. In addition, we assume that the reheating takes place instantaneously in which all energy density of the inflaton (minus the increase of the dark photon energy density due to the mass term acquired at the end of inflation) is transferred to the SM fields.
With these simplifications the equation for the mode function (3.15) after inflation in terms of cosmic time becomesü where we have dropped the polarization index in the mode function since the mode function for both polarizations are the same. From Eq. (3.24), we see that the evolution of the mode function and therefore the evolution of the energy density of the dark photon after inflation depends on the hierarchy among three scales k/a, m A and H. In the case that the mass scale is subdominant compared to the other two scales, then the equation is approximately the same as that of a massless vector field and the energy density after inflation evolves as radiation ρ A ∝ a −4 . However, if we assume that the mass scale is dominant then we get the following approximate solution and then in this regime ρ A ∝ 1/a 3 . For this to be true, the mass scale must be dominant over all relevant momenta, specially the momenta with significant contribution in the spectrum of the dark photons. As mentioned before, ρ A is related to the total energy accumulated from the gauge field fluctuations at the end of inflation ρ e given in Eq. (3.17). Since ρ e is constructed from the accumulated energy of all perturbations, we have to check the dominance of the mass scale over all of the modes which exhibit growth during inflation. Let us denote k min and k max as the first and the last modes which exhibit growth during inflation. Specifically, k min and k max are the modes that become tachyonic at the initial and end of inflation respectively (see the following section for more details). As we will see k min and k max are proportional to the scale of the horizon aH at the initial and end of inflation respectively with model dependent coefficients. Then we can roughly write k min ∼ e −N k max , where N is the duration of observed inflation. Therefore, when imposing the conditions above it is actually sufficient to impose the condition on the shortest mode which exhibits growth during inflation, i.e. m A k max /a. As a result, we define non-relativistic states as those satisfying both conditions m A H and m A k max /a so that ρ A ∝ a −3 . As a first order approximation, we assume an instantaneous transition from a −4 to a −3 behaviour at the time t = t NR when both non-relativistic conditions are met. Now the natural question is which one of the conditions m A ≥ H and m A ≥ k max /a happens later, guaranteeing the non-relativistic conditions. Let us parameterize the shortest growing scale as following where κ is a dimensionless function of the parameters of the problem which relates k max and the horizon scale at the end of inflation a e H e . First of all, note that comparing the two scales we have However, after inflation the universe is decelerating so aH is a decreasing function of time. As a result, if κ > 1 the above ratio is always bigger than unity which means that both non-relativistic conditions are satisfied when the condition m A ≥ k max /a is met. However, for parameters that κ < 1 this should be treated with care. More specifically the non-relativistic conditions are satisfied when a e a NR = R 1/4 , or a e a NR = depending on the fact that for which one a NR is larger, i.e. the condition is satisfied later. We can then imagine different scenarios depending on the values of R and κ. If R is very large such that R > κ 2 > 1, then both conditions are met right at the end of inflation so t NR happens right after reheating. This is a unique feature of our setup based on symmetry breaking which allows for R > 1 as well as efficient dark photon production during inflation. However, when R is not very large, for example R < κ 2 or R < 1, then t NR occurs sometime later than the time of reheating. Note that t NR must not be later than the time of matter-radiation equality for which we must have ρ A ∝ a −3 as a genuine dark matter candidate.
As mentioned before, to simplify the analysis we assume an instantaneous reheating for which the temperature T r is given by in which g * ,r = 106.75 is the number of relativistic degrees of freedom for energy density at the time of reheating. The energy density of the dark photon evolves like a −4 before the time t NR and evolves like a −3 after that. As a result the energy density of the dark photon today ρ A,0 will be where T 0 ∼ 10 −13 GeV is the CMB temperature today, T r is the reheating temperature, g s * ,0 = 3.9, g s * ,r = 106.75 are the relativistic degrees of freedom for entropy density today and at the time of reheating. Finally, ρ e is the energy density of dark photons at the end of inflation given by Eq. (3.17). Plugging the above values, the relic dimensionless energy density today of dark photon then is given by where the parameterization Eq. (3.18) has been used and we have used Eq. (3.29) to write H e in terms of reheating temperature. In each model of inflation with a vector field one can compute C, κ, R and T r . Then the general expression (3.31) gives the relic energy density of dark photons. An important constraint is that non-relativistic conditions must be met before the time of matterradiation equality. This imposes a lower lower bound on R as following where we have used T eq ∼ 10 −9 GeV as the matter-radiation equality temperature. As an estimation of lower bound on R, suppose we take the minimum reheat temperature allowed from Big Bang Nucleosynthesis, T r ∼ 10MeV and set κ ∼ 1. Then the above expression requires R 10 −15 . Of course, by taking higher value of the reheat temperature, say T r ∼ 10 12 GeV one can take lower values for R. Before closing this section let us examine the constraint on the model parameters from the CMB upper bound on the tension of cosmic strings. As mentioned before, the production of cosmic strings is a natural outcome of our setup with the tension µ given in Eq. (2.11). Using the definitions of m 2 A and R in Eqs. (2.9) and (2.10) and the formula for the reheat temperature Eq. (3.29) we can obtain a relation between µ, R and T r as follows ∼ 720 e 2 πg * ,r Gµ ∼ 2e 2 Gµ . In the final estimation we also used the assumption that e 2 < 1. This is because e is the effective gauge coupling at the end of inflation and in order to have the perturbative QFT of gauge field under control we require e < 1. We see that the bound Eq. (3.34) allows for a wide rang of the value of R. For example, taking T r = 10 12 GeV we only require R < 10 17 to satisfy the observational bound on the tension of strings. Increasing the reheat temperature the allowed values of R become smaller. Taking the extreme upper bound T r = 10 15 GeV, we require R < 10 5 . On the other hand, by reducing T r the upper bound on R becomes much more relaxed. For example, taking the lower bound T r ∼ 10MeV we only require R < 10 41 to satisfy the bound on the cosmic string tension.
In the following sections we will consider specific models within the general setup considered above. We first compute C and κ for each model. Then we need to look at the constraints from the backreaction conditions (3.19) and (3.22) as well as other possible constraints coming from the dynamics of perturbations to know the legitimate range of parameters. Finally we will compute the the present energy density of dark photons.

The specific model
As we mentioned before, during inflation the gauge field energy density decays quickly and this issue can be remedied if energy is pumped from the inflaton field to the gauge field through the couplings f (φ) and h(φ) in the action (2.1). Naturally, the conformal couplings f and h are functions of the inflaton field φ. However, during inflation the scale factor is related to the evolution of φ as in which V is given by Eq. (2.7). Therefore, we can imagine that the conformal couplings will be functions of a as well. Motivated from the past works on models of anisotropic inflation and the mechanism of magnetogenesis during inflation, we consider the following phenomenological ansatz for the conformal couplings where we have parametrized h e = γf e by means of a constant parameter γ without loss of generality. As mentioned before, the couplings of the form Eq. (4.3) are common in the context of magnetogenesis [29,[33][34][35][36][37][38][39][40][41][42][43], models of anisotropic inflation [24][25][26][27][28]52], inflation based on a pseudoscalar [29,30], multiple vector field models [67][68][69][70] and Schwinger mechanism during inflation [71][72][73]. In particular, the so-called anisotropic inflation corresponds to p 2 and γ = 0 but with a nonzero background value for the gauge field [24]. Moreover, as we will show, the case with p = 0 and q 1 corresponds to inflation based on a pseudoscalar field [30]. Of course, we can consider the case with p = 0 and γ = 0 and look for the effects of both terms simultaneously which is the subject of subsection 5.3.
As inflation ends, f and h reach their final values f e and h e , and stop evolving. We scale f e to be unity which leads to h e = γ. The kinetic term for the gauge field then takes the standard form and the parity violating term becomes a total derivative term so after inflation the dynamics of the gauge field is given by the Maxwell theory. An important criteria to take into account is that the effective gauge coupling is f −1 . In order to have a perturbative theory, we require f −1 < 1 throughout inflation, otherwise the theory becomes strongly interacting and the perturbative analysis cannot be trusted. This criteria prohibits negative values of p. So, in order to bypass the strong coupling problem we assume p > 0 throughout our analysis. Note that there is not such a constraint on the value of q.
Assuming the specific form of the couplings given in Eq. (4.3), the mode function Eq. (3.13) takes the following form where we have assumed that |q − p| 1 so that we can ignore the time dependence of the parameter ξ and treat it as a constant. As before, the special case of kinetic conformal coupling corresponds to ξ = 0 for which in the particular case of p = 2 we find the usual term 2/τ 2 in the mode function frequency which is the hallmark of the massless scale invariant perturbations in a fixed de Sitter background. In the case of p = 0 we have the axion-like coupling of the inflaton to the vector field similar to the setup of pseudoscalar inflation [30] where the time dependence of ξ is negligible to first order of slow-roll parameter. Note that such a model is considered before for example in Refs. [40,41], where they have set q = p and arrived at a similar equation of motion as Eq. (4.4). However, in that case the limit of p → 0 becomes nontrivial while with our definition we can obtain all limits consistently.
The solution for the mode function can be obtained from Eq. (4.4) assuming the Bunch-Davies initial condition as following in which W µ,ν (z) is the Whittaker function. By using the asymptotic behaviour W µ,ν (z) → z µ e −z/2 for z → ∞, one can show that v k,λ (τ ) approach the standard Bunch-Davies solutions at early times. Further, W µ,ν (z) → z 1/2−ν Γ(2ν)/Γ(ν + 1/2 − µ) for z → 0 so the super-horizon behaviour is given by First of all, we see the dependence on ξ is always through the combination λξ and, therefore, we assume that ξ > 0 without loss of generality. As we see, a nonzero value of ξ induces chirality in the spectrum of the produced dark photons so that one of the polarizations is enhanced significantly compared to the other. Furthermore, from the point of view of the equation of motion (and also from the properties of the Whittaker function) we see that the transformation p → 1 − p (or equivalently ν → −ν) does not change the properties of the mode function. Since we are interested in p ≥ 0 to avoid the strong coupling limit, it is then sufficient to study only the solution for p ≥ 1/2 (or equivalently ν ≥ 0) as the region 0 ≤ p ≤ 1/2 is the same as 1/2 ≤ p ≤ 1. Note also that the limit of axion-like coupling, i.e. p = 0 is the same as p = 1 (or equivalently ν = 1/2). However, we are ultimately interested in the energy density which unlike the mode function explicitly depends on p through the conformal factor f (φ) (see for example (3.14)). As a result, we need to consider the whole region p ≥ 0 and the true limit of axion-like coupling is p = 0.
We are interested in the regime where the mode function grows until the end of inflation. This is most efficient if the mode function undergoes tachyonic growth. To see this more clearly, we define a new variable x ≡ −kτ in term of which Eq. (4.4) takes the form We see that depending on the values of the parameters p and ξ, we may have completely different scenarios. If ξ 2 + p(p − 1) < 0, which is possible only for 0 ≤ p < 1, the frequencies for both polarizations are always positive and therefore there is no tachyonic regime. Note that this is different than the case ξ = 0 where we also do not have any tachyonic growth. On the other hand, tachyonic growth occurs if ξ 2 + p(p − 1) ≥ 0 which is possible for all p > 0. In this case, we have two different regimes during which tachyonic growth happens. First, we note that the condition ξ 2 + p(p − 1) ≥ 0 is satisfied for p ≥ 1 independent of the value of ξ. In this regime, which we call kinetic-like regime, the momenta which start becoming tachyonic right from the beginning of inflation at τ i until the end of inflation at τ e are given by *5 where we have defined x λ ≡ λξ + ξ 2 + p(p − 1). Note that x λ depend on the polarization for nonzero ξ. We will see that the energy density of the dark photon will be chiral in this regime and the chirality is controlled by ξ. Thus, for small ξ both polarizations play roles. The special case of ξ = 0 corresponds to the standard kinetic conformal coupling (as in anisotropic inflation) which is the reason why we called this regime kinetic-like. We will consider this limit in details in the following section. Note that in the kinetic-like case κ = x λ where κ is defined in Eq. (3.26). The second regime is defined by the condition ξ 2 + p(p − 1) ≥ 0 and 0 ≤ p < 1. In this regime, the parameter ξ is restricted as ξ 2 ≥ p(1 − p) and only the plus polarization becomes tachyonic. We thus reasonably call this regime chiral-like. The range of momenta that experience tachyonic growth in this regime are given by (see footnote *5) where we have definedx ± = ξ ± ξ 2 − p(1 − p). Note that up to here everything depends only on the properties of the equation of motion and therefore the results for the regions 0 ≤ p ≤ 1/2 and 1/2 ≤ p ≤ 1 would be similar. In this case κ =x + . We are interested in the energy density at the end of inflation given in Eq. (3.17), which following the parameterization (3.18) and after using Eq. (4.5), turns out to be (4.10) *5 We assume that k max is lower than the comoving wave number corresponding to the diffusion length of the dark photon, which is unknown and depends on the dark electric conductivity and thus on the nature of the interactions of the dark photon with other fields in the dark sector.
Here, the range of integration is x max = x λ , (kinetic-like, p ≥ 1) (4.11) for the kinetic-like regime and for the chiral-like regime. Note that in the former case we have summed over both polarizations as both polarizations are excited while in the latter case we only kept the dominant one which is λ = + in our convention. It can be shown that for p > 2 the dominant contribution of the integral (4.10) comes from its lower limit which is very small. As a result, it is sufficient to consider the leading contribution of the integrand in the limit of x → 0. In this case, Eq. (4.10) simplifies to which shows strong dependence on the number of e-folds N . However, for p ≤ 2 we need to compute the integral numerically. Before moving on to compute the relic energy density first we need to constrain the parameter space of the theory by demanding negligible backreactions on the background dynamics. The first backreaction constraint Eq. (3.19) can be written in the form where we have used the definition of the power spectrum of curvature perturbations and =φ 2 /2M 2 Pl H 2 is the slow-roll parameter. Note that since in the backreaction condition Eq. (3.19) we use ρ T and not ρ e we have to use the value of C during inflation when the gauge field is still massless, corresponding to R = 0.
The second backreaction parameter comes from the equation of motion of the inflaton field as discussed in Eq. (3.22). If we parametrize S φ ≡ (H/φ)H 4 Σ then the second condition can be written as 4π 2 3 PΣ 1 .  . The thick red curve shows the maximum limit at hand. We can see that for all interesting values we must have p 2.2 and ξ 4.6 to avoid backreaction. We have set the number of e-folds N = 50.
In the case that we have C ∼ Σ, the condition (4.16) is more restrictive than (4.14) by a factor of . Assuming that this is the case we have plotted the logarithm of the left hand side of Eq. (4.16) for different values of p and ξ in Fig. 2. This shows that we need to have p 2.2 and ξ 4.6 to respect the backreaction conditions. The former can be understood from the approximate result of Eq. (4.13) for p > 2 which has exponential dependence on N for p > 2.
Besides the above backreaction conditions we also need to make sure that the isocurvature perturbations do not violate the observational constraints. From our expression for ρ e in (3.14) and the parameterization in (3.17), the isocurvature perturbation can be estimated as where we have related the perturbation in the number of e-folds to the curvature perturbation R through the so-called δN formalism [74][75][76][77]. Note that we also have the contribution from the variation of R. However, the variations in R is related to the fluctuations of the waterfall field which are heavy and are suppressed, so the contribution from δR in the above isocurvature perturbations is suppressed. Interestingly, Eq. (4.18) indicates that the isocurvature and the curvature perturbations are fully correlated. As a result, the ratio of the dark matter isocurvature perturbations power spectrum to the curvature perturbations power spectrum is obtained to be where we have used the recent upper bound from Planck data for the fully correlated dark matter isocurvature perturbations [20]. The strong dependence on the number of e-folds only occurs for p > 2 for which from Eq. (4.13) we obtain β = (2p − 4) 2 yielding the upper bound p 2.01 . (4.20) As a result, the constraint on p from the isocurvature perturbations is stronger than what we already have obtained from the backreaction condition. Such a strong constraint on p motivates us to focus on the values p ≤ 2 in the rest of the paper. Note that the main reason, unlike [47], that we can look into this regime of p is that the dark photon is massless during inflation and acquires mass, which can be very large, only at the end of inflation through the dynamical symmetry breaking mechanism.

Different limits
In the following we first consider two specific limits: ξ = 0 and p = 0 separately. Then we consider cases where both ξ and p can be nonzero. We see that in all cases we can obtain enough relic energy density of vector dark matter in a wide range of dark photon mass.

Limit I: Kinetic conformal coupling
In this section, we consider the kinetic conformal coupling case which corresponds to which is equivalent to the case of γ = 0 (or ξ = 0) in our analysis in the previous section. To get a non-vanishing range for tachyonic momenta we must have p ≥ 1. For the vector dark matter scenario with which we are interested in this paper, the model Eq. (5.1) is recently studied in Refs. [45][46][47]. However, here we have a different and novel scenario of vector dark matter, thanks to the symmetry breaking mechanism which allows us to have much heavier dark photons. For ξ = 0, using the identity W (0, ν, z) = √ πz in which H (1) ν (x) is the Hankel function of the first kind, the mode function (4.5) can be written in a more familiar form where as before ν = p − 1/2. Correspondingly, the expression in Eq. (4.10) can also be written more compactly as where the integration region according to Eq. (4.11) is now given by  Note that by setting ξ = 0 all dependences on polarizations disappear, i.e. both polarizations contribute equally to the dark matter energy density.
In the case with p > 2, the dominant contribution to the integral in Eq. (5.3) comes from the lower limit x min ( 1) and therefore we can use the small argument limit of the Hankel functions to the first order of approximation. On the other hand, for 1 ≤ p ≤ 2, it is not clear if we can use the small argument limit of the Hankel functions. However, we have checked numerically that in computing the integral in Eq. (5.3) for all p ≥ 1 it is sufficient to use the small argument approximation of the Hankel function (5.5) Using the above approximation in Eq. (5.3) and then performing the integration, we find the following analytic result . (5.6) We can use this to compute the relic density from Eq. (3.31). We also need to compute the t NR which depends on R and κ = p(p − 1). For instance, assuming R > 1 the system becomes non-relativistic right at the end of inflation. However, for R < 1 it will take some time for the system to become non-relativistic which depends on the minimum value of R 1/4 and R 1/2 /κ. Fig. 3 shows the relic energy density today in this limit for 1 ≤ p 2.2 (note that from the isocurvature bound we must have p 2.01) for two cases R = 100 and R = 0.01. The red horizontal line is the observed dark matter density today. We can see that assuming 10 12 T r 10 13 GeV we can get all of the dark matter in form of the vector dark matter with appropriate choice of p. Also the dashed curve in Fig. 3 is obtained from the approximate analytic form of Eq.(5.6). As a last comment note that the case p < 2 is absent in the analysis of Ref. [47] since in their setup with an exponentially small value of R (R ∼ 10 −80 ), the produced dark photons would have very large momenta and will never become non-relativistic. In our setup with R > 1, larger momenta can become non-relativistic and we can have vector dark matter even for the case p < 2.

Limit II: Chiral coupling
In this section we consider the case known as the pseudoscalar inflation with chiral interaction with the vector field and is defined by which corresponds to the case with p = 0. Here we add the assumption that |q| 1 such that ξ defined in Eq. (4.4) is approximately constant in time. Further, note that this model is very similar to axion-like coupling of inflaton to the gauge field in which where f a is a dimensionful parameter which determines the strength of the interaction (the axion decay constant). In this case we can define [30] ξ ≡φ 2f a H , (5.9) which is approximately constant up to slow-roll parameters. With this definition, our phenomenological model also includes the results of the axion-like coupling case. Indeed, vector dark matter in the latter model with the coupling given in Eq. (5.8) was recently studied in Ref. [44]. However, in our symmetry breaking scenario we can have much heavier dark photons which makes our results different than those of Ref. [44]. In Ref. [44], the authors also took into account the effects of time dependence of ξ by numerical method. Here, we assume that ξ is constant as even with this assumption we have a new scenario thanks to the symmetry breaking mechanism.
Getting back to the model described by Eq. (5.7), we can use our solution Eq. (4.5) by simply setting p = 0. In this case we can use the identity W µ,1/2 (z) = e iarg[Γ(1−µ)] e −iπµ/2 iF 0 (iµ, iz/2) + G 0 (iµ, iz/2) in which F 0 and G 0 are the regular and irregular Coulomb functions respectively, to find where Θ ≡ arg[Γ(1 + iλξ)] determines a constant phase. The above result coincides with the result of Ref. [30] up to the phase factor e iΘ which is not important for our purpose. Contrary to the case of kinetic conformal coupling where we could use the small argument relation (5.5), here we cannot expand the Coulomb functions to find an analytic expression for the energy density of the dark photons. However, as it is shown in Ref. [30], we can instead use the WKB approximation to find where we have assumed that ξ > 0 as in previous section. As we see, the produced dark photons are completely chiral dominated by the plus polarization (with our convention with ξ > 0). Following Ref. [29] and using the approximate WKB solution for the mode function, we find that where we have neglected the contribution from the minus polarization as it is completely negligible in comparison to the exponentially enhanced plus polarization. Having Eq. (5.12) at hand, we can compute the relic energy density of dark matter from Eq. (3.31). As in the previous section, we consider the cases R > 1 and R < 1 separately. For R > 1 the vector field is non-relativistic right after inflation while for R < 1 it takes some time for the energy density to become non-relativistic which depends on the value of R and κ = 2ξ. In both cases the relic energy density is given by the general expression (3.31). The relic density for dark matter versus the chirality parameter ξ for both R > 1 and R < 1 are plotted in Fig. 4. The solid curves are computed numerically from Eq. (4.10) while the dashed curves are the approximate analytic result Eq. (5.12). We can see a reasonable consistency between the numeric results except for small ξ for which the WKB approximation breaks down. Finally, note that with 10 12 T r 10 13 GeV we can generate all of the observed dark matter from this model.

General case
After considering two special cases of ξ = 0 and p = 0 in subsections 5.1 and 5.2, in this subsection we consider the more general case where both ξ and p are nonzero. In this case we have an analytic solution for the mode function Eq. (4.5) in terms of the Whittaker function. However, we cannot perform the integrals in Eqs. (4.10) and (4.17) analytically so we have to perform the integrals numerically. The models studied in 5.1 and 5.2 then can be thought of as special cases of this general model. From Eq. (3.31), for each value of (p, ξ) we are interested in the corresponding reheating temperature T r such that all observed dark matter abundance can be generated in our model. Here, we set R = 100 such that the dark photons become non-relativistic right after inflation. Other values of R can be treated similarly taking into account the effect of t NR .
The kinetic-like scenario for which p ≥ 1 is shown in the left panel of Fig. 5 where the number over each contour is the logarithm of the corresponding reheating temperature in GeV units. We can see that a wide range of viable parameters can be used with temperature between around 10 12 ∼ 10 13 GeV. Also, as we have discussed earlier, the net chirality of dark photons is controlled by the parameter ξ. In the right panel of Fig. 5 we have plotted the parity parameter where ρ ± is the relic energy density of the plus and minus polarizations respectively. One can see from the right panel of Fig. 5 that upon increasing ξ the plus polarization dominates very easily. Finally, the chiral-like scenario is shown in Fig. 6 for which, as the other case, we have plotted the logarithm of the reheating temperature in GeV unit. Once again, all dark matter can be obtained from our model assuming temperature between around 10 12 − 10 13 GeV. Also in this case  Figure 6: Chiral-like case with 0 < p < 1 and R = 100. We have plotted the logarithm of the reheating temperature, log 10 (T r /GeV), needed to generate the observed dark matter density Ω DM = 0.27.
only the plus polarization is excited so the parameter δ defined in (5.13) is nearly equal to unity. However, note that there is a non-tachyonic region associated with the small values of ξ as denoted by the grey area in Fig. 6.

Summary and conclusions
The search for the nature of dark matter is an active area of research both in particle physics and in cosmology. The popular models of dark matter are mostly in the context of beyond SM particle physics such as the WIMPs scenarios which are still waiting for the direct or indirect verifications. On the other hand, inflation is a cornerstone of early universe cosmology, a working mechanism to generate large scale structure in the observed universe. It is an interesting question if dark matter can be produced naturally in the context of inflationary dynamics. In this work we have revisited the idea of "vector dark matter" in the context of inflation with symmetry breaking.
It has been shown that the longitudinal mode of dark gauge boson during inflation can play the roles of vector dark matter after inflation [8]. Recently, it was shown to be even possible to obtain the correct relic density for the vector dark matter from the transverse modes [44][45][46][47]. The idea is to pump energy from the inflaton to the dark photons sector through some appropriate couplings to prevent the gauge field energy density from being diluted during inflation. In Ref. [44] the model with the parity violating interaction term φFF has been studied and it is shown that to have the correct dark matter relic abundance the dark photon's mass is at the order of m A ≥ 10 −6 eV. On the other hand, in Ref. [47] the authors considered the kinetic conformal coupling model and have found much lighter mass for the dark photons m A ∼ 10 −21 eV. The difficulty with these scenarios is that in one hand one needs to assume a small mass during inflation to have an efficient dark photon production. On the other hand, one needs that the dark photon's mass to be large enough so that the energy density of the dark photon becomes non-relativistic before the time of matter-radiation equality.
In this paper, we have presented a new vector dark mater scenario which relaxes these difficulties. In this model, the gauge field is massless during almost all period of inflation and therefore the dark photons can be produced efficiently. Then the gauge field acquires mass through a dynamical symmetry breaking mechanism toward the end of inflation. This dynamical mass generation mechanism for the dark photons allows us to have large mass right after inflation. In this respect, we have significantly enlarged the parameter space of the models considered in Refs. [44] and [45][46][47] by relaxing the condition of having a small mass. Moreover, we have investigated a new setup containing both models of Refs. [44] and [45][46][47] simultaneously. We have shown that in this general model, we can obtain the correct relic abundance of the vector dark matter for a large space of the model parameters.
The symmetry breaking at the end of inflation allows us to have a dark photon much heavier than or comparable to the Hubble scale at the end of inflation, corresponding respectively to R 1 and R ∼ 1. With R not exponentially smaller than unity our setup predicts that we need a reheat temperature around and lower than 10 12 GeV in order to generate the observed dark matter abundance today. This is not a very high reheat temperature. As a result, our setup predicts that the tensor to scale ratio parameter r to be very small, typically r 10 −14 , which is too small to be detected in any foreseeable cosmological observations. Of course, choosing an exponentially small value of R, our model can saturate the current upper bound r 10 −2 .
With an instantaneous reheating scenario the mass of dark photon is related to the reheat temperature via m A ∼ √ RT 2 r /M P . So with T r = 10 12 GeV we obtain m A ∼ √ R × 10 6 GeV. For example with R = 10 2 we obtain m A ∼ 10 7 GeV while for R = 10 −2 we have m A ∼ 10 5 GeV. As a result, the dark photon in our setup is typically much heavier compared to the conventional WIMPs scenarios. However, with an exponentially small value of R, much lighter dark photons can also be generated as well.
Another generic prediction of our setup is that cosmic strings are produced at the end of inflation. Using the CMB upper bound on the tension of strings we have obtained the approximate upper bound Eq. (3.34) which allows for a wide range of the value of R.
In this work we have employed a heavy complex field, like the waterfall field of hybrid inflation, to undergo symmetry breaking and generate the mass for the dark gauge field via the Higgs mechanism. However, it is possible to consider a scenario in which the inflaton field itself is complex and is charged under the dark gauge field. If the inflationary potential has global minima as in Higgs potential, then the mechanism of vector dark matter production may occur in this setup as well. We would like to come back to this question elsewhere.
Finally we comment that to simplify the analysis we have neglected the dark photon production during the (p)reheating stages. Specifically, because of the couplings f (φ) 2 and h(φ) 2 the quanta of dark photons can in principle be generated via parametric resonance. To avoid this complication we have assumed that φ−dependence in the above couplings become very weak after inflation. Otherwise, one must take into account the additional energy transfer from the oscillating inflaton field to the dark photons during preheating as well. This can be particularly important for the dark photons lighter than H e , i.e. in models with R < 1. However, for heavy dark photons with R > 1 the particle production is not efficient and one can safely neglect its contribution in the observed dark matter relic density. Moreover, the details of the efficiency of the dark photon production via parametric resonance depend on the shape of the inflaton potential around its minimum. This is an interesting question which is beyond the scope of this work but we would like to come back to this question in future.