text:
# Functional relations in Stokes multipliers - Fun with 𝑥⁶+𝛼𝑥² potential-

## 1 Introduction

The eigenvalue problem of a one-body 1D Schrödinger operator is the most fundamental subject in quantum mechanics. Still, it provides vivid materials of research.

Besides a few exceptions where eigenvalues and wavefunctions are obtainable explicitly, one may employ several tools for analysis, e.g., the perturbation theory, the variational approach and so on. Among them, the exact WKB method - is unique in the sense that it provides non-perturbative information on the analytical structure of wavefunctions and spectral properties. We analytically continue $`x`$, original coordinate variable, to a complex number. The whole complex plane is divided into several sectors. In each sectors, there are two linearly independent solutions , as Schrödinger operator is the 2nd order differential operator. They are referred to as the fundamental system of solutions (FSS) in the sector . The relations among FSS in different sectors are central issue in the connection problem. The importance of the problem and consequently Stokes multipliers, in the WKB problem has been deeply recognized and emphasized in early ’80, especially in .

Recently, a remarkable link has been established among the spectral determinants of a 1D Schrödinger operator associated with the anharmonic oscillator, transfer matrices and $`𝐐`$ operators in CFT possessing $`U_q(\widehat{sl}(2))`$ -. Here the spectral determinants imply $`D(E)=_{E_j\mathrm{eigenvalue}}(1\frac{E}{E_j})`$ and its generalizations. A curious interplay between $`D(E)`$ and generalized Stokes multipliers is also found . In view of solvable models, a striking fact is that they share the same functional relations with transfer matrices in the fusion hierarchy possessing $`U_q(\widehat{sl}(2))`$ . This allows for applications of the strong machineries in the study of solvable models - to the studies of Stokes multipliers, spectral determinants and so on. Several results have been explicitly obtained for the anharmonic oscillator problem, and are extended to higher differential analogues .

In this note, we consider an anharmonic oscillator perturbed by a lower power potential term. It belongs to a class of potentials discussed generally in  with the exact resolution method. To be precise, we consider the eigenvalue problem,

$$(x,\alpha )\mathrm{\Psi }_k(x)=\left(\frac{d^2}{dx^2}+x^{2M}+\alpha x^{M1}\right)\mathrm{\Psi }_k(x)=E_k\mathrm{\Psi }_k(x).$$
(1)
Throughout this report we set $`\mathrm{}=1`$ and $`M>1`$.

The spectral problem concerning this Hamiltonian turns out to be in a category to which one can apply the tools in solvable models.

The sign of $`\alpha `$ seems to be crucial if one considers the operator (1) on the real axis.

$$\overline{)\mathrm{FIG}.1}$$

We will not expect much difference from the ”pure” anharmonic oscillator when $`\alpha >0`$, while we do expect change for $`\alpha <0`$ as the potential develops the double well.

It will be shown, however, that the negative $`\alpha `$ and the positive $`\alpha `$ problems are not separable when we discuss the global connection problem. Roughly speaking, the negative $`\alpha `$ problem is coupled to the positive $`\alpha `$ problem by crossing a border line of neighboring sectors and vice versa. See §2 for precise arguments. It may be then reasonable to consider a two-fold connection problem (crossing two adjacent lines) , or more generally, relations between sectors separated by even multiples of border lines. Some of the Stokes multipliers, in the generalized connection problem, possess expressions corresponding to the eigenvalues of the (fusion) transfer matrices of the 3 state Perk-Schulz (PS) model  of which underlying symmetry is $`U_q(\widehat{gl}(2|1))`$. Others can not be directly equated with the (fusion) transfer matrices but have relations with the the 3 state PS model as well. Thus we conclude that the perturbation $`\alpha x^{M1}`$ breaks the $`U_q(\widehat{sl}(2))`$ symmetry of the ”pure” anharmonic oscillator but it brings the new symmetry $`U_q(\widehat{gl}(2|1))`$. The deformation parameter $`q`$ is related to the exponent of the perturbation by $`q=\mathrm{exp}(i\frac{\pi }{M+1})`$. Through these findings, we can derive the nonlinear integral equations (NLIE) which characterize the energy levels of both the negative $`\alpha `$ problem and the positive $`\alpha `$ problem simultaneously.

The paper is organized as follows. In the next section, we will explore symmetries of solutions to (1). The precise definition of sectors is given. The connection problem is addressed in section 3. Certain components in fusion Stokes matrices are identified with eigenvalues of fusion transfer matrices associated to $`U_q(\widehat{gl}(2|1))`$. The spectral determinant is explicitly parameterized by FSS in a sector. The coupled NLIE are then derived in section4, which determine energy levels. We will perform analytical and numerical checks on the consistency of our result in section 5. Section 6 is devoted to summary and discussions on open problems.

## 2 Asymptotic Expansion and Symmetry of solutions

Let $`\varphi (x,\alpha ,E)`$ be an entire function of $`(x,\alpha ,E)`$ and a solution to $`(x,\alpha )\varphi (x,\alpha ,E)=E\varphi (x,\alpha ,E)`$.

The solution, which decays exponentially at $`x\mathrm{}`$, is of primary interest. By employing the argument in , we immediately find its asymptotic behavior,

$`\varphi (x,\alpha ,E)`$ $``$ $`x^{M/2\alpha /2}\mathrm{exp}({\displaystyle \frac{x^{M+1}}{M+1}}),`$ (2)
$`_x\varphi (x,\alpha ,E)`$ $``$ $`x^{M/2\alpha /2}\mathrm{exp}({\displaystyle \frac{x^{M+1}}{M+1}}).`$ (3)
The validity of the above expansion is not restricted to the real axis, but extends to the wedge, $`|\mathrm{arg}x|<\frac{3\pi }{2M+2}`$ .

The second order linear differential equation admits another independent solution. To specify it, or to deal with the global problem, it is convenient to extend $`x`$ to the complex plane as mentioned in introduction. Then, as in the case of $`\alpha =0`$, the solution exhibits a symmetry by rotating the complex $`x`$ plane by a specific angle.

The direct calculation proves the following.

###### Theorem 1

Let $`\varphi (x,\alpha ,E)`$ be the above solution and $`q=\mathrm{exp}(i\frac{\pi }{M+1})`$. Then $`\varphi (q^1x,q^{M+1}\alpha ,q^2E)`$ is also the solution to the differential equation, $`(x,\alpha )\varphi =E\varphi `$.

This is the desired second solution which grows exponentially on the positive real axis: $`x^{M/2+\alpha /2}\mathrm{exp}(\frac{x^{M+1}}{M+1})`$ for $`x\mathrm{}`$ We note that $`q^{M+1}\alpha =\alpha `$. This deserves an attention. As mentioned in introduction, the potential assumes the completely different structure for $`\alpha `$ positive and $`\alpha `$ negative on the real axis. The rotation in the complex $`x`$ plane by angle $`\frac{\pi }{M+1}`$, however, couples these two problems. Thus we shall treat the Hamiltonians with $`\pm \alpha `$ simultaneously. Similar pairing of differential equations is found for the positive and the negative angular momentum terms in a class of 3rd order differential equations.

This observation is crucial in our approach and can be generalized further. To state it, we prepare some notations.

Hereafter $`\alpha `$ always takes a non-negative real value. By $`^{(ϵ)}(x,\alpha )`$, we mean the Schrödinger operator,

$$\frac{d^2}{dx^2}+x^{2M}+ϵ\alpha x^{M1}$$

where $`ϵ=\pm 1`$.

Let $`𝒮_k`$ be a sector in the plane satisfying

$$|\mathrm{arg}x\frac{k\pi }{M+1}|\frac{\pi }{2M+2}.$$

$$\overline{)\mathrm{FIG}.2}$$

The FSS depends on the sector. We conveniently define

$$y_j^{(ϵ)}:=\frac{q^{j/2ϵ\alpha /2}}{\sqrt{2i}}\varphi (xq^j,ϵ\alpha ,q^{2j}E).$$

###### Theorem 2

For the $`^{(ϵ)}(x,\alpha )`$, the FSS in the sector $`𝒮_j`$ is given by $`(y_j^{(ϵ_j)},y_{j+1}^{(ϵ_{j+1})})`$ where $`ϵ_j=ϵ(1)^j`$.

For $`\alpha =0`$ case, this has been argued in . It is easily checked that $`y_j^{(ϵ_j)}`$ is the sub-dominant solution in $`𝒮_j`$; it tends to zero as $`x`$ tends to infinity along in any direction in the sector.

In the next section, we consider the global connection problem of these FSS in the complex $`x`$ plane.

## 3 Fusion Stokes multipliers, $`U_q(\widehat{gl}(2|1))`$ structure and Spectral Determinants

We introduce the Wronskian matrix

$$\mathrm{\Phi }_j^{(ϵ)}(x):=\left(\begin{array}{cc}y_j^{(ϵ)},& y_{j+1}^{(ϵ)}\\ _xy_j^{(ϵ)},& _xy_{j+1}^{(ϵ)}\end{array}\right).$$
(4)
and the Wronskian $`W_k^{(ϵ)}:=\text{det}\mathrm{\Phi }_k^{(ϵ)}(x)`$.

The linear dependence of the solution can be easily verified by evaluating the Wronskian at $`𝒮_{j+1/2}`$ using the asymptotic expansions (2) and (3) . The present normalization yields $`W_k^{(ϵ)}=1`$.

Let $`_{j,1}^{(ϵ)}`$ be the Stokes matrix connecting the Wronskian matrices $`\mathrm{\Phi }_j^{(ϵ)}(x)`$ and $`\mathrm{\Phi }_{j+1}^{(ϵ)}(x)`$,

$$\mathrm{\Phi }_j^{(ϵ)}(x)=\mathrm{\Phi }_{j+1}^{(ϵ)}(x)_{j,1}^{(ϵ)}.$$
(5)
It permits an explicit parameterization

$$_{j,1}^{(ϵ)}:=\left(\begin{array}{cc}\tau _j^{(ϵ)},& 1\\ 1,& 0\end{array}\right),$$
(6)
where $`\tau _j^{(ϵ)}`$ is referred to as the Stokes multiplier. We have two remarks. First, the (1,1) element is the function of $`\alpha `$ and $`q^{2j}E`$. We omit the dependency on $`\alpha `$. The dependency on $`E`$ is indicated by the index $`j`$. Second, the (2,1) element,$`1`$, is a consequence of the present normalization of the Wronskian.

To be more specific, we consider the operator $`^{(+)}(x,\alpha )`$ and start from the positive real axis (or more generally $`𝒮_0`$). The initial FSS is $`(y_0^{(+)},y_1^{()})`$.

$$\overline{)\mathrm{FIG}.3}$$

A linear relation follows from (5) between this FSS and $`(y_1^{()},y_2^{(+)})`$, the FSS at the neighboring sector $`𝒮_1`$,

$$y_0^{(+)}=\tau _0^{(+)}y_1^{()}y_2^{(+)}.$$
(7)
Similarly, the FSS in $`𝒮_2`$ is linked to the FSS in $`𝒮_1`$ by

$$y_1^{()}=\tau _0^{()}y_2^{(+)}y_3^{()}.$$
(8)
Judging from the upper indices which indicate the corresponding signs of $`\alpha `$, it may be natural to introduce a generalized Stokes matrix $`_{0,2}^{(+)}`$ connecting FSS $`(y_0^{(+)},y_1^{()})`$ and $`(y_2^{(+)},y_3^{()})`$. It is simply obtained by the matrix multiplication,

$$_{0,2}^{(+)}=_{1,1}^{()}_{0,1}^{(+)}=\left(\begin{array}{cc}\tau _1^{()}\tau _0^{(+)}1,& \tau _1^{()}\\ \tau _0^{(+)},& 1\end{array}\right).$$
(9)
Equations (7) and (8) yield $`\tau `$’s in terms of $`y`$’s. The (1,1) component in (9), hereafter denoted by $`T_{1,1}(E)`$, is then represented in terms of $`y`$ as,

