Backreaction of electromagnetic fields and the Schwinger effect in pseudoscalar inflation magnetogenesis

TL;DR

This paper investigates how electromagnetic fields generated during pseudoscalar inflation are affected by backreaction and the Schwinger effect, impacting magnetogenesis and the properties of resulting magnetic fields.

Contribution

It derives equations accounting for backreaction and Schwinger effect in pseudoscalar inflation, analyzing their influence on magnetic field generation.

Findings

Backreaction limits electromagnetic energy density during inflation.

Schwinger effect becomes significant when electric energy density exceeds a threshold.

Maximal magnetic field strength of 10^{-15} G with 1 pc correlation length at present.

Abstract

We study magnetogenesis in axionlike inflation driven by a pseudoscalar field $\phi$ coupled axially to the electromagnetic (EM) field $(\beta/M_{p})\phi F_{\mu\nu}\tilde{F}^{\mu\nu}$ with dimensionless coupling constant $\beta$. A set of equations for the inflaton field, scale factor, and expectation values of quadratic functions of the EM field is derived. These equations take into account the Schwinger effect and the backreaction of generated EM fields on the Universe expansion. It is found that the backreaction becomes important when the EM energy density reaches the value $\rho_{\rm EM}\sim (\sqrt{2\epsilon}/\beta)\rho_{\rm inf}$ ($\epsilon$ is the slow-roll parameter and $\rho_{\rm inf}$ is the energy density of the inflaton) slowing down the inflaton rolling and terminating magnetogenesis. The Schwinger effect becomes relevant when the electric energy density exceeds the value $\rho_{E}\sim \alpha_{\rm EM}^{-3} (\rho_{\rm tot}^{2}/M_{p}^{4})$, where $\rho_{\rm tot}=3H^{2}M_{p}^{2}$ is the total energy density and $\alpha_{\rm EM}$ is the EM coupling constant. For large $\beta$, produced charged particles could constitute a significant part of the Universe energy density even before the preheating stage. Numerically studying magnetogenesis in the $\alpha$-attractor model of inflation, we find that it is possible to generate helical magnetic fields with the maximal strength $10^{-15}\,{\rm G}$, however, only with the correlation length of order $1\,{\rm pc}$ at present.

Figures

Equations

$$g_{\mu \nu }=\mathrm{diag}(1,-a^2,-a^2,-a^2),\sqrt{-g}=a^3.$$

$$S=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi -V(\varphi )-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\frac{1}{4}I(\varphi )F_{\mu \nu }{\tilde{F}}^{\mu \nu }+{\mathcal{L}}_{\mathrm{charged}}(A_{\nu },\chi )\right],$$

$$I(\varphi )=\beta \frac{\varphi }{M_p},$$

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left[\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu }\varphi \right]+\frac{dV}{d\varphi }=-\frac{1}{4}{I}^{\prime }(\varphi )F_{\mu \nu }{\tilde{F}}^{\mu \nu },$$

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left[\sqrt{-g}{F}^{\mu \nu }\right]+I^{\prime }(\varphi ){\tilde{F}}^{\mu \nu }{\partial }_{\mu }\varphi =-j^{\nu },$$

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left[\sqrt{-g}{\tilde{F}}^{\mu \nu }\right]=0,$$

$$j^{\mu }=\frac{\partial {\mathcal{L}}_{\mathrm{charged}}(A_{\nu },\chi )}{\partial A_{\mu }}=(n,\mathbf{j})$$

