Anharmonic effects on the squeezing of axion perturbations

It is assumed in standard cosmology that the Universe underwent a period of inflation in its earliest phase, providing the seeds for structure formation through vacuum fluctuations of the inflaton scalar field. These fluctuations get stretched by the quasi-exponential expansion of the Universe and become squeezed. The aim of this paper is to deepen the understanding of the squeezing process, considering the effect of self-interactions. Axion-like particles can provide a useful setup to study this effect. Specifically we focus on the consequences that a non-trivial evolution of the background axion field has on the squeezing of the perturbations. We follow the evolution of the axion's fluctuation modes from the horizon exit during inflation to the radiation-dominated epoch. We compute Bogoliubov coefficients and squeezing parameters, which are linked to the axion particle number and isocurvature perturbation. We find that the quantum mechanical particle production and the squeezing of the perturbations are enhanced, if one accounts for anharmonic effects, i.e., the effect of higher order terms in the potential. This effect becomes particularly strong towards the hilltop of the potential.


Introduction
The current version of the standard cosmological model assumes that the Universe underwent a period of quasi-exponential expansion in its earliest stage, which solves several puzzles of the Big Bang cosmological scenario.This period of inflation turned out to be also an excellent mechanism to explain the origin of structure formation in the Universe.The exponential stretching of the fluctuations of the scalar field driving inflation, the so-called inflaton field, provides the seeds for the later gravitational collapse.In most inflationary models, these fluctuations have a quantum origin, however the CMB sky we observe today is classical.Therefore the question about the quantum or classical origin of the initial perturbations and how to probe it arise.
The classicalization of a quantum perturbation has been a subject of active research (see for example [1][2][3][4][5][6]).Actually, inflation itself provides an explanation for the "classicalization" of the originally quantum perturbations: they are squeezed due to the fast expansion of the Universe.A squeezed state is a special quantum state for which one variable has an arbitrarily small uncertainty, while its conjugate counterpart has a very big uncertainty, correspondingly.Since squeezed states are associated with a huge uncertainty in one variable, they are highly quantum mechanical states [7].However, within standard cosmological measurements, a squeezed state is indistinguishable from a classical phase-space distribution, if one considers a Gaussian state [8].Most of the previous works have focused on the study of nearly free fields.The aim of this paper is to deepen the understanding of the squeezing process considering the effect of self-interactions.
In general, to address the problem of the nature of the perturbations there is one further mechanism that takes place in the Early Universe that has to be taken into account: the decoherence during reheating.In fact, in order for a squeezed fluctuation to maintain its quantum nature into the later Universe, it is necessary that it does not undergo decoherence through the interactions with the cosmological plasma.The phenomenon of decoherence has been widely discussed in the literature; see for example [9][10][11][12][13].
From this perspective, axion-like particles [14,15] coupled to other matter fields only through gravitational interactions provide an ideal setup that can evade decoherence.In particular, axions can be produced non-thermally from an initial displacement of the axion field, which is called the misalignment mechanism.In this model axion-like particles are not required to directly interact with the thermal bath.Another interesting feature of axion-like particles is that their fluctuations around the classical axion background source axion isocurvature perturbations, that can be directly linked with observables.Their correlators constitute in principle a promising way to investigate the quantum origin of the cosmological perturbations.Indeed axions may evade decoherence and their bispectrum could provide a way to discern among the quantum or classical nature of the perturbations (see [16][17][18]).One further advantage in working with axions is that several constraints are already present in the literature, giving rise to strong bounds on the cosmological observables.Also for the purpose of studying interacting fields, axions are useful since they necessarily have self-interactions due to their periodic potential.The effects induced by higher-order terms in the potential is often referred to as "anharmonic effects," and axions provide a computable setup for incorporating the effect of self-interactions.Indeed in the late Universe the axion field settles down to the minimum of its potential and thus the self-interaction switches off.Axions thus reduce to effectively free fields in the asymptotic future, which allows one to apply Bogoliubov computations for studying the evolution of the fluctuations.
In this work we study the effect of self-interactions on the squeezing of axion fluc-tuations, focusing on how the background evolution modifies the evolution of the perturbations themselves.We focus on axion-like particles (instead of QCD axions [19][20][21]), and consider a scenario where the Peccei-Quinn symmetry is broken during inflation.In this way the axion field gets homogenized, giving rise to a background axion condensate.
We study the perturbations around this background.In previous works on axion squeezing [22,23], only the quadratic term has been considered.In this paper we go beyond these analyses, considering the effect of self-interactions.In particular we study the interaction between the perturbations and the zero-mode, i.e. the background axion field, neglecting the interactions among perturbations themselves.We find that the interaction with a nontrivial background gives already interesting effects.In particular, it strongly enhances the squeezing of the perturbations, which is directly linked with the number of particles created due to the expansion of the background.It has been known that anharmonic effects in the misalignment scenario give rise to the enhancement of the axion density [24,25], as well as the isocurvature perturbation [26][27][28].We revisit this effect using Bogoliubov calculations, and shed light on the anharmonic enhancement on the axion squeezing.The paper is organised as follows.In Section 2 we present our setup, and introduce the concepts of Bogoliubov coefficients and squeezing parameters.In Section 3 the evolution in time of the background field is analysed.In Section 4 the numerical results for the mode functions, the Bogoliubov coefficients and the squeezing parameters are reported.Their interpretation in terms of cosmological particle creation and squeezing are also presented.Moreover the dependence of the physical observables on the parameters of the system is studied.We conclude and discuss our results in Section 5.In Appendix A we give explanations on the squeezing parameters, while in Appendix B we provide analytic calculations for a quadratic axion potential.Finally, in Appendix C the cases of a noninstantaneous reheating and a quasi-de Sitter spacetime are considered.