$$T_{1,1}(E)=\tau _1^{()}\tau _0^{(+)}1=\frac{y_0^{(+)}}{y_2^{(+)}}+\frac{y_0^{(+)}y_3^{()}}{y_2^{(+)}y_1^{()}}+\frac{y_3^{()}}{y_1^{()}}.$$
(10)
The dependence of $`E`$ in the rhs is implicitly indicated by indices of $`y`$. We will comment on this representation in terms of a solvable model later.

There is another expression using both $`y`$’s and $`y`$’s. This form is of practical use in the following generalization. By applying the Cramer formula to (5), we immediately obtain

$$\tau _j^{(ϵ)}=\left|\begin{array}{cc}y_j^{(ϵ)},& y_{j+2}^{(ϵ)}\\ _xy_j^{(ϵ)},& _xy_{j+2}^{(ϵ)}\end{array}\right|.$$
(11)
Note that we use the fact that the Wronskian is normalized to be unity. The (1,1) entries in (9) is then given by

$$\left|\begin{array}{cc}y_0^{(+)},& y_3^{()}\\ _xy_0^{(+)},& _xy_3^{()}\end{array}\right|.$$

One can further generalize the above result. Naturally, the ”fusion” Stokes matrix $`_{j,2k}^{(+)}`$ is defined which relates FSS of $`𝒮_j`$ to $`𝒮_{j+2k}`$. Explicitly, it is given by

$$_{0,2k}^{(+)}=\left(\begin{array}{cc}\left|\begin{array}{cc}y_0^{(+)},& y_{2k+1}^{()}\\ _xy_0^{(+)},& _xy_{2k+1}^{()}\end{array}\right|,& \left|\begin{array}{cc}y_1^{()},& y_{2k+1}^{()}\\ _xy_1^{()},& _xy_{2k+1}^{()}\end{array}\right|\\ \left|\begin{array}{cc}y_0^{(+)},& y_{2k}^{(+)}\\ _xy_0^{(+)},& _xy_{2k}^{(+)}\end{array}\right|,& \left|\begin{array}{cc}y_1^{()},& y_{2k}^{(+)}\\ _xy_1^{()},& _xy_{2k}^{(+)}\end{array}\right|\end{array}\right),$$
(12)
for $`j=0`$. We can prove the above using the induction on $`k`$ most easily. Similar formula holds for $`_{0,2k}^{()}`$ by replacing all upper indices $`+`$.

We are now ready to relate an entry in a fusion Stokes matrix to the spectral determinant. Hereafter we assume $`M=2m1`$. It follows from the above argument that

$$\mathrm{\Phi }_{2m}^{(+)}=\mathrm{\Phi }_0^{(+)}(_{0,2m}^{(+)})^1.$$
(13)
$$\overline{)\mathrm{FIG}.4}$$

$`y_0^{(+)}`$( $`y_{2m}^{(+)}`$) stands for the subdominant solution on the positive (negative) real axis. They tend to zero asymptotically in their proper region, being appropriate basis for the eigenfunction. The eq.(13) tells, however, that $`y_{2m}^{(+)}`$ is combined to both $`y_0^{(+)}`$ and $`y_1^{()}`$ by rotating the complex plane by $`\pi `$,

$`y_{2m}^{(+)}`$ $`=`$ $`(c_1y_0^{(+)}+c_2y_1^{()})/\mathrm{det}_{0,2m}^{(+)}`$
$`c_1`$ $`=`$ $`(_{0,2m}^{(+)})_{1,1},c_2=(_{0,2m}^{(+)})_{2,1}.`$
This is an obstacle in constructing an eigenfunction defined on the whole real axis. The prescription is to demand that the coefficient of $`y_1^{()}`$ must vanish if $`E`$ is an eigenvalue. Consequently, it is proportional to the spectral determinant.

The coefficient is essentially equal to the (2,1) component of $`_{0,2m}^{(+)}`$, and it reads in terms of the original $`\varphi `$ function as

$$\frac{q^{\alpha (m+1)}}{2}\left|\begin{array}{cc}\varphi (x,\alpha ,E),& \varphi (x,\alpha ,E)\\ _x\varphi (x,\alpha ,E),& _x\varphi (x,\alpha ,E)\end{array}\right|.$$

The $`x`$ dependencies are spurious as the entities are products of Stokes multipliers which are obviously $`x`$ independent. We adopt the simplest choice $`x=0`$. The coefficient is now proportional to $`\varphi (0,\alpha ,E)_x\varphi (x,\alpha ,E)|_{x=0}`$. Thus we conclude that for an eigenvalue $`E_j^{(+)}`$ of $`^{(+)}(x,\alpha )`$,

$$\varphi (0,\alpha ,E_j^{(+)})=0\text{ or }_x\varphi (x,\alpha ,E_j^{(+)})|_{x=0}=0$$
(14)
must hold.

We can repeat the same argument starting from $`^{()}(x,\alpha )`$ on the positive real axis.

The above observation may lead to the identification

$`\varphi (0,ϵ\alpha ,E)`$ $``$ $`D_{}^{(ϵ)}(E):={\displaystyle \underset{j}{}}(1E_{,j}^{(ϵ)}),`$ (15)
$`_x\varphi (x,ϵ\alpha ,E)|_{x=0}`$ $``$ $`D_+^{(ϵ)}(E):={\displaystyle \underset{j}{}}(1E_{+,j}^{(ϵ)}).`$ (16)
The lower sign signifies the parity: the positive parity means a contribution from symmetric wave function. The product must be taken over eigenvalues with the corresponding parity. The total set eigenvalues $`\{E_j^{(ϵ)}\}`$ of $`^{(ϵ)}(x,\alpha )`$ consists of two subsets, $`\{E_j^{(ϵ)}\}`$ =$`\{E_{+,j}^{(ϵ)}\}\{E_{,j}^{(ϵ)}\}`$ and $`D^{(ϵ)}(E)=D_+^{(ϵ)}(E)D_{}^{(ϵ)}(E)`$.

We comment on the relation of the present result to an existing solved model. $`T_{1,1}(E)`$ in (10) can be represented, utilizing (16), as

$$T_{1,1}(E)=q^{\alpha 1}\frac{D_{}^{(+)}(E)}{D_{}^{(+)}(q^4E)}+q^{2\alpha }\frac{D_{}^{(+)}(E)D_{}^{()}(q^6E)}{D_{}^{(+)}(q^4E)D_{}^{()}(q^2E)}+q^{\alpha +1}\frac{D_{}^{()}(q^6E)}{D_{}^{()}(q^2E)},$$
(17)
where we safely choose $`x=0`$ in the rhs.

The above expression has similarity to the dressed vacuum form (DVF) of the (unfused) transfer matrix for the 3-state PS model with grading $`(+,,+)`$. The latter can be found in eq.(3.1) and (3.2) of . The spectral parameter $`v`$ corresponds to energy in the Schrödinger operator, precisely, $`E=\mathrm{exp}(\frac{2\pi v}{M+1})`$.

The spectral determinants have the following identification to the eigenvalues of Baxter’s $`Q`$ operators,

$`D_{}^{(+)}(E)`$ $`=`$ $`Q_2(v+{\displaystyle \frac{j2}{2}}i)`$
$`D_{}^{()}(E)`$ $`=`$ $`Q_1(v+{\displaystyle \frac{j1}{2}}i).`$
Those with positive parity may be identified with the second solutions of Baxter’s $`Q`$ operators.

The scalar factors (vacuum expectation values) $`f_a(x),g_a(x)`$ in  depend on the choice of the quantum space. We assume that the present quantum space space gives the simple scalars as in (17). In this sense, $`T_{1,1}(E)`$ exhibits the hidden $`U_q(\widehat{gl}(2|1))`$ symmetry behind the present Schrödinger operator, just as in the $`U_q(\widehat{sl}(2))`$ symmetry for $`\alpha =0`$ problem. This coincidence can be observed further. We have checked up to certain value of $`k`$ that the (1,1) element and the (2,2) element of $`_{0,2k}^{(+)}`$ coincide with DVF of symmetric fusion transfer matrices $`\mathrm{\Lambda }_k^{(1)}`$ and $`\mathrm{\Lambda }_{k1}^{(1)}`$ in , respectively. The interpretation of the (1,2) and the (2,1) element, in terms of fusion transfer matrices, is still an open problem.

One can adopt another description of $`T_{1,1}`$. The (2,1) component of (5) results

$$\tau _0^{(ϵ)}=\frac{y_0^{(ϵ)}}{y_1^{(ϵ)}}+\frac{y_2^{(ϵ)}}{y_1^{(ϵ)}}.$$
(18)
Proceeding as above, we arrive at,

$$T_{1,1}(E)=q^{\alpha +1}\frac{D_+^{(+)}(E)}{D_+^{(+)}(q^4E)}+q^{2\alpha }\frac{D_+^{(+)}(E)D_+^{()}(q^6E)}{D_+^{(+)}(q^4E)D_+^{()}(q^2E)}+q^{\alpha 1}\frac{D_+^{()}(q^6E)}{D_+^{()}(q^2E)}.$$
(19)
In the next section, we determine the energy levels by utilizing the above results.

## 4 Nonlinear Integral equations for eigenvalue problem

The Bethe ansatz equations follow from the pole-free property of $`T_{1,1}(E)`$ on the real $`E`$ axis,

$`q^{\alpha 1}{\displaystyle \frac{D_+^{()}(q^2E_{+,j}^{(+)})}{D_+^{()}(q^2E_{+,j}^{(+)})}}`$ $`=`$ $`q^{\alpha 1}{\displaystyle \frac{D_+^{(+)}(q^2E_{+,j}^{()})}{D_+^{(+)}(q^2E_{+,j}^{()})}}=1,`$ (20)
$`q^{\alpha +1}{\displaystyle \frac{D_{}^{()}(q^2E_{,j}^{(+)})}{D_{}^{()}(q^2E_{,j}^{(+)})}}`$ $`=`$ $`q^{\alpha +1}{\displaystyle \frac{D_{}^{(+)}(q^2E_{,j}^{()})}{D_{}^{(+)}(q^2E_{,j}^{()})}}=1.`$ (21)
To their analysis, we apply the strong machinery in solvable models, the method of nonlinear integral equations. Hereafter we shall confine ourselves to the case $`0\alpha M`$ where energies are nonnegative. The simple pattern of energy spectrum permits the following simple-mind choice of auxiliary functions,

$`a_ϵ^{}^{(ϵ)}(E)`$ $`:=`$ $`q^{ϵ\alpha ϵ^{}}{\displaystyle \frac{D_ϵ^{}^{(ϵ)}(q^2E)}{D_ϵ^{}^{(ϵ)}(q^2E)}}`$ (22)
$`A_ϵ^{}^{(ϵ)}(E)`$ $`:=`$ $`1+a_ϵ^{}^{(ϵ)}(E).`$ (23)
Thus

$$A_ϵ^{}^{(ϵ)}(E_{ϵ^{},j}^{(ϵ)})=0.$$
(24)
Remember that $`ϵ(=\pm 1)`$ represents the signature of the perturbation while $`ϵ^{}(=\pm 1)`$ denotes the parity.

In addition, we need some inputs about the asymptotic behaviors from the WKB method. Fortunately, they are already available  as the existence of lower power term does not alter them.

