text:
# Untitled Document

ON ZEROS OF SOLUTIONS OF LINEAR DIFFERENTIAL SYSTEMS.

by Alexei Grigoriev

Weizmann Institute of Science, Rehovot

The subject discussed in this text is the oscillatory behavior of solutions of holomorphic linear time-dependent differential systems of equations. In  a bound is provided for the number of times a trajectory of a polynomial vector field crosses a given algebraic hypersurface. This implies a bound also for linear time dependent vector fields with polynomial dependence on time. A different approach is taken here, giving better bound for some special cases.
This text is a preliminary report.
Besides of having an intrinsic interest, a related problem arises in estimating the number of zeros which the so called complete Abelian integrals may have, the latter problem being itself connected to the problem of estimating the number of limit cycles which can be born from polynomial perturbations of Hamiltonian (polynomial) vector fields on the plane (see, for example, ).
Remark. The bounds are given here for a real segment in the domain where the coefficients of the equations are holomorphic and real on the real axis, but they can be given in principle for any compact domain contained in the domain where the coefficients are holomorphic (). We also remark that the results on which the bounds appearing here are based, are not the best available (see discussion in ).
Let us be given a linear differential system of equations

$$\dot{x}=A(t)x$$

with coefficients holomorphic in some simply connected complex domain. Then the solution will be a vector valued holomorphic function in this domain. The problem in the large is to estimate the number of zeros the first component of the solution may have in that domain.
It is known that if the coefficients of the system belong to the class of polynomials of a given degree, a constructive bound on the number of zeros on some compact subset of the domain may be given which involves just the degree and the maximum modulus of the coefficients .
In the following, another way of obtaining bounds on the number of zeros of the solutions of such systems is discussed, which is the consequence of a simple but somewhat surprising fact.
In , one of the aims is to bound the complex oscillation of solutions of polynomial vector fields; the result for linear differential systems is obtained as a consequence. One writes the system as a polynomial differential system of equations with the coefficients and time treated as extra variables, the polynomial coefficients being integral. Then, by subsequent differentiating, one obtains the derivatives of the first variable written as polynomials with integer coefficients in all the variables of the system. It turns out that one can effectively bound the moment when the ideal generated by the latter polynomials stabilizes. This last derivative is then written as a linear combination of the former derivatives over the ring of polynomials in the variables of the system, which is a linear scalar differential equation for the first coordinate. Since the coefficients of the polynomials are integers, one uses a form of effective Nullstellensatz to bound the coefficients of this equation (on compact regions in the variable space). Since for linear scalar differential equations a bound on the number of zeros exists in terms of the maximum of the coefficients in the region (), a bound on the number of zeros is obtained.
One notes that the obtained scalar differential equation has a space of solutions of dimension much larger than the dimension of the original system, and thus has many solutions not corresponding (at least directly) to the solutions of the system.
In what we do here, we keep the integrality of the polynomial coefficients to be able to use the effective Nullstellensatz. As in , we can parameterize in this way all linear differential systems with coefficients being polynomials of a fixed degree. In the beginning, however, we assume only holomorphic dependence on the time and parameters so that not to obscure the main point.
So let there be given a linear differential system of order $`n`$ $`\dot{x}=A(t,p)x=(a_{ij})(t,p)(x)`$, where $`p`$ stands for parameters. We assume the dependence of the system matrix on $`t,p`$ to be holomorphic in some region $`U`$. By $`W`$ we denote the projection of the domain on the parameter space. We describe a procedure to obtain a scalar differential equation for the first system variable. Differentiating and substituting the relation we have for the first variable, we get:

$$x_1^{\prime \prime }=\underset{i=1}{\overset{n}{}}a_{1i}^{}x_i+a_{1i}x_i^{}=\underset{i=1}{\overset{n}{}}a_{1i}^{}x_i+a_{1i}\underset{i=1}{\overset{n}{}}a_{ij}x_j=\underset{i=1}{\overset{n}{}}(a_{1i}^{}+\underset{j=1}{\overset{n}{}}a_{1j}a_{ji})x_i.$$

Similarly we get the expressions for $`x_1^{(i)}`$ for any $`i`$. We write

$$x_1^{(i)}=a_j^{(i)}x_j=a^{(i)}x.$$

There is a minimal $`k`$, $`1kn`$, such that $`a^{(0)},a^{(1)},..,a^{(k)}`$ are linearly dependent. We then have a unique decomposition

$$a^{(k)}=A_{k1}a^{(k1)}+\mathrm{}+A_0a^{(0)},$$

$`A_i(U)`$, implying that the following equation must be satisfied by any solution for the $`x_1`$ variable of the initial system

$$x_1^{(k)}A_{k1}x_1^{(k1)}\mathrm{}A_0x_1=0.$$

In the sequel we call this differential equation simply ’the derived equation’. Notice that the determinant of a suitable $`k`$-minor of the matrix

$$(a^{(0)},a^{(1)},..,a^{(k1)})$$

is not identically zero, and that multiplying by it the equation, we get an equation with holomorphic coefficients. We thus call it the leading coefficient of the derived equation (corresponding to the minor chosen; of course in the case k=n there is only one such choice).
Claim 1.There exists an analytic set G such that for any $`pWG`$ any solution of the above equation is a solution for the $`x_1`$ variable of the initial system.
Proof. The leading coefficient of the derived equation is a holomorphic function in $`p,t`$, not identically zero. Thus the values of $`p`$ for which the leading coefficient vanishes identically in $`t`$, form an analytic subset of $`W`$, $`G`$. For $`pWG`$ for all but a discrete set of values of $`t`$ we have $`rank(a_j^{(i)}(t,\epsilon )=k`$.
Without limiting generality, we assume $`U`$ to be a polydisc. For any fixed $`p_0WG`$, there is a natural mapping from the space of solutions of the system to the space of solutions of the derived equation, which we can identify with the linear mapping

$$x_1(t_0,p_0),\mathrm{},x_n(t_0,p_0)x_1(t_0,p_0),\mathrm{},x_1^{(k)}(t_0,p_0)$$

where $`t_0`$ is such that the leading coefficient does not vanish on $`t_0,p_0`$, given by

$$x_1^{(i)}(t_0,p_0)=(a_j^{(i)}(t_0,p_0))(x_j(t_0,p_0)).$$

Since the leading coefficient does not vanish at $`t_0,p_0`$, $`rank(a_j^{(i)}(t,\epsilon ))=k`$, and the mapping is surjective. Hence $`p_0WG`$ every solution of the derived equation comes from a solution of the system.
Remark: if the dependence is algebraic rather than holomorphic, the set $`G`$ is algebraic.
Let us express the determinant of the $`i`$-th $`k`$-minor, the order of the derived equation being $`k`$, as $`b_{i0}(p)+b_{i1}(p)t+\mathrm{}+b_{ik_i}t^{k_i}`$. We call the ideal generated by all $`b_{ij}(p)`$ the degeneracy ideal, and the corresponding variety the degeneracy variety.
We now make a small digression, and explain a basic difficulty which might occur in estimating the number of zeros of the first component of a solution of the system through the derived above equation for this component. In  one finds the following theorem, derived from growth estimates on the solutions:
Theorem. Let

$$a(t)y^{(n)}+a_{n1}(t)y^{(n1)}+\mathrm{}+a_0(t)y=0$$

have coefficients holomorphic in a region $`U`$ containing a real segment $`K`$. Suppose that the modulus of all $`a_i(t)`$ is bounded in $`U`$ by $`A`$, and that

$$max_{tK}|a(t)|a.$$

Then the number of zeros on $`K`$ of any holomorphic in $`U`$ solution of the equation, real on the real axis, is not greater than

$$(A/a+n)^\mu ,$$

where $`\mu `$ depends only on the geometry of $`U`$ and $`K`$.
The main obstacle in applying this theorem to our case lies in the fact that, as parameters vary, $`max_{tK}|a(t,p)|`$ (which depends on the parameters) a priori may become arbitrarily small - the situation of singular perturbations - and the estimates explode.
Our next goal is to show that the scalar differential equation we derived for $`x_1`$ cannot be in fact singularly perturbed, at least when one parameter is involved. This is somewhat unexpected, and we explain why. Consider the following very simple $`\epsilon `$-dependent system:

$$\dot{x}=x+\epsilon y\dot{y}=x+y$$

One might expect that when $`\epsilon =0`$, the derived equation, which is of second order, will simply degenerate into the first order $`\dot{x}=x`$, implying that the derived equation is singularly perturbed. However, $`\ddot{x}=x+\epsilon y+\epsilon (x+y)=(1+\epsilon )x+2\epsilon y`$, and the derived equation is

$$\ddot{x}2\dot{x}+(1\epsilon )x=0.$$

This equation is not singularly perturbed. It has the property that for the exceptional value of the parameter, an additional solution appears which does not come from the original system, namely $`te^t`$. To explain why the singular perturbations cannot occur also in the general case, we proceed as follows. Until further notice the dimension of the parameter space is one, and the parameter will be denoted by $`\epsilon `$.
Definition. We call a meromorphic function $`f(t,\epsilon )`$ singularly perturbed at $`\epsilon _0`$, if when expressed as a quotient of two holomorphic functions $`f(t,\epsilon )=u(t,\epsilon )/v(t,\epsilon )`$, $`u(t,\epsilon )`$ is divisible by a higher or equal power of $`(\epsilon \epsilon _0)`$ than the one that divides $`v(t,\epsilon )`$.
Definition. A (scalar) differential equation with meromorphic coefficients depending on one parameter, normalized so that the leading coefficient is equal to 1, will be called singularly perturbed if at least one of its coefficients is singularly perturbed at some point in the parameter space.
We make the following obvious claim.
Lemma 1. Let us be given a differential equation with $`t,\epsilon `$-dependent coefficients meromorphic in a domain $`U`$, singularly perturbed at $`\epsilon _0`$. Then there is a point $`t_0`$, such that in the neighborhood of $`(t_0,\epsilon _0)U`$, the equation takes the form

$$(\epsilon \epsilon _0)^sy^{(n)}+a_{n1}(t,\epsilon )y^{(n1)}+\mathrm{}+a_0(t,\epsilon )y=0,$$

$`s>0`$, with coefficients holomorphic in the neighborhood of $`(t_0,\epsilon _0)`$, and at least one of them is not identically zero when $`\epsilon =\epsilon _0`$.
Theorem 1. Let a differential equation of order $`n`$ with coefficients which depend meromorphically on $`t,\epsilon `$ in some polydisc $`U\times V`$, be singularly perturbed for at least one parameter value $`\epsilon _0`$, and let it have l solutions $`y_1(t,\epsilon ),\mathrm{},y_l(t,\epsilon )`$ holomorphic in $`t`$ and $`\epsilon `$ in that polydisc. Then for each fixed parameter value $`\widehat{\epsilon }`$ from $`V`$,

$$dim_{}span_{}(y_1(t,\widehat{\epsilon }),\mathrm{},y_l(t,\widehat{\epsilon }))n1.$$

Proof. Using Lemma 1, we write the equation in the form given there in some smaller polydisc, and assume $`t_0=0,\epsilon _0=0`$:

$$\epsilon ^sy^{(n)}+a_{n1}(t,\epsilon )y^{(n1)}+\mathrm{}+a_0(t,\epsilon )y=0,$$

with $`a_k(0,0)0`$. We claim that if two holomorphic in $`t,\epsilon `$ solutions are such that they have equal $`y^{(n1)}(0,\epsilon ),\mathrm{},y^{(k+1)}(0,\epsilon ),y^{(k1)}(0,\epsilon ),\mathrm{},y(0,\epsilon )`$, they are in fact identical. Indeed, otherwise we get that there exists a nonzero holomorphic in $`t,\epsilon `$ solution $`y(t,\epsilon )`$ such that

$$y^{(n1)}(0,\epsilon )=\mathrm{}=y^{(k+1)}(0,\epsilon )=y^{(k1)}(0,\epsilon )=\mathrm{}=y(0,\epsilon )=0.$$

