Inflation with two-form field: the production of primordial black holes and gravitational waves

Antisymmetric tensor field (two-form field) is a ubiquitous component in string theory and generally couples to the scalar sector through its kinetic term. In this paper, we propose a cosmological scenario that the particle production of two-form field, which is triggered by the background motion of the coupled inflaton field, occurs at the intermediate stage of inflation and generates the sizable amount of primordial black holes as dark matter after inflation. We also compute the secondary gravitational waves sourced by the curvature perturbation and show that the resultant power spectra are testable with the future space-based laser interferometers.


Introduction
What is the nature of dark matter?Despite nearly a century of exploration, most of its properties are still wrapped in mystery.In fact, the potential mass scale of dark matter ranges over 90 orders of magnitudes from 10 −31 GeV to 10 60 GeV [1].For decades, dark matter searches have focused on particle models for dark matter, especially weakly interacting massive particles (WIMPs).Recently, however, the non-detection of WIMPs provides strong motivation for the search for alternative dark matter candidates.
Among various dark matter candidates, primordial black holes (PBHs) are well motivated, requiring no extensions to the standard model of particle physics, and have received a lot of attention owing to the recent detection of gravitational waves (GWs) from merging black hole binaries [2][3][4].PBHs can form by the gravitational collapse of local over-dense regions in the early universe [5][6][7].While the mass of PBHs can span many orders of magnitude in principle, the non-detection of PBHs constrains their abundance in many mass ranges.PBHs which formed with a mass ≲ 10 14 g would have evaporated by today, and could not contribute to the dark matter abundance today.For light PBHs (≲ 10 17 g), the evaporation of PBHs through Hawking radiation would affect the extragalactic/galactic γ-ray background [8][9][10][11].On the other hand, the abundance of non-evaporating heavy PBHs can be tested by their gravitational effects.The measurements of gravitational microlensing events limit the abundance of PBHs with mass range around 10 22 g − 10 35 g [12][13][14][15].Also, stellar-mass PBHs (≳ 10 33 g) would emit high energy photons via gas accretion and could modify the thermal history of universe [16,17].Hence, PBHs can be the dominant component of dark matter in the intermediate mass window, 10 17 g − 10 22 g 1 .
The abundance of PBHs is related to the amplitude of the primordial power spectrum of the curvature perturbation originating from cosmic inflation.In order for PBHs to compose a significant fraction of dark matter, it is required that a large amplitude of the curvature perturbation is produced on scales much smaller than the measurement scale of cosmic microwave background (CMB) anisotropies.In the simplest class of inflationary models, however, the predicted spectral shape is almost scale-invariant.From the measurement of CMB anisotropies, the amplitude of curvature perturbations is expected to be too small to predict a sizable amount of PBHs.To go beyond this naive expectation, many mechanisms for generating large curvature perturbation on small scales have been proposed, including the running mass model [26], axion inflation [27][28][29][30][31], a waterfall transition during hybrid inflation [32][33][34][35][36][37], a quartic action during inflation and a variable sound speed [38], an inflation coupled with the Gauss-Bonnet term [39], amongst many others.
Among these works, the mechanism of generating curvature perturbations via the particle production of coupled matter sectors has been intensively studied.Such a representative sector is the form fields.A vector field (dubbed gauge field or one-form field) and an antisymmetric tensor field (two-form field) are predicted by string theory, which are naturally coupled to the scalar sector in the low-energy effective action [40,41].In this framework, the form field can experience an instability caused by the motion of coupled scalar field and can source the coupled fluctuations on the super-horizon scales during the inflationary epoch, which finally leads to a rich phenomenology such as the generation of PBHs or the accompanied scale-dependent GWs.[28,31,[42][43][44][45][46][47][48][49].It is the aim of this paper to explore the particle production of the two-form field and its observable consequences.
In recent years, the phenomenology of particle production of two-form field has also been developed in parallel with studies of gauge field production.If the generation of the two-form field is on large scales, then a coherent mode of two-form field would break isotropy of universe due to its direction dependence [50][51][52].It would predict a statistically-anisotropic power spectrum [53], whose spectral anisotropies are different from the case of inflationary model with U (1) gauge field [54].On the other hand, the possibility of the generation of two-form field on small scales has been overlooked.In this paper, we extend our previous work on a U (1) gauge field [48] to the particle production of two-form field.We show that the two-form field can be amplified at an intermediate stage of inflation and can predict a sizable amount of PBHs as dark matter.
This paper is organized as follows.In Section 2, we build up the model of the 2-form field kinetically coupled to the inflaton.In Section 3, we calculate the curvature power spectrum sourced by two-form field.We estimate the PBH abundance in Section 4. In Section 5, we also discuss the generation of tensor perturbations sourced by two-form field and estimate the amount of secondary GWs after inflation.We finally summarize our study in Section 6.Throughout this paper, we will set the units ℏ = c = 1 unless otherwise specified. 1Although there have been other constraints discussed in this mass range such as the femtolensing events of γ-ray bursts [18], the dynamical capture of PBHs by stars [19][20][21][22], and the ignition of white dwarfs by PBHs [23].However, recent studies have revisited these constraints and claimed that this mass window remains open for PBHs as dark matter candidates [24,25].

Model Setup
In this paper, we consider the following inflationary model where an antisymmetric two-form field B µν is kinetically coupled to the inflaton: where is the antisymmetric field strength of the two-form field and I(φ) is a non-trivial kinetic function deriving from the compactified volume of extra dimensions in the framework of string theory [40,41].Hereafter, we assume that two-form field induces negligible backreaction to the background inflaton dynamics, which we will confirm in Section 3.For a gauge condition of two-form field, we can take with the gauge transformation δ g B µν = ∂ µ ξ ν − ∂ ν ξ µ and the redundant degree of freedom ξ µ → ξ µ + ∂ µ χ.One can show that B 0i is a non-dynamical variable and can be integrated out in the action.At a quadratic order in perturbation, B 0i could couple only to B ij , because B ij is assumed to have no background components.Due to the gauge conditions (2.2), however, those quadratic interactions between B 0i and B ij vanish and hence B 0i gives a higher-order contribution in perturbations.Therefore, hereafter we neglect B 0i in our analysis since its contribution becomes relevant only for the trispectrum.Throughout this paper, we assume that B µν does not have its homogeneous component that leads to the breaking of the isotropy of space [50][51][52].Hence, we consider the spatiallyflat FLRW metric, ds 2 = −dt 2 + a(t) 2 dx 2 = a(τ ) 2 (−dτ 2 + dx 2 ), and the perturbation of two-form field at zeroth-order level: B µν (t, x) = δB µν (t, x).We decompose δB ij in Fourier space as where ϵ ij is the antisymmetric tensor with k = (sin θ k cos φ k, sin θ k sin φ k, cos θ k) represented in terms of the polar coordinates.Then, it satisfies the following relationships: Then, the equation of motion (EoM) for δB k is given by where we ignored the inflaton perturbation δφ(t, x) = φ(t, x) − φ(t) in the kinetic function Ī ≡ I( φ(t)) because it gives a higher-order contribution in perturbation, and introduced a dimensionless conformal time variable, x ≡ −kτ .For a varying Ī, the EoM is modified and the two-form field can be produced.To avoid the strong coupling problem, we assume Ī decreases with time.Here, we define the logarithmic decay rate of the kinetic function with respect to the number of e-foldings dN = −Hdt.As we will see later, Ī decreases and the decay rate is positive during inflation, n(t) > 0.
To find useful analytic expressions, we first consider the case where n(t) is varying slowly enough to neglect its time-variation: n ≃ const.Then (2.6) is rewritten as With the Bunch-Davies initial condition, the solution is given by the Hankel function of the first kind.The "electromagnetic-like" components of the two-from field, E k ≡ Īδ Ḃk /a2 and M k ≡ k ĪδB k /a 3 , are computed as where we took the super-horizon limit in right proportional relations.The contribution to the electric energy density ρ E ≡ I 2 Ḃ2 ij /(4a 4 ) from each Fourier mode is proportional to |E k | 2 .Thus, the two-form field increases its energy density on the super-horizon scale if n > 1.Even if n(t) varies in time and the solution would not be given analytically, the particle production of two-form field occurs on the super-horizon regime for n ≳ 1 as we will see in Section 3. On the other hand, the gradient (magnetic) energy density is suppressed by a factor of x 2 in comparison with electric energy density.Therefore, hereafter we ignore the contribution from the gradient terms.
Regarding the kinetic function, we consider the following functional form [48] where B 1 , B 2 , and c 1 are model parameters which are all positive.In addition to a conventional exponential function, we also introduce another constant term expected to arise from a string-loop modification in powers of a dilaton-dependent coupling constant [55].We also assume this constant term plays a role of stabilizing the dilatonic field and controlling the generation of two-form field at around the end of inflation.For the infaton potential, we adopt the Starobinsky model consistent with the current (n s , r) constraint in CMB observation 2 and solve the background system of inflaton: where we have ignored the backreaction effect from the two-form field and we will discuss the validity of this in the next section.As discussed in Ref. [48], one can find approximated analytic expressions for the e-folds N (t) and the decay rate n(t) under the slow-roll approximation.However, they do not enable us to compute the power spectrum of the two-form field with sufficient accuracy, and we do not present them here.We set the following values of model parameters

Generation of scalar mode
In this section, we compute the scalar mode induced by the particle production of the twoform field via the interaction Lagrangian, L int = −I(φ) 2 H µνρ H µνρ /12.At leading order, it is written as where we ignored the contribution from the gradient term of two-form field because it is always suppressed on the super-horizon scale in comparison with the electric-like term (see (2.9) and (2.10)).We quantize the fluctuation of the inflaton field and decompose it in Fourier mode as δφ(t, x) = dk (2π) 3 δφ k e ik•x . (3. 2) The equation of motion for inflaton in the slow-roll regime is approximately given by ) The solution can be split into the vacuum mode and the sourced mode as By using the Green function method, the solution of sourced mode is given by with the retarded Green function G R (x, y) ≡ −Θ(y − x)(x 3 − y 3 )/(3xy) obtained by solving the homogeneous solution of (3.4) in de Sitter approximation.Since the two inflaton modes are statistically-independent, the power spectrum is written as Without loss of generality, we can choose k = (0, 0, 1).In this case, the following relation of the antisymmetric tensor holds: Combining (3.4), (3.6) and (3.7), we obtain where we approximate the lower bound of the time integral by |τ min | = min(1/p, 1/|k − p|), because the two-form field can grow only after the horizon crossing and we focus on the contribution from the super-horizon modes.
To perform the time and momentum integrals in (3.9) separately, we find an analytical formula for E k (τ ).The electric mode function of the two-form field is well approximated by the following Gaussian fitting function on the scales where the particle production is relevant:  where E peak (k) is a maximum amplitude of E k at τ = τ peak and σ characterizes the time period for which E k stays around the peak amplitude.We simply use the numerical result for E peak (k).τ peak corresponds to the time when n(t) falls below unity at N CMB − N ≈ 49 where the growth of the two-form field stops.σ is determined by n(t) which is governed by the background dynamics, and hence σ does not depend on k.Fig. 2 shows the time evolution of the electric mode function crossing the horizon at N k ≈ 25 , and compare it with the fitting solution.One can see that (3.10) fits well the numerical solution in the time domain.We also plot the time evolution of the dimensionless power spectrum of the electric field P E (k) ≡ k 3 |E k | 2 /(2π 2 ) in Fig. 3 and find that this fitting function agrees well with all the relevant modes at around the time when the particle production occurs significantly.
Using (3.10), we rewrite (3.9) in the late time limit as where p * ≡ p/k , |k − p| * ≡ |k − p|/k and we have used the slow-roll approximation Īφ / Ī ≃ n/(M Pl √ 2ϵ H ). We notice that the time integration of F can be extended to infinity because its integral receives its support almost around the peak τ = τ peak .We numerically integrate the time and momentum integrals in (3.11) and obtain P δφ,s .Now we convert the obtained P δφ,s into the power spectrum of the curvature perturbation using the relation on the flat-slicing, ζ = −Hδφ/ φ.Corresponding to (3.7), the curvature power spectrum is also separated into two contributions (3.13)The first term Pl ϵ H ) is the power spectrum of the vacuum mode normalized as P ζ,v (k CMB ) ≃ 2.1 × 10 −9 on CMB scale.On the other hand, the second term denotes that of the sourced mode which has a peak around k = k p ≃ 10 13 Mpc −1 , inheriting a similar feature as that seen in P E (k).Around its peak, the numerically computed P ζ,s (k) is well described by the following fitting function: In the plot of Fig. 4, we set In Fig. 4, one can confirm the agreement between the numerical and the fitting results.We will also use the same values of the parameters as (3.15) to estimate the abundance of PBHs in the next section.Before closing this section, we estimate the backreaction effect of the two-form field on our background dynamics.The electric energy density of the two-form field is included in the Friedmann equation and the equation of motion for inflaton as We can roughly estimate the magnitude of ⟨ρ E ⟩ as where ∆N = O(1) is the logarithmic width of the power spectrum at around the peak k = k peak .Therefore, the backreaction of the energy density of two-form field is negligible when the following conditions holds: ) where r v is a tensor-to-scalar ratio of vacuum modes.Therefore, (3.20) is automatically satisfied if (3.21) holds.The amount of r.h.s. in (3.21) is roughly 10 7 around the spectral peak, which is greater than l.h.s. by an order of magnitude (see Figure 3).Thus, we can safely ignore the backreaction of the two-form field at the background level.