$`\mathrm{ln}D_\pm ^{(ϵ)}(E)`$ $``$ $`{\displaystyle \frac{a_0}{2}}(E)^\mu ,|E|\mathrm{},|\mathrm{arg}(E)|<\pi `$ (25)
$`b_0E_j^{(ϵ)}`$ $``$ $`2\pi (j+{\displaystyle \frac{1}{2}}),j\mathrm{}`$ (26)
$`\mu `$ $`=`$ $`{\displaystyle \frac{M+1}{2M}},a_0={\displaystyle \frac{b_0}{2\mathrm{sin}\mu \pi }},`$
$`b_0`$ $`=`$ $`{\displaystyle \frac{\pi ^{1/2}\mathrm{\Gamma }(\frac{1}{2M})}{M\mathrm{\Gamma }(\frac{3}{2}+\frac{1}{2M})}}.`$
There might be several routes to reach nonlinear integral equations among $`a_{(ϵ^{})}^{(ϵ)}`$ and $`A_{(ϵ^{})}^{(ϵ)}`$. Here we choose the quickest way  which fully exploits the fact that zeros of $`D_\pm ^{(ϵ)}(E)`$ are on the positive real $`\theta `$ axis. In addition, we assume that there are no zeros of $`A_\pm ^{(ϵ)}(E)`$ inside the narrow strip including the positive real axis other than those from zeros of $`D_\pm ^{(ϵ)}(E)`$.

Apparently we have

$`\mathrm{log}a_ϵ^{}^{(ϵ)}(E)`$ $`=`$ $`{\displaystyle \frac{(ϵ\alpha ϵ^{})\pi }{M+1}}i+{\displaystyle \underset{j}{}}F({\displaystyle \frac{E}{E_{ϵ^{},j}^{(ϵ)}}})`$
$`F(E)`$ $`=`$ $`\mathrm{log}{\displaystyle \frac{1q^2E}{1q^2E}}.`$
The above assumption allows the representation of the summation part by an integral over contour $`𝒞_E`$ which surrounds the positive real axis counterclockwise,

$$\mathrm{log}a_ϵ^{}^{(ϵ)}(E)=\frac{(ϵ\alpha ϵ^{})\pi }{M+1}i+\frac{1}{2\pi i}_{𝒞_E}𝑑E^{}F(\frac{E}{E^{}})_E^{}\mathrm{log}A_ϵ^{}^{(ϵ)}(E^{}).$$
(27)
For convenience, we introduce a variable $`\theta `$  by

$$E=\mathrm{exp}(\theta /\mu )/\nu ^2\nu =(2M+2)^{\frac{1}{2\mu }}/\mathrm{\Gamma }(\frac{1}{2\mu })$$

which originates from the matching condition of the WKB result (25) and the $`Q`$ operator analysis .

Let $`𝔞_ϵ^{}^{(ϵ)}(\theta ),𝔄_ϵ^{}^{(ϵ)}(\theta )`$ be auxiliary functions defined in (22), (23) regarded as functions of $`\theta `$. Then eq.(27) reads,

$`\mathrm{log}𝔞_ϵ^{}^{(ϵ)}(\theta )`$ $`=`$ $`{\displaystyle \frac{(ϵ\alpha ϵ^{})\pi }{M+1}}i+{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{𝒞_\theta }}𝑑\theta ^{}G(\theta \theta ^{})_\theta ^{}\mathrm{log}𝔄_ϵ^{}^{(ϵ)}(\theta ^{})`$
$`G(\theta )`$ $`=`$ $`\mathrm{log}\left(q^2{\displaystyle \frac{\mathrm{sinh}(\frac{M\theta }{M+1}+i\frac{\pi }{M+1})}{\mathrm{sinh}(\frac{M\theta }{M+1}i\frac{\pi }{M+1})}}\right)`$
where $`𝒞_\theta `$ encircles the whole real axis counterclockwise.

For the reason which will be supplemented, we shall keep $`𝔞_ϵ^{}^{(ϵ)}(\theta )`$ in the lower half plane but use $`1/𝔞_ϵ^{}^{(ϵ)}(\theta )`$ in the upper half plane.

This requirement modifies the above expression as

$`\mathrm{log}𝔞_ϵ^{}^{(ϵ)}(\theta )`$ $`=`$ $`{\displaystyle \frac{(ϵ\alpha ϵ^{})\pi }{M+1}}i{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}_\theta G(\theta \theta ^{}+i0)\mathrm{log}𝔞_ϵ^{}^{(ϵ)}(\theta ^{}i0)`$ (28)
$`+{\displaystyle \frac{1}{\pi }}\mathrm{}\left({\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\theta ^{}_\theta G(\theta \theta ^{}+i0)\mathrm{log}𝔄_ϵ^{}^{(ϵ)}(\theta ^{}i0)\right),`$
where $`\theta `$ is assumed to possess small negative imaginary part. The property $`(a_ϵ^{}^{(ϵ)}(\theta ))^{}=\frac{1}{a_ϵ^{}^{(ϵ)}(\theta ^{})}`$ is employed in the above transformation.

We solve (28) in terms of $`\mathrm{log}𝔞_ϵ^{}^{(ϵ)}(\theta )`$ to reach the final expression of NLIE,

$`\mathrm{ln}𝔞_ϵ^{}^{(ϵ)}(\theta )`$ $`=`$ $`{\displaystyle \frac{i}{2}}b_0\nu ^{2\mu }\mathrm{e}^\theta +{\displaystyle \frac{\pi }{2}}i(ϵ^{}+ϵ{\displaystyle \frac{\alpha }{M}})`$ (29)
$`+`$ $`2i\mathrm{}\{{\displaystyle _{\mathrm{}}^{\mathrm{}}}K_1(\theta \theta ^{}+i0)\mathrm{ln}𝔄_ϵ^{}^{(ϵ)}(\theta ^{}i0)d\theta ^{}`$
$`+{\displaystyle _{\mathrm{}}^{\mathrm{}}}K_2(\theta \theta ^{}+i0)\mathrm{ln}𝔄_ϵ^{}^{(ϵ)}(\theta ^{}i0)d\theta ^{}\}.`$
The kernel functions read

$`K_1(\theta )`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{iw\theta }{\displaystyle \frac{\mathrm{sinh}^2\frac{\pi (M1)w}{2M}}{\mathrm{sinh}\pi w\mathrm{sinh}\frac{\pi w}{M}}}`$
$`K_2(\theta )`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{iw\theta }{\displaystyle \frac{\mathrm{sinh}\frac{\pi (M+1)w}{2M}\mathrm{sinh}\frac{\pi (M1)w}{2M}}{\mathrm{sinh}\pi w\mathrm{sinh}\frac{\pi w}{M}}}.`$
Few remarks are in order.

1. As a consequence of connection rules, the integral equations are coupled for auxiliary functions related to the positive and the negative coefficient of $`x^{M1}`$ .
2. On the other hand, equations with different parities are decoupled.
3. The constants are determined from the consistency by putting $`\theta \mathrm{}`$. The $`\alpha `$ dependence is only summarized in these constants.
4. The first term in the rhs is determined so that we recover the result from the WKB method (25) by dropping contributions of integrals. Clearly, $`𝔞_ϵ^{}^{(ϵ)}(\theta )`$ is bounded in the upper-half plane. This explains our choice of appropriate half planes for auxiliary functions.

The eigenvalues $`\{E_{j,\pm }^{(ϵ)}\}`$ are evaluated by

$$\mathrm{ln}𝔞_\pm ^{(ϵ)}(\theta _{j,\pm }^{(ϵ)})=(2j+1)\pi i,\text{ and }E_{j,\pm }^{(ϵ)}=\mathrm{exp}(\theta _{j,\pm }^{(ϵ)}/\mu )/\nu ^2.$$

More explicitly

$`{\displaystyle \frac{1}{2}}b_0\nu ^{2\mu }\mathrm{e}^{\theta _{j,ϵ^{}}^{(ϵ)}}=(2j+1ϵ^{}{\displaystyle \frac{1}{2}}+ϵ{\displaystyle \frac{\alpha }{2M}})\pi `$
$`+2\mathrm{}\{{\displaystyle _{\mathrm{}}^{\mathrm{}}}K_1(\theta _{j,ϵ^{}}^{(ϵ)}\theta ^{}+i0)\mathrm{ln}𝔄_ϵ^{}^{(ϵ)}(\theta ^{})𝑑\theta ^{}+{\displaystyle _{\mathrm{}}^{\mathrm{}}}K_2(\theta _{j,ϵ^{}}^{(ϵ)}\theta ^{}+i0)\mathrm{ln}𝔄_ϵ^{}^{(ϵ)}(\theta ^{})𝑑\theta ^{}\}`$
$`(j0).`$ (30)
We present examples of numerical solutions to (28) in Fig. 5. The real and the imaginary parts of $`\mathrm{ln}𝔄_+^{(\pm )}`$ are depicted for $`M=3,\alpha =1`$.

$$\overline{)\mathrm{FIG}.5}$$

## 5 benchmarks

We shall check the nonlinear integral equations analytically for limiting cases and numerically.

(1) $`\alpha =0`$ case

By putting, $`𝔞_\pm ^{(+)}(E)=𝔞_\pm ^{()}(E)`$ the coupled NLIE reduce to an identical integral equation. Immediately seen, the result coincides with the nonlinear integral equation in .

(2) $`\alpha =M`$ case

In this case, we have a duality in energy spectra; $`\{E_j^{(+)}\}`$ coincide with $`\{E_j^{()}\}`$ , except for $`E_0^{()}=0`$ in the latter. This degeneracy can be easily explained by the following representation of the Hamiltonians ,

$`^{()}(x,\alpha =M)`$ $`=`$ $`𝒟^{}𝒟`$
$`^{(+)}(x,\alpha =M)`$ $`=`$ $`𝒟𝒟^{}`$
$`𝒟`$ $`=`$ $`{\displaystyle \frac{1}{i}}{\displaystyle \frac{d}{dx}}ix^M.`$
Once an eigenvector $`^{()}(x,M)\psi _j^{()}=E_j^{()}\psi _j^{()}`$ is found, we can construct the eigenvector for $`^{(+)}(x,\alpha )`$ with the same energy by $`\psi _{j1}^{(+)}:=𝒟\psi _j^{()}`$. Only the exception is the $`j=0`$ case where $`𝒟\psi _{j=0}^{()}=0`$. It is interesting that the asymptotic form (2) from the WKB type argument is exact for all $`x`$ in this case. <sup>1</sup><sup>1</sup>1I thank V.V. Bazhanov, R.J. Baxter and B. Nienhuis for pointing out the explicit eigenfunction for $`\psi _0^{()}`$ right after my talk.

The above facts can be also verified from (29). Note that the rhs can be treated as the mod $`2\pi i`$ quantity. Then the choice $`\alpha =M`$ leads to the same coupled equations under identifications $`𝔞_+^{(+)}(\theta )𝔞_{}^{()}(\theta )`$, $`𝔞_{}^{(+)}(\theta )𝔞_+^{()}(\theta )`$. This explains the degeneracy of the spectra as it consists both from the negative and the positive parity contributions. The zero energy case must be treated more separately. By choosing $`j=0,ϵ=ϵ^{}=1,\alpha =M`$, we find the first term in lhs of (30) is null. So the first order approximation is $`\theta _{+,0}^{()}=\mathrm{}`$. Actually this is exact as we determine the constant terms so that NLIE are consistent in $`\theta \mathrm{}`$. See remark 3 after (29). This solution gives the missing energy 0.

Finally we present the preliminary numerical results for $`M=3`$.

Table 1 shows the results from the IMSL package (dsleig.f), the (naive) WKB method and those obtained by solving the nonlinear equations. The agreement is not yet precise enough (typically 3-4 digits). Some implement is still in need for the numerical accuracy. Nevertheless, the NLIE data already show much improvement from the (naive) WKB results.

By the (naive) WKB method, we mean a self-consistent determination of $`E_j^{(ϵ)}`$ by

