Superadiabatic basis in cosmological particle production: application to preheating

We review a commonly used definition of adiabatic vacuum, and will discuss its ambiguity in the next subsection. In this section, we will consider a scalar field with a time-dependent mass,

$\displaystyle\mathcal{L}=-\frac{1}{2}(\partial\chi)^{2}-\frac{1}{2}m^{2}(t)\chi^{2},$

(1)

and we Fourier expand the scalar field as

$\displaystyle\hat{\chi}(t,{\bf x})=\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}e^{{\rm i}{\bf k}\cdot{\bf x}}\left[\hat{a}_{\bf k}v_{k}(t)+\hat{a}_{-{\bf k}}^{\dagger}v^{*}_{k}(t)\right],$

(2)

where $v_{k}(t)$ is a mode function, $\hat{a}_{\bf k}$ and $\hat{a}_{\bf k}^{\dagger}$ are creation and annihilation operator, respectively, which satisfy $[a_{\bf k},a^{\dagger}_{{\bf k}^{\prime}}]=\delta^{3}({\bf k}-{\bf k}^{\prime})$. The mode function $v_{k}(t)$ should follow the normalization condition,

$$v_{k}\dot{v}_{k}^{*}-\dot{v}_{k}v_{k}^{*}=-{\rm i},$$

(3)

at any time $t$. The mode equation is given by

$$\ddot{v}_{k}+\omega_{k}^{2}v_{k}=0,$$

(4)

where $\omega_{k}^{2}=m^{2}(t)+k^{2}$.

The WKB method is commonly used in order to estimate the particle production rate, where we choose the positive and the negative basis function to be

$$f_{k}^{\pm}=\frac{1}{\sqrt{2\omega_{k}}}e^{\mp{\rm i}\int_{t_{0}}^{t}\omega_{k}(t^{\prime})dt^{\prime}},$$

(5)

where $t_{0}$ is a reference time often taken to be $t_{0}\to-\infty$. With these basis functions, we can formally find the mode function to be

$$v_{k}(t)=\alpha_{k}(t)f_{k}^{+}+\beta_{k}(t)f_{k}^{-}.$$

(6)

Here, the auxiliary functions $\alpha_{k}(t)$ and $\beta_{k}(t)$ satisfy

	$\displaystyle\dot{\alpha}_{k}(t)$	$\displaystyle=$	$\displaystyle\frac{\dot{\omega}_{k}}{2\omega_{k}}e^{2{\rm i}\int_{t_{0}}^{t}\omega_{k}dt^{\prime}}\beta_{k}(t),$		(7)
	$\displaystyle\dot{\beta}_{k}(t)$	$\displaystyle=$	$\displaystyle\frac{\dot{\omega}_{k}}{2\omega_{k}}e^{-2{\rm i}\int_{t_{0}}^{t}\omega_{k}dt^{\prime}}\alpha_{k}(t),$		(8)

so that the mode function (6) be the solution of the mode equation (4). Assuming $m^{2}(t)\to\text{const}$ as $t\to\pm\infty$, we are able to define “particles” in both the future and past infinity.⁴⁴4We will consider the backgrounds that are not asymptotically static. In this sense, it would be better to call thus defined “vacuum” as instantaneous adiabatic vacuum.

Within this formalism, the specification of the vacuum corresponds to the choice of the boundary conditions on $\alpha_{k}$ and $\beta_{k}$. The past adiabatic vacuum condition is

$$\alpha_{k}(-\infty)=1,\beta_{k}(-\infty)=0.$$

(9)

From the future observer, the number density of the particle produced by the time variation of the mass is given by $n_{k}=|\beta_{k}(\infty)|^{2}$. In most cosmological models, it is not always possible to assume $m(t)^{2}\to{\rm const}$ as $t\to-\infty$. In such a case, we need to impose conditions like (9) at some reference time $t_{0}$. As long as the total number density $n(t)=\int d^{3}k|\beta_{k}(t)|^{2}$ is finite, thus defined instantaneous vacuum can have overlap with the vacuum at different time $t$. In other words, if $n(t)$ were not finite, such a vacuum state is not well defined Parker:1969au ; Fulling:1979ac .

When we only discuss the particle number at the future and past infinity where the particle number is definite, the commonly used WKB method is sufficient to estimate the produced particle number. However, if we are interested in the particle number at an intermediate time, we should be careful: One may think that $n_{k}(t)=|\beta_{k}(t)|^{2}$ derived by solving (7) and (8) describes the particle number at a given time $t$. However, the “particle number” at intermediate time depends on the choice of the basis function $f_{k}^{\pm}$. As we will discuss below, the choice of the basis function is not unique, and in turn, the auxiliary functions $\alpha_{k}(t)$ and $\beta_{k}(t)$ are not uniquely determined.

Following Dabrowski:2014ica ; Dabrowski:2016tsx , we review more general definition of adiabatic vacuum. We use the following ansatz of the mode function