Estimate of primordial black hole mass function
In this section we will discuss the generation of PBHs in our model.We will consider the formation of PBHs from the collapse of large-amplitude scalar perturbations during the radiation dominated epoch.When a perturbation, which is initially super-horizon re-enters the horizon, it can quickly collapse to form a PBH if above some threshold amplitude.In order to estimate the abundance, we use the Press-Schechter-type (PS) formalism [57].A more accurate constraint on the power spectrum can be obtained by using peaks theory [58] (as well as modifications to peaks theory for application to PBHs [59,60]).However, using the Press-Schechter formalism is straightforwards to adapt to non-Gaussian distributionswhich needs to be accounted for since the abundance of PBHs depends strongly on any non-Gaussianity [33,[61][62][63][64][65][66][67][68][69][70].The PS approach provides constraints on the power spectrum which are accurate to O(10%) when compared to peaks theory calculations [71], which is sufficient for our purposes here.
If the amplitude of a density perturbation is large enough, it will collapse to form a PBH upon horizon entry.In order to determine whether a perturbation will collapse, the appropriate criterion to use is the compaction function [72][73][74] defined as where δM (x, R) is the mass excess within a sphere of areal radius R centred at location x, where we have used G = 1 in this section for convenience.We note that this is equivalent to the smoothed density contrast evaluated at horizon entry if the transfer function is neglected.On super-horizon scales, the time-dependence of δM and R cancel, making C a time-independent quantity.On super-horizon scales, at the centre of spherically symmetric peaks, the compaction function is related to the curvature perturbation ζ as [73] C where the prime denotes a spatial derivative (which here only includes the radial term due to the spherical symmetry).We have assumed spherical symmetry because the large, rare peaks, from which PBHs form, are expected to be close to spherically symmetric [58].The notation ζ ′ (r) then represents the value of ζ ′ at a distance r from the centre of a peak located at x.Typically, ζ is expected to follow a Gaussian distribution.However, in the model discussed here, the large contribution to the power spectrum P ζ,s is highly non-Gaussian and is expected to follow a χ-squared distribution, allowing us to write: where χ is a Gaussian variable -which will be helpful because the derivative of a Gaussiandistributed variable remains Gaussian.The compaction is then related to χ as where the component −rχ ′ (r) is what one would expect from smoothing the second derivative of 1 3 χ with a top-hat smoothing function of radius as below r [75].
where the subscript r indicates a smoothing at scale r, and the subscript 2 denotes the second spatial derivative.The factor of r 2 is included as the compaction is calculated as 1/r, whilst the smoothing function (below) includes a term 1/r 3 .The top hat smoothing function W (x, r) is given by In principle, χ(r) and χ 2,r are correlated Gaussian variables, and their respective PDFs should be integrated over to find the PDF of C(x, r).However, in the high peak limit relevant for PBHs, the PDF of χ(r) given a value of χ 2,r is well approximated by a Dirac-delta function.For typical PBH-forming perturbations, we can make the substitution χ(r m ) ≈ −0.756χ 2,r without significant loss in accuracy, where r m is the smoothing scale at which the compaction peaks [76] (see also appendix A for a brief discussion).This allows us to write the compaction as a function of a single variable, where in the second equality, for convenience we have rewritten χ 2 2,r = C 1 (to represent the linear component of the compaction).We here note that, since this follows a quadratic relationship, there is a maximum value for the compaction C max = 2/3 at C 1,max = 1.982.Perturbations with C 1 < C 1,max are referred to as type I perturbations, and perturbations with C 1 > C 1,max are referred to as type II perturbations [73].We will not consider further type II perturbations for 2 reasons: firstly, because the abundance of such perturbations is exponentially suppressed and has a negligible impact on PBH abundance; and secondly, the evolution of type II perturbations cannot be simulated simply using density perturbations and is thus not well understood.
C 1 follows a χ-squared distribution given by The variance is calculated from the power spectrum P ζ,s as a function of the smoothing scale r as The variance 2σ 2 is plotted in Fig. 5 as a function of r, for the power spectrum given in equation (3.14), with parameter choices as in equation (3.15).In order to determine the abundance of PBHs, we need to know the mass of PBH which will form from a perturbation, and this depends on both the scale r and amplitude C of the perturbation as [75,77,78] where M H (r) is the horizon mass of the unperturbed background when the horizon scale is equal to the perturbation scale, (aH) −1 = r, C c ≈ 0.5 is the critical value for PBH formation (with corresponding value for the linear component C 1,c ≈ 0.992).The value of K depends mildly on the specific profile shape of the perturbation collapsing, but we will here simply take K = 4, which is valid for typical profile shapes, and γ = 0.36 during radiation domination [75].Substituting equation (4.8) to write the PBH mass in terms of C 1 , we obtain Whilst there is no upper limit for the horizon mass which can form a PBH of given mass, there is a lower limit, given by Figure 5.The variance of the linear component of the compaction, ⟨C 2 1 ⟩, is plotted against the smoothing scale r, for the power spectrum given in equation (3.14).PBH formation peaks strongly at the scale where the variance of perturbations is largest.
Solving to write C 1 as a function of the mass, which will be useful later, gives where we have kept only the solution corresponding to type 1 perturbations.The horizon mass associated with a particular smoothing scale r is given by [71] M H (r) ≈ 0.5 g * 10.75 where g * is the number of relativistic degrees of freedom (changes in g * have a minimal effect due to the power of 1/6 and will be neglected here).This allows us to relate the smoothing scale r to a specific horizon mass M H -and therefore calculate the variance as a function of horizon mass, 2σ 2 = 2σ 2 (M H ). The fraction of the universe collapsing to form PBHs at the time of formation (we consider this as the time of horizon entry [74]) can be calculated by following a PS formalism by integrating over the range of perturbations which form PBHs: where the scale-dependence is encoded in the horizon mass M H . Assuming that the universe evolves strictly as radiation-dominated up until the time of matter-radiation equality, we can calculate the density parameter for PBHs at matter-radiation equality: where M eq = 2.8 × 10 17 M ⊙ is the horizon mass at matter-radiation equality (using the same parameter choices as [79]).The term (M eq /M H ) 1/2 is included to account for the redshift of PBH energy density (which evolves as matter) during radiation domination.Solving the integral numerically, we find that PBHs form with the correct abundance to make up the entirety of dark matter when the amplitude of the power spectrum is A ≃ 3.2 × 10 −4 , for the power spectrum given in equation (3.14).Finally, we define the mass function of PBHs as where f (M PBH ) is normalized to integrate to unity if PBHs make up the entirety of dark matter.Expressing the integral in equation (4.16) in terms of M PBH using equation (4.14) allows us to write the final expression PBH mass function The PBH mass function f (M PBH ) is shown in Figure 6, for which PBHs make up the entirety of dark matter, and which evades constraints on the PBH abundance.

Generation of tensor modes
In this section, we compute the generation of tensor modes h ij defined as in our model.There are apparently two significant contributions to the tensor mode.One is the perturbation of the two-form field itself, and the other is the curvature perturbation sourced by the two-form field.We will consider them in order.

Primordial tensor modes
We first calculate the primordial tensor power spectrum directly sourced by the two-form field.The EoM for h ij is given by where Π lm ij is a projection operator into transverse and traceless components.We decompose h ij into Fourier modes 4  The green region around M < O(10 17 )g is excluded by the extragalactic gamma ray [8] and the orange region O(10 22 )g < M is excluded by the Subaru/HSC [12].We note that the constraints plotted are the constraints on a monochromatic PBH mass function (whilst the black line the extended mass function predicted by the model presented here) -and so are not completely applicable to our mass function.However, the constraints for extended mass functions are similar to the constraints for monochromatic mass functions. (5.4) Regarding the calculation of the polarization tensors, we use the following identities Remarkably, they vanish in the integration over ϕ p.Therefore, as has been shown in the previous studies [51,53], the two-form field does not directly produce the primordial tensor modes at leading order5 .

Induced tensor modes
Next, we compute the tensor mode induced by the scalar modes after inflation by following a method developed in the previous works [80][81][82][83][84][85][86][87][88] (see e.g.[89] for detailed review).We take the conformal Newtonian gauge: where we neglected the vector perturbations and assumed that the two scalar perturbations Φ and Ψ satisfy the condition of no anisotropic pressure: Φ = Ψ.The density power spectrum has the peak at k ∼ 10 13 Mpc −1 in our setup, and such modes reenter the horizon during the radiation-dominated era, τ < τ eq .Then, the equation of motion for tensor mode is given by with H ≡ aH, and we obtain the solution of the induced tensor mode in the momentum space h s k,i (τ ): S k (τ ) = e s ij ( k) To evaluate (5.12), we decompose Ψ k (τ ) into the primordial field ψ k and the transfer function Ψ(kτ ) during the radiation-dominated era: Using ψ k = −2ζ k /3 at the radiation-dominated era, we can obtain the power spectrum of the tensor mode P h (k, τ ) as where I is given by the time integral of transfer functions (its explicit form is given in Appendix B).Then, we evaluate the current logarithmic energy density of GWs where ρ c = 3M 2 Pl H 2 0 is the critical density of present universe.In terms of the entropy conservation law, it is found to be [83] g * ,c 106.75 where Ω r,0 is the radiation density parameter at present, h ≃ 0.68 is the dimensionless Hubble parameter, τ c is the time when the production of GWs finished, and g * ,c is the relativistic degree of freedom at τ c .The over line on P h (k, τ c ) denotes the average of the oscillation of I in (5.15).The averaged product of two I's can be approximately decomposed as the following analytical formula: where K i is obtained as [86] K (more explanations on how to obtain this formula is written in Appendix B).Now, let us compute the 4-point correlation function ⟨ζ 4 ⟩ in (5.15).The explicit form of ζ with a momentum mode k is given by the sourced scalar mode (3.6) evaluated at the time of inflation end τ = τ end : where we have used (2.7) and the slow-roll approximation.It is important to notice that it cannot be merely replaced by ⟨ζ 2 ⟩ 2 since ζ is non-Gaussian: ζ ∼ δBδB.In particular, we need to calculate the 8-point correlation function of primordial electric component of two-form field in ⟨ζ 4 ⟩: As shown in the previous studies [45,[90][91][92], these non-Gaussian sources contribute to the generation of induced GWs via three distinct diagrams shown in Figure 7.In the left diagram (called "Reducible" diagram), the contraction is taken as and is therefore simply replaced by the sourced power spectrum of curvature perturbation P ζ .To perform this, let us calculate the contribution to the reducible diagram.Using the slow-roll approximation, it leads to where we have used G R (τ, τ i ) ≡ Θ(τ − τ i )(τ 3 − τ 3 i )/(3τ τ i ) and defined (5.24) The factor 8 comes from the symmetry.Therefore, we can rewrite (5.23) as the simplest form (5.25)This is the same expression of the power spectrum of induced GWs as that from the standard scenario in which the curvature perturbation is Gaussian [80][81][82][83][84][85][86][87][88].However, we notice that the relationship of abundance between PBH mass function and the curvature power spectrum is different because of the non-Gaussianity of ζ in our scenario.In the other two contributions ("Planar" and "Non-Planar" diagrams), δBs are contracted across more than two ζs, and hence we cannot replace the contracted ζ 2 with P ζ .The contractions of two-form fields in Planar diagram lead to where × ϵ op ( p + q − p 1 )ϵ op (− (p − p 1 ))ϵ qr (− (p + q − p 1 ))ϵ qr (− (k − p + p 1 )) . (5.27) And that in Non-Planar diagram is given by where ϵ n (k, p, q, p 1 ) ≡ ϵ ij ( p1 )ϵ ij ( p − p 1 )ϵ kl ( k − p + q + p 1 )ϵ kl (− (q + p 1 )) × ϵ op (− p1 )ϵ op ( q + p 1 )ϵ qr (− (k − p + q + p 1 ))ϵ qr (− (p − p 1 )) . (5.29) These contributions have been computed numerically in the previous work in the case of not two-form field but U(1) gauge field [45] or the local-type non-Gaussian field [92] and it was found that the Reducible diagram and Planar diagram had the same order of contributions to the spectrum, whereas that from Non-Planar diagram is suppressed compared to the other two diagrams 6 .Therefore, in this work, we numerically estimate the contributions only from the Reducible diagram and the Planar diagram.
The computation of the Reducible diagram is easily performed since we can factorize the two loops of B-B into the power spectrum of ζ.The formula of reducible diagram is written in the same way to that of induced GWs without non-Gaussianity [80][81][82][83][84][85][86][87][88] Ω GW (k)h 2 = 0.39 g * ,c 106.75 where we followed expressions in [48].On the other hand, the Planar diagram has three loops of momentum integrals which cannot be fuctorized into P ζ , so that the simple analytical method cannot be used.In  8.Both the Reducible and the Planar type GWs have similar contributions, while the Planar type has a peak at smaller frequency than Reducible type.This relation is consistent with the previous study on the axion-U (1) inflationary model [45] 7 .The GWs around this frequency is potentially observable with the future interferometer experiments like LISA, DECIGO, and BBO, while the strength of GWs depends on the uncertainty on the PBH formation rate and the momentum dependence of the factor F (3.12).

Conclusion
In this work, we have suggested an inflationary model where the inflaton is kinetically coupled to a two-form field, and developed a possibility of particle production of the two-form field on intermediate scales during inflation.Depending on the configuration of the kinetic function, Figure 8.The power spectrum of induced GWs is shown.We evaluate the Reducible-type GWs (dotted black line) in (5.30), the Planar-type GWs (dashed black line) in (C.4), and a sum of them (solid black line).We also show the future sensitivity curves of interferometers as orange lines; LISA (solid) [93], DECIGO (dotted), and BBO (dashed), respectively [94]. the two-form field is enhanced at a later stage of inflation and amplifies the coupled curvature perturbation at second-order level.Then, we showed that the sourced curvature perturbation can explain the formation of PBHs in sufficient abundance to act as dark matter.This production mechanism would be similar to that from the model of U (1) gauge field coupled with a dilatonic field [48] or an axionic field [44,45].Regarding tensor modes, however, the two-form field does not directly produce primordial tensor modes because of its antisymmetric feature 8 .Therefore, the absence of primordial GWs could be a unique prediction from twoform field model.On the other hand, our scenario predicts the associated induced GWs after inflation.Since the primordial curvature perturbation is sourced by the second-order of two-form field, three diagrammatic contributions appear in the spectrum of induced GWs: Reducible, Planar, and Non-Planar diagram.We computed the contributions from Reducible and Planar diagrams and found that they prodive a comparable amount of induced tensor spectrum with different frequency regime, which modifies the relationship of spectral peaks between PBH mass function and induced GWs from the standard prediction in Gaussian field model.We also showed that the resultant spectraum is potentially testable with future laser interferometers.
It is interesting to consider the possibility that this model can also explain a recent NanoGrav event [95].In our setup, however, the generation of sourced curvature perturbation on that scale would be challenging.The reason is as follows.To provide NanoGrav event, we would need the particle production of two-form field occurring on scales much greater than scales corresponding to a few number of e-foldings before inflation ends.However, on such scales the slow-roll condition is sufficiently met and hence the value of n (given by the speed of scalar field (2.7)) does not drastically change in time.This feature would make it hard to provide the peaky spectral shape in PBH mass spectrum and associated induced GWs testable with NanoGrav event.Therefore, we would need to extend this model to realize it and we leave further consideration of this to future study.

B Calculation of integral I
In this appendix, we briefly present the computation of integral I in the power spectrum of induced GWs.In terms of (5.11)- (5.14), the original expression of I is obtained as I(ν, u, x) = Since the GW spectrum observed today is sufficiently within the sub-horizon regime, we consider the late-time limit of I: kτ c → ∞.In this limit, K 1 (ν, u, kτ c ) and K 2 (ν, u, kτ c ) converge to constant and their asymptotic expressions (5.19) and (5.20) have been found [86].Then, we need to take an oscillation average with respect to kτ c to evaluate the amplitude of spectrum.The product of two I's can be decomposed as where we use the momentum conservation, k = k ′ , and sin(kτ c ) 2 = cos(kτ c ) 2 = 1/2 and sin(kτ c ) cos(kτ c ) = 0, where overlines represent the average over τ c .

C Calculation of the Planar diagram of induced GWs
In this appendix, we present the computation of Planar diagram of induced GWs (see Figure 7).By using the fitting function of E k (3.10) and the following relationship where r pq ≡ p + q and F (defined in (3.12)) contains the integration over τ .
Here, we investigate the efficient method to calculate the planar diagram.Although all diagrams contain three loops and require 3 × 3 = 9 integration of momenta, the Reducible diagram is easy to estimate because we can firstly perform the momentum loop of two-form field.Similarly, integration of triangle loops ζ-ζ-B in Planar diagram can be performed at first to derive the "reduced vertex" of hBB:  where P ζ,v (k p ) is a power spectrum of curvature fluctuations directly induced by vacuum mode fluctuations on the peak scale k p , which is slightly larger than that in CMB scale for the Starobinsky inflation model (2.12).The time integration F is 1 ∼ 1.1 in our parameter setup.

Figure 2 .
Figure 2. The time evolution of E k which crosses the horizon at 25 e-folds before the inflation end.We normalize it to be dimensionless and multiply √ n ∝ I φ /I to match the source term of the inflaton perturbation in (3.3).The solid blue line shows the exact numerical solution of (2.6).The dashed black line shows the fitting function (3.10) with σ 2 = 1.0 and |τ peak | ∼ 2 × 10 −18 Mpc.The vertical dotted lines denote the time when n(t) = 1, n max and 0.1 from right to left.

6  E /H 4 Figure 3 .
Figure3.The power spectrum of the two-form field P E (k) at several e-folds; N = N (t 1 ) ≃ 4 (solid red), N = N (t 1 ) − 0.5 (solid orange) and N = N (t 1 ) − 1 (solid blue), where t 1 denotes the time when the growth of the two-form field stops, n(t 1 ) = 1, and the spectral peak gets a maximum value.The black lines show our fitting function (3.10) evaluated at N = N (t 1 ) (dotted) and N = N (t 1 ) − 0.5 (dashed), and N = N (t 1 ) − 1 (dot-dashed).They agree well until the peak amplitude becomes around two orders of magnitude less than that at t = t 1 .

Figure 4 .
Figure 4. Power spectrum of the curvature perturbation P ζ = P ζ,v + P ζ,s against the wave number k[Mpc −1 ].The solid blue line represents the full result which numerically estimates (3.11), whereas the dashed blue line shows the contribution only from the vacuum mode.The black dashed line denotes the approximate fitting function (3.14) with the parameter set (3.15).

Figure 6 .
Figure 6.The mass spectrum of PBHs in our model when we choose the parameter set in(2.15).The green region around M < O(10 17 )g is excluded by the extragalactic gamma ray[8] and the orange region O(10 22 )g < M is excluded by the Subaru/HSC[12].We note that the constraints plotted are the constraints on a monochromatic PBH mass function (whilst the black line the extended mass function predicted by the model presented here) -and so are not completely applicable to our mass function.However, the constraints for extended mass functions are similar to the constraints for monochromatic mass functions.

Figure 7 .
Figure 7. Diagrams of the contractions of 2-form field to calculate the power spectrum of induced GWs.We label these diagrams as "Reducible" (left), "Planar" (center) and "Non-Planar" (right).The external solid line represents the gravitational wave perturbation h +/× .The intermediate dashed (wiggly) line represents the curvature perturbation ζ (the two-form field B).
this work, we investigate an efficient method to calculate the Planar type diagram and show the reduced formulae on the Planar-type GWs in Appendix.C. The basic idea is as follows: the Planar diagram has two ζ-ζ-B loops and we can firstly contract them as two effective vertices of h-B-B.Afterwards, we can calculate only one residual B-B loop in terms of the obtained vertices.To reduce the computational costs, we approximate the profile of two-form mode function as a Gaussian fitting function (3.10).We numerically estimate Ω GW (k) for Reducible diagram (5.30) (dotted black line) and Planar diagram (C.4) (dashed black line) shown in Fig.