Setup
We fix the metric to a flat FRW, ds 2 = a 2 (τ ) −dτ 2 + dx 2 .Here τ is the conformal time, linked with the cosmic time via dt = a(τ ) dτ .We consider the evolution of the Universe as simply given by a de Sitter (dS) phase of inflation, followed by a period of radiation domination (RD).We further assume, for simplicity, that the reheating phase is instantaneous1 .In the simple scenario here introduced, the Hubble rate is given by: Here H inf is the constant Hubble scale during the de Sitter phase and a reh is the value of the scale factor at the reheating.We consider the axion to have negligible effect on the cosmological evolution.It is assumed that the Peccei-Quinn symmetry is broken before inflation, and specifically that f > max(T reh , H inf ), where f is the symmetry breaking scale and T reh is the reheating temperature.This ensures that the symmetry is not restored during the reheating phase.For our assumption of instantaneous reheating, Let us now consider the axion potential, which can be written as where ϕ is the axion field.We assume that the axion mass is constant in time, and also that m < H inf ; this ensures the axion mass becomes comparable with the Hubble rate during radiation domination, after reheating has occurred.After the time when m ∼ H inf , the field starts oscillating around its minimum and the equations of motion becomes that of a damped oscillator.
We consider the axion field to be a homogeneous background plus vacuum fluctuations.The initial value of the axion field affects the subsequent evolution of the background.Specifically, if the starting point is near the hilltop of the potential then the field scans the full cosine potential as it rolls down.When deviations from the quadratic potential are considered, the perturbations are affected by the background evolution.These perturbation modes exit the horizon during inflation.After inflation has ended they reenter the horizon.Small-scale modes re-enter the horizon before the axion field begins to oscillate; hence they propagate unaffected by the oscillations.On the contrary, large-scale modes feel the oscillations and their evolution in time is consequently affected.Since we focus on the effect of the oscillations on the evolution of the perturbations, we will focus on large-scale modes that re-enter the horizon after the onset of the oscillations.

Canonical Quantization
In order to compute the squeezing parameters describing the system, we first need to calculate the Bogoliubov coefficients for the axion field.We work with the full potential for the background field evolution, while considering only the first order expansion of the potential for the perturbations.The equation of motion for the axion field is: Here the dot indicates a derivative with respect to the physical time t, and Latin letters are used to denote spatial indices.We split the axion field as where φ is the homogeneous classical background field, and δϕ is the small perturbation around it.We are going to expand the axion up to quadratic order in δϕ.The first derivative of the potential, up to linear order, can be written as The equation of motion for the homogeneous background field becomes while at first order in the perturbations one has Here we have moved to Fourier space, and we used k = |k|.The corresponding quadratic action for the perturbations is (2.9) Here the prime denotes a τ -derivative.One can now introduce a new variable χ = aδϕ.
Then the action can be rewritten as (2.10) The mixed term χχ ′ can be removed by integrating by parts and dropping surface terms as From this action we can construct the Hamiltonian as H = d 3 x (pχ ′ − L), where p = χ ′ and L is the Lagrangian.Writing in terms of the Fourier components, the Hamiltonian acquires the following form: where the effective mass m2 eff is given by Notice that the effective mass depends on time through the evolution of the background φ.Moreover it can become tachyonic, i.e. m 2 eff < 0, due to the cosine factor.We also introduce the time dependent frequency ω k , defined as In order to quantize the fields, we introduce time-dependent ladder operators as follows: respecting the canonical commutation relations 2 (2.17) All the other commutators vanish.One should note that the frequency ω k can become imaginary, hence we chose to normalize the fields with |ω k | in (2.15) and (2.16).

Bogoliubov Transformation
Let us introduce an 'initial' time τ 0 during inflation when the wave modes of interest are still deep inside the Hubble horizon.Then the time-dependent ladder operators at later times are related to those at τ 0 through a Bogoliubov transformation, where a 0 k = a k (τ 0 ).From the commutation relations for the ladder operators, one finds the following relation among the Bogoliubov coefficients: (2.19) The variables χ k and p k can be written in terms of the time-independent ladder operators directly, introducing the mode functions u k (τ ): ) The normalization condition for the mode functions is given by It arises from the commutations of the real space variables and the ladder operators.The mode functions can then be used to find an expression for the Bogoliubov coefficients α k and β k .Indeed, rewriting (2.15) and (2.16) one gets: With these definitions for the Bogoliubov coefficients, and using (2.22), we can compute a quantity which will be useful later: Having introduced the Bogoliubov coefficients, we need to interpret these quantities in terms of physical observables.Their evolution can be understood in terms of particle creation by a spacetime-dependent background [29][30][31][32].Consider two sets of ladder operators (a,a † ) and (b,b † ) and their corresponding vacua defined through a|0⟩ a = 0 and b|0⟩ b = 0.These sets can be linked via a constant Bogoliubov transformation: (2.27) The two vacua |0⟩ a and |0⟩ b do not necessarily coincide.The expectation value of the number operator of one basis with respect to the vacuum of the other basis is given by |B| 2 .In other words, an observer in the a-basis feels the vacuum as filled with b-particles.
In the cosmological case, where time-dependent Bogoliubov transformations has been introduced, the instantaneous vacuum defined by the time-dependent ladder operators (a k (τ ), a † k (τ )) is filled with particles associated with the initial-time operators (a 0 k , a † 0 k ).However, unlike in Minkowski space, in a general curved spacetime there is no way to uniquely define a vacuum state.This ambiguity can be overcome if the effective spacetime seen by the fluctuations approaches Minkowski both in the far past and in the far future, namely, if adiabatic vacua exist asymptotically such that the mode functions are given by plane waves.There the fluctuations can be interpreted as particles, and the Bogoliubov coefficients relating the ladder operators in the asymptotic past and future give the final particle number as To verify the above discussion explicitly, let us compute the expectation value of the Hamiltonian (2.12).Using (2.15) and (2.16), we can rewrite the Hamiltonian in terms of the ladder operators as (2.28) Since we are interested in the late time limit of (2.28), we can simply drop the absolute value for ω k and simplify the previous expression.In fact during radiation domination the scale factor follows a ′′ /a = 0, and after the onset of the oscillations the effective mass squared is positive, hence ω 2 k becomes positive.Taking the expectation value over the initial vacuum, i.e. the vacuum annihilated by the initial time-independent annihilation operator, a 0 k |0⟩ = 0 for all k, we get: The divergent vacuum energy can be removed by a normal ordering.The expectation value for the ladder operators can be computed from (2.18) as follows: where given in terms of the mode functions, as in (2.26).Notice that the particle interpretation is valid only when the perturbation mode does not feel the expansion, i.e. in the subhorizon limit.In this limit, ω k is positive and can be interpreted as the frequency of the perturbation mode.|β k | 2 can then be interpreted as the expectation value of the number of particles with comoving momentum k, created per comoving volume due to the timedependent background.The divergent factor δ 3 (0) arises because we are considering an infinite spatial volume, but the mean density of particles per comoving volume is finite.We therefore expect the norm of β k to approach a constant in the asymptotic future.