$$v_{k}=\frac{1}{\sqrt{2W_{k}}}e^{-{\rm i}\int_{t_{0}}^{t}W_{k}dt^{\prime}},$$

(10)

where the function $W_{k}(t)$ is an unspecified function. Substituting this ansatz to (4), we find the condition on $W_{k}(t)$

$$W_{k}^{2}(t)=\omega_{k}^{2}(t)-\left[\frac{\ddot{W}_{k}}{2W_{k}}-\frac{3}{4}\left(\frac{\dot{W}_{k}}{W_{k}}\right)^{2}\right].$$

(11)

Finding the solution $W_{k}(t)$ is equivalent to solve the mode equation exactly. However, it is impossible in general. Therefore, we will instead consider approximations: Let us assume adiabatic behavior of $\omega_{k}$ namely $|\dot{\omega}_{k}/\omega_{k}^{2}|\ll 1$, which may be violated at some region. Assuming adiabaticity, we are able to solve the equation by iteration as

$$(W_{k}^{(j+1)})^{2}=\omega_{k}^{2}-\left[\frac{\ddot{W}_{k}}{2W_{k}}-\frac{3}{4}\left(\frac{\dot{W}_{k}}{W_{k}}\right)^{2}\right]\Biggr{|}_{W_{k}=W_{k}^{(j)}},$$

(12)

where $j=0,1,2,\cdots$ with $W_{k}^{(0)}=\omega_{k}$ as the lowest solution. At the $(j+1)$-th order, one has to take the terms containing at most $2(j+1)$ time derivatives. For example, the first order is given by

$$(W_{k}^{(1)})^{2}=\omega_{k}^{2}-\left[\frac{\ddot{\omega}_{k}}{2\omega_{k}}-\frac{3}{4}\left(\frac{\dot{\omega}_{k}}{\omega_{k}}\right)^{2}\right].$$

(13)

This is known as adiabatic expansion, which is often used to regularize the UV divergence in time dependent background such as expanding Universe Zeldovich:1971mw ; Parker:1974qw . This can be thought of the generalization of the WKB approximation, which corresponds to the 0-th order in the adiabatic expansion. It is known that the adiabatic expansion is asymptotic series, namely the series does not converge. Therefore, we need to truncate the series at the optimal order. The optimal order can be estimated by the phase integral along the turning points as shown in appendix A. For derivation of the estimate, we refer Barry:1989zz ; Li:2019ves ; Sou:2021juh .

For any adiabatic order, we are able to construct the solution of the mode equation with the help of auxiliary functions as we showed in the previous section. Let us take the $j$-th order, and we formally write the mode function as

$$v_{k}(t)=\alpha_{k}(t)f^{+(j)}_{k}+\beta_{k}(t)f_{k}^{-(j)},$$

(14)

where

$$f_{k}^{\pm(j)}\equiv\frac{1}{\sqrt{2W_{k}^{(j)}}}\exp\left[\mp{\rm i}\int_{t_{0}}^{t}W_{k}^{(j)}dt^{\prime}\right]$$

(15)

are the $j$-th order basis functions and $\alpha_{k}(t)$ and $\beta_{k}(t)$ are unspecified auxiliary functions, which play the role of the Bogoliubov coefficients. We require the first time-derivative of $v_{k}$ to be

$$\dot{v}_{k}(t)=(-{\rm i}W_{k}^{(j)}+V_{k})\alpha_{k}f_{k}^{+(j)}+({\rm i}W_{k}^{(j)}+V_{k})\beta_{k}f_{k}^{-(j)},$$

(16)

where $V_{k}(t)$ is a real function of time. We can easily check that the normalization condition $v_{k}\dot{\bar{v}}_{k}-\dot{v}_{k}\bar{v}_{k}={\rm i}$ is satisfied if $|\alpha_{k}(t)|^{2}-|\beta_{k}(t)|^{2}=1$ independently of the real function $V_{k}(t)$. Therefore, $V_{k}(t)$ is another ambiguity of the definition of the basis function. The “natural choice” of this degree of freedom is proposed in Dabrowski:2016tsx as $V_{k}=-\frac{\dot{W}_{k}^{(j)}}{2W_{k}^{(j)}}$.

In order for $\alpha_{k}(t)$ and $\beta_{k}(t)$ to satisfy (16) and the mode equation, the evolution equations of $\alpha_{k}$ and $\beta_{k}$ should be

$$\left(\begin{array}[]{c}\dot{\alpha}_{k}(t)\\ \dot{\beta}_{k}(t)\end{array}\right)=\left(\begin{array}[]{cc}\delta_{k}&(\Delta_{k}+\delta_{k})e^{2{\rm i}\int_{t_{0}}^{t}W_{k}^{(j)}dt^{\prime}}\\ (\Delta_{k}-\delta_{k})e^{-2{\rm i}\int_{t_{0}}^{t}W_{k}^{(j)}dt^{\prime}}&-\delta_{k}\end{array}\right)\left(\begin{array}[]{c}\alpha_{k}(t)\\ \beta_{k}(t)\end{array}\right),$$