$$|p|𝑑x=_{x_0}^{x_0}\sqrt{E_j^{(ϵ)}x^6ϵ\alpha x^2}𝑑x=(j+\frac{1}{2})\pi $$
(31)
where $`E_j^{(ϵ)}x_0^6ϵ\alpha x_0^2=0`$. Particularly, for the value with asterisk, this method has subtlety. Immediately seen, $`E_0^{(ϵ)}=0,x_0=2^{1/4}`$ is a formal solution to (31) for $`j=0,ϵ=,\alpha =2`$. It however involves an isolated turning point of the 2nd order at the origin if $`E_0^{(ϵ)}=0`$, which spoils the simple application of the condition (31). The value with $``$ has similar difficultly. We however skip further discussion on the validity on the (naive) WKB method as it is out of the present subject.

Summarizing, we check the consistency of (29) in some limiting cases and by numerical methods.

## 6 Summary and Discussion

In this report, the eigenvalue problem has been addressed for the 1D quantum systems of which Hamiltonians include double well potentials. We have successfully derived the coupled NLIE which determine energy levels of the systems with potential terms of $`\pm \alpha x^{M1}+x^{2M}`$ at the same time.

The essence of our strategy is to utilize the following correspondences between 1D quantum mechanics and 1+1 D solvable models,

| energy | $``$ | spectral parameter |
| --- | --- | --- |
| Stokes multipliers | $``$ | transfer matrices |
| eigenfunctions or derivatives at $`x=0`$ | $``$ | vacuum expectation values of $`Q`$ operators |

We are then entitled to apply the strong machinery of the latter developed since Baxter’s revolution.

There are several open questions.

1. In this report we confine ourselves to the simplest case $`\alpha M`$. For $`\alpha >M`$, the existence of negative eigenvalues ruins the analyticity assumptions on auxiliary functions. Still, formal expressions of NLIE are possible which are similar to excited states TBA equations. The integration contour is, however, not so simple as described here. This is an apparent drawback in actual numerical investigations. The clever choice of auxiliary functions may be desired.
2. The understanding is lacking on the intrinsic reason why affine symmetry like $`U_q(\widehat{sl}(2))`$ or $`U_q(\widehat{gl}(2|1))`$ comes into play in this simple 1D quantum mechanical model.
3. This is somewhat related to the above, but is the most intriguing question. Where is the Yang-Baxter equation in the 1D Schrödinger operator problem?. Once this is known, the fusion hierarchy, useful in the present study, is a mere corollary of it.

We hope to answer these in the future publication.

## Acknowledgments

The author thanks V. V. Bazhanov, A. Kuniba, J.M. Maillard, B. Nienhuis P. Pearce, and C. Richard for comments, discussions and encouragement. He also would like to thank the organizers of ”Baxter’s revolution in mathematical physics”.

FIGURE CAPTIONS

Fig 1. An anharmonic oscillator perturbed by a positive or a negative perturbation term.

Fig 2. The complex plane is divided into sectors. $`𝒮_0`$ and $`𝒮_1`$ are indicated as examples.

Fig 3. The FSS in $`𝒮_0`$ and $`𝒮_1`$ are related by the matrix $`_{0,1}^{(+)}`$.

Fig 4. The connection of FSS on the negative and the positive real axis is accomplished by $`_{0,2m}^{(+)}`$.

Fig 5. Left: the real part of $`\mathrm{ln}𝔄_+^{(\pm )}`$, Right: the imaginary part of $`\mathrm{ln}𝔄_+^{(\pm )}`$.

TABLE CAPTIONS

Table 1. First two energy levels calculated by the IMSL library, the (naive) WKB method and the NLIE method. We choose $`M=3`$ and adopt various $`\alpha `$. See asterisk and $``$ for the text.


format: text

text:
# On the Theory of Fermionic Preheating

## I Introduction

During cosmic inflation, it is assumed that entropy and temperature associated with particles of matter are diluted to practically zero values together with the number density of particles. After inflation, the inflaton field $`\varphi `$ oscillates around the minimum of its effective potential $`V(\varphi )`$. The energy of inflaton oscillations is converted to the energy of newly created particles of matter. Eventually particles of different species are settled in a state of thermal equilibrium which marks the beginning of the conventional epoch of the hot Friedmann radiation domination. The actual process of particle creation from the classical background inflaton oscillations occurs very rapidly in the regime of parametric resonance  before the thermal equilibrium will be settled. Particles are created non-perturbatively in the out-of-equilibrium state. The theory of this process, preheating, is elaborated in details for the creation of bosons . For bosons (denoted $`\chi `$), the leading effect is the stimulated process of particle creation in the regime of parametric resonance, where the number density of created particles copiously increases with time as $`n_\chi \mathrm{exp}\mu 𝑑t`$. Soon, the backreaction of created $`\chi `$ particles becomes important, so that the self-consistent dynamics of interacting bose fields $`\varphi (t,\stackrel{}{x})`$ and $`\chi (t,\stackrel{}{x})`$, which can be treated classically, can be revealed with lattice simulations .

At the beginning of the preheating investigation, a study of fermion production did not look very interesting relative to the bosonic case. Indeed, the number density of fermions $`\psi `$ is bounded by Pauli blocking. Therefore, it was not expected that fermions will influence the dynamics of the inflaton and other scalar fields (despite a numerical study  which down-played an observation that fermion production from an oscillating scalar is different from the perturbative prediction).

However, it was understood with surprise  that creation of fermions from the coherently oscillating inflaton field occurs very differently than what the conventional perturbation theory of inflaton decay $`\varphi \overline{\psi }\psi `$ (say, due to a Yukawa coupling $`h\overline{\psi }\varphi \psi `$) would predict, and this may have many interesting cosmological applications. It turns out that the occupation number of fermions very quickly, within about ten inflaton oscillations, is saturated at the time-average value of about $`n_\psi 1/2`$. Moreover, in momentum space, fermions are excited within a non-degenerate ”Fermi sphere” of a large radius $`kq^{1/4}m`$, where $`m`$ is the frequency of inflaton oscillations, $`q`$ is the usual dimensionless parameter of parametric excitation, $`q=\frac{h^2\varphi _0^2}{m^2}`$. Ironically, it may be parametric excitations of fermions that are responsible for the most important observable signatures or observational constraints of preheating. Indeed, in some inflationary models there is significant production of gravitinos during the preheating stage . Gravitinos are cosmologically dangerous relics, for a mass of $`100`$ Gev their abundance relative to that of the relic photons cannot exceed the bound $`n_{3/2}/n_\gamma 10^{15}`$. Theory of gravitino production from preheating is rooted in the theory of spin 1/2 fermionic preheating. Another potentially important application of fermionic preheating is a possibility to produce superheavy fermions with a mass as large as $`10^{18}`$ Gev from inflatons of mass $`10^{13}`$ Gev, as was noticed in  and investigated in . Superheavy fermions may be interesting for the dark matter problem and for the problem of ultra-high energy cosmic rays. There are other interesting cosmological applications for the creation of fermions, e.g. the scenario of instant preheating  and the creation of massive fermions during inflation . Recently, fermionic preheating in hybrid inflation for some range of parameters was thoroughly studied . At present, it is hard to say how important fermionic preheating will be in the self-consistent non-linear dynamics of bose- and fermi-fields during preheating. Technical difficulty is to incorporate fermions into lattice simulations alongside with bosons. However, in  it was reported that even within the Hartree approximation backreaction of fermions catalyzes bosonic preheating. After all, for complete decay of inflatons one needs a ‘three-legs’ interaction, which is provided naturally by the Yukawa coupling with fermions.

In this paper we develop the theory of fermionic preheating in an expanding universe, following our short paper , with an emphasis on the analytic results. In particular, we will generalize for the fermionic case some of the methods which we earlier developed for bosonic preheating , in particular, the method of parabolic scattering, which works for large values of $`q`$, gives us an analytic formula for the occupation number of fermions $`n_k(t)`$ as function of time and momentum. We extend this method for production of superheavy fermions from a moving scalar field. We begin Sect. 2 of this paper with the Dirac equation and different VEVs for fermions interacting with a time-dependent background scalar field. In Sect. 3 the creation of fermions without expansion of the universe will be considered. We mostly will consider scalar fields oscillating in a potential $`V=\frac{m^2\varphi ^2}{2}`$, although the methods can easily be applied to others. We will give a semi-analytic treatment of the problem based on some earlier results. For the case $`q1`$ we develop the method of successive parabolic scatterings. In Sect. 4 we take into account expansion of the universe, when preheating of fermions acquires new qualitative features. We extend the theory to describe the production of superheavy fermions. Our formalism also includes the case when fermions are created not from inflaton oscillations, but from a single instance, when inflaton field crosses a certain level.

## II  Formalism: Fermions Coupling with Background Scalar in FRW metrics.

We will consider the creation of spin-$`\frac{1}{2}`$ Dirac fermions $`\psi `$ by a homogeneous, oscillating scalar field $`\varphi `$ in an expanding, flat FRW universe. Our strategy will be to solve the Heisenberg equation of motion for the quantum $`\psi `$-field in the presence of the classical backgrounds $`\varphi `$ and $`g_{\mu \nu }`$. Furthermore, we will assume the energy-momentum of the $`\varphi `$-field alone determines the expansion rate of the Universe. The matter action we use

$$S_\mathrm{M}[\varphi ,\psi ,e_\mu ^\alpha ]=d^4xe\left[\frac{1}{2}_\mu \varphi ^\mu \varphi V(\varphi )+i\overline{\psi }\overline{\gamma }^\mu \underset{\mu }{\overset{}{D}}\psi (m_\psi +h\varphi )\overline{\psi }\psi \right],$$
(1)
contains a simple Yukawa coupling between the scalar field with effective potential $`V(\varphi )`$ and fermion with bare-mass $`m_\psi `$. Here $`\overline{\gamma }^\mu `$ is a space-time dependent Dirac gamma matrix, $`e_\mu ^\alpha `$ is the vierbein with $`e`$ its determinant, and $`D_\mu =_\mu +\frac{1}{4}\gamma _{\alpha \beta }\omega _\mu ^{\alpha \beta }`$ is the spin-$`\frac{1}{2}`$ covariant derrivative with vierbein dependent spin-connection, $`\omega _\mu ^{\alpha \beta }`$. The unbarred gammas are standard, Minkowski space-time Dirac matrices and $`\gamma _{\alpha \beta }\gamma _{[\alpha }\gamma _{\beta ]}`$.

### A  Classical Background Fields

We consider flat FRW metrics $`ds^2=dt^2a(t)^2d\stackrel{}{x}^2`$ where $`a(t)`$ is the scale factor of the universe. The Hubble parameter $`H\frac{\dot{a}}{a}`$ is determined by the equation $`H^2=\frac{8\pi }{3M_{\mathrm{pl}}^2}\left(\frac{1}{2}(\dot{\varphi })^2+V(\varphi )\right)`$. The background homogeneous scalar field obeys the equation of motion $`\ddot{\varphi }+3H\dot{\varphi }+\frac{V}{\varphi }=0`$. We will be most interested in fermion creation while the $`\varphi `$-field oscillates about a minimum of its potential. For illustration of the methods, we will consider fermion creation in the context of chaotic inflation scenarios with potentials of the form $`V(\varphi )=\frac{1}{2}m_\varphi ^2\varphi ^2+\frac{\lambda }{4}\varphi ^4`$. After inflation, the $`\varphi `$-field oscillates quasi-periodically with a slowly decreasing amplitude due to the Hubble expansion. The specific value of the parameter $`m_\varphi `$ or $`\lambda `$ (or some combination thereof) will be dictated by cosmology. We mostly will consider quadratic potential with $`\lambda =0`$. In this paper we ignore backreaction of created particles.

### B  Quantum Dirac Field

Variation of the action (1) with respect to $`\overline{\psi }`$ leads to the general relativistic generalization of the Dirac equation

$$\left[i\overline{\gamma }^\mu D_\mu (m_\psi +h\varphi (t))\right]\psi (x)=0.$$
(2)
It can be shown that for an FRW space-time we have $`\overline{\gamma }^0=\gamma ^0`$ and $`\overline{\gamma }^i=\frac{1}{a(t)}\gamma ^i`$, where the $`\gamma ^\alpha `$’s are standard, Minkowski space-time Dirac matrices. Furthermore, spin connection leads to

