Dark Photon Dark Matter from Cosmic Strings and Gravitational Wave Background

Dark photon dark matter may be produced by the cosmic strings in association with the dark U(1) gauge symmetry breaking. We perform three-dimensional lattice simulations of the Abelian-Higgs model and follow the evolution of cosmic strings. In particular, we simulate the case of (very) light vector boson and find that such vector bosons are efficiently produced by the collapse of small loops while the production is inefficient in the case of heavy vector boson. We calculate the spectrum of the gravitational wave background produced by the cosmic string loops for the light vector boson case and find characteristic features in the spectrum, which can serve as a probe of the dark photon dark matter scenario. In particular, we find that the current ground-based detectors may be sensitive to such gravitational wave signals and also on-going/future pulsar timing observations give stringent constraint on the dark photon dark matter scenario.


Introduction
Dark photon is one of the candidates of dark matter (DM) of the Universe [1,2].The dark photon may have string-theoretic origins [3,4] and phenomenologically it is characterized by the mass and coupling through the kinetic mixing with the Standard Model photon.There are several known mechanisms for the dark photon production processes in order for the dark photon to be DM: gravitational production [5][6][7][8][9][10], gravitational thermal scattering [11,12], production through axion couplings [13][14][15], through Higgs couplings [16,17] and cosmic strings [18]. 1 In this paper we focus on the production through the cosmic strings.
If the dark photon, A µ , is associated with the dark U(1) gauge symmetry, it is natural that the finite dark photon mass arises from the Higgs mechanism, where the dark Higgs field Φ obtains a vacuum expectation value (VEV) and the dark U(1) symmetry is spontaneously broken.If the symmetry breaking happens after inflation, cosmic string networks appear in the Universe as topological defects [24].Once formed, the cosmic strings radiate dark photons if the dark photon is sufficiently light and the dark photon abundance produced in this way can be consistent with the present DM abundance [18].
One of the main purposes of our study is to confirm that the cosmic strings efficiently emit dark photons with the use of three-dimensional lattice simulation.On this point, one should note that there are already a series of studies for the case of global strings, i.e., the axion emission from cosmic strings associated with the global U(1) symmetry breaking .An important property of the global strings is that the average string number per Hubble volume grows logarithmically with time: ξ(t) = β + α log(m r /H), which we call the logarithmic violation of the scaling law [42-44, 46, 47, 49].On the other hand, Refs.[45,48] claim the conventional scaling law: ξ(t) = const ≃ 1.Although these results have been obtained for global strings, one may expect that the production of the longitudinal vector boson should be similar to the axion thanks to the Goldstone boson equivalence theorem if the momentum is high enough.We will reveal it by an explicit calculation.
Another main purpose of this paper is to estimate the gravitational wave (GW) background spectrum in the cosmic string scenario for the dark photon DM.The estimation of GW background from the cosmic strings, in particular from cosmic string loops, have been done in many works mostly in the case of local strings [51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66].Recently GWs from global string have been calculated in connection with the recent development of the global string simulation with the improved scaling properties [67,68].The inclusion of axion or vector boson (dark photon) emission drastically affects the resulting GW spectrum, since the most energy of string loops goes to the axion or dark photon emission while only tiny fraction of the loop energy is consumed as GW emission. 2Even more complexity arises when the dark photon is light but has finite mass.In such a case, loops longer than some threshold cannot emit dark photon while shorter ones can emit it.Thus the cosmic strings in this scenario will be a mixture of "local-type" and "global-type" strings.Such nontrivial structures should be imprinted in the shape of the GW spectrum.Ref. [18] already pointed out such a possibility, but detailed calculations have not been performed.We show explicitly the GW spectrum as a probe of the dark photon DM scenario.
This paper is organized as follows.In Sec. 2 we perform a three-dimensional lattice simulation of the local strings with light dark photon.In Sec. 3 we estimate the GW background from the cosmic string loops in the dark photon DM scenario.In Sec. 4 we give conclusions and discussion.

Simulation of Abelian-Higgs string networks
In this section, we numerically show the evolution of Abelian-Higgs string networks with both local-type and near-global-type cases.First, we focus on the scaling behavior for the cosmic string length.Then, we show the emission of dark photons from the string network.