Squeezing Parameters
Thanks to the constraint (2.19), the two complex Bogoliubov coefficients can be parameterised by three real variables, the so-called squeezing parameters, as: where r k is restricted to be non-negative.Inverting these relations yields the expressions for the squeezing parameters in terms of the Bogoliubov coefficients: (2.32) The process of particle creation explained in Section 2.3 can be equivalently described by means of the squeezing formalism.In particular, the cosmological particle creation amounts to squeezing the vacuum.The physical meaning of the squeezing parameters can be understood studying a quadrature pair in phase space of the physical system of interest.A squeezed state is given by a quadrature pair for which the uncertainty in one variable is small, while the uncertainty in the other variable is correspondingly large, in order to respect the Heisenberg uncertainty principle.This state is described by a rotated ellipse in phase space.The r k parameter is directly linked with the semi-axes of this ellipse, while the φ k parameter is related to its orientation.In particular, the r k parameter tells us how much the ellipse is squeezed along one direction while elongated along the other.The ϑ k angle is actually unphysical.A detailed description of the physical meaning of the squeezing parameters is given in Appendix A.
The squeezing parameters are also linked to cosmological observables.We have seen in Section 2.3 that |β k | 2 can be interpreted as the number of particles created by the expansion of spacetime.The parameter r k is given in terms of this net number of particles, as shown in (2.32).The parameter φ k is instead linked to the power spectrum of the field fluctuations.As explained in Appendix A, the initial state we are working with, i.e. the vacuum state, is rotationally invariant; hence only one phase is physically meaningful.From the expectation value of the field fluctuation: the power spectrum is written as, (2.35) In terms of the squeezing parameters, arg (α k β k ) corresponds to φ k , which is therefore a physical quantity.