Since $`y(t,\epsilon )`$ is not identically zero, $`y^{(k)}(0,\epsilon )`$ cannot be identically zero. Thus $`y^{(k)}(0,\epsilon )=\epsilon ^hg(\epsilon )`$, $`h0,g(0)0`$. Assume $`h>0`$. $`y(t,0)`$ is the solution of the differential equation with $`\epsilon `$ set to zero. Since $`a_k(0,0)0`$ and the initial conditions at $`\epsilon =0`$ are zero ($`h>0`$), we have $`y(t,0)=0`$ identically. This can only happen if $`y(t,\epsilon )=\epsilon y_1(t,\epsilon )`$, $`y_1(t,\epsilon )`$ being also a holomorphic solution of the system with the same zero initial conditions except for $`y_1^{(k)}(0,\epsilon )=\epsilon ^{h1}g(\epsilon )`$. We continue in the same way, until we get a holomorphic solution $`y_h(t,\epsilon )=y(t,\epsilon )/\epsilon ^h`$ with the initial conditions

$$y^{(n1)}(0,\epsilon )=\mathrm{}=y^{(k+1)}(0,\epsilon )=y^{(k1)}(0,\epsilon )=\mathrm{}=y(0,\epsilon )=0,$$

and $`y^{(k)}(o,\epsilon )=g(\epsilon )`$, $`g(0)0`$. However substituting $`y_h(t,\epsilon )`$ into the original equation and setting $`t=0`$, we get:

$$\epsilon ^sy^{(n)}(0,\epsilon )+a_k(0,\epsilon )y^{(k)}(0,\epsilon )=0,$$

which is a contradiction ($`a_k(0,0),y^{(k)}(0,0)0`$).
Thus in the smaller polydisc the holomorphic in $`t,\epsilon `$ solutions can be identified with sets of holomorphically dependent on $`\epsilon `$ sets of initial conditions $`y^{(n1)}(0,\epsilon ),\mathrm{},y^{(k+1)}(0,\epsilon ),y^{(k1)}(0,\epsilon ),\mathrm{},y(0,\epsilon )`$, omitting $`y^{(k)}(0,\epsilon )`$. Now we prove the theorem. Suppose we have holomorphic in $`U\times V`$ solutions of our singularly perturbed equation $`y_1(t,\epsilon ),..,y_l(t,\epsilon )`$. Over the field meromorphic in $`V`$ functions there can be at most $`n1`$ linearly independent such solutions since the latter are identified with $`𝒪(V)^{n1}`$. So we choose a maximal set $`S`$ of linearly independent over $`(V)`$ solutions and express all the other as their linear combination over $`(V)`$. It follows that all the other solutions lie in the space spanned by the solutions from $`S`$ for all but a discrete set of $`\epsilon `$-s. Thus for almost all values of $`\epsilon `$, the space of holomorphic in $`t,\epsilon `$ solutions is $`n1`$ dimensional, which implies that it is at most $`n1`$ dimensional for all values of $`\epsilon `$ by the semicontinuity of the dimension.
Theorem 2. Suppose that we are given a linear differential system with coefficients depending on one parameter, holomorphic both in $`t`$ and $`\epsilon `$ in a polydisc $`U\times V`$. Then the derived scalar differential equation for the variable $`x_1`$ cannot be singularly perturbed.
Proof. Returning to notations of proof of Claim 1, we take a point $`(t_0,\epsilon _0)UZ`$, in which the leading coefficient of the derived equation is not zero, and $`rank(a_j^{(i)}(t,\epsilon )=k`$.
For each $`m`$, there is a holomorphic in $`U`$ solution of the system defined by the initial conditions (constant in $`\epsilon `$):

$$x_1(t_0,\epsilon )=0\mathrm{}x_m(t_0,\epsilon )=1\mathrm{}x_n(t_0,\epsilon )=0.$$

We denote the 1-st component of this solution by $`x_1^m(t,\epsilon )`$. By discussion from the proof of Claim 1, since $`x_1^1(t,\epsilon _0),\mathrm{},x_1^n(t,\epsilon _0)`$ span the space of solutions of the derived equation at $`\epsilon _0`$, and since at this (generic) point the order of the derived equation is $`k`$, we have

$$dim_{}span_{}x_1^1(t,\epsilon _0),\mathrm{},x_1^n(t,\epsilon _0)=k.$$

However, if the derived equation is singularly perturbed even at only one point $`\epsilon _1V`$, from Theorem 1 it would follow that for every $`\epsilon _0V`$

$$dim_{}span_{}x_1^1(t,\epsilon _0),\mathrm{},x_1^n(t,\epsilon _0)k1,$$

giving a contradiction. Thus the derived equation cannot be singularly perturbed.
We explain now our goal in what follows. As it was already mentioned, we can parameterize the problem of general linear system of order $`n`$ with polynomial coefficients of degree $`d`$ as:

$$\dot{x}_i=\underset{j=1}{\overset{n}{}}\underset{k=0}{\overset{d}{}}a_{ijk}t^kx_ji=1,..,n$$

the parameters being $`a_{ijk}`$, thus making the system matrix to be polynomial in the time and parameters, with integer coefficients. Ideally we would like to bound the number of zeros which the first variable of this system might have, but the results available here are still not permitting such general bounds. Instead, we restrict ourselves to the case when the system matrix is polynomial in time and one additional parameter, then making several remarks about the general linear differential system with polynomial coefficients as above.
So let us have a linear differential system of order $`n`$, with coefficients being polynomials in $`t,\epsilon `$ of degree $`d`$, with integer coefficients of maximum modulus $`M`$. Let the order of the derived equation be $`k`$.
Lemma 2.  Let a nonsingularly perturbed rational function be written as

$$\frac{a_0(\epsilon )+a_1(\epsilon )t+\mathrm{}+a_k(\epsilon )t^k}{b_0(\epsilon )+b_1(\epsilon )t+\mathrm{}+b_l(\epsilon )t^l}.$$

Then each $`a_j`$ belongs to the ideal generated by $`b_0,b_1,..,b_l`$.
Proof. $`a_j`$ must vanish at points where all $`b_1,..,b_l`$ vanish (otherwise the fraction would be singularly perturbed). Write $`b_i(\epsilon )=c_i(\epsilon )b(\epsilon )`$, where $`b(\epsilon )`$ is the common factor. Then $`c_i(\epsilon )`$ have no common zero and therefore, since we are working over $``$, $`1=h_1(\epsilon )c_1(\epsilon )+..+h_l(\epsilon )c_l(\epsilon )`$, for some polynomials $`h_1,..,h_l`$. But $`a_j`$ must be divisible by $`b(\epsilon )`$: otherwise we get singular perturbations. Thus:

$$a_j(\epsilon )=d(\epsilon )b(\epsilon )=d(\epsilon )b(\epsilon )(h_1(\epsilon )c_1(\epsilon )+..+h_l(\epsilon )c_l(\epsilon ))=d(\epsilon )h_1(\epsilon )b_1(\epsilon )+..+d(\epsilon )h_l(\epsilon )b_l(\epsilon ).$$

We derive now some computational estimates we will need.
Lemma 3. Suppose we are given a non-singularly perturbed rational function of $`t,\epsilon `$,

$$\frac{a_0(\epsilon )+a_1(\epsilon )t+\mathrm{}+a_k(\epsilon )t^k}{b_0(\epsilon )+b_1(\epsilon )t+\mathrm{}+b_l(\epsilon )t^l},$$

with degrees of the numerator and the denominator being at most $`d`$, and the coefficients being integers of maximum modulus M (we mean by this possibly complex integers). Then one may write:

$$\frac{a_0(\epsilon )+a_1(\epsilon )t+\mathrm{}+a_k(\epsilon )t^k}{b_0(\epsilon )+b_1(\epsilon )t+\mathrm{}+b_l(\epsilon )t^l}=\frac{\underset{i=1}{\overset{k}{}}t^i\underset{j=1}{\overset{l}{}}h_{ij}(\epsilon )b_j(\epsilon )}{b_0(\epsilon )+b_1(\epsilon )t+\mathrm{}+b_l(\epsilon )t^l},$$

with $`h_{ij}(\epsilon )`$ being polynomials of degree at most $`2d1`$ with integer coefficients of maximum modulus $`(eM)^{Cd^3}`$ for some constant $`C`$.
Proof. From Lemma 2, $`a_i=h_{ij}b_j`$. By , since $`ideg(a_i)d`$, $`jdeg(b_j)d`$, we have that $`a_i=h_{ij}b_j,deg(h_{ij})2d1`$. Observe that this gives us an integer coefficient linear system of dimension at most $`2d(d+1)`$, maximum modulus of the coefficients not greater than $`M`$. By the integrality trick this gives us the estimate $`(2d(d+1))!M^{2d(d+1)}(eM)^{Cd^3}`$ on the coefficients of the polynomials $`h_{ij}`$.
Lemma 4. In the setting of Lemma 3, suppose that $`|\epsilon |`$ is bounded by $`E`$. Then the above fraction can always be written as

$$\frac{\underset{i=1}{\overset{d}{}}c_it^i}{_{i=1}^dd_it^i}$$

with $`c_i(eM)^{Cd^3}(E^{2d1}+1)`$ for some constant $`C`$ , with some $`d_i`$ equal to $`1`$ and all $`|d_i|`$ being not greater than 1.
Proof. For each $`\epsilon `$ in the disc, take the $`b_i`$ having the maximum modulus, and divide by it both the numerator and the denominator: use Lemma 3 (the constant here is different from the one in Lemma 3).
Lemma 5. Suppose we are given a non-singularly perturbed differential equation of order $`k`$

$$y^{(k)}+\frac{a_{k1}(t,\epsilon )}{b(t,\epsilon )}y^{(k1)}+\mathrm{}+\frac{a_0(t,\epsilon )}{b(t,\epsilon )}y=0,$$

all $`a_i(t,\epsilon )`$ and $`b(t,\epsilon )`$ being polynomials of degree at most $`d`$, having integer coefficients of maximum modulus $`M`$.
Let $`[R/2,R/2]`$ be a segment on the real axis containing $`[1,1]`$. Fix any real $`\epsilon `$ of modulus less than $`E`$. Let $`f`$ be a holomorphic solution of the equation for this parameter value, real on the real axis. Then the number of zeros $`f`$ can have on $`[R/2,R/2]`$, is bounded by

$$\left((Me)^{Cd^3}(E^{2d}+1)R^{d+1}+k\right)^\sigma ,$$

with $`\sigma `$ and C being some constants.
Proof. By Lemma 4, we can write the equation in the form

$$(t^s+\underset{is}{}d_it^i)y^{(k)}+(c_{k\mathrm{1\; 0}}+..+c_{k1d}t^d)y^{(k1)}+\mathrm{}+(c_{00}+\mathrm{}+c_{0d}t^d)y=0.$$

with $`|c_{ij}|`$ being bounded by the bound from Lemma 2, and all $`|d_i|1`$. If we could find now some point in $`[R/2,R/2]`$ where $`t^s+_{is}d_it^i`$ has some value bounded away from zero, we are done. This because then we may just apply the Theorem cited from .
So we write

$$|t^s+\underset{is}{}d_it^i||t^s+d_{s1}t^{s1}+\mathrm{}+d_0||t^{s+1}(d_{s+1}+d_{s+2}t+..)|$$

and notice that by Cartan lemma (), one may delete disks of total diameter $`h`$ from $``$, such that on the rest of $``$ the modulus of the monic polynomial $`t^s+d_{s1}t^{s1}+\mathrm{}+d_0`$ is $`(h/4e)^s`$. Thus for each $`h>0`$ there exists a point $`t_h`$ of modulus at most $`h/2`$, at which the value of the polynomial is $`(h/4e)^s`$. But then:

$$|t_h^s+\underset{is}{}d_it_h^i||t_h^s+d_{s1}t_h^{s1}+\mathrm{}+d_0||t_h^{s+1}(d_{s+1}+d_{s+2}t_h+..)|(h/4e)^sh^{s+1}(ds).$$

Taking for example $`h=1/(2(4e)^s(ds))`$ we get that there exists a point on the unit segment at which $`|t^s+_{is}d_it^i|`$ is bounded from below by

$$\frac{1}{(4e)^{s^2+s}2^{s+1}(ds)^s}\frac{1}{(4e)^{s^2+s}2^{s+1}d^s}$$