Model
Let us consider the Abelian-Higgs model in the hidden sector with the dark U(1) gauge symmetry given by the following Lagrangian, where Φ is the dark Higgs field, D µ = ∂ µ − ieA µ is the covariant derivative operator with A µ being the dark photon field, F µν = ∂ µ A ν − ∂ ν A µ is the field strength tensor, e, λ and v are respectively the gauge coupling constant, the scalar self coupling constant and the vacuum expectation value.Since |Φ| = v at the minimum of the potential, the U(1) symmetry is spontaneously broken and cosmic strings are formed in association with the symmetry braking.The equations of motion are derived as follows, Here and in what follows, we assume the flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric as the background spacetime and use the conformal time, τ , as time variable.In addition, we impose the temporal gauge, A 0 = 0, throughout this paper.The system of evolution equations can be written as where the prime denotes the derivative with respect to the conformal time, a is the scale factor, H = a ′ /a, E i ≡ F 0i and repeated indices are summed. 3Note that two of these equations (Eq.( 4) and ( 5)) determine the field dynamics and the Eq. ( 6) is the constraint equation corresponding to the Gauss's law.Here and in what follows, we assume the radiation-dominated universe, i.e. a ∝ τ and H = 1/τ and the energy density of the Abelian-Higgs system is negligible with respect to the background radiation density.Note that the dynamics depends only on the ratio λ/(2e 2 ) after redefining the field values.It corresponds to the square of the ratio between the mass of the radial component of the scalar field and the mass of the dark photon m A = √ 2ev around the symmetry-breaking vacuum.Specifically, the dark photon is much lighter than the radial mode for e 2 ≪ λ and one can expect that string configuration approaches to the global string one.In the following, we numerically show the scaling behavior and emission of dark photons in both local and near-global cases.In the following numerical analysis, we set λ = 2 throughout.

