Chemical-Potential-Assisted Particle Production in FRW Spacetimes

We analyze gravitational particle production assisted by chemical potential. By utilizing the uniformly smoothed Stokes-line method and Borel summation, we gain insight into the fine-grained history of enhanced particle production. Analytic/semi-analytic formulae describing the production amount, time and width are obtained for both spin-1 and spin-1/2 particles in various FRW spacetimes. Our work also serves as a concrete demonstration of the uniformly smoothed Stokes-line method applied to cosmology.


Introduction
The study of particle production in cosmology is, in a sense, the study of the origin of matter itself. Arguably, the most natural choice of the initial state of the universe is the one in which strict spatial homogeneity is maintained, if space is to retain its physical meaning at all. Such a perfectly homogeneous state must also be devoid of any ordinary matter, which we view as a collection of fundamental particles. In addition, on the classical level, this strict homogeneity will be kept throughout the cosmic evolution if the Lagrangian respects translational symmetry. This would naively suggest a dull universe with no matter but zero modes only. However, this is not the whole story. The dynamics of the background spacetime can bring the vacuum fluctuations to reality, breaking the spatial homogeneity on a quantum level. This spontaneous breaking of spatial translational invariance is accompanied with the appearance of perturbative modes in spacetime, which upon quantization, become real particles.

JHEP06(2021)129
Such is the scenario for the paradigm of inflation [1][2][3][4]. The initial condition of the inflationary universe is chosen as a Bunch-Davis (BD) vacuum, i.e., a coherent state of the inflaton field that is annihilated by the positive-frequency part of all fields with non-zero momenta. However, the expansion of the spacetime stretches the wavelength of different modes φ k and produces a squeezed state with occupation number |β(k)| 2 = 0. For φ being the inflaton or the graviton, this occupation number is exponentially large so that these particles decohere and become essentially classical waves, which source the primordial inhomogeneity for the later cosmic evolutions. For φ being heavier degrees of freedom, the occupation number is finite and typically suppressed by an exponentially small factor |β(k)| 2 ∝ e −2πm/H , where H 10 13 GeV is the Hubble parameter during inflation. These gravitationally produced particles are the very first matter emergent from the BD vacuum in the inflationary universe and their interactions leave characteristic non-Gaussian imprints on the Cosmic Microwave Background (CMB) as well as the Large Scale Structure (LSS). This recently thriving field known as cosmological collider physics [5][6][7][8][9][10] thus has an intimate relation with the phenomenon of particle production from the vacuum. For example, the signal strength is directly proportional to the square root of the production amount, S ∝ |β(k)|. Thus heavy particles with m H are extremely difficult to probe if they are produced in the purely gravitational way.
This problem motivates several proposals where the exponential suppression can be alleviated. One possibility is to consider special inflation models. For instance, in axionmonodromy inflation [11,12], time-dependent mass terms violate adiabaticity and lead to a dramatic amplification of particle number density and thus the size of non-Gaussianity [13]. In warm inflation [14,15], thermally produced heavy particles are Boltzmann-suppressed by an alternative temperature much higher than the Hubble scale, also giving rise to enhanced cosmological collider signals [16].
Alternatively, one can consider the interesting possibility of introducing a chemical potential. This can be naturally achieved via a rolling scalar field φ coupled to massive fields through operators of the form ∂ µ φJ µ , where J µ is a certain current made from massive fields. For massive spin-1 particles, choosing the Chern-Simons current J CS results in an enhanced production of vector boson [17][18][19][20][21]. For massive spin-1/2 fermions, choosing the chiral current J 5 leads to a natural amplification of fermion production rate [22][23][24][25]. Similar amplification is also found for charged scalar particles if we generalize the chemical potential operator to κ µ J µ , and reinterpret the chemical potential κ µ as a background gauge field creating Schwinger pairs [26][27][28]. These mechanisms of chemical-potentialassisted particle production typically break parity or rotational invariance, hence leaving sizable and distinctive signatures in cosmological observables.
In addition to physics during inflation, chemical potential can also play a role in the late universe. The aforementioned chemical potential introduced by a rolling Axion-Like Particle (ALP) generates a tachyonic instability in the gauge boson sector. This can, for example, efficiently convert the ALP to Dark Photon Dark Matter (DPDM) [29], and produce chiral gravitational waves [30][31][32]. In the fermion sector, chemical potential also sources the helicity asymmetry of fermion numbers, which can be important for baryogenesis [22,33].

JHEP06(2021)129
Thus, chemical-potential-assisted particle production is a generic phenomenon that appears in many different contexts and setups. A general and systematic investigation of particle production in the presence of chemical potential is necessary.
In this work, we re-derive Berry's uniformly smoothed Stokes-line method [34][35][36][37] and apply it to analyze the fine-grained production history of particles of mass m with chemical potential κ. By fine-grained production history, we mean the production amount, time and width (duration) exact to the leading order in the super-adiabatic expansion. We found that for a constant chemical potential κ = const, spin-1/2 fermions and spin-1 vector bosons share a similar production history, with a simple yet subtle replacement rule m 2 ↔ m 2 +κ 2 . We also derive analytic/semi-analytic formulae for the production history in five common FRW spacetimes.
We note that there are many past studies in the literature that utilize the Stokes phenomenon to study particle production. For example, it is applied to the Sauter-Schwinger effect [38][39][40][41][42][43], to Hawking radiation [44], to the adiabatic particle number in global de Sitter (dS) spacetime [45][46][47], to dark matter production at the end of inflation [48], to preheating [49,50], and to particle production triggered by vacuum decay [51]. An excellent review is recently given in [52]. However, we point out that most of them (with the exception of [40,47,51]) focus on the asymptotic production amount far away from the Stokes line. To our knowledge, the analytic calculation of production time and width presented in this work is a new ingredient in this area, with or without chemical potential. Understanding these fine production details is useful, for example, in the loop-level estimation of cosmological collider signals, or in estimation of backreaction time scales. We hope this work also serves to demonstrate the application of the uniformly smoothed Stokes-line method applied to cosmology. This paper is organized as follows. In section 2, we first discuss the generalities of chemical potential and define the model we work with. Then in section 3, we re-derive the uniformly smoothed Stokes phenomenon for both bosonic and fermionic systems and justify its validity for the case with significant particle production. In section 4, we move on to work out the production details for spin-1,1/2 particles in five common FRW spacetimes. At last, we summarize and give outlooks in section 5. For readers who wish to skip the detailed analysis to directly look up the results, we assemble our formulae into a checklist in appendix A.

Chemical potential: generalities and dS solutions
In this section, loosely following [53], we give some general discussions on chemical potential in cosmology. In thermodynamics, chemical potential is originally introduced by Gibbs to describe the change of internal energy of a system with respect to the change of particle numbers when the entropy and volume are held fixed. More formally in statistical mechanics, it is identified with the Lagrange multiplier κ of total particle number when the systems in a grand canonical ensemble are allowed to exchange particles with each other. Starting from the partition function for a grand canonical ensemble, Z = e −(H−κN )/T , we can straightforwardly generalize it to the field theory context as a path integral in the JHEP06(2021)129 phase space of a field φ, (2.1) where H is the Hamiltonian and N is the particle number operator associated with a certain symmetry (which may or may not be exact). Throughout this work, we will assume a spatially flat FRW spacetime background with ds 2 = −dt 2 + a(t) 2 dx 2 . Going into the field configuration space and write the chemical potential term as the integral over a local density, we have Z = Dφe i dtd 3 x √ |g|(L(φ,∂φ)+κJ 0 (φ,∂φ)) , (2.2) where N ≡ d 3 xa 3 J 0 . Now we can turn on the spacetime dependence of κ and interpret it as the zeroth component of a local vector field κ µ (t, x). Thus the general form of a chemical potential as a background field coupled to the matter field is Now if one inspects the effect of introducing such a chemical potential term into the matter Lagrangian, the result will depend on two aspects. First, in the absence of κ µ , if the matter current is derived from an exact symmetry, it is conserved as an operator identity, i.e., ∇ µ J µ = 0, where ∇ is a covariant derivative with respect to the metric g. Then one can always consistently gauge this symmetry by minimally coupling the conserved current to a vector potential. We are free to choose κ µ (x) to be such a background gauge field. If the chemical potential κ = κ µ dx µ is closed in the 1-form sense, d D κ = dκ + κ ∧ κ = 0, where D is compatible with the gauge field connection κ, the background gauge field then has no field strength and is gauge-equivalent to vacuum. Thus we can perform a gauge transformation to eliminate κ locally. Such is the case if κ = κ(t) is Abelian and spatially homogeneous. However, we point out that the elimination of the chemical potential term is not completely trivial in the scalar case, as we will see below, since it shifts the scalar mass in a quadratic way. The second possibility is that if the matter current is not built from an exact symmetry and is hence not conserved. This suggests that there is no consistent way of coupling it to a background gauge field and thus interpreting it as the chemical potential. In summary, we give the following necessary condition of a chemical potential term that has physical effects other than quadratically shifting the mass of the matter particle, In the following discussion, we provide several examples that satisfy the above criterion and have interesting particle production features. We will start the general discussions in flat FRW spacetime and retreat to exact dS spacetime when solving the Equations of Motion (EoMs). The general FRW EoMs will be discussed in later sections.