As $`s`$ can be from $`0`$ to $`d1`$ (it can be shown that the case $`s=d`$ is also covered by the obtained estimate), we are led to the conclusion that there will be a point in $`[R/2,R/2]`$ where

$$|t^s+\underset{is}{}d_it^i|\frac{1}{(4e)^{d^2+d}2^{d+1}d^d}.$$

In turn, $`|_{j=0}^dc_{ij}t^j|`$ cannot, when $`\epsilon `$ varies over disk of radius $`E`$ and $`t`$ over $`D_R`$, be larger than

$$\underset{j=0}{\overset{d}{}}R^j(eM)^{Cd^3}(E^{2d1}+1)(eM)^{Cd^3}(E^{2d1}+1)R^{d+1}.$$

By Theorem cited from , one deduces that for any given real $`|\epsilon |E`$, any holomorphic in $`t`$ solution of the given differential equation, real on the real axis, may have not more than

$$\left(\frac{(eM)^{Cd^3}(E^{2d1}+1)R^{d+1}}{1/(4e)^{d^2+d}2^{d+1}d^d}+k\right)^\sigma $$

zeros on $`[R/2,R/2]`$, $`\sigma `$ not depending on $`R`$ because of properties of $`\sigma (K,U)`$ from the Theorem. Simplifying, one gets the bound we stated.
To estimate the number of zeros a solution (more precisely, its $`x_1`$ component) of the initial linear system

$$\dot{x}=A(t,\epsilon )x$$

might have, we now have only to bound the maximum modulus of the (integer) coefficients of the derived equation. We will do the computation for an arbitrary number of parameters; to emphasize this we shall write $`p`$ instead of $`\epsilon `$. To remind the notations, we denote the dimension of the system by $`n`$, the degree of the polynomial coefficients by $`d`$, the maximum modulus of the coefficients of those polynomials by $`M`$, the number of parameters by $`q`$, and the order of the derived equation for $`x_1`$ by $`k`$.
Lemma 6. The following recurrent formulas hold for $`a^{(i)}`$, defined in the beginning of this text:

$$a^{(i)}=\dot{a}^{(i1)}+A(t,p)^Ta^{(i1)}.$$

Proof. Immediate.
The following is a computation:
Lemma 7. Let $`P_1`$, $`P_2`$ be two polynomials of degrees $`d_1,d_2`$ in s variables, with the modulus of their coefficients being not larger than $`M_1,M_2`$, respectively. Then the degree of $`P_1P_2`$ is $`d_1+d_2`$ and the coefficients are not larger in absolute value than $`(1+min(d_1,d_2))^sM_1M_2`$.
Lemma 8. The degrees and the maximum modulus of coefficients of the entries of $`a^{(i)}`$ are $`di`$ and $`n^i(d+(d+1)^{q+1}M)^i`$, respectively.
Proof. By writing recursive relations.
We now compute the maximum modulus and the degrees for the coefficients of the derived equation.
Lemma 9. Let

$$y^{(k)}+\frac{\alpha _{k1}(t,p)}{\beta (t,p)}y^{(k1)}+\mathrm{}+\frac{\alpha _0(t,p)}{\beta (t,p)}y=0$$

be the equation derived for $`x_1`$ from the original system. Then the degrees of both $`\alpha _i`$ and $`\beta `$ are bounded by $`k(k+1)d/2`$, and the maximum modulus for the (integer) coefficients of those polynomials is not larger than $`(eM(d+1))^{C(q+1)n^3}`$ for some constant $`C`$.
Proof. Again by recursive relations.
Theorem 2. Let the initial linear differential system have as coefficients polynomials in $`t,\epsilon `$ ($`\epsilon `$ one dimensional) with integer coefficients of maximum modulus $`M`$, and let $`[R/2,R/2]`$ be a segment on the real axis containing $`[1,1]`$. Fix any real $`\epsilon `$ of modulus less than $`E`$. Let $`f`$ be a holomorphic solution of the equation for this parameter value, real on the real axis. Then the number of zeros $`f`$ can have on $`[R/2,R/2]`$, is bounded by

$$\left((Me)^{Cn^9d^4}(E^{n(n+1)d}+1)R^{n(n+1)d/2+1}+n\right)^\sigma ,$$

for some constants $`C`$ and $`\sigma `$.
Proof. Just substitute the bounds stated in Lemma 9 with $`q=1`$ into Lemma 5 and simplify the expression.
Of course, the concrete form of the bound given above is irrelevant; what is important is the rate of its growth, which is simple exponential in $`n,d`$ compared with the tower of exponents which can be derived from  in this situation.
We now make some remarks on the general case when many parameters are allowed. First, it is clear that the absence of singular perturbations in the case of dependence on one parameter implies absence of singular perturbations along any given holomorphic curve in the parameter space, and the following is true.
Lemma 2a.  Let a rational function dependent on several parameters be written as

$$\frac{a_0(p)+a_1(p)t+\mathrm{}+a_k(p)t^k}{b_0(p)+b_1(p)t+\mathrm{}+b_l(p)t^l},$$

and let it be such that it is not singularly perturbed along any line in the parameter space. Then each $`a_j`$ vanishes on the variety where all $`b_0,b_1,..,b_l`$ vanish.
We cannot in general have the conclusion of Lemma 2 as the example of the nonsingularly perturbed fraction

$$\frac{ab}{a^2+b^2t}$$

shows. That is, in the derived equation for system with polynomial coefficients depending on many parameters we can only guarantee that the coefficients in the numerator vanish on the degeneracy variety, but not that they in fact belong to the degeneracy ideal. However, if by some additional reasoning for, say, the general linear system with polynomial coefficients

$$\dot{x}_i=\underset{j=1}{\overset{n}{}}\underset{k=0}{\overset{d}{}}a_{ijk}t^kx_ji=1,..,n,$$

it will turn out to be true, then the bound, though drastically grows because of the double exponential in $`n^2(d+1)`$ bounds in the effective Nullstellensatz (see ), has the asymptotic form (in terms of Theorem 2)

$$e^{e^{e^{P(n,d)}}}(ER)^{Cn^2d},$$

where $`P`$ is a polynomial. This is still a considerable improvement compared to  (of course if the assumption is valid).

References.

1. Novikov, D.; Yakovenko, S. Trajectories of polynomial vector fields and ascending chains of polynomial ideals. Ann. Inst. Fourier (Grenoble) 49 (1999), no. 2, 563–609.
2. Il’yashenko, Y.; Yakovenko, S. Counting real zeros of analytic functions satisfying linear ordinary differential equations. J. Differential Equations 126 (1996), no. 1, 87–105.
3. Yakovenko S. On functions and curves defined by ordinary differential equations , to appear in Proceedings of the Arnoldfest (Ed. by Bierstone, Khesin, Khovanskii, Marsden), Fields Institute Communications, 1998.
4. Il’yashenko, Y.; Yakovenko, S. Double exponential estimate for the number of zeros of complete abelian integrals and rational envelopes of linear ordinary differential equations with an irreducible monodromy group. Invent. Math. 121 (1995), no. 3, 613–650.
5. Shiffman, B. Degree bounds for the division problem in polynomial ideals. Michigan Math. J. 36 (1989), no. 2, 163–171.


format: text

text:
# Preheating of massive fermions after inflation: analytical results

## 1 Introduction

One of the key ingredients for our understanding of the early Universe is the mechanism of inflation , which constitutes a very elegant solution to several cosmological problems.

Despite of the simplicity of the general idea, the details of the physical processes which govern it are still somehow unclear and matter of intense work. The two main aims of these studies are

* to embed inflation in a context more motivated by particle physics (for a review see ) and
* to understand the reheating phase which converts the energy density that drives inflation into the matter and radiation that we presently see.

This second issue has been deeply influenced in the last decade by the possibility of particle creation through parametric resonance . The application of this phenomenon to creation of matter after inflation has been called *preheating* in the paper , since (with the exception of some very recent versions ) it is usually followed by a stage of (ordinary) perturbative reheating.

Preheating of bosons  is characterized by a very efficient and explosive creation, due to the coherent effect of the oscillations of the inflaton field. This allows significant production even when single particle decay is kinematically forbidden. It has been very recently noticed  that preheating of fermions can also be very efficient despite the production is in this case limited by Pauli blocking.

Parametric creation of spin $`1/2`$ fermions has been the subject of some works in the past. Pure gravitational production has been examined in refs. , while creation by an oscillating background field is instead considered in the works . References  report results for creation in a Minkowski space. Reference  studies the production of massless fermions after a $`\lambda \varphi ^4`$ inflation, exploiting the fact that this case can also be reconducted to a static one. In this work, production in a static Universe after chaotic inflation is also considered, and some conjectures on the effects of the expansion are made. Moreover, the full calculation of preheating of massive fermions after chaotic inflation in an expanding Universe has been performed numerically in ref. .

These last works had a great impact on the most recent studies. For example, their results turned also useful to the study of gravitinos production at preheating . This issue is particularly important, since gravitinos can easily overclose the Universe (if they are stable) or (if they decay) spoil the successful predictions of primordial nucleosynthesis through photodissociation of the light elements. Gravitinos can be thermally produced during the stage of reheating. To avoid this overproduction, the reheating temperature $`T_{\mathrm{RH}}`$ after inflation cannot be larger than $`(10^810^9)`$ GeV . However, it has been realized  that the non-thermal production of helicity $`\pm 1/2`$ gravitinos (whose equation of motion can be reconducted to the one of an ordinary spin $`1/2`$ Dirac particle) can easily be more efficient than the thermal one, and this in general leads to more stringent upper bounds on $`T_{\mathrm{RH}}`$. Several papers related to the works  have recently appeared .

Another important implication of preheating of fermions is constituted by leptogenesis, as the work  and the related papers  show. In this scheme , a leptonic asymmetry is first created from the decay of right-handed neutrinos, and then partially converted to baryon asymmetry through sphaleronic interactions. Since leptogenesis is very sensitive to both the mechanism of creation of the heavy neutrinos and to the neutrino mass matrices, it could constitute an interesting link between preheating and the experiments on neutrino oscillations.

Other phenomenological implications of these works appear in refs. , with preheating as a possible mechanism for creating superheavy relic particles responsible for the ultrahigh energy cosmic rays, and in ref. , where the possible impact of fermions produced *during* inflation on the microwave background anisotropies and on the large-structure surveys is considered. Finally, fermionic production can play an interesting role in hybrid inflationary models .

Due to the large number of these studies, it may be worth to reconsider the basic mechanism of fermionic preheating. In particular, it should be important to give an analytical confirmation to the results of production of massive fermions in the expanding Universe, which are so far known only numerically from the analysis . This is the aim of the present work.

In the next section we revise the basic formalism for preheating of fermions. We consider creation of very massive particles right after chaotic inflation. The coupling of the fermions to the oscillating inflaton gives them a time varying mass. As it is known, this can cause a non-adiabatical change of the frequency of the fermions and their consequent creation. In case of very massive fermions, the non-adiabaticity condition can be satisfied only when their total mass vanishes, and the production occurs at discrete intervals, until the inflaton oscillations become too small for the total fermionic mass to vanish.

In section 3 we derive analytical formulae for the spectra of the fermions after a generic production. Our derivation follows the one developed in ref.  for preheating of bosons. It exploits the fact that the production occurs in very short intervals around the zeros of the total fermionic mass: the calculation is made possible from the fact that the occupation number can be considered as constant outside these small regions, and that the expansion of the Universe can be neglected inside them. As a result, the only physical quantities relevant for the creation are the time derivative $`\varphi ^{}`$ of the inflaton field and the value of the scale factor $`a`$ at each production. The derivation of the “fermionic counterpart” of the formulae obtained for preheating of bosons in ref.  has also been done to a certain extent in the work , where the results of a single production during inflation is given. However, when one is interested in the successive productions, a more detailed study is necessary, as our analysis shows. As may be expected, the final results that we obtain closely resemble the ones of the work . What is most surprising is their excellent agreement with the numerical results, as some figures provided manifest.

In section 4 we consider the production in a non-expanding Universe. In this case our analytical formulae considerably simplify and agree with the ones of refs. . In particular, they show the presence of resonance bands which are anyhow limited by Pauli blocking.