(17)

where

	$\displaystyle\delta_{k}$	$\displaystyle\equiv$	$\displaystyle\frac{1}{2{\rm i}W_{k}^{(j)}}(\omega_{k}^{2}-(W^{(j)}_{k})^{2}+\dot{V}_{k}+V_{k}^{2}),$		(18)
	$\displaystyle\Delta_{k}$	$\displaystyle\equiv$	$\displaystyle\frac{\dot{W}_{k}^{(j)}}{2W_{k}^{(j)}}+V_{k}.$		(19)

Here, we have left $V_{k}$ in a general form, and the “standard” WKB formula discussed in the previous section corresponds to the choice $V_{k}(t)=0$ with $W_{k}^{(0)}=\omega_{k}(t)$. The natural choice $V_{k}=-\frac{\dot{W}_{k}^{(j)}}{2W_{k}^{(j)}}$ simplifies the expression because $\Delta_{k}=0$. The simplified differential equations for the Bogoliubov coefficients are

$$\left(\begin{array}[]{c}\dot{\alpha}_{k}(t)\\ \dot{\beta}_{k}(t)\end{array}\right)=\delta_{k}\left(\begin{array}[]{cc}1&e^{2{\rm i}\int_{t_{0}}^{t}W_{k}^{(j)}dt^{\prime}}\\ -e^{-2{\rm i}\int_{t_{0}}^{t}W_{k}^{(j)}dt^{\prime}}&-1\end{array}\right)\left(\begin{array}[]{c}\alpha_{k}(t)\\ \beta_{k}(t)\end{array}\right),$$

(20)

and one can easily confirm that the condition $|\alpha_{k}(t)|^{2}-|\beta_{k}(t)|^{2}=1$ is satisfied at any time $t$.

As explicitly seen from the discussion so far, the definition of the basis function or equivalently $(\alpha_{k}(t),\beta_{k}(t))$ is indeed ambiguous. We may interpret this ambiguity as follows. The choice of basis functions is the definition of energy quantum, which accordingly defines a “particle”. Since the time translation invariance is broken by the background field, we have no unique choice of the energy quanta defining particles. Nevertheless, the optimal choice of the “energy” quantum would give us the optimal definition of a “particle”. We may summarize the construction so far as follows. First, we need to define the energy quantum, namely, need to choose the adiabatic order of $W_{k}^{(j)}$. This is yet insufficient to define a particle. The ambiguity comes from the first order derivative $\dot{\chi}$. We have a real function degree of freedom $V_{k}(t)$, which determines the first order differential equation of the auxiliary functions $\alpha_{k}(t)$ and $\beta_{k}(t)$. Once we fix the real function $V_{k}(t)$, an “adiabatic particle” is defined. The adiabatic vacuum state corresponds to the boundary conditions for $\alpha_{k},\beta_{k}$.

We also note that in Dabrowski:2016tsx it is shown that there are several approaches to obtain adiabatic particle spectra, Klein-Gordon equation, Ermakov-Milne equation, Riccati equation, and spectral function approach. The authors of Dabrowski:2016tsx show the equivalence between them. In the following discussion, we will take the Klein-Gordon equation approach, which is more commonly used to discuss e.g. particle production.⁵⁵5For numerical simulations, Ermakov-Milne equation seems useful, and we review it in appendix B. Nevertheless, the equivalence shown in Dabrowski:2016tsx ensures any approach gives the same result.

In this section, we briefly review the Stokes phenomenon and its relation to particle production. We refer Froman1 for more details of Stokes phenomenon and WKB/phase integral methods (see also Hashiba:2021npn ). Instead of a general review of Stokes phenomena, we will explain how to evaluate the connection formula associated with Stokes phenomenon in the next section within simple examples. The explanation here is rather abstract.

The particle production caused by time-dependent backgrounds can be understood from asymptotic behavior of mode functions. In order to understand the asymptotic behavior, the WKB method or more generally the phase integral method is quite useful. For a given mode equation, the analytic property of solutions is characterized by $\omega_{k}^{2}(t)$. The behavior of the solution on complex $t$-plane significantly changes around turning points $t_{c}$ at which the frequency vanishes $\omega_{k}^{2}(t_{c})=0$. From turning points, Stokes lines emanate, on which ${\rm i}\int^{t}\omega_{k}dt$ is real. On the stokes line, the mode function $e^{\pm{\rm i}\int\omega_{k}dt}$ increases or decreases significantly. Such lines ends up at poles of $\omega_{k}$, different turning points or infinity. Thus, Stokes lines separate the complex plane into several pieces. When we consider analytic continuation of basis functions from one region to another separated by Stokes lines, the analytic continuation leads to non-trivial mixing between positive and negative frequency modes, which may give rise to nonvanishing $\beta_{k}$ even if we take $\beta_{k}=0$ namely an adiabatic vacuum condition at an initial time. See Fig. 1. The nontrivial mixing between positive and negative frequency basis functions is what we call Stokes phenomenon.⁶⁶6In Fig. 1, we describe the Stokes phenomenon by multiplication of matrices on $\alpha_{k},\beta_{k}$. This is equivalent to mixing of positive and negative frequency modes. Even if turning points do not appear on real $t$-axis, Stokes lines cross real $t$-axis and therefore Stokes phenomenon, or equivalently, “particle production” necessarily takes place.