$$T_{\mu \nu }=\frac{2}{\sqrt{-g}}\frac{\delta S}{\delta g^{\mu \nu }}={\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -g^{\lambda \rho }{F}_{\mu \lambda }{F}_{\nu \rho }-g_{\mu \nu }{\mathcal{L}}_0+T_{\mu \nu }^{(\mathrm{charged})},$$

$$\mathbf{E}=-\frac{1}{a}\dot{\mathbf{A}},\mathbf{B}=\frac{1}{a^2}\mathrm{rot}\mathbf{A},$$

$$F^{0i}=\frac{1}{a}{E}^i,F_{ij}=a^2{\epsilon }_{ijk}{B}^k,{\tilde{F}}^{0i}=\frac{1}{a}{B}^i,{\tilde{F}}_{ij}=-a^2{\epsilon }_{ijk}{E}^k,$$

$$\rho =\left[\frac{1}{2}{\dot{\varphi }}^2+V(\varphi )\right]+\frac{1}{2}⟨E^2+B^2⟩+{\rho }_{\chi }={\rho }_{\inf }+{\rho }_{\mathrm{EM}}+{\rho }_{\chi },$$

$$H^2=\frac{1}{3M_p^2}\left({\rho }_{\inf }+{\rho }_{\mathrm{EM}}+{\rho }_{\chi }\right),$$

$$\ddot{\varphi }+3H\dot{\varphi }+\frac{dV}{d\varphi }=I^{\prime }(\varphi )⟨\mathbf{E}\cdot \mathbf{B}⟩,$$

$$\dot{\mathbf{E}}+2H\mathbf{E}-\frac{1}{a}\mathrm{rot}\mathbf{B}+I^{\prime }(\varphi )\dot{\varphi }\mathbf{B}=-a\mathbf{j},$$

$$\dot{\mathbf{B}}+2H\mathbf{B}+\frac{1}{a}\mathrm{rot}\mathbf{E}=0,$$

$$\mathrm{div}\mathbf{B}=0,\mathrm{div}\mathbf{E}=0.$$

$$\mathbf{j}=\frac{1}{a}\sigma \mathbf{E}.$$

$${\mathcal{E}}^{(n)}$$

$${\mathcal{G}}^{(n)}$$

$${\mathcal{B}}^{(n)}$$

$$⟨\mathrm{rot}\mathbf{f}\cdot \mathbf{g}⟩-⟨\mathbf{f}\cdot \mathrm{rot}\mathbf{g}⟩=\mathrm{div}⟨[\mathbf{f}\times \mathbf{g}]⟩=0$$

$${\mathcal{B}}^{(-1)}\equiv \mathcal{H}=\frac{1}{a^3}⟨\mathbf{A}\cdot \mathrm{rot}\mathbf{A}⟩=\frac{1}{a}⟨\mathbf{A}\cdot \mathbf{B}⟩,$$

$$\displaystyle\dot{\mathscr{E}}^{(n)}+[(n+4)H+2\sigma]\mathscr{E}^{(n)}-2I^{\prime}(\phi)\dot{\phi}\mathscr{G}^{(n)}+2\mathscr{G}^{(n+1)}=\big{[}\dot{\mathscr{E}}^{(n)}\big{]}_{b},$$

$$\displaystyle\dot{\mathscr{G}}^{(n)}+[(n+4)H+\sigma]\mathscr{G}^{(n)}-\mathscr{E}^{(n+1)}+\mathscr{B}^{(n+1)}-I^{\prime}(\phi)\dot{\phi}\mathscr{B}^{(n)}=\big{[}\dot{\mathscr{G}}^{(n)}\big{]}_{b},$$

$$\displaystyle\dot{\mathscr{B}}^{(n)}+(n+4)H\mathscr{B}^{(n)}-2\mathscr{G}^{(n+1)}=\big{[}\dot{\mathscr{B}}^{(n)}\big{]}_{b}.$$

$${\dot{\rho }}_{\mathrm{EM}}+4H{\rho }_{\mathrm{EM}}+2\sigma {\rho }_E=I^{\prime }(\varphi )\dot{\varphi }{\mathcal{G}}^{(0)}+[{\dot{\rho }}_{\mathrm{EM}}{]}_b.$$

$${\dot{\rho }}_{\inf }+3H({\rho }_{\inf }+p_{\inf })=-I^{\prime }(\varphi )\dot{\varphi }{\mathcal{G}}^{(0)},$$

$$|{\mathcal{G}}^{(0)}|\sim \sqrt{{\mathcal{E}}^{(0)}{\mathcal{B}}^{(0)}}=2\sqrt{{\rho }_E{\rho }_B}\sim {\rho }_{\mathrm{EM}}.$$

$${\rho }_{\mathrm{EM}}\sim \frac{M_p}{\beta }|V^{\prime }(\varphi )|\sim \frac{\sqrt{2ϵ}}{\beta }{\rho }_{\inf },$$

$$ϵ=\frac{M_p^2}{2}{\left(\frac{V_{\varphi }^{\prime }}{V}\right)}^2$$

$$\hat{\mathbf{A}}(t,\mathbf{x})=\int \frac{d^3\mathbf{k}}{(2\pi {)}^{3/2}}\sum _{\lambda =\pm }\left\{{\mathbf{\epsilon }}_{\lambda }(\mathbf{k}){\hat{b}}_{\lambda ,\mathbf{k}}{A}_{\lambda }(t,\mathbf{k})e^{i\mathbf{k}\cdot \mathbf{x}}+{\mathbf{\epsilon }}_{\lambda }^{\ast }(\mathbf{k}){\hat{b}}_{\lambda ,\mathbf{k}}^{†}{A}_{\lambda }^{\ast }(t,\mathbf{k})e^{-i\mathbf{k}\cdot \mathbf{x}}\right\},$$

$$\mathbf{k}\cdot {\mathbf{\epsilon }}_{\lambda }(\mathbf{k})=0,{\mathbf{\epsilon }}_{\lambda }^{\ast }(\mathbf{k})={\mathbf{\epsilon }}_{-\lambda }(\mathbf{k}),[i\mathbf{k}\times {\mathbf{\epsilon }}_{\lambda }(\mathbf{k})]=\lambda k{\mathbf{\epsilon }}_{\lambda }(\mathbf{k}),{\mathbf{\epsilon }}_{\lambda }^{\ast }(\mathbf{k})\cdot {\mathbf{\epsilon }}_{{\lambda }^{\prime }}(\mathbf{k})={\delta }_{\lambda {\lambda }^{\prime }}.$$

$$[{\hat{b}}_{\lambda ,\mathbf{k}},{\hat{b}}_{{\lambda }^{\prime },{\mathbf{k}}^{\prime }}^{†}]={\delta }_{\lambda {\lambda }^{\prime }}{\delta }^{(3)}(\mathbf{k}-{\mathbf{k}}^{\prime }).$$

$${\ddot{A}}_{\lambda }(t,\mathbf{k})+H{\dot{A}}_{\lambda }(t,\mathbf{k})+\left[\frac{k^2}{a^2}-\lambda \frac{k}{a}2\xi H\right]A_{\lambda }(t,\mathbf{k})=0$$

Full text

magnetogenesis, pseudoscalar inflation, backreaction, Schwinger effect

I Introduction

Recently magnetic fields with very large coherence scale $\lambda_{B}\gtrsim 1\,$Mpc were detected in cosmic voids through the gamma-ray observations of distant blazars Neronov:2010 ; Tavecchio:2010 ; Taylor:2011 ; Dermer:2011 ; Caprini:2015 . Combining the corresponding results with the data of the cosmic microwave background (CMB) Planck:2015-pmf ; Sutton:2017 ; Jedamzik:2018 ; Paoletti:2018 ; Giovannini:2018b and ultra-high-energy cosmic rays Bray:2018 observations constrains the strength of these fields to $10^{-18}\lesssim B\lesssim 10^{-9}\,$G. It is well known also that magnetic fields are present on much smaller scales in the Universe Kronberg:1994 ; Grasso:2001 ; Widrow:2002 ; Giovannini:2004 ; Kandus:2011 ; Durrer:2013 ; Subramanian:2016 ; Giovannini:2018b including galaxies and clusters of galaxies. Observed intergalactic magnetic fields may be of either astrophysical or primordial origin and both magnetogenesis scenarios are currently under active consideration and study. Although astrophysical mechanisms based on the Biermann battery Biermann:1950 ; Pudritz:1989 ; Gnedin:2000 have been proposed to generate the “seed” magnetic fields and various types of dynamo can enhance them Zeldovich:1980book ; Lesch:1995 ; Kulsrud:1997 ; Colgate:2001 , it is problematic to use these mechanisms for the generation of magnetic fields with very large correlation length in cosmic voids indicating on the primordial origin of these magnetic fields. Clearly, this opens an intriguing possibility to obtain important information on the physical processes in the early Universe and poses the problem of how these fields were generated.

One of the natural mechanisms of primordial magnetic field generation is connected to the phase transitions in the early Universe. However, although they may lead to magnetic fields of the necessary strength Hogan:1983 ; Quashnock:1989 ; Vachaspati:1991 ; Cheng:1994 ; Sigl:1997 ; Ahonen:1998 , the coherence length is determined by the horizon size at the moment of phase transition and is much less than Mpc today. On the other hand, the inflationary magnetogenesis Turner:1988 ; Ratra:1992 can provide the necessary seeds for the observed magnetic fields and attain very large coherence length.

In order to generate electromagnetic (EM) fields during inflation, the conformal invariance of Maxwell’s action should be broken because fluctuations of the EM field are not enhanced in conformally flat inflationary background Parker:1968 . The conformal invariance can be broken by introducing the interaction of EM fields with scalar or pseudoscalar inflaton field or with the curvature scalar (see the pioneering works of Refs. Turner:1988 ; Ratra:1992 ; Garretson:1992 ; Dolgov:1993 ).

Inflationary magnetogenesis typically faces a backreaction problem which can put strong constraints on the parameter space of the model Demozzi:2009 . Usually in the study of inflationary magnetogenesis, one solves first the equations governing the evolution of the inflaton field and the scale factor and then determines generated EM fields. However, such an approach is possible only if generated fields are weak enough not to influence the background evolution. In the opposite case, one has to take into account generated EM fields self-consistently and solve all equations of motion simultaneously. Clearly, this makes the problem of inflationary magnetogenesis much more complicated. Moreover, the backreaction of EM fields can modify the inflaton dynamics, change the Universe expansion, and make an impact on the spectrum of primordial perturbations. In our previous work Sobol:2018 , we considered the case of kinetic coupling between the scalar inflaton field and the EM field introduced by Ratra Ratra:1992 and revisited many times in the literature Giovannini:2001 ; Bamba:2004 ; Martin:2008 ; Demozzi:2009 ; Kanno:2009 ; Ferreira:2013 ; Ferreira:2014 ; Vilchinskii:2017 ; Sharma:2017b ; Shtanov:2018 and showed that the backreaction of generated electric fields is indeed very important because it strongly suppresses magnetogenesis.

Because of the Schwinger effect Schwinger:1951 the generated electric fields can spontaneously create from vacuum the pairs of charged particles. The simplest case of pair creation in a constant and homogeneous background electric field in de Sitter space-time was investigated by many authors  Kobayashi:2014 ; Froeb:2014 ; Bavarsad:2016 ; Stahl:2016a ; Stahl:2016b ; Hayashinaka:2016a ; Hayashinaka:2016b ; Sharma:2017 ; Tangarife:2017 ; Bavarsad:2018 ; Hayashinaka:2018 ; Hayashinaka:thesis ; Stahl:2018 ; Geng:2018 ; Bavarsad:2018 , where different space-time dimensions and the various types of charge carriers were considered. However, the case of a constant electric energy density is unphysical because maintaining this regime would require the existence of ad hoc currents that could violate the second law of thermodynamics Giovannini:2018a . Nevertheless, the expressions for the Schwinger current derived in this approximation have the same functional dependence as in a time-dependent electric background in the strong-field regime Kitamoto:2018 .

It should be noted that the cosmological Schwinger effect contains interesting features Kobayashi:2014 ; Hayashinaka:2016a ; Hayashinaka:2016b ; Hayashinaka:2018 ; Hayashinaka:thesis that are absent in its flat-space counterpart such as (i) the infrared hyperconductivity in the bosonic case when the conductivity becomes very large in the limit of a small mass of charged particles and (ii) the negative conductivity in the weak field regime $eE\ll H^{2}$ which can, in principle, lead to the enhancement of the electric field Stahl:2018 . However, the definition of a finite part of the induced current is ambiguous and depends on the subtraction scheme. Thus, the negative conductivity is rather speculative Hayashinaka:2018 ; Hayashinaka:thesis and may be an artifact of the used subtraction scheme. Moreover, in Ref. Banyeres:2018 it was shown that the negative contribution to the conductivity is due to the nonlinear corrections to the Maxwell action and the logarithmic running of the coupling constant. Therefore, it cannot lead to any enhancement.

In this paper we consider the axial coupling of the EM field to a pseudoscalar inflaton field via the term $I(\phi)F_{\mu\nu}\tilde{F}^{\mu\nu}$ (here $\tilde{F}^{\mu\nu}$ is the dual EM tensor) taking into account the Schwinger effect and the backreaction of generated EM fields on the background evolution. In the case of the kinetic coupling model Kanno:2009 ; Sobol:2018 , the electric component dominates the EM energy density causing strong backreaction and spoiling magnetogenesis. Therefore, the kinetic coupling model is unfavorable in contrast to pseudoscalar inflation Durrer:2011 . Magnetogenesis in the latter model is promising because it generates helical magnetic fields Garretson:1992 ; Durrer:2011 ; Anber:2006 ; Anber:2010 ; Barnaby:2012 ; Caprini:2014 ; Anber:2015 ; Ng:2015 ; Adshead:2015 ; Adshead:2016 ; Shtanov:2019 which are more stable against dissipation and, in addition, their coherence length can be further increased due to the inverse cascade process Joyce:1997 ; Cornwall:1997 ; Boyarsky:2012 ; Sydorenko:2016 ; Boyarsky:2015 ; Gorbar:2016a ; Gorbar:2016b in the chiral primordial plasma. Since electric fields are also generated in this model, the pair production due to the Schwinger effect may become important and affect magnetogenesis. The same applies to the backreaction problem due to generated EM fields Garretson:1992 ; Durrer:2011 ; Anber:2006 ; Anber:2010 ; Barnaby:2012 ; Caprini:2014 ; Anber:2015 ; Ng:2015 ; Adshead:2015 ; Adshead:2016 ; Fujita:2015 ; Domcke:2018 ; Notari:2016 ; Figueroa:2018 ; Durrer:2019 . This means that both the Schwinger effect and backreaction of EM fields should be taken into account in the study of magnetogenesis in pseudoscalar inflation.

Since the pseudoscalar inflation model looks promising (pseudoscalar field may appear also as a Kalb-Ramond axion in string theory inspired models Basilakos:2019a ; Basilakos:2019b ), this provides the main motivation to extend our approach developed in Sobol:2018 to the case of pseudoscalar inflation magnetogenesis. Since it is quite costly computationally to solve directly the equations for the EM field for each mode, our approach in Ref. Sobol:2018 was connected with the derivation of a self-consistent set of equations for the electric field energy density, scale factor, and inflaton field taking into account the pair production due to the Schwinger effect. Perhaps not surprisingly, the implementation of this approach for magnetogenesis in the pseudoscalar inflation model turned out to be quite difficult. First of all, since the axial coupling involves the scalar product $\mathbf{E}\cdot\mathbf{B}$ of electric and magnetic fields, it is immediately clear that the function $\mathbf{E}\cdot\mathbf{B}$ should be necessarily taken into account in addition to the electric field energy density $\rho_{E}=\mathbf{E}^{2}/2$. Unfortunately, this strongly complicates the analysis because then the self-consistency of equations for the EM field requires the introduction of new quadratic functions of electric and magnetic fields with additional spatial curls. Therefore, the corresponding system of equations governing the evolution of the EM field is not closed, and one has to use certain physically reasonable approximations which make it possible to close the system.

This paper is organized as follows. We formulate a self-consistent system of equations that govern the joint evolution of the scale factor, inflaton field, and quadratic functions of the EM field in Sec. II, taking into account the backreaction of generated fields and the Schwinger effect. In Sec. III, we derive the explicit form of the boundary terms in the Maxwell equations which describe the quantum-to-classical transition for modes undergoing the tachyonic instability during inflation. Section IV is devoted to the determination of initial conditions for the EM field which should be imposed far from the end of inflation in order to study the system numerically. In Sec. V we analyze the conditions when the Schwinger effect becomes important and affects magnetogenesis. Section VI explains the approximation scheme which allows one to truncate the infinite system of equations and contains the numerical results in the $\alpha$-attractor inflationary model. The postinflationary evolution of generated magnetic fields is described in Sec. VII, where the present day values of the magnetic field strength and correlation length are determined. Finally, the summary of our work is given in Sec. VIII. Some technical details of the derivation of the initial conditions are given in Appendix A. We use in this paper the natural units and set $\hbar=c=1$.

II Evolution equations

We consider a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) Universe with metric