Background Evolution
To study the evolution in time of the Bogoliubov coefficients and then in turn of the squeezing parameters, we first analyse the background field dynamics by numerically solving the homogeneous equation of motion (2.6).It is convenient to introduce the background displacement angle θ = φ/f, ( and solve the equation in terms of this variable.In this way we do not need to specify the value of the symmetry breaking scale f .The equation of motion for the background angle in conformal time is given by: We expect the axion background to start oscillating after the Hubble rate becomes comparable with the axion mass.Solving the equations of motion numerically for different initial field values during inflation, we obtain the results reported in Figure 1 (left).Here we considered the Hubble rate in equation (2.1) with H inf = 10 8 GeV, and took the mass as m = 10 2 GeV.Then the third term in (3.2) becomes relevant at around ln(a/a reh ) = ln(H inf /m) 1/2 ≈ 7, as shown in Figure 1 (left).The different colours represent different initial angles as shown in the legend, while the vertical dashed line identifies when the mass of the axion field becomes comparable with the Hubble rate.If we choose as the initial field value θ 0 = 0.1 π, it is close enough to the minimum such that the axion potential is well approximated by a quadratic.Indeed, if we compute the background trajectory with a quadratic potential we find the same result.Looking at Figure 1 (left) we see that starting close to the minimum, the amplitude of the oscillations is reduced.A huge difference in both the amplitude of the oscillations and the time of the onset of the oscillations can be seen if we set the initial field value to be very close to the hilltop.This behaviour can be further understood by looking at the energy density of the background field.The energy density can be obtained as Given the solution for the background field evolution, the energy density can be directly computed.The time evolution of the energy density normalized by f 2 is reported in Figure 1 (right), which zooms into the RD epoch around the time when the oscillations begin.The energy density is constant until the background field starts oscillating; thereafter it decays as a −3 .The vertical dashed lines in Figure 1 (right) show when the oscillations begin; this is obtained as the time when the energy densities extrapolated from the asymptotic behaviours in the past (ρ = const.)and in the future (ρ ∝ a −3 ) cross each other.The value of the energy density in the asymptotic future in terms of the initial displacement angle is reported in Figure 2. It has be computed at the final time ln(a/a reh ) = 10.The final numerical value of ρ/f 2 itself is not relevant for our discussion.What is relevant is that the energy density becomes strongly enhanced towards the hilltop.This behavior can be understood by noticing that the field starts oscillating later as one approaches the hilltop.

Evolution of Axion Fluctuations 4.1 Mode Functions
In order to evaluate the evolution in time for the Bogoliubov coefficients and the squeezing parameters, we need firstly to solve the equation of motion for the mode functions, given by: If the potential has only the mass term, an analytical solution can be obtained, while if we work with the full cosine potential we have to rely on numerical computations.The analytical derivation for a mass term potential is reported in Appendix B.1.
The time evolution of the amplitude of the mode functions is reported in Figure 3.It has been numerically computed for two different values of the initial background field during inflation: θ 0 = 0.1 π and θ 0 = 0.999999 π.The initial condition for the mode function is set such that it starts as a positive-frequency WKB solution obeying the normalization condition (2.22), when the wave mode is deep inside the horizon during inflation.The parameters chosen for the evaluation are H inf = 10 8 GeV, m = 10 2 GeV and k/a reh = 10 2 GeV.These values have been chosen in order to consider a mode that leaves the horizon during inflation, and stays outside the horizon when the background field starts oscillating around the minimum, during radiation domination.The dashed vertical lines in Figure 3 indicate the time of the horizon exit, reheating and the onset of the oscillations of the background field.
Looking at Figure 3, we see that the mode function's amplitude stays constant until horizon exit, then grows outside the horizon until the background field starts to oscillate.This behaviour is more or less independent of the initial misalignment angle (and thus the lines overlap).Instead, at the onset of the oscillations the amplitude of the mode functions differ based on the value of the initial displacement angle.If the initial field is close to the minimum, then the approximation to the quadratic potential is enough to describe the system.The time evolution obtained from the analytical solution for a quadratic potential matches the time evolution for an initial field value θ 0 = 0.1π.However, if the initial field is close to the hilltop an extra enhancement in the amplitude of the mode functions is observed due to the anharmonic effects.

Bogoliubov Coefficients
Exploiting relations (2.24), (2.25) and (2.32), it is possible to evaluate the evolution in time of the Bogoliubov coefficients and in turn of the squeezing parameters.We start our analysis from the Bogoliubov coefficients, and specifically from |β k | 2 , which is directly linked with physical quantities, as explained in Section 2.3.The numerical results are reported in Figure 4, for the same values of parameters used in Figure 3 for the amplitude of the mode functions.The modulus square of the beta coefficient approaches zero in the asymptotic past, when the mode is far inside the horizon.Moreover it approaches a constant value in the asymptotic future, and specifically it becomes constant when the axion field starts oscillating.We further notice that the behaviour of |β k | 2 changes depending on the initial field value.The evolution in time is more or less independent of the initial condition while the axion field is frozen by the Hubble friction; however, when the potential term becomes comparable with the Hubble rate, the field dynamics right before the onset of the oscillations affects the asymptotic value of the Bogoliubov coefficients.We see in Figure 4 that, approaching the hilltop of the potential, the asymptotic value of |β k | 2 increases.On the other hand with θ 0 = 0.1 π, the evolution of |β k | 2 follows that for a quadratic potential, which is also analytically studied in Appendix B.1.
The spikes in the evolution of |β k | 2 seen in the plots, as well as the discontinuous jump at reheating, can be understood in terms of the frequency (2.14).In the asymptotic past, ω 2 k is dominated by k 2 .When the perturbation leaves the horizon during inflation, at k = aH inf , the frequency squared reduces to ω 2 k ≃ −k 2 .This implies that right before the horizon exit the frequency squared has turned from positive to negative.This gives a divergence in |β k | 2 , since it is defined in (2.26) with ω k in the denominator.The spike after reheating has a similar explanation.During RD, ω 2 k = k 2 + m 2 eff a 2 ; the k 2 term is always positive, while m 2 eff could be either positive or negative, depending on the value of the cosine.If the initial field is close to the hilltop, then m 2 eff is originally negative.With the parameter choice in the plot (m = k/a reh = 10 2 GeV), the tachyonic mass term  dominates the frequency right after reheating.This implies that the frequency squared flips its sign immediately after reheating, causing the spike in |β k | 2 .There is another spike at the onset of the oscillations, depending on the initial displacement angle.In fact, when the background field starts oscillating, m 2 eff goes from negative to positive.Since ω 2 k ≃ m 2 eff a 2 at that time, it crosses zero, making |β k | 2 diverge.The jump at reheating is instead caused by the discontinuity of the frequency ω k , which arises due to a ′′ /a suddenly vanishing upon instantaneous reheating.
However we remark that the spikes and jumps of |β k | 2 appear only while adiabaticity is violated and hence |β k | 2 cannot be interpreted as the particle number, as we will explain in the next subsection.The irregular behaviours are artifacts due to our definition of |β k | 2 , as well as the choice of the cosmological background, and hence should not have physical consequences.We actually show in Appendix C that a non-instantaneous reheating removes the jump of |β k | 2 , but does not alter the final value of |β k | 2 .

Adiabaticity Condition
We have seen that |β k | 2 approaches a constant when the background field starts oscillating.In order to understand why this is the case, let us analyse the so-called adiabaticity condition, which is given by: Here ω k is the effective frequency associated with the mode, defined in (2.14).If ω k is real and the adiabaticity condition holds, then the solution to the equation of motion for the perturbations should approach a plane wave solution.In particular, it implies that the WKB approximation holds so that the mode function can be written as When this is the case, the time-dependent Bogoliubov coefficients α k and β k reduce to the C-numbers A k and B k , up to a time-dependent phase.
Let us analytically assess the adiabaticity condition for a quadratic potential, i.e. m 2 eff = m 2 .In the dS epoch, the first adiabaticity parameter |ω ′ k /ω 2 k | is tiny before the mode exits the horizon, however when the mode is far outside the horizon then it settles to a constant value of 1/ √ 2. This can be shown by ignoring the axion mass which is much smaller than the inflationary Hubble scale, giving The second adiabaticity parameter ω ′′ k /ω 3 k exhibits a similar behaviour.The adiabaticity condition becomes satisfied again after inflation when the axion field starts oscillating.During radiation domination, the first adiabaticity parameter reduces to Given this result, there are two possible regimes: We restrict ourselves to wave modes that are outside the horizon at the onset of the oscillations, i.e. k < a osc m.Hence after the oscillation begins, namely when a > a osc and H < m, the adiabaticity parameter is given by the second line, which is smaller than unity.Therefore starting from that point the adiabaticity condition is always recovered.The behaviour long before a osc depends on the actual value of k.The second adiabaticity parameter is Considering again the two possible regimes, we have: Also in this case the condition is surely satisfied when a > a osc .
The above discussions apply similarly to the full cosine potential, for which we numerically checked with initial misalignment angles close to the hilltop.The behaviours of the adiabaticity parameters are the same while the axion is frozen by the Hubble friction.Soon after the axion starts to oscillate the adiabaticity parameters also oscillate, however, within a few e-folds the adiabaticity condition becomes satisfied.

Enhancement in the Number of Particles
In order to understand the enhancement in the beta coefficient due to the anharmonic effects, we show in Figure 5 (left, purple line) the value of |β| 2 in the asymptotic future as a function of the initial field value θ 0 = φ0 /f .(The pink line will be explained in Section 4.2.4.)The asymptotic value of the Bogoliubov coefficient significantly increases as one approaches the hilltop at θ 0 = π.This behaviour can be interpreted analysing the time at which the field starts oscillating.As we have seen in Section 3, increasing the initial field value delays the onset of the oscillations.Approaching the hilltop this delay becomes more and more evident.If we consider the field perturbation around a background value close to the minimum, the delay in the onset of the oscillation is too tiny to affect the evolution in time of δϕ k .However, if the field perturbation is considered around a background field value close to the maximum of the potential, then the evolution in time of the perturbation changes dramatically.In fact, close to the hilltop even a small difference in the initial position leads to a huge difference in the time of the onset of the oscillations.This implies that patches of the Universe that differ by a tiny variation of the initial misalignment angle will start to oscillate at very different times, sourcing huge fluctuations.This in turn reflects into the modulus square of the mode functions, that enters directly into |β k | 2 .In the hypothetical situation where the angle is set initially to π, the field does not start oscillating at all.Hence δϕ k is going to become infinite, thus making |β k | 2 divergent.
To better analyse the behaviour of |β k | 2 near the hilltop it is worth working in terms of π/(π − θ 0 ).In terms of this variable, the hilltop at θ 0 = π is shifted to ∞.In Figure 5 (right, purple line), |β k | 2 is reported in terms of this new variable.We find that the asymptotic value of |β k | 2 grows towards the hilltop3 as (π − θ 0 ) −2 .
We should remark that the above discussions do not apply arbitrarily close to the hilltop.In particular if the initial displacement from the hilltop, π − θ 0 , is comparable to or smaller than the angular quantum fluctuation δθ (whose amplitude is of H inf /f ), then the perturbative expansion around the classical background breaks down.In this regime the field can roll down to the other side of the hill in some patches of the universe, and give rise to axionic domain walls.In this work we do not consider such cases.

Parameter Dependence
Let us now study the dependence of |β k | 2 on the parameters H inf , m, k.The dependence on the parameters can be derived analytically for a quadratic potential.The complete derivation is reported in Appendix B.2, where we show that the asymptotic value of |β k | 2 for a quadratic potential is given by A physical interpretation of this result in terms of the energy density will be presented in Section 4.2.4.The dependence of |β k | 2 ∝ k −3 shows the scale invariance of the perturbation.
The expression (4.9) applies for initial field values close to the minimum, such that the potential is well approximated by a quadratic.To understand the dependence on the parameters for field values where anharmonic effects are relevant, we have to rely on numerical simulations.The results are shown in Figure 6, where the asymptotic value of |β k | 2 as a function of the initial misalignment angle is plotted.The different lines correspond to different choices for the Hubble rate, the axion mass, or the wave number; the horizontal dashed lines are the analytical value computed using (4.9).What we find is that |β k | 2 always depends on the same combination of parameters as given in (4.9), which sets the overall normalization for |β k | 2 .Since the dependence of |β k | 2 on k is the same, we can also see that the energy density remains scale independent when the anharmonic effects are taken into account.Moreover, we see that the shape of |β k | 2 as a function of θ 0 is identical for each choice of the parameters.This indicates that the anharmonic enhancement of |β k | 2 is determined only by the initial angle.
The final value for |β k | 2 also allows us to compute the number of axions today.From (2.29) (or (4.17)), the physical number density of axion particles that were quantum me- chanically produced in the early universe is where n k denotes the number density of particles with comoving momentum ∼ k.(It should be noted that the total number of axions is dominated by the zero mode axions classically produced from the vacuum misalignment.)Using the expression (4.9) for axions produced close to the potential minimum, and rewriting the redshift at instantaneous reheating in terms of the inflation scale as the number density today is obtained as (4.12)

Alternative Calculation for |β k | 2
In order to understand the physics underlying the expression (4.9), let us present an alternative derivation of the asymptotic value of |β k | 2 using the energy density of the axion field: Expanding this in terms of the field fluctuation δϕ, then the expectation value of the density operator in the initial vacuum defined above (2.29) is written as Here ρ is the energy density of the background field as given in (3.3), and we considered terms linear in δϕ to have vanishing expectation values.The quadratic-order density operator ρ 2 is given by whose vacuum expectation value is written in terms of the Bogoliubov coefficients as This is equivalent to the asymptotic value of the Hamiltonian density ⟨H⟩/a 4 V , cf. (2.29).
We can also compute the axion density by treating it as a purely classical quantity.For this, notice that the axion field during inflation follows a slow-roll attractor where the field velocity is set by the field value as ϕ ′ ∝ −dV /dϕ.Hence the asymptotic energy density can be considered to be uniquely determined by the axion field value at some initial time τ i during inflation, i.e., ρ = ρ(ϕ i ).Then the spatial inhomogeneity of the asymptotic density can be understood to originate from that of the initial field, ϕ i (x) = φi + δϕ i (x), where φi is the average field value over the entire universe at τ i , and δϕ i is the classical fluctuation around the average.(We are being somewhat sloppy and using the same notation as we used in (2.4) for splitting between the classical background and quantum fluctuations, but we stress that here everything is classical.)Expanding the asymptotic density in terms of δϕ i (x), and also taking a spatial average yields By identifying the spatially averaged quantities with the vacuum expectation values of the corresponding operators, we can equate (4.18) with (4.14).Moreover, the variance ⟨δϕ 2 ⟩ avg can be evaluated as the vacuum expectation value of the quadratic operator, as shown in (2.33).Hence the second term of (4.18) is rewritten as Here we ignored modes that are sub-horizon at the time τ i , since the contributions from such modes should be renormalized and absorbed into the cosmological constant.Additionally, for the mode function during inflation, we substituted its value in the superhorizon limit for an effectively massless field, Upon comparing (4.14) and (4.18), let us identify ρ with ρ( φi ), and further ignore the higher-order correlators.Then we can equate the second terms of each expression, namely, (4.17) and (4.19).Let us further suppose that this equality holds mode by mode, and also drop the 1/2 in (4.17).Then we find a relation between the asymptotic value of the Bogoliubov coefficient, and the second-order derivative of the background energy density in terms of the initial field value, Since the energy density after the axion started oscillating redshifts as ρ ∝ a −3 , the righthand side is actually time-independent.
If the initial field value is close to the potential minimum, i.e. | φi | ≪ f , then the axion potential is approximately quadratic and the field begins a harmonic oscillation when H ∼ m.Hence the asymptotic axion density can approximately be written as where upon moving to the far right-hand side we used a osc ∼ a reh (H inf /m) 1/2 .Plugging this into (4.20),we can estimate the Bogoliubov coefficient around the minimum as This matches with the exact result (4.9) up to a numerical factor of order unity. 4 We have numerically computed both sides of (4.20); these are displayed as functions of the initial angle5 in Figure 5, where the purple line shows the left-hand side, while the pink line shows the right-hand side.At θ 0 ≲ 1, where the potential can be approximated by a quadratic, the two results are in good agreement.In particular, both asymptote to the value of (4.9) as θ 0 → 0. Close to the hilltop, on the other hand, |β k | 2 is consistently larger than the right-hand side of (4.20) by an order of magnitude.We have also checked that the θ 0 -dependences of either sides of (4.20) are unaltered by changing the values of m, H inf , and k; hence the behaviours discussed above are also parameter-independent.
We should remark that our calculation incorporates axion self-interactions through the interaction between the homogeneous background field (zero mode) and the fluctuation mode, however we have ignored mode-mode interactions by expanding the action only up to quadratic order in the fluctuation.The relation (4.20) has also been derived by neglecting higher-order correlators.One would naively expect the mode-mode interactions to be suppressed by a hierarchical parameter choice such that m ≪ H inf ≪ f .However the uniform discrepancy between both sides of (4.20) at θ 0 ≳ 1 may be implying the importance of the interactions between the fluctuation modes.We leave a detailed investigation of this for the future.

Squeezing Parameters
We can now proceed to the analysis of the squeezing parameters.Through relations (2.32), one gets the evolution in time of the squeezing parameters.The result is shown in Figure 7 4 The numerical factor of (4.9) can also be obtained from (4.20) by using the exact result for the asymptotic axion density [34] (see their Eq.(A3), and combine with H = H inf (a reh /a) 2 ): for r k and φ k , which are the two parameters linked with physical quantities, as explained in Section 2.4 and in Appendix A. Their evolution is reported for four different values of the initial misalignment angle.The r k parameter is directly linked with the modulus square of the beta coefficient; we can therefore see an enhancement in its value towards the hilltop.r k gives a measure of the stretching of the ellipse describing the phase space distribution of two conjugate variables, as explained in Appendix A. If r k increases, as observed taking into account the anharmonic effects, then the ellipse gets stretched along one of its axis, and the state described by the two conjugate variables becomes more and more squeezed.We therefore conclude that anharmonic effects give extra squeezing at the onset of the oscillations.
The angle φ k is instead linked with the rotation of the ellipse in phase space.φ k is constant until the mode exits the horizon, then shifts to another constant value after the horizon exit.When the background field starts oscillating, and the adiabaticity condition starts to hold, the phase starts to rotate.This behaviour is visible in Figure 7.The sharp lines are present since the phase is restricted to the range −π ≤ φ k ≤ π.In phase space, the ellipse described by this phase is constantly rotating by this time dependent angle.

Conclusions
We have computed the time evolution of the Bogoliubov coefficients and the squeezing parameters for an axion-like field during inflation and radiation domination, keeping into account self-interactions.In particular we considered the non-trivial evolution of the background field.We have shown that the net number of particles created due to the expansion of the Universe, given by |β k | 2 , increases as the initial background field approaches the hilltop of the cosine potential.The r k parameter determining the amount of squeezing increases accordingly.Hence we conclude that the self-interaction due to the cosine potential, and specifically the interaction between the perturbation and the non-trivial background field dynamics, gives rise to extra post-inflation squeezing.
An analytical expression for |β k | 2 in terms of the relevant parameters of the system has also been found in the quadratic potential limit in (4.9).A full derivation has been carried out in terms of the analytical mode functions in Appendix B.1.The expression (4.9) has also been confirmed via numerical computations to set the overall normalization of |β k | 2 .The behaviour of the final value of |β k | 2 as a function of the initial displacement angle is not affected by the physical parameters H inf , m and k, provided the mode considered reenters the horizon after the background has started oscillating.Moreover the anharmonic enhancement of the squeezing is solely determined by the initial misalignment angle.
We have also derived a relation between the asymptotic values of the Bogoliubov coefficients, and the field derivatives of the background energy density, as shown in (4.20).This relation exactly holds for small misalignment angles where the axion potential is well approximated by a quadratic, however we also found that it breaks down towards the hilltop.Here we note that in this paper we evaluated the effect of axion self-interactions through the interaction between the homogeneous background field and the fluctuation, while we have ignored interactions among the fluctuation modes.The breakdown of (4.20) towards the hilltop may be indicating that mode-mode couplings can further enhance axion squeezing.We leave an investigation of this effect for future work.
A connection with observables has also been pointed out.The net number of particles created by the expansion of the Universe is given by |β k | 2 , which is directly linked also to the parameter r k and consequently to the amount of squeezing of the perturbations.The more particles are created, the more squeezed the state is.The other physically relevant squeezing parameter, i.e. the phase φ k , is related with the power spectrum of the field fluctuations, as shown in (2.35).
One interesting question left for future work is how we could observationally discern among the quantum or classical origin of the axion perturbations.The main problem in addressing this question is the fact that the two-point correlation functions of a squeezed state can be perfectly mimicked by a stochastic classical distribution.Despite this, a squeezed state is an intrinsically quantum state and in particular it is an entangled state.Hence one possible test of the quantum nature of the perturbations is a Bell type experiment on the Cosmic Microwave Background (CMB).In the case of the curvature perturbations this kind of experiment is however inconclusive, see [35] and references therein.Even if considering Leggett-Garg inequalities [36] instead of Bell inequalities, the curvature perturbations could be useless, since the conjugate momentum of the curvature perturbation ζ, i.e. its time derivative, is vanishing on super-horizon scales.Working with axion fields may avoid this latter problem.We have indeed shown that the phase φ k oscillates in the late time limit, making the ellipse in phase space describing the system to rotate.This suggests that the time derivative of the perturbation is non-negligible, and hence axions may provide an ideal target for a cosmological Leggett-Garg experiment.
Another way to probe the nature of the primordial fluctuations is to consider interacting theories and compute higher order correlators.One attempt in this regard have been performed in [16] (see also in a different context, e.g., [17,18]).In the recent work by Green and Porto [16] it has been pointed out how the quantum origin of the perturbations could be singled out by looking at the poles of the bispectrum, therefore a careful analysis of the axion bispectrum might be enlightening.In order to compute the bispectrum of an axion field, the effect of the non-trivial background computed in this work will have to be taken into account.
Although we have focused on axion fields throughout this paper, the study of the squeezing of the perturbations can also be important for addressing the evolution of gauge fields that could have been excited in the early Universe.Recently a paper has been published on this topic [37] and further analysis along this direction is an interesting expansion of our work.
Figure 8: Phase space distribution of a conjugate pair (x k , y k ).The green circle corresponds to a vacuum state, while the green ellipse is a squeezed state described by the parameter r k and rotation angle φ k .This picture has been adapted from [45]. of these two operators can be computed using the expansion e Rewriting the squeezing parameters in terms of the Bogoliubov coefficients using (2.32), we see that (A.4) matches with (2.18).These expressions also show that the modes +k and −k are mixed, hence a two-mode squeeze operator is needed to describe particle pairs in states ±k.
The rotation operator shifts the phase of the ladder operators as Ra 0 The vacuum state is invariant under this rotation [39,40], and hence the parameter ϑ k is unphysical.Physical observables thus depend only on the parameters r k and φ k .Using (2.32), this in turn implies that the relevant combinations of the Bogoliubov coefficients are |β k | 2 and arg (α k β k ).In particular, the squeezing parameter r k is related to the number of particles created.The expectation value of the number operator in the vacuum annihilated by a 0 k follows from (A.4) as, as we have also derived in (2.30), and V is the spatial volume.Thus cosmological particle creation amounts to squeezing the vacuum, as observed by Grishchuk and Sidorov [2].
In order to understand better the physical meaning of the squeezing parameters we need to move to the phase space representation for a quadrature pair.A squeezed state is described by a quadrature pair for which fluctuations in one variable are reduced below the symmetric quantum limit, while fluctuations in the other canonically conjugate quadrature analytical results presented coincide with the numerical results obtained for the full cosine potential in the limit of a small initial displacement angle, close to the minimum of the potential.

B.1 Mode Functions
The equation of motion for the mode function is given by (4.1), with m eff replaced by a constant and positive mass m, as here we focus on a quadratic potential.Let us solve this equation in the dS and RD epochs, assuming the two epochs to be connected by an instantaneous reheating.Here we set the time at reheating as τ reh = −1/a reh H inf .The equation of motion then reduces to the form, Supposing m 2 /H 2 inf < 9/4, the solution is given in terms of Hankel functions as, 3) The two constants c 1 and c 2 can be fixed by requiring that the mode function approaches a positive frequency solution in the asymptotic past, and also that it satisfies the normalization condition (2.22).This yields, up to an unphysical phase, where we have rewritten the solution in terms of the scale factor.
RD epoch In a radiation-dominated universe where H ∝ a −2 , the equation of motion can be rewritten as Solutions to this equation are the parabolic cylinder functions, The two constants c 3 and c 4 are fixed by imposing continuity of the solutions (B.4) and (B.7), as well as their time derivatives, at the time of reheating.