The adiabatic expansion is in general an asymptotic series, namely the series does not converge. Therefore, it is not mathematically well-defined. In order to have a meaningful solution, we may consider Borel transformation of adiabatic series, which gives well-defined convergent solutions. Such an analysis is known as exact WKB analysis. (See Taya:2020dco and references therein.) From exact WKB analysis viewpoints, the Stokes line corresponds to the place where the series are not Borel summable, and therefore the exact WKB solution is discontinuous on Stokes lines. Therefore, one needs to consider analytic continuations of solutions when crossing Stokes lines. The analytic continuation gives non-trivial mixing of positive and negative frequency solutions as mentioned above.

Berry applied partial Borel resummation technique to later terms of the adiabatic expansion and showed a smooth connection formula of $\beta_{k}$ around the Stokes line crossing Barry:1989zz , which turn out to be a good analytical formula describing the “time-dependent particle numbers” Dabrowski:2014ica ; Dabrowski:2016tsx . We summarize the formula in Appendix A without derivations, and refer to Berry’s original paper Barry:1989zz for the derivation. (See also recent papers Li:2019ves ; Sou:2021juh .)

Let us consider the following effective frequency

$$\omega_{k}^{2}(t)=k^{2}+m^{2}+\lambda\mu^{4}t^{2},$$

(28)

where $\mu$ is a mass parameter. This simple model captures the broad resonance regime of preheating and the moving D-brane model Kofman:1997yn ; Kofman:2004yc . Indeed, if we take $h(\phi)=\phi^{2}$ in (21) and approximate the field dynamics as $\phi\sim\mu^{2}t$ while ignoring expansion of the Universe, we find (28).

In order to estimate the amount of particle production, we need to discuss Stokes phenomenon in this setup. The turning points, at which $\omega_{k}^{2}=0$, can be found at

$$t_{c}={\rm i}\frac{\Omega_{k}}{\sqrt{\lambda}\mu^{2}},$$

(29)

and its conjugate $\bar{t}_{c}$, where $\Omega_{k}=\sqrt{k^{2}+m^{2}}$. Since $\omega_{k}(t)$ is a real function of $t$ on the real axis, the Stokes line crossing the real axis at $t=0$ connects two turning points. Note that three Stokes lines emanate from each turning point, and two of them do not cross the real axis and go to infinity. (See Fig. 2.)

Using the Dingle-Berry formula (58), we can estimate the time dependence of the particle number when we take the optimal adiabatic order. Let us explicitly show the quantities necessary to obtain an analytic formula in the following. The amount of the particle production is determined by the phase integral along the Stokes line connecting two turning points (61), which is given by

$$F_{k}^{(0)}={\rm i}\int_{t_{c}}^{\bar{t}_{c}}\omega_{k}dt=\frac{\pi\Omega_{k}^{2}}{2\sqrt{\lambda}\mu^{2}},$$

(30)

where we have chosen the phase of $\omega_{k}$ so that $F_{k}^{(0)}>0$. The Stokes line crosses real axis at $t=s_{c}=0$. Then, we can approximate $F_{k}(t)$ (60) as

	$\displaystyle F_{k}(t)$	$\displaystyle=$	$\displaystyle 2{\rm i}\left(\int_{t_{c}}^{s_{c}}+\int_{s_{c}}^{t}\right)\omega_{k}(t^{\prime})dt^{\prime}$		(31)
		$\displaystyle=$	$\displaystyle F_{k}^{(0)}+2{\rm i}\int_{s_{c}}^{t}\omega_{k}(t^{\prime})dt^{\prime}$
		$\displaystyle\sim$	$\displaystyle F_{k}^{(0)}+2{\rm i}\omega_{k}(s_{c})(t-s_{c})$
		$\displaystyle=$	$\displaystyle\frac{\pi\Omega_{k}^{2}}{2\sqrt{\lambda}\mu^{2}}+2{\rm i}\Omega_{k}t.$

This expression is valid around $t\sim s_{c}(=0)$. With this expression, we find $\sigma_{k}$ (59) as

$$\sigma_{k}(t)=\frac{{\rm Im}F_{k}(t)}{\sqrt{2{\rm Re}F_{k}(t)}}\sim\frac{2\lambda^{1/4}\mu t}{\sqrt{\pi}}.$$