[TABLE]

In order to study magnetogenesis in pseudoscalar inflation, we use the following action for the inflaton field $\phi$ and the EM 4-potential $A_{\mu}$:

[TABLE]

where $V(\phi)$ is the effective potential of the inflaton, $I(\phi)$ is a coupling function, and $\tilde{F}^{\mu\nu}=\frac{1}{2}\eta^{\mu\nu\lambda\rho}F_{\lambda\rho}$ is the dual EM tensor. In addition, $\eta^{\mu\nu\lambda\rho}=\varepsilon^{\mu\nu\lambda\rho}/\sqrt{-g}$, where $\varepsilon^{\mu\nu\lambda\rho}$ is the totally antisymmetric Levi-Civita symbol with $\varepsilon^{0123}=+1$ and $\mathcal{L}_{\rm charged}(A_{\nu},\chi)$ is a gauge-invariant Lagrangian of a charged field $\chi$ interacting with the EM field $A_{\mu}$. The charged field $\chi$ is a generic field which could stand for a boson or fermion field as well as for any collection of them. It can, in principle, interact also with the inflaton directly. This interaction is crucial during the preheating stage because it determines the production of charged particles by the quickly oscillating inflaton. However, particle production is not effective during the slow-roll inflation because the inflaton field changes adiabatically slowly. Since the preheating stage lies beyond the scope of our paper, we will neglect the interaction between $\chi$ and $\phi$.

Since $F_{\mu\nu}\tilde{F}^{\mu\nu}$ is a pseudoscalar, in order to preserve the parity symmetry the coupling function of the pseudoscalar inflaton field $\phi$ with $F_{\mu\nu}\tilde{F}^{\mu\nu}$ should be odd, $I(-\phi)=-I(\phi)$. In this paper, we consider the simplest axial coupling function linear in $\phi$

[TABLE]

where $M_{p}=(8\pi G)^{-1/2}\approx 2.4\times 10^{18}\,{\rm GeV}$ is the reduced Planck mass and $\beta$ is a dimensionless coupling constant.

The Euler-Lagrange equations corresponding to action (2) and the Bianchi identity have the form

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

is the electric four-current. The right-hand side of Eq. (4) could be neglected if the generated fields are weak. If they are not, then the backreaction of these fields on the background evolution should be taken into account and a self-consistent system of equations should be solved.

According to the first Friedmann equation, the expansion rate of the FLRW Universe is described by the 00 component of the stress-energy tensor. In the case under consideration, the latter is given by

[TABLE]

where $\mathcal{L}_{0}=\tfrac{1}{2}(\partial\phi)^{2}-V(\phi)-\tfrac{1}{4}F_{\mu\nu}F^{\mu\nu}$ and $T_{\mu\nu}^{({\rm charged})}$ describes the contribution due to the charged field $\chi$. Notice also that the axial coupling term does not contribute to the energy-momentum tensor because it does not depend on metric.

It is convenient to use the Coulomb gauge for the EM field $A_{\mu}=(0,\,\mathbf{A})$ where ${\rm div\,}\mathbf{A}=0$. Then electric and magnetic fields are defined as follows:

[TABLE]

where $a(t)$ is a scale factor of the FLRW background. The components of the EM tensors are expressed through the electric and magnetic fields as usual

[TABLE]

where $\varepsilon_{ijk}$ is the three-dimensional Levi-Civita symbol and the Latin indices $i,\,j,\,k$ on the right-hand side denote the components of 3-vectors. We would like to note that the electric and magnetic fields defined in Eqs. (9) and (10) are physical fields measured by a comoving observer.

We will assume that the inflaton field is spatially homogeneous (the question whether perturbations of metric and the inflaton can always be neglected in the study of magnetogenesis is important and deserves a separate investigation). Then the energy density reads as

[TABLE]

where $\rho_{\rm inf}$ and $\rho_{\chi}$ are the energy densities of the inflaton and charged particles produced by the Schwinger effect, respectively. Since EM fields are generated through the enhancement of quantum fluctuations, their energy density $\rho_{\rm EM}=\frac{1}{2}\langle E^{2}+B^{2}\rangle$ is given by the expectation value of the corresponding operator.

The Friedmann, Klein-Gordon-Fock (KGF), and Maxwell equations define the following closed system of dynamic equations:

[TABLE]

We express the current $\mathbf{j}$ of charged particles in terms of the generalized conductivity $\sigma$

[TABLE]

The electric and magnetic fields are vector quantities and, in addition to the absolute value, are characterized by the direction. Since these fields are generated from quantum fluctuations, they are chaotically oriented. Operating with such vector objects is inconvenient. Even more crucially, solving the system of equations (12)–(16) requires determining the dependence of electric and magnetic fields on spatial coordinates that makes numerical calculations very demanding. Therefore, we generalize our approach in Ref.Sobol:2018 to the case of the pseudoscalar inflation model and introduce the following expectation values of scalar quadratic functions of electric and magnetic fields:

[TABLE]

where $n=0,1,2,...$. Note that we consider quadratic correlators of electric and magnetic fields and spatial derivatives of the latter through curls because the Maxwell equations (14) and (15) contain spatial derivatives of electric and magnetic fields only in the form of curls. As we will see below, the functions with $n>1$ are considered because the evolution equations for the expectation values with $n=0$ and $n=1$ do not form a closed system. Finally, we would like to note that the position of curl acting on the first or the second multiplier does not matter in view of the identity

[TABLE]

because the expectation value $\langle[\mathbf{f}\times\mathbf{g}]\rangle$ does not depend on the spatial coordinates.

In fact, functions (18)–(20) for $n=0$ are related to the well-known characteristics of the EM field. Indeed, $\mathscr{E}^{(0)}/2=\rho_{E}$ and $\mathscr{B}^{(0)}/2=\rho_{B}$ are the electric and magnetic field energy densities, respectively. Further, $\mathscr{G}^{(0)}=-\langle\mathbf{E}\cdot\mathbf{B}\rangle$ is an important quantity which depends on helical properties of the EM field and determines its backreaction on the inflaton. Moreover, the quantity $\mathscr{B}^{(-1)}$ defines the magnetic helicity [we treat the quantity ${\rm rot}^{-1}\mathbf{B}$ in Eq. (20) for $n=-1$ simply as $a^{-2}\mathbf{A}$ in accord with Eq. (9)]

[TABLE]

which is crucial for the postinflationary evolution of the magnetic field.

The Maxwell equations (14) and (15) relate the time derivatives of the electric and magnetic fields to their spatial derivatives (curls). Therefore, the evolution equation for any function (18)–(20) always contains a function with one additional spatial curl. This results in an infinite chain of equations for the functions of different order. Indeed, by using Eqs. (14) and (15), we find

[TABLE]

We introduced also extra terms on the right-hand side of the above equations. The reason for their appearance is the following. Electromagnetic fields are generated from quantum fluctuations and can be treated classically only after the corresponding modes exit the horizon Lyth:2008 and undergo the tachyonic instability [see the corresponding discussion after Eq. (39)]. Therefore, the number of modes, which undergo enhancement, changes as the inflation goes on. This leads to additional (boundary) terms in the equations of motion (22)–(24).

Adding Eqs. (22) and (24) for $n=0$ and dividing them by 2, we obtain the equation governing the evolution of the EM energy density

[TABLE]

The term $2\sigma\rho_{E}$ on the left-hand side of Eq. (25) describes the dissipation of the EM energy density due to the Schwinger effect. It will be considered in more detail in Sec. V. The first term on the right-hand side corresponds to the energy density transfer from the inflaton to the EM field due to the axial coupling between them. Obviously, the energy transfer term enters the equation for the inflaton energy density evolution with the opposite sign [this equation can be obtained from Eq. (13) multiplying it by $\dot{\phi}$]

[TABLE]

where $\rho_{\rm inf}=\dot{\phi}^{2}/2+V(\phi)$ and $p_{\rm inf}=\dot{\phi}^{2}/2-V(\phi)$.