Spin-0
Consider a complex scalar field σ, the only non-trivial chemical potential term we can find at quadratic level is with J µ = i(σ * ∂ µ σ − σ∂ µ σ * ). If the original Lagrangian is U(1)symmetric, this current is conserved. Therefore, according to (2.4), we have to go for the JHEP06(2021)129 first possibility. 1 We can write the following action in the FRW background, where g µν ≡ a(τ ) 2 η µν using comoving coordinates. Absorbing this chemical potential term into the derivative term and shifting the mass term accordingly, we have As mentioned above, if κ µ =κ 0 (τ )δ 0 µ , the whole system can be considered as a charged scalar moving in the vacuum with a new mass There is no enhancement of particle production unless M 2 (τ ) becomes negative or its time dependence violates adiabatic condition.
If the first criterion is satisfied, then there is a non-zero field strength where prime denotes a derivative with respect to the conformal time τ . The enhancement of particle production is exactly the Schwinger effect of this electric field [28]. The EoM of σ in momentum space reads where κ ≡ κ i aê i . Clearly the second term in the square bracket breaks rotational symmetry and stands for the effect of the background electric field. It introduces an angular dependence in the effective mass of different modes. Those with lighter effective mass tend to get produced more easily, especially for the mode traveling at the same direction as κ. We shall see that this term linear in momentum is typical in the presence of chemical potentials. They represent the bias on the effective mass of different modes introduced by the chemical potential. For large enough |κ|, there is even a tachyonic instability and the pair creation rate becomes exponentially large. To be more quantitative, we set |κ| = const and limit ourselves to dS by taking a = − 1 Hτ . Then the solution to the EoM is given by wherek = k |k| ,κ ≡ κ H ,m ≡ m H and W is the Whittaker function that matches the BD initial condition at τ → −∞. The late-time expansion reveals an angular-dependent production amount In fact, even if the U(1) symmetry is explicitly broken with ∇ · J = 0, it can be shown that there is only a quadratic mass-shift and no enhancement particle production is present [53]. This is why (2.4) is only a necessary condition, rather than being sufficient. However, we note that a different opinion is given in [54], where it is argued that a scalar chemical potential can also bring isotropic enhancement to particle production.

JHEP06(2021)129
To gain an intuitive understanding of the enhancement, we can go to the large mass limit m H. Then the leading order particle number is Therefore, the direct consequence of a chemical potential is a linear and biased (rather than quadratic and un-biased) shift of the effective mass of the particle modes, making them easier or harder to produce. 2

Spin-1 2
In the massive spin-1 2 case, the second possibility of (2.4) can be satisfied by choosing the axial current J µ 5 = e µ aΨ γ a γ 5 Ψ. Hence a time-like chemical potential κ µ (τ ) ∝ δ 0 µ has no interpretation as a trivial background pure-gauge and can become physically relevant. To illustrate how it assists gravitational particle production, we and choose a Majorana fermion model [23] written in a Weyl basis Ψ = ψ ψ † . The Dirac fermion case can be obtained by combing two Majorana fermions with analogous behaviors. The action of a Majorana fermion with chemical potential reads This chemical potential term can arise from, for instance, a dimension-5 coupling to a rolling scalar, After choosing a tetrad e a µ = a(τ )δ a µ , e µ a = a(τ ) −1 δ µ a and κ µ = a(τ )κδ 0 µ with κ =φ Λ = const, the action simplifies to where for simplicity, we have rescaled ψ → a −3/2 ψ. In momentum space, we can perform a standard decomposition into helicity eigenmodes, where σ ·kh s (k) = sh s (k). The EoM reads

JHEP06(2021)129
This set of equations is exactly solvable in dS. With a BD initial condition, the mode functions take the form wherem ≡ m H andκ = κ H . Similar to the scalar case, performing an IR expansion at τ = 0 and matching the Bogoliubov coefficients, one can arrive at the production amount formula [22] When the mass is large and chemical potential is small, the leading order particle number again takes the form of a Boltzmann factor with linearly biased effective mass, For a positive chemical potential, the negative-helicity mode gets amplified whereas the positive-helicity mode is suppressed. However, when the chemical potential is larger than the mass scale, the enhancement in the negative-helicity particle production begins to saturate, This is essentially the Pauli blocking phenomenon generic to all fermionic systems. The exclusion principle forbids any mode being occupied more than once. Actually, fermion production with constant chemical potential in general FRW spacetimes can be understood in an elegant way. In terms of the physical time t, the EoM is essentially a two-state system evolving according to a Schrödinger equation (2.20) Without loss of generality, let us consider the negative helicity mode with s = − (the s = + helicity mode is obtained by k → −k). The unitary evolution governed by a Schrödinger equation (2.20) preserves the normalization |u − | 2 + |v − | 2 = 1. The BD initial condition selects the positive frequency mode u − in the early time limit t → −∞. If κ > 0, there can be a time when the physical wavelength k a of the mode is comparable to the chemical potential scale κ. In the language of quantum physics, we have an avoided crossing. Namely when the diagonal elements of the Hamiltonian vanish, its instantaneous eigenvalues approach each other, but a complete degeneracy is avoided due to the offdiagonal terms. If this process occurs adiabatically, according to the adiabatic theorem, the state smoothly maintains its positive-frequency trajectory and there is not much particle production. However, if κ is large, the avoided crossing becomes a non-adiabatic one, with JHEP06(2021)129 almost all positive-frequency crossed into negative-frequency part and thus nearly maximal particle production. Viewed in this way, Pauli exclusion principle is a built-in feature of fermions such that the evolution of (u, v) is unitary, as opposed to the symplectic evolution of bosons (σ,σ), which can be unbounded from above and exponentiating (e.g., if κ m in (2.10)).
Interestingly, the large-chemical potential limit in the spin-1/2 fermion case enjoys a universal behavior in general FRW spacetimes. To see this more explicitly, we can assume κ m and expand around the crossing time t * , where κ = k/a(t * ), This is none other than a Landau-Zener (LZ) model with η = m and γ = 2κH(t * ) [55,56]. The corresponding LZ parameter describing adiabaticity is z = η 2 γ = m 2 2κH(t * ) , where z 1 corresponds to adiabatic transitions and 0 < z 1 corresponds to diabatic transitions. See figure 1 for justification of approximation using an LZ transition in dS. Thus the crossing probability is given by the exponential factor κH(t * (k)) , (2.22) in agreement with the exact dS result (2.19). Furthermore, we obtain the production time t * as the solution to the equation κ = k a(t * ) . Naively, we would expect the production width to be the time scale at which the LZ resonance happen, i.e., 4η γ = 2m κH(t * ) . However, we shall see in section 3.2 that this is not the case. The actual particle production width is , which is typically shorter than the LZ time scale [36].