B.2 Asymptotic |β k | 2
In order to evaluate the Bogoliubov coefficient |β k | 2 in the asymptotic future, we introduce linear combinations of parabolic cylinder functions, and rewrite the mode function in the RD epoch (B.7) as Here A k and B k are time-independent coefficients.The functions F + and F − are chosen such that in the asymptotic future, they respectively approach the WKB solutions with positive and negative frequencies: (a reh H inf z) This gives a good approximation for (B.13) at κ 2 ≲ µ ≪ 1; here note that κ 2 ≲ µ is equivalent to k ≲ (aH) osc .Hence for wave modes that are outside the Hubble horizon when the axion begins to oscillate during the RD epoch, the asymptotic value of the Bogoliubov coefficient is obtained as . (B.17) C More Realistic Models for the Evolution of the Universe In the results derived in the main text, we assumed a pure dS spacetime followed by an RD spacetime, further assuming the reheating phase to be instantaneous.In this Appendix, we check that changing the evolution of the Hubble rate does not affect our main conclusions on the asymptotic value for |β k | 2 and the squeezing parameters.

C.1 Smooth Reheating
We first smooth the Hubble rate evolution between the dS and RD epochs to see if our results are affected by a non-instantaneous reheating phase.The expression we use is the following: We numerically solved the axion's equations of motion with the Hubble rate (C.1), and we found that it does not change the evolution of the background field nor the fluctuation's mode function.However it affects the evolution of β k as shown in Figure 9.The behaviour of |β k | 2 is modified near the reheating phase, but it coincides with the result from instantaneous reheating after a few e-folds.Hence the final value of |β k | 2 is not affected by this modification to the Hubble evolution.This is because the post-inflation squeezing is determined by the field dynamics near the onset of the oscillations; the Hubble rates in (C.1) and (2.1) differ for a few e-folds around reheating, however they coincide by the time the oscillations begin.