$$\frac{1}{4}\overline{\gamma }^\mu \gamma _{\alpha \beta }\omega _\mu ^{\alpha \beta }=\frac{3}{2}\left(\frac{\dot{a}}{a}\right)\gamma ^0.$$
(3)
Thus, in our case, the Dirac equation becomes

$$\left[i\gamma ^0_0+i\frac{1}{a}\stackrel{}{\gamma }\stackrel{}{}+i\frac{3}{2}\left(\frac{\dot{a}}{a}\right)\gamma ^0(m_\psi +h\varphi )\right]\psi (x)=0.$$
(4)
Only standard gamma matrices now appear. Similarly, $`\overline{\psi }=\psi ^{}\gamma ^0`$ is found to obey the conjugate of Equation (4).

In the Heisenberg representation of QFT, the Dirac equation (2) or (4) becomes the equation of motion for the $`\psi (x)`$ field operator. $`\psi (x)`$ may be decomposed into eigen-spinors

$$\psi (𝐱,t)=\underset{s=\pm }{}\frac{d^3k}{(2\pi )^3}\left(\widehat{a}_{k,s}𝐮_{k,s}(t)e^{+i𝐤𝐱}+\widehat{b}_{k,s}^{}𝐯_{k,s}(t)e^{i𝐤𝐱}\right),$$
(5)
where $`𝐮_{k,\pm }(t)`$ is a positive-frequency eigenspinor of the Dirac equation (4) with helicity $`\pm \frac{1}{2}`$ and $`𝐯_{k,\pm }(t)=𝒞\overline{𝐮}_{k,\pm }^T(t)`$ is its charge conjugate. Here $`𝒞=i\gamma _0\gamma _2`$ is the standard charge conjugation matrix. To construct these eigen-spinors we use the ansatz

$$𝐮_{k,\pm }(t)e^{+i𝐤𝐱}=a\left[i\gamma ^0_0i\frac{1}{a}\stackrel{}{\gamma }\stackrel{}{}i\frac{3}{2}\left(\frac{\dot{a}}{a}\right)\gamma ^0(m_\psi +h\varphi )\right]\chi _k(t)R_\pm (𝐤)e^{+i𝐤𝐱},$$
(6)
where $`R_\pm (𝐤)`$ are eigenvectors of the helicity operator $`𝐤𝚺`$ such that $`𝐤𝚺R_\pm (𝐤)=\pm 1`$ and $`\gamma ^0R_\pm (𝐤)=+1`$. The time dependence of the eigen-spinor is contained intirely in the mode function $`\chi _k(t)`$. Substituting this ansatz into equation (4) leads to the mode equation

$$\ddot{\chi }_k+4\frac{\dot{a}}{a}\dot{\chi }_k+\left[\frac{k^2}{a^2}+m_{\mathrm{eff}}^2i\frac{\dot{(aM_{\mathrm{eff}})}}{a}+\frac{9}{4}\left(\frac{\dot{a}}{a}\right)^2+\frac{3}{2}\frac{\ddot{a}}{a}\right]\chi _k=0,$$
(7)
where $`k^2=|𝐤|^2`$ and the effective mass of the fermions is

$$M_{\mathrm{eff}}(m_\psi +h\varphi ).$$
(8)
The damping term in this equation may be removed by defining a new mode function $`X_k(t)=a^2\chi _k(t)`$. The mode equation becomes

$$\ddot{X}_k+\left[\frac{k^2}{a^2}+M_{\mathrm{eff}}^2i\frac{(aM_{\mathrm{eff}})}{a}^.+\mathrm{\Delta }(a)\right]X_k=0,$$
(9)
where $`\mathrm{\Delta }(a)[\frac{1}{4}\left(\frac{\dot{a}}{a}\right)^2\frac{1}{2}\left(\frac{\ddot{a}}{a}\right)]`$. When $`at^n`$ such as in a matter or radiation dominated universe, $`\mathrm{\Delta }(a)t^2`$ and can be neglected soon after inflation.

Using the mode equation and ansatz, we find

$$𝐮_{k,\pm }(t)=a^1\left[i\dot{X}_ki\frac{1}{2}\left(\frac{\dot{a}}{a}\right)X_k+(\gamma 𝐤M_{\mathrm{eff}})X_k\right]R_\pm (𝐤).$$
(10)
Taking the charge conjugate of $`𝐮_{k,\pm }`$, we find

$$𝐯_{k,\pm }(t)=a^1\left[i\dot{X}_k^{}+i\frac{1}{2}\left(\frac{\dot{a}}{a}\right)X_k^{}(\gamma 𝐤+M_{\mathrm{eff}})X_k^{}\right]\overline{R}_\pm (𝐤),$$
(11)
where $`\overline{R}_\pm (𝐤)i\gamma ^2R_\pm ^{}(𝐤)`$ is an eigenvector of helicity such that $`𝐤𝚺\overline{R}_\pm (𝐤)=\pm 1`$ and $`\gamma ^0\overline{R}_\pm (𝐤)=1`$.

The energy-momentum tensor is obtained from the ($`\psi `$ and $`\overline{\psi }`$ symmetrized) matter action by variation with respect to the vierbein

$$T^{\mu \nu }=\frac{i}{2}\left[\overline{\psi }\overline{\gamma }_{(\mu }\underset{\nu )}{\overset{}{D}}\psi \overline{\psi }\underset{(\mu }{\overset{}{D}}\overline{\gamma }_{\nu )}\psi \right],$$
(12)
and the Hamiltonian operator is

$$_D=d^3x\left[i\psi ^{}\dot{\psi }\right].$$
(13)
In general, if this Hamiltonian is diagonal in the annihilation (and creation) operators $`\widehat{a}_{k,s}`$ ($`\widehat{a}_{k,s}^{}`$) and $`\widehat{b}_{k,s}`$ ($`\widehat{b}_{k,s}^{}`$) at $`t=0`$, it will not be for later times. This is the signature of particle creation due to the time dependent background. In order to determine the number of particles produced, we perform a Bogliubov transformation on the creation and annihilation operators so as to diagonalize the Hamiltonian at time $`t`$.

For the problem of particle creation quite often it is useful to represent the wave function in the adiabatic ( semi-classical, WKB) form

$$X_k(t)=\alpha _kN_+e^{i_0^t𝑑t\mathrm{\Omega }_k(t)}+\beta _kN_{}e^{+i_0^t𝑑t\mathrm{\Omega }_k(t)},$$
(14)
where $`N_\pm =(2\mathrm{\Omega }_k(\mathrm{\Omega }_k\pm M_{\mathrm{eff}}))^{1/2}`$ and the coefficients $`\alpha _k`$ and $`\beta _k`$ correspond to the coefficients of the Bogliubov transformation.

Once the Bogliubov transformation is done, we may write the comoving number density of particles $`n_k(t)=|\beta _k|^2`$ in a given spin state through the solutions of the mode equation (9)

$$n_k(t)=a\left(\frac{\mathrm{\Omega }_kM_{\mathrm{eff}}}{2\mathrm{\Omega }_k}\right)\left[|\dot{X}_k|^2+\mathrm{\Omega }_k^2|X_k|^22\mathrm{\Omega }_kIm(X_k\dot{X}_k^{})\right],$$
(15)
where $`M_{\mathrm{eff}}(m_\psi +h\varphi )`$ as before and $`\mathrm{\Omega }_k^2\frac{k^2}{a^2}+M_{\mathrm{eff}}^2`$. The energy density in these particles is then

$$\rho _\psi (t)=\frac{2}{a^3}\frac{d^3k}{(2\pi )^2}\mathrm{\Omega }_kn_k(t)=\frac{1}{a^3\pi }𝑑kk^2\mathrm{\Omega }_kn_k(t).$$
(16)
The normalization of the solutions $`X_k(t)`$ is such that $`X_k(t0^{})=N_+e^{i\mathrm{\Omega }_kt}`$ and $`n_k(0)=0`$. Thus, we find $`N_+=(2\mathrm{\Omega }_k(\mathrm{\Omega }_kM_{\mathrm{eff}}))^{1/2}`$. These are the so-called positive frequency initial conditions.

A comment about the regularization of fermion VEVs. In principle, fermionic VEVs like $`T^{\mu \nu }`$ require regularization , which in the presence of background metrics and scalars is rather non-trivial, see e.g.  and references therein. We will consider only the processes of particle creation, ignoring vacuum polarization. The creation of particles, which corresponds to the imaginary part of the effective action, has no formal divergencies and does not require regularization. Assuming the particle creation process dominates over vacuum polarization, we will not consider the issues of regularization.

## III  Fermionic Preheating without Expansion of the Universe

It is convenient to begin the investigation of fermionic preheating due to an oscillating scalar field with a simplified setting neglecting the expansion of the universe. This setting may have not only methodological advantages. Indeed, whenever the frequency of $`\varphi `$ oscillations is much greater than the rate of cosmic expansion $`H`$, it is sensible to neglect the time dependence of the scale factor in solving the mode equation (9). This can be the case, for example, when spontaneous symmetry breaking occurs rapidly leaving the $`\varphi `$-field oscillating about a new minimum where the effective mass $`m_\varphi `$ happens to be much larger than the Hubble parameter at the time of the transition. One such example is hybrid inflation scenarios which end with TeV scale energy densities. Another example is the inflationary model with the potential $`\lambda \varphi ^4`$. This theory possesses conformal properties: at the stage of inflaton oscillations equations for the fields by means of conformal transformations can be reduced to the equations in Minkowskii space-time, see e.g. .

In this paper we will mostly use a chaotic inflationary model with quadratic potential $`V(\varphi )=\frac{1}{2}m_\varphi ^2\varphi ^2`$. If we make the replacement $`a=1`$ in all the formulas of Section 2, all the effects of expansion will be removed. Background oscillations take the form of harmonic oscillations $`\varphi (t)=\varphi _0f(t)`$ with $`f(t)=\mathrm{cos}(m_\varphi t)`$ and $`\varphi _0`$ the time independent amplitude. It is convenient to define a new, dimensionless time variable $`\tau m_\varphi t`$. With this change of variables, the mode equation (9) may be written

$$X_k^{\prime \prime }+[\kappa ^2+(\stackrel{}{m}+\sqrt{q}f)^2i\sqrt{q}f^{}]X_k=0,$$
(17)
where we have introduced the dimensionless momentum $`\kappa \frac{k}{m_\varphi }`$, the dimensionless fermion mass $`\stackrel{}{m}\frac{m_\psi }{m_\varphi }`$, and the resonance parameter $`q\frac{h^2\varphi _0^2}{m_\varphi ^2}`$. These three parameters completely determine the strength of the effect. In fact, this form of the mode equation is valid not only for harmonic background oscillations, but for generic $`\varphi `$ oscillations in a general potential. For this, we identify $`m_\varphi `$ with the frequency of oscillation, $`\varphi _o`$ with its amplitude, and $`f(\tau )`$ with the periodic background oscillations normalized to unit amplitude. Note that the frequency will be amplitude-dependent for a general, non-quadratic potential.

### A  Parametric Excitation of Fermions

If individual inflatons at rest are decaying into light fermions in a process $`\varphi \overline{\psi }\psi `$, perturbative calculations give the rate of decay $`\mathrm{\Gamma }_{\varphi \psi \psi }\frac{h^2m}{8\pi }`$ and the spectrum of created fermions is sharply peaked around $`m/2`$ with the width $`\mathrm{\Gamma }_{\varphi \psi \psi }^1`$. In Figs. 2, 2, however, we plot the time-dependence and spectrum of occupation number for fermions created from the coherently oscillating inflaton field, as follows from numerical solution of Eq. (17) and (15).

The spectrum and time evolution are drastically different from what is expected from the perturbative calculations. We therefore can talk about a specific phenomena, the parametric excitation of fermions interacting with coherent background oscillations.