Spin-1
One can also easily find a non-conserved current in the Abelian massive vector boson case, namely the Chern-Simons current J µ The chemical potential term looks like For a choice of κ µ = θ (τ )δ 0 µ ≡ a(τ )κδ 0 µ , or κ ≡θ = const, an integration-by-part gives a time-dependent θ-term (a rolling axion), Thus the system is described by an axion electrodynamics Lagrangian [57] A massive vector boson has three degrees of freedom, two transverse modes and a longitudinal mode, from which the time-like component is solved using the constraint ∇ µ A µ = 0. Here we have chosen m = 15H, κ = 60H, or, after translating to LZ model parameters, η = 15H, γ = 120H 2 . The gray band corresponds to the naive time scale 4η γ = 0.5H −1 while the red band corresponds to the actual production width ∆t * = 2π γ = 0.23H −1 . As long as κ m, the light pink band is narrow enough so that the expansion of spacetime becomes irrelevant, and the instantaneous energy eigenvalues approach to that of a LZ model.

JHEP06(2021)129
Decomposing the spatial components into helicity eigenstates, we have The EoM then reads As a result, the transverse modes are affected by the chemical potential while the longitudinal mode is not. Focusing on the transverse modes, we can solve the EoMs analogously as (2.7), (2.28) The transverse particle production amount is then (2.29) In the large mass limit, the effect of chemical potential again simplifies to a linear bias over the effective mass,

JHEP06(2021)129
Here, without the protection of the exclusion principle, the production amount starts to become exponentially large if |κ| > µ. This is commonly recognized as a tachyonic instability in the study of axions. In this parameter regime, the backreaction to the rolling θ(τ ) must be taken into account.

The uniformly smoothed Stokes-line method
In this section, we derive the uniformly smoothed Stokes-line method for both spin-1 vector bosons and spin-1/2 Majorana fermions, providing a framework to analyze the histories of particle production with chemical potentials. Mathematically, the evolution of spin-1 vector bosons (2.27) and spin-1/2 Majorana fermions (2.20) behave as second order differential equations and the Schrödinger equation with two quantum states respectively, and such systems can experience the emergence of negative-frequency part starting from an initial positive-frequency solution, known as the Stokes phenomenon: where Ψ ± (τ ) are the positive/negative-frequency parts of either bosons or fermions, and the two time-dependent functions α(τ ) and β(τ ) can be regarded as the Bogoliubov coefficients, associated with the particle production. The Stokes phenomena in second order differential equations [35,40,47,48,58,59] and transitions between two quantum states [36,37] have been studied in many works with the assumption that the magnitude of the emergent part |β(τ )| is exponentially small, and all of these works point out that the singulant defined as the difference between the positive and negative phases accumulated from the complex turning point τ c satisfying ω(τ c ) = 0, are important to describe the details of the Stokes phenomena. To be specific, when the systems evolve near the Stokes line, defined as the line linking τ c and τ * c with ImF (τ ) = 0, the negative-frequency part starts to produce, and the production histories including the amounts and widths can be calculated from the singulant F (τ ), proved with the technique of optimally truncating the asymptotic series solution of (3.1).
In the following subsections, we first follow previous studies to apply the optimally truncated asymptotic series solution to calculate the particle production for vector bosons and fermions respectively, and this framework works properly when the particle production is exponentially small. We then analyze the situation when the exponential particle production is significantly enhanced by the chemical potential, so that the optimal truncation technique is not applicable, and the Borel summation technique should be applied to obtain the particle production.

Bosonic case
Consider the mode expansion for the transverse component of a massive vector boson,

JHEP06(2021)129
The mode function satisfies a second order EoM of the general form where we use the dimensionless variable 3 z = kτ , and we denote as the derivative with respect to z starting from here. The helicity label s is also omitted for simplicity. The asymptotic parameter λ, defined for later analysis of asymptotic series, is supposed to be large, and we can choose λ = m if the mass is the largest parameter in the problem. (3.4) has the same structure as the one-dimensional time-independent Schrödinger equation with a barrier, and therefore we adopt the method of analyzing waves near Stokes lines [35] to study the particle production. We can express the solution of (3.4) with the WKB form where C is a constant fixed by the initial and normalization conditions, z c is the complex turning point defined by w(z c ) = 0 and located at the lower-half complex plane, and the function W (z) satisfies In general, we cannot obtain the exact solution of W (z), but we can apply the iterative adiabatic expansion to approximate it [40,47] with W (0) (z) = w(z). Such an iterative relation can be used to derive the asymptotic series solution of W (z) boson's dS effective mass. On the other hand, it is well-known that the terms in the asymptotic series (3.8) keep increasing when n is sufficiently large [35,40,47,58] is Dingle's singulant variable [59]. So we can truncate the series sum at a suitable order n to approximate the solution (3.5).
The asymptotic series solution with an optimal truncation order n, , cannot fully represent the WKB solution (3.5), but we can choose to expand the exact solution with the super-adiabatic basis formed by f (n) and f * (n) , where z i in the value of z at initial time and g(z) can be viewed as the instantaneous positive-frequency solution. Now the vector field can be expanded in an alternative form using g(z): where the new annihilation operator acquires a time dependence through the Bogoliubov And the original vacuum annihilated by a s k now contains a spectrum of particles, (3.14) Thus our aim is to solve the time dependence of the Bogoliubov coefficients α(z), β(z). The solution satisfies the initial and normalization (Wronskian) conditions as

JHEP06(2021)129
As pointed out in [40,47], the constant Wronskian implies a degree of freedom of defining the derivative of f : , (3.16) where V (z) is an arbitrary real function. Choosing the time-dependent function V (z) decides the evolution of the Bogoliubov coefficients d dz and with the initial condition Therefore, the appropriate choice of V (z) should minimize the change of the Bogoliubov coefficients as we intends to minimize the difference between the basis function f (n) and the exact solution, and such a choice is V (z) = − W (n) (z) 2W (n) (z) , as suggested in [35,40] with different reasons. With this choice, the evolution of the Bogoliubov coefficients satisfies d dz (3.20) and the diagonal term can be removed by defining variables implying that S ± are the Stokes multipliers for the positive and negative modes respectively, so the evolution equation is simplified as (3.22)

JHEP06(2021)129
The term in the square bracket can be interpreted as the O(λ −2n−1 ) error when we approximate the EoM (3.4) with the 2n-th order partial sum of the asymptotic series, and this can be calculated with the asymptotic series of the positive-frequency part f + Substituting this series solution into the EoM In all the scenarios that we study in section 4, w 2 (z) has a simple root at z c , so we approximate w(z) as 26) and the solution of this recurrence relation Applying the asymptotic series of f (z) (3.23), The term in the square bracket of (3.22) can be calculated explicitly

JHEP06(2021)129
where (3.24) is applied to obtain the last line. Therefore, the evolution of S ± (3.22) is reduced to where we keep only the term with O(λ −2n−1 ). It is clear that |dS + /dF | |dS − /dF | as there is an exponential suppression e −F for the former, so we can solve (3.30) perturbatively starting from the initial values S i ± = S where the value of prefactor R n = Γ 2n + 7 6 Γ 2n + 11 6 (2n + 1)!Γ(2n + 1) , (3.32) and the functionΓ(−1 − 2n, −F ) is the continuous version of the incomplete Gamma function, defined as The incomplete Gamma function of (3.31) can oscillate dramatically for general n, and we can choose an optimal truncation order n such that the phase is stationary at F = ReF , the moment when the integral receives dominant contribution and the solution is n = Int ReF It is noteworthy that the optimal truncation order cannot be determined by the stationary-phase condition when ReF < 2, and we will consider such a situation in the later part of this subsection. For simplifying the following perturbative calculation, we set the initial condition as − ) = (1, 0), and the cases with general initial conditions can be obtained easily based on (3.31).