In section 5 we study the more interesting case of production in an expanding Universe. As indicated in ref.  and as confirmed by the numerical results of the work , the creation is now very different with respect to the previous case. The expansion removes the resonance bands and the production (almost) saturates a Fermi sphere up to a maximal momentum. Our analytical results confirm this behavior.

In section 5 we also calculate the total energy density $`\rho _X`$ of produced fermions, which may be the quantity of most physical relevance. We compare our results with the ones of ref. , where it is shown that the final value of $`\rho _X`$ (normalized to the inflaton energy density) scales linearly with the parameter $`q(g^2\varphi _0^2)/(4m_\varphi ^2)`$,111In this expression $`g`$ is the Yukawa coupling between the inflaton and the fermions, $`\varphi _0`$ the initial value of the inflaton, and $`m_\varphi `$ its mass. while it depends very weakly on the fermion bare mass $`m_X`$. However, these numerical results are valid only in a limited range for $`m_X`$, and it has been wondered if the density $`\rho _X`$ decreases at values of $`m_X`$ below this range. Our analytical results give a positive answer in this regard. To see this, a proper average of the analytical formulae must be done, exploiting the fact that the expansion of the Universe gives the production a stochastic character. In this way one can get a “mean” function that interpolates very well between the maxima and the minima of the spectra of produced particles. Again we derived it in close analogy with what is done in the bosonic case, and again the results that we get are in very good agreement with the numerical ones of ref.  in the region of validity of the latter.

All this analysis neglects the backreaction of the produced fermions on the evolution of the inflaton field and of the scale factor. Despite the difficulty of a more complete treatment, backreaction effects can be understood at least in the Hartree approximation. This was done numerically in ref. . In section 6 we see that the analytical formulae here provided allow to understand the effects of backreaction observed in the numerical simulations.

## 2 Production of fermions at preheating

In this section we revise the basic formalism for production of fermions at preheating. Our presentation follows the ones of refs.  with the correction of few typos.

We start from the Dirac equation (in conformal time $`\eta `$) for a fermionic field $`X`$ in the Friedmann-Robertson-Walker background:

$$\left(\frac{i}{a}\gamma ^\mu _\mu +i\frac{3}{2}H\gamma ^0m\right)X=0,$$
(1)
where $`a`$ is the scale factor of the Universe, $`H=a^{}/a^2`$ the Hubble rate,222Here and in the following prime denotes derivative which respect to $`\eta `$. and $`m`$ the total mass of the fermion.

Fermionic production is possible if this last quantity varies non adiabatically with time. This may happen during the coherent oscillations of the inflaton $`\varphi `$ at reheating in presence of a Yukawa interaction $`\varphi \overline{X}X`$, such that the total mass is given by

$$m(\eta )=m_X+g\varphi (\eta ).$$
(2)
We will always consider fermions with very high bare mass $`m_X`$. In this case, non adiabaticity can be achieved if the Yukawa interaction is sufficiently strong to make the total mass (2) vanish. Fermions are then created whenever the inflaton field crosses the value $`\varphi _{}m_X/g`$.

We will study the production after chaotic inflation, that is while the inflaton field coherently oscillates about the minimum of the potential

$$V=\frac{1}{2}m_\varphi ^2\varphi ^2,m_\varphi 10^{13}\mathrm{GeV}.$$
(3)
If one neglects the backreaction of the created particles (this effect will be considered in the last part of the work), then, after few oscillations, the inflaton evolves according to the formula

$$\varphi (t)\frac{M_p}{\sqrt{3\pi }}\frac{\mathrm{cos}\left(m_\varphi t\right)}{m_\varphi t},$$
(4)
where $`t`$ is the physical time. The presence of the $`t`$ at the denominator of eq. (4) shows the damp of the oscillations due to the expansion of the Universe. It thus follows that there exists a final time after which $`|\varphi |<m_X/g`$ and the total mass no longer vanishes, so that the production ends.

To proceed with our analysis, we redefine $`\chi =a^{3/2}X`$. Equation (1) acquires the more familiar form

$$(i/am)\chi =0.$$
(5)
We decompose

$$\chi (x)=\frac{d^3k}{(2\pi )^{3/2}}e^{i𝐤𝐱}\underset{r}{}\left[u_r(k,\eta )a_r(k)+v_r(k,\eta )b_r^{}(k)\right],$$
(6)
with333We choose    $`\gamma ^0=\left(\begin{array}{cc}\mathrm{𝟙}& 0\\ 0& \mathrm{𝟙}\end{array}\right),\gamma ^1=\left(\begin{array}{cc}0& i\sigma _2\\ i\sigma _2& 0\end{array}\right),\gamma ^2=\left(\begin{array}{cc}0& i\sigma _1\\ i\sigma _1& 0\end{array}\right),`$    $`\gamma ^3=\left(\begin{array}{cc}0& \mathrm{𝟙}\\ \mathrm{𝟙}& 0\end{array}\right)C=\left(\begin{array}{cc}0& \sigma _1\\ \sigma _1& 0\end{array}\right).`$

$$v_r(k)=C\overline{u}_r^T\left(k\right).$$
(7)
It is convenient to take the momentum along the third direction $`kk_z`$ and to define

$$u_r[\frac{u_+(\eta )}{\sqrt{2}}\psi _r,\frac{u_{}(\eta )}{\sqrt{2}}\psi _r]^T,v_r[\frac{v_+\left(\eta \right)}{\sqrt{2}}\psi _r,\frac{v_{}\left(\eta \right)}{\sqrt{2}}\psi _r]^T$$
(8)
with $`\psi _+=`$ $`\left(\begin{array}{c}1\\ 0\end{array}\right)`$ and $`\psi _{}=`$ $`\left(\begin{array}{c}0\\ 1\end{array}\right)`$ eigenvectors of the helicity operator $`\sigma 𝐯/|𝐯|`$. From the normalization adopted in eq. (8), one can impose

$`u_r^+u_s=v_r^+v_s=\delta _{rs},u_r^+v_s=0,|u_+|^2+|u_{}|^2=2.`$ (9)
With this choice, the Dirac equation (5) rewrites

$$u_\pm ^{}(\eta )=iku_{}(\eta )iamu_\pm (\eta ),$$
(10)
which can be decoupled into444We do not deal with the analogous equations for $`v_\pm `$, since from condition (7) it simply follows $`v_+(\eta )=u_{}^{}(\eta ),v_{}(\eta )=u_+^{}(\eta )`$.

$`u_\pm ^{\prime \prime }+\left[\omega _k^2\pm i(am)^{}\right]u_\pm =0,\omega ^2(\eta )=k^2+a^2m^2.`$ (11)
Under this decomposition, the hamiltonian of the system

$$H=\frac{1}{a}d^3𝐱\chi ^+(x)\left(i_\eta \right)\chi (x)$$
(12)
rewrites