It is instructive to compare the spectrum and evolution of fermionic occupation number with those of bosonic occupation number. The mode function of a quantum bose scalar field $`\widehat{\chi }`$ coupling as $`g^2\chi ^2\varphi ^2`$ with oscillating inflaton obeys the bosonic oscillator-like equation

$$X_{b,k}^{\prime \prime }+\left[\kappa ^2+\stackrel{}{m_b}^2+q_bf^2\right]X_{b,k}=0,$$
(18)
where $`q_b=\frac{g^2\varphi _0}{m^2}`$, $`\stackrel{}{m_b}=\frac{m_\chi }{m}`$. In Figs. 4 and 4 we plot the spectrum and time evolution of bosonic occupation number calculated from bosonic oscillator-like equations (18). We use the familiar model of bosonic resonance due to the self-interaction in $`\lambda \varphi ^4`$ inflation, which corresponds to $`g^2=3\lambda `$, $`q_b=3`$, $`m_b=0`$, $`f(\tau )`$ is given by oscillations in $`\lambda \varphi ^4`$ theory.

In the bosonic case, there are distinct resonance bands in which bosonic modes are exponentially unstable, $`n_k^b(t)e^{\mu _k\tau }`$. Particle creation takes place outside of the resonant bands as well, although there the occupation number is bounded and oscillates periodically. In the fermionic case $`n_k1`$ is always bounded by Pauli blocking and oscillating periodically with time. The occupation number in both cases is changing in time. It is therefore convenient to introduce an envelope function $`F_k`$ of the particle spectrum , which corresponds to $`n_k`$ averaged over short-time intervals (order of the background oscillation period). A bosonic envelope function cannot be defined for the resonant bands (where it is $`e^{\mu _k\tau }`$.) This zone corresponds to the gap between almost vertical lines in Fig. 4. The fermionic envelope function $`F_k`$, on the contrary, can be defined everywhere. Although it is always bounded by unity, their structure is reminiscent of the resonant-band structure: for some range of $`k`$, the ‘resonant’ band, $`n_k`$ is close to unity, while in other ranges, the stable bands, it is significantly smaller if not zero. In  different levels of $`F_k`$ were plotted on the parameter plane $`(\kappa ,q)`$, revealing a structure which reminiscent of the stability/instability chart of the bosonic parametric resonance. One of the most important results is that fermionic parametric excitations occurs very quickly, within about ten(s) background oscillations. Interestingly, fermionic ”resonant bands” are excited the last, while non-resonant intervals fill first.

Comparison of bosonic and fermionic parametric excitations is useful to understand some features of bosonic preheating. As we have seen, there is production of bosons outside of the resonance band. If we take into account the expansion of the universe, in the most interesting case of large $`q_b`$, the difference between resonant and non-resonant excitations of bosons will be erased, and the regime of stochastic resonant production of bosons will be settled down .

### B  Some Generic Analytic Results

Dynamics of the Fermi field coupling to the background homogeneous scalar can be revealed with the second-order oscillator-like equation (17) for the mode function $`X_k(t)`$. For periodic background oscillations $`f(t)=f(t+T)`$, $`T`$ is a period, some generic analytic results were derived a long time ago  in the context of particle creation in a periodic external electromagnetic field. In particular, the occupation number of created particles at instances $`t=N_sT`$, i.e. exactly after $`N_s`$ background oscillations, is given by expression

$$n_k(N_sT)=\frac{k^2}{2\mathrm{\Omega }_k^2}\frac{\mathrm{sin}^2N_sd_k}{\mathrm{sin}^2d_k}\left(\mathrm{𝐼𝑚}X_k^{(1)}(T)\right)^2,$$
(19)
where $`\mathrm{cos}d_k=\mathrm{𝑅𝑒}X_k^{(1)}(T)`$. To derive this result, one introduces two fundamental solutions of Eq. (17), $`X_k^{(1)}(t)`$ and $`X_k^{(2)}(t)`$; such that initially $`X_k^{(1)}(0)=1`$, $`\dot{X}_k^{(1)}(0)=0`$ and $`X_k^{(2)}(0)=0`$, $`\dot{X}_k^{(2)}(0)=1`$. Expression (19) involves only the value of the first fundamental solution $`X_k^{(1)}(T)`$ exactly after the first oscillation. It says that the occupation number of created particles after $`N_s`$ background oscillations is modulated with a certain frequency $`\nu _k`$ (which, as we will see, does not coincide with $`d_k`$). However, practical application of the generic formula (19) is rather limited, because it does not address the full time evolution of $`n_k(t)`$, and cannot strictly determine a period of modulation $`\pi /\nu _k`$.

To get an idea of how the occupation number of created particles $`n_k(t)`$ evolves with time, again let us look at Fig. 2 and further at Fig. 7 for different values of the parameters $`\kappa `$ and $`q`$. For small and moderate $`q`$ (but not too small $`\kappa `$) the occupation number exhibits high frequency (period $`<\frac{T}{2}`$) oscillations which are modulated by a long period behavior. For large $`q`$ number of fermions jumps in a step-like manner at instances when effective mass of the fermi field crosses zero, superposed by very high frequency oscillations around almost constant values, as depicted in Fig. 7. These jumps are modulated with a frequency $`\nu _k/2`$: steps in the first half of the cycle up are accumulated until $`n_k(t)`$ reaches its maximum $`F_k`$, and then steps down to zero in the second half of the cycle.

However, this picture of high frequency features superimposed over long-period modulation is not universal. For small $`\kappa `$, or for one of the most interesting cases of $`q1`$ and moderate $`\kappa `$, the occupation number of fermions jumps between zero and one within time interval much shorter than period of background oscillations, as depicted by the dotted curve in Fig. 2. There are interesting situations when fermions are created in an instant (single kick) process, where formula (19) is not applicable. Therefore, for different ranges of parameters we will shall develop different approaches.

### C  Semi-Analytic Theory for Averaged Occupation Number

Numerical curves for the mode functions suggests splitting of the time evolution into higher frequency features with the time-scale comparable or less than the period of background oscillations $`T`$, and low-frequency modulations with the period $`\pi /\nu _k`$ greater than $`T`$. In this case we can utilize generic result (19). However, it is convenient to use not the values $`n_k(N_sT)`$, but rather the smoothed occupation number $`\overline{n}_k(t)`$ which is $`n_k(t)`$ averaged over high frequency oscillations, $`\overline{n}_k(\tau )=\frac{1}{T}_\tau ^{(\tau +T)}𝑑\tau n_k(\tau )`$. Then we can write the smoothed occupation number of fermions in a factorized form

$$\overline{n}_k(\tau )=F_k\mathrm{sin}^2\nu _k\tau ,$$
(20)
where we introduce an envelope function $`F_k`$. The average occupation number of fermions evolves periodically with time. The spectrum of $`\overline{n}_k`$ can be characterized by the envelope function $`F_k`$ and the period of modulation $`\frac{\pi }{\nu _k}`$ which depends also on the parameter $`q`$.

Now we will utilize the result (19). We found that the envelope function can be extrapolated by the factor in the front of $`\mathrm{sin}^2N_sd_k`$ in (19) and given by the expression

$$F_k=\frac{1}{\mathrm{sin}^2\nu _kT}\frac{\kappa ^2}{2\mathrm{\Omega }_k^2}\left(\mathrm{𝐼𝑚}X_k^{(1)}(T)\right)^2.$$
(21)
Next is to determine the frequency $`\nu _k`$ of the $`\overline{n}_k`$ modulations. The value $`d_k`$ defined after (19) cannot be the right answer, because it would incorrectly predict that the peaks of $`F_k`$ are filled up first, while actually they are filled last. We tried the combination $`\nu _k=\pi /2d_k`$, because it corresponds to correct order of saturation of $`F_k`$, and it occurs to work well. Therefore, the modulation frequency $`\nu _k`$ is given by the relation $`\mathrm{cos}\nu _kT=\mathrm{𝑅𝑒}X_k^{(1)}(T)`$. Thus, to find $`F_k`$ and $`\nu _k`$, one need only calculate the complex value $`X_k^{(1)}(T)`$ after a single background oscillation, instead of performing a full numerical integration of Eq. (17).

We calculated $`X_k^{(1)}(T)`$ numerically for $`\frac{1}{2}m^2\varphi ^2`$ background model and constructed the envelope function $`F_k`$ plotted in Fig. 6 (similar graph was plotted in  for $`\frac{1}{4}\lambda \varphi ^4`$ model). In Fig. 2 we show, using (21), how the fermionic resonance bands are filled after $`10`$ background oscillations.

The function $`\nu _k`$ gives us the time scale for fermion excitation. In Fig. 6 we plot the period of modulation $`\frac{\pi }{\nu _k}`$ as a function of $`k`$. This function is peaked where $`F_k`$ is peaked, i.e. the peaks of the resonance curve are the last to fill.

### D  Method of Successive Scatterings for Fermions

It turns out that for large values of the parameter $`q`$ we can significantly advanced in calculations of $`n_k(t)`$ beyond the results of Section (III C). One can expect that $`q`$ may be much greater that one. Indeed, $`q=\frac{h^2\varphi _0^2}{m^2}`$. In the context of chaotic inflation with the potential $`V(\varphi )=\frac{1}{2}m^2\varphi ^2`$ it follows from the theory of cosmological perturbations that $`\frac{\varphi _0^2}{m^2}10^{12}`$. On the other hand one can admit that Yukawa coupling can be $`h10^6`$, which provides $`q1`$.

We will generalize for fermions the method of parabolic scatterings introduced for the bosonic resonance in reference . The method is based on the observation that, for $`q1`$, change of the particle number occurs only during a short time interval $`\tau _{}`$ near the zeros of the effective mass of the particles. This occurs because equation (15) for the number of particles in terms of the mode functions is an adiabatic invariant of the mode equation.

For large $`q`$, the adiabaticity condition $`\mathrm{\Omega }_k^{}<\mathrm{\Omega }_k^2`$ is violated only near the times $`\tau _{}`$ when the effective mass vanishes. This leads to a step-like evolution in the number of fermions. This is illustrated in figure (7) which shows the evolution of $`n_k(\tau )`$ and $`M_{\mathrm{eff}}(\tau )=h\varphi (\tau )+m_\psi `$

Near the times $`\tau _{}`$, we may approximate the effective mass in equation (17) by $`(\stackrel{}{m}+\sqrt{q}f)=\frac{1}{m_\varphi }(m_\psi +h\varphi (\tau ))\frac{h}{m_\varphi }\varphi (\tau _{})^{}(\tau \tau _{})+O((\tau \tau _{})^2)`$. We may also write this as

$$(\stackrel{}{m}+\sqrt{q}f)\sqrt{q}f^{}(\tau _{})(\tau \tau _{})=\pm (q\stackrel{2}{\stackrel{}{m}})(\tau \tau _{}),$$
(22)
where the sign on the left hand side depends on whether $`f(\tau )`$ is increasing ($`+`$) or decreasing ($``$) at $`\tau _{}`$. Thus, in the neighborhood of $`\tau _{}`$, the mode equation becomes

$$X_k^{\prime \prime }+[\kappa ^2+(q\stackrel{2}{\stackrel{}{m}})^2(\tau \tau _{})^2i\mathrm{sgn}(f^{})\sqrt{q\stackrel{2}{\stackrel{}{m}}}]X_k=0,$$
(23)
which is a Schroedinger equation for scattering off a negative parabolic potential centered at $`\tau _{}`$.

The method of parabolic scattering uses the exact solution of equation (23) to provide a connection formula between the adiabatic approximations of the full mode equation on either side of $`\tau _{}`$. Suppose the $`j`$-th zero of $`M_{\mathrm{eff}}`$ occurs at time $`\tau _j`$. For times between $`\tau _{j1}<\tau <\tau _j`$, the general solution of the mode equation (17) takes the adiabatic (or semi-classical) form