JHEP06(2021)129
Assuming that the optimal truncation is applicable with ReF ≥ 2, we can thus approximate the integrand as a Gaussian function around the point with stationary phase, and thus the Stokes multiplier reduces to an error function.
We are ready to calculate the first-order term of S + by solving where z * is the intersection between the Stokes line and the real z axis, and the phase integral along the Stokes line e  37) and the definition of S ± (3.21) implies that the normalization of the Bogoliubov coefficients preserves under the first-order perturbation. For the situations with ReF < 2, we expect that higher-order perturbations are required to solve for (3.30), and the perturbation theory may break down when e −ReF → 1, so we should analyze such situations carefully. In the cases with ReF < 2, we cannot choose an optimal truncation because of the failure of the stationary phase condition (3.34) and the magnitudes of the terms in the asymptotic series solution (3.23) b n increase, starting from the first term. Such a divergent series defined by (3.27) behaves like the generalized hypergeometric function 2 F 0 (a, b; ; F −1 ) which diverges everywhere from its original definition, but it can be defined meaningfully by applying the Borel summation: where K a (z) is the modified Bessel function. The function B(F ) is discontinuous when it crosses the Stokes line with ImF = 0, implying that the exact solution after crossing the Stokes line should depend on different set of linear combination: where the C 1 and C 2 are constants. After knowing the expression of f (z), the Bogoliubov coefficients can be solved by combining (3.11) and (3.16), and the results can be fully recorded by the singulant F and and thus the constants C 1 and C 2 are chosen such that the Bogoliubov coefficients and their derivatives are continuous at ImF = 0. To compare with the particle production in dS spacetime, we set the initial condition as For example, the |β(k)| in dS spacetime (2.29) cannot be fully described as an exponential factor in some parameter ranges, and thus the behavior of large |β(k)| depends on the details of scenarios. The universal property is that it approaches to the exponential form e −ReF when |β(k)| decreases, so we use the exponential form to describe the tendency of the production amount but not its exact value for the cases with small ReF .
On the other hand, we also compare the Stokes multiplier S num (F ) obtain from the numerical result (3.41) with the approximations utilizing the incomplete gamma function S Γ (F ) with n = 0 from the first-order perturbation (3.31) and the error function S Erf (F ) (3.35) respectively, as shown in figure 3. It is clear that only the imaginary part of β(k) remains non-zero after finishing the particle production, and different approximations have significant errors of describing the real part of β(k) which vanishes rapidly after crossing the Stokes line, but the error function can still describe the width of the production process. After knowing how the production amount and width depends on the singulant F (z), we can write down a simple form of β(z) which includes the tendencies of particle JHEP06(2021)129  production with vacuum initial condition when related parameters are changed: where we replace λw → w in (3.4) for simplicity since λ and w(z) always appear together in the final results, and the generalization to arbitrary initial conditions is straightforward based on the results (3.31) and (3.41).

Fermionic case
The fermionic case is logically similar to the bosonic case, but with important differences in the mathematical details. To be more specific, let us consider the action (2.13) for a Majorana fermion with chemical potential in an FRW background. We will use the Van der Waerden notation for two-component spinors and the conventions follow from [60].