$`H`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle }d^3k{\displaystyle \underset{r}{}}\{E_k\left(\eta \right)[a_r^{}\left(k\right)a_r\left(k\right)b_r\left(k\right)b_r^{}\left(k\right)]+`$ (13)
$`+F_k\left(\eta \right)b_r(k)a_r\left(k\right)+F_k^{}\left(\eta \right)a_r^{}\left(k\right)b_r^{}(k)\},`$
where (using the equations of motion)

$`E_k`$ $`=`$ $`kRe\left(u_+^{}u_{}\right)+am\left(1u_+^{}u_+\right),`$
$`F_k`$ $`=`$ $`{\displaystyle \frac{k}{2}}\left(u_+^2u_{}^2\right)+amu_+u_{},E_k^2+|F_k|^2=\omega ^2.`$ (14)
It is always possible to choose an initial configuration with the hamiltonian in the standard (diagonal) form. If we take, for example,

$$u_\pm \left(\eta _0\right)=\left(1\frac{am}{\omega }\right)^{1/2}e^{i\varphi },\varphi \mathrm{arbitrary}\mathrm{phase},$$
(15)
we have $`E_k\left(\eta _0\right)=\omega ,`$ $`F_k\left(\eta _0\right)=0`$.

However, the evolution equations (11) drive $`F_k`$ different from zero and a diagonal form for $`H`$ can be recovered only after a Bogolyubov transformation on the creation/annihilation operators

$`\widehat{a}(\eta ,k)`$ $``$ $`\alpha _k(\eta )a(k)+\beta _k(\eta )b^+(k),`$
$`\widehat{b}^+(\eta ,k)`$ $``$ $`\beta _k^{}(\eta )a(k)+\alpha _k^{}(\eta )b^+(k).`$ (16)
While it is immediate to see that the equal time anticommutation relations on both sets $`\{a^{(+)},b^{(+)}\}`$ and $`\{\widehat{a}^{(+)},\widehat{b}^{(+)}\}`$ enforce

$$|\alpha _k|^2+|\beta _k|^2=1,$$
(17)
the use of some algebra shows that a diagonalization of the hamiltonian is possible with the choice

$$\alpha =\beta \left(\frac{E+\omega }{F^{}}\right),|\beta |^2=\frac{|F|^2}{2\omega \left(E+\omega \right)}=\frac{\omega E}{2\omega }.$$
(18)
The time varying operators $`\widehat{a}^{(+)}`$ and $`\widehat{b}^{(+)}`$ are employed to define time dependent Fock spaces, each of them built from the zero (quasi)particle state at the time $`\eta `$:

$$\widehat{a}(\eta )|0_\eta =\widehat{b}(\eta )|0_\eta =0.$$
(19)
Let us consider the state $`|0_{\eta _0}|0`$ which describes a system with zero particles at the initial time $`\eta _0`$. This is no longer true for generic time $`\eta >\eta _0`$, since the occupation number at $`\eta `$ is defined in terms of the operators which diagonalize the hamiltonian at that moment. The particle density per physical volume $`V=a^3`$ at time $`\eta `$ is indeed given by555Of course the same result is achieved also for antifermions.

$`n(\eta )0|{\displaystyle \frac{N}{V}}|0=0|{\displaystyle \frac{1}{a^3}}{\displaystyle \underset{r=\pm 1}{}}{\displaystyle \frac{d^3𝐤}{\left(2\pi \right)^3}\widehat{a}_r^+(\eta ,r)\widehat{a}_r(\eta ,r)}|0=`$
$`={\displaystyle \frac{1}{\pi ^2a^3}}{\displaystyle 𝑑kk^2\left|\beta _k\right|^2},`$ (20)
which is different from zero whenever the hamiltonian is non diagonal in terms of the initial operators. The occupation number of created fermions is thus given by $`n_k=|\beta _k|^2`$ and the condition (17) ensures that the Pauli limit $`n_k<1`$ is respected.

## 3 Analytical evaluation of the occupation number

In this section we calculate analytically the evolution of the Bogolyubov coefficients (16) during the oscillations of the inflaton field after chaotic inflation.

In this derivation, we exploit the fact that, in the regime of very massive fermions we are interested in, the creation occurs only for very short intervals about the points $`\varphi _{}m_X/g`$ where the total fermionic mass (cfr. eq. (2)) vanishes. As the perfect agreement with the numerical results will confirm, this consideration allows one to treat the fermionic production with the same formalism adopted in the bosonic case . While far from the zeros of the total mass $`m`$ the Bogolyubov coefficients can be treated as constant, whenever $`\varphi `$ crosses $`\varphi _{}`$ a sudden variation occurs. Since the interval of production is very narrow, one can safely neglect the expansion of the Universe during the production and also linearize the function $`\varphi (\eta )\varphi _{}+\varphi ^{}\left(\eta _{}\right)\left(\eta \eta _{}\right)`$. As a consequence, the frequency $`\omega `$ defined in eq. (11) acquires the form

$$\omega ^2k^2+a^2\left(\eta _{}\right)\varphi ^2\left(\eta _{}\right)\left(\eta \eta _{}\right)^2,$$
(21)
and the whole calculation strongly resembles the one for scattering through a quadratic potential.

The use of this formalism is very well established in case of production of bosons . For what concerns fermions, it has been recently adopted by Chung et al.  for the study of the production during inflation. Fermionic production during inflation is possible only if the coefficient $`g`$ of the Yukawa interaction has opposite sign with respect to the inflaton field during inflation, or the total mass (2) would never vanish.666 Of course this constraint does not apply to our case, since the inflaton field changes sign after each half oscillation. If this is the case, it is possible to choose the value of $`g`$ such that the production occurs only once during inflation, while during reheating $`|\varphi |`$ is always too small for having creation. Our derivation is strongly inspired by this work. However, we are interested in couplings for which the production occurs several times during reheating, and it is not at all guaranteed a priori that an analytical approximation may work also in this case. The present analysis not only positively answers to this question, but also provides very simple and explicit formulae valid for arbitrary number of productions.

As may be expected, the final formulae are very similar to the ones found for production of bosons . What is most surprising is the almost perfect agreement between the analytical results and the numerical simulations, as we will show with some figures at the end of this section.

The first step for the derivation of the analytical formulae is to consider asymptotic solutions of eqs. (10) and (11). We look for solutions valid for $`\varphi `$ not very close to $`\varphi _{}`$, where the adiabaticity condition $`\omega ^{}\omega ^2`$ holds.

In this regime, the most general solutions of eqs. (10) and (11) are

$`u_+`$ $`=`$ $`A\left(1{\displaystyle \frac{am}{\omega }}\right)^{1/2}e^{i^\eta \omega _k𝑑\eta }+B\left(1+{\displaystyle \frac{am}{\omega }}\right)^{1/2}e^{i^\eta \omega _k𝑑\eta },`$
$`u_{}`$ $`=`$ $`B\left(1{\displaystyle \frac{am}{\omega }}\right)^{1/2}e^{i^\eta \omega _k𝑑\eta }+A\left(1+{\displaystyle \frac{am}{\omega }}\right)^{1/2}e^{i^\eta \omega _k𝑑\eta },`$ (22)
with $`|A|^2+|B|^2=1`$ following from the condition $`|u_+|^2+|u_{}|^2=2`$.

We can always put solutions of eqs. (10) and (11) into the form (22), with $`A`$ and $`B`$ in general functions of $`\eta `$. However, in most of the evolution (whenever the adiabaticity condition holds) it is a very good approximation to treat the coefficients $`A`$ and $`B`$ as constant.

Substituting the expressions (22) into eqs. (14), one finds

$`F=2\omega AB,E=\omega \left(|A|^2|B|^2\right),`$ (23)
from which it follows

$$|\beta |^2=\frac{|F|^2}{2\omega \left(E+\omega \right)}=|B|^2.$$
(24)
One can thus choose (up to an irrelevant global phase)

$$A\alpha B\beta ^{},$$
(25)
where $`\alpha `$ and $`\beta `$ are nothing but the Bogolyubov coefficients we are interested in.

Notice that the initial choice $`A\left(\eta _0\right)=1`$, $`B\left(\eta _0\right)=0`$ corresponds to the initial condition (15) and to the zero particles state chosen in the previous section.

As we have said, these coefficients undergo a sudden change whenever $`\varphi `$ crosses $`\varphi _0`$ and then they stabilize to new (almost) constant values. Our aim is to find the values at the end of the variation in terms of the ones prior to it.

We have not specified the lower limit of the integrals appearing in eqs. (22). For present convenience we choose it to be the time $`\eta _1`$ of the first production (that is when $`\varphi =\varphi _{}`$ for the first time).

Let us consider the evolution equation (11) near the point $`\eta _1`$. Since for high mass $`m_X`$ the fermionic production is limited to a very short interval, one can neglect the expansion of the Universe during it and write the equation for $`\varphi (\eta )`$ in a linearized form. We can thus write

$$am(\eta )a_1g\varphi _1^{}\left(\eta \eta _{\mathrm{\hspace{0.17em}1}}\right),a_1a\left(\eta _{\mathrm{\hspace{0.17em}1}}\right),\varphi _1\frac{d\varphi }{d\eta }|_{\eta _1}.$$
(26)
Following the notation of , we define

$$p\frac{k}{\sqrt{g|\varphi _1^{}|a_1}},\tau =\sqrt{g|\varphi _1^{}|a_1}\left(\eta \eta _1\right).$$
(27)
In terms of these new quantities, eqs. (11) rewrite (dot denotes derivative with respect to $`\tau `$)777For the production at the moment $`\eta _n`$, when the total mass vanishes for the n-th time, eq. (29) must be replaced by
$$\ddot{u}_\pm +\left(p^2\pm isign\left(\varphi _n^{}\right)+\tau ^2\right)u_\pm =0.$$
(28)  The effect of this replacement on the final results is reported below.

$$\ddot{u}_\pm +\left(p^2i+\tau ^2\right)u_\pm =0.$$
(29)
The point $`\eta _1`$ is thus mapped into the origin of $`\tau `$ and the region of asymptotic adiabaticity is at large $`|\tau |`$.

In the asymptotic solutions (22) we can see the behaviors

$`\left(1+{\displaystyle \frac{am}{\omega }}\right)^{1/2}`$ $``$ $`{\displaystyle \frac{p}{\sqrt{2}\tau }},`$
$`\left(1{\displaystyle \frac{am}{\omega }}\right)^{1/2}`$ $``$ $`\sqrt{2},`$
$`e^{\pm i^\eta \omega _k𝑑\eta }`$ $``$ $`\left({\displaystyle \frac{2\tau }{p}}\right)^{\pm ip^2/2}e^{\pm i\tau ^2/2}e^{\pm ip^2/4},`$ (30)
for $`\tau +\mathrm{}`$, and

$`\left(1+{\displaystyle \frac{am}{\omega }}\right)^{1/2}`$ $``$ $`\sqrt{2},`$
$`\left(1{\displaystyle \frac{am}{\omega }}\right)^{1/2}`$ $``$ $`{\displaystyle \frac{p}{\sqrt{2}\tau }},`$
$`e^{\pm i^\eta \omega _k𝑑\eta }`$ $``$ $`\left({\displaystyle \frac{p}{2\tau }}\right)^{\pm ip^2/2}e^{i\tau ^2/2}e^{ip^2/4},`$ (31)
for $`\tau \mathrm{}`$.

Equations (29) are solved  by parabolic cylinder functions $`D_\lambda (z)`$ . More precisely, the combination that matches the asymptotic solution (22) at $`\tau \mathrm{}`$ is888We deal only with the function $`u_+`$, since the study of $`u_{}`$ leads to the same results.

$`u_+(\tau )=A^{}\sqrt{2}\left({\displaystyle \frac{p}{\sqrt{2}}}\right)^{1+ip^2/2}e^{i(\frac{\pi }{4}\frac{p^2}{4})}e^{\pi p^2/8}D_{1ip^2/2}\left(\left(1+i\right)\tau \right)+`$
$`+B^{}\sqrt{2}\left({\displaystyle \frac{p}{\sqrt{2}}}\right)^{ip^2/2}e^{ip^2/4}e^{\pi p^2/8}D_{ip^2/2}((1i)\tau ).`$ (32)
In the above expression, $`A^{}`$ and $`B^{}`$ denote the values of the Bogolyubov coefficients before $`\varphi `$ crosses $`\varphi _{}`$, while the functions which multiply them are *exact* solutions of the *linearized* equation (29). The analytical approximation consists in considering them as solutions of the true evolution equation (11). For $`\tau +\mathrm{}`$ it is convenient to rewrite the solution (32) in terms of two different parabolic cylinder functions999This rewriting is always possible since eq. (29) has only two linearly independent solutions.

$`u_+(\tau )=[A^{}{\displaystyle \frac{\sqrt{\pi }pe^{\pi p^2/4}}{\mathrm{\Gamma }\left(1+ip^2/2\right)}}e^{i(\frac{\pi }{4}\frac{p^2}{2}+\frac{p^2}{2}\mathrm{ln}\frac{p^2}{2})}+B^{}e^{\pi p^2/2}]\times `$
$`\times \left\{\sqrt{2}\left({\displaystyle \frac{p}{\sqrt{2}}}\right)^{ip^2/2}e^{ip^2/4}e^{\pi p^2/8}D_{ip^2/2}((1i)\tau )\right\}+`$
$`+[A^{}e^{\pi p^2/2}+B^{}{\displaystyle \frac{\sqrt{\pi }pe^{\pi p^2/4}}{\mathrm{\Gamma }\left(1ip^2/2\right)}}e^{i(\frac{\pi }{4}\frac{p^2}{2}+\frac{p^2}{2}\mathrm{ln}\frac{p^2}{2})}]\times `$
$`\times \left\{\sqrt{2}\left({\displaystyle \frac{p}{\sqrt{2}}}\right)^{1+ip^2/2}e^{i(\frac{\pi }{4}\frac{p^2}{4})}e^{\pi p^2/8}D_{1ip^2/2}((1+i)\tau )\right\}.`$ (33)
In this new expression, the functions within curly brackets correspond to the asymptotic forms at $`\tau =+\mathrm{}`$ of the two terms of the solution (22). The coefficients in front of them give thus the new Bogolyubov coefficients in terms of the old ones.

All this derivation can be easily generalized when productions at successive zeros of the total mass $`m`$ are considered. The only important points are

* a difference in the values of the scale factor $`a`$ and of the derivative $`\varphi ^{}`$ at different $`\eta _i`$’s,
* a change of sign in the transfer matrix (the one which gives the new coefficients in terms of the old ones) whenever $`\varphi `$ crosses $`\varphi _{}`$ from below to above (cfr. the footnote just before eq. (29)), and
* the phase $`e^{\pm i{\scriptscriptstyle \omega _k𝑑\eta }}`$ which accumulates between $`\eta _1`$ and the $`\eta _i`$ considered.

Putting all this together, one has

$`\left(\begin{array}{c}\alpha _n\\ \beta _n^{}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}F_n& H_n\\ H_n^{}& F_n^{}\end{array}\right)\left(\begin{array}{c}\alpha _{n1}\\ \beta _{n1}^{}\end{array}\right)\mathrm{for}n\mathrm{odd},`$
$`H_n`$ $``$ $`H_n\mathrm{for}n\mathrm{even},`$ (34)
where $`a_n,\beta _n`$ are the values of the Bogolyubov coefficients after the $`n`$-th production, and where

$`F_n`$ $`=`$ $`\sqrt{1e^{\pi p_n^2}}e^{i(\frac{\pi }{4}\mathrm{arg}\mathrm{\Gamma }(ip_n^2/2)+\frac{p_n^2}{2}\mathrm{ln}(p_n^2/2)p_n^2/2)},`$
$`H_n`$ $`=`$ $`e^{\pi p_n^2/2}e^{2i_{\eta _1}^{\eta _n}\omega _k𝑑\eta },|F_n|^2+|H_n|^2=1.`$ (35)
We remind that $`p_n=k/\sqrt{g|\varphi _n^{}|a_n}`$.

If one starts with no fermions at the beginning, one may choose $`\alpha _0=1,\beta _0=0`$. Then, applying successive times the “transfer” matrix (34), one can get the spectrum of fermions produced after every $`\eta _n`$.

Of course our calculation reproduces the result

$$N_1=|\beta _1|^2=e^{\pi p_1^2}$$
(36)
reported in ref. .

We numerically integrated the evolution equations for $`u_\pm `$, $`\varphi `$, and the scale factor $`a`$.101010Our starting point is at $`\varphi (0)=0.28M_p`$, short after inflation, $`\varphi ^{}(0)=0.15M_pm_\varphi `$ (as follows from a numerical evaluation of the inflaton alone during inflation), and $`a(0)=1`$. The results obtained with the analytical expression (34) are always in very good agreement with the numerical ones. Just to give a couple of examples, we present here two cases at different regimes (we show them only with illustrative purpose, and the values of the parameters chosen have no particular importance). In figure 1 we present the spectrum of the fermions after two productions, that is after one complete oscillation of the inflaton field. In analogy with the bosonic case, we measure the strength of the coupling inflaton-fermions with the quantity $`qg^2\varphi _0^2/\left(4m_\varphi ^2\right)`$, where $`\varphi _00.28M_p`$ is the value of the inflaton at the beginning of reheating. In figure 1 we choose $`q=10^8`$, while we fix the bare fermion mass to be $`m_X=100m_\varphi `$. In figure 2 we show instead the resulting spectrum after $`7`$ productions in the case $`q=10^4`$, $`m_X=4m_\varphi `$.

## 4 Production in a non-expanding Universe

In the bosonic case, the study of the non-perturbative production in a non-expanding Universe has proven very useful in understanding the effects of the production. It is shown in ref.  that the bosonic wave function satisfies the Mathieu equations, whose solutions are characterized by resonance bands (in momentum space) of very “explosive” and efficient production. It is then shown that, due to the expansion of the Universe, modes of a given comoving momentum $`k`$ cross several resonance bands during the evolution. This gives the creation the stochastic character described in the work . In ref.  it is understood that an analogous behavior occurs also for fermions. The expansion is expected also in this case to spoil the clear picture of distinct resonance bands. This fact may help the transfer of energy to fermions, since the resonance bands in the fermionic case are anyhow limited by the Pauli principle. The expansion allows thus new modes to be occupied, and the production is no longer limited to the regions of resonance. In ref.  it is stated that the production should then almost completely fill the whole Fermi sphere up to a maximal momentum $`k_{\mathrm{max}}`$. This behavior is confirmed by the numerical results of the work . In this section we will see that our analytical formulae can reproduce the resonance bands, while in the next one we will discuss the effects of the expansion of the Universe.

Let us consider the matrices (we drop the index $`n`$ in the matrix elements since all the $`p_n`$’s have now the same value $`p`$)

$`M=\left(\begin{array}{cc}F& G\\ G& F^{}\end{array}\right),T_1=\left(\begin{array}{cc}e^{i\vartheta _1^2}& 0\\ 0& e^{i\vartheta _1^2}\end{array}\right),T_2=\left(\begin{array}{cc}e^{i\vartheta _2^3}& 0\\ 0& e^{i\vartheta _2^3}\end{array}\right),`$ (37)
with $`Ge^{\pi p^2/2}`$ and $`\vartheta _i^j_{\eta _i}^{\eta _j}\omega _k𝑑\eta `$ ($`F`$ was defined in the previous section).

Without the expansion of the Universe, the inflaton field has the periodic evolution

$$\varphi (\eta t)=\varphi _0\mathrm{cos}\left(m_\varphi \eta \right),$$
(38)
and all the $`\vartheta _i^j`$ (34) are hence sums of $`\vartheta _1^2`$ and $`\vartheta _2^3`$ (remember the $`\eta _i`$ are the moments at which the total fermionic mass vanishes).

After the generic $`n`$-th complete oscillation one thus has

$$\left(\begin{array}{c}\alpha _{2n}\\ \beta _{2n}^{}\end{array}\right)=\left(\begin{array}{cc}e^{i\vartheta _1^{2n+1}}& 0\\ 0& e^{i\vartheta _1^{2n+1}}\end{array}\right)T_2M^TT_1M\mathrm{}T_2M^TT_1M\left(\begin{array}{c}1\\ 0\end{array}\right)$$
(39)
with the combination $`\widehat{O}T_2M^TT_1M`$ repeated $`n`$ times.

One has thus simply to study the eigenvalue problem for $`\widehat{O}`$ (notice $`det\widehat{O}=1`$). This operator has the form

$`\widehat{O}`$ $`=`$ $`\left(\begin{array}{cc}A& B\\ B^{}& A^{}\end{array}\right),`$
$`A`$ $`=`$ $`F^2e^{i\left(\vartheta _1^2+\vartheta _2^3\right)}+G^2e^{i\left(\vartheta _1^2\vartheta _2^3\right)},`$
$`B`$ $`=`$ $`FGe^{i\left(\vartheta _1^2+\vartheta _2^3\right)}F^{}Ge^{i\left(\vartheta _1^2\vartheta _2^3\right)}`$ (40)
and its eigenvalues are $`\lambda _{1,2}=e^{\pm i\mathrm{\Lambda }}`$ with $`\mathrm{cos}\mathrm{\Lambda }=ReA`$.

Rewriting the initial condition $`\left(\begin{array}{c}1\\ 0\end{array}\right)`$ in terms of the eigenvectors of $`\widehat{O}`$ and substituting in formula (39), one gets the number of produced fermions

$$N_n=|\beta _n|^2=\frac{|B|^2}{1\left(\mathrm{Re}A\right)^2}\mathrm{sin}^2\left(n\mathrm{\Lambda }\right)=\frac{|B|^2}{\mathrm{sin}^2\mathrm{\Lambda }}\mathrm{sin}^2\left(n\mathrm{\Lambda }\right)$$
(41)
after the complete $`n`$-th oscillation.111111In ref.  it is said that it is possible to extract the average over one period of the occupation number from the knowledge of the solution of eq. (11) with initial conditions $`u_+\left(\eta _0\right)=1,u_+^{}\left(\eta _0\right)=0`$. The matching of this procedure with our formulae (22) and (40) gives the average time evolution (in our notation) $`\overline{N}(t)=|B|^2\mathrm{sin}^2[(\pi \mathrm{\Lambda })t/T]/\mathrm{sin}^2[\pi \mathrm{\Lambda }]`$, where $`T`$ is the period of one oscillation. Clearly, this result coincides with our eq. (41) at any given complete oscillation.

We notice the presence of the envelope function

$$\stackrel{~}{E}\frac{|B|^2}{1\left(ReA\right)^2}=\frac{1|A|^2}{1\left(ReA\right)^2}$$
(42)
which modulates the oscillating function $`\mathrm{sin}^2(n\mathrm{\Lambda })`$.

With increasing $`n`$ this last function oscillates very rapidly and can be at all effects averaged to $`1/2`$. One is thus left with the envelope function which shows the presence of resonance bands. The resonance bands occur where $`A`$ is real and $`\stackrel{~}{E}1`$. It is easy to understand that their width exponentially decreases with increasing momenta $`k`$. To see this, let us consider the behavior of $`\stackrel{~}{E}`$ at high momenta. In this regime, the function $`A`$ is given by

$$A(1e^{\pi p^2})e^{i\varphi _A},p=\frac{k}{\sqrt{g|\varphi _{}^{}|}},$$
(43)
where the phase $`\varphi _A`$ can be read from eq. (40). Near to the points where $`\mathrm{cos}\varphi _A=1`$ the envelope function behaves like121212A completely analogous behavior occurs where $`\mathrm{cos}\varphi _A=1`$.

$$\stackrel{~}{E}\frac{1(1e^{\pi p^2})^2}{1\left(1e^{\pi p^2}\right)^2\mathrm{cos}^2\varphi _A}\frac{1}{e^{\pi p^2}\left(1\mathrm{cos}\varphi _A\right)+1}.$$
(44)
The width of the band can be defined as the distance between the two successive points at which $`\stackrel{~}{E}=1/2`$. From the last expression it follows that the difference between the phases of $`A`$ in these two points is given by $`\mathrm{\Delta }\alpha =2\sqrt{2}e^{\pi p^2}`$. Since the most rapidly varying term which contributes to the phase of $`A`$ is $`\left(p^2/2\right)\mathrm{log}\left(p^2/2\right)`$, the width of the band can be thus estimated to be

$$\mathrm{\Delta }p\frac{2\sqrt{2}}{p\mathrm{log}\left(p^2/2\right)}e^{\pi p^2}.$$
(45)
We show in figure 3 the envelope of the produced fermions in a static Universe for the parameters $`q=10^6`$ and $`m_X=100m_\varphi `$. The peaks occur where $`A`$ is real and it is confirmed that their width decreases very rapidly at increasing momenta. Due to the fact that the last peaks plotted are indeed very sharp, the resolution of figure 3 does not allow to see their top at $`n_k=1`$.

## 5 Expansion taken into account

As stated in the previous section, the resonance bands disappear when the expansion of the Universe is taken into account. In this section we will show how the occupation number varies when a non-vanishing Hubble parameter is considered. As we have seen in section 2, eqs. (34) and (35) give a very good agreement with the numerical results. On the other hand, the presence of phases in eq. (35) makes the exact analytical treatment of the occupation number impossible. Now, the same observation made in the bosonic case  turns very useful also to us: the phases in eq. (35), when the expansion of the Universe is taken into account, are not correlated among themselves. As a result, the final spectra present several high frequency oscillations about some average function. The positions of the peaks of these oscillations depend on the details of the phases. However, the “mean” function can be understood in a surprisingly easy way. Our problem can be treated as one customary does when dealing with the “random walk”. Imagine one has to calculate the quantity

$$S|A_1+A_2+A_3+\mathrm{}+A_n|^2,$$
(46)
where the $`A_i`$ are complex numbers with random phases. The “random walk recipe” indicates that the best estimation of the above quantity is achieved by summing the squares of the terms $`A_i`$, since the mixed products average to zero. The chaoticity of the final spectra for $`n_k`$ suggests that this may also be true in our case, and this is confirmed by comparison with the numerical results.

With this method, eqs. (34) and (35) turn into the much simpler relations

$$\left(\begin{array}{c}|\alpha _n|^2\\ |\beta _n^{}|^2\end{array}\right)=\left(\begin{array}{cc}|F_n|^2& |H_n|^2\\ |H_n|^2& |F_n|^2\end{array}\right)\left(\begin{array}{c}|\alpha _{n1}|^2\\ |\beta _{n1}^{}|^2\end{array}\right),$$
(47)
where we remember

$$|F_n|^2=1e^{\pi p_n^2},|H_n|^2=e^{\pi p_n^2},|F_n|^2+|H_n|^2=1.$$
(48)
Assuming no fermions in the initial state, and applying $`n`$ times this formula, it is easy to see that the occupation number after the $`n`$-th production is given by

$$N_n(k)=\frac{1}{2}\frac{1}{2}\underset{i=1}{\overset{n}{}}\left(12e^{\pi p_i^2}\right).$$
(49)
A similar result holds for preheating of bosons, cfr.  where the idea of averaging on almost random phases is first introduced. In the bosonic case, one can exploit the fact that, due to the high efficiency of the production, the occupation number after the $`(n+1)`$-th creation is (almost) proportional to the occupation number after the $`n`$-th one:

$$N_{n+1}\left(k\right)\left(1+2e^{\pi \kappa _n^2}\right)N_n\left(k\right),$$
(50)
where the quantity $`\kappa _n`$ is analogous to our parameter $`p_n`$. This simplification is not possible in our case, since the Pauli principle forbids $`N_n`$ to be sufficiently high. However our final result, eq. (49), is also cast in a very simple and immediate form.

The validity of eq. (49) is confirmed by our numerical investigations, as we show here in one particular case. In figure 4 we compare the behavior of the “mean” function for the spectra with respect to the numerical one. We choose the physical parameters to be $`q=10^4,m_X=4m_\varphi `$, and we look at the results after the $`7`$-th production (this corresponds to the choices made in figure 2). We see that the “mean” function interpolates very well between the maxima and the minima of the numerical spectrum, and that it can indeed be considered as a very good approximation of the actual result.

This is confirmed by figure 5, where we plot (for the same values adopted in figure 4) the quantity $`n_kk^2`$ rather than the occupation number alone. This quantity is of more physical relevance when one is interested in the total energy transferred to the fermions, since the total number density of produced fermions is (apart from the dilution due to the expansion of the Universe, that will be considered only in the final result)

$$N_X=\frac{2}{\pi ^2}𝑑kk^2n_k.$$
(51)
As shown in the plot, the result for $`N_X`$ in the numerical and in the approximated case are in very good agreement.

After checking the validity of the approximation given by the “mean” function, we adopt it to understand how the production scales when different values of the parameters $`q`$ and $`m_X`$ are considered.

Equation (49) allows us to give an analytical estimate of the total amount of energy stored in the fermions after the $`n`$-th production, and in particular after that the whole process of non-perturbative production has been completed. Notice that all the dependence on the physical parameters is in the coefficients

$$z_i\frac{k}{\left(\pi ^{1/2}p_i\right)},$$
(52)
where $`k`$ is the comoving momentum, and the only part we have to determine are the numbers

$$z_i^2=\frac{2\sqrt{q}}{\pi }a\left(\eta _i\right)\frac{\left|\varphi ^{}\left(\eta _i\right)\right|}{\varphi _0}.$$
(53)
Here, and in what follows, we express the dimensionful quantities in units of $`m_\varphi `$, the inflaton mass, apart from the inflaton field $`\varphi `$ that is given in units of $`M_p`$.

We also introduce the ratio $`R2q^{1/2}/m_X`$. This quantity is the most relevant, since it determines the zeros of the total fermionic mass and thus the values of the $`z_i`$’s. Indeed, from eq. (2) we see that the total mass vanishes for $`\varphi _{}=\varphi _0/R`$. It is convenient to study the production in terms of the two independent parameters $`q`$ and $`R`$ (rather than $`q`$ and $`m_X`$) since, at fixed $`R`$, all the spectra are the same provided we rescale $`kq^{1/4}`$ (cfr. eq. (53)).

One can now proceed in two different ways, and we devote the next two subsections to each of them. First, one can evolve the equations for the inflaton field alone and find numerically the values $`z_i`$ for given $`q`$ and $`R`$. Inserting these values in eq. (49) one can get final values for the production which, as we have reported, are very close to the numerical ones. This method allows to get results which average the actual ones, and it has the advantage of being much more rapid than a full numerical evolution.

Alternatively, one can proceed with a full analytical study in order to understand the results given by the first semi-analytical method.

### 5.1 Semi-analytical results

In this subsection we evaluate the analytical formula (49), taking the coefficients $`z_i`$ from a numerical evolution of the inflaton field. As we have said, this method gives results which are very close to the numerical ones, due to the fact that the “mean” function interpolates very well the actual spectra of the produced fermions.

The first thing worth noticing is that, for each choice of $`q`$ and $`R`$, the maximum $`z_i`$ occurs at about half of the whole process of non-perturbative creation. It thus follows that fermions of maximal comoving momenta will be mainly produced at the half of the process. Our semi-analytical evaluations support this idea: the total energy stored in the created fermions increases slowly during the second part of the preheating. To see this, we show in figure 6 the numerical results for the quantity $`N_X=2𝑑kk^2n_k/\pi ^2`$ as function of the number of production. We fix the physical parameters to be $`q=10^6`$ and $`R=30`$, that is $`m_X=67m_\varphi `$. With this choice, the total mass vanishes $`10`$ times, and figure 5 shows the results after each step of this production. We see indeed that the final productions are less efficient than the previous ones.

From eq. (49) it is also possible to show the evolution of the spectra with the number of productions. We do this in figure 7. We choose the same parameters as in figure 6, and we show the results after each complete oscillation of the inflaton field (that is after each two productions). We observe that the production rapidly approaches a step function in the momentum space, i.e. there exists a maximum momentum below which Pauli blocking is saturated (notice that the value $`1/2`$ follows from the average understood in the “mean” function), and above which $`n_k0`$. The fact that the last productions do not contribute much to the total energy is also confirmed.131313The behavior at small $`k`$ is inessential, since this region does not significantly contribute to the total energy.

We now turn our attention to the total energy transferred to fermions after the whole preheating process is completed. We fix the parameter $`q`$ to the value $`10^6`$ and we investigate how the total integral $`N_X`$ changes with different values for the parameter $`R`$.141414As reported, the scaling of the final result with $`q`$ at fixed $`R`$ is simply understood from eq. (53). The results are shown in figure 8, for $`R`$ ranging from $`5`$ to $`10000`$. For the last value the total fermionic mass changes sign more than $`3700`$ times and a full numerical evaluation would appear very problematic. This can be done in our case, thanks to the analytical expression (49) found, and our results extend the validity region of the previous full numerical study .

In figure 8, the results of our semi-analytical method are also compared to the full analytical ones of the next subsection. This comparison will be discussed below.

From the scaling of $`N_X`$ with $`R`$ just reported, it is now easy to estimate the total energy transferred to fermions for generic values of $`q`$ and $`R`$. We are interested in comparing our results to the numerical ones of the work . To do so, we consider the ratio between the energy density given to fermions and the one in the inflaton field151515This ratio should be calculated at a time $`t_{\mathrm{end}}`$ at the end of preheating, when the total fermionic mass stops vanishing. Thus in the denominator of eq. (54) the comoving momentum $`k`$ should be replaced by the physical one $`p=k/a`$, with $`a`$ scale factor of the Universe at $`t_{\mathrm{end}}`$. Analogously, the value $`\varphi _0`$ of the beginning of reheating should be replaced by the one at $`t_{\mathrm{end}}`$. However, both the replacements cancel out in the ratio, since both $`\rho _\chi `$ and $`\rho _\varphi `$ redshift as energy densities of matter.

$$\frac{\rho _X}{\rho _\varphi }=\frac{2m_X}{\pi ^2}𝑑kk^2n_k\frac{2}{m_\varphi ^2\varphi _0^2}.$$
(54)
We present our results in figure 9. We report them in terms of $`q`$ and $`m_X`$, in order to have an immediate comparison to the analogous figure 1 of ref. .

In figure 9, for any fixed $`q`$, the greatest plotted value for $`m_X`$ corresponds to the choice $`R=5`$. We are not interested in extending this limit since we know that for greater $`m_X`$ (actually for values greater than the bound $`m_X\sqrt{q}/2`$) the production suddenly stops. The smallest value plotted for $`m_X`$ (at any fixed value $`q`$) corresponds instead to $`R=10000`$, that is to considering more than $`3700`$ productions in the numerical evaluation of eq. (49).

Our final values are in good agreement with the ones of figure 1 of ref.  in the regime of validity of the latter. The numerical results reported in that figure exhibit small fluctuations about an average function $`\rho _X\left(m_X\right)`$. Our results give this average function. This was expected, since the expression that we integrated, eq. (49), interpolates between the maxima and the minima of the numerical spectra.

The numerical results of ref.  have a smaller range of validity than the ones here presented. This occurs because the numerical evolution of that work is limited to the first $`20`$ oscillations of the inflaton field, and so fermionic production has not come to its end for small values of $`m_X`$. We confirm that at high values of $`m_X`$ (actually at small values of $`R`$ for any given $`q`$) the production depends very weakly on $`m_X`$. In addition, our results show a decrease of the energy transferred to fermions for smaller values of $`m_X`$. This behavior will be explained in details in the next subsection.

### 5.2 Analytical results

We want now to show that all the results presented in the previous section can be also achieved with a full analytical study of eq. (49).

First of all, we have to estimate the quantities $`z_i`$ given in eq. (53). To do this, it is more convenient to work in terms of the physical time $`t`$: after the first few oscillations, the inflaton evolution is very well approximated by the expression (remember $`t`$ is expressed in units of $`m_\varphi ^1`$, while $`\varphi `$ in units of $`M_p`$)

$$\varphi (t)\frac{1}{\sqrt{3\pi }}\frac{\mathrm{cos}\left(t\right)}{t}.$$
(55)
The scale factor of the Universe follows the “matter-domination” law, and it is well approximated by $`a(t)=t^{2/3}`$.

The values $`t_i`$ are determined by the condition of vanishing of the total mass of the fermions, that is, by making use of eq. (55),

$$\mathrm{cos}\left(t_i\right)=\frac{t_i}{RA},$$
(56)
where we remind $`R(2\sqrt{q})/m_X`$. The parameter $`A(\sqrt{3\pi }\varphi _0)^1`$ is of order one and will not play any special role in what follows. Notice that the last production occurs at $`tRA`$.

Hence, keeping only the dominant contribution to the derivative of $`\varphi `$ with respect to the physical time, we get the expression

$$z_i^2\frac{2\sqrt{q}}{\pi }R^{1/3}A^{4/3}\left(\frac{t_i}{RA}\right)^{1/3}\sqrt{1\left(\frac{t_i}{RA}\right)^2},$$
(57)
where we can assume $`t_ii\pi `$.

Equation (57) exhibits a very good agreement with the numerical evaluation of the same quantity. It also shows that the maximal value for $`z_i`$ is reached at $`t_i=AR/2`$, that is, at half of the whole process of non-perturbative creation. This was anticipated in the previous section, where we showed that the most of the fermionic production occurs in the first part of preheating.

Starting from eq. (57) we can also calculate the number density of produced particles

$$N_X(q,R)=\frac{2}{\pi ^2}𝑑kk^2N(k),$$
(58)
where $`N(k)`$ is obtained from eq. (49) with $`n=n_{\mathrm{max}}=(RA)/\pi `$.

For $`R`$ big enough, the product in eq. (49) can be written as the exponential of an integral. Thus, we obtain

$`N_X`$ $`=`$ $`{\displaystyle \frac{m_X}{4\pi ^3}}\left({\displaystyle \frac{2A^{4/3}}{\pi }}\right)^{3/2}q^{3/4}R^{1/2}\times `$
$`\times {\displaystyle _0^{\mathrm{}}}dtt^2\{1\mathrm{exp}\left[{\displaystyle \frac{RA}{\pi }}{\displaystyle _0^1}dy\mathrm{log}\right|12\mathrm{exp}({\displaystyle \frac{t^2}{y^{1/3}\sqrt{1y^2}}})\left|\right]\},`$
with the substitution $`t=k\sqrt{\pi /\left(2A^{4/3}q^{1/2}R^{1/3}\right)}`$.

The integral in $`dy`$ which appears in eq. (5.2) cannot be calculated analytically. Anyhow, we can approximate it by

$$_0^1𝑑y\mathrm{log}\left|12\mathrm{exp}\left(\frac{t^2}{y^{1/3}\sqrt{1y^2}}\right)\right|g\left(t\right)\mathrm{log}\left|12e^{1.5t^2}\right|,$$
(60)
where $`g(t)`$ is a function of order one that, for our purposes, can be approximated by a constant $`c`$ in the range $`0.5c1`$.

Hence, the integrand within curly brackets in eq. (5.2) rewrites

$$1\left|12e^{1.5t^2}\right|^{c\frac{RA}{\pi }}.$$
(61)
In the large-$`R`$ limit this function approximates a step function, which evaluates to one for

$$\sqrt{\frac{\pi \mathrm{log}2}{2RAc}}t\sqrt{\mathrm{log}\left(\frac{2RAc}{\pi \mathrm{log}2}\right)}$$
(62)
and to zero for the remaining values of $`t`$.

Since the quantity (61) is proportional to the occupation number $`n_k`$, our analytical calculation confirms the usual assumption that, after few oscillations, the fermionic production saturates the Fermi sphere up to a given maximum momentum $`k_{\mathrm{max}}`$. This is also shown in the previous figure 7. However, our derivation gives a different scaling for $`k_{\mathrm{max}}`$ with respect to the previous literature . Indeed, from eq. (62) it follows (apart from proportionality factors)

$$k_{\mathrm{max}}\frac{q^{1/3}}{m_X^{1/6}}\sqrt{\mathrm{log}\left(\frac{q^{1/2}}{m_X}\right)},$$
(63)
which is however quite close to the result given in . We will comment more about the origin of the scaling (63) in the conclusions.

From eqs. (5.2) and (62), we obtain the final expression for the number density of the fermions created during the whole process

$$N_X(q,R)=\frac{1}{3\pi ^2}\left(\frac{2A^{4/3}}{1.5\pi }\right)^{3/2}q^{3/4}R^{1/2}\left[\mathrm{log}\left(\frac{4Ac}{\pi \mathrm{log}2}\frac{q^{1/2}}{m_X}\right)\right]^{3/2}.$$
(64)
We can now go back to figure 8, where this last equation (called “analytical” in the figure) is compared to the result with the semi-analytical method of the previous subsection. In plotting eq. (64) we chose $`c=0.78`$ for the numerical factor involved. As we can see, the final results achieved with the two methods are in very good agreement with each other, thus confirming the validity of the formula (64).161616The small discrepancy between the two curves can be attributed to the fact that $`c`$ is not exactly constant.

Rewriting eq. (64) in terms of $`q`$ and $`m_X`$ we can draw some conclusions. First, apart from a logarithmic correction, the scaling of the total energy

$$\rho _Xm_XN_Xqm_X^{1/2}\left[\mathrm{log}\left(\frac{4Ac}{\pi \mathrm{log}2}\frac{q^{1/2}}{m_X}\right)\right]^{3/2}$$
(65)
is linear in $`q`$, as generally expected . The dependence of $`\rho _X`$ on $`m_X`$ requires some more care: the threshold value for $`m_X`$ is given by the condition $`R4`$, that is,171717The number $`4`$ comes from the fact that the value of the inflaton at its first minimum is $`\varphi 0.07M_p`$, while at beginning $`\varphi _0=0.28M_p`$.

$$\left(m_X\right)_{\mathrm{th}}\frac{\sqrt{q}}{2}.$$
(66)
For values of $`m_X`$ not much smaller than $`(m_X)_{\mathrm{th}}`$, the total energy depends very weakly on $`m_X`$, this result being in agreement with the numerical evaluation given in . On the other hand, for values of $`m_X`$ much smaller than $`(m_X)_{\mathrm{th}}`$, the factor $`m_X^{1/2}`$ starts to dominate, and we expect it to determine the scaling of the total energy when $`m_X0`$.

## 6 Backreaction

The results presented so far have been achieved neglecting the backreaction of the produced fermions on the evolution of the inflaton field and on the scale factor. This is a common approximation, since a more complete treatment (especially an analytical one) of the whole phenomenon is a very difficult task. However, the effects of the backreaction can be understood to a good degree of accuracy in the Hartree approximation.

For what concerns preheating of fermions, the Hartree approximation consists in taking into account the term

$$g\overline{X}X$$
(67)
into the evolution equations for the inflaton and the scale factor. The equation for the field $`\varphi `$ thus rewrites (in physical time)

$$\ddot{\varphi }+3H\dot{\varphi }+m_\varphi ^2\varphi +g\overline{X}X=0.$$
(68)
The study of this effect has been performed numerically in ref. , where it is shown that backreaction starts to be important for $`q>q_b10^810^{10}`$. Figure 7 of that work shows how the evolution of the inflaton field $`\varphi `$ is modified when backreaction is considered and $`q`$ is sufficiently high. First, one observes that the amplitude of the oscillations of the inflaton is very damped already after the first production. This effect is the most obvious one, since the term (67) takes into account the decay of the inflaton into fermion-antifermion pairs, while in its absence the equation for $`\varphi `$ considers only the damp due to the expansion of the Universe. The second feature that emerges from the evolution performed in ref.  is that at the beginning the field $`\varphi `$ does not oscillate about the minimum of the potential $`V=m_\varphi ^2\varphi ^2/2`$, but about the point $`\varphi _{}`$ where the total fermionic mass vanishes. Moreover, the frequency of these oscillations is higher than $`m_\varphi `$.

These last two effects are due to the change in the effective potential for $`\varphi `$ induced by the term (67), and disappear when the quantity $`\overline{X}X`$ is decreased by the expansion of the Universe. Their net effect is to render the whole mechanism of preheating more efficient, since the rise in the frequency of the oscillations of $`\varphi `$ increases the number of productions of fermions. In ref.  it is indeed shown that for $`q=10^8`$ the total production is about $`5\%`$ larger than the one without backreaction, while for $`q=10^{10}`$ the increase is about $`50\%`$.

We now briefly study the evolution of the inflaton $`\varphi `$ under eq. (68) by means of the analytical results presented above. We show that even a very approximate analysis confirms the numerical result that indicates in $`q10^810^{10}`$ the threshold above which backreaction should be considered.

To begin on, the term (67) needs to be normal ordered. Doing so, one gets

$$\overline{X}X=\frac{2}{\left(2\pi a\right)^3}d^3𝐤\left(1+\frac{ma}{\omega }|u_{}|^2\right).$$
(69)
In terms of the Bogolyubov coefficients, this quantity evaluates to

$$\overline{X}X=\frac{4}{\left(2\pi a\right)^3}d^3𝐤\left[|\beta |^2\frac{am}{\omega }+Re\left(\alpha \beta \frac{k}{\omega }e^{2i{\scriptscriptstyle \omega 𝑑\eta }}\right)\right].$$
(70)
We see that $`\overline{X}X`$ vanishes for $`\beta =0`$. This is obvious, since backreaction starts only after fermions are produced. Some approximations can render eq. (70) more manageable. First, we notice that the oscillating term in the exponential averages to very small values the integral of the second term in square brackets. Second, we see that $`k|am|`$ where $`\beta _k`$ is significantly different from zero. From both these considerations, the integrand in eq. (70) can be approximated (up to the sign of $`m`$) by the occupation number $`|\beta |^2`$, so that the whole effect is (approximatively) proportional to the number of produced fermions.

The numerical results of ref.  show that, for the values of $`q`$ for which the backreaction is to be considered, its effects can be seen already in the first oscillation of the inflaton field. Since we are only interested in estimating the order of magnitude of $`q_b`$, we thus concentrate on the first oscillation of $`\varphi `$, neglecting the expansion of the Universe in this short interval.

With all these approximations, eq. (68) rewrites

$$y^{\prime \prime }+y+10^{12}q^{5/4}\frac{m}{|m|}=0,$$
(71)
where we have rescaled $`y\varphi /\varphi _0`$ and we remind that the time is given in units of $`m_\varphi ^1`$.

This equation is very similar to the one obtained in the bosonic case . The last term changes sign each time $`m=0`$ and, when sufficiently high, forces the inflaton field to oscillate about the point $`\varphi _{}`$ at which the total fermionic mass vanishes. It is also responsible for the increase of the frequency of the oscillations. To see this, we assume that this term dominates over the second one of eq. (71) and we solve it right after the first fermionic production at the time $`t_1`$. In the absence of the second term, eq. (71) is obviously solved by a segment of parabola until the time $`t_2t_1+2|y^{}\left(t_1\right)|/\left(q^{5/4}10^{12}\right)`$ at which $`m`$ vanishes again and the last term of eq. (71) changes sign. As long as the second term of eq. (71) can be neglected, the inflaton evolution proceeds along segments of parabola among successive zeros of the mass $`m`$. The “time” duration of these segments is expected to be of the same order of the first one, since the successive fermionic productions balance the decrease of $`\overline{X}X`$ due to the expansion of the Universe.181818However, after the first part of the process, the production looses its efficiency and the expansion of the Universe dominates. As we have said, the term $`\overline{X}X`$ can then be neglected and the inflaton starts oscillating about the minimum of the tree level potential.

The period of these oscillations can thus be roughly estimated to be $`T2\left(t_2t_1\right)`$. We see that, for $`q10^9`$, this period is smaller than the one that the inflaton oscillations would have neglecting backreaction. Since this increase of the frequency is the main responsible for the higher fermionic production, the result $`q_b10^9`$ can be considered our estimate for the value of $`q`$ above which backreaction should be taken into account.

This result, although obtained with several approximations, is in agreement with the numerical one of ref. .

## 7 Conclusions

Preheating of fermions and bosons present several analogies. In both cases, especially when one is interested in the production of very massive particles, the creation occurs for very short intervals of time, during which the frequency of the particles varies non adiabatically and their occupation number cannot be defined. This strong similarity can be exploited to extend to the fermionic case the formalism developed in ref.  for the analytical study of preheating of bosons. This was first done in ref. , where the result for the particles created in the first production is reported. When one wants to consider the successive productions, a more detailed analysis is necessary. This has been the object of the first part of the present work. The formulae that we obtained, valid for an arbitrary number of productions, give a very good agreement with the numerical results, as some examples provided manifest.

Another strong similarity between the two phenomena is that in both cases the production would occur only through resonance bands (in momentum space) was not for the expansion of the Universe. For what concerns fermions, this feature was first studied in ref. , where production in a static Universe is considered. As it is well known, the expansion of the Universe removes the resonance bands and, as a result, this allows the fermionic production to saturate the whole Fermi sphere up to a maximal momentum $`k_{\mathrm{max}}`$. Our analytical results confirm both the presence of the resonance bands in the static case (where our formulae considerably simplify) and their disappearance when the expansion is considered.

Rather than the detailed form of the spectrum, a quantity which may be of more physical interest is the total energy of the produced particles. The analogy with the bosonic case turns useful also in this regard. In the final spectra one notices the presence of some maxima and minima, whose exact position is determined by the detailed knowledge of the phases of the Bogolyubov coefficients (adopted in the analytical derivation) after each production. However, due to the expansion of the Universe, these phases are effectively uncorrelated among themselves. This gives the production a stochastic character which is responsible for the disappearance of the resonance bands. Averaging over these phases, it is possible to get an analytical function for the mean occupation number. For what concerns preheating of bosons, the results are simplified by the fact that, in the region of very efficient production, the mean occupation number after a given production is proportional to the mean occupation number of the previous one. This assumption is not possible in the fermionic case, since the Pauli principle prevents the occupation number to exceed $`1`$. However, it is remarkable that also for fermions the final result can be cast in a very simple and readable form.

The analytical formula for the mean occupation number confirms the saturation of the Fermi sphere (times the factor $`1/2`$ which comes from the average) up to a momentum $`k_{\mathrm{max}}`$. In order to get this quantity, it is sufficient to know the derivative of the inflaton field and the value of the scale factor at the points where the productions occur. These values can be calculated numerically evolving the equations for the inflaton field alone (at least when backreaction is neglected), or analytically. The former possibility is much more rapid than a full numerical computation, since one does not have to consider the equations for the fermions. Moreover, its results agree very well with the full numerical ones . Thanks to the increase in the rapidity, this evaluation (that is called “semi-analytical” in the present work) allows to get results with a more extended range of validity with respect to the previous full numerical ones.

The results provided in this way can also be achieved with a full analytical study. From both these methods, we deduce that the ratio $`\rho _X/\rho _\varphi `$ scales as

$$\frac{\rho _X}{\rho _\varphi }qm_X^{1/2}\left[\mathrm{log}\frac{q^{1/2}}{m_X}\right]^{3/2},$$
(72)
up to $`\left(m_X\right)_{\mathrm{th}}\sqrt{q}/2`$.

For $`m_X`$ not too smaller than $`\left(m_X\right)_{\mathrm{th}}`$ this scaling is in good agreement with the one $`\rho _X/\rho _\varphi q`$ of the numerical work . However, we see that the density of produced fermions decreases at smaller values of $`m_X`$.

We would like to conclude our work showing that the scaling (72) can be also achieved from very immediate considerations. As we said, preheating of fermions saturates a Fermi sphere in momentum space up to a maximal momentum $`k_{\mathrm{max}}`$. From the analytical formula (47), we notice that at high momenta $`k`$ the occupation number is well approximated by

$$N_n(k)\underset{i=1}{\overset{n}{}}e^{k^2/z_i^2},$$
(73)
where we remember $`z_iq^{1/4}a^{1/2}\left(\eta _i\right)|\varphi ^{}\left(\eta _i\right)|^{1/2}`$.

In this last equation, we replace all the parameters $`z_i`$ with a mean value $`\overline{z}`$, so that $`N_nn\mathrm{exp}\left(k^2/\overline{z}^2\right)`$. The scaling of $`\overline{z}`$ with the physical parameters $`q`$ and $`m_X`$ follows from the scaling of all the $`z_i`$. The maximal momentum $`k_{\mathrm{max}}`$ is thus expected to scale as the quantity $`z_i(\mathrm{log}n)^{1/2}`$. Considering now the evolution of the inflaton field in physical time $`t`$, we notice that both the number $`n`$ of productions and the times $`t_i`$ at which they occur are proportional to the parameter $`R=q^{1/2}/(2m_X)`$. Moreover, we see that the $`z_i`$’s scale as

$$z_iq^{1/4}a_i\left[\frac{d\varphi }{dt}|_i\right]^{1/2}q^{1/4}t_i^{2/3}\frac{1}{t_i^{1/2}}q^{1/4}R^{1/6}.$$
(74)
We thus get $`k_{\mathrm{max}}q^{1/4}R^{1/6}[\mathrm{log}R]^{1/2}`$, from which (remember $`\rho _Xm_Xk_{\mathrm{max}}^3`$) the scaling (72) simply follows.

###### Acknowledgments.

We would like to thank Antonio Riotto for some suggestions and for very useful discussions.