(32)

Thus, for the optimal truncation order, $\beta_{k}(t)$ would behave as

$$\beta_{k}(t)\sim\frac{\rm i}{2}e^{-\frac{\pi\Omega_{k}^{2}}{2\sqrt{\lambda}\mu^{2}}}{\rm Erfc}\left(-\frac{2\lambda^{1/4}\mu t}{\sqrt{\pi}}\right),$$

(33)

and also the time-dependent particle number is given by

$$n_{k}\sim\frac{1}{4}e^{-\frac{\pi\Omega_{k}}{\sqrt{\lambda}\mu^{2}}}\left({\rm Erfc}\left(-\frac{2\lambda^{1/4}\mu t}{\sqrt{\pi}}\right)\right)^{2}.$$

(34)

We find that this simple formula actually captures the behavior of the particle number evolution at the optimal order with the natural choice $V_{k}=-\frac{\dot{W}_{k}^{(j)}}{2W_{k}^{(j)}}$. We note that the form of $\beta_{k}$ in (33) actually exhibits its non-perturbative property in the coupling constant $\lambda$. As we expected, $\beta_{k}$ decays exponentially as $k$ increases, and does not lead to UV divergences.

For illustration, we show the numerically derived time-dependent particle number at different adiabatic orders with either $V_{k}=-\frac{\dot{W}_{k}^{(j)}}{2W_{k}^{(j)}}$ or $V_{k}=0$ in Fig. 3 (for natural choice of $V_{k}$), and Fig. 4 (for $V_{k}=0$).

As we clearly see from the figures, the first adiabatic expansion with the natural choice of $V_{k}$ shows a smooth transition, which is in good agreement with the approximate formula (34). We also confirm that the asymptotic value of $n_{k}$ converges to the same, irrespective of the adiabatic order or the choice of $V_{k}$. Therefore, as long as we are interested in the final values of produced particle number, the different choices of either adiabatic orders $j$ or the choice of $V_{k}$ would not matter. Nevertheless, it would be appropriate to define a particle with $j=1$ and the natural choice of $V_{k}$ within this parameter set, which yields “stable” particle number. We note that, according to Berry’s estimate of the optimal order (62), our parameter choice is optimal at $j=1$, which is consistent with our result.¹³¹³13We should note that we find the numerical value $\frac{1}{2}(F^{(0)}_{k}-1)\sim 1.96$, which marginally indicates $j=1$ as the optimal order. This might be because the value of $F_{k}^{(0)}$ is not so large for our parameter choice, whereas Berry’s derivation of the optimal order assumes a large value of $F_{k}^{(0)}$. We may understand the small discrepancy between (34) and the first adiabatic order result in Fig. 3 in the same way.

Thus, we have confirmed that with an appropriate choice of the adiabatic order and $V_{k}$, we can have a proper adiabatic particle, which shows (1) smooth particle production around Stokes line crossing, (2) stable (non-oscillatory) behavior of particle number. Besides that, the proper adiabatic particle number or $\beta_{k}(t)$ is well approximated by a very simple formula (34). It is worth emphasizing that if we can have an analytic approximation for $\beta_{k}$, we are able to calculate physical quantities such as energy density (see (26)). Note that $\alpha_{k}$ should satisfy the normalization condition $|\alpha_{k}|^{2}-|\beta_{k}|^{2}=1$, and once $|\beta_{k}|$ is fixed, the absolute value of $\alpha_{k}$ can be determined. Therefore, we would be able to have approximate formulas for various physical quantities with a simple formula (34).

In the following, we will reconsider the cases where multiple particle production events are relevant and interference of events affect the time-dependent particle number.

We discuss the production of a scalar particle $\chi$ coupled to an oscillating scalar field $\phi$ in Minkowski background, which is a toy model of preheating after inflation. More explicitly we consider $h(\phi)=\phi(t)=A\cos(Mt)$, where $A$ is a dimensionful constant. Then, the effective frequency of $\chi$ for a mode with momentum $k$ is given by

$$\omega_{k}^{2}=k^{2}+m^{2}+\lambda A\cos(Mt).$$

(35)

Depending on the set of parameters, we have three different regime: For $\lambda A>k^{2}+m^{2}$, there are time intervals when the effective frequency becomes pure imaginary, namely a tachyonic instability regime. We will not discuss this regime here.¹⁴¹⁴14In Hashiba:2021npn , Stokes phenomenon in the tachyonic regime is discussed with in $R^{2}$ inflation model. In this case, the Stokes line appears along the real time axis. See also Dufaux:2006ee . For $\lambda A\ll k^{2}+m^{2}<M^{2}$, the modes in a certain momentum range experiences the so-called narrow resonance regime, which will be explained below. For $M^{2}\ll\lambda A<k^{2}+m^{2}$, we find the broad resonance regime. The broad resonance, in particular, cannot be understood from a perturbative scattering viewpoint, since the decay of $\phi$ to $\chi$ is kinematically forbidden.