JHEP06(2021)129
We begin by rewriting the mode expansion of ψ α in the Van der Waerden notation, Notice that it is sometimes customary to omit the zeroth Pauli matrix since it is an identity matrix in component form. However, for the sake of balancing the indices, we will keep them explicit here. The eigenvalue equation for the helicity basis can be written in a number of different equivalent forms: It is sometimes useful to choose an explicit component form of the helicity basis, withk pointing toward the (θ, φ) direction in spherical coordinates. Substituting the mode expansion into the equation of ψ α obtained from varying the action (2.13), we obtain the EoM of the mode functions, This EoM preserves the combination |u s | 2 + |v s | 2 , with the normalization constant determined by the canonical quantization condition {ψ α (τ, x), . After plugging in the mode expansion, this is reduced to a c-number equation Taking the trace and the determinant of the above equation yields thus fixing the normalization condition |u s | 2 + |v s | 2 = 1 separately for different helicities. As mentioned in section 2.2, the EoM (3.46) can be interpreted as describing the transition of a two-level system with a Hamiltonian The instantaneous eigenstates of H(τ ) are given by

JHEP06(2021)129
Therefore, an ansatz of the solution of (3.46) can be constructed as whereC s ,S s are slowly varying functions whose detailed form as a super-adiabatic basis will be computed later. If we choose |C s | 2 + |S s | 2 = 1, the coefficient functions α s , β s will satisfy the normalization |α s | 2 + |β s | 2 = 1, as required by the normalization of u s , v s and unitarity. Now we can insert the ansatz back into the mode expansion of ψ α , The time-dependent creation/annihilation operators are selected according to the instantaneous negative/positive frequency parts of ψ α . Therefore, we can regroup the terms according to the dynamical phase e ∓i Esdτ . First, we note the relation where η s (k) is a phase factor satisfying (3.56) where the new time-dependent annihilation operator is obtained as a Bogoliubov transformation [22,61], The anti-commutation relation is preserved: Therefore, the vacuum annihilated by the original operator b s k now contains a spectrum of particles with comoving number density

JHEP06(2021)129
where V is the comoving volume. Here rotational symmetry demands the isotropy of particle production spectrum, as β s only depends on the magnitude of the momentum. The structure of the EoM of fermion (3.46) is similar to the two-state systems in quantum mechanics, and thus we adopt the framework of analyzing the quantum transition histories of such systems [36]. To calculate the particle production, it is convenient to use the bra-ket notation. Similar to the case of boson, we rewrite the equation of motion with the asymptotic parameter λ where we use ψ to denote the two-component mode function, and z = kτ . For the ansatz (3.52) in the bra-ket notation the left hand side of (3.60) is whereas the right hand side is If |ψ α and |ψ β are the two exact solutions, the positive and negative modes evolve independently with constant α and β, implying that Similar to the case of boson, we approximate the solutions with the asymptotic series where |ψ Similar to the case of boson, we solve for c j (z) and d j (z) near the complex root z c of In all the scenarios that we study in section 4, z c is a first-order root, implying that . (3.69) With this approximation of θ (z), we obtain the recurrence relation of c j from (3.67) where A is defined by E ≈ A(z − z c ) 1/2 , and the solutions of c j (z) and d j (z) with c 0 = 1 and d 0 = 0 are and respectively, where (a) k = Γ(a + k)/Γ(a) is the Pochhammer symbol, and the singulant is Clearly c j (z) and d j (z) are divergent asymptotic series, and thus we truncate them at the order of n and denote the partial sums of (3.65) as |ψ (3.74) The values of matrix elements can be obtained by evaluating (3.64) with the truncated series and thus d dz where which has the same structure as vector boson (3.20). We remove the diagonal term by defining the Stokes multipliers and d dF By setting the vacuum initial condition (S i − ) = (1, 0), the first-order perturbation of S + implies that can be approximated as Similar to the bosonic case, the situations with ReF < 1 implies the failure of choosing an optimal truncation for the asymptotic series (3.65), and we may apply the Borel sum to evaluate such divergent series: where U (a, b, z) is the confluent hypergeometric function, and the last line is obtained from the relation between c j and d j (3.67). Since I(F ) and J(F ) are discontinuous at ImF = 0, we rewrite the exact solution (3.61) into two parts:  where α i and β i are the initial Bogoliubov coefficients, C 1 and C 2 are constants to let |ψ and its derivative continuous at ImF = 0. The time-dependent Bogoliubov coefficients defined with respect to the eigenstates |ψ  Similar to the bosonic case, the asymptotic parameter λ always appear with the instantaneous eigenvalue E(z), so we may set λH → H in (3.60). We summarize a simple form of β(z) which reflects the tendency of the particle production starting from vacuum initial condition where it is noteworthy that there is no additional prefactor −i compared to the bosonic case (3.42), and the generalization to arbitrary initial conditions can be easily done base on (3.81) and (3.86). Note that the production histories solely depend on the singulant F , and the definition of the singulant F in the fermion case is similar to the boson case, with

Analysis of particle production in various spacetimes
Armed with these powerful mathematical tools, we are now in a position to compute the fine-grained particle production histories of both massive vector bosons and Majorana fermions in various setups.
In this section, we have in mind that the chemical potential is provided as a external source by, for instance, a rolling scalar field. The backreaction to the external field that generates the chemical potential is also assumed to be negligible. In particular, we will assume the chemical potential κ is a constant in spacetime. The reason for such a choice is three-fold. First of all, this is indeed true in some cases. For example, the Hubble friction during inflation drives a rolling scalar to an attractor phase with constant speedφ, which corresponds to a constant κ when coupled to vectors or fermions. In a radiation/matterdominated universe, specifically chosen scalar potentials also give rise to constant rolling JHEP06(2021)129 speeds. Second, physically speaking, for any slowly-varying κ(τ ) = 0, a constant chemical potential is always a leading order approximation. As long as the typical time scale of κ(τ ) is longer than the particle production time scale, this approximation will be valid. Third, mathematically speaking, a constant κ leads to simple and analytical results that already contain lots of information in the general cases, which can always be dealt with using numerical methods.
We will focus on five familiar types of FRW spacetimes whose singulant integrals are exactly computable. The resulting production amount, time and width are given explicitly as analytical expressions. Some of these results are exact while others are approximate or empirical with percent-level error in most parameter regimes. To distinguish them from each other, we will use = when the result is exact. We use for results which are easily computable to any desired precision but which are shown with finite accuracy. And ≈ will be used for empirical results whose relative error is at percent-level.
Due to the replacement rule mentioned in section 3, we will only work out the spin-1 case with |κ| < m, and obtain the spin-1/2 results for all parameter regions by simple substitutions. Also because different helicities are related by a sign flip of κ, we will focus on the negative helicity state, whose production is enhanced if κ is positive. Throughout this section, we will be working in comoving coordinates and using conformal time rather than cosmic time. 4 Tilde variables will be used to define dimensionless parameters measured in units of a certain Hubble scale, e.g.,m ≡ m H ,κ ≡ κ H , etc. And we will typically expand quantities in powers ofκ m = κ m .

dS
In an exact dS spacetime, the scale factor has a time dependence a(τ ) = − 1 Hτ = e Ht . The EoM for spin-1 particles with chemical potential can be written in terms of a dimensionless variable z = −kτ = k aH , Notice that the variable z now runs from the right to the left, and the physical region is the positive real axis z > 0. After analytical continuation, one can define w(z) on the whole complex plane. The two roots of w 2 (z) = 0 lie at z c and z * c , with Notice that z c lies on the upper half complex plane since the original τ c is on the lower half complex plane and they differ by a sign. The absence of tachyonic instability requiresm > |κ|. Therefore, the two complex turning points lie in a symmetric fashion across the real axis. Starting from z c and z * c are two branch cuts that meet their ends at the pole at z = 0. In the z-domain, the singulant is evaluated as Notice the sign change compared to (3.10). The phase integral is exactly solvable: where the first line is purely imaginary for z lying on the positive real axis. Its behavior on the right-half z-plane as well as the Stokes lines are shown in figure 6.
• Production amount. The production amount is straightforwardly given by the formula with z i m taken on the real axis.

JHEP06(2021)129
• Production time. The crossing time z * lies on the real axis, and therefore satisfies a real-numbered equation: . into the equation (4.6) and collect the terms order-by-order inκ m , we are able to solve the coefficients iteratively, In this way, we obtain the particle production time as z * (m,κ) 0.6627m + 0.3435κ − 0.0102κ 2 m + 0.0064κ 3 m 2 + · · · . (4.9) Clearly, with a larger effective mass and a larger positive chemical potential, the production time becomes earlier.
• Production width. The derivative of the singulant is none other than the frequency itself: (4.10) Thus the production width in the z-domain is .

(4.11)
We can also translate the production width into the t-domain by ∆t * = ∆z * Hz * . Another useful measure of production width is the e-folding numbers during which the production is complete. In terms of a power series inκ/m, we have

JHEP06(2021)129
Now the alert readers may find an inconsistency here. If one compares the production amount (4.5) computed from the Stokes-line method with that of the exact result, namely (2.29) or (2.30), one finds that there is a mismatch of mass: the true result should contain the dS-corrected mass µ = m 2 − 1 4 instead ofm as given by the Stokes-line method. This mismatch implies that the results obtained by the naive application of Stokesline method are subjected to a relative error O(m −2 ), which can be important ifm is small.
To trace the origin of this 1/4 puzzle, let us go back to the super-adiabatic basis for the EoM (4.1), where S(z) is the Stokes multiplier. W (z) is solved order-by-order as specified in section 3.1.
The leading order reads W (0) (z) = w(z) and This leading order solution is sometimes called WKB approximation. The late-time behavior of the frequency function is Upon integration over z, the phase of the positive-frequency mode becomes linearly increasing with cosmic time t: where we have used z = −kτ = − k H e −Ht . Hence the late time behavior of the WKB basis is f (0) ∼ e ∓imt . This, however, corresponds to a wrong oscillation frequency for spin-1 vector particles. The correct frequency can be easily obtained by inspecting the late-time behavior of the EoM itself. Namely we can plug in the ansatz f ∼ z ∆−1 and expand (4.1) to leading order in z, This gives ∆ ± = 3 2 ± iµ and the correct IR behavior f ∼ z (1/2±iµ)t ∝ e ∓iµHt = e ∓imt . As a result, the leading order super-adiabatic solution does not capture the correct IR oscillation frequency, which is dictated by dS symmetries. Particles in dS are classified according to the unitary irreducible representation of the dS group [62,63] and a massive spin-S particle in the principal series has a conformal weight [10] ∆ (S) In other words, the geometry of dS modifies the effective mass of spinning particles in the IR (small z), and this fact is not taken into account by the naive WKB approximation (4.14).

JHEP06(2021)129
Fortunately, the advantage of the smoothed Stokes-line method is that the higher-order terms in the super-adiabatic basis can, and actually do, give essential corrections to the leading order solution. The first order correction to W is This function has two third-order poles at z c and z * c . The branch cuts brought by w(z) are still present and they connect the third-order poles to a simple pole at z = 0, The effect of this pole at the origin is exactly to give an O(m −2 ) correction to the IR oscillation frequency: Including the higher order corrections in the super-adiabatic series, we recover the correct IR oscillation frequency, Therefore, the 1/4 puzzle can be resolved by taking into account the full super-adiabatic basis and resumming the higher-order corrections to the mass. We can recover the correct production amount by redefining the singulant as an integral of W instead of W (0) = w. Technically, since this integral is ill-defined around the complex turning points where δW (n) (z) diverges, we need to manually impose a principal value prescription. We deform the integration contour to lie along the branch cut and tour along a semi-circle around the pole at origin (see figure 7 for illustration). Then by some arguments of complex analysis, the only non-zero contribution comes from the semi-circle at the origin, where the IR frequency correction is at work. Thus the corrected production amount can be written as half of the residue of W (z) at z = 0 + , subtracting half of the residue at z = ∞, For a detailed mathematical discussion of these arguments, we refer the readers to appendix B. In summary, the super-adiabatic corrections completely fix the mass mismatch and we only need to replacem → µ in (4.5), (4.9) and (4.12) to recover the true results. We list them below for the sake of clarity. • Production time (4.25) • Production width For spin-1/2 fermions, however, there is no such a problem, since the leading order WKB result already capture the correct IR behavior. The instantaneous eigenvalue of the Hamiltonian for its EoM (2.20) written in z-domain is Thus the mode functions behave as u, v ∼ e ±i E ± (z)dz ∼ e ∓i √m 2 +κ 2 t . This oscillation frequency indeed agrees with the late-time behavior of (2.20). Therefore, higher orders in the super-adiabatic basis do not offer any O(m −n ) corrections to the oscillation phase, hence to the production history. The production amount, time and width can simply be obtained by replacing µ → √m 2 +κ 2 in (4.24-4.26). Its Stokes lines are shown together with spin-1 particles in figure 8.
Finally, we plot the parameter dependence of production histories for both bosons and fermions in figure 9. From the plots, one can see that bosons and fermions share similar JHEP06(2021)129  . The production histories in exact dS. Left panel: production amount as a function of the dimensionless chemical potential for different particle masses. Middle panel: the z-domain production time dependence on chemical potential and mass. Right panel: production width measured in e-folding numbers. In all three plots, solid lines represent spin-1 bosons while dashed lines represent spin-1/2 fermions. Particles with different masses are distinguished by the colors of the lines according to the legend in the left panel. For bosons, we limit the range of chemical potential to be smaller than the mass, so that no tachyonic instability is induced. For fermions, the chemical potential is not restricted and we allow it to take arbitrarily large values. production histories when the chemical potential is small. However, they begin to depart from each other as |κ| becomes large, and in the end their large chemical-potential limits are drastically different: bosons enter the tachyonic regime, while fermions saturate and approach an asymptotic limit where z * ∝ |κ| and ∆N * ∝ |κ| −1/2 . Note that for bosons, both z * and ∆N * are monotonic functions of the chemical potential, whereas neither are monotonic for fermions.

Deviation from dS of the -type
The exact and rigid dS spacetime is maximally symmetric with a simple time dependence in the scale factor. The mode functions for free fields are exactly solvable, even though knowing the production time and width still requires the technique of smoothed Stokes phenomenon. However, dS only works as a leading order approximation to certain stages of cosmic evolution such as inflation or dark-energy dominated era. The actual evolution can deviate from that of dS in different ways, with distinctive impacts on the particle production history that we are after. In this subsection and the next, we will focus on two simplest ways to deform the dS geometry, namely with a constant parameter and with a constant η parameter. We will call these deviations the -type and the η-type.
The simplest kind of deformation is to introduce a (small) constant parameter, Integrating over τ yields a scale factor where H p ≡ H(τ p ) is the Hubble parameter evaluated at the time τ p when the scale factor is a(τ p ) = (1 − ) −1/ . Because the Hubble parameter is decreasing with conformal time as a power-law, modes that exit the horizon experience a slightly different gravitational background. This soft breaking of scale invariance, as we will see, manifests itself in the scale dependence of production history. Defining the dimensionless variable z = −kτ , the EoM of a massive spin-1 particle reads where we have denoted the scale-dependent dimensionless mass and chemical potential as The complex turning points lie at z c , z * c , with The phase integral is Here F 1 (a; b 1 , b 2 ; c; x, y) is the Appell hypergeometric function defined by the double series

JHEP06(2021)129
where (a) n = Γ(a+n) Γ(a) is the Pochhammer symbol. The Stokes lines are plotted in figure 10. Comparing to figure 8, one finds that with a non-zero , the distribution of turning points becomes asymmetric around the imaginary axis and the Stokes lines become more squeezed with a larger positive chemical potential.
• Production amount. For z i lying on the real axis, the above integral is real, suggesting no contribution to the particle production amount e −2Re F (z i ) , z i ∈ R. Therefore, the particle production amount only receives contribution at the lower end z c . Converting this to an integral along the whole Stokes line, we obtain For k Hp not far away from unity, a small dS-deformation parameter 1 can be used as an expansion parameter. After a further expansion in powers ofκ p and resummation, Therefore, the particle production amount is now scale dependent: where we have resummed the O( 0 ) super-adiabatic corrections to the mass and replaced The O( ) order, however, cannot be treated in the same way as in (4.23), because of the presence of a branch cut extending from the origin all the way to infinity. Therefore, considering the fact that the adiabatic parameter is w (z) w(z) 2 ∼m −2 p , we expect that the relative error of (4.36) is of order H 2 p m 2 , which is negligible if the mass m is large. For a positive , the production amount necessarily drops with scale k, because the effective mass in the Boltzmann factor is measured in units of the time-dependent Hubble parameter, which is decreasing during inflation.
• Production time. This can only be solved numerically from Im F (z * ) = 0. However, by an educated guess, we found an empirical formula that describes the O( ) contribution to z * very well, with an error of 2% on average (see figure 11 for example). Namely, Here the first line (z (0) * ) has the same form as the exact solution found in dS spacetime and the second line (z (1) * ) represents the leading-order slow-roll correction as an empirical formula, with an uncertainty due to the unresummed super-adiabatic corrections.
• Production width. The production width can also be expressed partially analytically. On the real axis, the singulant function accumulates no real parts and therefore JHEP06(2021)129 Re F (z * ) = Re F (z i ). Plugging in the expression for Im F (z * ), we obtain (4.39) When expressed in units of e-folding numbers and expanded into powers of κ m , we have the empirical expression (1) * ) is approximate, with the exception of the coefficient before the running term ln k m , which inherits its exactness from the first line.
The fermionic case is then easily obtained by applying the replacement µ p ,m p → m 2 p +κ 2 p to the above results. Since the zeroth order in is the same as the exact dS result with a Hubble constant H p , we focus on the O( ) corrections due to the deformation and plot them in figure 12. As shown in the plots, the behavior of bosons and fermions are again different for large chemical potentials. production amount excluding zeroth-order contribution as a function of the dimensionless chemical potential for different particle masses. Middle panel: the correction to z-domain production time dependence on chemical potential and mass. Right panel: the correction to production width measured in e-folding numbers. In all three plots, solid lines represent spin-1 bosons while dashed lines represent spin-1/2 fermions. Particles with different masses are distinguished by the colors of the lines according to the legend in the left panel. For bosons, we limit the range of chemical potential to be smaller than the mass, so that no tachyonic instability is induced. For fermions, the chemical potential is not restricted and we allow it to take arbitrarily large values.