In contrast to the case of the kinetic coupling model studied in Ref. Sobol:2018 , Eq. (25) for the evolution of the EM energy density and the KGF equation (13) do not form a closed system of equations because they contain a new quantity $\mathscr{G}^{(0)}=-\langle\mathbf{E}\cdot\mathbf{B}\rangle$. Equation (23) for $\mathscr{G}^{(0)}$, in turn, contains the higher-order functions $\mathscr{E}^{(1)}$ and $\mathscr{B}^{(1)}$. Clearly, this generates an infinite chain of equations and we need to impose some additional approximations in order to terminate it. Our approximation scheme will be discussed in Sec. VI. As to the boundary terms on the right-hand side of Eqs.(22)–(24), their explicit expressions will be considered in Sec. III.

Before closing this section, let us estimate the value of the EM energy density when the backreaction becomes relevant. For this purpose, we consider Eq. (13) with the coupling function (3) and compare the potential term $V^{\prime}(\phi)$ which determines the slow rolling of the inflaton with the term on the right-hand side of the KGF equation. Previous numerical studies of magnetogenesis in pseudoscalar inflation Fujita:2015 ; Figueroa:2018 found that the electric and magnetic fields make comparable contributions to the energy density, i.e., $\rho_{E}\sim\rho_{B}$ (this is in contrast to the kinetic coupling model where the electric component strongly dominates over the magnetic one, see, e.g., Refs. Martin:2008 ; Demozzi:2009 ; Kanno:2009 ; Vilchinskii:2017 ; Sobol:2018 ). Therefore, we can estimate the scalar product $\mathbf{E}\cdot\mathbf{B}$ as follows:

[TABLE]

Consequently, the backreaction becomes important when the EM energy density reaches the value

[TABLE]

where

[TABLE]

is the slow-roll parameter and we used the fact that $\rho_{\rm inf}\approx V(\phi)$ in the slow-roll regime $\epsilon\ll 1$. Since the typical value of the coupling constant, for which the backreaction can take place during inflation, equals $\beta\sim 10$ (see, e.g., Refs. Fujita:2015 ; Figueroa:2018 ), we conclude that the corresponding EM energy density is a few orders of magnitude less than that of the inflaton. This means that generated EM fields affect primarily only the evolution of the inflaton field while contribute to the Friedmann equation only indirectly, through the inflaton. It is worth noting that the same picture holds also for the kinetic coupling model Sobol:2018 ; Kanno:2009 .

III Boundary terms

In this section, we determine the boundary terms for the Maxwell equations (22)–(24), which describe the quantum-to-classical transition for the modes crossing the horizon. For this, we should consider the quantized EM field. There is, however, one complication connected with the pair production of charged particles by the classical electric field due to the Schwinger effect. The current of these particles backreacts on the EM field evolution making it, in general, nonlinear. However, in order to define the boundary terms, we will neglect the Schwinger effect. This is possible because the boundary terms are important during the early stages of inflation when the EM field is weak and each new mode crossing the horizon makes a notable contribution. Later, when the EM field is strong and a lot of modes are outside the horizon, the contribution due to new modes crossing the horizon is negligible. Vice versa, the Schwinger effect is negligible at the beginning of inflation and becomes important only close to the end of inflation when the generated electric field becomes strong. Therefore, the boundary term could be derived neglecting the Schwinger effect and the Schwinger current can be added phenomenologically in the final system of equations.

In the Coulomb gauge, the EM vector-potential operator can be decomposed over a set of creation and annihilation operators of two transverse polarizations

[TABLE]

where $\bm{\varepsilon}_{\lambda}(\mathbf{k})$ are two circular polarization 3-vectors, which satisfy the following conditions:

[TABLE]

The creation and annihilation operators have the standard commutation relations

[TABLE]

Substituting decomposition (30) into Eq. (14) we find the following equation governing the evolution of the mode function:

[TABLE]

or in conformal time $\eta=\int^{t}dt^{\prime}/a(t^{\prime})$,

[TABLE]

where we introduced the parameter

[TABLE]

In order to define the unique solution of Eqs. (33) or (34), we chose the initial conditions for the mode function as those in the Bunch-Davies vacuum Bunch:1978 (see the corresponding discussion in Sec.IV)

[TABLE]

Inspecting Eq. (33), we note that the tachyonic instability takes place for the modes with helicity $\lambda={\rm sign}(\xi)$ and $k/(aH)\lesssim 2|\xi|$ and one can expect the exponential amplification of these modes. As was mentioned above, the number of such modes changes in time. To derive the explicit expressions of the boundary terms we need to express Eqs. (18)–(20) in terms of the Fourier modes of the EM field. Using decomposition (30) and taking into account that only one polarization with $\lambda={\rm sign}(\xi)$ undergoes enhancement, we obtain

[TABLE]

where $k_{c}=2|\xi|aH$ is the critical momentum for which the instability occurs. The upper integration limit in Eqs. (37)–(39) is defined by $k_{c}$ because subhorizon modes have an oscillatory behavior and should be regarded as quantum fluctuations. Therefore, such modes do not contribute to the classical observables (18)–(20) and are excluded from the integration. Throughout the paper, we refer to a mode as subhorizon if its momentum is larger than the critical one $k>k_{c}$ [its wavelength is smaller than the effective horizon size $1/(2|\xi|H)$] and superhorizon in the opposite case.

In order to derive the boundary term, we note that Eq. (34) by using the replacements

[TABLE]

can be rewritten in the form of the Schrödinger equation

[TABLE]

For modes with momentum $k\ll k_{c}$, one can apply the Wentzel-Kramers-Brillouin (WKB) approximation and find the following approximate solution of Eq. (34) with the boundary conditions (36):

[TABLE]

where $\frac{p}{\hbar}=\sqrt{k^{2}-2aH|\xi|k}$. This solution is applicable when

[TABLE]

Here only the scale factor $a$ has been differentiated as the fastest varying quantity.

By using $k=2|\xi|aH(1+\delta)$, we obtain that the semiclassical solution (41) is applicable for

[TABLE]

Clearly, the semiclassical approximation fails when the left-hand and right-hand sides in the above inequality are approximately equal. A solution of the corresponding equation can be easily found in the two limiting cases $|\xi|\gg 1$ ($\delta_{c}\simeq[4|\xi|]^{-2/3}$) and $|\xi|\ll 1$ ($\delta_{c}\simeq[4|\xi|]^{-1/2}$) and then fit to the exact solution with accuracy better than 3% for all values of $\xi$ by the following function:

[TABLE]

As usual, the WKB method is inapplicable near the “classical turning point” $k=k_{c}$. Still we can approximate the value of $A_{\lambda}(t,k_{c})$ by its value at the boundary of validity of the semiclassical approximation at a slightly larger momentum $k_{c}(1+\delta_{c})$. Using Eqs. (41), (43), and (44), we have

[TABLE]

where $Q(|\xi|)=4|\xi|\delta_{c}\approx(4|\xi|)^{1/2}\left[1+(4|\xi|)^{2/3}\right]^{-1/4}$. Then the boundary terms equal

[TABLE]

Finally, let us comment on the magnetic helicity and its connection with the correlation length of the magnetic field. The energy density and helicity of the magnetic field are given by Durrer:2013 ; Subramanian:2016

[TABLE]

and its correlation length equals

[TABLE]

Comparing Eqs. (52) and (53), we derive the so-called realizability condition:

[TABLE]

This inequality is correct for any configuration of the magnetic field. However, in the case of maximally helical magnetic field, where only one polarization is enhanced, it approaches the equality

[TABLE]

which can be used to determine the correlation length of the magnetic field.

IV Initial conditions

It is well known that in order to resolve several cosmological problems, such as the flatness and horizon problems, one needs at least 50–60 $e$-foldings of inflation. Moreover, the primordial power spectrum inspected by the CMB observations is generated also at 50–60 $e$-foldings before the end of inflation. Therefore, although inflation can last more $e$-foldings, it is natural in numerical analysis to take the starting point of the simulation at 60 $e$-foldings before the end of inflation. Then, in order to find a solution to Eqs. (22)–(24), one should specify the initial conditions for the corresponding functions.

As we discussed above, modes deep inside the horizon oscillate and correspond to solution (36) in the Bunch-Davies vacuum. In this regime, they represent quantum fluctuations and do not contribute to the classical observables. Later due to the Universe expansion, the physical momentum of a given mode decreases and the square of the effective frequency changes its sign causing a tachyonic instability. The corresponding mode function stops to oscillate and can be interpreted as a classical field. Therefore, in order to impose the initial conditions, we should take into account all modes which underwent the tachyonic instability before the initial moment of time.

Let us consider how the parameter $\xi$ evolves during the slow-roll inflation in the absence of the backreaction of generated EM fields. Using the KGF equation (13) for the inflaton and the Friedmann equation (12) in the slow-roll regime, we find

[TABLE]

where $\epsilon_{1}$ is the first Hubble-flow function which coincides with the slow-roll parameter $\epsilon$ [see Eq. (29)] in the slow-roll regime Martin:2014 . The rate of time evolution of $\xi$ is given by

[TABLE]

where $\epsilon_{2}$ is the second Hubble-flow function from an infinite sequence of such quantities Martin:2014

[TABLE]

and $H_{i}$ is the initial value of the Hubble parameter, which is used as a reference point.

For the axionlike coupling (3), only the second term in Eq. (57) remains. Then we conclude that during inflation $\xi$ changes adiabatically slowly and can be considered as a constant for a few Hubble times. On the other hand, over long periods of time it slightly changes and its value is determined by the current values of $\dot{\phi}$ and $H$. In other words, it could be taken as a constant parameter in the Maxwell equations, therefore, the generation of EM fields is governed only by the value of $\xi$ at a given moment of time. We can also approximate the space-time metric by the de Sitter one assuming that the scale factor grows strictly exponentially

[TABLE]

Clearly, this approximation is valid with high accuracy in the slow-roll regime far from the end of inflation. Then Eq. (34) takes the form