Let us discuss the narrow resonance case $\lambda A\ll k^{2}+m^{2}<M^{2}$. For most momentum modes, the amount of particle production is quite small. Nevertheless, the modes in the resonance band around $4(k^{2}+m^{2})\sim M^{2}$ shows particle number exponentially growing in time. One of the reasons to reconsider the narrow resonance regime here is for the following question: Can we describe the narrow resonance in the same way as the broad resonance? It is well known that resonance in this model can be understood in terms of Mathieu equation: The mode equation of $v_{k}$ is rewritten as

$$v^{\prime\prime}_{k}+(A_{k}-2q\cos(2z))v_{k}=0,$$

(36)

where $A_{k}=4(k^{2}+m^{2})/M^{2}$, $q=2\lambda A/M^{2}$, and $z=Mt/2+\pi$. The first resonance band condition is $1-q<A_{k}<1+q$ Kofman:1997yn . The center of the resonance band $k^{2}+m^{2}\sim M^{2}/4$ seems similar to the possible kinematic configuration of $\phi\to\chi\chi$ decay channel. It is known that the mode function inside the resonance band grows as Kofman:1997yn

$$v_{k}(t)\propto\exp\left(\frac{1}{2}qz\right)\propto\exp\left(\frac{\lambda At}{2M}\right),$$

(37)

and therefore the particle number would grow as $n_{k}\propto\exp(\lambda At/M)$. Therefore, we interpret $\Gamma=\lambda A/M$ as the production rate of $\chi$. The narrow resonance can be understood from the perturbative decay of inflaton to $\chi$ with Bose enhancement effects (see e.g. Mukhanov:2005sc ).¹⁵¹⁵15We note however that the production rate is quite different from that given by a perturbative analysis. For clarification of this point, we give a perturbative analysis of $\chi$-production in appendix C.

On the other hand, in the broad resonance regime, produced particle number is expressed in terms of non-analytic functions of a coupling constant $\lambda$ as it can be understood from the simplified model in Sec. 4.1. If we try to describe both the narrow and broad resonance regime, it would mean that we need to treat perturbative and non-perturbative contributions simultaneously. We will show that this is possible, namely, we are able to describe both regime in a unified way. This is not so surprising since the both regime can be described by a single differential equation, Mathieu equation. Nevertheless, it is important to confirm that Stokes phenomenon captures both perturbative and non-perturbative effects.

Let us discuss the Stokes phenomenon in this model. Turning points at which $\omega_{k}^{2}(t_{c})=0$ are

$$t_{c}=\frac{1}{M}\left((2n-1)\pi\pm{\rm i}\log(c+\sqrt{c^{2}+1})\right),$$

(38)

where $c=(k^{2}+m^{2})/(\lambda A)$, and $n$ is an arbitrary integer. The appearance of $n$ represents the repeated particle production events, and we will take $n=1$ to be the first event. A turning point and its complex conjugate are connected by a Stokes line that crosses the real $t$ axis. We illustrate the Stokes line structure in this model in Fig. 5. The amount of the particle produced at each Stokes line crossing is determined by

$$F^{(0)}_{k}={\rm i}\int_{\bar{t}_{c}}^{t_{c}}dt\omega_{k}=-4{\rm i}\frac{\sqrt{\lambda A}}{M}\sqrt{c-1}E\left(\frac{\rm i}{2}\log(c+\sqrt{c^{2}+1})\bigl{|}-\frac{2}{c-1}\right),$$

(39)

where $E(a|x)=\int_{0}^{a}\sqrt{1-x^{2}\sin^{2}\theta}d\theta$ is the incomplete elliptic integral of the second kind. The interference effect among each Stokes line crossing is characterized by the phase shift

$$\theta=\int_{0}^{2\pi}\omega_{k}dt=\frac{2\sqrt{k^{2}+m^{2}}}{M\sqrt{c}}\left(\sqrt{c-1}E\left(-\frac{2}{c-1}\right)+\sqrt{c+1}E\left(\frac{2}{1+c}\right)\right),$$

(40)

where $E(x)=E(\frac{\pi}{2}|x)$ is the complete elliptic integral of the second kind. Note that for the large $c$ limit, we find $\theta\to\frac{2\pi\sqrt{k^{2}+m^{2}}}{M}$, and for $c\to 1$, $\theta\to 0$. Therefore, at the first resonance band $k^{2}+m^{2}\sim M^{2}/4$, the phase factor becomes $\theta\sim\pi$. Our approximate connection matrix (69) at $n$-th Stokes line crossing becomes