JHEP06(2021)129
Before ending this subsection, we note that the validity of our perturbative expansion is actually controlled by ln max{m p ,κ p } instead of just . Thus, for particles with extremely large mass or chemical potential, i.e.,m p ,κ p e 1/ , perturbation theory fails and one would have to rely on the full result (4.35) for production amount and numerically solve production time and width. For parameters chosen in figure 12, the corrections are small compared to the zeroth order dS results, suggesting the perturbative expansion is valid.

Deviation from dS of the η-type
The second type of deformation of dS is obtained from introducing a weak time dependence in (τ ) described by a non-zero η parameter, .
A further simplification appears if (τ ) η(τ ) for τ lying in the range of interest. We will call this type of deviation the η-type. It is reasonable to analyze particle production in such a scenario since for inflation, the Planck 2018 data [64] favors a smaller first slow-roll parameter compared to the second slow-roll parameter.
To study η-type deviation from dS, we suppose the scale factor can be expressed as a series with variable −H i τ , where H i has the dimension of Hubble. So the leading-order deviation can be approximated as with η i > 0 and τ < 0. With this scale factor, the Hubble parameter is and therefore the physical meaning of H i is the Hubble parameter at τ → −∞. On the other hand, the first and second slow-roll parameters are The asymptotic behaviors of the slow-roll parameters in the early-/late-time limit are is computable to the linear order in c(k) 1: where we have resummed the O(c(k) 0 ) super-adiabatic corrections to the mass and replacedm The O(c(k) 1 ) super-adiabatic corrections are complicated by the branch cut from z = 0 to z = −∞ and we are not able to analyze them. Therefore, as in the -deviation case, we have indicated in (4.54) the presence of an O(m −2 ) relative error.
• Production time. The equation Im F (z * ) = 0 can only be solved numerically. However, as in the -type deviation case, we have found a useful empirical formula to the first order in η i , whose relative error is less than 2% within most parameter regions (see JHEP06(2021)129  ) is approximate. 5 • Production width. Using the empirical formula (4.56), this can also be obtained approximately as a truncated power series, (4.57) 5 Note that in this subsection, we will denote X (m,n) as the O(c(k) m η n i ) correction to a certain quantity X. Since X (1,0) follows from an exact dS spacetime with a rescaled Hubble constant, it can always be trivially obtained by rescaling Hi in X (0,0) . In contrast, X (1,1) represents the leading order nontrivial η-corrections.

JHEP06(2021)129
In terms of e-folding numbers, the width is The results for fermions are again obtained via the corresponding replacement. The O(c(k)η i ) corrections to production histories in η-type deformed dS are plotted in figure 16. Comparing to figure 12, one can see that the overall behavior is similar to that of -type deformed dS. Indeed, if we choose ∼ c(k)η i , the corrections roughly match in size. This interesting fact will be discussed below. Nevertheless, we notice that there are important differences in the scale dependence of various production history parameters.
Comparison between -type and η-type. Gravitational particle production histories are crucially influenced by the Hubble parameter H(τ ), as it directly enters the expression of dimensionless mass and chemical potential. Two spacetimes with different H(τ ) are intrinsically different for any process that is non-local in time, including the smoothed Stokes phenomenon. Thus to test this intrinsic difference of particle production for the two types of dS deviations, we carefully select their Hubble parameters to be tangent to each other at time τ 0 , so that both the Hubble and the first slow-roll parameter are equal, For the purpose of demonstration, we will choose the following solution,  The Hubble parameter for the choice (4.61) is shown in figure 17. The production amount and width can be computed using formula given above. Their scale dependence is shown in figure 18. As shown in this figure, particles created near the tangent time τ 0 approximately have the same production amount and width, with mismatches of the same order as the higher-order errors in (4.54) and (4.58).
As a result, one can approximate the η-deformed dS by a tangential -deformed dS locally in time, and obtain the details of particles produced then, up to some higher-order errors. However, this can only be done mode-by-mode, since the scale dependence for these two types of deviations is different. Conversely, if we wish to observationally distinguish JHEP06(2021)129 the two scenarios, we can either accurately measure the production amount/width for a single mode up to some higher orders (e.g., O(c(k) 2 )), or probe the scale dependence by looking at different modes.