[TABLE]

Since the conformal time during the inflation stage has the negative sign, the second term in brackets in Eq. (60) for $\lambda={\rm sign}(\xi)$ is negative and, for a certain range of momenta $k/(aH)<2|\xi|$, it becomes larger than the first term and causes the tachyonic instability for the corresponding modes. This instability is absent for modes with another polarization, therefore, generated EM fields are helical. For $|\xi|\gtrsim 1$, the amplification is very effective and the corresponding fields can be considered as maximally helical, i.e., one can simply neglect the second circular polarization as we do below.

According to Ref. Anber:2010 , the solution of Eq. (60) can be found in terms of Coulomb wave functions Luke:book . However, it is more convenient to use the approximate expression for the growing mode for $k/(aH)<2|\xi|$ found in Ref. Durrer:2011

[TABLE]

where $W(|\xi|)=[e^{\pi|\xi|}\,{\rm sh\,}\pi|\xi|/(\pi|\xi|)]^{1/2}$, $z=2\sqrt{2|\xi|}\sqrt{k/(aH)}$, and $K_{\nu}(z)$ is the Macdonald function. Substituting these expressions into Eqs. (37)–(39), we obtain the initial conditions for the corresponding expectation values as functions of the parameter $\xi_{0}$ at the initial moment of time:

[TABLE]

[TABLE]

[TABLE]

These integrals can be expressed in terms of the Meijer $G$ function and the corresponding results are listed in Appendix A [see Eqs. (100)–(101)] together with their asymptotic expressions in the two limiting cases $|\xi_{0}|\ll 1$ and $|\xi_{0}|\gg 1$.

V Schwinger effect

In this section, we determine the Schwinger conductivity $\sigma$ and the energy density of charged particles $\rho_{\chi}$ created via the Schwinger process. For the Schwinger current, we will use the expressions derived in the minimal subtraction scheme in Refs. Kobayashi:2014 ; Hayashinaka:2016a . Since the general expressions are rather cumbersome, it is more convenient to use their asymptotics governed by the parameters

[TABLE]

where $m$ and $e$ are the mass and electric charge of created particles, respectively.

The typical value of the Hubble parameter, which can be fixed from the observations of the amplitude of primordial perturbations Planck:2018-infl , is $H\sim 10^{-5}\,M_{p}\sim 10^{13}\,{\rm GeV}$. Since the lightest charged spin-$1/2$ particle, the electron, has a mass $m\sim 10^{-3}\,{\rm GeV}$ and the lightest charged scalar particle, the pion, has a mass $m\sim 0.1\,{\rm GeV}$, we are interested in the small mass limit $M\ll 1$.

At the beginning of inflation, the parameter $\xi$ is usually small because it is proportional to the square root of the slow-roll parameter $\epsilon$, see Eq. (56). Thus, we may use Eq. (103) to estimate the initial value of the electric energy density. Then the parameter $L$ equals

[TABLE]

where $b=\exp(\gamma_{E})$ and $\gamma_{E}=0.577...$ is the Euler-Mascheroni constant. For a general model of inflation consistent with the CMB observations, $\epsilon$ at 50–60 $e$-foldings before the end of inflation can be expressed in terms of the tensor-to-scalar power ratio as $\epsilon=r/16\lesssim 10^{-3}$ (according to the latest CMB observations by Planck Collaboration and BICEP- Keck Array Planck:2018-infl , the upper limit on $r$ is $r<0.064$). Then, for the coupling constant $\beta\sim 10$, we have $|\xi_{0}|\sim 0.1$ and $L_{0}\sim 10^{-4}$.

On the other hand, close to the end of inflation and for the sufficiently large couplings that can cause the backreaction, the electric energy density is given by $\mathscr{E}^{(0)}/2=\rho_{E}\sim\sqrt{\epsilon}\beta^{-1}\rho_{\rm tot}=3\sqrt{\epsilon}\beta^{-1}H^{2}M_{p}^{2}$ [see estimate (28)]. Then $L\sim e\epsilon^{1/4}\beta^{-1/2}M_{p}/H\sim 10^{2}-10^{3}$ for $H\sim 10^{-5}\,M_{p}$ and $\beta\sim 10$. Therefore, the values of $L$ can vary in a wide range.

In the weak-field regime, $L\ll 1$, the fermionic Schwinger conductivity has the form Hayashinaka:2016a

[TABLE]

where $\alpha_{\rm EM}=e^{2}/(4\pi)\approx 1/137$ is the fine structure constant. It is easy to notice that the conductivity is negative in the low mass limit $m\ll H$ and could, in principle, enhance the electric field. Such a possibility was investigated in Ref. Stahl:2018 where it was found that the enhancement is negligible compared to the damping of the field due to the Universe expansion because the Schwinger conductivity is much less by the absolute value than the Hubble parameter for physically relevant particle masses. Moreover, the authors of Ref. Banyeres:2018 claim that the negative conductivity is spurious and corresponds to nonlinear corrections to the Maxwell action and the logarithmic running of the coupling constant. In any case, the ratio $\sigma/H$ is always small and the Schwinger effect is negligible in the small-field regime.

For bosons, the conductivity is positive and attains very large values in the low mass regime where it diverges like $1/M^{2}$ for $M\to 0$ Kobayashi:2014 ; Hayashinaka:2016b . In the literature this phenomenon is known as the infrared hyperconductivity. Since such a behavior is observed only in the extremely low field regime when $L\ll M\ll 1$, it is never realized in the problem under consideration because as we showed above $L$ cannot take extremely small values. For $L>M$, the conductivity quickly diminishes to values typical in the fermionic case (67) (see Ref. Hayashinaka:2016b ). Thus, the Schwinger effect does not play any role in the weak field regime.

The strong-field expressions ($L\gg 1$) for the Schwinger conductivity are similar for bosons and fermions and have the form

[TABLE]

where $g_{b}=1$ and $g_{f}=2$ are the number of spin degrees of freedom. Expressing the conductivity in terms of the energy density and neglecting the particles mass, we obtain

[TABLE]

Then we conclude that the Schwinger effect is important only in the strong field regime when $L>3\pi^{2}/\alpha_{\rm EM}\gg 1$. The corresponding condition for the electric energy density reads as

[TABLE]

Therefore, the strong-field expression (69) could be used all time. Although this expression is not correct in the weak field regime, this is irrelevant for the dynamics of EM fields.

Finally, we need also an equation governing the evolution of created particles. It is natural to require that the energy dissipated by EM fields is transferred into created particles. Then the energy density of produced particles can be described phenomenologically by the following equation:

[TABLE]

where the term on the right-hand side is exactly the same as in Eq. (25) for the electric energy density, however, has the opposite sign. Since we consider the particle mass to be smaller than the Hubble parameter, we treat the created particles as ultrarelativistic that gives the factor $4H$ in the above equation.

VI Numerical results

Having derived an infinite chain of Eqs. (22)–(24), we introduce in this section an approximation scheme which allows us to truncate this chain of equations and then numerically solve it in a certain inflationary model. The crucial physical input that allows us to perform this truncation is the assumption that electric and magnetic fields contribute equally to the energy density. This is a numerically well-justified approximation, especially for the small values of the coupling constant $\beta\lesssim 10$ (see, e.g., Fig. 6 in Ref. Figueroa:2018 ). The generalization of this approximation scheme to higher order quantities is straightforward. Indeed, we can assume that $\mathscr{E}^{(n)}=\mathscr{B}^{(n)}$ for a given $n$. This makes it possible to keep only 2 equations in the $n$th order. The first is the equation for $\mathscr{G}^{(n)}$, where $\mathscr{E}^{(n+1)}=\mathscr{B}^{(n+1)}$ is assumed and $\mathscr{B}^{(n)}$ is replaced by $\mathscr{E}^{(n)}$,

[TABLE]

The second equation for $(\mathscr{E}^{(n)}+\mathscr{B}^{(n)})/2\equiv\mathscr{E}^{(n)}$,

[TABLE]

is obtained as the sum of equations (22) and (24) divided by 2. Clearly, the two above equations together with the Friedmann equation (12) for the scale factor, the KGF equation (13) for the inflaton, and the equation for the helicity $\mathcal{H}$ form a closed system of evolution equations. Thus, in the approximation scheme of order $n$, we deal with $3n+5$ ordinary differential equations.

Below we provide numerical solutions of the system of equations discussed above for the inflationary model with concave and even potential. The latter property ensures the parity symmetry is preserved for the pseudoscalar inflaton field. The requirement for the inflaton potential to be concave follows from the recent CMB observations Planck:2018-infl . The best correspondence with observational data have the plateau models. The $T$ model of $\alpha$ attractors with hyperbolic tangent potential Kallosh:2013a ; Ferrara:2013 ; Kallosh:2013b

[TABLE]

is a typical representative of the plateau models. Here $\alpha$ is a free parameter. For given values of $\alpha$ and the number of $e$-foldings before the end of inflation when the pivot mode crosses the Hubble horizon $N_{\ast}$, one can determine the scalar spectral index $n_{s}$ and the tensor-to-scalar power ratio $r$

[TABLE]

Varying $\alpha$ and $N_{\ast}$ one can fill the region on the $(n_{s},\,r)$ diagram and compare the predictions of the model with the constraints set by the CMB observations of Planck Collaboration Planck:2018-infl . In particular, for $N_{\ast}=60$, we found that the range of $\alpha$ compatible with the 68% C.L. constraints is $\alpha<20$. For definiteness, we use $\alpha=1$ in the numerical analysis. For this value of the parameter $\alpha$, the potential has the form which is very similar to that of the well-known Starobinsky model, which is another typical plateau model.