$$X_k^j(\tau )=\alpha _k^jN_+e^{i_0^\tau 𝑑\tau \mathrm{\Omega }_k(\tau )}+\beta _k^jN_{}e^{+i_0^\tau 𝑑\tau \mathrm{\Omega }_k(\tau )},$$
(24)
where the coefficients $`\alpha _k^j`$ and $`\beta _k^j`$ are constant for $`\tau _{j1}<\tau <\tau _j`$. After the scattering at $`\tau _j`$, $`X_k(t)`$, within the interval $`\tau _j<\tau <\tau _{j+1}`$, again has the adiabatic form

$$X_k^{j+1}(\tau )=\alpha _k^{j+1}N_+e^{i_0^\tau 𝑑\tau \mathrm{\Omega }_k(\tau )}+\beta _k^{j+1}N_{}e^{+i_0^\tau 𝑑\tau \mathrm{\Omega }_k(\tau )},$$
(25)
with new coefficients $`\alpha _k^{j+1}`$ and $`\beta _k^{j+1}`$ that are constant for $`\tau _j<\tau <\tau _{j+1}`$. At $`\tau =0`$, our vacuum positive frequency condition requires $`\alpha _k^1=1`$ and $`\beta _k^1=0`$. Particle creation occurs when, after scattering at the times $`\tau _j`$, the initial positive frequency wave acquires a negative frequency part. The number density of produced particles with momentum $`k`$ is $`n_k^j=|\beta _k^j|^2`$ for times $`\tau _j<\tau <\tau _{j+1}`$. Furthermore, normalization requires $`|\alpha _k^j|^2+|\beta _k^j|^2=1`$ for all $`j`$.

The important observation is that the outgoing amplitudes ($`\alpha _k^{j+1}`$, $`\beta _k^{j+1}`$) can be expressed through the incoming amplitudes ($`\alpha _k^j`$, $`\beta _k^j`$) by means of the reflection $`R_k`$ and transmission $`D_k`$ amplitudes for scattering at $`t_j`$:

$$\left(\begin{array}{c}\alpha _k^{j+1}e^{i\theta _k^j}\\ \beta _k^{j+1}e^{+i\theta _k^j}\end{array}\right)=\left(\begin{array}{cc}\frac{1}{D_k}& \frac{R_k^{}}{D_k^{}}\\ \frac{R_k}{D_k}& \frac{1}{D_k^{}}\end{array}\right)\left(\begin{array}{c}\alpha _k^je^{i\theta _k^j}\\ \beta _k^je^{+i\theta _k^j}\end{array}\right).$$
(26)
Here $`\theta _k^j=\underset{0}{\overset{\tau _j}{}}𝑑\tau \mathrm{\Omega }(\tau )`$ is the phase accumulated by the moment $`\tau _j`$.

If we let $`x=(q\stackrel{2}{\stackrel{}{m}})^{1/4}\tau `$, equation (23) may be written

$$\frac{d^2X_k}{dx^2}+\left[\mathrm{\Lambda }_k^2i(1)^j+x^2\right]X_k=0,$$
(27)
where the parameter

$$\mathrm{\Lambda }_k^2\frac{\kappa ^2}{|hf^{}(\tau _{})|}=\frac{\kappa ^2}{\sqrt{q\stackrel{2}{\stackrel{}{m}}}}.$$
(28)
A general analytic solution of equation (27) is a linear combination of the parabolic cylinder functions $`W(\frac{(\mathrm{\Lambda }_k^2i(1)^j)}{2};\pm \sqrt{2}x)`$. The reflection $`R_k`$ and transmission $`D_k`$ amplitudes for scattering on the parabolic potential can be found from the asymptotic forms of these analytic solutions:

$$R_k=\frac{ie^{i\phi _k}}{\sqrt{1e^{\pi \mathrm{\Lambda }_k^2}}},D_k=\frac{e^{i\phi _k}}{\sqrt{1e^{\pi \mathrm{\Lambda }_k^2}}},$$
(29)
where the angle $`\phi _k`$ is $`\phi _k=\mathrm{arg}\mathrm{\Gamma }\left(\frac{1+i\mathrm{\Lambda }_k^2}{2}\right)+\frac{\mathrm{\Lambda }_k^2}{2}\left(1+\mathrm{ln}\frac{2}{\mathrm{\Lambda }_k^2}\right)`$. Note that the angle $`\phi `$ depends on the momentum $`k`$. Substituting (29) into (26), we can determine the change in induced in the $`\alpha _k`$ and $`\beta _k`$ coefficients by a single parabolic scattering in terms of the parameters of the parabolic potential and the phase $`\theta _k^j`$ only. Specifically, we find

$`\left(\begin{array}{c}\alpha _k^{j+1}\\ \beta _k^{j+1}\end{array}\right)=\left(\begin{array}{cc}\sqrt{1e^{\pi \mathrm{\Lambda }_k^2}}e^{i\phi _k}& (1)^je^{\frac{\pi }{2}\mathrm{\Lambda }_k^2+2i\theta _k^j}\\ (1)^je^{\frac{\pi }{2}\mathrm{\Lambda }_k^22i\theta _k^j}& \sqrt{1e^{\pi \mathrm{\Lambda }_k^2}}e^{i\phi _k}\end{array}\right)\left(\begin{array}{c}\alpha _k^j\\ \beta _k^j\end{array}\right).`$ (30)
It is now a simple matter to find the change in particle number after one scattering. From the normalization of $`\alpha _k`$ and $`\beta _k`$ we have the relation

$$\begin{array}{ccccccc}\hfill \left(\begin{array}{cc}\alpha _{}^{}{}_{k}{}^{j+1}& \beta _{}^{}{}_{k}{}^{j+1}\end{array}\right)& \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)& \left(\begin{array}{c}\alpha _k^{j+1}\\ \beta _k^{j+1}\end{array}\right)\hfill & =& |\alpha _k|^2|\beta _k|^2& =& 12n_k^{j+1}\end{array},$$
(31)
which, when applied to equation (30) gives the desired result

$`n_k^{j+1}`$ $`=`$ $`e^{\pi \mathrm{\Lambda }_k^2}+\left(12e^{\pi \mathrm{\Lambda }_k^2}\right)n_k^j`$ (32)
$``$ $`2(1)^je^{\frac{\pi }{2}\mathrm{\Lambda }_k^2}\sqrt{1e^{\pi \mathrm{\Lambda }_k^2}}\sqrt{n_k^j(1n_k^j)}\mathrm{sin}\theta _{tot}^j,`$ (33)
where the phase $`\theta _{tot}^j=2\theta _k^j\phi _k+\mathrm{arg}\beta _k^j\mathrm{arg}\alpha _k^j`$. Let us discuss this formula, which is the main result of our paper. If fermions are light, $`m_\psi =0`$, their occupations number is changing with time only when background field crosses zeros. Without expansion of the universe, the phases $`\theta `$ accumulated between successive zeros are equal. In this case one can try to proceed to find the solution of the matrix equation (26), as it was done for bosons . However, for bosons that was needed to find the stability/instability bands, which is not so interesting for fermions. In the case of massive fermions time intervals between successive zeros of $`M_{eff}`$ are not equal, and problem of finding matrix solution of eq. (26) became even more complicated. However, in the most interesting case of expanding universe the phases between successive zeros of $`M_{eff}`$ became random (see  for details). Then we can just put $`\theta _{tot}`$ as a random phase and use formula (33) as it is.

## IV  Parametric Excitation of Fermions with Expansion of the Universe

To address the problem of fermionic preheating after $`m^2\varphi ^2`$ chaotic inflation, we must now deal with the full mode equation (9), which no longer has a periodic time dependence in the complex frequency. Nevertheless, it is still convenient to work with the form (17) in the time variable $`\tau =m_\varphi t`$ where the parameters $`q`$ and $`\kappa ^2`$ are now understood to be time dependent:

$$X_k^{\prime \prime }+\left[\frac{\kappa ^2}{a^2}+M_{\mathrm{eff}}^2i\frac{(aM_{\mathrm{eff}})^{}}{a}+\mathrm{\Delta }(a)\right]X_k=0,$$
(34)
Check where  stands for the derivative in respect with $`\tau `$, $`\kappa `$ is a comoving momentum scaled in units of $`m_\varphi `$. There is often used another from of the fermionic mode equation written in terms of conformal time $`\eta =𝑑t/a`$ and mode function $`Y_k=a^{3/2}\chi _k`$

$$_\eta ^2Y_k+\left[\kappa ^2+M_{\mathrm{eff}}^2i_\eta (aM_{\mathrm{eff}})\right]Y_k=0.$$
(35)
This form of equation is useful, for instance, in conformal theory $`V(\varphi )=\lambda \varphi ^4`$ when the problem can be reduced to the problem in Minkowski space-time . In this case background field is oscillating periodically in respect with time $`\eta `$. In the case of quadratic inflaton potential background field is oscillating with the constant period in terms of physical time $`t`$ (or $`\tau `$). We therefore will use mode equation (34). The parameter $`q`$ is now understood to be time dependent. Specifically, we have $`q(\tau )\frac{h^2\mathrm{\Phi }^2(\tau )}{m_\varphi ^2}`$, scaled physical momentum $`p\frac{1}{a}\frac{k}{m_\varphi }`$ Here, $`\mathrm{\Phi }(\tau )`$ is the time dependent amplitude of inflaton oscillations. As is well known (c.f. ), oscillations of $`\varphi `$ in this model quickly approach the asymptotic solution $`\varphi (\tau )=\mathrm{\Phi }(\tau )\mathrm{cos}(\tau )`$, $`\mathrm{\Phi }(\tau )\frac{\varphi _0}{a^{3/2}}\frac{M_{\mathrm{pl}}}{\sqrt{3\pi }(1+\tau )}`$ where time is measured from the start of inflaton oscillations, $`\varphi _0\frac{M_{\mathrm{pl}}}{\sqrt{3\pi }}`$, and the scale factor averaged over several oscillations behaves as in a matter dominated universe: $`a(\tau )=(1+\tau )^{2/3}`$.

### A  Stochastic Parametric Excitations of Fermions

Let us first consider the case of light fermions with a large initial resonance parameter: $`h\varphi _0m_\varphi m_\psi `$. Fig. (9) shows a numerical solution of equation (17) in the absence of expansion with a resonance paramter $`q=10^6`$. As expected from for such a large $`q`$ parameter, we see a step-like change in occupation number with periodic modulation.

In Fig. (9), we show a numerical solution to the equation (34) in the presence of expansion with the initial parameters $`q_0=10^6`$. Here, the comoving number density of particles is plotted. We see that, as for the non-expanding case, evolution between the zeros of the effective mass (oe equally inflaton field), $`M_{\mathrm{eff}}\sqrt{q(\tau )}\mathrm{cos}(\tau )`$, is adiabatic. We can thus apply the results of Section III,D in particular the analytic formula (33), mutatis mutandis. The Pauli principle is obviously still obeyed and the typical step in particle number is still suppressed by a factor $`e^{\pi \mathrm{\Lambda }_k}`$ which is now time dependent. The most important qualitative change is that the accumulated phase $`\theta _{\mathrm{tot}}`$ is now uncorrelated between successive zeros of the effective mass (inflaton field). This occurs because the effective frequency $`\mathrm{\Omega }_k=\sqrt{p^2+q^2(\tau )\mathrm{cos}^2\tau }`$ is no longer periodic and the accumulated phase $`\theta =_j^{j+1}\mathrm{\Omega }_k\frac{2h\mathrm{\Phi }(\tau )}{m}+O(k^2)`$ changes substantially in magnitude within one inflaton oscillation, $`\delta \theta _k\frac{\sqrt{q}}{2N_s^2}`$, after the $`N_s`$-th oscillation . The result is that the $`\mathrm{sin}(\theta _{\mathrm{tot}})`$ term in (33) becomes a random variable. As is readily apparent in Fig. (9), this destroys the periodic modulation of $`n_k`$ and the parametric excitation of fermions becomes stochastic, as anticipated in .