Radiation-domination era
Now we turn to other completely different FRW backgrounds, namely those describing the post-inflationary evolution of the universe. Roughly speaking, the scale factor evolves as a(τ ) ∝ τ 2 3w+1 , where w is the equation of state of the dominating component. In this section and the next, we will consider radiation domination era with w = 1 3 and matter domination era with w = 0, respectively. The radiation-dominated universe has a scale factor linearly dependent on the conformal time, a(τ ) = c r τ , where c r > 0 and τ runs from 0 to +∞. The Hubble parameter is H(τ ) = 1 crτ 2 . At the origin lies the Big Bang Singularity, H(0 + ) → ∞. This singularity can be removed by continuously deforming the spacetime to other geometries such as that of inflationary [1][2][3][4], ekpyrotic [65,66], bouncing [67,68] and string gas cosmology [69,70]. These physical continuations provide a cutoff time τ i > 0 with H(τ i ) M p . As we are interested in particle production during the later stage of the radiation-dominated era, we will limit ourselves to the modes with kτ i 1 and assume that they have a vacuum initial condition at τ = τ i ≈ 0 + prepared by the earlier evolution history. 6 Thus for τ τ i , 6 If the initial condition is non-trivial with a non-zero particle number density, one can model this initial particle population by a non-zero βi. The phase of βi may be fixed if the initial particles are prepared coherently, whereas it is random for thermally prepared particles. In such cases, one can return to (3.31), (3.41), (3.81), (3.86), and add to β(τ ) a term proportional to βi and then compute the change of particle number by taking the square and performing an additional ensemble average over the phase of βi if it is thermally prepared.

JHEP06(2021)129
any particle mode with comoving momentum k has a time-dependent physical momentum k a(τ ) = k crτ . This scale is to be compared with H(τ ), m and κ. Defining z = kτ , the EoM now reads with scale-dependent effective chemical potential and mass Their physical meaning is the mass and chemical potential measured in units of Hubble parameter at horizon re-entry. The lower zero of w(z; k) lies at The singulant integral reads (4.65) The super-adiabatic corrections to W (z) yield no simple pole at the origin or at infinity: Therefore, unlike the case of dS in conformal time coordinates, there is no super-adiabatic mass correction and hence no need for resummation. The behavior of the singulant as well as the Stokes multiplier are shown in figure 19. The details of particle production can be easily obtained as follows.
• Production amount. Taking the imaginary part of F on the real axis, we have There are two interesting aspects in this formula. First, |β(k)| 2 is symmetric under the flip κ ↔ −κ (equivalent to a parity transformation that flips helicities), seemingly suggesting an equal enhancement for both helicities. However, this turns out to be superfluous as we will see from inspecting the production history below. Second, for a given mode k, the production amount does not seem to increase with mass monotonically. This phenomenon may be somewhat counterintuitive, as heavier particles naively should be more difficult to produce. This puzzle is resolved when one recalls that the radiation-dominated universe does not have a constant background temperature like dS. It effective "temperature" is likened to the decreasing Hubble parameter H(τ ) = 1 crτ 2 . Raising the mass may lead to two competing effects, one being increasing the difficulty of producing a real particle, the other being pushing the production time earlier, when the effective "temperature" is higher. We will see later that for the radiation-dominated universe, the first effect dominates the applicable range of our method and heavier particles come with a smaller production amount. However, in the matter-dominated universe, this is not the case. Finally, the production amount sharply drops to zero for k 2crm 3 m 2 −κ 2 , for which the "temperature" is too low to support any real particles. and admits a straightforward continuation to the unphysical region −∞ < z < 0. This mathematically continued EoM enjoys a Z 2 symmetry that is well-defined at the origin: 7 κ ↔ −κ, z ↔ −z. This is the cause of the apparent parity symmetry in |β(k)| 2 . In fact, the production time z * is also in the unphysical region for κ < 0. Thus for negative κ, both (4.67) and (4.68) must be understood in the analytically-continued sense. Namely, only if the initial condition at z = 0 + is prepared so as to match the solution of the analytically-continued EoM with Bunch-Davies initial condition at z → −∞, the Stokesline method results are valid. For κ > 0 and 0 < ∆z * 2 z * , these results agree with that of the usual vacuum initial condition at z = 0 + since the Stokes line is far right to the origin and particles do not get produced until a late time. Otherwise, a direct application of the Stokes-line method may be inaccurate. In that case, the particle production history depends on the actual initial condition set at z i = kτ i 1 by an earlier cosmic evolution, which is a physical continuation to the region z < 0.

JHEP06(2021)129
• Production width. The production width in z-domain is simply calculated as Interestingly, the production width in the z-domain for spin-1 boson does not depend on the chemical potential κ. The parameter region where our continuation interpretation matches that of the vacuum initial condition at z = 0 + is where the Ginzburg criterion is satisfied, These conditions actually limitsm < 4 π . The parameter region satisfying the Ginzburg criterion is shown in figure 21, and one can see that bothκ andm are bounded from above, i.e.,κ,m O(1). This suggests that the modes are still relativistic at horizon re-entry, and becomes non-relativistic only after the production time z * .
The fermion case is again obtained by a simple substitutionm → √m 2 +κ 2 . Here one useful check is to go to the large-chemical-potential limit with κ m. There the production amount of fermions reduces to κH(τ * (k)) , (4.71) which is exactly what we expect from the LZ model (2.22). The Ginzburg criterion must also be applied to fermions. Hence, unlike the previous scenarios in dS and its deviations, the chemical potential of fermions is bounded from above (as well as below) by the applicability of our method. We plot the production histories in figure 22. In the valid parameter region, we found that the production amount is monotonically decreasing with mass. Bosons are produced later with larger chemical potential, while the production time of fermions first increase and then decrease with chemical potential. The boson production width is independent of the chemical potential whereas that of fermions decreases with chemical potential, sinceκ enters the expression for effective mass. Due to the imposed constraint, the purple dashed line for fermions withm = 0.5 is absent in all three plots, since it would correspond to artificial initial conditions, as discussed in the main text above.

Matter-dominated era
Now we turn to the matter-dominated era with w = 0 and a scale factor quadratically dependent on the conformal time, a(τ ) = c m τ 2 . The Hubble parameter is H(τ ) = 2 cmτ 3 . Here c m > 0 and τ runs from 0 to +∞. The initial Big Bang singularity is understood to be removed by attaching a period of radiation domination era and some former primordial eras. A vacuum initial condition is still assumed, therefore we will still impose a Ginzburg criterion so that the Stokes-line method gives physical results.
Defining z = kτ , the EoM of a massive vector boson reads with scale-dependent chemical potential and mass The w(z, k) in the matter-dominated universe has four simple zeros on the complex z-plane, namely, ±z c and ±z * c with The z c -z * c pair lies in the right-half plane with positive real parts. Therefore, they are joined by a Stokes line that cross the real axis in the physical region (see the Stokes lines in figure 23 and figure 24). On the other hand, the Stokes line joining −z c and −z * c crosses the real axis in the unphysical region, which can only be understood in the aforementioned continuation sense. • Production time. The production time is solved from Im F (z * ) = 0. The result can be expressed as a power series