Finally, we should determine the strength of the potential which gives the correct value of the scalar power spectrum amplitude in Eq. (84) below,

[TABLE]

Thus, fixing $N_{\ast}$ and $\alpha$, we fully determine the inflaton potential. Further, for a certain value of the coupling constant $\beta$, we calculate $\xi_{0}$ and use it to determine the initial conditions given by Eqs. (106)–(108). For definiteness, we consider the case of fermion charged particles created due to the Schwinger effect and then calculate their electric conductivity by using Eq. (69). Then, we apply the $n$th order approximation scheme and compute numerically the time dependence of the inflaton, scale factor, and EM quadratic functions (18)–(20).

Only relatively small values of the coupling constant $\beta\lesssim 10$ are physically allowed, because for larger values the reheating stage is spoiled and large non-Gaussianities in primordial spectra are generated Fujita:2015 ; Figueroa:2018 . Nevertheless, we will apply our method also for larger values of $\beta$ in order to study its applicability and to show more clearly some characteristic features of the magnetogenesis process. Figure 1 shows the time dependence of the inflaton field [panel (a)] and the EM energy density [panel (b)] during the inflation stage calculated in the approximation scheme with $n=1$ for three nonzero values of the coupling constant. The curves for the smallest value $\beta=10$ are shown by the red dashed lines. In this case, the inflaton field evolves in the same way as in the absence of the coupling $\beta=0$ (shown by the blue solid line), meaning that the backreaction of the generated EM field on the inflaton evolution is negligible.

On the other hand, for larger values of the coupling constant $\beta=20$ and $\beta=40$ shown by the green dashed-dotted line and the purple dotted line, respectively, the time dependence of the energy density at a certain moment of time begins to deviate from the unperturbed solution and changes from fast exponential growth to almost constant until the preheating stage when it starts to decay. As we showed in Sec. II, the critical value of the EM energy density when the backreaction becomes relevant is of the order $\rho_{\rm EM}\sim(\sqrt{2\epsilon}/\beta)\rho_{\rm inf}$. In addition, for larger couplings, the initial values and the boundary terms are larger too leading to the generation of stronger fields. Consequently, the backreaction takes place earlier for larger coupling constant. These features are easily seen in Fig. 1(b). The inflaton field also changes its behavior in the backreaction regime evolving more slowly and reaching the potential minimum where it oscillates with small amplitude. Slowing down of the inflaton rolling stops the generation of EM fields which attain constant values.

Now, let us study how the results depend on the order of the approximation scheme $n$ and how different components of the energy density evolve with time. The corresponding results are shown in Fig. 2 for $\beta=10$ for which the backreaction is irrelevant. Panel (a) compares different approximation schemes. The blue solid line corresponds to the exact solution of the mode equation (34) with the Bunch-Davies initial condition (36) in the assumption that the backreaction of the generated fields and the Schwinger effect do not affect the background evolution. The purple dotted line, green dashed-dotted line, and red dashed line show the EM energy density calculated using the approximation schemes with $n=0$ (equipartition between electric and magnetic energy densities), $n=1$, and $n=2$, respectively. First of all, it is very important that during the whole inflation stage all three used approximation schemes give results very close to the solution obtained by the numerical integration of Eq. (34) without taking into account the backreaction and Schwinger effect, see the main plot in Fig. 2(a). The inset shows the time interval close to the end of inflation in more detail. Inspecting it, we conclude that higher approximation schemes are closer to the exact solution. We should note, however, that the convergence is asymptotic and higher order approximation schemes with $n\geq 3$ deviate from the exact solution at late times. Perhaps this can be explained by a large number of equations ($3n+5$) which we have to deal with when applying the scheme with large $n$ as the computational errors increase rapidly with $n$.

Figure 2(b) shows the different components of the energy density calculated in the approximation scheme with $n=2$. The red dashed line gives the magnetic energy density, the green dashed-dotted line shows the electric energy density, the purple dotted line corresponds to the energy density of charged particles created due to the Schwinger effect. Finally, the blue solid line gives the total energy density. We see that the EM energy density is almost equally distributed between its electric and magnetic parts. This was observed previously in numerical lattice simulations performed in Refs. Fujita:2015 ; Figueroa:2018 . At the beginning of inflation, the evolution of the EM field is determined by the competition of two processes: the Hubble damping due to the Universe expansion [the second term on the left-hand side of Eqs. (22)–(24)] and enhancement due to new modes crossing the horizon [the boundary terms on the right-hand side of Eqs. (22)–(24)]. Since for $\xi\sim 1$ the magnetic boundary term is larger than the electric one, see Eqs. (48) and (50), and the Hubble damping has the same form, the “equilibrium” value of the magnetic energy density is larger than the electric one. However, close to the end of inflation, the electric energy density starts growing faster and crosses the magnetic energy density curve, see Fig. 2(b). This growth can be explained by a faster rolling of the inflaton and by the presence of the term $-2I^{\prime}(\phi)\dot{\phi}\mathscr{G}^{(0)}$ in Eq. (22) for $n=0$ which has the negative sign and, thus, increases the time derivative. Note that such a term is absent in Eq. (24) for the magnetic counterpart. The second crossing of these curves is explained by a fast decrease of the electric energy density due to the Schwinger effect coming into play at the end of inflation. The energy density of charged particles is negligibly small compared to that of the EM field during almost the whole inflation stage. Only close to the end of inflation, the electric field becomes strong enough to satisfy condition (70) whose lower bound is shown in Fig. 2(b) by the thin black dashed line. After that, the Schwinger effect produces the charged particles more intensively and their energy density quickly becomes comparable to that of the EM field.

Further, let us consider larger values of the coupling constant. The plots similar to Fig. 2 are shown in Fig. 3 for $\beta=20$. According to panel (a), the backreaction becomes important close to the end of inflation. Consequently, the EM energy density suddenly stops growing and deviates from the result obtained neglecting the backreaction (blue solid line). As in the case of small $\beta$, higher order approximation schemes are in better correspondence with the exact solution, however, only for small $n$. Finally, Fig. 3(b) shows that the Schwinger effect comes into play a bit earlier than for $\beta=10$ and produces much more charged particles. Close to the end of inflation, their energy density becomes larger than that of the EM field and, due to the decrease of the inflaton energy density, it finally becomes the dominant part of the energy density of the Universe. Thus, the Schwinger effect helps to fill the Universe with charged particles even before the stage of inflaton fast oscillations. In the literature, this is known as the “Schwinger preheating” scenario Tangarife:2017 ; Sobol:2018 .

We also studied how the Schwinger effect influences the electric field. Figure 4 shows the time dependence of the electric energy density for three values of the electric charge of created particles. The blue solid line corresponds to the case without the Schwinger effect ($e=0$). Here, after entering the backreaction regime, the electric energy density remains almost constant with small oscillations. For nonzero electric charge, the Schwinger effect turns on near the end of inflation and suppresses the electric energy density compared to the case $e=0$. This suppression is stronger and occurs earlier for larger values of the electric charge. Indeed, the lower bound for the electric energy density (70), when the Schwinger effect becomes important, depends on the electric charge as $e^{-6}$. Therefore, increasing the charge twice should suppress the electric field by $2^{6}=64$, which is in nice correspondence with the results in Fig. 4.

Finally, we calculated the magnetic energy density and the magnetic correlation length at the end of inflation. These quantities, determined in three approximation schemes with $n=0$, 1, and 2, are given in Fig. 5 as functions of the coupling constant $\beta$. Panel (a) shows that the magnetic energy density at first grows with $\beta$ and then starting from $\beta\sim 10$ approaches a constant value $\rho_{B}\sim 10^{-12}M_{p}^{4}$, which slightly depends on the applied approximation scheme. The latter behavior can be explained by the fact that for $\beta>10$ the backreaction regime occurs prior to the end of inflation and the growth of the EM field terminates. Thus, for large $\beta$, the backreaction imposes a universal upper bound for the magnetic energy density generated during inflation which does not depend on $\beta$. The correlation length $\lambda_{B}$ always has the same order of magnitude, namely, $\lambda_{B}\sim 10^{6}M_{p}^{-1}$, which is of the same order as the horizon at the end of inflation $\lambda_{H}=1/H_{e}$. For large values of $\beta$, all approximation schemes give almost the same results recovering the universal behavior in the backreaction regime. However, for small $\beta$, different approximation schemes produce the results differing by a factor $2$ to $5$. The correlation scale (53) is proportional to the ratio of $\mathscr{B}^{(0)}$ to $\mathscr{B}^{(-1)}$ and changes in a small range. Thus, even relatively small differences in $\mathscr{B}^{(n)}$ in different approximation schemes [see, e.g., Figs. 2(a) and 3(a)] lead to visible deviations in $\lambda_{B}$.

VII Postinflationary evolution of the magnetic field

Another important issue in the study of magnetogenesis is to determine how magnetic fields generated during inflation evolve through the sequence of cosmological epochs until the present. After the end of inflation, the preheating takes place when the inflaton field oscillates in the potential minimum and decays into different particles. The EM field could be amplified during the preheating stage due to the interaction with the rapidly oscillating inflaton, e.g., through the mechanism of parametric resonance Fujita:2015 ; Adshead:2015 ; Adshead:2016 ; Deskins:2013 . However, our method is inapplicable at this stage for several reasons: (i) the boundary terms were derived for the inflationary background when new modes cross the horizon, while during the preheating stage they reenter the horizon; (ii) the expressions for the Schwinger current were derived for the de Sitter-like background where $H$ is almost constant; (iii) finally, our approximation scheme assuming $\mathscr{E}^{(n)}=\mathscr{B}^{(n)}$ would break down because the Schwinger effect suppresses the electric component compared to the magnetic one. Moreover, the backreaction regime, which is favorable for obtaining large magnetic fields, slows down the inflaton evolution decreasing its oscillation amplitude and reducing the efficiency of the parametric resonance. Another source of the magnetic field amplification could be the power exchange with the electric field due to the Faraday law Kobayashi:2019 . However, this mechanism requires electric fields much stronger than the magnetic ones, which is not the case in the problem under consideration. Therefore, we will assume in what follows that the magnetic field is not enhanced during the preheating stage.