Numerical simulation
One can follow the evolution of Abelian-Higgs string network by solving the system of equations (4), (5) under the constraint equation (6).In general, the dynamics of topological defects is highly nonlinear and hence we adopt three-dimensional lattice simulation.
The Abelian-Higgs lattice simulation has been performed in the literature [71,[80][81][82][83][84][85][86][87][88][89] and we follow the formulation based on the lattice gauge theory [82].In this formalism, the scalar field is placed at grid points and the gauge field is placed at links between two grid points.Discretized covariant derivative is defined so as to satisfy the exact gauge invariance (up to the machine precision) while it coincides with the continuous one within the second order accuracy with respect to the lattice spacing.Field values are updated following Eqs.( 4),( 5) by using the Leap-frog method.The Gauss's law ( 6) is also discretized and, in this formalism, it is exactly satisfied at each grid point at each step. 4s the initial condition, we set the random value for the phase component of the complex scalar field with the uniform probability density but the radial component is fixed to be |Φ| = v everywhere.Initial random fields are generated in the Fourier space with the cutoff scale k c /a = H, above which we set zero for the field value.The time derivative of the scalar field and each component of the gauge field is initially set to be zero.Regarding the identification of cosmic string in the lattice space, we follow the method in [90] using the gauge-invariant sum of phase differences.
In order to evaluate the scaling behavior of the network of cosmic strings, we introduce two quantities.One is a mean separation between neighboring cosmic strings which is conventionally adopted in the context of the Abelian-Higgs string simulation.It has a dimension of length and defined by where V box is the physical volume of the simulation box and ℓ tot is the total physical length of cosmic strings.When the string network follows the scaling regime, one can express it as m r d sep = am r τ + b with a and b being some numerical constants and m r ≡ λ/2v.The other is defined by the ratio of the string length to the Hubble length in one Hubble volume of the universe, which is conventionally used in the axion string simulation.It is dimensionless quantity defined by with t the cosmic time.In this case, ξ = O(1) (constant) when the string network follows the scaling regime.We start the simulation with the initial conformal time m r τ = 1 and end it at a∆x = m −1 r and aL = 2H −1 for physical and fat string cases respectively (see below), with ∆x and L being respectively the comoving lattice spacing and the comoving box size.The Setup of our simulation is summarized in Table 1 and we start the simulation with m r τ = 1. 5 Note that, in addition to the "physical" string case with λ and e being constant with time, we also consider the "fat" string case where λ and e are replaced by the time-dependent quantities, λa(t) −2 and ea(t) −1 respectively, in order to keep the string core size and the vector core size     larger than the spatial resolution of the lattice simulation.We have run 5 (10) simulations for the physical (fat) string case.Fig. 1 shows the time evolution of the mean string separation, d sep , for physical (left) and fat (right) string cases.Linear fitting lines are overlaid in these figures.The latter half of data points (m r τ = 33.6-64.0 and 133-256 for physical and fat string case respectively) are used for the linear fitting lines, am r τ +b with a and b being fitting parameters.In addition, we also use the former half of data points for the physical string case with e = 0.01.The linear fitting parameters obtained by our simulations are summarized in Table 1.The obtained values are consistent with previous studies, [86] for local string and [45] for global string, corresponding respectively to e = 1 and 0.01 in our case.One can clearly see that in the physical string case with e = 0.01, there is a transition around m r τ = 30.This implies that the string with e = 0.01 can be regarded as the near-global string in the early time and it finally behaves like the local string.Square points (lower points) in Fig. 2 and 3 show the time evolution of ξ for the physical and fat string cases respectively.Because of the constant shift from the linear function in d sep (b in the above expression), ξ cannot settle down to a constant value even in the local string case, at least in our limited simulation time.The constant shift can be a numerical artifact depending on the initial condition.Then, let us remove it following the analysis in [45].Here the string length density parameter is replaced by where t 0 and τ 0 represent the constant shift with respect to the cosmic time and conformal time respectively and the tilde denotes the comoving length/volume.τ 0 can be computed by the linear fitting of d sep (namely, τ 0 = b/a).The corrected value ξ is shown by the circle points (upper points) in Fig. 2 and 3.One can see that the corrected value ξ approaches a constant value while the original value ξ increases with time. 6We also fitted the string length parameter ξ by the log function, α log(m r /H) + β in both physical and fat string cases with e = 0.01.We adopted the former half of the data for physical string case because the late-time dynamics behaves like a local string one.The data can be well fitted by the log function and the fitting parameters are summarized in Table 2.Those parameters are near the values reported in [42,43,47] for the axion (global) string case.After the Hubble  parameter becomes smaller than the dark photon mass, the dark photon emission from the loop collapse is kinematically suppressed.Then, one can expect that, even in the case with e ≪ 1, the cosmic strings behave like local strings for m A > H and the scaling violation is no longer sustained.Thus, one can expect that ξ tends to be constant after that.However, it is still controversial whether the scaling violation indeed occurs or not.In particular, the error becomes large near the end of simulation, which makes it difficult to conclude which scenario is true.Thus, one needs larger-scale (longer-time) simulation with smaller gauge coupling constant to confirm whether the scaling violation occurs or not.We leave it for future work and in what follows we take ξ as a free parameter, which takes O(1) -O(10).
In the rest of this paper, we consider only the physical string.Figs. 5 and 6 show snapshots of time slices with log(m r /H) = 3.4, 4.6, 5.3 from left to left to right for e = 1 and 0.01 respectively.One can see more abundant small loops in the case with e = 0.01 corresponding to the near-global string case especially in the late time.It is due to the efficient vector boson emission as shown below.
The loop production efficiency is quantitatively shown in Figure 7, showing the number distribution of loops with logarithmic interval of loop size (ℓ) over the Hubble length (ℓ H = H −1 ).Those data are obtained from 10 independent simulations of physical strings with 1024 3 grids.The vertical dotted line corresponds to m A ℓ = 1 and thus the loop collapse is suppressed in the right side of this line.Note that the typical initial size of the loop decoupled from the network is roughly the Hubble length, i.e. ℓ ∼ ℓ H . Hence the number of small loops cannot increase efficiently for ℓ H > m −1 A as one can clearly see in the local string case with e = 1 (top panels).This is in contrast to the near-global string case, i.e. e = 0.01 (bottom panels), where plenty of small loops are observed.In that case, smaller loops are more abundant because the simulation box has more Hubble patches and thus more initial Hubble-size loops are produced at earlier times.Note that total number of strings decreases with time in any case because the number of independent Hubble patches in a simulation box decreases due to the cosmic expansion.
Let us consider the dark photon emission from the network of cosmic strings.In the Abelian-Higgs model, the Nambu-Goldstone mode (i.e. the phase component) of the complex scalar field and the longitudinal mode of the gauge field after the symmetry breaking are related through the Gauss constraint (6). 7Therefore, in order to evaluate the energy density of the longitudinal mode, we should also sum up the Goldstone boson term.More explicitly, the energy density of the total system given by To pick up the energy density of the "longitudinal" part ρ (L) A , we must include both the E 2 i term and the angular component of the Φ kinetic energy term.As a result, we obtain In order to screen the cosmic string and extract only the contribution of the free dark photon, the factor |Φ|/v is multiplied as in the axion string case [42,49].Namely, |Φ|/v ≃ 0 near the cosmic sting and |Φ|/v = 1 far from the string core.Note that the gradient term is highly contaminated by the cosmic string especially in the near-global string case (i.e. e ≪ 1), so the inclusion the gradient term may cause the overestimate of the free dark photon abundance.Instead, factor 2 is multiplied in order to take into account the gradient term or mass term both of which would be comparable to the kinetic term.By using the Gauss constraint (6), the energy density ( 11) is also rewritten as ρ The first term, which was originally the Goldstone kinetic term, is dominant for p ≫ m A with a typical physical momentum p.It is consistent with the intuition from the Goldstone equivalence theorem.To estimate the final relic abundance of the dark photon, it is convenient to calculate the number density, n A , rather than the energy density, since n A is conserved in a comoving volume.It can be calculated by with The energy spectrum dρ A /dk can be calculated by integrating the Fourier transformation of Im(Φ * Φ ′ ) and |Φ|E i with respect to the k-space solid angle.Dark photons emitted from the collapse of loops have typically Hubble-scale momenta, k/a ∼ H. On the other hand, the emission is kinematically suppressed for H < m A .Thus, a typical energy of dark photons is E A (k) ∼ H at the emission.Since the energy density of emitted dark photons is proportional to that of the sting network and if the network follows the scaling law, one obtains the number density of newly produced dark photons at some cosmic time t as ∆n A ∼ ρ str /H ∝ 1/t and hence a 3 ∆n A ∝ t 1/2 .Therefore, such newly produced dark photons account for the dominant part of the total abundance, implying that the final dark photon abundance is determined by the number of dark photons produced at H ∼ m A .When the dark photons are continuously produced from the collapse of string loops, the number density of longitudinal dark photons at each time slice can be roughly estimated as [18] n (L) where ĒA is the mean energy of produced dark photons.
Figure 8 shows the evolution (left) and spectrum (right) of the comoving number density of longitudinal (thick-solid lines) and transverse (thin-dotted lines) dark photons.One can see that the dark photon production is not efficient for e = 1.Persistent dark photon production is observed for e = 0.01 until log(m r /H) ∼ 5 but it stops afterward.This result can be well-fitted by Eq. ( 13) with µ = πv 2 and ĒA /H = 10 (dashed-black line).The spectrum is mildly red-tilted (∝ 1/k), which explains the relative large value of ĒA /H.Note that the transverse mode is also excited around H ∼ m A , and thus additional two polarization modes equally contribute to the total dark photon abundance.It is equivalent to multiplying n (L) A by a factor 3 to estimate the total n A .There are several ways to understand this.In the unitary gauge, the longitudinal boson coupling is enhanced by an amount E A /m A due to the kinetic term normalization compared with the transverse one, whose coupling is just given by e.When E A ∼ H ∼ m A , they are comparable.In the Goldstone picture, the Goldstone mode θ couples like ∼ (∂ µ θ/v) j µ while the transverse mode couples like ∼ eA (T) µ j µ where j µ = 2 Im(Φ * ∂ µ Φ) is the U(1) current which only exists at the string core.Thus it is seen that the effective Goldstone coupling is E A /(ev) ∼ E A /m A times that of the transverse mode. 8Therefore, the production of transverse vector boson is inefficient at early epoch, but eventually it is comparable to the longitudinal one.This is in contrast to the statement in Ref. [18], where it has been claimed that the transverse mode production is negligible.
After the production stops, the number density of dark photons divided by the entropy density, n A /s, is conserved if there is not any additional entropy production in the subsequent cosmological history.Then, one can obtain the final abundance of the dark photon where we have substituted ĒA /H = 10 obtained from the simulation (see Fig. 8).If there is no scaling violation, ξ is constant and ξ = 3 -5 from our simulation with e = 0.01.If one assumes the scaling violation, ξ is given by the following fitting formula supported by the simulation of global strings in Refs.[42,47], Our simulation is performed with only a limited time range.To determine more precise values of the final dark photon DM abundance, we need more refined numerical simulations using more grid points.It is left for future work.
As a remark, we assumed the radiation-dominated Universe around H = m A in deriving Eq. ( 14).If the reheating temperature T R is low enough, i.e.T R ≲ √ m A M Pl with M Pl being the reduced Planck scale, the Universe might be inflaton oscillation dominated at H = m A .In such a case, assuming that the inflaton oscillation behaves as non-relativistic matter, the final abundance should be multiplied by a factor T R / √ m A M Pl to Eq. ( 14), up to some corrections to numerical factor, which indicates that the abundance becomes independent of m A .

Gravitational wave background from cosmic strings
In this section we derive the spectrum of stochastic GW background from cosmic strings.In the dark photon DM scenario, a big difference from the case of local strings and pure global strings is that the dark photon is very light but still massive.The short loops that are created in the early Universe promptly decay through the dark photon emission and hence the GW emission efficiency is suppressed, while long loops that are created in relatively recent epoch cannot lose their energy through the dark photon emission.Thus the resulting GW spectrum will be significantly different from both the (pure) local string and global case.
from the string loops is estimated as P GW : P L : P T ∼ Gµ : 1 : A , where G denotes the Newton constant.In the opposite case ℓ > m −1 A we just assume that loops cannot emit a vector boson.See next section for estimation of the GW spectrum.

Energy loss of loops
The energy density of the cosmic string network in the scaling regime is where the ξ parameter is given by (15).The string tension may be given by µ = 2πv where Γ GW ∼ 50 and Γ vec ∼ 65 are numerical constants [24,27,68,[91][92][93].The step function ensures that the loop cannot emit vector boson if the inverse loop size is smaller than the vector boson mass.For small loops m A ℓ i < 1, the loop length at the time t with an initial length ℓ i is given by where For estimating the loop lifetime, the GW contribution is negligible unless the symmetry breaking v is very close to the Planck scale.Thus the loop lifetime is Up to a numerical factor of O(1), τ (ℓ i ) ∼ ℓ i in this case.Since ℓ i < t i , the lifetime is shorter than the Hubble time.On the other hand, for large loops m A ℓ i > 1, the GW emission is the only energy loss process 9 and we have Thus, for ℓ i > Γ GW Gµt i , the loops are long-lived, i.e., the lifetime is longer than the Hubble time.

Loop density
To maintain the scaling solution ( 16), the infinite strings cast their energy into string loops.An energy fraction C loop ∼ 10 % of the total string energy density is considered to go into the loops within one Hubble time.The loops lose their energy due to the emission of vector bosons or GWs.As far as v ≪ M Pl and the vector boson is much lighter than the inverse loop size m A ≲ ℓ −1 , the emission of (longitudinal) vector boson is much more efficient than that into GWs.Even the Hubble-size loops decay within one Hubble time by the emission of vector bosons.Still, a small fraction of loop energy goes into GWs and also notice that loops larger than the inverse vector boson mass cannot emit vector boson and hence GW emission is the dominant source of the energy loss.Below we estimate the stochastic GW background spectrum generated from the string loop decays. 10 Let us denote the number density of string loops with a loop length ℓ per logarithmic bin of ℓ produced at the cosmic time t as [dn ℓ /d ln ℓ(t)] prod .The total energy density of the string loops, averaged over the Hubble time, is where the subscript "prod" indicates the energy/number density of loops that are produced at the time t with the Hubble time interval and the factor t+τ (ℓ) τ (ℓ) is introduced taking account of the finite lifetime of the loop.If the loop lifetime τ (ℓ) is shorter than the Hubble time t, loops exist only in a short fractional time interval τ (ℓ)/t within one Hubble time.
The loops have some distributions on the initial loop length ℓ i .There are several studies on the initial loop size, but no definite conclusions have not been reached yet.Our simulation in Sec. 2 is also not enough to determine the loop distribution for very late Universe due to the limited dynamical range of the simulation.Instead, we consider typical possibilities for the initial distribution of α.One is the monochromatic distribution, where α 0 = 0.1 will be taken.The overall normalization is determined such that ρ loop (t) = C loop ρ str (t).The other is the flat distribution in the log spacing of ℓ in the interval m −1 r < ℓ < α 0 t [42] (which we call "log-flat distribution"), where the function Θ(t; x, y) is defined as 10 GWs are also produced from infinite strings [94][95][96], which we will not pursue in this paper.
Again, the overall normalization is determined such that ρ loop (t) = C loop ρ str (t).If the vector boson emission is efficient we have τ (ℓ) ∼ ℓ and if it is kinematically suppressed, the GW emission is dominant and we have τ (ℓ) ≫ ℓ.Thus the numerical factor in the parenthesis is always O(1).
After the formation of loops at t = t i , their number density decrease as a(t) −3 until they decay at t = t i + τ (ℓ).Thus the number density of loops of length ℓ measured at t that were produced at t = t i is where the initial loop size ℓ i will evolve according to (18).The total loop number density of length ℓ at the time t is then given by the integration over the production time t i : For the case of monochromatic distribution ( 23), the loop length ℓ i and the production time t i has one-to-one correspondence as ℓ i = α 0 t i .We obtain where For the case of log-flat distribution (24), by noting that the integral over t i is dominated by the lower edge t i ≃ ℓ i /α 0 , we obtain By noting that the factor τ (ℓ) τ (ℓ)+t i is at least O(0.1) in the monochromatic distribution with α 0 ∼ O(0.1), we can see that the monochromatic and log-flat distributions do not give much differences for the purpose of orders-of-magnitude estimation.Below we use the monochromatic distribution for simplicity.

Gravitational wave spectrum
The GW power radiated from each loop with the size ℓ into the GW frequency f is given by where the spectral function S(x) takes the form where q > 1 is assumed and the overall normalization is determined so that dxS(x) = Γ GW ∼ 50.The slope q depends on the mechanism of the GW emission.Such a power low tail in the high frequency part is produced by cusps or kinks on the loop.For cusps on the loop, we have S(x) ∝ x −4/3 [51] and for kinks S(x) ∝ x −5/3 [52,56,60].Since the cusp contribution is dominant, we take q = 4/3 in the following.
The GW energy density with the frequency between [f, f + df ] produced by the loops with length between [ℓ, ℓ + dℓ] per unit time at the time t is given by The present GW frequency f 0 is given by f 0 = f (t)a(t)/a 0 and the energy density scales as ρ GW ∝ a −4 .Thus the present GW spectrum is expressed as It is convenient to express it in terms of the density parameter Fig. 9 shows the present day gravitational wave spectrum for m A = (∞, 10 −5 , 10 −15 , 10 −20 , 10 −30 , 0) eV, respectively.We have taken v = 10 15 GeV in the top panel and v = 10 14 GeV in the bottom panel.For reference, v = 10 15 (10 14 ) GeV corresponds to Gµ = 4 × 10 −8 (4 × 10 −10 ) in terms of the frequently used parametrization Gµ for the case of local string.We have assumed monochromatic distribution for the loops.For illustrative purpose, we also show the spectrum with neglecting the log dependence by hand (i.e.ξ = 0.15 and µ = 2πv 2 ) in Fig. 10.We have used the fitting formula given in Ref. [97] for the relativistic degrees of freedom g * and g * s for arbitrary redshift.
Let us explain a bit more detail.First let us see Fig. 10 where the log dependence is dropped by hand.
• The case without vector boson emission ("heavy" line in Fig. 10) is just the same as the local string case.The flat part is roughly given by ∼ √ GµΩ rad , since the fraction of loop energy density, which was created at t = t i , to the total energy density is ∼ Gµ × [a(τ (ℓ))/a(t i )] ∼ √ Gµ when the loop decays at t = τ (ℓ) and all the energy of loops go to  The present day gravitational wave spectrum for m A = (∞, 10 −5 , 10 −15 , 10 −20 , 10 −30 , 0) eV, respectively.We have taken v = 10 15 GeV in the top panel and v = 10 14 GeV in the bottom panel.The present day gravitational wave spectrum for m A = (∞, 10 −5 , 10 −15 , 10 −20 , 10 −30 , 0) eV, respectively.We have taken v = 10 15 GeV in the top panel and v = 10 14 GeV in the bottom panel.We have neglected all the log dependence by hand for illustrative purpose.
the GW emission. 12The low frequency limit of the GW spectrum is dominated by the loops with lifetime longer than the present age of the Universe, which emit GWs around z ∼ 1.If αt eq < Γ GW Gµt 0 , with t eq being the cosmic time at the matter-radiation equality, all the loops that are created at the radiation dominated era completely decays until the present age.In this case, the low frequency spectrum is determined by the loops that are created at the matter domination era and it behaves as Ω GW ∝ f .In the opposite case, αt eq > Γ GW Gµt 0 , loops that are created in the radiation dominated era can also exist in the present Universe and it gives a contribution like Ω GW ∝ f 3/2 .For even lower frequency, loops that are created in the matter-dominated era contributes, which behaves as Ω GW ∝ f .The peak of the spectrum corresponds to GWs from loops that are decaying in the present Universe.Such loops produce f −1/3 tail towards higher frequency part until it crosses the flat part mentioned above.The peak frequency is roughly given by f ∼ (Γ GW Gµt 0 ) −1 and the peak spectrum is given by Ω GW ∼ Gµ for αt eq < Γ GW Gµt 0 and Ω GW ∼ Gµα/Γ GW Ω 3/4 rad for αt eq > Γ GW Gµt 0 .
• On the other hand, for the massless limit ("massless" line in Fig. 10), the string loops are short-lived, i.e., they decay within one Hubble time.The overall normalization of the GW spectrum at the flat part is given by ∼ (Gµ) 2 Ω rad , since the most energy is consumed by the vector boson emission and only the small fraction ∼ Gµ of each loop goes to the GW emission.The lower frequency part is proportional to f −1/3 , which is a tail of the spectrum (32) sourced by the loops of the size ∼ α 0 t 0 , which are produced around z ∼ 1.
• By assuming finite vector boson mass m A , there appear a break in the GW spectrum.
For larger loops that are produced at later epochs the vector boson emission is prohibited and hence the GW spectrum in the low-frequency limit is similar to those in the massive limit.On the other hand, for smaller loops the GW emission is strongly suppressed because of the vector boson emission.Thus the high frequency part is dominated by the f −1/3 tail of the GW emission from large loops as described by the function (32).Practically, small loops that are created in earlier epoch of the Universe decay into vector bosons and do not contribute to the present GW spectrum.It is hidden by the f −1/3 tail coming from large long-lived loops that are decaying around z ∼ 1.
• The overall picture is similar for the case including the log dependence (Fig. 9).Due to the log enhancement, the overall spectrum tilt is seen and the overall normalization also becomes significantly larger than the case without log dependence, since the string tension µ and ξ parameter gets an log enhancement factor that directly increases the GW amplitude.The "massless" line of Fig. 9 is consistent with that given in Ref. [68], where the case of global strings with massless axion has been analyzed.
An interesting feature seen in Fig. 11 is that the three lines (a)-(c) converge in highfrequency limit.It is understood as follows.As explained above, the overall normalization of the GW spectrum is proportional to √ Gµ for long-lived loops.The spectral break appears at the frequency f break that corresponds to the inverse of the loop size ℓ i = α 0 t i ∼ m −1 A , which decays at τ (ℓ i ) ∼ ℓ i /(Γ GW Gµ).The GW spectrum behaves as Ω GW ∝ f −1/3 for f > f break .The break frequency is estimated as f break ∼ ℓ −1 i [a(τ (ℓ i ))/a 0 ] ∝ m A /µ assuming that the corresponding loops are produced in the radiation-dominated era.On the other hand, the dark photon DM abundance is proportional to the combination µ √ m A (see Eq. ( 14)).By requiring the consistent dark photon DM abundance, we obtain f break ∝ µ −3/2 .Therefore, we have Ω GW (f break ) ∝ f −1/3 break .This power law dependence coincides with the f −1/3 tail (32).This is the reason why three lines (a)-(c) converge.It means that the GW amplitude at high frequency (e.g, in the DECIGO/BBO frequency range) is nearly independent of the dark photon mass if its abundance is consistent with the observed DM abundance.The parameter degeneracy is solved by combining with the low frequency observation, e.g. by SKA.
It is seen that v ≫ 10 12 GeV or m A ≪ 10 −5 eV is already excluded by the pulsar timing constraints if the log dependence of µ and ξ are taken into account.Also it should be noticed that the on-going ground-based experiments, LIGO/Virgo, may be able to detect GWs independently of the dark photon mass.However, one should note that there may be several orders-of-magnitude uncertainties on the GW spectrum, especially from the log dependence of the string tension and ξ parameter and also the loop abundance and such uncertainties significantly affect the resulting constraint on the symmetry breaking scale v and the dark photon mass m A .For example, if the GW spectrum is reduced by one order of magnitude due to these uncertainties, m A ∼ 10 −10 eV may be allowed, as it is clearly seen in the bottom panel of Fig. 11.Thus we need more accurate theoretical calculations to derive precise constraint.In any case, future space laser-interferometers such as LISA and DECIGO/BBO, ground-based detectors such as ET and also the radio wave observations of pulsars by SKA have good chances to detect GWs from cosmic strings in the dark photon DM scenario.

Conclusions and discussion
Dark photon is one of the DM candidates and it is naturally produced by the cosmic string dynamics in association with the dark U(1) symmetry breaking.For a consistent DM abundance, dark photon mass may be very light and it implies that the dark gauge coupling constant should be very small.We first presented a field-theoretic simulation of cosmic strings with such a small gauge coupling.We confirmed that, within a precision of lattice simulation, the string dynamics is consistent with the naive expectation that the vector boson emission (or the Goldstone boson emission) is efficient for small loops ℓ ≲ m −1 A and the emission is suppressed otherwise.Based on this observation, we calculated the stochastic GW background produced by the cosmic strings.We found that the dark photon DM scenario from cosmic strings predict sizable amount of stochastic GW background, which may already be excluded by pulsar timing observations for v ≫ 10 12 GeV or m A ≪ 10 −5 eV.The constraint on the string tension is stronger than the conventional local string case due to the logarithmic enhancement of the string tension and ξ parameter.However, we need more precise study on this issue considering the large scale separation between v and m A or the Hubble parameter, which is far beyond the reach of numerical simulation, since O(10) change in the GW spectrum would lead to orders-of-magnitude change to the constraint on v or m A .Some notes are in order.In calculating the GW spectrum, we made a simple assumption that loops longer than the inverse vector boson mass m −1 A cannot emit vector bosons.It may be a too simplified assumption since small scale irregularities may be present on long loops that allow to emit vector boson and long loops may lose their energy through the vector boson emission.Still we expect that the overall picture does not change drastically, since the energy of the order of µℓ remains for the loop length ℓ even if such irregularities are smoothed out by the vector boson emission, and eventually it should be compensated by the GW emission.Another general remark is that the f −1/3 tail of the GW spectrum represented by (32) is quite important in broad frequency range as seen in figures.This power law behavior assumes that the cusp on the loop dominates the GW emission.While it is natural to consider that cusps take important roles since cusps are necessarily formed periodically on the string loop according to the Nambu-Goto action [24], still it requires more investigation whether it remains true in the case of our interest.These are left for a future work.

Figure 1 :
Figure 1: Time evolution of the mean separation of cosmic strings for physical (left) and fat (right) string cases.The red and blue lines correspond to e = 1 and e = 0.01 respectively.Shaded region shows the 1-σ error.Data points are fitted by linear functions shown by thin dashed lines.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Square points (lower points) show time evolution of ξ for the physical strings with e = 1 (left) and e = 0.01 (right).Circle points (upper points) correspond to the corrected value ξ.Shaded region shows the 1-σ error.

Figure 7 :
Figure 7: The number distribution of loops with respect to the loop size for e = 1 (top), 0.1 (middle), 0.01 (bottom).In the horizontal axis, ℓ H ≡ H −1 is the Hubble length.Time slices correspond to log(m r /H) = 4.6, 4.9, 5.3 from left to right.The vertical dashed line corresponds to m A ℓ = 1.

Figure 8 :
Figure 8: Time evolution of the dark photon number density in a comoving box (top-left), its spectrum at log(m r /H) = 8.3 (m r τ = 64) (top-right) with λ = 2 (m r = v), e = 1 (red), 0.01 (blue) and the evolution of the spectrum of the longitudinal mode (bottom) where time evolves from left to right (m r τ = 1 -64).We ran 5 simulations and the shaded region represents 1-σ error.(The error bars are negligibly small in the left panel.)Thick solid lines and thin dotted lines correspond respectively to the longitudinal and transverse mode.In the top-left figure, The black-dashed line represents the analytic estimation given by Eq. (13) with ĒA /H = 10, µ = πv 2 with numerically obtained values of ξ.Magenta line represents the result with twice larger simulation box size.In the right panel, horizontal axis is normalized by the horizon scale and each vertical dashed line represents k/a = m A .

Figure 10 :
Figure 10:The present day gravitational wave spectrum for m A = (∞, 10 −5 , 10 −15 , 10 −20 , 10 −30 , 0) eV, respectively.We have taken v = 10 15 GeV in the top panel and v = 10 14 GeV in the bottom panel.We have neglected all the log dependence by hand for illustrative purpose.

Figure 11 :
Figure 11: (Top) GW spectrum from cosmic strings in the dark photon DM scenario, compared with current and future sensitivities of various experiments.Three theoretical lines correspond to benchmark points in which correct dark photon DM abundance is obtained: (a) (v, m A ) = (10 14 GeV, 10 −13 eV), (b) (10 13 GeV, 10 −9 eV), (c) (10 12 GeV, 10 −4 eV).(Bottom) The same as the top panel, but dropping the log dependence of µ and ξ by hand.

Table 1 :
Simulation setup and linear fitting parameters of the mean string separation in terms of the conformal time, defined by m r d sep = am r τ + b.

Table 2 :
Log fitting parameters of the string length parameter, defined by ξ = α log(m r /H)+ β.
2 log(m r /M ) where M ≡ max[H, m A ]. Within Hubble time, cosmic strings intersect with each other, forming loops which eventually decays emitting GWs or light vector bosons.The energy loss rate of a loop is given by