C.2 Single Field Inflation
We now study a quasi de Sitter inflationary background, by considering a single-field inflation model with a quadratic potential    We assume that the background dynamics is completely determined by the inflaton, with the axion field being only a spectator field.We also suppose instantaneous reheating, followed by a period of radiation domination.Therefore the Hubble rate is now given by: In contrast to the exact dS background invoked in the main text, now the Hubble rate varies in time during inflation, with the inflation scale set by the inflaton mass.Upon computing the evolution of the axion fluctuations in the above cosmological background, we have chosen m φ = 10 8 GeV (hence H inf ∼ 10 8 GeV), combined with m = 10 2 GeV and k/a reh = 10 2 GeV; these parameters have the same order of magnitude as those in Figure 9.We checked that the mode function present similar behaviours as in the case of an exact dS spacetime.The results for the time evolution of |β k | 2 are reported in Figure 10 for different initial misalignment angles.The general features of the evolution of |β k | 2 are similar to those in Figure 9 for an exact dS spacetime followed by an instantaneous reheating; |β k | 2 approaches a constant value after the onset of the oscillations, and its numerical value depends on the initial field.
From these analyses we conclude that our main conclusions for the axion squeezing are independent of the details of the inflation and reheating models.

Figure 1 :Figure 2 :
Figure 1: Time evolution of the background field (left) and of the background energy density (right) for different initial field values.The values of the initial angle are specified in the legend in the left plot.The parameters chosen are H inf = 10 8 GeV and m = 10 2 GeV.The grey dashed vertical line shows the time when the axion mass becomes equal to the Hubble rate.Dashed vertical lines in other colours denote the onset of the axion oscillation for the different initial displacement angles.