Further, during reheating, the plasma of created particles thermalizes and the Universe enters its hot radiation dominated phase. We will assume that at this stage the turbulent regime does not occur and the magnetic field and its correlation length evolve adiabatically, i.e., the corresponding comoving quantities remain constant:

[TABLE]

where the rescaling factor $a_{0}/a_{e}$ is determined below.

There are many models describing reheating in the literature, see Refs. Kofman:1997 ; Bassett:2006 ; Allahverdi:2010 ; Frolov:2010 ; Amin:2015 for a review. However, the lack of observational data does not allow one to discriminate between them. Nevertheless, this stage can be described phenomenologically by two parameters such as the temperature at the end of reheating $T_{\rm reh}$ and the effective equation of state $\bar{w}$ which is defined as follows Ringeval:2010 ; Martin:2014 :

[TABLE]

where $w(N)=p/\rho$ is the current value of the equation of state parameter and $N_{e}$ and $N_{\rm reh}$ are the numbers of $e$-foldings at the end of inflation and reheating, respectively. Using the energy conservation law

[TABLE]

we can relate the scale factors and energies at the end of inflation and reheating through $\bar{w}$ as

[TABLE]

The energy density at the end of inflation depends on the inflationary model and in the absence of the backreaction of EM fields can be expressed in terms of the inflaton potential

[TABLE]

Indeed, using Eq. (80) and the Hubble equation $H^{2}=\rho/(3M_{p}^{2})$, we find the following relation between pressure and the energy density:

[TABLE]

where $\epsilon_{1}$ is the first Hubble flow parameter, see Eq. (58). Taking into account that $\epsilon_{1}=1$ when inflation ends, we obtain $p_{e}=-\rho_{e}/3$ which implies Eq. (82). Moreover, Eq. (82) can be also used as an order-of-magnitude estimate of the energy density at the end of inflation even in the case when the backreaction occurs, because, as we see from the numerical results in, e.g., Fig. 3(b), the energy density is distributed between the inflaton, the EM field, and charged particles in comparable amounts.

Further, the CMB observations Planck:2018-infl provide constraints on the amplitude and spectral properties of the primordial perturbations which, in turn, can be expressed through the inflaton potential at the moment when the pivot scale $k_{\ast}$ exits the Hubble horizon (in what follows we mark all quantities at this moment by asterisk). In particular, the amplitude of the power spectrum of scalar perturbations equals

[TABLE]

The tensor-to-scalar power ratio depends on the first slow-roll parameter according to $r=16\epsilon_{\ast}$. Combining these two equations, we find the corresponding value of the inflaton potential

[TABLE]

Then the energy density at the end of inflation is given by

[TABLE]

where $A_{s}\simeq 2.1\times 10^{-9}$ is the spectral amplitude of scalar perturbations. The tensor-to-scalar power ratio is bounded from above by $r<0.064$ Planck:2018-infl and its typical value in the plateau models is of order $0.01$. We would like to note that $V_{e}/V_{\ast}\sim\mathcal{O}(0.1)$ in the majority of inflationary models except the case of the large-field inflation $V(\phi)\propto\phi^{n}$, which is strongly disfavored according to the results of the Planck mission Planck:2018-infl . Then the typical value of the energy density at the end of inflation is

[TABLE]

The energy density at the end of reheating can be expressed through the corresponding temperature as usual

[TABLE]

where $g_{\rm reh}$ is the effective number of relativistic degrees of freedom at that time. Using the entropy conservation relation, we can relate the scale factors at the end of reheating and the present time as follows:

[TABLE]

where $T_{0}=2.725\,{\rm K}=2.35\times 10^{-13}\,{\rm GeV}$ is the present CMB temperature and $g_{0}=2+3\times 2\times\frac{7}{8}\times\frac{4}{11}\approx 3.9$ is the effective number of relativistic degrees of freedom at present. Collecting all together, we obtain the rescaling factor

[TABLE]

which determines the comoving quantities (78) at the end of reheating.

Further, in the radiatively dominated epoch, the plasma is highly dense and conductive with extremely large kinetic and magnetic Reynolds numbers. In this turbulent regime the helical magnetic fields undergo the inverse cascade. The comoving helicity conservation and the equipartition of energy between the kinetic motion in plasma and the magnetic field imply Banerjee:2004 ; Durrer:2013 ; Subramanian:2016

[TABLE]

where $\eta$ is the conformal time. The spectrum of the magnetic field shifts to larger wavelength and lower magnitude and is self-similar all the time.

The inverse cascade takes place only in the turbulent regime which occurs when the Reynolds number is large. However, during the radiation domination epoch the turbulent phases alternate with the MHD viscous stages which include the periods of adiabatic evolution with $\tilde{\lambda}_{B}={\rm const}$ and the periods of rapid damping due to the free-streaming of the neutrinos or photons giving the major contribution to the viscosity. It was shown in Ref. Banerjee:2004 that the evolution during the viscous stage gives approximately such values for the comoving magnetic field and its correlation length as if the viscous stage were absent at all. In other words, for numerical estimates of the comoving quantities, it is possible to assume that the magnetic field from the end of reheating until recombination evolves in the inverse cascade regime.

Taking into account that $\eta\propto a\propto T^{-1}$ during the radiation dominated stage and $\eta\propto a^{1/2}\propto T^{-1/2}$ during the matter dominated epoch, we obtain the following comoving quantities at the moment of recombination:

[TABLE]

[TABLE]

where $T_{\rm eq}=6.54\times 10^{4}\Omega_{m0}h^{2}\,{\rm K}\simeq 9.16\times 10^{3}\,{\rm K}$ is the temperature at the matter-radiation equality and $T_{\rm rec}\simeq 4\times 10^{3}\,{\rm K}$ is the temperature at recombination.

If the correlation scale of the magnetic field becomes larger than the cosmic diffusion length $L_{\rm diff}\sim 1\,{\rm A.U.}$, then the magnetic field after recombination will undergo the free adiabatic evolution and Eqs. (92) and (93) give the present day values of the magnetic field and its correlation length. On the contrary, if the correlation length is smaller than $L_{\rm diff}$, the modes with shorter wavelengths than $L_{\rm diff}$ decay after recombination. In this case, in order to find the present day value of the magnetic field, we should know the spectral index $n_{B}$ which parametrizes the magnetic power spectrum $d\rho_{B}/d\,\ln k\propto k^{n_{B}}$. Assuming that $n_{B}>1$ (this is correct for magnetogenesis in pseudoscalar inflation, see Refs. Durrer:2011 ; Durrer:2013 ), we obtain the following present day value of the magnetic field and its coherence length for $\tilde{\lambda}_{B}(t_{\rm rec})<L_{\rm diff}$:

[TABLE]

We will see below that for realistic values of the inflaton energy density and reheating temperature the correlation length at recombination appears to be larger than the cosmic diffusion scale and rescaling (94) is not needed.

Finally, using Eqs. (78), (90), (92), and (93), we derive the general formulas for the present day magnetic field and its coherence length

[TABLE]

[TABLE]

In order to estimate the numerical values of the present day quantities (95)–(96), we consider two possible equations of state during reheating, namely, $\bar{w}=0$ if the preheating is described by decaying oscillations of the inflaton and $\bar{w}=1/3$ if the Schwinger effect is strong enough to transfer almost all energy density into ultrarelativistic charged particles or assuming that reheating is instantaneous. Choosing $10^{14}\,{\rm GeV}$ as a reference point for the reheating temperature, we obtain

[TABLE]

[TABLE]

Since the correlation scale is typically lower than $1\,{\rm Mpc}$, in order to compare our results with the data of gamma-ray observations Neronov:2010 ; Tavecchio:2010 ; Taylor:2011 ; Dermer:2011 ; Caprini:2015 , we find the following expression for the effective field $B_{\rm eff}=\tilde{B}\sqrt{\tilde{\lambda}_{B}/(1\,{\rm Mpc})}$ by using Eqs.(97) and (98):

[TABLE]

We determined in the previous section the numerical values of the magnetic energy density and correlation length at the end of inflation in the $\alpha$-attractor model. By using them we estimate now the present day values of the magnetic field and its correlation scale. The energy density at the end of inflation is of order $\rho_{e}\sim 10^{-10}M_{p}^{4}$. In the most optimistic case, the magnetic energy density is of order $\rho_{B}(t_{e})\sim 10^{-12}M_{p}^{4}$ and the correlation scale is always $\lambda_{B}(t_{e})\sim 10^{6}M_{p}^{-1}$. Then, for the equation of state with $\bar{w}=0$, we obtain $B_{0}\simeq 5\times 10^{-16}\,{\rm G}\left(\frac{T_{\rm reh}}{10^{14}\,{\rm GeV}}\right)^{1/3}$ and $\lambda_{B,0}\simeq 10^{-6}\,{\rm Mpc}\left(\frac{T_{\rm reh}}{10^{14}\,{\rm GeV}}\right)^{1/3}$. The effective magnitude of the magnetic field is equal to $B_{\rm eff}\simeq 5\times 10^{-19}\,{\rm G}\left(\frac{T_{\rm reh}}{10^{14}\,{\rm GeV}}\right)^{1/2}$. For $\bar{w}=1/3$, the present day quantities are $B_{0}\simeq 4\times 10^{-15}\,{\rm G}\left(\frac{T_{\rm reh}}{10^{14}\,{\rm GeV}}\right)^{-1/3}$, $\lambda_{B,0}\simeq 4\times 10^{-7}\,{\rm Mpc}\left(\frac{T_{\rm reh}}{10^{14}\,{\rm GeV}}\right)^{1/3}$, and $B_{\rm eff}\simeq 2.5\times 10^{-18}\,{\rm G}$.