$$C_{n}=\left(\begin{array}[]{cc}\sqrt{1+|E_{k,n}(t)|^{2}}&{\rm i}E_{k,n}(t)\\ -{\rm i}(E_{k,n}(t))^{*}&\sqrt{1+|E_{k,n}(t)|^{2}}\end{array}\right),$$

(41)

where

$$E_{k,n}(t)=\frac{1}{2}e^{2{\rm i}n\theta}e^{-F^{(0)}_{k}}{\rm Erfc}(-\sigma_{k,n}(t))$$

(42)

and

$$\sigma_{k,n}(t)\sim\frac{2\sqrt{k^{2}+m^{2}-\lambda A}}{\sqrt{2F_{k}^{(0)}}}(t-(2n-1)\pi).$$

(43)

Note that due to the periodicity of the background, the phase factor is simply given by the multiple of $e^{{\rm i}\theta}$. One can obtain $(\alpha_{k}(t),\beta_{k}(t))$ including particle production up to the $n$-th event as

$$\left(\begin{array}[]{c}\alpha_{k}(t)\\ \beta_{k}(t)\end{array}\right)=\prod_{j=1}^{n}C_{j}\left(\begin{array}[]{c}1\\ 0\end{array}\right),$$

(44)

where we have assumed the adiabatic vacuum condition $(\alpha_{k}(0),\beta_{k}(0))=(1,0)$. One may use the simpler approximation formula (63) if $n_{k}<1$, but it fails for the case $n_{k}>1$. The formula (44) is all we need for evaluation of the particle production. One can immediately obtain the particle number $n_{k}=|\beta_{k}|^{2}$ from the formula (44).

We give some general remarks about our approximate formula (41). The formula shows time-dependence in good agreement with numerically solved mode equation in various parameter regions. However, it fails to reproduce the numerical result exactly. There are several reasons for the discrepancy. Our conjectured formula (69) is based on the combination of Dingle-Berry’s error function formula and the connection formula of the parabolic cylinder function. This approximation would be applicable if each Stokes line crossing the real axis is well separated from others. If we try to improve the connection matrix formula, we need to perform more careful analysis of Stokes phenomenon using e.g. symmetry relations Froman1 ; Kutlin:2017dpe .¹⁶¹⁶16It would also be important to modify Berry’s derivation of the error function formula in order to find the time dependence of the Bogoliubov coefficients $\alpha_{k}(t),\beta_{k}(t)$. To our best knowledge, the generalization of Berry’s formula in the presence of multiple turning points is not known. Nevertheless, our simple formula (69) would be useful for practical purposes and reproduce correct behavior in most cases.

Let us show the comparison between $n_{k}=|\beta_{k}(t)|^{2}$ with our approximate formula (44) and that of numerical solutions. In evaluating particle number numerically, we used Emarkov-Milne equation reviewed in appendix B. In Fig. 6, we show the narrow resonance case. We see plateau regions for the numerical solution, but such regions do not appear in our approximate formula. This is because the width of the error function is too large and separated particle production events are smoothly connected. Nevertheless, we see the growth behavior consistent with the numerical solution, up to an $\mathcal{O}(1)$ overall factor. We show a simple modification in appendix D, which shows a better fitting. We also show the long time behavior of the narrow resonance regime in Fig. 7 (with a slightly different parameter set). We find the exponential growth rate expected from the instability band analysis is approximately reproduced. We would like to clarify the origin of the growth. As we mentioned earlier, the phase factor $\theta\sim\pi$ in this regime, and therefore, the phase factor appearing in the connection matrix does not give any interference effects. Effectively, after $n$ events, we just multiply the same matrix for $n$ times. One can easily confirm that such multiplication of matrices leads to exponential growth of $|\beta_{k}|$. This is the reason for the narrow resonance from our viewpoint. Such a viewpoint seems different from the instability analysis of Mathieu equation or Floque theory Amin:2014eta .

In Fig. 8, the broad resonance case is shown. In this case, we find a very good agreement between the numerical solution and our approximate formula. We also show a non-resonant case in Fig. 9. We find that our approximate formula looks an envelope of the numerical result. Thus, our simple formula is in good agreement with numerical evaluation of particle number in different parameter regions. However, we should emphasize a possible problem in our approximation. In the case of resonant particle production, the particle number grows exponentially, and therefore, small deviation at earlier time could be enhanced. Indeed, as shown in Fig. 7, the growth rate of our analytic approximation shows a slightly larger value than that of numerical result or expectation from the instability band analysis. This discrepancy results in large deviation due to the exponential dependence. As we quoted in appendix A, it is not so easy to improve this point. Therefore, one needs to be careful about such an issue when applying our analytic formulas to models of explosive particle production. A similar issue appears in the preheating in expanding universe. As we will see in the next section, it is quite difficult to determine the interference effects among number of events. Therefore, we are able to find at most the estimate with a stochastic interference effects, which is known from earlier work Kofman:1997yn . Independently of these issues, one has to take into account backreaction to the background field, which cannot be neglected if the particle number is exponentially large. Then, one needs to perform e.g. lattice simulations Khlebnikov:1996mc ; Khlebnikov:1996wr ; Prokopec:1996rr ; Khlebnikov:1996zt to take the backreaction into account. Even if one tries to use Stokes phenomenon analysis, the backreaction needs to be taken into account since it changes the background dynamics and accordingly the structure of Stokes lines. Then, the interference effects would become stochastic one. Therefore, a simple exponential growth would not be realized in reality.