Figure 3 :
Figure 3: Evolution in time of the amplitude of the mode functions for a wave mode k/a reh = 10 2 GeV, and for two different initial values of the background field.An increase in |u k | 2 is observed approaching the hilltop of the potential.Dashed vertical lines denote when the mode exits the horizon, reheating and the onset of the oscillations, respectively.The values chosen for the Hubble rate and the axion mass are H inf = 10 8 GeV and m = 10 2 GeV.

Figure 4 :
Figure 4: Evolution in time of the modulus square of the β k -coefficient for few different initial values of the background field.The parameters chosen are the same as Figure 3.When θ 0 approaches the hilltop of the potential the final value of |β k | 2 significantly increases.Dashed vertical lines denote when the mode exits the horizon, reheating and the onset of the oscillations, respectively.

Figure 5 :
Figure 5: Asymptotic value of |β k | 2 as functions of the initial displacement angle (left), and of π/(π − θ) (right).The parameters are chosen as H inf = 10 8 GeV, m = 10 2 GeV, and k/a reh = 10 2 GeV.The purple line refers to |β k | 2 computed through the mode functions.The pink line shows the second derivative of the axion density in terms of the initial field value, as given in the right-hand side of (4.20).

Figure 6 :
Figure 6: Asymptotic value of the Bogoliubov coefficient |β k | 2 as a function of the axion misalignment angle θ 0 .Colours denote different choices of the inflationary Hubble scale, axion mass, and the wave number upon reheating, as listed in the legend in units of GeV.Dashed horizontal lines denote the offset at θ 0 → 0, as given in (4.9).