Once the periodic modulation of $`n_k`$ is destroyed, the construction of Section III,C is no longer valid: the occupation number cannot be characterized by an amplitude and period. In fact, the spectrum of created particles is even simpler. Stochastic excitation allows a given comoving mode $`\kappa `$ to obtain any amplitude in the range $`0n_k1`$ if there are a sufficient number of parabolic scatterings, $`N_s`$. This gives us the picture of stochastically filling a (Fermi) sphere in the momentum space. Numerical calculations confirm this picture, see Fig. 10. Let us find its radius $`\kappa _s`$. A comoving mode will be excited if $`\mathrm{\Lambda }_k(\tau )1/\pi `$. Since $`\mathrm{\Lambda }_k(\tau )=\frac{\kappa ^2}{a^2\sqrt{q(\tau )}}`$, we have

$$\frac{\kappa _s^2}{a^2\sqrt{q(\tau )}}=a^{1/2}\frac{\kappa _s^2}{\sqrt{q_0}}1/\pi .$$
(36)
Therefore for the light fermions, comoving radius of excited modes is increasing with time as $`\kappa _sa^{1/4}`$. The radius of the sphere is scaled as $`\kappa _sq_0^{1/4}a^{1/4}`$.

Due to the expansion, the amplitude $`\mathrm{\Phi }`$ is decreasing and the parameter $`q(\tau )`$ drops. At some moment the Fermi sphere will quit expansing. Once $`q(\tau )`$ is order of unity, the excitation is no longer strong, and redshifts of fermion modes will be fast enough to prevent from parametric excitation. Eventually fermions will be produced in the perturbative regime. If fermionic mass is non-zero, the perturbative decay goes unless $`2m_\psi >m_\varphi `$.

### B  Analytic Results for Production of Supermassive Fermions

One of the most significant differences between bosonic and fermionic parametric amplification is the possibility to create superheavy fermions from much liter inflatons, if inflatons are oscillating coherently. As we seen, for large $`q`$, massive fermions are created not when inflaton field itself crosses zero, but when the combination $`M_{eff}=h\varphi +m_\psi `$ crosses zero. This is because at those instances even very heavy fermions are effectively massless, because their bare mass is “compensated” by large value of $`h\varphi `$, if $`\varphi `$ can have a large amplitude. This is exactly the case when $`\varphi `$ is the inflaton field. For instance, in the chaotic inflationary scenario amplitude of its oscillations immediately after inflation can be as large as $`0.1M_P`$.

In this case we again have a situation when mode function of heavy fermions has WKB form between zeros of $`M_{eff}`$, and can be described with parabolic scattering around this instance. Therfore our general formula (33) works for superheavy fermions as well. The only modification is that the phase between successive scatterings will be defined by the integral over time intervals between them. In formula (33) parameter $`\mathrm{\Lambda }_k(\tau )`$ will be

$$\mathrm{\Lambda }_k(\tau )=\frac{\pi \kappa ^2}{a^2\sqrt{\frac{q_0}{a^3}\stackrel{2}{\stackrel{}{m}}}}$$
(37)
In case of massive fermions, the criteria for excitation instead of (36) will be

$$\frac{\kappa _s^2}{a^2\sqrt{\frac{q_0}{a^3}\stackrel{2}{\stackrel{}{m}}}}1/\pi .$$
(38)
At the start of background oscillations for large $`q_0`$ a term $`\stackrel{2}{\stackrel{}{m}}`$ is negligible. Thus, as in the last section, we expect the Fermi sphere to fill with the radius $`\kappa _s=\sqrt{q_0}a^{1/4}`$. However, as $`q`$ drops, equation (38) solved for $`\kappa _s^2`$ reaches a maximum, $`\kappa _m`$. and then decreases. This maximum occurs when the scale factor is $`a_m=\left(\frac{q_0}{4\stackrel{2}{\stackrel{}{m}}}\right)^{1/3}`$. This gives $`\kappa _m^2=\frac{\sqrt{3}}{\pi }\left(\frac{q_0}{4\sqrt{\stackrel{}{m}}}\right)^{2/3}`$. Notice that $`k_m`$ scales as $`q_0^{1/3}`$. This is still not the final width of the band. This is because there are a number of background oscillations which occurs between the moment $`a_m`$ and the moment when denominator in (37) approaches zero. We put a number of oscillation before $`a_m`$ and between $`a_m`$ and the moment of zero denominator as approximately equal and equal to $`N_m`$. In this time after $`a_m`$ even the modes suppressed exponentially as $`e^{\mathrm{\Lambda }_k}`$ can achieve significant occupation numbers in their random stochastic walk. We found the amplification factor in the from of the exponent is $`\sqrt{N_m}`$. In contrary to the ligh fermions case, there is no perturbative end of the process. Excitation of superheavy fermions is abruptly terminated when $`\frac{q_0}{a^3}\stackrel{2}{\stackrel{}{m}}`$.

We found the spectrum of created superheavy particles after the process is terminated is given by formula

$$n_k=\frac{1}{2}\mathrm{exp}\left(2\frac{(\kappa \gamma \kappa _m)^2}{\kappa _m^2}\right)$$
(39)
where $`2\gamma ^2=\mathrm{ln}\left(\frac{q_0}{4\pi ^2\stackrel{2}{\stackrel{}{m}}}\right)`$. This formula is valid for $`N_m`$ greater than a few. In Fig. 12 we plot numerically calculated final spectrum of the supermassive fermions vs. an analytic formula 39, which are in a good agreement. In the case of superheavy particles in an expanding universe the maximum radius of the k-sphere is scaled with $`q_0`$ as $`\gamma k_mq_0^{1/3}\mathrm{ln}q_0`$, which is in close agreement with . However, one shall also check that at the moment when excitation of superheavy fermions terminates, their back reaction $`h\overline{\psi }\psi `$ to the dynamics of inflaton oscillation is still negligible.

In Fig. 12 we plot spectrum of superheavy fermions after the process is terminated. Dashed curves are obtained with analytic formula (39).

### C  Fermionic Production from Single Kick

As we seen formula (33) can be extended to the case of expanding universe. Consider the very first term in formula (33). As It corresponds to the creation of fermions after effective mass $`M_{eff}(t)`$ crosses zero for the first time. In this case we found

$`n_k=e^{\pi \mathrm{\Lambda }_k^2}`$ (40)
Occupation number of fermions generated from the single kick is plotted in Fig. 12. Two curves corresponding to the formula (40) and to the numerical solutions are shown and practically indistinguishable.

## V  Discussion

In this paper we developed a non-perturbative theory for fermion production by an oscillating inflaton field. As we have seen, the production of fermions can be characterized as parametric excitation. Even in the simple model of a Yukawa coupling between fermions and a background scalar field in an expanding universe oscillating around the minimum of its quadratic potential, the theory of fermionic parametric excitation is rich and leads to important results. Fermions are created very quickly, within about ten(s) oscillations, in out-of-equilibrium states. For large values of the resonance parameter $`q=\frac{h\varphi _0}{m_\varphi }`$, occupation number of light fermions ($`m_\psi <m_\varphi h\varphi _0`$) is changing in a step-like manner at instances $`t_j`$ when the inflaton amplitude $`\varphi (t_j)=0`$ passes through zero, $`j=1,2,3,..`$ We have developed the method of parabolic scatterings for fermions, based closely on a similar approach for the bosonic resonance . It is possible to derive a unified recursive formula, which relates the occupation number of fermions or bosons $`n_k^{j+1}`$ at the moment $`t_{j+1}`$ to the earlier value $`n_k^j`$ :

$`n_k^{j+1}`$ $`=`$ $`e^{\pi \mathrm{\Lambda }_k^2}+\left(1\pm 2e^{\pi \mathrm{\Lambda }_k^2}\right)n_k^j`$ (41)
$``$ $`2(1)^je^{\frac{\pi }{2}\mathrm{\Lambda }_k^2}\sqrt{1\pm e^{\pi \mathrm{\Lambda }_k^2}}\sqrt{n_k^j(1\pm n_k^j)}\mathrm{sin}\theta _{tot}^j,`$ (42)
For bosons one shall use an upper sign and neglect the $`(1)^j`$ terms , while for fermions one shall use the lower sign. For light fermions and bosons $`\mathrm{\Lambda }_k^2=\frac{k^2}{\sqrt{q}m_\varphi ^2}`$. For large $`q`$ the angle $`\theta _{tot}^j`$ can be treated as a random phase. As a result, formula (42) predicts the stochastic character of parametric excitation of both bosons and fermions. In the fermionic case it leads to the conclusion that in the momentum space fermions chaotically fill up a broad sphere of the radius $`q^{1/4}m`$. This formula also clearly shows the features of: spontaneous emission: for both bosons and fermions, the first oscillation leads to the spectrum $`n_k=e^{\pi \mathrm{\Lambda }_k}`$; stimulated emission for bosons with $`n_k1`$, we have $`(n_k^{j+1}n_k^j)n_k^j`$; Pauli Blocking for fermions, if $`n_k^j=1`$, the next value will always be $`n_k^{j+1}=1e^{\pi \mathrm{\Lambda }_k}`$ even in the stochastic case. This prevents the occupation number of fermions from exceeding $`1`$.

Formula (42) can be extended to the case of massive bosons and fermions. However, here important differences emerge. For bosons we will have $`\mathrm{\Lambda }_k=\frac{k^2+m_b^2}{\sqrt{q}m_\varphi }`$, where $`m_b`$ is the $`\chi `$-boson mass. It leads to the conclusion that the creation of superheavy, $`m_b>m_\varphi `$, bosons is exponentially suppressed. However, for fermions we have $`\mathrm{\Lambda }_k=\frac{k^2}{m_\varphi ^2sqrth^2\varphi _0^2m_\psi ^2}`$. Therefore even superheavy fermions with the mass as large as $`h\varphi _0`$ can be created in abundance from the coherent inflaton oscillations . This occurs because the effective mass of fermions is given by the algebraic combination $`h\varphi (t)+m_\psi `$ and the creation of fermions occurs when the effective mass goes to zero. The method of parabolic scatterings, applied for instances when this happens, leads to both formulae $`(\text{42})`$ with the corresponding $`\mathrm{\Lambda }_k`$.

There are situations where single instance (single kick) of particle creation may lead to interesting effects. An example is the scenario of instant preheating, which is especially important for inflationary models without minima of the inflaton potential . Another example is the interaction of the inflaton with superheavy fermions during the inflationary stage when the combination $`h\varphi (t)+m_\psi `$ can go through zero only once. The generic formula (42) embraces the case where bosons and fermions are created by a single kick. We shall put $`j=0`$, $`n_k^0=0`$, and then $`n_k^1=e^{\pi \mathrm{\Lambda }_k^2}`$. However, instead of the parameter $`q`$ we shall use another combination of coupling constants and $`\varphi _0`$ and $`m_\varphi `$. In this case the effect is defined by the velocity $`\dot{\varphi }_{}`$ at the moment $`t_{}`$, which is different for bosons and fermions. For bosons single instance creation of particles gives the spectrum

$`n_k=e^{\frac{\pi (k^2+m_{\chi ^2})}{g\dot{\varphi }_{}}}`$ (43)
where $`t_{}`$ corresponds to $`\varphi (t_{})=0`$. For fermions single instance creation gives

$`n_k=e^{\frac{\pi k^2}{h\dot{\varphi }_{}}}`$ (44)
where $`\varphi (t_{})+m_\psi /h=0`$.

We believe that the theory of fermionic preheating will be important ingredient of the realistic scenarios of the reheating of the universe after inflation.

Acknowledments. Collaboration and discussions with Andrei Linde and Alexei Starobinsky significantly influenced this work. We also thank Jügen Baacke for useful discussion. This work was supported by NSERC and CIAR.