It is also important to check the UV behavior of $\beta_{k}(t)$. As we showed, the absolute value of $\beta_{k}$ is characterized by $e^{-F_{k}^{(0)}}$. Let us see how $F_{k}^{(0)}$ behaves for high momentum region, which corresponds to $c\to\infty$. From the definition of the incomplete Elliptic integral of the second kind, we may approximate $F_{k}^{(0)}$ as

$$F_{k}^{(0)}\sim-4{\rm i}\frac{\sqrt{\lambda A}}{M}\sqrt{c}\int_{0}^{\frac{\rm i}{2}\log(2c)}d\varphi=\frac{2\sqrt{\lambda A}}{M}\sqrt{c}\log(2c).$$

(45)

Noting $c=(k^{2}+m^{2})/(\lambda A)$, we find exponential decay of $|\beta_{k}|$ for large $k$. Thus, the UV divergence does not appear e.g. in $\langle\rho_{\rm np}\rangle$ in (27). The UV finiteness is reasonable since the high momentum modes would not feel oscillation of mass effectively.

We have shown that the narrow and broad resonance regime are described in a unified way. In the broad resonance regime, the particle production rate at one Stokes line crossing is approximated by the model in Sec. 4.1, which is non-perturbative in the coupling constant $\lambda$.¹⁷¹⁷17Around the minimum of the potential $Mt_{0}=\pi/2+2n\pi$, one can expand $\cos Mt\sim-1+M^{2}(t-t_{0})^{2}/2$ and find the correspondence to the model in Sec. 4.1. Such an approximation is used in the analysis of broad resonance in preheating Kofman:1997yn . On the other hand, narrow resonance regime shows a perturbative behavior with respect to $\lambda$, since the growth rate is $\Gamma\propto\lambda$. The interpolation between perturbative and non-perturbative regimes is observed in the context of Schwinger effect Dunne:2005sx ; Basar:2015xna . The weak field limit or multi-photon regime would correspond to the narrow resonance and the Schwinger limit/strong field limit corresponds to the broad resonance. Clarifying the correspondence would be interesting from theoretical viewpoint, and we will leave it for future study.

In this section, we discuss preheating in the expanding Universe. We take in (21) $a(t)=(1+H_{\rm ini}t)^{2/3}$, $h(\phi)=\phi^{2}$, $\phi(t)=A\sin(Mt)/t$, where $H_{\rm ini}$ denotes the initial value of the Hubble parameter at the beginning of preheating stage, and $A$ is a dimensionless constant.¹⁸¹⁸18The coupling constant $\lambda$ in this case is dimensionless. For simplicity, we take $H_{\rm ini}=M$, and note that this setup does not describe the beginning of inflaton oscillation precisely, but is a good approximation in later time. As we will see, the most important dynamics is the oscillating part and the details of the rest would not affect much. We here neglect the terms including Hubble parameter and its time derivative in (24). The effective frequency of $\chi$ with momentum $k$ (24) is given by

$$\omega_{k}^{2}=\frac{k^{2}}{(1+Mt)^{4/3}}+m^{2}+\frac{\lambda A^{2}}{t^{2}}\sin^{2}(Mt).$$

(46)

Due to the expansion of the Universe, the amplitude of inflaton oscillation decays in time, and therefore, particle production is dominated by the early stage. We assume $A\gg 1$, and then the oscillating term dominates the time-dependence of (46).

In order to make an estimate of particle production rate, let us discuss the Stokes phenomenon in this model. When $A\gg 1$, the following approximation used in Hashiba:2021npn can be applied: Since the dominant term in (46) is the sinusoidal term, the turning point would be near $Mt\sim n\pi$, where $n$ is a positive integer. We assume the turning point to be $Mt_{n}=n\pi+\delta_{n}$. Unless $n\gg 1$, we expect $|\delta_{n}|\ll 1$, and therefore,

$$\frac{k^{2}}{(1+n\pi)^{4/3}}+m^{2}+\frac{\lambda A^{2}}{(n\pi)^{2}}\sin^{2}(\delta_{n})\sim 0,$$

(47)

or equivalently,

$$\sin^{2}\delta_{n}\sim-c_{n}^{2},$$

(48)

where

$$c_{n}^{2}\equiv\frac{(n\pi)^{2}}{\lambda A^{2}}\left(\frac{k^{2}}{(1+n\pi)^{4/3}}+m^{2}\right).$$

(49)