(4. 16 )
After the onset of the axion oscillations the effective mass asymptotes to the mass around the minimum, m eff ≃ m (≫ H, k/a), and the frequency to |ω k | ≃ am.In this epoch(4.16)takes the form of

Figure 7 :
Figure 7: Time evolution of the squeezing parameters r k and φ k for different values of the initial misalignment angle.The parameters chosen for the evaluation are H inf = 10 8 GeV, m = 10 2 GeV, and k/a reh = 10 2 GeV.The dashed vertical lines set the horizon exit (k = aH), reheating (a = a reh ) and the onset of the oscillations of the background field (m = H).
dS epoch With a constant Hubble rate H = H inf , the conformal time is related to the scale factor as τ = − 1 aH inf .(B.1)

1 2 . 2 + 4 √ 4 √ 2µD λ− 1 . (B. 13 ) 7 B
(B.11)This can be checked by noting that for an argument w with a sufficiently large amplitude and s ̸ = 0, 1, 2, • • • , the parabolic cylinder function is approximated byD s (w) ∼ e −w 2 4 w s − √ 2π Γ(−s) e iπs+ w 2 4 w −s−1 , (B.12) for π/4 < arg w < 5π/4.For |arg w| < 3π/4, the approximate expression is given only by the first term in the right hand side.Hence the amplitudes of the coefficients A k and B k in (B.10) match with those of the Bogoliubov coefficients α k and β k in the asymptotic future.By matching the mode functions (B.4) and (B.10), as well as their time derivatives u ′ k at a = a reh , the coefficient B k can be obtained as, ν+ iµ)D −λ + e −i 3π 2µ D −λ+1 iD −λ+1 D λ−1 + D −λ D λ + e −i πHere the arguments of the parabolic cylinder functions are e iπ/4 √ 2µ for D −λ and D −λ+1 , while e i3π/4 √ 2µ for D λ and D λ−1 .This expression for B k can by simplified by first expanding in small κ, and then further simplifying the leading term by considering a small µ.This gives,

7 1 2(
We use that the Hankel function of the first kind at w → 0 asymptotes to H Re(s) > 0. We also use D

Figure 10 :
Figure 10: Time evolution of |β k | 2 in a quasi de Sitter inflationary background, for different initial background field values.The parameters used are m φ = 10 8 GeV, m = 10 2 GeV and k/a reh = 10 2 GeV.