We study a realization of the inflationary scenario where the PecceiQuinn PQ symmetry is spontaneously broken during inflation facilitated by its nonminimal coupling to gravity. This results in effectively a twofield inflation the early stage is driven by an inflaton field with the PQ symmetry intact and the later stage is driven by the PQ scalar after its effective mass becomes tachyonic causing destabilization from the origin. The nonminimal coupling serves the dual purpose of restoring the PQ symmetry during early inflation and flattening the PQ potential posttachyonic shift allowing for continued slowroll. We analyze the inflationary background solutions and scalar perturbations which are amplified at small scales via significant isocurvature perturbations generated near the symmetrybreaking epoch. These perturbations lead to secondorder gravitational waves detectable by nextgeneration spacebased experiments.
['Cosmic inflation', 'Inflationary universe', 'Quantum cosmology', 'Non-standard theories of gravity', 'Gravitational waves', '319', '784', '1313', '1118', '678', 'High Energy Physics - Theory', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology']
Mesoinflationary PecceiQuinn Symmetry Breaking with Nonminimal Coupling
{'Meso-inflationary Peccei-Quinn symmetry breaking with non-minimal coupling': 'Yermek Aldabergenov, a, 1 Ding Ding, a, 2 Wei Lin, a, 3 Yidun Wan a,b, 4  \n- a State Key Laboratory of Surface Physics, Center for Astronomy and Astrophysics, Department of Physics, Center for Field Theory and Particle Physics, and Institute for Nanoelectronic devices and Quantum computing, Fudan University, 2005 Songhu Road, Shanghai 200433, China\n- b Shanghai Research Center for Quantum Sciences, 99 Xiupu Road, Shanghai 201315, China', 'Abstract': 'We study a realization of the inflationary scenario where the Peccei-Quinn (PQ) symmetry is spontaneously broken during inflation, facilitated by its non-minimal coupling to gravity. This results in effectively two-field inflation: the early stage is driven by an inflaton field with the PQ symmetry intact, and the later stage is driven by the PQ scalar after its effective mass becomes tachyonic, causing destabilization from the origin. The non-minimal coupling serves the dual purpose of restoring the PQ symmetry during early inflation and flattening the PQ potential post-tachyonic shift, allowing for continued slow roll. We analyze the inflationary background solutions and scalar perturbations, which are amplified at small scales via significant isocurvature perturbations generated near the symmetry-breaking epoch. These perturbations lead to second-order gravitational waves, detectable by next-generation space-based experiments.', '1 Introduction': 'A natural solution to the strong CP problem of QCD is provided by the existence of a pseudoNambu-Goldstone boson, or axion, which couples to the gluon topological density, originating from the spontaneous breaking of a global abelian symmetry known as the Peccei-Quinn (PQ) symmetry. Additionally, the axion is one of the leading candidates for dark matter [1-6], motivating cosmological and astrophysical studies in this direction.  \nAlong this line, we realize a meso-inflationary Peccei-Quinn Symmetry Breaking (PQSB) during the middle of observable inflation, facilitated by the non-minimal coupling of the PecceiQuinn field to gravity. This has two key consequences:  \n- 1. A two-stage inflation: the first stage driven by a PQ-neutral inflaton field with restored PQ symmetry, and the second stage driven by the PQ scalar after it becomes tachyonic, allowing for continued slow-roll inflation due to the flattening of the PQ potential by the non-minimal coupling.\n- 2. Amplification of scalar perturbations near the PQSB epoch, sourcing second-order gravitational waves (GW) potentially detectable by next-generation space-based experiments.  \nFurthermore, our framework adheres to general constraints on axion mass and decay constant, providing a broad perspective without adhering to specific axion models or production mechanisms. Although for demonstration we will use the Starobinsky model, our approach works for generic inflaton scalar potentials.  \nBackground and motivation . The origin of the axion is described by the Peccei-Quinn theory, where a complex scalar with an appropriate potential spontaneously breaks the PQ symmetry U (1) PQ at some high energy scale. Aside from the QCD axion, other axions and  \naxion-like particles are expected to arise in string theory due to topological properties of extra dimensions [7]. These axions and the respective Peccei-Quinn symmetries can play important roles in the evolution of the universe, for example affecting inflation, reheating, and dark matter abundance.  \nWhen it comes to the QCD axion, its mass can be constrained from combined astrophysical and QCD effects as [8]  \n10 -11 eV ≲ m a ≲ 10 -2 eV , (1)  \nwhich translates into constraints on the decay constant,  \n10 9 GeV ≲ f ≲ 10 17 GeV , (2)  \nby using the one-to-one relation [9]  \nm a = 5 . 7 × 10 -6 eV ( 10 12 GeV f ) . (3)  \nTighter constraints can be obtained when considering specific models and production mechanisms of axion dark matter. We will use the constraints (1) and (2) as guidelines, without referring to a particular axion model or production mechanism, to keep the framework as general as possible.  \nIn the context of inflationary cosmology, PQSB is usually assumed to occur either before observable inflation (pre-inflationary PQSB) or after inflation (post-inflationary PQSB). In the former case, cosmic strings that form during the PQ phase transition are diluted by inflation, and the misalignment mechanism is regarded as the main source of axion dark matter abundance. This scenario, however, can be problematic due to the presence of the (massless) axion during the horizon exit of the CMB scales. This presence of axion can potentially lead to large isocurvature perturbations. In the post-inflationary PQSB, the cosmic strings, and subsequently domain walls, form after inflation, and can be dangerous because they can overclose the universe. The decay of these cosmic defects also contributes to the axion dark matter abundance in this case, constraining the models. Therefore, it is important to study possible PQ dynamics during and after inflation to understand axion dark matter.  \nRecently, in addition to pre- and post-inflationary PQSB, some attention has gone to a third scenario, where the PQ symmetry is broken in the middle of (observable) inflation, first proposed in [10]. The authors of [11] and [12] studied the production mechanisms and constraints on axion dark matter in this scenario, while in [13, 14] the authors showed that the cosmic string-domain wall network can produce a large stochastic gravitational wave background, potentially relevant for pulsar timing array experiments, as well as primordial black holes (PBH) (see also [15-17] for axion domain walls and their cosmological signatures). A general proposal of how to achieve meso-inflationary PQSB is to consider a time-dependent (e.g. inflaton-dependent) PQ mass that can turn tachyonic at some point during inflation, triggering the symmetry breaking. A scenario where PQ phase transition happens during warm inflation was studied in [18].  \nIn our scenario, inflation is split into two stages: 1) The first stage is driven by a (PQneutral) inflaton field ϕ , while the PQ symmetry is restored by the non-minimal coupling (CMB scales exit the horizon at this stage). 2) The second stage starts after the critical point in the potential is reached, and the PQ scalar ρ acquires a tachyonic mass. Around the critical point, the classical force on the ρ -field gives way to quantum diffusion, which kicks ρ away from the origin. This triggers the PQSB and starts the second slow-roll stage, which is made possible thanks to the flattening of the PQ potential by the non-minimal coupling. The difference with  \nthe usual Higgs/Peccei-Quinn inflationary models is that in our case, the PQ-driven inflation happens around the symmetric point of the potential rather than at large field values. For the inflaton sector, we can, in principle, use any single-field model, but as a benchmark example, we will use the Starobinsky model. This scenario can be described as hybrid inflation [19-21] with a long waterfall regime, where the waterfall field is identified with the PQ radial scalar.  \nThe paper is organized as follows. Section 2 studies the case of single-field Peccei-Quinn inflation around the symmetric point and demonstrate the effect of the non-minimal coupling on the duration of inflation. Section 3 brings together the inflaton and non-minimally coupled Peccei-Quinn Lagrangians, and then studies the resulting two-field inflation, both with analytical estimates and numerically. At this point we introduce the method of stochastic inflation to study the transient stage where the symmetry breaking happens. Section 4 focuses on scalar perturbations and shows how the power spectrum is amplified at small scales. Section 5 concludes the paper and discusses about future works.', '2 Peccei-Quinn inflation around the symmetric point': 'Before studying the two-field dynamics of the inflaton ϕ and PQ scalar ρ , let us consider whether the simplest PQ potential is suitable for inflation around its symmetric point.', '2.1 Minimal coupling': 'The simplest PQ Lagrangian is given by  \n√ -g -1 L = 1 2 R -∂S∂S -1 4 λ ( f 2 -2 SS ) 2 , (4)  \nwhere S is the PQ complex scalar, which can be parametrized as S = 1 √ 2 ρe iσ . The angular mode σ shifts by a real constant under the U (1) PQ rotation and is identified with the QCD axion. The radial part ρ eventually acquires non-vanishing VEV ⟨ ρ ⟩ = f , spontaneously breaking the PQ symmetry.  \nAs usual, to study inflation, we take the FLRW metric ds 2 = -dt 2 + a 2 ( t ) δ ij dx i dx j and homogeneous background field ρ ( t ), and obtain from Einstein equations,  \n3 H 2 = 1 2 ˙ ρ 2 + V , ˙ H = -1 2 ˙ ρ 2 , (5)  \nwhere H ≡ ˙ a/a is the Hubble function. The scalar field (background) equation of motion (EOM) is  \n¨ ρ +3 H ˙ ρ + V ,ρ = 0 , (6)  \nwhere subscript ρ denotes the respective partial derivative. We define the potential and Hubble slow-roll parameters as  \nϵ V ≡ V 2 ,ρ 2 V 2 , η V ≡ V ,ρρ V , ϵ H ≡ -˙ H H 2 , η H ≡ ˙ ϵ H Hϵ H . (7)  \nThe usual slow-roll inflation is characterized by the conditions ϵ V , ϵ H , | η V | , | η H | ≪ 1, under which the potential and Hubble slow-roll parameters are related to each other as  \nϵ V ≃ ϵ H , η V ≃ 2 ϵ H -1 2 η H . (8)  \nInflation can also happen under ultra-slow-roll (USR) regime, where the slope of the potential in Eq. (6) nearly vanishes, rendering the first two terms proportional to each-other. In this  \ncase, by writing η H = 2¨ ρ/ ( H ˙ ρ ) + 2 ϵ H (where we used ϵ H = ˙ ρ 2 / (2 H 2 )), and using ¨ ρ = -3 H ˙ ρ , we get η H ≈ -6 assuming that ϵ H is small. Since during ultra-slow-roll V ,ρ is much smaller than 3 H ˙ ρ , we have the hierarchy ϵ V ≪ ϵ H . This can lead to overestimation of the number of efolds if in the equation  \n∆ N = ∫ ρ 2 ρ 1 dtH = ∫ ρ 2 ρ 1 dρH/ ˙ ρ = ∫ ρ 2 ρ 1 dρ √ 2 ϵ H , (9)  \nwe replace ϵ H with ϵ V in order to analytically estimate ∆ N . Nevertheless, analytical estimates based on ϵ V can still be useful for obtaining an upper bound  \n∆ N ≲ ∫ ρ 2 ρ 1 dρ √ 2 ϵ V . (10)  \nIf regular slow-roll holds between ρ 1 and ρ 2 , this upper bound is saturated. Since a sufficiently small initial value of ρ can lead to vanishing V ,ρ , USR is possible in the case of PQ inflation around the origin of ρ . We note that the parameters of the PQ potential in our scenario are not necessarily fixed by the CMB observations because the CMB scales exit the horizon during the ϕ -driven inflation, when ρ is assumed to be stabilized around zero. Rather, we can expect a certain combination of the inflation and PQ parameters to be fixed by CMB data.  \nLet us now analytically estimate the duration of inflation in the PQ model (4), which yields the following potential slow-roll parameters,  \nϵ V = 8 ρ 2 ( f 2 -ρ 2 ) 2 , η V = -4 f 2 -3 ρ 2 ( f 2 -ρ 2 ) 2 . (11)  \nWe immediately see that when ρ ≪ f and for f ≪ 1 (in Planck units), η V becomes large, while ϵ V can still be small if the initial value of ρ is small enough. We can use ϵ V from (11) to estimate the upper bound on the number of efolds (10). To get the conservative upper bound we can take the final value of ρ at its VEV, ρ 2 = f . The resulting ∆ N is  \n∆ N ≲ f 2 4 ( log f ρ 1 -1 2 ) . (12)  \nHere ρ 1 is the initial value of ρ , which is small, ρ 1 ≪ f . The estimate (12) tells us that the duration of inflation in this model is severely suppressed if we assume f ≪ 1. More concretely, if we take for example a (an optimistically) small value ρ 1 /f = 10 -10 , we get ∆ N ≲ 6 f 2 (we also confirm these estimates by numerically solving the EOM).  \nRemarkably, we find that a simple non-minimal coupling of the PQ field to scalar curvature ( ∼ ρ 2 R ) can dramatically increase the number of efolds of the Peccei-Quinn inflation, and at the same time restore U (1) PQ during the ϕ -driven first stage of inflation, without adding explicit ϕ -ρ interactions to the Lagrangian.', '2.2 Non-minimal coupling': 'The Lagrangian for the non-minimally coupled PQ field is  \n√ -g -1 L = 1 2 A ( SS ) R -∂S∂S -λ 4 ( f 2 -2 SS ) 2 . (13)  \nBy rescaling the metric as  \ng mn → A -1 g mn , (14)  \nwe can bring the action to the Einstein frame,  \n√ -g -1 L = 1 2 R -A -1 ∂S∂S -3 4 A -2 ∂A∂A -λ 4 A -2 ( f 2 -2 SS ) . (15)  \nAs usual, we take a linear (in SS ) form of the non-minimal coupling function,  \nA = µ 2 +2 ξSS = µ 2 + ξρ 2 , (16)  \nwhere µ is a real parameter which can be eliminated by requiring that A = M 2 P at the vacuum, or in Planck units,  \n⟨ A ⟩ = µ 2 + ξf 2 = 1 , (17)  \nso that µ 2 = 1 -ξf 2 . Then we can write  \nA = 1 -ξ ( f 2 -ρ 2 ) , (18)  \nand the Lagrangian becomes (ignoring the axion),  \n√ -g -1 L = 1 2 R -1 2 G ρρ ∂ρ∂ρ -λ ( f 2 -ρ 2 ) 2 4[1 -ξ ( f 2 -ρ 2 )] 2 , (19)  \nwhere the field-space metric of ρ , denoted G ρρ , is  \nG ρρ = 1 -ξf 2 +(1 + 6 ξ ) ξρ 2 [1 -ξ ( f 2 -ρ 2 )] 2 . (20)  \nThe canonical parametrization of ρ is given by the solution to  \nd ˜ ρ dρ = √ G ρρ , (21)  \nwhere ˜ ρ is the canonical scalar.  \nThe second derivative of the scalar potential of Eq. (19) at ρ = 0 is V ,ρρ ≃ -λf 2 / (1 -ξf 2 ) 3 . This shows that large negative ξf 2 flattens the potential while keeping the negative sign of V ,ρρ and the correct sign of the kinetic term of ρ around ρ = 0.  \nHaving taken into account the canonical parametrization (21), the potential slow-roll parameters become  \nϵ V = V 2 , ˜ ρ 2 V 2 = V 2 ,ρ 2 G ρρ V 2 = 8 ρ 2 ( f 2 -ρ 2 ) 2 ( A +6 ξ 2 ρ 2 ) ≃ 8 ρ 2 f 4 (1 -ξf 2 +6 ξ 2 ρ 2 ) . (22)  \nand  \nη V = V , ˜ ρ ˜ ρ V = V ,ρρ G ρρ V -∂ ρ G ρρ V ,ρ 2 G 2 ρρ V ≃ -4(1 -ξf 2 ) 2 f 2 (1 -ξf 2 +6 ξ 2 ρ 2 ) 2 , (23)  \nwhere the approximation ρ ≪ f simplifies the expressions, and A is given by (18).  \nFor | ξ | f 2 ≫ 1, ϵ V is suppressed compared to the minimally coupled PQ model. For η V in Eq. (23) the situation is more interesting. At ρ = 0 (and when 6 | ξ | ρ 2 ≪ f 2 ), it reduces to η V = -4 /f 2 , which is the same value as in the minimally coupled PQ model. Nonetheless, as ρ moves away from the origin, it can be quickly suppressed once 6 ξ 2 ρ 2 becomes much larger than | ξ | f 2 (this can happen while ρ/f is still small because we assume | ξ | ≫ f -2 ≫ 1). Therefore,  \nwe can anticipate a larger ∆ N in this case if the non-minimal coupling is sufficiently large. Indeed, we find that the number of efolds between some small ρ 1 and ρ 2 = f is  \n∆ N ≃ ∫ ˜ ρ 2 ˜ ρ 1 d ˜ ρ √ 2 ϵ V = ∫ ρ 2 ρ 1 dρ √ G ρρ 2 ϵ V ≃ f 2 4 ( log f ρ 1 -1 2 ) + 3 4 [ | ξ | f 2 -log(1 + | ξ | f 2 ) ] , (24)  \nwhere we used ρ 1 ≪ f and ξ = -| ξ | . As can be seen, the number of efolds is roughly proportional to | ξ | f 2 for large | ξ | f 2 . This opens up a possibility of meso-inflationary PQ symmetry breaking, where the symmetry is broken for ∆ N ≫ 1 efolds before the end of inflation.', '3 Two-field inflation with Peccei-Quinn phase': "We now consider effective two-field inflation where the usual inflaton scalar ϕ is responsible for the first ∆ N 1 efolds of inflation, while the non-minimally coupled PQ field ρ (either alone or in combination with ϕ ) drives the remaining ∆ N 2 efolds, such that ∆ N 1 +∆ N 2 = 50 ∼ 60. The Lagrangian reads  \n√ -g -1 L = 1 2 A ( SS ) R -1 2 ∂ϕ∂ϕ -∂S∂S -U 1 ( ϕ ) -U 2 ( SS ) , (25)  \nwhere U 1 and U 2 are an inflaton potential and a PQ potential, respectively (we denote Jordanframe potentials by U , and Einstein-frame potentials by V ). As a concrete example we take the Starobinsky potential [22] for the inflaton, and the usual sombrero-shaped potential for the PQ field, 1  \nU 1 = 3 4 M 2 ( 1 -e -√ 2 3 ϕ ) 2 , (26)  \nU 2 = λ 4 ( f 2 -2 SS ) 2 . (27)  \nWe do not consider any direct couplings between ϕ and S in the Jordan frame. If ϕ -S coupling exists in the Jordan frame, it can affect the inflationary history depending on the form of the coupling. For example, it is possible to restore the PQ symmetry during early inflation without the non-minimal coupling to gravity (or with a positive sign of the parameter ξ in (29)), by introducing a suitable ϕ -dependent effective mass for ρ in the Lagrangian. However, the simplest 'sombrero' potential is not flat enought to realize a ρ -driven inflationary phase for subplanckian decay constant f . This would require further modification of the scalar potential. The non-minimal coupling of ρ can naturally achieve both the restoration of U (1) PQ during early inflation, and flattening of the potential in ρ -direction. Then, for simplicity we assume that the ϕ -S couplings in the Jordan frame are absent or suppressed.  \nAfter the Weyl rescaling, from (25) we obtain the Einstein frame Lagrangian,  \n√ -g -1 L = 1 2 R -1 2 A -1 ∂ϕ∂ϕ -3 4 A -2 ∂A∂A -A -1 ∂S∂S -A -2 ( U 1 + U 2 ) , (28)  \nwhere ϕ -S interactions through the function A ( SS ) are now generated, both in kinetic and potential terms.  \nAs discussed in the previous section, we take  \nA = 1 -ξ ( f 2 -ρ 2 ) , (29)  \nsuch that at the vacuum, A = 1. This leads to the effective two-field Lagrangian (writing in terms of ρ = √ 2 | S | )  \n√ -g -1 L = 1 2 R -1 2 G ϕϕ ∂ϕ∂ϕ -1 2 G ρρ ∂ρ∂ρ -A -2 ( U 1 + U 2 ) , (30)  \nwhere the field-space metric reads  \nG ϕϕ = 1 1 -ξ ( f 2 -ρ 2 ) , G ρρ = 1 -ξf 2 +(1 + 6 ξ ) ξρ 2 [1 -ξ ( f 2 -ρ 2 )] 2 . (31)  \nKlein-Gordon equations with non-trivial field-space metric can be written with the help of covariant formalism,  \n□ Φ C +Γ C AB ∂ Φ A ∂ Φ B -G AC ∂ A V = 0 , (32)  \nwhere the indices A, B, C denote (real) scalar fields of the model, G AB is the inverse of the field-space metric, Γ C AB are the corresponding Christoffel symbols, and □ = ∇ µ ∇ µ ( ∂ without explicit indices denote space-time derivatives). For the Lagrangian (30) with the field-space metric (31), we get the background equations  \n¨ ϕ +3 H ˙ ϕ + G ϕϕ ( ∂ ρ G ϕϕ ˙ ϕ ˙ ρ + V ,ϕ ) = 0 , (33)  \n¨ ρ +3 H ˙ ρ + 1 2 G ρρ ( ∂ ρ G ρρ ˙ ρ 2 -∂ ρ G ϕϕ ˙ ϕ 2 +2 V ,ρ ) = 0 , (34)  \nand Friedmann equations,  \n3 H 2 = 1 2 ( G ϕϕ ˙ ϕ 2 + G ρρ ˙ ρ 2 ) + V , (35)  \n˙ H = -1 2 ( G ϕϕ ˙ ϕ 2 + G ρρ ˙ ρ 2 ) . (36)  \nFor future convenience, we write the Einstein-frame scalar potential V = A -2 ( U 1 + U 2 ) as  \nV = 3 M 2 (1 -x ) 2 + λ ( f 2 -ρ 2 ) 2 4[1 -ξ ( f 2 -ρ 2 )] 2 = 3 M 2 (1 -x ) 2 + ˆ λ (1 -ˆ ρ 2 ) 2 4[1 + ˆ ξ (1 -ˆ ρ 2 )] 2 , (37)  \nwhere x ≡ e -√ 2 / 3 ϕ , ξ < 0 (as assumed earlier), and  \nˆ ρ ≡ ρ/f , ˆ ξ ≡ | ξ | f 2 , ˆ λ ≡ λf 4 3 M 2 , (38)  \nsuch that ⟨ ˆ ρ ⟩ = 1 (not to confuse with ˜ ρ which we define as the canonical PQ scalar). The parameter ˆ ξ controls the flatness of the PQ potential, and ˆ λ shows the magnitude of the PQ potential relative to the Starobinsky potential (during early inflation when x and ˆ ρ are small).", '3.1 First stage': 'We define the first stage of inflation as starting from the horizon exit of CMB scales, when ρ is stabilized around zero and ϕ ≫ 1 (or x ≪ 1), and ending when the effective mass of the ρ -field vanishes, at some critical value ϕ c of the inflaton. After ϕ c , the mass-squared of ρ becomes negative, and quantum fluctuations destabilize the classical trajectory from ρ = 0, and the second stage begins. During the transition, slow roll may or may not be violated, depending on whether or not ϵ H becomes large, ϵ H ≥ 1, before we reach the critical point. From the potential (37), by setting V ,ρρ | ρ =0 = 0, we find the critical value,  \nx c ≡ e -√ 2 3 ϕ c = 1 -√ ˆ λ/ ˆ ξ . (39)  \nNevertheless, if the velocity of ϕ becomes non-negligible near ϕ c , the effective (canonical) mass of ρ around ρ = 0 gains the corresponding correction term thanks to the non-trivial field-space metric (see Eq. (34)),  \nm 2 ˜ ρ, eff = ∂ ρ ( G ρρ V ,ρ -1 2 G ρρ ∂ ρ G ϕϕ ˙ ϕ 2 ) | ρ =0 = 3 M 2 ˆ ξ f 2 A 2 0 [ (1 -x ) 2 -ˆ λ ˆ ξ -ϵ ϕ [(1 -x ) 2 + ˆ λ ] 2(3 -ϵ ϕ ) ] , (40)  \nwhere ϵ ϕ ≡ ϵ H | ˙ ρ =0 = G ϕϕ ˙ ϕ 2 / (2 H 2 ) as the ϕ -driven slow-roll parameter relevant for the first stage of inflation, G ϕϕ = G ρρ = 1 /A 0 (with A 0 = 1 + ˆ ξ ), and H 2 = V/ (3 -ϵ ϕ ) at ρ = ˙ ρ = 0. The effective critical point in time, t eff c , can be found by solving m 2 ˜ ρ, eff = 0, and the effective critical value of the inflaton is then given by ϕ eff c = ϕ ( t eff c ), which approaches ϕ c from (39) as ϵ ϕ → 0.  \nDuring early inflation, when x, ϵ ϕ → 0, we have  \nm 2 ˜ ρ, eff ≃ 3 M 2 f 2 A 2 0 ( ˆ ξ -ˆ λ ) , (41)  \nand for the stability of ρ in this regime (and also for the existence of the critical point (39)), we require ˆ ξ > ˆ λ . As ϕ approaches ϕ eff c , the effective mass-squared (40) decreases, and once it vanishes, quantum diffusion in ρ -direction becomes important. We will discuss this regime in the next section and focus on the first stage here.  \nDuring early first stage (at ρ, x → 0), the scale of inflation is given by the approximate potential  \nV ≃ 3 M 2 4 A 2 0 (1 + ˆ λ ) , (42)  \nwhere we can identify three different scenarios depending on the hierarchy of the parameters,  \n(a) ˆ λ ∼ 1 , (b) ˆ λ ≪ 1 , (c) ˆ λ ≫ 1 . (43)  \nIn case (a), both the Starobinsky potential and PQ potential are comparable in magnitude; in case (b), the first stage is dominated by the Starobinsky potential, V ∼ 3 4 M 2 /A 2 0 ; in case (c), it is dominated by the constant term of the PQ potential, V ∼ 1 4 λf 4 /A 2 0 . In each case, the CMB scale amplitude of scalar perturbations [25],  \nP ζ ≃ V 24 π 2 ϵ V ≈ 2 . 1 × 10 -9 , (44)  \ncan be used to fix one relevant parameter, or a combination of parameters. From the potential slow-roll parameters (at ρ = ˙ ρ = 0) at the horizon exit, we can obtain the scalar spectral index and tensor-to-scalar ratio,  \nn s ≃ 1 + 2 η V -6 ϵ V , r ≃ 16 ϵ V , (45)  \nwhich are constrained by observations as [25]  \nn s ≈ 0 . 9649 ± 0 . 0042 (68% CL) , r < 0 . 064 (95% CL) . (46)  \nNote that the slow-roll parameters ϵ V and η V during the first stage are modified compared to the regular Starobinsky inflation due to the non-canonical kinetic term of the inflaton ϕ . In  \nparticular, when ρ = 0, the coefficient of the kinetic term is G ϕϕ | ρ =0 = 1 /A 0 = 1 / (1 + ˆ ξ ). This means that the appropriate slow-roll parameters are given by  \nϵ V = A 0 V 2 ,ϕ 2 V 2 = 4 A 0 x 2 (1 -x ) 2 3[ ˆ λ +(1 -x ) 2 ] 2 , (47)  \nη V = A 0 V ,ϕϕ V = -4 A 0 x (1 -2 x ) 3[ ˆ λ +(1 -x ) 2 ] , (48)  \nand the scalar power spectrum during the first stage is  \nP ζ ≃ 3 M 2 [ ˆ λ +(1 -x ∗ ) 2 ] 3 128 A 3 0 π 2 x 2 ∗ (1 -x ∗ ) 2 , (49)  \nwhere x ∗ is the CMB value of x .  \nCase (a). We begin with the case free of strong hierarchy between the parameters 3 M 2 and λf 4 , i.e. where ˆ λ is close to O (1). Using the CMB value P ζ ≈ 2 . 1 × 10 -9 , from (49) we get  \nM 2 A 3 0 ≃ 8 . 84 × 10 -7 x 2 ∗ (1 -x ∗ ) 2 [ ˆ λ +(1 -x ∗ ) 2 ] 3 . (50)  \nBy combining the above expression and the definition ˆ λ ≡ λf 4 / (3 M 2 ), we can write  \nf 4 ≃ 2 . 65 × 10 -6 ˆ λA 3 0 x 2 ∗ (1 -x ∗ ) 2 λ [ ˆ λ +(1 -x ∗ ) 2 ] 3 . (51)  \nWe can then relate permitted values of n s and r to the axion decay constant f by using their slow-roll expressions (45),  \nn s ≃ 1 -8 A 0 x ∗ (1 -2 x ∗ ) 3[ ˆ λ +(1 -x ∗ ) 2 ] -8 A 0 x 2 ∗ (1 -x ∗ ) 2 [ ˆ λ +(1 -x ∗ ) 2 ] 2 , (52)  \nr ≃ 64 A 0 x 2 ∗ (1 -x ∗ ) 2 3[ ˆ λ +(1 -x ∗ ) 2 ] 2 , (53)  \nand drawing parametric plots of n s ( x ∗ ) versus f ( x ∗ ), and r ( x ∗ ) versus f ( x ∗ ) for the range 0 < x ∗ < 1. See Fig. 1 (left and center). Fig. 1 (right) also plots possible values of n s ( x ∗ ) against r ( x ∗ ). In these figures, we fix ˆ λ = 1, and take A 0 = 1 , 10 , 100 for instances ( A 0 = 1 corresponds to ˆ ξ = 0, i.e. the minimal coupling limit). The main takeaway from these plots is that even without studying the second stage, we can conclude that f is bounded from below as ≳ 10 -2 in case (a). In addition, from the n s -r plot of Fig. 1, we have learned that large A 0 (viz large ˆ ξ ) can reduce the tensor-to-scalar ratio significantly.  \nCase (b). We proceed to the small ˆ λ case where the potential during the first stage can be approximated by that of the Starobinsky,  \nV ≃ 3 M 2 4 A 2 0 (1 -x ) 2 . (54)  \nIn this case, from the amplitude (44), we can obtain (at the leading order in ˆ λ )  \nM ≃ 0 . 94 × 10 -3 A 3 / 2 0 x ∗ (1 -x ∗ ) 2 . (55)  \nFigure 1: Parametric plots of scalar tilt n s ( x ∗ ), tensor-to-scalar ratio r ( x ∗ ), and the axion decay constant f ( x ∗ ) in the Starobinsky-Peccei-Quinn hybrid inflation. Plots cover the whole range 0 < x ∗ < 1 ( x ∗ ≡ e -√ 2 / 3 ϕ ∗ corresponds to the horizon exit of the CMB scale) for A 0 = 1 , 10 , 100 and ˆ λ = 1. The vertical blue line in the left plot represents the values of ϕ ∗ close to 0 . 6, leading to large r of around 0 . 9, and can be ignored (analogous orange and green vertical lines are outside of the plot, corresponding to f ≳ 1).  \n<!-- image -->  \nThe number of efolds (9) during the first stage, i.e. between x ∗ and x c (the critical point) is  \n∆ N 1 ≃ 3 4 A 0 ( x c -x ∗ x c x ∗ +log x ∗ x c ) . (56)  \nThe Starobinsky inflation is of the large-field type, i.e. horizon exit of the CMB scales occurs when ϕ is large, or x ∗ ≪ 1, in which case we can approximate the above equation as  \n∆ N 1 ≃ 3 4 A 0 x ∗ , (57)  \nwhere x c is ignored because it gives only subleading contributions. After substituting x ∗ from (57) into (55), we obtain √  \nM ≈ 7 . 05 × 10 -4 A 0 ∆ N 1 . (58)  \nSimilarly, for n s and r we find  \nn s ≃ 1 -2 ∆ N 1 , r ≃ 12 A 0 ∆ N 2 1 , (59)  \nby using (52) and (53). The difference from the predictions of the usual Starobinsky model can be seen in the expression for r containing A 0 = 1+ ˆ ξ , which comes from the non-canonical kinetic term of ϕ . For large A 0 , tensor-to-scalar ratio is suppressed, in line with our results in Fig. 1. If the second stage is long, e.g. O (10) efolds, ∆ N 1 will be reduced from the usual 50 ∼ 60 efolds, and the value of n s will also decrease from standard prediction of the Starobinsky inflation. This can be addressed by modifying the Starobisnky potential 2 , or by considering a different potential that predicts a larger n s . One such example is the axion monodromy-type model with V ∼ [ (1 + ϕ 2 ) 1 / 4 -1 ] . Using this as the Jordan frame potential (26), we get approximately n s ≃ 1 -5 / (4∆ N 1 ). Nevertheless, we do not know a priori the value of ∆ N 1 , and so the best-fitting inflaton potential should be chosen depending on the results from the second stage of inflation, which will depend on the parameter choice. Modifying the  \npotential will affect the results quantitatively, but the qualitative behavior of the solutions are expected to hold, at least for large field inflation. For instance, the initial conditions for the ρ -driven second stage of inflation (derived from quantum diffusion near the symmetric point) can vary depending on the form of the inflaton potential. Yet, the slow-roll solution of the second stage is largely insensitive to these initial conditions due to the Hubble friction in Eq. (34). For concreteness, we will continue with the Starobinsky potential as a reference model, which will also allow us to use numerical methods when necessary.  \nWe define the end of the first stage as the time when x = x eff c , viz. when the effective ρ -mass vanishes, and the slow-roll regime may or may not hold at this point depending on the parameter choice. We can consider the expression  \nm 2 ˜ ρ, eff ∝ (1 -x ) 2 -ˆ λ/ ˆ ξ (1 -x ) 2 + ˆ λ -ϵ ϕ 2(3 -ϵ ϕ ) , (60)  \nobtained from (40) by dividing the expression in the square brackets by (1 -x ) 2 + ˆ λ . Upon slowroll violation, e.g. at ϵ ϕ = 1, we check the sign of (60): If it is negative (positive), the critical point is reached before (after) slow-roll violation. In case (b) ( ˆ λ ≪ 1), if ˆ λ, ˆ λ/ ˆ ξ ≪ (1 -x ) 2 , the mass-squared is still positive at ϵ ϕ = 1, so slow-roll is broken before the critical point. If ˆ λ/ ˆ ξ is close to (1 -x ) 2 , slow-roll may hold at t eff c .  \nCase (c). In this case, ˆ λ ≫ 1, i.e. the potential is dominated by the PQ constant term ∼ λf 4 , which suppresses the slow-roll parameters in comparison to case (b), as is clear from (47) and (48). The positivity of the effective ρ mass-squared (41) during the first stage requires ˆ ξ > ˆ λ , causing A 0 ∼ ˆ ξ ≳ ˆ λ for ˆ λ ≫ 1, i.e., A 0 cannot be small. Hence, both potential slow-roll parameters can be small only if x ≪ 1.  \nFrom the power spectrum (44) at the CMB scale, we can obtain the condition (for ˆ λ ≫ 1)  \n( λf 4 ) 3 A 3 0 M 4 ≃ 2 . 4 × 10 -5 x 2 ∗ (1 -x ∗ ) 2 , (61)  \nand the number of efolds becomes  \n∆ N 1 ≃ 3 ˆ λ 4 A 0 [ x c -x ∗ x c x ∗ +log x c (1 -x ∗ ) x ∗ (1 -x c ) ] . (62)  \nEquation (62) can re-express M in terms of ∆ N 1 (after using ˆ λ = λf 4 / (3 M 2 )) and casts (61) as  \nf ≃ 0 . 035 A 1 / 4 0 √ ∆ N 1 λ 1 / 4 √ x ∗ (1 -x ∗ ) [ x c -x ∗ x c x ∗ +log x c (1 -x ∗ ) x ∗ (1 -x c ) ] 1 2 . (63)  \nFor large ˆ λ , we can approximate n s as  \nn s ≃ 1 -8 A 0 ˆ λ x ∗ (1 -2 x ∗ ) , (64)  \nwhich, after eliminating ˆ λ by (62), can be written as  \nx ∗ (1 -2 x ∗ ) [ x c -x ∗ x c x ∗ +log x c (1 -x ∗ ) x ∗ (1 -x c ) ] ≃ 1 -n s 2 ∆ N 1 . (65)  \nBy combining this equation with (63) and using n s ≈ 0 . 9649 3 , we get  \nf ≃ 4 . 64 × 10 -3 A 1 / 4 0 λ 1 / 4 √ 1 -x ∗ 1 -2 x ∗ , (66)  \nwhere it becomes clear that f is bounded from below,  \nf ≳ 4 . 64 × 10 -3 A 1 / 4 0 λ 1 / 4 , (67)  \nassuming x ∗ < 1 / 2, as otherwise we will get a blue-tilted spectrum. Consistency with the astrophysical upper bound (2) ( f ≲ 0 . 1 in Planck units) requires that  \nA 0 = 1 + ˆ ξ ≲ 10 4 , (68)  \nfor λ = 0 . 1 or so. On the other hand, ˆ ξ is necessarily positive, rendering A 0 ≥ 1. Thus, f in (67) cannot be smaller than O (10 -3 ) in this case (we do not consider large λ because it would break perturbative regime).', '3.2 Transient stage: quantum diffusion and classicalization': "At the end of stage one, there is a brief transient stage where quantum effects become nonnegligible for the radial part of the PQ field ρ . To see this, we first recall that ρ -field's value is effectively zero during the first stage due to the U (1) PQ symmetry restored by the non-minimal coupling. As the ρ -field's effective mass becomes tachyonic, the ρ -field now sits at an unstable critical point, where the uncertainty principle dictates that quantum fluctuations will push it down to the stable minima at ρ = f . The usual classical equations of motion for the ρ -field are not sufficient to evolve the system to the subsequent stages. Instead, a framework including the effects of quantum backreaction is needed.  \nLet us make our statements more concrete. When the effective mass (and the slope of the potential) of the ρ -field vanishes at the end of stage one, the field enters into the USR regime 4 , where quantum diffusion becomes important. Consider first the classical EOM (34) for the ρ -field. At the end of stage one, the classical force ρ -field feels from the potential vanishes while ρ ∼ 0. Force terms including the factors ∂ ρ G ρρ = 0 and ∂ ρ G ϕϕ = 0 also vanish. Thus, the ρ -velocity is suppressed due to the Hubble friction  \n¨ ρ +3 H ˙ ρ = 0 ⇒ ˙ ρ ∼ e -3 Ht . (69)  \nThis suppression implies that classical evolution by itself favors the position in field space where ρ ∼ 0, possibly extending inflation for too long and producing curvature perturbations that are too large (after inflation ends in the ϕ -direction). Aside from velocity suppression, the common assumption (eg. in [27, 28]) of neglecting ¨ ρ may not be taken for granted, and in our case may not be valid, even around ρ = 0. In particular, for our two-field model, the  \nsecond (Hubble) slow-roll parameter η H is not close to -6 (viz, the typical single-field USR value) but instead  \nη H ∼ 2 ϵ H + O ( ¨ ϕ/ (3 H ˙ ϕ )) (70)  \nwhen ρ, ˙ ρ ∼ 0. If standard slow-roll still holds in the ϕ -direction, such that ¨ ϕ ≪ 3 H ˙ ϕ , the second slow-roll parameter satisfies | η H | ≪ 1 and stays small. Nonetheless, slow-roll can be temporarily violated for a certain parameter region.  \nThe non-standard evolution described above signals the need for a quantum description beyond the traditional slow-roll attractor description. In principle, the full description we need is quantum gravity because quantum diffusion occurs via backreaction of the small-scale modes on the (locally) homogeneous background field, which couples to geometry through the Friedmann equations. Fortunately, we expect the homogeneous part of the ρ -field to deviate only slightly from a classical state, provided the transient stage is brief (this briefness will be verified later). So, there is no need for a full quantum description. When this is indeed the case, the effective description provided by stochastic inflation [29,30] can be applied.  \nIn the framework of stochastic inflation, quantum diffusion occurs when small-scale modes collectively backreact on the (locally) homogeneous part of the field as they get stretched by (quasi) de Sitter expansion. The classical equations of motion for the background field are thus modified with noise terms arising from the quantum fluctuations of these small-scale modes. The homogeneous field becomes stochastic by inheriting the statistics of the quantum noise terms through Langevin-like equations. The statements can formally be made rigorous starting from canonical quantization and then obtaining an effective phase-space distribution of ρ that mimics some (but not all) of the features of its true quantum wavefunction.  \nFor our purpose, we wish to use the stochastic formalism to obtain the classical values of ρ and its conjugate momentum after the brief transient stage of quantum diffusion. These values serve as initial conditions for the subsequent classical evolution, which we call stage II. Multi-stage evolution containing a brief quantum diffusion stage in between two classical phases is common in mutli-field models, especially in hybrid inflation ( [28,31-34] for example). But because the classicalization process in general is not fully understood, the transition from a quantum diffusion stage (probabilistic) to the subsequent classical stage (deterministic) is difficult to implement, especially numerically. In some cases, one manually selects a set of classical values after the quantum diffusion phase based on order-of-magnitude estimates, without having to worry about stitching stochastic evolution to classical deterministic evolution. In other cases where quantum diffusion is treated more carefully with the stochastic formalism, one sometimes assumes an effectively instantaneous classicalizaton time for smoother numerical implementation but runs into the danger of a priori excluding the parameter space with prolonged quantum diffusion, which could affect primordial black hole abundance [35] if the ρ -field dominates the inflaton dynamics. In the following, we will use effective equations of motion for the quantum diffusion stage that are derived rigorously with canonical quantization and the stochastic framework. We will not ignore any time-derivative (or slow-roll) terms by default, as our model does not follow all the usual slow-roll assumptions. Lastly, we will employ a classicalization criteria, which helps determine the classicalization time, allowing for a smooth transition to the subsequent classical evolution.", '3.2.1 Quantum diffusion': "To describe quantum diffusion in the phase-space of the ρ -field we introduce the momentum field (we start with physical time in the Lagrangian (30) and then use d dt = H d dN ):  \nπ ρ ≡ 1 a 3 H 2 ∂ L ( t ( N )) ∂ρ ,N = G ρρ ρ ,N , π ρ = G ρρ π ρ = ρ ,N , (71)  \nwhere the subscript ' N ' denotes the derivative w.r.t. the efold number, and G ρρ raises the index (field metric is diagonal) and leads to the contravariant momentum. By this definition, from the ρ EOM (34), we obtain the equation system  \nρ ,N = π ρ , π ρ ,N = -( 3 -ϵ H + 1 2 G ρρ ∂ ρ G ρρ π ρ ) π ρ + G ρρ ( 1 2 ∂ ρ G ϕϕ ϕ 2 ,N -V ,ρ /H 2 ) . (72)  \nBy splitting ρ and π ρ into IR and UV modes (see Appendix A), we can obtain the stochastic equations  \nρ ,N = π ρ + ξ ρ , π ρ ,N = -( 3 -ϵ H + 1 2 G ρρ ∂ ρ G ρρ π ρ ) π ρ + G ρρ ( 1 2 ∂ ρ G ϕϕ ϕ 2 ,N -V ,ρ /H 2 ) + ξ π ρ , (73)  \nwhere ρ and π ρ should be understood as IR modes, in which the equations are in general nonlinear, while they have been linearized in the UV modes, which in turn generate the stochastic noise terms ξ ρ and ξ π ρ defined as  \nξ ρ ≡ -∫ d 3 k (2 π ) 3 / 2 ∂ N W ( k k s ) ( e -i kx ρ k ˆ a k +h . c . ) , (74)  \nξ π ρ ≡ -∫ d 3 k (2 π ) 3 / 2 ∂ N W ( k k s ) ( e -i kx π ρ k ˆ a k +h . c . ) , (75)  \nwhere W ( x ) = θ ( x -1) is a step function that selects out large-wavenumber modes. To avoid confusion, ξ and ˆ ξ without any indices always denote the non-minimal coupling, while ξ with subscript ρ or π ρ refers to the stochastic noise.  \nNear the origin ρ = 0, we can further linearize Eqs. (73) w.r.t. ρ and π ρ , but keep nonlinearities in ϕ and its derivative, treating them as a classical dynamical background (possibly non-slow-roll when ˙ ϕ is large),  \nρ ,N = π ρ + ξ ρ , π ρ ,N = -(3 -ϵ ϕ ) π ρ -H -2 m 2 ˜ ρ, eff ρ + ξ π ρ , (76)  \nwhere G ρρ ≃ G ϕϕ ≃ 1 /A 0 at the leading order in ρ , and ϵ ϕ is defined as  \nϵ ϕ ≡ ϵ H | ρ ,N =0 = 1 2 A -1 0 ϕ 2 ,N . (77)  \nThe effective mass term can be written as  \nm 2 ˜ ρ, eff H 2 = 4 f 2 ˆ ξ (3 -ϵ ϕ ) { (1 -x ) 2 -ˆ λ/ ˆ ξ (1 -x ) 2 -ˆ λ -ϵ ϕ 2(3 -ϵ ϕ ) } , (78)  \nwhere the evolution of ϕ is given by its classical EOM.  \nIf the noise is assumed to be Gaussian, it suffices to consider quadratic expectation values. From (76), we obtain their equations of motion:  \n⟨ ρ 2 ⟩ ,N = 2 Q + ⟨ ξ 2 ρ ⟩ , Q ,N = -m 2 ˜ ρ, eff H 2 ⟨ ρ 2 ⟩ -(3 -ϵ ϕ ) Q + ⟨ π ρ 2 ⟩ + ⟨ ξ ρ ξ π ρ ⟩ + ⟨ ξ π ρ ξ ρ ⟩ , ⟨ π ρ 2 ⟩ ,N = -2 m 2 ˜ ρ, eff H 2 Q -2(3 -ϵ ϕ ) ⟨ π ρ 2 ⟩ + ⟨ ξ 2 π ρ ⟩ , (79)  \nwhere Q ≡ 1 2 ( ⟨ ρπ ρ ⟩ + ⟨ π ρ ρ ⟩ ), and ⟨ ξ I ξ J ⟩ (with I, J = ρ, π ρ ) are quadratic (time-dependent) amplitudes of the noise correlators, defined as  \n⟨ ξ I ( N ) ξ J ( N ' ) ⟩ ≡ ⟨ ξ I ξ J ⟩ δ ( N -N ' ) . (80)  \nIn terms of the Fourier modes, the noise correlators are given by  \n⟨ ξ 2 ρ ⟩ = | ρ k s | 2 2 π 2 k 3 s | 1 -ϵ ϕ | , ⟨ ξ 2 π ρ ⟩ = | π ρ k s | 2 2 π 2 k 3 s | 1 -ϵ ϕ | , ⟨ ξ ρ ξ π ρ ⟩ + ⟨ ξ π ρ ξ ρ ⟩ = ρ k s ( π ρ k s ) ∗ +c . c . 2 π 2 k 3 s | 1 -ϵ ϕ | . (81)  \nThe noise amplitudes can be derived by using the Fourier mode solution ρ k . In the slow-roll limit, we can approximate it as 5  \nρ k ≃ √ A 0 e -ikτ a √ 2 k (1 + iaH/k ) , π ρ k = ρ k,N , (82)  \nwhere √ A 0 appears after the canonical normalization of ρ . This leads to the quadratic noise amplitudes,  \n⟨ ξ 2 ρ ⟩ ≃ A 0 H 2 4 π 2 , ⟨ ξ 2 π ρ ⟩ ≃ A 0 H 2 4 π 2 ( ϵ ϕ + s 2 ) 2 , ⟨ ξ ρ ξ π ρ ⟩ + ⟨ ξ π ρ ξ ρ ⟩ ≃ -A 0 H 2 2 π 2 ( ϵ ϕ + s 2 ) . (83)  \nSince s, ϵ ϕ ≪ 1 during slow-roll, ⟨ ξ 2 ρ ⟩ becomes the dominant noise term, and the rest can be ignored. In the cases where slow-roll is violated near N eff c , expressions (83) can be treated as order-of-magnitude estimates.  \nWe numerically solve the system (79) starting from N eff c , with vanishing initial conditions for ⟨ ρ 2 ⟩ , Q , and ⟨ π ρ 2 ⟩ . We then compute the duration ∆ N q of the quantum diffusion stage at the end of which, the solution 'classicalizes', viz., the stochastic noise term becomes negligible compared to the field momentum,  \n⟨ π ρ 2 ⟩ = κ 2 ⟨ ξ 2 ρ ⟩ , (84)  \nfor some constant κ ≫ 1 (for numerical solutions we will use κ = 100). After the solution √ ⟨ ρ 2 ⟩ is found, its value at the time N eff c +∆ N q (as well as its velocity) is taken as the initial condition for the second classical stage of inflation.  \nFigure 2: Solution to stochastic equations (79) (here π stands for π ρ ). The right plot shows how the average momentum becomes much larger than the noise term, by choosing for example κ = 100.  \n<!-- image -->  \nFigure 2 showcases a solution to stochastic equations (79) for the parameters ˆ λ = 1, M = 8 . 62 × 10 -5 , λ = 0 . 1 ( f is then derived as f ≈ 0 . 022). Here, we first run the numerical solution to the background equation (33) for ϕ (when ρ = 0) starting from N = 0 until the effective critical point N eff c ≈ 77 . 19 plus e.g. an additional efold, in order to obtain the ϕ -dependent coefficients in stochastic equations (79). The solution to (79) (Fig. 2, left) shows that ⟨ π ρ 2 ⟩ value starts to grow first, followed by Q and ⟨ ρ 2 ⟩ . The plot on the right indicates the 'classicalization' time when the momentum average crosses κ √ ⟨ ξ 2 ρ ⟩ , where we set κ = 100. This yields ∆ N q ≈ 0 . 024 in the current example. From this point, we proceed to the second classical stage of inflation.", '3.3 Second stage and numerical background solutions': "Here, we obtain full numerical trajectories for the background system (33)-(37). To summarize the strategy, we divide the inflationary solution into three parts: 1) the first classical solution running from N = 0 until the effective critical point N eff c ; 2) the quantum diffusion stage lasting ∆ N q efolds, as explained above; 3) the second classical stage taking place between N eff c +∆ N q and until the end of slow-roll, ϵ = 1. For the second classical stage, the ρ initial condition is provided by its quantum diffusion stage, while the ϕ initial condition can be found from its classical EOM (33), since the quantum diffusion of ϕ can be neglected due to its non-vanishing classical velocity (it in fact grows as we approach the critical point).  \nWe set λ = 0 . 1 as a reference value, while the inflaton mass parameter M can be fixed from the CMB scalar amplitude (44) 6 , and we are left with two free parameters f and ξ . As mentioned earlier, it is convenient to trade them in favor of the derived parameters  \nˆ ξ ≡ | ξ | f 2 , ˆ λ ≡ λf 4 3 M 2 . (85)  \nIn our first set of examples, we will take ˆ ξ = 30, with the expectation that it will lead to ∆ N 2 = 20 ∼ 30 efolds of PQ inflation, and vary ˆ λ around unity. The parameter ˆ λ is bounded from above, ˆ λ < ˆ ξ , as required by the positivity of the effective ρ mass (41) (i.e., restoration  \nFigure 3: Hubble function and the slow-roll parameters during the last 55 efolds of inflation. H is normalized by using Ω ≡ λf 4 /A 2 0 . Full parameter sets are shown in Table 1.  \n<!-- image -->  \nof the PQ symmetry) during the first stage of inflation. Figure 3 shows the Hubble function and the Hubble slow-roll parameters ϵ H and η for a range of ˆ λ H , during the last 55 efolds of inflation ( N = 0 corresponds to the horizon exit of CMB scales). Table 1 records the rest of the parameters used, as well as the initial conditions from the quantum diffusion stage. We find that as we increase ˆ λ , δN 2 gradually decreases until we get close to the upper bound ˆ λ = ˆ ξ = 30, at which it starts to rapidly increase. At ˆ λ ≈ 28 . 8, the second stage lasts more than 53 efolds. At the same time the duration of the quantum diffusion stage ∆ N q increases with ˆ λ , and at ˆ λ ≈ 28 . 8 it reaches ∼ 1 . 5 efolds.  \nFrom Fig. 3, we can see the step-like behavior of the Hubble function between the two slow-roll stages for small ˆ λ , while for large ˆ λ , the Hubble function is smooth and flat until the end of the second stage. By looking at the slow-roll parameters (Fig. 3, center and right), we can see spikes and a brief period of oscillations between the two slow-roll phases, where the oscillations are more violent for smaller values of ˆ λ , and slow-roll can be temporarily broken. Nevertheless, the choice ˆ λ ≈ 28 . 8 should be distinguished here because the spike in the slowroll parameters in this case does not indicate the transition to the second stage as we have defined it, but happens during the second stage (the seamless transition to the second stage happens near N = 0).  \nLet us elaborate on this behavior below by comparing the field trajectories for e.g. ˆ λ = 10 with ˆ λ = 28 . 787, as can be seen in Fig. 4, where the classical trajectories are shown in blue, and the short quantum diffusion stage in red. The plot on the left depicts the inflationary trajectory for ˆ λ = 10: the start of the second inflationary stage (right after the red portion of the trajectory) is followed by a short transition period where the ϕ -ρ trajectory is diagonal. Subsequently, we can see a few oscillations of ϕ and finally the ρ -driven inflation where ϕ = 0. In contrast, in the extreme case of ˆ λ = 28 . 787 (right-side plot), the transition period with the diagonal part of the trajectory is much longer and lasts from N ≈ 0 . 1 until the peak in the slow roll parameters shown in Fig. 3, which coincides with the onset of the ρ -driven stage. In summary, when ˆ λ is close to 30, the second stage can be further divided into the two-field part, and the effectively single-field, ρ -driven part. For smaller values of ˆ λ , the two-field part is short. In our examples, the ρ -mass around the true (Minkowski) vacuum varies between m ˜ ρ ∼ 10 -7 (for ˆ λ = 0 . 001) and m ˜ ρ ∼ 10 -6 (for ˆ λ ≈ 28 . 8), viz., a few orders below the inflaton mass given by M (the effective ρ -mass during the first stage is larger, between 10 -2 and 10 -5 , respectively, so ρ is strongly stabilized).  \nThe field-space trajectories for ˆ λ = 10 , 1 , 0 . 1 can reveal the origin of the oscillations in the slow-roll parameters. Figure 5 shows the trajectories as we decrease ˆ λ , around the transition to the second stage. The plot on the left shows few oscillations of ϕ , while ρ gradually grows until it reaches its VEV ⟨ ρ ⟩ = f . For ˆ λ = 1 (center plot), ϕ oscillates faster, but ρ still increases smoothly. For ˆ λ = 0 . 1 (right plot) the oscillations of ϕ further increase, but now ρ  \nTable 1: Parameters sets (for ˆ ξ = 30) and the ρ initial conditions (ˆ ρ ≡ ρ/f ) obtained from the stochastic equations. The last two rows show the duration of the quantum diffusion stage ∆ N q and the second classical stage of inflation ∆ N 2 .  \nFigure 4: Inflationary trajectories on the scalar potential for ˆ λ = 10 (left) and ˆ λ = 28 . 787 (right). Ω ≡ λf 4 /A 2 0 . Red portions represent the transient stage.  \n<!-- image -->  \nalso starts to oscillate around the origin (see Fig. 6), which results in a 'chaotic' behavior of the trajectory where ρ goes from positive to negative values before the trajectory stabilizes around ϕ = 0, and ρ finally reaches its (negative) VEV. Note that the negative values of ρ should be unphysical by the definition of ρ as the radial direction of the complex PQ field S . Nevertheless, these negative excursions of ρ can be understood as the angular change of S , i.e., the change in the axion value after the PQSB. This can be seen by parametrizing S as real and imaginary parts, S = S 1 + iS 2 . Since the EOM are U (1) PQ -symmetric, this yields identical equations for the scalars S 1 and S 2 , and the negative values in the trajectory of ρ can be understood as negative values of both S 1 and S 2 (we confirm this by numerically solving the system with scalars S 1 and S 2 for the parameter choice ˆ λ = 0 . 1 in Fig. 5).  \nTo see how ˆ ξ controls the duration ∆ N 2 of the second stage, we plot ∆ N 2 as a function of ˆ ξ for a given ˆ λ in Fig. 7. For a given ˆ λ , ˆ ξ is bounded as ˆ ξ > ˆ λ by the positivity of the ρ mass-squared (41) during the first stage, and f is varied between ∼ 10 -3 and ∼ 10 -2 (increases with ˆ λ ) in the given examples. Almost linear dependence can be seen when ˆ ξ ≫ ˆ λ , as expected from our analysis of the single field (non-minimally coupled) PQ inflation in Section 2.2 (see Eq. (24) for the corresponding number of efolds). Nonetheless, when ˆ ξ approaches ˆ λ from above, a multi-field effect is seen: ∆ N 2 undergoes steep growth, which is steeper for smaller values of ˆ λ . As a consequence, ∆ N 2 has a minimum for a given ˆ λ , and this minimum increases with ˆ λ . We have seen this sudden growth of ∆ N 2 in the case of ˆ λ ≈ 28 . 8 and ˆ ξ = 30, and the multi-field trajectory responsible for the large ∆ N 2 can be seen in Fig. 4 (right plot).  \n<!-- image -->  \n-  \n0.6  \n-  \n0.4  \n-  \n0.2 0.0  \n0.2  \n0.4  \n0.6  \n0.8  \nϕ  \n-  \n0.6  \n-  \n0.4  \n-  \n0.2 0.0 0.2 0.4 0.6 0.8  \nϕ  \nFigure 5: Field-space trajectories for ˆ λ = 10 , 1 , 0 . 1 (from left to right).Figure 6: Oscillations of ϕ and ρ at the start of the second stage, ˆ λ = 0 . 1.  \n<!-- image -->  \nFigure 7: ∆ N 2 as a function of ˆ ξ for a range of ˆ λ .  \n<!-- image -->  \nf  \n/  \nρ  \n0.5  \n0.4  \n0.3  \n0.2  \n0.1  \n0.0  \nf  \n/  \nρ  \n-  \n0.1  \n-  \n0.2  \n-  \n0.3  \n-  \n0.4  \n0.1  \n0.0", '4 Scalar perturbations': "Let us now turn to scalar perturbations and their superhorizon evolution due to multi-field dynamics. We will work with gauge-invariant perturbations  \nδ Φ A g = δ Φ A + ˙ Φ A ψ/H , (86)  \nwhere δ Φ A is the perturbation of the Φ A field, and ψ is the perturbation of the spatial part of the (spacetime) metric tensor. In the covariant formalism (see e.g. [36,37]), we can write the EOM for δ Φ A as (from now we will drop subscript ' g ')  \nD 2 t δ Φ A +3 HD t δ Φ A + ( k 2 a 2 δ A B + M A B ) δ Φ B = 0 , (87)  \nwhere the field-space covariant derivative along the time direction D t of a vector V A is  \nD t V A = ˙ V A + ˙ Φ B Γ A BC V C , (88)  \nand the effective mass matrix M A B is given by  \nM A B = ∇ A ∇ B V + R A CBD ˙ Φ C ˙ Φ D -a -3 D t ( a 3 ˙ Φ A ˙ Φ B /H ) . (89)  \nFrom the field basis we can switch to the adiabatic-isocurvature basis, which for the two-field case consists of an adiabatic mode Q σ and a single isocurvature mode Q s :  \nQ σ = G AB e A σ δ Φ B , Q s = G AB e A s δ Φ B , (90)  \nwhere the basis unit vectors can be written as  \ne A σ = ˙ Φ A √ G BC ˙ Φ B ˙ Φ C , e A s = D t e A σ √ G BC D t e B σ D t e C σ . (91)  \nIn terms of Q σ and Q s , the power spectra of curvature and isocurvature perturbations are  \nP ζ = k 3 2 π 2 H 2 | Q σ | 2 G AB ˙ Φ A ˙ Φ B , P S = k 3 2 π 2 H 2 | Q s | 2 G AB ˙ Φ A ˙ Φ B . (92)  \nFor the case ˆ ξ = 30 with different values of ˆ λ , Fig. 8 plots the superhorizon changes of the power spectrum by first considering the CMB scale k ∗ (we evolve the perturbation from around 5 efolds before the horizon exit until the end of inflation). For ˆ λ ≳ 1, our testing shows relatively stable behaviour of the k ∗ mode, but for small ˆ λ , significant amplification can be seen after rapid oscillations around the beginning of the second stage (numerical computation time for small ˆ λ also becomes increasingly longer, and prone to stiffness errors). This amplification is likely due to the spikes in the isocurvature effective mass caused by the oscillatory evolution of the fields: [36-39]  \nm 2 iso = G AC e C s e B s m A B +3 G AB D t e A σ D t e B σ , (93)  \nwhere  \nm A B = ∇ A ∇ B V + R A CBD ˙ Φ C ˙ Φ D . (94)  \nFor ˆ λ = 10 , 1 , 0 . 1 (and ˆ ξ = 30), we plot m 2 iso /H 2 in Fig. 9, where the narrow upward (positive) spikes are already evident for ˆ λ = 1. The spikes become increasingly higher and narrower, while for ˆ λ = 0 . 1, the isocurvature mass-squared develops downward (negative) spikes as well.  \nFigure 8: Superhorizon evolution of the power spectra P ζ and P S for the CMB mode k ∗ , normalized ( P ζ ≈ 2 . 1 × 10 -9 ) at the moment of horizon exit.  \n<!-- image -->  \nFigure 9: Effective isocurvature mass for ˆ λ = 10 , 1 , 0 . 1. The values of M are shown in Table 1.  \n<!-- image -->  \nRegarding the normalization of the CMB modes for small ˆ λ (Fig. 8) at ˆ ξ = 30, decreasing the parameter M (and hence f ) in order to decrease the amplitude at the end of inflation (to match P ζ ( k ∗ ) ≈ 2 . 1 × 10 -9 ) is not straightforward because of the non-trivial superhorizon evolution, unlike what one would expect from usual slow-roll models. For example, monotonically decreasing M can result in sudden changes (both positive and negative) in the final amplitude P ζ ( k ∗ , N end ), due to the highly oscillatory transient regime between the two slow-roll stages. Such high oscillations require more sophisticated numerical methods to deal with stiffness and lengthy computation times. On the other hand, because of the sensitivity of the final amplitude on the oscillatory phase, possible particle production (introducing additional friction term) or change in the shape of the scalar potential around the transient stage may significantly modify the superhorizon evolution of the perturbations. For ˆ λ ≳ O (1) and especially for larger values, the results and inflationary trajectories in general are less model dependent (in terms of the inflaton sector) due to the definition of ˆ λ ≡ λf 4 / (3 M 2 ) as the ratio between the PQ potential at ρ = 0 and the inflaton potential at the plateau region (large ˆ λ means that specific shape of the inflaton potential becomes less important).  \nFinally, we address the power spectrum and the amplification of perturbations at small scales. Figure 10 plots P ζ ( k ) for ˆ ξ = 30 and ˆ ξ = 20, and several values of ˆ λ 7 . The peaks in the power spectrum are located approximately at the critical values of k , viz., at the scales exiting the horizon at the critical time t eff c when the effective mass of ρ vanishes (and the amplification of each nearby k -mode happens right after the second classical stage begins). For ˆ λ < 1, we find even larger peaks in the power spectrum, and already at ˆ λ = 0 . 5, the peak reaches around P ζ ∼ 0 . 1 (we do not plot the whole power spectrum for these cases because of long numerical  \nFigure 10: Scalar power spectrum for a range of ˆ λ and ˆ ξ = 30 (left) and ˆ ξ = 20 (right).  \n<!-- image -->  \nFigure 11: Scalar-induced gravitational wave density spectrum for the scalar power spectra of Fig. 10 (with ˆ ξ = 30).  \n<!-- image -->  \ncomputation times). By computing the resulting PBH abundance (using the Press-Schechter formalism [40, 41]), we find overproduction of PBH for ˆ λ ∼ 0 . 5 or smaller, and exclude these values.  \nFor the power spectra shown in Fig. 10, while the PBH abundance is found to be negligible, scalar-induced gravitational waves can be significant. Scalar-induced GW, generated at the second order ( P h ∼ P 2 ζ ) can be computed for a given scalar power spectrum by using the method described in Refs. [42,43].Figure 11 depicts the scalar-induced GW density Ω GW ( h is the reduced Hubble parameter h = 0 . 67) as a function of frequency ν . It shows the stochastic GW background reaching the sensitivities of the upcoming spaceborne interferometers, such as LISA [44], TianQin [45], Taiji [46], and DECIGO [47]. The GW densities are shown for ˆ ξ = 30, corresponding to ∆ N 2 ∼ 22 efolds. For smaller values of ˆ ξ (shorter second stage), the curves move to higher frequencies (for ˆ ξ = 20 they are outside of the shown interferometer sensitivities). Likewise, larger ˆ ξ will prolong the second stage and decrease the peak frequency. We also confirm that first-order primordial GW are not amplified and can be neglected in comparison to scalar-induced GW at the amplified scales.", '5 Discussion': 'In this work, we considered a cosmological scenario where PQ symmetry breaking occurs in the middle of inflation due to the non-minimal coupling to gravity. The early stage of inflation is driven by an inflaton field ϕ , while the non-minimal coupling of the PQ field ρ ensures its stability at the symmetric point. As inflation goes on, the effective PQ mass turns tachyonic, and spontaneous symmetry breaking is triggered by the quantum diffusion  \nof ρ , which we analyze by using stochastic formalism. After a short period of time (typically ∆ N ≪ 1), the momentum/velocity of ρ becomes much larger than the stochastic noise, and classical equations of motion can be used again to describe the second slow-roll stage, with initial conditions provided by the stochastic equations. Since slow-roll resumes during the second stage, the inflationary trajectories are largely insensitive to these initial conditions. The main ingredient that makes inflation in the ρ -direction (i.e., the second stage) possible is the negative non-minimal coupling, ξ < 0, with large enough magnitude. More specifically, the duration of the second stage ∆ N 2 grows almost linearly with ˆ ξ ≡ | ξ | f 2 as shown in Fig. 7, aside from extreme cases where ˆ ξ is close to ˆ λ ≡ λf 4 / (3 M 2 ). In such cases, the ∆ N 2 grows extremely fast as ˆ ξ approaches ˆ λ from above (restoration of the PQ symmetry during the first stage requires ˆ ξ > ˆ λ ). Another important consequence of our scenario is the existence of the minimal duration ∆ N 2 of the second stage: for each given ˆ λ , ∆ N 2 as a function of ˆ ξ has a non-zero minimum proportional to the value of ˆ λ . In other words, by using the non-minimal coupling to restore PQ symmetry during the ϕ -driven inflation, we automatically get non-zero ∆ N 2 , viz meso-inflationary PQSB. For small ˆ λ , the minimal ∆ N 2 is small, but for ˆ λ ≳ 1, it grows proportionally.  \nThe parameter ˆ λ measures the ratio of the PQ potential to the inflaton potential at the slow-roll plateau, and can be used to outline three distinct parameter regimes: case (a) where ˆ λ ∼ O (1) and both potentials are close in magnitude; case (b) where ˆ λ ≪ 1 and the inflaton potential dominates; case (c) where ˆ λ ≫ 1 and the PQ potential dominates. In cases (a) and (c), the axion decay constant f is constrained from below by the observed amplitude of CMB scalar fluctuations, and in combination with experimental constraints, we have 10 -3 ≲ f ≲ 10 -1 in Planck units, i.e., a GUT-scale decay constant. In case (b), the decay constant can in principle be arbitrarily small. In transition between the two stages of inflation, slow-roll can be badly violated in case (b) due to large oscillations of both scalar fields but mildly in cases (a) and (c) (e.g., ϵ H can be small while η H ≥ 1) or absent.  \nIn all three cases, the brief transient stage of quantum diffusion (between stages I and II) plays a vital role as it sets the initial conditions for the second stage of our model (see e.g. Table 1). During the transient stage, the classical force from the potential is small enough that the quantum backreaction of the short-wavelength quantum modes becomes important for background evolution. Physically, quantum uncertainty will select out a specific set of ρ and ∂ N ρ values to break the symmetry, via wavefunction collapse. The entire evolution between stage I and stage II concerns a classical-to-quantum and then a quantum-to-classical process. A first-principle description of such evolution is challenging as the physics behind a wavefunction collapse remain elusive. Instead, we have adopted a strategy where we apply the stochastic formalism and approximate the collapsed field value with its statistical variance, evaluated when the stochastic noise is much smaller than the classical force. For the present work, this approximation suffices. But for more generic applications, say, in the case where non-Gaussianities are relevant, one may wish to adopt a more generic quantum backreaction framework, such as in [48,49].  \nBy studying scalar perturbations in this model, we find that the non-trivial field trajectories at the start of the second stage can lead to superhorizon changes of the power spectrum, including the CMB scales. In particular, the modes that exit the horizon around the effective critical point (when the ρ -mass vanishes) is amplified the most, which can be interesting from the perspective of PBHs and GW. We find that for small ˆ λ (which translates into sub-GUTscale decay constant), the superhorizon growth of the perturbations can become uncontrollably large, which can be used to exclude this parameter region by PBH overproduction. At the same time, the superhorizon evolution in this case becomes more sensitive to the field oscillation  \npatterns (at the start of the second stage), which depend on the shape of the scalar potential as well as friction terms caused by possible particle production, which can reduce the number of oscillations (we also find that for small ˆ λ , the Mathematica computation times can become extremely lengthy, and other numerical problems can arise). This sensitivity deserves more research, including more sophisticated numerical methods. In contrast, for ˆ λ ≳ 1, the field oscillations are mild, while the perturbations are relatively under control and less sensitive to the shape of the scalar potential (namely the inflaton potential). In this case, we find the amplification of the power spectrum at small scales up to P ζ ∼ 10 -2 (for ˆ ξ ∼ O (10)), yet with no significant PBH production. Nonetheless, these power spectra (Fig. 10) can generate large second-order GW, as shown in Fig. 11, potentially reaching the sensitivities of the upcoming spaceborne GW experiments.  \nLet us also comment on the values of the non-minimal coupling ξ . The duration of the second stage is roughly set by the value of ˆ ξ (aside from the extreme case ˆ ξ ≈ ˆ λ ). Hence, for meso-inflationary PQSB, where ∆ N 2 is between O (1) and O (10), we need ˆ ξ also in that range. Then, the magnitude of the original non-minimal coupling ξ is | ξ | = c/f 2 , where c sits between O (1) and O (10). As was shown, when ˆ λ ≳ 1, the CMB normalization leads to the GUT-scale decay constant f of around O (10 -3 ) ∼ O (10 -1 ), in which case | ξ | = O (10 2 ) ∼ O (10 7 ). For ˆ λ < 1, when f can be much smaller, the magnitude of ξ becomes even larger and may be problematic in view of effective field theory, see, e.g. [50-53]. This is an additional argument favoring large ˆ λ models, which bear GUT-scale decay constant. On the other hand, the nonminimal coupling of the PQ field may help solve the axion quality problem as discussed in [54], for the values | ξ | ≳ 2 × 10 3 , which are compatible with our scenario of meso-inflationary PQSB, even for large ˆ λ .  \nOur models can be considered as hybrid inflation with long waterfall regime. Nevertheless, in contrast to the conventional hybrid inflation with long waterfall [27,28,31-34,55], we focused on large field inflation (for the first stage), while the long waterfall regime in the PQ direction is realized by the flattening effect of the negative non-minimal coupling of the PQ field (without the need for super-Planckian axion decay constant), which at the same time restores U (1) PQ during the first stage. One can also see some similarities (including the methods used) with the model of [37], where the authors studied a scalar field coupled to R + R 2 -gravity, resulting in effectively a two-field inflation with a transient quantum diffusion phase.  \nEven though we used the Starobinsky model for the inflaton sector as an example, one can easily generalize it to any single-field inflationary model of interest, while the PQ potential can be considered as a general symmetry-breaking potential, not necessarily relevant for the strong CP problem of QCD.', 'Acknowledgements': 'DD is supported by China Postdoctoral Science Foundation (Certificate No.2023M730704). YW is supported by the General Program of Science and Technology of Shanghai No. 21ZR1406700, and Shanghai Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01). YW is grateful for the hospitality of the Perimeter Institute during his visit, where the main part of this work is done. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science.', 'A Canonical stochastic formalism': "Here we formulate the canonical (or phase-space) version standard stochastic formalism. The reason to repeat the exercise is mainly two-fold.  \nFirst, in the context of multi-field inflation with a brief quantum diffusion phase in between the prolonged classical evolution, one often uses the Lagrangian and slow-roll version of the Langevin equation. In this case, only ρ ' ( τ ) appears in the equations of motion (where τ is some yet unspecified parameterization of time). Schematically, when τ is chosen to be efold time, it may look something like  \nρ ' ( τ ) ≃ V ' ( ρ ) 3 H 2 + κξ . (95)  \nHere, ξ is a stochastic variable following a normalized Gaussian statistics and κ characterizes its (typically constant) amplitude.  \nIn this version, the field acceleration and its associated noise are often assumed to play a weaker role so as to be disregarded in the equations of motion. This might not be applicable, say, in cases where there is a large deceleration of ρ ' ( τ ) caused by Hubble friction, like in USR. The canonical momentum term and its noise cannot be neglected a priori .  \nSecond, a rigorous formulation of the stochastic formalism should start from canonical quantization of phase-space variables and then use the Wigner formalism to obtain an approximate phase-space distribution [30, 56] to represent stochastic uncertainty. While one could skip all this and heuristically add noise terms (with presumed statistics) to the equations of motion, one would not be able to know when or how to generalize these equations when they fail to approximate the true quantum evolution, especially when non-Gaussianites are not negligible. (We will not consider non-Gaussian aspects in this work. But the concern stands in general.)", 'A.1 Classical dynamics': "We will now recall key elements of the (canonical version) of the stochastic formalism. For our purpose, we start with the classical dynamics of a two-field model, with fields ρ and ϕ , as in the main text. But we will assume quantum diffusion is only non-negligible in the ρ -direction of the field space.  \nAs in the main text, one considers the action (in Einstein frame and M P = √ 1 / (8 πG N ) = 1)  \nS [ ρ, ϕ, g µν ] = ∫ dτd 3 x √ -g [ R 2 -1 2 G ρρ ∂ρ∂ρ -1 2 G ϕϕ ∂ϕ∂ϕ -V 1 ( ϕ ) -V 2 ( ρ ) ] . (96)  \nThe metric is assumed to be in an FLRW form with unspecified time parameterization τ  \nds 2 = g µν dx µ dx ν = -N 2 ( τ ) dτ 2 + a 2 ( τ ) δ ij dx i dx j . (97)  \nThe lapse function N ( τ ) is commonly chosen to be N ( τ ) = 1 , a ( τ ), and H -1 for cosmic time, conformal time, and efold number respectively. (In the appendix, the function N ( τ ) will be reserved for the lapse function and we will never commit to a specific time parameterization for τ , so that the efold number - also denoted as N in the main text - will never appear here.)  \nThe Ricci scalar and the square root of the metric determinant give  \nR = 6 ( a '' N ( τ ) 2 a + ( a ' ) 2 N ( τ ) 2 a 2 -a ' N ' ( τ ) aN ( τ ) 3 ) , √ -g = N ( τ ) a 3 , (98)  \nwhere a prime denotes a derivative with respect to the argument of the function. The purely gravitational part of the action can thus be simplified with integration-by-parts  \nS grav ≡ 1 2 ∫ dτd 3 x √ -gR = -3 ∫ dτd 3 x a ( a ' ) 2 N ( τ ) . (99)  \nThe canonical conjugates are thus  \nΠ a = -6 a N ( τ ) a ' , Π ρ = a 3 G ρρ N ( τ ) ρ ' , Π ϕ = a 3 G ϕϕ N ( τ ) ϕ ' , Π N = 0 . (100)  \nTo obtain the Hamiltonian, first recall that the lack of N ' ( τ ) in the Lagrangian gives rise to the primary constraint Π N = 0. This requires us to introduce a Lagrangian multiplier λ , giving us the Hamiltonian  \nH = ∫ d 3 x ( H g + H m ) H g = -N ( τ ) 12 a Π 2 a + λ Π N H m = 1 2 ( N ( τ ) a 3 G ρρ Π 2 ρ + N ( τ ) aG ρρ ∂ i ρ∂ j ρδ ij + N ( τ ) a 3 G ϕϕ Π 2 ϕ + N ( τ ) aG ϕϕ ∂ i ϕ∂ j ϕδ ij + N ( τ ) a 3 ( V 1 + V 2 ) ︸ ︷︷ ︸ N ( τ ) U ) ,  \nwhere we have defined the function U ≡ U ( ρ, ϕ, Π ϕ ) to separate the ρ -kinetic terms from the rest of the matter Hamiltonian (density).  \nFinally, the classical dynamics are given by the Hamilton's equations. The gravity part reads  \nΠ a = -6 a N ( τ ) a ' , λ = N ' ( τ ) Π ' N = 1 12 a Π 2 a -∂ H m ∂N ( τ ) = 0 (from constraint Π N = 0) Π ' a = -N ( τ ) 12 a 2 Π 2 a -∂ H m ∂a .  \nThe N ' ( τ )-equation reflects the time-reparameterization freedom due to the choice of the function λ ( τ ), and the second line is the familiar ' H 2 -equation' of the Friedmann equations.  \nThe matter part reads  \nΠ ρ = a 3 G ρρ N ( τ ) ρ ' , Π ϕ = a 3 G ϕϕ N ( τ ) ϕ ' Π ' ρ = N ( τ ) 2 a 3 G 2 ρρ ∂ ρ G ρρ Π 2 ρ + N ( τ )( aG ρρ ∂ i ∂ j ρδ ij -a 2 ∂ ρ G ρρ ∂ i ρ∂ j ρδ ij -∂ ρ U ) Π ' ϕ = N ( τ ) aG ϕϕ ∂ i ∂ j ϕδ ij -N ( τ ) a 3 ∂ ϕ ( V 1 + V 2 ) . (101)", 'A.2 Backreaction: stochastic corrections': "The fields ρ and ϕ of the Hamiltonian above can be canonically quantized. Nonetheless, since the field couples to the geometry via the Friedmann equations, quantizing the full system concerns quantum gravity and is beyond our present reach. Instead, we split the system into large and small scales separated by a time-dependent cut-off scale r σ = 1 /k σ = ( σaH ) -1 with constant σ ≪ 1. We quantize the small-scale modes with wavenumbers above k σ and view their backreaction on the large-scale modes effectively as noise.  \nMore concretely, the splitting is achieved with the help of a window function W ( x ) such that W ( x ) ∼ 1 for x ≫ 1 and W ( x ) ∼ 0 when x ≪ 1:  \nρ = ¯ ρ ( τ, x ) + ρ s ( τ, x ) , Π ρ = ¯ Π( τ, x ) + Π s ( τ, x ) ρ s ( τ, x ) = ∫ d 3 k (2 π ) 3 / 2 W ( k/k σ )˜ ρ k ( τ ) e -i k · x , Π s ( τ, x ) = ∫ d 3 k (2 π ) 3 / 2 W ( k/k σ ) ˜ Π k ( τ ) e -i k · x . (102)  \nNote that we have dropped the ρ -labels because in the following we will only be concerned with the stochastic dynamics of the ρ -field.  \nThe Fourier functions in (102) is linear in the mode functions of standard perturbation theory, which are typically kept at linear order in the equations of motion. (We will not consider non-Gaussianities in this work.) Following such standard analysis, we will assume that ˜ ρ k and ˜ Π k with | k | ≫ k σ solve the linear-order version of the equation set (101) when expanded around the background ¯ ρ and ¯ Π.  \nIn particular, this means that if we rewrite the full equations of motion (101) for ρ -field into the form (suppressing ϕ - and a -DOFs)  \nρ ' = f ( ρ, Π) , Π ' ρ = g ( ρ, Π) , (103)  \nthen ˜ ρ k and ˜ Π k with | k | ≫ k σ solve  \n˜ ρ ' k = ∂ ρ f ( ρ, Π) | ρ =¯ ρ, Π= ¯ Π ˜ ρ k + ∂ Π f ( ρ, Π) | ρ =¯ ρ, Π= ¯ Π ˜ Π k ˜ Π ' k = ∂ ρ g ( ρ, Π) | ρ =¯ ρ, Π= ¯ Π ˜ ρ k + ∂ Π g ( ρ, Π) | ρ =¯ ρ, Π= ¯ Π ˜ Π k , (104)  \nwhere spatial derivatives are understood in Fourier space, namely ∂ ρ ( ∂ i ρ )˜ ρ k → ik i ∂ ρ ( ρ )˜ ρ k and so on.  \nTo obtain the equations of motion for ¯ ρ and ¯ Π up to linear order in ˜ ρ k and ¯ Π k (or equivalently linear in ρ s and Π s ), we substitute the splitting (102) into the equations of motion (103) and use properties (104). The linear expansion of the LHS and RHS of (103) will partially cancel, leaving behind  \n¯ ρ ' + ∫ d 3 k (2 π ) 3 / 2 ∂ τ W ( k/k σ )˜ ρ k e -i k · x = f (¯ ρ, ¯ Π) ¯ Π ' + ∫ d 3 k (2 π ) 3 / 2 ∂ τ W ( k/k σ ) ˜ Π k e -i k · x = g (¯ ρ, ¯ Π) . (105)  \nA further simplification comes from the fact that one can ignore any spatial derivatives in f (¯ ρ, ¯ Π) and g (¯ ρ, ¯ Π). This is because σ ≪ 1 and the window function (more precisely 1 -W ( k/k σ )) ensures that ∂ i ¯ ρ and ∂ i ¯ Π are suppressed by σ when the derivatives are applied to the barred quantities.  \nOnce we apply canonical quantization to ˜ ρ k and ˜ Π k , they will follow a quantum distribution (wavefunction), which can be approximated by classical statistics in certain cases. A common  \nrealization is to use a Bunch-Davies solution for ˜ ρ k and ˜ Π k so that statistics only rely on two-point functions, which we can approximate with Gaussian noise terms ξ ρ and ξ Π for (105):  \n¯ ρ ' = f (¯ ρ, ¯ Π) + ξ ρ ¯ Π ' = g (¯ ρ, ¯ Π) + ξ Π ξ ρ = -∫ d 3 k (2 π ) 3 / 2 ∂ τ W ( k/k σ )˜ ρ k e -i k · x ξ Π = -∫ d 3 k (2 π ) 3 / 2 ∂ τ W ( k/k σ ) ˜ Π k e -i k · x , (106)  \nwhere ˜ ρ k and ˜ Π k are now understood as stochastic quantities.  \nFinally, the canonical version of the stochastic formalism makes it clear where the backreaction (or noise) comes from: Backreaction occurs when short-wavelength DOFs flow into the (long-wavelength) background as they get stretched via exponential expansion (viz, ξ ∼ ∂ τ W ( k/k σ )).", 'B Derivation of equations for the quadratic expectation values': 'In the previous appendix, we recalled the stochastic formalism using the field variable and its conjugate momentum. In the main text, to make the equations of motion simpler notationwise (as in (76)), we used instead the rescaled the momentum field (71). In this section, we derive the equations of motion for the quadratic expectation values of these rescaled variables.  \nWe start by linearizing the equations of motion of the ρ -field in (34). This is justified if one is interested in the stochastic evolution around the critical point, where ρ and π ρ are very small. Based on the result (106) of the previous section, the equations of motion of the ρ -field (with rescaled momentum π ρ ) can be written as (76), which one can put collectively into the matrix-notation equation  \n∂ τ ρ = A ( τ ) ρ + ξ ( x , τ ) , (107)  \nwhere in our case ρ = ( ρ, π ρ ) T and ξ ( x , τ ) = ( ξ ρ , ξ π ρ ) T . If τ is taken as the efold time, the matrix A reads  \nA = ( 0 1 -H -2 m 2 ˜ ρ, eff -(3 -ϵ H ) . ) (108)  \nThe solution to (107) with initial condition ρ 0 can be derived as [30]  \nρ ( τ ) = G ( τ, τ 0 ) ρ 0 + ∫ τ τ 0 ds G ( τ, s ) ξ ( s ) , (109)  \nwhere the Green matrix G ( τ, τ 0 ) can be defined as  \nG ( τ, τ 0 ) ≡ U ( τ ) U -1 ( τ 0 ) θ ( τ -τ 0 ) . (110)  \nThe matrix U in the definition of Green matrix can be constructed from two independent special solutions of classical EOM ∂ τ ρ = A ρ as  \nU ( τ ) = ( ρ (1) s ρ (2) s ) , (111)  \nin which each ρ s ( i ) represents one special solution of ∂ τ ρ = A ρ . We also note that based on the natural requirement ρ ( τ 0 ) = ρ 0 , we have G ( τ 0 , τ 0 ) = I and therefore θ (0) = 1 for the step-function. This is in contrast to other common choices such as θ (0) = 1 / 2 or 0. (In our context, the choice matters only for certain integrals involving delta-functions, which will be demonstrated below.)  \nFor Gaussian noise ⟨ 0 | ξ | 0 ⟩ = 0, the vacuum expectation value of ρ is exactly the one driven by a fully classical process. Therefore, to demonstrate the quantum noise effect, it is necessary to consider the variance of ρ defined as  \n∆ ≡ ⟨ 0 | ρρ † | 0 ⟩ - ⟨ 0 | ρ | 0 ⟩ ⟨ 0 | ρ † | 0 ⟩ = ∫ τ τ 0 ds G ( τ, s ) Z ( s ) G † ( τ, s ) . (112)  \nFor ρ = ( ρ 1 , ρ 2 ) T ≡ ( ρ, π ρ ) T , Gaussian noise implies ⟨ 0 | ξ ρ i ( x , τ 1 ) ξ ρ j ( x , τ 2 ) | 0 ⟩ ∝ δ ( τ 1 -τ 2 ) [30]. The matrix elements of matrix Z in the variance of ρ are defined through  \n⟨ 0 | ξ ρ i ( x , τ 1 ) ξ ρ j ( x , τ 2 ) | 0 ⟩ ≡ Z ρ i ρ j ( τ 1 ) δ ( τ 1 -τ 2 ) , (113)  \nwhich can be found once the mode functions are known.  \nTo derive the evolution of ∆ , we recall the definition of U and G . One finds that these matrices satisfy d U dτ = A ( τ ) U ( τ ) and ∂ G ( τ, τ 0 ) ∂τ = A ( τ ) G ( τ, τ 0 ) + I δ ( τ -τ 0 ). Therefore, taking the time-derivative of (112) one obtains  \n∂ τ ∆ = ∂ τ (∫ τ τ 0 ds G ( τ, s ) Z ( s ) G † ( τ, s ) ) = ∫ τ τ 0 ds ( ∂ τ G ( τ, s ) Z ( s ) G † ( τ, s ) + G ( τ, s ) Z ( s ) ∂ τ G † ( τ, s ) ) + ∫ τ + dτ τ ds G ( τ, s ) Z ( s ) G † ( τ, s ) /dτ = ∫ τ τ 0 ds ( A ( τ ) G ( τ, s ) Z ( s ) G † ( τ, s ) + G ( τ, s ) Z ( s ) G † ( τ, s ) A † ( τ ) ) + ∫ τ τ 0 δ ( τ -s )( Z ( s ) G † ( τ, s ) + G ( τ, s ) Z ( s )) + G ( τ, τ ) Z ( τ ) G † ( τ, τ ) = A ( τ ) ∆ ( τ ) + ∆ ( τ ) A † ( τ ) + Z ( τ ) . (114)  \nTo get the last equality, we used the identity associated with choosing θ (0) = 1  \n∫ τ τ 0 dsf ( s ) δ ( τ -s ) = ∫ τ τ 0 dsf ( s )( -∂ s θ ( τ -s )) = ∫ τ τ 0 ds∂ s f ( s ) θ ( τ -s ) -[ f ( τ ) θ ( τ -τ ) -f ( τ 0 ) θ ( τ -τ 0 )] = f ( τ ) -f ( τ 0 ) -[ f ( τ ) -f ( τ 0 )] = 0 .  \nUsing definition (112), the evolution equation (114) can be equivalently expressed as  \n∂ τ 〈 ΦΦ † 〉 = A ( τ ) 〈 ΦΦ † 〉 + 〈 ΦΦ † 〉 A † ( τ ) + Z ( τ ) . (115)  \nBy using our two-field equations of motion (76) along with (115), we obtain equations of motion (79) for the quadratic expectation values.', 'References': "- [1] J. Preskill, M. B. Wise, and F. Wilczek, 'Cosmology of the Invisible Axion,' Phys. Lett. B 120 (1983) 127-132.\n- [2] L. F. Abbott and P. Sikivie, 'A Cosmological Bound on the Invisible Axion,' Phys. Lett. B 120 (1983) 133-136.\n- [3] M. Dine and W. Fischler, 'The Not So Harmless Axion,' Phys. Lett. B 120 (1983) 137-141.\n- [4] A. S. Sakharov and M. Y. Khlopov, 'The Nonhomogeneity problem for the primordial axion field,' Phys. Atom. Nucl. 57 (1994) 485-487.\n- [5] A. S. Sakharov, D. D. Sokoloff, and M. Y. Khlopov, 'Large scale modulation of the distribution of coherent oscillations of a primordial axion field in the universe,' Phys. Atom. Nucl. 59 (1996) 1005-1010.\n- [6] M. Y. Khlopov, A. S. Sakharov, and D. D. Sokoloff, 'The nonlinear modulation of the density distribution in standard axionic CDM and its cosmological impact,' Nucl. Phys. B Proc. Suppl. 72 (1999) 105-109.\n- [7] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, 'String Axiverse,' Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th] .\n- [8] F. Chadha-Day, J. Ellis, and D. J. E. Marsh, 'Axion dark matter: What is it and why now?,' Sci. Adv. 8 no. 8, (2022) abj3618, arXiv:2105.01406 [hep-ph] .\n- [9] G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, 'The QCD axion, precisely,' JHEP 01 (2016) 034, arXiv:1511.02867 [hep-ph] .\n- [10] A. D. Linde, 'Axions in inflationary cosmology,' Phys. Lett. B 259 (1991) 38-47.\n- [11] M. Redi and A. Tesi, 'Meso-inflationary QCD axion,' Phys. Rev. D 107 no. 9, (2023) 095032, arXiv:2211.06421 [hep-ph] .\n- [12] K. Harigaya and L.-T. Wang, 'More axions from diluted domain walls,' arXiv:2211.08289 [hep-ph] .\n- [13] S. Ge, 'Stochastic gravitational wave background: birth from axionic string-wall death,' arXiv:2307.08185 [gr-qc] .\n- [14] S. Ge, J. Guo, and J. Liu, 'A New Mechanism for Primordial Black Hole Formation from QCD Axion,' arXiv:2309.01739 [hep-ph] .\n- [15] F. Takahashi and W. Yin, 'Kilobyte Cosmic Birefringence from ALP Domain Walls,' JCAP 04 (2021) 007, arXiv:2012.11576 [hep-ph] .\n- [16] D. Gonzalez, N. Kitajima, F. Takahashi, and W. Yin, 'Stability of domain wall network with initial inflationary fluctuations and its implications for cosmic birefringence,' Phys. Lett. B 843 (2023) 137990, arXiv:2211.06849 [hep-ph] .  \n- [17] N. Kitajima, J. Lee, F. Takahashi, and W. Yin, 'Stability of domain walls with inflationary fluctuations under potential bias, and gravitational wave signatures,' arXiv:2311.14590 [hep-ph] .\n- [18] J. a. G. Rosa and L. B. Ventura, 'Spontaneous breaking of the Peccei-Quinn symmetry during warm inflation,' arXiv:2105.05771 [hep-ph] .\n- [19] A. D. Linde, 'Hybrid inflation,' Phys. Rev. D 49 (1994) 748-754, arXiv:astro-ph/9307002 .\n- [20] J. Garcia-Bellido, A. D. Linde, and D. Wands, 'Density perturbations and black hole formation in hybrid inflation,' Phys. Rev. D 54 (1996) 6040-6058, arXiv:astro-ph/9605094 .\n- [21] J. Garcia-Bellido and D. Wands, 'The Spectrum of curvature perturbations from hybrid inflation,' Phys. Rev. D 54 (1996) 7181-7185, arXiv:astro-ph/9606047 .\n- [22] A. A. Starobinsky, 'A New Type of Isotropic Cosmological Models Without Singularity,' Phys. Lett. B 91 (1980) 99-102.\n- [23] E. McDonough, A. H. Guth, and D. I. Kaiser, 'Nonminimal Couplings and the Forgotten Field of Axion Inflation,' arXiv:2010.04179 [hep-th] .\n- [24] D. L. Lorenzoni, D. I. Kaiser, and E. McDonough, 'Natural Inflation with Exponentially Small Tensor-To-Scalar Ratio,' arXiv:2405.13881 [astro-ph.CO] .\n- [25] Planck Collaboration, Y. Akrami et al. , 'Planck 2018 results. X. Constraints on inflation,' Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO] .\n- [26] Y. Aldabergenov and S. V. Ketov, 'Primordial Black Holes from Volkov-Akulov-Starobinsky Supergravity,' Fortsch. Phys. 71 no. 6-7, (2023) 2300039, arXiv:2301.12750 [hep-th] .\n- [27] Y. Tada and M. Yamada, 'Primordial black hole formation in hybrid inflation,' Phys. Rev. D 107 no. 12, (2023) 123539, arXiv:2304.01249 [astro-ph.CO] .\n- [28] Y. Tada and M. Yamada, 'Stochastic dynamics of multi-waterfall hybrid inflation and formation of primordial black holes,' JCAP 11 (2023) 089, arXiv:2306.07324 [astro-ph.CO] .\n- [29] A. A. Starobinsky, 'STOCHASTIC DE SITTER (INFLATIONARY) STAGE IN THE EARLY UNIVERSE,' Lect. Notes Phys. 246 (1986) 107-126.\n- [30] V. Vennin, Stochastic inflation and primordial black holes . PhD thesis, U. Paris-Saclay, 6, 2020. arXiv:2009.08715 [astro-ph.CO] .\n- [31] S. Clesse, 'Hybrid inflation along waterfall trajectories,' Phys. Rev. D 83 (2011) 063518, arXiv:1006.4522 [gr-qc] .\n- [32] S. Clesse and J. Garc'ıa-Bellido, 'Massive Primordial Black Holes from Hybrid Inflation as Dark Matter and the seeds of Galaxies,' Phys. Rev. D 92 no. 2, (2015) 023524, arXiv:1501.07565 [astro-ph.CO] ."}