Higher values of the magnetic field are generated for larger reheating temperature $T_{\rm reh}$. However, the latter cannot exceed $T_{\rm max}\sim 10^{15}\,{\rm GeV}$ for which all energy density of the inflaton $\rho_{e}$ instantaneously transfers into radiation during reheating. In general, lower values of $T_{\rm reh}$ lead to lower $B_{0}$ except for $\bar{w}=1/3$. However, in this case, the correlation scale decreases with $T_{\rm reh}$ so that the effective strength of the magnetic field on the Mpc scale always remains of order $10^{-18}\,{\rm G}$. Therefore, even the most optimistic estimate can hardly explain the strength of the intergalactic magnetic fields from the distant blazar observations Neronov:2010 ; Tavecchio:2010 ; Taylor:2011 ; Dermer:2011 ; Caprini:2015 .

VIII Conclusion

In this work we investigated the influence of the backreaction of the EM field and the Schwinger effect on magnetogenesis in the pseudoscalar inflation model with the axial coupling of the EM field to the inflaton via the term $\beta(\phi/M_{p})F_{\mu\nu}\tilde{F}^{\mu\nu}$. In order to perform a self-consistent study of this system, we extended the framework proposed in Ref. Sobol:2018 , where magnetogenesis in the kinetic coupling model was studied. Its essence is a closed set of ordinary differential equations that describe the self-consistent evolution of classical observables in the form of quadratic functions of the electric and magnetic fields with an arbitrary number of the curls. The quantum origin of these observables is incorporated in the form of boundary terms in the corresponding evolution equations. These terms come from the contribution of EM modes which cross the horizon and undergo the quantum-to-classical transition transforming from the Bunch-Davies vacuum state to a state with large filling number.

The evolution equations for the expectation values of quadratic functions of the electric and magnetic fields contain terms with an extra power of the spatial curl. This results in an infinite chain of equations for the corresponding classical observables. Therefore, some additional physical arguments should be used in order to terminate this system of equations to a finite set of equations. Since the coupling to the inflaton is symmetric with respect to the electric and magnetic fields, the physically reasonable approximation could be to set $\langle\mathbf{E}\cdot({\rm rot})^{n}\mathbf{E}\rangle=\langle\mathbf{B}\cdot({\rm rot})^{n}\mathbf{B}\rangle$ for a certain $n$. For $n=0$, this corresponds to an approximate equipartition of the energy density between the electric and magnetic components as found in the previous numerical studies Fujita:2015 ; Figueroa:2018 for small coupling constant. Going to higher order in $n$ seems to be a natural straightforward generalization.

The backreaction of the generated EM fields occurs when their energy density exceeds the value $\rho_{\rm EM}\sim(\sqrt{2\epsilon}/\beta)\rho_{\rm inf}$ slowing down the inflaton rolling along the slope of the potential. This makes the magnetogenesis less effective and the energy density of the EM field terminates the exponential growth of electromagnetic fields which remain almost constant until the preheating stage.

The generated strong electric fields lead to the charged particles production via the Schwinger process. Although there is also a particle production by the time-dependent metric of the expanding Universe, this contribution can be neglected in the strong field regime $\rho_{E}\gg H^{4}/\alpha_{\rm EM}$. Analytic expressions for the Schwinger current can be derived only in the simple case of a constant electric field in de Sitter background Kobayashi:2014 ; Froeb:2014 ; Bavarsad:2016 ; Stahl:2016a ; Stahl:2016b ; Hayashinaka:2016a ; Hayashinaka:2016b ; Sharma:2017 ; Bavarsad:2018 ; Hayashinaka:2018 ; Hayashinaka:thesis . However, as was shown in Ref. Geng:2018 , the form of the current for a time-dependent electric field in the strong field regime is the same as in the case of a constant field. We analyzed the importance of the Schwinger effect and found that it can affect magnetogenesis only if the inequalities $\rho_{E}\gtrsim(9\pi^{3}/8)H^{4}/\alpha_{\rm EM}^{3}\gg H^{4}/\alpha_{\rm EM}$ are satisfied, i.e., only in the strong field regime.

For numerical analysis, we used the $\alpha$-attractor model, which is a typical representative of the plateau models of inflation. The strength of the potential was chosen from the requirement of generation of the correct scalar power spectrum amplitude. We considered only the last 60 $e$-foldings of inflation although it can last much more. Taking this into account we imposed the initial conditions for the EM observables. Although the previous numerical studies Fujita:2015 ; Figueroa:2018 excluded the range of coupling constants $\beta>10$ in view of large inflaton perturbations and non-Gaussianities in the primordial spectra, we consider also large values of $\beta$ in order to test our method and to demonstrate more clearly the features of the backreaction regime.

We applied three approximation schemes with $n=0$, 1, 2, and compared their results with those found via the exact integration of the mode equation (34), which does not take into account the backreaction and the Schwinger effect. For early times when the EM field is weak and does not cause the backreaction, all approximation schemes give the results which satisfactorily agree with the exact result. Close to the end of inflation, the backreaction takes place for $\beta>10$ and we observe the changes in the behavior of the EM energy density as predicted analytically. The Schwinger effect comes into play also only near the end of inflation. It suppresses the electric energy density and produces charged ultrarelativistic particles which can make the Universe radiation dominated even before the reheating stage. Indeed, for large values of the coupling constant $\beta\gtrsim 20$, the energy density of created particles becomes comparable with that of the inflaton at the end of inflation. This means that the Schwinger preheating should be also considered as an important complementary scenario to the standard one with the inflaton perturbative decay and parametric resonance.

Finally, we took into account the postinflationary evolution of the magnetic field and estimated its present day value. Because of the fact that the axial coupling to the inflaton enhances only one circular polarization of the EM modes, the generated magnetic fields have nontrivial helicity. Moreover, if the enhancement is significant, then the generated magnetic field is close to maximally helical because the second polarization may be neglected at all. The evolution of such a magnetic field is not described by the adiabatic decay due to the Universe expansion $B\propto a^{-2}$. In the turbulent plasma the magnetic field undergoes the inverse cascade when the magnetic energy is transferred from the short to large scales due to the helicity conservation. This helps to increase the correlation scale of the generated magnetic field and, as a result, to avoid the cosmic diffusion after the recombination stage. Taking the most optimistic estimate, the present day value of the magnetic field cannot exceed $10^{-15}$ G with the correlation length of 1 pc. This means that the effective magnetic field strength, which should be compared with the distant blazar observations, equals $B_{\rm eff}=B_{0}\sqrt{\lambda_{B,0}/(1\,{\rm Mpc})}\sim 10^{-18}\,{\rm G}$ and can hardly explain the observational data.

In this work we did not consider the inflaton and curvature perturbations generated by the EM field. Since the EM field affects the inflaton evolution, it can spoil the spectrum of the scalar and tensor perturbations probed by the CMB observations. This can also be used to constrain the parameter space of the model. However, the self-consistent inclusion of these perturbations into the system of equations would significantly complicate calculations and, therefore, it should be addressed elsewhere.

Acknowledgements.

The work of O. O. S. was supported by the ERC-AdG-2015 Grant No. 694896 and by the Swiss National Science Foundation Grant No. 200020B_182864. The work of S. I. V. was supported by the Swiss National Science Foundation Grant No. IZSEZ0_186551. S. I. V. and O. O. S. are grateful to Professor Mikhail Shaposhnikov for his kind hospitality at the Institute of Physics, École Polytechnique Fédérale de Lausanne, Switzerland, where a part of this work was done. S. I. V. is grateful to Professor Marc Vanderhaeghen and Dr. Vladimir Pascalutsa for their support and kind hospitality at the Institut für Kernphysik, Johannes Gutenberg-Universität Mainz, Germany, where a part of this work was done. The work of E.V.G. was supported partially by the Ukrainian State Foundation for Fundamental Research.

Appendix A Derivation of the initial conditions

In this appendix we present the explicit expressions for the initial values of the EM quadratic functions (18)–(20). They can be calculated in the approximation of an adiabatically slow change of the parameter $\xi$ defined in Eq. (35). Substituting the approximate solution (61) into Eqs. (37)–(39), we obtained the initial conditions (62)–(64) in the integral form as functions of $\xi_{0}$ which is the value of $\xi$ at the initial moment of time when the simulation starts.

The corresponding integrals can be calculated in terms of Meijer $G$ functions Luke:book . We have

[TABLE]

[TABLE]

[TABLE]

For further convenience, we derive also the asymptotic expressions in the limiting cases of small and large $|\xi_{0}|$. For $|\xi_{0}|\ll 1$, we can expand the Macdonald functions to the leading order for $z\to 0$ and then the integrals can be easily calculated

[TABLE]

where $b=\exp(\gamma_{E})$ and $\gamma_{E}=0.577...$ is the Euler-Mascheroni constant. In the opposite case $|\xi_{0}|\gg 1$, the upper bound of integration in Eqs. (62)–(64) can be extended to infinity. Then by using Eqs.(109)–(111), we obtain

[TABLE]

We present here also the integrals used in the derivation of Eqs. (106)–(107)

[TABLE]