JHEP06(2021)129
Although z * is always positive for |κ| <m, the Ginzburg criterion must still be imposed to match the vacuum initial condition.
• Production width. This can also be obtained as a power series,   Another piece of information attainable from the LZ model is an upper bound onκ due to the Ginzburg criterion: in agreement with figure 25.
As before, we plot the production histories for both bosons and fermions in figure 26. The left panel clearly demonstrates the seemingly counter-intuitive mass dependence mentioned previously. Here, the valid region covers the part where the production amount increases with the mass. This shows the fact that production amount decreases with mass in the radiation-domination era is just a coincidence, and that with a time-dependent Hub-JHEP06(2021)129 ble parameter, there are two opposite effects competing against each other, in which case the resulting mass dependence can be subtle. Another notable aspect is that the production time can be either earlier or later than horizon re-entry.

Summary and outlook
Ranging from cosmological collider physics in the primordial era to baryogenesis in the late universe, chemical potential plays an important role in the process of spontaneous creation of particles. In this paper, we focused on the impact of chemical potential on gravitational massive particle production. We first introduced the general form of chemical potential term and gave a necessary condition for its physical effects. After reviewing the chemical potential for particles with different spins, we extracted their essential features and likened the corresponding Bogoliubov coefficients to the coefficients of instantaneous positive/negative frequency solutions. Then the mathematical tools such as asymptotic series, Berry's smoothing techniques of Stokes-lines and Borel resummation were introduced to solve the coefficients. Having checked the applicability of this method at |β| 2 1, we obtain the recipe of particle production histories for both spin-1 bosons and spin-1/2 fermions, which are related by a simple replacement formula. At last, applying this recipe to cosmology, we gave a fine-grained analysis of chemical-potential-assisted particle production in five common FRW spacetimes. The production amount, time and width are obtained as analytic/semi-analytical expressions, each with characteristic dependences on chemical potential and mass.
In summary, our method demonstrates the application of uniformly smoothed Stokesline method to fine-grained particle production. In addition, our results serve as valuable theoretical data for future studies of chemical potential as well as general particle production.
Despite the heavy mathematical machinery and the detailed analysis in this current work, there are still many questions left unanswered which we hope to address in the future. We list a few of them as outlooks below.
• Starting with vacuum initial condition, the introduction of chemical potential invites the interesting possibility of significant particle production with |β| 2 ∼ 1, even with large masses. This mathematically corresponds to the failure of choosing an optimal truncation order n (ReF < 2 for vector bosons and ReF < 1 for fermions) determined by either the stationary phase condition of the first-order perturbation or the minimum term in the asymptotic series solution, and we proposed to use the Borel summation to evaluate the whole divergent asymptotic series, extending the workable parameter regions to ReF 0.5 for bosons and ReF 0.2 for fermions respectively. However, large errors in evaluating the particle production amount were still found when the particle production is too large and runs outside the mentioned workable regions, and part of the reason may be attributed to the failure of approximating w(z) or E(z) around its complex root z c . Although we are currently unable to fully resolve the problem when |β| 2 ∼ 1, it is interesting to note that the result with |β| 2 1, when naively extrapolated JHEP06(2021)129 to the |β| 2 ∼ 1 case, actually gives very accurate answers for the production amount (e.g., in dS and radiation domination era). A systematic method of calculating particle production which can link to the limit with |β| 2 ∼ 1, the tachynonic instability for bosons and the exact Landau-Zener model for fermions (2.21), may require new techniques, and we leave it for future works.
• Throughout the analysis of particle production in section 4, we have assumed a chemical potential constant in space and time. Although this can be justified as leading order approximations, the full understanding can only be acquired by introducing appropriate spacetime dependences according to different contexts.
• In dS and its two types of deviations, it is natural to assume a vacuum initial condition. However, in radiation domination era and matter domination era, quantum fields do not necessarily evolve from the vacuum. In fact, it is expected to have some initial particle population produced in earlier stages of the universe such as inflation or (p)reheating. Yet a non-vacuum initial condition is highly model-dependent. In this work, we choose vacuum initial condition because we focus more on a model-independent analysis of particle production due to a later effect of chemical potential. The treatment of other initial conditions is briefly described in section 4.4. However, it is worthwhile to note some interesting behaviors. If the initial particles are thermally prepared, the interference term in |β(z)| 2 is averaged out by taking the ensemble average over the phase of β i . Thus particle number generally increases due to chemical potential, as expected. In contrast, if the initial particles are coherently prepared with a common phase of β i , chemical potential can serve to produce or destroy particles, depending on the sign of the interference term. If the particle number decreases, one can understand it as the "decay" of particles with energy injection into the background chemical potential sector. It would be interesting to investigate these possibilities with concrete models in the future.
• The knowledge of the production time and width can be helpful in the estimation of signal strength in cosmological collider physics. As mentioned before, the O(|β|) oscillatory signatures on the cosmological collider originate from the interference between the positive frequency part and the negative frequency part, whose presence is controlled by the Stokes multiplier S(z). This fact can be useful when estimating the loop diagrams. At loop level, the momentum integral receives contribution from the UV region with z 1. Usually, this UV divergent part can be regularized by a momentum cutoff at Hubble scale, i.e., k a < H. Then the signal strength follows from dimensional analysis. This is convenient if the mass of the particle running the loop is close to Hubble scale [23]. However, if the particle is much heavier, the dimensionless parameter µ > 1 can enter in complicated ways. Adding chemical potential introduces yet another dimensionless parameterκ, thus invalidating the naive dimensional analysis [25]. However, with the knowledge of particle production history, the momentum cutoff can be posed more precisely at z * . This is because physically speaking, the particles that generate the signals do not get produced until their momentum drops below the production scale, k a < H(z * ± ∆z * ). This potentially offers a better way to estimation signal strength, which deserves further explorations.
• Aside from chemical potential, there are many other sophisticated mechanisms of cosmological particle production, to which the smoothed version of Stokes-line method can be applied. For example, parametric resonance is widely used in models of preheating [71][72][73], generation of primordial black holes [74,75] and primordial gravitational waves [76,77]. For a periodic effective frequency w(z), the turning points form periodic pairs on the complex z-plane, joined by periodic Stokes lines. Then the resonance condition can be viewed as the constructive interference of particle production amplitudes when crossing each Stokes line. In the literature, there are already preliminary attempts in this direction [50,52], but using the traditional Stokes-line method without uniform smoothing. This will be accurate if the Stokes lines are well-separated so that one can apply the "dilute gas" approximation, treating each crossing separately as sudden jumps in particle number. However, if the production widths are as wide as the separation between two neighboring Stokes lines, one may need to go to the fine-grained picture and perform the analysis using the smoothed version of Stokes-line method. It is interesting to compare this method with traditional ones such as the Floquet theory. (1 + c(k)) (0.6m + 0.3κ + · · · ) 1 − c(k)  Table 1. A checklist of results for spin-1 vector bosons. Here RD and MD stand for radiation domination era and matter domination era, respectively. m is the mass while κ is the chemical potential. The =, , ≈ symbols are used to indicate whether the result is exact, numerically exact, or empirical with 2% error. In the schematic expressions, to display the most salient features, we have omitted the Hubble parameter and blurred the difference between m H and µ = m 2 H 2 − 1 4 in dS. The widths in dS, -type dS and η-type dS are measured in e-folds, whereas the widths in RD and MD are measured in z-domain. The results for spin-1/2 fermion is obtained via a simple replacement rule m 2 → m 2 +κ 2 . Setting κ = 0 also gives the purely gravitational production results. Now consider a closed integration path C 1 ∪ C 2 ∪ C 3 ∪ C 4 that tours around the branch cut. The integral can be converted to a residue at z → ∞, Separating the full frequency into the zeroth order and higher orders, we have

JHEP06(2021)129
where we have used the fact that δW (n) (z) drops as at infinity and therefore does not contribute to the residue. Then we separate the closed contour into two parts, 2πiκ = Across the branch cut, the phase of W (z) jumps by π and its modulus remains continuous. Therefore, along the branch cut, the integral on C 1 and C 3 gives the same result. The arcs around z c and z * c on C 1 and C 3 , however, must be taken with a grain of salt. As these two singularities are of high orders, the integral there is ill-defined when the radius of the arc goes to zero. As a result, we need to manually impose a principal value prescription. By symmetry, the function H(θ) ≡ i If we choose θ m to be the angle of the meeting point of C 1 and C 3 , then the integrals along them give the same result:

JHEP06(2021)129
In addition, due to the sign flip across the branch cut, the integral along C 2 and C 4 are related by